Respond to more @smharper comments on #463

Move abstract interface block in surface_header. Make error messages more
informative. Add description in openmc-update-inputs epilog.

scripts/openmc-update-inputs

@@ -36,9 +36,10 @@ they will be moved to a new file with '.original' appended to their name.
 
 Formatting changes that will be made:
 
-geometry.xml:  Lattices containing 'outside' attributes/tags will be replaced
+geometry.xml: Lattices containing 'outside' attributes/tags will be replaced
   with lattices containing 'outer' attributes, and the appropriate
-  cells/universes will be added.
+  cells/universes will be added. Any 'surfaces' attributes/elements on a cell
+  will be renamed 'region'.
 """
 
 

src/input_xml.F90

@@ -1135,7 +1135,7 @@ contains
         call tokenize(region_spec, tokens)
 
         ! Use shunting-yard algorithm to determine RPN for surface algorithm
-        call generate_rpn(tokens, rpn)
+        call generate_rpn(c%id, tokens, rpn)
 
         ! Copy region spec and RPN form to cell arrays
         allocate(c % region(tokens%size()))
@@ -4801,7 +4801,8 @@ contains
 ! the infix notation.
 !===============================================================================
 
-  subroutine generate_rpn(tokens, output)
+  subroutine generate_rpn(cell_id, tokens, output)
+    integer, intent(in) :: cell_id
     type(VectorInt), intent(in) :: tokens    ! infix notation
     type(VectorInt), intent(inout) :: output ! RPN notation
 
@@ -4851,7 +4852,8 @@ contains
           ! If we run out of operators without finding a left parenthesis, it
           ! means there are mismatched parentheses.
           if (stack%size() == 0) then
-            call fatal_error('Mimatched parentheses in region specification')
+            call fatal_error('Mimatched parentheses in region specification &
+                 &for cell ' // trim(to_str(cell_id)) // '.')
           end if
 
           op = stack%data(stack%size())
@@ -4871,7 +4873,8 @@ contains
 
       ! If the operator is a parenthesis, it is mismatched
       if (op >= OP_RIGHT_PAREN) then
-        call fatal_error('Mimatched parentheses in region specification')
+        call fatal_error('Mimatched parentheses in region specification &
+             &for cell ' // trim(to_str(cell_id)) // '.')
       end if
 
       call output%push_back(op)

src/string.F90

@@ -76,6 +76,7 @@ contains
 
     integer :: i       ! current index
     integer :: i_start ! starting index of word
+    integer :: token
     character(len=len_trim(string)) :: string_
 
     ! Remove leading blanks
@@ -97,8 +98,22 @@ contains
         case ('(')
           call tokens%push_back(OP_LEFT_PAREN)
         case (')')
+          if (tokens%size() > 0) then
+            token = tokens%data(tokens%size())
+            if (token >= OP_UNION .and. token < OP_RIGHT_PAREN) then
+              call fatal_error("Right parentheses cannot follow an operator in &
+                   &region specification: " // trim(string))
+            end if
+          end if
           call tokens%push_back(OP_RIGHT_PAREN)
         case ('|')
+          if (tokens%size() > 0) then
+            token = tokens%data(tokens%size())
+            if (.not. (token < OP_UNION .or. token == OP_RIGHT_PAREN)) then
+              call fatal_error("Union cannot follow an operator in region &
+                   &specification: " // trim(string))
+            end if
+          end if
           call tokens%push_back(OP_UNION)
         case ('~')
           call tokens%push_back(OP_COMPLEMENT)

src/surface_header.F90

@@ -24,10 +24,45 @@ module surface_header
     procedure(iNormal),   deferred :: normal
   end type Surface
 
+  abstract interface
+    pure function iEvaluate(this, xyz) result(f)
+      import Surface
+      class(Surface), intent(in) :: this
+      real(8), intent(in) :: xyz(3)
+      real(8) :: f
+    end function iEvaluate
+
+    pure function iDistance(this, xyz, uvw, coincident) result(d)
+      import Surface
+      class(Surface), intent(in) :: this
+      real(8), intent(in) :: xyz(3)
+      real(8), intent(in) :: uvw(3)
+      logical, intent(in) :: coincident
+      real(8) :: d
+    end function iDistance
+
+    pure function iNormal(this, xyz) result(uvw)
+      import Surface
+      class(Surface), intent(in) :: this
+      real(8), intent(in) :: xyz(3)
+      real(8) :: uvw(3)
+    end function iNormal
+  end interface
+
+!===============================================================================
+! SURFACECONTAINER allows us to store an array of different types of surfaces
+!===============================================================================
+
   type :: SurfaceContainer
     class(Surface), allocatable :: obj
   end type SurfaceContainer
 
+!===============================================================================
+! All the derived types below are extensions of the abstract Surface type. They
+! inherent the reflect() and sense() type-bound procedures and must implement
+! evaluate(), distance(), and normal()
+!===============================================================================
+
   type, extends(Surface) :: SurfaceXPlane
     ! x = x0
     real(8) :: x0
@@ -148,31 +183,6 @@ module surface_header
     procedure :: normal => z_cone_normal
   end type SurfaceZCone
 
-  abstract interface
-    pure function iEvaluate(this, xyz) result(f)
-      import Surface
-      class(Surface), intent(in) :: this
-      real(8), intent(in) :: xyz(3)
-      real(8) :: f
-    end function iEvaluate
-
-    pure function iDistance(this, xyz, uvw, coincident) result(d)
-      import Surface
-      class(Surface), intent(in) :: this
-      real(8), intent(in) :: xyz(3)
-      real(8), intent(in) :: uvw(3)
-      logical, intent(in) :: coincident
-      real(8) :: d
-    end function iDistance
-
-    pure function iNormal(this, xyz) result(uvw)
-      import Surface
-      class(Surface), intent(in) :: this
-      real(8), intent(in) :: xyz(3)
-      real(8) :: uvw(3)
-    end function iNormal
-  end interface
-
 contains
 
 !===============================================================================
@@ -369,14 +379,14 @@ contains
     real(8) :: d
 
     real(8) :: f
-    real(8) :: tmp
+    real(8) :: projection
 
     f = this%A*xyz(1) + this%B*xyz(2) + this%C*xyz(3) - this%D
-    tmp = this%A*uvw(1) + this%B*uvw(2) + this%C*uvw(3)
-    if (coincident .or. abs(f) < FP_COINCIDENT .or. tmp == ZERO) then
+    projection = this%A*uvw(1) + this%B*uvw(2) + this%C*uvw(3)
+    if (coincident .or. abs(f) < FP_COINCIDENT .or. projection == ZERO) then
       d = INFINITY
     else
-      d = -f/tmp
+      d = -f/projection
       if (d < ZERO) d = INFINITY
     end if
   end function plane_distance