diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 1ac4d8a2..4a3065e2 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 }}; mfem=${{ matrix.mfem }} runs-on: ubuntu-latest strategy: fail-fast: false matrix: libmesh: [true, false] moab: [true, false] + mfem: [true, false] exclude: - moab: false libmesh: false + mfem: false steps: - name: Checkout @@ -81,6 +83,15 @@ jobs: cd libmesh git checkout v1.7.0 + - name: mfem Clone + if: ${{ matrix.mfem }} + shell: bash + run: | + cd ~ + git clone https://github.com/mfem/mfem.git + cd mfem + git checkout v4.9 + - name: OpenMP Environment Variables run: | echo "OMP_NUM_THREADS=1" >> $GITHUB_ENV @@ -97,6 +108,11 @@ jobs: # Enforce that we're using the debug build of libMesh echo "METHOD=dbg" >> $GITHUB_ENV + - name: mfem Environment Variables + if: ${{ matrix.mfem }} + run: | + echo "MFEM_SHA"=$(cd ~/mfem && git rev-parse HEAD) >> $GITHUB_ENV + - name: MOAB Cache if: ${{ matrix.moab }} id: moab-cache @@ -117,6 +133,16 @@ jobs: path: ~/LIBMESH key: libmesh-${{ runner.os }}-${{ env.cache-name }}-${{ env.LIBMESH_SHA }} + # - name: mfem Cache + # if: ${{ matrix.mfem }} + # id: mfem-cache + # uses: actions/cache@v3 + # env: + # cache-name: mfem-cache + # with: + # path: ~/mfem + # key: mfem-${{ runner.os }}-${{ env.cache-name }}-${{ env.MFEM_SHA }} + - if: ${{ matrix.moab && steps.moab-cache.outputs.cache-hit != 'true' }} name: Build MOAB run: | @@ -141,6 +167,23 @@ jobs: make -j4 sudo make install + - if: ${{ matrix.mfem && steps.mfem-cache.outputs.cache-hit != 'true' }} + name: Build mfem + shell: bash + run: | + # first let's get netcdf + sudo apt-get update + sudo apt-get install -y libnetcdf-dev libhdf5-dev + cd ~ + cd mfem + git checkout v4.9 + mkdir build + cd build + # manually entering the hdf5 dir since cmake can't find it on its own + cmake .. -DMFEM_USE_HYPRE=OFF -DMFEM_USE_NETCDF=ON -DCMAKE_INSTALL_PREFIX=../install -DCMAKE_BUILD_TYPE=Debug -DHDF5_DIR=/usr/lib/x86_64-linux-gnu/hdf5/serial -DBUILD_SHARED_LIBS=ON -DCMAKE_POSITION_INDEPENDENT_CODE=ON + make -j4 + sudo make install + - name: Build shell: bash run: | @@ -157,7 +200,13 @@ 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.mfem }}" = 'true' ]; then + if [ -n "$CMAKE_PREFIX_PATH" ]; then + CMAKE_PREFIX_PATH="$CMAKE_PREFIX_PATH;" + fi + CMAKE_PREFIX_PATH="${CMAKE_PREFIX_PATH}$HOME/mfem/install" + 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_MFEM=${{ matrix.mfem && 'ON' || 'OFF' }} make -j4 install - name: Test diff --git a/.gitignore b/.gitignore index 17800a19..ba5458c9 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,7 @@ bld -build +build* .vscode docs/_build docs/Doxyfile .vscode +*.jou diff --git a/CMakeLists.txt b/CMakeLists.txt index 0c185556..8a09030f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -43,6 +43,17 @@ if (NOT MOAB_USE_HDF5) endif() endif() +#=============================================================================== +# MFEM +#=============================================================================== +if (XDG_ENABLE_MFEM) + find_package(MFEM REQUIRED HINTS ${MFEM_DIR}) +if (NOT MFEM_FOUND) + message(FATAL_ERROR "MFEM package was not found") + endif() + message(STATUS "Found MFEM ${MFEM_VERSION} at ${MFEM_DIR}") +endif() + if (XDG_ENABLE_EMBREE) # find Embree for CPU ray tracing @@ -167,10 +178,11 @@ 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_MFEM) message(FATAL_ERROR "No mesh backend enabled. Enable at least one of:\n" " -DXDG_ENABLE_MOAB=ON\n" + " -DXDG_ENABLE_MFEM=ON\n" " -DXDG_ENABLE_LIBMESH=ON") endif() @@ -240,6 +252,13 @@ src/moab/metadata.cpp ) endif() +if (XDG_ENABLE_MFEM) +list(APPEND xdg_sources +# MFEM +src/mfem/mesh_manager.cpp +) +endif() + #=============================================================================== # RPATH information (from OpenMC) #=============================================================================== @@ -314,6 +333,10 @@ if (XDG_ENABLE_LIBMESH) target_compile_definitions(xdg PUBLIC XDG_ENABLE_LIBMESH) endif() +if (XDG_ENABLE_MFEM) + target_compile_definitions(xdg PUBLIC XDG_ENABLE_MFEM) +endif() + if (XDG_ENABLE_EMBREE) target_compile_definitions(xdg PUBLIC XDG_ENABLE_EMBREE) endif() @@ -369,6 +392,10 @@ if (XDG_ENABLE_MOAB) target_link_libraries(xdg PRIVATE MOAB) endif() +if (XDG_ENABLE_MFEM) + target_link_libraries(xdg PUBLIC mfem) +endif() + #================================================================= # Installation & Packaging #================================================================= diff --git a/include/xdg/bbox.h b/include/xdg/bbox.h index 170faea1..34b545ec 100644 --- a/include/xdg/bbox.h +++ b/include/xdg/bbox.h @@ -33,6 +33,15 @@ bool operator ==(const BoundingBox& other) { max_z == other.max_z; } +void update(const double* v) { + min_x = std::min(min_x, v[0]); + min_y = std::min(min_y, v[1]); + min_z = std::min(min_z, v[2]); + max_x = std::max(max_x, v[0]); + max_y = std::max(max_y, v[1]); + max_z = std::max(max_z, v[2]); +} + void update(const Vertex& v) { min_x = std::min(min_x, v.x); min_y = std::min(min_y, v.y); diff --git a/include/xdg/constants.h b/include/xdg/constants.h index caf43fd3..abb590b4 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, + MFEM }; // 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::MFEM, "MFEM"} }; static const std::map RT_LIB_TO_STR = diff --git a/include/xdg/libmesh/mesh_manager.h b/include/xdg/libmesh/mesh_manager.h index 5a869b7b..2356c7de 100644 --- a/include/xdg/libmesh/mesh_manager.h +++ b/include/xdg/libmesh/mesh_manager.h @@ -106,7 +106,7 @@ class LibMeshManager : public MeshManager { int num_vertices() const override; - std::vector get_volume_elements(MeshID volume) const; + std::vector get_volume_elements(MeshID volume) const override; std::vector get_surface_faces(MeshID surface) const override; diff --git a/include/xdg/mesh_managers.h b/include/xdg/mesh_managers.h index fe6912ce..46b4853c 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_MFEM +#include "xdg/mfem/mesh_manager.h" +#endif \ No newline at end of file diff --git a/include/xdg/mfem/mesh_manager.h b/include/xdg/mfem/mesh_manager.h new file mode 100644 index 00000000..b9704954 --- /dev/null +++ b/include/xdg/mfem/mesh_manager.h @@ -0,0 +1,266 @@ +#ifndef _XDG_MFEM_MESH_MANAGER +#define _XDG_MFEM_MESH_MANAGER + +#include + +#include "xdg/constants.h" +#include "xdg/element_face_accessor.h" +#include "xdg/mesh_manager_interface.h" +#include "xdg/error.h" + +#include "mfem/mfem.hpp" + +namespace xdg { +class MfemMeshManager : public MeshManager { +public: + + MfemMeshManager() {}; + + ~MfemMeshManager() override = default; + + // Backend methods + + void load_file(const std::string &filepath) override; + + void init() override; + + // Accessors + const mfem::Mesh* mesh() const { return mesh_.get(); } + mfem::Mesh* mesh() { return mesh_.get(); } + + // Interface methods + MeshLibrary mesh_library() const override { return MeshLibrary::MFEM; } + + 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; + } + } + + // I think we need to count the number of elements with the attribute "volume" + int num_volume_elements(MeshID volume) const override { + fatal_error("MfemMeshManager::num_volume_elements() not implemented yet"); + } + + int num_volume_elements() const override { + return mesh_->GetNE(); + } + + int num_boundary_elements() const { + return mesh_->GetNBE(); + } + + // count the number of faces with the attribute "volume" + int num_volume_faces(MeshID volume) const override { + fatal_error("MfemMeshManager::num_volume_faces() not implemented yet"); + } + + int num_surface_faces(MeshID surface) const override { + fatal_error("MfemMeshManager::num_surface_faces() not implemented yet"); + } + + // get all of the elements in this volume + virtual std::vector get_volume_elements(MeshID volume) const override; + + virtual std::vector get_surface_faces(MeshID surface) const override; + + virtual std::vector element_vertices(MeshID element) const override; + std::vector bdr_element_vertices(MeshID element) const; + + virtual std::vector face_vertices(MeshID element) const override; + + SurfaceFaceType get_surface_face_type(MeshID surface) const override; + + VolumeElementType get_volume_element_type(MeshID volume) const override; + + // The table works wonders for this + virtual MeshID adjacent_element(MeshID element, int face) const override; + + virtual MeshID get_boundary_face_element(MeshID face) const override { + fatal_error("MfemMeshManager::get_boundary_face_element() not implemented yet"); + } + + virtual Sense surface_sense(MeshID surface, MeshID volume) const override { + fatal_error("MfemMeshManager::surface_sense() not implemented yet"); + } + + virtual int num_vertices() const override { + return mesh_->GetNV(); + } + + virtual double element_volume(MeshID element) const override { + return mesh_->GetElementVolume(element); + } + + virtual Vertex vertex_coordinates(MeshID vertex_id) const override { + Vertex output; + const mfem::real_t* vertices = mesh_->GetVertex(vertex_id); + for (int i=0; iDimension(); i++) output[i] = vertices[i]; + return output; + } + + virtual std::vector element_connectivity(MeshID element) const override { + fatal_error("MfemMeshManager::element_connectivity() not implemented yet"); + } + + virtual std::vector face_connectivity(MeshID face) const override { + fatal_error("MfemMeshManager::face_connectivity() not implemented yet"); + } + + // Topology + + std::vector get_volume_surfaces(MeshID volume) const override; + + std::pair surface_senses(MeshID surface) const override; + + // Seems like it's only used to create the implicit complement + MeshID create_volume() override; + + void add_surface_to_volume(MeshID volume, MeshID surface, Sense sense, bool overwrite=false) override; + + // This will be largely a copy of the libmesh version + void determine_surface_senses(); + + // Metadata methods + void parse_metadata() override; + + // Accessors + const std::unique_ptr& mfem_mesh() const { + return mesh_; + } + + // Data members +private: + std::unique_ptr mesh_ {nullptr}; + + // For each volume of the mesh, keep a set of the interior element IDs + std::map> volume_to_element_map_; + + // For each sideset of the mesh, keep a set of the boundary element IDs + // This the libmesh equivalent of surface_map_ + std::map> sideset_to_element_map_; + + // map to keep track of each sideset held by a particular + // volume + std::map> volumes_to_sidesets_; + + // set to capture all of the valid volumes/attributes + // It's a set (not vector) to prevent double counting + // Note that this is somewhat redundant, since the base + // class has the volumes_ vector, which we copy the + // contents of this into. We keep BOTH, with attributes_ + // meant to resemble the TRUE volumes that the mfem mesh + // recogonises. volumes_ will be amended to include the + // implicit complement as well... + std::set attributes_; + + int num_interior_faces_; + int num_boundary_faces_; + + //! Mapping of surfaces to the volumes on either side. Volumes are ordered + //! based on their sense with respect to the surface triangles. We reuase + //! whichever ordering mfem decides on when the mesh is constructed. + std::unordered_map> surface_senses_; +}; + +struct MfemMeshElementFaceAccessor : public ElementFaceAccessor { + MfemMeshElementFaceAccessor(const MfemMeshManager* mesh_manager, MeshID element) : + ElementFaceAccessor(element), mesh_manager_(mesh_manager), element_(element) { + + // for each face, fetch the vertices of the face + auto& mesh = mesh_manager_->mfem_mesh(); + + // TODO: circumvent this method; fetch the table directly + mfem::Array ori; // don't care abour ori + mesh->GetElementFaces(element_, faces_, ori); + + // TODO: Fix hardcoded 4 faces + for (int f=0; f<4; f++) { + int face_no = faces_[f]; + + // pointer to the element object that defines this face + auto face_obj = mesh->GetFace(face_no); + face_obj->GetVertices(vertex_indices_); + + for (int v=0; vGetVertex( vertex_indices_[v] ); + for (int d=0; dSpaceDimension(); d++) + face_vertices_[f][v][d] = vertices[d]; + } + } + } + + // TODO: Is this correct? We're getting interior vertices here + // + // Clarification: we are concerned with the face here. The moab + // stores the vertices of the element, and picks the correct three + // that correspond to this face. Why not just get the face from + // the mesh itself? It exposes the vertices + std::vector face_vertices(int i) const override { + std::vector output; + + // we have already gathered the vertices for this face. + // copy them into the output array + std::copy(face_vertices_[i], face_vertices_[i+1], std::back_inserter(output)); + + // we need mesh_->GetFaceElementTransformations + auto& mesh = mesh_manager_->mfem_mesh(); + + // we get the face element trafo for the face we are + // currently talking about + int face_no = faces_[i]; + auto face_el_tx = mesh->GetFaceElementTransformations(face_no); + + // The mfem docu implies that FaceElementTransformations::Elem1No + // is the one that the normal vector is supposed to point out of. + // We use this info to switch some stuff around + + if ( face_el_tx->Elem2No == element_ ) + // This element is NOT the one that the normal vector points out of. + // switch two of the vertices around to make sure the cross product is good. + std::swap( output[0], output[1] ); + + return output; + } + + int num_faces() const override { + // first, look up the geom type + auto geom = mesh_manager_->get_volume_element_type(element_); + + switch (geom) { + case VolumeElementType::TET: return 4; + case VolumeElementType::HEX: return 6; + default: fatal_error(""); + } + } + + // data members + const MfemMeshManager* mesh_manager_; + + // 4 faces, 3 vertices each + Vertex face_vertices_[4][3]; + // indices for each of the faces on this element + mfem::Array faces_; + mfem::Array vertex_indices_; + + MeshID element_; +}; + +// helper functions to convert mfem's element types to xdg +VolumeElementType GetVolumeElementTypeFromMfem( mfem::Element::Type t ); +// SurfaceElementType GetSurfaceElementTypeFromMfem( mfem::Element::Type t ); + +} // namespace xdg + +#endif // include guard \ No newline at end of file diff --git a/include/xdg/moab/mesh_manager.h b/include/xdg/moab/mesh_manager.h index 1593ddd6..4062f296 100644 --- a/include/xdg/moab/mesh_manager.h +++ b/include/xdg/moab/mesh_manager.h @@ -38,7 +38,7 @@ class MOABMeshManager : public MeshManager { // Interface Methods MeshLibrary mesh_library() const override { return MeshLibrary::MOAB; } - void load_file(const std::string& filepath); + void load_file(const std::string& filepath) override; void init() override; diff --git a/include/xdg/moab/tag_conventions.h b/include/xdg/moab/tag_conventions.h index 5ad75a83..e42c6496 100644 --- a/include/xdg/moab/tag_conventions.h +++ b/include/xdg/moab/tag_conventions.h @@ -1,6 +1,6 @@ // Borrowed from MOAB #ifndef _XDG_MOABTAG_CONVENTIONS_H -#define XDG_MOAB_TAG_CONVENTIONS_H +#define _XDG_MOABTAG_CONVENTIONS_H namespace xdg { diff --git a/src/element_face_accessor.cpp b/src/element_face_accessor.cpp index d9f61ac2..2f990dc5 100644 --- a/src/element_face_accessor.cpp +++ b/src/element_face_accessor.cpp @@ -10,6 +10,9 @@ #endif #include "xdg/testing/mesh_mocks.h" +#ifdef XDG_ENABLE_MFEM +#include "xdg/mfem/mesh_manager.h" +#endif namespace xdg { @@ -26,6 +29,12 @@ std::shared_ptr ElementFaceAccessor::create(const MeshManag return std::make_shared(libmesh_mesh_manager, element); } #endif + #ifdef XDG_ENABLE_MFEM + if (mesh_manager->mesh_library() == MeshLibrary::MFEM) { + const MfemMeshManager* mfem_mesh_manager = dynamic_cast(mesh_manager); + return std::make_shared(mfem_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..6154323c 100644 --- a/src/embree/ray_tracer.cpp +++ b/src/embree/ray_tracer.cpp @@ -313,7 +313,7 @@ EmbreeRayTracer::ray_fire(SurfaceTreeID tree, const Direction& direction, const double dist_limit, HitOrientation orientation, - std::vector* const exclude_primitves) + std::vector* const exclude_primitives) { RTCScene scene = surface_volume_tree_to_scene_map_.at(tree); RTCDualRayHit rayhit; @@ -327,7 +327,7 @@ EmbreeRayTracer::ray_fire(SurfaceTreeID tree, rayhit.ray.mask = -1; // no mask rayhit.ray.volume_tree = tree; - if (exclude_primitves != nullptr) rayhit.ray.exclude_primitives = exclude_primitves; + if (exclude_primitives != nullptr) rayhit.ray.exclude_primitives = exclude_primitives; // fire the ray { @@ -342,7 +342,7 @@ EmbreeRayTracer::ray_fire(SurfaceTreeID tree, return {INFTY, ID_NONE}; else - if (exclude_primitves) exclude_primitves->push_back(rayhit.hit.primitive_ref->primitive_id); + if (exclude_primitives) exclude_primitives->push_back(rayhit.hit.primitive_ref->primitive_id); return {rayhit.ray.dtfar, rayhit.hit.surface}; } diff --git a/src/mfem/mesh_manager.cpp b/src/mfem/mesh_manager.cpp new file mode 100644 index 00000000..50c03fc0 --- /dev/null +++ b/src/mfem/mesh_manager.cpp @@ -0,0 +1,437 @@ +#include + +#include "xdg/mfem/mesh_manager.h" +#include "xdg/util/str_utils.h" + +namespace xdg { +void MfemMeshManager::load_file(const std::string &filepath) { + mesh_ = std::make_unique(filepath.c_str(), 1, 1); +} + +void MfemMeshManager::init() { + // Ensure the mesh is 3-dimensional + if (mesh_->Dimension() != 3) { + fatal_error("Mesh must be 3-dimensional"); + } + + // this is done in the mesh reader + // mesh_->FinalizeTopology(); + + // set the volumes/attributes... + // This set should have one entry per volume/attribute type + for (int i=0; iattributes.Size(); i++) { + attributes_.insert( mesh_->attributes[i] ); + } + + // Create a set for each volume attribute. Gather the IDs of all the + // interior elements with this characteristic + // TODO: This won't work with ParMesh + // + // While we are here, we can lazily populate the volume_element_id_map_ + // To do that, we need a vector to store the element IDs into. This + // is very lazy, and no better than just filling a vector with sequential + // ints anyway. + std::vector volume_element_ids(mesh_->GetNE()); + for (int i=0; iGetNE(); i++) { + int volume_id = mesh_->GetAttribute(i); + volume_to_element_map_[volume_id].insert(i); + volume_element_ids[i] = i; + } + + // Finish the BlockMapping for volume elements + volume_element_id_map_ = IDBlockMapping(volume_element_ids); + + // create a set for capturing all of the sideset IDs + // without repeats + std::set sideset_ids; + + // same for boundary attributes + for (int i=0; iGetNBE(); i++) { + int sideset = mesh_->GetBdrAttribute(i); + + sideset_to_element_map_[sideset].insert(i); + + // We also want to count the number of sidesets that each volume has. + // So, while we are looping over each boundary element, we look at + // which sideset it's a member of. We then look at + // its immediate neighbour on the interior of the mesh. We query + // this neighbour for which volume it's a member of, and register + // this sideset as a member of that volume. + // We want each sideset to be a member of exactly one volume, but + // that's probably too much to ask. + int elem_no, info; + mesh_->GetBdrElementAdjacentElement(i, elem_no, info); + + int volume = mesh_->GetAttribute(elem_no); + volumes_to_sidesets_[volume].insert(sideset); + + // we want to populate the surfaces_ array from the base class + sideset_ids.insert(sideset); + } + + /* + // we wanna check how many implicit boundaries we detect. So let's + // create a set to collect them + std::set< std::pair > implicit_bdr; + for (int f=0; fGetNumFaces(); f++) { + int e1, e2; + mesh_->GetFaceElements(f, &e1, &e2); + + // if the el index is -1, then just set the volume id to -1. + // this means the outside + int vol1 = (e1==-1) ? -1 : mesh_->GetAttribute(e1); + int vol2 = (e2==-1) ? -1 : mesh_->GetAttribute(e2); + + if (vol1!=vol2) { + // we have found implicit boundary. add to list! + if (vol1first << " and " << iter->second << "\n"; + } + */ + + // We've read in the mesh and counted all the attributes, i.e. a unique + // list of all the attributes we've seen. Let's copy the contents of + // attributes_ into volumes_, so the base class has access to the list + // of volume IDs + std::copy(attributes_.begin(), attributes_.end(), std::back_inserter(volumes_)); + + // ditto for surfaces + std::copy(sideset_ids.begin(), sideset_ids.end(), std::back_inserter(surfaces_)); + + // set these two attributes related to interior/boundary faces + num_interior_faces_ = mesh_->GetNumFaces(); + num_boundary_faces_ = mesh_->GetNBE(); + + determine_surface_senses(); + + create_implicit_complement(); +} + +// TODO: very slow, and could be done during init() +std::vector MfemMeshManager::get_volume_elements(MeshID volume) const { + std::vector output; + + // Copy the way that libmesh does it, which is if we call + // this function on the implicit complement, then return an + // empty vector. + // Note the deliberate use of attributes_ (which comes from + // the mfem mesh, naming convention intact) and not volumes_. + // This is because we would have modified volumes_ to include + // an implicit complement. More concretetly, attributes_ + // represents the TRUE volumes that the mfem mesh recognises. + if (attributes_.find(volume) == attributes_.end()) { + // check that we are looking at the implcit complement + if (std::find(volumes_.begin(), volumes_.end(), volume) == volumes_.end()) { + // This is an error now. It's not a true volume + // or the implicit complement + std::ostringstream output; + output << "Couldn't find volume " << volume << "\n"; + fatal_error(output.str()); + } + + // simply return the empty vector if this is the implicit + // complement + return output; + } + + // gather all the element IDs that have this attribute + // this method is absolutely criminal. Could be done at the start + // when we run over all the elements anyway... + for (int i=0; iGetNE(); i++) { + if ( mesh_->GetAttribute(i) == volume ) output.push_back(i); + } + + return output; +} + +// Should return all of the sidesets that are a part of this volume +std::vector MfemMeshManager::get_volume_surfaces(MeshID volume) const { + // walk the surface senses and return the surfaces that have this volume + // as an entry + std::vector surfaces; + for (const auto& [surface, senses] : surface_senses_) { + if (senses.first == volume || senses.second == volume) { + surfaces.push_back(surface); + } + } + return surfaces; +} + +std::vector MfemMeshManager::get_surface_faces(MeshID surface) const { + // get the set associated with this surface + const std::set& boundary_faces = sideset_to_element_map_.at(surface); + + // copy it into a vector + std::vector output(boundary_faces.begin(), boundary_faces.end()); + + // We want [0->mesh_->GetNumFaces() ) to represent interior faces. + // and we want [ mesh_->GetNumFaces(), mesh_->GetNumFaces() + mesh_->GetNBE() ) + // to represent the boundary faces. + // When we query the face vertices later, we need to take this + // mapping into account. All we do here is increase the MeshIDs by + // num_interior_faces_ to effect this mapping + std::transform( output.begin(), output.end(), output.begin(), [&](int in){ return in + num_interior_faces_; } ); + + return output; +} + +std::vector MfemMeshManager::face_vertices(MeshID element) const { + mfem::Array index_array; + + if (element >= num_interior_faces_) { + // we are actually talking about a boundary element here. this is + // us taking the mapping into account. see comments at the end of + // get_surfaces_faces + MeshID bdr_element = element - num_interior_faces_; + mesh_->GetBdrElementVertices(bdr_element, index_array); + + mfem::Element* bdr_el = mesh_->GetBdrElement(bdr_element); + int* vertices = bdr_el->GetVertices(); + } + + else { + // index_array gets populated with the indices of the vertices itself + mesh_->GetFaceVertices(element, index_array); + } + + // for (int i=0; iGetVertex( index_array[i] ); + + // for (int d=0; dSpaceDimension(); d++) output[i][d] = vertices[d]; + // } + + // copy these indices into std::vector and return it + std::vector output( index_array.GetData(), index_array.GetData() + index_array.Size() ); + return output; +} + +std::pair MfemMeshManager::surface_senses(MeshID surface) const { + + // TODO: make the second value one more than the largest volume ID we've seen + // i.e. since the only volume in the jezebel/brick is 1, the second id must be 2, + // to denote the implicit complement + return surface_senses_.at(surface); +} + +std::vector MfemMeshManager::element_vertices(MeshID element) const { + mfem::Array index_array; + + // ask the mesh for the vertices of this element + mesh_->GetElementVertices(element, index_array); + + std::vector output(index_array.Size()); + + for (int i=0; iGetVertex( index_array[i] ); + + for (int d=0; dSpaceDimension(); d++) output[i][d] = vertices[d]; + } + + return output; +} + +// I've written this extra function because the mesh manager needs to support one +// continuous list of elements for boundary and interior. So at some point we need +// to map them together. +// Update: not sure that's true. moab mesh manager reports the same number of elements +// on the jezebel as mesh->GetNE() +std::vector MfemMeshManager::bdr_element_vertices(MeshID element) const { + mfem::Array index_array; + + // ask the mesh for the vertices of this element + mesh_->GetBdrElementVertices(element, index_array); + + std::vector output(index_array.Size()); + + for (int i=0; iGetVertex( index_array[i] ); + + for (int d=0; dSpaceDimension(); d++) output[i][d] = vertices[d]; + } + + return output; +} + +MeshID MfemMeshManager::adjacent_element(MeshID element, int face) const { + // Would be nice if we had the face element accessor still available + // TODO: do something novel if we are already on the boundary + mfem::Array faces, ori; // don't care about ori + + mesh_->GetElementFaces(element, faces, ori); + + if (face==ID_NONE) return ID_NONE; + + // not quite. faces[face] is just the faceID. We need + // the element that it's connected to. + int e1, e2; + mesh_->GetFaceElements(faces[face], &e1, &e2); + + // e1 and e2 are now the element IDs of two elements + // that share this face. + assert(element == e1 or element == e2); + + if (element == e1) return e2; + else if (element == e2) return e1; + + fatal_error("Shouldn't reach this far!"); +} + +void MfemMeshManager::determine_surface_senses() { + for (auto &[surface_id, surface_faces] : sideset_to_element_map_) { + if (surface_faces.size() == 0) continue; + + int face_index = *surface_faces.begin(); + + // first, get the volume that this sideset is living on. + // we do it the dumb way + int elem_no, info; + mesh_->GetBdrElementAdjacentElement(face_index, elem_no, info); + + // TODO: does this return the same number for every element on + // sideset 3? + + int volume = mesh_->GetAttribute(elem_no); + + int face_no = mesh_->GetBdrElementFaceIndex(face_index); + // check if the connectivity is still there + int e1, e2; + mesh_->GetFaceElements(face_no, &e1, &e2); + + // if we have both elements nontrivial (i.e. != -1) then we ask the + // the mesh which is elem1 and which is elem2 + // The normal vector is supposed to point from the reverse sense + // to the forward sense + if (e1 != -1 and e2 !=-1) { + auto face_el_tx = mesh_->GetFaceElementTransformations(face_no); + + // check that face_el_tx has picked out the correct elements + assert( (face_el_tx->Elem1No == e1 or face_el_tx->Elem1No == e2) + and (face_el_tx->Elem2No == e1 or face_el_tx->Elem2No == e2) + ); + + // Elem1 is the reverse sense and Elem2 is the forwards sense, since + // by construction, this is the way the normal vectors are pointing. + // We can't put the element IDs in to the array, so we have to ask + // the mesh for their attributes + surface_senses_[surface_id] = { + mesh_->GetAttribute(face_el_tx->Elem1No), mesh_->GetAttribute(face_el_tx->Elem2No) + }; + } + + // We have a surface on a true boundary. The second element is simply + // the implcit complement + else { + surface_senses_[surface_id] = {volume, ID_NONE}; + } + } +} + +MeshID MfemMeshManager::create_volume() { + MeshID next_volume_id = *std::max_element(volumes_.begin(), volumes_.end()) + 1; + return next_volume_id; +} + +void MfemMeshManager::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_senses_[surface] = {volume, senses.second}; + } + + else { + if (!overwrite && senses.second != ID_NONE) { + fatal_error("Surface already has a reverse sense"); + } + surface_senses_[surface] = {senses.first, volume}; + } +} + +void MfemMeshManager::parse_metadata() { + auto& attr_sets = mesh_->attribute_sets; + + for (auto vol_name : attr_sets.GetAttributeSetNames()) { + const mfem::Array &attrs = attr_sets.GetAttributeSet(vol_name); + remove_substring(vol_name, "mat:"); + + // loop over every volume that is named the same thing + for (int i = 0; i < attrs.Size(); i++) { + const int vol_id = attrs[i]; + if (vol_name.empty()) + volume_metadata_[{vol_id, PropertyType::MATERIAL}] = VOID_MATERIAL; + else + volume_metadata_[{vol_id, PropertyType::MATERIAL}] = {PropertyType::MATERIAL, vol_name}; + + } + } + + // ditto for boundaries + auto& bdr_attr_sets = mesh_->bdr_attribute_sets; + for (auto bdr_name : bdr_attr_sets.GetAttributeSetNames()) { + const mfem::Array &attrs = bdr_attr_sets.GetAttributeSet(bdr_name); + remove_substring(bdr_name, "boundary:"); + + // loop over every surface that is named the same thing + for (int i = 0; i < attrs.Size(); i++) { + const int bdr_id = attrs[i]; + const auto key = std::make_pair(bdr_id, PropertyType::BOUNDARY_CONDITION); + + surface_metadata_[key] = {PropertyType::BOUNDARY_CONDITION, bdr_name}; + } + } +} + +// This is not quite correct. mfem does support mixed meshes. +// The intention of the caller is that the argument (surface) corresponds +// to the sideset. So we need to find a typical element from the sideset +// that is marked by the argument surface. For now, this will do +SurfaceFaceType MfemMeshManager::get_surface_face_type(MeshID surface) const { + mfem::Geometry::Type geom = mesh_->GetFaceGeometry(surface); + + switch(geom) { + case mfem::Geometry::TRIANGLE: return SurfaceFaceType::TRI; + case mfem::Geometry::TETRAHEDRON : return SurfaceFaceType::QUAD; + default: + fatal_error("Unsupported geom"); + } +} + +// Same problem as above. Volume is supposed to be a block_id, and this +// function is interpreting it as an element index. We should have a LOT more +// elements than blocks, so it's safe, but wrong. When dealing with a mixed +// mesh, it will fail tests +VolumeElementType MfemMeshManager::get_volume_element_type(MeshID volume) const { + mfem::Geometry::Type geom = mesh_->GetElementBaseGeometry(volume); + + switch (geom) { + case mfem::Geometry::TETRAHEDRON: return VolumeElementType::TET; + case mfem::Geometry::TRIANGLE: return VolumeElementType::HEX; + default: + fatal_error("Unsupported geom"); + } + + return VolumeElementType::TET; +} + + +// helper function to convert mfem's element types to xdg +VolumeElementType GetTypeFromMfem( mfem::Element::Type t ) { + switch (t) { + case mfem::Element::TETRAHEDRON: return VolumeElementType::TET; + case mfem::Element::HEXAHEDRON: return VolumeElementType::HEX; + default: + fatal_error("Unsupported element type\n"); + } +} + +} // namespace xdg diff --git a/src/xdg.cpp b/src/xdg.cpp index 75face29..3d1b500a 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_MFEM + if (mesh_lib == MeshLibrary::MFEM) 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_MFEM + msg += " MFEM"; + #endif fatal_error(msg); }; diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 23b95f1d..45c6b7fa 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -31,7 +31,11 @@ if (XDG_ENABLE_LIBMESH) list(APPEND TEST_NAMES test_libmesh) endif() -if (XDG_ENABLE_MOAB AND XDG_ENABLE_LIBMESH) +if (XDG_ENABLE_MFEM) + list(APPEND TEST_NAMES test_mfem) +endif() + +if (XDG_ENABLE_MOAB AND XDG_ENABLE_LIBMESH AND XDG_ENABLE_MFEM) list(APPEND TEST_NAMES test_mesh_library_cross_check) endif() @@ -43,6 +47,7 @@ if (XDG_ENABLE_MOAB AND XDG_BUILD_TOOLS) list(APPEND TEST_NAMES test_overlap_check) endif() + foreach(test ${TEST_NAMES}) add_executable(${test} test_main.cpp ${test}.cpp) target_link_libraries(${test} PRIVATE xdg fmt::fmt Catch2::Catch2) diff --git a/tests/test_mesh_library_cross_check.cpp b/tests/test_mesh_library_cross_check.cpp index 6ebb7ac2..1af17b38 100644 --- a/tests/test_mesh_library_cross_check.cpp +++ b/tests/test_mesh_library_cross_check.cpp @@ -60,20 +60,21 @@ class CrossCheck { TEST_CASE("Test MOAB-libMesh Cross-Check 1 Vol") { - auto harness = CrossCheck({{"jezebel.exo", MeshLibrary::LIBMESH}, {"jezebel.h5m", MeshLibrary::MOAB}}); + auto harness = CrossCheck({{"jezebel.exo", MeshLibrary::LIBMESH}, {"jezebel.h5m", MeshLibrary::MOAB}, {"jezebel.exo", MeshLibrary::MFEM}}); harness.transport(); harness.check(); } TEST_CASE("Test MOAB-libMesh Cross-Check 2 Vol") { - auto harness = CrossCheck({{"cyl-brick.exo", MeshLibrary::LIBMESH}, {"cyl-brick.h5m", MeshLibrary::MOAB}}); + auto harness = CrossCheck({{"cyl-brick.exo", MeshLibrary::LIBMESH}, {"cyl-brick.h5m", MeshLibrary::MOAB}, {"cyl-brick.exo", MeshLibrary::MFEM}}); harness.transport(); harness.check(); } TEST_CASE("Test MOAB-libMesh Cross-Check Pincell -- Implicit libMesh Boundaries") { + // Skip this one for mfem until we have implicit boundaries working auto harness = CrossCheck({{"pincell-implicit.exo", MeshLibrary::LIBMESH}, {"pincell.h5m", MeshLibrary::MOAB}}); harness.transport(); harness.check(); @@ -130,4 +131,92 @@ TEST_CASE("Test MOAB-libMesh Cross-Check Tallies -- Simple Cubes, Tet Mesh") REQUIRE_THAT(moab_tracks[j].second, Catch::Matchers::WithinAbs(libmesh_tracks[j].second, 1e-10)); } } -} \ No newline at end of file +} + +TEST_CASE("Test MOAB-libMesh Cross-Check Tallies -- JEZEBEL") +{ + auto xdg_moab = XDG::create(MeshLibrary::MOAB); + xdg_moab->mesh_manager()->load_file("jezebel.h5m"); + xdg_moab->mesh_manager()->init(); + xdg_moab->mesh_manager()->parse_metadata(); + xdg_moab->prepare_raytracer(); + + auto xdg_libmesh = XDG::create(MeshLibrary::LIBMESH); + xdg_libmesh->mesh_manager()->load_file("jezebel.exo"); + xdg_libmesh->mesh_manager()->init(); + xdg_libmesh->mesh_manager()->parse_metadata(); + xdg_libmesh->prepare_raytracer(); + + auto xdg_mfem = XDG::create(MeshLibrary::MFEM); + xdg_mfem->mesh_manager()->load_file("jezebel.exo"); + xdg_mfem->mesh_manager()->init(); + xdg_mfem->mesh_manager()->parse_metadata(); + xdg_mfem->prepare_raytracer(); + + // check that the global bounding box of the model and various model counts are the same + REQUIRE(xdg_moab->mesh_manager()->num_vertices() == xdg_libmesh->mesh_manager()->num_vertices()); + REQUIRE(xdg_moab->mesh_manager()->num_vertices() == xdg_mfem->mesh_manager()->num_vertices()); + + REQUIRE(xdg_moab->mesh_manager()->num_volume_elements() == xdg_libmesh->mesh_manager()->num_volume_elements()); + REQUIRE(xdg_moab->mesh_manager()->num_volume_elements() == xdg_mfem->mesh_manager()->num_volume_elements()); + + REQUIRE(xdg_moab->mesh_manager()->num_volumes() == xdg_libmesh->mesh_manager()->num_volumes()); + REQUIRE(xdg_moab->mesh_manager()->num_volumes() == xdg_mfem->mesh_manager()->num_volumes()); + + auto moab_bounding_box = xdg_moab->mesh_manager()->global_bounding_box(); + auto libmesh_bounding_box = xdg_libmesh->mesh_manager()->global_bounding_box(); + auto mfem_bounding_box = xdg_mfem->mesh_manager()->global_bounding_box(); + + REQUIRE_THAT(moab_bounding_box.min_x, Catch::Matchers::WithinAbs(libmesh_bounding_box.min_x, 1e-6)); + REQUIRE_THAT(moab_bounding_box.min_y, Catch::Matchers::WithinAbs(libmesh_bounding_box.min_y, 1e-6)); + REQUIRE_THAT(moab_bounding_box.min_z, Catch::Matchers::WithinAbs(libmesh_bounding_box.min_z, 1e-6)); + REQUIRE_THAT(moab_bounding_box.max_x, Catch::Matchers::WithinAbs(libmesh_bounding_box.max_x, 1e-6)); + REQUIRE_THAT(moab_bounding_box.max_y, Catch::Matchers::WithinAbs(libmesh_bounding_box.max_y, 1e-6)); + REQUIRE_THAT(moab_bounding_box.max_z, Catch::Matchers::WithinAbs(libmesh_bounding_box.max_z, 1e-6)); + + REQUIRE_THAT(moab_bounding_box.min_x, Catch::Matchers::WithinAbs(mfem_bounding_box.min_x, 1e-6)); + REQUIRE_THAT(moab_bounding_box.min_y, Catch::Matchers::WithinAbs(mfem_bounding_box.min_y, 1e-6)); + REQUIRE_THAT(moab_bounding_box.min_z, Catch::Matchers::WithinAbs(mfem_bounding_box.min_z, 1e-6)); + REQUIRE_THAT(moab_bounding_box.max_x, Catch::Matchers::WithinAbs(mfem_bounding_box.max_x, 1e-6)); + REQUIRE_THAT(moab_bounding_box.max_y, Catch::Matchers::WithinAbs(mfem_bounding_box.max_y, 1e-6)); + REQUIRE_THAT(moab_bounding_box.max_z, Catch::Matchers::WithinAbs(mfem_bounding_box.max_z, 1e-6)); + + // sample start and end locations within the bounding box of these models + int num_samples = 10000; + for (int i = 0; i < num_samples; i++) { + Position start = moab_bounding_box.sample_location(); + Position end = moab_bounding_box.sample_location(); + + auto moab_element = xdg_moab->find_element(start); + auto libmesh_element = xdg_libmesh->find_element(start); + auto mfem_element = xdg_mfem->find_element(start); + + if (libmesh_element == ID_NONE) { + // we want the others to be ID_NONE as well + REQUIRE(moab_element == ID_NONE); + REQUIRE(mfem_element == ID_NONE); + continue; + } + + REQUIRE(libmesh_element != ID_NONE); + REQUIRE(moab_element != ID_NONE); + REQUIRE(mfem_element != ID_NONE); + + // check element equivalence by index b/c IDs may be different depending on the library conventions + REQUIRE(xdg_moab->mesh_manager()->element_index(moab_element) == xdg_libmesh->mesh_manager()->element_index(libmesh_element)); + REQUIRE(xdg_moab->mesh_manager()->element_index(moab_element) == xdg_mfem->mesh_manager()->element_index(mfem_element)); + + auto moab_tracks = xdg_moab->segments(start, end); + auto libmesh_tracks = xdg_libmesh->segments(start, end); + auto mfem_tracks = xdg_mfem->segments(start, end); + + REQUIRE(moab_tracks.size() == libmesh_tracks.size()); + for (size_t j = 0; j < moab_tracks.size(); j++) { + REQUIRE(xdg_moab->mesh_manager()->element_index(moab_tracks[j].first) == xdg_libmesh->mesh_manager()->element_index(libmesh_tracks[j].first)); + REQUIRE(xdg_moab->mesh_manager()->element_index(moab_tracks[j].first) == xdg_mfem->mesh_manager()->element_index(mfem_tracks[j].first)); + + REQUIRE_THAT(moab_tracks[j].second, Catch::Matchers::WithinAbs(libmesh_tracks[j].second, 1e-10)); + REQUIRE_THAT(moab_tracks[j].second, Catch::Matchers::WithinAbs(mfem_tracks[j].second, 1e-10)); + } + } +} diff --git a/tests/test_mfem.cpp b/tests/test_mfem.cpp new file mode 100644 index 00000000..d0c58be2 --- /dev/null +++ b/tests/test_mfem.cpp @@ -0,0 +1,454 @@ +// stl includes +#include + + +// testing includes +#include +#include + +// xdg includes +#include "xdg/error.h" +#include "xdg/mesh_managers.h" +#include "xdg/xdg.h" +#include "util.h" + + +using namespace xdg; + +TEST_CASE("Test MFEM Initialization") +{ + std::unique_ptr mesh_manager = std::make_unique(); + + mesh_manager->load_file("cyl-brick.exo"); + mesh_manager->init(); + + REQUIRE(mesh_manager->num_volume_elements() == 16624); + + // property type +} + +// Read in the brick, read the element type and check it matches what +// we are expecting. +// The brick is meshed with tets +TEST_CASE("MFEM element types") +{ + std::unique_ptr mesh_manager = std::make_unique(); + mesh_manager->load_file("brick.exo"); + + mesh_manager->init(); + REQUIRE(mesh_manager->num_volume_elements() == 8790); + + // At time of writing, brick.exo does not have sidesets labelled, so we just check + // each of the elements + // for (int i=0; inum_boundary_elements(); i++) + // REQUIRE( mesh_manager->get_surface_element_type(i) == SurfaceElementType::TRI ); + +} + +// next, emulate the Find Element Method +TEST_CASE("TEST MFEM Find Element Method") +{ + std::shared_ptr xdg = XDG::create(MeshLibrary::MFEM); + + const auto& mesh_manager = xdg->mesh_manager(); + mesh_manager->load_file("jezebel.exo"); + mesh_manager->init(); + + size_t num_elements = mesh_manager->num_volume_elements(); + REQUIRE(num_elements == 10333); + + xdg->prepare_raytracer(); + + MeshID volume = 1; + MeshID element = xdg->find_element(volume, {0.0, 0.0, 100.0}); + REQUIRE(element == ID_NONE); // should not find an element since the point is outside the volume + + element = xdg->find_element(volume, {0.0, 0.0, 0.0}); + REQUIRE(element != ID_NONE); // should find an element + + auto next_element = xdg->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_CASE("TEST Ray Fire Brick") +{ + std::shared_ptr xdg = XDG::create(MeshLibrary::MFEM); + + const auto& mesh_manager = xdg->mesh_manager(); + mesh_manager->load_file("brick.exo"); + mesh_manager->init(); + xdg->prepare_raytracer(); + + MeshID volume = 1; + + Position origin {0.0, 0.0, 0.0}; + Direction direction {0.0, 0.0, 1.0}; + std::pair intersection; + + intersection = xdg->ray_fire(volume, origin, direction); + REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(5.0, 1e-6)); + + origin = {0.0, 0.0, 0.0}; + REQUIRE(xdg->point_in_volume(volume, origin)); +} + +TEST_CASE("Test Ray Fire Jezebel") +{ + std::shared_ptr xdg = XDG::create(MeshLibrary::MFEM); + xdg->mesh_manager()->mesh_library(); + REQUIRE(xdg->mesh_manager()->mesh_library() == MeshLibrary::MFEM); + const auto& mesh_manager = xdg->mesh_manager(); + mesh_manager->load_file("jezebel.exo"); + mesh_manager->init(); + xdg->prepare_raytracer(); + + MeshID volume = 1; + + // fire ray from the center of the cube + Position origin {0.0, 0.0, 0.0}; + Direction direction {0.0, 0.0, 1.0}; + + int n_rays {1000}; + + for (int i = 0; i < n_rays; i++) { + direction = rand_dir(); + std::pair intersection; + intersection = xdg->ray_fire(volume, origin, direction); + if (intersection.second == ID_NONE) + fatal_error("Ray did not intersect any geometry"); + if (intersection.first > 6.4) { + fatal_error("Ray intersected geometry at distance greater than 6.4 cm"); + } + REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(6.3849, 1e-1)); + } +} + +TEST_CASE("Test Cylinder-Brick Initialization") +{ + std::unique_ptr mesh_manager {std::make_unique()}; + + mesh_manager->load_file("cyl-brick.exo"); + + mesh_manager->init(); + + REQUIRE(mesh_manager->num_volumes() == 3); + + REQUIRE(mesh_manager->num_surfaces() == 12); + + // get an element from each volume and check its volume ID + auto vol1_elems = mesh_manager->get_volume_elements(1); + REQUIRE(!vol1_elems.empty()); + + auto vol2_elems = mesh_manager->get_volume_elements(2); + REQUIRE(!vol2_elems.empty()); + + mesh_manager->parse_metadata(); + + xdg::Property prop; + + // check the model's metadata + prop = mesh_manager->get_volume_property(1, PropertyType::MATERIAL); + REQUIRE(prop.type == PropertyType::MATERIAL); + REQUIRE(prop.value == "steel"); + + prop = mesh_manager->get_volume_property(2, PropertyType::MATERIAL); + REQUIRE(prop.type == PropertyType::MATERIAL); + REQUIRE(prop.value == "iron"); + + for (auto s : mesh_manager->surfaces()) { + prop = mesh_manager->get_surface_property(s, PropertyType::BOUNDARY_CONDITION); + std::cout << s << ", " << prop.value << std::endl; + REQUIRE(prop.type == PropertyType::BOUNDARY_CONDITION); + if (s == 3) { + REQUIRE(prop.value == "transmission"); + } else if (s == 4) { + REQUIRE(prop.value == "reflective"); + } else { + REQUIRE(prop.value == "vacuum"); + } + } +} + +TEST_CASE("Test Brick w/ Sidesets") +{ + std::unique_ptr mesh_manager {std::make_unique()}; + mesh_manager->load_file("brick-sidesets.exo"); + mesh_manager->init(); + + REQUIRE(mesh_manager->num_volumes() == 2); + REQUIRE(mesh_manager->num_surfaces() == 6); +} + +TEST_CASE("Test BVH Build Brick") +{ + std::shared_ptr mesh_manager = std::make_shared(); + + mesh_manager->load_file("brick.exo"); + mesh_manager->init(); + + REQUIRE(mesh_manager->num_volumes() == 2); + REQUIRE(mesh_manager->num_surfaces() == 1); + + std::unique_ptr ray_tracing_interface = std::make_unique(); + for (auto volume : mesh_manager->volumes()) { + ray_tracing_interface->register_volume(mesh_manager, volume); + } + + // volume elements will be detected on the mfem mesh, so three trees will be registered + REQUIRE(ray_tracing_interface->num_registered_trees() == 3); +} + + +TEST_CASE("Test BVH Build Brick w/ Sidesets") +{ + std::shared_ptr mesh_manager = std::make_shared(); + mesh_manager->load_file("brick-sidesets.exo"); + mesh_manager->init(); + + REQUIRE(mesh_manager->num_volumes() == 2); + REQUIRE(mesh_manager->num_surfaces() == 6); + + std::unique_ptr ray_tracing_interface = std::make_unique(); + + for (auto volume : mesh_manager->volumes()) { + ray_tracing_interface->register_volume(mesh_manager, volume); + } + // volume elements will be detected on the mfem mesh, so two trees will be registered + REQUIRE(ray_tracing_interface->num_registered_trees() == 3); +} + + +TEST_CASE("Test Ray Fire Brick") +{ + std::shared_ptr xdg = XDG::create(MeshLibrary::MFEM); + xdg->mesh_manager()->mesh_library(); + REQUIRE(xdg->mesh_manager()->mesh_library() == MeshLibrary::MFEM); + const auto& mesh_manager = xdg->mesh_manager(); + mesh_manager->load_file("brick.exo"); + mesh_manager->init(); + xdg->prepare_raytracer(); + + MeshID volume = 1; + + Position origin {0.0, 0.0, 0.0}; + Direction direction {0.0, 0.0, 1.0}; + std::pair intersection; + + intersection = xdg->ray_fire(volume, origin, direction); + // this cube is 10 cm on a side, so the ray should hit the surface at 5 cm + REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(5.0, 1e-6)); + + origin = {0.0, 0.0, 0.0}; + REQUIRE(xdg->point_in_volume(volume, origin)); +} + +TEST_CASE("Test Ray Fire Cylinder-Brick") +{ + std::shared_ptr xdg = XDG::create(MeshLibrary::MFEM); + xdg->mesh_manager()->mesh_library(); + REQUIRE(xdg->mesh_manager()->mesh_library() == MeshLibrary::MFEM); + const auto& mesh_manager = xdg->mesh_manager(); + mesh_manager->load_file("cyl-brick.exo"); + mesh_manager->init(); + xdg->prepare_raytracer(); + + MeshID volume = 2; + + // fire ray from the center of the cube + Position origin {0.0, 0.0, 10.0}; + Direction direction {0.0, 0.0, 1.0}; + std::pair intersection; + intersection = xdg->ray_fire(volume, origin, direction); + // this cube is 10 cm on a side, so the ray should hit the surface at 5 cm + REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(5.0, 1e-6)); + + // fire ray in the opposite direction + direction = {0.0, 0.0, -1.0}; + intersection = xdg->ray_fire(volume, origin, direction); + // this cube is 10 cm on a side, so the ray should hit the surface at 5 cm + REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(5.0, 1e-6)); + + // fire ray from the center of the cylinder in the negative z direction + volume = 1; + origin = {0.0, 0.0, 0.0}; + intersection = xdg->ray_fire(volume, origin, direction); + REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(5.0, 1e-6)); + + // fire a ray from the center of the cylinder in the positive z direction + direction = {0.0, 0.0, 1.0}; + intersection = xdg->ray_fire(volume, origin, direction); + REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(5.0, 1e-6)); + + // fire a ray from the center of the cylinder in the positive x direction + direction = {1.0, 0.0, 0.0}; + intersection = xdg->ray_fire(volume, origin, direction); + REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(5.0, 1e-3)); + + volume = 1; + origin = {0.0, 0.0, 0.0}; + REQUIRE(xdg->point_in_volume(volume, origin)); + + volume = 2; + origin = {0.0, 0.0, 10.0}; + REQUIRE(xdg->point_in_volume(volume, origin)); +} + +TEST_CASE("Test Volume Element Count Jezebel") +{ + std::shared_ptr xdg = XDG::create(MeshLibrary::MFEM); + xdg->mesh_manager()->mesh_library(); + REQUIRE(xdg->mesh_manager()->mesh_library() == MeshLibrary::MFEM); + const auto& mesh_manager = xdg->mesh_manager(); + mesh_manager->load_file("jezebel.exo"); + mesh_manager->init(); + xdg->prepare_raytracer(); + + MeshID volume = 1; + + auto elements = mesh_manager->get_volume_elements(volume); + REQUIRE(elements.size() == 10333); + REQUIRE(mesh_manager->num_volume_elements() == 10333); +} + +TEST_CASE("Test Point Location Jezebel") +{ + std::shared_ptr xdg = XDG::create(MeshLibrary::MFEM); + xdg->mesh_manager()->mesh_library(); + REQUIRE(xdg->mesh_manager()->mesh_library() == MeshLibrary::MFEM); + const auto& mesh_manager = xdg->mesh_manager(); + mesh_manager->load_file("jezebel.exo"); + mesh_manager->init(); + xdg->prepare_raytracer(); + + MeshID volume = 1; + + // fire ray from the center of the cube + Position origin {0.0, 0.0, 0.0}; + Direction direction {0.0, 0.0, 1.0}; + + // the origin of the problem should be in the volume + MeshID volume_id = xdg->find_volume(origin, direction); + REQUIRE(volume_id == volume); + + // a point outside of the sphere should be in the implicit complement + origin = {0.0, 0.0, 10.0}; + volume_id = xdg->find_volume(origin, direction); + REQUIRE(volume_id == xdg->mesh_manager()->implicit_complement()); +} + +TEST_CASE("Test Point Location Cylinder-Brick") +{ + std::shared_ptr xdg = XDG::create(MeshLibrary::MFEM); + xdg->mesh_manager()->mesh_library(); + REQUIRE(xdg->mesh_manager()->mesh_library() == MeshLibrary::MFEM); + const auto& mesh_manager = xdg->mesh_manager(); + mesh_manager->load_file("cyl-brick.exo"); + mesh_manager->init(); + xdg->prepare_raytracer(); + + REQUIRE(mesh_manager->num_volume_elements() == 16624); + + MeshID expected_volume = 1; + + // fire ray from the center of the cube + Position origin {0.0, 0.0, 0.0}; + Direction direction {0.0, 0.0, 1.0}; + + // test a point inside the cylinder + MeshID volume_id = xdg->find_volume(origin, direction); + REQUIRE(volume_id == expected_volume); + + // test a point inside the cube + expected_volume = 2; + origin = {0.0, 0.0, 10.0}; + volume_id = xdg->find_volume(origin, direction); + REQUIRE(volume_id == expected_volume); + + // a point outside of the sphere should be in the implicit complement + origin = {0.0, 0.0, 100.0}; + volume_id = xdg->find_volume(origin, direction); + REQUIRE(volume_id == xdg->mesh_manager()->implicit_complement()); +} + +TEST_CASE("Test Volume Element Count Cylinder-Brick") +{ + std::shared_ptr xdg = XDG::create(MeshLibrary::MFEM); + xdg->mesh_manager()->mesh_library(); + REQUIRE(xdg->mesh_manager()->mesh_library() == MeshLibrary::MFEM); + const auto& mesh_manager = xdg->mesh_manager(); + mesh_manager->load_file("cyl-brick.exo"); + mesh_manager->init(); + xdg->prepare_raytracer(); + + MeshID volume = 1; + auto elements = mesh_manager->get_volume_elements(volume); + REQUIRE(elements.size() == 7587); + + volume = 2; + elements = mesh_manager->get_volume_elements(volume); + REQUIRE(elements.size() == 9037); +} + +TEST_CASE("Test Find Element Brick") +{ + std::shared_ptr xdg = XDG::create(MeshLibrary::MFEM); + xdg->mesh_manager()->mesh_library(); + REQUIRE(xdg->mesh_manager()->mesh_library() == MeshLibrary::MFEM); + const auto& mesh_manager = xdg->mesh_manager(); + mesh_manager->load_file("brick.exo"); + mesh_manager->init(); + xdg->prepare_raytracer(); + + REQUIRE(mesh_manager->num_volume_elements() == 8790); + + MeshID volume = 1; + + MeshID element = xdg->find_element(volume, {0.0, 0.0, 0.0}); + REQUIRE(element != ID_NONE); + + element = xdg->find_element(volume, {0.0, 0.0, 100.0}); + REQUIRE(element == ID_NONE); +} + +TEST_CASE("Test Track Exiting Mesh Brick") +{ + std::shared_ptr xdg = XDG::create(MeshLibrary::MFEM); + xdg->mesh_manager()->mesh_library(); + REQUIRE(xdg->mesh_manager()->mesh_library() == MeshLibrary::MFEM); + const auto& mesh_manager = xdg->mesh_manager(); + mesh_manager->load_file("brick.exo"); + mesh_manager->init(); + xdg->prepare_raytracer(); + + MeshID volume = 1; + Position start {0.0, 0.0, -1000.0}; + Position end {0.0, 0.0, 1000.0}; + auto tracks = xdg->segments(volume, start, end); + + 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)); +} + +// TEST_CASE("Multiblock sidesets") +// { +// std::shared_ptr xdg = XDG::create(MeshLibrary::MFEM); +// REQUIRE(xdg->mesh_manager()->mesh_library() == MeshLibrary::MFEM); +// const auto& mesh_manager = xdg->mesh_manager(); +// mesh_manager->load_file("cube-w-multiblock-sideset.exo"); +// mesh_manager->init(); +// mesh_manager->parse_metadata(); + +// xdg->prepare_raytracer(); + +// MeshID volume = 1; +// Position start {0.0, 0.0, 0.0}; +// Position end {15.0, 20.0, 3.0}; +// auto tracks = xdg->segments(volume, start, end); + +// assert(tracks.size() > 0); +// } + diff --git a/tests/util.h b/tests/util.h index a9841086..9b4c538b 100644 --- a/tests/util.h +++ b/tests/util.h @@ -72,6 +72,13 @@ inline bool mesh_library_available(xdg::MeshLibrary mesh) { #else return false; #endif + + case xdg::MeshLibrary::MFEM: + #ifdef XDG_ENABLE_MFEM + return true; + #else + return false; + #endif } return false;