From 163b6e3de5a7235e1152926880e5474ec44af8fb Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Wed, 11 Oct 2017 16:58:00 -0400 Subject: [PATCH] Remove Fortran code for surface distance and sense --- src/geometry.F90 | 67 ++--- src/surface_header.C | 2 +- src/surface_header.F90 | 560 +---------------------------------------- 3 files changed, 23 insertions(+), 606 deletions(-) diff --git a/src/geometry.F90 b/src/geometry.F90 index 8a449a284b..873582dec9 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -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 diff --git a/src/surface_header.C b/src/surface_header.C index 6992c2e2ac..cf200e5e46 100644 --- a/src/surface_header.C +++ b/src/surface_header.C @@ -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; } } diff --git a/src/surface_header.F90 b/src/surface_header.F90 index ed4ef86e09..b0da37c1e0 100644 --- a/src/surface_header.F90 +++ b/src/surface_header.F90 @@ -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)