feat(discretiztion/backend/mfem): added boundary_trace and mass_diagonal backend

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 <memory>
+#include <utility>
+
+#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 <typename TermT>
+    concept BoundaryTraceTerm = blocks::GlobalScalarTerm<TermT> || (blocks::SpatialTerm<TermT> && blocks::SpatialTermTraits<TermT>::rank == 0 && blocks::SpatialTermTraits<TermT>::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<const DofSelection> 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<const DofSelection> rows, const DofMap& input, std::shared_ptr<const ::mfem::SparseMatrix> interpolation);
+
+        friend ScalarBoundaryTrace detail::build_boundary_trace(std::shared_ptr<const DofSelection>, const DofMap&, const ::mfem::FiniteElementSpace*);
+
+        std::shared_ptr<const DofSelection> m_rows;
+        const DofMap& m_binding;
+        std::shared_ptr<const ::mfem::SparseMatrix> 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 <BoundaryTraceTerm TermT, domain::schema::IsSchema SchemaT>
+    [[nodiscard]] ScalarBoundaryTrace make_boundary_trace(std::shared_ptr<const DofSelection> rows, const TermBinding<TermT, SchemaT>& input)
+    {
+        if constexpr (blocks::GlobalScalarTerm<TermT>)
+        {
+            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 <BoundaryTraceTerm TermT, domain::schema::IsSchema SchemaT>
+    ScalarBoundaryTrace make_boundary_trace(std::shared_ptr<const DofSelection>, TermBinding<TermT, SchemaT>&&) = delete;
+
+    /** @brief Const temporary bindings have the same lifetime restriction. */
+    template <BoundaryTraceTerm TermT, domain::schema::IsSchema SchemaT>
+    ScalarBoundaryTrace make_boundary_trace(std::shared_ptr<const DofSelection>, const TermBinding<TermT, SchemaT>&&) = delete;
+
+}

src/include/serif/discretization/backend/mfem/mass_diagonal.hpp

@@ -0,0 +1,28 @@
+#pragma once
+
+#include <expected>
+#include <string>
+#include <vector>
+
+#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::vector<double>, 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 <blocks::SpatialTerm Term, domain::schema::IsSchema Schema>
+    [[nodiscard]] std::expected<std::vector<double>, std::string> reference_mass_diagonal(const TermBinding<Term, Schema>& binding)
+    {
+        constexpr bool intrinsic_vector = blocks::SpatialTermTraits<Term>::family == functional_space::Family::hdiv;
+        return detail::reference_mass_diagonal(binding, intrinsic_vector);
+    }
+}