Factoring out a new SructuredMesh class. Plot method placement still needs to be resolved.

Updating unstructured mesh to allow use of indices only.

Going back to a single unstructure mesh class until we need the added abstraction.

include/openmc/mesh.h

@@ -69,6 +69,43 @@ public:
   //! \return Mesh bin
   virtual int get_bin(Position r) const = 0;
 
+  //! Get the number of mesh cells.
+  virtual int n_bins() const = 0;
+
+  //! Get the number of mesh cell surfaces.
+  virtual int n_surface_bins() const = 0;
+
+  //! Write mesh data to an HDF5 group
+  //
+  //! \param[in] group HDF5 group
+  virtual void to_hdf5(hid_t group) const = 0;
+
+  //! Find the mesh lines that intersect an axis-aligned slice plot
+  //
+  //! \param[in] plot_ll The lower-left coordinates of the slice plot.
+  //! \param[in] plot_ur The upper-right coordinates of the slice plot.
+  //! \return A pair of vectors indicating where the mesh lines lie along each
+  //!   of the plot's axes.  For example an xy-slice plot will get back a vector
+  //!   of x-coordinates and another of y-coordinates.  These vectors may be
+  //!   empty for low-dimensional meshes.
+  virtual std::pair<std::vector<double>, std::vector<double>>
+  plot(Position plot_ll, Position plot_ur) const = 0;
+
+  //! Get a label for the mesh bin
+  virtual std::string bin_label(int bin) const = 0;
+
+  // Data members
+  int id_ {-1};  //!< User-specified ID
+  int n_dimension_; //!< Number of dimensions
+};
+
+class StructuredMesh : public Mesh
+{
+public:
+  StructuredMesh() = default;
+  StructuredMesh(pugi::xml_node node) : Mesh {node} {};
+  virtual ~StructuredMesh() = default;
+
   //! Get bin given mesh indices
   //
   //! \param[in] Array of mesh indices
@@ -88,32 +125,21 @@ public:
   //! \param[out] ijk Mesh indices
   virtual void get_indices_from_bin(int bin, int* ijk) const = 0;
 
-  //! Get the number of mesh cells.
-  virtual int n_bins() const = 0;
+  //! Get a label for the mesh bin
+  std::string bin_label(int bin) const override {
+    std::vector<int> ijk(n_dimension_);
+    get_indices_from_bin(bin, ijk.data());
 
-  //! Get the number of mesh cell surfaces.
-  virtual int n_surface_bins() const = 0;
-
-  //! Find the mesh lines that intersect an axis-aligned slice plot
-  //
-  //! \param[in] plot_ll The lower-left coordinates of the slice plot.
-  //! \param[in] plot_ur The upper-right coordinates of the slice plot.
-  //! \return A pair of vectors indicating where the mesh lines lie along each
-  //!   of the plot's axes.  For example an xy-slice plot will get back a vector
-  //!   of x-coordinates and another of y-coordinates.  These vectors may be
-  //!   empty for low-dimensional meshes.
-  virtual std::pair<std::vector<double>, std::vector<double>>
-  plot(Position plot_ll, Position plot_ur) const = 0;
-
-  //! Write mesh data to an HDF5 group
-  //
-  //! \param[in] group HDF5 group
-  virtual void to_hdf5(hid_t group) const = 0;
+    if (n_dim > 2) {
+      return fmt::format("Mesh Index ({}, {}, {})", ijk[0], ijk[1], ijk[2]);
+    } else if (n_dim > 1) {
+      return fmt::format("Mesh Index ({}, {})", ijk[0], ijk[1]);
+    } else {
+      return fmt::format("Mesh Index ({})", ijk[0]) ;
+    }
+  }
 
   // Data members
-
-  int id_ {-1};  //!< User-specified ID
-  int n_dimension_; //!< Number of dimensions
   xt::xtensor<double, 1> lower_left_; //!< Lower-left coordinates of mesh
   xt::xtensor<double, 1> upper_right_; //!< Upper-right coordinates of mesh
 };
@@ -122,7 +148,7 @@ public:
 //! Tessellation of n-dimensional Euclidean space by congruent squares or cubes
 //==============================================================================
 
-class RegularMesh : public Mesh
+class RegularMesh : public StructuredMesh
 {
 public:
   // Constructors
@@ -188,7 +214,7 @@ private:
 };
 
 
-class RectilinearMesh : public Mesh
+class RectilinearMesh : public StructuredMesh
 {
 public:
   // Constructors
@@ -245,6 +271,7 @@ class UnstructuredMesh : public Mesh {
 public:
   UnstructuredMesh() { };
   UnstructuredMesh(pugi::xml_node);
+  ~UnstructuredMesh() = default;
 
   //! Determine which bins were crossed by a particle.
   //
@@ -316,14 +343,11 @@ intersect_track(const moab::CartVect& start,
   //! \return MOAB EntityHandle of tet
   moab::EntityHandle get_ent_handle_from_bin(int bin) const;
 
-  int get_bin_from_indices(const int* ijk) const override;
+  int get_bin_from_index(int idx) const;
 
-  void get_indices(Position r, int* ijk, bool* in_mesh) const override;
+  int get_index(Position r, bool* in_mesh) const;
 
-  void get_indices_from_bin(int bin, int* ijk) const override;
-
-  std::pair<std::vector<double>, std::vector<double>>
-  plot(Position plot_ll, Position plot_ur) const override;
+  int get_index_from_bin(int bin) const;
 
   //! Builds a KDTree for all tetrahedra in the mesh. All
   //! triangles representing 2D faces of the mesh are
@@ -333,6 +357,10 @@ intersect_track(const moab::CartVect& start,
   void build_kdtree(const moab::Range& all_tets);
 
 public:
+
+  std::pair<std::vector<double>, std::vector<double>>
+  plot(Position plot_ll, Position plot_ur) const override;
+
   //! Determine which surface bins were crossed by a particle.
   //
   //! \param[in] p Particle to check
@@ -356,6 +384,10 @@ public:
 
   int n_surface_bins() const override;
 
+  Position centroid(moab::EntityHandle tet) const;
+
+  std::string bin_label(int bin) const override;
+
   std::string filename_; //<! Path to unstructured mesh file
 
 private:

src/mesh.cpp

@@ -85,7 +85,7 @@ Mesh::Mesh(pugi::xml_node node) {
 //==============================================================================
 
 RegularMesh::RegularMesh(pugi::xml_node node)
-  : Mesh {node}
+  : StructuredMesh {node}
 {
   // Determine number of dimensions for mesh
   if (check_for_node(node, "dimension")) {
@@ -851,7 +851,7 @@ RegularMesh::count_sites(const Particle::Bank* bank, int64_t length,
 //==============================================================================
 
 RectilinearMesh::RectilinearMesh(pugi::xml_node node)
-  : Mesh {node}
+  : StructuredMesh {node}
 {
   n_dimension_ = 3;
 
@@ -1535,7 +1535,7 @@ UnstructuredMesh::UnstructuredMesh(pugi::xml_node node) : Mesh(node) {
   }
 
   // create MOAB instance
-  mbi_ = std::shared_ptr<moab::Interface>(new moab::Core());
+  mbi_ = std::unique_ptr<moab::Interface>(new moab::Core());
   // create meshset to load mesh into
   moab::ErrorCode rval = mbi_->create_meshset(moab::MESHSET_SET, meshset_);
   if (rval != moab::MB_SUCCESS) {
@@ -1855,7 +1855,6 @@ UnstructuredMesh::compute_barycentric_data(const moab::Range& all_tets) {
     a = a.transpose().inverse();
     baryc_data_.at(get_bin_from_ent_handle(tet)) = a;
   }
-
 }
 
 // TODO: write this function
@@ -1911,25 +1910,24 @@ UnstructuredMesh::point_in_tet(const moab::CartVect& r, moab::EntityHandle tet)
 }
 
 int
-UnstructuredMesh::get_bin_from_indices(const int* ijk) const {
-  if (ijk[0] >= n_bins()) {
+UnstructuredMesh::get_bin_from_index(int idx) const {
+  if (idx >= n_bins()) {
     std::stringstream s;
-    s << "Invalid bin: " << ijk[0];
+    s << "Invalid bin index: " << idx;
     fatal_error(s);
   }
-  int bin = ehs_[ijk[0]] - ehs_[0];
+  return ehs_[idx] - ehs_[0];
+}
+
+int
+UnstructuredMesh::get_index(Position r, bool* in_mesh) const {
+  int bin = get_bin(r);
+  *in_mesh = bin != -1;
   return bin;
 }
 
-void
-UnstructuredMesh::get_indices(Position r, int* ijk, bool* in_mesh) const {
-  int bin = get_bin(r);
-  ijk[0]= bin;
-  *in_mesh = bin != -1;
-}
-
-void UnstructuredMesh::get_indices_from_bin(int bin, int* ijk) const {
-  ijk[0] = bin;
+int UnstructuredMesh::get_index_from_bin(int bin) const {
+  return bin;
 }
 
 std::pair<std::vector<double>, std::vector<double>>
@@ -1972,6 +1970,42 @@ int UnstructuredMesh::n_surface_bins() const {
   return 2 * tris.size();
 }
 
+Position
+UnstructuredMesh::centroid(moab::EntityHandle tet) const {
+  moab::ErrorCode rval;
+
+  // look up the tet connectivity
+  std::vector<moab::EntityHandle> conn;
+  rval = mbi_->get_connectivity(&tet, 1, conn);
+  if (rval != moab::MB_SUCCESS) {
+    warning("Failed to get connectivity of a mesh element.");
+    return {};
+  }
+
+  // get the coordinates
+  std::vector<moab::CartVect> coords(conn.size());
+  rval = mbi_->get_coords(&conn.front(), conn.size(), coords[0].array());
+  if (rval != moab::MB_SUCCESS) {
+    warning("Failed to get the coordinates of a mesh element.");
+    return {};
+  }
+
+  // compute the centroid of the elements
+  moab::CartVect centroid(0.0);
+  for(const auto& coord : coords) {
+    centroid += coord;
+  }
+  centroid /= double(coords.size());
+
+  return {centroid[0], centroid[1], centroid[2]};
+}
+
+std::string
+UnstructuredMesh::bin_label(int bin) const {
+  std::stringstream out;
+  out << "Mesh Index (" << bin << ")";
+};
+
 #endif
 
 

src/tallies/filter_mesh.cpp

@@ -55,18 +55,7 @@ std::string
 MeshFilter::text_label(int bin) const
 {
   auto& mesh = *model::meshes[mesh_];
-  int n_dim = mesh.n_dimension_;
-
-  std::vector<int> ijk(n_dim);
-  mesh.get_indices_from_bin(bin, ijk.data());
-
-  if (n_dim > 2) {
-    return fmt::format("Mesh Index ({}, {}, {})", ijk[0], ijk[1], ijk[2]);
-  } else if (n_dim > 1) {
-    return fmt::format("Mesh Index ({}, {})", ijk[0], ijk[1]);
-  } else {
-    return fmt::format("Mesh Index ({})", ijk[0]) ;
-  }
+  return mesh.bin_label(bin);
 }
 
 void