Merge pull request #1873 from aprilnovak/mesh-multiplication

Allow constant scaling factor to be applied to libMesh unstructured mesh

docs/source/io_formats/tallies.rst

@@ -318,6 +318,10 @@ attributes/sub-elements:
   :dimension:
     The number of mesh cells in each direction. (For regular mesh only.)
 
+  :length_multiplier:
+    A multiplicative factor to apply to the mesh coordinates in all directions.
+    (For unstructured mesh only.)
+
   :lower_left:
     The lower-left corner of the structured mesh. If only two coordinates are
     given, it is assumed that the mesh is an x-y mesh. (For regular mesh only.)

include/openmc/mesh.h

@@ -331,6 +331,15 @@ public:
     true}; //!< Write tallies onto the unstructured mesh at the end of a run
   std::string filename_; //!< Path to unstructured mesh file
 
+protected:
+  //! Set the length multiplier to apply to each point in the mesh
+  void set_length_multiplier(const double length_multiplier);
+
+  // Data members
+  double length_multiplier_ {
+    1.0}; //!< Constant multiplication factor to apply to mesh coordinates
+  bool specified_length_multiplier_ {false};
+
 private:
   //! Setup method for the mesh. Builds data structures,
   //! sets up element mapping, creates bounding boxes, etc.
@@ -344,7 +353,7 @@ public:
   // Constructors
   MOABMesh() = default;
   MOABMesh(pugi::xml_node);
-  MOABMesh(const std::string& filename);
+  MOABMesh(const std::string& filename, double length_multiplier = 1.0);
   MOABMesh(std::shared_ptr<moab::Interface> external_mbi);
 
   // Overridden Methods
@@ -492,7 +501,7 @@ class LibMesh : public UnstructuredMesh {
 public:
   // Constructors
   LibMesh(pugi::xml_node node);
-  LibMesh(const std::string& filename);
+  LibMesh(const std::string& filename, double length_multiplier = 1.0);
 
   // Overridden Methods
   void bins_crossed(Position r0, Position r1, const Direction& u,

openmc/mesh.py

@@ -615,8 +615,8 @@ class UnstructuredMesh(MeshBase):
         Unique identifier for the mesh
     name : str
         Name of the mesh
-    size : int
-        Number of elements in the unstructured mesh
+    length_multiplier: float
+        Constant multiplier to apply to mesh coordinates
 
     Attributes
     ----------
@@ -626,11 +626,16 @@ class UnstructuredMesh(MeshBase):
         Name of the mesh
     filename : str
         Name of the file containing the unstructured mesh
+    length_multiplier: float
+        Multiplicative factor to apply to mesh coordinates
     library : str
         Mesh library used for the unstructured mesh tally
     output : bool
         Indicates whether or not automatic tally output should
         be generated for this mesh
+    specified_length_multiplier: bool
+        Indicates whether a non-unity length multiplier has been
+        applied to this mesh
     volumes : Iterable of float
         Volumes of the unstructured mesh elements
     total_volume : float
@@ -639,13 +644,16 @@ class UnstructuredMesh(MeshBase):
         An iterable of element centroid coordinates, e.g. [(0.0, 0.0, 0.0),
         (1.0, 1.0, 1.0), ...]
     """
-    def __init__(self, filename, library, mesh_id=None, name=''):
+    def __init__(self, filename, library, mesh_id=None, name='',
+                        length_multiplier=1.0):
         super().__init__(mesh_id, name)
         self.filename = filename
         self._volumes = None
         self._centroids = None
         self.library = library
         self._output = True
+        self._specified_length_multiplier = False
+        self.length_multiplier = length_multiplier
 
     @property
     def filename(self):
@@ -713,6 +721,20 @@ class UnstructuredMesh(MeshBase):
                       Iterable, Real)
         self._centroids = centroids
 
+    @property
+    def length_multiplier(self):
+        return self._length_multiplier
+
+    @length_multiplier.setter
+    def length_multiplier(self, length_multiplier):
+        cv.check_type("Unstructured mesh length multiplier", 
+                                length_multiplier,
+                                Real)
+        self._length_multiplier = length_multiplier
+
+        if (self._length_multiplier != 1.0):
+          self._specified_length_multiplier = True
+
     def __repr__(self):
         string = super().__repr__()
         string += '{: <16}=\t{}\n'.format('\tFilename', self.filename)
@@ -763,7 +785,7 @@ class UnstructuredMesh(MeshBase):
 
         for centroid in self.centroids:
             # create a point for each centroid
-            point_id = points.InsertNextPoint(centroid)
+            point_id = points.InsertNextPoint(centroid * self.length_multiplier)
             # create a cell of type "Vertex" for each point
             cell_id = vertices.InsertNextCell(cell_dim, (point_id,))
 
@@ -817,6 +839,9 @@ class UnstructuredMesh(MeshBase):
         mesh.centroids = np.reshape(centroids, (vol_data.shape[0], 3))
         mesh.size = mesh.volumes.size
 
+        if 'length_multiplier' in group:
+            mesh.length_multiplier = group['length_multiplier'][()]
+
         return mesh
 
     def to_xml_element(self):
@@ -836,6 +861,9 @@ class UnstructuredMesh(MeshBase):
         subelement = ET.SubElement(element, "filename")
         subelement.text = self.filename
 
+        if (self._specified_length_multiplier):
+          element.set("length_multiplier", str(self.length_multiplier))
+
         return element
 
     @classmethod
@@ -855,5 +883,6 @@ class UnstructuredMesh(MeshBase):
         mesh_id = int(get_text(elem, 'id'))
         filename = get_text(elem, 'filename')
         library = get_text(elem, 'library')
+        length_multiplier = float(get_text(elem, 'length_multiplier', 1.0))
 
-        return cls(filename, library, mesh_id)
+        return cls(filename, library, mesh_id, '', length_multiplier)

src/mesh.cpp

@@ -35,6 +35,7 @@
 #ifdef LIBMESH
 #include "libmesh/mesh_tools.h"
 #include "libmesh/numeric_vector.h"
+#include "libmesh/mesh_modification.h"
 #endif
 
 namespace openmc {
@@ -171,6 +172,12 @@ UnstructuredMesh::UnstructuredMesh(pugi::xml_node node) : Mesh(node)
     }
   }
 
+  // check if a length unit multiplier was specified
+  if (check_for_node(node, "length_multiplier")) {
+    length_multiplier_ = std::stod(get_node_value(node, "length_multiplier"));
+    specified_length_multiplier_ = true;
+  }
+
   // get the filename of the unstructured mesh to load
   if (check_for_node(node, "filename")) {
     filename_ = get_node_value(node, "filename");
@@ -218,9 +225,21 @@ void UnstructuredMesh::to_hdf5(hid_t group) const
 
   write_dataset(mesh_group, "volumes", tet_vols);
   write_dataset(mesh_group, "centroids", centroids);
+
+  if (specified_length_multiplier_)
+    write_dataset(mesh_group, "length_multiplier", length_multiplier_);
+
   close_group(mesh_group);
 }
 
+void UnstructuredMesh::set_length_multiplier(double length_multiplier)
+{
+  length_multiplier_ = length_multiplier;
+
+  if (length_multiplier_ != 1.0)
+    specified_length_multiplier_ = true;
+}
+
 void StructuredMesh::get_indices(Position r, int* ijk, bool* in_mesh) const
 {
   *in_mesh = true;
@@ -1585,9 +1604,10 @@ MOABMesh::MOABMesh(pugi::xml_node node) : UnstructuredMesh(node)
   initialize();
 }
 
-MOABMesh::MOABMesh(const std::string& filename)
+MOABMesh::MOABMesh(const std::string& filename, double length_multiplier)
 {
   filename_ = filename;
+  set_length_multiplier(length_multiplier);
   initialize();
 }
 
@@ -1634,6 +1654,34 @@ void MOABMesh::initialize()
     fatal_error("Failed to add tetrahedra to an entity set.");
   }
 
+  if (specified_length_multiplier_) {
+    // get the connectivity of all tets
+    moab::Range adj;
+    rval = mbi_->get_adjacencies(ehs_, 0, true, adj, moab::Interface::UNION);
+    if (rval != moab::MB_SUCCESS) {
+      fatal_error("Failed to get adjacent vertices of tetrahedra.");
+    }
+    // scale all vertex coords by multiplier (done individually so not all
+    // coordinates are in memory twice at once)
+    for (auto vert : adj) {
+      // retrieve coords
+      std::array<double, 3> coord;
+      rval = mbi_->get_coords(&vert, 1, coord.data());
+      if (rval != moab::MB_SUCCESS) {
+        fatal_error("Could not get coordinates of vertex.");
+      }
+      // scale coords
+      for (auto& c : coord) {
+        c *= length_multiplier_;
+      }
+      // set new coords
+      rval = mbi_->set_coords(&vert, 1, coord.data());
+      if (rval != moab::MB_SUCCESS) {
+        fatal_error("Failed to set new vertex coordinates");
+      }
+    }
+  }
+
   // build acceleration data structures
   compute_barycentric_data(ehs_);
   build_kdtree(ehs_);
@@ -2144,9 +2192,10 @@ LibMesh::LibMesh(pugi::xml_node node) : UnstructuredMesh(node)
   initialize();
 }
 
-LibMesh::LibMesh(const std::string& filename)
+LibMesh::LibMesh(const std::string& filename, double length_multiplier)
 {
   filename_ = filename;
+  set_length_multiplier(length_multiplier);
   initialize();
 }
 
@@ -2162,6 +2211,11 @@ void LibMesh::initialize()
 
   m_ = make_unique<libMesh::Mesh>(*settings::libmesh_comm, n_dimension_);
   m_->read(filename_);
+
+  if (specified_length_multiplier_) {
+    libMesh::MeshTools::Modification::scale(*m_, length_multiplier_);
+  }
+
   m_->prepare_for_use();
 
   // ensure that the loaded mesh is 3 dimensional