XDGMesh Object

.github/workflows/ci.yml

@@ -47,6 +47,7 @@ jobs:
         mpi: [n, y]
         omp: [n, y]
         dagmc: [n]
+        xdg: [n]
         libmesh: [n]
         event: [n]
 
@@ -64,6 +65,11 @@ jobs:
             python-version: "3.12"
             mpi: y
             omp: y
+          - xdg: y
+            python-version: "3.12"
+            libmesh: y
+            mpi: n
+            omp: y
           - libmesh: y
             python-version: "3.12"
             mpi: y
@@ -78,13 +84,15 @@ jobs:
             mpi: n
     name: "Python ${{ matrix.python-version }} (omp=${{ matrix.omp }},
       mpi=${{ matrix.mpi }}, dagmc=${{ matrix.dagmc }},
-      libmesh=${{ matrix.libmesh }}, event=${{ matrix.event }}"
+      xdg=${{ matrix.xdg }}, libmesh=${{ matrix.libmesh }},
+      event=${{ matrix.event }}"
 
     env:
       MPI: ${{ matrix.mpi }}
       PHDF5: ${{ matrix.mpi }}
       OMP: ${{ matrix.omp }}
       DAGMC: ${{ matrix.dagmc }}
+      XDG: ${{ matrix.xdg }}
       EVENT: ${{ matrix.event }}
       LIBMESH: ${{ matrix.libmesh }}
       NPY_DISABLE_CPU_FEATURES: "AVX512F AVX512_SKX"
@@ -146,6 +154,12 @@ jobs:
           sudo update-alternatives --set mpirun /usr/bin/mpirun.mpich
           sudo update-alternatives --set mpi-x86_64-linux-gnu /usr/include/x86_64-linux-gnu/mpich
 
+      - name: Optional apt dependencies for XDG
+        shell: bash
+        if: ${{ matrix.xdg == 'y' }}
+        run: |
+          sudo apt install -y libembree-dev
+
       - name: install
         shell: bash
         run: |

CMakeLists.txt

@@ -39,6 +39,7 @@ option(OPENMC_BUILD_TESTS     "Build tests"
 option(OPENMC_ENABLE_PROFILE  "Compile with profiling flags"                         OFF)
 option(OPENMC_ENABLE_COVERAGE "Compile with coverage analysis flags"                 OFF)
 option(OPENMC_USE_DAGMC       "Enable support for DAGMC (CAD) geometry"              OFF)
+option(OPENMC_USE_XDG         "Enable support for XDG discretized CAD geometry"      OFF)
 option(OPENMC_USE_LIBMESH     "Enable support for libMesh unstructured mesh tallies" OFF)
 option(OPENMC_USE_MPI         "Enable MPI"                                           OFF)
 option(OPENMC_USE_UWUW        "Enable UWUW"                                          OFF)
@@ -156,6 +157,15 @@ if(OPENMC_USE_DAGMC)
   endif()
 endif()
 
+#===============================================================================
+# XDG Mesh and Geometry Support
+#===============================================================================
+
+if(OPENMC_USE_XDG)
+  find_package(XDG REQUIRED PATH_SUFFIXES lib/cmake)
+  message(STATUS "Found XDG: ${XDG_DIR} (version ${XDG_VERSION})")
+endif()
+
 #===============================================================================
 # libMesh Unstructured Mesh Support
 #===============================================================================
@@ -428,6 +438,7 @@ list(APPEND libopenmc_SOURCES
   src/string_utils.cpp
   src/summary.cpp
   src/surface.cpp
+  src/xdg.cpp
   src/tallies/derivative.cpp
   src/tallies/filter.cpp
   src/tallies/filter_azimuthal.cpp
@@ -544,6 +555,11 @@ elseif(OPENMC_USE_UWUW)
   message(FATAL_ERROR "DAGMC must be enabled when UWUW is enabled.")
 endif()
 
+if(OPENMC_USE_XDG)
+  target_compile_definitions(libopenmc PRIVATE OPENMC_XDG_ENABLED)
+  target_link_libraries(libopenmc xdg::xdg)
+endif()
+
 if(OPENMC_USE_LIBMESH)
   target_compile_definitions(libopenmc PRIVATE OPENMC_LIBMESH_ENABLED)
   target_link_libraries(libopenmc PkgConfig::LIBMESH)

cmake/OpenMCConfig.cmake.in

@@ -1,12 +1,11 @@
 get_filename_component(OpenMC_CMAKE_DIR "${CMAKE_CURRENT_LIST_FILE}" DIRECTORY)
 
-# Compute the install prefix from this file's location
-get_filename_component(_OPENMC_PREFIX "${OpenMC_CMAKE_DIR}/../../.." ABSOLUTE)
+include(CMakeFindDependencyMacro)
 
-find_package(fmt CONFIG REQUIRED HINTS ${_OPENMC_PREFIX})
-find_package(pugixml CONFIG REQUIRED HINTS ${_OPENMC_PREFIX})
-if(@OPENMC_USE_DAGMC@)
-  find_package(DAGMC REQUIRED HINTS @DAGMC_DIR@)
+if(@OPENMC_USE_XDG@)
+  if(NOT TARGET xdg::xdg)
+    find_dependency(XDG CONFIG REQUIRED HINTS "@XDG_DIR@")
+  endif()
 endif()
 
 if(@OPENMC_USE_LIBMESH@)
@@ -16,18 +15,40 @@ if(@OPENMC_USE_LIBMESH@)
   pkg_check_modules(LIBMESH REQUIRED @LIBMESH_PC_FILE@>=1.7.0 IMPORTED_TARGET)
 endif()
 
-find_package(PNG)
+if(@OPENMC_USE_DAGMC@)
+  if(NOT TARGET dagmc-shared)
+    find_dependency(DAGMC REQUIRED HINTS "@DAGMC_DIR@")
+  endif()
+endif()
 
-if(NOT TARGET OpenMC::libopenmc)
-  include("${OpenMC_CMAKE_DIR}/OpenMCTargets.cmake")
+if(NOT TARGET fmt::fmt)
+  find_dependency(fmt CONFIG REQUIRED)
+endif()
+
+if(NOT TARGET pugixml::pugixml AND NOT TARGET pugixml)
+  find_dependency(pugixml CONFIG REQUIRED)
+endif()
+
+if(@PNG_FOUND@)
+  if(NOT TARGET PNG::PNG)
+    find_dependency(PNG REQUIRED)
+  endif()
 endif()
 
 if(@OPENMC_USE_MPI@)
-  find_package(MPI REQUIRED)
+  if(NOT TARGET MPI::MPI_CXX)
+    find_dependency(MPI REQUIRED)
+  endif()
 endif()
 
 if(@OPENMC_USE_OPENMP@)
-  find_package(OpenMP REQUIRED)
+  if(NOT TARGET OpenMP::OpenMP_CXX)
+    find_dependency(OpenMP REQUIRED)
+  endif()
+endif()
+
+if(NOT TARGET OpenMC::libopenmc)
+  include("${OpenMC_CMAKE_DIR}/OpenMCTargets.cmake")
 endif()
 
 if(@OPENMC_USE_UWUW@ AND NOT ${DAGMC_BUILD_UWUW})

include/openmc/mesh.h

@@ -47,6 +47,12 @@ namespace openmc {
 
 enum class ElementType { UNSUPPORTED = -1, LINEAR_TET, LINEAR_HEX };
 
+struct NextMeshCell {
+  double distance {INFTY};
+  int face_idx {-1};
+  std::array<int, 3> next_ijk;
+};
+
 //==============================================================================
 // Global variables
 //==============================================================================
@@ -771,6 +777,9 @@ class UnstructuredMesh : public Mesh {
   //! Get the library used for this unstructured mesh
   virtual std::string library() const = 0;
 
+  //! Get the mesh filename
+  virtual const std::string& filename() const { return filename_; }
+
   // Data members
   bool output_ {
     true}; //!< Write tallies onto the unstructured mesh at the end of a run
@@ -797,7 +806,36 @@ class UnstructuredMesh : public Mesh {
   //! \param[in] coords Coordinates of the tetrahedron
   //! \param[in] seed Random number generation seed
   //! \return Sampled position within the tetrahedron
-  Position sample_tet(std::array<Position, 4> coords, uint64_t* seed) const;
+  template<typename V>
+  Position sample_tet(span<V> coords, uint64_t* seed) const
+  {
+    // Uniform distribution
+    double s = prn(seed);
+    double t = prn(seed);
+    double u = prn(seed);
+
+    // From PyNE implementation of moab tet sampling C. Rocchini & P. Cignoni
+    // (2000) Generating Random Points in a Tetrahedron, Journal of Graphics
+    // Tools, 5:4, 9-12, DOI: 10.1080/10867651.2000.10487528
+    if (s + t > 1) {
+      s = 1.0 - s;
+      t = 1.0 - t;
+    }
+    if (s + t + u > 1) {
+      if (t + u > 1) {
+        double old_t = t;
+        t = 1.0 - u;
+        u = 1.0 - s - old_t;
+      } else if (t + u <= 1) {
+        double old_s = s;
+        s = 1.0 - t - u;
+        u = old_s + t + u - 1;
+      }
+    }
+    V result = s * (coords[1] - coords[0]) + t * (coords[2] - coords[0]) +
+               u * (coords[3] - coords[0]) + coords[0];
+    return {result[0], result[1], result[2]};
+  }
 
   // Data members
   double length_multiplier_ {

include/openmc/xdg.h

@@ -0,0 +1,113 @@
+#ifndef OPENMC_XDG_H
+#define OPENMC_XDG_H
+
+namespace openmc {
+extern "C" const bool XDG_ENABLED;
+}
+
+// always include the XML interface header
+#include "openmc/xml_interface.h"
+
+#ifdef OPENMC_XDG_ENABLED
+
+#include "xdg/xdg.h"
+
+#include "openmc/mesh.h"
+#include "openmc/position.h"
+namespace openmc {
+
+class XDGMesh : public UnstructuredMesh {
+
+public:
+  // Constructors
+  XDGMesh() = default;
+  XDGMesh(pugi::xml_node node);
+  XDGMesh(hid_t group);
+  XDGMesh(const std::string& filename, double length_multiplier = 1.0);
+  XDGMesh(std::shared_ptr<xdg::XDG> external_xdg);
+
+  static const std::string mesh_lib_type;
+
+  const std::shared_ptr<xdg::XDG>& xdg_instance() const { return xdg_; }
+
+  // Overridden Methods
+
+  //! Perform any preparation needed to support use in mesh filters
+  void prepare_for_point_location() override;
+
+  Position sample_element(int32_t bin, uint64_t* seed) const override;
+
+  void bins_crossed(Position r0, Position r1, const Direction& u,
+    vector<int>& bins, vector<double>& lengths) const override;
+
+  int get_bin(Position r) const override;
+
+  bool bin_is_valid(int bin) const { return bin >= 0 && bin < n_bins(); }
+
+  xdg::MeshID bin_to_mesh_id(int bin) const;
+
+  int mesh_id_to_bin(xdg::MeshID id) const;
+
+  int n_bins() const override;
+
+  int n_surface_bins() const override;
+
+  std::pair<vector<double>, vector<double>> plot(
+    Position plot_ll, Position plot_ur) const override;
+
+  std::string library() const override;
+
+  std::string mesh_library() const;
+
+  //! Add a score to the mesh instance
+  void add_score(const std::string& score) override {};
+
+  //! Remove all scores from the mesh instance
+  void remove_scores() override {};
+
+  //! Set data for a score
+  void set_score_data(const std::string& score, const vector<double>& values,
+    const vector<double>& std_dev) override {};
+
+  //! Write the mesh with any current tally data
+  void write(const std::string& base_filename) const override;
+
+  Position centroid(int bin) const override;
+
+  int n_vertices() const override;
+
+  Position vertex(int id) const override;
+
+  std::vector<int> connectivity(int id) const override;
+
+  //! Get the volume of a mesh bin
+  //
+  //! \param[in] bin Bin to return the volume for
+  //! \return Volume of the bin
+  double volume(int bin) const override;
+
+  //! Get the distance to the nearest boundary for a given position and
+  //! direction \param[in] g GeometryState object containing position and
+  //! direction \return NextMeshCell struct containing distance, face index, and
+  //! next indices
+  NextMeshCell distance_to_bin_boundary(GeometryState& g) const;
+
+  //! Get the distance to the nearest boundary for a given position and
+  //! direction \param[in] r Position to check \param[in] u Direction to check
+  //! \return Distance to the nearest boundary
+  NextMeshCell distance_to_bin_boundary(
+    int bin, const Position& r, const Direction& u) const;
+
+private:
+  void initialize() override;
+
+  std::shared_ptr<xdg::XDG> xdg_; //!< XDG instance
+  xdg::MeshLibrary mesh_library_ {
+    xdg::MeshLibrary::LIBMESH}; //!< Mesh library type
+};
+
+} // namespace openmc
+
+#endif // OPENMC_XDG_ENABLED
+
+#endif // OPENMC_XDG_H

openmc/__init__.py

@@ -34,6 +34,7 @@
 from openmc.search import *
 from openmc.polynomial import *
 from openmc.tracks import *
+from openmc.xdg import *
 from .config import *
 
 # Import a few names from the model module

openmc/filter.py

@@ -925,6 +925,8 @@ def mesh(self, mesh):
                 self.bins = list(range(len(mesh.volumes)))
             else:
                 self.bins = []
+        elif isinstance(mesh, openmc.XDGMesh):
+            self.bins = []
         else:
             self.bins = list(mesh.indices)
 

openmc/lib/__init__.py

@@ -40,6 +40,9 @@
 def _dagmc_enabled():
     return c_bool.in_dll(_dll, "DAGMC_ENABLED").value
 
+def _xdg_enabled():
+    return c_bool.in_dll(_dll, "XDG_ENABLED").value
+
 def _coord_levels():
     return c_int.in_dll(_dll, "n_coord_levels").value
 

openmc/mesh.py

@@ -337,6 +337,8 @@ def from_xml_element(cls, elem: ET.Element):
             an openmc mesh object
 
         """
+        from .xdg import XDGMesh
+
         mesh_type = get_text(elem, 'type')
 
         if mesh_type == 'regular' or mesh_type is None:
@@ -347,6 +349,8 @@ def from_xml_element(cls, elem: ET.Element):
             mesh = CylindricalMesh.from_xml_element(elem)
         elif mesh_type == 'spherical':
             mesh = SphericalMesh.from_xml_element(elem)
+        elif mesh_type == 'xdg':
+            mesh = XDGMesh.from_xml_element(elem)
         elif mesh_type == 'unstructured':
             mesh = UnstructuredMesh.from_xml_element(elem)
         else:
@@ -2653,7 +2657,7 @@ def library(self):
 
     @library.setter
     def library(self, lib: str):
-        cv.check_value('Unstructured mesh library', lib, ('moab', 'libmesh'))
+        cv.check_value('Unstructured mesh library', lib, ('xdg','moab', 'libmesh'))
         self._library = lib
 
     @property