tboudreaux/MeanField_SERiF
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

feat(discretization/domain): updated domain system to use unified mfem backend

Browse Source 
This means that the domain system never directly uses mfem code any more
This commit is contained in:
Emily Boudreaux
2026-10-14 11:01:33 -04:00
parent d2e43c82f1
commit 7f3ebbae45
2 changed files with 43 additions and 51 deletions
Download Patch File Download Diff File Expand all files Collapse all files

37
src/include/serif/discretization/domain/mesh/topology.hpp
View File

@@ -10,17 +10,6 @@
#include <cstddef>
#include <vector>
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<int> elementDomainIDs;
std::vector<FaceElements> faceElements;
std::vector<int> boundaryElementBoundaryIDs;
std::vector<int> 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;

57
src/lib/serif/discretization/domain/mesh/topology.cpp
View File

@@ -1,45 +1,38 @@
#include "serif/discretization/domain/mesh/topology.hpp"
#include <mfem.hpp>
#include <utility>
#include "serif/utils/error/errors.hpp"
namespace serif::discretization::domain::mesh {
const std::vector<int> 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<std::size_t>(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<SERiFErrorCode::INVALID_ARGUMENT>(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<SERiFErrorCode::OUT_OF_RANGE_ERROR>((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<SERiFErrorCode::OUT_OF_RANGE_ERROR>(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<int>(boundaryElementID));
}
}
}