load mesh objects from weight_windows.h5 file

include/openmc/mesh.h

@@ -132,6 +132,7 @@ class Mesh {
   // Constructors and destructor
   Mesh() = default;
   Mesh(pugi::xml_node node);
+  Mesh(hid_t group);
   virtual ~Mesh() = default;
 
   // Methods
@@ -258,6 +259,7 @@ class StructuredMesh : public Mesh {
 public:
   StructuredMesh() = default;
   StructuredMesh(pugi::xml_node node) : Mesh {node} {};
+  StructuredMesh(hid_t group) : Mesh {group} {};
   virtual ~StructuredMesh() = default;
 
   using MeshIndex = std::array<int, 3>;
@@ -423,6 +425,7 @@ class PeriodicStructuredMesh : public StructuredMesh {
 public:
   PeriodicStructuredMesh() = default;
   PeriodicStructuredMesh(pugi::xml_node node) : StructuredMesh {node} {};
+  PeriodicStructuredMesh(hid_t group) : StructuredMesh {group} {};
 
   Position local_coords(const Position& r) const override
   {
@@ -442,6 +445,7 @@ class RegularMesh : public StructuredMesh {
   // Constructors
   RegularMesh() = default;
   RegularMesh(pugi::xml_node node);
+  RegularMesh(hid_t group);
 
   // Overridden methods
   int get_index_in_direction(double r, int i) const override;
@@ -492,6 +496,7 @@ class RectilinearMesh : public StructuredMesh {
   // Constructors
   RectilinearMesh() = default;
   RectilinearMesh(pugi::xml_node node);
+  RectilinearMesh(hid_t group);
 
   // Overridden methods
   int get_index_in_direction(double r, int i) const override;
@@ -534,6 +539,7 @@ class CylindricalMesh : public PeriodicStructuredMesh {
   // Constructors
   CylindricalMesh() = default;
   CylindricalMesh(pugi::xml_node node);
+  CylindricalMesh(hid_t group);
 
   // Overridden methods
   virtual MeshIndex get_indices(Position r, bool& in_mesh) const override;
@@ -598,6 +604,7 @@ class SphericalMesh : public PeriodicStructuredMesh {
   // Constructors
   SphericalMesh() = default;
   SphericalMesh(pugi::xml_node node);
+  SphericalMesh(hid_t group);
 
   // Overridden methods
   virtual MeshIndex get_indices(Position r, bool& in_mesh) const override;
@@ -668,6 +675,7 @@ class UnstructuredMesh : public Mesh {
   // Constructors
   UnstructuredMesh() { n_dimension_ = 3; };
   UnstructuredMesh(pugi::xml_node node);
+  UnstructuredMesh(hid_t group);
 
   static const std::string mesh_type;
   virtual std::string get_mesh_type() const override;
@@ -774,6 +782,7 @@ class MOABMesh : public UnstructuredMesh {
   // Constructors
   MOABMesh() = default;
   MOABMesh(pugi::xml_node);
+  MOABMesh(hid_t group);
   MOABMesh(const std::string& filename, double length_multiplier = 1.0);
   MOABMesh(std::shared_ptr<moab::Interface> external_mbi);
 
@@ -943,6 +952,7 @@ class LibMesh : public UnstructuredMesh {
 public:
   // Constructors
   LibMesh(pugi::xml_node node);
+  LibMesh(hid_t group);
   LibMesh(const std::string& filename, double length_multiplier = 1.0);
   LibMesh(libMesh::MeshBase& input_mesh, double length_multiplier = 1.0);
 
@@ -1069,6 +1079,11 @@ class AdaptiveLibMesh : public LibMesh {
 //! \param[in] root XML node
 void read_meshes(pugi::xml_node root);
 
+//! Read meshes from a HDF5 file
+//
+//! \param[in] group HDF5 group
+void read_meshes(hid_t group);
+
 //! Write mesh data to an HDF5 group
 //
 //! \param[in] group HDF5 group

src/mesh.cpp

@@ -238,6 +238,17 @@ Mesh::Mesh(pugi::xml_node node)
     name_ = get_node_value(node, "name");
 }
 
+Mesh::Mesh(hid_t group)
+{
+  // Read mesh ID
+  read_attribute(group, "id", id_);
+
+  // Read mesh name
+  if (object_exists(group, "name")) {
+    read_dataset(group, "name", name_);
+  }
+}
+
 void Mesh::set_id(int32_t id)
 {
   assert(id >= 0 || id == C_NONE);
@@ -627,6 +638,46 @@ UnstructuredMesh::UnstructuredMesh(pugi::xml_node node) : Mesh(node)
   }
 }
 
+UnstructuredMesh::UnstructuredMesh(hid_t group) : Mesh(group)
+{
+  n_dimension_ = 3;
+
+  // check the mesh type
+  if (object_exists(group, "type")) {
+    std::string temp;
+    read_dataset(group, "type", temp);
+    if (temp != mesh_type) {
+      fatal_error(fmt::format("Invalid mesh type: {}", temp));
+    }
+  }
+
+  // check if a length unit multiplier was specified
+  if (object_exists(group, "length_multiplier")) {
+    read_dataset(group, "length_multiplier", length_multiplier_);
+  }
+
+  // get the filename of the unstructured mesh to load
+  if (object_exists(group, "filename")) {
+    read_dataset(group, "filename", filename_);
+    if (!file_exists(filename_)) {
+      fatal_error("Mesh file '" + filename_ + "' does not exist!");
+    }
+  } else {
+    fatal_error(fmt::format(
+      "No filename supplied for unstructured mesh with ID: {}", id_));
+  }
+
+  if (attribute_exists(group, "options")) {
+    read_attribute(group, "options", options_);
+  }
+
+  // check if mesh tally data should be written with
+  // statepoint files
+  if (attribute_exists(group, "output")) {
+    read_attribute(group, "output", output_);
+  }
+}
+
 void UnstructuredMesh::determine_bounds()
 {
   double xmin = INFTY;
@@ -1174,6 +1225,73 @@ RegularMesh::RegularMesh(pugi::xml_node node) : StructuredMesh {node}
   }
 }
 
+RegularMesh::RegularMesh(hid_t group) : StructuredMesh {group}
+{
+  // Determine number of dimensions for mesh
+  if (!object_exists(group, "dimension")) {
+    fatal_error("Must specify <dimension> on a regular mesh.");
+  }
+
+  xt::xtensor<int, 1> shape;
+  read_dataset(group, "dimension", shape);
+  int n = n_dimension_ = shape.size();
+  if (n != 1 && n != 2 && n != 3) {
+    fatal_error("Mesh must be one, two, or three dimensions.");
+  }
+  std::copy(shape.begin(), shape.end(), shape_.begin());
+
+  // Check that dimensions are all greater than zero
+  if (xt::any(shape <= 0)) {
+    fatal_error("All entries on the <dimension> element for a tally "
+                "mesh must be positive.");
+  }
+
+  // Check for lower-left coordinates
+  if (object_exists(group, "lower_left")) {
+    // Read mesh lower-left corner location
+    read_dataset(group, "lower_left", lower_left_);
+  } else {
+    fatal_error("Must specify <lower_left> on a mesh.");
+  }
+
+  // Make sure lower_left and dimension match
+  if (shape.size() != lower_left_.size()) {
+    fatal_error("Number of entries on <lower_left> must be the same "
+                "as the number of entries on <dimension>.");
+  }
+
+  if (object_exists(group, "upper_right")) {
+
+    read_dataset(group, "upper_right", upper_right_);
+
+    // Check to ensure width has same dimensions
+    auto n = upper_right_.size();
+    if (n != lower_left_.size()) {
+      fatal_error("Number of entries on <upper_right> must be the "
+                  "same as the number of entries on <lower_left>.");
+    }
+
+    // Check that upper-right is above lower-left
+    if (xt::any(upper_right_ < lower_left_)) {
+      fatal_error("The <upper_right> coordinates must be greater than "
+                  "the <lower_left> coordinates on a tally mesh.");
+    }
+
+    // Set width
+    width_ = xt::eval((upper_right_ - lower_left_) / shape);
+  } else {
+    fatal_error("Must specify either <upper_right> or <width> on a mesh.");
+  }
+
+  // Set material volumes
+  volume_frac_ = 1.0 / xt::prod(shape)();
+
+  element_volume_ = 1.0;
+  for (int i = 0; i < n_dimension_; i++) {
+    element_volume_ *= width_[i];
+  }
+}
+
 int RegularMesh::get_index_in_direction(double r, int i) const
 {
   return std::ceil((r - lower_left_[i]) / width_[i]);
@@ -1343,6 +1461,19 @@ RectilinearMesh::RectilinearMesh(pugi::xml_node node) : StructuredMesh {node}
   }
 }
 
+RectilinearMesh::RectilinearMesh(hid_t group) : StructuredMesh {group}
+{
+  n_dimension_ = 3;
+
+  read_dataset(group, "x_grid", grid_[0]);
+  read_dataset(group, "y_grid", grid_[1]);
+  read_dataset(group, "z_grid", grid_[2]);
+
+  if (int err = set_grid()) {
+    fatal_error(openmc_err_msg);
+  }
+}
+
 const std::string RectilinearMesh::mesh_type = "rectilinear";
 
 std::string RectilinearMesh::get_mesh_type() const
@@ -1478,6 +1609,19 @@ CylindricalMesh::CylindricalMesh(pugi::xml_node node)
   }
 }
 
+CylindricalMesh::CylindricalMesh(hid_t group) : PeriodicStructuredMesh {group}
+{
+  n_dimension_ = 3;
+  read_dataset(group, "r_grid", grid_[0]);
+  read_dataset(group, "phi_grid", grid_[1]);
+  read_dataset(group, "z_grid", grid_[2]);
+  read_dataset(group, "origin", origin_);
+
+  if (int err = set_grid()) {
+    fatal_error(openmc_err_msg);
+  }
+}
+
 const std::string CylindricalMesh::mesh_type = "cylindrical";
 
 std::string CylindricalMesh::get_mesh_type() const
@@ -1756,6 +1900,20 @@ SphericalMesh::SphericalMesh(pugi::xml_node node)
   }
 }
 
+SphericalMesh::SphericalMesh(hid_t group) : PeriodicStructuredMesh {group}
+{
+  n_dimension_ = 3;
+
+  read_dataset(group, "r_grid", grid_[0]);
+  read_dataset(group, "theta_grid", grid_[1]);
+  read_dataset(group, "phi_grid", grid_[2]);
+  read_dataset(group, "origin", origin_);
+
+  if (int err = set_grid()) {
+    fatal_error(openmc_err_msg);
+  }
+}
+
 const std::string SphericalMesh::mesh_type = "spherical";
 
 std::string SphericalMesh::get_mesh_type() const
@@ -2520,6 +2678,11 @@ MOABMesh::MOABMesh(pugi::xml_node node) : UnstructuredMesh(node)
   initialize();
 }
 
+MOABMesh::MOABMesh(hid_t group) : UnstructuredMesh(group)
+{
+  initialize();
+}
+
 MOABMesh::MOABMesh(const std::string& filename, double length_multiplier)
   : UnstructuredMesh()
 {
@@ -3228,6 +3391,15 @@ LibMesh::LibMesh(pugi::xml_node node) : UnstructuredMesh(node)
   initialize();
 }
 
+LibMesh::LibMesh(hid_t group) : UnstructuredMesh(group)
+{
+  // filename_ and length_multiplier_ will already be set by the
+  // UnstructuredMesh constructor
+  set_mesh_pointer_from_filename(filename_);
+  set_length_multiplier(length_multiplier_);
+  initialize();
+}
+
 // create the mesh from a pointer to a libMesh Mesh
 LibMesh::LibMesh(libMesh::MeshBase& input_mesh, double length_multiplier)
 {
@@ -3649,6 +3821,77 @@ void read_meshes(pugi::xml_node root)
   }
 }
 
+void read_meshes(hid_t group)
+{
+  std::unordered_set<int> mesh_ids;
+
+  std::array<int, 1> ids;
+  read_attribute(group, "ids", ids);
+
+  for (auto id : ids) {
+
+    // Check to make sure multiple meshes in the same file don't share IDs
+    if (contains(mesh_ids, id)) {
+      fatal_error(fmt::format("Two or more meshes use the same unique ID "
+                              "'{}' in the same HDF5 input file",
+        id));
+    }
+    mesh_ids.insert(id);
+
+    // If we've already read a mesh with the same ID in a *different* file,
+    // assume it is the same here
+    if (model::mesh_map.find(id) != model::mesh_map.end()) {
+      warning(fmt::format("Mesh with ID={} appears in multiple files.", id));
+      continue;
+    }
+
+    std::string name = fmt::format("mesh {}", id);
+    hid_t mesh_group = open_group(group, name.c_str());
+
+    std::string mesh_type;
+    if (object_exists(mesh_group, "type")) {
+      read_dataset(mesh_group, "type", mesh_type);
+    } else {
+      mesh_type = "regular";
+    }
+
+    // determine the mesh library to use
+    std::string mesh_lib;
+    if (object_exists(mesh_group, "library")) {
+      read_dataset(mesh_group, "library", mesh_lib);
+    }
+
+    // Read mesh and add to vector
+    if (mesh_type == RegularMesh::mesh_type) {
+      model::meshes.push_back(make_unique<RegularMesh>(mesh_group));
+    } else if (mesh_type == RectilinearMesh::mesh_type) {
+      model::meshes.push_back(make_unique<RectilinearMesh>(mesh_group));
+    } else if (mesh_type == CylindricalMesh::mesh_type) {
+      model::meshes.push_back(make_unique<CylindricalMesh>(mesh_group));
+    } else if (mesh_type == SphericalMesh::mesh_type) {
+      model::meshes.push_back(make_unique<SphericalMesh>(mesh_group));
+#ifdef OPENMC_DAGMC_ENABLED
+    } else if (mesh_type == UnstructuredMesh::mesh_type &&
+               mesh_lib == MOABMesh::mesh_lib_type) {
+      model::meshes.push_back(make_unique<MOABMesh>(mesh_group));
+#endif
+#ifdef OPENMC_LIBMESH_ENABLED
+    } else if (mesh_type == UnstructuredMesh::mesh_type &&
+               mesh_lib == LibMesh::mesh_lib_type) {
+      model::meshes.push_back(make_unique<LibMesh>(mesh_group));
+#endif
+    } else if (mesh_type == UnstructuredMesh::mesh_type) {
+      fatal_error("Unstructured mesh support is not enabled or the mesh "
+                  "library is invalid.");
+    } else {
+      fatal_error("Invalid mesh type: " + mesh_type);
+    }
+
+    // Map ID to position in vector
+    model::mesh_map[model::meshes.back()->id_] = model::meshes.size() - 1;
+  }
+}
+
 void meshes_to_hdf5(hid_t group)
 {
   // Write number of meshes