crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-29 06:35:48 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Remove Fortran code for surface distance and sense

Browse source
This commit is contained in:
Sterling Harper 2017-10-11 16:58:00 -04:00
parent 72cb25e677
commit 163b6e3de5
3 changed files with 23 additions and 606 deletions
Download patch file Download diff file Expand all files Collapse all files

67
src/geometry.F90
View file

@ -15,34 +15,34 @@ module geometry
implicit none
interface
function surface_sense_c(surf_ind, xyz, uvw) &
pure function surface_sense_c(surf_ind, xyz, uvw) &
bind(C, name='surface_sense') result(sense)
use ISO_C_BINDING
implicit none
integer(C_INT), value :: surf_ind;
real(C_DOUBLE) :: xyz(3);
real(C_DOUBLE) :: uvw(3);
logical(C_BOOL) :: sense;
integer(C_INT), intent(in), value :: surf_ind;
real(C_DOUBLE), intent(in) :: xyz(3);
real(C_DOUBLE), intent(in) :: uvw(3);
logical(C_BOOL) :: sense;
end function surface_sense_c
function surface_distance_c(surf_ind, xyz, uvw, coincident) &
pure function surface_distance_c(surf_ind, xyz, uvw, coincident) &
bind(C, name='surface_distance') result(d)
use ISO_C_BINDING
implicit none
integer(C_INT), value :: surf_ind;
real(C_DOUBLE) :: xyz(3);
real(C_DOUBLE) :: uvw(3);
logical(C_BOOL), value :: coincident;
real(C_DOUBLE) :: d;
integer(C_INT), intent(in), value :: surf_ind;
real(C_DOUBLE), intent(in) :: xyz(3);
real(C_DOUBLE), intent(in) :: uvw(3);
logical(C_BOOL), intent(in), value :: coincident;
real(C_DOUBLE) :: d;
end function surface_distance_c
subroutine surface_normal_c(surf_ind, xyz, uvw) &
pure subroutine surface_normal_c(surf_ind, xyz, uvw) &
bind(C, name='surface_normal')
use ISO_C_BINDING
implicit none
integer(C_INT), value :: surf_ind;
real(C_DOUBLE) :: xyz(3);
real(C_DOUBLE) :: uvw(3);
integer(C_INT), intent(in), value :: surf_ind;
real(C_DOUBLE), intent(in) :: xyz(3);
real(C_DOUBLE), intent(out) :: uvw(3);
end subroutine surface_normal_c
end interface
@ -62,8 +62,7 @@ contains
! expression using a stack, similar to how a RPN calculator would work.
!===============================================================================
!pure function cell_contains(c, p) result(in_cell)
function cell_contains(c, p) result(in_cell)
pure function cell_contains(c, p) result(in_cell)
type(Cell), intent(in) :: c
type(Particle), intent(in) :: p
logical :: in_cell
@ -75,8 +74,7 @@ contains
end if
end function cell_contains
!pure function simple_cell_contains(c, p) result(in_cell)
function simple_cell_contains(c, p) result(in_cell)
pure function simple_cell_contains(c, p) result(in_cell)
type(Cell), intent(in) :: c
type(Particle), intent(in) :: p
logical :: in_cell
@ -84,8 +82,6 @@ contains
integer :: i
integer :: token
logical :: actual_sense ! sense of particle wrt surface
logical :: alt_sense
real(8) :: uvw1(3), uvw2(3)
in_cell = .true.
do i = 1, size(c%rpn)
@ -101,25 +97,8 @@ contains
in_cell = .false.
exit
else
actual_sense = surfaces(abs(token))%obj%sense(&
actual_sense = surface_sense_c(abs(token)-1, &
p%coord(p%n_coord)%xyz, p%coord(p%n_coord)%uvw)
alt_sense = surface_sense_c(abs(token)-1, &
p%coord(p%n_coord)%xyz, p%coord(p%n_coord)%uvw)
if (actual_sense .neqv. alt_sense) then
write(*, *) surfaces(abs(token)) % obj % id
write(*, *) p % coord (p % n_coord) % xyz
write(*, *) p % coord (p % n_coord) % uvw
write(*, *) actual_sense, alt_sense
call fatal_error("")
end if
uvw1 = surfaces(abs(token))%obj%normal(p%coord(p%n_coord)%xyz)
call surface_normal_c(abs(token)-1, p%coord(p%n_coord)%xyz, uvw2)
if (.not. all(uvw1 - uvw2 == ZERO)) then
write(*, *) "Surface normals don't match"
write(*, *) uvw1
write(*, *) uvw2
call fatal_error("")
end if
if (actual_sense .neqv. (token > 0)) then
in_cell = .false.
exit
@ -167,7 +146,7 @@ contains
elseif (-token == p%surface) then
stack(i_stack) = .false.
else
actual_sense = surfaces(abs(token))%obj%sense(&
actual_sense = surface_sense_c(abs(token)-1, &
p%coord(p%n_coord)%xyz, p%coord(p%n_coord)%uvw)
stack(i_stack) = (actual_sense .eqv. (token > 0))
end if
@ -539,7 +518,6 @@ contains
type(Cell), pointer :: c
class(Surface), pointer :: surf
class(Lattice), pointer :: lat
real(8) :: d_alt
! inialize distance to infinity (huge)
dist = INFINITY
@ -573,13 +551,8 @@ contains
if (index_surf >= OP_UNION) cycle
! Calculate distance to surface
surf => surfaces(index_surf) % obj
d = surf % distance(p % coord(j) % xyz, p % coord(j) % uvw, coincident)
d_alt = surface_distance_c(index_surf-1, p % coord(j) % xyz, &
d = surface_distance_c(index_surf-1, p % coord(j) % xyz, &
p % coord(j) % uvw, logical(coincident, kind=C_BOOL))
if (d /= d_alt) then
write(*, *) d, d_alt
end if
! Check if calculated distance is new minimum
if (d < d_surf) then

2
src/surface_header.C
View file

@ -33,7 +33,7 @@ const char* get_coeff_str(pugi::xml_node surf_node)
return surf_node.child_value("coeffs");
} else {
std::cout << "ERROR: Found a surface with no coefficients" << std::endl;
return NULL;
return nullptr;
}
}

560
src/surface_header.F90
View file

@ -21,10 +21,8 @@ module surface_header
integer :: i_periodic = NONE ! Index of corresponding periodic surface
character(len=104) :: name = "" ! User-defined name
contains
procedure :: sense
procedure :: reflect
procedure(surface_evaluate_), deferred :: evaluate
procedure(surface_distance_), deferred :: distance
procedure(surface_normal_), deferred :: normal
end type Surface
@ -36,15 +34,6 @@ module surface_header
real(8) :: f
end function surface_evaluate_
pure function surface_distance_(this, xyz, uvw, coincident) result(d)
import Surface
class(Surface), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8), intent(in) :: uvw(3)
logical, intent(in) :: coincident
real(8) :: d
end function surface_distance_
pure function surface_normal_(this, xyz) result(uvw)
import Surface
class(Surface), intent(in) :: this
@ -63,8 +52,8 @@ module surface_header
!===============================================================================
! All the derived types below are extensions of the abstract Surface type. They
! inherent the reflect() and sense() type-bound procedures and must implement
! evaluate(), distance(), and normal()
! inherent the reflect() type-bound procedure and must implement evaluate(),
! and normal()
!===============================================================================
type, extends(Surface) :: SurfaceXPlane
@ -72,7 +61,6 @@ module surface_header
real(8) :: x0
contains
procedure :: evaluate => x_plane_evaluate
procedure :: distance => x_plane_distance
procedure :: normal => x_plane_normal
end type SurfaceXPlane
@ -81,7 +69,6 @@ module surface_header
real(8) :: y0
contains
procedure :: evaluate => y_plane_evaluate
procedure :: distance => y_plane_distance
procedure :: normal => y_plane_normal
end type SurfaceYPlane
@ -90,7 +77,6 @@ module surface_header
real(8) :: z0
contains
procedure :: evaluate => z_plane_evaluate
procedure :: distance => z_plane_distance
procedure :: normal => z_plane_normal
end type SurfaceZPlane
@ -102,7 +88,6 @@ module surface_header
real(8) :: D
contains
procedure :: evaluate => plane_evaluate
procedure :: distance => plane_distance
procedure :: normal => plane_normal
end type SurfacePlane
@ -113,7 +98,6 @@ module surface_header
real(8) :: r
contains
procedure :: evaluate => x_cylinder_evaluate
procedure :: distance => x_cylinder_distance
procedure :: normal => x_cylinder_normal
end type SurfaceXCylinder
@ -124,7 +108,6 @@ module surface_header
real(8) :: r
contains
procedure :: evaluate => y_cylinder_evaluate
procedure :: distance => y_cylinder_distance
procedure :: normal => y_cylinder_normal
end type SurfaceYCylinder
@ -135,7 +118,6 @@ module surface_header
real(8) :: r
contains
procedure :: evaluate => z_cylinder_evaluate
procedure :: distance => z_cylinder_distance
procedure :: normal => z_cylinder_normal
end type SurfaceZCylinder
@ -147,7 +129,6 @@ module surface_header
real(8) :: r
contains
procedure :: evaluate => sphere_evaluate
procedure :: distance => sphere_distance
procedure :: normal => sphere_normal
end type SurfaceSphere
@ -159,7 +140,6 @@ module surface_header
real(8) :: r2
contains
procedure :: evaluate => x_cone_evaluate
procedure :: distance => x_cone_distance
procedure :: normal => x_cone_normal
end type SurfaceXCone
@ -171,7 +151,6 @@ module surface_header
real(8) :: r2
contains
procedure :: evaluate => y_cone_evaluate
procedure :: distance => y_cone_distance
procedure :: normal => y_cone_normal
end type SurfaceYCone
@ -183,7 +162,6 @@ module surface_header
real(8) :: r2
contains
procedure :: evaluate => z_cone_evaluate
procedure :: distance => z_cone_distance
procedure :: normal => z_cone_normal
end type SurfaceZCone
@ -192,7 +170,6 @@ module surface_header
real(8) :: A, B, C, D, E, F, G, H, J, K
contains
procedure :: evaluate => quadric_evaluate
procedure :: distance => quadric_distance
procedure :: normal => quadric_normal
end type SurfaceQuadric
@ -205,36 +182,6 @@ module surface_header
contains
!===============================================================================
! SENSE determines whether a point is on the 'positive' or 'negative' side of a
! surface. This routine is crucial for determining what cell a particular point
! is in. The positive side is indicated by a returned value of .true. and the
! negative side is indicated by a returned value of .false.
!===============================================================================
pure function sense(this, xyz, uvw) result(s)
class(Surface), intent(in) :: this ! surface
real(8), intent(in) :: xyz(3)
real(8), intent(in) :: uvw(3)
logical :: s ! sense of particle
real(8) :: f ! surface function evaluated at point
! Evaluate the surface equation at the particle's coordinates to determine
! which side the particle is on
f = this%evaluate(xyz)
! Check which side of surface the point is on
if (abs(f) < FP_COINCIDENT) then
! Particle may be coincident with this surface. To determine the sense, we
! look at the direction of the particle relative to the surface normal (by
! default in the positive direction) via their dot product.
s = (dot_product(uvw, this%normal(xyz)) > ZERO)
else
s = (f > ZERO)
end if
end function sense
!===============================================================================
! REFLECT determines the direction a particle will travel if it is specularly
! reflected from the surface at a given position and direction. The position is
@ -275,24 +222,6 @@ contains
f = xyz(1) - this%x0
end function x_plane_evaluate
pure function x_plane_distance(this, xyz, uvw, coincident) result(d)
class(SurfaceXPlane), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8), intent(in) :: uvw(3)
logical, intent(in) :: coincident
real(8) :: d
real(8) :: f
f = this%x0 - xyz(1)
if (coincident .or. abs(f) < FP_COINCIDENT .or. uvw(1) == ZERO) then
d = INFINITY
else
d = f/uvw(1)
if (d < ZERO) d = INFINITY
end if
end function x_plane_distance
pure function x_plane_normal(this, xyz) result(uvw)
class(SurfaceXPlane), intent(in) :: this
real(8), intent(in) :: xyz(3)
@ -313,24 +242,6 @@ contains
f = xyz(2) - this%y0
end function y_plane_evaluate
pure function y_plane_distance(this, xyz, uvw, coincident) result(d)
class(SurfaceYPlane), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8), intent(in) :: uvw(3)
logical, intent(in) :: coincident
real(8) :: d
real(8) :: f
f = this%y0 - xyz(2)
if (coincident .or. abs(f) < FP_COINCIDENT .or. uvw(2) == ZERO) then
d = INFINITY
else
d = f/uvw(2)
if (d < ZERO) d = INFINITY
end if
end function y_plane_distance
pure function y_plane_normal(this, xyz) result(uvw)
class(SurfaceYPlane), intent(in) :: this
real(8), intent(in) :: xyz(3)
@ -353,24 +264,6 @@ contains
end function z_plane_evaluate
pure function z_plane_distance(this, xyz, uvw, coincident) result(d)
class(SurfaceZPlane), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8), intent(in) :: uvw(3)
logical, intent(in) :: coincident
real(8) :: d
real(8) :: f
f = this%z0 - xyz(3)
if (coincident .or. abs(f) < FP_COINCIDENT .or. uvw(3) == ZERO) then
d = INFINITY
else
d = f/uvw(3)
if (d < ZERO) d = INFINITY
end if
end function z_plane_distance
pure function z_plane_normal(this, xyz) result(uvw)
class(SurfaceZPlane), intent(in) :: this
real(8), intent(in) :: xyz(3)
@ -391,26 +284,6 @@ contains
f = this%A*xyz(1) + this%B*xyz(2) + this%C*xyz(3) - this%D
end function plane_evaluate
pure function plane_distance(this, xyz, uvw, coincident) result(d)
class(SurfacePlane), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8), intent(in) :: uvw(3)
logical, intent(in) :: coincident
real(8) :: d
real(8) :: f
real(8) :: projection
f = this%A*xyz(1) + this%B*xyz(2) + this%C*xyz(3) - this%D
projection = this%A*uvw(1) + this%B*uvw(2) + this%C*uvw(3)
if (coincident .or. abs(f) < FP_COINCIDENT .or. projection == ZERO) then
d = INFINITY
else
d = -f/projection
if (d < ZERO) d = INFINITY
end if
end function plane_distance
pure function plane_normal(this, xyz) result(uvw)
class(SurfacePlane), intent(in) :: this
real(8), intent(in) :: xyz(3)
@ -435,60 +308,6 @@ contains
f = y*y + z*z - this%r*this%r
end function x_cylinder_evaluate
pure function x_cylinder_distance(this, xyz, uvw, coincident) result(d)
class(SurfaceXCylinder), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8), intent(in) :: uvw(3)
logical, intent(in) :: coincident
real(8) :: d
real(8) :: y, z, k, a, c, quad
a = ONE - uvw(1)*uvw(1) ! v^2 + w^2
if (a == ZERO) then
d = INFINITY
else
y = xyz(2) - this%y0
z = xyz(3) - this%z0
k = y*uvw(2) + z*uvw(3)
c = y*y + z*z - this%r*this%r
quad = k*k - a*c
if (quad < ZERO) then
! no intersection with cylinder
d = INFINITY
elseif (coincident .or. abs(c) < FP_COINCIDENT) then
! particle is on the cylinder, thus one distance is positive/negative
! and the other is zero. The sign of k determines if we are facing in or
! out
if (k >= ZERO) then
d = INFINITY
else
d = (-k + sqrt(quad))/a
end if
elseif (c < ZERO) then
! particle is inside the cylinder, thus one distance must be negative
! and one must be positive. The positive distance will be the one with
! negative sign on sqrt(quad)
d = (-k + sqrt(quad))/a
else
! particle is outside the cylinder, thus both distances are either
! positive or negative. If positive, the smaller distance is the one
! with positive sign on sqrt(quad)
d = (-k - sqrt(quad))/a
if (d < ZERO) d = INFINITY
end if
end if
end function x_cylinder_distance
pure function x_cylinder_normal(this, xyz) result(uvw)
class(SurfaceXCylinder), intent(in) :: this
real(8), intent(in) :: xyz(3)
@ -515,60 +334,6 @@ contains
f = x*x + z*z - this%r*this%r
end function y_cylinder_evaluate
pure function y_cylinder_distance(this, xyz, uvw, coincident) result(d)
class(SurfaceYCylinder), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8), intent(in) :: uvw(3)
logical, intent(in) :: coincident
real(8) :: d
real(8) :: x, z, k, a, c, quad
a = ONE - uvw(2)*uvw(2) ! u^2 + w^2
if (a == ZERO) then
d = INFINITY
else
x = xyz(1) - this%x0
z = xyz(3) - this%z0
k = x*uvw(1) + z*uvw(3)
c = x*x + z*z - this%r*this%r
quad = k*k - a*c
if (quad < ZERO) then
! no intersection with cylinder
d = INFINITY
elseif (coincident .or. abs(c) < FP_COINCIDENT) then
! particle is on the cylinder, thus one distance is positive/negative
! and the other is zero. The sign of k determines if we are facing in or
! out
if (k >= ZERO) then
d = INFINITY
else
d = (-k + sqrt(quad))/a
end if
elseif (c < ZERO) then
! particle is inside the cylinder, thus one distance must be negative
! and one must be positive. The positive distance will be the one with
! negative sign on sqrt(quad)
d = (-k + sqrt(quad))/a
else
! particle is outside the cylinder, thus both distances are either
! positive or negative. If positive, the smaller distance is the one
! with positive sign on sqrt(quad)
d = (-k - sqrt(quad))/a
if (d < ZERO) d = INFINITY
end if
end if
end function y_cylinder_distance
pure function y_cylinder_normal(this, xyz) result(uvw)
class(SurfaceYCylinder), intent(in) :: this
real(8), intent(in) :: xyz(3)
@ -595,60 +360,6 @@ contains
f = x*x + y*y - this%r*this%r
end function z_cylinder_evaluate
pure function z_cylinder_distance(this, xyz, uvw, coincident) result(d)
class(SurfaceZCylinder), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8), intent(in) :: uvw(3)
logical, intent(in) :: coincident
real(8) :: d
real(8) :: x, y, k, a, c, quad
a = ONE - uvw(3)*uvw(3) ! u^2 + v^2
if (a == ZERO) then
d = INFINITY
else
x = xyz(1) - this%x0
y = xyz(2) - this%y0
k = x*uvw(1) + y*uvw(2)
c = x*x + y*y - this%r*this%r
quad = k*k - a*c
if (quad < ZERO) then
! no intersection with cylinder
d = INFINITY
elseif (coincident .or. abs(c) < FP_COINCIDENT) then
! particle is on the cylinder, thus one distance is positive/negative
! and the other is zero. The sign of k determines if we are facing in or
! out
if (k >= ZERO) then
d = INFINITY
else
d = (-k + sqrt(quad))/a
end if
elseif (c < ZERO) then
! particle is inside the cylinder, thus one distance must be negative
! and one must be positive. The positive distance will be the one with
! negative sign on sqrt(quad)
d = (-k + sqrt(quad))/a
else
! particle is outside the cylinder, thus both distances are either
! positive or negative. If positive, the smaller distance is the one
! with positive sign on sqrt(quad)
d = (-k - sqrt(quad))/a
if (d < ZERO) d = INFINITY
end if
end if
end function z_cylinder_distance
pure function z_cylinder_normal(this, xyz) result(uvw)
class(SurfaceZCylinder), intent(in) :: this
real(8), intent(in) :: xyz(3)
@ -676,55 +387,6 @@ contains
f = x*x + y*y + z*z - this%r*this%r
end function sphere_evaluate
pure function sphere_distance(this, xyz, uvw, coincident) result(d)
class(SurfaceSphere), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8), intent(in) :: uvw(3)
logical, intent(in) :: coincident
real(8) :: d
real(8) :: x, y, z, k, c, quad
x = xyz(1) - this%x0
y = xyz(2) - this%y0
z = xyz(3) - this%z0
k = x*uvw(1) + y*uvw(2) + z*uvw(3)
c = x*x + y*y + z*z - this%r*this%r
quad = k*k - c
if (quad < ZERO) then
! no intersection with sphere
d = INFINITY
elseif (coincident .or. abs(c) < FP_COINCIDENT) then
! particle is on the sphere, thus one distance is positive/negative and
! the other is zero. The sign of k determines if we are facing in or out
if (k >= ZERO) then
d = INFINITY
else
d = -k + sqrt(quad)
end if
elseif (c < ZERO) then
! particle is inside the sphere, thus one distance must be negative and
! one must be positive. The positive distance will be the one with
! negative sign on sqrt(quad)
d = -k + sqrt(quad)
else
! particle is outside the sphere, thus both distances are either positive
! or negative. If positive, the smaller distance is the one with positive
! sign on sqrt(quad)
d = -k - sqrt(quad)
if (d < ZERO) d = INFINITY
end if
end function sphere_distance
pure function sphere_normal(this, xyz) result(uvw)
class(SurfaceSphere), intent(in) :: this
real(8), intent(in) :: xyz(3)
@ -750,59 +412,6 @@ contains
f = y*y + z*z - this%r2*x*x
end function x_cone_evaluate
pure function x_cone_distance(this, xyz, uvw, coincident) result(d)
class(SurfaceXCone), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8), intent(in) :: uvw(3)
logical, intent(in) :: coincident
real(8) :: d
real(8) :: x, y, z, k, a, b, c, quad
x = xyz(1) - this%x0
y = xyz(2) - this%y0
z = xyz(3) - this%z0
a = uvw(2)*uvw(2) + uvw(3)*uvw(3) - this%r2*uvw(1)*uvw(1)
k = y*uvw(2) + z*uvw(3) - this%r2*x*uvw(1)
c = y*y + z*z - this%r2*x*x
quad = k*k - a*c
if (quad < ZERO) then
! no intersection with cone
d = INFINITY
elseif (coincident .or. abs(c) < FP_COINCIDENT) then
! particle is on the cone, thus one distance is positive/negative and the
! other is zero. The sign of k determines which distance is zero and which
! is not.
if (k >= ZERO) then
d = (-k - sqrt(quad))/a
else
d = (-k + sqrt(quad))/a
end if
else
! calculate both solutions to the quadratic
quad = sqrt(quad)
d = (-k - quad)/a
b = (-k + quad)/a
! determine the smallest positive solution
if (d < ZERO) then
if (b > ZERO) then
d = b
end if
else
if (b > ZERO) d = min(d, b)
end if
end if
! If the distance was negative, set boundary distance to infinity
if (d <= ZERO) d = INFINITY
end function x_cone_distance
pure function x_cone_normal(this, xyz) result(uvw)
class(SurfaceXCone), intent(in) :: this
real(8), intent(in) :: xyz(3)
@ -830,59 +439,6 @@ contains
f = x*x + z*z - this%r2*y*y
end function y_cone_evaluate
pure function y_cone_distance(this, xyz, uvw, coincident) result(d)
class(SurfaceYCone), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8), intent(in) :: uvw(3)
logical, intent(in) :: coincident
real(8) :: d
real(8) :: x, y, z, k, a, b, c, quad
x = xyz(1) - this%x0
y = xyz(2) - this%y0
z = xyz(3) - this%z0
a = uvw(1)*uvw(1) + uvw(3)*uvw(3) - this%r2*uvw(2)*uvw(2)
k = x*uvw(1) + z*uvw(3) - this%r2*y*uvw(2)
c = x*x + z*z - this%r2*y*y
quad = k*k - a*c
if (quad < ZERO) then
! no intersection with cone
d = INFINITY
elseif (coincident .or. abs(c) < FP_COINCIDENT) then
! particle is on the cone, thus one distance is positive/negative and the
! other is zero. The sign of k determines which distance is zero and which
! is not.
if (k >= ZERO) then
d = (-k - sqrt(quad))/a
else
d = (-k + sqrt(quad))/a
end if
else
! calculate both solutions to the quadratic
quad = sqrt(quad)
d = (-k - quad)/a
b = (-k + quad)/a
! determine the smallest positive solution
if (d < ZERO) then
if (b > ZERO) then
d = b
end if
else
if (b > ZERO) d = min(d, b)
end if
end if
! If the distance was negative, set boundary distance to infinity
if (d <= ZERO) d = INFINITY
end function y_cone_distance
pure function y_cone_normal(this, xyz) result(uvw)
class(SurfaceYCone), intent(in) :: this
real(8), intent(in) :: xyz(3)
@ -910,59 +466,6 @@ contains
f = x*x + y*y - this%r2*z*z
end function z_cone_evaluate
pure function z_cone_distance(this, xyz, uvw, coincident) result(d)
class(SurfaceZCone), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8), intent(in) :: uvw(3)
logical, intent(in) :: coincident
real(8) :: d
real(8) :: x, y, z, k, a, b, c, quad
x = xyz(1) - this%x0
y = xyz(2) - this%y0
z = xyz(3) - this%z0
a = uvw(1)*uvw(1) + uvw(2)*uvw(2) - this%r2*uvw(3)*uvw(3)
k = x*uvw(1) + y*uvw(2) - this%r2*z*uvw(3)
c = x*x + y*y - this%r2*z*z
quad = k*k - a*c
if (quad < ZERO) then
! no intersection with cone
d = INFINITY
elseif (coincident .or. abs(c) < FP_COINCIDENT) then
! particle is on the cone, thus one distance is positive/negative and the
! other is zero. The sign of k determines which distance is zero and which
! is not.
if (k >= ZERO) then
d = (-k - sqrt(quad))/a
else
d = (-k + sqrt(quad))/a
end if
else
! calculate both solutions to the quadratic
quad = sqrt(quad)
d = (-k - quad)/a
b = (-k + quad)/a
! determine the smallest positive solution
if (d < ZERO) then
if (b > ZERO) then
d = b
end if
else
if (b > ZERO) d = min(d, b)
end if
end if
! If the distance was negative, set boundary distance to infinity
if (d <= ZERO) d = INFINITY
end function z_cone_distance
pure function z_cone_normal(this, xyz) result(uvw)
class(SurfaceZCone), intent(in) :: this
real(8), intent(in) :: xyz(3)
@ -989,65 +492,6 @@ contains
end associate
end function quadric_evaluate
pure function quadric_distance(this, xyz, uvw, coincident) result(d)
class(SurfaceQuadric), intent(in) :: this
real(8), intent(in) :: xyz(3)
real(8), intent(in) :: uvw(3)
logical, intent(in) :: coincident
real(8) :: d
real(8) :: a, k, c
real(8) :: quad, b
associate (x => xyz(1), y => xyz(2), z => xyz(3), &
u => uvw(1), v => uvw(2), w => uvw(3))
a = this%A*u*u + this%B*v*v + this%C*w*w + this%D*u*v + this%E*v*w + &
this%F*u*w
k = (this%A*u*x + this%B*v*y + this%C*w*z + HALF*(this%D*(u*y + v*x) + &
this%E*(v*z + w*y) + this%F*(w*x + u*z) + this%G*u + this%H*v + &
this%J*w))
c = this%A*x*x + this%B*y*y + this%C*z*z + this%D*x*y + this%E*y*z + &
this%F*x*z + this%G*x + this%H*y + this%J*z + this%K
quad = k*k - a*c
if (quad < ZERO) then
! no intersection with surface
d = INFINITY
elseif (coincident .or. abs(c) < FP_COINCIDENT) then
! particle is on the surface, thus one distance is positive/negative and
! the other is zero. The sign of k determines which distance is zero and
! which is not.
if (k >= ZERO) then
d = (-k - sqrt(quad))/a
else
d = (-k + sqrt(quad))/a
end if
else
! calculate both solutions to the quadratic
quad = sqrt(quad)
d = (-k - quad)/a
b = (-k + quad)/a
! determine the smallest positive solution
if (d < ZERO) then
if (b > ZERO) then
d = b
end if
else
if (b > ZERO) d = min(d, b)
end if
end if
! If the distance was negative, set boundary distance to infinity
if (d <= ZERO) d = INFINITY
end associate
end function quadric_distance
pure function quadric_normal(this, xyz) result(uvw)
class(SurfaceQuadric), intent(in) :: this
real(8), intent(in) :: xyz(3)