From a0d3d1132813feffd0289ecdf10b17eca86e5e96 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Wed, 29 May 2019 11:57:13 -0400 Subject: [PATCH] Clean up RectilinearMesh implementation --- include/openmc/mesh.h | 20 ++++- openmc/mesh.py | 82 +++++++++++++++---- src/mesh.cpp | 127 ++++++++++++++++++----------- src/tallies/filter_mesh.cpp | 39 +++------ src/tallies/filter_meshsurface.cpp | 17 +--- 5 files changed, 175 insertions(+), 110 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index ecb7a59ca..12158339d 100644 --- a/include/openmc/mesh.h +++ b/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 - count_sites(const std::vector& bank, bool* outside) const; + xt::xarray count_sites(const std::vector& bank, + bool* outside) const; // Data members double volume_frac_; //!< Volume fraction of each mesh element - int n_dimension_; //!< Number of dimensions xt::xarray shape_; //!< Number of mesh elements in each dimension xt::xarray lower_left_; //!< Lower-left coordinates of mesh xt::xarray 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 diff --git a/openmc/mesh.py b/openmc/mesh.py index 0961b238e..ef3212b90 100644 --- a/openmc/mesh.py +++ b/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 ')) diff --git a/src/mesh.cpp b/src/mesh.cpp index 0e3198c63..2cb91b98f 100644 --- a/src/mesh.cpp +++ b/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& bank, RectilinearMesh::RectilinearMesh(pugi::xml_node node) : Mesh {node} { + n_dimension_ = 3; + grid_.resize(3); grid_[0] = get_node_array(node, "x_grid"); grid_[1] = get_node_array(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(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(model::meshes[index].get())->shape_.data(); - *n = dynamic_cast(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 shape = {static_cast(n)}; - auto& m = *dynamic_cast(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(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(model::meshes[index].get()); std::vector shape = {static_cast(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; diff --git a/src/tallies/filter_mesh.cpp b/src/tallies/filter_mesh.cpp index 9658e552e..5dce9c618 100644 --- a/src/tallies/filter_mesh.cpp +++ b/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(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(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(m_ptr)) { - n_bins_ = 1; - for (auto dim : m->shape_) n_bins_ *= dim; - } else if (auto* m = dynamic_cast(m_ptr)) { - n_bins_ = 1; - for (auto dim : m->shape_) n_bins_ *= dim; - } + n_bins_ = model::meshes[mesh_]->n_bins(); } //============================================================================== diff --git a/src/tallies/filter_meshsurface.cpp b/src/tallies/filter_meshsurface.cpp index aa6452013..edb0492e9 100644 --- a/src/tallies/filter_meshsurface.cpp +++ b/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(model::meshes[mesh_].get())) { - n_dim = mesh->n_dimension_; - } else if (auto* mesh = dynamic_cast(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(m_ptr)) { - n_bins_ = 4 * m->n_dimension_; - for (auto dim : m->shape_) n_bins_ *= dim; - } else if (auto* m = dynamic_cast(m_ptr)) { - n_bins_ = 4 * 3; - for (auto dim : m->shape_) n_bins_ *= dim; - } + n_bins_ = model::meshes[mesh_]->n_surface_bins(); } //==============================================================================