diff --git a/src/geometry.F90 b/src/geometry.F90 index 8639ede581..62afc3fb3b 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -30,6 +30,7 @@ contains integer :: i_surface ! index in surfaces array (with sign) logical :: specified_sense ! specified sense of surface in list logical :: actual_sense ! sense of particle wrt surface + class(Surface), pointer :: s SURFACE_LOOP: do i = 1, c % n_surfaces ! Lookup surface @@ -47,7 +48,8 @@ contains ! Determine the specified sense of the surface in the cell and the actual ! sense of the particle with respect to the surface - actual_sense = sense(p, surfaces(abs(i_surface))%obj) + s => surfaces(abs(i_surface))%obj + actual_sense = s%sense(p%coord(p%n_coord)%xyz, p%coord(p%n_coord)%uvw) specified_sense = (c % surfaces(i) > 0) ! Compare sense of point to specified sense @@ -794,40 +796,6 @@ contains end subroutine distance_to_boundary -!=============================================================================== -! 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. -!=============================================================================== - - recursive function sense(p, surf) result(s) - type(Particle), intent(inout) :: p - class(Surface), intent(in) :: surf ! surface - logical :: s ! sense of particle - - integer :: j - real(8) :: func ! surface function evaluated at point - - j = p%n_coord - - ! Evaluate the surface equation at the particle's coordinates to determine - ! which side the particle is on - func = surf%evaluate(p%coord(j)%xyz) - - ! Check which side of surface the point is on - if (abs(func) < FP_COINCIDENT) then - ! Particle may be coincident with this surface. Artifically move the - ! particle forward a tiny bit. - p%coord(j)%xyz = p%coord(j)%xyz + TINY_BIT * p%coord(j)%uvw - s = sense(p, surf) - elseif (func > 0) then - s = .true. - else - s = .false. - end if - - end function sense - !=============================================================================== ! NEIGHBOR_LISTS builds a list of neighboring cells to each surface to speed up ! searches when a cell boundary is crossed. diff --git a/src/surface_header.F90 b/src/surface_header.F90 index 43b6c58322..9ba0e0d01a 100644 --- a/src/surface_header.F90 +++ b/src/surface_header.F90 @@ -1,6 +1,6 @@ module surface_header - use constants, only: ONE, TWO, ZERO, INFINITY, FP_COINCIDENT + use constants, only: ONE, TWO, ZERO, INFINITY, FP_COINCIDENT, TINY_BIT implicit none @@ -17,6 +17,7 @@ module surface_header integer :: bc ! Boundary condition character(len=52) :: name = "" ! User-defined name contains + procedure :: sense procedure(iEvaluate), deferred :: evaluate procedure(iDistance), deferred :: distance procedure(iReflect), deferred :: reflect @@ -192,6 +193,38 @@ 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. +!=============================================================================== + + recursive function sense(this, xyz, uvw) result(s) + class(Surface), intent(in) :: this ! surface + real(8), intent(inout) :: 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. Artifically move the + ! particle forward a tiny bit. + xyz(:) = xyz + TINY_BIT * uvw + s = this%sense(xyz, uvw) + elseif (f > 0) then + s = .true. + else + s = .false. + end if + + end function sense + !=============================================================================== ! SurfaceXPlane Implementation !===============================================================================