mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 14:15:42 -04:00
load mesh objects from weight_windows.h5 file (#3598)
Co-authored-by: Patrick Shriwise <pshriwise@gmail.com>
This commit is contained in:
GuySten
GitHub
parent
3ac5d6f8f5
commit
c31032cf25
5 changed files with 579 additions and 74 deletions
|
|
@ -132,8 +132,14 @@ public:
|
|||
// Constructors and destructor
|
||||
Mesh() = default;
|
||||
Mesh(pugi::xml_node node);
|
||||
Mesh(hid_t group);
|
||||
virtual ~Mesh() = default;
|
||||
|
||||
// Factory method for creating meshes from either an XML node or HDF5 group
|
||||
template<typename T>
|
||||
static const std::unique_ptr<Mesh>& create(
|
||||
T dataset, const std::string& mesh_type, const std::string& mesh_library);
|
||||
|
||||
// Methods
|
||||
//! Perform any preparation needed to support point location within the mesh
|
||||
virtual void prepare_for_point_location() {};
|
||||
|
|
@ -258,6 +264,7 @@ class StructuredMesh : public Mesh {
|
|||
public:
|
||||
StructuredMesh() = default;
|
||||
StructuredMesh(pugi::xml_node node) : Mesh {node} {};
|
||||
StructuredMesh(hid_t group) : Mesh {group} {};
|
||||
virtual ~StructuredMesh() = default;
|
||||
|
||||
using MeshIndex = std::array<int, 3>;
|
||||
|
|
@ -423,6 +430,7 @@ class PeriodicStructuredMesh : public StructuredMesh {
|
|||
public:
|
||||
PeriodicStructuredMesh() = default;
|
||||
PeriodicStructuredMesh(pugi::xml_node node) : StructuredMesh {node} {};
|
||||
PeriodicStructuredMesh(hid_t group) : StructuredMesh {group} {};
|
||||
|
||||
Position local_coords(const Position& r) const override
|
||||
{
|
||||
|
|
@ -442,6 +450,7 @@ public:
|
|||
// Constructors
|
||||
RegularMesh() = default;
|
||||
RegularMesh(pugi::xml_node node);
|
||||
RegularMesh(hid_t group);
|
||||
|
||||
// Overridden methods
|
||||
int get_index_in_direction(double r, int i) const override;
|
||||
|
|
@ -481,6 +490,8 @@ public:
|
|||
//! Return the volume for a given mesh index
|
||||
double volume(const MeshIndex& ijk) const override;
|
||||
|
||||
int set_grid();
|
||||
|
||||
// Data members
|
||||
double volume_frac_; //!< Volume fraction of each mesh element
|
||||
double element_volume_; //!< Volume of each mesh element
|
||||
|
|
@ -492,6 +503,7 @@ public:
|
|||
// Constructors
|
||||
RectilinearMesh() = default;
|
||||
RectilinearMesh(pugi::xml_node node);
|
||||
RectilinearMesh(hid_t group);
|
||||
|
||||
// Overridden methods
|
||||
int get_index_in_direction(double r, int i) const override;
|
||||
|
|
@ -534,6 +546,7 @@ public:
|
|||
// Constructors
|
||||
CylindricalMesh() = default;
|
||||
CylindricalMesh(pugi::xml_node node);
|
||||
CylindricalMesh(hid_t group);
|
||||
|
||||
// Overridden methods
|
||||
virtual MeshIndex get_indices(Position r, bool& in_mesh) const override;
|
||||
|
|
@ -598,6 +611,7 @@ public:
|
|||
// Constructors
|
||||
SphericalMesh() = default;
|
||||
SphericalMesh(pugi::xml_node node);
|
||||
SphericalMesh(hid_t group);
|
||||
|
||||
// Overridden methods
|
||||
virtual MeshIndex get_indices(Position r, bool& in_mesh) const override;
|
||||
|
|
@ -668,6 +682,7 @@ public:
|
|||
// Constructors
|
||||
UnstructuredMesh() { n_dimension_ = 3; };
|
||||
UnstructuredMesh(pugi::xml_node node);
|
||||
UnstructuredMesh(hid_t group);
|
||||
|
||||
static const std::string mesh_type;
|
||||
virtual std::string get_mesh_type() const override;
|
||||
|
|
@ -774,6 +789,7 @@ public:
|
|||
// Constructors
|
||||
MOABMesh() = default;
|
||||
MOABMesh(pugi::xml_node);
|
||||
MOABMesh(hid_t group);
|
||||
MOABMesh(const std::string& filename, double length_multiplier = 1.0);
|
||||
MOABMesh(std::shared_ptr<moab::Interface> external_mbi);
|
||||
|
||||
|
|
@ -943,6 +959,7 @@ class LibMesh : public UnstructuredMesh {
|
|||
public:
|
||||
// Constructors
|
||||
LibMesh(pugi::xml_node node);
|
||||
LibMesh(hid_t group);
|
||||
LibMesh(const std::string& filename, double length_multiplier = 1.0);
|
||||
LibMesh(libMesh::MeshBase& input_mesh, double length_multiplier = 1.0);
|
||||
|
||||
|
|
@ -1069,6 +1086,11 @@ private:
|
|||
//! \param[in] root XML node
|
||||
void read_meshes(pugi::xml_node root);
|
||||
|
||||
//! Read meshes from an HDF5 file
|
||||
//
|
||||
//! \param[in] group HDF5 group ("meshes" group)
|
||||
void read_meshes(hid_t group);
|
||||
|
||||
//! Write mesh data to an HDF5 group
|
||||
//
|
||||
//! \param[in] group HDF5 group
|
||||
|
|
|
|||
369
src/mesh.cpp
369
src/mesh.cpp
|
|
@ -230,6 +230,42 @@ void MaterialVolumes::add_volume_unsafe(
|
|||
// Mesh implementation
|
||||
//==============================================================================
|
||||
|
||||
template<typename T>
|
||||
const std::unique_ptr<Mesh>& Mesh::create(
|
||||
T dataset, const std::string& mesh_type, const std::string& mesh_library)
|
||||
{
|
||||
// Determine mesh type. Add to model vector and map
|
||||
if (mesh_type == RegularMesh::mesh_type) {
|
||||
model::meshes.push_back(make_unique<RegularMesh>(dataset));
|
||||
} else if (mesh_type == RectilinearMesh::mesh_type) {
|
||||
model::meshes.push_back(make_unique<RectilinearMesh>(dataset));
|
||||
} else if (mesh_type == CylindricalMesh::mesh_type) {
|
||||
model::meshes.push_back(make_unique<CylindricalMesh>(dataset));
|
||||
} else if (mesh_type == SphericalMesh::mesh_type) {
|
||||
model::meshes.push_back(make_unique<SphericalMesh>(dataset));
|
||||
#ifdef OPENMC_DAGMC_ENABLED
|
||||
} else if (mesh_type == UnstructuredMesh::mesh_type &&
|
||||
mesh_library == MOABMesh::mesh_lib_type) {
|
||||
model::meshes.push_back(make_unique<MOABMesh>(dataset));
|
||||
#endif
|
||||
#ifdef OPENMC_LIBMESH_ENABLED
|
||||
} else if (mesh_type == UnstructuredMesh::mesh_type &&
|
||||
mesh_library == LibMesh::mesh_lib_type) {
|
||||
model::meshes.push_back(make_unique<LibMesh>(dataset));
|
||||
#endif
|
||||
} else if (mesh_type == UnstructuredMesh::mesh_type) {
|
||||
fatal_error("Unstructured mesh support is not enabled or the mesh "
|
||||
"library is invalid.");
|
||||
} else {
|
||||
fatal_error(fmt::format("Invalid mesh type: {}", mesh_type));
|
||||
}
|
||||
|
||||
// Map ID to position in vector
|
||||
model::mesh_map[model::meshes.back()->id_] = model::meshes.size() - 1;
|
||||
|
||||
return model::meshes.back();
|
||||
}
|
||||
|
||||
Mesh::Mesh(pugi::xml_node node)
|
||||
{
|
||||
// Read mesh id
|
||||
|
|
@ -238,6 +274,17 @@ Mesh::Mesh(pugi::xml_node node)
|
|||
name_ = get_node_value(node, "name");
|
||||
}
|
||||
|
||||
Mesh::Mesh(hid_t group)
|
||||
{
|
||||
// Read mesh ID
|
||||
read_attribute(group, "id", id_);
|
||||
|
||||
// Read mesh name
|
||||
if (object_exists(group, "name")) {
|
||||
read_dataset(group, "name", name_);
|
||||
}
|
||||
}
|
||||
|
||||
void Mesh::set_id(int32_t id)
|
||||
{
|
||||
assert(id >= 0 || id == C_NONE);
|
||||
|
|
@ -265,7 +312,13 @@ void Mesh::set_id(int32_t id)
|
|||
|
||||
// Update ID and entry in the mesh map
|
||||
id_ = id;
|
||||
model::mesh_map[id] = model::meshes.size() - 1;
|
||||
|
||||
// find the index of this mesh in the model::meshes vector
|
||||
// (search in reverse because this mesh was likely just added to the vector)
|
||||
auto it = std::find_if(model::meshes.rbegin(), model::meshes.rend(),
|
||||
[this](const std::unique_ptr<Mesh>& mesh) { return mesh.get() == this; });
|
||||
|
||||
model::mesh_map[id] = std::distance(model::meshes.begin(), it.base()) - 1;
|
||||
}
|
||||
|
||||
vector<double> Mesh::volumes() const
|
||||
|
|
@ -627,6 +680,46 @@ UnstructuredMesh::UnstructuredMesh(pugi::xml_node node) : Mesh(node)
|
|||
}
|
||||
}
|
||||
|
||||
UnstructuredMesh::UnstructuredMesh(hid_t group) : Mesh(group)
|
||||
{
|
||||
n_dimension_ = 3;
|
||||
|
||||
// check the mesh type
|
||||
if (object_exists(group, "type")) {
|
||||
std::string temp;
|
||||
read_dataset(group, "type", temp);
|
||||
if (temp != mesh_type) {
|
||||
fatal_error(fmt::format("Invalid mesh type: {}", temp));
|
||||
}
|
||||
}
|
||||
|
||||
// check if a length unit multiplier was specified
|
||||
if (object_exists(group, "length_multiplier")) {
|
||||
read_dataset(group, "length_multiplier", length_multiplier_);
|
||||
}
|
||||
|
||||
// get the filename of the unstructured mesh to load
|
||||
if (object_exists(group, "filename")) {
|
||||
read_dataset(group, "filename", filename_);
|
||||
if (!file_exists(filename_)) {
|
||||
fatal_error("Mesh file '" + filename_ + "' does not exist!");
|
||||
}
|
||||
} else {
|
||||
fatal_error(fmt::format(
|
||||
"No filename supplied for unstructured mesh with ID: {}", id_));
|
||||
}
|
||||
|
||||
if (attribute_exists(group, "options")) {
|
||||
read_attribute(group, "options", options_);
|
||||
}
|
||||
|
||||
// check if mesh tally data should be written with
|
||||
// statepoint files
|
||||
if (attribute_exists(group, "output")) {
|
||||
read_attribute(group, "output", output_);
|
||||
}
|
||||
}
|
||||
|
||||
void UnstructuredMesh::determine_bounds()
|
||||
{
|
||||
double xmin = INFTY;
|
||||
|
|
@ -1086,6 +1179,72 @@ void StructuredMesh::surface_bins_crossed(
|
|||
// RegularMesh implementation
|
||||
//==============================================================================
|
||||
|
||||
int RegularMesh::set_grid()
|
||||
{
|
||||
auto shape = xt::adapt(shape_, {n_dimension_});
|
||||
|
||||
// Check that dimensions are all greater than zero
|
||||
if (xt::any(shape <= 0)) {
|
||||
set_errmsg("All entries for a regular mesh dimensions "
|
||||
"must be positive.");
|
||||
return OPENMC_E_INVALID_ARGUMENT;
|
||||
}
|
||||
|
||||
// Make sure lower_left and dimension match
|
||||
if (lower_left_.size() != n_dimension_) {
|
||||
set_errmsg("Number of entries in lower_left must be the same "
|
||||
"as the regular mesh dimensions.");
|
||||
return OPENMC_E_INVALID_ARGUMENT;
|
||||
}
|
||||
if (width_.size() > 0) {
|
||||
|
||||
// Check to ensure width has same dimensions
|
||||
if (width_.size() != n_dimension_) {
|
||||
set_errmsg("Number of entries on width must be the same as "
|
||||
"the regular mesh dimensions.");
|
||||
return OPENMC_E_INVALID_ARGUMENT;
|
||||
}
|
||||
|
||||
// Check for negative widths
|
||||
if (xt::any(width_ < 0.0)) {
|
||||
set_errmsg("Cannot have a negative width on a regular mesh.");
|
||||
return OPENMC_E_INVALID_ARGUMENT;
|
||||
}
|
||||
|
||||
// Set width and upper right coordinate
|
||||
upper_right_ = xt::eval(lower_left_ + shape * width_);
|
||||
|
||||
} else if (upper_right_.size() > 0) {
|
||||
|
||||
// Check to ensure upper_right_ has same dimensions
|
||||
if (upper_right_.size() != n_dimension_) {
|
||||
set_errmsg("Number of entries on upper_right must be the "
|
||||
"same as the regular mesh dimensions.");
|
||||
return OPENMC_E_INVALID_ARGUMENT;
|
||||
}
|
||||
|
||||
// Check that upper-right is above lower-left
|
||||
if (xt::any(upper_right_ < lower_left_)) {
|
||||
set_errmsg(
|
||||
"The upper_right coordinates of a regular mesh must be greater than "
|
||||
"the lower_left coordinates.");
|
||||
return OPENMC_E_INVALID_ARGUMENT;
|
||||
}
|
||||
|
||||
// Set width
|
||||
width_ = xt::eval((upper_right_ - lower_left_) / shape);
|
||||
}
|
||||
|
||||
// Set material volumes
|
||||
volume_frac_ = 1.0 / xt::prod(shape)();
|
||||
|
||||
element_volume_ = 1.0;
|
||||
for (int i = 0; i < n_dimension_; i++) {
|
||||
element_volume_ *= width_[i];
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
|
||||
RegularMesh::RegularMesh(pugi::xml_node node) : StructuredMesh {node}
|
||||
{
|
||||
// Determine number of dimensions for mesh
|
||||
|
|
@ -1100,12 +1259,6 @@ RegularMesh::RegularMesh(pugi::xml_node node) : StructuredMesh {node}
|
|||
}
|
||||
std::copy(shape.begin(), shape.end(), shape_.begin());
|
||||
|
||||
// Check that dimensions are all greater than zero
|
||||
if (xt::any(shape <= 0)) {
|
||||
fatal_error("All entries on the <dimension> element for a tally "
|
||||
"mesh must be positive.");
|
||||
}
|
||||
|
||||
// Check for lower-left coordinates
|
||||
if (check_for_node(node, "lower_left")) {
|
||||
// Read mesh lower-left corner location
|
||||
|
|
@ -1114,12 +1267,6 @@ RegularMesh::RegularMesh(pugi::xml_node node) : StructuredMesh {node}
|
|||
fatal_error("Must specify <lower_left> on a mesh.");
|
||||
}
|
||||
|
||||
// Make sure lower_left and dimension match
|
||||
if (shape.size() != lower_left_.size()) {
|
||||
fatal_error("Number of entries on <lower_left> must be the same "
|
||||
"as the number of entries on <dimension>.");
|
||||
}
|
||||
|
||||
if (check_for_node(node, "width")) {
|
||||
// Make sure one of upper-right or width were specified
|
||||
if (check_for_node(node, "upper_right")) {
|
||||
|
|
@ -1128,49 +1275,52 @@ RegularMesh::RegularMesh(pugi::xml_node node) : StructuredMesh {node}
|
|||
|
||||
width_ = get_node_xarray<double>(node, "width");
|
||||
|
||||
// Check to ensure width has same dimensions
|
||||
auto n = width_.size();
|
||||
if (n != lower_left_.size()) {
|
||||
fatal_error("Number of entries on <width> must be the same as "
|
||||
"the number of entries on <lower_left>.");
|
||||
}
|
||||
|
||||
// Check for negative widths
|
||||
if (xt::any(width_ < 0.0)) {
|
||||
fatal_error("Cannot have a negative <width> on a tally mesh.");
|
||||
}
|
||||
|
||||
// Set width and upper right coordinate
|
||||
upper_right_ = xt::eval(lower_left_ + shape * width_);
|
||||
|
||||
} else if (check_for_node(node, "upper_right")) {
|
||||
|
||||
upper_right_ = get_node_xarray<double>(node, "upper_right");
|
||||
|
||||
// Check to ensure width has same dimensions
|
||||
auto n = upper_right_.size();
|
||||
if (n != lower_left_.size()) {
|
||||
fatal_error("Number of entries on <upper_right> must be the "
|
||||
"same as the number of entries on <lower_left>.");
|
||||
}
|
||||
|
||||
// Check that upper-right is above lower-left
|
||||
if (xt::any(upper_right_ < lower_left_)) {
|
||||
fatal_error("The <upper_right> coordinates must be greater than "
|
||||
"the <lower_left> coordinates on a tally mesh.");
|
||||
}
|
||||
|
||||
// Set width
|
||||
width_ = xt::eval((upper_right_ - lower_left_) / shape);
|
||||
} else {
|
||||
fatal_error("Must specify either <upper_right> or <width> on a mesh.");
|
||||
}
|
||||
|
||||
// Set material volumes
|
||||
volume_frac_ = 1.0 / xt::prod(shape)();
|
||||
if (int err = set_grid()) {
|
||||
fatal_error(openmc_err_msg);
|
||||
}
|
||||
}
|
||||
|
||||
element_volume_ = 1.0;
|
||||
for (int i = 0; i < n_dimension_; i++) {
|
||||
element_volume_ *= width_[i];
|
||||
RegularMesh::RegularMesh(hid_t group) : StructuredMesh {group}
|
||||
{
|
||||
// Determine number of dimensions for mesh
|
||||
if (!object_exists(group, "dimension")) {
|
||||
fatal_error("Must specify <dimension> on a regular mesh.");
|
||||
}
|
||||
|
||||
xt::xtensor<int, 1> shape;
|
||||
read_dataset(group, "dimension", shape);
|
||||
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 for lower-left coordinates
|
||||
if (object_exists(group, "lower_left")) {
|
||||
// Read mesh lower-left corner location
|
||||
read_dataset(group, "lower_left", lower_left_);
|
||||
} else {
|
||||
fatal_error("Must specify lower_left dataset on a mesh.");
|
||||
}
|
||||
|
||||
if (object_exists(group, "upper_right")) {
|
||||
|
||||
read_dataset(group, "upper_right", upper_right_);
|
||||
|
||||
} else {
|
||||
fatal_error("Must specify either upper_right dataset on a mesh.");
|
||||
}
|
||||
|
||||
if (int err = set_grid()) {
|
||||
fatal_error(openmc_err_msg);
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -1343,6 +1493,19 @@ RectilinearMesh::RectilinearMesh(pugi::xml_node node) : StructuredMesh {node}
|
|||
}
|
||||
}
|
||||
|
||||
RectilinearMesh::RectilinearMesh(hid_t group) : StructuredMesh {group}
|
||||
{
|
||||
n_dimension_ = 3;
|
||||
|
||||
read_dataset(group, "x_grid", grid_[0]);
|
||||
read_dataset(group, "y_grid", grid_[1]);
|
||||
read_dataset(group, "z_grid", grid_[2]);
|
||||
|
||||
if (int err = set_grid()) {
|
||||
fatal_error(openmc_err_msg);
|
||||
}
|
||||
}
|
||||
|
||||
const std::string RectilinearMesh::mesh_type = "rectilinear";
|
||||
|
||||
std::string RectilinearMesh::get_mesh_type() const
|
||||
|
|
@ -1478,6 +1641,19 @@ CylindricalMesh::CylindricalMesh(pugi::xml_node node)
|
|||
}
|
||||
}
|
||||
|
||||
CylindricalMesh::CylindricalMesh(hid_t group) : PeriodicStructuredMesh {group}
|
||||
{
|
||||
n_dimension_ = 3;
|
||||
read_dataset(group, "r_grid", grid_[0]);
|
||||
read_dataset(group, "phi_grid", grid_[1]);
|
||||
read_dataset(group, "z_grid", grid_[2]);
|
||||
read_dataset(group, "origin", origin_);
|
||||
|
||||
if (int err = set_grid()) {
|
||||
fatal_error(openmc_err_msg);
|
||||
}
|
||||
}
|
||||
|
||||
const std::string CylindricalMesh::mesh_type = "cylindrical";
|
||||
|
||||
std::string CylindricalMesh::get_mesh_type() const
|
||||
|
|
@ -1756,6 +1932,20 @@ SphericalMesh::SphericalMesh(pugi::xml_node node)
|
|||
}
|
||||
}
|
||||
|
||||
SphericalMesh::SphericalMesh(hid_t group) : PeriodicStructuredMesh {group}
|
||||
{
|
||||
n_dimension_ = 3;
|
||||
|
||||
read_dataset(group, "r_grid", grid_[0]);
|
||||
read_dataset(group, "theta_grid", grid_[1]);
|
||||
read_dataset(group, "phi_grid", grid_[2]);
|
||||
read_dataset(group, "origin", origin_);
|
||||
|
||||
if (int err = set_grid()) {
|
||||
fatal_error(openmc_err_msg);
|
||||
}
|
||||
}
|
||||
|
||||
const std::string SphericalMesh::mesh_type = "spherical";
|
||||
|
||||
std::string SphericalMesh::get_mesh_type() const
|
||||
|
|
@ -2520,6 +2710,11 @@ MOABMesh::MOABMesh(pugi::xml_node node) : UnstructuredMesh(node)
|
|||
initialize();
|
||||
}
|
||||
|
||||
MOABMesh::MOABMesh(hid_t group) : UnstructuredMesh(group)
|
||||
{
|
||||
initialize();
|
||||
}
|
||||
|
||||
MOABMesh::MOABMesh(const std::string& filename, double length_multiplier)
|
||||
: UnstructuredMesh()
|
||||
{
|
||||
|
|
@ -3228,6 +3423,15 @@ LibMesh::LibMesh(pugi::xml_node node) : UnstructuredMesh(node)
|
|||
initialize();
|
||||
}
|
||||
|
||||
LibMesh::LibMesh(hid_t group) : UnstructuredMesh(group)
|
||||
{
|
||||
// filename_ and length_multiplier_ will already be set by the
|
||||
// UnstructuredMesh constructor
|
||||
set_mesh_pointer_from_filename(filename_);
|
||||
set_length_multiplier(length_multiplier_);
|
||||
initialize();
|
||||
}
|
||||
|
||||
// create the mesh from a pointer to a libMesh Mesh
|
||||
LibMesh::LibMesh(libMesh::MeshBase& input_mesh, double length_multiplier)
|
||||
{
|
||||
|
|
@ -3618,34 +3822,51 @@ void read_meshes(pugi::xml_node root)
|
|||
mesh_lib = get_node_value(node, "library", true, true);
|
||||
}
|
||||
|
||||
// Read mesh and add to vector
|
||||
if (mesh_type == RegularMesh::mesh_type) {
|
||||
model::meshes.push_back(make_unique<RegularMesh>(node));
|
||||
} else if (mesh_type == RectilinearMesh::mesh_type) {
|
||||
model::meshes.push_back(make_unique<RectilinearMesh>(node));
|
||||
} else if (mesh_type == CylindricalMesh::mesh_type) {
|
||||
model::meshes.push_back(make_unique<CylindricalMesh>(node));
|
||||
} else if (mesh_type == SphericalMesh::mesh_type) {
|
||||
model::meshes.push_back(make_unique<SphericalMesh>(node));
|
||||
#ifdef OPENMC_DAGMC_ENABLED
|
||||
} else if (mesh_type == UnstructuredMesh::mesh_type &&
|
||||
mesh_lib == MOABMesh::mesh_lib_type) {
|
||||
model::meshes.push_back(make_unique<MOABMesh>(node));
|
||||
#endif
|
||||
#ifdef OPENMC_LIBMESH_ENABLED
|
||||
} else if (mesh_type == UnstructuredMesh::mesh_type &&
|
||||
mesh_lib == LibMesh::mesh_lib_type) {
|
||||
model::meshes.push_back(make_unique<LibMesh>(node));
|
||||
#endif
|
||||
} else if (mesh_type == UnstructuredMesh::mesh_type) {
|
||||
fatal_error("Unstructured mesh support is not enabled or the mesh "
|
||||
"library is invalid.");
|
||||
} else {
|
||||
fatal_error("Invalid mesh type: " + mesh_type);
|
||||
Mesh::create(node, mesh_type, mesh_lib);
|
||||
}
|
||||
}
|
||||
|
||||
void read_meshes(hid_t group)
|
||||
{
|
||||
std::unordered_set<int> mesh_ids;
|
||||
|
||||
std::vector<int> ids;
|
||||
read_attribute(group, "ids", ids);
|
||||
|
||||
for (auto id : ids) {
|
||||
|
||||
// Check to make sure multiple meshes in the same file don't share IDs
|
||||
if (contains(mesh_ids, id)) {
|
||||
fatal_error(fmt::format("Two or more meshes use the same unique ID "
|
||||
"'{}' in the same HDF5 input file",
|
||||
id));
|
||||
}
|
||||
mesh_ids.insert(id);
|
||||
|
||||
// If we've already read a mesh with the same ID in a *different* file,
|
||||
// assume it is the same here
|
||||
if (model::mesh_map.find(id) != model::mesh_map.end()) {
|
||||
warning(fmt::format("Mesh with ID={} appears in multiple files.", id));
|
||||
continue;
|
||||
}
|
||||
|
||||
// Map ID to position in vector
|
||||
model::mesh_map[model::meshes.back()->id_] = model::meshes.size() - 1;
|
||||
std::string name = fmt::format("mesh {}", id);
|
||||
hid_t mesh_group = open_group(group, name.c_str());
|
||||
|
||||
std::string mesh_type;
|
||||
if (object_exists(mesh_group, "type")) {
|
||||
read_dataset(mesh_group, "type", mesh_type);
|
||||
} else {
|
||||
mesh_type = "regular";
|
||||
}
|
||||
|
||||
// determine the mesh library to use
|
||||
std::string mesh_lib;
|
||||
if (object_exists(mesh_group, "library")) {
|
||||
read_dataset(mesh_group, "library", mesh_lib);
|
||||
}
|
||||
|
||||
Mesh::create(mesh_group, mesh_type, mesh_lib);
|
||||
}
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -1328,6 +1328,10 @@ extern "C" int openmc_weight_windows_import(const char* filename)
|
|||
|
||||
hid_t weight_windows_group = open_group(ww_file, "weight_windows");
|
||||
|
||||
hid_t mesh_group = open_group(ww_file, "meshes");
|
||||
|
||||
read_meshes(mesh_group);
|
||||
|
||||
std::vector<std::string> names = group_names(weight_windows_group);
|
||||
|
||||
for (const auto& name : names) {
|
||||
|
|
|
|||
|
|
@ -5,6 +5,7 @@ set(TEST_NAMES
|
|||
test_interpolate
|
||||
test_math
|
||||
test_mcpl_stat_sum
|
||||
test_mesh
|
||||
# Add additional unit test files here
|
||||
)
|
||||
|
||||
|
|
|
|||
257
tests/cpp_unit_tests/test_mesh.cpp
Normal file
257
tests/cpp_unit_tests/test_mesh.cpp
Normal file
|
|
@ -0,0 +1,257 @@
|
|||
#include <cstdio>
|
||||
#include <iostream>
|
||||
#include <string>
|
||||
|
||||
#include <catch2/catch_test_macros.hpp>
|
||||
#include <pugixml.hpp>
|
||||
|
||||
#include "openmc/hdf5_interface.h"
|
||||
#include "openmc/mesh.h"
|
||||
|
||||
using namespace openmc;
|
||||
|
||||
TEST_CASE("Test mesh hdf5 roundtrip - regular")
|
||||
{
|
||||
// The XML data as a string
|
||||
std::string xml_string = R"(
|
||||
<mesh id="1">
|
||||
<dimension>3 4 5</dimension>
|
||||
<lower_left>-2 -3 -5</lower_left>
|
||||
<upper_right>2 3 5</upper_right>
|
||||
</mesh>
|
||||
)";
|
||||
|
||||
// Create a pugixml document object
|
||||
pugi::xml_document doc;
|
||||
|
||||
// Load the XML from the string
|
||||
pugi::xml_parse_result result = doc.load_string(xml_string.c_str());
|
||||
|
||||
pugi::xml_node root = doc.child("mesh");
|
||||
|
||||
auto mesh = RegularMesh(root);
|
||||
|
||||
hid_t file_id = file_open("mesh.h5", 'w');
|
||||
|
||||
mesh.to_hdf5(file_id);
|
||||
|
||||
file_close(file_id);
|
||||
|
||||
hid_t file_id2 = file_open("mesh.h5", 'r');
|
||||
|
||||
hid_t group = open_group(file_id2, "mesh 1");
|
||||
|
||||
auto mesh2 = RegularMesh(group);
|
||||
|
||||
file_close(file_id2);
|
||||
|
||||
remove("mesh.h5");
|
||||
|
||||
REQUIRE(mesh2.shape_ == mesh.shape_);
|
||||
|
||||
REQUIRE(mesh2.lower_left() == mesh.lower_left());
|
||||
|
||||
REQUIRE(mesh2.upper_right() == mesh.upper_right());
|
||||
}
|
||||
|
||||
TEST_CASE("Test mesh hdf5 roundtrip - rectilinear")
|
||||
{
|
||||
// The XML data as a string
|
||||
std::string xml_string = R"(
|
||||
<mesh id="1" type="rectilinear">
|
||||
<x_grid>0.0 1.0 5.0 10.0</x_grid>
|
||||
<y_grid>-10.0 -5.0 0.0</y_grid>
|
||||
<z_grid>-100.0 0.0 100.0</z_grid>
|
||||
</mesh>
|
||||
)";
|
||||
|
||||
// Create a pugixml document object
|
||||
pugi::xml_document doc;
|
||||
|
||||
// Load the XML from the string
|
||||
pugi::xml_parse_result result = doc.load_string(xml_string.c_str());
|
||||
|
||||
pugi::xml_node root = doc.child("mesh");
|
||||
|
||||
auto mesh = RectilinearMesh(root);
|
||||
|
||||
hid_t file_id = file_open("mesh.h5", 'w');
|
||||
|
||||
mesh.to_hdf5(file_id);
|
||||
|
||||
file_close(file_id);
|
||||
|
||||
hid_t file_id2 = file_open("mesh.h5", 'r');
|
||||
|
||||
hid_t group = open_group(file_id2, "mesh 1");
|
||||
|
||||
auto mesh2 = RectilinearMesh(group);
|
||||
|
||||
file_close(file_id2);
|
||||
|
||||
remove("mesh.h5");
|
||||
|
||||
REQUIRE(mesh2.shape_ == mesh.shape_);
|
||||
|
||||
REQUIRE(mesh2.grid_ == mesh.grid_);
|
||||
}
|
||||
|
||||
TEST_CASE("Test mesh hdf5 roundtrip - cylindrical")
|
||||
{
|
||||
// The XML data as a string
|
||||
std::string xml_string = R"(
|
||||
<mesh id="1" type="cylindrical">
|
||||
<r_grid>0.1 0.2 0.5 1.0</r_grid>
|
||||
<phi_grid>0.0 6.283185307179586</phi_grid>
|
||||
<z_grid>0.1 0.2 0.4 0.6 1.0</z_grid>
|
||||
<origin>0 0 0</origin>
|
||||
</mesh>
|
||||
)";
|
||||
|
||||
// Create a pugixml document object
|
||||
pugi::xml_document doc;
|
||||
|
||||
// Load the XML from the string
|
||||
pugi::xml_parse_result result = doc.load_string(xml_string.c_str());
|
||||
|
||||
pugi::xml_node root = doc.child("mesh");
|
||||
|
||||
auto mesh = CylindricalMesh(root);
|
||||
|
||||
hid_t file_id = file_open("mesh.h5", 'w');
|
||||
|
||||
mesh.to_hdf5(file_id);
|
||||
|
||||
file_close(file_id);
|
||||
|
||||
hid_t file_id2 = file_open("mesh.h5", 'r');
|
||||
|
||||
hid_t group = open_group(file_id2, "mesh 1");
|
||||
|
||||
auto mesh2 = CylindricalMesh(group);
|
||||
|
||||
file_close(file_id2);
|
||||
|
||||
remove("mesh.h5");
|
||||
|
||||
REQUIRE(mesh2.shape_ == mesh.shape_);
|
||||
|
||||
REQUIRE(mesh2.grid_ == mesh.grid_);
|
||||
}
|
||||
|
||||
TEST_CASE("Test mesh hdf5 roundtrip - spherical")
|
||||
{
|
||||
// The XML data as a string
|
||||
std::string xml_string = R"(
|
||||
<mesh id="1" type="spherical">
|
||||
<r_grid>0.1 0.2 0.5 1.0</r_grid>
|
||||
<theta_grid>0.0 3.141592653589793</theta_grid>
|
||||
<phi_grid>0.0 6.283185307179586</phi_grid>
|
||||
<origin>0.0 0.0 0.0</origin>
|
||||
</mesh>'
|
||||
)";
|
||||
|
||||
// Create a pugixml document object
|
||||
pugi::xml_document doc;
|
||||
|
||||
// Load the XML from the string
|
||||
pugi::xml_parse_result result = doc.load_string(xml_string.c_str());
|
||||
|
||||
pugi::xml_node root = doc.child("mesh");
|
||||
|
||||
auto mesh = SphericalMesh(root);
|
||||
|
||||
hid_t file_id = file_open("mesh.h5", 'w');
|
||||
|
||||
mesh.to_hdf5(file_id);
|
||||
|
||||
file_close(file_id);
|
||||
|
||||
hid_t file_id2 = file_open("mesh.h5", 'r');
|
||||
|
||||
hid_t group = open_group(file_id2, "mesh 1");
|
||||
|
||||
auto mesh2 = SphericalMesh(group);
|
||||
|
||||
file_close(file_id2);
|
||||
|
||||
remove("mesh.h5");
|
||||
|
||||
REQUIRE(mesh2.shape_ == mesh.shape_);
|
||||
|
||||
REQUIRE(mesh2.grid_ == mesh.grid_);
|
||||
}
|
||||
|
||||
TEST_CASE("Test multiple meshes HDF5 roundtrip - spherical")
|
||||
{
|
||||
// The XML data as a string
|
||||
std::string xml_string = R"(
|
||||
<meshes>
|
||||
<mesh id="1" type="spherical">
|
||||
<r_grid>0.1 0.2 0.5 1.0</r_grid>
|
||||
<theta_grid>0.0 3.141592653589793</theta_grid>
|
||||
<phi_grid>0.0 6.283185307179586</phi_grid>
|
||||
<origin>0.0 0.0 0.0</origin>
|
||||
</mesh>
|
||||
<mesh id="2">
|
||||
<dimension>3 4 5</dimension>
|
||||
<lower_left>-2 -3 -5</lower_left>
|
||||
<upper_right>2 3 5</upper_right>
|
||||
</mesh>
|
||||
</meshes>
|
||||
)";
|
||||
|
||||
// Create a pugixml document object
|
||||
pugi::xml_document doc;
|
||||
|
||||
// Load the XML from the string
|
||||
pugi::xml_parse_result result = doc.load_string(xml_string.c_str());
|
||||
|
||||
pugi::xml_node root = doc.child("meshes");
|
||||
|
||||
read_meshes(root);
|
||||
|
||||
const auto spherical_mesh_xml =
|
||||
dynamic_cast<SphericalMesh*>(model::meshes[0].get());
|
||||
const auto regular_mesh_xml =
|
||||
dynamic_cast<RegularMesh*>(model::meshes[1].get());
|
||||
|
||||
hid_t file_id = file_open("meshes.h5", 'w');
|
||||
|
||||
hid_t root_group = create_group(file_id, "root");
|
||||
|
||||
open_group(file_id, "root");
|
||||
|
||||
meshes_to_hdf5(root_group);
|
||||
|
||||
close_group(root_group);
|
||||
|
||||
file_close(file_id);
|
||||
|
||||
hid_t file_id2 = file_open("meshes.h5", 'r');
|
||||
|
||||
hid_t root_group_read = open_group(file_id2, "root");
|
||||
|
||||
hid_t mesh_group_read = open_group(root_group_read, "meshes");
|
||||
|
||||
read_meshes(mesh_group_read);
|
||||
|
||||
// increment mesh IDs to avoid collision during read
|
||||
for (auto& mesh : model::meshes) {
|
||||
mesh->set_id(mesh->id() + 10);
|
||||
}
|
||||
|
||||
const auto spherical_mesh_hdf5 = dynamic_cast<SphericalMesh*>(
|
||||
model::meshes[model::mesh_map[spherical_mesh_xml->id_]].get());
|
||||
const auto regular_mesh_hdf5 = dynamic_cast<RegularMesh*>(
|
||||
model::meshes[model::mesh_map[regular_mesh_xml->id_]].get());
|
||||
|
||||
remove("meshes.h5");
|
||||
|
||||
REQUIRE(spherical_mesh_hdf5->shape_ == spherical_mesh_xml->shape_);
|
||||
REQUIRE(spherical_mesh_hdf5->grid_ == spherical_mesh_xml->grid_);
|
||||
|
||||
REQUIRE(regular_mesh_hdf5->shape_ == regular_mesh_xml->shape_);
|
||||
REQUIRE(regular_mesh_hdf5->lower_left() == regular_mesh_xml->lower_left());
|
||||
REQUIRE(regular_mesh_hdf5->upper_right() == regular_mesh_xml->upper_right());
|
||||
}
|
||||
Loading…
Add table
Add a link
Reference in a new issue