From 9708a6a9c38a352d66b8aa181399986b000d3243 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sat, 22 Nov 2025 22:51:03 -0600 Subject: [PATCH 01/19] Minimal linkage to MFEM --- CMakeLists.txt | 26 +++++++++++ include/xdg/constants.h | 6 ++- include/xdg/mesh_managers.h | 4 ++ include/xdg/mfem/mesh_manager.h | 76 +++++++++++++++++++++++++++++++++ src/mfem/mesh_manager.cpp | 20 +++++++++ 5 files changed, 130 insertions(+), 2 deletions(-) create mode 100644 include/xdg/mfem/mesh_manager.h create mode 100644 src/mfem/mesh_manager.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 0c185556..146e7e20 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 @@ -240,6 +251,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 +332,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 +391,10 @@ if (XDG_ENABLE_MOAB) target_link_libraries(xdg PRIVATE MOAB) endif() +if (XDG_ENABLE_MFEM) + target_link_libraries(xdg mfem) +endif() + #================================================================= # Installation & Packaging #================================================================= 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/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..7f43278a --- /dev/null +++ b/include/xdg/mfem/mesh_manager.h @@ -0,0 +1,76 @@ +#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; + } + } + + 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_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"); + } + + BoundingBox element_bounding_box(MeshID element) const override { + mesh_.Get + } + + // Data members +private: + std::unique_ptr mesh_ {nullptr}; +}; + +} + +#endif // include guard \ No newline at end of file diff --git a/src/mfem/mesh_manager.cpp b/src/mfem/mesh_manager.cpp new file mode 100644 index 00000000..809500be --- /dev/null +++ b/src/mfem/mesh_manager.cpp @@ -0,0 +1,20 @@ +#include + +#include "xdg/mfem/mesh_manager.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"); + } + + // Finalize the mesh setup + mesh_->FinalizeTopology(); +} + +} // namespace xdg \ No newline at end of file From 7e46c85dfca60aa570bb6300c0c0a39ca3e63449 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 25 Nov 2025 12:21:25 -0600 Subject: [PATCH 02/19] Adding element bounding box computation based on element vertices --- include/xdg/bbox.h | 9 +++++++++ include/xdg/mfem/mesh_manager.h | 14 +++++++++++++- 2 files changed, 22 insertions(+), 1 deletion(-) 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/mfem/mesh_manager.h b/include/xdg/mfem/mesh_manager.h index 7f43278a..9130feed 100644 --- a/include/xdg/mfem/mesh_manager.h +++ b/include/xdg/mfem/mesh_manager.h @@ -63,7 +63,19 @@ class MfemMeshManager : public MeshManager { } BoundingBox element_bounding_box(MeshID element) const override { - mesh_.Get + auto elem = mesh_->GetElement(element); + if (!elem) { + fatal_error(fmt::format("MfemMeshManager::element_bounding_box(): invalid element ID {}", element)); + } + auto bbox = BoundingBox { + INFTY, INFTY, INFTY, + -INFTY, -INFTY, -INFTY + }; + for(int i = 0; i < elem->GetNVertices(); i++) { + const mfem::real_t* v = mesh_->GetVertex(*(i + elem->GetVertices())); + bbox.update(v); + } + return bbox; } // Data members From b7b258ec67f460146dc09001881964bd49e24aa9 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 25 Nov 2025 15:59:41 -0600 Subject: [PATCH 03/19] Add MFEM manger instantiation and file load test --- include/xdg/mfem/mesh_manager.h | 80 ++++++++++++++++++++++++++++++--- tests/CMakeLists.txt | 5 +++ tests/test_mfem.cpp | 23 ++++++++++ 3 files changed, 102 insertions(+), 6 deletions(-) create mode 100644 tests/test_mfem.cpp diff --git a/include/xdg/mfem/mesh_manager.h b/include/xdg/mfem/mesh_manager.h index 9130feed..a0501df5 100644 --- a/include/xdg/mfem/mesh_manager.h +++ b/include/xdg/mfem/mesh_manager.h @@ -13,7 +13,8 @@ namespace xdg { class MfemMeshManager : public MeshManager { public: - MfemMeshManager(); + + MfemMeshManager() {}; ~MfemMeshManager() override = default; @@ -62,11 +63,43 @@ class MfemMeshManager : public MeshManager { fatal_error("MfemMeshManager::num_surface_faces() not implemented yet"); } - BoundingBox element_bounding_box(MeshID element) const override { - auto elem = mesh_->GetElement(element); - if (!elem) { - fatal_error(fmt::format("MfemMeshManager::element_bounding_box(): invalid element ID {}", element)); - } + virtual std::vector get_volume_elements(MeshID volume) const override { + fatal_error("MfemMeshManager::get_volume_elements() not implemented yet"); + } + + virtual std::vector get_surface_faces(MeshID surface) const override { + fatal_error("MfemMeshManager::get_surface_faces() not implemented yet"); + } + + virtual std::vector element_vertices(MeshID element) const override { + fatal_error("MfemMeshManager::element_vertices() not implemented yet"); + } + + virtual std::array face_vertices(MeshID element) const override { + fatal_error("MfemMeshManager::face_vertices() not implemented yet"); + } + + virtual std::vector get_surface_vertices(MeshID surface) const override { + fatal_error("MfemMeshManager::get_surface_vertices() not implemented yet"); + } + + virtual MeshID adjacent_element(MeshID element, int face) const override { + fatal_error("MfemMeshManager::adjacent_element() not implemented yet"); + } + + virtual Sense surface_sense(MeshID surface, MeshID volume) const override { + fatal_error("MfemMeshManager::surface_sense() not implemented yet"); + } + + virtual std::pair, std::vector> get_surface_mesh(MeshID surface) const override { + fatal_error("MfemMeshManager::get_surface_mesh() not implemented yet"); + } + + virtual SurfaceElementType get_surface_element_type(MeshID element) const override { + fatal_error("MfemMeshManager::get_surface_element_type() not implemented yet"); + } + + BoundingBox _mfem_element_bounding_box(mfem::Element* elem) const { auto bbox = BoundingBox { INFTY, INFTY, INFTY, -INFTY, -INFTY, -INFTY @@ -78,6 +111,41 @@ class MfemMeshManager : public MeshManager { return bbox; } + BoundingBox element_bounding_box(MeshID element) const override { + auto elem = mesh_->GetElement(element); + if (!elem) { + fatal_error(fmt::format("MfemMeshManager::element_bounding_box(): invalid element ID {}", element)); + } + return _mfem_element_bounding_box(elem); + } + + // Topology + std::vector get_volume_surfaces(MeshID volume) const override { + fatal_error("MfemMeshManager::get_volume_surfaces() not implemented yet"); + } + + std::pair surface_senses(MeshID surface) const override { + fatal_error("MfemMeshManager::surface_senses() not implemented yet"); + } + + MeshID create_volume() override { + fatal_error("MfemMeshManager::create_volume() not implemented yet"); + } + + void add_surface_to_volume(MeshID volume, MeshID surface, Sense sense, bool overwrite=false) override { + fatal_error("MfemMeshManager::add_surface_to_volume() not implemented yet"); + } + + // Metadata methods + void parse_metadata() override { + fatal_error("MfemMeshManager::parse_metadata() not implemented yet"); + } + + // Accessors + const std::unique_ptr& mfem_mesh() const { + return mesh_; + } + // Data members private: std::unique_ptr mesh_ {nullptr}; diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 23b95f1d..4af15d5d 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -31,6 +31,10 @@ if (XDG_ENABLE_LIBMESH) list(APPEND TEST_NAMES test_libmesh) endif() +if (XDG_ENABLE_MFEM) + list(APPEND TEST_NAMES test_mfem) +endif() + if (XDG_ENABLE_MOAB AND XDG_ENABLE_LIBMESH) 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_mfem.cpp b/tests/test_mfem.cpp new file mode 100644 index 00000000..c65af493 --- /dev/null +++ b/tests/test_mfem.cpp @@ -0,0 +1,23 @@ +// stl includes +#include + + +// testing includes +#include + +// xdg includes +#include "xdg/error.h" +#include "xdg/mesh_managers.h" + + +using namespace xdg; + +TEST_CASE("Test MFEM Initialization") +{ + std::unique_ptr mesh_manager = std::make_unique(); + + mesh_manager->load_file("jezebel.exo"); + mesh_manager->init(); + + REQUIRE(mesh_manager->num_volume_elements() == 10333); +} \ No newline at end of file From a1aaec94a264b5f82dc6e9a6bb0727f3fab65c80 Mon Sep 17 00:00:00 2001 From: Sean Baccas Date: Thu, 16 Apr 2026 15:29:01 +0100 Subject: [PATCH 04/19] Rebase on main and make minimal changes to get it compiling again --- include/xdg/mfem/mesh_manager.h | 60 ++++++++++++++++++------------ include/xdg/moab/mesh_manager.h | 2 +- include/xdg/moab/tag_conventions.h | 2 +- 3 files changed, 38 insertions(+), 26 deletions(-) diff --git a/include/xdg/mfem/mesh_manager.h b/include/xdg/mfem/mesh_manager.h index a0501df5..81a89830 100644 --- a/include/xdg/mfem/mesh_manager.h +++ b/include/xdg/mfem/mesh_manager.h @@ -79,10 +79,6 @@ class MfemMeshManager : public MeshManager { fatal_error("MfemMeshManager::face_vertices() not implemented yet"); } - virtual std::vector get_surface_vertices(MeshID surface) const override { - fatal_error("MfemMeshManager::get_surface_vertices() not implemented yet"); - } - virtual MeshID adjacent_element(MeshID element, int face) const override { fatal_error("MfemMeshManager::adjacent_element() not implemented yet"); } @@ -91,34 +87,50 @@ class MfemMeshManager : public MeshManager { fatal_error("MfemMeshManager::surface_sense() not implemented yet"); } - virtual std::pair, std::vector> get_surface_mesh(MeshID surface) const override { - fatal_error("MfemMeshManager::get_surface_mesh() not implemented yet"); - } - virtual SurfaceElementType get_surface_element_type(MeshID element) const override { fatal_error("MfemMeshManager::get_surface_element_type() not implemented yet"); } - BoundingBox _mfem_element_bounding_box(mfem::Element* elem) const { - auto bbox = BoundingBox { - INFTY, INFTY, INFTY, - -INFTY, -INFTY, -INFTY - }; - for(int i = 0; i < elem->GetNVertices(); i++) { - const mfem::real_t* v = mesh_->GetVertex(*(i + elem->GetVertices())); - bbox.update(v); - } - return bbox; + virtual int num_vertices() const override { + fatal_error("MfemMeshManager::num_vertices() not implemented yet"); } - BoundingBox element_bounding_box(MeshID element) const override { - auto elem = mesh_->GetElement(element); - if (!elem) { - fatal_error(fmt::format("MfemMeshManager::element_bounding_box(): invalid element ID {}", element)); - } - return _mfem_element_bounding_box(elem); + virtual double element_volume(MeshID element) const override { + fatal_error("MfemMeshManager::element_volume() not implemented yet"); + } + + virtual Vertex vertex_coordinates(MeshID vertex_id) const override { + fatal_error("MfemMeshManager::vertex_coordinates() not implemented yet"); + } + + 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"); + } + + // BoundingBox _mfem_element_bounding_box(mfem::Element* elem) const { + // auto bbox = BoundingBox { + // INFTY, INFTY, INFTY, + // -INFTY, -INFTY, -INFTY + // }; + // for(int i = 0; i < elem->GetNVertices(); i++) { + // const mfem::real_t* v = mesh_->GetVertex(*(i + elem->GetVertices())); + // bbox.update(v); + // } + // return bbox; + // } + + // BoundingBox element_bounding_box(MeshID element) const { + // auto elem = mesh_->GetElement(element); + // if (!elem) { + // fatal_error(fmt::format("MfemMeshManager::element_bounding_box(): invalid element ID {}", element)); + // } + // return _mfem_element_bounding_box(elem); + // } + // Topology std::vector get_volume_surfaces(MeshID volume) const override { fatal_error("MfemMeshManager::get_volume_surfaces() not implemented yet"); 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 { From 2b37d6cf34d8216a632ed68744e93040d32a19d1 Mon Sep 17 00:00:00 2001 From: Sean Baccas Date: Fri, 15 May 2026 12:42:40 +0100 Subject: [PATCH 05/19] Work in progress. Haven't tried next element near the boundary yet --- include/xdg/mfem/mesh_manager.h | 173 ++++++++++++++++++++--------- src/element_face_accessor.cpp | 9 ++ src/mfem/mesh_manager.cpp | 191 +++++++++++++++++++++++++++++++- src/xdg.cpp | 6 + tests/test_mfem.cpp | 52 ++++++++- tools/CMakeLists.txt | 1 + tools/mfem_tool.cpp | 36 ++++++ 7 files changed, 413 insertions(+), 55 deletions(-) create mode 100644 tools/mfem_tool.cpp diff --git a/include/xdg/mfem/mesh_manager.h b/include/xdg/mfem/mesh_manager.h index 81a89830..30a1dbc8 100644 --- a/include/xdg/mfem/mesh_manager.h +++ b/include/xdg/mfem/mesh_manager.h @@ -31,12 +31,13 @@ class MfemMeshManager : public MeshManager { // Interface methods MeshLibrary mesh_library() const override { return MeshLibrary::MFEM; } + // This info might not be available in mfem int num_volumes() const override { - return volumes_.size(); + return mesh_->attribute_sets.GetAttributeSetNames().size(); } int num_surfaces() const override { - return surfaces_.size(); + return mesh_->bdr_attribute_sets.GetAttributeSetNames().size(); } int num_ents_of_dimension(int dim) const override { @@ -47,6 +48,7 @@ class MfemMeshManager : public MeshManager { } } + // 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"); } @@ -55,6 +57,11 @@ class MfemMeshManager : public MeshManager { 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"); } @@ -63,44 +70,41 @@ class MfemMeshManager : public MeshManager { fatal_error("MfemMeshManager::num_surface_faces() not implemented yet"); } - virtual std::vector get_volume_elements(MeshID volume) const override { - fatal_error("MfemMeshManager::get_volume_elements() 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 { - fatal_error("MfemMeshManager::get_surface_faces() not implemented yet"); - } + virtual std::vector get_surface_faces(MeshID surface) const override; - virtual std::vector element_vertices(MeshID element) const override { - fatal_error("MfemMeshManager::element_vertices() not implemented yet"); - } + // see Mesh::GetElementVertices + virtual std::vector element_vertices(MeshID element) const override; + std::vector bdr_element_vertices(MeshID element) const; - virtual std::array face_vertices(MeshID element) const override { - fatal_error("MfemMeshManager::face_vertices() not implemented yet"); - } + // this one is very easy - Mesh::GetFaceVertices returns the coords of face i at the elment level + virtual std::array face_vertices(MeshID element) const override; - virtual MeshID adjacent_element(MeshID element, int face) const override { - fatal_error("MfemMeshManager::adjacent_element() not implemented yet"); - } + // The table works wonders for this + virtual MeshID adjacent_element(MeshID element, int face) const override; virtual Sense surface_sense(MeshID surface, MeshID volume) const override { fatal_error("MfemMeshManager::surface_sense() not implemented yet"); } - virtual SurfaceElementType get_surface_element_type(MeshID element) const override { - fatal_error("MfemMeshManager::get_surface_element_type() not implemented yet"); - } + // mesh_->GetElement(0)->GetGeometryType() + virtual SurfaceElementType get_surface_element_type(MeshID element) const override; virtual int num_vertices() const override { - fatal_error("MfemMeshManager::num_vertices() not implemented yet"); + return mesh_->GetNV(); } virtual double element_volume(MeshID element) const override { - fatal_error("MfemMeshManager::element_volume() not implemented yet"); + return mesh_->GetElementVolume(element); } virtual Vertex vertex_coordinates(MeshID vertex_id) const override { - fatal_error("MfemMeshManager::vertex_coordinates() not implemented yet"); + 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 { @@ -111,35 +115,13 @@ class MfemMeshManager : public MeshManager { fatal_error("MfemMeshManager::face_connectivity() not implemented yet"); } - // BoundingBox _mfem_element_bounding_box(mfem::Element* elem) const { - // auto bbox = BoundingBox { - // INFTY, INFTY, INFTY, - // -INFTY, -INFTY, -INFTY - // }; - // for(int i = 0; i < elem->GetNVertices(); i++) { - // const mfem::real_t* v = mesh_->GetVertex(*(i + elem->GetVertices())); - // bbox.update(v); - // } - // return bbox; - // } - - // BoundingBox element_bounding_box(MeshID element) const { - // auto elem = mesh_->GetElement(element); - // if (!elem) { - // fatal_error(fmt::format("MfemMeshManager::element_bounding_box(): invalid element ID {}", element)); - // } - // return _mfem_element_bounding_box(elem); - // } - // Topology - std::vector get_volume_surfaces(MeshID volume) const override { - fatal_error("MfemMeshManager::get_volume_surfaces() not implemented yet"); - } - std::pair surface_senses(MeshID surface) const override { - fatal_error("MfemMeshManager::surface_senses() not implemented yet"); - } + 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 { fatal_error("MfemMeshManager::create_volume() not implemented yet"); } @@ -161,8 +143,99 @@ class MfemMeshManager : public MeshManager { // 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 + 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 + std::set attributes_; +}; + +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); + + mfem::Array vertex_indices; + 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::array face_vertices(int i) const override { + std::array verts; + + // we have already gathered the vertices for this face. + // copy them into the output array + std::copy(face_vertices_[i], face_vertices_[i+1], verts.begin()); + + // 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( verts[0], verts[1] ); + + return verts; + } + + // 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_; + 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/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/mfem/mesh_manager.cpp b/src/mfem/mesh_manager.cpp index 809500be..24803f4c 100644 --- a/src/mfem/mesh_manager.cpp +++ b/src/mfem/mesh_manager.cpp @@ -4,7 +4,7 @@ namespace xdg { void MfemMeshManager::load_file(const std::string &filepath) { -mesh_ = std::make_unique(filepath.c_str(), 1, 1); + mesh_ = std::make_unique(filepath.c_str(), 1, 1); } void MfemMeshManager::init() { @@ -13,8 +13,193 @@ void MfemMeshManager::init() { fatal_error("Mesh must be 3-dimensional"); } - // Finalize the mesh setup - mesh_->FinalizeTopology(); + // 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 + for (int i=0; iGetNE(); i++) { + int volume_id = mesh_->GetAttribute(i); + volume_to_element_map_[volume_id].insert(i); + } + + // 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'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_)); +} + +// TODO: very slow, and could be done during init() +std::vector MfemMeshManager::get_volume_elements(MeshID volume) const { + if (attributes_.find(volume) == attributes_.end()) { + std::ostringstream output; + output << "Couldn't find volume " << volume << "\n"; + fatal_error(output.str()); + } + + std::vector 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; +} + +SurfaceElementType MfemMeshManager::get_surface_element_type(MeshID element) const { + auto mfem_element_type = mesh_->GetBdrElement(element)->GetType(); + return GetSurfaceElementTypeFromMfem(mfem_element_type); +} + +// Should return all of the sidesets that are a part of this volume +std::vector MfemMeshManager::get_volume_surfaces(MeshID volume) const { + // get the set associated with this volume + const std::set& sidesets = volumes_to_sidesets_.at(volume); + + // create a vector from this set + std::vector output(sidesets.begin(), sidesets.end()); + + return output; +} + +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()); + + return output; +} + +std::array MfemMeshManager::face_vertices(MeshID element) const { + std::array output; + + // create an mfem array to be passed into Mesh::GetFaceVertices. + // this gets populated with the indices of the vertices itself + mfem::Array index_array; + mesh_->GetFaceVertices(element, index_array); + + for (int i=0; iGetVertex( index_array[i] ); + + for (int d=0; dSpaceDimension(); d++) output[i][d] = vertices[d]; + } + + return output; +} + +std::pair MfemMeshManager::surface_senses(MeshID surface) const { + // I am trying to get the raytracer preparation routines working with + // the jezebel, so just return {-1, 1}. i.e. implicit_complement, interior_volume. + // Even though we haven't created implicit_complement yet. + warning("MfemMeshManager::surface_senses() is hardcoded for jezebel"); + return {-1,1}; +} + +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); + + // face is in range [0,3). So we just need the one + // that the caller asked for + return faces[face]; +} + +// 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"); + } +} + +// this second function is somewhat redundant. The mfem enum captures all +// of the possible geometries, in all possible dimensions... +SurfaceElementType GetSurfaceElementTypeFromMfem( mfem::Element::Type t ) { + switch (t) { + case mfem::Element::TRIANGLE: return SurfaceElementType::TRI; + case mfem::Element::QUADRILATERAL: return SurfaceElementType::QUAD; + default: + fatal_error("Unsupported element type\n"); + } +} + + } // namespace xdg \ No newline at end of file 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/test_mfem.cpp b/tests/test_mfem.cpp index c65af493..10fcd5df 100644 --- a/tests/test_mfem.cpp +++ b/tests/test_mfem.cpp @@ -8,6 +8,8 @@ // xdg includes #include "xdg/error.h" #include "xdg/mesh_managers.h" +#include "xdg/xdg.h" +#include "util.h" using namespace xdg; @@ -16,8 +18,54 @@ 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 MOAB Find Element Method") +{ + std::shared_ptr xdg = XDG::create(MeshLibrary::MFEM); // let it pick whichever raytracer + + const auto& mesh_manager = xdg->mesh_manager(); mesh_manager->load_file("jezebel.exo"); mesh_manager->init(); - REQUIRE(mesh_manager->num_volume_elements() == 10333); -} \ No newline at end of file + 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); +} diff --git a/tools/CMakeLists.txt b/tools/CMakeLists.txt index 50a0c650..79949eec 100644 --- a/tools/CMakeLists.txt +++ b/tools/CMakeLists.txt @@ -7,6 +7,7 @@ point_in_volume overlap_check walk_elements tally_segments +mfem_tool ) #=============================================================================== diff --git a/tools/mfem_tool.cpp b/tools/mfem_tool.cpp new file mode 100644 index 00000000..512cf935 --- /dev/null +++ b/tools/mfem_tool.cpp @@ -0,0 +1,36 @@ +#include +#include +#include + +#include "xdg/error.h" +#include "xdg/mesh_manager_interface.h" +#include "xdg/mesh_managers.h" +#include "xdg/vec3da.h" +#include "xdg/xdg.h" + +#include "argparse/argparse.hpp" + +#include "particle_sim.h" + +using namespace xdg; + +int main(int argc, char** argv) { + + std::unique_ptr mesh_manager = std::make_unique(); + argparse::ArgumentParser args("MFEM debugging tool", "1.0", argparse::default_arguments::help); + + args.add_argument("filename").help("Path to the input file"); + + try { + args.parse_args(argc, argv); + } + catch (const std::runtime_error& err) { + std::cout << err.what() << std::endl; + std::cout << args; + exit(0); + } + + mesh_manager->load_file(args.get("filename")); + mesh_manager->init(); + +} From cceeea7d3a8d04e34f84e39ddff670d81567ddf6 Mon Sep 17 00:00:00 2001 From: Sean Baccas Date: Wed, 20 May 2026 12:20:36 +0100 Subject: [PATCH 06/19] Minor additions to get ray_fire working --- include/xdg/mfem/mesh_manager.h | 3 +++ src/embree/ray_tracer.cpp | 6 +++--- src/mfem/mesh_manager.cpp | 37 +++++++++++++++++++++++++++------ tests/test_mfem.cpp | 25 +++++++++++++++++++++- 4 files changed, 61 insertions(+), 10 deletions(-) diff --git a/include/xdg/mfem/mesh_manager.h b/include/xdg/mfem/mesh_manager.h index 30a1dbc8..b1bba370 100644 --- a/include/xdg/mfem/mesh_manager.h +++ b/include/xdg/mfem/mesh_manager.h @@ -157,6 +157,9 @@ class MfemMeshManager : public MeshManager { // set to capture all of the valid volumes/attributes // It's a set (not vector) to prevent double counting std::set attributes_; + + int num_interior_faces_; + int num_boundary_faces_; }; struct MfemMeshElementFaceAccessor : public ElementFaceAccessor { 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 index 24803f4c..564f7b7e 100644 --- a/src/mfem/mesh_manager.cpp +++ b/src/mfem/mesh_manager.cpp @@ -56,6 +56,10 @@ void MfemMeshManager::init() { // 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_)); + + // set these two attributes related to interior/boundary faces + num_interior_faces_ = mesh_->GetNumFaces(); + num_boundary_faces_ = mesh_->GetNBE(); } // TODO: very slow, and could be done during init() @@ -101,16 +105,37 @@ std::vector MfemMeshManager::get_surface_faces(MeshID surface) const { // 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::array MfemMeshManager::face_vertices(MeshID element) const { std::array output; - - // create an mfem array to be passed into Mesh::GetFaceVertices. - // this gets populated with the indices of the vertices itself mfem::Array index_array; - mesh_->GetFaceVertices(element, 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 { + // create an mfem array to be passed into Mesh::GetFaceVertices. + // this gets populated with the indices of the vertices itself + mesh_->GetFaceVertices(element, index_array); + } for (int i=0; iGetVertex( index_array[i] ); @@ -125,8 +150,8 @@ std::pair MfemMeshManager::surface_senses(MeshID surface) const { // I am trying to get the raytracer preparation routines working with // the jezebel, so just return {-1, 1}. i.e. implicit_complement, interior_volume. // Even though we haven't created implicit_complement yet. - warning("MfemMeshManager::surface_senses() is hardcoded for jezebel"); - return {-1,1}; + warning("MfemMeshManager::surface_senses() is hardcoded for single-volume meshes"); + return {1,2}; } std::vector MfemMeshManager::element_vertices(MeshID element) const { diff --git a/tests/test_mfem.cpp b/tests/test_mfem.cpp index 10fcd5df..2e4498cd 100644 --- a/tests/test_mfem.cpp +++ b/tests/test_mfem.cpp @@ -4,6 +4,7 @@ // testing includes #include +#include // xdg includes #include "xdg/error.h" @@ -47,7 +48,7 @@ TEST_CASE("MFEM element types") // next, emulate the Find Element Method TEST_CASE("TEST MOAB Find Element Method") { - std::shared_ptr xdg = XDG::create(MeshLibrary::MFEM); // let it pick whichever raytracer + std::shared_ptr xdg = XDG::create(MeshLibrary::MFEM); const auto& mesh_manager = xdg->mesh_manager(); mesh_manager->load_file("jezebel.exo"); @@ -69,3 +70,25 @@ TEST_CASE("TEST MOAB Find Element Method") 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)); +} From 181f8b97967e24952fb2cbf214521aa3d4b51a10 Mon Sep 17 00:00:00 2001 From: Sean Baccas Date: Fri, 26 Jun 2026 11:33:42 +0100 Subject: [PATCH 07/19] Adding these minor changes before rebase --- src/mfem/mesh_manager.cpp | 4 ++++ tests/test_mfem.cpp | 2 +- 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/src/mfem/mesh_manager.cpp b/src/mfem/mesh_manager.cpp index 564f7b7e..92cfc0e2 100644 --- a/src/mfem/mesh_manager.cpp +++ b/src/mfem/mesh_manager.cpp @@ -151,6 +151,10 @@ std::pair MfemMeshManager::surface_senses(MeshID surface) const { // the jezebel, so just return {-1, 1}. i.e. implicit_complement, interior_volume. // Even though we haven't created implicit_complement yet. warning("MfemMeshManager::surface_senses() is hardcoded for single-volume meshes"); + + // 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 {1,2}; } diff --git a/tests/test_mfem.cpp b/tests/test_mfem.cpp index 2e4498cd..32561999 100644 --- a/tests/test_mfem.cpp +++ b/tests/test_mfem.cpp @@ -46,7 +46,7 @@ TEST_CASE("MFEM element types") } // next, emulate the Find Element Method -TEST_CASE("TEST MOAB Find Element Method") +TEST_CASE("TEST MFEM Find Element Method") { std::shared_ptr xdg = XDG::create(MeshLibrary::MFEM); From a7440a6bf9de95a0ea52740e27ce0384202ddf99 Mon Sep 17 00:00:00 2001 From: Sean Baccas Date: Fri, 26 Jun 2026 12:13:22 +0100 Subject: [PATCH 08/19] Minor updates after rebase --- CMakeLists.txt | 2 +- include/xdg/mfem/mesh_manager.h | 4 ++++ tools/CMakeLists.txt | 1 - tools/mfem_tool.cpp | 36 --------------------------------- 4 files changed, 5 insertions(+), 38 deletions(-) delete mode 100644 tools/mfem_tool.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 146e7e20..2d35494e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -392,7 +392,7 @@ if (XDG_ENABLE_MOAB) endif() if (XDG_ENABLE_MFEM) - target_link_libraries(xdg mfem) + target_link_libraries(xdg PUBLIC mfem) endif() #================================================================= diff --git a/include/xdg/mfem/mesh_manager.h b/include/xdg/mfem/mesh_manager.h index b1bba370..ac7ea8bf 100644 --- a/include/xdg/mfem/mesh_manager.h +++ b/include/xdg/mfem/mesh_manager.h @@ -85,6 +85,10 @@ class MfemMeshManager : public MeshManager { // 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"); } diff --git a/tools/CMakeLists.txt b/tools/CMakeLists.txt index 79949eec..50a0c650 100644 --- a/tools/CMakeLists.txt +++ b/tools/CMakeLists.txt @@ -7,7 +7,6 @@ point_in_volume overlap_check walk_elements tally_segments -mfem_tool ) #=============================================================================== diff --git a/tools/mfem_tool.cpp b/tools/mfem_tool.cpp deleted file mode 100644 index 512cf935..00000000 --- a/tools/mfem_tool.cpp +++ /dev/null @@ -1,36 +0,0 @@ -#include -#include -#include - -#include "xdg/error.h" -#include "xdg/mesh_manager_interface.h" -#include "xdg/mesh_managers.h" -#include "xdg/vec3da.h" -#include "xdg/xdg.h" - -#include "argparse/argparse.hpp" - -#include "particle_sim.h" - -using namespace xdg; - -int main(int argc, char** argv) { - - std::unique_ptr mesh_manager = std::make_unique(); - argparse::ArgumentParser args("MFEM debugging tool", "1.0", argparse::default_arguments::help); - - args.add_argument("filename").help("Path to the input file"); - - try { - args.parse_args(argc, argv); - } - catch (const std::runtime_error& err) { - std::cout << err.what() << std::endl; - std::cout << args; - exit(0); - } - - mesh_manager->load_file(args.get("filename")); - mesh_manager->init(); - -} From 59861d63c68a97d9a0fce76a8c32ce5e0415e495 Mon Sep 17 00:00:00 2001 From: Sean Baccas Date: Wed, 8 Jul 2026 12:23:08 +0100 Subject: [PATCH 09/19] Getting surface senses working --- .gitignore | 3 +- include/xdg/libmesh/mesh_manager.h | 2 +- include/xdg/mfem/mesh_manager.h | 31 ++++--- src/mfem/mesh_manager.cpp | 141 +++++++++++++++++++++++++---- tests/test_mfem.cpp | 80 ++++++++++++++++ tests/util.h | 7 ++ 6 files changed, 231 insertions(+), 33 deletions(-) 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/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/mfem/mesh_manager.h b/include/xdg/mfem/mesh_manager.h index ac7ea8bf..fcb05eda 100644 --- a/include/xdg/mfem/mesh_manager.h +++ b/include/xdg/mfem/mesh_manager.h @@ -31,13 +31,12 @@ class MfemMeshManager : public MeshManager { // Interface methods MeshLibrary mesh_library() const override { return MeshLibrary::MFEM; } - // This info might not be available in mfem int num_volumes() const override { - return mesh_->attribute_sets.GetAttributeSetNames().size(); + return volumes_.size(); } int num_surfaces() const override { - return mesh_->bdr_attribute_sets.GetAttributeSetNames().size(); + return surfaces_.size(); } int num_ents_of_dimension(int dim) const override { @@ -75,11 +74,9 @@ class MfemMeshManager : public MeshManager { virtual std::vector get_surface_faces(MeshID surface) const override; - // see Mesh::GetElementVertices virtual std::vector element_vertices(MeshID element) const override; std::vector bdr_element_vertices(MeshID element) const; - // this one is very easy - Mesh::GetFaceVertices returns the coords of face i at the elment level virtual std::array face_vertices(MeshID element) const override; // The table works wonders for this @@ -93,7 +90,6 @@ class MfemMeshManager : public MeshManager { fatal_error("MfemMeshManager::surface_sense() not implemented yet"); } - // mesh_->GetElement(0)->GetGeometryType() virtual SurfaceElementType get_surface_element_type(MeshID element) const override; virtual int num_vertices() const override { @@ -126,13 +122,12 @@ class MfemMeshManager : public MeshManager { std::pair surface_senses(MeshID surface) const override; // Seems like it's only used to create the implicit complement - MeshID create_volume() override { - fatal_error("MfemMeshManager::create_volume() not implemented yet"); - } + MeshID create_volume() override; - void add_surface_to_volume(MeshID volume, MeshID surface, Sense sense, bool overwrite=false) override { - fatal_error("MfemMeshManager::add_surface_to_volume() not implemented yet"); - } + 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 { @@ -152,6 +147,7 @@ class MfemMeshManager : public MeshManager { 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 @@ -160,10 +156,21 @@ class MfemMeshManager : public MeshManager { // 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 { diff --git a/src/mfem/mesh_manager.cpp b/src/mfem/mesh_manager.cpp index 92cfc0e2..f47b8c9b 100644 --- a/src/mfem/mesh_manager.cpp +++ b/src/mfem/mesh_manager.cpp @@ -30,6 +30,10 @@ void MfemMeshManager::init() { volume_to_element_map_[volume_id].insert(i); } + // 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); @@ -49,6 +53,9 @@ void MfemMeshManager::init() { 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've read in the mesh and counted all the attributes, i.e. a unique @@ -57,21 +64,45 @@ void MfemMeshManager::init() { // 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 calling + // 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()) { - std::ostringstream output; - output << "Couldn't find volume " << volume << "\n"; - fatal_error(output.str()); + // 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; } - std::vector 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... @@ -89,13 +120,15 @@ SurfaceElementType MfemMeshManager::get_surface_element_type(MeshID element) con // Should return all of the sidesets that are a part of this volume std::vector MfemMeshManager::get_volume_surfaces(MeshID volume) const { - // get the set associated with this volume - const std::set& sidesets = volumes_to_sidesets_.at(volume); - - // create a vector from this set - std::vector output(sidesets.begin(), sidesets.end()); - - return output; + // 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 { @@ -132,8 +165,7 @@ std::array MfemMeshManager::face_vertices(MeshID element) const { } else { - // create an mfem array to be passed into Mesh::GetFaceVertices. - // this gets populated with the indices of the vertices itself + // index_array gets populated with the indices of the vertices itself mesh_->GetFaceVertices(element, index_array); } @@ -147,15 +179,11 @@ std::array MfemMeshManager::face_vertices(MeshID element) const { } std::pair MfemMeshManager::surface_senses(MeshID surface) const { - // I am trying to get the raytracer preparation routines working with - // the jezebel, so just return {-1, 1}. i.e. implicit_complement, interior_volume. - // Even though we haven't created implicit_complement yet. - warning("MfemMeshManager::surface_senses() is hardcoded for single-volume meshes"); // 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 {1,2}; + return surface_senses_.at(surface); } std::vector MfemMeshManager::element_vertices(MeshID element) const { @@ -209,6 +237,81 @@ MeshID MfemMeshManager::adjacent_element(MeshID element, int face) const { return faces[face]; } +// TODO: Mesh::GetFaceElements or Mesh::GetFaceInformation are what you need +// if +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}; + } +} + + // helper function to convert mfem's element types to xdg VolumeElementType GetTypeFromMfem( mfem::Element::Type t ) { switch (t) { diff --git a/tests/test_mfem.cpp b/tests/test_mfem.cpp index 32561999..d1a0eb8a 100644 --- a/tests/test_mfem.cpp +++ b/tests/test_mfem.cpp @@ -92,3 +92,83 @@ TEST_CASE("TEST Ray Fire Brick") 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()); + + // MFEM does not capture the right metadata for all this stuff!!! + // + // + // 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"); + // } + // } +} 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; From 0baa5b9d976d2df17a064799907d93165a643831 Mon Sep 17 00:00:00 2001 From: Sean Baccas Date: Wed, 8 Jul 2026 17:12:56 +0100 Subject: [PATCH 10/19] Get walk_elements working --- src/mfem/mesh_manager.cpp | 18 ++- tests/test_mfem.cpp | 264 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 279 insertions(+), 3 deletions(-) diff --git a/src/mfem/mesh_manager.cpp b/src/mfem/mesh_manager.cpp index f47b8c9b..e35ae077 100644 --- a/src/mfem/mesh_manager.cpp +++ b/src/mfem/mesh_manager.cpp @@ -232,9 +232,21 @@ MeshID MfemMeshManager::adjacent_element(MeshID element, int face) const { mesh_->GetElementFaces(element, faces, ori); - // face is in range [0,3). So we just need the one - // that the caller asked for - return faces[face]; + 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!"); } // TODO: Mesh::GetFaceElements or Mesh::GetFaceInformation are what you need diff --git a/tests/test_mfem.cpp b/tests/test_mfem.cpp index d1a0eb8a..9bfba757 100644 --- a/tests/test_mfem.cpp +++ b/tests/test_mfem.cpp @@ -172,3 +172,267 @@ TEST_CASE("Test Cylinder-Brick Initialization") // } // } } + +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)); +} + From 2d0eaa0badb30a8e0c41b53f48e152b611bcf892 Mon Sep 17 00:00:00 2001 From: Sean Baccas Date: Mon, 20 Jul 2026 13:51:55 +0100 Subject: [PATCH 11/19] Add metadata parseing --- include/xdg/mfem/mesh_manager.h | 4 +-- src/mfem/mesh_manager.cpp | 35 ++++++++++++++++++++++ tests/test_mfem.cpp | 53 ++++++++++++++++----------------- 3 files changed, 61 insertions(+), 31 deletions(-) diff --git a/include/xdg/mfem/mesh_manager.h b/include/xdg/mfem/mesh_manager.h index fcb05eda..a10e4348 100644 --- a/include/xdg/mfem/mesh_manager.h +++ b/include/xdg/mfem/mesh_manager.h @@ -130,9 +130,7 @@ class MfemMeshManager : public MeshManager { void determine_surface_senses(); // Metadata methods - void parse_metadata() override { - fatal_error("MfemMeshManager::parse_metadata() not implemented yet"); - } + void parse_metadata() override; // Accessors const std::unique_ptr& mfem_mesh() const { diff --git a/src/mfem/mesh_manager.cpp b/src/mfem/mesh_manager.cpp index e35ae077..4040b5ec 100644 --- a/src/mfem/mesh_manager.cpp +++ b/src/mfem/mesh_manager.cpp @@ -1,6 +1,7 @@ #include #include "xdg/mfem/mesh_manager.h" +#include "xdg/util/str_utils.h" namespace xdg { void MfemMeshManager::load_file(const std::string &filepath) { @@ -323,6 +324,40 @@ void MfemMeshManager::add_surface_to_volume(MeshID volume, MeshID surface, Sense } } +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}; + } + } +} + // helper function to convert mfem's element types to xdg VolumeElementType GetTypeFromMfem( mfem::Element::Type t ) { diff --git a/tests/test_mfem.cpp b/tests/test_mfem.cpp index 9bfba757..77996c75 100644 --- a/tests/test_mfem.cpp +++ b/tests/test_mfem.cpp @@ -143,34 +143,31 @@ TEST_CASE("Test Cylinder-Brick Initialization") auto vol2_elems = mesh_manager->get_volume_elements(2); REQUIRE(!vol2_elems.empty()); - // MFEM does not capture the right metadata for all this stuff!!! - // - // - // 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"); - // } - // } + 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") From 26e2e0b978b18a3130b279a034d5f66e008aad68 Mon Sep 17 00:00:00 2001 From: Sean Baccas Date: Thu, 6 Aug 2026 12:40:11 +0100 Subject: [PATCH 12/19] Add in cross checks --- src/mfem/mesh_manager.cpp | 43 ++++++++++- tests/test_mesh_library_cross_check.cpp | 95 ++++++++++++++++++++++++- tests/test_mfem.cpp | 19 +++++ 3 files changed, 151 insertions(+), 6 deletions(-) diff --git a/src/mfem/mesh_manager.cpp b/src/mfem/mesh_manager.cpp index 4040b5ec..c1352383 100644 --- a/src/mfem/mesh_manager.cpp +++ b/src/mfem/mesh_manager.cpp @@ -26,11 +26,21 @@ void MfemMeshManager::init() { // 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; @@ -58,6 +68,35 @@ void MfemMeshManager::init() { // 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 @@ -81,7 +120,7 @@ void MfemMeshManager::init() { std::vector MfemMeshManager::get_volume_elements(MeshID volume) const { std::vector output; - // Copy the way that libmesh does it, which is if we calling + // 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 @@ -250,8 +289,6 @@ MeshID MfemMeshManager::adjacent_element(MeshID element, int face) const { fatal_error("Shouldn't reach this far!"); } -// TODO: Mesh::GetFaceElements or Mesh::GetFaceInformation are what you need -// if void MfemMeshManager::determine_surface_senses() { for (auto &[surface_id, surface_faces] : sideset_to_element_map_) { if (surface_faces.size() == 0) continue; 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 index 77996c75..be209f27 100644 --- a/tests/test_mfem.cpp +++ b/tests/test_mfem.cpp @@ -433,3 +433,22 @@ TEST_CASE("Test Track Exiting Mesh Brick") 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); +// } + From c2b1d8e5ba595725376f9aa287690f2dabf8a3d9 Mon Sep 17 00:00:00 2001 From: Sean Baccas Date: Thu, 6 Aug 2026 17:35:20 +0100 Subject: [PATCH 13/19] Bringing branch up to speed with changes relating to hex/tet support --- include/xdg/mfem/mesh_manager.h | 41 +++++++++++++-------- src/mfem/mesh_manager.cpp | 64 ++++++++++++++++++++------------- tests/test_mfem.cpp | 4 +-- 3 files changed, 69 insertions(+), 40 deletions(-) diff --git a/include/xdg/mfem/mesh_manager.h b/include/xdg/mfem/mesh_manager.h index a10e4348..b9704954 100644 --- a/include/xdg/mfem/mesh_manager.h +++ b/include/xdg/mfem/mesh_manager.h @@ -77,7 +77,11 @@ class MfemMeshManager : public MeshManager { virtual std::vector element_vertices(MeshID element) const override; std::vector bdr_element_vertices(MeshID element) const; - virtual std::array face_vertices(MeshID element) const override; + 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; @@ -90,8 +94,6 @@ class MfemMeshManager : public MeshManager { fatal_error("MfemMeshManager::surface_sense() not implemented yet"); } - virtual SurfaceElementType get_surface_element_type(MeshID element) const override; - virtual int num_vertices() const override { return mesh_->GetNV(); } @@ -188,13 +190,11 @@ struct MfemMeshElementFaceAccessor : public ElementFaceAccessor { // pointer to the element object that defines this face auto face_obj = mesh->GetFace(face_no); + face_obj->GetVertices(vertex_indices_); - mfem::Array vertex_indices; - face_obj->GetVertices(vertex_indices); - - for (int v=0; vGetVertex( vertex_indices[v] ); + const double* vertices = mesh->GetVertex( vertex_indices_[v] ); for (int d=0; dSpaceDimension(); d++) face_vertices_[f][v][d] = vertices[d]; } @@ -207,12 +207,12 @@ struct MfemMeshElementFaceAccessor : public ElementFaceAccessor { // 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::array face_vertices(int i) const override { - std::array verts; + 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], verts.begin()); + std::copy(face_vertices_[i], face_vertices_[i+1], std::back_inserter(output)); // we need mesh_->GetFaceElementTransformations auto& mesh = mesh_manager_->mfem_mesh(); @@ -229,9 +229,20 @@ struct MfemMeshElementFaceAccessor : public ElementFaceAccessor { 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( verts[0], verts[1] ); + 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_); - return verts; + switch (geom) { + case VolumeElementType::TET: return 4; + case VolumeElementType::HEX: return 6; + default: fatal_error(""); + } } // data members @@ -241,12 +252,14 @@ struct MfemMeshElementFaceAccessor : public ElementFaceAccessor { 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 ); +// SurfaceElementType GetSurfaceElementTypeFromMfem( mfem::Element::Type t ); } // namespace xdg diff --git a/src/mfem/mesh_manager.cpp b/src/mfem/mesh_manager.cpp index c1352383..50c03fc0 100644 --- a/src/mfem/mesh_manager.cpp +++ b/src/mfem/mesh_manager.cpp @@ -153,11 +153,6 @@ std::vector MfemMeshManager::get_volume_elements(MeshID volume) const { return output; } -SurfaceElementType MfemMeshManager::get_surface_element_type(MeshID element) const { - auto mfem_element_type = mesh_->GetBdrElement(element)->GetType(); - return GetSurfaceElementTypeFromMfem(mfem_element_type); -} - // 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 @@ -189,8 +184,7 @@ std::vector MfemMeshManager::get_surface_faces(MeshID surface) const { return output; } -std::array MfemMeshManager::face_vertices(MeshID element) const { - std::array output; +std::vector MfemMeshManager::face_vertices(MeshID element) const { mfem::Array index_array; if (element >= num_interior_faces_) { @@ -209,12 +203,14 @@ std::array MfemMeshManager::face_vertices(MeshID element) const { mesh_->GetFaceVertices(element, index_array); } - for (int i=0; iGetVertex( index_array[i] ); + // for (int i=0; iGetVertex( index_array[i] ); - for (int d=0; dSpaceDimension(); d++) output[i][d] = vertices[d]; - } + // 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; } @@ -395,27 +391,47 @@ void MfemMeshManager::parse_metadata() { } } +// 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); -// 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; + switch(geom) { + case mfem::Geometry::TRIANGLE: return SurfaceFaceType::TRI; + case mfem::Geometry::TETRAHEDRON : return SurfaceFaceType::QUAD; default: - fatal_error("Unsupported element type\n"); + 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; } -// this second function is somewhat redundant. The mfem enum captures all -// of the possible geometries, in all possible dimensions... -SurfaceElementType GetSurfaceElementTypeFromMfem( mfem::Element::Type t ) { + +// helper function to convert mfem's element types to xdg +VolumeElementType GetTypeFromMfem( mfem::Element::Type t ) { switch (t) { - case mfem::Element::TRIANGLE: return SurfaceElementType::TRI; - case mfem::Element::QUADRILATERAL: return SurfaceElementType::QUAD; + case mfem::Element::TETRAHEDRON: return VolumeElementType::TET; + case mfem::Element::HEXAHEDRON: return VolumeElementType::HEX; default: fatal_error("Unsupported element type\n"); } } - -} // namespace xdg \ No newline at end of file +} // namespace xdg diff --git a/tests/test_mfem.cpp b/tests/test_mfem.cpp index be209f27..d0c58be2 100644 --- a/tests/test_mfem.cpp +++ b/tests/test_mfem.cpp @@ -40,8 +40,8 @@ TEST_CASE("MFEM element types") // 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 ); + // for (int i=0; inum_boundary_elements(); i++) + // REQUIRE( mesh_manager->get_surface_element_type(i) == SurfaceElementType::TRI ); } From c7fd1db7c8bf6d3f79b5ab14cbf50305ff5afe73 Mon Sep 17 00:00:00 2001 From: Sean Baccas Date: Tue, 18 Aug 2026 15:55:22 +0100 Subject: [PATCH 14/19] See if the ci will build mfem --- .github/workflows/ci.yml | 43 +++++++++++++++++++++++++++++++++++++++- 1 file changed, 42 insertions(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 1ac4d8a2..47a0e630 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -10,13 +10,14 @@ 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 @@ -81,6 +82,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 +107,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 +132,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.LIBMESH_SHA }} + - if: ${{ matrix.moab && steps.moab-cache.outputs.cache-hit != 'true' }} name: Build MOAB run: | @@ -141,6 +166,22 @@ 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 + cmake .. -DMFEM_USE_HYPRE=OFF -DMFEM_USE_NETCDF=ON -DCMAKE_INSTALL_PREFIX=../install -DCMAKE_BUILD_TYPE=Debug + make -j4 + sudo make install + - name: Build shell: bash run: | From d88272416ebfa2d5de335662577c899318a33cf2 Mon Sep 17 00:00:00 2001 From: Sean Baccas Date: Tue, 18 Aug 2026 16:09:14 +0100 Subject: [PATCH 15/19] Try and find the hdf5 dir --- .github/workflows/ci.yml | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 47a0e630..193378e8 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -173,12 +173,14 @@ jobs: # first let's get netcdf sudo apt-get update sudo apt-get install -y libnetcdf-dev libhdf5-dev + HDF5_DIR="$(dirname "$(dpkg -L libhdf5-dev | grep '/hdf5-config.cmake$' | head -1)")" + echo "HDF5_DIR=$HDF5_DIR" cd ~ cd mfem git checkout v4.9 mkdir build cd build - cmake .. -DMFEM_USE_HYPRE=OFF -DMFEM_USE_NETCDF=ON -DCMAKE_INSTALL_PREFIX=../install -DCMAKE_BUILD_TYPE=Debug + 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 make -j4 sudo make install From 7d3cb52e2eca47a6de3cb6845be7cfcebdcb1001 Mon Sep 17 00:00:00 2001 From: Sean Baccas Date: Wed, 19 Aug 2026 10:58:26 +0100 Subject: [PATCH 16/19] Finish building with mfem --- .github/workflows/ci.yml | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 193378e8..368418cf 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -173,14 +173,13 @@ jobs: # first let's get netcdf sudo apt-get update sudo apt-get install -y libnetcdf-dev libhdf5-dev - HDF5_DIR="$(dirname "$(dpkg -L libhdf5-dev | grep '/hdf5-config.cmake$' | head -1)")" - echo "HDF5_DIR=$HDF5_DIR" cd ~ cd mfem git checkout v4.9 mkdir build cd build - 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 + # 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 @@ -200,7 +199,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 From f40b139c37e52cb568307466cc9aa7f7aaf4736e Mon Sep 17 00:00:00 2001 From: Sean Baccas Date: Wed, 19 Aug 2026 13:22:10 +0100 Subject: [PATCH 17/19] Comment out mfem cache job. Hopefully forces it to rebuild mfem --- .github/workflows/ci.yml | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 368418cf..9996e3f2 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -132,15 +132,15 @@ 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.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 From e9a31fea4d64c4b05e84c762e2f3f71f72eb4044 Mon Sep 17 00:00:00 2001 From: Sean Baccas Date: Wed, 19 Aug 2026 15:18:28 +0100 Subject: [PATCH 18/19] Modify tests --- .github/workflows/ci.yml | 1 + tests/CMakeLists.txt | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 9996e3f2..4a3065e2 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -21,6 +21,7 @@ jobs: exclude: - moab: false libmesh: false + mfem: false steps: - name: Checkout diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 4af15d5d..f6d98ee9 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -35,7 +35,7 @@ if (XDG_ENABLE_MFEM) list(APPEND TEST_NAMES test_mfem) endif() -if (XDG_ENABLE_MOAB AND XDG_ENABLE_LIBMESH) +if (XDG_ENABLE_MOAB AND XDG_ENABLE_LIBMESH and XDG_ENABLE_MFEM) list(APPEND TEST_NAMES test_mesh_library_cross_check) endif() From cb486d52ef75d2f66178c94a2569ebfcb3ed8c59 Mon Sep 17 00:00:00 2001 From: Sean Baccas Date: Wed, 2 Sep 2026 15:30:04 +0100 Subject: [PATCH 19/19] Fix CI --- CMakeLists.txt | 3 ++- tests/CMakeLists.txt | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 2d35494e..8a09030f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -178,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() diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index f6d98ee9..45c6b7fa 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -35,7 +35,7 @@ if (XDG_ENABLE_MFEM) list(APPEND TEST_NAMES test_mfem) endif() -if (XDG_ENABLE_MOAB AND XDG_ENABLE_LIBMESH and XDG_ENABLE_MFEM) +if (XDG_ENABLE_MOAB AND XDG_ENABLE_LIBMESH AND XDG_ENABLE_MFEM) list(APPEND TEST_NAMES test_mesh_library_cross_check) endif()