Merge pull request #1012 from smharper/c_cell

Move most geometry data to C++

CMakeLists.txt

@@ -432,9 +432,12 @@ set(LIBOPENMC_FORTRAN_SRC
   src/tallies/trigger_header.F90
 )
 set(LIBOPENMC_CXX_SRC
+  src/cell.cpp
   src/initialize.cpp
   src/finalize.cpp
+  src/geometry_aux.cpp
   src/hdf5_interface.cpp
+  src/lattice.cpp
   src/math_functions.cpp
   src/message_passing.cpp
   src/plot.cpp

docs/source/devguide/styleguide.rst

@@ -219,8 +219,8 @@ Curly braces
 
 For a function definition, the opening and closing braces should each be on
 their own lines.  This helps distinguish function code from the argument list.
-If the entire function fits on one line, then the braces can be on the same
-line. e.g.:
+If the entire function fits on one or two lines, then the braces can be on the
+same line. e.g.:
 
 .. code-block:: C++
 
@@ -238,6 +238,9 @@ line. e.g.:
 
     int return_one() {return 1;}
 
+    int return_one()
+    {return 1;}
+
 For a conditional, the opening brace should be on the same line as the end of
 the conditional statement. If there is a following ``else if`` or ``else``
 statement, the closing brace should be on the same line as that following

src/api.F90

@@ -207,7 +207,7 @@ contains
 
     if (found) then
       if (rtype == 1) then
-        id = cells(p % coord(p % n_coord) % cell) % id
+        id = cells(p % coord(p % n_coord) % cell) % id()
       elseif (rtype == 2) then
         if (p % material == MATERIAL_VOID) then
           id = 0

src/cell.cpp

@@ -0,0 +1,523 @@
+#include "cell.h"
+
+#include <cmath>
+#include <limits>
+#include <sstream>
+#include <string>
+
+#include "constants.h"
+#include "error.h"
+#include "hdf5_interface.h"
+#include "lattice.h"
+#include "surface.h"
+#include "xml_interface.h"
+
+
+namespace openmc {
+
+//==============================================================================
+// Constants
+//==============================================================================
+
+constexpr int32_t OP_LEFT_PAREN   {std::numeric_limits<int32_t>::max()};
+constexpr int32_t OP_RIGHT_PAREN  {std::numeric_limits<int32_t>::max() - 1};
+constexpr int32_t OP_COMPLEMENT   {std::numeric_limits<int32_t>::max() - 2};
+constexpr int32_t OP_INTERSECTION {std::numeric_limits<int32_t>::max() - 3};
+constexpr int32_t OP_UNION        {std::numeric_limits<int32_t>::max() - 4};
+
+extern "C" double FP_PRECISION;
+
+//==============================================================================
+// Global variables
+//==============================================================================
+
+int32_t n_cells {0};
+
+std::vector<Cell*> global_cells;
+std::unordered_map<int32_t, int32_t> cell_map;
+
+std::vector<Universe*> global_universes;
+std::unordered_map<int32_t, int32_t> universe_map;
+
+//==============================================================================
+//! Convert region specification string to integer tokens.
+//!
+//! The characters (, ), |, and ~ count as separate tokens since they represent
+//! operators.
+//==============================================================================
+
+std::vector<int32_t>
+tokenize(const std::string region_spec) {
+  // Check for an empty region_spec first.
+  std::vector<int32_t> tokens;
+  if (region_spec.empty()) {
+    return tokens;
+  }
+
+  // Parse all halfspaces and operators except for intersection (whitespace).
+  for (int i = 0; i < region_spec.size(); ) {
+    if (region_spec[i] == '(') {
+      tokens.push_back(OP_LEFT_PAREN);
+      i++;
+
+    } else if (region_spec[i] == ')') {
+      tokens.push_back(OP_RIGHT_PAREN);
+      i++;
+
+    } else if (region_spec[i] == '|') {
+      tokens.push_back(OP_UNION);
+      i++;
+
+    } else if (region_spec[i] == '~') {
+      tokens.push_back(OP_COMPLEMENT);
+      i++;
+
+    } else if (region_spec[i] == '-' || region_spec[i] == '+'
+               || std::isdigit(region_spec[i])) {
+      // This is the start of a halfspace specification.  Iterate j until we
+      // find the end, then push-back everything between i and j.
+      int j = i + 1;
+      while (j < region_spec.size() && std::isdigit(region_spec[j])) {j++;}
+      tokens.push_back(std::stoi(region_spec.substr(i, j-i)));
+      i = j;
+
+    } else if (std::isspace(region_spec[i])) {
+      i++;
+
+    } else {
+      std::stringstream err_msg;
+      err_msg << "Region specification contains invalid character, \""
+              << region_spec[i] << "\"";
+      fatal_error(err_msg);
+    }
+  }
+
+  // Add in intersection operators where a missing operator is needed.
+  int i = 0;
+  while (i < tokens.size()-1) {
+    bool left_compat {(tokens[i] < OP_UNION) || (tokens[i] == OP_RIGHT_PAREN)};
+    bool right_compat {(tokens[i+1] < OP_UNION)
+                       || (tokens[i+1] == OP_LEFT_PAREN)
+                       || (tokens[i+1] == OP_COMPLEMENT)};
+    if (left_compat && right_compat) {
+      tokens.insert(tokens.begin()+i+1, OP_INTERSECTION);
+    }
+    i++;
+  }
+
+  return tokens;
+}
+
+//==============================================================================
+//! Convert infix region specification to Reverse Polish Notation (RPN)
+//!
+//! This function uses the shunting-yard algorithm.
+//==============================================================================
+
+std::vector<int32_t>
+generate_rpn(int32_t cell_id, std::vector<int32_t> infix)
+{
+  std::vector<int32_t> rpn;
+  std::vector<int32_t> stack;
+
+  for (int32_t token : infix) {
+    if (token < OP_UNION) {
+      // If token is not an operator, add it to output
+      rpn.push_back(token);
+
+    } else if (token < OP_RIGHT_PAREN) {
+      // Regular operators union, intersection, complement
+      while (stack.size() > 0) {
+        int32_t op = stack.back();
+
+        if (op < OP_RIGHT_PAREN &&
+             ((token == OP_COMPLEMENT && token < op) ||
+             (token != OP_COMPLEMENT && token <= op))) {
+          // 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.
+          rpn.push_back(op);
+          stack.pop_back();
+        } else {
+          break;
+        }
+      }
+
+      stack.push_back(token);
+
+    } else if (token == OP_LEFT_PAREN) {
+      // If the token is a left parenthesis, push it onto the stack
+      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.
+      for (auto it = stack.rbegin(); *it != OP_LEFT_PAREN; it++) {
+        // If we run out of operators without finding a left parenthesis, it
+        // means there are mismatched parentheses.
+        if (it == stack.rend()) {
+          std::stringstream err_msg;
+          err_msg << "Mismatched parentheses in region specification for cell "
+                  << cell_id;
+          fatal_error(err_msg);
+        }
+
+        rpn.push_back(stack.back());
+        stack.pop_back();
+      }
+
+      // Pop the left parenthesis.
+      stack.pop_back();
+    }
+  }
+
+  while (stack.size() > 0) {
+    int32_t op = stack.back();
+
+    // If the operator is a parenthesis it is mismatched.
+    if (op >= OP_RIGHT_PAREN) {
+      std::stringstream err_msg;
+      err_msg << "Mismatched parentheses in region specification for cell "
+              << cell_id;
+      fatal_error(err_msg);
+    }
+
+    rpn.push_back(stack.back());
+    stack.pop_back();
+  }
+
+  return rpn;
+}
+
+//==============================================================================
+// Cell implementation
+//==============================================================================
+
+Cell::Cell(pugi::xml_node cell_node)
+{
+  if (check_for_node(cell_node, "id")) {
+    id = stoi(get_node_value(cell_node, "id"));
+  } else {
+    fatal_error("Must specify id of cell in geometry XML file.");
+  }
+
+  //TODO: don't automatically lowercase cell and surface names
+  if (check_for_node(cell_node, "name")) {
+    name = get_node_value(cell_node, "name");
+  }
+
+  if (check_for_node(cell_node, "universe")) {
+    universe = stoi(get_node_value(cell_node, "universe"));
+  } else {
+    universe = 0;
+  }
+
+  if (check_for_node(cell_node, "fill")) {
+    fill = stoi(get_node_value(cell_node, "fill"));
+  } else {
+    fill = C_NONE;
+  }
+
+  if (check_for_node(cell_node, "material")) {
+    //TODO: read material ids.
+    material.push_back(C_NONE+1);
+    material.shrink_to_fit();
+  } else {
+    material.push_back(C_NONE);
+    material.shrink_to_fit();
+  }
+
+  // Make sure that either material or fill was specified.
+  if ((material[0] == C_NONE) && (fill == C_NONE)) {
+    std::stringstream err_msg;
+    err_msg << "Neither material nor fill was specified for cell " << id;
+    fatal_error(err_msg);
+  }
+
+  // Make sure that material and fill haven't been specified simultaneously.
+  if ((material[0] != C_NONE) && (fill != C_NONE)) {
+    std::stringstream err_msg;
+    err_msg << "Cell " << id << " has both a material and a fill specified; "
+            << "only one can be specified per cell";
+    fatal_error(err_msg);
+  }
+
+  // Read the region specification.
+  std::string region_spec;
+  if (check_for_node(cell_node, "region")) {
+    region_spec = get_node_value(cell_node, "region");
+  }
+
+  // Get a tokenized representation of the region specification.
+  region = tokenize(region_spec);
+  region.shrink_to_fit();
+
+  // Convert user IDs to surface indices.
+  for (auto &r : region) {
+    if (r < OP_UNION) {
+      r = copysign(surface_map[abs(r)] + 1, r);
+    }
+  }
+
+  // Convert the infix region spec to RPN.
+  rpn = generate_rpn(id, region);
+  rpn.shrink_to_fit();
+
+  // Check if this is a simple cell.
+  simple = true;
+  for (int32_t token : rpn) {
+    if ((token == OP_COMPLEMENT) || (token == OP_UNION)) {
+      simple = false;
+      break;
+    }
+  }
+}
+
+//==============================================================================
+
+bool
+Cell::contains(const double xyz[3], const double uvw[3],
+               int32_t on_surface) const
+{
+  if (simple) {
+    return contains_simple(xyz, uvw, on_surface);
+  } else {
+    return contains_complex(xyz, uvw, on_surface);
+  }
+}
+
+//==============================================================================
+
+std::pair<double, int32_t>
+Cell::distance(const double xyz[3], const double uvw[3],
+               int32_t on_surface) const
+{
+  double min_dist {INFTY};
+  int32_t i_surf {std::numeric_limits<int32_t>::max()};
+
+  for (int32_t token : rpn) {
+    // Ignore this token if it corresponds to an operator rather than a region.
+    if (token >= OP_UNION) {continue;}
+
+    // Calculate the distance to this surface.
+    // Note the off-by-one indexing
+    bool coincident {token == on_surface};
+    double d {surfaces_c[abs(token)-1]->distance(xyz, uvw, coincident)};
+
+    // Check if this distance is the new minimum.
+    if (d < min_dist) {
+      if (std::abs(d - min_dist) / min_dist >= FP_PRECISION) {
+        min_dist = d;
+        i_surf = -token;
+      }
+    }
+  }
+
+  return {min_dist, i_surf};
+}
+
+//==============================================================================
+
+void
+Cell::to_hdf5(hid_t cell_group) const
+{
+  if (!name.empty()) {
+    write_string(cell_group, "name", name, false);
+  }
+
+  //TODO: Fix the off-by-one indexing.
+  write_int(cell_group, 0, nullptr, "universe",
+            &global_universes[universe-1]->id, false);
+
+  // Write the region specification.
+  if (!region.empty()) {
+    std::stringstream region_spec {};
+    for (int32_t token : region) {
+      if (token == OP_LEFT_PAREN) {
+        region_spec << " (";
+      } else if (token == OP_RIGHT_PAREN) {
+        region_spec << " )";
+      } else if (token == OP_COMPLEMENT) {
+        region_spec << " ~";
+      } else if (token == OP_INTERSECTION) {
+      } else if (token == OP_UNION) {
+        region_spec << " |";
+      } else {
+        // Note the off-by-one indexing
+        region_spec << " " << copysign(surfaces_c[abs(token)-1]->id, token);
+      }
+    }
+    write_string(cell_group, "region", region_spec.str(), false);
+  }
+}
+
+//==============================================================================
+
+bool
+Cell::contains_simple(const double xyz[3], const double uvw[3],
+                      int32_t on_surface) const
+{
+  for (int32_t token : rpn) {
+    if (token < OP_UNION) {
+      // 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().
+      if (token == on_surface) {
+      } else if (-token == on_surface) {
+        return false;
+      } else {
+        // Note the off-by-one indexing
+        bool sense = surfaces_c[abs(token)-1]->sense(xyz, uvw);
+        if (sense != (token > 0)) {return false;}
+      }
+    }
+  }
+  return true;
+}
+
+//==============================================================================
+
+bool
+Cell::contains_complex(const double xyz[3], const double uvw[3],
+                       int32_t on_surface) const
+{
+  // Make a stack of booleans.  We don't know how big it needs to be, but we do
+  // know that rpn.size() is an upper-bound.
+  bool stack[rpn.size()];
+  int i_stack = -1;
+
+  for (int32_t token : rpn) {
+    // 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.
+    if (token == OP_UNION) {
+      stack[i_stack-1] = stack[i_stack-1] || stack[i_stack];
+      i_stack --;
+    } else if (token == OP_INTERSECTION) {
+      stack[i_stack-1] = stack[i_stack-1] && stack[i_stack];
+      i_stack --;
+    } else if (token == OP_COMPLEMENT) {
+      stack[i_stack] = !stack[i_stack];
+    } else {
+      // 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 ++;
+      if (token == on_surface) {
+        stack[i_stack] = true;
+      } else if (-token == on_surface) {
+        stack[i_stack] = false;
+      } else {
+        // Note the off-by-one indexing
+        bool sense = surfaces_c[abs(token)-1]->sense(xyz, uvw);;
+        stack[i_stack] = (sense == (token > 0));
+      }
+    }
+  }
+
+  if (i_stack == 0) {
+    // The one remaining bool on the stack indicates whether the particle is
+    // in the cell.
+    return stack[i_stack];
+  } else {
+    // This case occurs if there is no region specification since i_stack will
+    // still be -1.
+    return true;
+  }
+}
+
+//==============================================================================
+// Non-method functions
+//==============================================================================
+
+extern "C" void
+read_cells(pugi::xml_node *node)
+{
+  // Count the number of cells.
+  for (pugi::xml_node cell_node: node->children("cell")) {n_cells++;}
+  if (n_cells == 0) {
+    fatal_error("No cells found in geometry.xml!");
+  }
+
+  // Allocate the vector of Cells.
+  global_cells.reserve(n_cells);
+
+  // Loop over XML cell elements and populate the array.
+  for (pugi::xml_node cell_node: node->children("cell")) {
+    global_cells.push_back(new Cell(cell_node));
+  }
+
+  // Populate the Universe vector and map.
+  for (int i = 0; i < global_cells.size(); i++) {
+    int32_t uid = global_cells[i]->universe;
+    auto it = universe_map.find(uid);
+    if (it == universe_map.end()) {
+      global_universes.push_back(new Universe());
+      global_universes.back()->id = uid;
+      global_universes.back()->cells.push_back(i);
+      universe_map[uid] = global_universes.size() - 1;
+    } else {
+      global_universes[it->second]->cells.push_back(i);
+    }
+  }
+}
+
+//==============================================================================
+// Fortran compatibility functions
+//==============================================================================
+
+extern "C" {
+  Cell* cell_pointer(int32_t cell_ind) {return global_cells[cell_ind];}
+
+  int32_t cell_id(Cell *c) {return c->id;}
+
+  void cell_set_id(Cell *c, int32_t id) {c->id = id;}
+
+  int cell_type(Cell *c) {return c->type;}
+
+  void cell_set_type(Cell *c, int type) {c->type = type;}
+
+  int32_t cell_universe(Cell *c) {return c->universe;}
+
+  void cell_set_universe(Cell *c, int32_t universe) {c->universe = universe;}
+
+  int32_t cell_fill(Cell *c) {return c->fill;}
+
+  int32_t* cell_fill_ptr(Cell *c) {return &c->fill;}
+
+  int32_t cell_n_instances(Cell *c) {return c->n_instances;}
+
+  bool cell_simple(Cell *c) {return c->simple;}
+
+  bool cell_contains(Cell *c, double xyz[3], double uvw[3], int32_t on_surface)
+  {return c->contains(xyz, uvw, on_surface);}
+
+  void cell_distance(Cell *c, double xyz[3], double uvw[3], int32_t on_surface,
+                     double *min_dist, int32_t *i_surf)
+  {
+    std::pair<double, int32_t> out = c->distance(xyz, uvw, on_surface);
+    *min_dist = out.first;
+    *i_surf = out.second;
+  }
+
+  int32_t cell_offset(Cell *c, int map) {return c->offset[map];}
+
+  void cell_to_hdf5(Cell *c, hid_t group) {c->to_hdf5(group);}
+
+  void extend_cells_c(int32_t n)
+  {
+    global_cells.reserve(global_cells.size() + n);
+    for (int32_t i = 0; i < n; i++) {
+      global_cells.push_back(new Cell());
+    }
+    n_cells = global_cells.size();
+  }
+}
+
+
+} // namespace openmc

src/cell.h

@@ -0,0 +1,118 @@
+#ifndef CELL_H
+#define CELL_H
+
+#include <cstdint>
+#include <string>
+#include <unordered_map>
+#include <vector>
+
+#include "hdf5.h"
+#include "pugixml/pugixml.hpp"
+
+
+namespace openmc {
+
+//==============================================================================
+// Constants
+//==============================================================================
+
+extern "C" int FILL_MATERIAL;
+extern "C" int FILL_UNIVERSE;
+extern "C" int FILL_LATTICE;
+
+//==============================================================================
+// Global variables
+//==============================================================================
+
+extern "C" int32_t n_cells;
+
+class Cell;
+extern std::vector<Cell*> global_cells;
+extern std::unordered_map<int32_t, int32_t> cell_map;
+
+class Universe;
+extern std::vector<Universe*> global_universes;
+extern std::unordered_map<int32_t, int32_t> universe_map;
+
+//==============================================================================
+//! A geometry primitive that fills all space and contains cells.
+//==============================================================================
+
+class Universe
+{
+public:
+  int32_t id;                  //!< Unique ID
+  std::vector<int32_t> cells;  //!< Cells within this universe
+  //double x0, y0, z0;           //!< Translation coordinates.
+};
+
+//==============================================================================
+//! A geometry primitive that links surfaces, universes, and materials
+//==============================================================================
+
+class Cell 
+{
+public:
+  int32_t id;                //!< Unique ID
+  std::string name;          //!< User-defined name
+  int type;                  //!< Material, universe, or lattice
+  int32_t universe;          //!< Universe # this cell is in
+  int32_t fill;              //!< Universe # filling this cell
+  int32_t n_instances{0};    //!< Number of instances of this cell
+
+  //! \brief Material(s) within this cell.
+  //!
+  //! May be multiple materials for distribcell.  C_NONE signifies a universe.
+  std::vector<int32_t> material;
+
+  //! Definition of spatial region as Boolean expression of half-spaces
+  std::vector<std::int32_t> region;
+  //! Reverse Polish notation for region expression
+  std::vector<std::int32_t> rpn;
+  bool simple;  //!< Does the region contain only intersections?
+
+  std::vector<int32_t> offset;  //!< Distribcell offset table
+
+  Cell() {};
+
+  explicit Cell(pugi::xml_node cell_node);
+
+  //! \brief Determine 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.
+  //!
+  //! The function is split into two cases, one for simple cells (those
+  //! involving only the intersection of half-spaces) and one for complex cells.
+  //! Simple cells can be evaluated with short circuit evaluation, i.e., as soon
+  //! as we know that one half-space is not satisfied, we can exit. This
+  //! provides a performance benefit for the common case. In
+  //! contains_complex, we evaluate the RPN expression using a stack, similar to
+  //! how a RPN calculator would work.
+  //! @param xyz[3] The 3D Cartesian coordinate to check.
+  //! @param uvw[3] A direction used to "break ties" the coordinates are very
+  //!   close to a surface.
+  //! @param on_surface The signed index of a surface that the coordinate is
+  //!   known to be on.  This index takes precedence over surface sense
+  //!   calculations.
+  bool
+  contains(const double xyz[3], const double uvw[3], int32_t on_surface) const;
+
+  //! Find the oncoming boundary of this cell.
+  std::pair<double, int32_t>
+  distance(const double xyz[3], const double uvw[3], int32_t on_surface) const;
+
+  //! \brief Write cell information to an HDF5 group.
+  //! @param group_id An HDF5 group id.
+  void to_hdf5(hid_t group_id) const;
+
+protected:
+  bool contains_simple(const double xyz[3], const double uvw[3],
+                       int32_t on_surface) const;
+  bool contains_complex(const double xyz[3], const double uvw[3],
+                        int32_t on_surface) const;
+};
+
+} // namespace openmc
+#endif // CELL_H

src/constants.F90

@@ -113,6 +113,9 @@ module constants
        FILL_MATERIAL = 1, & ! Cell with a specified material
        FILL_UNIVERSE = 2, & ! Cell filled by a separate universe
        FILL_LATTICE  = 3    ! Cell filled with a lattice
+  integer(C_INT), bind(C, name='FILL_MATERIAL') :: FILL_MATERIAL_C = FILL_MATERIAL
+  integer(C_INT), bind(C, name='FILL_UNIVERSE') :: FILL_UNIVERSE_C = FILL_UNIVERSE
+  integer(C_INT), bind(C, name='FILL_LATTICE')  :: FILL_LATTICE_C  = FILL_LATTICE
 
   ! Void material
   integer, parameter :: MATERIAL_VOID = -1
@@ -418,6 +421,7 @@ module constants
 
   ! indicates that an array index hasn't been set
   integer, parameter :: NONE = 0
+  integer, parameter :: C_NONE = -1
 
   ! Codes for read errors -- better hope these numbers are never used in an
   ! input file!

src/constants.h

@@ -0,0 +1,16 @@
+#ifndef CONSTANTS_H
+#define CONSTANTS_H
+
+#include <limits>
+
+
+namespace openmc{
+
+extern "C" double FP_COINCIDENT;
+extern "C" double FP_PRECISION;
+constexpr double INFTY {std::numeric_limits<double>::max()};
+constexpr int C_NONE {-1};
+
+} // namespace openmc
+
+#endif // CONSTANTS_H

src/geometry.F90

@@ -14,125 +14,36 @@ module geometry
 
   implicit none
 
+  interface
+    function cell_contains_c(cell_ptr, xyz, uvw, on_surface) &
+         bind(C, name="cell_contains") result(in_cell)
+      import C_PTR, C_DOUBLE, C_INT32_T, C_BOOL
+      type(C_PTR),        intent(in), value :: cell_ptr
+      real(C_DOUBLE),     intent(in)        :: xyz(3)
+      real(C_DOUBLE),     intent(in)        :: uvw(3)
+      integer(C_INT32_T), intent(in), value :: on_surface
+      logical(C_BOOL)                       :: in_cell
+    end function cell_contains_c
+
+    function count_universe_instances(search_univ, target_univ_id) bind(C) &
+         result(count)
+      import C_INT32_T, C_INT
+      integer(C_INT32_T), intent(in), value :: search_univ
+      integer(C_INT32_T), intent(in), value :: target_univ_id
+      integer(C_INT)                        :: count
+    end function
+  end interface
+
 contains
 
-!===============================================================================
-! 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.
-!
-! The function is split into two cases, one for simple cells (those involving
-! only the intersection of half-spaces) and one for complex cells. Simple cells
-! can be evaluated with short circuit evaluation, i.e., as soon as we know that
-! one half-space is not satisfied, we can exit. This provides a performance
-! benefit for the common case. In complex_cell_contains, we evaluate the RPN
-! expression using a stack, similar to how a RPN calculator would work.
-!===============================================================================
-
-  pure function cell_contains(c, p) result(in_cell)
+  function cell_contains(c, p) result(in_cell)
     type(Cell), intent(in) :: c
     type(Particle), intent(in) :: p
     logical :: in_cell
-
-    if (c%simple) then
-      in_cell = simple_cell_contains(c, p)
-    else
-      in_cell = complex_cell_contains(c, p)
-    end if
+    in_cell = cell_contains_c(c%ptr, p%coord(p%n_coord)%xyz, &
+                              p%coord(p%n_coord)%uvw, p%surface)
   end function cell_contains
 
-  pure function simple_cell_contains(c, p) result(in_cell)
-    type(Cell), intent(in) :: c
-    type(Particle), intent(in) :: p
-    logical :: in_cell
-
-    integer :: i
-    integer :: token
-    logical :: actual_sense    ! sense of particle wrt surface
-
-    in_cell = .true.
-    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().
-        if (token == p%surface) then
-          cycle
-        elseif (-token == p%surface) then
-          in_cell = .false.
-          exit
-        else
-          actual_sense = surfaces(abs(token)) % sense( &
-               p%coord(p%n_coord)%xyz, p%coord(p%n_coord)%uvw)
-          if (actual_sense .neqv. (token > 0)) then
-            in_cell = .false.
-            exit
-          end if
-        end if
-      end if
-    end do
-  end function simple_cell_contains
-
-  pure function complex_cell_contains(c, p) result(in_cell)
-    type(Cell), intent(in) :: c
-    type(Particle), intent(in) :: p
-    logical :: in_cell
-
-    integer :: i
-    integer :: token
-    integer :: i_stack
-    logical :: actual_sense    ! sense of particle wrt surface
-    logical :: stack(size(c%rpn))
-
-    i_stack = 0
-    do i = 1, size(c%rpn)
-      token = c%rpn(i)
-
-      ! 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.
-      select case (token)
-      case (OP_UNION)
-        stack(i_stack - 1) = stack(i_stack - 1) .or. stack(i_stack)
-        i_stack = i_stack - 1
-      case (OP_INTERSECTION)
-        stack(i_stack - 1) = stack(i_stack - 1) .and. stack(i_stack)
-        i_stack = i_stack - 1
-      case (OP_COMPLEMENT)
-        stack(i_stack) = .not. stack(i_stack)
-      case default
-        ! 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)) % sense( &
-               p%coord(p%n_coord)%xyz, p%coord(p%n_coord)%uvw)
-          stack(i_stack) = (actual_sense .eqv. (token > 0))
-        end if
-      end select
-
-    end do
-
-    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 region specification since i_stack will
-      ! still be zero.
-      in_cell = .true.
-    end if
-  end function complex_cell_contains
-
 !===============================================================================
 ! CHECK_CELL_OVERLAP checks for overlapping cells at the current particle's
 ! position using cell_contains and the LocalCoord's built up by find_cell
@@ -166,8 +77,8 @@ contains
           ! 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: " &
-                 &// trim(to_str(cells(index_cell) % id)) // ", " &
-                 &// trim(to_str(cells(p % coord(j) % cell) % id)) &
+                 &// trim(to_str(cells(index_cell) % id())) // ", " &
+                 &// trim(to_str(cells(p % coord(j) % cell) % id())) &
                  &// " on universe " // trim(to_str(univ % id)))
           end if
 
@@ -223,7 +134,7 @@ contains
       if (use_search_cells) then
         i_cell = search_cells(i)
         ! check to make sure search cell is in same universe
-        if (cells(i_cell) % universe /= i_universe) cycle
+        if (cells(i_cell) % universe() /= i_universe) cycle
       else
         i_cell = universes(i_universe) % cells(i)
       end if
@@ -236,7 +147,7 @@ contains
         ! Show cell information on trace
         if (verbosity >= 10 .or. trace) then
           call write_message("    Entering cell " // trim(to_str(&
-               cells(i_cell) % id)))
+               cells(i_cell) % id())))
         end if
 
         found = .true.
@@ -246,7 +157,7 @@ contains
 
     if (found) then
       associate(c => cells(i_cell))
-        CELL_TYPE: if (c % type == FILL_MATERIAL) then
+        CELL_TYPE: if (c % type() == FILL_MATERIAL) then
           ! ======================================================================
           ! AT LOWEST UNIVERSE, TERMINATE SEARCH
 
@@ -262,20 +173,20 @@ contains
             distribcell_index = c % distribcell_index
             offset = 0
             do k = 1, p % n_coord
-              if (cells(p % coord(k) % cell) % type == FILL_UNIVERSE) then
-                offset = offset + cells(p % coord(k) % cell) % &
-                     offset(distribcell_index)
-              elseif (cells(p % coord(k) % cell) % type == FILL_LATTICE) then
+              if (cells(p % coord(k) % cell) % type() == FILL_UNIVERSE) then
+                offset = offset + cells(p % coord(k) % cell) &
+                     % offset(distribcell_index-1)
+              elseif (cells(p % coord(k) % cell) % type() == FILL_LATTICE) then
                 if (lattices(p % coord(k + 1) % lattice) % obj &
                      % are_valid_indices([&
                      p % coord(k + 1) % lattice_x, &
                      p % coord(k + 1) % lattice_y, &
                      p % coord(k + 1) % lattice_z])) then
-                  offset = offset + lattices(p % coord(k + 1) % lattice) % obj % &
-                       offset(distribcell_index, &
-                       p % coord(k + 1) % lattice_x, &
+                  offset = offset + lattices(p % coord(k + 1) % lattice) % obj &
+                       % offset(distribcell_index - 1, &
+                       [p % coord(k + 1) % lattice_x, &
                        p % coord(k + 1) % lattice_y, &
-                       p % coord(k + 1) % lattice_z)
+                       p % coord(k + 1) % lattice_z])
                 end if
               end if
             end do
@@ -300,7 +211,7 @@ contains
             p % sqrtkT = c % sqrtkT(1)
           end if
 
-        elseif (c % type == FILL_UNIVERSE) then CELL_TYPE
+        elseif (c % type() == FILL_UNIVERSE) then CELL_TYPE
           ! ======================================================================
           ! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL
 
@@ -311,7 +222,7 @@ contains
           ! Move particle to next level and set universe
           j = j + 1
           p % n_coord = j
-          p % coord(j) % universe = c % fill
+          p % coord(j) % universe = c % fill() + 1
 
           ! Apply translation
           if (allocated(c % translation)) then
@@ -328,11 +239,11 @@ contains
           call find_cell(p, found)
           j = p % n_coord
 
-        elseif (c % type == FILL_LATTICE) then CELL_TYPE
+        elseif (c % type() == FILL_LATTICE) then CELL_TYPE
           ! ======================================================================
           ! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL
 
-          associate (lat => lattices(c % fill) % obj)
+          associate (lat => lattices(c % fill() + 1) % obj)
             ! Determine lattice indices
             i_xyz = lat % get_indices(p % coord(j) % xyz + TINY_BIT * p % coord(j) % uvw)
 
@@ -341,7 +252,7 @@ contains
             p % coord(j + 1) % uvw = p % coord(j) % uvw
 
             ! set particle lattice indices
-            p % coord(j + 1) % lattice   = c % fill
+            p % coord(j + 1) % lattice   = c % fill() + 1
             p % coord(j + 1) % lattice_x = i_xyz(1)
             p % coord(j + 1) % lattice_y = i_xyz(2)
             p % coord(j + 1) % lattice_z = i_xyz(3)
@@ -350,18 +261,18 @@ contains
             if (lat % are_valid_indices(i_xyz)) then
               ! Particle is inside the lattice.
               p % coord(j + 1) % universe = &
-                   lat % universes(i_xyz(1), i_xyz(2), i_xyz(3))
+                   lat % get([i_xyz(1), i_xyz(2), i_xyz(3)]) + 1
 
             else
               ! Particle is outside the lattice.
-              if (lat % outer == NO_OUTER_UNIVERSE) then
+              if (lat % outer() == NO_OUTER_UNIVERSE) then
                 call warning("Particle " // trim(to_str(p %id)) &
-                     // " is outside lattice " // trim(to_str(lat % id)) &
+                     // " is outside lattice " // trim(to_str(lat % id())) &
                      // " but the lattice has no defined outer universe.")
                 found = .false.
                 return
               else
-                p % coord(j + 1) % universe = lat % outer
+                p % coord(j + 1) % universe = lat % outer() + 1
               end if
             end if
           end associate
@@ -396,7 +307,7 @@ contains
     lat => lattices(p % coord(j) % lattice) % obj
 
     if (verbosity >= 10 .or. trace) then
-      call write_message("    Crossing lattice " // trim(to_str(lat % id)) &
+      call write_message("    Crossing lattice " // trim(to_str(lat % id())) &
            &// ". Current position (" // trim(to_str(p % coord(j) % lattice_x)) &
            &// "," // trim(to_str(p % coord(j) % lattice_y)) // "," &
            &// trim(to_str(p % coord(j) % lattice_z)) // ")")
@@ -428,7 +339,8 @@ contains
     else OUTSIDE_LAT
 
       ! Find cell in next lattice element
-      p % coord(j) % universe = lat % universes(i_xyz(1), i_xyz(2), i_xyz(3))
+      p % coord(j) % universe = &
+           lat % get([i_xyz(1), i_xyz(2), i_xyz(3)]) + 1
 
       call find_cell(p, found)
       if (.not. found) then
@@ -465,26 +377,15 @@ contains
     integer,        intent(out)   :: lattice_translation(3)
     integer,        intent(out)   :: next_level
 
-    integer :: i                  ! index for surface in cell
     integer :: j
-    integer :: index_surf         ! index in surfaces array (with sign)
     integer :: i_xyz(3)           ! lattice indices
     integer :: level_surf_cross   ! surface crossed on current level
     integer :: level_lat_trans(3) ! lattice translation on current level
-    real(8) :: x,y,z              ! particle coordinates
     real(8) :: xyz_t(3)           ! local particle coordinates
-    real(8) :: beta, gama         ! skewed particle coordiantes
-    real(8) :: u,v,w              ! particle directions
-    real(8) :: beta_dir           ! skewed particle direction
-    real(8) :: gama_dir           ! skewed particle direction
-    real(8) :: edge               ! distance to oncoming edge
-    real(8) :: d                  ! evaluated distance
     real(8) :: d_lat              ! distance to lattice boundary
     real(8) :: d_surf             ! distance to surface
-    real(8) :: x0,y0,z0           ! coefficients for surface
     real(8) :: xyz_cross(3)       ! coordinates at projected surface crossing
     real(8) :: surf_uvw(3)        ! surface normal direction
-    logical :: coincident         ! is particle on surface?
     type(Cell),       pointer :: c
     class(Lattice),   pointer :: lat
 
@@ -502,35 +403,11 @@ contains
       ! get pointer to cell on this level
       c => cells(p % coord(j) % cell)
 
-      ! copy directional cosines
-      u = p % coord(j) % uvw(1)
-      v = p % coord(j) % uvw(2)
-      w = p % coord(j) % uvw(3)
-
       ! =======================================================================
       ! FIND MINIMUM DISTANCE TO SURFACE IN THIS CELL
 
-      SURFACE_LOOP: do i = 1, size(c % region)
-        index_surf = c % region(i)
-        coincident = (index_surf == p % surface)
-
-        ! ignore this token if it corresponds to an operator rather than a
-        ! region.
-        index_surf = abs(index_surf)
-        if (index_surf >= OP_UNION) cycle
-
-        ! Calculate distance to surface
-        d = surfaces(index_surf) % distance(p % coord(j) % xyz, &
-             p % coord(j) % uvw, logical(coincident, kind=C_BOOL))
-
-        ! Check if calculated distance is new minimum
-        if (d < d_surf) then
-          if (abs(d - d_surf)/d_surf >= FP_PRECISION) then
-            d_surf = d
-            level_surf_cross = -c % region(i)
-          end if
-        end if
-      end do SURFACE_LOOP
+      call c % distance(p % coord(j) % xyz, p % coord(j) % uvw, p % surface, &
+                        d_surf, level_surf_cross)
 
       ! =======================================================================
       ! FIND MINIMUM DISTANCE TO LATTICE SURFACES
@@ -538,185 +415,22 @@ contains
       LAT_COORD: if (p % coord(j) % lattice /= NONE) then
         lat => lattices(p % coord(j) % lattice) % obj
 
+        i_xyz(1) = p % coord(j) % lattice_x
+        i_xyz(2) = p % coord(j) % lattice_y
+        i_xyz(3) = p % coord(j) % lattice_z
+
         LAT_TYPE: select type(lat)
 
         type is (RectLattice)
-          ! copy local coordinates
-          x = p % coord(j) % xyz(1)
-          y = p % coord(j) % xyz(2)
-          z = p % coord(j) % xyz(3)
-
-          ! determine oncoming edge
-          x0 = sign(lat % pitch(1) * HALF, u)
-          y0 = sign(lat % pitch(2) * HALF, v)
-
-          ! left and right sides
-          if (abs(x - x0) < FP_PRECISION) then
-            d = INFINITY
-          elseif (u == ZERO) then
-            d = INFINITY
-          else
-            d = (x0 - x)/u
-          end if
-
-          d_lat = d
-          if (u > 0) then
-            level_lat_trans(:) = [1, 0, 0]
-          else
-            level_lat_trans(:) = [-1, 0, 0]
-          end if
-
-          ! front and back sides
-          if (abs(y - y0) < FP_PRECISION) then
-            d = INFINITY
-          elseif (v == ZERO) then
-            d = INFINITY
-          else
-            d = (y0 - y)/v
-          end if
-
-          if (d < d_lat) then
-            d_lat = d
-            if (v > 0) then
-              level_lat_trans(:) = [0, 1, 0]
-            else
-              level_lat_trans(:) = [0, -1, 0]
-            end if
-          end if
-
-          if (lat % is_3d) then
-            z0 = sign(lat % pitch(3) * HALF, w)
-
-            ! top and bottom sides
-            if (abs(z - z0) < FP_PRECISION) then
-              d = INFINITY
-            elseif (w == ZERO) then
-              d = INFINITY
-            else
-              d = (z0 - z)/w
-            end if
-
-            if (d < d_lat) then
-              d_lat = d
-              if (w > 0) then
-                level_lat_trans(:) = [0, 0, 1]
-              else
-                level_lat_trans(:) = [0, 0, -1]
-              end if
-            end if
-          end if
+          call lat % distance(p % coord(j) % xyz, p % coord(j) % uvw, &
+                              i_xyz, d_lat, level_lat_trans)
 
         type is (HexLattice) LAT_TYPE
-          ! Copy local coordinates.
-          z = p % coord(j) % xyz(3)
-          i_xyz(1) = p % coord(j) % lattice_x
-          i_xyz(2) = p % coord(j) % lattice_y
-          i_xyz(3) = p % coord(j) % lattice_z
-
-          ! Compute velocities along the hexagonal axes.
-          beta_dir = u*sqrt(THREE)/TWO + v/TWO
-          gama_dir = u*sqrt(THREE)/TWO - v/TWO
-
-          ! Note that hexagonal lattice distance calculations are performed
-          ! using the particle's coordinates relative to the neighbor lattice
-          ! cells, not relative to the particle's current cell.  This is done
-          ! because there is significant disagreement between neighboring cells
-          ! on where the lattice boundary is due to the worse finite precision
-          ! of hex lattices.
-
-          ! Upper right and lower left sides.
-          edge = -sign(lat % pitch(1)/TWO, beta_dir)  ! Oncoming edge
-          if (beta_dir > ZERO) then
-            xyz_t = lat % get_local_xyz(p % coord(j - 1) % xyz, i_xyz+[1, 0, 0])
-          else
-            xyz_t = lat % get_local_xyz(p % coord(j - 1) % xyz, i_xyz+[-1, 0, 0])
-          end if
-          beta = xyz_t(1)*sqrt(THREE)/TWO + xyz_t(2)/TWO
-          if (abs(beta - edge) < FP_PRECISION) then
-            d = INFINITY
-          else if (beta_dir == ZERO) then
-            d = INFINITY
-          else
-            d = (edge - beta)/beta_dir
-          end if
-
-          d_lat = d
-          if (beta_dir > 0) then
-            level_lat_trans(:) = [1, 0, 0]
-          else
-            level_lat_trans(:) = [-1, 0, 0]
-          end if
-
-          ! Lower right and upper left sides.
-          edge = -sign(lat % pitch(1)/TWO, gama_dir)  ! Oncoming edge
-          if (gama_dir > ZERO) then
-            xyz_t = lat % get_local_xyz(p % coord(j - 1) % xyz, i_xyz+[1, -1, 0])
-          else
-            xyz_t = lat % get_local_xyz(p % coord(j - 1) % xyz, i_xyz+[-1, 1, 0])
-          end if
-          gama = xyz_t(1)*sqrt(THREE)/TWO - xyz_t(2)/TWO
-          if (abs(gama - edge) < FP_PRECISION) then
-            d = INFINITY
-          else if (gama_dir == ZERO) then
-            d = INFINITY
-          else
-            d = (edge - gama)/gama_dir
-          end if
-
-          if (d < d_lat) then
-            d_lat = d
-            if (gama_dir > 0) then
-              level_lat_trans(:) = [1, -1, 0]
-            else
-              level_lat_trans(:) = [-1, 1, 0]
-            end if
-          end if
-
-          ! Upper and lower sides.
-          edge = -sign(lat % pitch(1)/TWO, v)  ! Oncoming edge
-          if (v > ZERO) then
-            xyz_t = lat % get_local_xyz(p % coord(j - 1) % xyz, i_xyz+[0, 1, 0])
-          else
-            xyz_t = lat % get_local_xyz(p % coord(j - 1) % xyz, i_xyz+[0, -1, 0])
-          end if
-          if (abs(xyz_t(2) - edge) < FP_PRECISION) then
-            d = INFINITY
-          else if (v == ZERO) then
-            d = INFINITY
-          else
-            d = (edge - xyz_t(2))/v
-          end if
-
-          if (d < d_lat) then
-            d_lat = d
-            if (v > 0) then
-              level_lat_trans(:) = [0, 1, 0]
-            else
-              level_lat_trans(:) = [0, -1, 0]
-            end if
-          end if
-
-          ! Top and bottom sides.
-          if (lat % is_3d) then
-            z0 = sign(lat % pitch(2) * HALF, w)
-
-            if (abs(z - z0) < FP_PRECISION) then
-              d = INFINITY
-            elseif (w == ZERO) then
-              d = INFINITY
-            else
-              d = (z0 - z)/w
-            end if
-
-            if (d < d_lat) then
-              d_lat = d
-              if (w > 0) then
-                level_lat_trans(:) = [0, 0, 1]
-              else
-                level_lat_trans(:) = [0, 0, -1]
-              end if
-            end if
-          end if
+          xyz_t(1) = p % coord(j-1) % xyz(1)
+          xyz_t(2) = p % coord(j-1) % xyz(2)
+          xyz_t(3) = p % coord(j) % xyz(3)
+          call lat % distance(xyz_t, p % coord(j) % uvw, &
+                              i_xyz, d_lat, level_lat_trans)
         end select LAT_TYPE
 
         if (d_lat < ZERO) then
@@ -738,7 +452,7 @@ contains
           ! positive half-space were given in the region specification. Thus, we
           ! have to explicitly check which half-space the particle would be
           ! traveling into if the surface is crossed
-          if (.not. c % simple) then
+          if (.not. c % simple()) then
             xyz_cross(:) = p % coord(j) % xyz + d_surf*p % coord(j) % uvw
             call surfaces(abs(level_surf_cross)) % normal(xyz_cross, surf_uvw)
             if (dot_product(p % coord(j) % uvw, surf_uvw) > ZERO) then
@@ -820,388 +534,4 @@ contains
 
   end subroutine neighbor_lists
 
-!===============================================================================
-! CALC_OFFSETS calculates and stores the offsets in all fill cells. This
-! routine is called once upon initialization.
-!===============================================================================
-
-  subroutine calc_offsets(univ_id, map, univ, counts, found)
-
-    integer, intent(in)        :: univ_id         ! target universe ID
-    integer, intent(in)        :: map          ! map index in vector of maps
-    type(Universe), intent(in) :: univ         ! universe searching in
-    integer, intent(inout)     :: counts(:,:)  ! target count
-    logical, intent(inout)     :: found(:,:)   ! target found
-
-    integer :: i                          ! index over cells
-    integer :: j, k, m                    ! indices in lattice
-    integer :: n                          ! number of cells to search
-    integer :: offset                     ! total offset for a given cell
-    integer :: cell_index                 ! index in cells array
-    type(Cell),     pointer :: c          ! pointer to current cell
-    type(Universe), pointer :: next_univ  ! next universe to cycle through
-    class(Lattice), pointer :: lat        ! pointer to current lattice
-
-    n = size(univ % cells)
-    offset = 0
-
-    do i = 1, n
-
-      cell_index = univ % cells(i)
-
-      ! get pointer to cell
-      c => cells(cell_index)
-
-      ! ====================================================================
-      ! AT LOWEST UNIVERSE, TERMINATE SEARCH
-      if (c % type == FILL_MATERIAL) then
-
-      ! ====================================================================
-      ! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL
-      elseif (c % type == FILL_UNIVERSE) then
-        ! Set offset for the cell on this level
-        c % offset(map) = offset
-
-        ! Count contents of this cell
-        next_univ => universes(c % fill)
-        offset = offset + count_target(next_univ, counts, found, univ_id, map)
-
-        ! Move into the next universe
-        next_univ => universes(c % fill)
-        c => cells(cell_index)
-
-      ! ====================================================================
-      ! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL
-      elseif (c % type == FILL_LATTICE) then
-
-        ! Set current lattice
-        lat => lattices(c % fill) % obj
-
-        select type (lat)
-
-        type is (RectLattice)
-
-          ! Loop over lattice coordinates
-          do m = 1, lat % n_cells(3)
-            do k = 1, lat % n_cells(2)
-              do j = 1, lat % n_cells(1)
-                lat % offset(map, j, k, m) = offset
-                next_univ => universes(lat % universes(j, k, m))
-                offset = offset + &
-                     count_target(next_univ, counts, found, univ_id, map)
-              end do
-            end do
-          end do
-
-        type is (HexLattice)
-
-          ! Loop over lattice coordinates
-          do m = 1, lat % n_axial
-            do k = 1, 2*lat % n_rings - 1
-              do j = 1, 2*lat % n_rings - 1
-                ! This array location is never used
-                if (j + k < lat % n_rings + 1) then
-                  cycle
-                ! This array location is never used
-                else if (j + k > 3*lat % n_rings - 1) then
-                  cycle
-                else
-                  lat % offset(map, j, k, m) = offset
-                  next_univ => universes(lat % universes(j, k, m))
-                  offset = offset + &
-                       count_target(next_univ, counts, found, univ_id, map)
-                end if
-              end do
-            end do
-          end do
-        end select
-
-      end if
-    end do
-
-  end subroutine calc_offsets
-
-!===============================================================================
-! COUNT_TARGET recursively totals the numbers of occurances of a given
-! universe ID beginning with the universe given.
-!===============================================================================
-
-  recursive function count_target(univ, counts, found, univ_id, map) result(count)
-
-    type(Universe), intent(in) :: univ         ! universe to search through
-    integer, intent(inout)     :: counts(:,:)  ! target count
-    logical, intent(inout)     :: found(:,:)   ! target found
-    integer, intent(in)        :: univ_id         ! target universe ID
-    integer, intent(in)        :: map          ! current map
-
-    integer :: i                           ! index over cells
-    integer :: j, k, m                     ! indices in lattice
-    integer :: n                           ! number of cells to search
-    integer :: cell_index                  ! index in cells array
-    integer :: count                       ! number of times target located
-    type(Cell),     pointer :: c           ! pointer to current cell
-    type(Universe), pointer :: next_univ   ! next univ to loop through
-    class(Lattice), pointer :: lat         ! pointer to current lattice
-
-    ! Don't research places already checked
-    if (found(universe_dict % get(univ % id), map)) then
-      count = counts(universe_dict % get(univ % id), map)
-      return
-    end if
-
-    ! If this is the target, it can't contain itself.
-    ! Count = 1, then quit
-    if (univ % id == univ_id) then
-      count = 1
-      counts(universe_dict % get(univ % id), map) = 1
-      found(universe_dict % get(univ % id), map) = .true.
-      return
-    end if
-
-    count = 0
-    n = size(univ % cells)
-
-    do i = 1, n
-
-      cell_index = univ % cells(i)
-
-      ! get pointer to cell
-      c => cells(cell_index)
-
-      ! ====================================================================
-      ! AT LOWEST UNIVERSE, TERMINATE SEARCH
-      if (c % type == FILL_MATERIAL) then
-
-      ! ====================================================================
-      ! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL
-      elseif (c % type == FILL_UNIVERSE) then
-
-        next_univ => universes(c % fill)
-
-        ! Found target - stop since target cannot contain itself
-        if (next_univ % id == univ_id) then
-          count = count + 1
-          return
-        end if
-
-        count = count + count_target(next_univ, counts, found, univ_id, map)
-        c => cells(cell_index)
-
-      ! ====================================================================
-      ! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL
-      elseif (c % type == FILL_LATTICE) then
-
-        ! Set current lattice
-        lat => lattices(c % fill) % obj
-
-        select type (lat)
-
-        type is (RectLattice)
-
-          ! Loop over lattice coordinates
-          do m = 1, lat % n_cells(3)
-            do k = 1, lat % n_cells(2)
-              do j = 1, lat % n_cells(1)
-                next_univ => universes(lat % universes(j, k, m))
-
-                ! Found target - stop since target cannot contain itself
-                if (next_univ % id == univ_id) then
-                  count = count + 1
-                  cycle
-                end if
-
-                count = count + &
-                     count_target(next_univ, counts, found, univ_id, map)
-
-              end do
-            end do
-          end do
-
-          type is (HexLattice)
-
-            ! Loop over lattice coordinates
-            do m = 1, lat % n_axial
-              do k = 1, 2*lat % n_rings - 1
-                do j = 1, 2*lat % n_rings - 1
-                  ! This array location is never used
-                  if (j + k < lat % n_rings + 1) then
-                    cycle
-                  ! This array location is never used
-                  else if (j + k > 3*lat % n_rings - 1) then
-                    cycle
-                  else
-                    next_univ => universes(lat % universes(j, k, m))
-
-                    ! Found target - stop since target cannot contain itself
-                    if (next_univ % id == univ_id) then
-                      count = count + 1
-                      cycle
-                    end if
-
-                    count = count + &
-                         count_target(next_univ, counts, found, univ_id, map)
-                  end if
-                end do
-              end do
-            end do
-
-          end select
-
-      end if
-    end do
-
-    counts(universe_dict % get(univ % id), map) = count
-    found(universe_dict % get(univ % id), map) = .true.
-
-  end function count_target
-
-!===============================================================================
-! COUNT_INSTANCE recursively totals the number of occurrences of all cells
-! beginning with the universe given.
-!===============================================================================
-
-  recursive subroutine count_instance(univ)
-
-    type(Universe), intent(in) :: univ  ! universe to search through
-
-    integer :: i        ! index over cells
-    integer :: j, k, m  ! indices in lattice
-    integer :: n        ! number of cells to search
-
-    n = size(univ % cells)
-
-    do i = 1, n
-      associate (c => cells(univ % cells(i)))
-        c % instances = c % instances + 1
-
-        ! ====================================================================
-        ! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL
-        if (c % type == FILL_UNIVERSE) then
-
-          call count_instance(universes(c % fill))
-
-        ! ====================================================================
-        ! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL
-        elseif (c % type == FILL_LATTICE) then
-
-          ! Set current lattice
-          associate (lat => lattices(c % fill) % obj)
-
-            select type (lat)
-            type is (RectLattice)
-
-              ! Loop over lattice coordinates
-              do m = 1, lat % n_cells(3)
-                do k = 1, lat % n_cells(2)
-                  do j = 1, lat % n_cells(1)
-                    call count_instance(universes(lat % universes(j, k, m)))
-                  end do
-                end do
-              end do
-
-            type is (HexLattice)
-
-              ! Loop over lattice coordinates
-              do m = 1, lat % n_axial
-                do k = 1, 2*lat % n_rings - 1
-                  do j = 1, 2*lat % n_rings - 1
-                    ! This array location is never used
-                    if (j + k < lat % n_rings + 1) then
-                      cycle
-                      ! This array location is never used
-                    else if (j + k > 3*lat % n_rings - 1) then
-                      cycle
-                    else
-                      call count_instance(universes(lat % universes(j, k, m)))
-                    end if
-                  end do
-                end do
-              end do
-
-            end select
-          end associate
-        end if
-      end associate
-    end do
-
-  end subroutine count_instance
-
-!===============================================================================
-! MAXIMUM_LEVELS determines the maximum number of nested coordinate levels in
-! the geometry
-!===============================================================================
-
-  recursive function maximum_levels(univ) result(levels)
-
-    type(Universe), intent(in) :: univ  ! universe to search through
-    integer :: levels                   ! maximum number of levels for this universe
-
-    integer :: i                          ! index over cells
-    integer :: j, k, m                    ! indices in lattice
-    integer :: levels_below               ! max levels below this universe
-    type(Cell),     pointer :: c          ! pointer to current cell
-    type(Universe), pointer :: next_univ  ! next universe to loop through
-    class(Lattice), pointer :: lat        ! pointer to current lattice
-
-    levels_below = 0
-    do i = 1, size(univ % cells)
-      c => cells(univ % cells(i))
-
-      ! ====================================================================
-      ! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL
-      if (c % type == FILL_UNIVERSE) then
-
-        next_univ => universes(c % fill)
-        levels_below = max(levels_below, maximum_levels(next_univ))
-
-      ! ====================================================================
-      ! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL
-      elseif (c % type == FILL_LATTICE) then
-
-        ! Set current lattice
-        lat => lattices(c % fill) % obj
-
-        select type (lat)
-
-        type is (RectLattice)
-
-          ! Loop over lattice coordinates
-          do m = 1, lat % n_cells(3)
-            do k = 1, lat % n_cells(2)
-              do j = 1, lat % n_cells(1)
-                next_univ => universes(lat % universes(j, k, m))
-                levels_below = max(levels_below, maximum_levels(next_univ))
-              end do
-            end do
-          end do
-
-        type is (HexLattice)
-
-          ! Loop over lattice coordinates
-          do m = 1, lat % n_axial
-            do k = 1, 2*lat % n_rings - 1
-              do j = 1, 2*lat % n_rings - 1
-                ! This array location is never used
-                if (j + k < lat % n_rings + 1) then
-                  cycle
-                ! This array location is never used
-                else if (j + k > 3*lat % n_rings - 1) then
-                  cycle
-                else
-                  next_univ => universes(lat % universes(j, k, m))
-                  levels_below = max(levels_below, maximum_levels(next_univ))
-                end if
-              end do
-            end do
-          end do
-
-        end select
-
-      end if
-    end do
-
-    levels = 1 + levels_below
-
-  end function maximum_levels
-
 end module geometry

src/geometry_aux.cpp

@@ -0,0 +1,341 @@
+#include "geometry_aux.h"
+
+#include <algorithm>  // for std::max
+#include <sstream>
+#include <unordered_set>
+
+#include "cell.h"
+#include "constants.h"
+#include "error.h"
+#include "lattice.h"
+
+
+namespace openmc {
+
+//==============================================================================
+
+void
+adjust_indices_c()
+{
+  // Adjust material/fill idices.
+  for (Cell *c : global_cells) {
+    if (c->material[0] == C_NONE) {
+      int32_t id = c->fill;
+      auto search_univ = universe_map.find(id);
+      auto search_lat = lattice_map.find(id);
+      if (search_univ != universe_map.end()) {
+        c->type = FILL_UNIVERSE;
+        c->fill = search_univ->second;
+      } else if (search_lat != lattice_map.end()) {
+        c->type = FILL_LATTICE;
+        c->fill = search_lat->second;
+      } else {
+        std::stringstream err_msg;
+        err_msg << "Specified fill " << id << " on cell " << c->id
+                << " is neither a universe nor a lattice.";
+        fatal_error(err_msg);
+      }
+    } else {
+      //TODO: materials
+      c->type = FILL_MATERIAL;
+    }
+  }
+
+  // Change cell.universe values from IDs to indices.
+  for (Cell *c : global_cells) {
+    auto search = universe_map.find(c->universe);
+    if (search != universe_map.end()) {
+      //TODO: Remove this off-by-one indexing.
+      c->universe = search->second + 1;
+    } else {
+      std::stringstream err_msg;
+      err_msg << "Could not find universe " << c->universe
+              << " specified on cell " << c->id;
+      fatal_error(err_msg);
+    }
+  }
+
+  // Change all lattice universe values from IDs to indices.
+  for (Lattice *l : lattices_c) {
+    l->adjust_indices();
+  }
+}
+
+//==============================================================================
+
+int32_t
+find_root_universe()
+{
+  // Find all the universes listed as a cell fill.
+  std::unordered_set<int32_t> fill_univ_ids;
+  for (Cell *c : global_cells) {
+    fill_univ_ids.insert(c->fill);
+  }
+
+  // Find all the universes contained in a lattice.
+  for (Lattice *lat : lattices_c) {
+    for (auto it = lat->begin(); it != lat->end(); ++it) {
+      fill_univ_ids.insert(*it);
+    }
+    if (lat->outer != NO_OUTER_UNIVERSE) {
+      fill_univ_ids.insert(lat->outer);
+    }
+  }
+
+  // Figure out which universe is not in the set.  This is the root universe.
+  bool root_found {false};
+  int32_t root_univ;
+  for (int32_t i = 0; i < global_universes.size(); i++) {
+    auto search = fill_univ_ids.find(global_universes[i]->id);
+    if (search == fill_univ_ids.end()) {
+      if (root_found) {
+        fatal_error("Two or more universes are not used as fill universes, so "
+                    "it is not possible to distinguish which one is the root "
+                    "universe.");
+      } else {
+        root_found = true;
+        root_univ = i;
+      }
+    }
+  }
+  if (!root_found) fatal_error("Could not find a root universe.  Make sure "
+       "there are no circular dependencies in the geometry.");
+
+  return root_univ;
+}
+
+//==============================================================================
+
+void
+allocate_offset_tables(int n_maps)
+{
+  for (Cell *c : global_cells) {
+    if (c->type != FILL_MATERIAL) {
+      c->offset.resize(n_maps, C_NONE);
+    }
+  }
+
+  for (Lattice *lat : lattices_c) {
+    lat->allocate_offset_table(n_maps);
+  }
+}
+
+//==============================================================================
+
+void
+count_cell_instances(int32_t univ_indx)
+{
+  for (int32_t cell_indx : global_universes[univ_indx]->cells) {
+    Cell &c = *global_cells[cell_indx];
+    ++c.n_instances;
+
+    if (c.type == FILL_UNIVERSE) {
+      // This cell contains another universe.  Recurse into that universe.
+      count_cell_instances(c.fill);
+
+    } else if (c.type == FILL_LATTICE) {
+      // This cell contains a lattice.  Recurse into the lattice universes.
+      Lattice &lat = *lattices_c[c.fill];
+      for (auto it = lat.begin(); it != lat.end(); ++it) {
+        count_cell_instances(*it);
+      }
+    }
+  }
+}
+
+//==============================================================================
+
+int
+count_universe_instances(int32_t search_univ, int32_t target_univ_id)
+{
+  //  If this is the target, it can't contain itself.
+  if (global_universes[search_univ]->id == target_univ_id) {
+    return 1;
+  }
+
+  int count {0};
+  for (int32_t cell_indx : global_universes[search_univ]->cells) {
+    Cell &c = *global_cells[cell_indx];
+
+    if (c.type == FILL_UNIVERSE) {
+      int32_t next_univ = c.fill;
+      count += count_universe_instances(next_univ, target_univ_id);
+
+    } else if (c.type == FILL_LATTICE) {
+      Lattice &lat = *lattices_c[c.fill];
+      for (auto it = lat.begin(); it != lat.end(); ++it) {
+        int32_t next_univ = *it;
+        count += count_universe_instances(next_univ, target_univ_id);
+      }
+    }
+  }
+
+  return count;
+}
+
+//==============================================================================
+
+void
+fill_offset_tables(int32_t target_univ_id, int map)
+{
+  for (Universe *univ : global_universes) {
+    int32_t offset {0};  // TODO: is this a bug?  It matches F90 implementation.
+    for (int32_t cell_indx : univ->cells) {
+      Cell &c = *global_cells[cell_indx];
+
+      if (c.type == FILL_UNIVERSE) {
+        c.offset[map] = offset;
+        int32_t search_univ = c.fill;
+        offset += count_universe_instances(search_univ, target_univ_id);
+
+      } else if (c.type == FILL_LATTICE) {
+        Lattice &lat = *lattices_c[c.fill];
+        offset = lat.fill_offset_table(offset, target_univ_id, map);
+      }
+    }
+  }
+}
+
+//==============================================================================
+
+std::string
+distribcell_path_inner(int32_t target_cell, int32_t map, int32_t target_offset,
+                       const Universe &search_univ, int32_t offset)
+{
+  std::stringstream path;
+
+  path << "u" << search_univ.id << "->";
+
+  // Check to see if this universe directly contains the target cell.  If so,
+  // write to the path and return.
+  for (int32_t cell_indx : search_univ.cells) {
+    if ((cell_indx == target_cell) && (offset == target_offset)) {
+      Cell &c = *global_cells[cell_indx];
+      path << "c" << c.id;
+      return path.str();
+    }
+  }
+
+  // The target must be further down the geometry tree and contained in a fill
+  // cell or lattice cell in this universe.  Find which cell contains the
+  // target.
+  std::vector<std::int32_t>::const_reverse_iterator cell_it
+       {search_univ.cells.crbegin()};
+  for (; cell_it != search_univ.cells.crend(); ++cell_it) {
+    Cell &c = *global_cells[*cell_it];
+
+    // Material cells don't contain other cells so ignore them.
+    if (c.type != FILL_MATERIAL) {
+      int32_t temp_offset;
+      if (c.type == FILL_UNIVERSE) {
+        temp_offset = offset + c.offset[map];
+      } else {
+        Lattice &lat = *lattices_c[c.fill];
+        int32_t indx = lat.universes.size()*map + lat.begin().indx;
+        temp_offset = offset + lat.offsets[indx];
+      }
+
+      // The desired cell is the first cell that gives an offset smaller or
+      // equal to the target offset.
+      if (temp_offset <= target_offset) break;
+    }
+  }
+
+  // Add the cell to the path string.
+  Cell &c = *global_cells[*cell_it];
+  path << "c" << c.id << "->";
+
+  if (c.type == FILL_UNIVERSE) {
+    // Recurse into the fill cell.
+    offset += c.offset[map];
+    path << distribcell_path_inner(target_cell, map, target_offset,
+                                   *global_universes[c.fill], offset);
+    return path.str();
+  } else {
+    // Recurse into the lattice cell.
+    Lattice &lat = *lattices_c[c.fill];
+    path << "l" << lat.id;
+    for (ReverseLatticeIter it = lat.rbegin(); it != lat.rend(); ++it) {
+      int32_t indx = lat.universes.size()*map + it.indx;
+      int32_t temp_offset = offset + lat.offsets[indx];
+      if (temp_offset <= target_offset) {
+        offset = temp_offset;
+        path << "(" << lat.index_to_string(it.indx) << ")->";
+        path << distribcell_path_inner(target_cell, map, target_offset,
+                                       *global_universes[*it], offset);
+        return path.str();
+      }
+    }
+  }
+}
+
+//==============================================================================
+
+int
+distribcell_path_len(int32_t target_cell, int32_t map, int32_t target_offset,
+                     int32_t root_univ)
+{
+  Universe &root = *global_universes[root_univ];
+  std::string path_ {distribcell_path_inner(target_cell, map, target_offset,
+                                            root, 0)};
+  return path_.size() + 1;
+}
+
+//==============================================================================
+
+void
+distribcell_path(int32_t target_cell, int32_t map, int32_t target_offset,
+                 int32_t root_univ, char *path)
+{
+  Universe &root = *global_universes[root_univ];
+  std::string path_ {distribcell_path_inner(target_cell, map, target_offset,
+                                            root, 0)};
+  path_.copy(path, path_.size());
+  path[path_.size()] = '\0';
+}
+
+//==============================================================================
+
+int
+maximum_levels(int32_t univ)
+{
+  int levels_below {0};
+
+  for (int32_t cell_indx : global_universes[univ]->cells) {
+    Cell &c = *global_cells[cell_indx];
+    if (c.type == FILL_UNIVERSE) {
+      int32_t next_univ = c.fill;
+      levels_below = std::max(levels_below, maximum_levels(next_univ));
+    } else if (c.type == FILL_LATTICE) {
+      Lattice &lat = *lattices_c[c.fill];
+      for (auto it = lat.begin(); it != lat.end(); ++it) {
+        int32_t next_univ = *it;
+        levels_below = std::max(levels_below, maximum_levels(next_univ));
+      }
+    }
+  }
+
+  ++levels_below;
+  return levels_below;
+}
+
+//==============================================================================
+
+void
+free_memory_geometry_c()
+{
+  for (Cell *c : global_cells) {delete c;}
+  global_cells.clear();
+  cell_map.clear();
+  n_cells = 0;
+
+  for (Universe *u : global_universes) {delete u;}
+  global_universes.clear();
+  universe_map.clear();
+
+  for (Lattice *lat : lattices_c) {delete lat;}
+  lattices_c.clear();
+  lattice_map.clear();
+}
+
+} // namespace openmc

src/geometry_aux.h

@@ -0,0 +1,112 @@
+//! \file geometry_aux.h
+//! Auxilary functions for geometry initialization and general data handling.
+
+#ifndef GEOMETRY_AUX_H
+#define GEOMETRY_AUX_H
+
+#include <cstdint>
+
+
+namespace openmc {
+
+//==============================================================================
+//! Replace Universe, Lattice, and Material IDs with indices.
+//==============================================================================
+
+extern "C" void adjust_indices_c();
+
+//==============================================================================
+//! Figure out which Universe is the root universe.
+//!
+//! This function looks for a universe that is not listed in a Cell::fill or in
+//! a Lattice.
+//! @return The index of the root universe.
+//==============================================================================
+
+extern "C" int32_t find_root_universe();
+
+//==============================================================================
+//! Allocate storage in Lattice and Cell objects for distribcell offset tables.
+//==============================================================================
+
+extern "C" void allocate_offset_tables(int n_maps);
+
+//==============================================================================
+//! Recursively search through the geometry and count cell instances.
+//!
+//! This function will update the Cell::n_instances value for each cell in the
+//! geometry.
+//! @param univ_indx The index of the universe to begin searching from (probably
+//!   the root universe).
+//==============================================================================
+
+extern "C" void count_cell_instances(int32_t univ_indx);
+
+//==============================================================================
+//! Recursively search through universes and count universe instances.
+//! @param search_univ The index of the universe to begin searching from.
+//! @param target_univ_id The ID of the universe to be counted.
+//! @return The number of instances of target_univ_id in the geometry tree under
+//!   search_univ.
+//==============================================================================
+
+extern "C" int
+count_universe_instances(int32_t search_univ, int32_t target_univ_id);
+
+//==============================================================================
+//! Populate Cell and Lattice distribcell offset tables.
+//! @param target_univ_id The ID of the universe to be counted.
+//! @param map The index of the distribcell map that defines the offsets for the
+//!   target universe.
+//==============================================================================
+
+extern "C" void fill_offset_tables(int32_t target_univ_id, int map);
+
+//==============================================================================
+//! Find the length necessary for a string to contain a distribcell path.
+//! @param target_cell The index of the Cell in the global Cell array.
+//! @param map The index of the distribcell mapping corresponding to the target
+//!   cell.
+//! @param target_offset An instance number for a distributed cell.
+//! @param root_univ The index of the root Universe in the global Universe
+//!   array.
+//! @return The size of a character array needed to fit the distribcell path.
+//==============================================================================
+
+extern "C" int
+distribcell_path_len(int32_t target_cell, int32_t map, int32_t target_offset,
+                     int32_t root_univ);
+
+//==============================================================================
+//! Build a character array representing the path to a distribcell instance.
+//! @param target_cell The index of the Cell in the global Cell array.
+//! @param map The index of the distribcell mapping corresponding to the target
+//!   cell.
+//! @param target_offset An instance number for a distributed cell.
+//! @param root_univ The index of the root Universe in the global Universe
+//!   array.
+//! @param[out] path The unique traversal through the geometry tree that leads
+//!   to the desired instance of the target cell.
+//==============================================================================
+
+extern "C" void
+distribcell_path(int32_t target_cell, int32_t map, int32_t target_offset,
+                 int32_t root_univ, char *path);
+
+//==============================================================================
+//! Determine the maximum number of nested coordinate levels in the geometry.
+//! @param univ The index of the universe to begin seraching from (probably the
+//!   root universe).
+//! @return The number of coordinate levels.
+//==============================================================================
+
+extern "C" int maximum_levels(int32_t univ);
+
+//==============================================================================
+//! Deallocates global vectors and maps for cells, universes, and lattices.
+//==============================================================================
+
+extern "C" void free_memory_geometry_c();
+
+} // namespace openmc
+#endif // GEOMETRY_AUX_H

src/geometry_header.F90

@@ -6,6 +6,7 @@ module geometry_header
   use constants,       only: HALF, TWO, THREE, INFINITY, K_BOLTZMANN, &
                              MATERIAL_VOID, NONE
   use dict_header,     only: DictCharInt, DictIntInt
+  use hdf5_interface,  only: HID_T
   use material_header, only: Material, materials, material_dict, n_materials
   use nuclide_header
   use sab_header
@@ -14,13 +15,192 @@ module geometry_header
 
   implicit none
 
+  interface
+    function cell_pointer_c(cell_ind) bind(C, name='cell_pointer') result(ptr)
+      import C_PTR, C_INT32_T
+      integer(C_INT32_T), intent(in), value :: cell_ind
+      type(C_PTR)                           :: ptr
+    end function cell_pointer_c
+
+    function cell_id_c(cell_ptr) bind(C, name='cell_id') result(id)
+      import C_PTR, C_INT32_T
+      type(C_PTR), intent(in), value :: cell_ptr
+      integer(C_INT32_T)             :: id
+    end function cell_id_c
+
+    subroutine cell_set_id_c(cell_ptr, id) bind(C, name='cell_set_id')
+      import C_PTR, C_INT32_T
+      type(C_PTR),        intent(in), value :: cell_ptr
+      integer(C_INT32_T), intent(in), value :: id
+    end subroutine cell_set_id_c
+
+    function cell_type_c(cell_ptr) bind(C, name='cell_type') result(type)
+      import C_PTR, C_INT
+      type(C_PTR), intent(in), value :: cell_ptr
+      integer(C_INT)                 :: type
+    end function cell_type_c
+
+    subroutine cell_set_type_c(cell_ptr, type) bind(C, name='cell_set_type')
+      import C_PTR, C_INT
+      type(C_PTR),    intent(in), value :: cell_ptr
+      integer(C_INT), intent(in), value :: type
+    end subroutine cell_set_type_c
+
+    function cell_universe_c(cell_ptr) bind(C, name='cell_universe') &
+         result(universe)
+      import C_PTR, C_INT32_T
+      type(C_PTR), intent(in), value :: cell_ptr
+      integer(C_INT32_T)             :: universe
+    end function cell_universe_c
+
+    subroutine cell_set_universe_c(cell_ptr, universe) &
+         bind(C, name='cell_set_universe')
+      import C_PTR, C_INT32_T
+      type(C_PTR),        intent(in), value :: cell_ptr
+      integer(C_INT32_T), intent(in), value :: universe
+    end subroutine cell_set_universe_c
+
+    function cell_fill_c(cell_ptr) bind(C, name="cell_fill") result(fill)
+      import C_PTR, C_INT32_T
+      type(C_PTR), intent(in), value :: cell_ptr
+      integer(C_INT32_T)             :: fill
+    end function cell_fill_c
+
+    function cell_fill_ptr(cell_ptr) bind(C) result(fill_ptr)
+      import C_PTR
+      type(C_PTR), intent(in), value :: cell_ptr
+      type(C_PTR)                    :: fill_ptr
+    end function cell_fill_ptr
+
+    function cell_n_instances_c(cell_ptr) bind(C, name='cell_n_instances') &
+         result(n_instances)
+      import C_PTR, C_INT32_T
+      type(C_PTR), intent(in), value :: cell_ptr
+      integer(C_INT32_T)             :: n_instances
+    end function cell_n_instances_c
+
+    function cell_simple_c(cell_ptr) bind(C, name='cell_simple') result(simple)
+      import C_PTR, C_BOOL
+      type(C_PTR), intent(in), value :: cell_ptr
+      logical(C_BOOL)                :: simple
+    end function cell_simple_c
+
+    subroutine cell_distance_c(cell_ptr, xyz, uvw, on_surface, min_dist, &
+                               i_surf) bind(C, name="cell_distance")
+      import C_PTR, C_INT32_T, C_DOUBLE
+      type(C_PTR),        intent(in), value :: cell_ptr
+      real(C_DOUBLE),     intent(in)        :: xyz(3)
+      real(C_DOUBLE),     intent(in)        :: uvw(3)
+      integer(C_INT32_T), intent(in), value :: on_surface
+      real(C_DOUBLE),     intent(out)       :: min_dist
+      integer(C_INT32_T), intent(out)       :: i_surf
+    end subroutine cell_distance_c
+
+    function cell_offset_c(cell_ptr, map) bind(C, name="cell_offset") &
+         result(offset)
+      import C_PTR, C_INT, C_INT32_T
+      type(C_PTR),    intent(in), value :: cell_ptr
+      integer(C_INT), intent(in), value :: map
+      integer(C_INT32_T)                :: offset
+    end function cell_offset_c
+
+    subroutine cell_to_hdf5_c(cell_ptr, group) bind(C, name='cell_to_hdf5')
+      import HID_T, C_PTR
+      type(C_PTR),    intent(in), value :: cell_ptr
+      integer(HID_T), intent(in), value :: group
+    end subroutine cell_to_hdf5_c
+
+    function lattice_pointer_c(lat_ind) bind(C, name='lattice_pointer') &
+         result(ptr)
+      import C_PTR, C_INT32_T
+      integer(C_INT32_T), intent(in), value :: lat_ind
+      type(C_PTR)                           :: ptr
+    end function lattice_pointer_c
+
+    function lattice_id_c(lat_ptr) bind(C, name='lattice_id') result(id)
+      import C_PTR, C_INT32_T
+      type(C_PTR), intent(in), value :: lat_ptr
+      integer(C_INT32_T)             :: id
+    end function lattice_id_c
+
+    function lattice_are_valid_indices_c(lat_ptr, i_xyz) &
+         bind(C, name='lattice_are_valid_indices') result (is_valid)
+      import C_PTR, C_INT, C_BOOL
+      type(C_PTR),    intent(in), value :: lat_ptr
+      integer(C_INT), intent(in)        :: i_xyz(3)
+      logical(C_BOOL)                   :: is_valid
+    end function lattice_are_valid_indices_c
+
+    subroutine lattice_distance_c(lat_ptr, xyz, uvw, i_xyz, d, lattice_trans) &
+         bind(C, name='lattice_distance')
+      import C_PTR, C_INT, C_DOUBLE
+      type(C_PTR),    intent(in), value :: lat_ptr
+      real(C_DOUBLE), intent(in)        :: xyz(3)
+      real(C_DOUBLE), intent(in)        :: uvw(3)
+      integer(C_INT), intent(in)        :: i_xyz(3)
+      real(C_DOUBLE), intent(out)       :: d
+      integer(C_INT), intent(out)       :: lattice_trans(3)
+    end subroutine lattice_distance_c
+
+    subroutine lattice_get_indices_c(lat_ptr, xyz, i_xyz) &
+         bind(C, name='lattice_get_indices')
+      import C_PTR, C_INT, C_DOUBLE
+      type(C_PTR),    intent(in), value :: lat_ptr
+      real(C_DOUBLE), intent(in)        :: xyz(3)
+      integer(C_INT), intent(out)       :: i_xyz(3)
+    end subroutine lattice_get_indices_c
+
+    subroutine lattice_get_local_xyz_c(lat_ptr, global_xyz, i_xyz, local_xyz) &
+         bind(C, name='lattice_get_local_xyz')
+      import C_PTR, C_INT, C_DOUBLE
+      type(C_PTR),    intent(in), value :: lat_ptr
+      real(C_DOUBLE), intent(in)        :: global_xyz(3)
+      integer(C_INT), intent(in)        :: i_xyz(3)
+      real(C_DOUBLE), intent(out)       :: local_xyz(3)
+    end subroutine lattice_get_local_xyz_c
+
+    subroutine lattice_to_hdf5_c(lat_ptr, group) bind(C, name='lattice_to_hdf5')
+      import HID_T, C_PTR
+      type(C_PTR),    intent(in), value :: lat_ptr
+      integer(HID_T), intent(in), value :: group
+    end subroutine lattice_to_hdf5_c
+
+    function lattice_offset_c(lat_ptr, map, i_xyz) &
+         bind(C, name='lattice_offset') result(offset)
+      import C_PTR, C_INT, C_INT32_T
+      type(C_PTR),    intent(in), value :: lat_ptr
+      integer(C_INT), intent(in), value :: map
+      integer(C_INT), intent(in)        :: i_xyz(3)
+      integer(C_INT32_T)                :: offset
+    end function lattice_offset_c
+
+    function lattice_outer_c(lat_ptr) bind(C, name='lattice_outer') &
+         result(outer)
+      import C_PTR, C_INT32_T
+      type(C_PTR), intent(in), value :: lat_ptr
+      integer(C_INT32_T)             :: outer
+    end function lattice_outer_c
+
+    function lattice_universe_c(lat_ptr, i_xyz) &
+         bind(C, name='lattice_universe') result(univ)
+      import C_PTR, C_INT32_t, C_INT
+      type(C_PTR),    intent(in), value :: lat_ptr
+      integer(C_INT), intent(in)        :: i_xyz(3)
+      integer(C_INT32_T)                :: univ
+    end function lattice_universe_c
+
+    subroutine extend_cells_c(n) bind(C)
+      import C_INT32_t
+      integer(C_INT32_T), intent(in), value :: n
+    end subroutine extend_cells_c
+  end interface
+
 !===============================================================================
 ! UNIVERSE defines a geometry that fills all phase space
 !===============================================================================
 
   type Universe
     integer :: id                     ! Unique ID
-    integer :: type                   ! Type
     integer, allocatable :: cells(:)  ! List of cells within
     real(8) :: x0                     ! Translation in x-coordinate
     real(8) :: y0                     ! Translation in y-coordinate
@@ -32,68 +212,24 @@ module geometry_header
 !===============================================================================
 
   type, abstract :: Lattice
-    integer              :: id               ! Universe number for lattice
-    character(len=104)   :: name = ""        ! User-defined name
-    real(8), allocatable :: pitch(:)         ! Pitch along each axis
-    integer, allocatable :: universes(:,:,:) ! Specified universes
-    integer              :: outside          ! Material to fill area outside
-    integer              :: outer            ! universe to tile outside the lat
-    logical              :: is_3d            ! Lattice has cells on z axis
-    integer, allocatable :: offset(:,:,:,:)  ! Distribcell offsets
+    type(C_PTR) :: ptr
   contains
-    procedure(lattice_are_valid_indices_), deferred :: are_valid_indices
-    procedure(lattice_get_indices_),       deferred :: get_indices
-    procedure(lattice_get_local_xyz_),     deferred :: get_local_xyz
+    procedure :: id => lattice_id
+    procedure :: are_valid_indices => lattice_are_valid_indices
+    procedure :: distance => lattice_distance
+    procedure :: get => lattice_get
+    procedure :: get_indices => lattice_get_indices
+    procedure :: get_local_xyz => lattice_get_local_xyz
+    procedure :: offset => lattice_offset
+    procedure :: outer => lattice_outer
+    procedure :: to_hdf5 => lattice_to_hdf5
   end type Lattice
 
-  abstract interface
-
-!===============================================================================
-! ARE_VALID_INDICES returns .true. if the given lattice indices fit within the
-! bounds of the lattice.  Returns false otherwise.
-
-    function lattice_are_valid_indices_(this, i_xyz) result(is_valid)
-      import Lattice
-      class(Lattice), intent(in) :: this
-      integer,        intent(in) :: i_xyz(3)
-      logical                    :: is_valid
-    end function lattice_are_valid_indices_
-
-!===============================================================================
-! GET_INDICES returns the indices in a lattice for the given global xyz.
-
-    function lattice_get_indices_(this, global_xyz) result(i_xyz)
-      import Lattice
-      class(Lattice), intent(in) :: this
-      real(8),        intent(in) :: global_xyz(3)
-      integer                    :: i_xyz(3)
-    end function lattice_get_indices_
-
-!===============================================================================
-! GET_LOCAL_XYZ returns the translated local version of the given global xyz.
-
-    function lattice_get_local_xyz_(this, global_xyz, i_xyz) result(local_xyz)
-      import Lattice
-      class(Lattice), intent(in) :: this
-      real(8),        intent(in) :: global_xyz(3)
-      integer,        intent(in) :: i_xyz(3)
-      real(8)                    :: local_xyz(3)
-    end function lattice_get_local_xyz_
-  end interface
-
 !===============================================================================
 ! RECTLATTICE extends LATTICE for rectilinear arrays.
 !===============================================================================
 
   type, extends(Lattice) :: RectLattice
-    integer              :: n_cells(3)     ! Number of cells along each axis
-    real(8), allocatable :: lower_left(:)  ! Global lower-left corner of lat
-
-    contains
-
-    procedure :: are_valid_indices => valid_inds_rect
-    procedure :: get_indices => get_inds_rect
-    procedure :: get_local_xyz => get_local_rect
   end type RectLattice
 
 !===============================================================================
@@ -101,15 +237,6 @@ module geometry_header
 !===============================================================================
 
   type, extends(Lattice) :: HexLattice
-    integer              :: n_rings   ! Number of radial ring cell positoins
-    integer              :: n_axial   ! Number of axial cell positions
-    real(8), allocatable :: center(:) ! Global center of lattice
-
-    contains
-
-    procedure :: are_valid_indices => valid_inds_hex
-    procedure :: get_indices => get_inds_hex
-    procedure :: get_local_xyz => get_local_hex
   end type HexLattice
 
 !===============================================================================
@@ -125,25 +252,13 @@ module geometry_header
 !===============================================================================
 
   type Cell
-    integer :: id                          ! Unique ID
-    character(len=104) :: name = ""        ! User-defined name
-    integer :: type                        ! Type of cell (normal, universe,
-                                           !  lattice)
-    integer :: universe                    ! universe # this cell is in
-    integer :: fill                        ! universe # filling this cell
-    integer :: instances                   ! number of instances of this cell in
-                                           !  the geom
+    type(C_PTR) :: ptr
+
     integer, allocatable :: material(:)    ! Material within cell.  Multiple
                                            !  materials for distribcell
                                            !  instances.  0 signifies a universe
-    integer, allocatable :: offset(:)      ! Distribcell offset for tally
-                                           !  counter
     integer, allocatable :: region(:)      ! Definition of spatial region as
                                            !  Boolean expression of half-spaces
-    integer, allocatable :: rpn(:)         ! Reverse Polish notation for region
-                                           !  expression
-    logical :: simple                      ! Is the region simple (intersections
-                                           !  only)
     integer :: distribcell_index           ! Index corresponding to this cell in
                                            !  distribcell arrays
     real(8), allocatable :: sqrtkT(:)      ! Square root of k_Boltzmann *
@@ -154,6 +269,22 @@ module geometry_header
     real(8), allocatable :: translation(:)
     real(8), allocatable :: rotation(:)
     real(8), allocatable :: rotation_matrix(:,:)
+
+  contains
+
+    procedure :: id => cell_id
+    procedure :: set_id => cell_set_id
+    procedure :: type => cell_type
+    procedure :: set_type => cell_set_type
+    procedure :: universe => cell_universe
+    procedure :: set_universe => cell_set_universe
+    procedure :: fill => cell_fill
+    procedure :: n_instances => cell_n_instances
+    procedure :: simple => cell_simple
+    procedure :: distance => cell_distance
+    procedure :: offset => cell_offset
+    procedure :: to_hdf5 => cell_to_hdf5
+
   end type Cell
 
   ! array index of the root universe
@@ -161,7 +292,6 @@ module geometry_header
 
   integer(C_INT32_T), bind(C) :: n_cells     ! # of cells
   integer(C_INT32_T), bind(C) :: n_universes ! # of universes
-  integer(C_INT32_T), bind(C) :: n_lattices  ! # of lattices
 
   type(Cell),             allocatable, target :: cells(:)
   type(Universe),         allocatable, target :: universes(:)
@@ -174,172 +304,150 @@ module geometry_header
 
 contains
 
-!===============================================================================
+  function lattice_id(this) result(id)
+    class(Lattice), intent(in) :: this
+    integer(C_INT32_T)         :: id
+    id = lattice_id_c(this % ptr)
+  end function lattice_id
 
-  function valid_inds_rect(this, i_xyz) result(is_valid)
-    class(RectLattice), intent(in) :: this
-    integer,            intent(in) :: i_xyz(3)
-    logical                        :: is_valid
+  function lattice_are_valid_indices(this, i_xyz) result (is_valid)
+    class(Lattice), intent(in) :: this
+    integer(C_INT), intent(in) :: i_xyz(3)
+    logical(C_BOOL)            :: is_valid
+    is_valid = lattice_are_valid_indices_c(this % ptr, i_xyz)
+  end function lattice_are_valid_indices
 
-    is_valid = all(i_xyz > 0 .and. i_xyz <= this % n_cells)
-  end function valid_inds_rect
+  subroutine lattice_distance(this, xyz, uvw, i_xyz, d, lattice_trans)
+    class(Lattice), intent(in)  :: this
+    real(C_DOUBLE), intent(in)  :: xyz(3)
+    real(C_DOUBLE), intent(in)  :: uvw(3)
+    integer(C_INT), intent(in)  :: i_xyz(3)
+    real(C_DOUBLE), intent(out) :: d
+    integer(C_INT), intent(out) :: lattice_trans(3)
+    call lattice_distance_c(this % ptr, xyz, uvw, i_xyz, d, lattice_trans)
+  end subroutine lattice_distance
+
+  function lattice_get(this, i_xyz) result(univ)
+    class(Lattice), intent(in)  :: this
+    integer(C_INT), intent(in)  :: i_xyz(3)
+    integer(C_INT32_T)          :: univ
+    univ = lattice_universe_c(this % ptr, i_xyz)
+  end function lattice_get
+
+  function lattice_get_indices(this, xyz) result(i_xyz)
+    class(Lattice), intent(in)  :: this
+    real(C_DOUBLE), intent(in)  :: xyz(3)
+    integer(C_INT)              :: i_xyz(3)
+    call lattice_get_indices_c(this % ptr, xyz, i_xyz)
+  end function lattice_get_indices
+
+  function lattice_get_local_xyz(this, global_xyz, i_xyz) &
+       result(local_xyz)
+    class(Lattice), intent(in) :: this
+    real(C_DOUBLE), intent(in) :: global_xyz(3)
+    integer(C_INT), intent(in) :: i_xyz(3)
+    real(C_DOUBLE)             :: local_xyz(3)
+    call lattice_get_local_xyz_c(this % ptr, global_xyz, i_xyz, local_xyz)
+  end function lattice_get_local_xyz
+
+  function lattice_offset(this, map, i_xyz) result(offset)
+    class(Lattice), intent(in) :: this
+    integer(C_INT), intent(in) :: map
+    integer(C_INT), intent(in) :: i_xyz(3)
+    integer(C_INT32_T)         :: offset
+    offset = lattice_offset_c(this % ptr, map, i_xyz)
+  end function lattice_offset
+
+  function lattice_outer(this) result(outer)
+    class(Lattice), intent(in) :: this
+    integer(C_INT32_T)         :: outer
+    outer = lattice_outer_c(this % ptr)
+  end function lattice_outer
+
+  subroutine lattice_to_hdf5(this, group)
+    class(Lattice), intent(in) :: this
+    integer(HID_T), intent(in) :: group
+    call lattice_to_hdf5_c(this % ptr, group)
+  end subroutine lattice_to_hdf5
 
 !===============================================================================
 
-  function valid_inds_hex(this, i_xyz) result(is_valid)
-    class(HexLattice), intent(in) :: this
-    integer,           intent(in) :: i_xyz(3)
-    logical                       :: is_valid
+  function cell_id(this) result(id)
+    class(Cell), intent(in) :: this
+    integer(C_INT32_T)      :: id
+    id = cell_id_c(this % ptr)
+  end function cell_id
 
-    is_valid = (all(i_xyz > 0) .and. &
-               &i_xyz(1) < 2*this % n_rings .and. &
-               &i_xyz(2) < 2*this % n_rings .and. &
-               &i_xyz(1) + i_xyz(2) > this % n_rings .and. &
-               &i_xyz(1) + i_xyz(2) < 3*this % n_rings .and. &
-               &i_xyz(3) <= this % n_axial)
-  end function valid_inds_hex
+  subroutine cell_set_id(this, id)
+    class(Cell),        intent(in) :: this
+    integer(C_INT32_T), intent(in) :: id
+    call cell_set_id_c(this % ptr, id)
+  end subroutine cell_set_id
 
-!===============================================================================
+  function cell_type(this) result(type)
+    class(Cell), intent(in) :: this
+    integer(C_INT)          :: type
+    type = cell_type_c(this % ptr)
+  end function cell_type
 
-  function get_inds_rect(this, global_xyz) result(i_xyz)
-    class(RectLattice), intent(in) :: this
-    real(8),            intent(in) :: global_xyz(3)
-    integer                        :: i_xyz(3)
+  subroutine cell_set_type(this, type)
+    class(Cell),    intent(in) :: this
+    integer(C_INT), intent(in) :: type
+    call cell_set_type_c(this % ptr, type)
+  end subroutine cell_set_type
 
-    real(8) :: xyz(3)  ! global_xyz alias
+  function cell_universe(this) result(universe)
+    class(Cell), intent(in) :: this
+    integer(C_INT32_T)      :: universe
+    universe = cell_universe_c(this % ptr)
+  end function cell_universe
 
-    xyz = global_xyz
+  subroutine cell_set_universe(this, universe)
+    class(Cell),        intent(in) :: this
+    integer(C_INT32_T), intent(in) :: universe
+    call cell_set_universe_c(this % ptr, universe)
+  end subroutine cell_set_universe
 
-    i_xyz(1) = ceiling((xyz(1) - this % lower_left(1))/this % pitch(1))
-    i_xyz(2) = ceiling((xyz(2) - this % lower_left(2))/this % pitch(2))
-    if (this % is_3d) then
-      i_xyz(3) = ceiling((xyz(3) - this % lower_left(3))/this % pitch(3))
-    else
-      i_xyz(3) = 1
-    end if
-  end function get_inds_rect
+  function cell_fill(this) result(fill)
+    class(Cell), intent(in) :: this
+    integer(C_INT32_T)      :: fill
+    fill = cell_fill_c(this % ptr)
+  end function cell_fill
 
-!===============================================================================
+  function cell_n_instances(this) result(n_instances)
+    class(Cell), intent(in) :: this
+    integer(C_INT32_T)      :: n_instances
+    n_instances = cell_n_instances_c(this % ptr)
+  end function cell_n_instances
 
-  function get_inds_hex(this, global_xyz) result(i_xyz)
-    class(HexLattice), intent(in) :: this
-    real(8),           intent(in) :: global_xyz(3)
-    integer                       :: i_xyz(3)
+  function cell_simple(this) result(simple)
+    class(Cell), intent(in) :: this
+    logical(C_BOOL)         :: simple
+    simple = cell_simple_c(this % ptr)
+  end function cell_simple
 
-    real(8) :: xyz(3)    ! global xyz relative to the center
-    real(8) :: alpha     ! Skewed coord axis
-    real(8) :: xyz_t(3)  ! Local xyz
-    real(8) :: d, d_min  ! Squared distance from cell centers
-    integer :: i, j, k   ! Iterators
-    integer :: k_min     ! Minimum distance index
+  subroutine cell_distance(this, xyz, uvw, on_surface, min_dist, i_surf)
+    class(Cell),        intent(in)  :: this
+    real(C_DOUBLE),     intent(in)  :: xyz(3)
+    real(C_DOUBLE),     intent(in)  :: uvw(3)
+    integer(C_INT32_T), intent(in)  :: on_surface
+    real(C_DOUBLE),     intent(out) :: min_dist
+    integer(C_INT32_T), intent(out) :: i_surf
+    call cell_distance_c(this % ptr, xyz, uvw, on_surface, min_dist, i_surf)
+  end subroutine cell_distance
 
-    xyz(1) = global_xyz(1) - this % center(1)
-    xyz(2) = global_xyz(2) - this % center(2)
+  function cell_offset(this, map) result(offset)
+    class(Cell), intent(in)    :: this
+    integer(C_INT), intent(in) :: map
+    integer(C_INT32_T)         :: offset
+    offset = cell_offset_c(this % ptr, map)
+  end function cell_offset
 
-    ! Index z direction.
-    if (this % is_3d) then
-      xyz(3) = global_xyz(3) - this % center(3)
-      i_xyz(3) = ceiling(xyz(3)/this % pitch(2) + HALF*this % n_axial)
-    else
-      xyz(3) = global_xyz(3)
-      i_xyz(3) = 1
-    end if
-
-    ! Convert coordinates into skewed bases.  The (x, alpha) basis is used to
-    ! find the index of the global coordinates to within 4 cells.
-    alpha = xyz(2) - xyz(1) / sqrt(THREE)
-    i_xyz(1) = floor(xyz(1) / (sqrt(THREE) / TWO * this % pitch(1)))
-    i_xyz(2) = floor(alpha / this % pitch(1))
-
-    ! Add offset to indices (the center cell is (i_x, i_alpha) = (0, 0) but
-    ! the array is offset so that the indices never go below 1).
-    i_xyz(1) = i_xyz(1) + this % n_rings
-    i_xyz(2) = i_xyz(2) + this % n_rings
-
-    ! Calculate the (squared) distance between the particle and the centers of
-    ! the four possible cells.  Regular hexagonal tiles form a centroidal
-    ! Voronoi tessellation so the global xyz should be in the hexagonal cell
-    ! that it is closest to the center of.  This method is used over a
-    ! method that uses the remainders of the floor divisions above because it
-    ! provides better finite precision performance.  Squared distances are
-    ! used becasue they are more computationally efficient than normal
-    ! distances.
-    k = 1
-    d_min = INFINITY
-    do i = 0, 1
-      do j = 0, 1
-        xyz_t = this % get_local_xyz(global_xyz, i_xyz + [j, i, 0])
-        d = xyz_t(1)**2 + xyz_t(2)**2
-        if (d < d_min) then
-          d_min = d
-          k_min = k
-        end if
-        k = k + 1
-      end do
-    end do
-
-    ! Select the minimum squared distance which corresponds to the cell the
-    ! coordinates are in.
-    if (k_min == 2) then
-      i_xyz(1) = i_xyz(1) + 1
-    else if (k_min == 3) then
-      i_xyz(2) = i_xyz(2) + 1
-    else if (k_min == 4) then
-      i_xyz(1) = i_xyz(1) + 1
-      i_xyz(2) = i_xyz(2) + 1
-    end if
-  end function get_inds_hex
-
-!===============================================================================
-
-  function get_local_rect(this, global_xyz, i_xyz) result(local_xyz)
-    class(RectLattice), intent(in) :: this
-    real(8),            intent(in) :: global_xyz(3)
-    integer,            intent(in) :: i_xyz(3)
-    real(8)                        :: local_xyz(3)
-
-    real(8) :: xyz(3)  ! global_xyz alias
-
-    xyz = global_xyz
-
-    local_xyz(1) = xyz(1) - (this % lower_left(1) + &
-         (i_xyz(1) - HALF)*this % pitch(1))
-    local_xyz(2) = xyz(2) - (this % lower_left(2) + &
-         (i_xyz(2) - HALF)*this % pitch(2))
-    if (this % is_3d) then
-      local_xyz(3) = xyz(3) - (this % lower_left(3) + &
-           (i_xyz(3) - HALF)*this % pitch(3))
-    else
-      local_xyz(3) = xyz(3)
-    end if
-  end function get_local_rect
-
-!===============================================================================
-
-  function get_local_hex(this, global_xyz, i_xyz) result(local_xyz)
-    class(HexLattice), intent(in) :: this
-    real(8),           intent(in) :: global_xyz(3)
-    integer,           intent(in) :: i_xyz(3)
-    real(8)                       :: local_xyz(3)
-
-    real(8) :: xyz(3)  ! global_xyz alias
-
-    xyz = global_xyz
-
-    ! x_l = x_g - (center + pitch_x*cos(30)*index_x)
-    local_xyz(1) = xyz(1) - (this % center(1) + &
-         sqrt(THREE) / TWO * (i_xyz(1) - this % n_rings) * this % pitch(1))
-    ! y_l = y_g - (center + pitch_x*index_x + pitch_y*sin(30)*index_y)
-    local_xyz(2) = xyz(2) - (this % center(2) + &
-         (i_xyz(2) - this % n_rings) * this % pitch(1) + &
-         (i_xyz(1) - this % n_rings) * this % pitch(1) / TWO)
-    if (this % is_3d) then
-      local_xyz(3) = xyz(3) - this % center(3) &
-           + (HALF*this % n_axial - i_xyz(3) + HALF) * this % pitch(2)
-    else
-      local_xyz(3) = xyz(3)
-    end if
-  end function get_local_hex
+  subroutine cell_to_hdf5(this, group)
+    class(Cell),    intent(in) :: this
+    integer(HID_T), intent(in) :: group
+    call cell_to_hdf5_c(this % ptr, group)
+  end subroutine cell_to_hdf5
 
 !===============================================================================
 ! GET_TEMPERATURES returns a list of temperatures that each nuclide/S(a,b) table
@@ -408,10 +516,15 @@ contains
 !===============================================================================
 
   subroutine free_memory_geometry()
+    interface
+      subroutine free_memory_geometry_c() bind(C)
+      end subroutine free_memory_geometry_c
+    end interface
+
+    call free_memory_geometry_c()
 
     n_cells = 0
     n_universes = 0
-    n_lattices = 0
 
     if (allocated(cells)) deallocate(cells)
     if (allocated(universes)) deallocate(universes)
@@ -432,6 +545,7 @@ contains
     integer(C_INT32_T), value, intent(in) :: n
     integer(C_INT32_T), optional, intent(out) :: index_start
     integer(C_INT32_T), optional, intent(out) :: index_end
+    integer(C_INT32_T) :: i
     integer(C_INT) :: err
 
     type(Cell), allocatable :: temp(:) ! temporary cells array
@@ -453,7 +567,12 @@ contains
     ! Return indices in cells array
     if (present(index_start)) index_start = n_cells + 1
     if (present(index_end)) index_end = n_cells + n
-    n_cells = n_cells + n
+
+    ! Extend the C++ cells array and get pointers to the C++ objects
+    call extend_cells_c(n)
+    do i = n_cells - n, n_cells
+      cells(i) % ptr = cell_pointer_c(i - 1)
+    end do
 
     err = 0
   end function openmc_extend_cells
@@ -490,14 +609,14 @@ contains
     err = 0
     if (index >= 1 .and. index <= size(cells)) then
       associate (c => cells(index))
-        type = c % type
+        type = c % type()
         select case (type)
         case (FILL_MATERIAL)
           n = size(c % material)
           indices = C_LOC(c % material(1))
         case (FILL_UNIVERSE, FILL_LATTICE)
           n = 1
-          indices = C_LOC(c % fill)
+          indices = cell_fill_ptr(c % ptr)
         end select
       end associate
     else
@@ -514,7 +633,7 @@ contains
     integer(C_INT) :: err
 
     if (index >= 1 .and. index <= size(cells)) then
-      id = cells(index) % id
+      id = cells(index) % id()
       err = 0
     else
       err = E_OUT_OF_BOUNDS
@@ -541,7 +660,7 @@ contains
           if (allocated(c % material)) deallocate(c % material)
           allocate(c % material(n))
 
-          c % type = FILL_MATERIAL
+          call c % set_type(FILL_MATERIAL)
           do i = 1, n
             j = indices(i)
             if ((j >= 1 .and. j <= n_materials) .or. j == MATERIAL_VOID) then
@@ -553,9 +672,9 @@ contains
             end if
           end do
         case (FILL_UNIVERSE)
-          c % type = FILL_UNIVERSE
+          call c % set_type(FILL_UNIVERSE)
         case (FILL_LATTICE)
-          c % type = FILL_LATTICE
+          call c % set_type(FILL_LATTICE)
         end select
       end associate
     else
@@ -573,7 +692,7 @@ contains
     integer(C_INT) :: err
 
     if (index >= 1 .and. index <= n_cells) then
-      cells(index) % id = id
+      call cells(index) % set_id(id)
       call cell_dict % set(id, index)
       err = 0
     else
@@ -609,7 +728,7 @@ contains
     if (index >= 1 .and. index <= size(cells)) then
 
       ! error if the cell is filled with another universe
-      if (cells(index) % fill /= NONE) then
+      if (cells(index) % fill() /= C_NONE) then
         err = E_GEOMETRY
         call set_errmsg("Cannot set temperature on a cell filled &
              &with a universe.")

src/hdf5_interface.h

@@ -85,12 +85,21 @@ extern "C" void write_tally_results(hid_t group_id, hsize_t n_filter, hsize_t n_
                                     const double* results);
 
 template<std::size_t array_len> void
-write_double_1D(hid_t group_id, char const *name,
-                std::array<double, array_len> &buffer)
+write_int(hid_t group_id, char const *name,
+          const std::array<int, array_len> &buffer, bool indep)
+{
+  hsize_t dims[1] {array_len};
+  write_dataset(group_id, 1, dims, name, H5T_NATIVE_INT, buffer.data(), indep);
+}
+
+
+template<std::size_t array_len> void
+write_double(hid_t group_id, char const *name,
+             const std::array<double, array_len> &buffer, bool indep)
 {
   hsize_t dims[1] {array_len};
   write_dataset(group_id, 1, dims, name, H5T_NATIVE_DOUBLE,
-                buffer.data(), false);
+                buffer.data(), indep);
 }
 
 } // namespace openmc

src/lattice.cpp

@@ -0,0 +1,946 @@
+#include "lattice.h"
+
+#include <cmath>
+#include <sstream>
+#include <vector>
+
+#include "cell.h"
+#include "error.h"
+#include "geometry_aux.h"
+#include "hdf5_interface.h"
+#include "string_utils.h"
+#include "xml_interface.h"
+
+
+namespace openmc {
+
+//==============================================================================
+// Global variables
+//==============================================================================
+
+std::vector<Lattice*> lattices_c;
+
+std::unordered_map<int32_t, int32_t> lattice_map;
+
+//==============================================================================
+// Lattice implementation
+//==============================================================================
+
+Lattice::Lattice(pugi::xml_node lat_node)
+{
+  if (check_for_node(lat_node, "id")) {
+    id = stoi(get_node_value(lat_node, "id"));
+  } else {
+    fatal_error("Must specify id of lattice in geometry XML file.");
+  }
+
+  if (check_for_node(lat_node, "name")) {
+    name = get_node_value(lat_node, "name");
+  }
+
+  if (check_for_node(lat_node, "outer")) {
+    outer = stoi(get_node_value(lat_node, "outer"));
+  }
+}
+
+//==============================================================================
+
+LatticeIter Lattice::begin()
+{return LatticeIter(*this, 0);}
+
+LatticeIter Lattice::end()
+{return LatticeIter(*this, universes.size());}
+
+ReverseLatticeIter Lattice::rbegin()
+{return ReverseLatticeIter(*this, universes.size()-1);}
+
+ReverseLatticeIter Lattice::rend()
+{return ReverseLatticeIter(*this, -1);}
+
+//==============================================================================
+
+void
+Lattice::adjust_indices()
+{
+  // Adjust the indices for the universes array.
+  for (LatticeIter it = begin(); it != end(); ++it) {
+    int uid = *it;
+    auto search = universe_map.find(uid);
+    if (search != universe_map.end()) {
+      *it = search->second;
+    } else {
+      std::stringstream err_msg;
+      err_msg << "Invalid universe number " << uid << " specified on "
+           "lattice " << id;
+      fatal_error(err_msg);
+    }
+  }
+
+  // Adjust the index for the outer universe.
+  if (outer != NO_OUTER_UNIVERSE) {
+    auto search = universe_map.find(outer);
+    if (search != universe_map.end()) {
+      outer = search->second;
+    } else {
+      std::stringstream err_msg;
+      err_msg << "Invalid universe number " << outer << " specified on "
+           "lattice " << id;
+      fatal_error(err_msg);
+    }
+  }
+}
+
+//==============================================================================
+
+int32_t
+Lattice::fill_offset_table(int32_t offset, int32_t target_univ_id, int map)
+{
+  for (LatticeIter it = begin(); it != end(); ++it) {
+    offsets[map * universes.size() + it.indx] = offset;
+    offset += count_universe_instances(*it, target_univ_id);
+  }
+  return offset;
+}
+
+//==============================================================================
+
+void
+Lattice::to_hdf5(hid_t lattices_group) const
+{
+  // Make a group for the lattice.
+  std::string group_name {"lattice "};
+  group_name += std::to_string(id);
+  hid_t lat_group = create_group(lattices_group, group_name);
+
+  // Write the name and outer universe.
+  if (!name.empty()) {
+    write_string(lat_group, "name", name, false);
+  }
+
+  if (outer != NO_OUTER_UNIVERSE) {
+    int32_t outer_id = global_universes[outer]->id;
+    write_int(lat_group, 0, nullptr, "outer", &outer_id, false);
+  } else {
+    write_int(lat_group, 0, nullptr, "outer", &outer, false);
+  }
+
+  // Call subclass-overriden function to fill in other details.
+  to_hdf5_inner(lat_group);
+
+  close_group(lat_group);
+}
+
+//==============================================================================
+// RectLattice implementation
+//==============================================================================
+
+RectLattice::RectLattice(pugi::xml_node lat_node)
+  : Lattice {lat_node}
+{
+  // Read the number of lattice cells in each dimension.
+  std::string dimension_str {get_node_value(lat_node, "dimension")};
+  std::vector<std::string> dimension_words {split(dimension_str)};
+  if (dimension_words.size() == 2) {
+    n_cells[0] = stoi(dimension_words[0]);
+    n_cells[1] = stoi(dimension_words[1]);
+    n_cells[2] = 1;
+    is_3d = false;
+  } else if (dimension_words.size() == 3) {
+    n_cells[0] = stoi(dimension_words[0]);
+    n_cells[1] = stoi(dimension_words[1]);
+    n_cells[2] = stoi(dimension_words[2]);
+    is_3d = true;
+  } else {
+    fatal_error("Rectangular lattice must be two or three dimensions.");
+  }
+
+  // Read the lattice lower-left location.
+  std::string ll_str {get_node_value(lat_node, "lower_left")};
+  std::vector<std::string> ll_words {split(ll_str)};
+  if (ll_words.size() != dimension_words.size()) {
+    fatal_error("Number of entries on <lower_left> must be the same as the "
+                "number of entries on <dimension>.");
+  }
+  lower_left[0] = stod(ll_words[0]);
+  lower_left[1] = stod(ll_words[1]);
+  if (is_3d) {lower_left[2] = stod(ll_words[2]);}
+
+  // Read the lattice pitches.
+  std::string pitch_str {get_node_value(lat_node, "pitch")};
+  std::vector<std::string> pitch_words {split(pitch_str)};
+  if (pitch_words.size() != dimension_words.size()) {
+    fatal_error("Number of entries on <pitch> must be the same as the "
+                "number of entries on <dimension>.");
+  }
+  pitch[0] = stod(pitch_words[0]);
+  pitch[1] = stod(pitch_words[1]);
+  if (is_3d) {pitch[2] = stod(pitch_words[2]);}
+
+  // Read the universes and make sure the correct number was specified.
+  std::string univ_str {get_node_value(lat_node, "universes")};
+  std::vector<std::string> univ_words {split(univ_str)};
+  if (univ_words.size() != nx*ny*nz) {
+    std::stringstream err_msg;
+    err_msg << "Expected " << nx*ny*nz
+            << " universes for a rectangular lattice of size "
+            << nx << "x" << ny << "x" << nz << " but " << univ_words.size()
+            << " were specified.";
+    fatal_error(err_msg);
+  }
+
+  // Parse the universes.
+  universes.resize(nx*ny*nz, C_NONE);
+  for (int iz = 0; iz < nz; iz++) {
+    for (int iy = ny-1; iy > -1; iy--) {
+      for (int ix = 0; ix < nx; ix++) {
+        int indx1 = nx*ny*iz + nx*(ny-iy-1) + ix;
+        int indx2 = nx*ny*iz + nx*iy + ix;
+        universes[indx1] = stoi(univ_words[indx2]);
+      }
+    }
+  }
+}
+
+//==============================================================================
+
+int32_t&
+RectLattice::operator[](const int i_xyz[3])
+{
+  int indx = nx*ny*i_xyz[2] + nx*i_xyz[1] + i_xyz[0];
+  return universes[indx];
+}
+
+//==============================================================================
+
+bool
+RectLattice::are_valid_indices(const int i_xyz[3]) const
+{
+  return (   (i_xyz[0] >= 0) && (i_xyz[0] < n_cells[0])
+          && (i_xyz[1] >= 0) && (i_xyz[1] < n_cells[1])
+          && (i_xyz[2] >= 0) && (i_xyz[2] < n_cells[2]));
+}
+
+//==============================================================================
+
+std::pair<double, std::array<int, 3>>
+RectLattice::distance(const double xyz[3], const double uvw[3],
+                      const int i_xyz[3]) const
+{
+  // Get short aliases to the coordinates.
+  double x {xyz[0]};
+  double y {xyz[1]};
+  double z {xyz[2]};
+  double u {uvw[0]};
+  double v {uvw[1]};
+
+  // Determine the oncoming edge.
+  double x0 {copysign(0.5 * pitch[0], u)};
+  double y0 {copysign(0.5 * pitch[1], v)};
+
+  // Left and right sides
+  double d {INFTY};
+  std::array<int, 3> lattice_trans;
+  if ((std::abs(x - x0) > FP_PRECISION) && u != 0) {
+    d = (x0 - x) / u;
+    if (u > 0) {
+      lattice_trans = {1, 0, 0};
+    } else {
+      lattice_trans = {-1, 0, 0};
+    }
+  }
+
+  // Front and back sides
+  if ((std::abs(y - y0) > FP_PRECISION) && v != 0) {
+    double this_d = (y0 - y) / v;
+    if (this_d < d) {
+      d = this_d;
+      if (v > 0) {
+        lattice_trans = {0, 1, 0};
+      } else {
+        lattice_trans = {0, -1, 0};
+      }
+    }
+  }
+
+  // Top and bottom sides
+  if (is_3d) {
+    double w {uvw[2]};
+    double z0 {copysign(0.5 * pitch[2], w)};
+    if ((std::abs(z - z0) > FP_PRECISION) && w != 0) {
+      double this_d = (z0 - z) / w;
+      if (this_d < d) {
+        d = this_d;
+        if (w > 0) {
+          lattice_trans = {0, 0, 1};
+        } else {
+          lattice_trans = {0, 0, -1};
+        }
+      }
+    }
+  }
+
+  return {d, lattice_trans};
+}
+
+//==============================================================================
+
+std::array<int, 3>
+RectLattice::get_indices(const double xyz[3]) const
+{
+  int ix {static_cast<int>(std::ceil((xyz[0] - lower_left[0]) / pitch[0]))-1};
+  int iy {static_cast<int>(std::ceil((xyz[1] - lower_left[1]) / pitch[1]))-1};
+  int iz;
+  if (is_3d) {
+    iz = static_cast<int>(std::ceil((xyz[2] - lower_left[2]) / pitch[2]))-1;
+  } else {
+    iz = 0;
+  }
+  return {ix, iy, iz};
+}
+
+//==============================================================================
+
+std::array<double, 3>
+RectLattice::get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const
+{
+  std::array<double, 3> local_xyz;
+  local_xyz[0] = global_xyz[0] - (lower_left[0] + (i_xyz[0] + 0.5)*pitch[0]);
+  local_xyz[1] = global_xyz[1] - (lower_left[1] + (i_xyz[1] + 0.5)*pitch[1]);
+  if (is_3d) {
+    local_xyz[2] = global_xyz[2] - (lower_left[2] + (i_xyz[2] + 0.5)*pitch[2]);
+  } else {
+    local_xyz[2] = global_xyz[2];
+  }
+  return local_xyz;
+}
+
+//==============================================================================
+
+int32_t&
+RectLattice::offset(int map, const int i_xyz[3])
+{
+  return offsets[nx*ny*nz*map + nx*ny*i_xyz[2] + nx*i_xyz[1] + i_xyz[0]];
+}
+
+//==============================================================================
+
+std::string
+RectLattice::index_to_string(int indx) const
+{
+  int iz {indx / (nx * ny)};
+  int iy {(indx - nx * ny * iz) / nx};
+  int ix {indx - nx * ny * iz - nx * iy};
+  std::string out {std::to_string(ix)};
+  out += ',';
+  out += std::to_string(iy);
+  if (is_3d) {
+    out += ',';
+    out += std::to_string(iz);
+  }
+  return out;
+}
+
+//==============================================================================
+
+void
+RectLattice::to_hdf5_inner(hid_t lat_group) const
+{
+  // Write basic lattice information.
+  write_string(lat_group, "type", "rectangular", false);
+  if (is_3d) {
+    write_double(lat_group, "pitch", pitch, false);
+    write_double(lat_group, "lower_left", lower_left, false);
+    write_int(lat_group, "dimension", n_cells, false);
+  } else {
+    std::array<double, 2> pitch_short {{pitch[0], pitch[1]}};
+    write_double(lat_group, "pitch", pitch_short, false);
+    std::array<double, 2> ll_short {{lower_left[0], lower_left[1]}};
+    write_double(lat_group, "lower_left", ll_short, false);
+    std::array<int, 2> nc_short {{n_cells[0], n_cells[1]}};
+    write_int(lat_group, "dimension", nc_short, false);
+  }
+
+  // Write the universe ids.  The convention here is to switch the ordering on
+  // the y-axis to match the way universes are input in a text file.
+  if (is_3d) {
+    hsize_t nx {static_cast<hsize_t>(n_cells[0])};
+    hsize_t ny {static_cast<hsize_t>(n_cells[1])};
+    hsize_t nz {static_cast<hsize_t>(n_cells[2])};
+    int out[nx*ny*nz];
+
+    for (int m = 0; m < nz; m++) {
+      for (int k = 0; k < ny; k++) {
+        for (int j = 0; j < nx; j++) {
+          int indx1 = nx*ny*m + nx*k + j;
+          int indx2 = nx*ny*m + nx*(ny-k-1) + j;
+          out[indx2] = global_universes[universes[indx1]]->id;
+        }
+      }
+    }
+
+    hsize_t dims[3] {nz, ny, nx};
+    write_int(lat_group, 3, dims, "universes", out, false);
+
+  } else {
+    hsize_t nx {static_cast<hsize_t>(n_cells[0])};
+    hsize_t ny {static_cast<hsize_t>(n_cells[1])};
+    int out[nx*ny];
+
+    for (int k = 0; k < ny; k++) {
+      for (int j = 0; j < nx; j++) {
+        int indx1 = nx*k + j;
+        int indx2 = nx*(ny-k-1) + j;
+        out[indx2] = global_universes[universes[indx1]]->id;
+      }
+    }
+
+    hsize_t dims[3] {1, ny, nx};
+    write_int(lat_group, 3, dims, "universes", out, false);
+  }
+}
+
+//==============================================================================
+// HexLattice implementation
+//==============================================================================
+
+HexLattice::HexLattice(pugi::xml_node lat_node)
+  : Lattice {lat_node}
+{
+  // Read the number of lattice cells in each dimension.
+  n_rings = stoi(get_node_value(lat_node, "n_rings"));
+  if (check_for_node(lat_node, "n_axial")) {
+    n_axial = stoi(get_node_value(lat_node, "n_axial"));
+    is_3d = true;
+  } else {
+    n_axial = 1;
+    is_3d = false;
+  }
+
+  // Read the lattice center.
+  std::string center_str {get_node_value(lat_node, "center")};
+  std::vector<std::string> center_words {split(center_str)};
+  if (is_3d && (center_words.size() != 3)) {
+    fatal_error("A hexagonal lattice with <n_axial> must have <center> "
+                "specified by 3 numbers.");
+  } else if (!is_3d && center_words.size() != 2) {
+    fatal_error("A hexagonal lattice without <n_axial> must have <center> "
+                "specified by 2 numbers.");
+  }
+  center[0] = stod(center_words[0]);
+  center[1] = stod(center_words[1]);
+  if (is_3d) {center[2] = stod(center_words[2]);}
+
+  // Read the lattice pitches.
+  std::string pitch_str {get_node_value(lat_node, "pitch")};
+  std::vector<std::string> pitch_words {split(pitch_str)};
+  if (is_3d && (pitch_words.size() != 2)) {
+    fatal_error("A hexagonal lattice with <n_axial> must have <pitch> "
+                "specified by 2 numbers.");
+  } else if (!is_3d && (pitch_words.size() != 1)) {
+    fatal_error("A hexagonal lattice without <n_axial> must have <pitch> "
+                "specified by 1 number.");
+  }
+  pitch[0] = stod(pitch_words[0]);
+  if (is_3d) {pitch[1] = stod(pitch_words[1]);}
+
+  // Read the universes and make sure the correct number was specified.
+  int n_univ = (3*n_rings*n_rings - 3*n_rings + 1) * n_axial;
+  std::string univ_str {get_node_value(lat_node, "universes")};
+  std::vector<std::string> univ_words {split(univ_str)};
+  if (univ_words.size() != n_univ) {
+    std::stringstream err_msg;
+    err_msg << "Expected " << n_univ
+            << " universes for a hexagonal lattice with " << n_rings
+            << " rings and " << n_axial << " axial levels" << " but "
+            << univ_words.size() << " were specified.";
+    fatal_error(err_msg);
+  }
+
+  // Parse the universes.
+  // Universes in hexagonal lattices are stored in a manner that represents
+  // a skewed coordinate system: (x, alpha) rather than (x, y).  There is
+  // no obvious, direct relationship between the order of universes in the
+  // input and the order that they will be stored in the skewed array so
+  // the following code walks a set of index values across the skewed array
+  // in a manner that matches the input order.  Note that i_x = 0, i_a = 0
+  // corresponds to the center of the hexagonal lattice.
+
+  universes.resize((2*n_rings-1) * (2*n_rings-1) * n_axial, C_NONE);
+  int input_index = 0;
+  for (int m = 0; m < n_axial; m++) {
+    // Initialize lattice indecies.
+    int i_x = 1;
+    int i_a = n_rings - 1;
+
+    // Map upper triangular region of hexagonal lattice which is found in the
+    // first n_rings-1 rows of the input.
+    for (int k = 0; k < n_rings-1; k++) {
+      // Walk the index to lower-left neighbor of last row start.
+      i_x -= 1;
+
+      // Iterate over the input columns.
+      for (int j = 0; j < k+1; j++) {
+        int indx = (2*n_rings-1)*(2*n_rings-1) * m
+                    + (2*n_rings-1) * (i_a+n_rings-1)
+                    + (i_x+n_rings-1);
+        universes[indx] = stoi(univ_words[input_index]);
+        input_index++;
+        // Walk the index to the right neighbor (which is not adjacent).
+        i_x += 2;
+        i_a -= 1;
+      }
+
+      // Return the lattice index to the start of the current row.
+      i_x -= 2 * (k+1);
+      i_a += (k+1);
+    }
+
+    // Map the middle square region of the hexagonal lattice which is found in
+    // the next 2*n_rings-1 rows of the input.
+    for (int k = 0; k < 2*n_rings-1; k++) {
+      if ((k % 2) == 0) {
+        // Walk the index to the lower-left neighbor of the last row start.
+        i_x -= 1;
+      } else {
+        // Walk the index to the lower-right neighbor of the last row start.
+        i_x += 1;
+        i_a -= 1;
+      }
+
+      // Iterate over the input columns.
+      for (int j = 0; j < n_rings - (k % 2); j++) {
+        int indx = (2*n_rings-1)*(2*n_rings-1) * m
+                    + (2*n_rings-1) * (i_a+n_rings-1)
+                    + (i_x+n_rings-1);
+        universes[indx] = stoi(univ_words[input_index]);
+        input_index++;
+        // Walk the index to the right neighbor (which is not adjacent).
+        i_x += 2;
+        i_a -= 1;
+      }
+
+      // Return the lattice index to the start of the current row.
+      i_x -= 2*(n_rings - (k % 2));
+      i_a += n_rings - (k % 2);
+    }
+
+    // Map the lower triangular region of the hexagonal lattice.
+    for (int k = 0; k < n_rings-1; k++) {
+      // Walk the index to the lower-right neighbor of the last row start.
+      i_x += 1;
+      i_a -= 1;
+
+      // Iterate over the input columns.
+      for (int j = 0; j < n_rings-k-1; j++) {
+        int indx = (2*n_rings-1)*(2*n_rings-1) * m
+                    + (2*n_rings-1) * (i_a+n_rings-1)
+                    + (i_x+n_rings-1);
+        universes[indx] = stoi(univ_words[input_index]);
+        input_index++;
+        // Walk the index to the right neighbor (which is not adjacent).
+        i_x += 2;
+        i_a -= 1;
+      }
+
+      // Return lattice index to start of current row.
+      i_x -= 2*(n_rings - k - 1);
+      i_a += n_rings - k - 1;
+    }
+  }
+}
+
+//==============================================================================
+
+int32_t&
+HexLattice::operator[](const int i_xyz[3])
+{
+  int indx = (2*n_rings-1)*(2*n_rings-1) * i_xyz[2]
+              + (2*n_rings-1) * i_xyz[1]
+              + i_xyz[0];
+  return universes[indx];
+}
+
+//==============================================================================
+
+LatticeIter HexLattice::begin()
+{return LatticeIter(*this, n_rings-1);}
+
+ReverseLatticeIter HexLattice::rbegin()
+{return ReverseLatticeIter(*this, universes.size()-n_rings);}
+
+//==============================================================================
+
+bool
+HexLattice::are_valid_indices(const int i_xyz[3]) const
+{
+  return ((i_xyz[0] >= 0) && (i_xyz[1] >= 0) && (i_xyz[2] >= 0)
+          && (i_xyz[0] < 2*n_rings-1) && (i_xyz[1] < 2*n_rings-1)
+          && (i_xyz[0] + i_xyz[1] > n_rings-2)
+          && (i_xyz[0] + i_xyz[1] < 3*n_rings-2)
+          && (i_xyz[2] < n_axial));
+}
+
+//==============================================================================
+
+std::pair<double, std::array<int, 3>>
+HexLattice::distance(const double xyz[3], const double uvw[3],
+                     const int i_xyz[3]) const
+{
+  // Compute the direction on the hexagonal basis.
+  double beta_dir = uvw[0] * std::sqrt(3.0) / 2.0 + uvw[1] / 2.0;
+  double gamma_dir = uvw[0] * std::sqrt(3.0) / 2.0 - uvw[1] / 2.0;
+
+  // Note that hexagonal lattice distance calculations are performed
+  // using the particle's coordinates relative to the neighbor lattice
+  // cells, not relative to the particle's current cell.  This is done
+  // because there is significant disagreement between neighboring cells
+  // on where the lattice boundary is due to finite precision issues.
+
+  // Upper-right and lower-left sides.
+  double d {INFTY};
+  std::array<int, 3> lattice_trans;
+  double edge = -copysign(0.5*pitch[0], beta_dir);  // Oncoming edge
+  std::array<double, 3> xyz_t;
+  if (beta_dir > 0) {
+    const int i_xyz_t[3] {i_xyz[0]+1, i_xyz[1], i_xyz[2]};
+    xyz_t = get_local_xyz(xyz, i_xyz_t);
+  } else {
+    const int i_xyz_t[3] {i_xyz[0]-1, i_xyz[1], i_xyz[2]};
+    xyz_t = get_local_xyz(xyz, i_xyz_t);
+  }
+  double beta = xyz_t[0] * std::sqrt(3.0) / 2.0 + xyz_t[1] / 2.0;
+  if ((std::abs(beta - edge) > FP_PRECISION) && beta_dir != 0) {
+    d = (edge - beta) / beta_dir;
+    if (beta_dir > 0) {
+      lattice_trans = {1, 0, 0};
+    } else {
+      lattice_trans = {-1, 0, 0};
+    }
+  }
+
+  // Lower-right and upper-left sides.
+  edge = -copysign(0.5*pitch[0], gamma_dir);
+  if (gamma_dir > 0) {
+    const int i_xyz_t[3] {i_xyz[0]+1, i_xyz[1]-1, i_xyz[2]};
+    xyz_t = get_local_xyz(xyz, i_xyz_t);
+  } else {
+    const int i_xyz_t[3] {i_xyz[0]-1, i_xyz[1]+1, i_xyz[2]};
+    xyz_t = get_local_xyz(xyz, i_xyz_t);
+  }
+  double gamma = xyz_t[0] * std::sqrt(3.0) / 2.0 - xyz_t[1] / 2.0;
+  if ((std::abs(gamma - edge) > FP_PRECISION) && gamma_dir != 0) {
+    double this_d = (edge - gamma) / gamma_dir;
+    if (this_d < d) {
+      if (gamma_dir > 0) {
+        lattice_trans = {1, -1, 0};
+      } else {
+        lattice_trans = {-1, 1, 0};
+      }
+      d = this_d;
+    }
+  }
+
+  // Upper and lower sides.
+  edge = -copysign(0.5*pitch[0], uvw[1]);
+  if (uvw[1] > 0) {
+    const int i_xyz_t[3] {i_xyz[0], i_xyz[1]+1, i_xyz[2]};
+    xyz_t = get_local_xyz(xyz, i_xyz_t);
+  } else {
+    const int i_xyz_t[3] {i_xyz[0], i_xyz[1]-1, i_xyz[2]};
+    xyz_t = get_local_xyz(xyz, i_xyz_t);
+  }
+  if ((std::abs(xyz_t[1] - edge) > FP_PRECISION) && uvw[1] != 0) {
+    double this_d = (edge - xyz_t[1]) / uvw[1];
+    if (this_d < d) {
+      if (uvw[1] > 0) {
+        lattice_trans = {0, 1, 0};
+      } else {
+        lattice_trans = {0, -1, 0};
+      }
+      d = this_d;
+    }
+  }
+
+  // Top and bottom sides
+  if (is_3d) {
+    double z {xyz[2]};
+    double w {uvw[2]};
+    double z0 {copysign(0.5 * pitch[1], w)};
+    if ((std::abs(z - z0) > FP_PRECISION) && w != 0) {
+      double this_d = (z0 - z) / w;
+      if (this_d < d) {
+        d = this_d;
+        if (w > 0) {
+          lattice_trans = {0, 0, 1};
+        } else {
+          lattice_trans = {0, 0, -1};
+        }
+        d = this_d;
+      }
+    }
+  }
+
+  return {d, lattice_trans};
+}
+
+//==============================================================================
+
+std::array<int, 3>
+HexLattice::get_indices(const double xyz[3]) const
+{
+  // Offset the xyz by the lattice center.
+  double xyz_o[3] {xyz[0] - center[0], xyz[1] - center[1], xyz[2]};
+  if (is_3d) {xyz_o[2] -= center[2];}
+
+  // Index the z direction.
+  std::array<int, 3> out;
+  if (is_3d) {
+    out[2] = static_cast<int>(std::ceil(xyz_o[2] / pitch[1] + 0.5 * n_axial))-1;
+  } else {
+    out[2] = 0;
+  }
+
+  // Convert coordinates into skewed bases.  The (x, alpha) basis is used to
+  // find the index of the global coordinates to within 4 cells.
+  double alpha = xyz_o[1] - xyz_o[0] / std::sqrt(3.0);
+  out[0] = static_cast<int>(std::floor(xyz_o[0]
+                                       / (0.5*std::sqrt(3.0) * pitch[0])));
+  out[1] = static_cast<int>(std::floor(alpha / pitch[0]));
+
+  // Add offset to indices (the center cell is (i_x, i_alpha) = (0, 0) but
+  // the array is offset so that the indices never go below 0).
+  out[0] += n_rings-1;
+  out[1] += n_rings-1;
+
+  // Calculate the (squared) distance between the particle and the centers of
+  // the four possible cells.  Regular hexagonal tiles form a Voronoi
+  // tessellation so the xyz should be in the hexagonal cell that it is closest
+  // to the center of.  This method is used over a method that uses the
+  // remainders of the floor divisions above because it provides better finite
+  // precision performance.  Squared distances are used becasue they are more
+  // computationally efficient than normal distances.
+  int k {1};
+  int k_min {1};
+  double d_min {INFTY};
+  for (int i = 0; i < 2; i++) {
+    for (int j = 0; j < 2; j++) {
+      int i_xyz[3] {out[0] + j, out[1] + i, 0};
+      std::array<double, 3> xyz_t = get_local_xyz(xyz, i_xyz);
+      double d = xyz_t[0]*xyz_t[0] + xyz_t[1]*xyz_t[1];
+      if (d < d_min) {
+        d_min = d;
+        k_min = k;
+      }
+      k++;
+    }
+  }
+
+  // Select the minimum squared distance which corresponds to the cell the
+  // coordinates are in.
+  if (k_min == 2) {
+    out[0] += 1;
+  } else if (k_min == 3) {
+    out[1] += 1;
+  } else if (k_min == 4) {
+    out[0] += 1;
+    out[1] += 1;
+  }
+
+  return out;
+}
+
+//==============================================================================
+
+std::array<double, 3>
+HexLattice::get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const
+{
+  std::array<double, 3> local_xyz;
+
+  // x_l = x_g - (center + pitch_x*cos(30)*index_x)
+  local_xyz[0] = global_xyz[0] - (center[0]
+       + std::sqrt(3.0)/2.0 * (i_xyz[0] - n_rings + 1) * pitch[0]);
+  // y_l = y_g - (center + pitch_x*index_x + pitch_y*sin(30)*index_y)
+  local_xyz[1] = global_xyz[1] - (center[1]
+       + (i_xyz[1] - n_rings + 1) * pitch[0]
+       + (i_xyz[0] - n_rings + 1) * pitch[0] / 2.0);
+  if (is_3d) {
+    local_xyz[2] = global_xyz[2] - center[2]
+         + (0.5 * n_axial - i_xyz[2] - 0.5) * pitch[1];
+  } else {
+    local_xyz[2] = global_xyz[2];
+  }
+
+  return local_xyz;
+}
+
+//==============================================================================
+
+bool
+HexLattice::is_valid_index(int indx) const
+{
+  int nx {2*n_rings - 1};
+  int ny {2*n_rings - 1};
+  int nz {n_axial};
+  int iz = indx / (nx * ny);
+  int iy = (indx - nx*ny*iz) / nx;
+  int ix = indx - nx*ny*iz - nx*iy;
+  int i_xyz[3] {ix, iy, iz};
+  return are_valid_indices(i_xyz);
+}
+
+//==============================================================================
+
+int32_t&
+HexLattice::offset(int map, const int i_xyz[3])
+{
+  int nx {2*n_rings - 1};
+  int ny {2*n_rings - 1};
+  int nz {n_axial};
+  return offsets[nx*ny*nz*map + nx*ny*i_xyz[2] + nx*i_xyz[1] + i_xyz[0]];
+}
+
+//==============================================================================
+
+std::string
+HexLattice::index_to_string(int indx) const
+{
+  int nx {2*n_rings - 1};
+  int ny {2*n_rings - 1};
+  int iz {indx / (nx * ny)};
+  int iy {(indx - nx * ny * iz) / nx};
+  int ix {indx - nx * ny * iz - nx * iy};
+  std::string out {std::to_string(ix - n_rings + 1)};
+  out += ',';
+  out += std::to_string(iy - n_rings + 1);
+  if (is_3d) {
+    out += ',';
+    out += std::to_string(iz);
+  }
+  return out;
+}
+
+//==============================================================================
+
+void
+HexLattice::to_hdf5_inner(hid_t lat_group) const
+{
+  // Write basic lattice information.
+  write_string(lat_group, "type", "hexagonal", false);
+  write_int(lat_group, 0, nullptr, "n_rings", &n_rings, false);
+  write_int(lat_group, 0, nullptr, "n_axial", &n_axial, false);
+  if (is_3d) {
+    write_double(lat_group, "pitch", pitch, false);
+    write_double(lat_group, "center", center, false);
+  } else {
+    std::array<double, 1> pitch_short {{pitch[0]}};
+    write_double(lat_group, "pitch", pitch_short, false);
+    std::array<double, 2> center_short {{center[0], center[1]}};
+    write_double(lat_group, "center", center_short, false);
+  }
+
+  // Write the universe ids.
+  hsize_t nx {static_cast<hsize_t>(2*n_rings - 1)};
+  hsize_t ny {static_cast<hsize_t>(2*n_rings - 1)};
+  hsize_t nz {static_cast<hsize_t>(n_axial)};
+  int out[nx*ny*nz];
+
+  for (int m = 0; m < nz; m++) {
+    for (int k = 0; k < ny; k++) {
+      for (int j = 0; j < nx; j++) {
+        int indx = nx*ny*m + nx*k + j;
+        if (j + k < n_rings - 1) {
+          // This array position is never used; put a -1 to indicate this.
+          out[indx] = -1;
+        } else if (j + k > 3*n_rings - 3) {
+          // This array position is never used; put a -1 to indicate this.
+          out[indx] = -1;
+        } else {
+          out[indx] = global_universes[universes[indx]]->id;
+        }
+      }
+    }
+  }
+
+  hsize_t dims[3] {nz, ny, nx};
+  write_int(lat_group, 3, dims, "universes", out, false);
+}
+
+//==============================================================================
+// Non-method functions
+//==============================================================================
+
+extern "C" void
+read_lattices(pugi::xml_node *node)
+{
+  for (pugi::xml_node lat_node : node->children("lattice")) {
+    lattices_c.push_back(new RectLattice(lat_node));
+  }
+  for (pugi::xml_node lat_node : node->children("hex_lattice")) {
+    lattices_c.push_back(new HexLattice(lat_node));
+  }
+
+  // Fill the lattice map.
+  for (int i_lat = 0; i_lat < lattices_c.size(); i_lat++) {
+    int id = lattices_c[i_lat]->id;
+    auto in_map = lattice_map.find(id);
+    if (in_map == lattice_map.end()) {
+      lattice_map[id] = i_lat;
+    } else {
+      std::stringstream err_msg;
+      err_msg << "Two or more lattices use the same unique ID: " << id;
+      fatal_error(err_msg);
+    }
+  }
+}
+
+//==============================================================================
+// Fortran compatibility functions
+//==============================================================================
+
+extern "C" {
+  Lattice* lattice_pointer(int lat_ind) {return lattices_c[lat_ind];}
+
+  int32_t lattice_id(Lattice *lat) {return lat->id;}
+
+  bool lattice_are_valid_indices(Lattice *lat, const int i_xyz[3])
+  {return lat->are_valid_indices(i_xyz);}
+
+  void lattice_distance(Lattice *lat, const double xyz[3], const double uvw[3],
+                        const int i_xyz[3], double *d, int lattice_trans[3])
+  {
+    std::pair<double, std::array<int, 3>> ld {lat->distance(xyz, uvw, i_xyz)};
+    *d = ld.first;
+    lattice_trans[0] = ld.second[0];
+    lattice_trans[1] = ld.second[1];
+    lattice_trans[2] = ld.second[2];
+  }
+
+  void lattice_get_indices(Lattice *lat, const double xyz[3], int i_xyz[3])
+  {
+    std::array<int, 3> inds = lat->get_indices(xyz);
+    i_xyz[0] = inds[0];
+    i_xyz[1] = inds[1];
+    i_xyz[2] = inds[2];
+  }
+
+  void lattice_get_local_xyz(Lattice *lat, const double global_xyz[3],
+                             const int i_xyz[3], double local_xyz[3])
+  {
+    std::array<double, 3> xyz = lat->get_local_xyz(global_xyz, i_xyz);
+    local_xyz[0] = xyz[0];
+    local_xyz[1] = xyz[1];
+    local_xyz[2] = xyz[2];
+  }
+
+  int32_t lattice_offset(Lattice *lat, int map, const int i_xyz[3])
+  {return lat->offset(map, i_xyz);}
+
+  int32_t lattice_outer(Lattice *lat) {return lat->outer;}
+
+  void lattice_to_hdf5(Lattice *lat, hid_t group) {lat->to_hdf5(group);}
+
+  int32_t lattice_universe(Lattice *lat, const int i_xyz[3])
+  {return (*lat)[i_xyz];}
+}
+
+} // namespace openmc

src/lattice.h

@@ -0,0 +1,259 @@
+#ifndef LATTICE_H
+#define LATTICE_H
+
+#include <array>
+#include <cstdint>
+#include <string>
+#include <unordered_map>
+#include <vector>
+
+#include "constants.h"
+#include "hdf5.h"
+#include "pugixml/pugixml.hpp"
+
+
+namespace openmc {
+
+//==============================================================================
+// Module constants
+//==============================================================================
+
+constexpr int32_t NO_OUTER_UNIVERSE{-1};
+
+//==============================================================================
+// Global variables
+//==============================================================================
+
+class Lattice;
+extern std::vector<Lattice*> lattices_c;
+
+extern std::unordered_map<int32_t, int32_t> lattice_map;
+
+//==============================================================================
+//! \class Lattice
+//! \brief Abstract type for ordered array of universes.
+//==============================================================================
+
+class LatticeIter;
+class ReverseLatticeIter;
+
+class Lattice
+{
+public:
+  int32_t id;                         //!< Universe ID number
+  std::string name;                   //!< User-defined name
+  std::vector<int32_t> universes;     //!< Universes filling each lattice tile
+  int32_t outer {NO_OUTER_UNIVERSE};  //!< Universe tiled outside the lattice
+  std::vector<int32_t> offsets;       //!< Distribcell offset table
+
+  explicit Lattice(pugi::xml_node lat_node);
+
+  virtual ~Lattice() {}
+
+  virtual int32_t& operator[](const int i_xyz[3]) = 0;
+
+  virtual LatticeIter begin();
+  LatticeIter end();
+
+  virtual ReverseLatticeIter rbegin();
+  ReverseLatticeIter rend();
+
+  //! Convert internal universe values from IDs to indices using universe_map.
+  void adjust_indices();
+
+  //! Allocate offset table for distribcell.
+  void allocate_offset_table(int n_maps)
+  {offsets.resize(n_maps * universes.size(), C_NONE);}
+
+  //! Populate the distribcell offset tables.
+  int32_t fill_offset_table(int32_t offset, int32_t target_univ_id, int map);
+
+  //! \brief Check lattice indices.
+  //! @param i_xyz[3] The indices for a lattice tile.
+  //! @return true if the given indices fit within the lattice bounds.  False
+  //!   otherwise.
+  virtual bool are_valid_indices(const int i_xyz[3]) const = 0;
+
+  //! \brief Find the next lattice surface crossing
+  //! @param xyz[3] A 3D Cartesian coordinate.
+  //! @param uvw[3] A 3D Cartesian direction.
+  //! @param i_xyz[3] The indices for a lattice tile.
+  //! @return The distance to the next crossing and an array indicating how the
+  //!   lattice indices would change after crossing that boundary.
+  virtual std::pair<double, std::array<int, 3>>
+  distance(const double xyz[3], const double uvw[3], const int i_xyz[3]) const
+  = 0;
+
+  //! \brief Find the lattice tile indices for a given point.
+  //! @param xyz[3] A 3D Cartesian coordinate.
+  //! @return An array containing the indices of a lattice tile.
+  virtual std::array<int, 3> get_indices(const double xyz[3]) const = 0;
+
+  //! \brief Get coordinates local to a lattice tile.
+  //! @param global_xyz[3] A 3D Cartesian coordinate.
+  //! @param i_xyz[3] The indices for a lattice tile.
+  //! @return Local 3D Cartesian coordinates.
+  virtual std::array<double, 3>
+  get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const = 0;
+
+  //! \brief Check flattened lattice index.
+  //! @param indx The index for a lattice tile.
+  //! @return true if the given index fit within the lattice bounds.  False
+  //!   otherwise.
+  virtual bool is_valid_index(int indx) const
+  {return (indx >= 0) && (indx < universes.size());}
+
+  //! \brief Get the distribcell offset for a lattice tile.
+  //! @param The map index for the target cell.
+  //! @param i_xyz[3] The indices for a lattice tile.
+  //! @return Distribcell offset i.e. the largest instance number for the target
+  //!  cell found in the geometry tree under this lattice tile.
+  virtual int32_t& offset(int map, const int i_xyz[3]) = 0;
+
+  //! \brief Convert an array index to a useful human-readable string.
+  //! @param indx The index for a lattice tile.
+  //! @return A string representing the lattice tile.
+  virtual std::string index_to_string(int indx) const = 0;
+
+  //! \brief Write lattice information to an HDF5 group.
+  //! @param group_id An HDF5 group id.
+  void to_hdf5(hid_t group_id) const;
+
+protected:
+  bool is_3d;  //!< Has divisions along the z-axis?
+
+  virtual void to_hdf5_inner(hid_t group_id) const = 0;
+};
+
+//==============================================================================
+//! An iterator over lattice universes.
+//==============================================================================
+
+class LatticeIter
+{
+public:
+  int indx;  //!< An index to a Lattice universes or offsets array.
+
+  LatticeIter(Lattice &lat_, int indx_)
+    : lat(lat_),
+      indx(indx_)
+  {}
+
+  bool operator==(const LatticeIter &rhs) {return (indx == rhs.indx);}
+
+  bool operator!=(const LatticeIter &rhs) {return !(*this == rhs);}
+
+  int32_t& operator*() {return lat.universes[indx];}
+
+  LatticeIter& operator++()
+  {
+    while (indx < lat.universes.size()) {
+      ++indx;
+      if (lat.is_valid_index(indx)) return *this;
+    }
+    indx = lat.universes.size();
+    return *this;
+  }
+
+protected:
+  Lattice &lat;
+};
+
+//==============================================================================
+//! A reverse iterator over lattice universes.
+//==============================================================================
+
+class ReverseLatticeIter : public LatticeIter
+{
+public:
+  ReverseLatticeIter(Lattice &lat_, int indx_)
+    : LatticeIter {lat_, indx_}
+  {}
+
+  ReverseLatticeIter& operator++()
+  {
+    while (indx > -1) {
+      --indx;
+      if (lat.is_valid_index(indx)) return *this;
+    }
+    indx = -1;
+    return *this;
+  }
+};
+
+//==============================================================================
+
+class RectLattice : public Lattice
+{
+public:
+  explicit RectLattice(pugi::xml_node lat_node);
+
+  int32_t& operator[](const int i_xyz[3]);
+
+  bool are_valid_indices(const int i_xyz[3]) const;
+
+  std::pair<double, std::array<int, 3>>
+  distance(const double xyz[3], const double uvw[3], const int i_xyz[3]) const;
+
+  std::array<int, 3> get_indices(const double xyz[3]) const;
+
+  std::array<double, 3>
+  get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const;
+
+  int32_t& offset(int map, const int i_xyz[3]);
+
+  std::string index_to_string(int indx) const;
+
+  void to_hdf5_inner(hid_t group_id) const;
+
+private:
+  std::array<int, 3> n_cells;       //!< Number of cells along each axis
+  std::array<double, 3> lower_left; //!< Global lower-left corner of the lattice
+  std::array<double, 3> pitch;      //!< Lattice tile width along each axis
+
+  // Convenience aliases
+  int &nx {n_cells[0]};
+  int &ny {n_cells[1]};
+  int &nz {n_cells[2]};
+};
+
+//==============================================================================
+
+class HexLattice : public Lattice
+{
+public:
+  explicit HexLattice(pugi::xml_node lat_node);
+
+  int32_t& operator[](const int i_xyz[3]);
+
+  LatticeIter begin();
+
+  ReverseLatticeIter rbegin();
+
+  bool are_valid_indices(const int i_xyz[3]) const;
+
+  std::pair<double, std::array<int, 3>>
+  distance(const double xyz[3], const double uvw[3], const int i_xyz[3]) const;
+
+  std::array<int, 3> get_indices(const double xyz[3]) const;
+
+  std::array<double, 3>
+  get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const;
+
+  bool is_valid_index(int indx) const;
+
+  int32_t& offset(int map, const int i_xyz[3]);
+
+  std::string index_to_string(int indx) const;
+
+  void to_hdf5_inner(hid_t group_id) const;
+
+private:
+  int n_rings;                   //!< Number of radial tile positions
+  int n_axial;                   //!< Number of axial tile positions
+  std::array<double, 3> center;  //!< Global center of lattice
+  std::array<double, 2> pitch;   //!< Lattice tile width and height
+};
+
+} //  namespace openmc
+#endif // LATTICE_H

src/nuclide_header.F90

@@ -10,7 +10,6 @@ module nuclide_header
   use endf_header,            only: Function1D, Polynomial, Tabulated1D
   use error
   use hdf5_interface
-  use list_header,            only: ListInt
   use math,                   only: faddeeva, w_derivative, &
                                     broaden_wmp_polynomials
   use multipole_header,       only: FORM_RM, FORM_MLBW, MP_EA, RM_RT, RM_RA, &

src/output.F90

@@ -234,7 +234,7 @@ contains
       ! Print cell for this level
       if (p % coord(i) % cell /= NONE) then
         c => cells(p % coord(i) % cell)
-        write(ou,*) '    Cell             = ' // trim(to_str(c % id))
+        write(ou,*) '    Cell             = ' // trim(to_str(c % id()))
       end if
 
       ! Print universe for this level
@@ -246,7 +246,7 @@ contains
       ! Print information on lattice
       if (p % coord(i) % lattice /= NONE) then
         l => lattices(p % coord(i) % lattice) % obj
-        write(ou,*) '    Lattice          = ' // trim(to_str(l % id))
+        write(ou,*) '    Lattice          = ' // trim(to_str(l % id()))
         write(ou,*) '    Lattice position = (' // trim(to_str(&
              p % coord(i) % lattice_x)) // ',' // trim(to_str(&
              p % coord(i) % lattice_y)) // ')'
@@ -258,7 +258,7 @@ contains
     end do
 
     ! Print surface
-    if (p % surface /= NONE) then
+    if (p % surface /= ERROR_INT) then
       write(ou,*) '  Surface = ' // to_str(sign(surfaces(i)%id(), p % surface))
     end if
 
@@ -633,7 +633,7 @@ contains
     write(ou,100) 'Cell ID','No. Overlap Checks'
 
     do i = 1, n_cells
-      write(ou,101) cells(i) % id, overlap_check_cnt(i)
+      write(ou,101) cells(i) % id(), overlap_check_cnt(i)
       if (overlap_check_cnt(i) < 10) num_sparse = num_sparse + 1
     end do
     write(ou,*)
@@ -643,7 +643,7 @@ contains
     do i = 1, n_cells
       if (overlap_check_cnt(i) < 10) then
         j = j + 1
-        write(ou,'(1X,A8)', advance='no') trim(to_str(cells(i) % id))
+        write(ou,'(1X,A8)', advance='no') trim(to_str(cells(i) % id()))
         if (modulo(j,8) == 0) write(ou,*)
       end if
     end do

src/plot.F90

@@ -85,7 +85,7 @@ contains
       if (pl % color_by == PLOT_COLOR_MATS) then
         ! Assign color based on material
         associate (c => cells(p % coord(j) % cell))
-          if (c % type == FILL_UNIVERSE) then
+          if (c % type() == FILL_UNIVERSE) then
             ! If we stopped on a middle universe level, treat as if not found
             rgb = pl % not_found % rgb
             id = -1
@@ -101,7 +101,7 @@ contains
       else if (pl % color_by == PLOT_COLOR_CELLS) then
         ! Assign color based on cell
         rgb = pl % colors(p % coord(j) % cell) % rgb
-        id = cells(p % coord(j) % cell) % id
+        id = cells(p % coord(j) % cell) % id()
       else
         rgb = 0
         id = -1

src/simulation_header.F90

@@ -72,9 +72,6 @@ module simulation_header
 
   logical :: trace
 
-  ! Number of distribcell maps
-  integer :: n_maps
-
 !$omp threadprivate(trace, thread_id, current_work)
 
 contains

src/string_utils.h

@@ -0,0 +1,35 @@
+#ifndef STRING_UTILS_H
+#define STRING_UTILS_H
+
+#include <string>
+#include <vector>
+
+
+namespace openmc {
+
+std::vector<std::string>
+split(const std::string &in)
+{
+  std::vector<std::string> out;
+
+  for (int i = 0; i < in.size(); ) {
+    // Increment i until we find a non-whitespace character.
+    if (std::isspace(in[i])) {
+      i++;
+
+    } else {
+      // Find the next whitespace character at j.
+      int j = i + 1;
+      while (j < in.size() && std::isspace(in[j]) == 0) {j++;}
+
+      // Push-back everything between i and j.
+      out.push_back(in.substr(i, j-i));
+      i = j + 1; // j is whitespace so leapfrog to j+1
+    }
+  }
+
+  return out;
+}
+
+} // namespace openmc
+#endif // STRING_UTILS_H

src/summary.F90

@@ -15,7 +15,6 @@ module summary
   use surface_header
   use string,          only: to_str
   use tally_header,    only: TallyObject
-  use tally_filter_distribcell, only: find_offset
 
   implicit none
   private
@@ -160,8 +159,7 @@ contains
   subroutine write_geometry(file_id)
     integer(HID_T), intent(in) :: file_id
 
-    integer :: i, j, k, m
-    integer, allocatable :: lattice_universes(:,:,:)
+    integer :: i, j, cell_fill
     integer, allocatable :: cell_materials(:)
     integer, allocatable :: cell_ids(:)
     real(8), allocatable :: cell_temperatures(:)
@@ -169,8 +167,7 @@ contains
     integer(HID_T) :: cells_group, cell_group
     integer(HID_T) :: surfaces_group
     integer(HID_T) :: universes_group, univ_group
-    integer(HID_T) :: lattices_group, lattice_group
-    character(:), allocatable :: region_spec
+    integer(HID_T) :: lattices_group
     type(Cell),     pointer :: c
     class(Lattice), pointer :: lat
 
@@ -179,7 +176,7 @@ contains
     call write_attribute(geom_group, "n_cells", n_cells)
     call write_attribute(geom_group, "n_surfaces", n_surfaces)
     call write_attribute(geom_group, "n_universes", n_universes)
-    call write_attribute(geom_group, "n_lattices", n_lattices)
+    call write_attribute(geom_group, "n_lattices", size(lattices))
 
     ! ==========================================================================
     ! WRITE INFORMATION ON CELLS
@@ -190,16 +187,12 @@ contains
     ! Write information on each cell
     CELL_LOOP: do i = 1, n_cells
       c => cells(i)
-      cell_group = create_group(cells_group, "cell " // trim(to_str(c%id)))
+      cell_group = create_group(cells_group, "cell " // trim(to_str(c%id())))
 
-      ! Write name for this cell
-      call write_dataset(cell_group, "name", c%name)
-
-      ! Write universe for this cell
-      call write_dataset(cell_group, "universe", universes(c%universe)%id)
+      call c % to_hdf5(cell_group)
 
       ! Write information on what fills this cell
-      select case (c%type)
+      select case (c%type())
       case (FILL_MATERIAL)
         call write_dataset(cell_group, "fill_type", "material")
 
@@ -231,12 +224,7 @@ contains
 
       case (FILL_UNIVERSE)
         call write_dataset(cell_group, "fill_type", "universe")
-        call write_dataset(cell_group, "fill", universes(c%fill)%id)
-        if (allocated(c%offset)) then
-          if (size(c%offset) > 0) then
-            call write_dataset(cell_group, "offset", c%offset)
-          end if
-        end if
+        call write_dataset(cell_group, "fill", universes(c%fill()+1)%id)
 
         if (allocated(c%translation)) then
           call write_dataset(cell_group, "translation", c%translation)
@@ -247,30 +235,12 @@ contains
 
       case (FILL_LATTICE)
         call write_dataset(cell_group, "fill_type", "lattice")
-        call write_dataset(cell_group, "lattice", lattices(c%fill)%obj%id)
+        ! Do not access the 'lattices' array with 'c % fill() + 1' directly; it
+        ! causes a segfault in GCC 7.3.0
+        cell_fill = c % fill() + 1
+        call write_dataset(cell_group, "lattice", lattices(cell_fill)%obj%id())
       end select
 
-      ! Write list of bounding surfaces
-      region_spec = ""
-      do j = 1, size(c%region)
-        k = c%region(j)
-        select case(k)
-        case (OP_LEFT_PAREN)
-          region_spec = trim(region_spec) // " ("
-        case (OP_RIGHT_PAREN)
-          region_spec = trim(region_spec) // " )"
-        case (OP_COMPLEMENT)
-          region_spec = trim(region_spec) // " ~"
-        case (OP_INTERSECTION)
-        case (OP_UNION)
-          region_spec = trim(region_spec) // " |"
-        case default
-          region_spec = trim(region_spec) // " " // to_str(&
-               sign(surfaces(abs(k))%id(), k))
-        end select
-      end do
-      call write_dataset(cell_group, "region", adjustl(region_spec))
-
       call close_group(cell_group)
     end do CELL_LOOP
 
@@ -305,7 +275,7 @@ contains
         if (size(u % cells) > 0) then
           allocate(cell_ids(size(u % cells)))
           do j = 1, size(u % cells)
-            cell_ids(j) = cells(u % cells(j)) % id
+            cell_ids(j) = cells(u % cells(j)) % id()
           end do
           call write_dataset(univ_group, "cells", cell_ids)
           deallocate(cell_ids)
@@ -320,87 +290,12 @@ contains
     ! ==========================================================================
     ! WRITE INFORMATION ON LATTICES
 
-    ! Create lattices group (nothing directly written here) then close
     lattices_group = create_group(geom_group, "lattices")
 
-    ! Write information on each lattice
-    LATTICE_LOOP: do i = 1, n_lattices
+    do i = 1, size(lattices)
       lat => lattices(i)%obj
-      lattice_group = create_group(lattices_group, "lattice " // trim(to_str(lat%id)))
-
-      ! Write name, pitch, and outer universe
-      call write_dataset(lattice_group, "name", lat%name)
-      call write_dataset(lattice_group, "pitch", lat%pitch)
-      if (lat % outer > 0) then
-        call write_dataset(lattice_group, "outer", universes(lat % outer) % id)
-      else
-        call write_dataset(lattice_group, "outer", lat % outer)
-      end if
-
-      select type (lat)
-      type is (RectLattice)
-        ! Write lattice type.
-        call write_dataset(lattice_group, "type", "rectangular")
-
-        ! Write lattice dimensions, lower left corner, and pitch
-        if (lat % is_3d) then
-          call write_dataset(lattice_group, "dimension", lat % n_cells)
-          call write_dataset(lattice_group, "lower_left", lat % lower_left)
-        else
-          call write_dataset(lattice_group, "dimension", lat % n_cells(1:2))
-          call write_dataset(lattice_group, "lower_left", lat % lower_left)
-        end if
-
-        ! Write lattice universes.
-        allocate(lattice_universes(lat%n_cells(1), lat%n_cells(2), &
-             &lat%n_cells(3)))
-        do j = 1, lat%n_cells(1)
-          do k = 0, lat%n_cells(2) - 1
-            do m = 1, lat%n_cells(3)
-              lattice_universes(j, k+1, m) = &
-                   universes(lat%universes(j, lat%n_cells(2) - k, m))%id
-            end do
-          end do
-        end do
-
-      type is (HexLattice)
-        ! Write lattice type.
-        call write_dataset(lattice_group, "type", "hexagonal")
-
-        ! Write number of lattice cells.
-        call write_dataset(lattice_group, "n_rings", lat%n_rings)
-        call write_dataset(lattice_group, "n_axial", lat%n_axial)
-
-        ! Write lattice center
-        call write_dataset(lattice_group, "center", lat%center)
-
-        ! Write lattice universes.
-        allocate(lattice_universes(2*lat%n_rings - 1, 2*lat%n_rings - 1, &
-             &lat%n_axial))
-        do m = 1, lat%n_axial
-          do k = 1, 2*lat%n_rings - 1
-            do j = 1, 2*lat%n_rings - 1
-              if (j + k < lat%n_rings + 1) then
-                ! This array position is never used; put a -1 to indicate this
-                lattice_universes(j,k,m) = -1
-                cycle
-              else if (j + k > 3*lat%n_rings - 1) then
-                ! This array position is never used; put a -1 to indicate this
-                lattice_universes(j,k,m) = -1
-                cycle
-              end if
-              lattice_universes(j,k,m) = universes(lat%universes(j,k,m))%id
-            end do
-          end do
-        end do
-      end select
-
-      ! Write lattice universes
-      call write_dataset(lattice_group, "universes", lattice_universes)
-      deallocate(lattice_universes)
-
-      call close_group(lattice_group)
-    end do LATTICE_LOOP
+      call lat % to_hdf5(lattices_group)
+    end do
 
     call close_group(lattices_group)
     call close_group(geom_group)

src/surface.cpp

@@ -11,12 +11,25 @@
 
 namespace openmc {
 
+//==============================================================================
+// Module constant definitions
+//==============================================================================
+
+extern "C" const int BC_TRANSMIT {0};
+extern "C" const int BC_VACUUM {1};
+extern "C" const int BC_REFLECT {2};
+extern "C" const int BC_PERIODIC {3};
+
 //==============================================================================
 // Global variables
 //==============================================================================
 
 int32_t n_surfaces;
 
+Surface **surfaces_c;
+
+std::map<int, int> surface_map;
+
 //==============================================================================
 // Helper functions for reading the "coeffs" node of an XML surface element
 //==============================================================================
@@ -305,7 +318,7 @@ void SurfaceXPlane::to_hdf5_inner(hid_t group_id) const
 {
   write_string(group_id, "type", "x-plane", false);
   std::array<double, 1> coeffs {{x0}};
-  write_double_1D(group_id, "coefficients", coeffs);
+  write_double(group_id, "coefficients", coeffs, false);
 }
 
 bool SurfaceXPlane::periodic_translate(PeriodicSurface *other, double xyz[3],
@@ -370,7 +383,7 @@ void SurfaceYPlane::to_hdf5_inner(hid_t group_id) const
 {
   write_string(group_id, "type", "y-plane", false);
   std::array<double, 1> coeffs {{y0}};
-  write_double_1D(group_id, "coefficients", coeffs);
+  write_double(group_id, "coefficients", coeffs, false);
 }
 
 bool SurfaceYPlane::periodic_translate(PeriodicSurface *other, double xyz[3],
@@ -436,7 +449,7 @@ void SurfaceZPlane::to_hdf5_inner(hid_t group_id) const
 {
   write_string(group_id, "type", "z-plane", false);
   std::array<double, 1> coeffs {{z0}};
-  write_double_1D(group_id, "coefficients", coeffs);
+  write_double(group_id, "coefficients", coeffs, false);
 }
 
 bool SurfaceZPlane::periodic_translate(PeriodicSurface *other, double xyz[3],
@@ -497,7 +510,7 @@ void SurfacePlane::to_hdf5_inner(hid_t group_id) const
 {
   write_string(group_id, "type", "plane", false);
   std::array<double, 4> coeffs {{A, B, C, D}};
-  write_double_1D(group_id, "coefficients", coeffs);
+  write_double(group_id, "coefficients", coeffs, false);
 }
 
 bool SurfacePlane::periodic_translate(PeriodicSurface *other, double xyz[3],
@@ -629,7 +642,7 @@ void SurfaceXCylinder::to_hdf5_inner(hid_t group_id) const
 {
   write_string(group_id, "type", "x-cylinder", false);
   std::array<double, 3> coeffs {{y0, z0, r}};
-  write_double_1D(group_id, "coefficients", coeffs);
+  write_double(group_id, "coefficients", coeffs, false);
 }
 
 //==============================================================================
@@ -663,7 +676,7 @@ void SurfaceYCylinder::to_hdf5_inner(hid_t group_id) const
 {
   write_string(group_id, "type", "y-cylinder", false);
   std::array<double, 3> coeffs {{x0, z0, r}};
-  write_double_1D(group_id, "coefficients", coeffs);
+  write_double(group_id, "coefficients", coeffs, false);
 }
 
 //==============================================================================
@@ -697,7 +710,7 @@ void SurfaceZCylinder::to_hdf5_inner(hid_t group_id) const
 {
   write_string(group_id, "type", "z-cylinder", false);
   std::array<double, 3> coeffs {{x0, y0, r}};
-  write_double_1D(group_id, "coefficients", coeffs);
+  write_double(group_id, "coefficients", coeffs, false);
 }
 
 //==============================================================================
@@ -768,7 +781,7 @@ void SurfaceSphere::to_hdf5_inner(hid_t group_id) const
 {
   write_string(group_id, "type", "sphere", false);
   std::array<double, 4> coeffs {{x0, y0, z0, r}};
-  write_double_1D(group_id, "coefficients", coeffs);
+  write_double(group_id, "coefficients", coeffs, false);
 }
 
 //==============================================================================
@@ -885,7 +898,7 @@ void SurfaceXCone::to_hdf5_inner(hid_t group_id) const
 {
   write_string(group_id, "type", "x-cone", false);
   std::array<double, 4> coeffs {{x0, y0, z0, r_sq}};
-  write_double_1D(group_id, "coefficients", coeffs);
+  write_double(group_id, "coefficients", coeffs, false);
 }
 
 //==============================================================================
@@ -919,7 +932,7 @@ void SurfaceYCone::to_hdf5_inner(hid_t group_id) const
 {
   write_string(group_id, "type", "y-cone", false);
   std::array<double, 4> coeffs {{x0, y0, z0, r_sq}};
-  write_double_1D(group_id, "coefficients", coeffs);
+  write_double(group_id, "coefficients", coeffs, false);
 }
 
 //==============================================================================
@@ -953,7 +966,7 @@ void SurfaceZCone::to_hdf5_inner(hid_t group_id) const
 {
   write_string(group_id, "type", "z-cone", false);
   std::array<double, 4> coeffs {{x0, y0, z0, r_sq}};
-  write_double_1D(group_id, "coefficients", coeffs);
+  write_double(group_id, "coefficients", coeffs, false);
 }
 
 //==============================================================================
@@ -1047,7 +1060,7 @@ void SurfaceQuadric::to_hdf5_inner(hid_t group_id) const
 {
   write_string(group_id, "type", "quadric", false);
   std::array<double, 10> coeffs {{A, B, C, D, E, F, G, H, J, K}};
-  write_double_1D(group_id, "coefficients", coeffs);
+  write_double(group_id, "coefficients", coeffs, false);
 }
 
 //==============================================================================
@@ -1116,12 +1129,12 @@ read_surfaces(pugi::xml_node *node)
     }
   }
 
-  // Fill the surface dictionary.
+  // Fill the surface map.
   for (int i_surf = 0; i_surf < n_surfaces; i_surf++) {
     int id = surfaces_c[i_surf]->id;
-    auto in_dict = surface_dict.find(id);
-    if (in_dict == surface_dict.end()) {
-      surface_dict[id] = i_surf;
+    auto in_map = surface_map.find(id);
+    if (in_map == surface_map.end()) {
+      surface_map[id] = i_surf;
     } else {
       std::stringstream err_msg;
       err_msg << "Two or more surfaces use the same unique ID: " << id;
@@ -1211,7 +1224,7 @@ read_surfaces(pugi::xml_node *node)
           }
         } else {
           // Convert the surface id to an index.
-          surf->i_periodic = surface_dict[surf->i_periodic];
+          surf->i_periodic = surface_map[surf->i_periodic];
         }
       } else {
         // This is a SurfacePlane.  We won't try to find it's partner if the
@@ -1223,7 +1236,7 @@ read_surfaces(pugi::xml_node *node)
           fatal_error(err_msg);
         } else {
           // Convert the surface id to an index.
-          surf->i_periodic = surface_dict[surf->i_periodic];
+          surf->i_periodic = surface_map[surf->i_periodic];
         }
       }
 
@@ -1242,43 +1255,36 @@ read_surfaces(pugi::xml_node *node)
 // Fortran compatibility functions
 //==============================================================================
 
-extern "C" Surface* surface_pointer(int surf_ind) {return surfaces_c[surf_ind];}
+extern "C" {
+  Surface* surface_pointer(int surf_ind) {return surfaces_c[surf_ind];}
 
-extern "C" int surface_id(Surface *surf) {return surf->id;}
+  int surface_id(Surface *surf) {return surf->id;}
 
-extern "C" int surface_bc(Surface *surf) {return surf->bc;}
+  int surface_bc(Surface *surf) {return surf->bc;}
 
-extern "C" bool surface_sense(Surface *surf, double xyz[3], double uvw[3])
-{return surf->sense(xyz, uvw);}
+  void surface_reflect(Surface *surf, double xyz[3], double uvw[3])
+  {surf->reflect(xyz, uvw);}
 
-extern "C" void surface_reflect(Surface *surf, double xyz[3], double uvw[3])
-{surf->reflect(xyz, uvw);}
+  void surface_normal(Surface *surf, double xyz[3], double uvw[3])
+  {return surf->normal(xyz, uvw);}
 
-extern "C" double
-surface_distance(Surface *surf, double xyz[3], double uvw[3], bool coincident)
-{return surf->distance(xyz, uvw, coincident);}
+  void surface_to_hdf5(Surface *surf, hid_t group) {surf->to_hdf5(group);}
 
-extern "C" void surface_normal(Surface *surf, double xyz[3], double uvw[3])
-{return surf->normal(xyz, uvw);}
+  int surface_i_periodic(PeriodicSurface *surf) {return surf->i_periodic;}
 
-extern "C" void surface_to_hdf5(Surface *surf, hid_t group)
-{surf->to_hdf5(group);}
+  bool
+  surface_periodic(PeriodicSurface *surf, PeriodicSurface *other, double xyz[3],
+                   double uvw[3])
+  {return surf->periodic_translate(other, xyz, uvw);}
 
-extern "C" int surface_i_periodic(PeriodicSurface *surf)
-{return surf->i_periodic;}
-
-extern "C" bool
-surface_periodic(PeriodicSurface *surf, PeriodicSurface *other, double xyz[3],
-                 double uvw[3])
-{return surf->periodic_translate(other, xyz, uvw);}
-
-extern "C" void free_memory_surfaces_c()
-{
-  for (int i = 0; i < n_surfaces; i++) {delete surfaces_c[i];}
-  delete surfaces_c;
-  surfaces_c = nullptr;
-  n_surfaces = 0;
-  surface_dict.clear();
+  void free_memory_surfaces_c()
+  {
+    for (int i = 0; i < n_surfaces; i++) {delete surfaces_c[i];}
+    delete surfaces_c;
+    surfaces_c = nullptr;
+    n_surfaces = 0;
+    surface_map.clear();
+  }
 }
 
 } // namespace openmc

src/surface.h

@@ -8,25 +8,19 @@
 #include "hdf5.h"
 #include "pugixml/pugixml.hpp"
 
+#include "constants.h"
+
 
 namespace openmc {
 
 //==============================================================================
-// Module constants
+// Module constant declarations (defined in .cpp)
 //==============================================================================
 
-extern "C" const int BC_TRANSMIT {0};
-extern "C" const int BC_VACUUM {1};
-extern "C" const int BC_REFLECT {2};
-extern "C" const int BC_PERIODIC {3};
-
-//==============================================================================
-// Constants that should eventually be moved out of this file
-//==============================================================================
-
-extern "C" double FP_COINCIDENT;
-constexpr double INFTY {std::numeric_limits<double>::max()};
-constexpr int C_NONE {-1};
+extern "C" const int BC_TRANSMIT;
+extern "C" const int BC_VACUUM;
+extern "C" const int BC_REFLECT;
+extern "C" const int BC_PERIODIC;
 
 //==============================================================================
 // Global variables
@@ -35,9 +29,9 @@ constexpr int C_NONE {-1};
 extern "C" int32_t n_surfaces;
 
 class Surface;
-Surface **surfaces_c;
+extern Surface **surfaces_c;
 
-std::map<int, int> surface_dict;
+extern std::map<int, int> surface_map;
 
 //==============================================================================
 //! Coordinates for an axis-aligned cube that bounds a geometric object.
@@ -64,7 +58,7 @@ public:
   //int neighbor_pos[],        //!< List of cells on positive side
   //    neighbor_neg[];        //!< List of cells on negative side
   int bc;                    //!< Boundary condition
-  std::string name{""};      //!< User-defined name
+  std::string name;          //!< User-defined name
 
   explicit Surface(pugi::xml_node surf_node);
 
@@ -107,6 +101,7 @@ public:
 
   //! Write all information needed to reconstruct the surface to an HDF5 group.
   //! @param group_id An HDF5 group id.
+  //TODO: this probably needs to include i_periodic for PeriodicSurface
   void to_hdf5(hid_t group_id) const;
 
 protected:
@@ -383,19 +378,21 @@ public:
 // Fortran compatibility functions
 //==============================================================================
 
-extern "C" Surface* surface_pointer(int surf_ind);
-extern "C" int surface_id(Surface *surf);
-extern "C" int surface_bc(Surface *surf);
-extern "C" bool surface_sense(Surface *surf, double xyz[3], double uvw[3]);
-extern "C" void surface_reflect(Surface *surf, double xyz[3], double uvw[3]);
-extern "C" double surface_distance(Surface *surf, double xyz[3], double uvw[3],
-                                   bool coincident);
-extern "C" void surface_normal(Surface *surf, double xyz[3], double uvw[3]);
-extern "C" void surface_to_hdf5(Surface *surf, hid_t group);
-extern "C" int surface_i_periodic(PeriodicSurface *surf);
-extern "C" bool surface_periodic(PeriodicSurface *surf, PeriodicSurface *other,
-                                 double xyz[3], double uvw[3]);
-extern "C" void free_memory_surfaces_c();
+extern "C" {
+    Surface* surface_pointer(int surf_ind);
+    int surface_id(Surface *surf);
+    int surface_bc(Surface *surf);
+    bool surface_sense(Surface *surf, double xyz[3], double uvw[3]);
+    void surface_reflect(Surface *surf, double xyz[3], double uvw[3]);
+    double surface_distance(Surface *surf, double xyz[3], double uvw[3],
+                            bool coincident);
+    void surface_normal(Surface *surf, double xyz[3], double uvw[3]);
+    void surface_to_hdf5(Surface *surf, hid_t group);
+    int surface_i_periodic(PeriodicSurface *surf);
+    bool surface_periodic(PeriodicSurface *surf, PeriodicSurface *other,
+                          double xyz[3], double uvw[3]);
+    void free_memory_surfaces_c();
+}
 
 } // namespace openmc
 #endif // SURFACE_H

src/surface_header.F90

@@ -30,16 +30,6 @@ module surface_header
       integer(C_INT)                 :: bc
     end function surface_bc_c
 
-    pure function surface_sense_c(surf_ptr, xyz, uvw) &
-         bind(C, name='surface_sense') result(sense)
-      use ISO_C_BINDING
-      implicit none
-      type(C_PTR),    intent(in), value :: surf_ptr
-      real(C_DOUBLE), intent(in)        :: xyz(3)
-      real(C_DOUBLE), intent(in)        :: uvw(3)
-      logical(C_BOOL)                   :: sense
-    end function surface_sense_c
-
     pure subroutine surface_reflect_c(surf_ptr, xyz, uvw) &
          bind(C, name='surface_reflect')
       use ISO_C_BINDING
@@ -49,17 +39,6 @@ module surface_header
       real(C_DOUBLE), intent(inout)     :: uvw(3);
     end subroutine surface_reflect_c
 
-    pure function surface_distance_c(surf_ptr, xyz, uvw, coincident) &
-         bind(C, name='surface_distance') result(d)
-      use ISO_C_BINDING
-      implicit none
-      type(C_PTR),     intent(in), value :: surf_ptr
-      real(C_DOUBLE),  intent(in)        :: xyz(3);
-      real(C_DOUBLE),  intent(in)        :: uvw(3);
-      logical(C_BOOL), intent(in), value :: coincident;
-      real(C_DOUBLE)                     :: d;
-    end function surface_distance_c
-
     pure subroutine surface_normal_c(surf_ptr, xyz, uvw) &
          bind(C, name='surface_normal')
       use ISO_C_BINDING
@@ -115,9 +94,7 @@ module surface_header
 
     procedure :: id => surface_id
     procedure :: bc => surface_bc
-    procedure :: sense => surface_sense
     procedure :: reflect => surface_reflect
-    procedure :: distance => surface_distance
     procedure :: normal => surface_normal
     procedure :: to_hdf5 => surface_to_hdf5
     procedure :: i_periodic => surface_i_periodic
@@ -146,14 +123,6 @@ contains
     bc = surface_bc_c(this % ptr)
   end function surface_bc
 
-  pure function surface_sense(this, xyz, uvw) result(sense)
-    class(Surface), intent(in) :: this
-    real(C_DOUBLE), intent(in) :: xyz(3)
-    real(C_DOUBLE), intent(in) :: uvw(3)
-    logical(C_BOOL)            :: sense
-    sense = surface_sense_c(this % ptr, xyz, uvw)
-  end function surface_sense
-
   pure subroutine surface_reflect(this, xyz, uvw)
     class(Surface), intent(in)    :: this
     real(C_DOUBLE), intent(in)    :: xyz(3);
@@ -161,15 +130,6 @@ contains
     call surface_reflect_c(this % ptr, xyz, uvw)
   end subroutine surface_reflect
 
-  pure function surface_distance(this, xyz, uvw, coincident) result(d)
-    class(Surface),  intent(in) :: this
-    real(C_DOUBLE),  intent(in) :: xyz(3);
-    real(C_DOUBLE),  intent(in) :: uvw(3);
-    logical(C_BOOL), intent(in) :: coincident;
-    real(C_DOUBLE)              :: d;
-    d = surface_distance_c(this % ptr, xyz, uvw, coincident)
-  end function surface_distance
-
   pure subroutine surface_normal(this, xyz, uvw)
     class(Surface), intent(in)  :: this
     real(C_DOUBLE), intent(in)  :: xyz(3);

src/tallies/tally_filter_cell.F90

@@ -80,7 +80,7 @@ contains
 
     allocate(cell_ids(size(this % cells)))
     do i = 1, size(this % cells)
-      cell_ids(i) = cells(this % cells(i)) % id
+      cell_ids(i) = cells(this % cells(i)) % id()
     end do
     call write_dataset(filter_group, "bins", cell_ids)
   end subroutine to_statepoint_cell
@@ -114,7 +114,7 @@ contains
     integer,           intent(in) :: bin
     character(MAX_LINE_LEN)       :: label
 
-    label = "Cell " // to_str(cells(this % cells(bin)) % id)
+    label = "Cell " // to_str(cells(this % cells(bin)) % id())
   end function text_label_cell
 
 !===============================================================================

src/tallies/tally_filter_cellborn.F90

@@ -74,7 +74,7 @@ contains
     call write_dataset(filter_group, "n_bins", this % n_bins)
     allocate(cell_ids(size(this % cells)))
     do i = 1, size(this % cells)
-      cell_ids(i) = cells(this % cells(i)) % id
+      cell_ids(i) = cells(this % cells(i)) % id()
     end do
     call write_dataset(filter_group, "bins", cell_ids)
   end subroutine to_statepoint_cellborn
@@ -108,7 +108,7 @@ contains
     integer,               intent(in) :: bin
     character(MAX_LINE_LEN)           :: label
 
-    label = "Birth Cell " // to_str(cells(this % cells(bin)) % id)
+    label = "Birth Cell " // to_str(cells(this % cells(bin)) % id())
   end function text_label_cellborn
 
 end module tally_filter_cellborn

src/tallies/tally_filter_cellfrom.F90

@@ -64,7 +64,7 @@ contains
 
     allocate(cell_ids(size(this % cells)))
     do i = 1, size(this % cells)
-      cell_ids(i) = cells(this % cells(i)) % id
+      cell_ids(i) = cells(this % cells(i)) % id()
     end do
     call write_dataset(filter_group, "bins", cell_ids)
   end subroutine to_statepoint_cell_from
@@ -74,7 +74,7 @@ contains
     integer,               intent(in) :: bin
     character(MAX_LINE_LEN)           :: label
 
-    label = "Cell from " // to_str(cells(this % cells(bin)) % id)
+    label = "Cell from " // to_str(cells(this % cells(bin)) % id())
   end function text_label_cell_from
 
 end module tally_filter_cellfrom

src/tallies/tally_filter_distribcell.F90

@@ -14,7 +14,6 @@ module tally_filter_distribcell
 
   implicit none
   private
-  public :: find_offset
 
 !===============================================================================
 ! DISTRIBCELLFILTER specifies which distributed geometric cells tally events
@@ -58,20 +57,20 @@ contains
       distribcell_index = cells(this % cell) % distribcell_index
       offset = 0
       do i = 1, p % n_coord
-        if (cells(p % coord(i) % cell) % type == FILL_UNIVERSE) then
-          offset = offset + cells(p % coord(i) % cell) % &
-               offset(distribcell_index)
-        elseif (cells(p % coord(i) % cell) % type == FILL_LATTICE) then
+        if (cells(p % coord(i) % cell) % type() == FILL_UNIVERSE) then
+          offset = offset + cells(p % coord(i) % cell) &
+               % offset(distribcell_index-1)
+        elseif (cells(p % coord(i) % cell) % type() == FILL_LATTICE) then
           if (lattices(p % coord(i + 1) % lattice) % obj &
                % are_valid_indices([&
                p % coord(i + 1) % lattice_x, &
                p % coord(i + 1) % lattice_y, &
                p % coord(i + 1) % lattice_z])) then
-            offset = offset + lattices(p % coord(i + 1) % lattice) % obj % &
-                 offset(distribcell_index, &
-                 p % coord(i + 1) % lattice_x, &
+            offset = offset + lattices(p % coord(i + 1) % lattice) % obj &
+                 % offset(distribcell_index - 1, &
+                 [p % coord(i + 1) % lattice_x, &
                  p % coord(i + 1) % lattice_y, &
-                 p % coord(i + 1) % lattice_z)
+                 p % coord(i + 1) % lattice_z])
           end if
         end if
         if (this % cell == p % coord(i) % cell) then
@@ -88,7 +87,7 @@ contains
 
     call write_dataset(filter_group, "type", "distribcell")
     call write_dataset(filter_group, "n_bins", this % n_bins)
-    call write_dataset(filter_group, "bins", cells(this % cell) % id)
+    call write_dataset(filter_group, "bins", cells(this % cell) % id())
   end subroutine to_statepoint_distribcell
 
   subroutine initialize_distribcell(this)
@@ -102,7 +101,7 @@ contains
     val = cell_dict % get(id)
     if (val /= EMPTY) then
       this % cell = val
-      this % n_bins = cells(this % cell) % instances
+      this % n_bins = cells(this % cell) % n_instances()
     else
       call fatal_error("Could not find cell " // trim(to_str(id)) &
            &// " specified on tally filter.")
@@ -114,13 +113,8 @@ contains
     integer,                  intent(in) :: bin
     character(MAX_LINE_LEN)              :: label
 
-    integer                 :: offset
-    type(Universe), pointer :: univ
-
-    univ => universes(root_universe)
-    offset = 0
     label = ''
-    call find_offset(this % cell, univ, bin-1, offset, label)
+    call find_offset(this % cell, bin-1, label)
     label = "Distributed Cell " // label
   end function text_label_distribcell
 
@@ -130,279 +124,46 @@ contains
 ! the target cell with the given offset
 !===============================================================================
 
-  recursive subroutine find_offset(i_cell, univ, target_offset, offset, path)
-
+  subroutine find_offset(i_cell, target_offset, path)
     integer, intent(in) :: i_cell         ! The target cell index
-    type(Universe), intent(in) :: univ  ! Universe to begin search
-    integer, intent(in) :: target_offset        ! Target offset
-    integer, intent(inout) :: offset    ! Current offset
-    character(*), intent(inout) :: path ! Path to offset
+    integer, intent(in) :: target_offset  ! Target offset
+    character(*), intent(inout) :: path   ! Path to offset
 
     integer :: map                  ! Index in maps vector
-    integer :: i, j                 ! Index over cells
-    integer :: k, l, m              ! Indices in lattice
-    integer :: old_k, old_l, old_m  ! Previous indices in lattice
-    integer :: n_x, n_y, n_z        ! Lattice cell array dimensions
-    integer :: n                    ! Number of cells to search
-    integer :: cell_index           ! Index in cells array
-    integer :: lat_offset           ! Offset from lattice
-    integer :: temp_offset          ! Looped sum of offsets
-    integer :: i_univ               ! index in universes array
-    logical :: this_cell = .false.  ! Advance in this cell?
-    logical :: later_cell = .false. ! Fill cells after this one?
-    type(Cell), pointer :: c           ! Pointer to current cell
-    type(Universe), pointer :: next_univ  ! Next universe to loop through
-    class(Lattice), pointer :: lat        ! Pointer to current lattice
+    integer :: i                    ! Index over cells
+
+    integer(C_INT) :: path_len
+    character(kind=C_CHAR), allocatable, target :: path_c(:)
+
+    interface
+      function distribcell_path_len(target_cell, map, target_offset, root_univ)&
+           bind(C) result(len)
+        import C_INT32_T, C_INT
+        integer(C_INT32_T), intent(in), value :: target_cell, map, &
+                                                 target_offset, root_univ
+        integer(C_INT) :: len
+      end function distribcell_path_len
+
+      subroutine distribcell_path(target_cell, map, target_offset, root_univ, &
+                                  path) bind(C)
+        import C_INT32_T, C_CHAR
+        integer(C_INT32_T), intent(in), value :: target_cell, map, &
+                                                 target_offset, root_univ
+        character(kind=C_CHAR), intent(out)   :: path(*)
+      end subroutine distribcell_path
+    end interface
 
     ! Get the distribcell index for this cell
     map = cells(i_cell) % distribcell_index
 
-    n = size(univ % cells)
-
-    ! Write to the geometry stack
-    i_univ = universe_dict % get(univ % id)
-    if (i_univ == root_universe) then
-      path = trim(path) // "u" // to_str(univ%id)
-    else
-      path = trim(path) // "->u" // to_str(univ%id)
-    end if
-
-    ! Look through all cells in this universe
-    do i = 1, n
-      ! If the cell matches the goal and the offset matches final, write to the
-      ! geometry stack
-      if (univ % cells(i) == i_cell .and. offset == target_offset) then
-        c => cells(univ % cells(i))
-        path = trim(path) // "->c" // to_str(c % id)
-        return
-      end if
-    end do
-
-    ! Find the fill cell or lattice cell that we need to enter
-    do i = 1, n
-
-      later_cell = .false.
-
-      c => cells(univ % cells(i))
-
-      this_cell = .false.
-
-      ! If we got here, we still think the target is in this universe
-      ! or further down, but it's not this exact cell.
-      ! Compare offset to next cell to see if we should enter this cell
-      if (i /= n) then
-
-        do j = i+1, n
-
-          c => cells(univ % cells(j))
-
-          ! Skip normal cells which do not have offsets
-          if (c % type == FILL_MATERIAL) cycle
-
-          ! Break loop once we've found the next cell with an offset
-          exit
-        end do
-
-        ! Ensure we didn't just end the loop by iteration
-        if (c % type /= FILL_MATERIAL) then
-
-          ! There are more cells in this universe that it could be in
-          later_cell = .true.
-
-          ! Two cases, lattice or fill cell
-          if (c % type == FILL_UNIVERSE) then
-            temp_offset = c % offset(map)
-
-          ! Get the offset of the first lattice location
-          else
-            lat => lattices(c % fill) % obj
-            temp_offset = lat % offset(map, 1, 1, 1)
-          end if
-
-          ! If the final offset is in the range of offset - temp_offset+offset
-          ! then the goal is in this cell
-          if (target_offset < temp_offset + offset) then
-            this_cell = .true.
-          end if
-        end if
-      end if
-
-      if (n == 1 .and. c % type /= FILL_MATERIAL) then
-        this_cell = .true.
-      end if
-
-      if (.not. later_cell) then
-        this_cell = .true.
-      end if
-
-      ! Get pointer to THIS cell because target must be in this cell
-      if (this_cell) then
-
-        cell_index = univ % cells(i)
-        c => cells(cell_index)
-
-        path = trim(path) // "->c" // to_str(c%id)
-
-        ! ====================================================================
-        ! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL
-        if (c % type == FILL_UNIVERSE) then
-
-          ! Enter this cell to update the current offset
-          offset = c % offset(map) + offset
-
-          next_univ => universes(c % fill)
-          call find_offset(i_cell, next_univ, target_offset, offset, path)
-          return
-
-        ! ====================================================================
-        ! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL
-        elseif (c % type == FILL_LATTICE) then
-
-          ! Set current lattice
-          lat => lattices(c % fill) % obj
-
-          select type (lat)
-
-          ! ==================================================================
-          ! RECTANGULAR LATTICES
-          type is (RectLattice)
-
-            ! Write to the geometry stack
-            path = trim(path) // "->l" // to_str(lat%id)
-
-            n_x = lat % n_cells(1)
-            n_y = lat % n_cells(2)
-            n_z = lat % n_cells(3)
-            old_m = 1
-            old_l = 1
-            old_k = 1
-
-            ! Loop over lattice coordinates
-            do m = 1, n_z
-              do l = 1, n_y
-                do k = 1, n_x
-
-                  if (target_offset >= lat % offset(map, k, l, m) + offset) then
-                    if (k == n_x .and. l == n_y .and. m == n_z) then
-                      ! This is last lattice cell, so target must be here
-                      lat_offset = lat % offset(map, k, l, m)
-                      offset = offset + lat_offset
-                      next_univ => universes(lat % universes(k, l, m))
-                      if (lat % is_3d) then
-                        path = trim(path) // "(" // trim(to_str(k-1)) // &
-                             "," // trim(to_str(l-1)) // "," // &
-                             trim(to_str(m-1)) // ")"
-                      else
-                        path = trim(path) // "(" // trim(to_str(k-1)) // &
-                             "," // trim(to_str(l-1)) // ")"
-                      end if
-                      call find_offset(i_cell, next_univ, target_offset, offset, path)
-                      return
-                    else
-                      old_m = m
-                      old_l = l
-                      old_k = k
-                      cycle
-                    end if
-                  else
-                    ! Target is at this lattice position
-                    lat_offset = lat % offset(map, old_k, old_l, old_m)
-                    offset = offset + lat_offset
-                    next_univ => universes(lat % universes(old_k, old_l, old_m))
-                    if (lat % is_3d) then
-                      path = trim(path) // "(" // trim(to_str(old_k-1)) // &
-                           "," // trim(to_str(old_l-1)) // "," // &
-                           trim(to_str(old_m-1)) // ")"
-                    else
-                      path = trim(path) // "(" // trim(to_str(old_k-1)) // &
-                           "," // trim(to_str(old_l-1)) // ")"
-                    end if
-                    call find_offset(i_cell, next_univ, target_offset, offset, path)
-                    return
-                  end if
-
-                end do
-              end do
-            end do
-
-          ! ==================================================================
-          ! HEXAGONAL LATTICES
-          type is (HexLattice)
-
-            ! Write to the geometry stack
-            path = trim(path) // "->l" // to_str(lat%id)
-
-            n_z = lat % n_axial
-            n_y = 2 * lat % n_rings - 1
-            n_x = 2 * lat % n_rings - 1
-            old_m = 1
-            old_l = 1
-            old_k = 1
-
-            ! Loop over lattice coordinates
-            do m = 1, n_z
-              do l = 1, n_y
-                do k = 1, n_x
-
-                  ! This array position is never used
-                  if (k + l < lat % n_rings + 1) then
-                    cycle
-                  ! This array position is never used
-                  else if (k + l > 3*lat % n_rings - 1) then
-                    cycle
-                  end if
-
-                  if (target_offset >= lat % offset(map, k, l, m) + offset) then
-                    if (k == lat % n_rings .and. l == n_y .and. m == n_z) then
-                      ! This is last lattice cell, so target must be here
-                      lat_offset = lat % offset(map, k, l, m)
-                      offset = offset + lat_offset
-                      next_univ => universes(lat % universes(k, l, m))
-                      if (lat % is_3d) then
-                        path = trim(path) // "(" // &
-                             trim(to_str(k - lat % n_rings)) // "," // &
-                             trim(to_str(l - lat % n_rings)) // "," // &
-                             trim(to_str(m - 1)) // ")"
-                      else
-                        path = trim(path) // "(" // &
-                             trim(to_str(k - lat % n_rings)) // "," // &
-                             trim(to_str(l - lat % n_rings)) // ")"
-                      end if
-                      call find_offset(i_cell, next_univ, target_offset, offset, path)
-                      return
-                    else
-                      old_m = m
-                      old_l = l
-                      old_k = k
-                      cycle
-                    end if
-                  else
-                    ! Target is at this lattice position
-                    lat_offset = lat % offset(map, old_k, old_l, old_m)
-                    offset = offset + lat_offset
-                    next_univ => universes(lat % universes(old_k, old_l, old_m))
-                    if (lat % is_3d) then
-                      path = trim(path) // "(" // &
-                           trim(to_str(old_k - lat % n_rings)) // "," // &
-                           trim(to_str(old_l - lat % n_rings)) // "," // &
-                           trim(to_str(old_m - 1)) // ")"
-                    else
-                      path = trim(path) // "(" // &
-                           trim(to_str(old_k - lat % n_rings)) // "," // &
-                           trim(to_str(old_l - lat % n_rings)) // ")"
-                    end if
-                    call find_offset(i_cell, next_univ, target_offset, offset, path)
-                    return
-                  end if
-
-                end do
-              end do
-            end do
-
-          end select
-
-        end if
-      end if
+    path_len = distribcell_path_len(i_cell-1, map-1, target_offset, &
+                                    root_universe-1)
+    allocate(path_c(path_len))
+    call distribcell_path(i_cell-1, map-1, target_offset, root_universe-1, &
+                          path_c)
+    do i = 1, min(path_len, MAX_LINE_LEN)
+      if (path_c(i) == C_NULL_CHAR) exit
+      path(i:i) = path_c(i)
     end do
   end subroutine find_offset
 

src/tracking.F90

@@ -174,7 +174,7 @@ contains
         p % coord(p % n_coord) % cell = NONE
         if (any(lattice_translation /= 0)) then
           ! Particle crosses lattice boundary
-          p % surface = NONE
+          p % surface = ERROR_INT
           call cross_lattice(p, lattice_translation)
           p % event = EVENT_LATTICE
         else
@@ -205,7 +205,7 @@ contains
              call score_surface_tally(p, active_meshsurf_tallies)
 
         ! Clear surface component
-        p % surface = NONE
+        p % surface = ERROR_INT
 
         if (run_CE) then
           call collision(p)
@@ -475,7 +475,7 @@ contains
     ! COULDN'T FIND PARTICLE IN NEIGHBORING CELLS, SEARCH ALL CELLS
 
     ! Remove lower coordinate levels and assignment of surface
-    p % surface = NONE
+    p % surface = ERROR_INT
     p % n_coord = 1
     call find_cell(p, found)
 

src/volume_calc.F90

@@ -205,7 +205,8 @@ contains
       elseif (this % domain_type == FILTER_CELL) THEN
         do level = 1, p % n_coord
           do i_domain = 1, size(this % domain_id)
-            if (cells(p % coord(level) % cell) % id == this % domain_id(i_domain)) then
+            if (cells(p % coord(level) % cell) % id() &
+                 == this % domain_id(i_domain)) then
               i_material = p % material
               call check_hit(i_domain, i_material, indices, hits, n_mat)
             end if

src/xml_interface.cpp

@@ -2,9 +2,7 @@
 
 #include <algorithm>  // for std::transform
 #include <sstream>
-#include <string>
 
-#include "pugixml/pugixml.hpp"
 #include "error.h"
 
 

src/xml_interface.h

@@ -2,6 +2,7 @@
 #define XML_INTERFACE_H
 
 #include <string>
+#include <vector>
 
 #include "pugixml/pugixml.hpp"