diff --git a/docs/source/io_formats/tallies.rst b/docs/source/io_formats/tallies.rst
index 7f95f29850..3cde6fc854 100644
--- a/docs/source/io_formats/tallies.rst
+++ b/docs/source/io_formats/tallies.rst
@@ -312,8 +312,8 @@ element with the tag name ````. This element has the following
attributes/sub-elements:
:type:
- The type of mesh. This can be either "regular", "rectilinear", or
- "unstructured".
+ The type of mesh. This can be either "regular", "rectilinear",
+ "cylindrical", "spherical", or "unstructured".
:dimension:
The number of mesh cells in each direction. (For regular mesh only.)
@@ -340,7 +340,16 @@ attributes/sub-elements:
The mesh divisions along the y-axis. (For rectilinear mesh only.)
:z_grid:
- The mesh divisions along the z-axis. (For rectilinear mesh only.)
+ The mesh divisions along the z-axis. (For rectilinear and cylindrical meshes only.)
+
+ :r_grid:
+ The mesh divisions along the r-axis. (For cylindrical and spherical meshes only.)
+
+ :phi_grid:
+ The mesh divisions along the phi-axis. (For cylindrical and spherical meshes only.)
+
+ :theta_grid:
+ The mesh divisions along the theta-axis. (For spherical mesh only.)
:library:
The mesh library used to represent an unstructured mesh. This can be either
diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h
index d34ad2087b..ae1e648b38 100644
--- a/include/openmc/mesh.h
+++ b/include/openmc/mesh.h
@@ -123,6 +123,9 @@ public:
//! \param[in] bin Mesh bin to generate a label for
virtual std::string bin_label(int bin) const = 0;
+ //! Return the mesh type
+ virtual std::string get_mesh_type() const = 0;
+
// Data members
int id_ {-1}; //!< User-specified ID
int n_dimension_; //!< Number of dimensions
@@ -134,6 +137,22 @@ public:
StructuredMesh(pugi::xml_node node) : Mesh {node} {};
virtual ~StructuredMesh() = default;
+ using MeshIndex = std::array;
+
+ struct MeshDistance {
+ MeshDistance() = default;
+ MeshDistance(int _index, bool _max_surface, double _distance)
+ : next_index {_index}, max_surface {_max_surface}, distance {_distance}
+ {}
+ int next_index {-1};
+ double distance {INFTY};
+ bool max_surface {true};
+ bool operator<(const MeshDistance& o) const
+ {
+ return distance < o.distance;
+ }
+ };
+
int get_bin(Position r) const override;
int n_bins() const override;
@@ -143,6 +162,19 @@ public:
void bins_crossed(Position r0, Position r1, const Direction& u,
vector& bins, vector& lengths) const override;
+ void surface_bins_crossed(Position r0, Position r1, const Direction& u,
+ vector& bins) const override;
+
+ //! Determine which cell or surface bins were crossed by a particle
+ //
+ //! \param[in] r0 Previous position of the particle
+ //! \param[in] r1 Current position of the particle
+ //! \param[in] u Particle direction
+ //! \param[in] tally Functor that eventually stores the tally data
+ template
+ void raytrace_mesh(
+ Position r0, Position r1, const Direction& u, T tally) const;
+
//! Count number of bank sites in each mesh bin / energy bin
//
//! \param[in] Pointer to bank sites
@@ -155,20 +187,20 @@ public:
//
//! \param[in] Array of mesh indices
//! \return Mesh bin
- virtual int get_bin_from_indices(const int* ijk) const;
+ virtual int get_bin_from_indices(const MeshIndex& ijk) const;
//! Get mesh indices given a position
//
//! \param[in] r Position to get indices for
- //! \param[out] ijk Array of mesh indices
//! \param[out] in_mesh Whether position is in mesh
- virtual void get_indices(Position r, int* ijk, bool* in_mesh) const;
+ //! \return Array of mesh indices
+ virtual MeshIndex get_indices(Position r, bool& in_mesh) const;
//! Get mesh indices corresponding to a mesh bin
//
//! \param[in] bin Mesh bin
- //! \param[out] ijk Mesh indices
- virtual void get_indices_from_bin(int bin, int* ijk) const;
+ //! \return ijk Mesh indices
+ virtual MeshIndex get_indices_from_bin(int bin) const;
//! Get mesh index in a particular direction
//!
@@ -176,38 +208,34 @@ public:
//! \param[in] i Direction index
virtual int get_index_in_direction(double r, int i) const = 0;
- //! Check where a line segment intersects the mesh and if it intersects at all
- //
- //! \param[in,out] r0 In: starting position, out: intersection point
- //! \param[in] r1 Ending position
- //! \param[out] ijk Indices of the mesh bin containing the intersection point
- //! \return Whether the line segment connecting r0 and r1 intersects mesh
- virtual bool intersects(Position& r0, Position r1, int* ijk) const;
-
- //! Get the coordinate for the mesh grid boundary in the positive direction
+ //! Get the closest distance from the coordinate r to the grid surface
+ //! in i direction that bounds mesh cell ijk and that is larger than l
+ //! The coordinate r does not have to be inside the mesh cell ijk. In
+ //! curved coordinates, multiple crossings of the same surface can happen,
+ //! these are selected by the parameter l
//!
//! \param[in] ijk Array of mesh indices
- //! \param[in] i Direction index
- virtual double positive_grid_boundary(int* ijk, int i) const = 0;
-
- //! Get the coordinate for the mesh grid boundary in the negative direction
- //!
- //! \param[in] ijk Array of mesh indices
- //! \param[in] i Direction index
- virtual double negative_grid_boundary(int* ijk, int i) const = 0;
+ //! \param[in] i direction index of grid surface
+ //! \param[in] r0 position, from where to calculate the distance
+ //! \param[in] u direction of flight. actual position is r0 + l * u
+ //! \param[in] l actual chord length
+ //! \return MeshDistance struct with closest distance, next cell index in
+ //! i-direction and min/max surface indicator
+ virtual MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i,
+ const Position& r0, const Direction& u, double l) const = 0;
//! Get a label for the mesh bin
std::string bin_label(int bin) const override;
+ //! Get shape as xt::xtensor
+ xt::xtensor get_x_shape() const;
+
// Data members
xt::xtensor lower_left_; //!< Lower-left coordinates of mesh
xt::xtensor upper_right_; //!< Upper-right coordinates of mesh
- xt::xtensor shape_; //!< Number of mesh elements in each dimension
+ std::array shape_; //!< Number of mesh elements in each dimension
protected:
- virtual bool intersects_1d(Position& r0, Position r1, int* ijk) const;
- virtual bool intersects_2d(Position& r0, Position r1, int* ijk) const;
- virtual bool intersects_3d(Position& r0, Position r1, int* ijk) const;
};
//==============================================================================
@@ -221,14 +249,14 @@ public:
RegularMesh(pugi::xml_node node);
// Overridden methods
- void surface_bins_crossed(Position r0, Position r1, const Direction& u,
- vector& bins) const override;
-
int get_index_in_direction(double r, int i) const override;
- double positive_grid_boundary(int* ijk, int i) const override;
+ virtual std::string get_mesh_type() const override;
- double negative_grid_boundary(int* ijk, int i) const override;
+ static const std::string mesh_type;
+
+ MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i,
+ const Position& r0, const Direction& u, double l) const override;
std::pair, vector> plot(
Position plot_ll, Position plot_ur) const override;
@@ -236,6 +264,17 @@ public:
void to_hdf5(hid_t group) const override;
// New methods
+ //! Get the coordinate for the mesh grid boundary in the positive direction
+ //!
+ //! \param[in] ijk Array of mesh indices
+ //! \param[in] i Direction index
+ double positive_grid_boundary(const MeshIndex& ijk, int i) const;
+
+ //! Get the coordinate for the mesh grid boundary in the negative direction
+ //!
+ //! \param[in] ijk Array of mesh indices
+ //! \param[in] i Direction index
+ double negative_grid_boundary(const MeshIndex& ijk, int i) const;
//! Count number of bank sites in each mesh bin / energy bin
//
@@ -257,25 +296,151 @@ public:
RectilinearMesh(pugi::xml_node node);
// Overridden methods
- void surface_bins_crossed(Position r0, Position r1, const Direction& u,
- vector& bins) const override;
-
int get_index_in_direction(double r, int i) const override;
- double positive_grid_boundary(int* ijk, int i) const override;
+ virtual std::string get_mesh_type() const override;
- double negative_grid_boundary(int* ijk, int i) const override;
+ static const std::string mesh_type;
+
+ MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i,
+ const Position& r0, const Direction& u, double l) const override;
std::pair, vector> plot(
Position plot_ll, Position plot_ur) const override;
void to_hdf5(hid_t group) const override;
- vector> grid_;
+ // New methods
+ //! Get the coordinate for the mesh grid boundary in the positive direction
+ //!
+ //! \param[in] ijk Array of mesh indices
+ //! \param[in] i Direction index
+ double positive_grid_boundary(const MeshIndex& ijk, int i) const;
+
+ //! Get the coordinate for the mesh grid boundary in the negative direction
+ //!
+ //! \param[in] ijk Array of mesh indices
+ //! \param[in] i Direction index
+ double negative_grid_boundary(const MeshIndex& ijk, int i) const;
+
+ array, 3> grid_;
int set_grid();
};
+class CylindricalMesh : public StructuredMesh {
+public:
+ // Constructors
+ CylindricalMesh() = default;
+ CylindricalMesh(pugi::xml_node node);
+
+ // Overridden methods
+ virtual MeshIndex get_indices(Position r, bool& in_mesh) const override;
+
+ int get_index_in_direction(double r, int i) const override;
+
+ virtual std::string get_mesh_type() const override;
+
+ static const std::string mesh_type;
+
+ MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i,
+ const Position& r0, const Direction& u, double l) const override;
+
+ std::pair, vector> plot(
+ Position plot_ll, Position plot_ur) const override;
+
+ void to_hdf5(hid_t group) const override;
+
+ array, 3> grid_;
+
+ int set_grid();
+
+private:
+ double find_r_crossing(
+ const Position& r, const Direction& u, double l, int shell) const;
+ double find_phi_crossing(
+ const Position& r, const Direction& u, double l, int shell) const;
+ StructuredMesh::MeshDistance find_z_crossing(
+ const Position& r, const Direction& u, double l, int shell) const;
+
+ bool full_phi_ {false};
+
+ constexpr inline int sanitize_angular_index(int idx, bool full, int N) const
+ {
+ if ((idx > 0) and (idx <= N)) {
+ return idx;
+ } else if (full) {
+ return (idx + N - 1) % N + 1;
+ } else {
+ return 0;
+ }
+ }
+
+ inline int sanitize_phi(int idx) const
+ {
+ return sanitize_angular_index(idx, full_phi_, shape_[1]);
+ }
+};
+
+class SphericalMesh : public StructuredMesh {
+public:
+ // Constructors
+ SphericalMesh() = default;
+ SphericalMesh(pugi::xml_node node);
+
+ // Overridden methods
+ virtual MeshIndex get_indices(Position r, bool& in_mesh) const override;
+
+ int get_index_in_direction(double r, int i) const override;
+
+ virtual std::string get_mesh_type() const override;
+
+ static const std::string mesh_type;
+
+ MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i,
+ const Position& r0, const Direction& u, double l) const override;
+
+ std::pair, vector> plot(
+ Position plot_ll, Position plot_ur) const override;
+
+ void to_hdf5(hid_t group) const override;
+
+ array, 3> grid_;
+
+ int set_grid();
+
+private:
+ double find_r_crossing(
+ const Position& r, const Direction& u, double l, int shell) const;
+ double find_theta_crossing(
+ const Position& r, const Direction& u, double l, int shell) const;
+ double find_phi_crossing(
+ const Position& r, const Direction& u, double l, int shell) const;
+
+ bool full_theta_ {false};
+ bool full_phi_ {false};
+
+ constexpr inline int sanitize_angular_index(int idx, bool full, int N) const
+ {
+ if ((idx > 0) and (idx <= N)) {
+ return idx;
+ } else if (full) {
+ return (idx + N - 1) % N + 1;
+ } else {
+ return 0;
+ }
+ }
+
+ inline int sanitize_theta(int idx) const
+ {
+ return sanitize_angular_index(idx, full_theta_, shape_[1]);
+ }
+ inline int sanitize_phi(int idx) const
+ {
+ return sanitize_angular_index(idx, full_phi_, shape_[2]);
+ }
+};
+
// Abstract class for unstructured meshes
class UnstructuredMesh : public Mesh {
@@ -285,6 +450,9 @@ public:
UnstructuredMesh(pugi::xml_node node);
UnstructuredMesh(const std::string& filename);
+ static const std::string mesh_type;
+ virtual std::string get_mesh_type() const override;
+
// Overridden Methods
void surface_bins_crossed(Position r0, Position r1, const Direction& u,
vector& bins) const override;
@@ -356,6 +524,8 @@ public:
MOABMesh(const std::string& filename, double length_multiplier = 1.0);
MOABMesh(std::shared_ptr external_mbi);
+ static const std::string mesh_lib_type;
+
// Overridden Methods
void bins_crossed(Position r0, Position r1, const Direction& u,
@@ -503,6 +673,8 @@ public:
LibMesh(pugi::xml_node node);
LibMesh(const std::string& filename, double length_multiplier = 1.0);
+ static const std::string mesh_lib_type;
+
// Overridden Methods
void bins_crossed(Position r0, Position r1, const Direction& u,
vector& bins, vector& lengths) const override;
diff --git a/openmc/lib/mesh.py b/openmc/lib/mesh.py
index dc528573ac..43f66d7a12 100644
--- a/openmc/lib/mesh.py
+++ b/openmc/lib/mesh.py
@@ -11,7 +11,7 @@ from . import _dll
from .core import _FortranObjectWithID
from .error import _error_handler
-__all__ = ['RegularMesh', 'RectilinearMesh', 'meshes']
+__all__ = ['RegularMesh', 'RectilinearMesh', 'CylindricalMesh', 'SphericalMesh', 'meshes']
# Mesh functions
_dll.openmc_extend_meshes.argtypes = [c_int32, c_char_p, POINTER(c_int32),
@@ -56,6 +56,26 @@ _dll.openmc_regular_mesh_set_params.argtypes = [
_dll.openmc_regular_mesh_set_params.restype = c_int
_dll.openmc_regular_mesh_set_params.errcheck = _error_handler
+_dll.openmc_cylindrical_mesh_get_grid.argtypes = [c_int32,
+ POINTER(POINTER(c_double)), POINTER(c_int), POINTER(POINTER(c_double)),
+ POINTER(c_int), POINTER(POINTER(c_double)), POINTER(c_int)]
+_dll.openmc_cylindrical_mesh_get_grid.restype = c_int
+_dll.openmc_cylindrical_mesh_get_grid.errcheck = _error_handler
+_dll.openmc_cylindrical_mesh_set_grid.argtypes = [c_int32, POINTER(c_double),
+ c_int, POINTER(c_double), c_int, POINTER(c_double), c_int]
+_dll.openmc_cylindrical_mesh_set_grid.restype = c_int
+_dll.openmc_cylindrical_mesh_set_grid.errcheck = _error_handler
+
+_dll.openmc_spherical_mesh_get_grid.argtypes = [c_int32,
+ POINTER(POINTER(c_double)), POINTER(c_int), POINTER(POINTER(c_double)),
+ POINTER(c_int), POINTER(POINTER(c_double)), POINTER(c_int)]
+_dll.openmc_spherical_mesh_get_grid.restype = c_int
+_dll.openmc_spherical_mesh_get_grid.errcheck = _error_handler
+_dll.openmc_spherical_mesh_set_grid.argtypes = [c_int32, POINTER(c_double),
+ c_int, POINTER(c_double), c_int, POINTER(c_double), c_int]
+_dll.openmc_spherical_mesh_set_grid.restype = c_int
+_dll.openmc_spherical_mesh_set_grid.errcheck = _error_handler
+
class Mesh(_FortranObjectWithID):
"""Base class to represent mesh objects
@@ -275,9 +295,184 @@ class RectilinearMesh(Mesh):
ny, z_grid, nz)
+class CylindricalMesh(Mesh):
+ """CylindricalMesh stored internally.
+
+ This class exposes a mesh that is stored internally in the OpenMC
+ library. To obtain a view of a mesh with a given ID, use the
+ :data:`openmc.lib.meshes` mapping.
+
+ Parameters
+ ----------
+ index : int
+ Index in the `meshes` array.
+
+ Attributes
+ ----------
+ id : int
+ ID of the mesh
+ dimension : iterable of int
+ The number of mesh cells in each direction.
+ lower_left : numpy.ndarray
+ The lower-left corner of the structured mesh.
+ upper_right : numpy.ndarray
+ The upper-right corner of the structrued mesh.
+ width : numpy.ndarray
+ The width of mesh cells in each direction.
+
+ """
+ mesh_type = 'cylindrical'
+
+ def __init__(self, uid=None, new=True, index=None):
+ super().__init__(uid, new, index)
+
+ @property
+ def lower_left(self):
+ return self._get_parameters()[0]
+
+ @property
+ def upper_right(self):
+ return self._get_parameters()[1]
+
+ @property
+ def dimension(self):
+ return self._get_parameters()[2]
+
+ @property
+ def width(self):
+ return self._get_parameters()[3]
+
+ def _get_parameters(self):
+ gx = POINTER(c_double)()
+ nx = c_int()
+ gy = POINTER(c_double)()
+ ny = c_int()
+ gz = POINTER(c_double)()
+ nz = c_int()
+ # Call C API to get grid parameters
+ _dll.openmc_cylindrical_mesh_get_grid(self._index, gx, nx, gy, ny, gz,
+ nz)
+
+ # Convert grid parameters to Numpy arrays
+ grid_x = as_array(gx, (nx.value,))
+ grid_y = as_array(gy, (ny.value,))
+ grid_z = as_array(gz, (nz.value,))
+
+ # Calculate lower_left, upper_right, width, and dimension from grid
+ lower_left = np.array((grid_x[0], grid_y[0], grid_z[0]))
+ upper_right = np.array((grid_x[-1], grid_y[-1], grid_z[-1]))
+ dimension = np.array((nx.value - 1, ny.value - 1, nz.value - 1))
+ width = np.zeros(list(dimension) + [3])
+
+ for i, diff_x in enumerate(np.diff(grid_x)):
+ for j, diff_y in enumerate(np.diff(grid_y)):
+ for k, diff_z in enumerate(np.diff(grid_z)):
+ width[i, j, k, :] = diff_x, diff_y, diff_z
+
+ return (lower_left, upper_right, dimension, width)
+
+ def set_grid(self, x_grid, y_grid, z_grid):
+ nx = len(x_grid)
+ x_grid = (c_double*nx)(*x_grid)
+ ny = len(y_grid)
+ y_grid = (c_double*ny)(*y_grid)
+ nz = len(z_grid)
+ z_grid = (c_double*nz)(*z_grid)
+ _dll.openmc_cylindrical_mesh_set_grid(self._index, x_grid, nx, y_grid,
+ ny, z_grid, nz)
+
+class SphericalMesh(Mesh):
+ """SphericalMesh stored internally.
+
+ This class exposes a mesh that is stored internally in the OpenMC
+ library. To obtain a view of a mesh with a given ID, use the
+ :data:`openmc.lib.meshes` mapping.
+
+ Parameters
+ ----------
+ index : int
+ Index in the `meshes` array.
+
+ Attributes
+ ----------
+ id : int
+ ID of the mesh
+ dimension : iterable of int
+ The number of mesh cells in each direction.
+ lower_left : numpy.ndarray
+ The lower-left corner of the structured mesh.
+ upper_right : numpy.ndarray
+ The upper-right corner of the structrued mesh.
+ width : numpy.ndarray
+ The width of mesh cells in each direction.
+
+ """
+ mesh_type = 'spherical'
+
+ def __init__(self, uid=None, new=True, index=None):
+ super().__init__(uid, new, index)
+
+ @property
+ def lower_left(self):
+ return self._get_parameters()[0]
+
+ @property
+ def upper_right(self):
+ return self._get_parameters()[1]
+
+ @property
+ def dimension(self):
+ return self._get_parameters()[2]
+
+ @property
+ def width(self):
+ return self._get_parameters()[3]
+
+ def _get_parameters(self):
+ gx = POINTER(c_double)()
+ nx = c_int()
+ gy = POINTER(c_double)()
+ ny = c_int()
+ gz = POINTER(c_double)()
+ nz = c_int()
+ # Call C API to get grid parameters
+ _dll.openmc_spherical_mesh_get_grid(self._index, gx, nx, gy, ny, gz,
+ nz)
+
+ # Convert grid parameters to Numpy arrays
+ grid_x = as_array(gx, (nx.value,))
+ grid_y = as_array(gy, (ny.value,))
+ grid_z = as_array(gz, (nz.value,))
+
+ # Calculate lower_left, upper_right, width, and dimension from grid
+ lower_left = np.array((grid_x[0], grid_y[0], grid_z[0]))
+ upper_right = np.array((grid_x[-1], grid_y[-1], grid_z[-1]))
+ dimension = np.array((nx.value - 1, ny.value - 1, nz.value - 1))
+ width = np.zeros(list(dimension) + [3])
+
+ for i, diff_x in enumerate(np.diff(grid_x)):
+ for j, diff_y in enumerate(np.diff(grid_y)):
+ for k, diff_z in enumerate(np.diff(grid_z)):
+ width[i, j, k, :] = diff_x, diff_y, diff_z
+
+ return (lower_left, upper_right, dimension, width)
+
+ def set_grid(self, x_grid, y_grid, z_grid):
+ nx = len(x_grid)
+ x_grid = (c_double*nx)(*x_grid)
+ ny = len(y_grid)
+ y_grid = (c_double*ny)(*y_grid)
+ nz = len(z_grid)
+ z_grid = (c_double*nz)(*z_grid)
+ _dll.openmc_spherical_mesh_set_grid(self._index, x_grid, nx, y_grid,
+ ny, z_grid, nz)
+
+
_MESH_TYPE_MAP = {
'regular': RegularMesh,
- 'rectilinear': RectilinearMesh
+ 'rectilinear': RectilinearMesh,
+ 'cylindrical': CylindricalMesh,
+ 'spherical': SphericalMesh
}
diff --git a/openmc/mesh.py b/openmc/mesh.py
index 9aacdac283..3cd208e89d 100644
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@@ -79,6 +79,10 @@ class MeshBase(IDManagerMixin, ABC):
return RegularMesh.from_hdf5(group)
elif mesh_type == 'rectilinear':
return RectilinearMesh.from_hdf5(group)
+ elif mesh_type == 'cylindrical':
+ return CylindricalMesh.from_hdf5(group)
+ elif mesh_type == 'spherical':
+ return SphericalMesh.from_hdf5(group)
elif mesh_type == 'unstructured':
return UnstructuredMesh.from_hdf5(group)
else:
@@ -624,6 +628,334 @@ class RectilinearMesh(MeshBase):
return element
+class CylindricalMesh(MeshBase):
+ """A 3D cylindrical 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
+ dimension : Iterable of int
+ The number of mesh cells in each direction.
+ n_dimension : int
+ Number of mesh dimensions (always 3 for a CylindricalMesh).
+ r_grid : Iterable of float
+ Mesh boundary points along the r-axis.
+ Requirement is r >= 0.
+ phi_grid : Iterable of float
+ Mesh boundary points along the phi-axis.
+ The default value is [0, 360], i.e. the full phi range.
+ z_grid : Iterable of float
+ Mesh boundary points along the z-axis.
+ indices : Iterable of tuple
+ An iterable 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._r_grid = None
+ self._phi_grid = [0, 360]
+ self._z_grid = None
+
+ @property
+ def dimension(self):
+ return (len(self.r_grid) - 1,
+ len(self.phi_grid) - 1,
+ len(self.z_grid) - 1)
+
+ @property
+ def n_dimension(self):
+ return 3
+
+ @property
+ def r_grid(self):
+ return self._r_grid
+
+ @property
+ def phi_grid(self):
+ return self._phi_grid
+
+ @property
+ def z_grid(self):
+ return self._z_grid
+
+ @property
+ def indices(self):
+ nr, np, nz = self.dimension
+ np = len(self.phi_grid) - 1
+ nz = len(self.z_grid) - 1
+ return ((r, p, z)
+ for z in range(1, nz + 1)
+ for p in range(1, np + 1)
+ for r in range(1, nr + 1))
+
+ @r_grid.setter
+ def r_grid(self, grid):
+ cv.check_type('mesh r_grid', grid, Iterable, Real)
+ self._r_grid = grid
+
+ @phi_grid.setter
+ def phi_grid(self, grid):
+ cv.check_type('mesh phi_grid', grid, Iterable, Real)
+ self._phi_grid = np.array(grid)
+
+ @z_grid.setter
+ def z_grid(self, grid):
+ cv.check_type('mesh z_grid', grid, Iterable, Real)
+ self._z_grid = np.array(grid)
+
+ def __repr__(self):
+ fmt = '{0: <16}{1}{2}\n'
+ string = super().__repr__()
+ string += fmt.format('\tDimensions', '=\t', self.n_dimension)
+ r_grid_str = str(self._r_grid) if self._r_grid is None else len(self._r_grid)
+ string += fmt.format('\tN R pnts:', '=\t', r_grid_str)
+ if self._r_grid is not None:
+ string += fmt.format('\tR Min:', '=\t', self._r_grid[0])
+ string += fmt.format('\tR Max:', '=\t', self._r_grid[-1])
+ phi_grid_str = str(self._phi_grid) if self._phi_grid is None else len(self._phi_grid)
+ string += fmt.format('\tN Phi pnts:', '=\t', phi_grid_str)
+ if self._phi_grid is not None:
+ string += fmt.format('\tPhi Min:', '=\t', self._phi_grid[0])
+ string += fmt.format('\tPhi Max:', '=\t', self._phi_grid[-1])
+ z_grid_str = str(self._z_grid) if self._z_grid is None else len(self._z_grid)
+ string += fmt.format('\tN Z pnts:', '=\t', z_grid_str)
+ if self._z_grid is not None:
+ string += fmt.format('\tZ Min:', '=\t', self._z_grid[0])
+ string += fmt.format('\tZ Max:', '=\t', self._z_grid[-1])
+ return string
+
+ @classmethod
+ def from_hdf5(cls, group):
+ mesh_id = int(group.name.split('/')[-1].lstrip('mesh '))
+
+ # Read and assign mesh properties
+ mesh = cls(mesh_id)
+ mesh.r_grid = group['r_grid'][()]
+ mesh.phi_grid = 180 / np.pi * group['phi_grid'][()]
+ mesh.z_grid = group['z_grid'][()]
+
+ return mesh
+
+ def to_xml_element(self):
+ """Return XML representation of the mesh
+
+ Returns
+ -------
+ element : xml.etree.ElementTree.Element
+ XML element containing mesh data
+
+ """
+
+ element = ET.Element("mesh")
+ element.set("id", str(self._id))
+ element.set("type", "cylindrical")
+
+ subelement = ET.SubElement(element, "r_grid")
+ subelement.text = ' '.join(map(str, self.r_grid))
+
+ subelement = ET.SubElement(element, "phi_grid")
+ subelement.text = ' '.join(map(str, self.phi_grid))
+
+ subelement = ET.SubElement(element, "z_grid")
+ subelement.text = ' '.join(map(str, self.z_grid))
+
+ return element
+
+ def calc_mesh_volumes(self):
+ """Return Volumes for every mesh cell
+
+ Returns
+ -------
+ volumes : Iterable of float
+ Volumes
+
+ """
+
+ V_r = np.diff(np.array(self.r_grid)**2 / 2)
+ V_p = np.diff(np.array(self.phi_grid) * np.pi / 180.0)
+ V_z = np.diff(np.array(self.z_grid))
+
+ return np.multiply.outer(np.outer(V_r, V_p), V_z)
+
+
+class SphericalMesh(MeshBase):
+ """A 3D spherical 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
+ dimension : Iterable of int
+ The number of mesh cells in each direction.
+ n_dimension : int
+ Number of mesh dimensions (always 3 for a SphericalMesh).
+ r_grid : Iterable of float
+ Mesh boundary points along the r-axis.
+ Requirement is r >= 0.
+ theta_grid : Iterable of float
+ Mesh boundary points along the theta-axis in degrees.
+ The default value is [0, 180], i.e. the full theta range.
+ phi_grid : Iterable of float
+ Mesh boundary points along the phi-axis in degrees.
+ The default value is [0, 360], i.e. the full phi range.
+ indices : Iterable of tuple
+ An iterable 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._r_grid = None
+ self._theta_grid = [0, 180]
+ self._phi_grid = [0, 360]
+
+ @property
+ def dimension(self):
+ return (len(self.r_grid) - 1,
+ len(self.theta_grid) - 1,
+ len(self.phi_grid) - 1)
+
+ @property
+ def n_dimension(self):
+ return 3
+
+ @property
+ def r_grid(self):
+ return self._r_grid
+
+ @property
+ def theta_grid(self):
+ return self._theta_grid
+
+ @property
+ def phi_grid(self):
+ return self._phi_grid
+
+ @property
+ def indices(self):
+ nr, nt, np = self.dimension
+ nt = len(self.theta_grid) - 1
+ np = len(self.phi_grid) - 1
+ return ((r, t, p)
+ for p in range(1, np + 1)
+ for t in range(1, nt + 1)
+ for r in range(1, nr + 1))
+
+ @r_grid.setter
+ def r_grid(self, grid):
+ cv.check_type('mesh r_grid', grid, Iterable, Real)
+ self._r_grid = grid
+
+ @theta_grid.setter
+ def theta_grid(self, grid):
+ cv.check_type('mesh theta_grid', grid, Iterable, Real)
+ self._theta_grid = np.array(grid)
+
+ @phi_grid.setter
+ def phi_grid(self, grid):
+ cv.check_type('mesh phi_grid', grid, Iterable, Real)
+ self._phi_grid = np.array(grid)
+
+ def __repr__(self):
+ fmt = '{0: <16}{1}{2}\n'
+ string = super().__repr__()
+ string += fmt.format('\tDimensions', '=\t', self.n_dimension)
+ r_grid_str = str(self._r_grid) if not self._r_grid else len(self._r_grid)
+ string += fmt.format('\tN R pnts:', '=\t', r_grid_str)
+ if self._r_grid:
+ string += fmt.format('\tR Min:', '=\t', self._r_grid[0])
+ string += fmt.format('\tR Max:', '=\t', self._r_grid[-1])
+ theta_grid_str = str(self._theta_grid) if not self._theta_grid else len(self._theta_grid)
+ string += fmt.format('\tN Theta pnts:', '=\t', theta_grid_str)
+ if self._theta_grid:
+ string += fmt.format('\tTheta Min:', '=\t', self._theta_grid[0])
+ string += fmt.format('\tTheta Max:', '=\t', self._theta_grid[-1])
+ phi_grid_str = str(self._phi_grid) if not self._phi_grid else len(self._phi_grid)
+ string += fmt.format('\tN Phi pnts:', '=\t', phi_grid_str)
+ if self._phi_grid:
+ string += fmt.format('\tPhi Min:', '=\t', self._phi_grid[0])
+ string += fmt.format('\tPhi Max:', '=\t', self._phi_grid[-1])
+ return string
+
+ @classmethod
+ def from_hdf5(cls, group):
+ mesh_id = int(group.name.split('/')[-1].lstrip('mesh '))
+
+ # Read and assign mesh properties
+ mesh = cls(mesh_id)
+ mesh.r_grid = group['r_grid'][()]
+ mesh.theta_grid = 180 / np.pi * group['theta_grid'][()]
+ mesh.phi_grid = 180 / np.pi * group['phi_grid'][()]
+
+ return mesh
+
+ def to_xml_element(self):
+ """Return XML representation of the mesh
+
+ Returns
+ -------
+ element : xml.etree.ElementTree.Element
+ XML element containing mesh data
+
+ """
+
+ element = ET.Element("mesh")
+ element.set("id", str(self._id))
+ element.set("type", "spherical")
+
+ subelement = ET.SubElement(element, "r_grid")
+ subelement.text = ' '.join(map(str, self.r_grid))
+
+ subelement = ET.SubElement(element, "theta_grid")
+ subelement.text = ' '.join(map(str, self.theta_grid))
+
+ subelement = ET.SubElement(element, "phi_grid")
+ subelement.text = ' '.join(map(str, self.phi_grid))
+
+ return element
+
+ def calc_mesh_volumes(self):
+ """Return Volumes for every mesh cell
+
+ Returns
+ -------
+ volumes : Iterable of float
+ Volumes
+
+ """
+
+ V_r = np.diff(np.array(self.r_grid)**3 / 3)
+ V_t = np.diff(-np.cos(np.pi * np.array(self.theta_grid) / 180.0))
+ V_p = np.diff(np.array(self.phi_grid) * np.pi / 180)
+
+ return np.multiply.outer(np.outer(V_r, V_t), V_p)
+
+
+
class UnstructuredMesh(MeshBase):
"""A 3D unstructured mesh
diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp
index b35281f4ad..5584210a1f 100644
--- a/src/eigenvalue.cpp
+++ b/src/eigenvalue.cpp
@@ -577,7 +577,7 @@ void ufs_count_sites()
#ifdef OPENMC_MPI
// Send source fraction to all processors
- int n_bins = xt::prod(simulation::ufs_mesh->shape_)();
+ int n_bins = simulation::ufs_mesh->n_bins();
MPI_Bcast(
simulation::source_frac.data(), n_bins, MPI_DOUBLE, 0, mpi::intracomm);
#endif
diff --git a/src/mesh.cpp b/src/mesh.cpp
index cabdd6b3e3..6654684acf 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -33,9 +33,9 @@
#include "openmc/xml_interface.h"
#ifdef LIBMESH
+#include "libmesh/mesh_modification.h"
#include "libmesh/mesh_tools.h"
#include "libmesh/numeric_vector.h"
-#include "libmesh/mesh_modification.h"
#endif
namespace openmc {
@@ -145,8 +145,7 @@ void Mesh::set_id(int32_t id)
std::string StructuredMesh::bin_label(int bin) const
{
- vector ijk(n_dimension_);
- get_indices_from_bin(bin, ijk.data());
+ MeshIndex ijk = get_indices_from_bin(bin);
if (n_dimension_ > 2) {
return fmt::format("Mesh Index ({}, {}, {})", ijk[0], ijk[1], ijk[2]);
@@ -157,6 +156,11 @@ std::string StructuredMesh::bin_label(int bin) const
}
}
+xt::xtensor StructuredMesh::get_x_shape() const
+{
+ return xt::adapt(shape_, {n_dimension_});
+}
+
//==============================================================================
// Unstructured Mesh implementation
//==============================================================================
@@ -167,7 +171,7 @@ UnstructuredMesh::UnstructuredMesh(pugi::xml_node node) : Mesh(node)
// check the mesh type
if (check_for_node(node, "type")) {
auto temp = get_node_value(node, "type", true, true);
- if (temp != "unstructured") {
+ if (temp != mesh_type) {
fatal_error(fmt::format("Invalid mesh type: {}", temp));
}
}
@@ -196,6 +200,13 @@ UnstructuredMesh::UnstructuredMesh(pugi::xml_node node) : Mesh(node)
}
}
+const std::string UnstructuredMesh::mesh_type = "unstructured";
+
+std::string UnstructuredMesh::get_mesh_type() const
+{
+ return mesh_type;
+}
+
void UnstructuredMesh::surface_bins_crossed(
Position r0, Position r1, const Direction& u, vector& bins) const
{
@@ -211,7 +222,7 @@ void UnstructuredMesh::to_hdf5(hid_t group) const
{
hid_t mesh_group = create_group(group, fmt::format("mesh {}", id_));
- write_dataset(mesh_group, "type", "unstructured");
+ write_dataset(mesh_group, "type", mesh_type);
write_dataset(mesh_group, "filename", filename_);
write_dataset(mesh_group, "library", this->library());
// write volume of each element
@@ -240,18 +251,21 @@ void UnstructuredMesh::set_length_multiplier(double length_multiplier)
specified_length_multiplier_ = true;
}
-void StructuredMesh::get_indices(Position r, int* ijk, bool* in_mesh) const
+StructuredMesh::MeshIndex StructuredMesh::get_indices(
+ Position r, bool& in_mesh) const
{
- *in_mesh = true;
+ MeshIndex ijk;
+ in_mesh = true;
for (int i = 0; i < n_dimension_; ++i) {
ijk[i] = get_index_in_direction(r[i], i);
if (ijk[i] < 1 || ijk[i] > shape_[i])
- *in_mesh = false;
+ in_mesh = false;
}
+ return ijk;
}
-int StructuredMesh::get_bin_from_indices(const int* ijk) const
+int StructuredMesh::get_bin_from_indices(const MeshIndex& ijk) const
{
switch (n_dimension_) {
case 1:
@@ -265,8 +279,9 @@ int StructuredMesh::get_bin_from_indices(const int* ijk) const
}
}
-void StructuredMesh::get_indices_from_bin(int bin, int* ijk) const
+StructuredMesh::MeshIndex StructuredMesh::get_indices_from_bin(int bin) const
{
+ MeshIndex ijk;
if (n_dimension_ == 1) {
ijk[0] = bin + 1;
} else if (n_dimension_ == 2) {
@@ -277,24 +292,25 @@ void StructuredMesh::get_indices_from_bin(int bin, int* ijk) const
ijk[1] = (bin % (shape_[0] * shape_[1])) / shape_[0] + 1;
ijk[2] = bin / (shape_[0] * shape_[1]) + 1;
}
+ return ijk;
}
int StructuredMesh::get_bin(Position r) const
{
// Determine indices
- vector ijk(n_dimension_);
bool in_mesh;
- get_indices(r, ijk.data(), &in_mesh);
+ MeshIndex ijk = get_indices(r, in_mesh);
if (!in_mesh)
return -1;
// Convert indices to bin
- return get_bin_from_indices(ijk.data());
+ return get_bin_from_indices(ijk);
}
int StructuredMesh::n_bins() const
{
- return xt::prod(shape_)();
+ return std::accumulate(
+ shape_.begin(), shape_.begin() + n_dimension_, 1, std::multiplies<>());
}
int StructuredMesh::n_surface_bins() const
@@ -357,360 +373,178 @@ xt::xtensor StructuredMesh::count_sites(
return counts;
}
-bool StructuredMesh::intersects(Position& r0, Position r1, int* ijk) const
+// raytrace through the mesh. The template class T will do the tallying.
+// A modern optimizing compiler can recognize the noop method of T and eleminate
+// that call entirely.
+template
+void StructuredMesh::raytrace_mesh(
+ Position r0, Position r1, const Direction& u, T tally) const
{
- switch (n_dimension_) {
- case 1:
- return intersects_1d(r0, r1, ijk);
- case 2:
- return intersects_2d(r0, r1, ijk);
- case 3:
- return intersects_3d(r0, r1, ijk);
- default:
- throw std::runtime_error {"Invalid number of mesh dimensions."};
- }
-}
-bool StructuredMesh::intersects_1d(Position& r0, Position r1, int* ijk) const
-{
- // Copy coordinates of starting point
- double x0 = r0.x;
- double y0 = r0.y;
- double z0 = r0.z;
-
- // Copy coordinates of ending point
- double x1 = r1.x;
- double y1 = r1.y;
- double z1 = r1.z;
-
- // Copy coordinates of mesh lower_left and upper_right
- double xm0 = lower_left_[0];
- double xm1 = upper_right_[0];
-
- double min_dist = INFTY;
-
- // Check if line intersects left surface -- calculate the intersection point
- // (y,z)
- if ((x0 < xm0 && x1 > xm0) || (x0 > xm0 && x1 < xm0)) {
- double yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0);
- double zi = z0 + (xm0 - x0) * (z1 - z0) / (x1 - x0);
- if (check_intersection_point(xm0, x0, yi, yi, zi, zi, r0, min_dist)) {
- ijk[0] = 1;
- }
- }
-
- // Check if line intersects right surface -- calculate the intersection point
- // (y,z)
- if ((x0 < xm1 && x1 > xm1) || (x0 > xm1 && x1 < xm1)) {
- double yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0);
- double zi = z0 + (xm1 - x0) * (z1 - z0) / (x1 - x0);
- if (check_intersection_point(xm1, x0, yi, yi, zi, zi, r0, min_dist)) {
- ijk[0] = shape_[0];
- }
- }
-
- return min_dist < INFTY;
-}
-
-bool StructuredMesh::intersects_2d(Position& r0, Position r1, int* ijk) const
-{
- // Copy coordinates of starting point
- double x0 = r0.x;
- double y0 = r0.y;
- double z0 = r0.z;
-
- // Copy coordinates of ending point
- double x1 = r1.x;
- double y1 = r1.y;
- double z1 = r1.z;
-
- // Copy coordinates of mesh lower_left
- double xm0 = lower_left_[0];
- double ym0 = lower_left_[1];
-
- // Copy coordinates of mesh upper_right
- double xm1 = upper_right_[0];
- double ym1 = upper_right_[1];
-
- double min_dist = INFTY;
-
- // Check if line intersects left surface -- calculate the intersection point
- // (y,z)
- if ((x0 < xm0 && x1 > xm0) || (x0 > xm0 && x1 < xm0)) {
- double yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0);
- double zi = z0 + (xm0 - x0) * (z1 - z0) / (x1 - x0);
- if (yi >= ym0 && yi < ym1) {
- if (check_intersection_point(xm0, x0, yi, y0, zi, zi, r0, min_dist)) {
- ijk[0] = 1;
- ijk[1] = get_index_in_direction(yi, 1);
- }
- }
- }
-
- // Check if line intersects back surface -- calculate the intersection point
- // (x,z)
- if ((y0 < ym0 && y1 > ym0) || (y0 > ym0 && y1 < ym0)) {
- double xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0);
- double zi = z0 + (ym0 - y0) * (z1 - z0) / (y1 - y0);
- if (xi >= xm0 && xi < xm1) {
- if (check_intersection_point(xi, x0, ym0, y0, zi, zi, r0, min_dist)) {
- ijk[0] = get_index_in_direction(xi, 0);
- ijk[1] = 1;
- }
- }
- }
-
- // Check if line intersects right surface -- calculate the intersection point
- // (y,z)
- if ((x0 < xm1 && x1 > xm1) || (x0 > xm1 && x1 < xm1)) {
- double yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0);
- double zi = z0 + (xm1 - x0) * (z1 - z0) / (x1 - x0);
- if (yi >= ym0 && yi < ym1) {
- if (check_intersection_point(xm1, x0, yi, y0, zi, zi, r0, min_dist)) {
- ijk[0] = shape_[0];
- ijk[1] = get_index_in_direction(yi, 1);
- }
- }
- }
-
- // Check if line intersects front surface -- calculate the intersection point
- // (x,z)
- if ((y0 < ym1 && y1 > ym1) || (y0 > ym1 && y1 < ym1)) {
- double xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0);
- double zi = z0 + (ym1 - y0) * (z1 - z0) / (y1 - y0);
- if (xi >= xm0 && xi < xm1) {
- if (check_intersection_point(xi, x0, ym1, y0, zi, zi, r0, min_dist)) {
- ijk[0] = get_index_in_direction(xi, 0);
- ijk[1] = shape_[1];
- }
- }
- }
-
- return min_dist < INFTY;
-}
-
-bool StructuredMesh::intersects_3d(Position& r0, Position r1, int* ijk) const
-{
- // Copy coordinates of starting point
- double x0 = r0.x;
- double y0 = r0.y;
- double z0 = r0.z;
-
- // Copy coordinates of ending point
- double x1 = r1.x;
- double y1 = r1.y;
- double z1 = r1.z;
-
- // Copy coordinates of mesh lower_left
- double xm0 = lower_left_[0];
- double ym0 = lower_left_[1];
- double zm0 = lower_left_[2];
-
- // Copy coordinates of mesh upper_right
- double xm1 = upper_right_[0];
- double ym1 = upper_right_[1];
- double zm1 = upper_right_[2];
-
- double min_dist = INFTY;
-
- // Check if line intersects left surface -- calculate the intersection point
- // (y,z)
- if ((x0 < xm0 && x1 > xm0) || (x0 > xm0 && x1 < xm0)) {
- double yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0);
- double zi = z0 + (xm0 - x0) * (z1 - z0) / (x1 - x0);
- if (yi >= ym0 && yi < ym1 && zi >= zm0 && zi < zm1) {
- if (check_intersection_point(xm0, x0, yi, y0, zi, z0, r0, min_dist)) {
- ijk[0] = 1;
- ijk[1] = get_index_in_direction(yi, 1);
- ijk[2] = get_index_in_direction(zi, 2);
- }
- }
- }
-
- // Check if line intersects back surface -- calculate the intersection point
- // (x,z)
- if ((y0 < ym0 && y1 > ym0) || (y0 > ym0 && y1 < ym0)) {
- double xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0);
- double zi = z0 + (ym0 - y0) * (z1 - z0) / (y1 - y0);
- if (xi >= xm0 && xi < xm1 && zi >= zm0 && zi < zm1) {
- if (check_intersection_point(xi, x0, ym0, y0, zi, z0, r0, min_dist)) {
- ijk[0] = get_index_in_direction(xi, 0);
- ijk[1] = 1;
- ijk[2] = get_index_in_direction(zi, 2);
- }
- }
- }
-
- // Check if line intersects bottom surface -- calculate the intersection
- // point (x,y)
- if ((z0 < zm0 && z1 > zm0) || (z0 > zm0 && z1 < zm0)) {
- double xi = x0 + (zm0 - z0) * (x1 - x0) / (z1 - z0);
- double yi = y0 + (zm0 - z0) * (y1 - y0) / (z1 - z0);
- if (xi >= xm0 && xi < xm1 && yi >= ym0 && yi < ym1) {
- if (check_intersection_point(xi, x0, yi, y0, zm0, z0, r0, min_dist)) {
- ijk[0] = get_index_in_direction(xi, 0);
- ijk[1] = get_index_in_direction(yi, 1);
- ijk[2] = 1;
- }
- }
- }
-
- // Check if line intersects right surface -- calculate the intersection point
- // (y,z)
- if ((x0 < xm1 && x1 > xm1) || (x0 > xm1 && x1 < xm1)) {
- double yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0);
- double zi = z0 + (xm1 - x0) * (z1 - z0) / (x1 - x0);
- if (yi >= ym0 && yi < ym1 && zi >= zm0 && zi < zm1) {
- if (check_intersection_point(xm1, x0, yi, y0, zi, z0, r0, min_dist)) {
- ijk[0] = shape_[0];
- ijk[1] = get_index_in_direction(yi, 1);
- ijk[2] = get_index_in_direction(zi, 2);
- }
- }
- }
-
- // Check if line intersects front surface -- calculate the intersection point
- // (x,z)
- if ((y0 < ym1 && y1 > ym1) || (y0 > ym1 && y1 < ym1)) {
- double xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0);
- double zi = z0 + (ym1 - y0) * (z1 - z0) / (y1 - y0);
- if (xi >= xm0 && xi < xm1 && zi >= zm0 && zi < zm1) {
- if (check_intersection_point(xi, x0, ym1, y0, zi, z0, r0, min_dist)) {
- ijk[0] = get_index_in_direction(xi, 0);
- ijk[1] = shape_[1];
- ijk[2] = get_index_in_direction(zi, 2);
- }
- }
- }
-
- // Check if line intersects top surface -- calculate the intersection point
- // (x,y)
- if ((z0 < zm1 && z1 > zm1) || (z0 > zm1 && z1 < zm1)) {
- double xi = x0 + (zm1 - z0) * (x1 - x0) / (z1 - z0);
- double yi = y0 + (zm1 - z0) * (y1 - y0) / (z1 - z0);
- if (xi >= xm0 && xi < xm1 && yi >= ym0 && yi < ym1) {
- if (check_intersection_point(xi, x0, yi, y0, zm1, z0, r0, min_dist)) {
- ijk[0] = get_index_in_direction(xi, 0);
- ijk[1] = get_index_in_direction(yi, 1);
- ijk[2] = shape_[2];
- }
- }
- }
-
- return min_dist < INFTY;
-}
-
-void StructuredMesh::bins_crossed(Position r0, Position r1, const Direction& u,
- vector& bins, vector& lengths) const
-{
- // ========================================================================
- // Determine where the track intersects the mesh and if it intersects at all.
+ // TODO: when c++-17 is available, use "if constexpr ()" to compile-time
+ // enable/disable tally calls for now, T template type needs to provide both
+ // surface and track methods, which might be empty. modern optimizing
+ // compilers will (hopefully) eleminate the complete code (including
+ // calculation of parameters) but for the future: be explicit
// Compute the length of the entire track.
double total_distance = (r1 - r0).norm();
if (total_distance == 0.0)
return;
- // While determining if this track intersects the mesh, offset the starting
- // and ending coords by a bit. This avoid finite-precision errors that can
- // occur when the mesh surfaces coincide with lattice or geometric surfaces.
- Position last_r = r0 + TINY_BIT * u;
- Position r = r1 - TINY_BIT * u;
+ const int n = n_dimension_;
- // Determine the mesh indices for the starting and ending coords. Here, we
- // use arrays for ijk0 and ijk1 instead of vector because we obtain a
- // small performance improvement by forcing this data to live on the stack,
- // rather than on the heap. We know the maximum length is 3, and by
- // ensuring that all loops are only indexed up to n_dimension, we will not
- // access any non-initialized values. The same concept is used throughout.
- int n = n_dimension_;
- int ijk0[3], ijk1[3];
- bool start_in_mesh;
- get_indices(last_r, ijk0, &start_in_mesh);
- bool end_in_mesh;
- get_indices(r, ijk1, &end_in_mesh);
+ // Flag if position is inside the mesh
+ bool in_mesh;
- // Reset coordinates and check for a mesh intersection if necessary.
- if (start_in_mesh) {
- // The initial coords lie in the mesh, use those coords for tallying.
- last_r = r0;
- } else {
- // The initial coords do not lie in the mesh. Check to see if the particle
- // eventually intersects the mesh and compute the relevant coords and
- // indices.
- if (!intersects(last_r, r, ijk0))
- return;
- }
- r = r1;
+ // Position is r = r0 + u * traveled_distance, start at r0
+ double traveled_distance {0.0};
- // The TINY_BIT offsets above mean that the preceding logic cannot always find
- // the correct ijk0 and ijk1 indices. For tracks shorter than 2*TINY_BIT, just
- // assume the track lies in only one mesh bin. These tracks are very short so
- // any error caused by this assumption will be small. It is important that
- // ijk0 values are used rather than ijk1 because the previous logic guarantees
- // ijk0 is a valid mesh bin.
+ // Calculate index of current cell. Offset the position a tiny bit in
+ // direction of flight
+ MeshIndex ijk = get_indices(r0 + TINY_BIT * u, in_mesh);
+
+ // if track is very short, assume that it is completely inside one cell.
+ // Only the current cell will score and no surfaces
if (total_distance < 2 * TINY_BIT) {
- for (int i = 0; i < n; ++i)
- ijk1[i] = ijk0[i];
+ if (in_mesh) {
+ tally.track(ijk, 1.0);
+ }
+ return;
}
- // ========================================================================
- // Find which mesh cells are traversed and the length of each traversal.
+ // Calculate initial distances to next surfaces in all three dimensions
+ std::array distances;
+ for (int k = 0; k < n; ++k) {
+ distances[k] = distance_to_grid_boundary(ijk, k, r0, u, 0.0);
+ }
+ // Loop until r = r1 is eventually reached
while (true) {
- if (std::equal(ijk0, ijk0 + n, ijk1)) {
- // The track ends in this cell. Use the particle end location rather
- // than the mesh surface and stop iterating.
- double distance = (last_r - r).norm();
- bins.push_back(get_bin_from_indices(ijk0));
- lengths.push_back(distance / total_distance);
- break;
- }
- // The track exits this cell. Determine the distance to each mesh surface.
- double d[3];
- for (int k = 0; k < n; ++k) {
- if (std::fabs(u[k]) < FP_PRECISION) {
- d[k] = INFTY;
- } else if (u[k] > 0) {
- double xyz_cross = positive_grid_boundary(ijk0, k);
- d[k] = (xyz_cross - last_r[k]) / u[k];
- } else {
- double xyz_cross = negative_grid_boundary(ijk0, k);
- d[k] = (xyz_cross - last_r[k]) / u[k];
+ if (in_mesh) {
+
+ // find surface with minimal distance to current position
+ const auto k = std::min_element(distances.begin(), distances.end()) -
+ distances.begin();
+
+ // Tally track length delta since last step
+ tally.track(ijk,
+ (std::min(distances[k].distance, total_distance) - traveled_distance) /
+ total_distance);
+
+ // update position and leave, if we have reached end position
+ traveled_distance = distances[k].distance;
+ if (traveled_distance >= total_distance)
+ return;
+
+ // If we have not reached r1, we have hit a surface. Tally outward current
+ tally.surface(ijk, k, distances[k].max_surface, false);
+
+ // Update cell and calculate distance to next surface in k-direction.
+ // The two other directions are still valid!
+ ijk[k] = distances[k].next_index;
+ distances[k] =
+ distance_to_grid_boundary(ijk, k, r0, u, traveled_distance);
+
+ // Check if we have left the interior of the mesh
+ in_mesh = ((ijk[k] >= 1) && (ijk[k] <= shape_[k]));
+
+ // If we are still inside the mesh, tally inward current for the next cell
+ if (in_mesh)
+ tally.surface(ijk, k, !distances[k].max_surface, true);
+
+ } else { // not inside mesh
+
+ // For all directions outside the mesh, find the distance that we need to
+ // travel to reach the next surface. Use the largest distance, as only
+ // this will cross all outer surfaces.
+ int k_max {0};
+ for (int k = 0; k < n; ++k) {
+ if ((ijk[k] < 1 || ijk[k] > shape_[k]) &&
+ (distances[k].distance > traveled_distance)) {
+ traveled_distance = distances[k].distance;
+ k_max = k;
+ }
}
- }
- // Pick the closest mesh surface and append this traversal to the output.
- auto j = std::min_element(d, d + n) - d;
- double distance = d[j];
- bins.push_back(get_bin_from_indices(ijk0));
- lengths.push_back(distance / total_distance);
+ // If r1 is not inside the mesh, exit here
+ if (traveled_distance >= total_distance)
+ return;
- // Translate to the oncoming mesh surface.
- last_r += distance * u;
-
- // Increment the indices into the next mesh cell.
- if (u[j] > 0.0) {
- ++ijk0[j];
- } else {
- --ijk0[j];
- }
-
- // If the next indices are invalid, then the track has left the mesh and
- // we are done.
- bool in_mesh = true;
- for (int i = 0; i < n; ++i) {
- if (ijk0[i] < 1 || ijk0[i] > shape_[i]) {
- in_mesh = false;
- break;
+ // Calculate the new cell index and update all distances to next surfaces.
+ ijk = get_indices(r0 + (traveled_distance + TINY_BIT) * u, in_mesh);
+ for (int k = 0; k < n; ++k) {
+ distances[k] =
+ distance_to_grid_boundary(ijk, k, r0, u, traveled_distance);
}
+
+ // If inside the mesh, Tally inward current
+ if (in_mesh)
+ tally.surface(ijk, k_max, !distances[k_max].max_surface, true);
}
- if (!in_mesh)
- break;
}
}
+void StructuredMesh::bins_crossed(Position r0, Position r1, const Direction& u,
+ vector& bins, vector& lengths) const
+{
+
+ // Helper tally class.
+ // stores a pointer to the mesh class and references to bins and lengths
+ // parameters. Performs the actual tally through the track method.
+ struct TrackAggregator {
+ TrackAggregator(
+ const StructuredMesh* _mesh, vector& _bins, vector& _lengths)
+ : mesh(_mesh), bins(_bins), lengths(_lengths)
+ {}
+ void surface(const MeshIndex& ijk, int k, bool max, bool inward) const {}
+ void track(const MeshIndex& ijk, double l) const
+ {
+ bins.push_back(mesh->get_bin_from_indices(ijk));
+ lengths.push_back(l);
+ }
+
+ const StructuredMesh* mesh;
+ vector& bins;
+ vector& lengths;
+ };
+
+ // Perform the mesh raytrace with the helper class.
+ raytrace_mesh(r0, r1, u, TrackAggregator(this, bins, lengths));
+}
+
+void StructuredMesh::surface_bins_crossed(
+ Position r0, Position r1, const Direction& u, vector& bins) const
+{
+
+ // Helper tally class.
+ // stores a pointer to the mesh class and a reference to the bins parameter.
+ // Performs the actual tally through the surface method.
+ struct SurfaceAggregator {
+ SurfaceAggregator(const StructuredMesh* _mesh, vector& _bins)
+ : mesh(_mesh), bins(_bins)
+ {}
+ void surface(const MeshIndex& ijk, int k, bool max, bool inward) const
+ {
+ int i_bin =
+ 4 * mesh->n_dimension_ * mesh->get_bin_from_indices(ijk) + 4 * k;
+ if (max)
+ i_bin += 2;
+ if (inward)
+ i_bin += 1;
+ bins.push_back(i_bin);
+ }
+ void track(const MeshIndex& idx, double l) const {}
+
+ const StructuredMesh* mesh;
+ vector& bins;
+ };
+
+ // Perform the mesh raytrace with the helper class.
+ raytrace_mesh(r0, r1, u, SurfaceAggregator(this, bins));
+}
+
//==============================================================================
// RegularMesh implementation
//==============================================================================
@@ -722,14 +556,15 @@ RegularMesh::RegularMesh(pugi::xml_node node) : StructuredMesh {node}
fatal_error("Must specify on a regular mesh.");
}
- shape_ = get_node_xarray(node, "dimension");
- int n = n_dimension_ = shape_.size();
+ xt::xtensor shape = get_node_xarray(node, "dimension");
+ 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)) {
+ if (xt::any(shape <= 0)) {
fatal_error("All entries on the element for a tally "
"mesh must be positive.");
}
@@ -743,7 +578,7 @@ RegularMesh::RegularMesh(pugi::xml_node node) : StructuredMesh {node}
}
// Make sure lower_left and dimension match
- if (shape_.size() != lower_left_.size()) {
+ if (shape.size() != lower_left_.size()) {
fatal_error("Number of entries on must be the same "
"as the number of entries on .");
}
@@ -769,7 +604,7 @@ RegularMesh::RegularMesh(pugi::xml_node node) : StructuredMesh {node}
}
// Set width and upper right coordinate
- upper_right_ = xt::eval(lower_left_ + shape_ * width_);
+ upper_right_ = xt::eval(lower_left_ + shape * width_);
} else if (check_for_node(node, "upper_right")) {
upper_right_ = get_node_xarray(node, "upper_right");
@@ -788,13 +623,13 @@ RegularMesh::RegularMesh(pugi::xml_node node) : StructuredMesh {node}
}
// Set width
- width_ = xt::eval((upper_right_ - lower_left_) / shape_);
+ width_ = xt::eval((upper_right_ - lower_left_) / shape);
} else {
fatal_error("Must specify either or on a mesh.");
}
// Set volume fraction
- volume_frac_ = 1.0 / xt::prod(shape_)();
+ volume_frac_ = 1.0 / xt::prod(shape)();
}
int RegularMesh::get_index_in_direction(double r, int i) const
@@ -802,158 +637,41 @@ int RegularMesh::get_index_in_direction(double r, int i) const
return std::ceil((r - lower_left_[i]) / width_[i]);
}
-double RegularMesh::positive_grid_boundary(int* ijk, int i) const
+const std::string RegularMesh::mesh_type = "regular";
+
+std::string RegularMesh::get_mesh_type() const
+{
+ return mesh_type;
+}
+
+double RegularMesh::positive_grid_boundary(const MeshIndex& ijk, int i) const
{
return lower_left_[i] + ijk[i] * width_[i];
}
-double RegularMesh::negative_grid_boundary(int* ijk, int i) const
+double RegularMesh::negative_grid_boundary(const MeshIndex& ijk, int i) const
{
return lower_left_[i] + (ijk[i] - 1) * width_[i];
}
-void RegularMesh::surface_bins_crossed(
- Position r0, Position r1, const Direction& u, vector& bins) const
+StructuredMesh::MeshDistance RegularMesh::distance_to_grid_boundary(
+ const MeshIndex& ijk, int i, const Position& r0, const Direction& u,
+ double l) const
{
- // Determine indices for starting and ending location.
- int n = n_dimension_;
- vector ijk0(n), ijk1(n);
- bool start_in_mesh;
- get_indices(r0, ijk0.data(), &start_in_mesh);
- bool end_in_mesh;
- get_indices(r1, ijk1.data(), &end_in_mesh);
+ MeshDistance d;
+ d.next_index = ijk[i];
+ if (std::abs(u[i]) < FP_PRECISION)
+ return d;
- // Check if the track intersects any part of the mesh.
- if (!start_in_mesh) {
- Position r0_copy = r0;
- vector ijk0_copy(ijk0);
- if (!intersects(r0_copy, r1, ijk0_copy.data()))
- return;
- }
-
- // ========================================================================
- // Find which mesh surfaces are crossed.
-
- // Calculate number of surface crossings
- int n_cross = 0;
- for (int i = 0; i < n; ++i)
- n_cross += std::abs(ijk1[i] - ijk0[i]);
- if (n_cross == 0)
- return;
-
- // Bounding coordinates
- Position xyz_cross;
- for (int i = 0; i < n; ++i) {
- if (u[i] > 0.0) {
- xyz_cross[i] = positive_grid_boundary(ijk0.data(), i);
- } else {
- xyz_cross[i] = negative_grid_boundary(ijk0.data(), i);
- }
- }
-
- for (int j = 0; j < n_cross; ++j) {
- // Set the distances to infinity
- Position d {INFTY, INFTY, INFTY};
-
- // Determine closest bounding surface. We need to treat
- // special case where the cosine of the angle is zero since this would
- // result in a divide-by-zero.
- double distance = INFTY;
- for (int i = 0; i < n; ++i) {
- if (u[i] == 0) {
- d[i] = INFTY;
- } else {
- d[i] = (xyz_cross[i] - r0[i]) / u[i];
- }
- distance = std::min(distance, d[i]);
- }
-
- // Loop over the dimensions
- for (int i = 0; i < n; ++i) {
- // Check whether distance is the shortest distance
- if (distance == d[i]) {
-
- // Check whether the current indices are within the mesh bounds
- bool in_mesh = true;
- for (int j = 0; j < n; ++j) {
- if (ijk0[j] < 1 || ijk0[j] > shape_[j]) {
- in_mesh = false;
- break;
- }
- }
-
- // Check whether particle is moving in positive i direction
- if (u[i] > 0) {
-
- // Outward current on i max surface
- if (in_mesh) {
- int i_surf = 4 * i + 3;
- int i_mesh = get_bin_from_indices(ijk0.data());
- int i_bin = 4 * n * i_mesh + i_surf - 1;
-
- bins.push_back(i_bin);
- }
-
- // Advance position
- ++ijk0[i];
- xyz_cross[i] += width_[i];
- in_mesh = true;
- for (int j = 0; j < n; ++j) {
- if (ijk0[j] < 1 || ijk0[j] > shape_[j]) {
- in_mesh = false;
- break;
- }
- }
-
- // If the particle crossed the surface, tally the inward current on
- // i min surface
- if (in_mesh) {
- int i_surf = 4 * i + 2;
- int i_mesh = get_bin_from_indices(ijk0.data());
- int i_bin = 4 * n * i_mesh + i_surf - 1;
-
- bins.push_back(i_bin);
- }
-
- } else {
- // The particle is moving in the negative i direction
-
- // Outward current on i min surface
- if (in_mesh) {
- int i_surf = 4 * i + 1;
- int i_mesh = get_bin_from_indices(ijk0.data());
- int i_bin = 4 * n * i_mesh + i_surf - 1;
-
- bins.push_back(i_bin);
- }
-
- // Advance position
- --ijk0[i];
- xyz_cross[i] -= width_[i];
- in_mesh = true;
- for (int j = 0; j < n; ++j) {
- if (ijk0[j] < 1 || ijk0[j] > shape_[j]) {
- in_mesh = false;
- break;
- }
- }
-
- // If the particle crossed the surface, tally the inward current on
- // i max surface
- if (in_mesh) {
- int i_surf = 4 * i + 4;
- int i_mesh = get_bin_from_indices(ijk0.data());
- int i_bin = 4 * n * i_mesh + i_surf - 1;
-
- bins.push_back(i_bin);
- }
- }
- }
- }
-
- // Calculate new coordinates
- r0 += distance * u;
+ d.max_surface = (u[i] > 0);
+ if (d.max_surface && (ijk[i] <= shape_[i])) {
+ d.next_index++;
+ d.distance = (positive_grid_boundary(ijk, i) - r0[i]) / u[i];
+ } else if (!d.max_surface && (ijk[i] >= 1)) {
+ d.next_index--;
+ d.distance = (negative_grid_boundary(ijk, i) - r0[i]) / u[i];
}
+ return d;
}
std::pair, vector> RegularMesh::plot(
@@ -1002,7 +720,7 @@ void RegularMesh::to_hdf5(hid_t group) const
hid_t mesh_group = create_group(group, "mesh " + std::to_string(id_));
write_dataset(mesh_group, "type", "regular");
- write_dataset(mesh_group, "dimension", shape_);
+ write_dataset(mesh_group, "dimension", get_x_shape());
write_dataset(mesh_group, "lower_left", lower_left_);
write_dataset(mesh_group, "upper_right", upper_right_);
write_dataset(mesh_group, "width", width_);
@@ -1073,7 +791,6 @@ RectilinearMesh::RectilinearMesh(pugi::xml_node node) : StructuredMesh {node}
{
n_dimension_ = 3;
- grid_.resize(n_dimension_);
grid_[0] = get_node_array(node, "x_grid");
grid_[1] = get_node_array(node, "y_grid");
grid_[2] = get_node_array(node, "z_grid");
@@ -1083,16 +800,45 @@ RectilinearMesh::RectilinearMesh(pugi::xml_node node) : StructuredMesh {node}
}
}
-double RectilinearMesh::positive_grid_boundary(int* ijk, int i) const
+const std::string RectilinearMesh::mesh_type = "rectilinear";
+
+std::string RectilinearMesh::get_mesh_type() const
+{
+ return mesh_type;
+}
+
+double RectilinearMesh::positive_grid_boundary(
+ const MeshIndex& ijk, int i) const
{
return grid_[i][ijk[i]];
}
-double RectilinearMesh::negative_grid_boundary(int* ijk, int i) const
+double RectilinearMesh::negative_grid_boundary(
+ const MeshIndex& ijk, int i) const
{
return grid_[i][ijk[i] - 1];
}
+StructuredMesh::MeshDistance RectilinearMesh::distance_to_grid_boundary(
+ const MeshIndex& ijk, int i, const Position& r0, const Direction& u,
+ double l) const
+{
+ MeshDistance d;
+ d.next_index = ijk[i];
+ if (std::abs(u[i]) < FP_PRECISION)
+ return d;
+
+ d.max_surface = (u[i] > 0);
+ if (d.max_surface && (ijk[i] <= shape_[i])) {
+ d.next_index++;
+ d.distance = (positive_grid_boundary(ijk, i) - r0[i]) / u[i];
+ } else if (!d.max_surface && (ijk[i] > 0)) {
+ d.next_index--;
+ d.distance = (negative_grid_boundary(ijk, i) - r0[i]) / u[i];
+ }
+ return d;
+}
+
int RectilinearMesh::set_grid()
{
shape_ = {static_cast(grid_[0].size()) - 1,
@@ -1105,12 +851,11 @@ int RectilinearMesh::set_grid()
"must each have at least 2 points");
return OPENMC_E_INVALID_ARGUMENT;
}
- for (int i = 1; i < g.size(); ++i) {
- if (g[i] <= g[i - 1]) {
- set_errmsg("Values in for x-, y-, and z- grids for "
- "rectilinear meshes must be sorted and unique.");
- return OPENMC_E_INVALID_ARGUMENT;
- }
+ if (std::adjacent_find(g.begin(), g.end(), std::greater_equal<>()) !=
+ g.end()) {
+ set_errmsg("Values in for x-, y-, and z- grids for "
+ "rectilinear meshes must be sorted and unique.");
+ return OPENMC_E_INVALID_ARGUMENT;
}
}
@@ -1120,169 +865,6 @@ int RectilinearMesh::set_grid()
return 0;
}
-void RectilinearMesh::surface_bins_crossed(
- Position r0, Position r1, const Direction& u, vector& bins) const
-{
- // Determine indices for starting and ending location.
- int ijk0[3], ijk1[3];
- bool start_in_mesh;
- get_indices(r0, ijk0, &start_in_mesh);
- bool end_in_mesh;
- get_indices(r1, ijk1, &end_in_mesh);
-
- // If the starting coordinates do not lie in the mesh, compute the coords and
- // mesh indices of the first intersection, and add the bin for this first
- // intersection. Return if the particle does not intersect the mesh at all.
- if (!start_in_mesh) {
- // Compute the incoming intersection coordinates and indices.
- if (!intersects(r0, r1, ijk0))
- return;
-
- // Determine which surface the particle entered.
- double min_dist = INFTY;
- int i_surf;
- for (int i = 0; i < 3; ++i) {
- if (u[i] > 0.0 && ijk0[i] == 1) {
- double d = std::abs(r0[i] - grid_[i][0]);
- if (d < min_dist) {
- min_dist = d;
- i_surf = 4 * i + 2;
- }
- } else if (u[i] < 0.0 && ijk0[i] == shape_[i]) {
- double d = std::abs(r0[i] - grid_[i][shape_[i]]);
- if (d < min_dist) {
- min_dist = d;
- i_surf = 4 * i + 4;
- }
- } // u[i] == 0 intentionally skipped
- }
-
- // Add the incoming current bin.
- int i_mesh = get_bin_from_indices(ijk0);
- int i_bin = 4 * 3 * i_mesh + i_surf - 1;
- bins.push_back(i_bin);
- }
-
- // If the ending coordinates do not lie in the mesh, compute the coords and
- // mesh indices of the last intersection, and add the bin for this last
- // intersection.
- if (!end_in_mesh) {
- // Compute the outgoing intersection coordinates and indices.
- intersects(r1, r0, ijk1);
-
- // Determine which surface the particle exited.
- double min_dist = INFTY;
- int i_surf;
- for (int i = 0; i < 3; ++i) {
- if (u[i] > 0.0 && ijk1[i] == shape_[i]) {
- double d = std::abs(r1[i] - grid_[i][shape_[i]]);
- if (d < min_dist) {
- min_dist = d;
- i_surf = 4 * i + 3;
- }
- } else if (u[i] < 0.0 && ijk1[i] == 1) {
- double d = std::abs(r1[i] - grid_[i][0]);
- if (d < min_dist) {
- min_dist = d;
- i_surf = 4 * i + 1;
- }
- } // u[i] == 0 intentionally skipped
- }
-
- // Add the outgoing current bin.
- int i_mesh = get_bin_from_indices(ijk1);
- int i_bin = 4 * 3 * i_mesh + i_surf - 1;
- bins.push_back(i_bin);
- }
-
- // ========================================================================
- // Find which mesh surfaces are crossed.
-
- // Calculate number of surface crossings
- int n_cross = 0;
- for (int i = 0; i < 3; ++i)
- n_cross += std::abs(ijk1[i] - ijk0[i]);
- if (n_cross == 0)
- return;
-
- // Bounding coordinates
- Position xyz_cross;
- for (int i = 0; i < 3; ++i) {
- if (u[i] > 0.0) {
- xyz_cross[i] = positive_grid_boundary(ijk0, i);
- } else {
- xyz_cross[i] = negative_grid_boundary(ijk0, i);
- }
- }
-
- for (int j = 0; j < n_cross; ++j) {
- // Set the distances to infinity
- Position d {INFTY, INFTY, INFTY};
-
- // Determine closest bounding surface. We need to treat
- // special case where the cosine of the angle is zero since this would
- // result in a divide-by-zero.
- double distance = INFTY;
- for (int i = 0; i < 3; ++i) {
- if (u[i] == 0) {
- d[i] = INFTY;
- } else {
- d[i] = (xyz_cross[i] - r0[i]) / u[i];
- }
- distance = std::min(distance, d[i]);
- }
-
- // Loop over the dimensions
- for (int i = 0; i < 3; ++i) {
- // Check whether distance is the shortest distance
- if (distance == d[i]) {
-
- // Check whether particle is moving in positive i direction
- if (u[i] > 0) {
-
- // Outward current on i max surface
- int i_surf = 4 * i + 3;
- int i_mesh = get_bin_from_indices(ijk0);
- int i_bin = 4 * 3 * i_mesh + i_surf - 1;
- bins.push_back(i_bin);
-
- // Advance position
- ++ijk0[i];
- xyz_cross[i] = grid_[i][ijk0[i]];
-
- // Inward current on i min surface
- i_surf = 4 * i + 2;
- i_mesh = get_bin_from_indices(ijk0);
- i_bin = 4 * 3 * i_mesh + i_surf - 1;
- bins.push_back(i_bin);
-
- } else {
- // The particle is moving in the negative i direction
-
- // Outward current on i min surface
- int i_surf = 4 * i + 1;
- int i_mesh = get_bin_from_indices(ijk0);
- int i_bin = 4 * 3 * i_mesh + i_surf - 1;
- bins.push_back(i_bin);
-
- // Advance position
- --ijk0[i];
- xyz_cross[i] = grid_[i][ijk0[i] - 1];
-
- // Inward current on i min surface
- i_surf = 4 * i + 4;
- i_mesh = get_bin_from_indices(ijk0);
- i_bin = 4 * 3 * i_mesh + i_surf - 1;
- bins.push_back(i_bin);
- }
- }
- }
-
- // Calculate new coordinates
- r0 += distance * u;
- }
-}
-
int RectilinearMesh::get_index_in_direction(double r, int i) const
{
return lower_bound_index(grid_[i].begin(), grid_[i].end(), r) + 1;
@@ -1330,6 +912,513 @@ void RectilinearMesh::to_hdf5(hid_t group) const
close_group(mesh_group);
}
+//==============================================================================
+// CylindricalMesh implementation
+//==============================================================================
+
+CylindricalMesh::CylindricalMesh(pugi::xml_node node) : StructuredMesh {node}
+{
+ n_dimension_ = 3;
+
+ grid_[0] = get_node_array(node, "r_grid");
+ grid_[1] = get_node_array(node, "phi_grid");
+ grid_[2] = get_node_array(node, "z_grid");
+
+ if (int err = set_grid()) {
+ fatal_error(openmc_err_msg);
+ }
+}
+
+const std::string CylindricalMesh::mesh_type = "cylindrical";
+
+std::string CylindricalMesh::get_mesh_type() const
+{
+ return mesh_type;
+}
+
+StructuredMesh::MeshIndex CylindricalMesh::get_indices(
+ Position r, bool& in_mesh) const
+{
+ Position mapped_r;
+
+ mapped_r[0] = std::hypot(r.x, r.y);
+ mapped_r[2] = r[2];
+ if (mapped_r[0] < FP_PRECISION) {
+ mapped_r[1] = 0.0;
+ } else {
+ mapped_r[1] = std::atan2(r.y, r.x);
+ if (mapped_r[1] < 0)
+ mapped_r[1] += 2 * M_PI;
+ }
+
+ MeshIndex idx = StructuredMesh::get_indices(mapped_r, in_mesh);
+
+ idx[1] = sanitize_phi(idx[1]);
+
+ return idx;
+}
+
+double CylindricalMesh::find_r_crossing(
+ const Position& r, const Direction& u, double l, int shell) const
+{
+
+ if ((shell < 0) || (shell >= shape_[0]))
+ return INFTY;
+
+ // solve r.x^2 + r.y^2 == r0^2
+ // x^2 + 2*s*u*x + s^2*u^2 + s^2*v^2+2*s*v*y + y^2 -r0^2 = 0
+ // s^2 * (u^2 + v^2) + 2*s*(u*x+v*y) + x^2+y^2-r0^2 = 0
+
+ const double r0 = grid_[0][shell];
+ const double denominator = u.x * u.x + u.y * u.y;
+
+ // Direction of flight is in z-direction. Will never intersect r.
+ if (std::abs(denominator) < FP_PRECISION)
+ return INFTY;
+
+ // inverse of dominator to help the compiler to speed things up
+ const double inv_denominator = 1.0 / denominator;
+
+ const double p = (u.x * r.x + u.y * r.y) * inv_denominator;
+ double D = p * p + (r0 * r0 - r.x * r.x - r.y * r.y) * inv_denominator;
+
+ if (D < 0.0)
+ return INFTY;
+
+ D = std::sqrt(D);
+
+ // the solution -p - D is always smaller as -p + D : Check this one first
+ if (-p - D > l)
+ return -p - D;
+ if (-p + D > l)
+ return -p + D;
+
+ return INFTY;
+}
+
+double CylindricalMesh::find_phi_crossing(
+ const Position& r, const Direction& u, double l, int shell) const
+{
+ // Phi grid is [0, 360], thus there is no real surface to cross
+ if (full_phi_ && (shape_[1] == 1))
+ return INFTY;
+
+ shell = sanitize_phi(shell);
+
+ const double p0 = grid_[1][shell];
+
+ // solve y(s)/x(s) = tan(p0) = sin(p0)/cos(p0)
+ // => x(s) * cos(p0) = y(s) * sin(p0)
+ // => (y + s * v) * cos(p0) = (x + s * u) * sin(p0)
+ // = s * (v * cos(p0) - u * sin(p0)) = - (y * cos(p0) - x * sin(p0))
+
+ const double c0 = std::cos(p0);
+ const double s0 = std::sin(p0);
+
+ const double denominator = (u.x * s0 - u.y * c0);
+
+ // Check if direction of flight is not parallel to phi surface
+ if (std::abs(denominator) > FP_PRECISION) {
+ const double s = -(r.x * s0 - r.y * c0) / denominator;
+ // Check if solution is in positive direction of flight and crosses the
+ // correct phi surface (not -phi)
+ if ((s > l) && ((c0 * (r.x + s * u.x) + s0 * (r.y + s * u.y)) > 0.0))
+ return s;
+ }
+
+ return INFTY;
+}
+
+StructuredMesh::MeshDistance CylindricalMesh::find_z_crossing(
+ const Position& r, const Direction& u, double l, int shell) const
+{
+ MeshDistance d;
+ d.next_index = shell;
+
+ // Direction of flight is within xy-plane. Will never intersect z.
+ if (std::abs(u.z) < FP_PRECISION)
+ return d;
+
+ d.max_surface = (u.z > 0.0);
+ if (d.max_surface && (shell <= shape_[2])) {
+ d.next_index += 1;
+ d.distance = (grid_[2][shell] - r.z) / u.z;
+ } else if (!d.max_surface && (shell > 0)) {
+ d.next_index -= 1;
+ d.distance = (grid_[2][shell - 1] - r.z) / u.z;
+ }
+ return d;
+}
+
+StructuredMesh::MeshDistance CylindricalMesh::distance_to_grid_boundary(
+ const MeshIndex& ijk, int i, const Position& r0, const Direction& u,
+ double l) const
+{
+
+ if (i == 0) {
+
+ return std::min(
+ MeshDistance(ijk[i] + 1, true, find_r_crossing(r0, u, l, ijk[i])),
+ MeshDistance(ijk[i] - 1, false, find_r_crossing(r0, u, l, ijk[i] - 1)));
+
+ } else if (i == 1) {
+
+ return std::min(MeshDistance(sanitize_phi(ijk[i] + 1), true,
+ find_phi_crossing(r0, u, l, ijk[i])),
+ MeshDistance(sanitize_phi(ijk[i] - 1), false,
+ find_phi_crossing(r0, u, l, ijk[i] - 1)));
+
+ } else {
+ return find_z_crossing(r0, u, l, ijk[i]);
+ }
+}
+
+int CylindricalMesh::set_grid()
+{
+ shape_ = {static_cast(grid_[0].size()) - 1,
+ static_cast(grid_[1].size()) - 1,
+ static_cast(grid_[2].size()) - 1};
+
+ for (const auto& g : grid_) {
+ if (g.size() < 2) {
+ set_errmsg("r-, phi-, and z- grids for cylindrical meshes "
+ "must each have at least 2 points");
+ return OPENMC_E_INVALID_ARGUMENT;
+ }
+ if (std::adjacent_find(g.begin(), g.end(), std::greater_equal<>()) !=
+ g.end()) {
+ set_errmsg("Values in for r-, phi-, and z- grids for "
+ "cylindrical meshes must be sorted and unique.");
+ return OPENMC_E_INVALID_ARGUMENT;
+ }
+ }
+ if (grid_[0].front() < 0.0) {
+ set_errmsg("r-grid for "
+ "cylindrical meshes must start at r >= 0.");
+ return OPENMC_E_INVALID_ARGUMENT;
+ }
+ if (grid_[1].front() < 0.0) {
+ set_errmsg("phi-grid for "
+ "cylindrical meshes must start at phi >= 0.");
+ return OPENMC_E_INVALID_ARGUMENT;
+ }
+ if (grid_[1].back() > 360) {
+ set_errmsg("phi-grids for "
+ "cylindrical meshes must end with theta <= 360 degree.");
+
+ return OPENMC_E_INVALID_ARGUMENT;
+ }
+
+ full_phi_ = (grid_[1].front() == 0.0) && (grid_[1].back() == 360.0);
+
+ // Transform phi-grid from degrees to radians
+ std::transform(grid_[1].begin(), grid_[1].end(), grid_[1].begin(),
+ [](double d) { return M_PI * d / 180.0; });
+
+ lower_left_ = {grid_[0].front(), grid_[1].front(), grid_[2].front()};
+ upper_right_ = {grid_[0].back(), grid_[1].back(), grid_[2].back()};
+
+ return 0;
+}
+
+int CylindricalMesh::get_index_in_direction(double r, int i) const
+{
+ return lower_bound_index(grid_[i].begin(), grid_[i].end(), r) + 1;
+}
+
+std::pair, vector> CylindricalMesh::plot(
+ Position plot_ll, Position plot_ur) const
+{
+ fatal_error("Plot of cylindrical Mesh not implemented");
+
+ // Figure out which axes lie in the plane of the plot.
+ array, 2> axis_lines;
+ return {axis_lines[0], axis_lines[1]};
+}
+
+void CylindricalMesh::to_hdf5(hid_t group) const
+{
+ hid_t mesh_group = create_group(group, "mesh " + std::to_string(id_));
+
+ write_dataset(mesh_group, "type", "cylindrical");
+ write_dataset(mesh_group, "r_grid", grid_[0]);
+ write_dataset(mesh_group, "phi_grid", grid_[1]);
+ write_dataset(mesh_group, "z_grid", grid_[2]);
+
+ close_group(mesh_group);
+}
+
+//==============================================================================
+// SphericalMesh implementation
+//==============================================================================
+
+SphericalMesh::SphericalMesh(pugi::xml_node node) : StructuredMesh {node}
+{
+ n_dimension_ = 3;
+
+ grid_[0] = get_node_array(node, "r_grid");
+ grid_[1] = get_node_array(node, "theta_grid");
+ grid_[2] = get_node_array(node, "phi_grid");
+
+ if (int err = set_grid()) {
+ fatal_error(openmc_err_msg);
+ }
+}
+
+const std::string SphericalMesh::mesh_type = "spherical";
+
+std::string SphericalMesh::get_mesh_type() const
+{
+ return mesh_type;
+}
+
+StructuredMesh::MeshIndex SphericalMesh::get_indices(
+ Position r, bool& in_mesh) const
+{
+ Position mapped_r;
+
+ mapped_r[0] = r.norm();
+ if (mapped_r[0] < FP_PRECISION) {
+ mapped_r[1] = 0.0;
+ mapped_r[2] = 0.0;
+ } else {
+ mapped_r[1] = std::acos(r.z / mapped_r.x);
+ mapped_r[2] = std::atan2(r.y, r.x);
+ if (mapped_r[2] < 0)
+ mapped_r[2] += 2 * M_PI;
+ }
+
+ MeshIndex idx = StructuredMesh::get_indices(mapped_r, in_mesh);
+
+ idx[1] = sanitize_theta(idx[1]);
+ idx[2] = sanitize_phi(idx[2]);
+
+ return idx;
+}
+
+double SphericalMesh::find_r_crossing(
+ const Position& r, const Direction& u, double l, int shell) const
+{
+ if ((shell < 0) || (shell >= shape_[0]))
+ return INFTY;
+
+ // solve |r+s*u| = r0
+ // |r+s*u| = |r| + 2*s*r*u + s^2 (|u|==1 !)
+ const double r0 = grid_[0][shell];
+ const double p = r.dot(u);
+ double D = p * p - r.dot(r) + r0 * r0;
+
+ if (D >= 0.0) {
+ D = std::sqrt(D);
+ // the solution -p - D is always smaller as -p + D : Check this one first
+ if (-p - D > l)
+ return -p - D;
+ if (-p + D > l)
+ return -p + D;
+ }
+
+ return INFTY;
+}
+
+double SphericalMesh::find_theta_crossing(
+ const Position& r, const Direction& u, double l, int shell) const
+{
+ // Theta grid is [0, 180], thus there is no real surface to cross
+ if (full_theta_ && (shape_[1] == 1))
+ return INFTY;
+
+ shell = sanitize_theta(shell);
+
+ // solving z(s) = cos/theta) * r(s) with r(s) = r+s*u
+ // yields
+ // a*s^2 + 2*b*s + c == 0 with
+ // a = cos(theta)^2 - u.z * u.z
+ // b = r*u * cos(theta)^2 - u.z * r.z
+ // c = r*r * cos(theta)^2 - r.z^2
+
+ const double cos_t = std::cos(grid_[1][shell]);
+ const bool sgn = std::signbit(cos_t);
+ const double cos_t_2 = cos_t * cos_t;
+
+ const double a = cos_t_2 - u.z * u.z;
+ const double b = r.dot(u) * cos_t_2 - r.z * u.z;
+ const double c = r.dot(r) * cos_t_2 - r.z * r.z;
+
+ // if factor of s^2 is zero, direction of flight is parallel to theta surface
+ if (std::abs(a) < FP_PRECISION) {
+ // if b vanishes, direction of flight is within theta surface and crossing
+ // is not possible
+ if (std::abs(b) < FP_PRECISION)
+ return INFTY;
+
+ const double s = -0.5 * c / b;
+ // Check if solution is in positive direction of flight and has correct sign
+ if ((s > l) && (std::signbit(r.z + s * u.z) == sgn))
+ return s;
+
+ // no crossing is possible
+ return INFTY;
+ }
+
+ const double p = b / a;
+ double D = p * p - c / a;
+
+ if (D < 0.0)
+ return INFTY;
+
+ D = std::sqrt(D);
+
+ // the solution -p-D is always smaller as -p+D : Check this one first
+ double s = -p - D;
+ // Check if solution is in positive direction of flight and has correct sign
+ if ((s > l) && (std::signbit(r.z + s * u.z) == sgn))
+ return s;
+
+ s = -p + D;
+ // Check if solution is in positive direction of flight and has correct sign
+ if ((s > l) && (std::signbit(r.z + s * u.z) == sgn))
+ return s;
+
+ return INFTY;
+}
+
+double SphericalMesh::find_phi_crossing(
+ const Position& r, const Direction& u, double l, int shell) const
+{
+ // Phi grid is [0, 360], thus there is no real surface to cross
+ if (full_phi_ && (shape_[2] == 1))
+ return INFTY;
+
+ shell = sanitize_phi(shell);
+
+ const double p0 = grid_[2][shell];
+
+ // solve y(s)/x(s) = tan(p0) = sin(p0)/cos(p0)
+ // => x(s) * cos(p0) = y(s) * sin(p0)
+ // => (y + s * v) * cos(p0) = (x + s * u) * sin(p0)
+ // = s * (v * cos(p0) - u * sin(p0)) = - (y * cos(p0) - x * sin(p0))
+
+ const double c0 = std::cos(p0);
+ const double s0 = std::sin(p0);
+
+ const double denominator = (u.x * s0 - u.y * c0);
+
+ // Check if direction of flight is not parallel to phi surface
+ if (std::abs(denominator) > FP_PRECISION) {
+ const double s = -(r.x * s0 - r.y * c0) / denominator;
+ // Check if solution is in positive direction of flight and crosses the
+ // correct phi surface (not -phi)
+ if ((s > l) && ((c0 * (r.x + s * u.x) + s0 * (r.y + s * u.y)) > 0.0))
+ return s;
+ }
+
+ return INFTY;
+}
+
+StructuredMesh::MeshDistance SphericalMesh::distance_to_grid_boundary(
+ const MeshIndex& ijk, int i, const Position& r0, const Direction& u,
+ double l) const
+{
+
+ if (i == 0) {
+
+ return std::min(
+ MeshDistance(ijk[i] + 1, true, find_r_crossing(r0, u, l, ijk[i])),
+ MeshDistance(ijk[i] - 1, false, find_r_crossing(r0, u, l, ijk[i] - 1)));
+
+ } else if (i == 1) {
+
+ return std::min(MeshDistance(sanitize_theta(ijk[i] + 1), true,
+ find_theta_crossing(r0, u, l, ijk[i])),
+ MeshDistance(sanitize_theta(ijk[i] - 1), false,
+ find_theta_crossing(r0, u, l, ijk[i] - 1)));
+
+ } else {
+
+ return std::min(MeshDistance(sanitize_phi(ijk[i] + 1), true,
+ find_phi_crossing(r0, u, l, ijk[i])),
+ MeshDistance(sanitize_phi(ijk[i] - 1), false,
+ find_phi_crossing(r0, u, l, ijk[i] - 1)));
+ }
+}
+
+int SphericalMesh::set_grid()
+{
+ shape_ = {static_cast(grid_[0].size()) - 1,
+ static_cast(grid_[1].size()) - 1,
+ static_cast(grid_[2].size()) - 1};
+
+ for (const auto& g : grid_) {
+ if (g.size() < 2) {
+ set_errmsg("x-, y-, and z- grids for spherical meshes "
+ "must each have at least 2 points");
+ return OPENMC_E_INVALID_ARGUMENT;
+ }
+ if (std::adjacent_find(g.begin(), g.end(), std::greater_equal<>()) !=
+ g.end()) {
+ set_errmsg("Values in for r-, theta-, and phi- grids for "
+ "spherical meshes must be sorted and unique.");
+ return OPENMC_E_INVALID_ARGUMENT;
+ }
+ if (g.front() < 0.0) {
+ set_errmsg("r-, theta-, and phi- grids for "
+ "spherical meshes must start at v >= 0.");
+ return OPENMC_E_INVALID_ARGUMENT;
+ }
+ }
+ if (grid_[1].back() > 180) {
+ set_errmsg("theta-grids for "
+ "spherical meshes must end with theta <= 180 degree.");
+
+ return OPENMC_E_INVALID_ARGUMENT;
+ }
+ if (grid_[2].back() > 360) {
+ set_errmsg("phi-grids for "
+ "spherical meshes must end with phi <= 180 degree.");
+ return OPENMC_E_INVALID_ARGUMENT;
+ }
+
+ full_theta_ = (grid_[1].front() == 0.0) && (grid_[1].back() == 180.0);
+ full_phi_ = (grid_[2].front() == 0.0) && (grid_[2].back() == 360.0);
+
+ // Transform theta- and phi-grid from degrees to radians
+ for (int i = 1; i < 3; i++)
+ std::transform(grid_[i].begin(), grid_[i].end(), grid_[i].begin(),
+ [](double d) { return M_PI * d / 180.0; });
+
+ lower_left_ = {grid_[0].front(), grid_[1].front(), grid_[2].front()};
+ upper_right_ = {grid_[0].back(), grid_[1].back(), grid_[2].back()};
+
+ return 0;
+}
+
+int SphericalMesh::get_index_in_direction(double r, int i) const
+{
+ return lower_bound_index(grid_[i].begin(), grid_[i].end(), r) + 1;
+}
+
+std::pair, vector> SphericalMesh::plot(
+ Position plot_ll, Position plot_ur) const
+{
+ fatal_error("Plot of spherical Mesh not implemented");
+
+ // Figure out which axes lie in the plane of the plot.
+ array, 2> axis_lines;
+ return {axis_lines[0], axis_lines[1]};
+}
+
+void SphericalMesh::to_hdf5(hid_t group) const
+{
+ hid_t mesh_group = create_group(group, "mesh " + std::to_string(id_));
+
+ write_dataset(mesh_group, "type", SphericalMesh::mesh_type);
+ write_dataset(mesh_group, "r_grid", grid_[0]);
+ write_dataset(mesh_group, "theta_grid", grid_[1]);
+ write_dataset(mesh_group, "phi_grid", grid_[2]);
+
+ close_group(mesh_group);
+}
+
//==============================================================================
// Helper functions for the C API
//==============================================================================
@@ -1357,6 +1446,13 @@ int check_mesh_type(int32_t index)
return 0;
}
+template
+bool is_mesh_type(int32_t index)
+{
+ T* mesh = dynamic_cast(model::meshes[index].get());
+ return mesh;
+}
+
//==============================================================================
// C API functions
//==============================================================================
@@ -1367,16 +1463,8 @@ extern "C" int openmc_mesh_get_type(int32_t index, char* type)
if (int err = check_mesh(index))
return err;
- RegularMesh* mesh = dynamic_cast(model::meshes[index].get());
- if (mesh) {
- std::strcpy(type, "regular");
- } else {
- RectilinearMesh* mesh =
- dynamic_cast(model::meshes[index].get());
- if (mesh) {
- std::strcpy(type, "rectilinear");
- }
- }
+ std::strcpy(type, model::meshes[index].get()->get_mesh_type().c_str());
+
return 0;
}
@@ -1389,10 +1477,14 @@ extern "C" int openmc_extend_meshes(
std::string mesh_type;
for (int i = 0; i < n; ++i) {
- if (std::strcmp(type, "regular") == 0) {
+ if (RegularMesh::mesh_type == type) {
model::meshes.push_back(make_unique());
- } else if (std::strcmp(type, "rectilinear") == 0) {
+ } else if (RectilinearMesh::mesh_type == type) {
model::meshes.push_back(make_unique());
+ } else if (CylindricalMesh::mesh_type == type) {
+ model::meshes.push_back(make_unique());
+ } else if (SphericalMesh::mesh_type == type) {
+ model::meshes.push_back(make_unique());
} else {
throw std::runtime_error {"Unknown mesh type: " + std::string(type)};
}
@@ -1412,14 +1504,14 @@ extern "C" int openmc_add_unstructured_mesh(
bool valid_lib = false;
#ifdef DAGMC
- if (lib_name == "moab") {
+ if (lib_name == MOABMesh::mesh_lib_type) {
model::meshes.push_back(std::move(make_unique(mesh_file)));
valid_lib = true;
}
#endif
#ifdef LIBMESH
- if (lib_name == "libmesh") {
+ if (lib_name == LibMesh::mesh_lib_type) {
model::meshes.push_back(std::move(make_unique(mesh_file)));
valid_lib = true;
}
@@ -1491,9 +1583,8 @@ extern "C" int openmc_regular_mesh_set_dimension(
RegularMesh* mesh = dynamic_cast(model::meshes[index].get());
// Copy dimension
- vector shape = {static_cast(n)};
- mesh->shape_ = xt::adapt(dims, n, xt::no_ownership(), shape);
- mesh->n_dimension_ = mesh->shape_.size();
+ mesh->n_dimension_ = n;
+ std::copy(dims, dims + n, mesh->shape_.begin());
return 0;
}
@@ -1529,15 +1620,15 @@ extern "C" int openmc_regular_mesh_set_params(
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->width_ = (m->upper_right_ - m->lower_left_) / m->get_x_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->upper_right_ = m->lower_left_ + m->get_x_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->lower_left_ = m->upper_right_ - m->get_x_shape() * m->width_;
} else {
set_errmsg("At least two parameters must be specified.");
return OPENMC_E_INVALID_ARGUMENT;
@@ -1546,14 +1637,45 @@ extern "C" int openmc_regular_mesh_set_params(
return 0;
}
-//! Get the rectilinear mesh grid
-extern "C" int openmc_rectilinear_mesh_get_grid(int32_t index, double** grid_x,
+//! Set the mesh parameters for rectilinear, cylindrical and spharical meshes
+template
+int openmc_structured_mesh_set_grid_impl(int32_t index, const double* grid_x,
+ const int nx, const double* grid_y, const int ny, const double* grid_z,
+ const int nz)
+{
+ if (int err = check_mesh_type(index))
+ return err;
+
+ C* m = dynamic_cast(model::meshes[index].get());
+
+ m->n_dimension_ = 3;
+
+ m->grid_[0].reserve(nx);
+ m->grid_[1].reserve(ny);
+ m->grid_[2].reserve(nz);
+
+ for (int i = 0; i < nx; i++) {
+ m->grid_[0].push_back(grid_x[i]);
+ }
+ for (int i = 0; i < ny; i++) {
+ m->grid_[1].push_back(grid_y[i]);
+ }
+ for (int i = 0; i < nz; i++) {
+ m->grid_[2].push_back(grid_z[i]);
+ }
+
+ int err = m->set_grid();
+ return err;
+}
+
+//! Get the mesh parameters for rectilinear, cylindrical and spherical meshes
+template
+int openmc_structured_mesh_get_grid_impl(int32_t index, double** grid_x,
int* nx, double** grid_y, int* ny, double** grid_z, int* nz)
{
- if (int err = check_mesh_type(index))
+ if (int err = check_mesh_type(index))
return err;
- RectilinearMesh* m =
- dynamic_cast(model::meshes[index].get());
+ C* m = dynamic_cast(model::meshes[index].get());
if (m->lower_left_.dimension() == 0) {
set_errmsg("Mesh parameters have not been set.");
@@ -1570,35 +1692,63 @@ extern "C" int openmc_rectilinear_mesh_get_grid(int32_t index, double** grid_x,
return 0;
}
+//! Get the rectilinear mesh grid
+extern "C" int openmc_rectilinear_mesh_get_grid(int32_t index, double** grid_x,
+ int* nx, double** grid_y, int* ny, double** grid_z, int* nz)
+{
+ return openmc_structured_mesh_get_grid_impl(
+ index, grid_x, nx, grid_y, ny, grid_z, nz);
+}
+
//! Set the rectilienar mesh parameters
extern "C" int openmc_rectilinear_mesh_set_grid(int32_t index,
const double* grid_x, const int nx, const double* grid_y, const int ny,
const double* grid_z, const int nz)
{
- if (int err = check_mesh_type(index))
- return err;
- RectilinearMesh* m =
- dynamic_cast(model::meshes[index].get());
+ return openmc_structured_mesh_set_grid_impl(
+ index, grid_x, nx, grid_y, ny, grid_z, nz);
+}
- m->n_dimension_ = 3;
- m->grid_.resize(m->n_dimension_);
+//! Get the cylindrical mesh grid
+extern "C" int openmc_cylindrical_mesh_get_grid(int32_t index, double** grid_x,
+ int* nx, double** grid_y, int* ny, double** grid_z, int* nz)
+{
+ return openmc_structured_mesh_get_grid_impl(
+ index, grid_x, nx, grid_y, ny, grid_z, nz);
+}
- for (int i = 0; i < nx; i++) {
- m->grid_[0].push_back(grid_x[i]);
- }
- for (int i = 0; i < ny; i++) {
- m->grid_[1].push_back(grid_y[i]);
- }
- for (int i = 0; i < nz; i++) {
- m->grid_[2].push_back(grid_z[i]);
- }
+//! Set the cylindrical mesh parameters
+extern "C" int openmc_cylindrical_mesh_set_grid(int32_t index,
+ const double* grid_x, const int nx, const double* grid_y, const int ny,
+ const double* grid_z, const int nz)
+{
+ return openmc_structured_mesh_set_grid_impl(
+ index, grid_x, nx, grid_y, ny, grid_z, nz);
+}
- int err = m->set_grid();
- return err;
+//! Get the spherical mesh grid
+extern "C" int openmc_spherical_mesh_get_grid(int32_t index, double** grid_x,
+ int* nx, double** grid_y, int* ny, double** grid_z, int* nz)
+{
+
+ return openmc_structured_mesh_get_grid_impl(
+ index, grid_x, nx, grid_y, ny, grid_z, nz);
+ ;
+}
+
+//! Set the spherical mesh parameters
+extern "C" int openmc_spherical_mesh_set_grid(int32_t index,
+ const double* grid_x, const int nx, const double* grid_y, const int ny,
+ const double* grid_z, const int nz)
+{
+ return openmc_structured_mesh_set_grid_impl(
+ index, grid_x, nx, grid_y, ny, grid_z, nz);
}
#ifdef DAGMC
+const std::string MOABMesh::mesh_lib_type = "moab";
+
MOABMesh::MOABMesh(pugi::xml_node node) : UnstructuredMesh(node)
{
initialize();
@@ -1868,7 +2018,7 @@ double MOABMesh::volume(int bin) const
std::string MOABMesh::library() const
{
- return "moab";
+ return mesh_lib_type;
}
double MOABMesh::tet_volume(moab::EntityHandle tet) const
@@ -2187,6 +2337,8 @@ void MOABMesh::write(const std::string& base_filename) const
#ifdef LIBMESH
+const std::string LibMesh::mesh_lib_type = "libmesh";
+
LibMesh::LibMesh(pugi::xml_node node) : UnstructuredMesh(node)
{
initialize();
@@ -2261,7 +2413,7 @@ Position LibMesh::centroid(int bin) const
std::string LibMesh::library() const
{
- return "libmesh";
+ return mesh_lib_type;
}
int LibMesh::n_bins() const
@@ -2438,19 +2590,25 @@ void read_meshes(pugi::xml_node root)
}
// Read mesh and add to vector
- if (mesh_type == "regular") {
+ if (mesh_type == RegularMesh::mesh_type) {
model::meshes.push_back(make_unique(node));
- } else if (mesh_type == "rectilinear") {
+ } else if (mesh_type == RectilinearMesh::mesh_type) {
model::meshes.push_back(make_unique(node));
+ } else if (mesh_type == CylindricalMesh::mesh_type) {
+ model::meshes.push_back(make_unique(node));
+ } else if (mesh_type == SphericalMesh::mesh_type) {
+ model::meshes.push_back(make_unique(node));
#ifdef DAGMC
- } else if (mesh_type == "unstructured" && mesh_lib == "moab") {
+ } else if (mesh_type == UnstructuredMesh::mesh_type &&
+ mesh_lib == MOABMesh::mesh_lib_type) {
model::meshes.push_back(make_unique(node));
#endif
#ifdef LIBMESH
- } else if (mesh_type == "unstructured" && mesh_lib == "libmesh") {
+ } else if (mesh_type == UnstructuredMesh::mesh_type &&
+ mesh_lib == LibMesh::mesh_lib_type) {
model::meshes.push_back(make_unique(node));
#endif
- } else if (mesh_type == "unstructured") {
+ } else if (mesh_type == UnstructuredMesh::mesh_type) {
fatal_error("Unstructured mesh support is not enabled or the mesh "
"library is invalid.");
} else {
diff --git a/src/relaxng/tallies.rnc b/src/relaxng/tallies.rnc
index 3d504a66c8..4e46c45a17 100644
--- a/src/relaxng/tallies.rnc
+++ b/src/relaxng/tallies.rnc
@@ -24,6 +24,24 @@ element tallies {
attribute y_grid { list { xsd:double+ } }) &
(element z_grid { list { xsd:double+ } } |
attribute z_grid { list { xsd:double+ } })
+ ) | (
+ (element type { ( "cylindrical" ) } |
+ attribute type { ( "cylindrical" ) }) &
+ (element r_grid { list { xsd:double+ } } |
+ attribute r_grid { list { xsd:double+ } }) &
+ (element phi_grid { list { xsd:double+ } } |
+ attribute phi_grid { list { xsd:double+ } }) &
+ (element z_grid { list { xsd:double+ } } |
+ attribute z_grid { list { xsd:double+ } })
+ ) | (
+ (element type { ( "spherical" ) } |
+ attribute type { ( "spherical" ) }) &
+ (element r_grid { list { xsd:double+ } } |
+ attribute r_grid { list { xsd:double+ } }) &
+ (element theta_grid { list { xsd:double+ } } |
+ attribute theta_grid { list { xsd:double+ } }) &
+ (element phi_grid { list { xsd:double+ } } |
+ attribute phi_grid { list { xsd:double+ } })
)
)
}* &
diff --git a/tests/regression_tests/filter_mesh/inputs_true.dat b/tests/regression_tests/filter_mesh/inputs_true.dat
index 37716c02e4..e188768517 100644
--- a/tests/regression_tests/filter_mesh/inputs_true.dat
+++ b/tests/regression_tests/filter_mesh/inputs_true.dat
@@ -53,36 +53,58 @@
-7.5 -6.617647058823529 -5.735294117647059 -4.852941176470589 -3.9705882352941178 -3.0882352941176467 -2.2058823529411766 -1.3235294117647065 -0.4411764705882355 0.4411764705882355 1.3235294117647065 2.2058823529411757 3.0882352941176467 3.9705882352941178 4.852941176470587 5.735294117647058 6.617647058823529 7.5
1.0 1.223224374241637 1.4962778697388448 1.8302835609029084 2.2388474634702153 2.7386127875258306 3.3499379133114306 4.09772570775871 5.012437964687018 6.131336292779302 7.500000000000001
+
+ 0.0 0.4411764705882353 0.8823529411764706 1.3235294117647058 1.7647058823529411 2.2058823529411766 2.6470588235294117 3.0882352941176467 3.5294117647058822 3.9705882352941178 4.411764705882353 4.852941176470588 5.294117647058823 5.735294117647059 6.1764705882352935 6.617647058823529 7.0588235294117645 7.5
+ 0.0 20.0 40.0 60.0 80.0 100.0 120.0 140.0 160.0 180.0 200.0 220.0 240.0 260.0 280.0 300.0 320.0 340.0 360.0
+ -7.5 -6.5625 -5.625 -4.6875 -3.75 -2.8125 -1.875 -0.9375 0.0 0.9375 1.875 2.8125 3.75 4.6875 5.625 6.5625 7.5
+
+
+ 0.0 0.4411764705882353 0.8823529411764706 1.3235294117647058 1.7647058823529411 2.2058823529411766 2.6470588235294117 3.0882352941176467 3.5294117647058822 3.9705882352941178 4.411764705882353 4.852941176470588 5.294117647058823 5.735294117647059 6.1764705882352935 6.617647058823529 7.0588235294117645 7.5
+ 0.0 22.5 45.0 67.5 90.0 112.5 135.0 157.5 180.0
+ 0.0 20.0 40.0 60.0 80.0 100.0 120.0 140.0 160.0 180.0 200.0 220.0 240.0 260.0 280.0 300.0 320.0 340.0 360.0
+
1
-
+
1
2
-
+
2
3
-
+
3
4
-
+
4
+
+ 5
+
+
+ 5
+
+
+ 6
+
+
+ 6
+
1
total
- 5
+ 7
current
@@ -90,7 +112,7 @@
total
- 6
+ 8
current
@@ -98,7 +120,7 @@
total
- 7
+ 9
current
@@ -106,7 +128,23 @@
total
- 8
+ 10
+ current
+
+
+ 5
+ total
+
+
+ 11
+ current
+
+
+ 6
+ total
+
+
+ 12
current
diff --git a/tests/regression_tests/filter_mesh/results_true.dat b/tests/regression_tests/filter_mesh/results_true.dat
index 332c8b36cf..97735130e1 100644
--- a/tests/regression_tests/filter_mesh/results_true.dat
+++ b/tests/regression_tests/filter_mesh/results_true.dat
@@ -1 +1 @@
-97cf099a24099d2af794752e145841a9ef9f24980c53ee71e3fdd400d4249a7e6f3f52f63132fe2fe8d1973604a72d87910ab13bc3721f7d0bec533b87756847
\ No newline at end of file
+6e02d01fc115c54e4c60477a9f963149bc83c82e6a7bdb2d29125b131546150c59de5ac3aae5b9469b89b343ecf2574f6c949b918ccb43231e8666ee2b7c1b7c
\ No newline at end of file
diff --git a/tests/regression_tests/filter_mesh/test.py b/tests/regression_tests/filter_mesh/test.py
index c8aa871a83..eddf0a9913 100644
--- a/tests/regression_tests/filter_mesh/test.py
+++ b/tests/regression_tests/filter_mesh/test.py
@@ -51,18 +51,32 @@ def model():
recti_mesh.y_grid = np.linspace(-7.5, 7.5, 18)
recti_mesh.z_grid = np.logspace(0, np.log10(7.5), 11)
+ cyl_mesh = openmc.CylindricalMesh()
+ cyl_mesh.r_grid = np.linspace(0, 7.5, 18)
+ cyl_mesh.phi_grid = np.linspace(0, 360, 19)
+ cyl_mesh.z_grid = np.linspace(-7.5, 7.5, 17)
+
+ sph_mesh = openmc.SphericalMesh()
+ sph_mesh.r_grid = np.linspace(0, 7.5, 18)
+ sph_mesh.theta_grid = np.linspace(0, 180, 9)
+ sph_mesh.phi_grid = np.linspace(0, 360, 19)
+
# Create filters
reg_filters = [
openmc.MeshFilter(mesh_1d),
openmc.MeshFilter(mesh_2d),
openmc.MeshFilter(mesh_3d),
- openmc.MeshFilter(recti_mesh)
+ openmc.MeshFilter(recti_mesh),
+ openmc.MeshFilter(cyl_mesh),
+ openmc.MeshFilter(sph_mesh)
]
surf_filters = [
openmc.MeshSurfaceFilter(mesh_1d),
openmc.MeshSurfaceFilter(mesh_2d),
openmc.MeshSurfaceFilter(mesh_3d),
- openmc.MeshSurfaceFilter(recti_mesh)
+ openmc.MeshSurfaceFilter(recti_mesh),
+ openmc.MeshSurfaceFilter(cyl_mesh),
+ openmc.MeshSurfaceFilter(sph_mesh)
]
# Create tallies
diff --git a/tests/regression_tests/mgxs_library_mesh/results_true.dat b/tests/regression_tests/mgxs_library_mesh/results_true.dat
index a8775984f7..7a6b4abe7c 100644
--- a/tests/regression_tests/mgxs_library_mesh/results_true.dat
+++ b/tests/regression_tests/mgxs_library_mesh/results_true.dat
@@ -181,11 +181,11 @@
mesh 1 group in nuclide mean std. dev.
x y surf
3 1 1 x-max in 1 total 0.1816 0.006129
-2 1 1 x-max out 1 total 0.2718 0.095077
+2 1 1 x-max out 1 total 0.1730 0.004382
1 1 1 x-min in 1 total 0.0000 0.000000
0 1 1 x-min out 1 total 0.0000 0.000000
7 1 1 y-max in 1 total 0.1790 0.011921
-6 1 1 y-max out 1 total 0.2290 0.034821
+6 1 1 y-max out 1 total 0.1802 0.016599
5 1 1 y-min in 1 total 0.0000 0.000000
4 1 1 y-min out 1 total 0.0000 0.000000
19 1 2 x-max in 1 total 0.1872 0.012447
@@ -194,23 +194,23 @@
16 1 2 x-min out 1 total 0.0000 0.000000
23 1 2 y-max in 1 total 0.0000 0.000000
22 1 2 y-max out 1 total 0.0000 0.000000
-21 1 2 y-min in 1 total 0.2290 0.034821
+21 1 2 y-min in 1 total 0.1802 0.016599
20 1 2 y-min out 1 total 0.1790 0.011921
11 2 1 x-max in 1 total 0.0000 0.000000
10 2 1 x-max out 1 total 0.0000 0.000000
-9 2 1 x-min in 1 total 0.2718 0.095077
+9 2 1 x-min in 1 total 0.1730 0.004382
8 2 1 x-min out 1 total 0.1816 0.006129
15 2 1 y-max in 1 total 0.1778 0.009484
-14 2 1 y-max out 1 total 0.2326 0.042782
+14 2 1 y-max out 1 total 0.1822 0.010077
13 2 1 y-min in 1 total 0.0000 0.000000
12 2 1 y-min out 1 total 0.0000 0.000000
27 2 2 x-max in 1 total 0.0000 0.000000
-26 2 2 x-max out 1 total 0.0260 0.026000
+26 2 2 x-max out 1 total 0.0000 0.000000
25 2 2 x-min in 1 total 0.1952 0.015948
24 2 2 x-min out 1 total 0.1872 0.012447
31 2 2 y-max in 1 total 0.0000 0.000000
-30 2 2 y-max out 1 total 0.0244 0.024400
-29 2 2 y-min in 1 total 0.2326 0.042782
+30 2 2 y-max out 1 total 0.0000 0.000000
+29 2 2 y-min in 1 total 0.1822 0.010077
28 2 2 y-min out 1 total 0.1778 0.009484
mesh 1 group in nuclide mean std. dev.
x y z
diff --git a/tests/unit_tests/test_lib.py b/tests/unit_tests/test_lib.py
index ccff33ac39..3572e9863a 100644
--- a/tests/unit_tests/test_lib.py
+++ b/tests/unit_tests/test_lib.py
@@ -580,6 +580,68 @@ def test_rectilinear_mesh(lib_init):
msf = openmc.lib.MeshSurfaceFilter(mesh)
assert msf.mesh == mesh
+def test_cylindrical_mesh(lib_init):
+ deg2rad = lambda deg: deg*np.pi/180
+ mesh = openmc.lib.CylindricalMesh()
+ x_grid = [0., 5., 10.]
+ y_grid = [0., 10., 20.]
+ z_grid = [10., 20., 30.]
+ mesh.set_grid(x_grid, y_grid, z_grid)
+ assert np.all(mesh.lower_left == (0., 0., 10.))
+ assert np.all(mesh.upper_right == (10., deg2rad(20.), 30.))
+ assert np.all(mesh.dimension == (2, 2, 2))
+ for i, diff_x in enumerate(np.diff(x_grid)):
+ for j, diff_y in enumerate(np.diff(y_grid)):
+ for k, diff_z in enumerate(np.diff(z_grid)):
+ assert np.all(mesh.width[i, j, k, :] == (5, deg2rad(10), 10))
+
+ with pytest.raises(exc.AllocationError):
+ mesh2 = openmc.lib.CylindricalMesh(mesh.id)
+
+ meshes = openmc.lib.meshes
+ assert isinstance(meshes, Mapping)
+ assert len(meshes) == 3
+
+ mesh = meshes[mesh.id]
+ assert isinstance(mesh, openmc.lib.CylindricalMesh)
+
+ mf = openmc.lib.MeshFilter(mesh)
+ assert mf.mesh == mesh
+
+ msf = openmc.lib.MeshSurfaceFilter(mesh)
+ assert msf.mesh == mesh
+
+def test_spherical_mesh(lib_init):
+ deg2rad = lambda deg: deg*np.pi/180
+ mesh = openmc.lib.SphericalMesh()
+ x_grid = [0., 5., 10.]
+ y_grid = [0., 10., 20.]
+ z_grid = [10., 20., 30.]
+ mesh.set_grid(x_grid, y_grid, z_grid)
+ assert np.all(mesh.lower_left == (0., 0., deg2rad(10.)))
+ assert np.all(mesh.upper_right == (10., deg2rad(20.), deg2rad(30.)))
+ assert np.all(mesh.dimension == (2, 2, 2))
+ for i, diff_x in enumerate(np.diff(x_grid)):
+ for j, diff_y in enumerate(np.diff(y_grid)):
+ for k, diff_z in enumerate(np.diff(z_grid)):
+ assert np.all(abs(mesh.width[i, j, k, :] - (5, deg2rad(10), deg2rad(10))) < 1e-16)
+
+ with pytest.raises(exc.AllocationError):
+ mesh2 = openmc.lib.SphericalMesh(mesh.id)
+
+ meshes = openmc.lib.meshes
+ assert isinstance(meshes, Mapping)
+ assert len(meshes) == 4
+
+ mesh = meshes[mesh.id]
+ assert isinstance(mesh, openmc.lib.SphericalMesh)
+
+ mf = openmc.lib.MeshFilter(mesh)
+ assert mf.mesh == mesh
+
+ msf = openmc.lib.MeshSurfaceFilter(mesh)
+ assert msf.mesh == mesh
+
def test_restart(lib_init, mpi_intracomm):
# Finalize and re-init to make internal state consistent with XML.