diff --git a/openmc/summary.py b/openmc/summary.py index 2ae7464846..c2c02da18a 100644 --- a/openmc/summary.py +++ b/openmc/summary.py @@ -213,7 +213,7 @@ class Summary(object): fill = self._f['geometry/cells'][key]['lattice'].value if 'surfaces' in self._f['geometry/cells'][key].keys(): - surfaces = self._f['geometry/cells'][key]['surfaces'][...] + surfaces = self._f['geometry/cells'][key]['surfaces'].value.decode() else: surfaces = [] @@ -241,12 +241,15 @@ class Summary(object): self._cell_fills[index] = (fill_type, fill) # Iterate over all Surfaces and add them to the Cell - for surface_halfspace in surfaces: - - halfspace = np.sign(surface_halfspace) - surface_id = abs(surface_halfspace) - surface = self.get_surface_by_id(surface_id) - cell.add_surface(surface, halfspace) + for token in surfaces.split(): + try: + surface_halfspace = int(token) + halfspace = np.sign(surface_halfspace) + surface_id = abs(surface_halfspace) + surface = self.get_surface_by_id(surface_id) + cell.add_surface(surface, halfspace) + except ValueError: + pass # Add the Cell to the global dictionary of all Cells self.cells[index] = cell diff --git a/src/geometry.F90 b/src/geometry.F90 index 2cde01e023..9e40155e52 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -17,56 +17,77 @@ module geometry contains !=============================================================================== -! SIMPLE_CELL_CONTAINS determines whether a given the current coordinates of the -! particle are inside a cell defined as the intersection of a series of surfaces +! CELL_CONTAINS determines if a cell contains the particle at a given +! location. The bounds of the cell are detemined by a logical expression +! involving surface half-spaces. At initialization, the expression was converted +! to RPN notation. In cell_contains, we evaluate the RPN expression using a +! stack, similar to how a RPN calculator would work. !=============================================================================== - function simple_cell_contains(c, p) result(in_cell) + pure function cell_contains(c, p) result(in_cell) type(Cell), intent(in) :: c - type(Particle), intent(inout) :: p + type(Particle), intent(in) :: p logical :: in_cell - integer :: i ! index of surfaces in cell - integer :: i_surface ! index in surfaces array (with sign) - logical :: specified_sense ! specified sense of surface in list + integer :: i + integer :: token + logical :: b1, b2 + integer :: i_stack logical :: actual_sense ! sense of particle wrt surface - class(Surface), pointer :: s + logical :: stack(size(c%rpn)) - SURFACE_LOOP: do i = 1, c % n_surfaces - ! Lookup surface - i_surface = c % surfaces(i) - - ! Check if the particle is currently on the specified surface - if (i_surface == p % surface) then - ! Particle is heading into the cell - cycle - elseif (i_surface == -p % surface) then - ! Particle is heading out of the cell - in_cell = .false. - return + i_stack = 0 + do i = 1, size(c%rpn) + token = c%rpn(i) + if (token < OP_UNION) then + ! If the token is not an operator, evaluate the sense of particle with + ! respect to the surface and see if the token matches the sense. If the + ! particle's surface attribute is set and matches the token, that + ! overrides the determination based on sense(). + i_stack = i_stack + 1 + if (token == p%surface) then + stack(i_stack) = .true. + elseif (-token == p%surface) then + stack(i_stack) = .false. + else + actual_sense = surfaces(abs(token))%obj%sense(& + p%coord(p%n_coord)%xyz, p%coord(p%n_coord)%uvw) + stack(i_stack) = (actual_sense .eqv. (token > 0)) + end if + else + ! If the token is a binary operator (intersection/union), apply it to + ! the last two items on the stack. If the token is a unary operator + ! (complement), apply it to the last item on the stack. + b1 = stack(i_stack) + select case (token) + case (OP_UNION) + b2 = stack(i_stack - 1) + stack(i_stack - 1) = b1 .or. b2 + i_stack = i_stack - 1 + case (OP_INTERSECTION) + b2 = stack(i_stack - 1) + stack(i_stack - 1) = b1 .and. b2 + i_stack = i_stack - 1 + case (OP_COMPLEMENT) + stack(i_stack) = .not. b1 + end select end if + end do - ! Determine the specified sense of the surface in the cell and the actual - ! sense of the particle with respect to the surface - s => surfaces(abs(i_surface))%obj - actual_sense = s%sense(p%coord(p%n_coord)%xyz, p%coord(p%n_coord)%uvw) - specified_sense = (c % surfaces(i) > 0) - - ! Compare sense of point to specified sense - if (actual_sense .neqv. specified_sense) then - in_cell = .false. - return - end if - end do SURFACE_LOOP - - ! If we've reached here, then the sense matched on every surface or there - ! are no surfaces. - in_cell = .true. - end function simple_cell_contains + if (i_stack == 1) then + ! The one remaining logical on the stack indicates whether the particle is + ! in the cell. + in_cell = stack(i_stack) + else + ! This case occurs if there is no surface specification since i_stack will + ! still be zero. + in_cell = .true. + end if + end function cell_contains !=============================================================================== ! CHECK_CELL_OVERLAP checks for overlapping cells at the current particle's -! position using simple_cell_contains and the LocalCoord's built up by find_cell +! position using cell_contains and the LocalCoord's built up by find_cell !=============================================================================== subroutine check_cell_overlap(p) @@ -93,7 +114,7 @@ contains index_cell = univ % cells(i) c => cells(index_cell) - if (simple_cell_contains(c, p)) then + if (cell_contains(c, p)) then ! the particle should only be contained in one cell per level if (index_cell /= p % coord(j) % cell) then call fatal_error("Overlapping cells detected: " & @@ -162,7 +183,7 @@ contains c => cells(index_cell) ! Move on to the next cell if the particle is not inside this cell - if (.not. simple_cell_contains(c, p)) cycle + if (.not. cell_contains(c, p)) cycle ! Set cell on this level p % coord(j) % cell = index_cell @@ -556,7 +577,7 @@ contains ! ======================================================================= ! FIND MINIMUM DISTANCE TO SURFACE IN THIS CELL - SURFACE_LOOP: do i = 1, cl % n_surfaces + SURFACE_LOOP: do i = 1, size(cl%surfaces) ! check for operators index_surf = cl%surfaces(i) coincident = (index_surf == p % surface) @@ -823,10 +844,11 @@ contains c => cells(i) ! loop over each surface specification - do j = 1, c % n_surfaces + do j = 1, size(c%surfaces) i_surface = c % surfaces(j) positive = (i_surface > 0) i_surface = abs(i_surface) + if (i_surface >= OP_UNION) cycle if (positive) then count_positive(i_surface) = count_positive(i_surface) + 1 else @@ -853,10 +875,11 @@ contains c => cells(i) ! loop over each surface specification - do j = 1, c % n_surfaces + do j = 1, size(c%surfaces) i_surface = c % surfaces(j) positive = (i_surface > 0) i_surface = abs(i_surface) + if (i_surface >= OP_UNION) cycle if (positive) then count_positive(i_surface) = count_positive(i_surface) + 1 diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index 13d806e8cf..c8c99783f3 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -128,14 +128,13 @@ module geometry_header ! the geom integer :: material ! Material within cell (0 for ! universe) - integer :: n_surfaces ! Number of surfaces within integer, allocatable :: offset (:) ! Distribcell offset for tally ! counter - integer, allocatable :: & - & surfaces(:) ! List of surfaces bounding cell - ! -- note that parentheses, union, - ! etc operators will be listed here - ! too + integer, allocatable :: surfaces(:) ! List of surfaces bounding cell -- + ! note that parentheses, union, etc + ! operators will be listed here too + integer, allocatable :: rpn(:) ! Reverse Polish notation for surface + ! expression ! Rotation matrix and translation vector real(8), allocatable :: translation(:) diff --git a/src/initialize.F90 b/src/initialize.F90 index 3671a40af7..d46fd8d48a 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -571,7 +571,7 @@ contains ! ADJUST SURFACE LIST FOR EACH CELL c => cells(i) - do j = 1, c%n_surfaces + do j = 1, size(c%surfaces) id = c%surfaces(j) if (id < OP_UNION) then if (surface_dict%has_key(abs(id))) then @@ -584,6 +584,14 @@ contains end if end do + do j = 1, size(c%rpn) + id = c%rpn(j) + if (id < OP_UNION) then + i_array = surface_dict%get_key(abs(id)) + c%rpn(j) = sign(i_array, id) + end if + end do + ! ======================================================================= ! ADJUST UNIVERSE INDEX FOR EACH CELL diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 308bdc66b0..156d33ce46 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -13,8 +13,9 @@ module input_xml use plot_header use random_lcg, only: prn use surface_header + use stl_vector, only: VectorInt use string, only: to_lower, to_str, str_to_int, str_to_real, & - starts_with, ends_with + starts_with, ends_with, tokenize use tally_header, only: TallyObject, TallyFilter use tally_initialize, only: add_tallies use xml_interface @@ -993,6 +994,7 @@ contains logical :: boundary_exists character(MAX_LINE_LEN) :: filename character(MAX_WORD_LEN) :: word + character(MAX_LINE_LEN) :: surface_spec type(Cell), pointer :: c class(Surface), pointer :: s class(Lattice), pointer :: lat @@ -1004,6 +1006,8 @@ contains type(NodeList), pointer :: node_surf_list => null() type(NodeList), pointer :: node_rlat_list => null() type(NodeList), pointer :: node_hlat_list => null() + type(VectorInt) :: tokens + type(VectorInt) :: rpn ! Display output message call write_message("Reading geometry XML file...", 5) @@ -1116,16 +1120,26 @@ contains ! Allocate array for surfaces and copy if (check_for_node(node_cell, "surfaces")) then - n = get_arraysize_integer(node_cell, "surfaces") - else - n = 0 - end if - c % n_surfaces = n + call get_node_value(node_cell, "surfaces", surface_spec) + if (len_trim(surface_spec) > 0) then + ! Create surfaces array from string + call tokenize(surface_spec, tokens) - if (n > 0) then - allocate(c % surfaces(n)) - call get_node_array(node_cell, "surfaces", c % surfaces) + ! Use shunting-yard algorithm to determine RPN for surface algorithm + call generate_rpn(tokens, rpn) + + ! Copy surface spec and RPN form to cell arrays + allocate(c % surfaces(tokens%size())) + allocate(c % rpn(rpn%size())) + c % surfaces(:) = tokens%data(1:tokens%size()) + c % rpn(:) = rpn%data(1:rpn%size()) + + call tokens%clear() + call rpn%clear() + end if end if + if (.not. allocated(c%surfaces)) allocate(c%surfaces(0)) + if (.not. allocated(c%rpn)) allocate(c%rpn(0)) ! Rotation matrix if (check_for_node(node_cell, "rotation")) then @@ -4766,4 +4780,88 @@ contains end subroutine expand_natural_element +!=============================================================================== +! GENERATE_RPN implements the shunting-yard algorithm to generate a Reverse +! polish notation (RPN) expression for the surface specification of a cell given +! the infix notation. +!=============================================================================== + + subroutine generate_rpn(tokens, output) + type(VectorInt), intent(in) :: tokens ! infix notation + type(VectorInt), intent(inout) :: output ! RPN notation + + integer :: i + integer :: token + integer :: op + type(VectorInt) :: stack + + do i = 1, tokens%size() + token = tokens%data(i) + + if (token < OP_UNION) then + ! If token is not an operator, add it to output + call output%push_back(token) + + elseif (token < OP_RIGHT_PAREN) then + ! Regular operators union, intersection, complement + do while (stack%size() > 0) + op = stack%data(stack%size()) + + if (op < OP_RIGHT_PAREN .and. & + ((token == OP_COMPLEMENT .and. token < op) .or. & + (token /= OP_COMPLEMENT .and. token <= op))) then + ! While there is an operator, op, on top of the stack, if the token + ! is left-associative and its precedence is less than or equal to + ! that of op or if the token is right-associative and its precedence + ! is less than that of op, move op to the output queue and push the + ! token on to the stack. Note that only complement is + ! right-associative. + call output%push_back(op) + call stack%pop_back() + else + exit + end if + end do + + call stack%push_back(token) + + elseif (token == OP_LEFT_PAREN) then + ! If the token is a left parenthesis, push it onto the stack + call stack%push_back(token) + + else + ! If the token is a right parenthesis, move operators from the stack to + ! the output queue until reaching the left parenthesis. + do + ! 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 surface specification') + end if + + op = stack%data(stack%size()) + if (op == OP_LEFT_PAREN) exit + call output%push_back(op) + call stack%pop_back() + end do + + ! Pop the left parenthesis. + call stack%pop_back() + end if + end do + + ! While there are operators on the stack, move them to the output queue + do while (stack%size() > 0) + op = stack%data(stack%size()) + + ! If the operator is a parenthesis, it is mismatched + if (op >= OP_RIGHT_PAREN) then + call fatal_error('Mimatched parentheses in surface specification') + end if + + call output%push_back(op) + call stack%pop_back() + end do + end subroutine generate_rpn + end module input_xml diff --git a/src/relaxng/geometry.rnc b/src/relaxng/geometry.rnc index 82f1f15b32..07a02af9ea 100644 --- a/src/relaxng/geometry.rnc +++ b/src/relaxng/geometry.rnc @@ -9,7 +9,7 @@ element geometry { (element material { ( xsd:int | "void" ) } | attribute material { ( xsd:int | "void" ) }) ) & - (element surfaces { list { xsd:int* } } | attribute surfaces { list { xsd:int* } })? & + (element surfaces { xsd:string } | attribute surfaces { xsd:string })? & (element rotation { list { xsd:double+ } } | attribute rotation { list { xsd:double+ } })? & (element translation { list { xsd:double+ } } | attribute translation { list { xsd:double+ } })? }* diff --git a/src/relaxng/geometry.rng b/src/relaxng/geometry.rng index 9bd573b346..6a68aa721e 100644 --- a/src/relaxng/geometry.rng +++ b/src/relaxng/geometry.rng @@ -63,18 +63,10 @@ - - - - - + - - - - - + diff --git a/src/string.F90 b/src/string.F90 index 5f57bb3855..2d690dded5 100644 --- a/src/string.F90 +++ b/src/string.F90 @@ -1,8 +1,10 @@ module string - use constants, only: MAX_WORDS, MAX_LINE_LEN, ERROR_INT, ERROR_REAL + use constants, only: MAX_WORDS, MAX_LINE_LEN, ERROR_INT, ERROR_REAL, & + OP_LEFT_PAREN, OP_RIGHT_PAREN, OP_COMPLEMENT, OP_INTERSECTION, OP_UNION use error, only: fatal_error, warning use global, only: master + use stl_vector, only: VectorInt implicit none @@ -63,63 +65,81 @@ contains end subroutine split_string !=============================================================================== -! SPLIT_STRING_WL takes a string that includes logical expressions for a list of -! bounding surfaces in a cell and splits it into separate words. The characters -! (, ), :, and # count as separate words since they represent operators. -! -! Arguments: -! string = input line -! words = array of words -! n = number of words +! TOKENIZE takes a string that includes logical expressions for a list of +! bounding surfaces in a cell and splits it into separate tokens. The characters +! (, ), ^, and ~ count as separate tokens since they represent operators. !=============================================================================== - subroutine split_string_wl(string, words, n) + subroutine tokenize(string, tokens) + character(*), intent(in) :: string + type(VectorInt), intent(inout) :: tokens - character(*), intent(in) :: string - character(*), intent(out) :: words(MAX_WORDS) - integer, intent(out) :: n + integer :: i ! current index + integer :: i_start ! starting index of word + character(len=len_trim(string)) :: string_ - character(1) :: chr ! current character - integer :: i ! current index - integer :: i_start ! starting index of word - integer :: i_end ! ending index of word + ! Remove leading blanks + string_ = adjustl(string) i_start = 0 - i_end = 0 - n = 0 - do i = 1, len_trim(string) - chr = string(i:i) - + i = 1 + do while (i <= len_trim(string_)) ! Check for special characters - if (index('():#', chr) > 0) then + if (index('()^~ ', string_(i:i)) > 0) then + ! If the special character appears immediately after a non-operator, + ! create a token with the surface half-space if (i_start > 0) then - i_end = i - 1 - n = n + 1 - words(n) = string(i_start:i_end) + call tokens%push_back(int(str_to_int(& + string_(i_start:i - 1)), 4)) end if - n = n + 1 - words(n) = chr + + select case (string_(i:i)) + case ('(') + call tokens%push_back(OP_LEFT_PAREN) + case (')') + call tokens%push_back(OP_RIGHT_PAREN) + case ('^') + call tokens%push_back(OP_UNION) + case ('~') + call tokens%push_back(OP_COMPLEMENT) + case (' ') + ! Find next non-space character + do while (string_(i+1:i+1) == ' ') + i = i + 1 + end do + + ! If previous token is not an operator and next token is not a left + ! parenthese or union operator, that implies that the whitespace is to + ! be interpreted as an intersection operator + if (i_start > 0) then + if (index(')^', string_(i+1:i+1)) == 0) then + call tokens%push_back(OP_INTERSECTION) + end if + end if + end select + i_start = 0 - i_end = 0 - cycle + else + ! Check for invalid characters + if (index('-0123456789', string_(i:i)) == 0) then + call fatal_error("Invalid character '" // string_(i:i) // "' in & + &surface specification.") + end if + + ! If we haven't yet reached the start of a word, start a new word + if (i_start == 0) i_start = i end if - if ((i_start == 0) .and. (chr /= ' ')) then - i_start = i - end if - if (i_start > 0) then - if (chr == ' ') i_end = i - 1 - if (i == len_trim(string)) i_end = i - if (i_end > 0) then - n = n + 1 - words(n) = string(i_start:i_end) - ! reset indices - i_start = 0 - i_end = 0 - end if - end if + i = i + 1 end do - end subroutine split_string_wl + + ! If we've reached the end and we're still in a word, create a token from it + ! and add it to the list + if (i_start > 0) then + call tokens%push_back(int(str_to_int(& + string_(i_start:len_trim(string_))), 4)) + end if + end subroutine tokenize !=============================================================================== ! CONCATENATE takes an array of words and concatenates them together in one diff --git a/src/summary.F90 b/src/summary.F90 index 587f00b6b7..49fe802305 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -107,13 +107,13 @@ contains integer :: i, j, k, m integer, allocatable :: lattice_universes(:,:,:) - integer, allocatable :: surface_ids(:) integer(HID_T) :: geom_group integer(HID_T) :: cells_group, cell_group integer(HID_T) :: surfaces_group, surface_group integer(HID_T) :: universes_group, univ_group integer(HID_T) :: lattices_group, lattice_group real(8), allocatable :: coeffs(:) + character(MAX_LINE_LEN) :: surface_spec type(Cell), pointer :: c class(Surface), pointer :: s type(Universe), pointer :: u @@ -176,15 +176,26 @@ contains end select ! Write list of bounding surfaces - if (c%n_surfaces > 0) then - allocate(surface_ids(c%n_surfaces)) - do j = 1, c%n_surfaces - k = c%surfaces(j) - surface_ids(j) = sign(surfaces(abs(k))%obj%id, k) - end do - call write_dataset(cell_group, "surfaces", surface_ids) - deallocate(surface_ids) - end if + surface_spec = "" + do j = 1, size(c%surfaces) + k = c%surfaces(j) + if (k < OP_UNION) then + surface_spec = trim(surface_spec) // " " // to_str(& + sign(surfaces(abs(k))%obj%id, k)) + else + select case(k) + case (OP_LEFT_PAREN) + surface_spec = trim(surface_spec) // " (" + case (OP_RIGHT_PAREN) + surface_spec = trim(surface_spec) // " (" + case (OP_COMPLEMENT) + surface_spec = trim(surface_spec) // " ~" + case (OP_UNION) + surface_spec = trim(surface_spec) // " ^" + end select + end if + end do + call write_dataset(cell_group, "surfaces", adjustl(surface_spec)) call close_group(cell_group) end do CELL_LOOP