diff --git a/Geometry/Cells/Tests/simpleCell_test.f90 b/Geometry/Cells/Tests/simpleCell_test.f90 index ca547b201..abd94a7c9 100644 --- a/Geometry/Cells/Tests/simpleCell_test.f90 +++ b/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 diff --git a/Geometry/Cells/Tests/unionCell_test.f90 b/Geometry/Cells/Tests/unionCell_test.f90 index f5d860215..7491fdbab 100644 --- a/Geometry/Cells/Tests/unionCell_test.f90 +++ b/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 diff --git a/Geometry/Cells/cell_inter.f90 b/Geometry/Cells/cell_inter.f90 index 9a51a2e2b..05f49dc30 100644 --- a/Geometry/Cells/cell_inter.f90 +++ b/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 diff --git a/Geometry/Cells/simpleCell_class.f90 b/Geometry/Cells/simpleCell_class.f90 index a20333c53..bc6648d01 100644 --- a/Geometry/Cells/simpleCell_class.f90 +++ b/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 diff --git a/Geometry/Cells/unionCell_class.f90 b/Geometry/Cells/unionCell_class.f90 index 0e3c70345..72b2bf633 100644 --- a/Geometry/Cells/unionCell_class.f90 +++ b/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 diff --git a/Geometry/Surfaces/QuadSurfaces/aPlane_class.f90 b/Geometry/Surfaces/QuadSurfaces/aPlane_class.f90 index 4ee31b0b3..490d93ffb 100644 --- a/Geometry/Surfaces/QuadSurfaces/aPlane_class.f90 +++ b/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 diff --git a/Geometry/Surfaces/QuadSurfaces/cone_class.f90 b/Geometry/Surfaces/QuadSurfaces/cone_class.f90 index 8286b8987..8811e4be6 100644 --- a/Geometry/Surfaces/QuadSurfaces/cone_class.f90 +++ b/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 diff --git a/Geometry/Surfaces/QuadSurfaces/cylinder_class.f90 b/Geometry/Surfaces/QuadSurfaces/cylinder_class.f90 index ed5017c5d..bc2acb47d 100644 --- a/Geometry/Surfaces/QuadSurfaces/cylinder_class.f90 +++ b/Geometry/Surfaces/QuadSurfaces/cylinder_class.f90 @@ -2,7 +2,7 @@ module cylinder_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 @@ -55,6 +55,7 @@ module cylinder_class procedure :: evaluate procedure :: distance procedure :: going + procedure :: normal procedure :: kill ! Local procedures @@ -275,7 +276,7 @@ end function distance !! See surface_inter for details !! pure function going(self, r, u) result(halfspace) - class(cylinder), intent(in) :: self + class(cylinder), intent(in) :: self real(defReal), dimension(3), intent(in) :: r real(defReal), dimension(3), intent(in) :: u logical(defBool) :: halfspace @@ -287,6 +288,26 @@ pure function going(self, r, u) result(halfspace) halfspace = dot_product(rp , up ) >= ZERO end function going + + !! + !! Return the normal corresponding to the cylinder rim + !! + pure function normal(self, r, u) result(n) + class(cylinder), intent(in) :: self + real(defReal), dimension(3), intent(in) :: r + real(defReal), dimension(3), intent(in) :: u + real(defReal), dimension(3) :: n + integer(shortInt), dimension(2) :: p + + n = ZERO + + ! Define for short-hand + p = self % plane + + n(p) = r(p) - self % origin(p) + n = n / norm2(n) + + end function normal !! !! Return to uninitialised state diff --git a/Geometry/Surfaces/QuadSurfaces/plane_class.f90 b/Geometry/Surfaces/QuadSurfaces/plane_class.f90 index dfb913868..2b3623e06 100644 --- a/Geometry/Surfaces/QuadSurfaces/plane_class.f90 +++ b/Geometry/Surfaces/QuadSurfaces/plane_class.f90 @@ -2,7 +2,7 @@ module plane_class use numPrecision use universalVariables, only : X_AXIS, Y_AXIS, Z_AXIS, INF - use genericProcedures, only : fatalError, dotProduct, numToChar + use genericProcedures, only : fatalError, numToChar use dictionary_class, only : dictionary use quadSurface_inter, only : quadSurface use surface_inter, only : kill_super => kill @@ -38,6 +38,7 @@ module plane_class procedure :: evaluate procedure :: distance procedure :: going + procedure :: normal procedure :: kill end type plane @@ -106,7 +107,7 @@ end subroutine init !! !! See surface_inter for details !! - !! Always returns infinate box (even when aligned with some axis) + !! Always returns infinite box (even when aligned with some axis) !! pure function boundingBox(self) result(aabb) class(plane), intent(in) :: self @@ -115,7 +116,6 @@ pure function boundingBox(self) result(aabb) aabb(1:3) = -INF aabb(4:6) = INF - end function boundingBox !! @@ -128,7 +128,7 @@ pure function evaluate(self, r) result(c) real(defReal), dimension(3), intent(in) :: r real(defReal) :: c - c = dotProduct(r, self % norm) - self % offset + c = dot_product(r, self % norm) - self % offset end function evaluate @@ -151,7 +151,7 @@ pure function distance(self, r, u) result(d) real(defReal) :: d real(defReal) :: k, c - k = dotProduct(u, self % norm) + k = dot_product(u, self % norm) c = self % evaluate(r) if ( k == ZERO .or. abs(c) < self % surfTol()) then ! Parallel or at the surface @@ -180,7 +180,7 @@ pure function going(self, r, u) result(halfspace) logical(defBool) :: halfspace real(defReal) :: proj - proj = dotProduct(u, self % norm) + proj = dot_product(u, self % norm) halfspace = proj > ZERO ! Special case of parallel direction @@ -190,6 +190,19 @@ pure function going(self, r, u) result(halfspace) end if end function going + + !! + !! Return the surface normal, already computed + !! + pure function normal(self, r, u) result(n) + class(plane), intent(in) :: self + real(defReal), dimension(3), intent(in) :: r + real(defReal), dimension(3), intent(in) :: u + real(defReal), dimension(3) :: n + + n = self % norm + + end function normal !! !! Return to uninitialised state diff --git a/Geometry/Surfaces/QuadSurfaces/sphere_class.f90 b/Geometry/Surfaces/QuadSurfaces/sphere_class.f90 index 853b6113b..cb519da03 100644 --- a/Geometry/Surfaces/QuadSurfaces/sphere_class.f90 +++ b/Geometry/Surfaces/QuadSurfaces/sphere_class.f90 @@ -2,7 +2,7 @@ module sphere_class use numPrecision use universalVariables, only : INF, SURF_TOL - use genericProcedures, only : fatalError, dotProduct, numToChar + use genericProcedures, only : fatalError, numToChar use dictionary_class, only : dictionary use surface_inter, only : kill_super => kill use quadSurface_inter, only : quadSurface @@ -48,6 +48,7 @@ module sphere_class procedure :: evaluate procedure :: distance procedure :: going + procedure :: normal procedure :: kill end type sphere @@ -135,7 +136,7 @@ pure function evaluate(self, r) result(c) diff = r - self % origin - c = dotProduct(diff, diff) - self % r_sq + c = dot_product(diff, diff) - self % r_sq end function evaluate @@ -158,7 +159,7 @@ pure function distance(self, r, u) result(d) ! Calculate quadratic components c = self % evaluate(r) - k = dotProduct(r - self % origin, u) + k = dot_product(r - self % origin, u) delta = k*k - c ! Technically delta/4 ! Calculate the distance @@ -193,9 +194,24 @@ pure function going(self, r, u) result(halfspace) real(defReal), dimension(3), intent(in) :: u logical(defBool) :: halfspace - halfspace = dotProduct(r - self % origin, u) >= ZERO + halfspace = dot_product(r - self % origin, u) >= ZERO end function going + + !! + !! Return the normal corresponding to the sphere surface + !! + pure function normal(self, r, u) result(n) + class(sphere), intent(in) :: self + real(defReal), dimension(3), intent(in) :: r + real(defReal), dimension(3), intent(in) :: u + real(defReal), dimension(3) :: n + + n = r - self % origin + n = n / norm2(n) + + end function normal + !! !! Return to uninitialised state diff --git a/Geometry/Surfaces/Tests/aPlane_test.f90 b/Geometry/Surfaces/Tests/aPlane_test.f90 index 2ba3a783e..d6a398d9c 100644 --- a/Geometry/Surfaces/Tests/aPlane_test.f90 +++ b/Geometry/Surfaces/Tests/aPlane_test.f90 @@ -283,5 +283,23 @@ subroutine testDistance(this) @assertEqual(INF, this % surf % distance(r, u2)) end subroutine testDistance + + !! + !! Test producing the normal vector + !! +@Test(cases=[1, 2, 3]) + subroutine testNormal(this) + class(test_aPlane), intent(inout) :: this + real(defReal), dimension(3) :: n + real(defReal), dimension(3) :: r, u + + r = [99.92_defReal, -6.0_defReal, 4.0_defReal] + u = [100, 200, 400] + + n = this % surf % normal(r, u) + + @assertEqual(ONE, n(this % axis)) + + end subroutine testNormal end module aPlane_test diff --git a/Geometry/Surfaces/Tests/box_test.f90 b/Geometry/Surfaces/Tests/box_test.f90 index 274c8993d..9ccb97594 100644 --- a/Geometry/Surfaces/Tests/box_test.f90 +++ b/Geometry/Surfaces/Tests/box_test.f90 @@ -258,6 +258,46 @@ subroutine testDistance() @assertEqual(ref, surf % distance(r, -u), TOL*ref) end subroutine testDistance + + !! + !! Test normal calculation + !! +@Test + subroutine testNormal() + real(defReal), dimension(3) :: r, u, n + real(defReal), parameter :: TOL = 1.0E-7 + + u = [0.3, 0.2, 5.0] + + ! Test on a plane + r = [TWO, TWO, ONE] + + n = surf % normal(r, u) + + @assertEqual(ONE, n(1), TOL) + @assertEqual(ZERO, n(2), TOL) + @assertEqual(ZERO, n(3), TOL) + + + ! Test on an edge + r = [ONE, 4.0_defReal, -2.0_defReal] + + n = surf % normal(r, u) + + @assertEqual(ZERO, n(1), TOL) + @assertEqual(ONE / sqrt(TWO), n(2), TOL) + @assertEqual(-ONE / sqrt(TWO), n(3), TOL) + + ! Test on a corner + r = [ZERO, ZERO, -TWO] + + n = surf % normal(r, u) + + @assertEqual(-ONE / sqrt(3.0_defReal), n(1), TOL) + @assertEqual(-ONE / sqrt(3.0_defReal), n(2), TOL) + @assertEqual(-ONE / sqrt(3.0_defReal), n(3), TOL) + + end subroutine testNormal !! !! Test Edge Cases @@ -273,7 +313,7 @@ subroutine testEdgeCases() ! ** Corner ! * Particle is almost at the corner - ! Either it is outside or can escape with a short movment in next step + ! Either it is outside or can escape with a short movement in next step ! ! Currently does not work for Y-position in plane (r(2) == ZERO) ! See `distance` documentation in box_class @@ -281,10 +321,11 @@ subroutine testEdgeCases() ! eps = 5.0_defReal * epsilon(eps) r = [2.0_defReal-eps, eps, 4.0_defReal-eps] - u = [HALF, ZERO, -ONE] + !u = [HALF, ZERO, -ONE] + u = [ONE/sqrt(3.0), ONE/sqrt(3.0), ONE/sqrt(3.0)] u = u/norm2(u) hs = surf % halfspace(r, u) - if (.not.hs) then ! Perform small movment + if (.not.hs) then ! Perform small movement d = surf % distance(r, u) @assertTrue( abs(d) < 1.0E-6) @assertTrue(surf % halfspace(r + d*u, u)) @@ -294,24 +335,24 @@ subroutine testEdgeCases() u = [-ONE, ZERO, HALF] u = u/norm2(u) hs = surf % halfspace(r, u) - if (.not.hs) then ! Perform small movment + if (.not.hs) then ! Perform small movement d = surf % distance(r, u) @assertTrue( abs(d) < 1.0E-6) @assertTrue(surf % halfspace(r + d*u, u)) end if - ! Try asymertic corner + ! Try asymmetric corner r = [2.0_defReal-TWO*eps, eps, 4.0_defReal-eps] u = [HALF, ZERO, -ONE] u = u/norm2(u) hs = surf % halfspace(r, u) - if (.not.hs) then ! Perform small movment + if (.not.hs) then ! Perform small movement d = surf % distance(r, u) @assertTrue( abs(d) < 1.0E-6) @assertTrue(surf % halfspace(r + d*u, u)) end if - ! Asymetric corner position + ! Asymmetric corner position ! Try other direction r = [2.0_defReal-eps, eps, 4.0_defReal-TWO*eps] u = [-ONE, ZERO, HALF] diff --git a/Geometry/Surfaces/Tests/cone_test.f90 b/Geometry/Surfaces/Tests/cone_test.f90 index 53d7b1078..c71dc6a17 100644 --- a/Geometry/Surfaces/Tests/cone_test.f90 +++ b/Geometry/Surfaces/Tests/cone_test.f90 @@ -70,6 +70,7 @@ end function toString !! Vertex 1.0, 1.0, 1.0 !! Tangent 1.0 !! ID 52 + !! Base 11.0 !! function newTestCase(dir) result(tst) type(dirParam), intent(in) :: dir @@ -477,5 +478,46 @@ subroutine testDistance(this) @assertEqual(ref, this % surf % distance(r, u), TOL * ref) end subroutine testDistance + + !! + !! Test normal calculation + !! +@Test(cases=[1, 2, 3]) + subroutine testNormal(this) + class(test_cone), intent(inout) :: this + integer(shortInt) :: a, p1, p2 + real(defReal), dimension(3) :: r, u, n + real(defReal), parameter :: TOL = 1.0E-7 + + ! Get axis and different planar directions + a = this % axis + p1 = this % plane(1) + p2 = this % plane(2) + + u = [0.3, 0.2, 5.0] + + ! Test on cylindrical surface/base + r(p1) = ONE + r(p2) = ONE + r(a) = ONE + 10.0_defReal + + n = this % surf % normal(r, u) + + @assertEqual(ZERO, n(p1), TOL) + @assertEqual(ZERO, n(p2), TOL) + @assertEqual(ONE, n(a), TOL) + + ! Test on the slope + r(p1) = ONE + ZERO + r(p2) = ONE + TWO + r(a) = ONE + TWO + + n = this % surf % normal(r, u) + + @assertEqual(ZERO, n(p1), TOL) + @assertEqual(ONE/sqrt(TWO), n(p2), TOL) + @assertEqual(-ONE/sqrt(TWO), n(a), TOL) + + end subroutine testNormal end module cone_test diff --git a/Geometry/Surfaces/Tests/cylinder_test.f90 b/Geometry/Surfaces/Tests/cylinder_test.f90 index ea1880241..bbae5de12 100644 --- a/Geometry/Surfaces/Tests/cylinder_test.f90 +++ b/Geometry/Surfaces/Tests/cylinder_test.f90 @@ -339,5 +339,34 @@ subroutine testDistance(this) @assertEqual(ref, this % surf % distance(r, -u), TOL * ref) end subroutine testDistance + + !! + !! Test producing the normal vector + !! +@Test(cases=[1, 2, 3]) + subroutine testNormal(this) + class(test_cylinder), intent(inout) :: this + real(defReal), dimension(3) :: n + real(defReal), dimension(3) :: r, u + integer(shortInt) :: a, p1, p2 + + ! Set axis and plane axis indices + a = this % axis + p1 = this % plane(1) + p2 = this % plane(2) + + r(a) = TWO + r(p1) = 3.0_defReal + r(p2) = 1.0_defReal + + u = [1, 0, 0] + + n = this % surf % normal(r, u) + + @assertEqual(n(a), ZERO) + @assertEqual(n(p1), ONE / sqrt(TWO)) + @assertEqual(n(p2), -ONE / sqrt(TWO)) + + end subroutine testNormal end module cylinder_test diff --git a/Geometry/Surfaces/Tests/plane_test.f90 b/Geometry/Surfaces/Tests/plane_test.f90 index 41e184467..f9fb0b043 100644 --- a/Geometry/Surfaces/Tests/plane_test.f90 +++ b/Geometry/Surfaces/Tests/plane_test.f90 @@ -4,7 +4,7 @@ module plane_test use universalVariables, only : SURF_TOL, VACUUM_BC, REFLECTIVE_BC, INF use dictionary_class, only : dictionary use dictParser_func, only : charToDict - use plane_class, only : plane + use plane_class, only : plane use funit implicit none @@ -179,5 +179,25 @@ subroutine testDistance() @assertEqual(INF, surf % distance(r, u2)) end subroutine testDistance + + !! + !! Test producing the normal vector + !! +@Test + subroutine testNormal() + real(defReal), dimension(3) :: n + real(defReal), dimension(3) :: r, u + real(defReal), parameter :: TOL = 1.0E-7 + + r = [99.92, -6.0, 4.0] + u = [100, 200, 400] + + n = surf % normal(r, u) + + @assertEqual(ONE/sqrt(3.0), n(1), TOL) + @assertEqual(ONE/sqrt(3.0), n(2), TOL) + @assertEqual(ONE/sqrt(3.0), n(3), TOL) + + end subroutine testNormal end module plane_test diff --git a/Geometry/Surfaces/Tests/sphere_test.f90 b/Geometry/Surfaces/Tests/sphere_test.f90 index 34576351e..98b04e58e 100644 --- a/Geometry/Surfaces/Tests/sphere_test.f90 +++ b/Geometry/Surfaces/Tests/sphere_test.f90 @@ -181,5 +181,30 @@ subroutine testDistance() @assertEqual(ref, surf % distance(r, -u)) end subroutine testDistance + + !! + !! Test producing the normal vector + !! +@Test + subroutine testNormal() + real(defReal), dimension(3) :: n + real(defReal), dimension(3) :: r, u + + r = [3.0_defReal, TWO, ONE] + u = [100, 200, 400] + n = surf % normal(r, u) + + @assertEqual(ONE, n(1)) + @assertEqual(ZERO, n(2)) + @assertEqual(ZERO, n(3)) + + r = [ONE, 3.0_defReal, ZERO] + n = surf % normal(r, u) + + @assertEqual(ZERO, n(1)) + @assertEqual(ONE/sqrt(TWO), n(2)) + @assertEqual(-ONE/sqrt(TWO), n(3)) + + end subroutine testNormal end module sphere_test diff --git a/Geometry/Surfaces/Tests/squareCylinder_test.f90 b/Geometry/Surfaces/Tests/squareCylinder_test.f90 index 42cb749cf..7ad1005f8 100644 --- a/Geometry/Surfaces/Tests/squareCylinder_test.f90 +++ b/Geometry/Surfaces/Tests/squareCylinder_test.f90 @@ -10,7 +10,6 @@ module squareCylinder_test !! !! Test parameter wrapper around AN INTEGER (bit of boilerplate) !! - !! @testParameter(constructor=newParam) type, extends (AbstractTestParameter) :: dirParam integer(shortInt) :: dir @@ -327,10 +326,6 @@ subroutine testHalfspace(this) @assertTrue(this % surf % halfspace(r + eps*u, u2)) @assertFalse(this % surf % halfspace(r - eps*u, u2)) - ! At the surface in diffrent quadrant - - - end subroutine testHalfspace !! @@ -416,6 +411,47 @@ subroutine testDistance(this) @assertEqual(ref, this % surf % distance(r, -u), ref * TOL) end subroutine testDistance + + !! + !! Test normal calculation + !! +@Test(cases=[1, 2, 3]) + subroutine testNormal(this) + class(test_squareCylinder), intent(inout) :: this + integer(shortInt) :: ax, p1, p2 + real(defReal), dimension(3) :: r, u, n + + ! Get axis and different planar directions + ax = this % axis + p1 = this % plane(1) + p2 = this % plane(2) + + u = [0.3, 0.2, 5.0] + + ! Test on a plane + r(p1) = ONE + TWO + r(p2) = TWO + r(ax) = TWO + + n = this % surf % normal(r, u) + + @assertEqual(n(p1), ONE) + @assertEqual(n(p2), ZERO) + @assertEqual(n(ax), ZERO) + + + ! Test on an edge + r(p1) = ONE - TWO + r(p2) = TWO + 3.0_defReal + r(ax) = TWO + + n = this % surf % normal(r, u) + + @assertEqual(n(p1), -ONE / sqrt(TWO)) + @assertEqual(n(p2), ONE / sqrt(TWO)) + @assertEqual(n(ax), ZERO) + + end subroutine testNormal !! !! Test Edge Cases @@ -431,7 +467,7 @@ subroutine testEdgeCases(this) real(defReal) :: eps, d logical(defBool) :: hs - ! Get axis and diffrent planar directions + ! Get axis and different planar directions ax = this % axis p1 = this % plane(1) p2 = this % plane(2) @@ -444,7 +480,7 @@ subroutine testEdgeCases(this) u([ax, p1, p2]) = [ZERO, TWO, -ONE] u = u/norm2(u) hs = this % surf % halfspace(r, u) - if (.not.hs) then ! Perform small movment + if (.not.hs) then ! Perform small movement d = this % surf % distance(r, u) @assertTrue( abs(d) < 1.0E-6) @assertTrue(this % surf % halfspace(r + d*u, u)) @@ -454,19 +490,19 @@ subroutine testEdgeCases(this) u([ax, p1, p2]) = [ZERO, -ONE, TWO] u = u/norm2(u) hs = this % surf % halfspace(r, u) - if (.not.hs) then ! Perform small movment + if (.not.hs) then ! Perform small movement d = this % surf % distance(r, u) @assertTrue( abs(d) < 1.0E-6) @assertTrue(this % surf % halfspace(r + d*u, u)) end if - ! Try asymertic corner + ! Try asymmetric corner ! Point is not exactly at the diagonal r([ax, p1, p2]) = [ZERO, -1.0_defReal+eps, -1.0_defReal+eps*TWO] u([ax, p1, p2]) = [ZERO, TWO, -ONE] u = u/norm2(u) hs = this % surf % halfspace(r, u) - if (.not.hs) then ! Perform small movment + if (.not.hs) then ! Perform small movement d = this % surf % distance(r, u) @assertTrue( abs(d) < 1.0E-6) @assertTrue(this % surf % halfspace(r + d*u, u)) @@ -477,7 +513,7 @@ subroutine testEdgeCases(this) u([ax, p1, p2]) = [ZERO, -ONE, TWO] u = u/norm2(u) hs = this % surf % halfspace(r, u) - if (.not.hs) then ! Perform small movment + if (.not.hs) then ! Perform small movement d = this % surf % distance(r, u) @assertTrue( abs(d) < 1.0E-6) @assertTrue(this % surf % halfspace(r + d*u, u)) diff --git a/Geometry/Surfaces/Tests/truncCylinder_test.f90 b/Geometry/Surfaces/Tests/truncCylinder_test.f90 index 1a7216f20..d1d42ea97 100644 --- a/Geometry/Surfaces/Tests/truncCylinder_test.f90 +++ b/Geometry/Surfaces/Tests/truncCylinder_test.f90 @@ -485,6 +485,78 @@ subroutine testDistance(this) @assertEqual(ref, this % surf % distance(r, u), ref * TOL) end subroutine testDistance + + !! + !! Test normal calculation + !! +@Test(cases=[1, 2, 3]) + subroutine testNormal(this) + class(test_truncCylinder), intent(inout) :: this + integer(shortInt) :: ax, p1, p2 + real(defReal), dimension(3) :: r, u, n + real(defReal), parameter :: TOL = 1.0E-7 + + ! Get axis and diffrent planar directions + ax = this % axis + p1 = this % plane(1) + p2 = this % plane(2) + + u = [0.3, 0.2, 5.0] + + ! Test on cylindrical surface + r(p1) = ONE + TWO + r(p2) = TWO + r(ax) = TWO + + n = this % surf % normal(r, u) + + @assertEqual(ONE, n(p1), TOL) + @assertEqual(ZERO, n(p2), TOL) + @assertEqual(ZERO, n(ax), TOL) + + r(p1) = ONE - sqrt(TWO) + r(p2) = TWO + sqrt(TWO) + r(ax) = TWO + + n = this % surf % normal(r, u) + + @assertEqual(-ONE / sqrt(TWO), n(p1), TOL) + @assertEqual(ONE / sqrt(TWO), n(p2), TOL) + @assertEqual(ZERO, n(ax), TOL) + + ! Test on a circle surface + r(p1) = ONE + ONE + r(p2) = TWO + HALF + r(ax) = TWO + 1.5_defReal + + n = this % surf % normal(r, u) + + @assertEqual(ZERO, n(p1), TOL) + @assertEqual(ZERO, n(p2), TOL) + @assertEqual(ONE, n(ax), TOL) + + r(p1) = ONE + 1.1_defReal + r(p2) = TWO - 0.3_defReal + r(ax) = TWO - 1.5_defReal + + n = this % surf % normal(r, u) + + @assertEqual(ZERO, n(p1), TOL) + @assertEqual(ZERO, n(p2), TOL) + @assertEqual(-ONE, n(ax), TOL) + + ! Test on the edge + r(p1) = ONE + TWO + r(p2) = TWO + r(ax) = HALF + + n = this % surf % normal(r, u) + + @assertEqual(ONE / sqrt(TWO), n(p1), TOL) + @assertEqual(ZERO, n(p2), TOL) + @assertEqual(-ONE / sqrt(TWO), n(ax), TOL) + + end subroutine testNormal !! !! Test Edge Cases diff --git a/Geometry/Surfaces/box_class.f90 b/Geometry/Surfaces/box_class.f90 index 221903c9b..f8e1eb995 100644 --- a/Geometry/Surfaces/box_class.f90 +++ b/Geometry/Surfaces/box_class.f90 @@ -9,6 +9,8 @@ module box_class implicit none private + real(defReal), parameter :: EDGE_TOL = 1.0E-6 + !! !! Axis Aligned Box !! @@ -48,6 +50,7 @@ module box_class procedure :: evaluate procedure :: distance procedure :: going + procedure :: normal procedure :: kill procedure :: setBC procedure :: explicitBC @@ -155,12 +158,12 @@ end function evaluate !! Otherwise either near or far will be the closest distance depending on whether !! the particle is in or outside the box !! - !! NOTE: The algorithim as implemented may have some problems in the corners. If particle + !! NOTE: The algorithm as implemented may have some problems in the corners. If particle !! is EXACTLY (to the bit) in the corner, or it goes through corner while being PERFECTLY !! (again to the bit) in one of the planes !! !! The reason for this is that a ray contained EXACTLY in one of the boundary planes - !! may be assigned with an empty section (instead of infinate length). + !! may be assigned with an empty section (instead of infinite length). !! pure function distance(self, r, u) result(d) class(box), intent(in) :: self @@ -241,13 +244,13 @@ end function distance !! Works by: !! 1) Determine a plane in which direction is the closest !! 2) Use normal for this plane and project distance - !! 3) Determinie halfspace based on sign of the projection + !! 3) Determine halfspace based on sign of the projection !! !! Note: - !! For parallel direction halfspace is asigned by the sign of `evaluate` result. + !! For parallel direction halfspace is assigned by the sign of `evaluate` result. !! pure function going(self, r, u) result(halfspace) - class(box), intent(in) :: self + class(box), intent(in) :: self real(defReal), dimension(3), intent(in) :: r real(defReal), dimension(3), intent(in) :: u logical(defBool) :: halfspace @@ -274,6 +277,47 @@ pure function going(self, r, u) result(halfspace) end function going + !! + !! Produces the normal for the box. Makes sure to check edges and corners. + !! Does so by checking each cardinal direction and comparing. + !! + pure function normal(self, r, u) result(n) + class(box), 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) :: rl, rAbs + real(defReal) :: maxDist + integer(shortInt) :: i + + rl = r - self % origin + rAbs = abs(rl) - self % halfwidth + + ! Check each dimension, with a tolerance to detect corners/edges + n = ZERO + maxDist = ZERO + + do i = 1, 3 + + ! Intersecting an edge/corner + if (abs(rAbs(i) - maxDist) < EDGE_TOL) then + n(i) = sign(ONE, rl(i)) + + ! Intersecting a regular plane + elseif (rAbs(i) > maxDist) then + maxDist = rAbs(i) + n = ZERO + n(i) = sign(ONE, rl(i)) + end if + + end do + + ! Normalise + n = n / norm2(n) + + end function normal + + !! !! Return to uninitialised state !! @@ -396,6 +440,8 @@ subroutine transformBC(self, r, u) ! Calculate halfwidth reduced by the surface_tolerance ! Necessary to capture particles at the boundary + ! SHOULD THE SURFTOL BE RELATIVE? MIGHT RISK LOSING PARTICLES IN SMALL GEOMETRIES IF NOT + !a_bar = self % halfwidth *(ONE - self % surfTol()) a_bar = self % halfwidth - self % surfTol() ! Calculate distance (in # of transformations) in each direction diff --git a/Geometry/Surfaces/squareCylinder_class.f90 b/Geometry/Surfaces/squareCylinder_class.f90 index 5e41d95ec..a4059eb54 100644 --- a/Geometry/Surfaces/squareCylinder_class.f90 +++ b/Geometry/Surfaces/squareCylinder_class.f90 @@ -9,6 +9,8 @@ module squareCylinder_class implicit none private + real(defReal), parameter :: EDGE_TOL = 1.0E-6 + !! !! Square cylinder aligned with x,y or z axis !! @@ -64,6 +66,7 @@ module squareCylinder_class procedure :: evaluate procedure :: distance procedure :: going + procedure :: normal procedure :: kill procedure :: setBC procedure :: explicitBC @@ -288,7 +291,7 @@ end function distance !! Works by: !! 1) Determine a plane in which direction is the closest !! 2) Use normal for this plane and project distance - !! 3) Determinie halfspace based on sign of the projection + !! 3) Determine halfspace based on sign of the projection !! !! Note: !! For parallel direction halfspace is asigned by the sign of `evaluate` result. @@ -312,7 +315,7 @@ pure function going(self, r, u) result(halfspace) ! Projection of direction on the normal proj = ul(maxCom) * sign(ONE, rl(maxCom)) - + halfspace = proj > ZERO ! Parallel direction @@ -322,6 +325,43 @@ pure function going(self, r, u) result(halfspace) end if end function going + + !! + !! Produces the normal for the square cylinder. Makes sure to check edges. + !! Does so by checking each cardinal direction and comparing. + !! + pure function normal(self, r, u) result(n) + class(squareCylinder), 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(2) :: rl, rAbs + integer(shortInt), dimension(2) :: p + integer(shortInt) :: i + + p = self % plane + + rl = r(self % plane) - self % origin + rAbs = abs(rl) - self % halfwidth + + ! Check each dimension, with a tolerance to detect edges + n = ZERO + + ! Intersecting an edge + if (abs(rAbs(1) - rAbs(2)) < EDGE_TOL) then + n(p(1)) = sign(ONE, rl(1)) + n(p(2)) = sign(ONE, rl(2)) + + ! Intersecting a plane + else + i = maxloc(rAbs,1) + n(p(i)) = sign(ONE, rl(i)) + end if + + ! Normalise + n = n / norm2(n) + + end function normal !! !! Return to uninitialised state @@ -458,7 +498,7 @@ subroutine transformBC(self, r, u) Ri = ceiling(abs(r(self % plane) - self % origin) / a_bar) / 2 ! Loop over directions & number of transformations - ! Becouse of the mix of 2D and 3D vectors to get right component use: + ! Because of the mix of 2D and 3D vectors to get right component use: ! i -> for 2D vectors ! ax -> for 3D vectors (r & u) axis : do i = 1, 2 diff --git a/Geometry/Surfaces/surface_inter.f90 b/Geometry/Surfaces/surface_inter.f90 index e0b28b2c2..46d193bba 100644 --- a/Geometry/Surfaces/surface_inter.f90 +++ b/Geometry/Surfaces/surface_inter.f90 @@ -42,6 +42,7 @@ module surface_inter !! evaluate -> Return remainder of the surface equation c = F(r) !! distance -> Return distance to the surface !! going -> Determine to which halfspace particle is currently going + !! normal -> Return the normal of the surface on which the particle is sat !! explicitBC -> Apply explicit BCs !! transformBC -> Apply transform BCs !! @@ -67,6 +68,7 @@ module surface_inter procedure(evaluate), deferred :: evaluate procedure(distance), deferred :: distance procedure(going), deferred :: going + procedure(normal), deferred :: normal procedure :: explicitBC procedure :: transformBC end type surface @@ -104,7 +106,7 @@ end subroutine init !! Return axis-aligned bounding box !! !! Note: - !! If bounding box is infinate in any axis its smalles value is -INF and highest INF + !! If bounding box is infinite in any axis its smallest value is -INF and highest INF !! !! Args: !! None @@ -180,6 +182,25 @@ pure function going(self, r, u) result(halfspace) real(defReal), dimension(3), intent(in) :: u logical(defBool) :: halfspace end function going + + !! + !! Return normal to the surface at the given point + !! + !! Args: + !! r [in] -> Position of the particle + !! u [in] -> DIrection of the particle. Assume norm2(u) = 1.0 + !! + !! Result: + !! Normal vector at the point, normalised + !! in +ve direction, returns INF + !! + pure function normal(self, r, u) result(n) + import :: surface, defReal + class(surface), intent(in) :: self + real(defReal), dimension(3), intent(in) :: r + real(defReal), dimension(3), intent(in) :: u + real(defReal), dimension(3) :: n + end function normal end interface contains diff --git a/Geometry/Surfaces/truncCylinder_class.f90 b/Geometry/Surfaces/truncCylinder_class.f90 index c45a56f08..a08ede803 100644 --- a/Geometry/Surfaces/truncCylinder_class.f90 +++ b/Geometry/Surfaces/truncCylinder_class.f90 @@ -9,6 +9,8 @@ module truncCylinder_class implicit none private + real(defReal), parameter :: EDGE_TOL = 1.0E-6 + !! !! Finite length cylinder aligned with one of the co-ord axis (x,y or z) !! @@ -62,6 +64,7 @@ module truncCylinder_class procedure :: evaluate procedure :: distance procedure :: going + procedure :: normal procedure :: kill procedure :: setBC procedure :: explicitBC @@ -226,7 +229,7 @@ pure function distance(self, r, u) result(d) a = ONE - u(ax)**2 delta = k*k - a*c1 ! Technically delta/4 - ! Find closes & furthest distance + ! Find closest & furthest distance if (delta <= ZERO .or. a == ZERO) then ! No intersection far = INF near = sign(INF, c1) ! If ray is parallel inside the cylinder it must be fully contained @@ -338,6 +341,49 @@ pure function going(self, r, u) result(halfspace) end if end function going + + !! + !! Produces the normal for the truncated cylinder. Makes sure to check edges. + !! Does so by checking each cardinal direction and comparing. + !! + pure function normal(self, r, u) result(n) + class(truncCylinder), 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(2) :: rp + integer(shortInt), dimension(2) :: p + real(defReal) :: dCirc0, dCirc, dCyl + integer(shortInt) :: a + + n = ZERO + + ! Get plane and axis indexes into shorter variables (for clarity) + p = self % plane + a = self % axis + + ! Check distances to see whether intersecting the circle or cylinder or both + ! Cylinder distance + rp = r(p) - self % origin(p) + dCyl = sum(rp**2) - self % r * self % r + + ! Circle distance + dCirc0 = r(a) - self % origin(a) + dCirc = abs(dCirc0) - self % a + + ! Check for overlap + if (abs(dCyl - dCirc) < EDGE_TOL) then + n(p) = rp / norm2(rp) + n(a) = sign(ONE, dCirc0) + elseif (abs(dCyl) < abs(dCirc)) then + n(p) = rp + else + n(a) = sign(ONE, dCirc0) + end if + + n = n / norm2(n) + + end function normal !! !! Return to uninitialised state diff --git a/Geometry/Tests/geometryStd_iTest.f90 b/Geometry/Tests/geometryStd_iTest.f90 index 75def8e56..d9fbbac99 100644 --- a/Geometry/Tests/geometryStd_iTest.f90 +++ b/Geometry/Tests/geometryStd_iTest.f90 @@ -131,7 +131,7 @@ subroutine test_lattice_geom() u = [ZERO, -ONE, ZERO] call coords % init(r, u) - ! Collosion movement + ! Collision movement maxDist = 1.0_defReal call geom % moveGlobal(coords, maxDist, event) @@ -161,7 +161,7 @@ subroutine test_lattice_geom() @assertEqual(idx, coords % matIdx) @assertEqual(0.26_defReal, maxDist, TOL) - !*** Normal Movment (easy case) + !*** Normal Movement (easy case) r = [-0.63_defReal, -0.63_defReal, 0.0_defReal] u = [ZERO, -ONE, ZERO] call coords % init(r, u) diff --git a/Geometry/Universes/Tests/cellUniverse_test.f90 b/Geometry/Universes/Tests/cellUniverse_test.f90 index 875f71403..3d90fbf2d 100644 --- a/Geometry/Universes/Tests/cellUniverse_test.f90 +++ b/Geometry/Universes/Tests/cellUniverse_test.f90 @@ -187,7 +187,6 @@ subroutine test_distance() @assertEqual(ref, d, TOL * ref) @assertEqual(surfs % getIdx(1), surfIdx ) - ! ** In local cell 2 distance to surface 2 pos % r = [7.0_defReal, 0.0_defReal, 0.0_defReal] pos % dir = [-ONE, ZERO, ZERO] @@ -260,6 +259,27 @@ subroutine test_cellOffset() end subroutine test_cellOffset + !! + !! Test getting a surface normal + !! +@Test + subroutine test_normal() + type(coord) :: pos + integer(shortInt) :: idx + real(defReal), dimension(3) :: n + + ! Cross from cell 1 to cell 2 + pos % r = [0.0_defReal, 2.0_defReal, 0.0_defReal] + pos % dir = [ONE, ZERO, ZERO] + pos % cellIdx = cells % getIdx(1) + pos % localId = 1 + + idx = surfs % getIdx(1) + n = uni % getNormal(idx, pos) + + @assertEqual([0, 1, 0], n) + + end subroutine test_normal end module cellUniverse_test diff --git a/Geometry/Universes/Tests/latUniverse_test.f90 b/Geometry/Universes/Tests/latUniverse_test.f90 index 6384dac7f..acf900cba 100644 --- a/Geometry/Universes/Tests/latUniverse_test.f90 +++ b/Geometry/Universes/Tests/latUniverse_test.f90 @@ -9,7 +9,7 @@ module latUniverse_test use coord_class, only : coord use surfaceShelf_class, only : surfaceShelf use cellShelf_class, only : cellShelf - use latUniverse_class, only : latUniverse + use latUniverse_class use funit implicit none @@ -421,5 +421,42 @@ subroutine test_cellOffset() end subroutine test_cellOffset + !! + !! Test getting a normal + !! +@Test + subroutine test_normal() + type(coord) :: pos + real(defReal), dimension(3) :: n + + ! Get normals at each edge + pos % r = [-1.0_defReal, 0.0_defReal, -0.5_defReal] + pos % dir = [-ONE, ONE, -ONE] + pos % dir = pos % dir / norm2(pos % dir) + + n = uni1 % getNormal(X_MIN, pos) + @assertEqual([-1, 0, 0], n) + + n = uni1 % getNormal(X_MAX, pos) + @assertEqual([1, 0, 0], n) + + n = uni1 % getNormal(Y_MIN, pos) + @assertEqual([0, -1, 0], n) + + n = uni1 % getNormal(Y_MAX, pos) + @assertEqual([0, 1, 0], n) + + n = uni1 % getNormal(Z_MIN, pos) + @assertEqual([0, 0, -1], n) + + n = uni1 % getNormal(Z_MAX, pos) + @assertEqual([0, 0, 1], n) + + pos % r = [3.0_defReal, 0.0_defReal, 0.0_defReal] + n = uni1 % getNormal(OUTLINE_SURF, pos) + @assertEqual([1, 0, 0], n) + + end subroutine test_normal + end module latUniverse_test diff --git a/Geometry/Universes/Tests/pinUniverse_test.f90 b/Geometry/Universes/Tests/pinUniverse_test.f90 index 603d76a62..f6cd274c6 100644 --- a/Geometry/Universes/Tests/pinUniverse_test.f90 +++ b/Geometry/Universes/Tests/pinUniverse_test.f90 @@ -278,5 +278,39 @@ subroutine test_edgeCases() end subroutine test_edgeCases + !! + !! Test getting a normal + !! +@Test + subroutine test_normal() + type(coord) :: pos + integer(shortInt) :: idx + real(defReal), dimension(3) :: n + + pos % r = [0.0_defReal, -1.5_defReal, 0.0_defReal] + pos % dir = [ZERO, ONE, ZERO] + pos % localId = 1 + + ! Surface should be innermost annulus + idx = MOVING_IN + n = uni % getNormal(idx, pos) + @assertEqual([0, -1, 0], n) + + ! Surface should be second innermost annulus - results will be identical + pos % r = [0.0_defReal, -1.8_defReal, 0.0_defReal] + pos % localId = 2 + idx = MOVING_OUT + n = uni % getNormal(idx, pos) + @assertEqual([0, -1, 0], n) + + ! What about outermost? + pos % r = [9999.9_defReal, 0.0_defReal, 0.0_defReal] + pos % dir = [ZERO, ONE, ZERO] + pos % localId = 3 + idx = MOVING_OUT + n = uni % getNormal(idx, pos) + @assertEqual([1, 0, 0], n) + + end subroutine test_normal end module pinUniverse_test diff --git a/Geometry/Universes/Tests/rootUniverse_test.f90 b/Geometry/Universes/Tests/rootUniverse_test.f90 index e3af614cb..dae73ad5b 100644 --- a/Geometry/Universes/Tests/rootUniverse_test.f90 +++ b/Geometry/Universes/Tests/rootUniverse_test.f90 @@ -192,6 +192,25 @@ subroutine test_cellOffset() end subroutine test_cellOffset + !! + !! Test getting normal + !! +@Test + subroutine test_normal() + type(coord) :: pos + integer(shortInt) :: idx + real(defReal), dimension(3) :: n + + ! Cross into outside + pos % r = [2.0_defReal, 0.0_defReal, 0.0_defReal] + pos % dir = [ONE, ZERO, ZERO] + + idx = surfs % getIdx(1) + + n = uni % getNormal(idx, pos) + @assertEqual([1, 0, 0], n) + + end subroutine test_normal end module rootUniverse_test diff --git a/Geometry/Universes/cellUniverse_class.f90 b/Geometry/Universes/cellUniverse_class.f90 index ce528ab0f..d54928c12 100644 --- a/Geometry/Universes/cellUniverse_class.f90 +++ b/Geometry/Universes/cellUniverse_class.f90 @@ -61,6 +61,7 @@ module cellUniverse_class procedure :: distance procedure :: cross procedure :: cellOffset + procedure :: getNormal end type cellUniverse contains @@ -203,6 +204,34 @@ function cellOffset(self, coords) result (offset) offset = ZERO end function cellOffset + + !! + !! Return normal for the given surfIdx at a given point + !! + !! See universe_inter for details. + !! + function getNormal(self, surfIdx, coords) result (normal) + class(cellUniverse), intent(in) :: self + integer(shortInt), intent(in) :: surfIdx + type(coord), intent(in) :: coords + real(defReal), dimension(3) :: normal + integer(shortInt) :: cIdx + character(100), parameter :: Here = 'getNormal (cellUniverse_class.f90)' + + ! Ensure that the current cell of coords has surfIdx as a component + cIdx = coords % localID + + ! Ensure that there are this many cells + if (cIdx > size(self % cells) .or. cIdx < 1) call fatalError(Here, & + 'cIdx is invalid: '//numToChar(cIdx)) + + ! Ensure that the surfIdx is valid + if (surfIdx < 1) call fatalError(Here, & + 'surfIdx is invalid: '//numToChar(surfIdx)) + + normal = self % cells(cIdx) % ptr % getNormal(surfIdx, coords % r, coords % dir) + + end function getNormal !! !! Return to uninitialised state diff --git a/Geometry/Universes/latUniverse_class.f90 b/Geometry/Universes/latUniverse_class.f90 index ca950e9cd..c08b0d8d3 100644 --- a/Geometry/Universes/latUniverse_class.f90 +++ b/Geometry/Universes/latUniverse_class.f90 @@ -17,8 +17,9 @@ module latUniverse_class ! Parameters ! Note X/Y/Z MIN/MAX MUST have the value they have (or generated by a function) - integer(shortInt), parameter :: X_MIN = -1, X_MAX = -2, Y_MIN = -3, Y_MAX = -4, Z_MIN = -5, & - Z_MAX = -6, OUTLINE_SURF = -7 + ! Public so they can be accessed by unit tests + integer(shortInt), public, parameter :: X_MIN = -1, X_MAX = -2, Y_MIN = -3, Y_MAX = -4, & + Z_MIN = -5, Z_MAX = -6, OUTLINE_SURF = -7 !! !! 2D or 3D Cartesian lattice with constant pitch @@ -90,6 +91,7 @@ module latUniverse_class procedure :: distance procedure :: cross procedure :: cellOffset + procedure :: getNormal end type latUniverse contains @@ -341,6 +343,48 @@ function cellOffset(self, coords) result (offset) end if end function cellOffset + + !! + !! Return normal for the given surfIdx at a given point + !! + !! See universe_inter for details. + !! + function getNormal(self, surfIdx, coords) result (normal) + class(latUniverse), intent(in) :: self + integer(shortInt), intent(in) :: surfIdx + type(coord), intent(in) :: coords + real(defReal), dimension(3) :: normal + character(100), parameter :: Here = 'getNormal (latUniverse_class.f90)' + + ! Convert surfIdx to appropriate axis + select case(surfIdx) + case(X_MIN) + normal = [-1, 0, 0] + + case(X_MAX) + normal = [1, 0, 0] + + case(Y_MIN) + normal = [0, -1, 0] + + case(Y_MAX) + normal = [0, 1, 0] + + case(Z_MIN) + normal = [0, 0, -1] + + case(Z_MAX) + normal = [0, 0, 1] + + case(OUTLINE_SURF) + normal = self % outline % normal(coords % r, coords % dir) + + case default + call fatalError(Here, 'Unrecognised surfIdx: '//numToChar(surfIdx)) + + end select + + end function getNormal !! !! Return to uninitialised state diff --git a/Geometry/Universes/pinUniverse_class.f90 b/Geometry/Universes/pinUniverse_class.f90 index e22fae1cc..f120f9057 100644 --- a/Geometry/Universes/pinUniverse_class.f90 +++ b/Geometry/Universes/pinUniverse_class.f90 @@ -58,6 +58,7 @@ module pinUniverse_class procedure :: distance procedure :: cross procedure :: cellOffset + procedure :: getNormal end type pinUniverse contains @@ -234,8 +235,8 @@ end subroutine distance !! subroutine cross(self, coords, surfIdx) class(pinUniverse), intent(inout) :: self - type(coord), intent(inout) :: coords - integer(shortInt), intent(in) :: surfIdx + type(coord), intent(inout) :: coords + integer(shortInt), intent(in) :: surfIdx character(100), parameter :: Here = 'cross (pinUniverse_class.f90)' if (surfIdx == MOVING_IN) then @@ -265,6 +266,38 @@ function cellOffset(self, coords) result (offset) offset = ZERO end function cellOffset + + !! + !! Return normal for the given surfIdx at a given point + !! + !! See universe_inter for details. + !! + function getNormal(self, surfIdx, coords) result (normal) + class(pinUniverse), intent(in) :: self + integer(shortInt), intent(in) :: surfIdx + type(coord), intent(in) :: coords + real(defReal), dimension(3) :: normal + integer(shortInt) :: cIdx + character(100), parameter :: Here = 'getNormal (pinUniverse_class.f90)' + + ! Local ID and surfIdx should be sufficient to identify which annulus + ! is being crossed and therefore which normal to produce + cIdx = coords % localID + + if (surfIdx == MOVING_IN) then + cIdx = cIdx + elseif (surfIdx == MOVING_OUT) then + cIdx = cIdx - 1 + else + call fatalError(Here, 'Unrecognised surfIdx: '//numToChar(surfIdx)) + end if + + if (cIdx > size(self % annuli) .or. cIdx < 1) call fatalError(Here,& + 'Invalid cIdx requested: '// numToChar(cIdx)) + + normal = self % annuli(cIdx) % normal(coords % r, coords % dir) + + end function getNormal !! !! Return to uninitialised state diff --git a/Geometry/Universes/rootUniverse_class.f90 b/Geometry/Universes/rootUniverse_class.f90 index 54f6d76be..2d9051b1e 100644 --- a/Geometry/Universes/rootUniverse_class.f90 +++ b/Geometry/Universes/rootUniverse_class.f90 @@ -56,6 +56,7 @@ module rootUniverse_class procedure :: distance procedure :: cross procedure :: cellOffset + procedure :: getNormal ! Subclass procedures procedure :: border @@ -193,6 +194,26 @@ function cellOffset(self, coords) result (offset) offset = ZERO end function cellOffset + + !! + !! Return normal for the given surfIdx at a given point + !! + !! See universe_inter for details. + !! + function getNormal(self, surfIdx, coords) result (normal) + class(rootUniverse), intent(in) :: self + integer(shortInt), intent(in) :: surfIdx + type(coord), intent(in) :: coords + real(defReal), dimension(3) :: normal + character(100), parameter :: Here = 'getNormal (rootUniverse_class.f90)' + + ! Ensure surfIdx is the border + if (surfIdx /= self % border()) call fatalError(Here, & + 'Provided surfIdx is not the border. SurfIdx: '//numToChar(surfIdx)) + + normal = self % surf % normal(coords % r, coords % dir) + + end function getNormal !! !! Return to uninitialised state diff --git a/Geometry/Universes/universe_inter.f90 b/Geometry/Universes/universe_inter.f90 index c903a50bd..967306838 100644 --- a/Geometry/Universes/universe_inter.f90 +++ b/Geometry/Universes/universe_inter.f90 @@ -81,6 +81,7 @@ module universe_inter procedure(distance), deferred :: distance procedure(cross), deferred :: cross procedure(cellOffset), deferred :: cellOffset + procedure(getNormal), deferred :: getNormal end type universe abstract interface @@ -208,6 +209,25 @@ function cellOffset(self, coords) result (offset) type(coord), intent(in) :: coords real(defReal), dimension(3) :: offset end function cellOffset + + !! + !! Return normal for the given surfIdx at a given point + !! + !! Args: + !! surfIdx [in] -> Idx of surface being crossed + !! coords [in] -> Coordinates placed in the universe (after transformations and with + !! local ID set). + !! + !! Result: + !! Surface normal (3D position vector). + !! + function getNormal(self, surfIdx, coords) result (normal) + import :: universe, shortInt, coord, defReal + class(universe), intent(in) :: self + integer(shortInt), intent(in) :: surfIdx + type(coord), intent(in) :: coords + real(defReal), dimension(3) :: normal + end function getNormal end interface diff --git a/Geometry/coord_class.f90 b/Geometry/coord_class.f90 index d7900d9a7..2fc11667f 100644 --- a/Geometry/coord_class.f90 +++ b/Geometry/coord_class.f90 @@ -1,7 +1,7 @@ module coord_class use numPrecision - use universalVariables, only : HARDCODED_MAX_NEST, FP_REL_TOL + use universalVariables, only : HARDCODED_MAX_NEST, FP_REL_TOL, UNDEF_MAT use genericProcedures, only : rotateVector, fatalError, numToChar implicit none @@ -81,11 +81,12 @@ module coord_class !! cell -> Return cellIdx at the lowest level !! assignPosition -> Set position and take ABOVE GEOMETRY !! assignDirection -> Set direction and do not change coordList state + !! assignDirectionBottom -> Set direction in the lowest universe level and do not change coordList state !! type, public :: coordList integer(shortInt) :: nesting = 0 type(coord), dimension(HARDCODED_MAX_NEST) :: lvl - integer(shortInt) :: matIdx = -3 + integer(shortInt) :: matIdx = UNDEF_MAT integer(shortInt) :: uniqueId = -3 contains ! Build procedures @@ -107,6 +108,7 @@ module coord_class procedure :: cell procedure :: assignPosition procedure :: assignDirection + procedure :: assignDirectionLevel end type coordList @@ -205,7 +207,7 @@ elemental subroutine kill(self) class(coordList), intent(inout) :: self self % nesting = 0 - self % matIdx = -3 + self % matIdx = UNDEF_MAT self % uniqueID = -3 ! Kill coordinates @@ -243,7 +245,7 @@ elemental function isAbove(self) result(isIt) class(coordList), intent(in) :: self logical(defBool) :: isIt - isIt = (self % matIdx < 0) .and. (self % uniqueID < 0) .and. (self % nesting == 1) + isIt = (self % matIdx == UNDEF_MAT) .and. (self % uniqueID < 0) .and. (self % nesting == 1) end function isAbove @@ -294,7 +296,7 @@ elemental subroutine takeAboveGeom(self) class(coordList), intent(inout) :: self self % nesting = 1 - self % matIdx = -3 + self % matIdx = UNDEF_MAT self % uniqueID = -3 end subroutine takeAboveGeom @@ -466,4 +468,40 @@ pure subroutine assignDirection(self, u) end subroutine assignDirection + !! + !! Assign new direction at given universe level, propagating changes to the top. + !! Sets the level of the coordList to that provided + !! + !! Args: + !! u [in] -> New normalised direction at given level + !! level [in] -> Level at which the direction is set + !! + !! NOTE: + !! Does not check if u is normalised! + !! + subroutine assignDirectionLevel(self, u, level) + class(coordList), intent(inout) :: self + real(defReal), dimension(3), intent(in) :: u + integer(shortInt), intent(in) :: level + integer(shortInt) :: i + character(100),parameter :: Here = 'assignDirectionLevel (coord_class.f90)' + + if (level > self % nesting) call fatalError(Here,'Input level exceeded current nesting') + + ! Assign new direction in lowest level + self % lvl(level) % dir = u + + ! Propage the change to higher levels + do i = level - 1, 1, -1 + if(self % lvl(i+1) % isRotated) then + self % lvl(i) % dir = matmul(transpose(self % lvl(i+1) % rotMat), self % lvl(i+1) % dir) + else + self % lvl(i) % dir = self % lvl(i+1) % dir + end if + end do + + self % nesting = level + + end subroutine assignDirectionLevel + end module coord_class diff --git a/Geometry/geometryStd_class.f90 b/Geometry/geometryStd_class.f90 index d53e33020..a9eb790d9 100644 --- a/Geometry/geometryStd_class.f90 +++ b/Geometry/geometryStd_class.f90 @@ -60,6 +60,7 @@ module geometryStd_class procedure :: move_noCache procedure :: move_withCache procedure :: moveGlobal + procedure :: moveNoBC procedure :: teleport procedure :: activeMats @@ -356,6 +357,95 @@ subroutine moveGlobal(self, coords, maxDist, event) call self % placeCoord(coords) end subroutine moveGlobal + + !! + !! Given coordinates placed in the geometry move point through the geometry. + !! Does not apply boundary conditions. Also provides normal of struck surfaces. + !! Used for ray plotting. + !! + !! See geometry_inter for details + !! + subroutine moveNoBC(self, coords, maxDist, event, n) + class(geometryStd), intent(in) :: self + type(coordList), intent(inout) :: coords + real(defReal), intent(inout) :: maxDist + integer(shortInt), intent(out) :: event + real(defReal), dimension(3), intent(out) :: n + integer(shortInt) :: surfIdx, level + real(defReal) :: dist + class(surface), pointer :: surf + class(universe), pointer :: uni + type(coordList) :: coordsTemp + character(100), parameter :: Here = 'moveNoBC (geometryStd_class.f90)' + + n = ZERO + + ! Find distance to the next surface + call self % closestDist(dist, surfIdx, level, coords) + + ! Moves within cell + ! Note: no normal should be produced since there is no + ! collision with a surface + if (maxDist <= dist .or. level == 0) then + + if (level > 1) then + call coords % moveLocal(maxDist, coords % nesting) + + ! Catch in case moving through the outside region + else + call coords % moveGlobal(dist) + + ! Place back in geometry + call self % placeCoord(coords) + end if + maxDist = maxDist ! Left for explicitness. Compiler will not stand it anyway + event = COLL_EV + + else if (surfIdx == self % geom % borderIdx .and. level == 1) then ! Hits domain boundary + ! Move global to the boundary + call coords % moveGlobal(dist) + event = BOUNDARY_EV + maxDist = dist + + ! Place back in geometry + call self % placeCoord(coords) + + ! Produce normal at domain boundary + surf => self % geom % surfs % getPtr(self % geom % borderIdx) + n = surf % normal(coords % lvl(1) % r, coords % lvl(1) % dir) + + else ! Crosses to different local cell + ! Move to boundary at hit level + call coords % moveLocal(dist, level) + event = CROSS_EV + maxDist = dist + + ! Get universe and cross to the next cell + uni => self % geom % unis % getPtr_fast(coords % lvl(level) % uniIdx) + call uni % cross(coords % lvl(level), surfIdx) + + ! Get material + call self % diveToMat(coords, level) + + ! Obtain surface normal in the local reference frame + ! If surfIdx is positive, this is an actual surface + if (surfIdx > 0) then + surf => self % geom % surfs % getPtr(surfIdx) + n = surf % normal(coords % lvl(coords % nesting) % r, coords % lvl(coords % nesting) % dir) + + ! If negative, this is an implied surface from a structured universe, e.g., lattice + else + n = uni % getNormal(surfIdx, coords % lvl(level)) + end if + + ! Rotate normal to the top universe reference frame + coordsTemp = coords + call coordsTemp % assignDirectionLevel(n, level) + n = coordsTemp % lvl(1) % dir + + end if + + end subroutine moveNoBC !! !! Move a particle in the top level without stopping diff --git a/Geometry/geometry_inter.f90 b/Geometry/geometry_inter.f90 index 973007549..73051eb3f 100644 --- a/Geometry/geometry_inter.f90 +++ b/Geometry/geometry_inter.f90 @@ -1,7 +1,7 @@ module geometry_inter use numPrecision - use universalVariables, only : X_AXIS, Y_AXIS, Z_AXIS, HARDCODED_MAX_NEST + use universalVariables, only : X_AXIS, Y_AXIS, Z_AXIS, HARDCODED_MAX_NEST, INFINITY, COLL_EV use genericProcedures, only : fatalError use dictionary_class, only : dictionary use charMap_class, only : charMap @@ -41,12 +41,15 @@ module geometry_inter procedure(move_noCache), deferred :: move_noCache procedure(move_withCache), deferred :: move_withCache procedure(moveGlobal), deferred :: moveGlobal + procedure(moveNoBC), deferred :: moveNoBC procedure(teleport), deferred :: teleport procedure(activeMats), deferred :: activeMats ! Common procedures procedure :: slicePlot - procedure :: voxelplot + procedure :: voxelPlot + procedure :: rayPlot + procedure :: phongTrace end type geometry abstract interface @@ -225,6 +228,41 @@ subroutine moveGlobal(self, coords, maxDist, event) integer(shortInt), intent(out) :: event end subroutine moveGlobal + !! + !! Given coordinates placed in the geometry move point through the geometry + !! + !! Move by up to maxDis stopping at domain boundary or until matIdx or uniqueID changes + !! When particle hits boundary, boundary conditions are NOT applied before returning. + !! Also supplies the normal of the surface struck. Normal is not defined if a surface is + !! not struck, i.e., event is not either BOUNDARY_EV or CROSS_EV + !! + !! Following events can be returned: + !! COLL_EV -> Particle moved by entire maxDist. Collision happens + !! BOUNDARY_EV -> Particle hit domain boundary + !! CROSS_EV -> Partilce crossed to a region with different matIdx or uniqueID + !! LOST_EV -> Something gone wrong in tracking and particle is lost + !! + !! Args: + !! coords [inout] -> Coordinate list of the particle to be moved through the geometry + !! maxDict [inout] -> Maximum distance to move the position. If movment is stopped + !! prematurely (e.g. hitting boundary), maxDist is set to the distance the particle has + !! moved by. + !! event [out] -> Event flag that specifies what finished the movement. + !! n [out] -> Normal vector of the surface struck at the end of the movement + !! + !! Errors: + !! If coords is not placed in the geometry behaviour is unspecified + !! If maxDist < 0.0 behaviour is unspecified + !! + subroutine moveNoBC(self, coords, maxDist, event, n) + import :: geometry, coordList, defReal, shortInt + class(geometry), intent(in) :: self + type(coordList), intent(inout) :: coords + real(defReal), intent(inout) :: maxDist + integer(shortInt), intent(out) :: event + real(defReal), dimension(3), intent(out) :: n + end subroutine moveNoBC + !! !! Move a particle in the top level without stopping !! @@ -425,5 +463,146 @@ subroutine voxelPlot(self, img, centre, what, width) end subroutine voxelPlot + !! + !! Generates images by ray tracing using Phong's approximation, as described in + !! Gavin Ridley's SNA 2024 paper. + !! + subroutine rayPlot(self, brightness, matIDs, camera, light, M, mats, fov, diffuse) + class(geometry), intent(in) :: self + real(defReal), dimension(:,:), intent(inout) :: brightness + integer(shortInt), dimension(:,:), intent(inout) :: matIDs + real(defReal), dimension(3), intent(in) :: camera + real(defReal), dimension(3), intent(in) :: light + real(defReal), dimension(3,3), intent(in) :: M + integer(shortInt), dimension(:), intent(in) :: mats + real(defReal), intent(in) :: fov + real(defReal), intent(in) :: diffuse + integer(shortInt) :: iv, nh, nv + integer(shortInt), save :: ih, matIdx + real(defReal), save :: bright + real(defReal) :: focalDist, dx, dy + real(defReal), dimension(3), save :: vec + type(coordlist), save :: ray + !$omp threadprivate(vec, ih, ray, matIdx, bright) + + ! The results is apparently insensitive to this value, as recommended by Ridley + focalDist = 10 + + nh = size(brightness, 1) + nv = size(brightness, 2) + dx = 2 * focalDist * tan(HALF * fov) + dy = real(nv, defReal) * dx / nh + + ! Loop over vertical pixels + !$omp parallel do + do iv = 1, nv + do ih = 1, nh + + vec = [focalDist, -dx * HALF + dx * real((ih - 1), defReal) / nh, -dy * HALF + dy * real((iv - 1), defReal) / nv] + vec = vec / norm2(vec) + + ! Place the ray in the right position and direction + call ray % init(camera, matmul(M,vec)) + + ! Get local material and its illumination + call self % phongTrace(ray, matIdx, bright, mats, diffuse, light) + matIDs(ih, iv) = matIdx + brightness(ih, iv) = bright + + end do + end do + !$omp end parallel do + + + end subroutine rayPlot + + !! + !! Procedure for tracing from a camera until an opaque object is hit. + !! Then calculates the luminous contribution from a light source. + !! + subroutine phongTrace(self, ray, matIdx, bright, mats, diffuse, light) + class(geometry), intent(in) :: self + type(coordlist), intent(inout) :: ray + integer(shortInt), intent(out) :: matIdx + real(defReal), intent(out) :: bright + integer(shortInt), dimension(:), intent(in) :: mats + real(defReal), intent(in) :: diffuse + real(defReal), dimension(3), intent(in) :: light + logical(defBool) :: lightTrace + real(defReal) :: dist, dotP, dNudge + integer(shortInt) :: event, matIdx0 + real(defReal), dimension(3) :: dir, normal0, normal + + lightTrace = .true. + + call self % placeCoord(ray) + + matIdx = ray % matIdx + + ! Trace until an opaque material is hit + do while (any(mats == matIdx)) + + dist = INFINITY + call self % moveNoBC(ray, dist, event, normal0) + + matIdx = ray % matIdx + + ! If distance is infinite or particle collided, not pointing towards anything. + ! Hence, short circuit the trace + if ((dist == INFINITY) .or. (event == COLL_EV)) then + lightTrace = .false. + exit + end if + + end do + + bright = ZERO + + ! Trace to light source to calculate brightness at the struck point + if (lightTrace) then + + ! Make sure normal is oriented correctly + dotP = dot_product(ray % lvl(1) % dir, normal0) + if (dotP > ZERO) normal0 = -normal0 + + ! Flip ray and nudge it backwards to get it out of the opaque material + ! This nudge may cause problems! Any better solutions? + dNudge = 1.0E-10 + call ray % assignDirection(-ray % lvl(1) % dir) + call self % moveNoBC(ray, dNudge, event, normal) + + ! Reorient ray towards light + dir = light - ray % lvl(1) % r + dir = dir / norm2(dir) + call ray % assignDirection(dir) + call self % placeCoord(ray) + + ! Compute product betwen normal and light direction + dotP = max(dot_product(normal0, dir), ZERO) + + ! Does the ray fly directly to the light? + matIdx0 = ray % matIdx + do while (any(mats == matIdx0)) + + ! Find maximum flight distance + dist = norm2(ray % lvl(1) % r - light) + call self % moveNoBC(ray, dist, event, normal) + + matIdx0 = ray % matIdx + + ! The ray flew straight to the light + if (event == COLL_EV) then + bright = (ONE - diffuse) * dotP + exit + end if + + end do + + ! Add diffuse/ambient contribution + bright = bright + diffuse + + end if + + end subroutine phongTrace end module geometry_inter diff --git a/InputFiles/Benchmarks/BEAVRS/BEAVRS2D b/InputFiles/Benchmarks/BEAVRS/BEAVRS2D index 0ace8d16c..899a59c71 100644 --- a/InputFiles/Benchmarks/BEAVRS/BEAVRS2D +++ b/InputFiles/Benchmarks/BEAVRS/BEAVRS2D @@ -114,7 +114,8 @@ geometry { thickGrid {type simpleCell; id 55; surfaces (90 -91); filltype mat; material Inconel;} thinGrid {type simpleCell; id 56; surfaces (-92 93); filltype mat; material Inconel;} - pressureVessel { type simpleCell; id 7; surfaces (-1 2); filltype mat; material CarbonSteel;} + ! Does not use surface 1 to define it as that bounds the geometry + pressureVessel { type simpleCell; id 7; surfaces (2); filltype mat; material CarbonSteel;} RPVLiner { type simpleCell; id 8; surfaces (-2 3); filltype mat; material SS304;} outerWater1 {type simpleCell; id 9; surfaces (-3 4 ); filltype mat; material Water;} diff --git a/InputFiles/Benchmarks/BEAVRS/BEAVRS b/InputFiles/Benchmarks/BEAVRS/BEAVRS_ARO similarity index 97% rename from InputFiles/Benchmarks/BEAVRS/BEAVRS rename to InputFiles/Benchmarks/BEAVRS/BEAVRS_ARO index 52e5cb575..dd3ed25c5 100644 --- a/InputFiles/Benchmarks/BEAVRS/BEAVRS +++ b/InputFiles/Benchmarks/BEAVRS/BEAVRS_ARO @@ -3,14 +3,12 @@ !! This is not a fully faithful replica of BEAVRS at HZP: !! all rods are fully withdrawn. !! In reality, some rods are partially inserted. -!! TODO: add rods in the appropriate assemblies at the -!! correct heights! !! type eigenPhysicsPackage; -pop 10000000; +pop 1000000; active 50; -inactive 200; +inactive 150; XSdata ce; dataType ce; @@ -216,7 +214,8 @@ geometry { thickGrid {type simpleCell; id 555; surfaces (90 ); filltype mat; material Inconel;} thinGrid {type simpleCell; id 556; surfaces (92 ); filltype mat; material Zircaloy;} - pressureVessel { type simpleCell; id 7; surfaces (-1 2); filltype mat; material CarbonSteel;} + // Don't need to bound PV by 1 since it is the bounding surface of the geometry. + pressureVessel { type simpleCell; id 7; surfaces ( 2); filltype mat; material CarbonSteel;} RPVLiner { type simpleCell; id 8; surfaces (-2 3); filltype mat; material SS304;} outerWater1 {type simpleCell; id 9; surfaces (-3 4 ); filltype mat; material Water;} diff --git a/InputFiles/Benchmarks/BEAVRS/BEAVRS_HZP b/InputFiles/Benchmarks/BEAVRS/BEAVRS_HZP new file mode 100644 index 000000000..5a740ce28 --- /dev/null +++ b/InputFiles/Benchmarks/BEAVRS/BEAVRS_HZP @@ -0,0 +1,2410 @@ +!! +!! 3D BEAVRS benchmark +!! At HZP, the D bank of rods is partially inserted up +!! to 115 steps inserted / 113 steps withdrawn. +!! A step corresponds to an increment of 1.58193cm +!! +type eigenPhysicsPackage; + +pop 1000000; +active 50; +inactive 250; +XSdata ce; +dataType ce; + +collisionOperator { neutronCE {type neutronCEstd;}} + +transportOperator { + !type transportOperatorDT; + type transportOperatorHT; cache 1; + } + +inactiveTally { + shannon { + type shannonEntropyClerk; + map {type multiMap; + maps (xax yax zax); + xax { type spaceMap; grid lin; min -161.2773; max 161.2773; N 15; axis x;} + yax { type spaceMap; grid lin; min -161.2773; max 161.2773; N 15; axis y;} + zax { type spaceMap; grid lin; min 36.748; max 402.508; N 15; axis z;} + } + cycles 200; + } + +} + +activeTally { + pinFissRadial { type collisionClerk; response (fission); fission { type macroResponse; MT -6;} + map {type multiMap; maps (xax yax); + xax {type spaceMap; axis x; grid lin; N 255; min -161.2773; max 161.2773; } + yax {type spaceMap; axis y; grid lin; N 255; min -161.2773; max 161.2773; } + } + } + assemblyFissRadial { type collisionClerk; response (fission); fission { type macroResponse; MT -6;} + map {type multiMap; maps (xax yax); + xax {type spaceMap; axis x; grid lin; N 15; min -161.2773; max 161.2773; } + yax {type spaceMap; axis y; grid lin; N 15; min -161.2773; max 161.2773; } + } + } + fissionAxial { type collisionClerk; response (fission); fission { type macroResponse; MT -6;} + map {type spaceMap; axis z; grid lin; N 60; min 36.748; max 402.508;} + } + fissionYZ { type collisionClerk; response (fission); fission { type macroResponse; MT -6;} + map {type multiMap; maps (yax zax); + yax {type spaceMap; axis y; grid lin; N 255; min -161.2773; max 161.2773; } + zax {type spaceMap; axis z; grid lin; N 60; min 36.748; max 402.508;} + } + } +} + +geometry { + type geometryStd; + boundary ( 0 0 0 0 0 0); + graph {type shrunk;} + + surfaces { + + // thickness specifications for RPV and RPV liner + outerRPV { id 1; type zTruncCylinder; radius 241.3; origin (0.0 0.0 230.0); halfwidth 230; } + innerRPV { id 2; type zCylinder; radius 219.710; origin (0.0 0.0 0.0); } + innerRPVLiner { id 3; type zCylinder; radius 219.150; origin (0.0 0.0 0.0); } + + // thickness specifications for neutron shield + outerBoundNS { id 4; type zCylinder; radius 201.630; origin (0.0 0.0 0.0); } + innerBoundNS { id 5; type zCylinder; radius 194.84; origin (0.0 0.0 0.0); } + + // thickness specifications for core barrel + outerCoreBarrel { id 6; type zCylinder; radius 193.675; origin (0.0 0.0 0.0); } + innerCoreBarrel { id 7; type zCylinder; radius 187.96; origin (0.0 0.0 0.0); } + + // four planes that intersect to bound the Neutron shield panel + P1 { id 8; type plane; coeffs (-0.48480962025 0.87461970714 0.0 0.0);} + P2 { id 9; type plane; coeffs (-0.87461970714 0.48480962025 0.0 0.0);} + P3 { id 10; type plane; coeffs (-0.87461970714 -0.48480962025 0.0 0.0);} + P4 { id 11; type plane; coeffs (-0.48480962025 -0.87461970714 0.0 0.0);} + + // bounding widths for baffle on various sides + // right & left refers to the side of the reactor that it is on + // close/away refers to its location in relation to the LATTICE it is a part of + // (NOT the reactor itself) + rightClose { id 50; type plane; coeffs (1.0 0.0 0.0 8.36662);} + rightAway { id 51; type plane; coeffs (1.0 0.0 0.0 10.58912);} + leftClose { id 52; type plane; coeffs (-1.0 0.0 0.0 8.36662);} + leftAway { id 53; type plane; coeffs (-1.0 0.0 0.0 10.58912);} + bottomClose { id 54; type plane; coeffs (0.0 -1.0 0.0 8.36662);} + bottomAway { id 55; type plane; coeffs (0.0 -1.0 0.0 10.58912);} + topClose { id 56; type plane; coeffs (0.0 1.0 0.0 8.36662);} + topAway { id 57; type plane; coeffs (0.0 1.0 0.0 10.58912);} + + // thickness specifications for grid with thickness of 0.0198cm (Inconel) + pinThickGridInner { id 90; type zSquareCylinder; origin (0.0 0.0 0.0); halfwidth (0.61015 0.61015 0.0); } + pinThickGridOuter { id 91; type zSquareCylinder; origin (0.0 0.0 0.0); halfwidth (0.62992 0.62992 0.0); } + + // thickness specifications for grid with thickness of 0.0194cm (Zircaloy) + pinThinGridInner { id 92; type zSquareCylinder; origin (0.0 0.0 0.0); halfwidth (0.61049 0.61049 0.0); } + pinThinGridOuter { id 93; type zSquareCylinder; origin (0.0 0.0 0.0); halfwidth (0.62992 0.62992 0.0); } + + // inner and outer surfaces of assembly sleeves (both SS and Zircaloy) + assemblySleeveInner { id 94; type zSquareCylinder; origin (0.0 0.0 0.0); halfwidth (10.70864 10.70864 0.0); } + assemblySleeveOuter { id 95; type zSquareCylinder; origin (0.0 0.0 0.0); halfwidth (10.74798 10.74798 0.0); } + + + // Axial planes across core height + // Names are based on axial heights + plane460 { id 100; type plane; coeffs (0.0 0.0 1.0 460.0 ); } + plane431p876 { id 101; type plane; coeffs (0.0 0.0 1.0 431.876); } + plane423p049 { id 102; type plane; coeffs (0.0 0.0 1.0 423.049); } + plane421p532 { id 103; type plane; coeffs (0.0 0.0 1.0 421.532); } + plane419p704 { id 104; type plane; coeffs (0.0 0.0 1.0 419.704); } + plane417p164 { id 105; type plane; coeffs (0.0 0.0 1.0 417.164); } + plane415p164 { id 106; type plane; coeffs (0.0 0.0 1.0 415.164); } + plane411p806 { id 107; type plane; coeffs (0.0 0.0 1.0 411.806); } + plane403p778 { id 108; type plane; coeffs (0.0 0.0 1.0 403.778); } + plane402p508 { id 109; type plane; coeffs (0.0 0.0 1.0 402.508); } + plane401p238 { id 110; type plane; coeffs (0.0 0.0 1.0 401.238); } + plane364p725 { id 111; type plane; coeffs (0.0 0.0 1.0 364.725); } + plane359p01 { id 112; type plane; coeffs (0.0 0.0 1.0 359.01 ); } + plane312p528 { id 113; type plane; coeffs (0.0 0.0 1.0 312.528); } + plane306p813 { id 114; type plane; coeffs (0.0 0.0 1.0 306.813); } + plane260p331 { id 115; type plane; coeffs (0.0 0.0 1.0 260.331); } + plane254p616 { id 116; type plane; coeffs (0.0 0.0 1.0 254.616); } + plane208p134 { id 117; type plane; coeffs (0.0 0.0 1.0 208.134); } + plane202p419 { id 118; type plane; coeffs (0.0 0.0 1.0 202.419); } + plane155p937 { id 119; type plane; coeffs (0.0 0.0 1.0 155.937); } + plane150p222 { id 120; type plane; coeffs (0.0 0.0 1.0 150.222); } + plane143p428 { id 121; type plane; coeffs (0.0 0.0 1.0 143.428); } + plane103p74 { id 122; type plane; coeffs (0.0 0.0 1.0 103.74 ); } + plane98p025 { id 123; type plane; coeffs (0.0 0.0 1.0 98.025 ); } + plane41p828 { id 124; type plane; coeffs (0.0 0.0 1.0 41.828 ); } + plane40p558 { id 125; type plane; coeffs (0.0 0.0 1.0 40.558 ); } + plane40p52 { id 126; type plane; coeffs (0.0 0.0 1.0 40.52 ); } + plane39p958 { id 127; type plane; coeffs (0.0 0.0 1.0 39.958 ); } + plane38p66 { id 128; type plane; coeffs (0.0 0.0 1.0 38.66 ); } + plane37p1621 { id 129; type plane; coeffs (0.0 0.0 1.0 37.1621); } + plane36p748 { id 130; type plane; coeffs (0.0 0.0 1.0 36.748 ); } + plane35 { id 131; type plane; coeffs (0.0 0.0 1.0 35.0 ); } + plane20 { id 132; type plane; coeffs (0.0 0.0 1.0 20.0 ); } + plane0 { id 133; type plane; coeffs (0.0 0.0 1.0 0.0 ); } + + planeSteelBottom { id 134; type plane; coeffs (0.0 0.0 1.0 400.638); } + planeCRLowerBottom { id 135; type plane; coeffs (0.0 0.0 1.0 402.508); } // same as 109 on withdrawal + planeCRUpperBottom { id 136; type plane; coeffs (0.0 0.0 1.0 504.108); } // out of core on withdrawal + + plane322p1861 {id 137; type plane; coeffs (0 0 1 322.18609);} + plane220p586 {id 138; type plane; coeffs (0 0 1 220.5861);} + plane218p716 {id 139; type plane; coeffs (0 0 1 218.71609);} + + } + + cells { + + // assembly wrappers and surrounding water at various heights + wrapper1 {type simpleCell; id 2001; surfaces (94 -95 129 -126); filltype mat; material SS304;} + wrapper2 {type simpleCell; id 2002; surfaces (94 -95 123 -122); filltype mat; material Zircaloy;} + wrapper3 {type simpleCell; id 2003; surfaces (94 -95 120 -119); filltype mat; material Zircaloy;} + wrapper4 {type simpleCell; id 2004; surfaces (94 -95 118 -117); filltype mat; material Zircaloy;} + wrapper5 {type simpleCell; id 2005; surfaces (94 -95 116 -115); filltype mat; material Zircaloy;} + wrapper6 {type simpleCell; id 2006; surfaces (94 -95 114 -113); filltype mat; material Zircaloy;} + wrapper7 {type simpleCell; id 2007; surfaces (94 -95 112 -111); filltype mat; material Zircaloy;} + wrapper8 {type simpleCell; id 2008; surfaces (94 -95 107 -106); filltype mat; material SS304;} + + assemWater0 {type simpleCell; id 2009; surfaces (94 -129); filltype mat; material Water;} + assemWater1 {type simpleCell; id 2010; surfaces (94 -95 126 -123); filltype mat; material Water;} + assemWater2 {type simpleCell; id 2011; surfaces (94 -95 122 -120); filltype mat; material Water;} + assemWater3 {type simpleCell; id 2012; surfaces (94 -95 119 -118); filltype mat; material Water;} + assemWater4 {type simpleCell; id 2013; surfaces (94 -95 117 -116); filltype mat; material Water;} + assemWater5 {type simpleCell; id 2014; surfaces (94 -95 115 -114); filltype mat; material Water;} + assemWater6 {type simpleCell; id 2015; surfaces (94 -95 113 -112); filltype mat; material Water;} + assemWater7 {type simpleCell; id 2016; surfaces (94 -95 111 -107); filltype mat; material Water;} + assemWater8 {type simpleCell; id 2017; surfaces (94 -95 106); filltype mat; material Water;} + assemWaterEx {type simpleCell; id 2018; surfaces (95); filltype mat; material Water;} + + // assemblies inside the wrappers + assem1424 {type simpleCell; id 2019; surfaces (-94); filltype uni; universe 1424;} + assem1416 {type simpleCell; id 2020; surfaces (-94); filltype uni; universe 1416;} + assem1431 {type simpleCell; id 2021; surfaces (-94); filltype uni; universe 1431;} + assem60316 {type simpleCell; id 2022; surfaces (-94); filltype uni; universe 60316;} + assem603112 {type simpleCell; id 2023; surfaces (-94); filltype uni; universe 603112;} + assem60313 {type simpleCell; id 2024; surfaces (-94); filltype uni; universe 60313;} + assem60319 {type simpleCell; id 2025; surfaces (-94); filltype uni; universe 60319;} + assem15315 {type simpleCell; id 2026; surfaces (-94); filltype uni; universe 15315;} + assem15317 {type simpleCell; id 2027; surfaces (-94); filltype uni; universe 15317;} + assem15311 {type simpleCell; id 2028; surfaces (-94); filltype uni; universe 15311;} + assem153111 {type simpleCell; id 2029; surfaces (-94); filltype uni; universe 153111;} + assem1631 {type simpleCell; id 2030; surfaces (-94); filltype uni; universe 1631;} + assem2031 {type simpleCell; id 2031; surfaces (-94); filltype uni; universe 2031;} + assem1224 {type simpleCell; id 2032; surfaces (-94); filltype uni; universe 1224;} + assem1624 {type simpleCell; id 2033; surfaces (-94); filltype uni; universe 1624;} + + // unrodded assemblies in wrappers + assem2424 {type simpleCell; id 2034; surfaces (-94); filltype uni; universe 2424;} + assem2416 {type simpleCell; id 2035; surfaces (-94); filltype uni; universe 2416;} + assem2431 {type simpleCell; id 2036; surfaces (-94); filltype uni; universe 2431;} + assem70316 {type simpleCell; id 2037; surfaces (-94); filltype uni; universe 70316;} + assem703112 {type simpleCell; id 2038; surfaces (-94); filltype uni; universe 703112;} + assem70313 {type simpleCell; id 2039; surfaces (-94); filltype uni; universe 70313;} + assem70319 {type simpleCell; id 2040; surfaces (-94); filltype uni; universe 70319;} + assem25315 {type simpleCell; id 2041; surfaces (-94); filltype uni; universe 25315;} + assem25317 {type simpleCell; id 2042; surfaces (-94); filltype uni; universe 25317;} + assem25311 {type simpleCell; id 2043; surfaces (-94); filltype uni; universe 25311;} + assem253111 {type simpleCell; id 2044; surfaces (-94); filltype uni; universe 253111;} + assem2631 {type simpleCell; id 2045; surfaces (-94); filltype uni; universe 2631;} + assem3031 {type simpleCell; id 2046; surfaces (-94); filltype uni; universe 3031;} + assem2224 {type simpleCell; id 2047; surfaces (-94); filltype uni; universe 2224;} + assem2624 {type simpleCell; id 2048; surfaces (-94); filltype uni; universe 2624;} + + // rodded assemblies in wrappers + assem1425 {type simpleCell; id 2049; surfaces (-94); filltype uni; universe 1425;} + assem1417 {type simpleCell; id 2050; surfaces (-94); filltype uni; universe 1417;} + + // pin grids - thick at the top and bottom, thin in fuelled region + thickGrid {type simpleCell; id 555; surfaces (90 ); filltype mat; material Inconel;} + thinGrid {type simpleCell; id 556; surfaces (92 ); filltype mat; material Zircaloy;} + + // Don't need to bound PV by 1 since it is the bounding surface of the geometry. + pressureVessel { type simpleCell; id 7; surfaces (2); filltype mat; material CarbonSteel;} + RPVLiner { type simpleCell; id 8; surfaces (-2 3); filltype mat; material SS304;} + outerWater1 {type simpleCell; id 9; surfaces (-3 4 ); filltype mat; material Water;} + + // Neutron shields + NS1 { type simpleCell; id 10; surfaces (-4 5 -8 9); filltype mat; material SS304;} + NS2 { type simpleCell; id 11; surfaces (-4 5 8 -9); filltype mat; material SS304;} + NS3 { type simpleCell; id 12; surfaces (-4 5 -10 11); filltype mat; material SS304;} + NS4 { type simpleCell; id 13; surfaces (-4 5 10 -11); filltype mat; material SS304;} + + // Water in the arc between neutron shields + outerWaterSeg1 {type simpleCell; id 14; surfaces (-4 5 ); filltype mat; material Water;} + outerWater2 {type simpleCell; id 15; surfaces (-5 6 ); filltype mat; material Water;} + + // Outer core + coreBarrel { type simpleCell; id 16; surfaces (-6 7); filltype mat; material SS304;} + core {type simpleCell; id 17; surfaces (-5); filltype uni; universe 9999;} + + + // Gridded pins + + // THICK GRID + // 2.4% in grid + grid24Thick {type simpleCell; id 253; surfaces (-90); filltype uni; universe 24000;} + // guide tube in grid + gridGTThick {type simpleCell; id 254; surfaces (-90); filltype uni; universe 12000;} + // 3.1% in grid + grid31Thick {type simpleCell; id 255; surfaces (-90); filltype uni; universe 31000;} + // 1.6 % in grid + grid16Thick {type simpleCell; id 256; surfaces (-90); filltype uni; universe 16000;} + // instrumentation tube in grid + gridITThick {type simpleCell; id 257; surfaces (-90); filltype uni; universe 14000;} + // empty GT at dashpot in grid + gridGTDPThick {type simpleCell; id 258; surfaces (-90); filltype uni; universe 1010;} + // pin upper fuel plenum in grid + gridFPPThick {type simpleCell; id 259; surfaces (-90); filltype uni; universe 1008;} + // stainless steel in guide tube in grid + gridSSGThick {type simpleCell; id 260; surfaces (-90); filltype uni; universe 1023;} + // stainless steel in dash pot in grid + gridSSDPThick {type simpleCell; id 261; surfaces (-90); filltype uni; universe 1024;} + // BP plenum in grid + gridBPPThick {type simpleCell; id 262; surfaces (-90); filltype uni; universe 1012;} + // Lower rodded GT in grid + gridLRGTThick {type simpleCell; id 263; surfaces (-90); filltype uni; universe 1014;} + // Upper rodded GT in grid + gridURGTThick {type simpleCell; id 463; surfaces (-90); filltype uni; universe 1013;} + + // THIN GRID + // 2.4% in grid + grid24Thin {type simpleCell; id 264; surfaces (-92); filltype uni; universe 24000;} + // guide tube in grid + gridGTThin {type simpleCell; id 265; surfaces (-92); filltype uni; universe 12000;} + // burnable poison in grid + gridBPThin {type simpleCell; id 266; surfaces (-92); filltype uni; universe 1000;} + // 3.1% in grid + grid31Thin {type simpleCell; id 267; surfaces (-92); filltype uni; universe 31000;} + // 1.6 % in grid + grid16Thin {type simpleCell; id 268; surfaces (-92); filltype uni; universe 16000;} + // instrumentation tube in grid + gridITThin {type simpleCell; id 269; surfaces (-92); filltype uni; universe 14000;} + // empty GT at dashpot in grid + gridGTDPThin {type simpleCell; id 270; surfaces (-92); filltype uni; universe 1010;} + // pin upper fuel plenum in grid + gridFPPThin {type simpleCell; id 271; surfaces (-92); filltype uni; universe 1008;} + // stainless steel in guide tube in grid + gridSSGTThin {type simpleCell; id 272; surfaces (-92); filltype uni; universe 1023;} + // stainless steel in dash pot in grid + gridSSDPThin {type simpleCell; id 273; surfaces (-92); filltype uni; universe 1024;} + // BP plenum in grid + gridBPPThin {type simpleCell; id 274; surfaces (-92); filltype uni; universe 1012;} + // Lower rodded GT in grid (not used when rods fully withdrawn) + gridLRGTThin {type simpleCell; id 275; surfaces (-92); filltype uni; universe 1014;} + // Upper rodded GT in grid (not used when rods fully withdrawn) + gridURGTThin {type simpleCell; id 475; surfaces (-92); filltype uni; universe 1013;} + + + // 3.1% enriched pins, axial layering + 31FP460 {type simpleCell; id 100; surfaces ( 101); filltype uni; universe 1001;} + 31FP431 {type simpleCell; id 101; surfaces (-101 102); filltype uni; universe 1003;} + 31FP423 {type simpleCell; id 102; surfaces (-102 104); filltype uni; universe 1001;} + 31FP419 {type simpleCell; id 103; surfaces (-104 105); filltype uni; universe 1006;} + 31FP417 {type simpleCell; id 104; surfaces (-105 106); filltype uni; universe 1008;} + 31FP415 {type simpleCell; id 105; surfaces (-106 107); filltype uni; universe 1017;} + 31FP411 {type simpleCell; id 106; surfaces (-107 109); filltype uni; universe 1008;} + 31FP402 {type simpleCell; id 107; surfaces (-109 111); filltype uni; universe 31000;} + 31FP364 {type simpleCell; id 108; surfaces (-111 112); filltype uni; universe 9231;} + 31FP359 {type simpleCell; id 109; surfaces (-112 113); filltype uni; universe 31000;} + 31FP312 {type simpleCell; id 110; surfaces (-113 114); filltype uni; universe 9231;} + 31FP306 {type simpleCell; id 111; surfaces (-114 115); filltype uni; universe 31000;} + 31FP260 {type simpleCell; id 112; surfaces (-115 116); filltype uni; universe 9231;} + 31FP254 {type simpleCell; id 113; surfaces (-116 117); filltype uni; universe 31000;} + 31FP208 {type simpleCell; id 114; surfaces (-117 118); filltype uni; universe 9231;} + 31FP202 {type simpleCell; id 115; surfaces (-118 119); filltype uni; universe 31000;} + 31FP155 {type simpleCell; id 116; surfaces (-119 120); filltype uni; universe 9231;} + 31FP150 {type simpleCell; id 117; surfaces (-120 122); filltype uni; universe 31000;} + 31FP103 {type simpleCell; id 118; surfaces (-122 123); filltype uni; universe 9231;} + 31FP98 {type simpleCell; id 119; surfaces (-123 126); filltype uni; universe 31000;} + 31FP4052 {type simpleCell; id 120; surfaces (-126 129); filltype uni; universe 1131;} + 31FP37 {type simpleCell; id 121; surfaces (-129 130); filltype uni; universe 31000;} + 31FP36 {type simpleCell; id 122; surfaces (-130 131); filltype uni; universe 1006;} + 31FP35 {type simpleCell; id 123; surfaces (-131 132); filltype uni; universe 1003;} + 31FP20 {type simpleCell; id 124; surfaces (-132 ); filltype uni; universe 1001;} + + + //2.4% enriched pins, axial layering + 24FP460 {type simpleCell; id 125; surfaces ( 101); filltype uni; universe 1001;} + 24FP431 {type simpleCell; id 126; surfaces (-101 102); filltype uni; universe 1003;} + 24FP423 {type simpleCell; id 127; surfaces (-102 104); filltype uni; universe 1001;} + 24FP419 {type simpleCell; id 128; surfaces (-104 105); filltype uni; universe 1006;} + 24FP417 {type simpleCell; id 129; surfaces (-105 106); filltype uni; universe 1008;} + 24FP415 {type simpleCell; id 130; surfaces (-106 107); filltype uni; universe 1017;} + 24FP411 {type simpleCell; id 131; surfaces (-107 109); filltype uni; universe 1008;} + 24FP402 {type simpleCell; id 132; surfaces (-109 111); filltype uni; universe 24000;} + 24FP364 {type simpleCell; id 133; surfaces (-111 112); filltype uni; universe 9224;} + 24FP359 {type simpleCell; id 134; surfaces (-112 113); filltype uni; universe 24000;} + 24FP312 {type simpleCell; id 135; surfaces (-113 114); filltype uni; universe 9224;} + 24FP306 {type simpleCell; id 136; surfaces (-114 115); filltype uni; universe 24000;} + 24FP260 {type simpleCell; id 137; surfaces (-115 116); filltype uni; universe 9224;} + 24FP254 {type simpleCell; id 138; surfaces (-116 117); filltype uni; universe 24000;} + 24FP208 {type simpleCell; id 139; surfaces (-117 118); filltype uni; universe 9224;} + 24FP202 {type simpleCell; id 140; surfaces (-118 119); filltype uni; universe 24000;} + 24FP155 {type simpleCell; id 141; surfaces (-119 120); filltype uni; universe 9224;} + 24FP150 {type simpleCell; id 142; surfaces (-120 122); filltype uni; universe 24000;} + 24FP103 {type simpleCell; id 143; surfaces (-122 123); filltype uni; universe 9224;} + 24FP98 {type simpleCell; id 144; surfaces (-123 126); filltype uni; universe 24000;} + 24FP4052 {type simpleCell; id 145; surfaces (-126 129); filltype uni; universe 1124;} + 24FP37 {type simpleCell; id 146; surfaces (-129 130); filltype uni; universe 24000;} + 24FP36 {type simpleCell; id 147; surfaces (-130 131); filltype uni; universe 1006;} + 24FP35 {type simpleCell; id 148; surfaces (-131 132); filltype uni; universe 1003;} + 24FP20 {type simpleCell; id 149; surfaces (-132); filltype uni; universe 1001;} + + //1.6% enriched pins, axial layering + 16FP460 {type simpleCell; id 150; surfaces ( 101); filltype uni; universe 1001;} + 16FP431 {type simpleCell; id 151; surfaces (-101 102); filltype uni; universe 1003;} + 16FP423 {type simpleCell; id 152; surfaces (-102 104); filltype uni; universe 1001;} + 16FP419 {type simpleCell; id 153; surfaces (-104 105); filltype uni; universe 1006;} + 16FP417 {type simpleCell; id 154; surfaces (-105 106); filltype uni; universe 1008;} + 16FP415 {type simpleCell; id 155; surfaces (-106 107); filltype uni; universe 1017;} + 16FP411 {type simpleCell; id 156; surfaces (-107 109); filltype uni; universe 1008;} + 16FP402 {type simpleCell; id 157; surfaces (-109 111); filltype uni; universe 16000;} + 16FP364 {type simpleCell; id 158; surfaces (-111 112); filltype uni; universe 9216;} + 16FP359 {type simpleCell; id 159; surfaces (-112 113); filltype uni; universe 16000;} + 16FP312 {type simpleCell; id 160; surfaces (-113 114); filltype uni; universe 9216;} + 16FP306 {type simpleCell; id 161; surfaces (-114 115); filltype uni; universe 16000;} + 16FP260 {type simpleCell; id 162; surfaces (-115 116); filltype uni; universe 9216;} + 16FP254 {type simpleCell; id 163; surfaces (-116 117); filltype uni; universe 16000;} + 16FP208 {type simpleCell; id 164; surfaces (-117 118); filltype uni; universe 9216;} + 16FP202 {type simpleCell; id 165; surfaces (-118 119); filltype uni; universe 16000;} + 16FP155 {type simpleCell; id 166; surfaces (-119 120); filltype uni; universe 9216;} + 16FP150 {type simpleCell; id 167; surfaces (-120 122); filltype uni; universe 16000;} + 16FP103 {type simpleCell; id 168; surfaces (-122 123); filltype uni; universe 9216;} + 16FP98 {type simpleCell; id 169; surfaces (-123 126); filltype uni; universe 16000;} + 16FP4052 {type simpleCell; id 170; surfaces (-126 129); filltype uni; universe 1116;} + 16FP37 {type simpleCell; id 171; surfaces (-129 130); filltype uni; universe 16000;} + 16FP36 {type simpleCell; id 172; surfaces (-130 131); filltype uni; universe 1006;} + 16FP35 {type simpleCell; id 173; surfaces (-131 132); filltype uni; universe 1003;} + 16FP20 {type simpleCell; id 174; surfaces (-132); filltype uni; universe 1001;} + + + //guide tube, with CR, axial layering + GC460 {type simpleCell; id 175; surfaces ( 101); filltype uni; universe 1014;} + GC431 {type simpleCell; id 176; surfaces (-101 106); filltype uni; universe 1014;} // Nozzle/support plate BW? Replaced just with rod + GC415 {type simpleCell; id 177; surfaces (-106 107); filltype uni; universe 1018;} // Rodded w/ grid + GC411 {type simpleCell; id 178; surfaces (-107 109); filltype uni; universe 1014;} + GC402 {type simpleCell; id 179; surfaces (-109 134); filltype uni; universe 1015;} + GC400 {type simpleCell; id 180; surfaces (-134 111); filltype uni; universe 12000;} + GC364 {type simpleCell; id 181; surfaces (-111 112); filltype uni; universe 9112;} + GC359 {type simpleCell; id 182; surfaces (-112 113); filltype uni; universe 12000;} + GC312 {type simpleCell; id 183; surfaces (-113 114); filltype uni; universe 9112;} + GC306 {type simpleCell; id 184; surfaces (-114 115); filltype uni; universe 12000;} + GC260 {type simpleCell; id 185; surfaces (-115 116); filltype uni; universe 9112;} + GC254 {type simpleCell; id 186; surfaces (-116 117); filltype uni; universe 12000;} + GC208 {type simpleCell; id 187; surfaces (-117 118); filltype uni; universe 9112;} + GC202 {type simpleCell; id 188; surfaces (-118 119); filltype uni; universe 12000;} + GC155 {type simpleCell; id 189; surfaces (-119 120); filltype uni; universe 9112;} + GC150 {type simpleCell; id 190; surfaces (-120 122); filltype uni; universe 12000;} + GC103 {type simpleCell; id 191; surfaces (-122 123); filltype uni; universe 9112;} + GC98 {type simpleCell; id 192; surfaces (-123 126); filltype uni; universe 1010;} + GC4052 {type simpleCell; id 193; surfaces (-126 127); filltype uni; universe 12000;} + GC39 {type simpleCell; id 194; surfaces (-127 129); filltype uni; universe 1019;} + GC37 {type simpleCell; id 195; surfaces (-129 131); filltype uni; universe 1010;} + GC35 {type simpleCell; id 196; surfaces (-131 132); filltype uni; universe 1005;} + GC20 {type simpleCell; id 197; surfaces (-132); filltype uni; universe 1001;} + + + // instrumentation tube, axial layering + IT460 {type simpleCell; id 198; surfaces ( 102); filltype uni; universe 1001;} + IT423 {type simpleCell; id 199; surfaces (-102 106); filltype uni; universe 14000;} + IT415 {type simpleCell; id 200; surfaces (-106 107); filltype uni; universe 1114;} + IT411 {type simpleCell; id 201; surfaces (-107 111); filltype uni; universe 14000;} + IT364 {type simpleCell; id 202; surfaces (-111 112); filltype uni; universe 9114;} + IT359 {type simpleCell; id 203; surfaces (-112 113); filltype uni; universe 14000;} + IT312 {type simpleCell; id 204; surfaces (-113 114); filltype uni; universe 9114;} + IT306 {type simpleCell; id 205; surfaces (-114 115); filltype uni; universe 14000;} + IT260 {type simpleCell; id 206; surfaces (-115 116); filltype uni; universe 9114;} + IT254 {type simpleCell; id 207; surfaces (-116 117); filltype uni; universe 14000;} + IT208 {type simpleCell; id 208; surfaces (-117 118); filltype uni; universe 9114;} + IT202 {type simpleCell; id 209; surfaces (-118 119); filltype uni; universe 14000;} + IT155 {type simpleCell; id 210; surfaces (-119 120); filltype uni; universe 9114;} + IT150 {type simpleCell; id 211; surfaces (-120 122); filltype uni; universe 14000;} + IT103 {type simpleCell; id 212; surfaces (-122 123); filltype uni; universe 9114;} + IT98 {type simpleCell; id 213; surfaces (-123 126); filltype uni; universe 14000;} + IT4052 {type simpleCell; id 214; surfaces (-126 129); filltype uni; universe 1114;} + IT37 {type simpleCell; id 215; surfaces (-129 131); filltype uni; universe 14000;} + IT35 {type simpleCell; id 216; surfaces (-131 132); filltype uni; universe 1005;} + IT20 {type simpleCell; id 217; surfaces (-132 ); filltype uni; universe 1011;} + + + // burnable absorber, axial layering + BA460 {type simpleCell; id 218; surfaces ( 101); filltype uni; universe 1001;} + BA431 {type simpleCell; id 219; surfaces (-101 102); filltype uni; universe 1002;} + BA423 {type simpleCell; id 220; surfaces (-102 103); filltype uni; universe 1023;} + BA421 {type simpleCell; id 230; surfaces (-103 106); filltype uni; universe 1012;} + BA415 {type simpleCell; id 231; surfaces (-106 107); filltype uni; universe 1027;} + BA411 {type simpleCell; id 232; surfaces (-107 110); filltype uni; universe 1012;} + BA401 {type simpleCell; id 233; surfaces (-110 111); filltype uni; universe 1000;} + BA364 {type simpleCell; id 234; surfaces (-111 112); filltype uni; universe 1110;} + BA359 {type simpleCell; id 235; surfaces (-112 113); filltype uni; universe 1000;} + BA312 {type simpleCell; id 236; surfaces (-113 114); filltype uni; universe 1110;} + BA306 {type simpleCell; id 237; surfaces (-114 115); filltype uni; universe 1000;} + BA260 {type simpleCell; id 238; surfaces (-115 116); filltype uni; universe 1110;} + BA254 {type simpleCell; id 239; surfaces (-116 117); filltype uni; universe 1000;} + BA208 {type simpleCell; id 240; surfaces (-117 118); filltype uni; universe 1110;} + BA202 {type simpleCell; id 241; surfaces (-118 119); filltype uni; universe 1000;} + BA155 {type simpleCell; id 242; surfaces (-119 120); filltype uni; universe 1110;} + BA150 {type simpleCell; id 243; surfaces (-120 122); filltype uni; universe 1000;} + BA103 {type simpleCell; id 244; surfaces (-122 123); filltype uni; universe 1110;} + BA98 {type simpleCell; id 245; surfaces (-123 125); filltype uni; universe 1000;} + BA4055 {type simpleCell; id 246; surfaces (-125 126); filltype uni; universe 1023;} + BA4052 {type simpleCell; id 247; surfaces (-126 127); filltype uni; universe 1021;} + BA39 {type simpleCell; id 248; surfaces (-127 128); filltype uni; universe 1025;} + BA38 {type simpleCell; id 249; surfaces (-128 129); filltype uni; universe 1019;} + BA37 {type simpleCell; id 250; surfaces (-129 131); filltype uni; universe 1010;} + BA35 {type simpleCell; id 251; surfaces (-131 132); filltype uni; universe 1005;} + BA20 {type simpleCell; id 252; surfaces (-132 ); filltype uni; universe 1001;} + + //guide tube, no CR, axial layering + GT460 {type simpleCell; id 353; surfaces ( 101); filltype uni; universe 1001;} + GT431 {type simpleCell; id 354; surfaces (-101 102); filltype uni; universe 1005;} + GT423 {type simpleCell; id 355; surfaces (-102 106); filltype uni; universe 12000;} + GT415 {type simpleCell; id 356; surfaces (-106 107); filltype uni; universe 1112;} + GT411 {type simpleCell; id 357; surfaces (-107 111); filltype uni; universe 12000;} + GT364 {type simpleCell; id 358; surfaces (-111 112); filltype uni; universe 9112;} + GT359 {type simpleCell; id 359; surfaces (-112 113); filltype uni; universe 12000;} + GT312 {type simpleCell; id 360; surfaces (-113 114); filltype uni; universe 9112;} + GT306 {type simpleCell; id 361; surfaces (-114 115); filltype uni; universe 12000;} + GT260 {type simpleCell; id 362; surfaces (-115 116); filltype uni; universe 9112;} + GT254 {type simpleCell; id 363; surfaces (-116 117); filltype uni; universe 12000;} + GT208 {type simpleCell; id 364; surfaces (-117 118); filltype uni; universe 9112;} + GT202 {type simpleCell; id 365; surfaces (-118 119); filltype uni; universe 12000;} + GT155 {type simpleCell; id 366; surfaces (-119 120); filltype uni; universe 9112;} + GT150 {type simpleCell; id 367; surfaces (-120 122); filltype uni; universe 12000;} + GT103 {type simpleCell; id 368; surfaces (-122 123); filltype uni; universe 9112;} + GT98 {type simpleCell; id 369; surfaces (-123 126); filltype uni; universe 1010;} + GT4052 {type simpleCell; id 370; surfaces (-126 127); filltype uni; universe 12000;} + GT39 {type simpleCell; id 371; surfaces (-127 129); filltype uni; universe 1019;} + GT37 {type simpleCell; id 372; surfaces (-129 131); filltype uni; universe 1010;} + GT35 {type simpleCell; id 373; surfaces (-131 132); filltype uni; universe 1005;} + GT20 {type simpleCell; id 374; surfaces (-132); filltype uni; universe 1001;} + + //guide tube, partially inserted CR, axial layering + // Note: universe is basically a combo of fully and partially inserted rods + GP460 {type simpleCell; id 375; surfaces ( 101); filltype uni; universe 1013;} + GP431 {type simpleCell; id 376; surfaces (-101 102); filltype uni; universe 1013;} + GP423 {type simpleCell; id 377; surfaces (-102 106); filltype uni; universe 1013;} + GP415 {type simpleCell; id 378; surfaces (-106 107); filltype uni; universe 1133;} + GP411 {type simpleCell; id 379; surfaces (-107 111); filltype uni; universe 1013;} + GP364 {type simpleCell; id 380; surfaces (-111 112); filltype uni; universe 9233;} + GP359 {type simpleCell; id 381; surfaces (-112 137); filltype uni; universe 1013;} + GP322 {type simpleCell; id 382; surfaces (-137 113); filltype uni; universe 1014;} + GP312 {type simpleCell; id 383; surfaces (-113 114); filltype uni; universe 9232;} + GP306 {type simpleCell; id 384; surfaces (-114 115); filltype uni; universe 1014;} + GP260 {type simpleCell; id 385; surfaces (-115 116); filltype uni; universe 9232;} + GP254 {type simpleCell; id 386; surfaces (-116 138); filltype uni; universe 1014;} + GP220 {type simpleCell; id 387; surfaces (-138 139); filltype uni; universe 1023;} + GP219 {type simpleCell; id 388; surfaces (-139 117); filltype uni; universe 12000;} + GP208 {type simpleCell; id 389; surfaces (-117 118); filltype uni; universe 9112;} + GP202 {type simpleCell; id 390; surfaces (-118 119); filltype uni; universe 12000;} + GP155 {type simpleCell; id 391; surfaces (-119 120); filltype uni; universe 9112;} + GP150 {type simpleCell; id 392; surfaces (-120 122); filltype uni; universe 12000;} + GP103 {type simpleCell; id 393; surfaces (-122 123); filltype uni; universe 9112;} + GP98 {type simpleCell; id 394; surfaces (-123 126); filltype uni; universe 1010;} + GP4052 {type simpleCell; id 395; surfaces (-126 127); filltype uni; universe 12000;} + GP39 {type simpleCell; id 396; surfaces (-127 129); filltype uni; universe 1019;} + GP37 {type simpleCell; id 397; surfaces (-129 131); filltype uni; universe 1010;} + GP35 {type simpleCell; id 398; surfaces (-131 132); filltype uni; universe 1005;} + GP20 {type simpleCell; id 399; surfaces (-132); filltype uni; universe 1001;} + + // control rod, axial layering + // Used (probably with some modification) only when fully inserted + //CR460 {type simpleCell; id 448; surfaces (-120 121); filltype uni; universe 1002;} + //CR415 {type simpleCell; id 449; surfaces (-126 403); filltype uni; universe 1013;} + //CR403 {type simpleCell; id 450; surfaces (-403 402); filltype uni; universe 1015;} + //CR402 {type simpleCell; id 451; surfaces (-402 1430); filltype uni; universe 1013;} + //CR143 {type simpleCell; id 452; surfaces (-1430 41); filltype uni; universe 1014;} + //CR41 {type simpleCell; id 453; surfaces (-41 143); filltype uni; universe 1002;} + //CR39 {type simpleCell; id 454; surfaces (-143 147); filltype uni; universe 1001;} + //CR35 {type simpleCell; id 455; surfaces (-147 148); filltype uni; universe 1005;} + //CR20 {type simpleCell; id 456; surfaces (-148 149); filltype uni; universe 1001;} + + + + outsideLeftBaffle { type simpleCell; id 52; surfaces (-50); filltype mat; material Water;} + leftBaffle { type simpleCell; id 53; surfaces (50 -51); filltype mat; material SS304;} + insideLeftBaffle { type simpleCell; id 54; surfaces (51); filltype mat; material Water;} + + outsideRightBaffle { type simpleCell; id 55; surfaces (-52); filltype mat; material Water;} + RightBaffle { type simpleCell; id 56; surfaces (52 -53); filltype mat; material SS304;} + insideRightBaffle { type simpleCell; id 57; surfaces (53); filltype mat; material Water;} + + outsideTopBaffle { type simpleCell; id 58; surfaces (-54); filltype mat; material Water;} + TopBaffle { type simpleCell; id 59; surfaces (54 -55); filltype mat; material SS304;} + insideTopBaffle { type simpleCell; id 60; surfaces (55); filltype mat; material Water;} + + outsideBottomBaffle { type simpleCell; id 61; surfaces (-56); filltype mat; material Water;} + BottomBaffle { type simpleCell; id 62; surfaces (56 -57); filltype mat; material SS304;} + insideBottomBaffle { type simpleCell; id 63; surfaces (57); filltype mat; material Water;} + + topLeftCornerBaffle1 { type simpleCell; id 64; surfaces (52 -53 -57); filltype mat; material SS304;} + topLeftCornerBaffle2 { type simpleCell; id 65; surfaces (56 -57 -52); filltype mat; material SS304;} + topLeftCornerGap1 { type simpleCell; id 66; surfaces (57); filltype mat; material Water;} + topLeftCornerGap2 { type simpleCell; id 67; surfaces (53); filltype mat; material Water;} + topLeftMajorGap { type simpleCell; id 68; surfaces (-56 -52); filltype mat; material Water;} + + topRightCornerBaffle1 { type simpleCell; id 69; surfaces (-57 50 -51); filltype mat; material SS304;} + topRightCornerBaffle2 { type simpleCell; id 70; surfaces (-50 56 -57); filltype mat; material SS304;} + topRightCornerGap1 { type simpleCell; id 71; surfaces (57); filltype mat; material Water;} + topRightCornerGap2 { type simpleCell; id 72; surfaces (51); filltype mat; material Water;} + topRightMajorGap { type simpleCell; id 73; surfaces (-56 -50); filltype mat; material Water;} + + bottomLeftCornerBaffle1 { type simpleCell; id 74; surfaces (-55 52 -53); filltype mat; material SS304;} + bottomLeftCornerBaffle2 { type simpleCell; id 75; surfaces (-55 54 -52); filltype mat; material SS304;} + bottomLeftCornerGap1 { type simpleCell; id 76; surfaces (55); filltype mat; material Water;} + bottomLeftCornerGap2 { type simpleCell; id 77; surfaces (53); filltype mat; material Water;} + bottomLeftMajorGap { type simpleCell; id 78; surfaces (-54 -52); filltype mat; material Water;} + + bottomRightCornerBaffle1 { type simpleCell; id 79; surfaces (-51 50 -55); filltype mat; material SS304;} + bottomRightCornerBaffle2 { type simpleCell; id 80; surfaces (-55 54 -50); filltype mat; material SS304;} + bottomRightCornerGap1 { type simpleCell; id 81; surfaces (51); filltype mat; material Water;} + bottomRightCornerGap2 { type simpleCell; id 82; surfaces (55); filltype mat; material Water;} + bottomRightMajorGap { type simpleCell; id 83; surfaces (-50 -54); filltype mat; material Water;} + + + TLSG1 { type simpleCell; id 84; surfaces (-56 -52); filltype mat; material Water;} + TLSG2 { type simpleCell; id 85; surfaces (-56 52); filltype mat; material Water;} + TLSG3 { type simpleCell; id 86; surfaces (56 -52); filltype mat; material Water;} + topLeftSquare { type simpleCell; id 87; surfaces (56 52); filltype mat; material SS304;} + + TRSG1 { type simpleCell; id 88; surfaces (-56 50); filltype mat; material Water;} + TRSG2 { type simpleCell; id 89; surfaces (-56 -50); filltype mat; material Water;} + TRSG3 { type simpleCell; id 90; surfaces (56 -50); filltype mat; material Water;} + topRightSquare { type simpleCell; id 91; surfaces (56 50); filltype mat; material SS304;} + + BLSG1 { type simpleCell; id 92; surfaces (54 -52); filltype mat; material Water;} + BLSG2 { type simpleCell; id 93; surfaces (-54 52); filltype mat; material Water;} + BLSG3 { type simpleCell; id 94; surfaces (-54 -52); filltype mat; material Water;} + bottomLeftSquare { type simpleCell; id 95; surfaces (54 52); filltype mat; material SS304;} + + BRSG1 { type simpleCell; id 96; surfaces (-54 50); filltype mat; material Water;} + BRSG2 { type simpleCell; id 97; surfaces (54 -50); filltype mat; material Water;} + BRSG3 { type simpleCell; id 98; surfaces (-54 -50); filltype mat; material Water;} + bottomRightSquare { type simpleCell; id 99; surfaces (54 50); filltype mat; material SS304;} + } + + universes { + root { id 1; type rootUniverse; border 1; fill u<8888>; } + + // Pin universes + + //Burnable poison + pinBPaboveDP { id 1000; type pinUniverse; radii (0.21400 0.23051 0.24130 0.42672 0.43688 0.48387 0.56134 0.60198 0.0); + fills (Air SS304 Helium BorosilicateGlass Helium SS304 Water Zircaloy Water);} + pinBPPlenumGeometry { id 1012; type pinUniverse; radii ( 0.21400 0.23051 0.43688 0.48387 0.50419 0.54610 0.0); + fills (Air SS304 Helium SS304 Water Zircaloy Water);} + + //guide tubes + pinGTaboveDP { id 12000; type pinUniverse; radii (0.56134 0.60198 0.0 ); fills (Water Zircaloy Water);} + pinGTatDP { id 1010; type pinUniverse; radii (0.50419 0.54610 0.0); fills (Water Zircaloy Water);} + + //INST Tube + pinIT { id 14000; type pinUniverse; radii (0.43688 0.48387 0.56134 0.60198 0.0 ); + fills (Air Zircaloy Water Zircaloy Water);} + pinBareInstrumentThimble { id 1011; type pinUniverse; radii (0.43688 0.48387 0.0); fills (Air Zircaloy Water);} + + // Fuel pins + pin16 { id 16000; type pinUniverse; radii (0.39218 0.40005 0.45720 0.0); + fills (UO2-16 Helium Zircaloy Water);} + pin24 { id 24000; type pinUniverse; radii (0.39218 0.40005 0.45720 0.0); + fills (UO2-24 Helium Zircaloy Water);} + pin31 { id 31000; type pinUniverse; radii (0.39218 0.40005 0.45720 0.0); + fills (UO2-31 Helium Zircaloy Water);} + // Higher enrichments not used + //pin32 { id 32000; type pinUniverse; radii (0.39218 0.40005 0.45720 0.0); + // fills (UO2-32 Helium Zircaloy Water);} + //pin34 { id 34000; type pinUniverse; radii (0.39218 0.40005 0.45720 0.0); + // fills (UO2-34 Helium Zircaloy Water);} + + pinWater { id 1001; type pinUniverse; radii ( 0.0); fills (Water);} + + + // Solid pins, assumed radius to be that of a fuel pin (0.45720) + pinNozzle_SupportSteel { id 1003; type pinUniverse; radii ( 0.45720 0.0); fills (SupportPlateSS Water);} + pinSupportPlateBW { id 1005; type pinUniverse; radii ( 0.45720 0.0); fills (SupportPlateBW Water);} + pinZircaloy { id 1006; type pinUniverse; radii ( 0.45720 0.0); fills (Zircaloy Water);} + + + SSinDashPot { id 1024; type pinUniverse; radii (0.50419 0.54610 0.0); fills (SS304 Zircaloy Water);} + SSinGuideTube { id 1023; type pinUniverse; radii ( 0.56134 0.60198 0.0); fills (SS304 Zircaloy Water);} + SSnoGuideTube { id 1002; type pinUniverse; radii ( 0.56134 0.0); fills (SS304 Water);} + + + pinUpperFuelPlenum { id 1008; type pinUniverse; radii ( 0.06459 0.40005 0.45720 0.0); + fills (Inconel Helium Zircaloy Water);} + + // Control rod pins + pinControlRodUpper { id 1013; type pinUniverse; radii ( 0.37338 0.38608 0.48387 0.56134 0.60198 0.0); + fills (B4C Helium SS304 Water Zircaloy Water);} + pinControlRodLower { id 1014; type pinUniverse; radii ( 0.38227 0.38608 0.48387 0.56134 0.60198 0.0); + fills (Ag-In-Cd Helium SS304 Water Zircaloy Water);} + pinControlRodSpacer { id 1015; type pinUniverse; radii ( 0.37845 0.38608 0.48387 0.56134 0.60198 0.0); + fills (SS304 Helium SS304 Water Zircaloy Water);} + pinControlRodPlenum { id 1016; type pinUniverse; radii ( 0.06459 0.38608 0.48387 0.56134 0.60198 0.0); + fills (Inconel Helium SS304 Water Zircaloy Water);} + + // pins that have grids + fuelRodPlenumWithGridThick { + id 1017; + type cellUniverse; + cells ( 259 555);} + + GTRodThick { + id 1018; + type cellUniverse; + cells (263 555);} + + dashPotGuideTubeGridThick { + id 1019; + type cellUniverse; + cells ( 258 555);} + + dashPotGuideTubeGridThin { + id 1020; + type cellUniverse; + cells ( 270 556);} + + SSinGuideTubeThick { + id 1021; + type cellUniverse; + cells ( 260 555);} + + SSinGuideTubeThin { + id 1022; + type cellUniverse; + cells ( 272 556);} + + SSinDashPotThick { + id 1025; + type cellUniverse; + cells ( 261 555);} + + SSinDashPotThin { + id 1026; + type cellUniverse; + cells ( 273 556);} + + BPPlenumThick { + id 1027; + type cellUniverse; + cells ( 262 555);} + + BPPlenumThin { + id 1028; + type cellUniverse; + cells ( 274 556);} + + BPaboveDPThin { + id 1110; + type cellUniverse; + cells (266 556);} + + GTThick { + id 1112; + type cellUniverse; + cells (254 555);} + + ITThick { + id 1114; + type cellUniverse; + cells (257 555);} + + pin16Thick { + id 1116; + type cellUniverse; + cells (256 555);} + + pin24Thick { + id 1124; + type cellUniverse; + cells (253 555);} + + pin31Thick { + id 1131; + type cellUniverse; + cells (255 555);} + + LowerRodGTThick { + id 1132; + type cellUniverse; + cells (263 555);} + + UpperRodGTThick { + id 1133; + type cellUniverse; + cells (463 555);} + + BPThin { // Is this necessary given 1110, BPaboveDPThin??? + id 9110; + type cellUniverse; + cells (266 556);} + + GTThin { + id 9112; + type cellUniverse; + cells (265 556);} + + ITThin { + id 9114; + type cellUniverse; + cells (269 556);} + + pin16Thin { + id 9216; + type cellUniverse; + cells (268 556);} + + pin24Thin { + id 9224; + type cellUniverse; + cells (264 556);} + + pin31Thin { + id 9231; + type cellUniverse; + cells (267 556);} + + LowerRodGTThin { + id 9232; + type cellUniverse; + cells (275 556);} + + UpperRodGTThin { + id 9233; + type cellUniverse; + cells (475 556);} + + // Axial stacks of universes to make up full pins + + // 3.1 % + fuelPin31 { + id 31; + type cellUniverse; + cells ( 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124);} + + // 2.4 % + fuelPin24 { + id 24; + type cellUniverse; + cells ( 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149);} + + //1.6 % + fuelPin16 { + id 16; + type cellUniverse; + cells ( 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174);} + + //burnable absorber + BP { + id 10; + type cellUniverse; + cells (218 219 220 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252);} + + // guide tube, with CR + GuideTubeRodded { + id 12; + type cellUniverse; + cells (175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197);} + + //control rod, fully retracted + GuideTubeEmpty { + id 13; + type cellUniverse; + cells (353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374);} + + //instr. tube + instrumentTube { + id 14; + type cellUniverse; + cells (198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217);} + + + // control rod, partially inserted + GuideTubePartial { + id 15; + type cellUniverse; + cells (375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399); + + } + + + + // Lattices w/o grid + // Names represent AE + A0E24 { + id 1424; + type latUniverse; + origin (0.0 0.0 0.0); + pitch (1.26 1.26 0.0); + shape (17 17 0); + padMat Water; + map ( + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 24 24 24 12 24 24 12 24 24 12 24 24 24 24 24 + 24 24 24 12 24 24 24 24 24 24 24 24 24 12 24 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 12 24 24 12 24 24 12 24 24 12 24 24 12 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 12 24 24 12 24 24 14 24 24 12 24 24 12 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 12 24 24 12 24 24 12 24 24 12 24 24 12 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 24 12 24 24 24 24 24 24 24 24 24 12 24 24 24 + 24 24 24 24 24 12 24 24 12 24 24 12 24 24 24 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 ); } + + // assembly with sleeves at different heights + A0E24Sleeve { + id 14240; + type cellUniverse; + cells ( + 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 + 2016 2017 2018 2019 + );} + + + A0E16 { + id 1416; + type latUniverse; + origin (0.0 0.0 0.0); + pitch (1.26 1.26 0.0); + shape (17 17 0); + padMat Water; + map ( + 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 + 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 + 16 16 16 16 16 12 16 16 12 16 16 12 16 16 16 16 16 + 16 16 16 12 16 16 16 16 16 16 16 16 16 12 16 16 16 + 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 + 16 16 12 16 16 12 16 16 12 16 16 12 16 16 12 16 16 + 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 + 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 + 16 16 12 16 16 12 16 16 14 16 16 12 16 16 12 16 16 + 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 + 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 + 16 16 12 16 16 12 16 16 12 16 16 12 16 16 12 16 16 + 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 + 16 16 16 12 16 16 16 16 16 16 16 16 16 12 16 16 16 + 16 16 16 16 16 12 16 16 12 16 16 12 16 16 16 16 16 + 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 + 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16); } + + A0E16Sleeve { + id 14160; + type cellUniverse; + cells ( + 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 + 2016 2017 2018 2020 + ); + } + + A0E31 { + id 1431; + type latUniverse; + origin (0.0 0.0 0.0); + pitch (1.26 1.26 0.0); + shape (17 17 0); + padMat Water; + map ( + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 12 31 31 12 31 31 12 31 31 31 31 31 + 31 31 31 12 31 31 31 31 31 31 31 31 31 12 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 12 31 31 12 31 31 12 31 31 12 31 31 12 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 12 31 31 12 31 31 14 31 31 12 31 31 12 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 12 31 31 12 31 31 12 31 31 12 31 31 12 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 12 31 31 31 31 31 31 31 31 31 12 31 31 31 + 31 31 31 31 31 12 31 31 12 31 31 12 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 ); } + + A0E31Sleeve { + id 14310; + type cellUniverse; + cells ( + 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 + 2016 2017 2018 2021 + ); + } + + + A6BE31B { + id 60316; + type latUniverse; + origin (0.0 0.0 0.0); + pitch (1.26 1.26 0.0); + shape (17 17 0); + padMat Water; + map ( + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 10 31 31 12 31 31 10 31 31 31 31 31 + 31 31 31 10 31 31 31 31 31 31 31 31 31 10 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 10 31 31 12 31 31 12 31 31 12 31 31 10 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 12 31 31 12 31 31 14 31 31 12 31 31 12 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 12 31 31 12 31 31 12 31 31 12 31 31 12 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 12 31 31 31 31 31 31 31 31 31 12 31 31 31 + 31 31 31 31 31 12 31 31 12 31 31 12 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 ); } + + + A6BE31BSleeve { + id 603160; + type cellUniverse; + cells ( + 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 + 2016 2017 2018 2022 + ); + } + + A6BE31T { + id 603112; + type latUniverse; + origin (0.0 0.0 0.0); + pitch (1.26 1.26 0.0); + shape (17 17 0); + padMat Water; + map ( + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 12 31 31 12 31 31 12 31 31 31 31 31 + 31 31 31 12 31 31 31 31 31 31 31 31 31 12 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 12 31 31 12 31 31 12 31 31 12 31 31 12 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 12 31 31 12 31 31 14 31 31 12 31 31 12 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 10 31 31 12 31 31 12 31 31 12 31 31 10 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 10 31 31 31 31 31 31 31 31 31 10 31 31 31 + 31 31 31 31 31 10 31 31 12 31 31 10 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 ); } + + A6BE31TSleeve { + id 6031120; + type cellUniverse; + cells ( + 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 + 2016 2017 2018 2023 + ); + } + + A6BE31R { + id 60313; + type latUniverse; + origin (0.0 0.0 0.0); + pitch (1.26 1.26 0.0); + shape (17 17 0); + padMat Water; + map ( + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 10 31 31 12 31 31 12 31 31 31 31 31 + 31 31 31 10 31 31 31 31 31 31 31 31 31 12 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 10 31 31 12 31 31 12 31 31 12 31 31 12 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 12 31 31 12 31 31 14 31 31 12 31 31 12 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 10 31 31 12 31 31 12 31 31 12 31 31 12 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 10 31 31 31 31 31 31 31 31 31 12 31 31 31 + 31 31 31 31 31 10 31 31 12 31 31 12 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 ); } + + A6BE31RSleeve { + id 603130; + type cellUniverse; + cells ( + 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 + 2016 2017 2018 2024 + ); + } + + A6BE31L { + id 60319; + type latUniverse; + origin (0.0 0.0 0.0); + pitch (1.26 1.26 0.0); + shape (17 17 0); + padMat Water; + map ( + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 12 31 31 12 31 31 10 31 31 31 31 31 + 31 31 31 12 31 31 31 31 31 31 31 31 31 10 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 12 31 31 12 31 31 12 31 31 12 31 31 10 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 12 31 31 12 31 31 14 31 31 12 31 31 12 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 12 31 31 12 31 31 12 31 31 12 31 31 10 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 12 31 31 31 31 31 31 31 31 31 10 31 31 31 + 31 31 31 31 31 12 31 31 12 31 31 10 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 ); } + + A6BE31LSleeve { + id 603190; + type cellUniverse; + cells ( + 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 + 2016 2017 2018 2025 + ); + } + + A15BE31BR { + id 15315; + type latUniverse; + origin (0.0 0.0 0.0); + pitch (1.26 1.26 0.0); + shape (17 17 0); + padMat Water; + map ( + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 10 31 31 10 31 31 10 31 31 31 31 31 + 31 31 31 10 31 31 31 31 31 31 31 31 31 12 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 10 31 31 10 31 31 10 31 31 10 31 31 12 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 10 31 31 10 31 31 14 31 31 10 31 31 12 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 10 31 31 10 31 31 10 31 31 10 31 31 12 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 12 31 31 31 31 31 31 31 31 31 12 31 31 31 + 31 31 31 31 31 12 31 31 12 31 31 12 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 ); } + + A15BE31BRSleeve { + id 153150; + type cellUniverse; + cells ( + 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 + 2016 2017 2018 2026 + ); + } + + A15BE31BL { + id 15317; + type latUniverse; + origin (0.0 0.0 0.0); + pitch (1.26 1.26 0.0); + shape (17 17 0); + padMat Water; + map ( + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 10 31 31 10 31 31 10 31 31 31 31 31 + 31 31 31 12 31 31 31 31 31 31 31 31 31 10 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 12 31 31 10 31 31 10 31 31 10 31 31 10 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 12 31 31 10 31 31 14 31 31 10 31 31 10 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 12 31 31 10 31 31 10 31 31 10 31 31 10 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 12 31 31 31 31 31 31 31 31 31 12 31 31 31 + 31 31 31 31 31 12 31 31 12 31 31 12 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 ); } + + A15BE31BLSleeve { + id 153170; + type cellUniverse; + cells ( + 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 + 2016 2017 2018 2027 + ); + } + + A15BE31TR { + id 15311; + type latUniverse; + origin (0.0 0.0 0.0); + pitch (1.26 1.26 0.0); + shape (17 17 0); + padMat Water; + map ( + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 12 31 31 12 31 31 12 31 31 31 31 31 + 31 31 31 12 31 31 31 31 31 31 31 31 31 12 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 10 31 31 10 31 31 10 31 31 10 31 31 12 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 10 31 31 10 31 31 14 31 31 10 31 31 12 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 10 31 31 10 31 31 10 31 31 10 31 31 12 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 10 31 31 31 31 31 31 31 31 31 12 31 31 31 + 31 31 31 31 31 10 31 31 10 31 31 10 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 ); } + + A15BE31TRSleeve { + id 153110; + type cellUniverse; + cells ( + 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 + 2016 2017 2018 2028 + ); + } + + A15BE31TL { + id 153111; + type latUniverse; + origin (0.0 0.0 0.0); + pitch (1.26 1.26 0.0); + shape (17 17 0); + padMat Water; + map ( + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 12 31 31 12 31 31 12 31 31 31 31 31 + 31 31 31 12 31 31 31 31 31 31 31 31 31 12 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 12 31 31 10 31 31 10 31 31 10 31 31 10 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 12 31 31 10 31 31 14 31 31 10 31 31 10 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 12 31 31 10 31 31 10 31 31 10 31 31 10 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 12 31 31 31 31 31 31 31 31 31 10 31 31 31 + 31 31 31 31 31 10 31 31 10 31 31 10 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 ); } + + A15BE31TLSleeve { + id 1531110; + type cellUniverse; + cells ( + 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 + 2016 2017 2018 2029 + ); + } + + A16BE31 { + id 1631; + type latUniverse; + origin (0.0 0.0 0.0); + pitch (1.26 1.26 0.0); + shape (17 17 0); + padMat Water; + map ( + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 10 31 31 10 31 31 10 31 31 31 31 31 + 31 31 31 10 31 31 31 31 31 31 31 31 31 10 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 10 31 31 12 31 31 12 31 31 12 31 31 10 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 10 31 31 12 31 31 14 31 31 12 31 31 10 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 10 31 31 12 31 31 12 31 31 12 31 31 10 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 10 31 31 31 31 31 31 31 31 31 10 31 31 31 + 31 31 31 31 31 10 31 31 10 31 31 10 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 ); } + + A16BE31Sleeve { + id 16310; + type cellUniverse; + cells ( + 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 + 2016 2017 2018 2030 + ); + } + + A20BE31 { + id 2031; + type latUniverse; + origin (0.0 0.0 0.0); + pitch (1.26 1.26 0.0); + shape (17 17 0); + padMat Water; + map ( + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 10 31 31 10 31 31 10 31 31 31 31 31 + 31 31 31 10 31 31 31 31 31 31 31 31 31 10 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 10 31 31 10 31 31 12 31 31 10 31 31 10 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 10 31 31 12 31 31 14 31 31 12 31 31 10 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 10 31 31 10 31 31 12 31 31 10 31 31 10 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 10 31 31 31 31 31 31 31 31 31 10 31 31 31 + 31 31 31 31 31 10 31 31 10 31 31 10 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 ); } + + A20BE31Sleeve { + id 20310; + type cellUniverse; + cells ( + 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 + 2016 2017 2018 2031 + ); + } + + A12BE24 { + id 1224; + type latUniverse; + origin (0.0 0.0 0.0); + pitch (1.26 1.26 0.0); + shape (17 17 0); + padMat Water; + map ( + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 24 24 24 10 24 24 12 24 24 10 24 24 24 24 24 + 24 24 24 10 24 24 24 24 24 24 24 24 24 10 24 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 10 24 24 12 24 24 12 24 24 12 24 24 10 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 12 24 24 12 24 24 14 24 24 12 24 24 12 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 10 24 24 12 24 24 12 24 24 12 24 24 10 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 24 10 24 24 24 24 24 24 24 24 24 10 24 24 24 + 24 24 24 24 24 10 24 24 12 24 24 10 24 24 24 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 ); } + + A12BE24Sleeve { + id 12240; + type cellUniverse; + cells ( + 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 + 2016 2017 2018 2032 + ); + } + + A16BE24 { + id 1624; + type latUniverse; + origin (0.0 0.0 0.0); + pitch (1.26 1.26 0.0); + shape (17 17 0); + padMat Water; + map ( + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 24 24 24 10 24 24 10 24 24 10 24 24 24 24 24 + 24 24 24 10 24 24 24 24 24 24 24 24 24 10 24 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 10 24 24 12 24 24 12 24 24 12 24 24 10 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 10 24 24 12 24 24 14 24 24 12 24 24 10 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 10 24 24 12 24 24 12 24 24 12 24 24 10 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 24 10 24 24 24 24 24 24 24 24 24 10 24 24 24 + 24 24 24 24 24 10 24 24 10 24 24 10 24 24 24 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 ); } + + A16BE24Sleeve { + id 16240; + type cellUniverse; + cells ( + 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 + 2016 2017 2018 2033 + ); + } + + // Unrodded assemblies + A0E24U { + id 2424; + type latUniverse; + origin (0.0 0.0 0.0); + pitch (1.26 1.26 0.0); + shape (17 17 0); + padMat Water; + map ( + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 24 24 24 13 24 24 13 24 24 13 24 24 24 24 24 + 24 24 24 13 24 24 24 24 24 24 24 24 24 13 24 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 13 24 24 13 24 24 13 24 24 13 24 24 13 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 13 24 24 13 24 24 14 24 24 13 24 24 13 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 13 24 24 13 24 24 13 24 24 13 24 24 13 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 24 13 24 24 24 24 24 24 24 24 24 13 24 24 24 + 24 24 24 24 24 13 24 24 13 24 24 13 24 24 24 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 ); } + + A0E24USleeve { + id 24240; + type cellUniverse; + cells ( + 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 + 2016 2017 2018 2034 + );} + + A0E16U { + id 2416; + type latUniverse; + origin (0.0 0.0 0.0); + pitch (1.26 1.26 0.0); + shape (17 17 0); + padMat Water; + map ( + 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 + 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 + 16 16 16 16 16 13 16 16 13 16 16 13 16 16 16 16 16 + 16 16 16 13 16 16 16 16 16 16 16 16 16 13 16 16 16 + 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 + 16 16 13 16 16 13 16 16 13 16 16 13 16 16 13 16 16 + 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 + 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 + 16 16 13 16 16 13 16 16 14 16 16 13 16 16 13 16 16 + 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 + 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 + 16 16 13 16 16 13 16 16 13 16 16 13 16 16 13 16 16 + 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 + 16 16 16 13 16 16 16 16 16 16 16 16 16 13 16 16 16 + 16 16 16 16 16 13 16 16 13 16 16 13 16 16 16 16 16 + 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 + 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16); } + + // sleeved + A0E16USleeve { + id 24160; + type cellUniverse; + cells ( + 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 + 2016 2017 2018 2035 + ); + } + + A0E31U { + id 2431; + type latUniverse; + origin (0.0 0.0 0.0); + pitch (1.26 1.26 0.0); + shape (17 17 0); + padMat Water; + map ( + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 13 31 31 13 31 31 13 31 31 31 31 31 + 31 31 31 13 31 31 31 31 31 31 31 31 31 13 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 13 31 31 13 31 31 13 31 31 13 31 31 13 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 13 31 31 13 31 31 14 31 31 13 31 31 13 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 13 31 31 13 31 31 13 31 31 13 31 31 13 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 13 31 31 31 31 31 31 31 31 31 13 31 31 31 + 31 31 31 31 31 13 31 31 13 31 31 13 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 ); } + + A0E31USleeve { + id 24310; + type cellUniverse; + cells ( + 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 + 2016 2017 2018 2036 + ); + } + + + A6BE31BU { + id 70316; + type latUniverse; + origin (0.0 0.0 0.0); + pitch (1.26 1.26 0.0); + shape (17 17 0); + padMat Water; + map ( + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 10 31 31 13 31 31 10 31 31 31 31 31 + 31 31 31 10 31 31 31 31 31 31 31 31 31 10 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 10 31 31 13 31 31 13 31 31 13 31 31 10 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 13 31 31 13 31 31 14 31 31 13 31 31 13 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 13 31 31 13 31 31 13 31 31 13 31 31 13 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 13 31 31 31 31 31 31 31 31 31 13 31 31 31 + 31 31 31 31 31 13 31 31 13 31 31 13 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 ); } + + + A6BE31BUSleeve { + id 703160; + type cellUniverse; + cells ( + 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 + 2016 2017 2018 2037 + ); + } + + A6BE31TU { + id 703112; + type latUniverse; + origin (0.0 0.0 0.0); + pitch (1.26 1.26 0.0); + shape (17 17 0); + padMat Water; + map ( + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 13 31 31 13 31 31 13 31 31 31 31 31 + 31 31 31 13 31 31 31 31 31 31 31 31 31 13 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 13 31 31 13 31 31 13 31 31 13 31 31 13 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 13 31 31 13 31 31 14 31 31 13 31 31 13 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 10 31 31 13 31 31 13 31 31 13 31 31 10 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 10 31 31 31 31 31 31 31 31 31 10 31 31 31 + 31 31 31 31 31 10 31 31 13 31 31 10 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 ); } + + A6BE31TUSleeve { + id 7031120; + type cellUniverse; + cells ( + 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 + 2016 2017 2018 2038 + ); + } + + A6BE31RU { + id 70313; + type latUniverse; + origin (0.0 0.0 0.0); + pitch (1.26 1.26 0.0); + shape (17 17 0); + padMat Water; + map ( + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 10 31 31 13 31 31 13 31 31 31 31 31 + 31 31 31 10 31 31 31 31 31 31 31 31 31 13 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 10 31 31 13 31 31 13 31 31 13 31 31 13 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 13 31 31 13 31 31 14 31 31 13 31 31 13 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 10 31 31 13 31 31 13 31 31 13 31 31 13 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 10 31 31 31 31 31 31 31 31 31 13 31 31 31 + 31 31 31 31 31 10 31 31 13 31 31 13 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 ); } + + A6BE31RUSleeve { + id 703130; + type cellUniverse; + cells ( + 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 + 2016 2017 2018 2039 + ); + } + + A6BE31LU { + id 70319; + type latUniverse; + origin (0.0 0.0 0.0); + pitch (1.26 1.26 0.0); + shape (17 17 0); + padMat Water; + map ( + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 13 31 31 13 31 31 10 31 31 31 31 31 + 31 31 31 13 31 31 31 31 31 31 31 31 31 10 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 13 31 31 13 31 31 13 31 31 13 31 31 10 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 13 31 31 13 31 31 14 31 31 13 31 31 13 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 13 31 31 13 31 31 13 31 31 13 31 31 10 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 13 31 31 31 31 31 31 31 31 31 10 31 31 31 + 31 31 31 31 31 13 31 31 13 31 31 10 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 ); } + + A6BE31LUSleeve { + id 703190; + type cellUniverse; + cells ( + 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 + 2016 2017 2018 2040 + ); + } + + A15BE31BRU { + id 25315; + type latUniverse; + origin (0.0 0.0 0.0); + pitch (1.26 1.26 0.0); + shape (17 17 0); + padMat Water; + map ( + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 10 31 31 10 31 31 10 31 31 31 31 31 + 31 31 31 10 31 31 31 31 31 31 31 31 31 13 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 10 31 31 10 31 31 10 31 31 10 31 31 13 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 10 31 31 10 31 31 14 31 31 10 31 31 13 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 10 31 31 10 31 31 10 31 31 10 31 31 13 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 13 31 31 31 31 31 31 31 31 31 13 31 31 31 + 31 31 31 31 31 13 31 31 13 31 31 13 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 ); } + + A15BE31BRUSleeve { + id 253150; + type cellUniverse; + cells ( + 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 + 2016 2017 2018 2041 + ); + } + + A15BE31BLU { + id 25317; + type latUniverse; + origin (0.0 0.0 0.0); + pitch (1.26 1.26 0.0); + shape (17 17 0); + padMat Water; + map ( + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 10 31 31 10 31 31 10 31 31 31 31 31 + 31 31 31 13 31 31 31 31 31 31 31 31 31 10 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 13 31 31 10 31 31 10 31 31 10 31 31 10 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 13 31 31 10 31 31 14 31 31 10 31 31 10 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 13 31 31 10 31 31 10 31 31 10 31 31 10 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 13 31 31 31 31 31 31 31 31 31 13 31 31 31 + 31 31 31 31 31 13 31 31 13 31 31 13 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 ); } + + A15BE31BLUSleeve { + id 253170; + type cellUniverse; + cells ( + 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 + 2016 2017 2018 2042 + ); + } + + A15BE31TRU { + id 25311; + type latUniverse; + origin (0.0 0.0 0.0); + pitch (1.26 1.26 0.0); + shape (17 17 0); + padMat Water; + map ( + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 13 31 31 13 31 31 13 31 31 31 31 31 + 31 31 31 13 31 31 31 31 31 31 31 31 31 13 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 10 31 31 10 31 31 10 31 31 10 31 31 13 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 10 31 31 10 31 31 14 31 31 10 31 31 13 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 10 31 31 10 31 31 10 31 31 10 31 31 13 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 10 31 31 31 31 31 31 31 31 31 13 31 31 31 + 31 31 31 31 31 10 31 31 10 31 31 10 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 ); } + + A15BE31TRUSleeve { + id 253110; + type cellUniverse; + cells ( + 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 + 2016 2017 2018 2043 + ); + } + + A15BE31TLU { + id 253111; + type latUniverse; + origin (0.0 0.0 0.0); + pitch (1.26 1.26 0.0); + shape (17 17 0); + padMat Water; + map ( + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 13 31 31 13 31 31 13 31 31 31 31 31 + 31 31 31 13 31 31 31 31 31 31 31 31 31 13 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 13 31 31 10 31 31 10 31 31 10 31 31 10 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 13 31 31 10 31 31 14 31 31 10 31 31 10 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 13 31 31 10 31 31 10 31 31 10 31 31 10 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 13 31 31 31 31 31 31 31 31 31 10 31 31 31 + 31 31 31 31 31 10 31 31 10 31 31 10 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 ); } + + A15BE31TLUSleeve { + id 2531110; + type cellUniverse; + cells ( + 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 + 2016 2017 2018 2044 + ); + } + + A16BE31U { + id 2631; + type latUniverse; + origin (0.0 0.0 0.0); + pitch (1.26 1.26 0.0); + shape (17 17 0); + padMat Water; + map ( + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 10 31 31 10 31 31 10 31 31 31 31 31 + 31 31 31 10 31 31 31 31 31 31 31 31 31 10 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 10 31 31 13 31 31 13 31 31 13 31 31 10 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 10 31 31 13 31 31 14 31 31 13 31 31 10 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 10 31 31 13 31 31 13 31 31 13 31 31 10 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 10 31 31 31 31 31 31 31 31 31 10 31 31 31 + 31 31 31 31 31 10 31 31 10 31 31 10 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 ); } + + A16BE31USleeve { + id 26310; + type cellUniverse; + cells ( + 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 + 2016 2017 2018 2045 + ); + } + + A20BE31U { + id 3031; + type latUniverse; + origin (0.0 0.0 0.0); + pitch (1.26 1.26 0.0); + shape (17 17 0); + padMat Water; + map ( + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 10 31 31 10 31 31 10 31 31 31 31 31 + 31 31 31 10 31 31 31 31 31 31 31 31 31 10 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 10 31 31 10 31 31 13 31 31 10 31 31 10 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 10 31 31 13 31 31 14 31 31 13 31 31 10 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 10 31 31 10 31 31 13 31 31 10 31 31 10 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 10 31 31 31 31 31 31 31 31 31 10 31 31 31 + 31 31 31 31 31 10 31 31 10 31 31 10 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 + 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 31 ); } + + A20BE31USleeve { + id 30310; + type cellUniverse; + cells ( + 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 + 2016 2017 2018 2046 + ); + } + + A12BE24U { + id 2224; + type latUniverse; + origin (0.0 0.0 0.0); + pitch (1.26 1.26 0.0); + shape (17 17 0); + padMat Water; + map ( + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 24 24 24 10 24 24 13 24 24 10 24 24 24 24 24 + 24 24 24 10 24 24 24 24 24 24 24 24 24 10 24 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 10 24 24 13 24 24 13 24 24 13 24 24 10 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 13 24 24 13 24 24 14 24 24 13 24 24 13 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 10 24 24 13 24 24 13 24 24 13 24 24 10 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 24 10 24 24 24 24 24 24 24 24 24 10 24 24 24 + 24 24 24 24 24 10 24 24 13 24 24 10 24 24 24 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 ); } + + A12BE24USleeve { + id 22240; + type cellUniverse; + cells ( + 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 + 2016 2017 2018 2047 + ); + } + + A16BE24U { + id 2624; + type latUniverse; + origin (0.0 0.0 0.0); + pitch (1.26 1.26 0.0); + shape (17 17 0); + padMat Water; + map ( + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 24 24 24 10 24 24 10 24 24 10 24 24 24 24 24 + 24 24 24 10 24 24 24 24 24 24 24 24 24 10 24 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 10 24 24 13 24 24 13 24 24 13 24 24 10 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 10 24 24 13 24 24 14 24 24 13 24 24 10 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 10 24 24 13 24 24 13 24 24 13 24 24 10 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 24 10 24 24 24 24 24 24 24 24 24 10 24 24 24 + 24 24 24 24 24 10 24 24 10 24 24 10 24 24 24 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 ); } + + A16BE24USleeve { + id 26240; + type cellUniverse; + cells ( + 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 + 2016 2017 2018 2048 + ); + } + + A0E24DBank { + id 1425; + type latUniverse; + origin (0.0 0.0 0.0); + pitch (1.26 1.26 0.0); + shape (17 17 0); + padMat Water; + map ( + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 24 24 24 15 24 24 15 24 24 15 24 24 24 24 24 + 24 24 24 15 24 24 24 24 24 24 24 24 24 15 24 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 15 24 24 15 24 24 15 24 24 15 24 24 15 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 15 24 24 15 24 24 14 24 24 15 24 24 15 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 15 24 24 15 24 24 15 24 24 15 24 24 15 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 24 15 24 24 24 24 24 24 24 24 24 15 24 24 24 + 24 24 24 24 24 15 24 24 15 24 24 15 24 24 24 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 + 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 24 ); } + + // assembly with sleeves at different heights + A0E24SleeveDBank { + id 14250; + type cellUniverse; + cells ( + 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 + 2016 2017 2018 2049 + );} + + A0E16DBank { + id 1417; + type latUniverse; + origin (0.0 0.0 0.0); + pitch (1.26 1.26 0.0); + shape (17 17 0); + padMat Water; + map ( + 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 + 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 + 16 16 16 16 16 15 16 16 15 16 16 15 16 16 16 16 16 + 16 16 16 15 16 16 16 16 16 16 16 16 16 15 16 16 16 + 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 + 16 16 15 16 16 15 16 16 15 16 16 15 16 16 15 16 16 + 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 + 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 + 16 16 15 16 16 15 16 16 14 16 16 15 16 16 15 16 16 + 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 + 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 + 16 16 15 16 16 15 16 16 15 16 16 15 16 16 15 16 16 + 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 + 16 16 16 15 16 16 16 16 16 16 16 16 16 15 16 16 16 + 16 16 16 16 16 15 16 16 15 16 16 15 16 16 16 16 16 + 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 + 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16); } + + A0E16SleeveDBank { + id 14170; + type cellUniverse; + cells ( + 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 + 2016 2017 2018 2050 + ); + } + + + leftBaffleUni { + id 5222; + type cellUniverse; + cells (52 53 54);} + + + rightBaffleUni { + id 5223; + type cellUniverse; + cells (55 56 57);} + + topBaffleUni { + id 5224; + type cellUniverse; + cells (58 59 60);} + + bottomBaffleUni { + id 5225; + type cellUniverse; + cells (61 62 63);} + + + topLeft { + id 5226; + type cellUniverse; + cells (64 65 66 67 68);} + + topRight { + id 5227; + type cellUniverse; + cells ( 69 70 71 72 73);} + + BottomLeft { + id 5228; + type cellUniverse; + cells ( 74 75 76 77 78);} + + BottomRight { + id 5229; + type cellUniverse; + cells ( 79 80 81 82 83);} + + + SQTL { + id 1500; + type cellUniverse; + cells (84 85 86 87);} + + + SQTR { + id 1600; + type cellUniverse; + cells (88 89 90 91);} + + SQBL { + id 1700; + type cellUniverse; + cells (92 93 94 95);} + + SQBR { + id 1800; + type cellUniverse; + cells (96 97 98 99);} + + + latCore { + id 9999; + type latUniverse; + origin (0.0 0.0 0.0); + pitch (21.50364 21.50364 0.0); + shape (17 17 0); + padMat Water; + map ( + 1001 1001 1001 1001 1800 5224 5224 5224 5224 5224 5224 5224 1700 1001 1001 1001 1001 + 1001 1001 1800 5224 5229 24310 7031120 24310 7031120 24310 7031120 24310 5228 5224 1700 1001 1001 + 1001 1800 5229 24310 14310 26310 14160 30310 14160 30310 14160 26310 14310 24310 5228 1700 1001 + 1001 5222 24310 2531110 26240 14160 26240 14160 26240 14160 26240 14160 26240 25311 24310 5223 1001 + 1800 5229 14310 26240 14250 26240 24160 22240 14160 22240 24160 26240 14250 26240 14310 5228 1700 + 5222 24310 26310 14160 26240 24160 22240 24160 22240 24160 22240 24160 26240 14160 26310 24310 5223 + 5222 703190 14160 26240 24160 22240 14160 22240 14160 22240 14160 22240 24160 26240 14160 703130 5223 + 5222 24310 30310 14160 22240 24160 22240 24160 26240 24160 22240 24160 22240 14160 30310 24310 5223 + 5222 703190 14160 26240 14160 22240 14160 26240 14170 26240 14160 22240 14160 26240 14160 703130 5223 + 5222 24310 30310 14160 22240 24160 22240 24160 26240 24160 22240 24160 22240 14160 30310 24310 5223 + 5222 703190 14160 26240 24160 22240 14160 22240 14160 22240 14160 22240 24160 26240 14160 703130 5223 + 5222 24310 26310 14160 26240 24160 22240 24160 22240 24160 22240 24160 26240 14160 26310 24310 5223 + 1600 5227 14310 26240 14250 26240 24160 22240 14160 22240 24160 26240 14250 26240 14310 5226 1500 + 1001 5222 24310 253170 26240 14160 26240 14160 26240 14160 26240 14160 26240 253150 24310 5223 1001 + 1001 1600 5227 24310 14310 26310 14160 30310 14160 30310 14160 26310 14310 24310 5226 1500 1001 + 1001 1001 1600 5225 5227 24310 703160 24310 703160 24310 703160 24310 5226 5225 1500 1001 1001 + 1001 1001 1001 1001 1600 5225 5225 5225 5225 5225 5225 5225 1500 1001 1001 1001 1001 ); } + +! Note partial rodded assemblies end with a 1, i.e., 14170 and 14250 + + coreAndStructures { + id 8888; + type cellUniverse; + cells (7 8 9 10 11 12 13 14 15 16 17);} + + + + } +} + + +viz { + bmpZ { + type bmp; + output imgXY; + what material; + centre (-17.13 240.69 167.74); + width (50 50); + axis z; + res (1000 1000); + } + bmpYZ { + type bmp; + output imgYZ; + what material; + centre (0.0 0.0 232.0); + width (100.0 200.0); + axis x; + res (1000 2000); + } +} + + +nuclearData { + handles { + ce {type aceNeutronDatabase; aceLibrary $SCONE_ACE; ures 0; } + } + + materials { + + // Note that commented nuclide densities are included in the specification + // but are not available in the JEFF-3.11 data library + + Air { + temp 566; + composition { + 18036.06 7.8730E-09; + 18038.06 1.4844E-09; + 18040.06 2.3506E-06; + 6012.06 6.7539E-08; + //6013.06 7.5658E-10; + 7014.06 1.9680E-04; + 7015.06 7.2354E-07; + 8016.06 5.2866E-05; + 8017.06 2.0084E-08; + //8018.06 1.0601E-07; + } + } + + SS304 { + temp 566; + composition { + 24050.06 7.6778E-04; + 24052.06 1.4806E-02; + 24053.06 1.6789E-03; + 24054.06 4.1791E-04; + 26054.06 3.4620E-03; + 26056.06 5.4345E-02; + 26057.06 1.2551E-03; + 26058.06 1.6703E-04; + 25055.06 1.7604E-03; + 28058.06 5.6089E-03; + 28060.06 2.1605E-03; + 28061.06 9.3917E-05; + 28062.06 2.9945E-04; + 28064.06 7.6261E-05; + 14028.06 9.5281E-04; + 14029.06 4.8381E-05; + 14030.06 3.1893E-05; } + } + + Helium { + temp 566; + composition { + 2003.06 4.8089E-10; + 2004.06 2.4044E-04; } + } + + BorosilicateGlass { + temp 566; + composition { + 13027.06 1.7352E-03; + 5010.06 9.6506E-04; + 5011.06 3.9189E-03; + 8016.06 4.6514E-02; + 8017.06 1.7671E-05; + //8018.06 9.3268E-05; + 14028.06 1.6926E-02; + 14029.06 8.5944E-04; + 14030.06 5.6654E-04; } + } + + Water { + temp 566; + moder {1001.06 (lwj3.11 lwj3.09); } + composition { + 5010.06 7.9714E-06; + 5011.06 3.2247E-05; + 1001.06 4.9456E-02; + 1002.06 7.7035E-06; + 8016.06 2.4673E-02; + 8017.06 9.3734E-06; + //8018.06 4.9474E-05; + } + } + + Zircaloy { + temp 566; + composition { + 24050.06 3.2962E-06; + 24052.06 6.3564E-05; + 24053.06 7.2076E-06; + 24054.06 1.7941E-06; + 26054.06 8.6698E-06; + 26056.06 1.3610E-04; + 26057.06 3.1431E-06; + 26058.06 4.1829E-07; + 8016.06 3.0744E-04; + 8017.06 1.1680E-07; + //8018.03 6.1648E-07; + 50112.06 4.6735E-06; + 50114.06 3.1799E-06; + 50115.06 1.6381E-06; + 50116.06 7.0055E-05; + 50117.06 3.7003E-05; + 50118.06 1.1669E-04; + 50119.06 4.1387E-05; + 50120.06 1.5697E-04; + 50122.06 2.2308E-05; + 50124.06 2.7897E-05; + 40090.06 2.1828E-02; + 40091.06 4.7601E-03; + 40092.06 7.2759E-03; + 40094.06 7.3734E-03; + 40096.06 1.1879E-03; } + } + + Inconel{ + temp 566; + composition { + 24050.06 7.8239E-04; + 24052.06 1.5088E-02; + 24053.06 1.7108E-03; + 24054.06 4.2586E-04; + 26054.06 1.4797E-03; + 26056.06 2.3229E-02; + 26057.06 5.3645E-04; + 26058.06 7.1392E-05; + 25055.06 7.8201E-04; + 28058.06 2.9320E-02; + 28060.06 1.1294E-02; + 28061.06 4.9094E-04; + 28062.06 1.5653E-03; + 28064.06 3.9864E-04; + 14028.06 5.6757E-04; + 14029.06 2.8820E-05; + 14030.06 1.8998E-05; } + } + + B4C{ + temp 566; + composition { + 5010.06 1.5206E-02; + 5011.06 6.1514E-02; + 6012.06 1.8972E-02; + //6013.06 2.1252E-04; + } + } + + Ag-In-Cd{ + temp 566; + composition { + 47107.06 2.3523E-02; + 47109.06 2.1854E-02; + 48106.06 3.3882E-05; + 48108.06 2.4166E-05; + 48110.06 3.3936E-04; + 48111.06 3.4821E-04; + 48112.06 6.5611E-04; + 48113.06 3.3275E-04; + 48114.06 7.8252E-04; + 48116.06 2.0443E-04; + 49113.06 3.4219E-04; + 49115.06 7.6511E-03; } + } + + UO2-16 { + temp 566; + tms 1; + composition { + 8016.03 4.5897E-02; + 8017.03 1.7436E-05; + //8018.03 9.2032E-05; + 92234.03 3.0131E-06; + 92235.03 3.7503E-04; + 92238.03 2.2625E-02;} + } + + UO2-24 { + temp 566; + tms 1; + composition { + 8016.03 4.5830E-02; + 8017.03 1.7411E-05; + //8018.03 9.1898E-05; + 92234.03 4.4842E-06; + 92235.03 5.5814E-04; + 92238.03 2.2407E-02;} + } + + UO2-31 { + temp 566; + tms 1; + composition { + 8016.03 4.5853E-02; + 8017.03 1.7420E-05; + //8018.03 9.1942E-05; + 92234.03 5.7987E-06; + 92235.03 7.2175E-04; + 92238.03 2.2253E-02;} + } + + UO2-32 { + temp 566; + tms 1; + composition { + 8016.03 4.6029E-02; + 8017.03 1.7487E-05; + //8018.03 9.2296E-05; + 92234.03 5.9959E-06; + 92235.03 7.4630E-04; + 92238.03 2.2317E-02; + } + } + + UO2-34 { + temp 566; + tms 1; + composition { + 8016.03 4.6110E-02; + 8017.03 1.7517E-05; + //8018.03 9.2459E-05; + 92234.03 6.4018E-06; + 92235.03 7.9681E-04; + 92238.03 2.2307E-02;} + } + + // vanadium51 was stated twice in carbonsteel below + // in the beavrs pdf - typo? + CarbonSteel { + temp 566; + composition { + 13027.06 4.3523E-05; + 5010.06 2.5833E-06; + 5011.06 1.0450E-05; + 6012.06 1.0442E-03; + //6013.06 1.1697E-05 ; + 20040.06 1.7043E-05; + 20042.06 1.1375E-07; + 20043.06 2.3734E-08; + 20044.06 3.6673E-07; + 20046.06 7.0322E-10; + 20048.06 3.2875E-08; + 24050.06 1.3738E-05; + 24052.06 2.6493E-04; + 24053.06 3.0041E-05; + 24054.06 7.4778E-06; + 29063.06 1.0223E-04; + 29065.06 4.5608E-05; + 26054.06 4.7437E-03; + 26056.06 7.4465E-02; + 26057.06 1.7197E-03; + 26058.06 2.2886E-04; + 25055.06 6.4126E-04; + 42100.06 2.9814E-05; + 42092.06 4.4822E-05; + 42094.06 2.8110E-05; + 42095.06 4.8567E-05; + 42096.06 5.1015E-05; + 42097.06 2.9319E-05; + 42098.06 7.4327E-05; + 41093.06 5.0559E-06; + 28058.06 4.0862E-04; + 28060.06 1.5740E-04; + 28061.06 6.8420E-06; + 28062.06 2.1815E-05; + 28064.06 5.5557E-06; + 15031.06 3.7913E-05; + 16032.06 3.4808E-05; + 16033.06 2.7420E-07; + 16034.06 1.5368E-06; + 16036.06 5.3398E-09; + 14028.06 6.1702E-04; + 14029.06 3.1330E-05; + 14030.06 2.0653E-05; + 22046.06 1.2144E-06; + 22047.06 1.0952E-06; + 22048.06 1.0851E-05; + 22049.06 7.9634E-07; + 22050.06 7.6249E-07; + //23050.06 1.1526E-07; + 23051.06 4.5989E-05; + } + } + + SupportPlateSS { + temp 566; + composition { + 24050.06 3.5223E-04; + 24052.06 6.7924E-03; + 24053.06 7.7020E-04; + 24054.06 1.9172E-04; + 26054.06 1.5882E-03; + 26056.06 2.4931E-02; + 26057.06 5.7578E-04; + 26058.06 7.6625E-05; + 25055.06 8.0762E-04; + 28058.06 2.5731E-03; + 28060.06 9.9117E-04; + 28061.06 4.3085E-05; + 28062.06 1.3738E-04; + 28064.06 3.4985E-05; + 14028.06 4.3711E-04; + 14029.06 2.2195E-05; + 14030.06 1.4631E-05;} + } + + SupportPlateBW { + temp 566; + moder {1001.06 (lwj3.11 lwj3.09); } + composition { + 5010.06 1.0559E-05; + 5011.06 4.2716E-05; + 1001.06 6.5512E-02; + 1002.06 1.0204E-05; + 8016.06 3.2683E-02; + 8017.06 1.2416E-05; + //8018.06 6.5535E-05; + } + } + + +} +} diff --git a/SharedModules/universalVariables.f90 b/SharedModules/universalVariables.f90 index fdb8a2d57..90d54e732 100644 --- a/SharedModules/universalVariables.f90 +++ b/SharedModules/universalVariables.f90 @@ -7,7 +7,7 @@ module universalVariables ! *** DON't CHANGE THIS. HARDCODED IS FINE ! CHANGE THIS: NUMBER MUST BE CALCULATED DURING INITIAL GEOMETRY PROCESSING ! Problematic for separating modules! - integer(shortInt), parameter, public :: HARDCODED_MAX_NEST = 8 + integer(shortInt), parameter, public :: HARDCODED_MAX_NEST = 12 integer(shortInt), parameter, public :: MAX_OUTGOING_PARTICLES = 5 ! CHANGE THIS: NUMBER WILL DEPEND ON SYSTEM ARCHITECTURE diff --git a/UserInterface/Graphics/imgBmp_func.f90 b/UserInterface/Graphics/imgBmp_func.f90 index 0bbcecf37..d76131416 100644 --- a/UserInterface/Graphics/imgBmp_func.f90 +++ b/UserInterface/Graphics/imgBmp_func.f90 @@ -3,11 +3,11 @@ !! !! Only 24-bit Color Bitmaps (without any compression or color palette data) can be created. !! -!! It is not ideal due to large size of the BMP files. However, the purpose of theese functions -!! is to provide some image creataion capability without external dependencies. +!! It is not ideal due to large size of the BMP files. However, the purpose of these functions +!! is to provide some image creation capability without external dependencies. !! !! Implementation should work even when the Fortran compiler prints data in big-endian format. -!! (I guess it depends on the compiler+system combination. I do wander what gfortran would +!! (I guess it depends on the compiler+system combination. I do wonder what gfortran would !! do for the bi-endian CPU) !! !! Has been based on: diff --git a/Visualisation/visualiser_class.f90 b/Visualisation/visualiser_class.f90 index a0cdcd768..a5303dda3 100644 --- a/Visualisation/visualiser_class.f90 +++ b/Visualisation/visualiser_class.f90 @@ -2,13 +2,13 @@ module visualiser_class use numPrecision use universalVariables - use genericProcedures, only : fatalError, numToChar + use genericProcedures, only : fatalError, numToChar, crossProduct use hashFunctions_func, only : knuthHash, FNV_1 use imgBmp_func, only : imgBmp_toFile use commandLineUI, only : getInputFile use dictionary_class, only : dictionary use geometry_inter, only : geometry - use materialMenu_mod, only : mm_colourMap => colourMap + use materialMenu_mod, only : mm_colourMap => colourMap, mm_nameMap => nameMap use outputVTK_class implicit none @@ -52,6 +52,7 @@ module visualiser_class procedure :: kill procedure, private :: makeVTK procedure, private :: makeBmpImg + procedure, private :: makeRayPlot end type contains @@ -121,9 +122,12 @@ subroutine makeViz(self) case('bmp') call self % makeBmpImg(tempDict) + + case('ray') + call self % makeRayPlot(tempDict) case default - call fatalError(here, 'Unrecognised visualisation - presently only accept vtk') + call fatalError(here, 'Unrecognised visualisation - presently only accept vtk, bmp, and ray') end select @@ -326,6 +330,191 @@ subroutine makeBmpImg(self, dict) call imgBmp_toFile(img, outputFile) end subroutine makeBmpImg + + !! + !! Generate a ray traced image of the geometry using Phong's approximation + !! + !! Follows the papers from Gavin Ridley and Brian Nease at SNA+MC 2024. + !! + !! The user can input the point *towards which* the camera is pointed, the point from which + !! the target is viewed, the horizontal field of view (defaults to recommended 70 degrees), + !! the point from which the light originates, and the number of pixels in 2D. + !! Additionally, the user may specify which materials are transparent. Otherwise the plot may + !! be quite boring. + !! + !! Optionally, the user can specify which direction is 'up' to determine the orientation of + !! the camera, the relative amount of diffuse light, the field of view, and the colour offset. + !! + !! Args: + !! dict [in] -> Dictionary with settings + !! + !! Sample dictionary input: + !! ray_plot { + !! type ray; + !! output img; // Name of output file without extension + !! centre (0.0 0.0 0.0); // Coordinates of the centre/target of the plot + !! camera (0.0 0.0 0.0); // Coordinates of the camera + !! res (300 300); // Resolution of the image + !! #light (0.0 0.0 0.0);# // Coordinates of the light source, defaults to camera position + !! #up (0.0 0.0 1.0);# // Which way is 'up'? Determines the rotation of the camera + !! #diffuse 0.3;# // Fraction of light which is diffuse rather than direct + !! #fov 70;# // Field-of-view in the horizontal axis in degrees + !! #offset 978; # // Parameter to 'randomize' the colour map + !! #transparent (mat1, mat2, ... );# // Names of transparent materials + !! } + !! + !! + subroutine makeRayPlot(self, dict) + class(visualiser), intent(inout) :: self + class(dictionary), intent(in) :: dict + real(defReal) :: fov, diffuse + real(defReal), dimension(3) :: centre, up, camera, light, cv, ch, d + real(defReal), dimension(:), allocatable :: temp + integer(shortInt), dimension(:), allocatable :: tempInt + character(nameLen), dimension(:), allocatable :: matNames + integer(shortInt), dimension(:), allocatable :: mats + integer(shortInt), dimension(:,:), allocatable :: img, matIDs + real(defReal), dimension(:,:), allocatable :: lum + real(defReal), dimension(3,3) :: M + integer(shortInt) :: i, offset + character(10) :: time + character(8) :: date + character(nameLen) :: outputFile + character(100), parameter :: Here = 'makeRayPlot (visualiser_class.f90)' + + ! Get name of the output file + call dict % get(outputFile, 'output') + outputFile = trim(outputFile) // '.bmp' + + ! Central/target point + call dict % get(temp, 'centre') + + if (size(temp) /= 3) then + call fatalError(Here, "'centre' must have size 3. Has: "//numToChar(size(temp))) + end if + + centre = temp + deallocate(temp) + + ! Camera origin + call dict % get(temp, 'camera') + + if (size(temp) /= 3) then + call fatalError(Here, "'camera' must have size 3. Has: "//numToChar(size(temp))) + end if + camera = temp + deallocate(temp) + + ! Get light location or default to camera location + if (dict % isPresent('light')) then + call dict % get(temp, 'light') + + if (size(temp) /= 3) then + call fatalError(Here, "'light' must have size 3. Has: "//numToChar(size(temp))) + end if + light = temp + deallocate(temp) + + else + light = camera + end if + + ! The up direction, which sets the camera rotation + if (dict % isPresent('up')) then + call dict % get(temp, 'up') + + if (size(temp) /= 3) then + call fatalError(Here, "'up' must have size 3. Has: "//numToChar(size(temp))) + end if + up = temp + up = up / norm2(up) + deallocate(temp) + + else + up = [0, 0, 1] + end if + + ! Optional field of view in degrees + ! Convert to radians + call dict % getOrDefault(fov, 'fov', 70.0_defReal) + fov = fov * PI / 180 + + ! Optional fraction of diffuse light + call dict % getOrDefault(diffuse, 'diffuse', 0.3_defReal) + if ((diffuse > 1) .or. (diffuse < 0)) call fatalError(Here,'The fraction of diffuse light must be between 0 and 1') + + ! Create governing vectors for the plot + d = (centre - camera) + d = d/norm2(d) + + ! Ensure that up is not colinear with the view direction + if (all(abs(crossProduct(up, d)) < 1E-6)) call fatalError(Here,"View direction is co-linear with 'up'.") + + cv = crossProduct(d, up) + cv = cv / norm2(cv) + + ch = crossProduct(cv, d) + ch = ch / norm2(ch) + + ! Create coordinate matrix + M(:,1) = d + M(:,2) = cv + M(:,3) = ch + + ! Resolution + call dict % get(tempInt, 'res') + + if (size(tempInt) /= 2) then + call fatalError(Here, "'res' must have size 2. Has: "//numToChar(size(tempInt))) + else if (any(tempInt <= 0)) then + call fatalError(Here, "Resolution must be +ve. There is 0 or -ve entry!") + end if + + allocate(matIDs(tempInt(1), tempInt(2))) + allocate(lum(tempInt(1), tempInt(2))) + + ! Transparent materials + ! Defaults to an array containing only VOID and OUTSIDE + if (dict % isPresent('transparent')) then + call dict % get(matNames, 'transparent') + allocate(mats(size(matNames) + 2)) + do i = 1, size(matNames) + mats(i) = mm_nameMap % get(matNames(i)) + end do + mats(i) = OUTSIDE_MAT + mats(i+1) = VOID_MAT + else + allocate(mats(2)) + mats(1) = OUTSIDE_MAT + mats(2) = VOID_MAT + end if + + ! Colourmap offset + ! If not given select randomly + if (dict % isPresent('offset')) then + call dict % get(offset, 'offset') + + else + call date_and_time(date, time) + call FNV_1(date // time, offset) + + end if + + ! Img contains luminosity values, matIDs identifies which materials were hit + call self % geom % rayPlot(lum, matIDs, camera, light, M, mats, fov, diffuse) + + ! Translate to an image. + ! Obtain material colours and scale by luminosity + matIDs = materialColour(matIDs, offset) + img = matIDs + + ! Scale image by brightness + img = brightnessScale(img, lum) + + ! Print image + call imgBmp_toFile(img, outputFile) + + end subroutine makeRayPlot !! !! Terminates visualiser @@ -394,5 +583,28 @@ elemental function uniqueIDColour(uniqueID) result(colour) end function uniqueIDColour + !! + !! Scales an integer image by the provided value of brightness + !! + elemental function brightnessScale(img0, brightness) result(img) + integer(shortInt), intent(in) :: img0 + real(defReal), intent(in) :: brightness + integer(shortInt) :: img + integer(shortInt) :: r, g, b + + ! Extract RGB components (assuming 0x00RRGGBB) + r = ishft(iand(img0, Z'00FF0000'), -16) + g = ishft(iand(img0, Z'0000FF00'), -8) + b = iand(img0, Z'000000FF') + + ! Scale and clamp to 0–255 + r = max(0, min(255, int(r * brightness))) + g = max(0, min(255, int(g * brightness))) + b = max(0, min(255, int(b * brightness))) + + ! Reassemble RGB integer + img = ior(ior(ishft(r, 16), ishft(g, 8)), b) + + end function brightnessScale end module visualiser_class diff --git a/docs/User Manual.rst b/docs/User Manual.rst index 56cc5eb9f..90d487a44 100644 --- a/docs/User Manual.rst +++ b/docs/User Manual.rst @@ -704,6 +704,8 @@ The possible types of files that the geometry is plotted in are: vtk ### +This produces a VTK output of the 3D geometry. + * corner: (x y z) array with the corner of the geometry [cm] * width: (x y z) array with the width of the mesh in each direction [cm] * vox: (x y z) array with the number of voxels requested in each direction @@ -718,6 +720,8 @@ Example: :: bmp ### +This produces a 2D slice plot of the geometry. + * centre: (x y z) array with the coordinates of the center of the plot [cm] * axis: ``x``, ``y`` or ``z``, it's the axis normal to the 2D plot * width (*optional*, default = whole geometry): (y z), (x z) or (x y) array @@ -736,6 +740,39 @@ Example: :: plotBMP { type bmp; axis z; centre (0.0 0.0 0.0); width (50 10); res (1000 200); output geomZ; what material; } +ray +### + +This produce a ray-traced plot of the geometry, coloured by material, with lighting. + +* centre: (x y z) array with the coordinates of the center of the plot [cm] +* camera: (x y z) array with the coordinates of the camera, which points towards centre [cm] +* res: (x y) array with the resolution of the plot in each direction +* output: name of the output file, with extension ``.bmp`` +* light (*optional*, default = camera): (x y z) array with the coordinates of the light source [cm]. + Defaults to the position of the camera. +* up (*optional*, default = (0 0 1)): (x y z) array which defines which direction is up [-]. Used + to change the orientation of the camera. +* diffuse (*optional*, default = 0.3): real value between 0 and 1 which defines the fraction of light + contributed by ambient sources, rather than directly from the light source [-]. +* fov (*optional*, default = 70): real value between 0 and 180 which defines the horizontal field of view + of the camera. [degrees] +* offset (*optional*, default = random): An integer (positive or negative) that + shifts the sequence of colours assigned to materials. Allows to change colours + from the default sequence in a parametric way. +* transparent (*optional*, default = void): A list of material names which are transparent to the rays. + It is advisable to include a few materials here to have an interesting plot. Void is always included as + a transparent material. + +Example: :: + + + ray_plot { + type ray; output myRayPlot; centre (0.0 0.0 0.0); camera (-10.0 0.0 20.0); + res (300 300); light (10.0 0.0 20.0); up (0.0 0.0 1.0); diffuse 0.2; + fov 70; offset 978; transparent (coolant, air, outerSteel ); + } + .. note:: SCONE can be run to visualise geometry without actually doing transport, by including ``--plot`` when running the application. In this case the visualiser