Implement bounding_box operation for meshes (#3119)

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)

include/openmc/bounding_box.h

@@ -0,0 +1,61 @@
+#ifndef OPENMC_BOUNDING_BOX_H
+#define OPENMC_BOUNDING_BOX_H
+
+#include <algorithm> // 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

include/openmc/cell.h

@@ -12,6 +12,7 @@
 #include "pugixml.hpp"
 #include <gsl/gsl-lite.hpp>
 
+#include "openmc/bounding_box.h"
 #include "openmc/constants.h"
 #include "openmc/memory.h" // for unique_ptr
 #include "openmc/neighbor_list.h"

include/openmc/mesh.h

@@ -11,6 +11,7 @@
 #include "xtensor/xtensor.hpp"
 #include <gsl/gsl-lite.hpp>
 
+#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<MaterialVolume> 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<double, 1> lower_left_;  //!< Lower-left coordinates of mesh
+  xt::xtensor<double, 1> 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<double, 1> lower_left_;  //!< Lower-left coordinates of mesh
-  xt::xtensor<double, 1> upper_right_; //!< Upper-right coordinates of mesh
   std::array<int, 3> 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.

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<int, int> surface_map;
 extern vector<unique_ptr<Surface>> 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.
 //==============================================================================

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 {

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]

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'

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<Position, 4> 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

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);

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)