Add skeleton of C++ lattices

CMakeLists.txt

@@ -439,6 +439,8 @@ set(LIBOPENMC_CXX_SRC
   src/cell.h
   src/error.h
   src/hdf5_interface.h
+  src/lattice.cpp
+  src/lattice.h
   src/random_lcg.cpp
   src/random_lcg.h
   src/surface.cpp

src/cell.cpp

@@ -301,23 +301,25 @@ Cell::to_hdf5(hid_t cell_group) const
   //write_int(cell_group, "universe", universe);
 
   // Write the region specification.
-  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);
+  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());
   }
-  write_string(cell_group, "region", region_spec.str());
 
 //  close_group(cell_group);
 }

src/cell.h

@@ -25,6 +25,19 @@ std::vector<Cell> cells_c;
 
 std::map<int, int> cell_dict;
 
+//==============================================================================
+//! A geometry primitive that fills all space and contains cells.
+//==============================================================================
+
+class Universe
+{
+  public:
+  int32_t id;                  //! Unique ID
+  int32_t type;
+  std::vector<int32_t> cells;  //! Cells within this universe
+  double x0, y0, z0;           //! Translation coordinates.
+};
+
 //==============================================================================
 //! A geometry primitive that links surfaces, universes, and materials
 //==============================================================================

src/geometry.F90

@@ -260,7 +260,7 @@ contains
               ! Particle is outside the lattice.
               if (lat % outer == NO_OUTER_UNIVERSE) then
                 call p % mark_as_lost("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.")
                 return
               else
@@ -299,7 +299,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)) // ")")

src/geometry_header.F90

@@ -79,6 +79,29 @@ module geometry_header
       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)
+      use ISO_C_BINDING
+      implicit none
+      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)
+      use ISO_C_BINDING
+      implicit none
+      type(C_PTR), intent(in), value :: lat_ptr
+      integer(C_INT32_T)             :: id
+    end function lattice_id_c
+
+    subroutine lattice_to_hdf5_c(lat_ptr, group) bind(C, name='lattice_to_hdf5')
+      use ISO_C_BINDING
+      use hdf5
+      implicit none
+      type(C_PTR),    intent(in), value :: lat_ptr
+      integer(HID_T), intent(in), value :: group
+    end subroutine lattice_to_hdf5_c
   end interface
 
 !===============================================================================
@@ -99,8 +122,9 @@ module geometry_header
 !===============================================================================
 
   type, abstract :: Lattice
-    integer              :: id               ! Universe number for lattice
-    character(len=104)   :: name = ""        ! User-defined name
+    type(C_PTR) :: ptr
+
+    !integer              :: id               ! Universe number for lattice
     real(8), allocatable :: pitch(:)         ! Pitch along each axis
     integer, allocatable :: universes(:,:,:) ! Specified universes
     integer              :: outside          ! Material to fill area outside
@@ -108,6 +132,10 @@ module geometry_header
     logical              :: is_3d            ! Lattice has cells on z axis
     integer, allocatable :: offset(:,:,:,:)  ! Distribcell offsets
   contains
+
+    procedure :: id => lattice_id
+    procedure :: to_hdf5 => lattice_to_hdf5
+
     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
@@ -247,6 +275,18 @@ 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
+
+  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_rect(this, i_xyz) result(is_valid)

src/input_xml.F90

@@ -60,6 +60,12 @@ module input_xml
       implicit none
       type(C_PTR) :: node_ptr
     end subroutine read_cells
+
+    subroutine read_lattices(node_ptr) bind(C)
+      use ISO_C_BINDING
+      implicit none
+      type(C_PTR) :: node_ptr
+    end subroutine read_lattices
   end interface
 
 contains
@@ -1240,6 +1246,8 @@ contains
     ! ==========================================================================
     ! READ LATTICES FROM GEOMETRY.XML
 
+    call read_lattices(root % ptr)
+
     ! Get pointer to list of XML <lattice>
     call get_node_list(root, "lattice", node_rlat_list)
     call get_node_list(root, "hex_lattice", node_hlat_list)
@@ -1253,30 +1261,13 @@ contains
     RECT_LATTICES: do i = 1, n_rlats
       allocate(RectLattice::lattices(i) % obj)
       lat => lattices(i) % obj
+      lat % ptr = lattice_pointer_c(i - 1)
       select type(lat)
       type is (RectLattice)
 
       ! Get pointer to i-th lattice
       node_lat = node_rlat_list(i)
 
-      ! ID of lattice
-      if (check_for_node(node_lat, "id")) then
-        call get_node_value(node_lat, "id", lat % id)
-      else
-        call fatal_error("Must specify id of lattice in geometry XML file.")
-      end if
-
-      ! Check to make sure 'id' hasn't been used
-      if (lattice_dict % has(lat % id)) then
-        call fatal_error("Two or more lattices use the same unique ID: " &
-             // to_str(lat % id))
-      end if
-
-      ! Copy lattice name
-      if (check_for_node(node_lat, "name")) then
-        call get_node_value(node_lat, "name", lat % name)
-      end if
-
       ! Read number of lattice cells in each dimension
       n = node_word_count(node_lat, "dimension")
       if (n == 2) then
@@ -1341,7 +1332,7 @@ contains
       n = node_word_count(node_lat, "universes")
       if (n /= n_x*n_y*n_z) then
         call fatal_error("Number of universes on <universes> does not match &
-             &size of lattice " // trim(to_str(lat % id)) // ".")
+             &size of lattice " // trim(to_str(lat % id())) // ".")
       end if
 
       allocate(temp_int_array(n))
@@ -1377,7 +1368,7 @@ contains
       end if
 
       ! Add lattice to dictionary
-      call lattice_dict % set(lat % id, i)
+      call lattice_dict % set(lat % id(), i)
 
       end select
     end do RECT_LATTICES
@@ -1385,30 +1376,13 @@ contains
     HEX_LATTICES: do i = 1, n_hlats
       allocate(HexLattice::lattices(n_rlats + i) % obj)
       lat => lattices(n_rlats + i) % obj
+      lat % ptr = lattice_pointer_c(n_rlats + i - 1)
       select type (lat)
       type is (HexLattice)
 
       ! Get pointer to i-th lattice
       node_lat = node_hlat_list(i)
 
-      ! ID of lattice
-      if (check_for_node(node_lat, "id")) then
-        call get_node_value(node_lat, "id", lat % id)
-      else
-        call fatal_error("Must specify id of lattice in geometry XML file.")
-      end if
-
-      ! Check to make sure 'id' hasn't been used
-      if (lattice_dict % has(lat % id)) then
-        call fatal_error("Two or more lattices use the same unique ID: " &
-             // to_str(lat % id))
-      end if
-
-      ! Copy lattice name
-      if (check_for_node(node_lat, "name")) then
-        call get_node_value(node_lat, "name", lat % name)
-      end if
-
       ! Read number of lattice cells in each dimension
       call get_node_value(node_lat, "n_rings", lat % n_rings)
       if (check_for_node(node_lat, "n_axial")) then
@@ -1454,7 +1428,7 @@ contains
       n = node_word_count(node_lat, "universes")
       if (n /= (3*n_rings**2 - 3*n_rings + 1)*n_z) then
         call fatal_error("Number of universes on <universes> does not match &
-             &size of lattice " // trim(to_str(lat % id)) // ".")
+             &size of lattice " // trim(to_str(lat % id())) // ".")
       end if
 
       allocate(temp_int_array(n))
@@ -1564,7 +1538,7 @@ contains
       end if
 
       ! Add lattice to dictionary
-      call lattice_dict % set(lat % id, n_rlats + i)
+      call lattice_dict % set(lat % id(), n_rlats + i)
 
       end select
     end do HEX_LATTICES
@@ -4353,7 +4327,7 @@ contains
               else
                 call fatal_error("Invalid universe number " &
                      &// trim(to_str(id)) // " specified on lattice " &
-                     &// trim(to_str(lat % id)))
+                     &// trim(to_str(lat % id())))
               end if
             end do
           end do
@@ -4374,7 +4348,7 @@ contains
               else
                 call fatal_error("Invalid universe number " &
                      &// trim(to_str(id)) // " specified on lattice " &
-                     &// trim(to_str(lat % id)))
+                     &// trim(to_str(lat % id())))
               end if
             end do
           end do
@@ -4388,7 +4362,7 @@ contains
         else
           call fatal_error("Invalid universe number " &
                &// trim(to_str(lat % outer)) &
-               &// " specified on lattice " // trim(to_str(lat % id)))
+               &// " specified on lattice " // trim(to_str(lat % id())))
         end if
       end if
 

src/lattice.cpp

@@ -0,0 +1,109 @@
+#include "lattice.h"
+
+#include <sstream>
+#include <vector>
+
+#include "error.h"
+#include "hdf5_interface.h"
+#include "xml_interface.h"
+
+//TODO: remove this include
+#include <iostream>
+
+
+namespace openmc {
+
+//==============================================================================
+// Global variables
+//==============================================================================
+
+// Braces force n_lattices to be defined here, not just declared.
+//extern "C" {int32_t n_lattices {0};}
+
+//Lattice **lattices_c;
+std::vector<Lattice*> lattices_c;
+
+std::map<int32_t, int32_t> lattice_dict;
+
+//==============================================================================
+// 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");
+  }
+}
+
+void
+Lattice::to_hdf5(hid_t lat_group) const
+{
+  if (!name.empty()) {
+    write_string(lat_group, "name", name);
+  }
+}
+
+//==============================================================================
+// RectLattice implementation
+//==============================================================================
+
+RectLattice::RectLattice(pugi::xml_node lat_node)
+  : Lattice {lat_node}
+{
+}
+
+//==============================================================================
+// HexLattice implementation
+//==============================================================================
+
+HexLattice::HexLattice(pugi::xml_node lat_node)
+  : Lattice {lat_node}
+{
+}
+
+//==============================================================================
+
+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 dictionary.
+  for (int i_lat = 0; i_lat < lattices_c.size(); i_lat++) {
+    int id = lattices_c[i_lat]->id;
+    auto in_dict = lattice_dict.find(id);
+    if (in_dict == lattice_dict.end()) {
+      lattice_dict[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;}
+
+  void lattice_to_hdf5(Lattice *lat, hid_t group) {lat->to_hdf5(group);}
+}
+
+} // namespace openmc

src/lattice.h

@@ -0,0 +1,77 @@
+#ifndef LATTICE_H
+#define LATTICE_H
+
+#include <cstdint>
+#include <map>
+#include <string>
+#include <vector>
+
+#include "hdf5.h"
+#include "pugixml/pugixml.hpp"
+
+
+namespace openmc {
+
+//==============================================================================
+// Global variables
+//==============================================================================
+
+//extern "C" int32_t n_lattice;
+
+class Lattice;
+//extern Lattice **lattices_c;
+extern std::vector<Lattice*> lattices_c;
+
+extern std::map<int32_t, int32_t> lattice_dict;
+
+//==============================================================================
+//! Abstract type for ordered array of universes
+//==============================================================================
+
+class Lattice
+{
+public:
+  int32_t id;                     //! Universe ID number
+  std::string name;               //! User-defined name
+  //std::vector<double> pitch;      //! Pitch along each basis
+  //std::vector<int32_t> universes; //! Universes filling each lattice tile
+  //int32_t outer;                  //! Universe tiled outside the lattice
+  //std::vector<int32_t> offset;    //! Distribcell offsets
+
+  explicit Lattice(pugi::xml_node lat_node);
+
+  virtual ~Lattice() {}
+
+  //virtual bool are_valid_indices(const int i_xyz[3]) const = 0;
+
+  //virtual void get_indices(const double global_xyz[3], int i_xyz[3]) const = 0;
+
+  //virtual void get_local_xyz(const double global_xyz[3], const int i_xyz[3],
+  //                           double local_xyz[3]) const = 0;
+
+  //! Write all information needed to reconstruct the lattice to an HDF5 group.
+  //! @param group_id An HDF5 group id.
+  void to_hdf5(hid_t group_id) const;
+};
+
+//==============================================================================
+//==============================================================================
+
+class RectLattice : public Lattice
+{
+public:
+  explicit RectLattice(pugi::xml_node lat_node);
+
+  virtual ~RectLattice() {}
+};
+
+class HexLattice : public Lattice
+{
+public:
+  explicit HexLattice(pugi::xml_node lat_node);
+
+  virtual ~HexLattice() {}
+};
+
+} //  namespace openmc
+#endif // LATTICE_H

src/output.F90

@@ -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)) // ')'

src/summary.F90

@@ -200,7 +200,7 @@ contains
 
       case (FILL_LATTICE)
         call write_dataset(cell_group, "fill_type", "lattice")
-        call write_dataset(cell_group, "lattice", lattices(c%fill)%obj%id)
+        call write_dataset(cell_group, "lattice", lattices(c%fill)%obj%id())
       end select
 
       call close_group(cell_group)
@@ -258,10 +258,11 @@ contains
     ! Write information on each lattice
     LATTICE_LOOP: do i = 1, n_lattices
       lat => lattices(i)%obj
-      lattice_group = create_group(lattices_group, "lattice " // trim(to_str(lat%id)))
+      lattice_group = create_group(lattices_group, "lattice " // trim(to_str(lat%id())))
+
+      call lat % to_hdf5(lattice_group)
 
       ! 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)

src/tallies/tally_filter_distribcell.F90

@@ -270,7 +270,7 @@ contains
           type is (RectLattice)
 
             ! Write to the geometry stack
-            path = trim(path) // "->l" // to_str(lat%id)
+            path = trim(path) // "->l" // to_str(lat%id())
 
             n_x = lat % n_cells(1)
             n_y = lat % n_cells(2)
@@ -332,7 +332,7 @@ contains
           type is (HexLattice)
 
             ! Write to the geometry stack
-            path = trim(path) // "->l" // to_str(lat%id)
+            path = trim(path) // "->l" // to_str(lat%id())
 
             n_z = lat % n_axial
             n_y = 2 * lat % n_rings - 1