diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 1ac4d8a2..97285702 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,18 @@ 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 + cd ~ + git clone https://github.com/sandialabs/seacas.git + + - name: OpenMP Environment Variables run: | echo "OMP_NUM_THREADS=1" >> $GITHUB_ENV @@ -97,6 +111,16 @@ jobs: # 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: 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 @@ -117,6 +141,26 @@ 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 }} + + - 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: | @@ -141,6 +185,43 @@ jobs: make -j4 sudo make install + - 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 + # 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: ${{ 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 \ + -DSEACASExodus_DIR=$HOME/SEACAS/lib/cmake/SEACASExodus \ + -DOmega_h_USE_SEACASExodus=ON + make -j4 + make install + - name: Build shell: bash run: | @@ -157,7 +238,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;$HOME/SEACAS" + 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 diff --git a/.github/workflows/openmc-test.yml b/.github/workflows/openmc-test.yml index f232e34e..234a8737 100644 --- a/.github/workflows/openmc-test.yml +++ b/.github/workflows/openmc-test.yml @@ -50,10 +50,22 @@ 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 + 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 @@ -73,6 +85,24 @@ 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 }} + + - 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 @@ -95,13 +125,54 @@ 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 + 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 \ + -DSEACASExodus_DIR=$HOME/SEACAS/lib/cmake/SEACASExodus \ + -DOmega_h_USE_SEACASExodus=ON + 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;$HOME/SEACAS" \ + -DXDG_ENABLE_MOAB=ON \ + -DXDG_ENABLE_LIBMESH=ON \ + -DXDG_ENABLE_OMEGA_H=ON make -j4 all install - name: Build OpenMC 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/CMakeLists.txt b/CMakeLists.txt index 0c185556..8d9ad0e9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -4,6 +4,7 @@ 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) 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" 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) @@ -43,6 +44,20 @@ if (NOT MOAB_USE_HDF5) endif() endif() +#=============================================================================== +# Omega_h backend +#=============================================================================== +if (XDG_ENABLE_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 @@ -148,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 OFF CACHE BOOL "Generate the fmt install target") @@ -165,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) +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() # Ensure at least one ray tracing backend is enabled @@ -240,6 +254,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) #=============================================================================== @@ -318,6 +338,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) @@ -358,6 +384,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 index 4953a2a5..40feda29 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,25 +44,59 @@ 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 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 + +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 .. \ + -DCMAKE_INSTALL_PREFIX=/XDG_TEST_SYSTEM/omega_h_install_dir \ + -DCMAKE_BUILD_TYPE=Release \ + -DOmega_h_USE_MPI=OFF \ + -DBUILD_TESTING=OFF \ + -DSEACASExodus_DIR=${SEACAS_INSTALL_PATH}/lib/cmake/SEACASExodus\ + -DOmega_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 + 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) + -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 ENV XDG_INSTALL_PATH=/XDG_TEST_SYSTEM/xdg_install_dir diff --git a/include/xdg/constants.h b/include/xdg/constants.h index caf43fd3..dacd61ce 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 = @@ -70,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/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/include/xdg/omega_h/mesh_manager.h b/include/xdg/omega_h/mesh_manager.h new file mode 100644 index 00000000..bd7b1f50 --- /dev/null +++ b/include/xdg/omega_h/mesh_manager.h @@ -0,0 +1,166 @@ +#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_mesh.hpp" + +namespace xdg { + +//! \brief Mesh manager backed by the Omega_h simplex mesh library. +class OmegaHMeshManager : public MeshManager { +public: + OmegaHMeshManager(); + + ~OmegaHMeshManager() override = default; + + 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: + 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_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::vector face_vertices(MeshID face) const override; + + //! \brief Vertices of a local face of a tetrahedral element + std::vector element_face_vertices(MeshID element, int local_face) const; + + SurfaceFaceType get_surface_face_type(MeshID surface_element_id) const override { + // As Omega_h simplex meshes always use triangular surface elements + return SurfaceFaceType::TRI; + } + + VolumeElementType get_volume_element_type(MeshID volume) const override { + // Omega_h simplex meshes always use tetrahedral volume elements + return VolumeElementType::TET; + } + + MeshID adjacent_element(MeshID element, int face) const override; + + MeshID create_volume() override { return next_volume_id(); }; + + 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; + + // 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(); + + //! \brief Fetch Omega_h's derived adjacencies and coordinate array exactly + //! once, single-threaded, and store them as members. + void cache_derived_arrays(); + + std::unique_ptr mesh_; + + //! Cached derived arrays/adjacencies + 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_; + + 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::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_; +}; + +} // namespace xdg + +#endif // XDG_OMEGA_H_MESH_MANAGER_H \ No newline at end of file diff --git a/src/element_face_accessor.cpp b/src/element_face_accessor.cpp index d9f61ac2..02f7c8dc 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_mocks.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) { if (const auto* tri_tet_mesh = dynamic_cast(mesh_manager)) { diff --git a/src/embree/ray_tracer.cpp b/src/embree/ray_tracer.cpp index d276212c..51c91bb2 100644 --- a/src/embree/ray_tracer.cpp +++ b/src/embree/ray_tracer.cpp @@ -56,6 +56,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 new file mode 100644 index 00000000..78840a90 --- /dev/null +++ b/src/omega_h/mesh_manager.cpp @@ -0,0 +1,399 @@ +#include "xdg/omega_h/mesh_manager.h" + +#include +#include +#include +#include +#include + +#include "xdg/error.h" +#include "xdg/geometry/measure.h" + +#include "Omega_h_file.hpp" + +namespace xdg { + +// Topological entity dimensions in Omega_h +static constexpr int OMEGA_H_FACE = 2; +static constexpr int OMEGA_H_REGION = 3; + +// 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; + +static constexpr int kLocalFaceVerts[FACES_PER_TET][VERTS_PER_TRI] = { + {0, 2, 1}, {0, 1, 3}, {1, 2, 3}, {2, 0, 3}}; + +// Constructor + +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. +Omega_h::Library &shared_omega_h_library() { + static Omega_h::Library lib(nullptr, nullptr); + return lib; +} +} // namespace + +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::init() { + // XDG operates on 3-dimensional volume meshes + if (mesh()->dim() != 3) { + fatal_error("Mesh must be 3-dimensional"); + } + + // in Omega_h the regions are 3D simplices or aka elements + num_elements_ = mesh_->nregions(); + + // 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. + if (has_classification()) { + discover_geometry(); + determine_surface_senses(); + } else { + discover_single_volume(); + } + + // create an implicit complement to bound the model + create_implicit_complement(); + + // map ID spaces into indices for ordered access by downstream applications + map_id_spaces(); +} + +void OmegaHMeshManager::cache_derived_arrays() { + 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) && + mesh_->has_tag(OMEGA_H_FACE, CLASS_DIM_TAG) && + mesh_->has_tag(OMEGA_H_FACE, CLASS_ID_TAG); +} + +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()); +} + +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::determine_surface_senses() { + + 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}; + } + + 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 f_base = face * VERTS_PER_TRI; + std::array natural_order = { + 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); + for (Omega_h::LO k = begin; k < end; ++k) { + 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; + } + } + } +} + +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); + + std::vector vertex_ids(mesh_->nverts()); + std::iota(vertex_ids.begin(), vertex_ids.end(), 0); + vertex_id_map_ = IDBlockMapping(vertex_ids); +} + +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 mesh_->nverts(); } + +std::vector +OmegaHMeshManager::get_volume_elements(MeshID volume) const { + 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 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 elements; +} + +std::vector OmegaHMeshManager::get_surface_faces(MeshID surface) const { + 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)) { + 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 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 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 faces; +} + +std::vector +OmegaHMeshManager::element_connectivity(MeshID element) const { + // tetrahedra store four vertices per element in element->vertex order + 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 { + // triangles store three vertices per face in face->vertex order + 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 + 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) const { + // 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)}; +} + +std::vector OmegaHMeshManager::element_vertices(MeshID element) const { + 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::vector OmegaHMeshManager::face_vertices(MeshID face) const { + return face_connectivity(face); +} + +std::vector +OmegaHMeshManager::element_face_vertices(MeshID element, int local_face) const { + auto connectivity = element_connectivity(element); + std::vector vertices(VERTS_PER_TRI); + 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 + 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 + 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)); + if (neighbor != element) { + return neighbor; + } + } + return ID_NONE; // boundary face, no neighbor +} + +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 { + // walk the surface senses and return the surfaces bounding this volume + std::vector result; + 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); + return volume == senses.first ? Sense::FORWARD : Sense::REVERSE; +} + +} // namespace xdg \ No newline at end of file diff --git a/src/xdg.cpp b/src/xdg.cpp index 75face29..15c5c5bb 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); }; diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 23b95f1d..d80bdbde 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -43,6 +43,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) @@ -55,6 +59,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..0e46f2ec --- /dev/null +++ b/tests/test_omega_h.cpp @@ -0,0 +1,323 @@ +#include +#include +#include +#include + +#include +#include +#include + +#include "util.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; + + + +TEST_CASE("Test Omega_h Initialization") { + std::unique_ptr mesh_manager = + std::make_unique(); + + mesh_manager->load_file("pincell-implicit.exo"); + mesh_manager->init(); + + 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); + + // 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("pincell-implicit.exo"); + 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("pincell-implicit.exo"); + mesh_manager->init(); + + // Omega_h simplex meshes use triangular surface elements throughout + for (const auto surface : mesh_manager->surfaces()) { + REQUIRE(mesh_manager->get_surface_face_type(surface) == + SurfaceFaceType::TRI); + } +} + +TEST_CASE("Omega_h Connectivity") { + std::unique_ptr mesh_manager = + std::make_unique(); + mesh_manager->load_file("pincell-implicit.exo"); + 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 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]); + } + + // 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("pincell-implicit.exo"); + mesh_manager->init(); + + constexpr unsigned 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) { + constexpr auto rt_backend = TestType::value; + + std::shared_ptr mesh_manager = + std::make_shared(); + 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 + auto rti = create_raytracer(rt_backend); + + 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); + } + + // 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); + } +} + +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("pincell-implicit.exo"); + mm->init(); + + xdg->prepare_raytracer(); + MeshID volume = 6; + + Position origin{0.0, 0.0, 0.0}; + Direction dir{1.0, 0.0, 0.0}; + + auto hit = xdg->ray_fire(volume, origin, dir); + REQUIRE(hit.second != ID_NONE); + 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(22.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) { + // 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)) { + 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(); + xdg->prepare_raytracer(); + + 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 + 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") +{ + // 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-sidesets.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(); + 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("pincell-implicit.exo"); + mesh_manager->init(); + xdg->prepare_raytracer(); + + 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 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(52.5, 1e-6)); +} \ No newline at end of file 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; }