From 64bcb0d98a9e7611de9a670f0585105995382e97 Mon Sep 17 00:00:00 2001 From: magnoxemo Date: Tue, 2 Jun 2026 12:20:22 -0500 Subject: [PATCH 01/26] adding omega_h install path to the build system --- CMakeLists.txt | 11 +++++++++++ include/xdg/omega_h/mesh_manager.h | 28 ++++++++++++++++++++++++++++ 2 files changed, 39 insertions(+) create mode 100644 include/xdg/omega_h/mesh_manager.h diff --git a/CMakeLists.txt b/CMakeLists.txt index fb762521..b2574df9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -4,6 +4,7 @@ project(openmc C CXX) option(XDG_ENABLE_MOAB "Enable support for the MOAB mesh library" ON) option(XDG_ENABLE_MFEM "Enable support for the MFEM mesh library" OFF) option(XDG_ENABLE_LIBMESH "Enable support for the libMesh mesh library" OFF) +option(XDG_ENABLE_OMEGA_H "Enable support for the Omega_h mesh library" ON) option(XDG_LINK_MPI "Link with MPI (for dependency compatibility)" OFF) option(XDG_ENABLE_EMBREE "Enable support for the Embree ray tracing library" ON) option(XDG_ENABLE_GPRT "Enable support for the GPRT ray tracing library" OFF) @@ -43,6 +44,16 @@ if (NOT MOAB_USE_HDF5) endif() endif() + +#=============================================================================== +# Omega_h backend +#=============================================================================== +if (XDG_ENABLE_OMEGA_H) + message(STATUS "OMEGA_H_INSTALL_PATH=$ENV{OMEGA_H_INSTALL_PATH}") + find_package(Omega_h REQUIRED HINTS "$ENV{OMEGA_H_INSTALL_PATH}/lib/cmake/Omega_h") +endif() + + if (XDG_ENABLE_EMBREE) # find Embree for CPU ray tracing diff --git a/include/xdg/omega_h/mesh_manager.h b/include/xdg/omega_h/mesh_manager.h new file mode 100644 index 00000000..86737d41 --- /dev/null +++ b/include/xdg/omega_h/mesh_manager.h @@ -0,0 +1,28 @@ +#ifndef XDG_OMEGA_H_MESH_MANAGER_H +#define XDG_OMEGA_H_MESH_MANAGER_H + + +#include "mesh_manager_interface.h" + +namespace xdg{ + +class OmegaHMeshManager:public MeshManager{ + +public: + OmegaHMeshManager(); + ~OmegaHMeshManager()=default; + + void load_file(const std::string& file_path) override; + + void init() override; + + int num_volumes() const override; + int num_surfaces() const override; + + int num_ents_of_dimension(int dim) const; + +}; + + +} +#endif //XDG_MESH_MANAGER_H From b2665149f9a32e8b2ac0719cf87e24b8ebbd83f8 Mon Sep 17 00:00:00 2001 From: magnoxemo Date: Sat, 6 Jun 2026 20:30:41 -0500 Subject: [PATCH 02/26] add omega_h to the build system --- CMakeLists.txt | 21 +++++++++++-- Dockerfile | 0 include/xdg/omega_h/mesh_manager.h | 49 +++++++++++++++++++++++++++++- src/omega_h/mesh_manager.cpp | 2 ++ 4 files changed, 68 insertions(+), 4 deletions(-) create mode 100644 Dockerfile create mode 100644 src/omega_h/mesh_manager.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index b2574df9..f27dca41 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -44,7 +44,6 @@ if (NOT MOAB_USE_HDF5) endif() endif() - #=============================================================================== # Omega_h backend #=============================================================================== @@ -53,7 +52,6 @@ if (XDG_ENABLE_OMEGA_H) find_package(Omega_h REQUIRED HINTS "$ENV{OMEGA_H_INSTALL_PATH}/lib/cmake/Omega_h") endif() - if (XDG_ENABLE_EMBREE) # find Embree for CPU ray tracing @@ -178,7 +176,7 @@ endif() # Ensure at least one mesh backend is enabled -if (NOT XDG_ENABLE_MOAB AND NOT XDG_ENABLE_LIBMESH) +if (NOT XDG_ENABLE_MOAB AND NOT XDG_ENABLE_LIBMESH AND NOT XDG_ENABLE_OMEGA_H) message(FATAL_ERROR "No mesh backend enabled. Enable at least one of:\n" " -DXDG_ENABLE_MOAB=ON\n" @@ -247,6 +245,12 @@ src/moab/metadata.cpp ) endif() +if (XDG_ENABLE_OMEGA_H) + list(APPEND xdg_sources + src/omega_h/mesh_manager.cpp + ) +endif () + #=============================================================================== # RPATH information (from OpenMC) #=============================================================================== @@ -325,6 +329,12 @@ if (XDG_ENABLE_EMBREE) target_compile_definitions(xdg PUBLIC XDG_ENABLE_EMBREE) endif() + +if (XDG_ENABLE_OMEGA_H) + target_compile_definitions(xdg PUBLIC XDG_ENABLE_OMEGA_H) +endif() + + if (XDG_ENABLE_GPRT) target_compile_definitions(xdg PUBLIC XDG_ENABLE_GPRT) target_link_options(xdg PRIVATE -Wl,--unresolved-symbols=ignore-in-shared-libs) @@ -365,6 +375,11 @@ if(XDG_ENABLE_LIBMESH) target_link_libraries(xdg PRIVATE PkgConfig::LIBMESH) endif() +if(XDG_ENABLE_OMEGA_H) + target_link_libraries(xdg PRIVATE Omega_h::omega_h) +endif() + + # this provides the ability to link MPI when needed for compatibility with other # libraries (e.g., libMesh) if (XDG_LINK_MPI) diff --git a/Dockerfile b/Dockerfile new file mode 100644 index 00000000..e69de29b diff --git a/include/xdg/omega_h/mesh_manager.h b/include/xdg/omega_h/mesh_manager.h index 86737d41..2ff4aadb 100644 --- a/include/xdg/omega_h/mesh_manager.h +++ b/include/xdg/omega_h/mesh_manager.h @@ -2,7 +2,8 @@ #define XDG_OMEGA_H_MESH_MANAGER_H -#include "mesh_manager_interface.h" +#include "xdg/mesh_manager_interface.h" + namespace xdg{ @@ -21,6 +22,52 @@ class OmegaHMeshManager:public MeshManager{ int num_ents_of_dimension(int dim) const; + + MeshID create_volume() override; + + void add_surface_to_volume(MeshID volume, MeshID surface, Sense sense, bool overwrite=false) override; + + // Mesh + int num_volume_elements(MeshID volume) const override; + + int num_volume_elements() const override; + + int num_volume_faces(MeshID volume) const override; + + int num_surface_faces(MeshID surface) const override; + + int num_vertices() const override; + + std::vector get_volume_elements(MeshID volume) const override; + + std::vector get_surface_faces(MeshID surface) const override; + + std::vector element_connectivity(MeshID element) const override; + + std::vector face_connectivity(MeshID face) const override; + + Vertex vertex_coordinates(MeshID vertex) const override; + + std::vector element_vertices(MeshID element) const override; + + std::array face_vertices(MeshID face) const override; + + SurfaceElementType get_surface_element_type(MeshID surface) const override; + + MeshID adjacent_element(MeshID element, int face) const override; + + double element_volume(MeshID element) const override; + + // Topology + std::pair surface_senses(MeshID surface) const override; + + std::vector get_volume_surfaces(MeshID volume) const override; + + Sense surface_sense(MeshID surface, MeshID volume) const override; + + // Metadata + void parse_metadata() override; + }; diff --git a/src/omega_h/mesh_manager.cpp b/src/omega_h/mesh_manager.cpp new file mode 100644 index 00000000..85dda9ea --- /dev/null +++ b/src/omega_h/mesh_manager.cpp @@ -0,0 +1,2 @@ +#include "xdg/omega_h/mesh_manager.h" + From 8067920a75fbb5a4346f35563a178bcceccd4d45 Mon Sep 17 00:00:00 2001 From: magnoxemo Date: Sat, 6 Jun 2026 21:11:45 -0500 Subject: [PATCH 03/26] override the pure virtual functions --- include/xdg/omega_h/mesh_manager.h | 101 ++++++++++++++--------------- 1 file changed, 50 insertions(+), 51 deletions(-) diff --git a/include/xdg/omega_h/mesh_manager.h b/include/xdg/omega_h/mesh_manager.h index 2ff4aadb..c5d3f810 100644 --- a/include/xdg/omega_h/mesh_manager.h +++ b/include/xdg/omega_h/mesh_manager.h @@ -1,75 +1,74 @@ #ifndef XDG_OMEGA_H_MESH_MANAGER_H #define XDG_OMEGA_H_MESH_MANAGER_H +#include +#include +#include #include "xdg/mesh_manager_interface.h" +#include "Omega_h_library.hpp" +#include "Omega_h_map.hpp" +#include "Omega_h_mesh.hpp" +#include "Omega_h_profile.hpp" +#include "Omega_h_quality.hpp" +#include "Omega_h_timer.hpp" -namespace xdg{ - -class OmegaHMeshManager:public MeshManager{ +namespace xdg { +class OmegaHMeshManager : public MeshManager { public: - OmegaHMeshManager(); - ~OmegaHMeshManager()=default; - - void load_file(const std::string& file_path) override; - - void init() override; - - int num_volumes() const override; - int num_surfaces() const override; - - int num_ents_of_dimension(int dim) const; - - - MeshID create_volume() override; - - void add_surface_to_volume(MeshID volume, MeshID surface, Sense sense, bool overwrite=false) override; - - // Mesh - int num_volume_elements(MeshID volume) const override; - - int num_volume_elements() const override; - - int num_volume_faces(MeshID volume) const override; - - int num_surface_faces(MeshID surface) const override; - - int num_vertices() const override; - - std::vector get_volume_elements(MeshID volume) const override; - - std::vector get_surface_faces(MeshID surface) const override; - - std::vector element_connectivity(MeshID element) const override; + OmegaHMeshManager(); + ~OmegaHMeshManager() = default; - std::vector face_connectivity(MeshID face) const override; + void load_file(const std::string &file_path) override; - Vertex vertex_coordinates(MeshID vertex) const override; + void init() override; - std::vector element_vertices(MeshID element) const override; + int num_volumes() const override; + int num_surfaces() const override; + int num_ents_of_dimension(int dim) const; - std::array face_vertices(MeshID face) const override; + MeshID create_volume() override; + void add_surface_to_volume(MeshID volume, MeshID surface, Sense sense, + bool overwrite = false) override; - SurfaceElementType get_surface_element_type(MeshID surface) const override; + int num_volume_elements(MeshID volume) const override; + int num_volume_elements() const override; + int num_volume_faces(MeshID volume) const override; + int num_surface_faces(MeshID surface) const override; + int num_vertices() const override; - MeshID adjacent_element(MeshID element, int face) const override; + std::vector get_volume_elements(MeshID volume) const override; + std::vector get_surface_faces(MeshID surface) const override; + std::vector element_connectivity(MeshID element) const override; + std::vector face_connectivity(MeshID face) const override; + Vertex vertex_coordinates(MeshID vertex) const override; + std::vector element_vertices(MeshID element) const override; + std::array face_vertices(MeshID face) const override; + SurfaceElementType get_surface_element_type(MeshID surface) const override; + MeshID adjacent_element(MeshID element, int face) const override; + double element_volume(MeshID element) const override; - double element_volume(MeshID element) const override; + std::pair surface_senses(MeshID surface) const override; + std::vector get_volume_surfaces(MeshID volume) const override; + Sense surface_sense(MeshID surface, MeshID volume) const override; - // Topology - std::pair surface_senses(MeshID surface) const override; + void parse_metadata() {}; - std::vector get_volume_surfaces(MeshID volume) const override; + Omega_h::Mesh *mesh() { return mesh_.get(); } + const Omega_h::Mesh *mesh() const { return mesh_.get(); } - Sense surface_sense(MeshID surface, MeshID volume) const override; +private: + std::unique_ptr library_; + std::unique_ptr mesh_; - // Metadata - void parse_metadata() override; + std::map> surface_senses_; + MeshID next_volume_id_{1}; + MeshID next_surface_id_{1}; }; +} // namespace xdg -} -#endif //XDG_MESH_MANAGER_H +#endif // XDG_OMEGA_H_MESH_MANAGER_H \ No newline at end of file From da32c5b2054aef9bfd2982487b1bd0d13d2b28f4 Mon Sep 17 00:00:00 2001 From: magnoxemo Date: Sat, 6 Jun 2026 23:46:22 -0500 Subject: [PATCH 04/26] update CI to include omega_h --- .github/workflows/ci.yml | 60 +++++++++++++++++++++++++++++++++++++--- CMakeLists.txt | 21 ++++++++------ 2 files changed, 69 insertions(+), 12 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 1ac4d8a2..be49a94c 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -10,16 +10,18 @@ on: jobs: main: name: | - libMesh=${{ matrix.libmesh }};MOAB=${{ matrix.moab }} + libMesh=${{ matrix.libmesh }};MOAB=${{ matrix.moab }};Omega_h=${{ matrix.omega_h }} runs-on: ubuntu-latest strategy: fail-fast: false matrix: libmesh: [true, false] moab: [true, false] + omega_h: [true, false] exclude: - moab: false libmesh: false + omega_h: false steps: - name: Checkout @@ -81,6 +83,15 @@ jobs: cd libmesh git checkout v1.7.0 + - name: Omega_h Clone + if: ${{ matrix.omega_h }} + shell: bash + run: | + cd ~ + git clone https://github.com/SCOREC/omega_h.git + cd omega_h + git checkout master + - name: OpenMP Environment Variables run: | echo "OMP_NUM_THREADS=1" >> $GITHUB_ENV @@ -94,9 +105,13 @@ jobs: if: ${{ matrix.libmesh }} run: | echo "LIBMESH_SHA"=$(cd ~/libmesh && git rev-parse HEAD) >> $GITHUB_ENV - # Enforce that we're using the debug build of libMesh echo "METHOD=dbg" >> $GITHUB_ENV + - name: Omega_h Environment Variables + if: ${{ matrix.omega_h }} + run: | + echo "OMEGA_H_SHA"=$(cd ~/omega_h && git rev-parse HEAD) >> $GITHUB_ENV + - name: MOAB Cache if: ${{ matrix.moab }} id: moab-cache @@ -117,6 +132,16 @@ jobs: path: ~/LIBMESH key: libmesh-${{ runner.os }}-${{ env.cache-name }}-${{ env.LIBMESH_SHA }} + - name: Omega_h Cache + if: ${{ matrix.omega_h }} + id: omega-h-cache + uses: actions/cache@v3 + env: + cache-name: omega-h-cache + with: + path: ~/OMEGA_H + key: omega-h-${{ runner.os }}-${{ env.cache-name }}-${{ env.OMEGA_H_SHA }} + - if: ${{ matrix.moab && steps.moab-cache.outputs.cache-hit != 'true' }} name: Build MOAB run: | @@ -141,6 +166,21 @@ jobs: make -j4 sudo make install + - if: ${{ matrix.omega_h && steps.omega-h-cache.outputs.cache-hit != 'true' }} + name: Build Omega_h + shell: bash + run: | + cd ~/omega_h + mkdir build + cd build + cmake .. \ + -DCMAKE_INSTALL_PREFIX=$HOME/OMEGA_H \ + -DCMAKE_BUILD_TYPE=Release \ + -DOmega_h_USE_MPI=OFF \ + -DBUILD_TESTING=OFF + make -j4 + make install + - name: Build shell: bash run: | @@ -157,7 +197,19 @@ jobs: fi CMAKE_PREFIX_PATH="${CMAKE_PREFIX_PATH}$HOME/MOAB" fi - cmake .. -DCMAKE_PREFIX_PATH="$CMAKE_PREFIX_PATH" -DCMAKE_INSTALL_PREFIX=$HOME/opt -DXDG_ENABLE_MOAB=${{ matrix.moab && 'ON' || 'OFF' }} -DXDG_ENABLE_LIBMESH=${{ matrix.libmesh && 'ON' || 'OFF' }} + if [ "${{ matrix.omega_h }}" = 'true' ]; then + if [ -n "$CMAKE_PREFIX_PATH" ]; then + CMAKE_PREFIX_PATH="$CMAKE_PREFIX_PATH;" + fi + CMAKE_PREFIX_PATH="${CMAKE_PREFIX_PATH}$HOME/OMEGA_H" + fi + cmake .. \ + -DCMAKE_PREFIX_PATH="$CMAKE_PREFIX_PATH" \ + -DCMAKE_INSTALL_PREFIX=$HOME/opt \ + -DXDG_ENABLE_MOAB=${{ matrix.moab && 'ON' || 'OFF' }} \ + -DXDG_ENABLE_LIBMESH=${{ matrix.libmesh && 'ON' || 'OFF' }} \ + -DXDG_ENABLE_OMEGA_H=${{ matrix.omega_h && 'ON' || 'OFF' }} \ + -DOMEGA_H_INSTALL_PATH=$HOME/OMEGA_H make -j4 install - name: Test @@ -173,4 +225,4 @@ jobs: - name: Setup tmate session if: ${{ failure() }} - uses: mxschmitt/action-tmate@v3 + uses: mxschmitt/action-tmate@v3 \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt index f27dca41..2e2dc60a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,5 +1,5 @@ cmake_minimum_required(VERSION 3.10 FATAL_ERROR) -project(openmc C CXX) +project(xdg C CXX) option(XDG_ENABLE_MOAB "Enable support for the MOAB mesh library" ON) option(XDG_ENABLE_MFEM "Enable support for the MFEM mesh library" OFF) @@ -15,7 +15,7 @@ option(XDG_BUILD_TOOLS "Enable tools and miniapps" ON set(XDG_VERSION_MAJOR 0) set(XDG_VERSION_MINOR 14) set(XDG_VERSION_RELEASE 1) -set(XDG_VERSION ${OPENMC_VERSION_MAJOR}.${OPENMC_VERSION_MINOR}.${OPENMC_VERSION_RELEASE}) +set(XDG_VERSION ${XDG_VERSION_MAJOR}.${XDG_VERSION_MINOR}.${XDG_VERSION_RELEASE}) if(NOT CMAKE_BUILD_TYPE) set(CMAKE_BUILD_TYPE Release CACHE STRING "Choose build type" FORCE) @@ -48,10 +48,16 @@ endif() # Omega_h backend #=============================================================================== if (XDG_ENABLE_OMEGA_H) - message(STATUS "OMEGA_H_INSTALL_PATH=$ENV{OMEGA_H_INSTALL_PATH}") - find_package(Omega_h REQUIRED HINTS "$ENV{OMEGA_H_INSTALL_PATH}/lib/cmake/Omega_h") + if (NOT OMEGA_H_INSTALL_PATH) + set(OMEGA_H_INSTALL_PATH "$ENV{OMEGA_H_INSTALL_PATH}") + endif() + message(STATUS "OMEGA_H_INSTALL_PATH=${OMEGA_H_INSTALL_PATH}") + find_package(Omega_h REQUIRED HINTS "${OMEGA_H_INSTALL_PATH}/lib/cmake/Omega_h") endif() +#=============================================================================== +# Embree +#=============================================================================== if (XDG_ENABLE_EMBREE) # find Embree for CPU ray tracing @@ -131,7 +137,7 @@ if(GIT_FOUND AND EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/.git") RESULT_VARIABLE GIT_SUBMOD_RESULT) if(NOT GIT_SUBMOD_RESULT EQUAL 0) message(FATAL_ERROR - "git submodule update --init failed with erro code ${GIT_SUBMOD_RESULT}, " + "git submodule update --init failed with error code ${GIT_SUBMOD_RESULT}, " "please checkout submodule at ${VENDOR_PATH} manually and try again.") endif() endforeach() @@ -157,7 +163,6 @@ if (XDG_BUILD_TESTS) message(STATUS "Using Catch2 ${Catch2_VERSION} at ${Catch2_DIR}") endif() - # fmt if (NOT fmt_FOUND) set(FMT_INSTALL ON CACHE BOOL "Generate the fmt install target") @@ -174,13 +179,13 @@ if(XDG_BUILD_TOOLS) add_subdirectory(${CMAKE_SOURCE_DIR}/vendor/indicators) endif() - # Ensure at least one mesh backend is enabled if (NOT XDG_ENABLE_MOAB AND NOT XDG_ENABLE_LIBMESH AND NOT XDG_ENABLE_OMEGA_H) message(FATAL_ERROR "No mesh backend enabled. Enable at least one of:\n" " -DXDG_ENABLE_MOAB=ON\n" - " -DXDG_ENABLE_LIBMESH=ON") + " -DXDG_ENABLE_LIBMESH=ON" + " -DXDG_ENABLE_OMEGA_H=ON") endif() # Ensure at least one ray tracing backend is enabled From b980783942f99cc437984cf8ceeb74e16d1d3896 Mon Sep 17 00:00:00 2001 From: magnoxemo Date: Mon, 8 Jun 2026 10:39:38 -0500 Subject: [PATCH 05/26] update openmc testing ci to include omega_h as well --- .github/workflows/ci.yml | 1 + .github/workflows/openmc-test.yml | 39 ++++++++++++++++++++++++++++++- CMakeLists.txt | 2 +- 3 files changed, 40 insertions(+), 2 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index be49a94c..3c001842 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -105,6 +105,7 @@ jobs: if: ${{ matrix.libmesh }} run: | echo "LIBMESH_SHA"=$(cd ~/libmesh && git rev-parse HEAD) >> $GITHUB_ENV + # Enforce that we're using the debug build of libMesh echo "METHOD=dbg" >> $GITHUB_ENV - name: Omega_h Environment Variables diff --git a/.github/workflows/openmc-test.yml b/.github/workflows/openmc-test.yml index f232e34e..359ddcb8 100644 --- a/.github/workflows/openmc-test.yml +++ b/.github/workflows/openmc-test.yml @@ -50,10 +50,19 @@ jobs: cd libmesh git checkout v1.7.0 + - name: Omega_h Clone + shell: bash + run: | + cd ~ + git clone https://github.com/SCOREC/omega_h.git + cd omega_h + git checkout master + - name: Cache Variables run: | echo "MOAB_SHA"=$(cd ~/moab && git rev-parse HEAD) >> $GITHUB_ENV echo "LIBMESH_SHA"=$(cd ~/libmesh && git rev-parse HEAD) >> $GITHUB_ENV + echo "OMEGA_H_SHA"=$(cd ~/omega_h && git rev-parse HEAD) >> $GITHUB_ENV - name: MOAB Cache id: moab-cache @@ -73,6 +82,15 @@ jobs: path: ~/LIBMESH key: libmesh-${{ runner.os }}-${{ env.cache-name }}-${{ env.LIBMESH_SHA }} + - name: Omega_h Cache + id: omega-h-cache + uses: actions/cache@v3 + env: + cache-name: omega-h-cache + with: + path: ~/OMEGA_H + key: omega-h-${{ runner.os }}-${{ env.cache-name }}-${{ env.OMEGA_H_SHA }} + - if: ${{ steps.moab-cache.outputs.cache-hit != 'true' }} name: Build MOAB shell: bash @@ -95,13 +113,32 @@ jobs: ../configure --prefix=$HOME/LIBMESH --enable-exodus --disable-netcdf4 --disable-eigen --disable-lapack --disable-mpi --disable-metaphysicl make -j4 make install + - if: ${{ steps.omega-h-cache.outputs.cache-hit != 'true' }} + name: Build Omega_h + shell: bash + run: | + cd ~/omega_h + mkdir build + cd build + cmake .. \ + -DCMAKE_INSTALL_PREFIX=$HOME/OMEGA_H \ + -DCMAKE_BUILD_TYPE=Release \ + -DOmega_h_USE_MPI=OFF \ + -DBUILD_TESTING=OFF + make -j4 + make install - name: Build XDG shell: bash run: | mkdir build cd build - cmake .. -DCMAKE_INSTALL_PREFIX=$HOME/opt -DCMAKE_PREFIX_PATH="$HOME/MOAB;$HOME/LIBMESH" -DXDG_ENABLE_MOAB=ON -DXDG_ENABLE_LIBMESH=ON + cmake .. \ + -DCMAKE_INSTALL_PREFIX=$HOME/opt \ + -DCMAKE_PREFIX_PATH="$HOME/MOAB;$HOME/LIBMESH;$HOME/OMEGA_H" \ + -DXDG_ENABLE_MOAB=ON \ + -DXDG_ENABLE_LIBMESH=ON \ + -DXDG_ENABLE_OMEGA_H=ON \ make -j4 all install - name: Build OpenMC diff --git a/CMakeLists.txt b/CMakeLists.txt index 2e2dc60a..560e7ef2 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -184,7 +184,7 @@ if (NOT XDG_ENABLE_MOAB AND NOT XDG_ENABLE_LIBMESH AND NOT XDG_ENABLE_OMEGA_H) message(FATAL_ERROR "No mesh backend enabled. Enable at least one of:\n" " -DXDG_ENABLE_MOAB=ON\n" - " -DXDG_ENABLE_LIBMESH=ON" + " -DXDG_ENABLE_LIBMESH=ON\n" " -DXDG_ENABLE_OMEGA_H=ON") endif() From e954124199a3c99e5583e257ecc4ca8891b3eb07 Mon Sep 17 00:00:00 2001 From: magnoxemo Date: Mon, 8 Jun 2026 11:15:01 -0500 Subject: [PATCH 06/26] concentrating on changes that are necessary only for omega_h support --- .github/workflows/openmc-test.yml | 2 +- CMakeLists.txt | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/.github/workflows/openmc-test.yml b/.github/workflows/openmc-test.yml index 359ddcb8..aafea7d1 100644 --- a/.github/workflows/openmc-test.yml +++ b/.github/workflows/openmc-test.yml @@ -138,7 +138,7 @@ jobs: -DCMAKE_PREFIX_PATH="$HOME/MOAB;$HOME/LIBMESH;$HOME/OMEGA_H" \ -DXDG_ENABLE_MOAB=ON \ -DXDG_ENABLE_LIBMESH=ON \ - -DXDG_ENABLE_OMEGA_H=ON \ + -DXDG_ENABLE_OMEGA_H=ON make -j4 all install - name: Build OpenMC diff --git a/CMakeLists.txt b/CMakeLists.txt index 560e7ef2..1af551d0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,10 +1,10 @@ cmake_minimum_required(VERSION 3.10 FATAL_ERROR) -project(xdg C CXX) +project(openmc C CXX) option(XDG_ENABLE_MOAB "Enable support for the MOAB mesh library" ON) option(XDG_ENABLE_MFEM "Enable support for the MFEM mesh library" OFF) option(XDG_ENABLE_LIBMESH "Enable support for the libMesh mesh library" OFF) -option(XDG_ENABLE_OMEGA_H "Enable support for the Omega_h mesh library" ON) +option(XDG_ENABLE_OMEGA_H "Enable support for the Omega_h mesh library" OFF) option(XDG_LINK_MPI "Link with MPI (for dependency compatibility)" OFF) option(XDG_ENABLE_EMBREE "Enable support for the Embree ray tracing library" ON) option(XDG_ENABLE_GPRT "Enable support for the GPRT ray tracing library" OFF) @@ -15,7 +15,7 @@ option(XDG_BUILD_TOOLS "Enable tools and miniapps" ON set(XDG_VERSION_MAJOR 0) set(XDG_VERSION_MINOR 14) set(XDG_VERSION_RELEASE 1) -set(XDG_VERSION ${XDG_VERSION_MAJOR}.${XDG_VERSION_MINOR}.${XDG_VERSION_RELEASE}) +set(XDG_VERSION ${OPENMC_VERSION_MAJOR}.${OPENMC_VERSION_MINOR}.${OPENMC_VERSION_RELEASE}) if(NOT CMAKE_BUILD_TYPE) set(CMAKE_BUILD_TYPE Release CACHE STRING "Choose build type" FORCE) From a8c949105e93902a7b1d37bef93e80dfd84d6706 Mon Sep 17 00:00:00 2001 From: magnoxemo Date: Mon, 8 Jun 2026 11:18:22 -0500 Subject: [PATCH 07/26] keeping my cicd docker file locally --- Dockerfile | 0 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 Dockerfile diff --git a/Dockerfile b/Dockerfile deleted file mode 100644 index e69de29b..00000000 From afc7b2c03b3f8443e2ddc57b9f44208d8568013f Mon Sep 17 00:00:00 2001 From: magnoxemo Date: Mon, 15 Jun 2026 10:08:10 -0500 Subject: [PATCH 08/26] still need to implement parse metadata --- include/xdg/omega_h/mesh_manager.h | 2 +- src/omega_h/mesh_manager.cpp | 346 +++++++++++++++++++++++++++++ 2 files changed, 347 insertions(+), 1 deletion(-) diff --git a/include/xdg/omega_h/mesh_manager.h b/include/xdg/omega_h/mesh_manager.h index c5d3f810..1349a80a 100644 --- a/include/xdg/omega_h/mesh_manager.h +++ b/include/xdg/omega_h/mesh_manager.h @@ -54,7 +54,7 @@ class OmegaHMeshManager : public MeshManager { std::vector get_volume_surfaces(MeshID volume) const override; Sense surface_sense(MeshID surface, MeshID volume) const override; - void parse_metadata() {}; + void parse_metadata(); Omega_h::Mesh *mesh() { return mesh_.get(); } const Omega_h::Mesh *mesh() const { return mesh_.get(); } diff --git a/src/omega_h/mesh_manager.cpp b/src/omega_h/mesh_manager.cpp index 85dda9ea..e6cdfee3 100644 --- a/src/omega_h/mesh_manager.cpp +++ b/src/omega_h/mesh_manager.cpp @@ -1,2 +1,348 @@ #include "xdg/omega_h/mesh_manager.h" +#include +#include +#include +#include +#include +#include +#include + +#include "xdg/error.h" +#include "xdg/geometry/measure.h" +#include "xdg/moab/tag_conventions.h" +#include "xdg/util/str_utils.h" + +#include "Omega_h_file.hpp" +#include "Omega_h_for.hpp" +#include "Omega_h_mesh.hpp" + +namespace xdg { + +// Omega_h dimension constants +static constexpr int OMEGA_H_VERTICES = 0; +static constexpr int OMEGA_H_EDGE = 1; +static constexpr int OMEGA_H_FACE = 2; +static constexpr int OMEGA_H_REGION = 3; + +// Tag names used to identify geometric classification stored on the mesh +// These match the convention written by Omega_h's gmsh/vtk readers when +// class_id / class_dim tags are present. + +static constexpr const char *CLASS_DIM_TAG = "class_dim"; +static constexpr const char *CLASS_ID_TAG = "class_id"; + +OmegaHMeshManager::OmegaHMeshManager() { + library_ = std::make_unique(nullptr, nullptr); + mesh_ = std::make_unique(library_.get()); +} + +void OmegaHMeshManager::load_file(const std::string &file_path) { + *mesh_ = Omega_h::binary::read(file_path, library_.get()); +} + +void OmegaHMeshManager::init() { + if (mesh()->dim() != 3) + fatal_error("Mesh must be 3-dimensional"); + + bool has_class = mesh_->has_tag(OMEGA_H_REGION, CLASS_DIM_TAG) && + mesh_->has_tag(OMEGA_H_REGION, CLASS_ID_TAG); + + if (has_class) { + + auto region_class_dim = + mesh_->get_array(OMEGA_H_REGION, CLASS_DIM_TAG); + auto region_class_id = + mesh_->get_array(OMEGA_H_REGION, CLASS_ID_TAG); + + std::set vol_set; + for (Omega_h::LO i = 0; i < mesh_->nregions(); ++i) { + if (region_class_dim.get(i) == 3) + vol_set.insert(static_cast(region_class_id.get(i))); + } + for (auto v : vol_set) + volumes_.push_back(v); + + if (mesh_->has_tag(OMEGA_H_FACE, CLASS_DIM_TAG) && + mesh_->has_tag(OMEGA_H_FACE, CLASS_ID_TAG)) { + auto face_class_dim = + mesh_->get_array(OMEGA_H_FACE, CLASS_DIM_TAG); + auto face_class_id = + mesh_->get_array(OMEGA_H_FACE, CLASS_ID_TAG); + + std::set surf_set; + for (Omega_h::LO i = 0; i < mesh_->nfaces(); ++i) { + if (face_class_dim.get(i) == 2) + surf_set.insert(static_cast(face_class_id.get(i))); + } + for (auto s : surf_set) + surfaces_.push_back(s); + } + + if (!surfaces_.empty()) { + auto face_to_region = mesh_->ask_up(OMEGA_H_FACE, OMEGA_H_REGION); + auto face_class_dim = + mesh_->get_array(OMEGA_H_FACE, CLASS_DIM_TAG); + auto face_class_id = + mesh_->get_array(OMEGA_H_FACE, CLASS_ID_TAG); + auto region_class_id = + mesh_->get_array(OMEGA_H_REGION, CLASS_ID_TAG); + + for (auto surf_id : surfaces_) { + surface_senses_[surf_id] = {ID_NONE, ID_NONE}; + } + + for (Omega_h::LO fi = 0; fi < mesh_->nfaces(); ++fi) { + if (face_class_dim.get(fi) != 2) + continue; + MeshID surf_id = static_cast(face_class_id.get(fi)); + + Omega_h::LO begin = face_to_region.a2ab.get(fi); + Omega_h::LO end = face_to_region.a2ab.get(fi + 1); + for (Omega_h::LO k = begin; k < end; ++k) { + Omega_h::LO ri = face_to_region.ab2b.get(k); + MeshID vol_id = static_cast(region_class_id.get(ri)); + auto &senses = surface_senses_[surf_id]; + if (senses.first == ID_NONE) + senses.first = vol_id; // forward + else if (senses.second == ID_NONE && senses.first != vol_id) + senses.second = vol_id; // reverse + } + } + } + } else { + MeshID vol_id = create_volume(); + MeshID surf_id = next_surface_id_++; + surfaces_.push_back(surf_id); + surface_senses_[surf_id] = {vol_id, ID_NONE}; + add_surface_to_volume(vol_id, surf_id, Sense::FORWARD); + } + + create_implicit_complement(); +} + +int OmegaHMeshManager::num_volumes() const { + return static_cast(mesh_->nelems()); +} + +int OmegaHMeshManager::num_surfaces() const { + return static_cast(surfaces_.size()); +} + +int OmegaHMeshManager::num_ents_of_dimension(int dim) const { + return static_cast(mesh_->nents(dim)); +} + +MeshID OmegaHMeshManager::create_volume() { + MeshID vol_id = next_volume_id_++; + volumes_.push_back(vol_id); + return vol_id; +} + +void OmegaHMeshManager::add_surface_to_volume(MeshID volume, MeshID surface, + Sense sense, bool overwrite) { + auto it = surface_senses_.find(surface); + if (it == surface_senses_.end()) { + surface_senses_[surface] = {ID_NONE, ID_NONE}; + it = surface_senses_.find(surface); + } + + auto &senses = it->second; + if (sense == Sense::FORWARD) { + if (senses.first != ID_NONE && !overwrite) + fatal_error("Surface {} forward sense already set", surface); + senses.first = volume; + } else if (sense == Sense::REVERSE) { + if (senses.second != ID_NONE && !overwrite) + fatal_error("Surface {} reverse sense already set", surface); + senses.second = volume; + } else { + fatal_error("Invalid sense for surface {}", surface); + } +} + +int OmegaHMeshManager::num_volume_elements(MeshID volume) const { + return static_cast(get_volume_elements(volume).size()); +} + +int OmegaHMeshManager::num_volume_elements() const { + // in omega_h number of regions = number of elements + // are maybe same. weird + return static_cast(mesh_->nregions()); +} + +int OmegaHMeshManager::num_volume_faces(MeshID volume) const { + int total = 0; + for (auto surf : get_volume_surfaces(volume)) + total += num_surface_faces(surf); + return total; +} + +int OmegaHMeshManager::num_surface_faces(MeshID surface) const { + return static_cast(get_surface_faces(surface).size()); +} + +int OmegaHMeshManager::num_vertices() const { + return static_cast(mesh_->nverts()); +} + +std::vector +OmegaHMeshManager::get_volume_elements(MeshID volume) const { + // Elements classified on a given volume have class_dim==3 and + // class_id==volume + std::vector result; + + auto class_dim = mesh_->get_array(OMEGA_H_REGION, CLASS_DIM_TAG); + auto class_id = mesh_->get_array(OMEGA_H_REGION, CLASS_ID_TAG); + for (Omega_h::LO i = 0; i < mesh_->nregions(); ++i) { + if (class_dim.get(i) == 3 && + static_cast(class_id.get(i)) == volume) { + result.push_back(static_cast(i)); + } + } + + return result; +} + +std::vector OmegaHMeshManager::get_surface_faces(MeshID surface) const { + std::vector result; + if (!mesh_->has_tag(OMEGA_H_FACE, CLASS_ID_TAG)) { + + auto face_to_region = mesh_->ask_up(OMEGA_H_FACE, OMEGA_H_REGION); + for (Omega_h::LO fi = 0; fi < mesh_->nfaces(); ++fi) { + Omega_h::LO count = + face_to_region.a2ab.get(fi + 1) - face_to_region.a2ab.get(fi); + if (count == 1) // exposed face + result.push_back(static_cast(fi)); + } + return result; + } + + auto class_dim = mesh_->get_array(OMEGA_H_FACE, CLASS_DIM_TAG); + auto class_id = mesh_->get_array(OMEGA_H_FACE, CLASS_ID_TAG); + for (Omega_h::LO fi = 0; fi < mesh_->nfaces(); ++fi) { + if (class_dim.get(fi) == 2 && + static_cast(class_id.get(fi)) == surface) + result.push_back(static_cast(fi)); + } + return result; +} + +std::vector +OmegaHMeshManager::element_connectivity(MeshID element) const { + // tet vertex indices: 4 verts per tet stored in ask_elem_verts() + auto ev2v = mesh_->ask_elem_verts(); // LOs, size = 4 * nregions + constexpr int verts_per_tet = 4; + std::vector conn(verts_per_tet); + Omega_h::LO base = static_cast(element) * verts_per_tet; + for (int i = 0; i < verts_per_tet; ++i) + conn[i] = static_cast(ev2v.get(base + i)); + return conn; +} + +std::vector OmegaHMeshManager::face_connectivity(MeshID face) const { + // tri vertex indices: 3 verts per triangle stored in ask_verts_of(2) + auto fv2v = mesh_->ask_verts_of(OMEGA_H_FACE); // LOs, size = 3 * n_faces + constexpr int verts_per_tri = 3; + std::vector conn(verts_per_tri); + Omega_h::LO base = static_cast(face) * verts_per_tri; + for (int i = 0; i < verts_per_tri; ++i) + conn[i] = static_cast(fv2v.get(base + i)); + return conn; +} + +Vertex OmegaHMeshManager::vertex_coordinates(MeshID vertex_id) const { + // coords() returns a flat array: [x0,y0,z0, x1,y1,z1, ...] + auto coords = mesh_->coords(); + Omega_h::LO base = static_cast(vertex_id) * 3; + return Vertex(coords.get(base), coords.get(base + 1), coords.get(base + 2)); +} + +std::vector OmegaHMeshManager::element_vertices(MeshID element) const { + auto conn = element_connectivity(element); + std::vector verts(conn.size()); + for (std::size_t i = 0; i < conn.size(); ++i) + verts[i] = vertex_coordinates(conn[i]); + return verts; +} + +std::array OmegaHMeshManager::face_vertices(MeshID face) const { + auto fv2v = mesh_->ask_verts_of(OMEGA_H_FACE); + constexpr int verts_per_tri = 3; + Omega_h::LO base = static_cast(face) * verts_per_tri; + auto coords = mesh_->coords(); + std::array verts; + for (int i = 0; i < verts_per_tri; ++i) { + Omega_h::LO vid = fv2v.get(base + i); + Omega_h::LO vbase = vid * 3; + verts[i] = + Vertex(coords.get(vbase), coords.get(vbase + 1), coords.get(vbase + 2)); + } + return verts; +} + +SurfaceElementType +OmegaHMeshManager::get_surface_element_type(MeshID /*surface*/) const { + // Omega_h simplex meshes always use triangular surface elements + return SurfaceElementType::TRI; +} + +MeshID OmegaHMeshManager::adjacent_element(MeshID element, int face) const { + // ask_up(face→region) gives us the regions adjacent to each face. + // First find the face index for the given local face of the element, + // then look at its upward adjacency. + auto er2f = + mesh_->ask_down(OMEGA_H_REGION, OMEGA_H_FACE); // region→face downward + constexpr int faces_per_tet = 4; + Omega_h::LO base = static_cast(element) * faces_per_tet; + Omega_h::LO fi = er2f.ab2b.get(base + face); // global face index + + auto f2r = mesh_->ask_up(OMEGA_H_FACE, OMEGA_H_REGION); + Omega_h::LO begin = f2r.a2ab.get(fi); + Omega_h::LO end = f2r.a2ab.get(fi + 1); + + for (Omega_h::LO k = begin; k < end; ++k) { + Omega_h::LO ri = f2r.ab2b.get(k); + if (static_cast(ri) != element) + return static_cast(ri); + } + return ID_NONE; // boundary face — no neighbor +} + +double OmegaHMeshManager::element_volume(MeshID element) const { + auto verts = element_vertices(element); + std::array tet_verts; + for (int i = 0; i < 4; ++i) + tet_verts[i] = verts[i]; + return tetrahedron_volume(tet_verts); +} + +std::pair +OmegaHMeshManager::surface_senses(MeshID surface) const { + auto it = surface_senses_.find(surface); + if (it == surface_senses_.end()) + return {ID_NONE, ID_NONE}; + return it->second; +} + +std::vector +OmegaHMeshManager::get_volume_surfaces(MeshID volume) const { + std::vector result; + for (const auto &[surf, senses] : surface_senses_) { + if (senses.first == volume || senses.second == volume) + result.push_back(surf); + } + return result; +} + +Sense OmegaHMeshManager::surface_sense(MeshID surface, MeshID volume) const { + auto senses = surface_senses(surface); + if (senses.first == volume) + return Sense::FORWARD; + if (senses.second == volume) + return Sense::REVERSE; + fatal_error("Volume {} is not a parent of surface {}", volume, surface); + return Sense::UNSET; +} + +} // namespace xdg \ No newline at end of file From 638dd41d99caeccc51d75c0bb6372fc404d2d0b0 Mon Sep 17 00:00:00 2001 From: Ebny Walid Ahammed <69362074+magnoxemo@users.noreply.github.com> Date: Mon, 15 Jun 2026 15:22:58 -0500 Subject: [PATCH 09/26] implementing tests --- Dockerfile | 67 +++++ include/xdg/constants.h | 6 +- include/xdg/omega_h/mesh_manager.h | 132 +++++++-- src/omega_h/mesh_manager.cpp | 455 +++++++++++++++-------------- tests/CMakeLists.txt | 8 + tests/test_omega_h.cpp | 324 ++++++++++++++++++++ 6 files changed, 748 insertions(+), 244 deletions(-) create mode 100644 Dockerfile create mode 100644 tests/test_omega_h.cpp diff --git a/Dockerfile b/Dockerfile new file mode 100644 index 00000000..4953a2a5 --- /dev/null +++ b/Dockerfile @@ -0,0 +1,67 @@ +FROM ubuntu:24.04 + +ARG BUILD_JOBS=1 + +# installing dependencies +RUN apt-get update --yes +RUN apt-get install --yes \ + git \ + make \ + autoconf \ + automake \ + libtool \ + flex \ + bison \ + cmake \ + g++ \ + gfortran \ + libhdf5-dev \ + libopenblas-dev \ + mpich \ + libmpich-dev \ + libtirpc-dev \ + curl \ + nano \ + libembree-dev \ + pkg-config + +# compile libmesh from source and install it +RUN git clone --recurse-submodules https://github.com/libMesh/libmesh.git /XDG_TEST_SYSTEM/libmesh +WORKDIR /XDG_TEST_SYSTEM/libmesh/build +RUN ../configure --prefix=/XDG_TEST_SYSTEM/libmesh_install_dir +RUN make -j $(BUILD_JOBS) +RUN make install + +ENV LIBMESH_INSTALL_PATH=/XDG_TEST_SYSTEM/libmesh_install_dir + +# build MOAB from source +RUN git clone --recurse-submodules https://bitbucket.org/fathomteam/moab.git /XDG_TEST_SYSTEM/moab +WORKDIR /XDG_TEST_SYSTEM/moab/build +RUN cmake .. \ + -DCMAKE_INSTALL_PREFIX=/XDG_TEST_SYSTEM/moab_install_dir \ + -DENABLE_HDF5=ON \ + -DHDF5_ROOT=/usr \ + -DBLAS_LIBRARIES=/usr/lib/x86_64-linux-gnu/libopenblas.so \ + -DBUILD_SHARED_LIBS=ON +RUN make -j $(BUILD_JOBS) +RUN make install + +ENV MOAB_INSTALL_PATH=/XDG_TEST_SYSTEM/moab_install_dir + + +# build XDG +RUN git clone --recurse-submodules https://github.com/xdg-org/xdg.git /XDG_TEST_SYSTEM/xdg +WORKDIR /XDG_TEST_SYSTEM/xdg/build +RUN cmake .. \ + -DCMAKE_INSTALL_PREFIX=/XDG_TEST_SYSTEM/xdg_install_dir \ + -DCMAKE_BUILD_TYPE=Release \ + -DCMAKE_C_COMPILER=mpicc \ + -DCMAKE_CXX_COMPILER=mpicxx \ + -DXDG_ENABLE_MOAB=ON \ + -DMOAB_DIR=${MOAB_INSTALL_PATH} \ + -DXDG_ENABLE_LIBMESH=ON \ + -DLIBMESH_DIR=${LIBMESH_INSTALL_PATH} +RUN make -j $(BUILD_JOBS) +RUN make install + +ENV XDG_INSTALL_PATH=/XDG_TEST_SYSTEM/xdg_install_dir diff --git a/include/xdg/constants.h b/include/xdg/constants.h index 38df9db1..a668e8f3 100644 --- a/include/xdg/constants.h +++ b/include/xdg/constants.h @@ -48,7 +48,8 @@ enum class Sense { enum class MeshLibrary { MOCK = 0, // mock testing interface MOAB, - LIBMESH + LIBMESH, + OMEGA_H }; // Ray Tracing library identifier @@ -61,7 +62,8 @@ static const std::map MESH_LIB_TO_STR = { {MeshLibrary::MOCK, "MOCK"}, {MeshLibrary::MOAB, "MOAB"}, - {MeshLibrary::LIBMESH, "LIBMESH"} + {MeshLibrary::LIBMESH, "LIBMESH"}, + {MeshLibrary::OMEGA_H,"OMEGA_H"} }; static const std::map RT_LIB_TO_STR = diff --git a/include/xdg/omega_h/mesh_manager.h b/include/xdg/omega_h/mesh_manager.h index 1349a80a..5e2185ca 100644 --- a/include/xdg/omega_h/mesh_manager.h +++ b/include/xdg/omega_h/mesh_manager.h @@ -1,72 +1,164 @@ #ifndef XDG_OMEGA_H_MESH_MANAGER_H #define XDG_OMEGA_H_MESH_MANAGER_H +#include #include #include #include +#include +#include "xdg/constants.h" +#include "xdg/element_face_accessor.h" #include "xdg/mesh_manager_interface.h" #include "Omega_h_library.hpp" -#include "Omega_h_map.hpp" #include "Omega_h_mesh.hpp" -#include "Omega_h_profile.hpp" -#include "Omega_h_quality.hpp" -#include "Omega_h_timer.hpp" namespace xdg { +//! \brief Mesh manager backed by the Omega_h simplex mesh library. +//! +//! Omega_h tetrahedral meshes carry the geometric model they were generated +//! from as "class_dim"/"class_id" classification tags on each entity. These +//! tags map mesh entities onto the volumes (3D model entities) and surfaces +//! (2D model entities) of the geometry, which is all the topology XDG needs. +//! When a mesh has no classification, the whole mesh is treated as a single +//! volume bounded by its exposed faces. class OmegaHMeshManager : public MeshManager { public: OmegaHMeshManager(); - ~OmegaHMeshManager() = default; + + OmegaHMeshManager(const Omega_h::Mesh* mesh); + + ~OmegaHMeshManager() override = default; + + // Interface methods + MeshLibrary mesh_library() const override { return MeshLibrary::OMEGA_H; } void load_file(const std::string &file_path) override; void init() override; - int num_volumes() const override; - int num_surfaces() const override; - int num_ents_of_dimension(int dim) const; + void parse_metadata() override {} - MeshID create_volume() override; - void add_surface_to_volume(MeshID volume, MeshID surface, Sense sense, - bool overwrite = false) override; + int num_volumes() const override { return volumes_.size(); } + + int num_surfaces() const override { return surfaces_.size(); } + + int num_ents_of_dimension(int dim) const override { + switch (dim) { + case 3: + return num_volumes(); + case 2: + return num_surfaces(); + default: + return 0; + } + } + + int num_volume_elements(MeshID volume) const override { + return get_volume_elements(volume).size(); + } + + int num_volume_elements() const override { return num_elements_; } + + int num_volume_faces(MeshID volume) const override { + int count = 0; + for (auto surface : get_volume_surfaces(volume)) { + count += num_surface_faces(surface); + } + return count; + } + + int num_surface_faces(MeshID surface) const override { + return get_surface_faces(surface).size(); + } - int num_volume_elements(MeshID volume) const override; - int num_volume_elements() const override; - int num_volume_faces(MeshID volume) const override; - int num_surface_faces(MeshID surface) const override; int num_vertices() const override; std::vector get_volume_elements(MeshID volume) const override; + std::vector get_surface_faces(MeshID surface) const override; + std::vector element_connectivity(MeshID element) const override; + std::vector face_connectivity(MeshID face) const override; + + MeshID get_boundary_face_element(MeshID face) const override; + Vertex vertex_coordinates(MeshID vertex) const override; + std::vector element_vertices(MeshID element) const override; + std::array face_vertices(MeshID face) const override; - SurfaceElementType get_surface_element_type(MeshID surface) const override; + + //! \brief Vertices of a local face of a tetrahedral element + std::array element_face_vertices(MeshID element, + int local_face) const; + + SurfaceElementType + get_surface_element_type(MeshID /*surface*/) const override { + // Omega_h simplex meshes always use triangular surface elements + return SurfaceElementType::TRI; + } + MeshID adjacent_element(MeshID element, int face) const override; + double element_volume(MeshID element) const override; + MeshID create_volume() override; + + void add_surface_to_volume(MeshID volume, MeshID surface, Sense sense, + bool overwrite = false) override; + std::pair surface_senses(MeshID surface) const override; + std::vector get_volume_surfaces(MeshID volume) const override; - Sense surface_sense(MeshID surface, MeshID volume) const override; - void parse_metadata(); + Sense surface_sense(MeshID surface, MeshID volume) const override; + // Accessors Omega_h::Mesh *mesh() { return mesh_.get(); } const Omega_h::Mesh *mesh() const { return mesh_.get(); } private: + //! \brief Whether the mesh carries volume and surface classification tags + bool has_classification() const; + + //! \brief Populate volumes_ and surfaces_ from the classification tags + void discover_geometry(); + + //! \brief Treat the whole mesh as a single volume bounded by its exposed + //! faces + void discover_single_volume(); + + //! \brief Assign the parent volumes (forward/reverse senses) of each surface + void determine_surface_senses(); + + //! \brief Map element and vertex ID spaces into contiguous index spaces + void map_id_spaces(); + std::unique_ptr library_; std::unique_ptr mesh_; + //! Mapping of surfaces to the volumes on either side. Volumes are ordered + //! based on their sense with respect to the surface triangles std::map> surface_senses_; - MeshID next_volume_id_{1}; - MeshID next_surface_id_{1}; + int32_t num_elements_{-1}; +}; + +//! \brief Face-vertex accessor for Omega_h tetrahedral elements +struct OmegaHElementFaceAccessor : public ElementFaceAccessor { + OmegaHElementFaceAccessor(const OmegaHMeshManager *mesh_manager, + MeshID element) + : ElementFaceAccessor(element), mesh_manager_(mesh_manager) {} + + std::array face_vertices(int i) const override { + return mesh_manager_->element_face_vertices(element_, i); + } + + const OmegaHMeshManager *mesh_manager_; }; } // namespace xdg diff --git a/src/omega_h/mesh_manager.cpp b/src/omega_h/mesh_manager.cpp index e6cdfee3..fb697217 100644 --- a/src/omega_h/mesh_manager.cpp +++ b/src/omega_h/mesh_manager.cpp @@ -1,37 +1,37 @@ #include "xdg/omega_h/mesh_manager.h" #include -#include +#include #include -#include -#include #include #include #include "xdg/error.h" #include "xdg/geometry/measure.h" -#include "xdg/moab/tag_conventions.h" -#include "xdg/util/str_utils.h" #include "Omega_h_file.hpp" -#include "Omega_h_for.hpp" -#include "Omega_h_mesh.hpp" namespace xdg { -// Omega_h dimension constants -static constexpr int OMEGA_H_VERTICES = 0; -static constexpr int OMEGA_H_EDGE = 1; +// Topological entity dimensions in Omega_h static constexpr int OMEGA_H_FACE = 2; static constexpr int OMEGA_H_REGION = 3; -// Tag names used to identify geometric classification stored on the mesh -// These match the convention written by Omega_h's gmsh/vtk readers when -// class_id / class_dim tags are present. +// Geometric model dimension a mesh entity may be classified on +static constexpr int CLASS_DIM_SURFACE = 2; +static constexpr int CLASS_DIM_VOLUME = 3; +// Classification tag names written by Omega_h static constexpr const char *CLASS_DIM_TAG = "class_dim"; static constexpr const char *CLASS_ID_TAG = "class_id"; +static constexpr int VERTS_PER_TET = 4; +static constexpr int VERTS_PER_TRI = 3; +static constexpr int FACES_PER_TET = 4; + +// Constructor + +OmegaHMeshManager::OmegaHMeshManager(const Omega_h::Mesh* mesh) {} OmegaHMeshManager::OmegaHMeshManager() { library_ = std::make_unique(nullptr, nullptr); mesh_ = std::make_unique(library_.get()); @@ -42,307 +42,318 @@ void OmegaHMeshManager::load_file(const std::string &file_path) { } void OmegaHMeshManager::init() { - if (mesh()->dim() != 3) + // XDG operates on 3-dimensional volume meshes + if (mesh()->dim() != 3) { fatal_error("Mesh must be 3-dimensional"); + } - bool has_class = mesh_->has_tag(OMEGA_H_REGION, CLASS_DIM_TAG) && - mesh_->has_tag(OMEGA_H_REGION, CLASS_ID_TAG); - - if (has_class) { - - auto region_class_dim = - mesh_->get_array(OMEGA_H_REGION, CLASS_DIM_TAG); - auto region_class_id = - mesh_->get_array(OMEGA_H_REGION, CLASS_ID_TAG); - - std::set vol_set; - for (Omega_h::LO i = 0; i < mesh_->nregions(); ++i) { - if (region_class_dim.get(i) == 3) - vol_set.insert(static_cast(region_class_id.get(i))); - } - for (auto v : vol_set) - volumes_.push_back(v); - - if (mesh_->has_tag(OMEGA_H_FACE, CLASS_DIM_TAG) && - mesh_->has_tag(OMEGA_H_FACE, CLASS_ID_TAG)) { - auto face_class_dim = - mesh_->get_array(OMEGA_H_FACE, CLASS_DIM_TAG); - auto face_class_id = - mesh_->get_array(OMEGA_H_FACE, CLASS_ID_TAG); - - std::set surf_set; - for (Omega_h::LO i = 0; i < mesh_->nfaces(); ++i) { - if (face_class_dim.get(i) == 2) - surf_set.insert(static_cast(face_class_id.get(i))); - } - for (auto s : surf_set) - surfaces_.push_back(s); - } - - if (!surfaces_.empty()) { - auto face_to_region = mesh_->ask_up(OMEGA_H_FACE, OMEGA_H_REGION); - auto face_class_dim = - mesh_->get_array(OMEGA_H_FACE, CLASS_DIM_TAG); - auto face_class_id = - mesh_->get_array(OMEGA_H_FACE, CLASS_ID_TAG); - auto region_class_id = - mesh_->get_array(OMEGA_H_REGION, CLASS_ID_TAG); - - for (auto surf_id : surfaces_) { - surface_senses_[surf_id] = {ID_NONE, ID_NONE}; - } + // in Omega_h the regions (3D simplices) are the volume elements + num_elements_ = mesh_->nregions(); - for (Omega_h::LO fi = 0; fi < mesh_->nfaces(); ++fi) { - if (face_class_dim.get(fi) != 2) - continue; - MeshID surf_id = static_cast(face_class_id.get(fi)); - - Omega_h::LO begin = face_to_region.a2ab.get(fi); - Omega_h::LO end = face_to_region.a2ab.get(fi + 1); - for (Omega_h::LO k = begin; k < end; ++k) { - Omega_h::LO ri = face_to_region.ab2b.get(k); - MeshID vol_id = static_cast(region_class_id.get(ri)); - auto &senses = surface_senses_[surf_id]; - if (senses.first == ID_NONE) - senses.first = vol_id; // forward - else if (senses.second == ID_NONE && senses.first != vol_id) - senses.second = vol_id; // reverse - } - } - } + // a classified mesh defines its volumes and surfaces directly through the + // class_dim/class_id tags. Otherwise treat the entire mesh as a single volume + // bounded by its exposed faces. + if (has_classification()) { + discover_geometry(); + determine_surface_senses(); } else { - MeshID vol_id = create_volume(); - MeshID surf_id = next_surface_id_++; - surfaces_.push_back(surf_id); - surface_senses_[surf_id] = {vol_id, ID_NONE}; - add_surface_to_volume(vol_id, surf_id, Sense::FORWARD); + discover_single_volume(); } + // create an implicit complement to bound the model create_implicit_complement(); -} -int OmegaHMeshManager::num_volumes() const { - return static_cast(mesh_->nelems()); + // map ID spaces into indices for ordered access by downstream applications + map_id_spaces(); } -int OmegaHMeshManager::num_surfaces() const { - return static_cast(surfaces_.size()); +bool OmegaHMeshManager::has_classification() const { + return mesh_->has_tag(OMEGA_H_REGION, CLASS_DIM_TAG) && + mesh_->has_tag(OMEGA_H_REGION, CLASS_ID_TAG) && + mesh_->has_tag(OMEGA_H_FACE, CLASS_DIM_TAG) && + mesh_->has_tag(OMEGA_H_FACE, CLASS_ID_TAG); } -int OmegaHMeshManager::num_ents_of_dimension(int dim) const { - return static_cast(mesh_->nents(dim)); +void OmegaHMeshManager::discover_geometry() { + // volumes are the unique class IDs of regions classified on a 3D model entity + auto region_class_dim = + mesh_->get_array(OMEGA_H_REGION, CLASS_DIM_TAG); + auto region_class_id = + mesh_->get_array(OMEGA_H_REGION, CLASS_ID_TAG); + + std::set volume_ids; + for (Omega_h::LO region = 0; region < mesh_->nregions(); ++region) { + if (region_class_dim.get(region) == CLASS_DIM_VOLUME) { + volume_ids.insert(static_cast(region_class_id.get(region))); + } + } + volumes_.assign(volume_ids.begin(), volume_ids.end()); + + // surfaces are the unique class IDs of faces classified on a 2D model entity + auto face_class_dim = + mesh_->get_array(OMEGA_H_FACE, CLASS_DIM_TAG); + auto face_class_id = + mesh_->get_array(OMEGA_H_FACE, CLASS_ID_TAG); + + std::set surface_ids; + for (Omega_h::LO face = 0; face < mesh_->nfaces(); ++face) { + if (face_class_dim.get(face) == CLASS_DIM_SURFACE) { + surface_ids.insert(static_cast(face_class_id.get(face))); + } + } + surfaces_.assign(surface_ids.begin(), surface_ids.end()); } -MeshID OmegaHMeshManager::create_volume() { - MeshID vol_id = next_volume_id_++; - volumes_.push_back(vol_id); - return vol_id; +void OmegaHMeshManager::discover_single_volume() { + // one volume bounded by a single surface made up of every exposed face + MeshID volume = create_volume(); + volumes_.push_back(volume); + + MeshID surface = next_surface_id(); + surfaces_.push_back(surface); + surface_senses_[surface] = {volume, ID_NONE}; } -void OmegaHMeshManager::add_surface_to_volume(MeshID volume, MeshID surface, - Sense sense, bool overwrite) { - auto it = surface_senses_.find(surface); - if (it == surface_senses_.end()) { +void OmegaHMeshManager::determine_surface_senses() { + // the parent volumes of a surface are the volumes of the regions on either + // side of its faces. The first region encountered defines the forward sense + // and the opposing region, if any, defines the reverse sense. + auto face_to_region = mesh_->ask_up(OMEGA_H_FACE, OMEGA_H_REGION); + auto face_class_dim = + mesh_->get_array(OMEGA_H_FACE, CLASS_DIM_TAG); + auto face_class_id = + mesh_->get_array(OMEGA_H_FACE, CLASS_ID_TAG); + auto region_class_id = + mesh_->get_array(OMEGA_H_REGION, CLASS_ID_TAG); + + for (auto surface : surfaces_) { surface_senses_[surface] = {ID_NONE, ID_NONE}; - it = surface_senses_.find(surface); } - auto &senses = it->second; - if (sense == Sense::FORWARD) { - if (senses.first != ID_NONE && !overwrite) - fatal_error("Surface {} forward sense already set", surface); - senses.first = volume; - } else if (sense == Sense::REVERSE) { - if (senses.second != ID_NONE && !overwrite) - fatal_error("Surface {} reverse sense already set", surface); - senses.second = volume; - } else { - fatal_error("Invalid sense for surface {}", surface); + for (Omega_h::LO face = 0; face < mesh_->nfaces(); ++face) { + if (face_class_dim.get(face) != CLASS_DIM_SURFACE) { + continue; + } + MeshID surface = static_cast(face_class_id.get(face)); + auto &senses = surface_senses_[surface]; + + Omega_h::LO begin = face_to_region.a2ab.get(face); + Omega_h::LO end = face_to_region.a2ab.get(face + 1); + for (Omega_h::LO k = begin; k < end; ++k) { + MeshID volume = + static_cast(region_class_id.get(face_to_region.ab2b.get(k))); + if (senses.first == ID_NONE) { + senses.first = volume; // forward sense + } else if (senses.second == ID_NONE && senses.first != volume) { + senses.second = volume; // reverse sense + } + } } } -int OmegaHMeshManager::num_volume_elements(MeshID volume) const { - return static_cast(get_volume_elements(volume).size()); -} +void OmegaHMeshManager::map_id_spaces() { + // Omega_h stores entities in a contiguous, zero-based index space, so element + // and vertex IDs are identical to their indices + std::vector element_ids(mesh_->nregions()); + std::iota(element_ids.begin(), element_ids.end(), 0); + volume_element_id_map_ = IDBlockMapping(element_ids); -int OmegaHMeshManager::num_volume_elements() const { - // in omega_h number of regions = number of elements - // are maybe same. weird - return static_cast(mesh_->nregions()); + std::vector vertex_ids(mesh_->nverts()); + std::iota(vertex_ids.begin(), vertex_ids.end(), 0); + vertex_id_map_ = IDBlockMapping(vertex_ids); } -int OmegaHMeshManager::num_volume_faces(MeshID volume) const { - int total = 0; - for (auto surf : get_volume_surfaces(volume)) - total += num_surface_faces(surf); - return total; -} +MeshID OmegaHMeshManager::create_volume() { return next_volume_id(); } -int OmegaHMeshManager::num_surface_faces(MeshID surface) const { - return static_cast(get_surface_faces(surface).size()); +void OmegaHMeshManager::add_surface_to_volume(MeshID volume, MeshID surface, + Sense sense, bool overwrite) { + auto senses = surface_senses(surface); + if (sense == Sense::FORWARD) { + if (!overwrite && senses.first != ID_NONE) { + fatal_error("Surface {} already has a forward sense", surface); + } + surface_senses_[surface] = {volume, senses.second}; + } else { + if (!overwrite && senses.second != ID_NONE) { + fatal_error("Surface {} already has a reverse sense", surface); + } + surface_senses_[surface] = {senses.first, volume}; + } } -int OmegaHMeshManager::num_vertices() const { - return static_cast(mesh_->nverts()); -} +int OmegaHMeshManager::num_vertices() const { return mesh_->nverts(); } std::vector OmegaHMeshManager::get_volume_elements(MeshID volume) const { - // Elements classified on a given volume have class_dim==3 and - // class_id==volume - std::vector result; + std::vector elements; + + // without classification every region belongs to the single volume + if (!mesh_->has_tag(OMEGA_H_REGION, CLASS_ID_TAG)) { + elements.resize(mesh_->nregions()); + std::iota(elements.begin(), elements.end(), 0); + return elements; + } + // otherwise gather the regions classified on this volume auto class_dim = mesh_->get_array(OMEGA_H_REGION, CLASS_DIM_TAG); auto class_id = mesh_->get_array(OMEGA_H_REGION, CLASS_ID_TAG); - for (Omega_h::LO i = 0; i < mesh_->nregions(); ++i) { - if (class_dim.get(i) == 3 && - static_cast(class_id.get(i)) == volume) { - result.push_back(static_cast(i)); + for (Omega_h::LO region = 0; region < mesh_->nregions(); ++region) { + if (class_dim.get(region) == CLASS_DIM_VOLUME && + static_cast(class_id.get(region)) == volume) { + elements.push_back(static_cast(region)); } } - - return result; + return elements; } std::vector OmegaHMeshManager::get_surface_faces(MeshID surface) const { - std::vector result; - if (!mesh_->has_tag(OMEGA_H_FACE, CLASS_ID_TAG)) { + std::vector faces; + // without classification the surface is the set of exposed boundary faces + // (faces adjacent to exactly one region) + if (!mesh_->has_tag(OMEGA_H_FACE, CLASS_ID_TAG)) { auto face_to_region = mesh_->ask_up(OMEGA_H_FACE, OMEGA_H_REGION); - for (Omega_h::LO fi = 0; fi < mesh_->nfaces(); ++fi) { - Omega_h::LO count = - face_to_region.a2ab.get(fi + 1) - face_to_region.a2ab.get(fi); - if (count == 1) // exposed face - result.push_back(static_cast(fi)); + for (Omega_h::LO face = 0; face < mesh_->nfaces(); ++face) { + Omega_h::LO n_adjacent = + face_to_region.a2ab.get(face + 1) - face_to_region.a2ab.get(face); + if (n_adjacent == 1) { + faces.push_back(static_cast(face)); + } } - return result; + return faces; } + // otherwise gather the faces classified on this surface auto class_dim = mesh_->get_array(OMEGA_H_FACE, CLASS_DIM_TAG); auto class_id = mesh_->get_array(OMEGA_H_FACE, CLASS_ID_TAG); - for (Omega_h::LO fi = 0; fi < mesh_->nfaces(); ++fi) { - if (class_dim.get(fi) == 2 && - static_cast(class_id.get(fi)) == surface) - result.push_back(static_cast(fi)); + for (Omega_h::LO face = 0; face < mesh_->nfaces(); ++face) { + if (class_dim.get(face) == CLASS_DIM_SURFACE && + static_cast(class_id.get(face)) == surface) { + faces.push_back(static_cast(face)); + } } - return result; + return faces; } std::vector OmegaHMeshManager::element_connectivity(MeshID element) const { - // tet vertex indices: 4 verts per tet stored in ask_elem_verts() - auto ev2v = mesh_->ask_elem_verts(); // LOs, size = 4 * nregions - constexpr int verts_per_tet = 4; - std::vector conn(verts_per_tet); - Omega_h::LO base = static_cast(element) * verts_per_tet; - for (int i = 0; i < verts_per_tet; ++i) - conn[i] = static_cast(ev2v.get(base + i)); - return conn; + // tetrahedra store four vertices per element in element->vertex order + auto elem_verts = mesh_->ask_elem_verts(); + Omega_h::LO base = static_cast(element) * VERTS_PER_TET; + std::vector connectivity(VERTS_PER_TET); + for (int i = 0; i < VERTS_PER_TET; ++i) { + connectivity[i] = static_cast(elem_verts.get(base + i)); + } + return connectivity; } std::vector OmegaHMeshManager::face_connectivity(MeshID face) const { - // tri vertex indices: 3 verts per triangle stored in ask_verts_of(2) - auto fv2v = mesh_->ask_verts_of(OMEGA_H_FACE); // LOs, size = 3 * n_faces - constexpr int verts_per_tri = 3; - std::vector conn(verts_per_tri); - Omega_h::LO base = static_cast(face) * verts_per_tri; - for (int i = 0; i < verts_per_tri; ++i) - conn[i] = static_cast(fv2v.get(base + i)); - return conn; + // triangles store three vertices per face in face->vertex order + auto face_verts = mesh_->ask_verts_of(OMEGA_H_FACE); + Omega_h::LO base = static_cast(face) * VERTS_PER_TRI; + std::vector connectivity(VERTS_PER_TRI); + for (int i = 0; i < VERTS_PER_TRI; ++i) { + connectivity[i] = static_cast(face_verts.get(base + i)); + } + return connectivity; +} + +MeshID OmegaHMeshManager::get_boundary_face_element(MeshID face) const { + // the owning element of a boundary face is its single adjacent region + auto face_to_region = mesh_->ask_up(OMEGA_H_FACE, OMEGA_H_REGION); + Omega_h::LO begin = face_to_region.a2ab.get(face); + Omega_h::LO end = face_to_region.a2ab.get(face + 1); + if (begin == end) { + return ID_NONE; + } + return static_cast(face_to_region.ab2b.get(begin)); } -Vertex OmegaHMeshManager::vertex_coordinates(MeshID vertex_id) const { - // coords() returns a flat array: [x0,y0,z0, x1,y1,z1, ...] +Vertex OmegaHMeshManager::vertex_coordinates(MeshID vertex) const { + // coords() is a flat array laid out as [x0, y0, z0, x1, y1, z1, ...] auto coords = mesh_->coords(); - Omega_h::LO base = static_cast(vertex_id) * 3; - return Vertex(coords.get(base), coords.get(base + 1), coords.get(base + 2)); + Omega_h::LO base = static_cast(vertex) * 3; + return {coords.get(base), coords.get(base + 1), coords.get(base + 2)}; } std::vector OmegaHMeshManager::element_vertices(MeshID element) const { - auto conn = element_connectivity(element); - std::vector verts(conn.size()); - for (std::size_t i = 0; i < conn.size(); ++i) - verts[i] = vertex_coordinates(conn[i]); - return verts; + auto connectivity = element_connectivity(element); + std::vector vertices; + vertices.reserve(connectivity.size()); + for (auto vertex : connectivity) { + vertices.push_back(vertex_coordinates(vertex)); + } + return vertices; } std::array OmegaHMeshManager::face_vertices(MeshID face) const { - auto fv2v = mesh_->ask_verts_of(OMEGA_H_FACE); - constexpr int verts_per_tri = 3; - Omega_h::LO base = static_cast(face) * verts_per_tri; - auto coords = mesh_->coords(); - std::array verts; - for (int i = 0; i < verts_per_tri; ++i) { - Omega_h::LO vid = fv2v.get(base + i); - Omega_h::LO vbase = vid * 3; - verts[i] = - Vertex(coords.get(vbase), coords.get(vbase + 1), coords.get(vbase + 2)); + auto connectivity = face_connectivity(face); + std::array vertices; + for (int i = 0; i < VERTS_PER_TRI; ++i) { + vertices[i] = vertex_coordinates(connectivity[i]); } - return verts; + return vertices; } -SurfaceElementType -OmegaHMeshManager::get_surface_element_type(MeshID /*surface*/) const { - // Omega_h simplex meshes always use triangular surface elements - return SurfaceElementType::TRI; +std::array +OmegaHMeshManager::element_face_vertices(MeshID element, int local_face) const { + // resolve the global face index for the requested local face, then return its + // vertices in their natural (face->vertex) orientation + auto region_to_face = mesh_->ask_down(OMEGA_H_REGION, OMEGA_H_FACE); + Omega_h::LO base = static_cast(element) * FACES_PER_TET; + Omega_h::LO global_face = region_to_face.ab2b.get(base + local_face); + return face_vertices(static_cast(global_face)); } MeshID OmegaHMeshManager::adjacent_element(MeshID element, int face) const { - // ask_up(face→region) gives us the regions adjacent to each face. - // First find the face index for the given local face of the element, - // then look at its upward adjacency. - auto er2f = - mesh_->ask_down(OMEGA_H_REGION, OMEGA_H_FACE); // region→face downward - constexpr int faces_per_tet = 4; - Omega_h::LO base = static_cast(element) * faces_per_tet; - Omega_h::LO fi = er2f.ab2b.get(base + face); // global face index - - auto f2r = mesh_->ask_up(OMEGA_H_FACE, OMEGA_H_REGION); - Omega_h::LO begin = f2r.a2ab.get(fi); - Omega_h::LO end = f2r.a2ab.get(fi + 1); - + // resolve the global face index for the requested local face of the element + auto region_to_face = mesh_->ask_down(OMEGA_H_REGION, OMEGA_H_FACE); + Omega_h::LO base = static_cast(element) * FACES_PER_TET; + Omega_h::LO global_face = region_to_face.ab2b.get(base + face); + + // the neighbor is the other region sharing that face, if any + auto face_to_region = mesh_->ask_up(OMEGA_H_FACE, OMEGA_H_REGION); + Omega_h::LO begin = face_to_region.a2ab.get(global_face); + Omega_h::LO end = face_to_region.a2ab.get(global_face + 1); for (Omega_h::LO k = begin; k < end; ++k) { - Omega_h::LO ri = f2r.ab2b.get(k); - if (static_cast(ri) != element) - return static_cast(ri); + MeshID neighbor = static_cast(face_to_region.ab2b.get(k)); + if (neighbor != element) { + return neighbor; + } } - return ID_NONE; // boundary face — no neighbor + return ID_NONE; // boundary face, no neighbor } double OmegaHMeshManager::element_volume(MeshID element) const { - auto verts = element_vertices(element); - std::array tet_verts; - for (int i = 0; i < 4; ++i) - tet_verts[i] = verts[i]; - return tetrahedron_volume(tet_verts); + auto vertices = element_vertices(element); + std::array tet{vertices[0], vertices[1], vertices[2], vertices[3]}; + return tetrahedron_volume(tet); } std::pair OmegaHMeshManager::surface_senses(MeshID surface) const { auto it = surface_senses_.find(surface); - if (it == surface_senses_.end()) + if (it == surface_senses_.end()) { return {ID_NONE, ID_NONE}; + } return it->second; } std::vector OmegaHMeshManager::get_volume_surfaces(MeshID volume) const { + // walk the surface senses and return the surfaces bounding this volume std::vector result; - for (const auto &[surf, senses] : surface_senses_) { - if (senses.first == volume || senses.second == volume) - result.push_back(surf); + for (const auto &[surface, senses] : surface_senses_) { + if (senses.first == volume || senses.second == volume) { + result.push_back(surface); + } } return result; } Sense OmegaHMeshManager::surface_sense(MeshID surface, MeshID volume) const { auto senses = surface_senses(surface); - if (senses.first == volume) - return Sense::FORWARD; - if (senses.second == volume) - return Sense::REVERSE; - fatal_error("Volume {} is not a parent of surface {}", volume, surface); - return Sense::UNSET; + return volume == senses.first ? Sense::FORWARD : Sense::REVERSE; } } // namespace xdg \ No newline at end of file diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 21691467..0026cba8 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -42,6 +42,10 @@ if (XDG_ENABLE_MOAB AND XDG_BUILD_TOOLS) list(APPEND TEST_NAMES test_overlap_check) endif() +if (XDG_ENABLE_OMEGA_H) + list(APPEND TEST_NAMES test_omega_h) +endif () + foreach(test ${TEST_NAMES}) add_executable(${test} test_main.cpp ${test}.cpp) target_link_libraries(${test} PRIVATE xdg fmt::fmt Catch2::Catch2) @@ -54,6 +58,10 @@ foreach(test ${TEST_NAMES}) if (XDG_ENABLE_MOAB) target_link_libraries(${test} PRIVATE MOAB) endif() + if (XDG_ENABLE_OMEGA_H) + target_link_libraries(${test} PRIVATE Omega_h::omega_h) + endif () + set_target_properties(${test} PROPERTIES BUILD_RPATH "$") catch_discover_tests(${test} diff --git a/tests/test_omega_h.cpp b/tests/test_omega_h.cpp new file mode 100644 index 00000000..516547de --- /dev/null +++ b/tests/test_omega_h.cpp @@ -0,0 +1,324 @@ +// stl includes +#include +#include +#include +#include + +// testing includes +#include +#include +#include + +// xdg includes +#include "util.h" +#include "xdg/constants.h" +#include "xdg/error.h" +#include "xdg/mesh_manager_interface.h" +#include "xdg/omega_h/mesh_manager.h" +#include "xdg/xdg.h" + +using namespace xdg; +using namespace xdg::test; + +// These tests use "cube.osh", an Omega_h binary mesh of a single cubic volume, +// 10 units on a side and centered at the origin (the same geometry used by the +// MOAB "cube.h5m" and libMesh "brick.exo" models). Its classification therefore +// describes one model volume bounded by six planar model surfaces. If the +// converted mesh in your test data differs, update the expected counts and +// ray-fire distances below to match. + +TEST_CASE("Test Omega_h Initialization") { + std::unique_ptr mesh_manager = + std::make_unique(); + + mesh_manager->load_file("cube.osh"); + mesh_manager->init(); + + // a cube is a single classified volume bounded by six surfaces + REQUIRE(mesh_manager->num_volumes() == 1); + REQUIRE(mesh_manager->num_surfaces() == 6); + + // init() already builds the implicit complement, so it should be counted + REQUIRE(mesh_manager->implicit_complement() != ID_NONE); + REQUIRE(mesh_manager->num_volumes() == 2); + + // num_ents_of_dimension should agree with the volume and surface counts + REQUIRE(mesh_manager->num_ents_of_dimension(3) == + mesh_manager->num_volumes()); + REQUIRE(mesh_manager->num_ents_of_dimension(2) == + mesh_manager->num_surfaces()); + + // every surface is bounded by the cube volume on its forward side + MeshID volume = mesh_manager->volumes().front(); + for (auto surface : mesh_manager->get_volume_surfaces(volume)) { + auto senses = mesh_manager->surface_senses(surface); + REQUIRE((senses.first == volume || senses.second == volume)); + } +} + +TEST_CASE("Omega_h Volume Elements") { + std::unique_ptr mesh_manager = + std::make_unique(); + mesh_manager->load_file("cube.osh"); + mesh_manager->init(); + + // unlike the MOAB surface meshes, an Omega_h mesh is volumetric: the cube + // volume must contain tetrahedral elements + MeshID volume = mesh_manager->volumes().front(); + auto elements = mesh_manager->get_volume_elements(volume); + REQUIRE(!elements.empty()); + REQUIRE(mesh_manager->num_volume_elements(volume) == + static_cast(elements.size())); + + // the implicit complement holds no elements + REQUIRE(mesh_manager->num_volume_elements( + mesh_manager->implicit_complement()) == 0); + + // the global element count is the sum over all volumes and equals the number + // of regions (tetrahedra) in the mesh + int total = 0; + for (auto v : mesh_manager->volumes()) { + total += mesh_manager->num_volume_elements(v); + } + REQUIRE(mesh_manager->num_volume_elements() == total); +} + +TEST_CASE("Omega_h Element Types") { + std::shared_ptr xdg = XDG::create(MeshLibrary::OMEGA_H); + REQUIRE(xdg->mesh_manager()->mesh_library() == MeshLibrary::OMEGA_H); + const auto &mesh_manager = xdg->mesh_manager(); + mesh_manager->load_file("cube.osh"); + mesh_manager->init(); + + // Omega_h simplex meshes use triangular surface elements throughout + for (const auto surface : mesh_manager->surfaces()) { + REQUIRE(mesh_manager->get_surface_element_type(surface) == + SurfaceElementType::TRI); + } +} + +TEST_CASE("Omega_h Connectivity") { + std::unique_ptr mesh_manager = + std::make_unique(); + mesh_manager->load_file("cube.osh"); + mesh_manager->init(); + + auto coords_match = [&](MeshID vertex, const Vertex &v) { + auto expected = mesh_manager->vertex_coordinates(vertex); + REQUIRE_THAT(v[0], Catch::Matchers::WithinAbs(expected[0], 1e-12)); + REQUIRE_THAT(v[1], Catch::Matchers::WithinAbs(expected[1], 1e-12)); + REQUIRE_THAT(v[2], Catch::Matchers::WithinAbs(expected[2], 1e-12)); + }; + + // tetrahedra have four vertices; their coordinates must round-trip + MeshID volume = mesh_manager->volumes().front(); + for (auto element : mesh_manager->get_volume_elements(volume)) { + auto conn = mesh_manager->element_connectivity(element); + REQUIRE(conn.size() == 4); + + auto verts = mesh_manager->element_vertices(element); + REQUIRE(verts.size() == conn.size()); + for (std::size_t i = 0; i < conn.size(); ++i) { + coords_match(conn[i], verts[i]); + } + + // a non-degenerate tetrahedron has positive volume + REQUIRE(std::abs(mesh_manager->element_volume(element)) > 0.0); + } + + // triangular faces have three vertices; their coordinates must round-trip + for (auto surface : mesh_manager->surfaces()) { + for (auto face : mesh_manager->get_surface_faces(surface)) { + auto conn = mesh_manager->face_connectivity(face); + REQUIRE(conn.size() == 3); + + auto verts = mesh_manager->face_vertices(face); + for (int i = 0; i < 3; ++i) { + coords_match(conn[i], verts[i]); + } + + // a boundary face is owned by exactly one element of the mesh + REQUIRE(mesh_manager->get_boundary_face_element(face) != ID_NONE); + } + } +} + +TEST_CASE("Omega_h Adjacency") { + std::unique_ptr mesh_manager = + std::make_unique(); + mesh_manager->load_file("cube.osh"); + mesh_manager->init(); + + constexpr int faces_per_tet = 4; + MeshID volume = mesh_manager->volumes().front(); + + for (auto element : mesh_manager->get_volume_elements(volume)) { + for (int face = 0; face < faces_per_tet; ++face) { + MeshID neighbor = mesh_manager->adjacent_element(element, face); + // a neighbor is either another element or a boundary (ID_NONE), but never + // the element itself + REQUIRE(neighbor != element); + if (neighbor != ID_NONE) { + // adjacency is symmetric: the neighbor must list this element back + bool reciprocal = false; + for (int j = 0; j < faces_per_tet; ++j) { + if (mesh_manager->adjacent_element(neighbor, j) == element) { + reciprocal = true; + break; + } + } + REQUIRE(reciprocal); + } + } + } +} + +TEMPLATE_TEST_CASE("Test BVH Build Omega_h", "[omega_h][bvh]", Embree_Raytracer, + GPRT_Raytracer) { + std::shared_ptr mesh_manager = + std::make_shared(); + mesh_manager->load_file("cube.osh"); + mesh_manager->init(); + + REQUIRE(mesh_manager->num_volumes() == 2); + REQUIRE(mesh_manager->num_surfaces() == 6); + + constexpr auto rt_backend = TestType::value; + + DYNAMIC_SECTION(fmt::format("Backend = {}", rt_backend)) { + check_ray_tracer_supported( + rt_backend); // skip if backend not enabled at configuration time + auto rti = create_raytracer(rt_backend); + + for (const auto &volume : mesh_manager->volumes()) { + rti->register_volume(mesh_manager, volume); + } + // at least one acceleration structure is registered per non-empty volume + REQUIRE(rti->num_registered_trees() > 0); + } +} + +TEMPLATE_TEST_CASE("Test Ray Fire Omega_h (all built backends)", + "[ray_tracer][omega_h]", Embree_Raytracer, GPRT_Raytracer) { + constexpr auto rt_backend = TestType::value; + + DYNAMIC_SECTION(fmt::format("Backend = {}", rt_backend)) { + check_ray_tracer_supported( + rt_backend); // skip if backend not enabled at configuration time + auto xdg = XDG::create(MeshLibrary::OMEGA_H, rt_backend); + REQUIRE(xdg->mesh_manager()->mesh_library() == MeshLibrary::OMEGA_H); + + const auto &mm = xdg->mesh_manager(); + mm->load_file("cube.osh"); + mm->init(); + xdg->prepare_raytracer(); + + MeshID volume = mm->volumes()[0]; + + Position origin{0.0, 0.0, 0.0}; + Direction dir{1.0, 0.0, 0.0}; + + // the cube is 10 units on a side and centered at the origin + auto hit = xdg->ray_fire(volume, origin, dir); + REQUIRE(hit.second != ID_NONE); + REQUIRE_THAT(hit.first, Catch::Matchers::WithinAbs(5.0, 1e-6)); + + origin = {3.0, 0.0, 0.0}; + hit = xdg->ray_fire(volume, origin, dir); + REQUIRE_THAT(hit.first, Catch::Matchers::WithinAbs(2.0, 1e-6)); + + origin = {0.0, 0.0, 0.0}; + REQUIRE(xdg->point_in_volume(volume, origin)); + } +} + +TEMPLATE_TEST_CASE("Test Omega_h Find Element Method", "[omega_h][elements]", + Embree_Raytracer) { + constexpr auto rt_backend = TestType::value; + + DYNAMIC_SECTION(fmt::format("Backend = {}", rt_backend)) { + check_ray_tracer_supported( + rt_backend); // skip if backend not enabled at configuration time + std::shared_ptr xdg = XDG::create(MeshLibrary::OMEGA_H, rt_backend); + REQUIRE(xdg->mesh_manager()->mesh_library() == MeshLibrary::OMEGA_H); + const auto &mesh_manager = xdg->mesh_manager(); + mesh_manager->load_file("cube.osh"); + mesh_manager->init(); + xdg->prepare_raytracer(); + + MeshID volume = mesh_manager->volumes().front(); + + MeshID element = xdg->find_element(volume, {0.0, 0.0, 100.0}); + REQUIRE(element == ID_NONE); // point lies outside the cube + + element = xdg->find_element(volume, {0.0, 0.0, 0.0}); + REQUIRE(element != ID_NONE); // point lies inside the cube + + // test the next_element method + auto next_element = + mesh_manager->next_element(element, {0.0, 0.0, 0.0}, {0.0, 0.0, 1.0}); + REQUIRE(next_element.first != ID_NONE); + REQUIRE(next_element.second != INFTY); + + // test the walk_elements method across the full width of the cube + auto walk_elements = mesh_manager->walk_elements(element, {0.0, 0.0, 0.0}, + {0.0, 0.0, 1.0}, 100.0); + double distance = + std::accumulate(walk_elements.begin(), walk_elements.end(), 0.0, + [](double total, const auto &segment) { + return total + segment.second; + }); + REQUIRE(distance > 0.0); + REQUIRE(distance <= 100.0); + for (const auto &segment : walk_elements) { + REQUIRE(segment.first != ID_NONE); + REQUIRE(segment.second >= 0.0); + } + } +} + +TEST_CASE("Omega_h Element ID and Index Mapping") { + std::unique_ptr mesh_manager = + std::make_unique(); + REQUIRE(mesh_manager->mesh_library() == MeshLibrary::OMEGA_H); + mesh_manager->load_file("cube.osh"); + mesh_manager->init(); + + // Omega_h stores entities in a contiguous, zero-based index space, so IDs and + // indices are identical for both elements and vertices + size_t num_elements = mesh_manager->num_volume_elements(); + REQUIRE(num_elements > 0); + for (size_t idx = 0; idx < num_elements; ++idx) { + MeshID element_id = mesh_manager->element_id(idx); + REQUIRE(element_id == static_cast(idx)); + REQUIRE(mesh_manager->element_index(element_id) == static_cast(idx)); + } + + size_t num_vertices = mesh_manager->num_vertices(); + REQUIRE(num_vertices > 0); + for (size_t idx = 0; idx < num_vertices; ++idx) { + MeshID vertex_id = mesh_manager->vertex_id(idx); + REQUIRE(vertex_id == static_cast(idx)); + REQUIRE(mesh_manager->vertex_index(vertex_id) == static_cast(idx)); + } +} + +TEST_CASE("Test Track Exiting Mesh Omega_h") { + std::shared_ptr xdg = XDG::create(MeshLibrary::OMEGA_H); + REQUIRE(xdg->mesh_manager()->mesh_library() == MeshLibrary::OMEGA_H); + const auto &mesh_manager = xdg->mesh_manager(); + mesh_manager->load_file("cube.osh"); + mesh_manager->init(); + xdg->prepare_raytracer(); + + MeshID volume = mesh_manager->volumes().front(); + Position start{0.0, 0.0, -1000.0}; + Position end{0.0, 0.0, 1000.0}; + auto tracks = xdg->segments(volume, start, end); + + // the accumulated track length through the cube equals its 10 unit extent + double length = std::accumulate( + tracks.begin(), tracks.end(), 0.0, + [](double sum, const auto &track) { return sum + track.second; }); + REQUIRE_THAT(length, Catch::Matchers::WithinAbs(10.0, 1e-6)); +} \ No newline at end of file From b805e6ae68a84db96a2447c08021f198cbfabceb Mon Sep 17 00:00:00 2001 From: magnoxemo Date: Sun, 21 Jun 2026 13:56:41 -0500 Subject: [PATCH 10/26] compile omega_h with exodus support --- .github/workflows/ci.yml | 3 ++- .gitignore | 2 ++ Dockerfile | 9 +++++++++ src/omega_h/mesh_manager.cpp | 12 +++++------- tests/test_omega_h.cpp | 11 ++++------- 5 files changed, 22 insertions(+), 15 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 3c001842..f4653c76 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -178,7 +178,8 @@ jobs: -DCMAKE_INSTALL_PREFIX=$HOME/OMEGA_H \ -DCMAKE_BUILD_TYPE=Release \ -DOmega_h_USE_MPI=OFF \ - -DBUILD_TESTING=OFF + -DBUILD_TESTING=OFF \ + -DOmega_h_USE_SEACASExodus=ON make -j4 make install diff --git a/.gitignore b/.gitignore index 17800a19..2b30f634 100644 --- a/.gitignore +++ b/.gitignore @@ -4,3 +4,5 @@ build docs/_build docs/Doxyfile .vscode + +.idea/* \ No newline at end of file diff --git a/Dockerfile b/Dockerfile index 4953a2a5..0208ceef 100644 --- a/Dockerfile +++ b/Dockerfile @@ -49,6 +49,15 @@ RUN make install ENV MOAB_INSTALL_PATH=/XDG_TEST_SYSTEM/moab_install_dir +RUN git clone https://github.com/SCOREC/omega_h.git /XDG_TEST_SYSTEM/omega_h +WORKDIR /XDG_TEST_SYSTEM/omega_h/build +RUN cmake .. \ + -DCMAKE_INSTALL_PREFIX=/XDG_TEST_SYSTEM/omega_h_install_dir \ + -DCMAKE_BUILD_TYPE=Release \ + -DOmega_h_USE_MPI=OFF \ + -DBUILD_TESTING=OFF \ + -Omega_h_USE_SEACASExodus=ON + # build XDG RUN git clone --recurse-submodules https://github.com/xdg-org/xdg.git /XDG_TEST_SYSTEM/xdg WORKDIR /XDG_TEST_SYSTEM/xdg/build diff --git a/src/omega_h/mesh_manager.cpp b/src/omega_h/mesh_manager.cpp index fb697217..e8bf0a2e 100644 --- a/src/omega_h/mesh_manager.cpp +++ b/src/omega_h/mesh_manager.cpp @@ -23,7 +23,7 @@ static constexpr int CLASS_DIM_VOLUME = 3; // Classification tag names written by Omega_h static constexpr const char *CLASS_DIM_TAG = "class_dim"; -static constexpr const char *CLASS_ID_TAG = "class_id"; +static constexpr const char *CLASS_ID_TAG = "class_id"; static constexpr int VERTS_PER_TET = 4; static constexpr int VERTS_PER_TRI = 3; @@ -47,11 +47,11 @@ void OmegaHMeshManager::init() { fatal_error("Mesh must be 3-dimensional"); } - // in Omega_h the regions (3D simplices) are the volume elements + // in Omega_h the regions are 3D simplices or aka elements num_elements_ = mesh_->nregions(); // a classified mesh defines its volumes and surfaces directly through the - // class_dim/class_id tags. Otherwise treat the entire mesh as a single volume + // class_dim/class_id tags. Otherwise, treat the entire mesh as a single volume // bounded by its exposed faces. if (has_classification()) { discover_geometry(); @@ -76,10 +76,8 @@ bool OmegaHMeshManager::has_classification() const { void OmegaHMeshManager::discover_geometry() { // volumes are the unique class IDs of regions classified on a 3D model entity - auto region_class_dim = - mesh_->get_array(OMEGA_H_REGION, CLASS_DIM_TAG); - auto region_class_id = - mesh_->get_array(OMEGA_H_REGION, CLASS_ID_TAG); + auto region_class_dim = mesh_->get_array(OMEGA_H_REGION, CLASS_DIM_TAG); + auto region_class_id = mesh_->get_array(OMEGA_H_REGION, CLASS_ID_TAG); std::set volume_ids; for (Omega_h::LO region = 0; region < mesh_->nregions(); ++region) { diff --git a/tests/test_omega_h.cpp b/tests/test_omega_h.cpp index 516547de..8f3f53a7 100644 --- a/tests/test_omega_h.cpp +++ b/tests/test_omega_h.cpp @@ -28,10 +28,9 @@ using namespace xdg::test; // ray-fire distances below to match. TEST_CASE("Test Omega_h Initialization") { - std::unique_ptr mesh_manager = - std::make_unique(); + std::unique_ptr mesh_manager = std::make_unique(); - mesh_manager->load_file("cube.osh"); + mesh_manager->load_file("brick.exo"); mesh_manager->init(); // a cube is a single classified volume bounded by six surfaces @@ -43,10 +42,8 @@ TEST_CASE("Test Omega_h Initialization") { REQUIRE(mesh_manager->num_volumes() == 2); // num_ents_of_dimension should agree with the volume and surface counts - REQUIRE(mesh_manager->num_ents_of_dimension(3) == - mesh_manager->num_volumes()); - REQUIRE(mesh_manager->num_ents_of_dimension(2) == - mesh_manager->num_surfaces()); + REQUIRE(mesh_manager->num_ents_of_dimension(3) == mesh_manager->num_volumes()); + REQUIRE(mesh_manager->num_ents_of_dimension(2) == mesh_manager->num_surfaces()); // every surface is bounded by the cube volume on its forward side MeshID volume = mesh_manager->volumes().front(); From b92d36275553e0406599d155dbecc2fc15673c6f Mon Sep 17 00:00:00 2001 From: magnoxemo Date: Sun, 21 Jun 2026 14:02:29 -0500 Subject: [PATCH 11/26] add omega_h to the docker file --- Dockerfile | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/Dockerfile b/Dockerfile index 0208ceef..d57b1f58 100644 --- a/Dockerfile +++ b/Dockerfile @@ -56,7 +56,13 @@ RUN cmake .. \ -DCMAKE_BUILD_TYPE=Release \ -DOmega_h_USE_MPI=OFF \ -DBUILD_TESTING=OFF \ - -Omega_h_USE_SEACASExodus=ON + -Omega_h_USE_SEACASExodus=ON \ + +RUN make -j $(BUILD_JOBS) +RUN make install +ENV OMEGA_H_INSTALL_PATH=/XDG_TEST_SYSTEM/omega_h_install_dir + + # build XDG RUN git clone --recurse-submodules https://github.com/xdg-org/xdg.git /XDG_TEST_SYSTEM/xdg @@ -69,7 +75,10 @@ RUN cmake .. \ -DXDG_ENABLE_MOAB=ON \ -DMOAB_DIR=${MOAB_INSTALL_PATH} \ -DXDG_ENABLE_LIBMESH=ON \ - -DLIBMESH_DIR=${LIBMESH_INSTALL_PATH} + -DLIBMESH_DIR=${LIBMESH_INSTALL_PATH} \ + -DXDG_ENABLE_OMEGA_H=ON \ + -DOMEGA_H_dir=${OMEGA_H_INSTALL_PATH} \ + RUN make -j $(BUILD_JOBS) RUN make install From 2b855a8d5726f86b295df108f4594245129e3c8d Mon Sep 17 00:00:00 2001 From: magnoxemo Date: Wed, 24 Jun 2026 10:41:43 -0500 Subject: [PATCH 12/26] fix testing mesh file name --- include/xdg/omega_h/mesh_manager.h | 5 ++--- src/omega_h/mesh_manager.cpp | 8 +++++++- tests/test_omega_h.cpp | 25 +++++++++---------------- 3 files changed, 18 insertions(+), 20 deletions(-) diff --git a/include/xdg/omega_h/mesh_manager.h b/include/xdg/omega_h/mesh_manager.h index 5e2185ca..1d8dd919 100644 --- a/include/xdg/omega_h/mesh_manager.h +++ b/include/xdg/omega_h/mesh_manager.h @@ -93,11 +93,10 @@ class OmegaHMeshManager : public MeshManager { std::array face_vertices(MeshID face) const override; //! \brief Vertices of a local face of a tetrahedral element - std::array element_face_vertices(MeshID element, - int local_face) const; + std::array element_face_vertices(MeshID element, int local_face) const; SurfaceElementType - get_surface_element_type(MeshID /*surface*/) const override { + get_surface_element_type(MeshID surface_element_id) const override { // Omega_h simplex meshes always use triangular surface elements return SurfaceElementType::TRI; } diff --git a/src/omega_h/mesh_manager.cpp b/src/omega_h/mesh_manager.cpp index e8bf0a2e..7d8dbd9c 100644 --- a/src/omega_h/mesh_manager.cpp +++ b/src/omega_h/mesh_manager.cpp @@ -31,7 +31,13 @@ static constexpr int FACES_PER_TET = 4; // Constructor -OmegaHMeshManager::OmegaHMeshManager(const Omega_h::Mesh* mesh) {} +OmegaHMeshManager::OmegaHMeshManager(const Omega_h::Mesh* mesh) { + + // I will rewrite this method later + library_ = std::make_unique(nullptr, nullptr); + mesh_ = std::make_unique(library_.get()); + +} OmegaHMeshManager::OmegaHMeshManager() { library_ = std::make_unique(nullptr, nullptr); mesh_ = std::make_unique(library_.get()); diff --git a/tests/test_omega_h.cpp b/tests/test_omega_h.cpp index 8f3f53a7..fa9b40e9 100644 --- a/tests/test_omega_h.cpp +++ b/tests/test_omega_h.cpp @@ -1,15 +1,12 @@ -// stl includes #include #include #include #include -// testing includes #include #include #include -// xdg includes #include "util.h" #include "xdg/constants.h" #include "xdg/error.h" @@ -20,10 +17,8 @@ using namespace xdg; using namespace xdg::test; -// These tests use "cube.osh", an Omega_h binary mesh of a single cubic volume, -// 10 units on a side and centered at the origin (the same geometry used by the -// MOAB "cube.h5m" and libMesh "brick.exo" models). Its classification therefore -// describes one model volume bounded by six planar model surfaces. If the +// These tests use "brick.exo". +// Its classification therefore describes one model volume bounded by six planar model surfaces. If the // converted mesh in your test data differs, update the expected counts and // ray-fire distances below to match. @@ -56,7 +51,7 @@ TEST_CASE("Test Omega_h Initialization") { TEST_CASE("Omega_h Volume Elements") { std::unique_ptr mesh_manager = std::make_unique(); - mesh_manager->load_file("cube.osh"); + mesh_manager->load_file("brick.exo"); mesh_manager->init(); // unlike the MOAB surface meshes, an Omega_h mesh is volumetric: the cube @@ -84,7 +79,7 @@ TEST_CASE("Omega_h Element Types") { std::shared_ptr xdg = XDG::create(MeshLibrary::OMEGA_H); REQUIRE(xdg->mesh_manager()->mesh_library() == MeshLibrary::OMEGA_H); const auto &mesh_manager = xdg->mesh_manager(); - mesh_manager->load_file("cube.osh"); + mesh_manager->load_file("brick.exo"); mesh_manager->init(); // Omega_h simplex meshes use triangular surface elements throughout @@ -97,7 +92,7 @@ TEST_CASE("Omega_h Element Types") { TEST_CASE("Omega_h Connectivity") { std::unique_ptr mesh_manager = std::make_unique(); - mesh_manager->load_file("cube.osh"); + mesh_manager->load_file("brick.exo"); mesh_manager->init(); auto coords_match = [&](MeshID vertex, const Vertex &v) { @@ -141,9 +136,8 @@ TEST_CASE("Omega_h Connectivity") { } TEST_CASE("Omega_h Adjacency") { - std::unique_ptr mesh_manager = - std::make_unique(); - mesh_manager->load_file("cube.osh"); + std::unique_ptr mesh_manager = std::make_unique(); + mesh_manager->load_file("brick.exo"); mesh_manager->init(); constexpr int faces_per_tet = 4; @@ -174,7 +168,7 @@ TEMPLATE_TEST_CASE("Test BVH Build Omega_h", "[omega_h][bvh]", Embree_Raytracer, GPRT_Raytracer) { std::shared_ptr mesh_manager = std::make_shared(); - mesh_manager->load_file("cube.osh"); + mesh_manager->load_file("brick.exo"); mesh_manager->init(); REQUIRE(mesh_manager->num_volumes() == 2); @@ -275,8 +269,7 @@ TEMPLATE_TEST_CASE("Test Omega_h Find Element Method", "[omega_h][elements]", } TEST_CASE("Omega_h Element ID and Index Mapping") { - std::unique_ptr mesh_manager = - std::make_unique(); + std::unique_ptr mesh_manager = std::make_unique(); REQUIRE(mesh_manager->mesh_library() == MeshLibrary::OMEGA_H); mesh_manager->load_file("cube.osh"); mesh_manager->init(); From e86a8a266540a938e0972e8059a2053a238ad0a6 Mon Sep 17 00:00:00 2001 From: magnoxemo Date: Wed, 24 Jun 2026 11:35:45 -0500 Subject: [PATCH 13/26] update mesh file name --- tests/test_omega_h.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/test_omega_h.cpp b/tests/test_omega_h.cpp index fa9b40e9..58446a2b 100644 --- a/tests/test_omega_h.cpp +++ b/tests/test_omega_h.cpp @@ -200,7 +200,7 @@ TEMPLATE_TEST_CASE("Test Ray Fire Omega_h (all built backends)", REQUIRE(xdg->mesh_manager()->mesh_library() == MeshLibrary::OMEGA_H); const auto &mm = xdg->mesh_manager(); - mm->load_file("cube.osh"); + mm->load_file("brick.exo"); mm->init(); xdg->prepare_raytracer(); @@ -233,7 +233,7 @@ TEMPLATE_TEST_CASE("Test Omega_h Find Element Method", "[omega_h][elements]", std::shared_ptr xdg = XDG::create(MeshLibrary::OMEGA_H, rt_backend); REQUIRE(xdg->mesh_manager()->mesh_library() == MeshLibrary::OMEGA_H); const auto &mesh_manager = xdg->mesh_manager(); - mesh_manager->load_file("cube.osh"); + mesh_manager->load_file("brick.exo"); mesh_manager->init(); xdg->prepare_raytracer(); @@ -271,7 +271,7 @@ TEMPLATE_TEST_CASE("Test Omega_h Find Element Method", "[omega_h][elements]", TEST_CASE("Omega_h Element ID and Index Mapping") { std::unique_ptr mesh_manager = std::make_unique(); REQUIRE(mesh_manager->mesh_library() == MeshLibrary::OMEGA_H); - mesh_manager->load_file("cube.osh"); + mesh_manager->load_file("brick.exo"); mesh_manager->init(); // Omega_h stores entities in a contiguous, zero-based index space, so IDs and @@ -297,7 +297,7 @@ TEST_CASE("Test Track Exiting Mesh Omega_h") { std::shared_ptr xdg = XDG::create(MeshLibrary::OMEGA_H); REQUIRE(xdg->mesh_manager()->mesh_library() == MeshLibrary::OMEGA_H); const auto &mesh_manager = xdg->mesh_manager(); - mesh_manager->load_file("cube.osh"); + mesh_manager->load_file("brick.exo"); mesh_manager->init(); xdg->prepare_raytracer(); From 2bfa6ccce357e61773d5919225c67f10de0e22e1 Mon Sep 17 00:00:00 2001 From: magnoxemo Date: Thu, 25 Jun 2026 12:35:56 -0500 Subject: [PATCH 14/26] install seacas in docker file to support exodus mesh in omega_h --- Dockerfile | 34 +++++++++++++++++++++++++++------- 1 file changed, 27 insertions(+), 7 deletions(-) diff --git a/Dockerfile b/Dockerfile index d57b1f58..a0c3a47e 100644 --- a/Dockerfile +++ b/Dockerfile @@ -23,13 +23,14 @@ RUN apt-get install --yes \ curl \ nano \ libembree-dev \ - pkg-config + pkg-config \ + libnetcdf-dev # compile libmesh from source and install it RUN git clone --recurse-submodules https://github.com/libMesh/libmesh.git /XDG_TEST_SYSTEM/libmesh WORKDIR /XDG_TEST_SYSTEM/libmesh/build RUN ../configure --prefix=/XDG_TEST_SYSTEM/libmesh_install_dir -RUN make -j $(BUILD_JOBS) +RUN make -j ${BUILD_JOBS} RUN make install ENV LIBMESH_INSTALL_PATH=/XDG_TEST_SYSTEM/libmesh_install_dir @@ -43,12 +44,31 @@ RUN cmake .. \ -DHDF5_ROOT=/usr \ -DBLAS_LIBRARIES=/usr/lib/x86_64-linux-gnu/libopenblas.so \ -DBUILD_SHARED_LIBS=ON -RUN make -j $(BUILD_JOBS) +RUN make -j${BUILD_JOBS} RUN make install ENV MOAB_INSTALL_PATH=/XDG_TEST_SYSTEM/moab_install_dir +# installing SEACAS to enable exodus support +RUN git clone --depth=1 https://github.com/sandialabs/seacas.git /XDG_TEST_SYSTEM/seacas +WORKDIR /XDG_TEST_SYSTEM/seacas/build + +RUN cmake .. \ + -DCMAKE_INSTALL_PREFIX=/XDG_TEST_SYSTEM/seacas_install_dir \ + -DCMAKE_BUILD_TYPE=Release \ + -DSeacas_ENABLE_ALL_PACKAGES=OFF \ + -DSeacas_ENABLE_SEACASExodus=ON \ + -DTPL_ENABLE_MPI=OFF \ + -DTPL_ENABLE_Netcdf=ON \ + -DTPL_Netcdf_INCLUDE_DIRS=/usr/include \ + -DTPL_Netcdf_LIBRARIES=/usr/lib/x86_64-linux-gnu/libnetcdf.so + +RUN make -j${BUILD_JOBS} +RUN make install +ENV SEACAS_INSTALL_PATH=/XDG_TEST_SYSTEM/seacas_install_dir + +# build omega_h RUN git clone https://github.com/SCOREC/omega_h.git /XDG_TEST_SYSTEM/omega_h WORKDIR /XDG_TEST_SYSTEM/omega_h/build RUN cmake .. \ @@ -56,9 +76,9 @@ RUN cmake .. \ -DCMAKE_BUILD_TYPE=Release \ -DOmega_h_USE_MPI=OFF \ -DBUILD_TESTING=OFF \ - -Omega_h_USE_SEACASExodus=ON \ + -DOmega_h_USE_SEACASExodus=ON -RUN make -j $(BUILD_JOBS) +RUN make -j ${BUILD_JOBS} RUN make install ENV OMEGA_H_INSTALL_PATH=/XDG_TEST_SYSTEM/omega_h_install_dir @@ -77,9 +97,9 @@ RUN cmake .. \ -DXDG_ENABLE_LIBMESH=ON \ -DLIBMESH_DIR=${LIBMESH_INSTALL_PATH} \ -DXDG_ENABLE_OMEGA_H=ON \ - -DOMEGA_H_dir=${OMEGA_H_INSTALL_PATH} \ + -DOMEGA_H_dir=${OMEGA_H_INSTALL_PATH} -RUN make -j $(BUILD_JOBS) +RUN make -j${BUILD_JOBS} RUN make install ENV XDG_INSTALL_PATH=/XDG_TEST_SYSTEM/xdg_install_dir From 929b883843496cd93ead47163b1b4dc3e2f423d9 Mon Sep 17 00:00:00 2001 From: magnoxemo Date: Thu, 25 Jun 2026 13:47:32 -0500 Subject: [PATCH 15/26] setting xdg branch in docker file to this branch for now --- Dockerfile | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/Dockerfile b/Dockerfile index a0c3a47e..a66c1523 100644 --- a/Dockerfile +++ b/Dockerfile @@ -53,16 +53,14 @@ ENV MOAB_INSTALL_PATH=/XDG_TEST_SYSTEM/moab_install_dir RUN git clone --depth=1 https://github.com/sandialabs/seacas.git /XDG_TEST_SYSTEM/seacas WORKDIR /XDG_TEST_SYSTEM/seacas/build - RUN cmake .. \ -DCMAKE_INSTALL_PREFIX=/XDG_TEST_SYSTEM/seacas_install_dir \ -DCMAKE_BUILD_TYPE=Release \ + -DBUILD_SHARED_LIBS=ON \ -DSeacas_ENABLE_ALL_PACKAGES=OFF \ -DSeacas_ENABLE_SEACASExodus=ON \ -DTPL_ENABLE_MPI=OFF \ - -DTPL_ENABLE_Netcdf=ON \ - -DTPL_Netcdf_INCLUDE_DIRS=/usr/include \ - -DTPL_Netcdf_LIBRARIES=/usr/lib/x86_64-linux-gnu/libnetcdf.so + -DTPL_ENABLE_Netcdf=ON RUN make -j${BUILD_JOBS} RUN make install @@ -76,7 +74,8 @@ RUN cmake .. \ -DCMAKE_BUILD_TYPE=Release \ -DOmega_h_USE_MPI=OFF \ -DBUILD_TESTING=OFF \ - -DOmega_h_USE_SEACASExodus=ON + -DSEACASExodus_DIR=${SEACAS_INSTALL_PATH}/lib/cmake/SEACASExodus\ + -DOmega_h_USE_SEACASExodus=ON RUN make -j ${BUILD_JOBS} RUN make install @@ -85,19 +84,20 @@ ENV OMEGA_H_INSTALL_PATH=/XDG_TEST_SYSTEM/omega_h_install_dir # build XDG -RUN git clone --recurse-submodules https://github.com/xdg-org/xdg.git /XDG_TEST_SYSTEM/xdg +# currently building docker image from this branch. As omega_h isn't include the cmake build system +# in upstream yet. I will change this later. +RUN git clone -b omega_h --recurse-submodules https://github.com/magnoxemo/xdg.git /XDG_TEST_SYSTEM/xdg + WORKDIR /XDG_TEST_SYSTEM/xdg/build RUN cmake .. \ - -DCMAKE_INSTALL_PREFIX=/XDG_TEST_SYSTEM/xdg_install_dir \ - -DCMAKE_BUILD_TYPE=Release \ - -DCMAKE_C_COMPILER=mpicc \ - -DCMAKE_CXX_COMPILER=mpicxx \ - -DXDG_ENABLE_MOAB=ON \ - -DMOAB_DIR=${MOAB_INSTALL_PATH} \ - -DXDG_ENABLE_LIBMESH=ON \ - -DLIBMESH_DIR=${LIBMESH_INSTALL_PATH} \ - -DXDG_ENABLE_OMEGA_H=ON \ - -DOMEGA_H_dir=${OMEGA_H_INSTALL_PATH} + -DCMAKE_INSTALL_PREFIX=/XDG_TEST_SYSTEM/xdg_install_dir \ + -DCMAKE_BUILD_TYPE=Release \ + -DCMAKE_C_COMPILER=mpicc \ + -DCMAKE_CXX_COMPILER=mpicxx \ + -DXDG_ENABLE_OMEGA_H=ON \ + -DXDG_ENABLE_MOAB=ON \ + -DXDG_ENABLE_LIBMESH=ON \ + -DCMAKE_PREFIX_PATH="${OMEGA_H_INSTALL_PATH};${SEACAS_INSTALL_PATH};${MOAB_INSTALL_PATH};${LIBMESH_INSTALL_PATH}" RUN make -j${BUILD_JOBS} RUN make install From 7e8a4797a0d8120688d30cc5a43af0bc0725018f Mon Sep 17 00:00:00 2001 From: magnoxemo Date: Sat, 27 Jun 2026 18:07:48 -0500 Subject: [PATCH 16/26] update main CI to include seacas --- .github/workflows/ci.yml | 20 ++++++++++++++++++++ Dockerfile | 2 +- 2 files changed, 21 insertions(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index f4653c76..9e3f61e9 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -91,6 +91,9 @@ jobs: git clone https://github.com/SCOREC/omega_h.git cd omega_h git checkout master + cd ~ + git clone https://github.com/sandialabs/seacas.git + - name: OpenMP Environment Variables run: | @@ -171,6 +174,22 @@ jobs: name: Build Omega_h shell: bash run: | + # before building omega_h we need to compile seacas on which omega_h depends + # for interacting with exodus mesh API + cd ~/seacas + mkdir build + cd build + cmake .. \ + -DCMAKE_INSTALL_PREFIX=$HOME/SEACAS \ + -DCMAKE_BUILD_TYPE=Release \ + -DBUILD_SHARED_LIBS=ON \ + -DSeacas_ENABLE_ALL_PACKAGES=OFF \ + -DSeacas_ENABLE_SEACASExodus=ON \ + -DTPL_ENABLE_MPI=OFF \ + -DTPL_ENABLE_Netcdf=ON + make -j4 + make install + # now let's install omega_h cd ~/omega_h mkdir build cd build @@ -179,6 +198,7 @@ jobs: -DCMAKE_BUILD_TYPE=Release \ -DOmega_h_USE_MPI=OFF \ -DBUILD_TESTING=OFF \ + -DSEACASExodus_DIR=$HOME/SEACAS/lib/cmake/SEACASExodus \ -DOmega_h_USE_SEACASExodus=ON make -j4 make install diff --git a/Dockerfile b/Dockerfile index a66c1523..0d115533 100644 --- a/Dockerfile +++ b/Dockerfile @@ -50,7 +50,7 @@ RUN make install ENV MOAB_INSTALL_PATH=/XDG_TEST_SYSTEM/moab_install_dir # installing SEACAS to enable exodus support -RUN git clone --depth=1 https://github.com/sandialabs/seacas.git /XDG_TEST_SYSTEM/seacas +RUN git clone https://github.com/sandialabs/seacas.git /XDG_TEST_SYSTEM/seacas WORKDIR /XDG_TEST_SYSTEM/seacas/build RUN cmake .. \ From be2d8492de59884a8a6d586e10f0566a2ab9c5b6 Mon Sep 17 00:00:00 2001 From: magnoxemo Date: Sat, 27 Jun 2026 18:36:51 -0500 Subject: [PATCH 17/26] update xdg headers and lib pointers to allows using omega_h there is problem with exodus lib version which I will address later --- include/xdg/mesh_managers.h | 4 ++++ src/omega_h/mesh_manager.cpp | 8 ++++++-- src/xdg.cpp | 6 ++++++ 3 files changed, 16 insertions(+), 2 deletions(-) diff --git a/include/xdg/mesh_managers.h b/include/xdg/mesh_managers.h index fe6912ce..8c3a5c82 100644 --- a/include/xdg/mesh_managers.h +++ b/include/xdg/mesh_managers.h @@ -6,3 +6,7 @@ #ifdef XDG_ENABLE_LIBMESH #include "xdg/libmesh/mesh_manager.h" #endif + +#ifdef XDG_ENABLE_OMEGA_H +#include "xdg/omega_h/mesh_manager.h" +#endif \ No newline at end of file diff --git a/src/omega_h/mesh_manager.cpp b/src/omega_h/mesh_manager.cpp index 7d8dbd9c..97c7640c 100644 --- a/src/omega_h/mesh_manager.cpp +++ b/src/omega_h/mesh_manager.cpp @@ -43,8 +43,12 @@ OmegaHMeshManager::OmegaHMeshManager() { mesh_ = std::make_unique(library_.get()); } -void OmegaHMeshManager::load_file(const std::string &file_path) { - *mesh_ = Omega_h::binary::read(file_path, library_.get()); + +void OmegaHMeshManager::load_file(const std::string &file_path) +{ + const int exodus_file = Omega_h::exodus::open(file_path); + Omega_h::exodus::read_mesh(exodus_file, mesh_.get()); + Omega_h::exodus::close(exodus_file); } void OmegaHMeshManager::init() { diff --git a/src/xdg.cpp b/src/xdg.cpp index cc738008..3900ac47 100644 --- a/src/xdg.cpp +++ b/src/xdg.cpp @@ -64,6 +64,9 @@ std::shared_ptr XDG::create(MeshLibrary mesh_lib, RTLibrary ray_tracing_lib #ifdef XDG_ENABLE_LIBMESH if (mesh_lib == MeshLibrary::LIBMESH) return std::make_shared(); #endif + #ifdef XDG_ENABLE_OMEGA_H + if (mesh_lib == MeshLibrary::OMEGA_H) return std::make_shared(); + #endif // If no supported mesh library throw an error std::string msg = fmt::format("Invalid mesh library '{}'. Supported:", MESH_LIB_TO_STR.at(mesh_lib)); @@ -73,6 +76,9 @@ std::shared_ptr XDG::create(MeshLibrary mesh_lib, RTLibrary ray_tracing_lib #ifdef XDG_ENABLE_LIBMESH msg += " LIBMESH"; #endif + #ifdef XDG_ENABLE_OMEGA_H + msg += " OMEGA_H"; + #endif fatal_error(msg); }; From 5e003ca797f5552fab1c9ace6d7af4e8ade43710 Mon Sep 17 00:00:00 2001 From: magnoxemo Date: Mon, 13 Jul 2026 10:01:23 -0500 Subject: [PATCH 18/26] change mesh file against which we will be testing I will have to change the gold params for this mesh file as well --- include/xdg/constants.h | 2 +- src/xdg.cpp | 2 +- tests/test_omega_h.cpp | 79 ++++++++++++++++++++++------------------- tests/util.h | 20 ++++++++--- 4 files changed, 59 insertions(+), 44 deletions(-) diff --git a/include/xdg/constants.h b/include/xdg/constants.h index a668e8f3..4c01d6a8 100644 --- a/include/xdg/constants.h +++ b/include/xdg/constants.h @@ -72,7 +72,7 @@ static const std::map RT_LIB_TO_STR = {RTLibrary::GPRT, "GPRT"} }; -// Mesh identifer type +// Mesh identifier type using MeshID = int32_t; using MeshIndex = int32_t; diff --git a/src/xdg.cpp b/src/xdg.cpp index 3900ac47..ea35c2f2 100644 --- a/src/xdg.cpp +++ b/src/xdg.cpp @@ -77,7 +77,7 @@ std::shared_ptr XDG::create(MeshLibrary mesh_lib, RTLibrary ray_tracing_lib msg += " LIBMESH"; #endif #ifdef XDG_ENABLE_OMEGA_H - msg += " OMEGA_H"; + msg += " OMEGA_H"; #endif fatal_error(msg); }; diff --git a/tests/test_omega_h.cpp b/tests/test_omega_h.cpp index 58446a2b..eaaa2e11 100644 --- a/tests/test_omega_h.cpp +++ b/tests/test_omega_h.cpp @@ -8,7 +8,6 @@ #include #include "util.h" -#include "xdg/constants.h" #include "xdg/error.h" #include "xdg/mesh_manager_interface.h" #include "xdg/omega_h/mesh_manager.h" @@ -18,27 +17,28 @@ using namespace xdg; using namespace xdg::test; // These tests use "brick.exo". -// Its classification therefore describes one model volume bounded by six planar model surfaces. If the -// converted mesh in your test data differs, update the expected counts and -// ray-fire distances below to match. +// Its classification therefore describes one model volume bounded by six planar +// model surfaces. If the converted mesh in your test data differs, update the +// expected counts and ray-fire distances below to match. TEST_CASE("Test Omega_h Initialization") { - std::unique_ptr mesh_manager = std::make_unique(); + std::unique_ptr mesh_manager = + std::make_unique(); - mesh_manager->load_file("brick.exo"); + mesh_manager->load_file("pincell-implicit.exo"); mesh_manager->init(); - // a cube is a single classified volume bounded by six surfaces - REQUIRE(mesh_manager->num_volumes() == 1); - REQUIRE(mesh_manager->num_surfaces() == 6); + REQUIRE(mesh_manager->num_volumes() == 5); + REQUIRE(mesh_manager->num_surfaces() == 13); // init() already builds the implicit complement, so it should be counted REQUIRE(mesh_manager->implicit_complement() != ID_NONE); - REQUIRE(mesh_manager->num_volumes() == 2); // num_ents_of_dimension should agree with the volume and surface counts - REQUIRE(mesh_manager->num_ents_of_dimension(3) == mesh_manager->num_volumes()); - REQUIRE(mesh_manager->num_ents_of_dimension(2) == mesh_manager->num_surfaces()); + REQUIRE(mesh_manager->num_ents_of_dimension(3) == + mesh_manager->num_volumes()); + REQUIRE(mesh_manager->num_ents_of_dimension(2) == + mesh_manager->num_surfaces()); // every surface is bounded by the cube volume on its forward side MeshID volume = mesh_manager->volumes().front(); @@ -51,7 +51,7 @@ TEST_CASE("Test Omega_h Initialization") { TEST_CASE("Omega_h Volume Elements") { std::unique_ptr mesh_manager = std::make_unique(); - mesh_manager->load_file("brick.exo"); + mesh_manager->load_file("pincell-implicit.exo"); mesh_manager->init(); // unlike the MOAB surface meshes, an Omega_h mesh is volumetric: the cube @@ -79,7 +79,7 @@ TEST_CASE("Omega_h Element Types") { std::shared_ptr xdg = XDG::create(MeshLibrary::OMEGA_H); REQUIRE(xdg->mesh_manager()->mesh_library() == MeshLibrary::OMEGA_H); const auto &mesh_manager = xdg->mesh_manager(); - mesh_manager->load_file("brick.exo"); + mesh_manager->load_file("pincell-implicit.exo"); mesh_manager->init(); // Omega_h simplex meshes use triangular surface elements throughout @@ -92,7 +92,7 @@ TEST_CASE("Omega_h Element Types") { TEST_CASE("Omega_h Connectivity") { std::unique_ptr mesh_manager = std::make_unique(); - mesh_manager->load_file("brick.exo"); + mesh_manager->load_file("pincell-implicit.exo"); mesh_manager->init(); auto coords_match = [&](MeshID vertex, const Vertex &v) { @@ -136,11 +136,12 @@ TEST_CASE("Omega_h Connectivity") { } TEST_CASE("Omega_h Adjacency") { - std::unique_ptr mesh_manager = std::make_unique(); - mesh_manager->load_file("brick.exo"); + std::unique_ptr mesh_manager = + std::make_unique(); + mesh_manager->load_file("pincell-implicit.exo"); mesh_manager->init(); - constexpr int faces_per_tet = 4; + constexpr unsigned int faces_per_tet = 4; MeshID volume = mesh_manager->volumes().front(); for (auto element : mesh_manager->get_volume_elements(volume)) { @@ -164,17 +165,17 @@ TEST_CASE("Omega_h Adjacency") { } } -TEMPLATE_TEST_CASE("Test BVH Build Omega_h", "[omega_h][bvh]", Embree_Raytracer, - GPRT_Raytracer) { +TEMPLATE_TEST_CASE("Test BVH Build Omega_h", "[omega_h][bvh]", + Embree_Raytracer) { + constexpr auto rt_backend = TestType::value; + std::shared_ptr mesh_manager = std::make_shared(); - mesh_manager->load_file("brick.exo"); + mesh_manager->load_file("pincell-implicit.exo"); mesh_manager->init(); - REQUIRE(mesh_manager->num_volumes() == 2); - REQUIRE(mesh_manager->num_surfaces() == 6); - - constexpr auto rt_backend = TestType::value; + REQUIRE(mesh_manager->num_volumes() == 5); + REQUIRE(mesh_manager->num_surfaces() == 13); DYNAMIC_SECTION(fmt::format("Backend = {}", rt_backend)) { check_ray_tracer_supported( @@ -184,8 +185,7 @@ TEMPLATE_TEST_CASE("Test BVH Build Omega_h", "[omega_h][bvh]", Embree_Raytracer, for (const auto &volume : mesh_manager->volumes()) { rti->register_volume(mesh_manager, volume); } - // at least one acceleration structure is registered per non-empty volume - REQUIRE(rti->num_registered_trees() > 0); + REQUIRE(rti->num_registered_trees() == 2); } } @@ -200,10 +200,10 @@ TEMPLATE_TEST_CASE("Test Ray Fire Omega_h (all built backends)", REQUIRE(xdg->mesh_manager()->mesh_library() == MeshLibrary::OMEGA_H); const auto &mm = xdg->mesh_manager(); - mm->load_file("brick.exo"); + mm->load_file("pincell-implicit.exo"); mm->init(); - xdg->prepare_raytracer(); + xdg->prepare_raytracer(); MeshID volume = mm->volumes()[0]; Position origin{0.0, 0.0, 0.0}; @@ -225,6 +225,7 @@ TEMPLATE_TEST_CASE("Test Ray Fire Omega_h (all built backends)", TEMPLATE_TEST_CASE("Test Omega_h Find Element Method", "[omega_h][elements]", Embree_Raytracer) { + // Gold values for this test needs to be fixed. I will do that later. constexpr auto rt_backend = TestType::value; DYNAMIC_SECTION(fmt::format("Backend = {}", rt_backend)) { @@ -233,14 +234,15 @@ TEMPLATE_TEST_CASE("Test Omega_h Find Element Method", "[omega_h][elements]", std::shared_ptr xdg = XDG::create(MeshLibrary::OMEGA_H, rt_backend); REQUIRE(xdg->mesh_manager()->mesh_library() == MeshLibrary::OMEGA_H); const auto &mesh_manager = xdg->mesh_manager(); - mesh_manager->load_file("brick.exo"); + mesh_manager->load_file("pincell-implicit.exo"); mesh_manager->init(); xdg->prepare_raytracer(); MeshID volume = mesh_manager->volumes().front(); - MeshID element = xdg->find_element(volume, {0.0, 0.0, 100.0}); - REQUIRE(element == ID_NONE); // point lies outside the cube + MeshID element = xdg->find_element( + volume, {0.0, 0.0, 100.0}); // I will need to fix this one as well + REQUIRE(element == ID_NONE); // point lies outside the cube element = xdg->find_element(volume, {0.0, 0.0, 0.0}); REQUIRE(element != ID_NONE); // point lies inside the cube @@ -269,14 +271,17 @@ TEMPLATE_TEST_CASE("Test Omega_h Find Element Method", "[omega_h][elements]", } TEST_CASE("Omega_h Element ID and Index Mapping") { - std::unique_ptr mesh_manager = std::make_unique(); + // Gold values for this test needs to be fixed. I will do that later. + std::unique_ptr mesh_manager = + std::make_unique(); REQUIRE(mesh_manager->mesh_library() == MeshLibrary::OMEGA_H); mesh_manager->load_file("brick.exo"); mesh_manager->init(); - // Omega_h stores entities in a contiguous, zero-based index space, so IDs and - // indices are identical for both elements and vertices - size_t num_elements = mesh_manager->num_volume_elements(); + // Omega_h stores entities in a contiguous, zero-based index space, so IDs + and + // indices are identical for both elements and vertices + size_t num_elements = mesh_manager->num_volume_elements(); REQUIRE(num_elements > 0); for (size_t idx = 0; idx < num_elements; ++idx) { MeshID element_id = mesh_manager->element_id(idx); @@ -297,7 +302,7 @@ TEST_CASE("Test Track Exiting Mesh Omega_h") { std::shared_ptr xdg = XDG::create(MeshLibrary::OMEGA_H); REQUIRE(xdg->mesh_manager()->mesh_library() == MeshLibrary::OMEGA_H); const auto &mesh_manager = xdg->mesh_manager(); - mesh_manager->load_file("brick.exo"); + mesh_manager->load_file("pincell-implicit.exo"); mesh_manager->init(); xdg->prepare_raytracer(); diff --git a/tests/util.h b/tests/util.h index a9841086..b6fbbc39 100644 --- a/tests/util.h +++ b/tests/util.h @@ -67,11 +67,17 @@ inline bool mesh_library_available(xdg::MeshLibrary mesh) { #endif case xdg::MeshLibrary::LIBMESH: - #ifdef XDG_ENABLE_LIBMESH - return true; - #else - return false; - #endif + #ifdef XDG_ENABLE_LIBMESH + return true; + #else + return false; + #endif + case xdg::MeshLibrary::OMEGA_H: + #ifdef XDG_ENABLE_OMEGA_H + return true; + #else + return false; + #endif } return false; @@ -102,6 +108,10 @@ create_mesh_manager(xdg::MeshLibrary mesh) { if (mesh == xdg::MeshLibrary::LIBMESH) return std::make_unique(); #endif + #ifdef XDG_ENABLE_OMEGA_H + if (mesh == xdg::MeshLibrary::OMEGA_H) + return std::make_unique(); + #endif return nullptr; } From 05072883e7fbff0b9a1298937f8ce38bf2000b80 Mon Sep 17 00:00:00 2001 From: magnoxemo Date: Thu, 23 Jul 2026 11:10:31 -0500 Subject: [PATCH 19/26] update test mesh I need to fix the rti->register_volume() method --- include/xdg/omega_h/mesh_manager.h | 3 ++- src/embree/ray_tracer.cpp | 1 + src/omega_h/mesh_manager.cpp | 13 +++++---- tests/test_omega_h.cpp | 42 ++++++++++++++---------------- 4 files changed, 29 insertions(+), 30 deletions(-) diff --git a/include/xdg/omega_h/mesh_manager.h b/include/xdg/omega_h/mesh_manager.h index 1d8dd919..4d0f5a7d 100644 --- a/include/xdg/omega_h/mesh_manager.h +++ b/include/xdg/omega_h/mesh_manager.h @@ -39,7 +39,8 @@ class OmegaHMeshManager : public MeshManager { void init() override; - void parse_metadata() override {} + // TODO: I will have to implement this one. + void parse_metadata() override {}; int num_volumes() const override { return volumes_.size(); } diff --git a/src/embree/ray_tracer.cpp b/src/embree/ray_tracer.cpp index a7484ea7..5c37499f 100644 --- a/src/embree/ray_tracer.cpp +++ b/src/embree/ray_tracer.cpp @@ -53,6 +53,7 @@ EmbreeRayTracer::create_surface_tree(const std::shared_ptr& mesh_ma SurfaceTreeID tree = next_surface_tree_id(); surface_trees_.push_back(tree); auto volume_scene = this->create_embree_scene(); + auto volume_surfaces = mesh_manager->get_volume_surfaces(volume_id); // allocate total storage for all the primtives in a volume diff --git a/src/omega_h/mesh_manager.cpp b/src/omega_h/mesh_manager.cpp index 97c7640c..a02a9a7d 100644 --- a/src/omega_h/mesh_manager.cpp +++ b/src/omega_h/mesh_manager.cpp @@ -126,13 +126,10 @@ void OmegaHMeshManager::determine_surface_senses() { // the parent volumes of a surface are the volumes of the regions on either // side of its faces. The first region encountered defines the forward sense // and the opposing region, if any, defines the reverse sense. - auto face_to_region = mesh_->ask_up(OMEGA_H_FACE, OMEGA_H_REGION); - auto face_class_dim = - mesh_->get_array(OMEGA_H_FACE, CLASS_DIM_TAG); - auto face_class_id = - mesh_->get_array(OMEGA_H_FACE, CLASS_ID_TAG); - auto region_class_id = - mesh_->get_array(OMEGA_H_REGION, CLASS_ID_TAG); + auto face_to_region = mesh_->ask_up(OMEGA_H_FACE, OMEGA_H_REGION); + auto face_class_dim = mesh_->get_array(OMEGA_H_FACE, CLASS_DIM_TAG); + auto face_class_id = mesh_->get_array(OMEGA_H_FACE, CLASS_ID_TAG); + auto region_class_id = mesh_->get_array(OMEGA_H_REGION, CLASS_ID_TAG); for (auto surface : surfaces_) { surface_senses_[surface] = {ID_NONE, ID_NONE}; @@ -364,4 +361,6 @@ Sense OmegaHMeshManager::surface_sense(MeshID surface, MeshID volume) const { return volume == senses.first ? Sense::FORWARD : Sense::REVERSE; } + + } // namespace xdg \ No newline at end of file diff --git a/tests/test_omega_h.cpp b/tests/test_omega_h.cpp index eaaa2e11..50932c3e 100644 --- a/tests/test_omega_h.cpp +++ b/tests/test_omega_h.cpp @@ -16,10 +16,7 @@ using namespace xdg; using namespace xdg::test; -// These tests use "brick.exo". -// Its classification therefore describes one model volume bounded by six planar -// model surfaces. If the converted mesh in your test data differs, update the -// expected counts and ray-fire distances below to match. + TEST_CASE("Test Omega_h Initialization") { std::unique_ptr mesh_manager = @@ -35,10 +32,8 @@ TEST_CASE("Test Omega_h Initialization") { REQUIRE(mesh_manager->implicit_complement() != ID_NONE); // num_ents_of_dimension should agree with the volume and surface counts - REQUIRE(mesh_manager->num_ents_of_dimension(3) == - mesh_manager->num_volumes()); - REQUIRE(mesh_manager->num_ents_of_dimension(2) == - mesh_manager->num_surfaces()); + REQUIRE(mesh_manager->num_ents_of_dimension(3) == mesh_manager->num_volumes()); + REQUIRE(mesh_manager->num_ents_of_dimension(2) == mesh_manager->num_surfaces()); // every surface is bounded by the cube volume on its forward side MeshID volume = mesh_manager->volumes().front(); @@ -59,17 +54,16 @@ TEST_CASE("Omega_h Volume Elements") { MeshID volume = mesh_manager->volumes().front(); auto elements = mesh_manager->get_volume_elements(volume); REQUIRE(!elements.empty()); - REQUIRE(mesh_manager->num_volume_elements(volume) == - static_cast(elements.size())); + REQUIRE(mesh_manager->num_volume_elements(volume) == static_cast(elements.size())); // the implicit complement holds no elements - REQUIRE(mesh_manager->num_volume_elements( - mesh_manager->implicit_complement()) == 0); + REQUIRE(mesh_manager->num_volume_elements( mesh_manager->implicit_complement()) == 0); // the global element count is the sum over all volumes and equals the number // of regions (tetrahedra) in the mesh int total = 0; - for (auto v : mesh_manager->volumes()) { + for (auto v : mesh_manager->volumes()) + { total += mesh_manager->num_volume_elements(v); } REQUIRE(mesh_manager->num_volume_elements() == total); @@ -174,17 +168,23 @@ TEMPLATE_TEST_CASE("Test BVH Build Omega_h", "[omega_h][bvh]", mesh_manager->load_file("pincell-implicit.exo"); mesh_manager->init(); + REQUIRE(mesh_manager->num_volume_elements() == 93170) ; REQUIRE(mesh_manager->num_volumes() == 5); REQUIRE(mesh_manager->num_surfaces() == 13); DYNAMIC_SECTION(fmt::format("Backend = {}", rt_backend)) { - check_ray_tracer_supported( - rt_backend); // skip if backend not enabled at configuration time + check_ray_tracer_supported( rt_backend); // skip if backend not enabled at configuration time auto rti = create_raytracer(rt_backend); - for (const auto &volume : mesh_manager->volumes()) { + auto volume_surfaces = mesh_manager->get_volume_surfaces(mesh_manager->volumes().front()); + std::cout<<"number of surfaces in this volume is "<volumes()){ rti->register_volume(mesh_manager, volume); } + REQUIRE(rti->num_registered_trees() == 2); } } @@ -240,8 +240,7 @@ TEMPLATE_TEST_CASE("Test Omega_h Find Element Method", "[omega_h][elements]", MeshID volume = mesh_manager->volumes().front(); - MeshID element = xdg->find_element( - volume, {0.0, 0.0, 100.0}); // I will need to fix this one as well + MeshID element = xdg->find_element( volume, {0.0, 0.0, 100.0}); // I will need to fix this one as well REQUIRE(element == ID_NONE); // point lies outside the cube element = xdg->find_element(volume, {0.0, 0.0, 0.0}); @@ -278,10 +277,9 @@ TEST_CASE("Omega_h Element ID and Index Mapping") { mesh_manager->load_file("brick.exo"); mesh_manager->init(); - // Omega_h stores entities in a contiguous, zero-based index space, so IDs - and - // indices are identical for both elements and vertices - size_t num_elements = mesh_manager->num_volume_elements(); + // Omega_h stores entities in a contiguous, zero-based index space, so IDs and + // indices are identical for both elements and vertices + size_t num_elements = mesh_manager->num_volume_elements(); REQUIRE(num_elements > 0); for (size_t idx = 0; idx < num_elements; ++idx) { MeshID element_id = mesh_manager->element_id(idx); From cc0db2dc060d6de8a346397168b6e4b6937c0b03 Mon Sep 17 00:00:00 2001 From: magnoxemo Date: Mon, 27 Jul 2026 11:15:32 -0500 Subject: [PATCH 20/26] trying to fix ci by providing seacas prefix path --- .github/workflows/ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 9e3f61e9..09e22aa1 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -223,7 +223,7 @@ jobs: if [ -n "$CMAKE_PREFIX_PATH" ]; then CMAKE_PREFIX_PATH="$CMAKE_PREFIX_PATH;" fi - CMAKE_PREFIX_PATH="${CMAKE_PREFIX_PATH}$HOME/OMEGA_H" + CMAKE_PREFIX_PATH="${CMAKE_PREFIX_PATH}$HOME/OMEGA_H; $HOME/SEACAS" fi cmake .. \ -DCMAKE_PREFIX_PATH="$CMAKE_PREFIX_PATH" \ From 634c3a1000ded639424ec278380d2358b2de7001 Mon Sep 17 00:00:00 2001 From: Ebny Walid Ahammed <69362074+magnoxemo@users.noreply.github.com> Date: Sun, 2 Aug 2026 17:02:49 -0500 Subject: [PATCH 21/26] Fix CMAKE_PREFIX_PATH formatting in CI workflow --- .github/workflows/ci.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 09e22aa1..383d5034 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -223,7 +223,7 @@ jobs: if [ -n "$CMAKE_PREFIX_PATH" ]; then CMAKE_PREFIX_PATH="$CMAKE_PREFIX_PATH;" fi - CMAKE_PREFIX_PATH="${CMAKE_PREFIX_PATH}$HOME/OMEGA_H; $HOME/SEACAS" + CMAKE_PREFIX_PATH="${CMAKE_PREFIX_PATH}$HOME/OMEGA_H;$HOME/SEACAS" fi cmake .. \ -DCMAKE_PREFIX_PATH="$CMAKE_PREFIX_PATH" \ @@ -247,4 +247,4 @@ jobs: - name: Setup tmate session if: ${{ failure() }} - uses: mxschmitt/action-tmate@v3 \ No newline at end of file + uses: mxschmitt/action-tmate@v3 From 7a6d1cf336d67298048a731102476474baa53f70 Mon Sep 17 00:00:00 2001 From: magnoxemo Date: Sun, 2 Aug 2026 17:07:02 -0500 Subject: [PATCH 22/26] regold test values --- include/xdg/omega_h/mesh_manager.h | 21 +++- src/element_face_accessor.cpp | 10 ++ src/omega_h/mesh_manager.cpp | 165 +++++++++++++++++++++-------- tests/test_omega_h.cpp | 32 ++++-- 4 files changed, 173 insertions(+), 55 deletions(-) diff --git a/include/xdg/omega_h/mesh_manager.h b/include/xdg/omega_h/mesh_manager.h index 4d0f5a7d..702e428d 100644 --- a/include/xdg/omega_h/mesh_manager.h +++ b/include/xdg/omega_h/mesh_manager.h @@ -138,9 +138,28 @@ class OmegaHMeshManager : public MeshManager { //! \brief Map element and vertex ID spaces into contiguous index spaces void map_id_spaces(); - std::unique_ptr library_; + //! \brief Fetch Omega_h's derived adjacencies and coordinate array exactly + //! once, single-threaded, and store them as members. Omega_h::Read<>/Adj + //! objects are reference-counted handles (Omega_h::SharedAlloc) whose + //! copy/destroy operations are NOT thread-safe. Every accessor below reads + //! only from these cached members instead of calling mesh_->ask_*()/coords() + //! itself, so no Read<>/Adj handle is ever copied or destroyed again after + //! init() -- which is what let concurrent Embree bounds-callback threads + //! race on Omega_h's internal refcounts and cause a double-free. + void cache_derived_arrays(); + std::unique_ptr mesh_; + //! Cached derived arrays/adjacencies (see cache_derived_arrays()). Never + //! touch mesh_->ask_*()/coords() anywhere outside of that function -- read + //! from these members instead so no Omega_h::Read<>/Adj handle is copied or + //! destroyed from more than one thread. + Omega_h::Reals coords_; + Omega_h::LOs elem_verts_; + Omega_h::LOs face_verts_; + Omega_h::Adj region_to_face_; + Omega_h::Adj face_to_region_; + //! Mapping of surfaces to the volumes on either side. Volumes are ordered //! based on their sense with respect to the surface triangles std::map> surface_senses_; diff --git a/src/element_face_accessor.cpp b/src/element_face_accessor.cpp index 4b965285..17e66090 100644 --- a/src/element_face_accessor.cpp +++ b/src/element_face_accessor.cpp @@ -9,6 +9,10 @@ #include "xdg/libmesh/mesh_manager.h" #endif +#ifdef XDG_ENABLE_OMEGA_H +#include "xdg/omega_h/mesh_manager.h" +#endif + #include "xdg/testing/mesh_mock.h" namespace xdg { @@ -26,6 +30,12 @@ std::shared_ptr ElementFaceAccessor::create(const MeshManag return std::make_shared(libmesh_mesh_manager, element); } #endif + #ifdef XDG_ENABLE_OMEGA_H + if (mesh_manager->mesh_library() == MeshLibrary::OMEGA_H) { + const OmegaHMeshManager* omega_h_mesh_manager = dynamic_cast(mesh_manager); + return std::make_shared(omega_h_mesh_manager, element); + } + #endif // for testing if (mesh_manager->mesh_library() == MeshLibrary::MOCK) { const MeshMock* mock_mesh_manager = dynamic_cast(mesh_manager); diff --git a/src/omega_h/mesh_manager.cpp b/src/omega_h/mesh_manager.cpp index a02a9a7d..2ade43e7 100644 --- a/src/omega_h/mesh_manager.cpp +++ b/src/omega_h/mesh_manager.cpp @@ -29,18 +29,37 @@ static constexpr int VERTS_PER_TET = 4; static constexpr int VERTS_PER_TRI = 3; static constexpr int FACES_PER_TET = 4; +// Omega_h's canonical tet local-face -> local-vertex template (elem_dim=3, +// bdry_dim=2 case of simplex_down_template() in Omega_h_simplex.hpp). Indexing +// an element's own vertex connectivity with this table always yields a +// winding whose normal points outward from that specific element -- the +// element's own natural vertex order alone does not, since a face's global +// (face->vertex) order is only outward-facing for one of its (up to two) +// adjacent elements. +static constexpr int kLocalFaceVerts[FACES_PER_TET][VERTS_PER_TRI] = { + {0, 2, 1}, {0, 1, 3}, {1, 2, 3}, {2, 0, 3}}; + // Constructor -OmegaHMeshManager::OmegaHMeshManager(const Omega_h::Mesh* mesh) { +namespace { +// Omega_h::Library wraps process-wide MPI/Kokkos init and finalize calls and +// must only be constructed once per process -- Kokkos aborts (and corrupts +// its allocator state) if initialized more than once. Share a single, +// lazily-constructed instance across every OmegaHMeshManager rather than +// creating one per manager, since Catch2 constructs a fresh manager in +// nearly every TEST_CASE. +Omega_h::Library& shared_omega_h_library() { + static Omega_h::Library lib(nullptr, nullptr); + return lib; +} +} // namespace +OmegaHMeshManager::OmegaHMeshManager(const Omega_h::Mesh* mesh) { // I will rewrite this method later - library_ = std::make_unique(nullptr, nullptr); - mesh_ = std::make_unique(library_.get()); - + mesh_ = std::make_unique(&shared_omega_h_library()); } OmegaHMeshManager::OmegaHMeshManager() { - library_ = std::make_unique(nullptr, nullptr); - mesh_ = std::make_unique(library_.get()); + mesh_ = std::make_unique(&shared_omega_h_library()); } @@ -60,6 +79,11 @@ void OmegaHMeshManager::init() { // in Omega_h the regions are 3D simplices or aka elements num_elements_ = mesh_->nregions(); + // Fetch every derived array/adjacency this class needs exactly once, here, + // single-threaded. See the member declarations in the header for why: it + // is not safe to call mesh_->ask_*()/coords() again after this point. + cache_derived_arrays(); + // a classified mesh defines its volumes and surfaces directly through the // class_dim/class_id tags. Otherwise, treat the entire mesh as a single volume // bounded by its exposed faces. @@ -77,6 +101,27 @@ void OmegaHMeshManager::init() { map_id_spaces(); } +void OmegaHMeshManager::cache_derived_arrays() { + // Omega_h::Read<>/Adj objects are reference-counted handles + // (Omega_h::SharedAlloc) whose copy/destroy operations are NOT + // thread-safe. Ray tracer BVH builds (e.g. Embree with RTC_BUILD_QUALITY_HIGH) + // invoke element_connectivity()/face_connectivity()/vertex_coordinates()/ + // element_face_vertices()/adjacent_element() from many worker threads in + // parallel via their bounds/intersect callbacks. Those methods used to call + // mesh_->ask_elem_verts()/ask_verts_of()/ask_down()/ask_up()/coords() + // directly, which copies Omega_h's internally cached shared array into a + // new Read<>/Adj handle on every single call -- racing the refcount of the + // underlying SharedAlloc across threads and causing a double-free once two + // threads both drop it to zero. Fetching each array exactly once here and + // reading only from these members afterward means no Read<>/Adj handle is + // ever copied or destroyed from more than one thread. + coords_ = mesh_->coords(); + elem_verts_ = mesh_->ask_elem_verts(); + face_verts_ = mesh_->ask_verts_of(OMEGA_H_FACE); + region_to_face_ = mesh_->ask_down(OMEGA_H_REGION, OMEGA_H_FACE); + face_to_region_ = mesh_->ask_up(OMEGA_H_FACE, OMEGA_H_REGION); +} + bool OmegaHMeshManager::has_classification() const { return mesh_->has_tag(OMEGA_H_REGION, CLASS_DIM_TAG) && mesh_->has_tag(OMEGA_H_REGION, CLASS_ID_TAG) && @@ -123,10 +168,16 @@ void OmegaHMeshManager::discover_single_volume() { } void OmegaHMeshManager::determine_surface_senses() { - // the parent volumes of a surface are the volumes of the regions on either - // side of its faces. The first region encountered defines the forward sense - // and the opposing region, if any, defines the reverse sense. - auto face_to_region = mesh_->ask_up(OMEGA_H_FACE, OMEGA_H_REGION); + // The parent volumes of a surface are the volumes of the regions on either + // side of its faces. Which region is "forward" vs "reverse" isn't just a + // matter of adjacency order: face_vertices()/face_normal() (used by the ray + // tracer) always return a face's fixed global winding, which only points + // outward from ONE of its (up to two) adjacent regions -- the other sees + // an inward-pointing normal that the ray tracer flips based on this + // forward/reverse assignment (see EmbreeRayTracer::register_surface() and + // TriangleIntersectionFunc()). So the region whose canonical outward local + // face winding (see kLocalFaceVerts) matches the face's actual global + // winding must be recorded as "forward"; the other, if any, as "reverse". auto face_class_dim = mesh_->get_array(OMEGA_H_FACE, CLASS_DIM_TAG); auto face_class_id = mesh_->get_array(OMEGA_H_FACE, CLASS_ID_TAG); auto region_class_id = mesh_->get_array(OMEGA_H_REGION, CLASS_ID_TAG); @@ -142,15 +193,50 @@ void OmegaHMeshManager::determine_surface_senses() { MeshID surface = static_cast(face_class_id.get(face)); auto &senses = surface_senses_[surface]; - Omega_h::LO begin = face_to_region.a2ab.get(face); - Omega_h::LO end = face_to_region.a2ab.get(face + 1); + Omega_h::LO fbase = face * VERTS_PER_TRI; + std::array natural_order = { + face_verts_.get(fbase + 0), face_verts_.get(fbase + 1), + face_verts_.get(fbase + 2)}; + + Omega_h::LO begin = face_to_region_.a2ab.get(face); + Omega_h::LO end = face_to_region_.a2ab.get(face + 1); for (Omega_h::LO k = begin; k < end; ++k) { - MeshID volume = - static_cast(region_class_id.get(face_to_region.ab2b.get(k))); - if (senses.first == ID_NONE) { - senses.first = volume; // forward sense - } else if (senses.second == ID_NONE && senses.first != volume) { - senses.second = volume; // reverse sense + Omega_h::LO region = face_to_region_.ab2b.get(k); + MeshID volume = static_cast(region_class_id.get(region)); + + // find which of the region's 4 local faces this global face is + Omega_h::LO rbase = region * FACES_PER_TET; + int local_face = -1; + for (int lf = 0; lf < FACES_PER_TET; ++lf) { + if (region_to_face_.ab2b.get(rbase + lf) == face) { + local_face = lf; + break; + } + } + + Omega_h::LO ebase = region * VERTS_PER_TET; + std::array outward_order; + for (int i = 0; i < VERTS_PER_TRI; ++i) { + outward_order[i] = elem_verts_.get(ebase + kLocalFaceVerts[local_face][i]); + } + + // the two possible windings of 3 vertices are cyclic rotations of + // either `outward_order` or its reverse; find where outward_order[0] + // falls in natural_order and compare the next element to tell them apart + int rot = 0; + while (natural_order[rot] != outward_order[0]) ++rot; + bool is_outward = natural_order[(rot + 1) % VERTS_PER_TRI] == outward_order[1]; + + // Only assign each slot once: a well-formed surface's every face + // agrees on which side is outward, so this is a no-op after the first + // face. Some meshes carry degenerate faces classified onto a surface + // that borders more than two distinct volumes (not a true 2-manifold + // interface); for those, keep whichever pair of volumes was recorded + // first rather than letting later faces keep overwriting the senses. + if (is_outward) { + if (senses.first == ID_NONE) senses.first = volume; + } else { + if (senses.second == ID_NONE && senses.first != volume) senses.second = volume; } } } @@ -217,10 +303,9 @@ std::vector OmegaHMeshManager::get_surface_faces(MeshID surface) const { // without classification the surface is the set of exposed boundary faces // (faces adjacent to exactly one region) if (!mesh_->has_tag(OMEGA_H_FACE, CLASS_ID_TAG)) { - auto face_to_region = mesh_->ask_up(OMEGA_H_FACE, OMEGA_H_REGION); for (Omega_h::LO face = 0; face < mesh_->nfaces(); ++face) { Omega_h::LO n_adjacent = - face_to_region.a2ab.get(face + 1) - face_to_region.a2ab.get(face); + face_to_region_.a2ab.get(face + 1) - face_to_region_.a2ab.get(face); if (n_adjacent == 1) { faces.push_back(static_cast(face)); } @@ -243,42 +328,38 @@ std::vector OmegaHMeshManager::get_surface_faces(MeshID surface) const { std::vector OmegaHMeshManager::element_connectivity(MeshID element) const { // tetrahedra store four vertices per element in element->vertex order - auto elem_verts = mesh_->ask_elem_verts(); Omega_h::LO base = static_cast(element) * VERTS_PER_TET; std::vector connectivity(VERTS_PER_TET); for (int i = 0; i < VERTS_PER_TET; ++i) { - connectivity[i] = static_cast(elem_verts.get(base + i)); + connectivity[i] = static_cast(elem_verts_.get(base + i)); } return connectivity; } std::vector OmegaHMeshManager::face_connectivity(MeshID face) const { // triangles store three vertices per face in face->vertex order - auto face_verts = mesh_->ask_verts_of(OMEGA_H_FACE); Omega_h::LO base = static_cast(face) * VERTS_PER_TRI; std::vector connectivity(VERTS_PER_TRI); for (int i = 0; i < VERTS_PER_TRI; ++i) { - connectivity[i] = static_cast(face_verts.get(base + i)); + connectivity[i] = static_cast(face_verts_.get(base + i)); } return connectivity; } MeshID OmegaHMeshManager::get_boundary_face_element(MeshID face) const { // the owning element of a boundary face is its single adjacent region - auto face_to_region = mesh_->ask_up(OMEGA_H_FACE, OMEGA_H_REGION); - Omega_h::LO begin = face_to_region.a2ab.get(face); - Omega_h::LO end = face_to_region.a2ab.get(face + 1); + Omega_h::LO begin = face_to_region_.a2ab.get(face); + Omega_h::LO end = face_to_region_.a2ab.get(face + 1); if (begin == end) { return ID_NONE; } - return static_cast(face_to_region.ab2b.get(begin)); + return static_cast(face_to_region_.ab2b.get(begin)); } Vertex OmegaHMeshManager::vertex_coordinates(MeshID vertex) const { - // coords() is a flat array laid out as [x0, y0, z0, x1, y1, z1, ...] - auto coords = mesh_->coords(); + // coords_ is a flat array laid out as [x0, y0, z0, x1, y1, z1, ...] Omega_h::LO base = static_cast(vertex) * 3; - return {coords.get(base), coords.get(base + 1), coords.get(base + 2)}; + return {coords_.get(base), coords_.get(base + 1), coords_.get(base + 2)}; } std::vector OmegaHMeshManager::element_vertices(MeshID element) const { @@ -302,26 +383,24 @@ std::array OmegaHMeshManager::face_vertices(MeshID face) const { std::array OmegaHMeshManager::element_face_vertices(MeshID element, int local_face) const { - // resolve the global face index for the requested local face, then return its - // vertices in their natural (face->vertex) orientation - auto region_to_face = mesh_->ask_down(OMEGA_H_REGION, OMEGA_H_FACE); - Omega_h::LO base = static_cast(element) * FACES_PER_TET; - Omega_h::LO global_face = region_to_face.ab2b.get(base + local_face); - return face_vertices(static_cast(global_face)); + auto connectivity = element_connectivity(element); + std::array vertices; + for (int i = 0; i < VERTS_PER_TRI; ++i) { + vertices[i] = vertex_coordinates(connectivity[kLocalFaceVerts[local_face][i]]); + } + return vertices; } MeshID OmegaHMeshManager::adjacent_element(MeshID element, int face) const { // resolve the global face index for the requested local face of the element - auto region_to_face = mesh_->ask_down(OMEGA_H_REGION, OMEGA_H_FACE); Omega_h::LO base = static_cast(element) * FACES_PER_TET; - Omega_h::LO global_face = region_to_face.ab2b.get(base + face); + Omega_h::LO global_face = region_to_face_.ab2b.get(base + face); // the neighbor is the other region sharing that face, if any - auto face_to_region = mesh_->ask_up(OMEGA_H_FACE, OMEGA_H_REGION); - Omega_h::LO begin = face_to_region.a2ab.get(global_face); - Omega_h::LO end = face_to_region.a2ab.get(global_face + 1); + Omega_h::LO begin = face_to_region_.a2ab.get(global_face); + Omega_h::LO end = face_to_region_.a2ab.get(global_face + 1); for (Omega_h::LO k = begin; k < end; ++k) { - MeshID neighbor = static_cast(face_to_region.ab2b.get(k)); + MeshID neighbor = static_cast(face_to_region_.ab2b.get(k)); if (neighbor != element) { return neighbor; } diff --git a/tests/test_omega_h.cpp b/tests/test_omega_h.cpp index 50932c3e..36cff6df 100644 --- a/tests/test_omega_h.cpp +++ b/tests/test_omega_h.cpp @@ -179,13 +179,14 @@ TEMPLATE_TEST_CASE("Test BVH Build Omega_h", "[omega_h][bvh]", auto volume_surfaces = mesh_manager->get_volume_surfaces(mesh_manager->volumes().front()); std::cout<<"number of surfaces in this volume is "<volumes()){ rti->register_volume(mesh_manager, volume); } - REQUIRE(rti->num_registered_trees() == 2); + // one surface tree per volume, plus one element tree per volume that + // actually contains elements (the implicit complement has none, so it + // contributes a surface tree only): 5 surface trees + 4 element trees + REQUIRE(rti->num_registered_trees() == 9); } } @@ -204,19 +205,23 @@ TEMPLATE_TEST_CASE("Test Ray Fire Omega_h (all built backends)", mm->init(); xdg->prepare_raytracer(); - MeshID volume = mm->volumes()[0]; + // Volume 6 (the moderator region) is bounded by a clean, fully-classified + // set of surfaces (7-12), unlike volumes 1/2 whose only registered + // boundary is a degenerate surface (class_id -1) spanning several + // unrelated volume interfaces in this mesh's classification data. Its + // outer boundary is an axis-aligned box extending to +/-25 on each axis. + MeshID volume = 6; Position origin{0.0, 0.0, 0.0}; Direction dir{1.0, 0.0, 0.0}; - // the cube is 10 units on a side and centered at the origin auto hit = xdg->ray_fire(volume, origin, dir); REQUIRE(hit.second != ID_NONE); - REQUIRE_THAT(hit.first, Catch::Matchers::WithinAbs(5.0, 1e-6)); + REQUIRE_THAT(hit.first, Catch::Matchers::WithinAbs(25.0, 1e-6)); origin = {3.0, 0.0, 0.0}; hit = xdg->ray_fire(volume, origin, dir); - REQUIRE_THAT(hit.first, Catch::Matchers::WithinAbs(2.0, 1e-6)); + REQUIRE_THAT(hit.first, Catch::Matchers::WithinAbs(22.0, 1e-6)); origin = {0.0, 0.0, 0.0}; REQUIRE(xdg->point_in_volume(volume, origin)); @@ -274,7 +279,7 @@ TEST_CASE("Omega_h Element ID and Index Mapping") { std::unique_ptr mesh_manager = std::make_unique(); REQUIRE(mesh_manager->mesh_library() == MeshLibrary::OMEGA_H); - mesh_manager->load_file("brick.exo"); + mesh_manager->load_file("brick-sidesets.exo"); mesh_manager->init(); // Omega_h stores entities in a contiguous, zero-based index space, so IDs and @@ -304,14 +309,19 @@ TEST_CASE("Test Track Exiting Mesh Omega_h") { mesh_manager->init(); xdg->prepare_raytracer(); - MeshID volume = mesh_manager->volumes().front(); + // volume 6 (the moderator) is only used to locate the starting element; + // segments() then walks element adjacency across every material the ray + // actually crosses (see volume 1/2's note in the Ray Fire test above for + // why they aren't suitable here) + MeshID volume = 6; Position start{0.0, 0.0, -1000.0}; Position end{0.0, 0.0, 1000.0}; auto tracks = xdg->segments(volume, start, end); - // the accumulated track length through the cube equals its 10 unit extent + // the accumulated track length equals the combined axial extent of every + // material the ray crosses through the pincell assembly double length = std::accumulate( tracks.begin(), tracks.end(), 0.0, [](double sum, const auto &track) { return sum + track.second; }); - REQUIRE_THAT(length, Catch::Matchers::WithinAbs(10.0, 1e-6)); + REQUIRE_THAT(length, Catch::Matchers::WithinAbs(52.5, 1e-6)); } \ No newline at end of file From 7b36ae7534562e766b16716a47711bd3e5bed708 Mon Sep 17 00:00:00 2001 From: magnoxemo Date: Tue, 4 Aug 2026 16:12:37 -0500 Subject: [PATCH 23/26] cleaning up --- include/xdg/omega_h/mesh_manager.h | 34 ++-------- src/omega_h/mesh_manager.cpp | 103 ++++++++++------------------- tests/test_omega_h.cpp | 31 ++++----- 3 files changed, 52 insertions(+), 116 deletions(-) diff --git a/include/xdg/omega_h/mesh_manager.h b/include/xdg/omega_h/mesh_manager.h index 702e428d..f27d43a5 100644 --- a/include/xdg/omega_h/mesh_manager.h +++ b/include/xdg/omega_h/mesh_manager.h @@ -17,35 +17,21 @@ namespace xdg { //! \brief Mesh manager backed by the Omega_h simplex mesh library. -//! -//! Omega_h tetrahedral meshes carry the geometric model they were generated -//! from as "class_dim"/"class_id" classification tags on each entity. These -//! tags map mesh entities onto the volumes (3D model entities) and surfaces -//! (2D model entities) of the geometry, which is all the topology XDG needs. -//! When a mesh has no classification, the whole mesh is treated as a single -//! volume bounded by its exposed faces. class OmegaHMeshManager : public MeshManager { public: OmegaHMeshManager(); - OmegaHMeshManager(const Omega_h::Mesh* mesh); - ~OmegaHMeshManager() override = default; - // Interface methods MeshLibrary mesh_library() const override { return MeshLibrary::OMEGA_H; } - void load_file(const std::string &file_path) override; - void init() override; // TODO: I will have to implement this one. void parse_metadata() override {}; int num_volumes() const override { return volumes_.size(); } - int num_surfaces() const override { return surfaces_.size(); } - int num_ents_of_dimension(int dim) const override { switch (dim) { case 3: @@ -94,11 +80,12 @@ class OmegaHMeshManager : public MeshManager { std::array face_vertices(MeshID face) const override; //! \brief Vertices of a local face of a tetrahedral element - std::array element_face_vertices(MeshID element, int local_face) const; + std::array element_face_vertices(MeshID element, + int local_face) const; SurfaceElementType get_surface_element_type(MeshID surface_element_id) const override { - // Omega_h simplex meshes always use triangular surface elements + // As Omega_h simplex meshes always use triangular surface elements return SurfaceElementType::TRI; } @@ -106,7 +93,7 @@ class OmegaHMeshManager : public MeshManager { double element_volume(MeshID element) const override; - MeshID create_volume() override; + MeshID create_volume() override { return next_volume_id(); }; void add_surface_to_volume(MeshID volume, MeshID surface, Sense sense, bool overwrite = false) override; @@ -139,21 +126,12 @@ class OmegaHMeshManager : public MeshManager { void map_id_spaces(); //! \brief Fetch Omega_h's derived adjacencies and coordinate array exactly - //! once, single-threaded, and store them as members. Omega_h::Read<>/Adj - //! objects are reference-counted handles (Omega_h::SharedAlloc) whose - //! copy/destroy operations are NOT thread-safe. Every accessor below reads - //! only from these cached members instead of calling mesh_->ask_*()/coords() - //! itself, so no Read<>/Adj handle is ever copied or destroyed again after - //! init() -- which is what let concurrent Embree bounds-callback threads - //! race on Omega_h's internal refcounts and cause a double-free. + //! once, single-threaded, and store them as members. void cache_derived_arrays(); std::unique_ptr mesh_; - //! Cached derived arrays/adjacencies (see cache_derived_arrays()). Never - //! touch mesh_->ask_*()/coords() anywhere outside of that function -- read - //! from these members instead so no Omega_h::Read<>/Adj handle is copied or - //! destroyed from more than one thread. + //! Cached derived arrays/adjacencies Omega_h::Reals coords_; Omega_h::LOs elem_verts_; Omega_h::LOs face_verts_; diff --git a/src/omega_h/mesh_manager.cpp b/src/omega_h/mesh_manager.cpp index 2ade43e7..8ec678fc 100644 --- a/src/omega_h/mesh_manager.cpp +++ b/src/omega_h/mesh_manager.cpp @@ -23,19 +23,12 @@ static constexpr int CLASS_DIM_VOLUME = 3; // Classification tag names written by Omega_h static constexpr const char *CLASS_DIM_TAG = "class_dim"; -static constexpr const char *CLASS_ID_TAG = "class_id"; +static constexpr const char *CLASS_ID_TAG = "class_id"; static constexpr int VERTS_PER_TET = 4; static constexpr int VERTS_PER_TRI = 3; static constexpr int FACES_PER_TET = 4; -// Omega_h's canonical tet local-face -> local-vertex template (elem_dim=3, -// bdry_dim=2 case of simplex_down_template() in Omega_h_simplex.hpp). Indexing -// an element's own vertex connectivity with this table always yields a -// winding whose normal points outward from that specific element -- the -// element's own natural vertex order alone does not, since a face's global -// (face->vertex) order is only outward-facing for one of its (up to two) -// adjacent elements. static constexpr int kLocalFaceVerts[FACES_PER_TET][VERTS_PER_TRI] = { {0, 2, 1}, {0, 1, 3}, {1, 2, 3}, {2, 0, 3}}; @@ -44,30 +37,21 @@ static constexpr int kLocalFaceVerts[FACES_PER_TET][VERTS_PER_TRI] = { namespace { // Omega_h::Library wraps process-wide MPI/Kokkos init and finalize calls and // must only be constructed once per process -- Kokkos aborts (and corrupts -// its allocator state) if initialized more than once. Share a single, -// lazily-constructed instance across every OmegaHMeshManager rather than -// creating one per manager, since Catch2 constructs a fresh manager in -// nearly every TEST_CASE. -Omega_h::Library& shared_omega_h_library() { +// its allocator state) if initialized more than once. +Omega_h::Library &shared_omega_h_library() { static Omega_h::Library lib(nullptr, nullptr); return lib; } } // namespace -OmegaHMeshManager::OmegaHMeshManager(const Omega_h::Mesh* mesh) { - // I will rewrite this method later - mesh_ = std::make_unique(&shared_omega_h_library()); -} OmegaHMeshManager::OmegaHMeshManager() { mesh_ = std::make_unique(&shared_omega_h_library()); } - -void OmegaHMeshManager::load_file(const std::string &file_path) -{ - const int exodus_file = Omega_h::exodus::open(file_path); - Omega_h::exodus::read_mesh(exodus_file, mesh_.get()); - Omega_h::exodus::close(exodus_file); +void OmegaHMeshManager::load_file(const std::string &file_path) { + const int exodus_file = Omega_h::exodus::open(file_path); + Omega_h::exodus::read_mesh(exodus_file, mesh_.get()); + Omega_h::exodus::close(exodus_file); } void OmegaHMeshManager::init() { @@ -79,14 +63,13 @@ void OmegaHMeshManager::init() { // in Omega_h the regions are 3D simplices or aka elements num_elements_ = mesh_->nregions(); - // Fetch every derived array/adjacency this class needs exactly once, here, - // single-threaded. See the member declarations in the header for why: it - // is not safe to call mesh_->ask_*()/coords() again after this point. + // We should cache all necessary derived objects as rediscovering those thing + // during run time causes performance penalties. cache_derived_arrays(); // a classified mesh defines its volumes and surfaces directly through the - // class_dim/class_id tags. Otherwise, treat the entire mesh as a single volume - // bounded by its exposed faces. + // class_dim/class_id tags. Otherwise, treat the entire mesh as a single + // volume bounded by its exposed faces. if (has_classification()) { discover_geometry(); determine_surface_senses(); @@ -102,19 +85,6 @@ void OmegaHMeshManager::init() { } void OmegaHMeshManager::cache_derived_arrays() { - // Omega_h::Read<>/Adj objects are reference-counted handles - // (Omega_h::SharedAlloc) whose copy/destroy operations are NOT - // thread-safe. Ray tracer BVH builds (e.g. Embree with RTC_BUILD_QUALITY_HIGH) - // invoke element_connectivity()/face_connectivity()/vertex_coordinates()/ - // element_face_vertices()/adjacent_element() from many worker threads in - // parallel via their bounds/intersect callbacks. Those methods used to call - // mesh_->ask_elem_verts()/ask_verts_of()/ask_down()/ask_up()/coords() - // directly, which copies Omega_h's internally cached shared array into a - // new Read<>/Adj handle on every single call -- racing the refcount of the - // underlying SharedAlloc across threads and causing a double-free once two - // threads both drop it to zero. Fetching each array exactly once here and - // reading only from these members afterward means no Read<>/Adj handle is - // ever copied or destroyed from more than one thread. coords_ = mesh_->coords(); elem_verts_ = mesh_->ask_elem_verts(); face_verts_ = mesh_->ask_verts_of(OMEGA_H_FACE); @@ -131,8 +101,10 @@ bool OmegaHMeshManager::has_classification() const { void OmegaHMeshManager::discover_geometry() { // volumes are the unique class IDs of regions classified on a 3D model entity - auto region_class_dim = mesh_->get_array(OMEGA_H_REGION, CLASS_DIM_TAG); - auto region_class_id = mesh_->get_array(OMEGA_H_REGION, CLASS_ID_TAG); + auto region_class_dim = + mesh_->get_array(OMEGA_H_REGION, CLASS_DIM_TAG); + auto region_class_id = + mesh_->get_array(OMEGA_H_REGION, CLASS_ID_TAG); std::set volume_ids; for (Omega_h::LO region = 0; region < mesh_->nregions(); ++region) { @@ -150,9 +122,8 @@ void OmegaHMeshManager::discover_geometry() { std::set surface_ids; for (Omega_h::LO face = 0; face < mesh_->nfaces(); ++face) { - if (face_class_dim.get(face) == CLASS_DIM_SURFACE) { + if (face_class_dim.get(face) == CLASS_DIM_SURFACE) surface_ids.insert(static_cast(face_class_id.get(face))); - } } surfaces_.assign(surface_ids.begin(), surface_ids.end()); } @@ -168,16 +139,7 @@ void OmegaHMeshManager::discover_single_volume() { } void OmegaHMeshManager::determine_surface_senses() { - // The parent volumes of a surface are the volumes of the regions on either - // side of its faces. Which region is "forward" vs "reverse" isn't just a - // matter of adjacency order: face_vertices()/face_normal() (used by the ray - // tracer) always return a face's fixed global winding, which only points - // outward from ONE of its (up to two) adjacent regions -- the other sees - // an inward-pointing normal that the ray tracer flips based on this - // forward/reverse assignment (see EmbreeRayTracer::register_surface() and - // TriangleIntersectionFunc()). So the region whose canonical outward local - // face winding (see kLocalFaceVerts) matches the face's actual global - // winding must be recorded as "forward"; the other, if any, as "reverse". + auto face_class_dim = mesh_->get_array(OMEGA_H_FACE, CLASS_DIM_TAG); auto face_class_id = mesh_->get_array(OMEGA_H_FACE, CLASS_ID_TAG); auto region_class_id = mesh_->get_array(OMEGA_H_REGION, CLASS_ID_TAG); @@ -187,16 +149,17 @@ void OmegaHMeshManager::determine_surface_senses() { } for (Omega_h::LO face = 0; face < mesh_->nfaces(); ++face) { - if (face_class_dim.get(face) != CLASS_DIM_SURFACE) { + if (face_class_dim.get(face) != CLASS_DIM_SURFACE) continue; - } + MeshID surface = static_cast(face_class_id.get(face)); auto &senses = surface_senses_[surface]; - Omega_h::LO fbase = face * VERTS_PER_TRI; + Omega_h::LO f_base = face * VERTS_PER_TRI; std::array natural_order = { - face_verts_.get(fbase + 0), face_verts_.get(fbase + 1), - face_verts_.get(fbase + 2)}; + face_verts_.get(f_base + 0), face_verts_.get(f_base + 1), + face_verts_.get(f_base + 2) + }; Omega_h::LO begin = face_to_region_.a2ab.get(face); Omega_h::LO end = face_to_region_.a2ab.get(face + 1); @@ -217,14 +180,17 @@ void OmegaHMeshManager::determine_surface_senses() { Omega_h::LO ebase = region * VERTS_PER_TET; std::array outward_order; for (int i = 0; i < VERTS_PER_TRI; ++i) { - outward_order[i] = elem_verts_.get(ebase + kLocalFaceVerts[local_face][i]); + outward_order[i] = + elem_verts_.get(ebase + kLocalFaceVerts[local_face][i]); } // the two possible windings of 3 vertices are cyclic rotations of // either `outward_order` or its reverse; find where outward_order[0] // falls in natural_order and compare the next element to tell them apart int rot = 0; - while (natural_order[rot] != outward_order[0]) ++rot; + while (natural_order[rot] != outward_order[0]) { + ++rot; + } bool is_outward = natural_order[(rot + 1) % VERTS_PER_TRI] == outward_order[1]; // Only assign each slot once: a well-formed surface's every face @@ -234,9 +200,11 @@ void OmegaHMeshManager::determine_surface_senses() { // interface); for those, keep whichever pair of volumes was recorded // first rather than letting later faces keep overwriting the senses. if (is_outward) { - if (senses.first == ID_NONE) senses.first = volume; + if (senses.first == ID_NONE) + senses.first = volume; } else { - if (senses.second == ID_NONE && senses.first != volume) senses.second = volume; + if (senses.second == ID_NONE && senses.first != volume) + senses.second = volume; } } } @@ -254,8 +222,6 @@ void OmegaHMeshManager::map_id_spaces() { vertex_id_map_ = IDBlockMapping(vertex_ids); } -MeshID OmegaHMeshManager::create_volume() { return next_volume_id(); } - void OmegaHMeshManager::add_surface_to_volume(MeshID volume, MeshID surface, Sense sense, bool overwrite) { auto senses = surface_senses(surface); @@ -386,7 +352,8 @@ OmegaHMeshManager::element_face_vertices(MeshID element, int local_face) const { auto connectivity = element_connectivity(element); std::array vertices; for (int i = 0; i < VERTS_PER_TRI; ++i) { - vertices[i] = vertex_coordinates(connectivity[kLocalFaceVerts[local_face][i]]); + vertices[i] = + vertex_coordinates(connectivity[kLocalFaceVerts[local_face][i]]); } return vertices; } @@ -440,6 +407,4 @@ Sense OmegaHMeshManager::surface_sense(MeshID surface, MeshID volume) const { return volume == senses.first ? Sense::FORWARD : Sense::REVERSE; } - - } // namespace xdg \ No newline at end of file diff --git a/tests/test_omega_h.cpp b/tests/test_omega_h.cpp index 36cff6df..aaaa020f 100644 --- a/tests/test_omega_h.cpp +++ b/tests/test_omega_h.cpp @@ -141,14 +141,16 @@ TEST_CASE("Omega_h Adjacency") { for (auto element : mesh_manager->get_volume_elements(volume)) { for (int face = 0; face < faces_per_tet; ++face) { MeshID neighbor = mesh_manager->adjacent_element(element, face); - // a neighbor is either another element or a boundary (ID_NONE), but never - // the element itself + // a neighbor is either another element or a boundary (ID_NONE), but never the element itself REQUIRE(neighbor != element); - if (neighbor != ID_NONE) { + if (neighbor != ID_NONE) + { // adjacency is symmetric: the neighbor must list this element back bool reciprocal = false; - for (int j = 0; j < faces_per_tet; ++j) { - if (mesh_manager->adjacent_element(neighbor, j) == element) { + for (int j = 0; j < faces_per_tet; ++j) + { + if (mesh_manager->adjacent_element(neighbor, j) == element) + { reciprocal = true; break; } @@ -195,8 +197,7 @@ TEMPLATE_TEST_CASE("Test Ray Fire Omega_h (all built backends)", constexpr auto rt_backend = TestType::value; DYNAMIC_SECTION(fmt::format("Backend = {}", rt_backend)) { - check_ray_tracer_supported( - rt_backend); // skip if backend not enabled at configuration time + check_ray_tracer_supported( rt_backend); // skip if backend not enabled at configuration time auto xdg = XDG::create(MeshLibrary::OMEGA_H, rt_backend); REQUIRE(xdg->mesh_manager()->mesh_library() == MeshLibrary::OMEGA_H); @@ -205,11 +206,6 @@ TEMPLATE_TEST_CASE("Test Ray Fire Omega_h (all built backends)", mm->init(); xdg->prepare_raytracer(); - // Volume 6 (the moderator region) is bounded by a clean, fully-classified - // set of surfaces (7-12), unlike volumes 1/2 whose only registered - // boundary is a degenerate surface (class_id -1) spanning several - // unrelated volume interfaces in this mesh's classification data. Its - // outer boundary is an axis-aligned box extending to +/-25 on each axis. MeshID volume = 6; Position origin{0.0, 0.0, 0.0}; @@ -234,10 +230,10 @@ TEMPLATE_TEST_CASE("Test Omega_h Find Element Method", "[omega_h][elements]", constexpr auto rt_backend = TestType::value; DYNAMIC_SECTION(fmt::format("Backend = {}", rt_backend)) { - check_ray_tracer_supported( - rt_backend); // skip if backend not enabled at configuration time + check_ray_tracer_supported(rt_backend); // skip if backend not enabled at configuration time std::shared_ptr xdg = XDG::create(MeshLibrary::OMEGA_H, rt_backend); REQUIRE(xdg->mesh_manager()->mesh_library() == MeshLibrary::OMEGA_H); + const auto &mesh_manager = xdg->mesh_manager(); mesh_manager->load_file("pincell-implicit.exo"); mesh_manager->init(); @@ -274,7 +270,8 @@ TEMPLATE_TEST_CASE("Test Omega_h Find Element Method", "[omega_h][elements]", } } -TEST_CASE("Omega_h Element ID and Index Mapping") { +TEST_CASE("Omega_h Element ID and Index Mapping") +{ // Gold values for this test needs to be fixed. I will do that later. std::unique_ptr mesh_manager = std::make_unique(); @@ -309,10 +306,6 @@ TEST_CASE("Test Track Exiting Mesh Omega_h") { mesh_manager->init(); xdg->prepare_raytracer(); - // volume 6 (the moderator) is only used to locate the starting element; - // segments() then walks element adjacency across every material the ray - // actually crosses (see volume 1/2's note in the Ray Fire test above for - // why they aren't suitable here) MeshID volume = 6; Position start{0.0, 0.0, -1000.0}; Position end{0.0, 0.0, 1000.0}; From 461210d2e834197e78fcda2c4330b1f20e4a0456 Mon Sep 17 00:00:00 2001 From: magnoxemo Date: Tue, 4 Aug 2026 16:15:19 -0500 Subject: [PATCH 24/26] don't build docker from my own branch --- Dockerfile | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/Dockerfile b/Dockerfile index 0d115533..40feda29 100644 --- a/Dockerfile +++ b/Dockerfile @@ -83,10 +83,7 @@ ENV OMEGA_H_INSTALL_PATH=/XDG_TEST_SYSTEM/omega_h_install_dir -# build XDG -# currently building docker image from this branch. As omega_h isn't include the cmake build system -# in upstream yet. I will change this later. -RUN git clone -b omega_h --recurse-submodules https://github.com/magnoxemo/xdg.git /XDG_TEST_SYSTEM/xdg +RUN git clone --recurse-submodules https://github.com/xdg-org/xdg.git /XDG_TEST_SYSTEM/xdg WORKDIR /XDG_TEST_SYSTEM/xdg/build RUN cmake .. \ From 280fc0000f72084dbc7030e089f1661a53708933 Mon Sep 17 00:00:00 2001 From: magnoxemo Date: Sat, 8 Aug 2026 14:47:22 -0500 Subject: [PATCH 25/26] rebase and adapt to API change issue --- include/xdg/omega_h/mesh_manager.h | 21 ++++++++++++--------- src/element_face_accessor.cpp | 1 - src/omega_h/mesh_manager.cpp | 19 ++++--------------- tests/test_omega_h.cpp | 9 ++++++--- 4 files changed, 22 insertions(+), 28 deletions(-) diff --git a/include/xdg/omega_h/mesh_manager.h b/include/xdg/omega_h/mesh_manager.h index f27d43a5..bd7b1f50 100644 --- a/include/xdg/omega_h/mesh_manager.h +++ b/include/xdg/omega_h/mesh_manager.h @@ -77,21 +77,22 @@ class OmegaHMeshManager : public MeshManager { std::vector element_vertices(MeshID element) const override; - std::array face_vertices(MeshID face) const override; + std::vector face_vertices(MeshID face) const override; //! \brief Vertices of a local face of a tetrahedral element - std::array element_face_vertices(MeshID element, - int local_face) const; + std::vector element_face_vertices(MeshID element, int local_face) const; - SurfaceElementType - get_surface_element_type(MeshID surface_element_id) const override { + SurfaceFaceType get_surface_face_type(MeshID surface_element_id) const override { // As Omega_h simplex meshes always use triangular surface elements - return SurfaceElementType::TRI; + return SurfaceFaceType::TRI; } - MeshID adjacent_element(MeshID element, int face) const override; + VolumeElementType get_volume_element_type(MeshID volume) const override { + // Omega_h simplex meshes always use tetrahedral volume elements + return VolumeElementType::TET; + } - double element_volume(MeshID element) const override; + MeshID adjacent_element(MeshID element, int face) const override; MeshID create_volume() override { return next_volume_id(); }; @@ -151,10 +152,12 @@ struct OmegaHElementFaceAccessor : public ElementFaceAccessor { MeshID element) : ElementFaceAccessor(element), mesh_manager_(mesh_manager) {} - std::array face_vertices(int i) const override { + std::vector face_vertices(int i) const override { return mesh_manager_->element_face_vertices(element_, i); } + int num_faces() const override { return 4; } + const OmegaHMeshManager *mesh_manager_; }; diff --git a/src/element_face_accessor.cpp b/src/element_face_accessor.cpp index ca1096b5..02f7c8dc 100644 --- a/src/element_face_accessor.cpp +++ b/src/element_face_accessor.cpp @@ -13,7 +13,6 @@ #include "xdg/omega_h/mesh_manager.h" #endif -#include "xdg/testing/mesh_mock.h" #include "xdg/testing/mesh_mocks.h" namespace xdg { diff --git a/src/omega_h/mesh_manager.cpp b/src/omega_h/mesh_manager.cpp index 8ec678fc..78840a90 100644 --- a/src/omega_h/mesh_manager.cpp +++ b/src/omega_h/mesh_manager.cpp @@ -338,19 +338,14 @@ std::vector OmegaHMeshManager::element_vertices(MeshID element) const { return vertices; } -std::array OmegaHMeshManager::face_vertices(MeshID face) const { - auto connectivity = face_connectivity(face); - std::array vertices; - for (int i = 0; i < VERTS_PER_TRI; ++i) { - vertices[i] = vertex_coordinates(connectivity[i]); - } - return vertices; +std::vector OmegaHMeshManager::face_vertices(MeshID face) const { + return face_connectivity(face); } -std::array +std::vector OmegaHMeshManager::element_face_vertices(MeshID element, int local_face) const { auto connectivity = element_connectivity(element); - std::array vertices; + std::vector vertices(VERTS_PER_TRI); for (int i = 0; i < VERTS_PER_TRI; ++i) { vertices[i] = vertex_coordinates(connectivity[kLocalFaceVerts[local_face][i]]); @@ -375,12 +370,6 @@ MeshID OmegaHMeshManager::adjacent_element(MeshID element, int face) const { return ID_NONE; // boundary face, no neighbor } -double OmegaHMeshManager::element_volume(MeshID element) const { - auto vertices = element_vertices(element); - std::array tet{vertices[0], vertices[1], vertices[2], vertices[3]}; - return tetrahedron_volume(tet); -} - std::pair OmegaHMeshManager::surface_senses(MeshID surface) const { auto it = surface_senses_.find(surface); diff --git a/tests/test_omega_h.cpp b/tests/test_omega_h.cpp index aaaa020f..0e46f2ec 100644 --- a/tests/test_omega_h.cpp +++ b/tests/test_omega_h.cpp @@ -78,8 +78,8 @@ TEST_CASE("Omega_h Element Types") { // Omega_h simplex meshes use triangular surface elements throughout for (const auto surface : mesh_manager->surfaces()) { - REQUIRE(mesh_manager->get_surface_element_type(surface) == - SurfaceElementType::TRI); + REQUIRE(mesh_manager->get_surface_face_type(surface) == + SurfaceFaceType::TRI); } } @@ -118,7 +118,10 @@ TEST_CASE("Omega_h Connectivity") { auto conn = mesh_manager->face_connectivity(face); REQUIRE(conn.size() == 3); - auto verts = mesh_manager->face_vertices(face); + auto vertex_ids = mesh_manager->face_vertices(face); + REQUIRE(vertex_ids == conn); + + auto verts = mesh_manager->face_vertex_coordinates(face); for (int i = 0; i < 3; ++i) { coords_match(conn[i], verts[i]); } From 4da8529fd219ba768c384e1a2749024deda08a97 Mon Sep 17 00:00:00 2001 From: magnoxemo Date: Wed, 12 Aug 2026 16:19:47 -0500 Subject: [PATCH 26/26] cache seacas in CI as well --- .github/workflows/ci.yml | 31 ++++++++++++++++++++----- .github/workflows/openmc-test.yml | 38 +++++++++++++++++++++++++++++-- src/omega_h/mesh_manager.cpp | 5 +--- 3 files changed, 62 insertions(+), 12 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 383d5034..97285702 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -116,6 +116,11 @@ jobs: run: | echo "OMEGA_H_SHA"=$(cd ~/omega_h && git rev-parse HEAD) >> $GITHUB_ENV + - name: Seacas Environment Variables + if: ${{ matrix.omega_h }} + run: | + echo "SEACAS_SHA"=$(cd ~/seacas && git rev-parse HEAD) >> $GITHUB_ENV + - name: MOAB Cache if: ${{ matrix.moab }} id: moab-cache @@ -146,6 +151,16 @@ jobs: path: ~/OMEGA_H key: omega-h-${{ runner.os }}-${{ env.cache-name }}-${{ env.OMEGA_H_SHA }} + - name: Seacas Cache + if: ${{ matrix.omega_h }} + id: seacas-cache + uses: actions/cache@v3 + env: + cache-name: seacas-cache + with: + path: ~/SEACAS + key: seacas-${{ runner.os }}-${{ env.cache-name }}-${{ env.SEACAS_SHA }} + - if: ${{ matrix.moab && steps.moab-cache.outputs.cache-hit != 'true' }} name: Build MOAB run: | @@ -170,14 +185,14 @@ jobs: make -j4 sudo make install - - if: ${{ matrix.omega_h && steps.omega-h-cache.outputs.cache-hit != 'true' }} - name: Build Omega_h + - if: ${{ matrix.omega_h && steps.seacas-cache.outputs.cache-hit != 'true' }} + name: Build Seacas shell: bash run: | - # before building omega_h we need to compile seacas on which omega_h depends + # before building omega_h we need to compile seacas on which omega_h depends # for interacting with exodus mesh API cd ~/seacas - mkdir build + mkdir build cd build cmake .. \ -DCMAKE_INSTALL_PREFIX=$HOME/SEACAS \ @@ -187,9 +202,13 @@ jobs: -DSeacas_ENABLE_SEACASExodus=ON \ -DTPL_ENABLE_MPI=OFF \ -DTPL_ENABLE_Netcdf=ON - make -j4 + make -j4 make install - # now let's install omega_h + + - if: ${{ matrix.omega_h && steps.omega-h-cache.outputs.cache-hit != 'true' }} + name: Build Omega_h + shell: bash + run: | cd ~/omega_h mkdir build cd build diff --git a/.github/workflows/openmc-test.yml b/.github/workflows/openmc-test.yml index aafea7d1..234a8737 100644 --- a/.github/workflows/openmc-test.yml +++ b/.github/workflows/openmc-test.yml @@ -57,12 +57,15 @@ jobs: git clone https://github.com/SCOREC/omega_h.git cd omega_h git checkout master + cd ~ + git clone https://github.com/sandialabs/seacas.git - name: Cache Variables run: | echo "MOAB_SHA"=$(cd ~/moab && git rev-parse HEAD) >> $GITHUB_ENV echo "LIBMESH_SHA"=$(cd ~/libmesh && git rev-parse HEAD) >> $GITHUB_ENV echo "OMEGA_H_SHA"=$(cd ~/omega_h && git rev-parse HEAD) >> $GITHUB_ENV + echo "SEACAS_SHA"=$(cd ~/seacas && git rev-parse HEAD) >> $GITHUB_ENV - name: MOAB Cache id: moab-cache @@ -91,6 +94,15 @@ jobs: path: ~/OMEGA_H key: omega-h-${{ runner.os }}-${{ env.cache-name }}-${{ env.OMEGA_H_SHA }} + - name: Seacas Cache + id: seacas-cache + uses: actions/cache@v3 + env: + cache-name: seacas-cache + with: + path: ~/SEACAS + key: seacas-${{ runner.os }}-${{ env.cache-name }}-${{ env.SEACAS_SHA }} + - if: ${{ steps.moab-cache.outputs.cache-hit != 'true' }} name: Build MOAB shell: bash @@ -113,6 +125,26 @@ jobs: ../configure --prefix=$HOME/LIBMESH --enable-exodus --disable-netcdf4 --disable-eigen --disable-lapack --disable-mpi --disable-metaphysicl make -j4 make install + - if: ${{ steps.seacas-cache.outputs.cache-hit != 'true' }} + name: Build Seacas + shell: bash + run: | + # before building omega_h we need to compile seacas on which omega_h depends + # for interacting with exodus mesh API + cd ~/seacas + mkdir build + cd build + cmake .. \ + -DCMAKE_INSTALL_PREFIX=$HOME/SEACAS \ + -DCMAKE_BUILD_TYPE=Release \ + -DBUILD_SHARED_LIBS=ON \ + -DSeacas_ENABLE_ALL_PACKAGES=OFF \ + -DSeacas_ENABLE_SEACASExodus=ON \ + -DTPL_ENABLE_MPI=OFF \ + -DTPL_ENABLE_Netcdf=ON + make -j4 + make install + - if: ${{ steps.omega-h-cache.outputs.cache-hit != 'true' }} name: Build Omega_h shell: bash @@ -124,7 +156,9 @@ jobs: -DCMAKE_INSTALL_PREFIX=$HOME/OMEGA_H \ -DCMAKE_BUILD_TYPE=Release \ -DOmega_h_USE_MPI=OFF \ - -DBUILD_TESTING=OFF + -DBUILD_TESTING=OFF \ + -DSEACASExodus_DIR=$HOME/SEACAS/lib/cmake/SEACASExodus \ + -DOmega_h_USE_SEACASExodus=ON make -j4 make install @@ -135,7 +169,7 @@ jobs: cd build cmake .. \ -DCMAKE_INSTALL_PREFIX=$HOME/opt \ - -DCMAKE_PREFIX_PATH="$HOME/MOAB;$HOME/LIBMESH;$HOME/OMEGA_H" \ + -DCMAKE_PREFIX_PATH="$HOME/MOAB;$HOME/LIBMESH;$HOME/OMEGA_H;$HOME/SEACAS" \ -DXDG_ENABLE_MOAB=ON \ -DXDG_ENABLE_LIBMESH=ON \ -DXDG_ENABLE_OMEGA_H=ON diff --git a/src/omega_h/mesh_manager.cpp b/src/omega_h/mesh_manager.cpp index 78840a90..ea5e0e41 100644 --- a/src/omega_h/mesh_manager.cpp +++ b/src/omega_h/mesh_manager.cpp @@ -195,10 +195,7 @@ void OmegaHMeshManager::determine_surface_senses() { // Only assign each slot once: a well-formed surface's every face // agrees on which side is outward, so this is a no-op after the first - // face. Some meshes carry degenerate faces classified onto a surface - // that borders more than two distinct volumes (not a true 2-manifold - // interface); for those, keep whichever pair of volumes was recorded - // first rather than letting later faces keep overwriting the senses. + // face. if (is_outward) { if (senses.first == ID_NONE) senses.first = volume;