mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 22:26:08 -04:00
Add Cylindrical & Spherical Mesh for Tally Filter
New method for ray-tracing the mesh, used for both MeshFilter & MeshSurfaceFilter. For all Structured Meshes! meshes negative/positive_grid_boundary to non-virtual function for speed MeshIndex typedef now a std::array instead of int[3] -> Get rid of pointers (see e.g. core quidelines I.13) MeshDistance class to stroe all relevant information
This commit is contained in:
Olaf Schumann
parent
e24069e439
commit
d723dd7f07
4 changed files with 1469 additions and 687 deletions
|
|
@ -134,6 +134,18 @@ public:
|
|||
StructuredMesh(pugi::xml_node node) : Mesh {node} {};
|
||||
virtual ~StructuredMesh() = default;
|
||||
|
||||
using MeshIndex = std::array<int, 3>;
|
||||
struct MeshDistance {
|
||||
MeshDistance() {};
|
||||
MeshDistance(int _index, bool _maxSurface, double _distance): nextIndex{_index}, maxSurface(_maxSurface), distance{_distance} {};
|
||||
int nextIndex { -1 };
|
||||
double distance { INFTY };
|
||||
bool maxSurface { 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 +155,18 @@ public:
|
|||
void bins_crossed(Position r0, Position r1, const Direction& u,
|
||||
vector<int>& bins, vector<double>& lengths) const override;
|
||||
|
||||
void surface_bins_crossed(Position r0, Position r1, const Direction& u,
|
||||
vector<int>& bins) const override;
|
||||
|
||||
//! Determine the 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 <class T> 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 +179,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 ijk 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;
|
||||
virtual MeshIndex get_indices_from_bin(int bin) const;
|
||||
|
||||
//! Get mesh index in a particular direction
|
||||
//!
|
||||
|
|
@ -176,25 +200,11 @@ 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 to coordinate r to the grid for the mesh grid boundary in the negative direction
|
||||
//!
|
||||
//! \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
|
||||
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;
|
||||
|
|
@ -205,9 +215,7 @@ public:
|
|||
xt::xtensor<int, 1> 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 +229,11 @@ public:
|
|||
RegularMesh(pugi::xml_node node);
|
||||
|
||||
// Overridden methods
|
||||
void surface_bins_crossed(Position r0, Position r1, const Direction& u,
|
||||
vector<int>& bins) const override;
|
||||
|
||||
int get_index_in_direction(double r, int i) const override;
|
||||
|
||||
double positive_grid_boundary(int* ijk, int i) const override;
|
||||
|
||||
double negative_grid_boundary(int* ijk, int i) const override;
|
||||
|
||||
MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i, const Position& r0, const Direction& u, double l) const override;
|
||||
|
||||
std::pair<vector<double>, vector<double>> plot(
|
||||
Position plot_ll, Position plot_ur) const override;
|
||||
|
|
@ -236,6 +241,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,14 +273,46 @@ public:
|
|||
RectilinearMesh(pugi::xml_node node);
|
||||
|
||||
// Overridden methods
|
||||
void surface_bins_crossed(Position r0, Position r1, const Direction& u,
|
||||
vector<int>& bins) const override;
|
||||
int get_index_in_direction(double r, int i) const override;
|
||||
|
||||
MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i, const Position& r0, const Direction& u, double l) const override;
|
||||
|
||||
std::pair<vector<double>, vector<double>> plot(
|
||||
Position plot_ll, Position plot_ur) const override;
|
||||
|
||||
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;
|
||||
|
||||
|
||||
vector<vector<double>> 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;
|
||||
|
||||
double positive_grid_boundary(int* ijk, int i) const override;
|
||||
|
||||
double negative_grid_boundary(int* ijk, int i) const override;
|
||||
MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i, const Position& r0, const Direction& u, double l) const override;
|
||||
|
||||
std::pair<vector<double>, vector<double>> plot(
|
||||
Position plot_ll, Position plot_ur) const override;
|
||||
|
|
@ -274,6 +322,78 @@ public:
|
|||
vector<vector<double>> grid_;
|
||||
|
||||
int set_grid();
|
||||
|
||||
protected:
|
||||
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;
|
||||
|
||||
MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i, const Position& r0, const Direction& u, double l) const override;
|
||||
|
||||
std::pair<vector<double>, vector<double>> plot(
|
||||
Position plot_ll, Position plot_ur) const override;
|
||||
|
||||
void to_hdf5(hid_t group) const override;
|
||||
|
||||
vector<vector<double>> grid_;
|
||||
|
||||
int set_grid();
|
||||
|
||||
protected:
|
||||
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
|
||||
|
|
|
|||
|
|
@ -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': RectilinearMesh,
|
||||
'spherical': RectilinearMesh
|
||||
}
|
||||
|
||||
|
||||
|
|
|
|||
337
openmc/mesh.py
337
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,339 @@ 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.
|
||||
phi_grid : Iterable of float
|
||||
Mesh boundary points along the phi-axis.
|
||||
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 = len(self.r_grid) - 1
|
||||
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])
|
||||
t_grid_str = str(self._phi_grid) if self._phi_grid is None else len(self._phi_grid)
|
||||
string += fmt.format('\tN Phi pnts:', '=\t', t_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])
|
||||
p_grid_str = str(self._z_grid) if self._z_grid is None else len(self._z_grid)
|
||||
string += fmt.format('\tN Z pnts:', '=\t', p_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['p_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, "p_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 = int_r int_phi int_phi r dr dphi dz
|
||||
|
||||
V_r = np.array(self.r_grid)**2 / 3
|
||||
V_r = V_r[1:] - V_r[:-1]
|
||||
V_p = -np.cos(np.pi * np.array(self.phi_grid) / 180.0)
|
||||
V_p = V_p[1:] - V_p[:-1]
|
||||
V_z = np.array(self.z_grid) * np.pi / 180
|
||||
V_z = V_z[1:] - V_z[:-1]
|
||||
|
||||
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 RectilinearMesh).
|
||||
r_grid : Iterable of float
|
||||
Mesh boundary points along the r-axis.
|
||||
theta_grid : Iterable of float
|
||||
Mesh boundary points along the theta-axis.
|
||||
phi_grid : Iterable of float
|
||||
Mesh boundary points along the phi-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._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 = len(self.r_grid) - 1
|
||||
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])
|
||||
t_grid_str = str(self._theta_grid) if not self._theta_grid else len(self._theta_grid)
|
||||
string += fmt.format('\tN Theta pnts:', '=\t', t_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])
|
||||
p_grid_str = str(self._phi_grid) if not self._phi_grid else len(self._phi_grid)
|
||||
string += fmt.format('\tN Phi pnts:', '=\t', p_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['t_grid'][()]
|
||||
mesh.phi_grid = 180 / np.pi * group['p_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, "t_grid")
|
||||
subelement.text = ' '.join(map(str, self.theta_grid))
|
||||
|
||||
subelement = ET.SubElement(element, "p_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 = int_r int_theta int_phi r^2 dr sin(theta) dtheta dphi
|
||||
|
||||
V_r = np.array(self.r_grid)**3 / 3
|
||||
V_r = V_r[1:] - V_r[:-1]
|
||||
V_t = -np.cos(np.pi * np.array(self.theta_grid) / 180.0)
|
||||
V_t = V_t[1:] - V_t[:-1]
|
||||
V_p = np.array(self.phi_grid) * np.pi / 180
|
||||
V_p = V_p[1:] - V_p[:-1]
|
||||
|
||||
return np.multiply.outer(np.outer(V_r, V_t), V_p)
|
||||
|
||||
|
||||
|
||||
|
||||
class UnstructuredMesh(MeshBase):
|
||||
"""A 3D unstructured mesh
|
||||
|
|
|
|||
1432
src/mesh.cpp
1432
src/mesh.cpp
File diff suppressed because it is too large
Load diff
Loading…
Add table
Add a link
Reference in a new issue