Ray plot

Geometry/Cells/Tests/simpleCell_test.f90

@@ -129,5 +129,32 @@ subroutine test_distance()
     @assertEqual(idx_ref, idx)
 
   end subroutine test_distance
+
+  !!
+  !! Test getting a normal
+  !!
+@Test
+  subroutine test_normal()
+    real(defReal), dimension(3) :: r, u
+    integer(shortInt)           :: idx
+    real(defReal), dimension(3) :: n
+
+    ! Y-Plane normal
+    r = [0.3_defReal, 0.4_defReal, 0.0_defReal]
+    u = [-ONE, ZERO, ZERO]
+    idx = surfs % getIdx(99)
+
+    n = cell % getNormal(idx, r, u)
+    @assertEqual([0, 1, 0], n)
+
+    ! Sphere normal
+    r = [0.0_defReal, 0.0_defReal, -2.0_defReal]
+    u = [-ONE, ZERO, ZERO]
+    idx = surfs % getIdx(13)
+
+    n = cell % getNormal(idx, r, u)
+    @assertEqual([0, 0, -1], n)
+
+  end subroutine test_normal
 
 end module simpleCell_test

Geometry/Cells/Tests/unionCell_test.f90

@@ -165,5 +165,32 @@ subroutine test_distance()
     @assertEqual(idx_ref, idx)
 
   end subroutine test_distance
+
+  !!
+  !! Test getting a normal
+  !!
+@Test
+  subroutine test_normal()
+    real(defReal), dimension(3) :: r, u
+    integer(shortInt)           :: idx
+    real(defReal), dimension(3) :: n
+
+    ! Sphere normal
+    r = [0.0_defReal, 0.0_defReal, -2.0_defReal]
+    u = [-ONE, ZERO, ZERO]
+    idx = surfs % getIdx(13)
+
+    n = easyCell % getNormal(idx, r, u)
+    @assertEqual([0, 0, -1], n)
+
+    ! Y-Plane normal
+    r = [0.1_defReal, 5.1_defReal, 0.0_defReal]
+    u = [-ONE, ZERO, ZERO]
+    idx = surfs % getIdx(5)
+
+    n = complexCell % getNormal(idx, r, u)
+    @assertEqual([0, 1, 0], n)
+
+  end subroutine test_normal
 
 end module unionCell_test

Geometry/Cells/cell_inter.f90

@@ -48,12 +48,13 @@ module cell_inter
     private
     integer(shortInt) :: cellId = 0
   contains
-    procedure(init), deferred     :: init
-    procedure(inside), deferred   :: inside
-    procedure(distance), deferred :: distance
-    procedure, non_overridable    :: setId
-    procedure, non_overridable    :: id
-    procedure                     :: kill
+    procedure(init), deferred      :: init
+    procedure(inside), deferred    :: inside
+    procedure(distance), deferred  :: distance
+    procedure(getNormal), deferred :: getNormal
+    procedure, non_overridable     :: setId
+    procedure, non_overridable     :: id
+    procedure                      :: kill
   end type cell
 
   abstract interface
@@ -117,6 +118,31 @@ pure subroutine distance(self, d, surfIdx, r, u)
       real(defReal), dimension(3), intent(in) :: r
       real(defReal), dimension(3), intent(in) :: u
     end subroutine distance
+
+    !!
+    !! Return normal of provided surface at provided point
+    !!
+    !! Args:
+    !!   surfIdx [in] -> Index of the surface to be crossed. If the surface is not defined on
+    !!     the surface shelf its value should be -ve.
+    !!   r [in]        -> Position
+    !!   u [in]        -> normalised direction (norm2(u) = 1.0)
+    !!
+    !! Output:
+    !!   normal -> outward normal from the given surface
+    !!
+    !! Errors:
+    !!   Throw a fatal error if the surfIdx is not a part of the cell
+    !!
+    function getNormal(self, surfIdx, r, u) result(normal)
+      import :: cell, defReal, shortInt
+      class(cell), intent(in)                 :: self
+      integer(shortInt), intent(in)           :: surfIdx
+      real(defReal), dimension(3), intent(in) :: r
+      real(defReal), dimension(3), intent(in) :: u
+      real(defReal), dimension(3)             :: normal
+
+    end function getNormal
 
   end interface
 

Geometry/Cells/simpleCell_class.f90

@@ -39,6 +39,7 @@ module simpleCell_class
     procedure :: init
     procedure :: inside
     procedure :: distance
+    procedure :: getNormal
     procedure :: kill
   end type simpleCell
 
@@ -138,6 +139,35 @@ pure subroutine distance(self, d, surfIdx, r, u)
     end do
 
   end subroutine distance
+
+  !!
+  !! Return normal of surfIdx at given co-ordinate
+  !!
+  !! See cell_inter for details
+  !!
+  function getNormal(self, surfIdx, r, u) result(normal)
+    class(simpleCell), intent(in)           :: self
+    integer(shortInt), intent(in)           :: surfIdx
+    real(defReal), dimension(3), intent(in) :: r
+    real(defReal), dimension(3), intent(in) :: u
+    real(defReal), dimension(3)             :: normal
+    integer(shortInt)                       :: i
+    character(100), parameter :: Here = 'getNormal (simpleCell_class.f90)'
+
+    ! To avoid compiler warnings
+    normal = ZERO
+
+    do i = 1, size(self % surfaces)
+      if (abs(self % surfaces(i) % surfIdx) == surfIdx) then
+        normal = self % surfaces(i) % ptr % normal(r, u)
+        return
+      end if
+    end do
+
+    ! Matching surface wasn't found
+    call fatalError(Here,'Provided surfIdx is not a member of the cell: '//numToChar(surfIdx))
+
+  end function getNormal
 
   !!
   !! Return to uninitialised state

Geometry/Cells/unionCell_class.f90

@@ -70,6 +70,7 @@ module unionCell_class
     procedure :: insideComplex
     procedure :: shortCircuit
     procedure :: distance
+    procedure :: getNormal
     procedure :: kill
   end type unionCell
 
@@ -459,6 +460,35 @@ pure subroutine distance(self, d, surfIdx, r, u)
     end do
 
   end subroutine distance
+
+  !!
+  !! Return normal of surfIdx at given co-ordinate
+  !!
+  !! See cell_inter for details
+  !!
+  function getNormal(self, surfIdx, r, u) result(normal)
+    class(unionCell), intent(in)            :: self
+    integer(shortInt), intent(in)           :: surfIdx
+    real(defReal), dimension(3), intent(in) :: r
+    real(defReal), dimension(3), intent(in) :: u
+    real(defReal), dimension(3)             :: normal
+    integer(shortInt)                       :: i
+    character(100), parameter :: Here = 'getNormal (unionCell_class.f90)'
+
+    ! To avoid compiler warnings
+    normal = ZERO
+
+    do i = 1, size(self % surfaces)
+      if (abs(self % surfaces(i) % surfIdx) == surfIdx) then
+        normal = self % surfaces(i) % ptr % normal(r, u)
+        return
+      end if
+    end do
+
+    ! Matching surface wasn't found
+    call fatalError(Here,'Provided surfIdx is not a member of the cell: '//numToChar(surfIdx))
+
+  end function getNormal
 
   !!
   !! Return to uninitialised state

Geometry/Surfaces/QuadSurfaces/aPlane_class.f90

@@ -41,6 +41,7 @@ module aPlane_class
     procedure :: evaluate
     procedure :: distance
     procedure :: going
+    procedure :: normal
     procedure :: kill
   end type aPlane
 
@@ -119,7 +120,7 @@ subroutine init(self, dict)
   end subroutine init
 
   !!
-  !! Return axix-align bounding box for the surface
+  !! Return axis-aligned bounding box for the surface
   !!
   !! See surface_inter for details
   !!
@@ -189,7 +190,7 @@ end function distance
   !!   For parallel direction halfspace is asigned by the sign of `evaluate` result.
   !!
   pure function going(self, r, u) result(halfspace)
-    class(aPlane), intent(in)              :: self
+    class(aPlane), intent(in)               :: self
     real(defReal), dimension(3), intent(in) :: r
     real(defReal), dimension(3), intent(in) :: u
     logical(defBool)                        :: halfspace
@@ -199,12 +200,26 @@ pure function going(self, r, u) result(halfspace)
     halfspace = ua > ZERO
 
     ! Special case of parallel direction
-    ! Partilce stays in its current halfspace
+    ! Particle stays in its current halfspace
     if (ua == ZERO) then
       halfspace = (r(self % axis) - self % a0) >= ZERO
     end if
 
   end function going
+
+  !!
+  !! Provides the normal for the plane
+  !!
+  pure function normal(self, r, u) result(n)
+    class(aPlane), intent(in)               :: self
+    real(defReal), dimension(3), intent(in) :: r
+    real(defReal), dimension(3), intent(in) :: u
+    real(defReal), dimension(3)             :: n
+
+    n = ZERO
+    n(self % axis) = ONE
+
+  end function normal
 
   !!
   !! Return to uninitialised state

Geometry/Surfaces/QuadSurfaces/cone_class.f90

@@ -2,7 +2,7 @@ module cone_class
 
   use numPrecision
   use universalVariables, only : SURF_TOL, INF, X_AXIS, Y_AXIS, Z_AXIS
-  use genericProcedures,  only : fatalError, numToChar, dotProduct
+  use genericProcedures,  only : fatalError, numToChar
   use dictionary_class,   only : dictionary
   use quadSurface_inter,  only : quadSurface
   use surface_inter,      only : kill_super => kill
@@ -75,6 +75,7 @@ module cone_class
     procedure :: evaluate
     procedure :: distance
     procedure :: going
+    procedure :: normal
     procedure :: kill
 
     ! Local procedures
@@ -396,40 +397,55 @@ pure function going(self, r, u) result(halfspace)
     real(defReal), dimension(3), intent(in) :: r
     real(defReal), dimension(3), intent(in) :: u
     logical(defBool)                        :: halfspace
-    real(defReal), dimension(3)             :: diff, norm
+    real(defReal), dimension(3)             :: n
     real(defReal)                           :: proj
 
+    n = self % normal(r, u)
+
+    proj = dot_product(n,u)
+
+    ! Determine halfspace
+    halfspace = proj > ZERO
+
+    ! Parallel direction
+    ! Need to use position to determine halfspace
+    if (proj == ZERO) then
+      halfspace = self % evaluate(r) >= ZERO
+    end if
+
+  end function going
+
+  !!
+  !! Produces the normal vector
+  !!
+  pure function normal(self, r, u) result(n)
+    class(cone), intent(in)                 :: self
+    real(defReal), dimension(3), intent(in) :: r
+    real(defReal), dimension(3), intent(in) :: u
+    real(defReal), dimension(3)             :: n
+    real(defReal), dimension(3)             :: diff
+
     diff = r - self % vertex
 
     ! Check the location of the particle, i.e., base or cone surface, to calculate
     ! the normal
     if (abs(diff(self % axis) - self % hMin) < self % surfTol()) then
-      norm(self % axis)  = -ONE
-      norm(self % plane) = ZERO
+      n(self % axis)  = -ONE
+      n(self % plane) = ZERO
 
     elseif (abs(diff(self % axis) - self % hMax) < self % surfTol()) then
-      norm(self % axis)  = ONE
-      norm(self % plane) = ZERO
+      n(self % axis)  = ONE
+      n(self % plane) = ZERO
 
     else
-      norm(self % plane) = diff(self % plane)
-      norm(self % axis)  = - self % tanSquare * diff(self % axis)
+      n(self % plane) = diff(self % plane)
+      n(self % axis)  = - self % tanSquare * diff(self % axis)
 
     end if
 
-    norm = norm/norm2(norm)
-    proj = dot_product(norm,u)
-
-    ! Determine halfspace
-    halfspace = proj > ZERO
-
-    ! Parallel direction
-    ! Need to use position to determine halfspace
-    if (proj == ZERO) then
-      halfspace = self % evaluate(r) >= ZERO
-    end if
+    n = n / norm2(n)
 
-  end function going
+  end function normal
 
   !!
   !! Return to uninitialised state