From 7f3ebbae45c59923471b42560b634f42dd9733fa Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Wed, 14 Oct 2026 11:01:33 -0400 Subject: [PATCH] feat(discretization/domain): updated domain system to use unified mfem backend This means that the domain system never directly uses mfem code any more --- .../discretization/domain/mesh/topology.hpp | 37 ++++++------ .../discretization/domain/mesh/topology.cpp | 57 ++++++++----------- 2 files changed, 43 insertions(+), 51 deletions(-) diff --git a/src/include/serif/discretization/domain/mesh/topology.hpp b/src/include/serif/discretization/domain/mesh/topology.hpp index d016b5c2..3ad3e8fb 100644 --- a/src/include/serif/discretization/domain/mesh/topology.hpp +++ b/src/include/serif/discretization/domain/mesh/topology.hpp @@ -10,17 +10,6 @@ #include #include -namespace mfem { - /** - * @brief Forward declaration of the mfem::Mesh class to avoid including the entire MFEM header. - * - * @todo We should probably consider doing this more often and in some central location since ideally we would avoid including MFEM headers in our own headers. Due to how - * heavily templated our code is that wont always be possible, but there are many places where we could. This has the potential to dramatically reduce compile times - * since our compile time is primarily limited by parsing through the MFEM headers at the moment. - */ - class Mesh; -} - // MeshTopology namespace serif::discretization::domain::mesh { /** @@ -32,6 +21,20 @@ namespace serif::discretization::domain::mesh { int secondElementID{-1}; }; + /** + * @brief Backend-independent material labels and face connectivity; + * + * @todo Implemented this with std::vector when I moved all the mfem code into backend/mfem; however, + * std::vector has all sorts of limitations and problems. We may want to consider using the + * types.hpp vector alias here instead? Or some other vector implementation. + */ + struct MeshTopologyData { + std::vector elementDomainIDs; + std::vector faceElements; + std::vector boundaryElementBoundaryIDs; + std::vector boundaryElementFaceIDs; + }; + /** * @class MeshTopology * @brief Provides a convenient interface for querying the topological relationships between elements, faces, and boundaries in a mesh. @@ -40,10 +43,10 @@ namespace serif::discretization::domain::mesh { public: /** - * @brief Constructs a MeshTopology object from an MFEM mesh. Note that MeshTopology does not own the mesh; - * @param mesh The MFEM mesh from which to extract topological information. + * @brief Owns material labels and connectivity independently of the mesh backend. + * @param data The topology data from which to derive element and boundary adjacency. */ - explicit MeshTopology(const mfem::Mesh &mesh); + explicit MeshTopology(MeshTopologyData data); /** * @brief Returns the number of elements in the mesh. @@ -68,7 +71,7 @@ namespace serif::discretization::domain::mesh { * @param elementID The ID of the element for which to retrieve the domain ID. * @return The domain ID of the specified element. * - * @note The domain ID is the MFEM element attribute (often called the material attribute) of a volume element, i.e. which DomainID the mesh claims this element belongs to. + * @note The domain ID is the material label of a volume element, i.e. which DomainID the mesh claims this element belongs to. */ [[nodiscard]] int element_domain_id(int elementID) const; @@ -77,7 +80,7 @@ namespace serif::discretization::domain::mesh { * @param boundaryElementID The ID of the boundary element for which to retrieve the boundary ID. * @return The boundary ID of the specified boundary element. * - * @note The boundary ID is the MFEM boundary attribute of a boundary element, i.e. which BoundaryID the mesh claims this boundary element belongs to. + * @note The boundary ID is the boundary label of a boundary element, i.e. which BoundaryID the mesh claims this boundary element belongs to. */ [[nodiscard]] int boundary_element_boundary_id(int boundaryElementID) const; @@ -100,7 +103,7 @@ namespace serif::discretization::domain::mesh { /** * @brief Returns the face a boundary element sits on. * @param boundaryElementID The ID of the boundary element for which to retrieve the face. - * @return The ID of the face the boundary element sits on, or -1 if MFEM reports none. + * @return The ID of the face the boundary element sits on, or -1 if no face was supplied. */ [[nodiscard]] int boundary_element_face(int boundaryElementID) const; diff --git a/src/lib/serif/discretization/domain/mesh/topology.cpp b/src/lib/serif/discretization/domain/mesh/topology.cpp index fa2f90f1..8e07951b 100644 --- a/src/lib/serif/discretization/domain/mesh/topology.cpp +++ b/src/lib/serif/discretization/domain/mesh/topology.cpp @@ -1,45 +1,38 @@ #include "serif/discretization/domain/mesh/topology.hpp" -#include +#include +#include "serif/utils/error/errors.hpp" namespace serif::discretization::domain::mesh { const std::vector MeshTopology::s_noElements{}; namespace { + /** @brief Convert a nonnegative topology coordinate to a container index. */ [[nodiscard]] std::size_t as_index(const int id) noexcept { return static_cast(id); } + /** @brief Test a signed topology coordinate before indexing storage. */ [[nodiscard]] bool in_range(const int id, const std::size_t size) noexcept { return id >= 0 && as_index(id) < size; } } - MeshTopology::MeshTopology(const mfem::Mesh &mesh) - : m_elementDomainIDs(as_index(mesh.GetNE()), 0), - m_faceElements(as_index(mesh.GetNumFaces())), - m_elementNeighbors(as_index(mesh.GetNE())), - m_boundaryElementBoundaryIDs(as_index(mesh.GetNBE()), 0), - m_boundaryElementFaceIDs(as_index(mesh.GetNBE()), -1), - m_boundaryElementsByFace(as_index(mesh.GetNumFaces())) { - - // Which domain each volume element claims to belong to. - for (int elementID = 0; elementID < mesh.GetNE(); ++elementID) { - m_elementDomainIDs[as_index(elementID)] = mesh.GetAttribute(elementID); - } - - // Face connectivity. - for (int faceID = 0; faceID < mesh.GetNumFaces(); ++faceID) { - int firstElementID = -1; - int secondElementID = -1; - - mesh.GetFaceElements(faceID, &firstElementID, &secondElementID); - - m_faceElements[as_index(faceID)] = FaceElements{ - .firstElementID = firstElementID, - .secondElementID = secondElementID - }; + MeshTopology::MeshTopology(MeshTopologyData data) + : m_elementDomainIDs(std::move(data.elementDomainIDs)), + m_faceElements(std::move(data.faceElements)), + m_elementNeighbors(m_elementDomainIDs.size()), + m_boundaryElementBoundaryIDs(std::move(data.boundaryElementBoundaryIDs)), + m_boundaryElementFaceIDs(std::move(data.boundaryElementFaceIDs)), + m_boundaryElementsByFace(m_faceElements.size()) { + using namespace utils::errors; + serif_ensure(m_boundaryElementBoundaryIDs.size() == m_boundaryElementFaceIDs.size(), "Each topology boundary element requires both a boundary ID and a face ID."); + // Face connectivity determines adjacency independently of backend numbering conventions. + for (const FaceElements& face : m_faceElements) { + const int firstElementID = face.firstElementID; + const int secondElementID = face.secondElementID; + serif_ensure((firstElementID == -1 || in_range(firstElementID, m_elementDomainIDs.size())) && (secondElementID == -1 || in_range(secondElementID, m_elementDomainIDs.size())), "Topology face elements must be valid element indices or -1 for the exterior; received {} and {} for {} elements.", firstElementID, secondElementID, m_elementDomainIDs.size()); if (firstElementID < 0 || secondElementID < 0) { continue; } @@ -48,16 +41,12 @@ namespace serif::discretization::domain::mesh { m_elementNeighbors[as_index(secondElementID)].push_back(firstElementID); } - // Boundary elements - for (int boundaryElementID = 0; boundaryElementID < mesh.GetNBE(); ++boundaryElementID) { - m_boundaryElementBoundaryIDs[as_index(boundaryElementID)] = mesh.GetBdrAttribute(boundaryElementID); - - const int faceID = mesh.GetBdrElementFaceIndex(boundaryElementID); - - m_boundaryElementFaceIDs[as_index(boundaryElementID)] = faceID; - + // A face can carry several boundary labels; retain all of them for schema validation. + for (std::size_t boundaryElementID = 0; boundaryElementID < m_boundaryElementFaceIDs.size(); ++boundaryElementID) { + const int faceID = m_boundaryElementFaceIDs[boundaryElementID]; + serif_ensure(faceID == -1 || in_range(faceID, m_faceElements.size()), "Topology boundary face {} must be a valid face index or -1 for an unresolved face.", faceID); if (in_range(faceID, m_boundaryElementsByFace.size())) { - m_boundaryElementsByFace[as_index(faceID)].push_back(boundaryElementID); + m_boundaryElementsByFace[as_index(faceID)].push_back(static_cast(boundaryElementID)); } } }