From 54dcd9f3f9ff903a22e297ce604fab53507f0d4b Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Mon, 5 Oct 2026 13:22:26 -0400 Subject: [PATCH] feat(discretiztion/backend/mfem): added boundary_trace and mass_diagonal backend --- .../backend/mfem/boundary_trace.hpp | 95 +++++++++++++++++++ .../backend/mfem/mass_diagonal.hpp | 28 ++++++ 2 files changed, 123 insertions(+) create mode 100644 src/include/serif/discretization/backend/mfem/boundary_trace.hpp create mode 100644 src/include/serif/discretization/backend/mfem/mass_diagonal.hpp diff --git a/src/include/serif/discretization/backend/mfem/boundary_trace.hpp b/src/include/serif/discretization/backend/mfem/boundary_trace.hpp new file mode 100644 index 00000000..50092e0f --- /dev/null +++ b/src/include/serif/discretization/backend/mfem/boundary_trace.hpp @@ -0,0 +1,95 @@ +#pragma once + +/** + * @file boundary_trace.hpp + * @author Emily M. Boudreaux + * + * @brief File containing the definition of the restricted dof space for a scalar boundary trace. Note in FEM a trace is a restriction on a function to the boundary of a domain. Therefore, + * This file contains the definitions for getting the DOFs at the boundary of a scalar field. + */ + +#include +#include + +#include "mfem/linalg/sparsemat.hpp" +#include "serif/discretization/backend/mfem/dof_selection.hpp" + +namespace serif::discretization::backend::mfem +{ + /** + * @brief Concept Limiting the input term to either a global scalar term or a spatial term of rank 0 in the H1 family. Outside of H1, the trace is a vector (e.g. the normal trace + * for a H(div) field or the tangential trace for a H(curl) field). Further, for higher + * rank H1 fields the trace is a tensor. + **/ + template + concept BoundaryTraceTerm = blocks::GlobalScalarTerm || (blocks::SpatialTerm && blocks::SpatialTermTraits::rank == 0 && blocks::SpatialTermTraits::family == functional_space::Family::h1); + + /** + * @todo We have this forward definition here to avoid a circular dependency; however, a better refactor should be made here to avoid needing this forward declaration. + */ + class ScalarBoundaryTrace; + + namespace detail + { + /** @brief Construct one interpolation operator in the selected residual-node ordering. */ + [[nodiscard]] ScalarBoundaryTrace build_boundary_trace(std::shared_ptr rows, const DofMap& input, const ::mfem::FiniteElementSpace* input_space); + } + + /** + * @brief An immutable scalar interpolation operator from supported coefficients to selected residual nodes. + * @note Both term bindings must outlive the trace. The sparse matrix and selection are shared immutable data; + * the same matrix supplies the forward evaluation and its exact coordinate transpose without quadrature weights. + */ + class ScalarBoundaryTrace final + { + public: + /** @brief Return the borrowed input coefficient binding. */ + [[nodiscard]] const DofMap& binding() const noexcept; + + /** @brief Return the owned selection defining the interpolation row ordering. */ + [[nodiscard]] const DofSelection& rows() const noexcept; + + /** @brief Return the interpolation matrix after checking its borrowed discretization. */ + [[nodiscard]] const ::mfem::SparseMatrix& interpolation() const; + + /** @brief Reject stale mesh topology, finite-element bindings, or domain and boundary labels. */ + void validate() const; + + private: + /** @brief Publish the completed matrix together with the bindings defining its coordinates. */ + ScalarBoundaryTrace(std::shared_ptr rows, const DofMap& input, std::shared_ptr interpolation); + + friend ScalarBoundaryTrace detail::build_boundary_trace(std::shared_ptr, const DofMap&, const ::mfem::FiniteElementSpace*); + + std::shared_ptr m_rows; + const DofMap& m_binding; + std::shared_ptr m_interpolation; + }; + + /** + * @brief Sample a continuous scalar field at the selected nodal H1 residual points. + * @note Input and residual spaces may have different orders and numbering but must share a mesh. + * Discontinuous fields require an explicit physical trace convention and are intentionally excluded. + */ + template + [[nodiscard]] ScalarBoundaryTrace make_boundary_trace(std::shared_ptr rows, const TermBinding& input) + { + if constexpr (blocks::GlobalScalarTerm) + { + return detail::build_boundary_trace(std::move(rows), input, nullptr); + } + else + { + return detail::build_boundary_trace(std::move(rows), input, &input.space()); + } + } + + /** @brief A trace borrows its input binding, so temporaries cannot supply its coefficient ordering. */ + template + ScalarBoundaryTrace make_boundary_trace(std::shared_ptr, TermBinding&&) = delete; + + /** @brief Const temporary bindings have the same lifetime restriction. */ + template + ScalarBoundaryTrace make_boundary_trace(std::shared_ptr, const TermBinding&&) = delete; + +} diff --git a/src/include/serif/discretization/backend/mfem/mass_diagonal.hpp b/src/include/serif/discretization/backend/mfem/mass_diagonal.hpp new file mode 100644 index 00000000..e8866be5 --- /dev/null +++ b/src/include/serif/discretization/backend/mfem/mass_diagonal.hpp @@ -0,0 +1,28 @@ +#pragma once + +#include +#include +#include + +#include "serif/discretization/backend/mfem/term_binding.hpp" + +namespace serif::discretization::backend::mfem +{ + namespace detail + { + /** @brief Assemble a supported reference mass diagonal for scalar copies or intrinsic H(div) vectors. */ + [[nodiscard]] std::expected, std::string> reference_mass_diagonal(const DofMap& binding, bool intrinsic_vector); + } + + /** + * @brief Measure the supported coefficients in their reference finite-element geometry. + * @return The exact mass-matrix diagonal in the binding's reduced ordering, or an invalid-metric diagnostic. + * @note The binding must remain current. Coordinates belong to the reference mesh, not a trial deformation. + */ + template + [[nodiscard]] std::expected, std::string> reference_mass_diagonal(const TermBinding& binding) + { + constexpr bool intrinsic_vector = blocks::SpatialTermTraits::family == functional_space::Family::hdiv; + return detail::reference_mass_diagonal(binding, intrinsic_vector); + } +}