From 39a2d46e2636707d11f951fb4f206c2155c77198 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 18 Aug 2024 23:09:32 -0500 Subject: [PATCH] Implement bounding_box operation for meshes (#3119) --- cmake/Modules/FindLIBMESH.cmake | 2 +- include/openmc/bounding_box.h | 61 +++++++++++++++++++++++++++++++++ include/openmc/cell.h | 1 + include/openmc/mesh.h | 52 +++++++++++++++++++++++++--- include/openmc/surface.h | 50 +-------------------------- include/openmc/universe.h | 1 + openmc/bounding_box.py | 3 +- openmc/lib/mesh.py | 30 ++++++++++++++++ src/mesh.cpp | 60 +++++++++++++++++++++++++++++--- src/particle.cpp | 2 +- tests/unit_tests/test_lib.py | 21 ++++++++++++ 11 files changed, 223 insertions(+), 60 deletions(-) create mode 100644 include/openmc/bounding_box.h diff --git a/cmake/Modules/FindLIBMESH.cmake b/cmake/Modules/FindLIBMESH.cmake index df5b5b9d8..048dfc2a8 100644 --- a/cmake/Modules/FindLIBMESH.cmake +++ b/cmake/Modules/FindLIBMESH.cmake @@ -14,7 +14,7 @@ if(DEFINED ENV{METHOD}) message(STATUS "Using environment variable METHOD to determine libMesh build: ${LIBMESH_PC_FILE}") endif() -include(FindPkgConfig) +find_package(PkgConfig REQUIRED) set(ENV{PKG_CONFIG_PATH} "$ENV{PKG_CONFIG_PATH}:${LIBMESH_PC}") set(PKG_CONFIG_USE_CMAKE_PREFIX_PATH True) pkg_check_modules(LIBMESH REQUIRED ${LIBMESH_PC_FILE}>=1.7.0 IMPORTED_TARGET) diff --git a/include/openmc/bounding_box.h b/include/openmc/bounding_box.h new file mode 100644 index 000000000..40f603583 --- /dev/null +++ b/include/openmc/bounding_box.h @@ -0,0 +1,61 @@ +#ifndef OPENMC_BOUNDING_BOX_H +#define OPENMC_BOUNDING_BOX_H + +#include // for min, max + +#include "openmc/constants.h" + +namespace openmc { + +//============================================================================== +//! Coordinates for an axis-aligned cuboid that bounds a geometric object. +//============================================================================== + +struct BoundingBox { + double xmin = -INFTY; + double xmax = INFTY; + double ymin = -INFTY; + double ymax = INFTY; + double zmin = -INFTY; + double zmax = INFTY; + + inline BoundingBox operator&(const BoundingBox& other) + { + BoundingBox result = *this; + return result &= other; + } + + inline BoundingBox operator|(const BoundingBox& other) + { + BoundingBox result = *this; + return result |= other; + } + + // intersect operator + inline BoundingBox& operator&=(const BoundingBox& other) + { + xmin = std::max(xmin, other.xmin); + xmax = std::min(xmax, other.xmax); + ymin = std::max(ymin, other.ymin); + ymax = std::min(ymax, other.ymax); + zmin = std::max(zmin, other.zmin); + zmax = std::min(zmax, other.zmax); + return *this; + } + + // union operator + inline BoundingBox& operator|=(const BoundingBox& other) + { + xmin = std::min(xmin, other.xmin); + xmax = std::max(xmax, other.xmax); + ymin = std::min(ymin, other.ymin); + ymax = std::max(ymax, other.ymax); + zmin = std::min(zmin, other.zmin); + zmax = std::max(zmax, other.zmax); + return *this; + } +}; + +} // namespace openmc + +#endif diff --git a/include/openmc/cell.h b/include/openmc/cell.h index d78614f05..70843140b 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -12,6 +12,7 @@ #include "pugixml.hpp" #include +#include "openmc/bounding_box.h" #include "openmc/constants.h" #include "openmc/memory.h" // for unique_ptr #include "openmc/neighbor_list.h" diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 3917e6368..6f3f2b0a4 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -11,6 +11,7 @@ #include "xtensor/xtensor.hpp" #include +#include "openmc/bounding_box.h" #include "openmc/error.h" #include "openmc/memory.h" // for unique_ptr #include "openmc/particle.h" @@ -185,9 +186,24 @@ public: vector material_volumes( int n_sample, int bin, uint64_t* seed) const; + //! Determine bounding box of mesh + // + //! \return Bounding box of mesh + BoundingBox bounding_box() const + { + auto ll = this->lower_left(); + auto ur = this->upper_right(); + return {ll.x, ur.x, ll.y, ur.y, ll.z, ur.z}; + } + + virtual Position lower_left() const = 0; + virtual Position upper_right() const = 0; + // Data members - int id_ {-1}; //!< User-specified ID - int n_dimension_ {-1}; //!< Number of dimensions + xt::xtensor lower_left_; //!< Lower-left coordinates of mesh + xt::xtensor upper_right_; //!< Upper-right coordinates of mesh + int id_ {-1}; //!< User-specified ID + int n_dimension_ {-1}; //!< Number of dimensions }; class StructuredMesh : public Mesh { @@ -325,14 +341,30 @@ public: return this->volume(get_indices_from_bin(bin)); } + Position lower_left() const override + { + int n = lower_left_.size(); + Position ll {lower_left_[0], 0.0, 0.0}; + ll.y = (n >= 2) ? lower_left_[1] : -INFTY; + ll.z = (n == 3) ? lower_left_[2] : -INFTY; + return ll; + }; + + Position upper_right() const override + { + int n = upper_right_.size(); + Position ur {upper_right_[0], 0.0, 0.0}; + ur.y = (n >= 2) ? upper_right_[1] : INFTY; + ur.z = (n == 3) ? upper_right_[2] : INFTY; + return ur; + }; + //! Get the volume of a specified element //! \param[in] ijk Mesh index to return the volume for //! \return Volume of the bin virtual double volume(const MeshIndex& ijk) const = 0; // Data members - xt::xtensor lower_left_; //!< Lower-left coordinates of mesh - xt::xtensor upper_right_; //!< Upper-right coordinates of mesh std::array shape_; //!< Number of mesh elements in each dimension protected: @@ -655,6 +687,15 @@ public: ElementType element_type(int bin) const; + Position lower_left() const override + { + return {lower_left_[0], lower_left_[1], lower_left_[2]}; + } + Position upper_right() const override + { + return {upper_right_[0], upper_right_[1], upper_right_[2]}; + } + protected: //! Set the length multiplier to apply to each point in the mesh void set_length_multiplier(const double length_multiplier); @@ -672,6 +713,9 @@ protected: -1.0}; //!< Multiplicative factor applied to mesh coordinates std::string options_; //!< Options for search data structures + //! Determine lower-left and upper-right bounds of mesh + void determine_bounds(); + private: //! Setup method for the mesh. Builds data structures, //! sets up element mapping, creates bounding boxes, etc. diff --git a/include/openmc/surface.h b/include/openmc/surface.h index 350775123..af235301c 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -9,6 +9,7 @@ #include "pugixml.hpp" #include "openmc/boundary_condition.h" +#include "openmc/bounding_box.h" #include "openmc/constants.h" #include "openmc/memory.h" // for unique_ptr #include "openmc/particle.h" @@ -28,55 +29,6 @@ extern std::unordered_map surface_map; extern vector> surfaces; } // namespace model -//============================================================================== -//! Coordinates for an axis-aligned cuboid that bounds a geometric object. -//============================================================================== - -struct BoundingBox { - double xmin = -INFTY; - double xmax = INFTY; - double ymin = -INFTY; - double ymax = INFTY; - double zmin = -INFTY; - double zmax = INFTY; - - inline BoundingBox operator&(const BoundingBox& other) - { - BoundingBox result = *this; - return result &= other; - } - - inline BoundingBox operator|(const BoundingBox& other) - { - BoundingBox result = *this; - return result |= other; - } - - // intersect operator - inline BoundingBox& operator&=(const BoundingBox& other) - { - xmin = std::max(xmin, other.xmin); - xmax = std::min(xmax, other.xmax); - ymin = std::max(ymin, other.ymin); - ymax = std::min(ymax, other.ymax); - zmin = std::max(zmin, other.zmin); - zmax = std::min(zmax, other.zmax); - return *this; - } - - // union operator - inline BoundingBox& operator|=(const BoundingBox& other) - { - xmin = std::min(xmin, other.xmin); - xmax = std::max(xmax, other.xmax); - ymin = std::min(ymin, other.ymin); - ymax = std::max(ymax, other.ymax); - zmin = std::min(zmin, other.zmin); - zmax = std::max(zmax, other.zmax); - return *this; - } -}; - //============================================================================== //! A geometry primitive used to define regions of 3D space. //============================================================================== diff --git a/include/openmc/universe.h b/include/openmc/universe.h index 26f33cb38..9fea06bcc 100644 --- a/include/openmc/universe.h +++ b/include/openmc/universe.h @@ -1,6 +1,7 @@ #ifndef OPENMC_UNIVERSE_H #define OPENMC_UNIVERSE_H +#include "openmc/bounding_box.h" #include "openmc/cell.h" namespace openmc { diff --git a/openmc/bounding_box.py b/openmc/bounding_box.py index 5f3d3a6cf..f0dc06a4a 100644 --- a/openmc/bounding_box.py +++ b/openmc/bounding_box.py @@ -42,7 +42,8 @@ class BoundingBox: def __repr__(self) -> str: return "BoundingBox(lower_left={}, upper_right={})".format( - tuple(self.lower_left), tuple(self.upper_right)) + tuple(float(x) for x in self.lower_left), + tuple(float(x) for x in self.upper_right)) def __getitem__(self, key) -> np.ndarray: return self._bounds[key] diff --git a/openmc/lib/mesh.py b/openmc/lib/mesh.py index 78566d449..16bec0198 100644 --- a/openmc/lib/mesh.py +++ b/openmc/lib/mesh.py @@ -2,6 +2,7 @@ from collections.abc import Mapping, Sequence from ctypes import (c_int, c_int32, c_char_p, c_double, POINTER, Structure, create_string_buffer, c_uint64, c_size_t) from random import getrandbits +import sys from weakref import WeakValueDictionary import numpy as np @@ -13,6 +14,7 @@ from .core import _FortranObjectWithID from .error import _error_handler from .material import Material from .plot import _Position +from ..bounding_box import BoundingBox __all__ = [ 'Mesh', 'RegularMesh', 'RectilinearMesh', 'CylindricalMesh', @@ -44,6 +46,10 @@ _dll.openmc_mesh_get_n_elements.errcheck = _error_handler _dll.openmc_mesh_get_volumes.argtypes = [c_int32, POINTER(c_double)] _dll.openmc_mesh_get_volumes.restype = c_int _dll.openmc_mesh_get_volumes.errcheck = _error_handler +_dll.openmc_mesh_bounding_box.argtypes = [ + c_int32, POINTER(c_double), POINTER(c_double)] +_dll.openmc_mesh_bounding_box.restype = c_int +_dll.openmc_mesh_bounding_box.errcheck = _error_handler _dll.openmc_mesh_material_volumes.argtypes = [ c_int32, c_int, c_int, c_int, POINTER(_MaterialVolume), POINTER(c_int), POINTER(c_uint64)] @@ -166,6 +172,22 @@ class Mesh(_FortranObjectWithID): self._index, volumes.ctypes.data_as(POINTER(c_double))) return volumes + @property + def bounding_box(self) -> BoundingBox: + inf = sys.float_info.max + ll = np.zeros(3) + ur = np.zeros(3) + _dll.openmc_mesh_bounding_box( + self._index, + ll.ctypes.data_as(POINTER(c_double)), + ur.ctypes.data_as(POINTER(c_double)) + ) + ll[ll == inf] = np.inf + ur[ur == inf] = np.inf + ll[ll == -inf] = -np.inf + ur[ur == -inf] = -np.inf + return BoundingBox(ll, ur) + def material_volumes( self, n_samples: int = 10_000, @@ -292,6 +314,8 @@ class RegularMesh(Mesh): Total number of mesh elements. volumes : numpy.ndarray Volume of each mesh element in [cm^3] + bounding_box : openmc.BoundingBox + Axis-aligned bounding box of the mesh """ mesh_type = 'regular' @@ -378,6 +402,8 @@ class RectilinearMesh(Mesh): Total number of mesh elements. volumes : numpy.ndarray Volume of each mesh element in [cm^3] + bounding_box : openmc.BoundingBox + Axis-aligned bounding box of the mesh """ mesh_type = 'rectilinear' @@ -481,6 +507,8 @@ class CylindricalMesh(Mesh): Total number of mesh elements. volumes : numpy.ndarray Volume of each mesh element in [cm^3] + bounding_box : openmc.BoundingBox + Axis-aligned bounding box of the mesh """ mesh_type = 'cylindrical' @@ -584,6 +612,8 @@ class SphericalMesh(Mesh): Total number of mesh elements. volumes : numpy.ndarray Volume of each mesh element in [cm^3] + bounding_box : openmc.BoundingBox + Axis-aligned bounding box of the mesh """ mesh_type = 'spherical' diff --git a/src/mesh.cpp b/src/mesh.cpp index c57dd5dcc..0f8b4b14e 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -319,6 +319,28 @@ UnstructuredMesh::UnstructuredMesh(pugi::xml_node node) : Mesh(node) } } +void UnstructuredMesh::determine_bounds() +{ + double xmin = INFTY; + double ymin = INFTY; + double zmin = INFTY; + double xmax = -INFTY; + double ymax = -INFTY; + double zmax = -INFTY; + int n = this->n_vertices(); + for (int i = 0; i < n; ++i) { + auto v = this->vertex(i); + xmin = std::min(v.x, xmin); + ymin = std::min(v.y, ymin); + zmin = std::min(v.z, zmin); + xmax = std::max(v.x, xmax); + ymax = std::max(v.y, ymax); + zmax = std::max(v.z, zmax); + } + lower_left_ = {xmin, ymin, zmin}; + upper_right_ = {xmax, ymax, zmax}; +} + Position UnstructuredMesh::sample_tet( std::array coords, uint64_t* seed) const { @@ -1372,8 +1394,10 @@ int CylindricalMesh::set_grid() full_phi_ = (grid_[1].front() == 0.0) && (grid_[1].back() == 2.0 * PI); - lower_left_ = {grid_[0].front(), grid_[1].front(), grid_[2].front()}; - upper_right_ = {grid_[0].back(), grid_[1].back(), grid_[2].back()}; + lower_left_ = {origin_[0] - grid_[0].back(), origin_[1] - grid_[0].back(), + origin_[2] + grid_[2].front()}; + upper_right_ = {origin_[0] + grid_[0].back(), origin_[1] + grid_[0].back(), + origin_[2] + grid_[2].back()}; return 0; } @@ -1687,8 +1711,9 @@ int SphericalMesh::set_grid() full_theta_ = (grid_[1].front() == 0.0) && (grid_[1].back() == PI); full_phi_ = (grid_[2].front() == 0.0) && (grid_[2].back() == 2 * PI); - lower_left_ = {grid_[0].front(), grid_[1].front(), grid_[2].front()}; - upper_right_ = {grid_[0].back(), grid_[1].back(), grid_[2].back()}; + double r = grid_[0].back(); + lower_left_ = {origin_[0] - r, origin_[1] - r, origin_[2] - r}; + upper_right_ = {origin_[0] + r, origin_[1] + r, origin_[2] + r}; return 0; } @@ -1899,6 +1924,26 @@ extern "C" int openmc_mesh_get_volumes(int32_t index, double* volumes) return 0; } +//! Get the bounding box of a mesh +extern "C" int openmc_mesh_bounding_box(int32_t index, double* ll, double* ur) +{ + if (int err = check_mesh(index)) + return err; + + BoundingBox bbox = model::meshes[index]->bounding_box(); + + // set lower left corner values + ll[0] = bbox.xmin; + ll[1] = bbox.ymin; + ll[2] = bbox.zmin; + + // set upper right corner values + ur[0] = bbox.xmax; + ur[1] = bbox.ymax; + ur[2] = bbox.zmax; + return 0; +} + extern "C" int openmc_mesh_material_volumes(int32_t index, int n_sample, int bin, int result_size, void* result, int* hits, uint64_t* seed) { @@ -2265,6 +2310,9 @@ void MOABMesh::initialize() } } } + + // Determine bounds of mesh + this->determine_bounds(); } void MOABMesh::prepare_for_tallies() @@ -2952,6 +3000,10 @@ void LibMesh::initialize() // bounding box for the mesh for quick rejection checks bbox_ = libMesh::MeshTools::create_bounding_box(*m_); + libMesh::Point ll = bbox_.min(); + libMesh::Point ur = bbox_.max(); + lower_left_ = {ll(0), ll(1), ll(2)}; + upper_right_ = {ur(0), ur(1), ur(2)}; } // Sample position within a tet for LibMesh type tets diff --git a/src/particle.cpp b/src/particle.cpp index 7d30e26bd..64c50c943 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -296,7 +296,7 @@ void Particle::event_cross_surface() if (surf->surf_source_ && surf->bc_) { add_surf_source_to_bank(*this, *surf); } - cross_surface(*surf); + this->cross_surface(*surf); // If no BC, add particle to surface source after crossing surface if (surf->surf_source_ && !surf->bc_) { add_surf_source_to_bank(*this, *surf); diff --git a/tests/unit_tests/test_lib.py b/tests/unit_tests/test_lib.py index e3c7ce3b7..64c16c238 100644 --- a/tests/unit_tests/test_lib.py +++ b/tests/unit_tests/test_lib.py @@ -570,6 +570,12 @@ def test_regular_mesh(lib_init): np.testing.assert_allclose(mesh.volumes, 1.0) + # bounding box + mesh.set_parameters(lower_left=ll, upper_right=ur) + bbox = mesh.bounding_box + np.testing.assert_allclose(bbox.lower_left, ll) + np.testing.assert_allclose(bbox.upper_right, ur) + meshes = openmc.lib.meshes assert isinstance(meshes, Mapping) assert len(meshes) == 1 @@ -650,6 +656,11 @@ def test_rectilinear_mesh(lib_init): np.testing.assert_allclose(mesh.volumes, 1000.0) + # bounding box + bbox = mesh.bounding_box + np.testing.assert_allclose(bbox.lower_left, (-10., 0., 10.)) + np.testing.assert_allclose(bbox.upper_right, (10., 20., 30.)) + with pytest.raises(exc.AllocationError): mesh2 = openmc.lib.RectilinearMesh(mesh.id) @@ -697,6 +708,11 @@ def test_cylindrical_mesh(lib_init): np.testing.assert_allclose(mesh.volumes[::2], 10/360 * pi * 5**2 * 10) np.testing.assert_allclose(mesh.volumes[1::2], 10/360 * pi * (10**2 - 5**2) * 10) + # bounding box + bbox = mesh.bounding_box + np.testing.assert_allclose(bbox.lower_left, (-10., -10., 10.)) + np.testing.assert_allclose(bbox.upper_right, (10., 10., 30.)) + with pytest.raises(exc.AllocationError): mesh2 = openmc.lib.CylindricalMesh(mesh.id) @@ -750,6 +766,11 @@ def test_spherical_mesh(lib_init): np.testing.assert_allclose(mesh.volumes[2::4], f * 5**3 * dtheta(10., 20.)) np.testing.assert_allclose(mesh.volumes[3::4], f * (10**3 - 5**3) * dtheta(10., 20.)) + # bounding box + bbox = mesh.bounding_box + np.testing.assert_allclose(bbox.lower_left, (-10., -10., -10.)) + np.testing.assert_allclose(bbox.upper_right, (10., 10., 10.)) + with pytest.raises(exc.AllocationError): mesh2 = openmc.lib.SphericalMesh(mesh.id)