Clean up RectilinearMesh implementation

include/openmc/mesh.h

@@ -79,6 +79,12 @@ 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 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
@@ -87,6 +93,7 @@ public:
   // Data members
 
   int id_ {-1};  //!< User-specified ID
+  int n_dimension_; //!< Number of dimensions
 };
 
 //==============================================================================
@@ -116,6 +123,10 @@ public:
 
   void get_indices_from_bin(int bin, int* ijk) const override;
 
+  int n_bins() const override;
+
+  int n_surface_bins() const override;
+
   void to_hdf5(hid_t group) const override;
 
   // New methods
@@ -133,13 +144,12 @@ public:
   //! \param[in] bank Array of bank sites
   //! \param[out] Whether any bank sites are outside the mesh
   //! \return Array indicating number of sites in each mesh/energy bin
-  virtual xt::xarray<double>
-  count_sites(const std::vector<Particle::Bank>& bank, bool* outside) const;
+  xt::xarray<double> count_sites(const std::vector<Particle::Bank>& bank,
+    bool* outside) const;
 
   // Data members
 
   double volume_frac_; //!< Volume fraction of each mesh element
-  int n_dimension_; //!< Number of dimensions
   xt::xarray<int> shape_; //!< Number of mesh elements in each dimension
   xt::xarray<double> lower_left_; //!< Lower-left coordinates of mesh
   xt::xarray<double> upper_right_; //!< Upper-right coordinates of mesh
@@ -173,6 +183,10 @@ public:
 
   void get_indices_from_bin(int bin, int* ijk) const override;
 
+  int n_bins() const override;
+
+  int n_surface_bins() const override;
+
   void to_hdf5(hid_t group) const override;
 
   // New methods

openmc/mesh.py

@@ -77,7 +77,7 @@ class MeshBase(IDManagerMixin, metaclass=ABCMeta):
 
 
 class Mesh(MeshBase):
-    """A structured Cartesian mesh in one, two, or three dimensions
+    """A regular Cartesian mesh in one, two, or three dimensions
 
     Parameters
     ----------
@@ -92,10 +92,10 @@ class Mesh(MeshBase):
         Unique identifier for the mesh
     name : str
         Name of the mesh
-    type : str
-        Type of the mesh
     dimension : Iterable of int
         The number of mesh cells in each direction.
+    n_dimension : int
+        Number of mesh dimensions.
     lower_left : Iterable of float
         The lower-left corner of the structured mesh. If only two coordinate are
         given, it is assumed that the mesh is an x-y mesh.
@@ -114,16 +114,11 @@ class Mesh(MeshBase):
         super().__init__(mesh_id, name)
 
         # Initialize Mesh class attributes
-        self._type = 'regular'
         self._dimension = None
         self._lower_left = None
         self._upper_right = None
         self._width = None
 
-    @property
-    def type(self):
-        return self._type
-
     @property
     def dimension(self):
         return self._dimension
@@ -166,14 +161,6 @@ class Mesh(MeshBase):
             nx, = self.dimension
             return ((x,) for x in range(1, nx + 1))
 
-    @type.setter
-    def type(self, meshtype):
-        cv.check_type('type for mesh ID="{0}"'.format(self._id),
-                      meshtype, str)
-        cv.check_value('type for mesh ID="{0}"'.format(self._id),
-                       meshtype, ['regular'])
-        self._type = meshtype
-
     @dimension.setter
     def dimension(self, dimension):
         cv.check_type('mesh dimension', dimension, Iterable, Integral)
@@ -215,7 +202,6 @@ class Mesh(MeshBase):
 
         # Read and assign mesh properties
         mesh = cls(mesh_id)
-        mesh.type = group['type'][()].decode()
         mesh.dimension = group['dimension'][()]
         mesh.lower_left = group['lower_left'][()]
         mesh.upper_right = group['upper_right'][()]
@@ -269,7 +255,7 @@ class Mesh(MeshBase):
 
         element = ET.Element("mesh")
         element.set("id", str(self._id))
-        element.set("type", self._type)
+        element.set("type", "regular")
 
         subelement = ET.SubElement(element, "dimension")
         subelement.text = ' '.join(map(str, self._dimension))
@@ -418,13 +404,58 @@ class Mesh(MeshBase):
 
 
 class RectilinearMesh(MeshBase):
+    """A 3D rectilinear Cartesian mesh
+
+    Parameters
+    ----------
+    mesh_id : int
+        Unique identifier for the mesh
+    name : str
+        Name of the mesh
+
+    Attributes
+    ----------
+    id : int
+        Unique identifier for the mesh
+    name : str
+        Name of the mesh
+    n_dimension : int
+        Number of mesh dimensions (always 3 for a RectilinearMesh).
+    x_grid : Iterable of float
+        Mesh boundary points along the x-axis.
+    y_grid : Iterable of float
+        Mesh boundary points along the y-axis.
+    z_grid : Iterable of float
+        Mesh boundary points along the z-axis.
+    indices : list of tuple
+        A list of mesh indices for each mesh element, e.g. [(1, 1, 1), (2, 1,
+        1), ...]
+
+    """
+
     def __init__(self, mesh_id=None, name=''):
         super().__init__(mesh_id, name)
 
+        self._x_grid = None
+        self._y_grid = None
+        self._z_grid = None
+
     @property
     def n_dimension(self):
         return 3
 
+    @property
+    def x_grid(self):
+        return self._x_grid
+
+    @property
+    def y_grid(self):
+        return self._y_grid
+
+    @property
+    def z_grid(self):
+        return self._z_grid
+
     @property
     def indices(self):
         nx = len(self.x_grid) - 1
@@ -435,6 +466,21 @@ class RectilinearMesh(MeshBase):
                 for y in range(1, ny + 1)
                 for x in range(1, nx + 1))
 
+    @x_grid.setter
+    def x_grid(self, grid):
+        cv.check_type('mesh x_grid', grid, Iterable, Real)
+        self._x_grid = grid
+
+    @y_grid.setter
+    def y_grid(self, grid):
+        cv.check_type('mesh y_grid', grid, Iterable, Real)
+        self._y_grid = grid
+
+    @z_grid.setter
+    def z_grid(self, grid):
+        cv.check_type('mesh z_grid', grid, Iterable, Real)
+        self._z_grid = grid
+
     @classmethod
     def from_hdf5(cls, group):
         mesh_id = int(group.name.split('/')[-1].lstrip('mesh '))

src/mesh.cpp

@@ -228,6 +228,18 @@ void RegularMesh::get_indices_from_bin(int bin, int* ijk) const
   }
 }
 
+int RegularMesh::n_bins() const
+{
+  int n_bins = 1;
+  for (auto dim : shape_) n_bins *= dim;
+  return n_bins;
+}
+
+int RegularMesh::n_surface_bins() const
+{
+  return 4 * n_dimension_ * n_bins();
+}
+
 bool RegularMesh::intersects(Position& r0, Position r1, int* ijk) const
 {
   switch(n_dimension_) {
@@ -788,6 +800,8 @@ RegularMesh::count_sites(const std::vector<Particle::Bank>& bank,
 RectilinearMesh::RectilinearMesh(pugi::xml_node node)
   : Mesh {node}
 {
+  n_dimension_ = 3;
+
   grid_.resize(3);
   grid_[0] = get_node_array<double>(node, "x_grid");
   grid_[1] = get_node_array<double>(node, "y_grid");
@@ -1099,6 +1113,18 @@ void RectilinearMesh::get_indices_from_bin(int bin, int* ijk) const
   ijk[2] = bin / (shape_[0] * shape_[1]) + 1;
 }
 
+int RectilinearMesh::n_bins() const
+{
+  int n_bins = 1;
+  for (auto dim : shape_) n_bins *= dim;
+  return n_bins;
+}
+
+int RectilinearMesh::n_surface_bins() const
+{
+  return 4 * n_dimension_ * n_bins();
+}
+
 void RectilinearMesh::to_hdf5(hid_t group) const
 {
   hid_t mesh_group = create_group(group, "mesh " + std::to_string(id_));
@@ -1222,6 +1248,32 @@ bool RectilinearMesh::intersects(Position& r0, Position r1, int* ijk) const
   return min_dist < INFTY;
 }
 
+//==============================================================================
+// Helper functions for the C API
+//==============================================================================
+
+int
+check_mesh(int32_t index)
+{
+  if (index < 0 || index >= model::meshes.size()) {
+    set_errmsg("Index in meshes array is out of bounds.");
+    return OPENMC_E_OUT_OF_BOUNDS;
+  }
+  return 0;
+}
+
+int
+check_regular_mesh(int32_t index, RegularMesh** mesh)
+{
+  if (int err = check_mesh(index)) return err;
+  *mesh = dynamic_cast<RegularMesh*>(model::meshes[index].get());
+  if (!*mesh) {
+    set_errmsg("This function is only valid for regular meshes.");
+    return OPENMC_E_INVALID_TYPE;
+  }
+  return 0;
+}
+
 //==============================================================================
 // C API functions
 //==============================================================================
@@ -1256,10 +1308,7 @@ openmc_get_mesh_index(int32_t id, int32_t* index)
 extern "C" int
 openmc_mesh_get_id(int32_t index, int32_t* id)
 {
-  if (index < 0 || index >= model::meshes.size()) {
-    set_errmsg("Index in meshes array is out of bounds.");
-    return OPENMC_E_OUT_OF_BOUNDS;
-  }
+  if (int err = check_mesh(index)) return err;
   *id = model::meshes[index]->id_;
   return 0;
 }
@@ -1268,10 +1317,7 @@ openmc_mesh_get_id(int32_t index, int32_t* id)
 extern "C" int
 openmc_mesh_set_id(int32_t index, int32_t id)
 {
-  if (index < 0 || index >= model::meshes.size()) {
-    set_errmsg("Index in meshes array is out of bounds.");
-    return OPENMC_E_OUT_OF_BOUNDS;
-  }
+  if (int err = check_mesh(index)) return err;
   model::meshes[index]->id_ = id;
   model::mesh_map[id] = index;
   return 0;
@@ -1281,12 +1327,10 @@ openmc_mesh_set_id(int32_t index, int32_t id)
 extern "C" int
 openmc_mesh_get_dimension(int32_t index, int** dims, int* n)
 {
-  if (index < 0 || index >= model::meshes.size()) {
-    set_errmsg("Index in meshes array is out of bounds.");
-    return OPENMC_E_OUT_OF_BOUNDS;
-  }
-  *dims = dynamic_cast<RegularMesh*>(model::meshes[index].get())->shape_.data();
-  *n = dynamic_cast<RegularMesh*>(model::meshes[index].get())->n_dimension_;
+  RegularMesh* mesh;
+  if (int err = check_regular_mesh(index, &mesh)) return err;
+  *dims = mesh->shape_.data();
+  *n = mesh->n_dimension_;
   return 0;
 }
 
@@ -1294,16 +1338,13 @@ openmc_mesh_get_dimension(int32_t index, int** dims, int* n)
 extern "C" int
 openmc_mesh_set_dimension(int32_t index, int n, const int* dims)
 {
-  if (index < 0 || index >= model::meshes.size()) {
-    set_errmsg("Index in meshes array is out of bounds.");
-    return OPENMC_E_OUT_OF_BOUNDS;
-  }
+  RegularMesh* mesh;
+  if (int err = check_regular_mesh(index, &mesh)) return err;
 
   // Copy dimension
   std::vector<std::size_t> shape = {static_cast<std::size_t>(n)};
-  auto& m = *dynamic_cast<RegularMesh*>(model::meshes[index].get());
-  m.shape_ = xt::adapt(dims, n, xt::no_ownership(), shape);
-  m.n_dimension_ = m.shape_.size();
+  mesh->shape_ = xt::adapt(dims, n, xt::no_ownership(), shape);
+  mesh->n_dimension_ = mesh->shape_.size();
 
   return 0;
 }
@@ -1312,22 +1353,18 @@ openmc_mesh_set_dimension(int32_t index, int n, const int* dims)
 extern "C" int
 openmc_mesh_get_params(int32_t index, double** ll, double** ur, double** width, int* n)
 {
-  if (index < 0 || index >= model::meshes.size()) {
-    set_errmsg("Index in meshes array is out of bounds.");
-    return OPENMC_E_OUT_OF_BOUNDS;
-  }
+  RegularMesh* m;
+  if (int err = check_regular_mesh(index, &m)) return err;
 
-  //auto& m = model::meshes[index];
-  auto& m = *dynamic_cast<RegularMesh*>(model::meshes[index].get());
-  if (m.lower_left_.dimension() == 0) {
+  if (m->lower_left_.dimension() == 0) {
     set_errmsg("Mesh parameters have not been set.");
     return OPENMC_E_ALLOCATE;
   }
 
-  *ll = m.lower_left_.data();
-  *ur = m.upper_right_.data();
-  *width = m.width_.data();
-  *n = m.n_dimension_;
+  *ll = m->lower_left_.data();
+  *ur = m->upper_right_.data();
+  *width = m->width_.data();
+  *n = m->n_dimension_;
   return 0;
 }
 
@@ -1336,26 +1373,22 @@ extern "C" int
 openmc_mesh_set_params(int32_t index, int n, const double* ll, const double* ur,
                        const double* width)
 {
-  if (index < 0 || index >= model::meshes.size()) {
-    set_errmsg("Index in meshes array is out of bounds.");
-    return OPENMC_E_OUT_OF_BOUNDS;
-  }
+  RegularMesh* m;
+  if (int err = check_regular_mesh(index, &m)) return err;
 
-  //auto& m = model::meshes[index];
-  auto& m = *dynamic_cast<RegularMesh*>(model::meshes[index].get());
   std::vector<std::size_t> shape = {static_cast<std::size_t>(n)};
   if (ll && ur) {
-    m.lower_left_ = xt::adapt(ll, n, xt::no_ownership(), shape);
-    m.upper_right_ = xt::adapt(ur, n, xt::no_ownership(), shape);
-    m.width_ = (m.upper_right_ - m.lower_left_) / m.shape_;
+    m->lower_left_ = xt::adapt(ll, n, xt::no_ownership(), shape);
+    m->upper_right_ = xt::adapt(ur, n, xt::no_ownership(), shape);
+    m->width_ = (m->upper_right_ - m->lower_left_) / m->shape_;
   } else if (ll && width) {
-    m.lower_left_ = xt::adapt(ll, n, xt::no_ownership(), shape);
-    m.width_ = xt::adapt(width, n, xt::no_ownership(), shape);
-    m.upper_right_ = m.lower_left_ + m.shape_ * m.width_;
+    m->lower_left_ = xt::adapt(ll, n, xt::no_ownership(), shape);
+    m->width_ = xt::adapt(width, n, xt::no_ownership(), shape);
+    m->upper_right_ = m->lower_left_ + m->shape_ * m->width_;
   } else if (ur && width) {
-    m.upper_right_ = xt::adapt(ur, n, xt::no_ownership(), shape);
-    m.width_ = xt::adapt(width, n, xt::no_ownership(), shape);
-    m.lower_left_ = m.upper_right_ - m.shape_ * m.width_;
+    m->upper_right_ = xt::adapt(ur, n, xt::no_ownership(), shape);
+    m->width_ = xt::adapt(width, n, xt::no_ownership(), shape);
+    m->lower_left_ = m->upper_right_ - m->shape_ * m->width_;
   } else {
     set_errmsg("At least two parameters must be specified.");
     return OPENMC_E_INVALID_ARGUMENT;

src/tallies/filter_mesh.cpp

@@ -55,43 +55,26 @@ MeshFilter::to_statepoint(hid_t filter_group) const
 std::string
 MeshFilter::text_label(int bin) const
 {
-  if (auto* mesh = dynamic_cast<RegularMesh*>(model::meshes[mesh_].get())) {
-    int n_dim = mesh->n_dimension_;
+  auto& mesh = *model::meshes[mesh_];
+  int n_dim = mesh.n_dimension_;
 
-    int ijk[n_dim];
-    mesh->get_indices_from_bin(bin, ijk);
+  int ijk[n_dim];
+  mesh.get_indices_from_bin(bin, ijk);
 
-    std::stringstream out;
-    out << "Mesh Index (" << ijk[0];
-    if (n_dim > 1) out << ", " << ijk[1];
-    if (n_dim > 2) out << ", " << ijk[2];
-    out << ")";
+  std::stringstream out;
+  out << "Mesh Index (" << ijk[0];
+  if (n_dim > 1) out << ", " << ijk[1];
+  if (n_dim > 2) out << ", " << ijk[2];
+  out << ")";
 
-    return out.str();
-
-  } else if (auto* mesh = dynamic_cast<RectilinearMesh*>(model::meshes[mesh_].get())) {
-    int ijk[3];
-    mesh->get_indices_from_bin(bin, ijk);
-
-    std::stringstream out;
-    out << "Mesh Index (" << ijk[0] << ", " << ijk[1] << ", " << ijk[2] << ")";
-
-    return out.str();
-  }
+  return out.str();
 }
 
 void
 MeshFilter::set_mesh(int32_t mesh)
 {
   mesh_ = mesh;
-  auto* m_ptr = model::meshes[mesh_].get();
-  if (auto* m = dynamic_cast<RegularMesh*>(m_ptr)) {
-    n_bins_ = 1;
-    for (auto dim : m->shape_) n_bins_ *= dim;
-  } else if (auto* m = dynamic_cast<RectilinearMesh*>(m_ptr)) {
-    n_bins_ = 1;
-    for (auto dim : m->shape_) n_bins_ *= dim;
-  }
+  n_bins_ = model::meshes[mesh_]->n_bins();
 }
 
 //==============================================================================

src/tallies/filter_meshsurface.cpp

@@ -18,12 +18,8 @@ MeshSurfaceFilter::get_all_bins(const Particle* p, int estimator,
 std::string
 MeshSurfaceFilter::text_label(int bin) const
 {
-  int n_dim;
-  if (auto* mesh = dynamic_cast<RegularMesh*>(model::meshes[mesh_].get())) {
-    n_dim = mesh->n_dimension_;
-  } else if (auto* mesh = dynamic_cast<RectilinearMesh*>(model::meshes[mesh_].get())) {
-    n_dim = 3;
-  }
+  auto& mesh = *model::meshes[mesh_];
+  int n_dim = mesh.n_dimension_;
 
   // Get flattend mesh index and surface index.
   int i_mesh = bin / (4 * n_dim);
@@ -79,14 +75,7 @@ void
 MeshSurfaceFilter::set_mesh(int32_t mesh)
 {
   mesh_ = mesh;
-  auto* m_ptr = model::meshes[mesh_].get();
-  if (auto* m = dynamic_cast<RegularMesh*>(m_ptr)) {
-    n_bins_ = 4 * m->n_dimension_;
-    for (auto dim : m->shape_) n_bins_ *= dim;
-  } else if (auto* m = dynamic_cast<RectilinearMesh*>(m_ptr)) {
-    n_bins_ = 4 * 3;
-    for (auto dim : m->shape_) n_bins_ *= dim;
-  }
+  n_bins_ = model::meshes[mesh_]->n_surface_bins();
 }
 
 //==============================================================================