diff --git a/src/include/serif/discretization/backend/mfem/basis_evaluation.hpp b/src/include/serif/discretization/backend/mfem/basis_evaluation.hpp new file mode 100644 index 00000000..8bb26ee0 --- /dev/null +++ b/src/include/serif/discretization/backend/mfem/basis_evaluation.hpp @@ -0,0 +1,231 @@ +#pragma once + +/** + * @file basis_evaluation.hpp + * @author Emily M. Boudreaux + * @brief Basis evaluation, pointwise mapping, and transpose application using dense or tensor-product representations. + */ +#include + +#include "mfem/fem/fe/fe_base.hpp" +#include "mfem/fem/intrules.hpp" + +#include "serif/discretization/backend/mfem/tensor_basis.hpp" + +namespace serif::discretization::backend::mfem { + /** + * @brief Reusable scratch storage for basis contractions, reducing repeated allocations. + */ + struct BasisWorkspace { + TensorWorkspace tensor; + ::mfem::Vector reference_values; + }; + + /** + * @class BasisEvaluation + * @brief Represents an element-local basis evaluation operator, supporting dense or tensor contractions, pointwise mappings, and transpose application. + * + * Consider a scalar field represented on one finite element: + * + * \f[ + * \rho_h(\xi) = \sum_{j=1}^{n} c_j \phi_j(\xi). + * \f] + * + * Here, \f$c_{j}\f$ is a scalar coefficient, \f$\phi_{j}\f$ is its associated basis function, + * and \f$\xi\f$ is the canonical reference coordinate. The field is expressed in + * these coordinates through the element mapping. + * + * For a quadrature rule with \f$N_{q}\f$ points, define the evaluation matrix \f$\mathbf{E}\f$ by + * + * \f[ + * E_{qj} = \phi_j(\xi_q). + * \f] + * + * \f$\mathbf{E}\f$ therefore has \f$N_{q}\f$ rows and \f$n\f$ columns. The field value at one point is + * + * \f[ + * \rho_h(\xi_q) = \sum_{j=1}^{n} E_{qj} c_j. + * \f] + * + * Collecting the samples into a vector gives + * + * \f[ + * \boldsymbol{\rho}_h = \mathbf{E}\mathbf{c}. + * \f] + * + * This class represents that linear operation without necessarily storing + * the full matrix \f$\mathbf{E}\f$: it can instead apply equivalent tensor contractions. + * + * @par Example + * Evaluate a constant scalar field on the reference cube, then compute + * its integral against each basis function. + * + * @code{.cpp} + * namespace backend = serif::discretization::backend::mfem; + * + * // First-order nodal H1 basis on the reference cube. + * ::mfem::H1_HexahedronElement element(1); + * const ::mfem::IntegrationRule& rule = + * ::mfem::IntRules.Get(::mfem::Geometry::CUBE, 2); + * + * backend::TensorBasisCache cache; + * backend::BasisWorkspace workspace; + * backend::BasisEvaluation basis; + * + * // Prepare E: one scalar basis value per quadrature point and DOF. + * basis.sample( + * element, + * rule, + * backend::ReferenceOperation::value, + * 1, // One copy of the scalar basis. + * cache, + * backend::ContractionMode::automatic + * ); + * + * // For this nodal basis, equal coefficients represent a constant field. + * ::mfem::Vector coefficients(element.GetDof()); + * coefficients = 2.0; + * + * // Compute values = E * coefficients. + * // Every quadrature-point value is approximately 2.0. + * ::mfem::Vector values; + * basis.evaluate(coefficients, values, workspace); + * + * // Supply quadrature weights explicitly: integrate does not add them. + * ::mfem::Vector weighted_values(rule.GetNPoints()); + * for (int q = 0; q < rule.GetNPoints(); ++q) { + * weighted_values[q] = rule.IntPoint(q).weight * values[q]; + * } + * + * // Compute element_integrals = E^T * weighted_values. + * ::mfem::Vector element_integrals; + * basis.integrate(weighted_values, element_integrals, workspace); + * @endcode + * + * Each output entry approximates + * \f[ + * \int_{\hat K} \phi_j(\xi)\rho_h(\xi)\,d\xi. + * \f] + * These are element-DOF integral contributions, not reconstructed field + * coefficients. Integration over a physical element additionally requires + * the appropriate volume Jacobian in the quadrature weights. + */ + class BasisEvaluation final { + public: + /** + * @brief Prepare basis values or derivatives at quadrature points using a dense matrix or tensor-product representation. + * @param element Finite element to sample. + * @param rule Integration rule to use for sampling. + * @param operation Reference operation to perform. + * @param copies Number of independent copies of the finite-element basis; for example, three scalar basis copies for a vector H1 field. + * @param cache Tensor basis cache to use. + * @param mode Contraction mode to use. + */ + void sample( + const ::mfem::FiniteElement& element, + const ::mfem::IntegrationRule& rule, + ReferenceOperation operation, + int copies, + TensorBasisCache& cache, + ContractionMode mode + ); + + /** + * @brief Set the basis evaluation to a constant basis with the given number of points. + * @param points Number of quadrature points to use for the constant basis. + */ + void set_constant( + int points + ); + + /** + * @brief Allocate storage for a mapped basis using the reference's coefficient ordering and quadrature points. + * @param reference Basis evaluation supplying the dimensions and any shared tensor representation. + * @param components Number of output components per quadrature point. + * + * @note Populate the destination at every quadrature point using map_point before evaluation or transpose application. + */ + void reset_from( + const BasisEvaluation& reference, + int components + ); + + /** + * @brief Apply a pointwise linear map to the reference basis or its geometry variation. + * @tparam MapT Type of the mapping function to apply. + * @param point Quadrature point index to apply the mapping at. + * @param reference Reference basis evaluation to map from. + * @param map Mapping function to apply. + * + * @pre The mapping callback must be linear in its input and fill every output component. + */ + template + void map_point( + int point, + const BasisEvaluation& reference, + MapT&& map + ); + + /** + * @brief Apply the basis to native coefficients; no integration weight is introduced. + * @param coefficients Input coefficients to apply the basis to. + * @param values Output values after applying the basis. + * @param workspace Workspace to use for intermediate computations. + */ + void evaluate( + const ::mfem::Vector& coefficients, + ::mfem::Vector& values, + BasisWorkspace& workspace + ) const; + + /** + * @brief Apply the coordinate transpose to already weighted quadrature duals. + * @param[in] values Quadrature data with any required integration weights already applied. + * @param[out] coefficients Element-DOF contributions computed as E^T times values. + * @param workspace Workspace to use for intermediate computations. + */ + void integrate( + const ::mfem::Vector& values, + ::mfem::Vector& coefficients, + BasisWorkspace& workspace + ) const; + + /// Number of flattened quadrature samples. + [[nodiscard]] int height() const noexcept; + + /// Number of native element coefficients. + [[nodiscard]] int width() const noexcept; + + /// Report whether contractions use verified Cartesian factors. + [[nodiscard]] bool uses_tensor() const noexcept; + + private: + int m_points{0}; ///< Number of quadrature points in the rule + int m_components{0}; ///< Number of output components per point + int m_width{0}; ///< Number of native element coefficients + std::shared_ptr m_tensor; ///< Canonical tensor representation, when the tensor path is selected + ::mfem::DenseMatrix m_dense; ///< Full evaluation matrix, including applied mappings, when the dense path is selected + ::mfem::DenseMatrix m_mapping; ///< Stacked pointwise maps from canonical tensor samples to current output components + }; +} // namespace serif::discretization::backend::mfem + +template +void serif::discretization::backend::mfem::BasisEvaluation::map_point(const int point, const BasisEvaluation& reference, MapT&& map) { + const ::mfem::DenseMatrix& source = reference.m_tensor ? reference.m_mapping : reference.m_dense; + + ::mfem::DenseMatrix& destination = m_tensor ? m_mapping : m_dense; + ::mfem::Vector input(reference.m_components); + ::mfem::Vector output(m_components); + + for (int column = 0; column < source.Width(); ++column) { + for (int axis = 0; axis < reference.m_components; ++axis) { + input[axis] = source(point * reference.m_components + axis, column); + } + + map(input, output); + + for (int axis = 0; axis < m_components; ++axis) { + destination(point * m_components + axis, column) = output[axis]; + } + } +} diff --git a/src/include/serif/discretization/backend/mfem/coefficient_tracking.hpp b/src/include/serif/discretization/backend/mfem/coefficient_tracking.hpp new file mode 100644 index 00000000..b821a6b8 --- /dev/null +++ b/src/include/serif/discretization/backend/mfem/coefficient_tracking.hpp @@ -0,0 +1,117 @@ +#pragma once + +/** + * @file coefficient_tracking.hpp + * @brief Snapshot and comparison policy for MFEM coefficient views. + */ + +#include + +#include "mfem/linalg/vector.hpp" + +namespace serif::discretization::backend::mfem { + /** + * @struct CoefficientTracking + * @brief Track coefficient values using an owned snapshot and the identity of their binding. + * + * A coefficient vector describes numerical values, while its binding + * determines how those values correspond to discretization degrees of freedom. + * Both pieces of information matter when deciding whether cached work can + * be reused. + * + * This policy provides three operations to the dependency system: + * + * - capture copies the current coefficient values into an owned snapshot. + * - matches compares current values against a previously captured snapshot. + * - identity identifies the persistent binding associated with the values. + * + * The vector's allocation is not part of its identity. A different vector + * containing equal coefficients can represent the same input when it uses + * the same binding. Conversely, equal coefficients associated with a + * different binding must not automatically reuse the previous calculation. + * + * This policy supplies comparisons; the dependency system decides which + * computations must run. + * + * @note Value comparison uses floating-point equality without a tolerance. + * It does not check whether the coefficients are finite or admissible; + * those checks belong to input validation. + */ + struct CoefficientTracking { + /** + * @brief Copy the current coefficient values into an independent snapshot. + * + * Subsequent changes to the source vector do not modify the returned + * snapshot. The entire supplied vector is copied; if the view contains + * only selected nodal coefficients, only those selected values are captured. + * + * @tparam ViewT View type exposing its coefficients through a values member. + * @param[in] view View of the current coefficients. + * @return An owned copy of the coefficient vector. + */ + template + [[nodiscard]] static ::mfem::Vector capture(const ViewT& view); + + /** + * @brief Check whether the current coefficients equal a captured snapshot. + * + * The vectors must have the same size and equal corresponding entries. + * Two empty vectors match. Their storage addresses do not affect the + * comparison. + * + * This method compares values only. Binding identity is supplied + * separately by identity and checked by the dependency system. + * + * @tparam ViewT View type exposing its coefficients through a values member. + * @param[in] captured Previously captured coefficient values. + * @param[in] view View of the current coefficients. + * @return True when the sizes and all corresponding values are equal. + */ + template + [[nodiscard]] bool matches( + const ::mfem::Vector& captured, + const ViewT& view + ) const; + + /** + * @brief Identify the persistent binding associated with the coefficients. + * + * The returned address identifies the binding object, rather than the + * temporary view or the coefficient vector's storage. Creating another + * view of the same binding therefore preserves this identity. + * + * @tparam ViewT View type exposing its persistent binding through a binding member. + * @param[in] view View referring to the coefficient binding. + * @return The address of the binding object. + * + * @pre The binding member refers to an object whose lifetime and address + * remain valid while the dependency system relies on this identity. + */ + template + [[nodiscard]] static const void* identity(const ViewT& view); + }; +} // namespace serif::discretization::backend::mfem + +template +::mfem::Vector serif::discretization::backend::mfem::CoefficientTracking::capture(const ViewT& view) { + return ::mfem::Vector(view.values); +} + +template +bool serif::discretization::backend::mfem::CoefficientTracking::matches( + const ::mfem::Vector& captured, + const ViewT& view +) const { + return captured.Size() == view.values.Size() && + (captured.Size() == 0 || + std::equal( + captured.GetData(), + captured.GetData() + captured.Size(), + view.values.GetData() + )); +} + +template +[[nodiscard]] const void* serif::discretization::backend::mfem::CoefficientTracking::identity(const ViewT& view) { + return &view.binding; +} diff --git a/src/include/serif/discretization/backend/mfem/dof_map.hpp b/src/include/serif/discretization/backend/mfem/dof_map.hpp new file mode 100644 index 00000000..e7cb36f8 --- /dev/null +++ b/src/include/serif/discretization/backend/mfem/dof_map.hpp @@ -0,0 +1,306 @@ +#pragma once +/** + * @file dof_map.hpp + * @author Emily M. Boudreaux + * @brief This file contains the @ref DofMap which allows for the full + * dof vector to be projected into a reduced space of dofs on supported elements for some field Generally, + * while much of the code here may be exercised, you will rarely see a DofMap be directly invoked. Rather, + * these are generally implicitly constructed through their child class term_binding. + */ + +#include +#include +#include + +#include "mfem/fem/fespace.hpp" +#include "mfem/fem/pfespace.hpp" +#include "serif/utils/error/errors.hpp" + +namespace serif::discretization::backend::mfem { + /** + * @brief A mapping from a full MFEM finite-element space to a reduced set of DOFs that + * are supporter by the given element support. That is to say that if one has some set of + * dofs D (perhaps it looks like [0, 1, 2, 3, 4] as a hypothetical) and dofs 0, 2, and 4 are + * associated to domain A while 1 and 3 are associated to domain B. Then if one were to + * construct a DofMap with element support [1, 0, 1, 0, 1] (where 1 indicates that the element + * is supported and 0 indicates that it is not) then the reduced dof map would be [0, 2, 4] and + * the full to reduced map would be [-1, -1, 0, -1, 1]. + * + * @note This is a key place which will need to be updated if we are to support, as I believe we + * must eventually, non confirming meshes. Further, this current implementation + * does not support multi-rank parallelism in a smart manner. + */ + class DofMap { + protected: + DofMap( + ::mfem::FiniteElementSpace& space, + const ::mfem::Array& element_support + ) : + m_space(&space), + m_mesh(*space.GetMesh()), + m_space_sequence(space.GetSequence()), + m_mesh_sequence(space.GetMesh()->GetSequence()), + m_attributes(space.GetNE()), + m_element_support(element_support), + m_full_to_reduced(space.GetVSize()) { + using namespace utils::errors; + +#ifdef MFEM_USE_MPI + serif_ensure(dynamic_cast<::mfem::ParFiniteElementSpace*>(&space) == nullptr,"DofMap does not yet assemble parallel true DOFs."); +#endif + serif_ensure(space.Conforming() && space.GetTrueVSize() == space.GetVSize(), "Term bindings currently require conforming serial spaces."); + serif_ensure(space.GetMesh()->Dimension() == 3 && space.GetMesh()->SpaceDimension() == 3,"Term bindings currently require a three-dimensional mesh."); + serif_ensure(element_support.Size() == space.GetNE(), "Element support does not match the space."); + + m_full_to_reduced = -1; + ::mfem::Array dofs; + + for (int element = 0; element < space.GetNE(); ++element) { + m_attributes[element] = space.GetMesh()->GetAttribute(element); + + if (element_support[element] == 0) { + continue; + } + + space.GetElementVDofs(element, dofs); + for (const int dof : dofs) { + m_full_to_reduced[dof < 0 ? -1 - dof : dof] = 0; + } + } + + for (int dof = 0; dof < space.GetVSize(); ++dof) { + if (m_full_to_reduced[dof] == 0) { + m_full_to_reduced[dof] = m_reduced_to_full.Size(); + m_reduced_to_full.Append(dof); + } + } + } + + DofMap( + ::mfem::Mesh& mesh, + const ::mfem::Array& element_support + ) : + m_space(nullptr), + m_mesh(mesh), + m_space_sequence(-1), + m_mesh_sequence(mesh.GetSequence()), + m_attributes(mesh.GetNE()), + m_element_support(element_support), + m_reduced_to_full(1), + m_full_to_reduced(1) { + using namespace utils::errors; + + serif_ensure(mesh.Dimension() == 3 && mesh.SpaceDimension() == 3, "Term bindings currently require a three-dimensional mesh."); + serif_ensure(element_support.Size() == mesh.GetNE(), "Element support does not match the mesh."); + + m_reduced_to_full[0] = m_full_to_reduced[0] = 0; + for (int element = 0; element < mesh.GetNE(); ++element) { + m_attributes[element] = mesh.GetAttribute(element); + } + } + + ~DofMap() = default; + + public: + DofMap(const DofMap&) = delete; + DofMap& operator=(const DofMap&) = delete; + DofMap(DofMap&&) = delete; + DofMap& operator=(DofMap&&) = delete; + + /** + * @brief Return the borrowed finite-element space associated with this DOF map. + * @return A reference to the finite-element space in m_space + */ + [[nodiscard]] ::mfem::FiniteElementSpace& space() const { + using namespace utils::errors; + serif_ensure(m_space != nullptr, "Global scalars do not have a finite-element space."); + return *m_space; + } + + /** + * @brief Return the borrowed mesh associated with this DOF map. + * @return A reference to the mesh in m_mesh + */ + [[nodiscard]] ::mfem::Mesh& mesh() const noexcept { return m_mesh; } + + /** + * @brief Return the number of full-space DOFs, including those not supported by the element support. + * @return The number of full-space DOFs (not reduced by element support) + */ + [[nodiscard]] int full_size() const noexcept { return m_full_to_reduced.Size(); } + + /** + * @brief Return the number of supported DOFs, which is the size of the reduced coefficient vector. + * @return The number of supported DOFs (reduced by element support) + */ + [[nodiscard]] int size() const noexcept { + return m_reduced_to_full.Size(); + } + + /** + * @brief Located a full space dof in the reduced space. E.g if I ask where is dof at index x in the full space this will tell me what index it is in the reduced space. If the dof is not supported this will return std::nullopt. + * @param full_dof index of the dof in the full space + * @return An optional index of the dof in the reduced space, or std::nullopt if the dof is not supported. + */ + [[nodiscard]] std::optional reduced_dof(const int full_dof) const { + using namespace utils::errors; + serif_ensure(full_dof >= 0 && full_dof < full_size(), "Full DOF index is out of range."); + const int reduced = m_full_to_reduced[full_dof]; + return reduced < 0 ? std::nullopt : std::optional{reduced}; + } + + /** + * @brief Check whether an element is supported by this dof map + * @param element element index + * @return true if element is supported by this dof map, false otherwise + */ + [[nodiscard]] bool supports(const int element) const { + using namespace utils::errors; + + serif_ensure(element >= 0 && element < m_element_support.Size(), "Element index is out of range."); + return m_element_support[element] != 0; + } + + /** + * @brief Validate the dof map + */ + void validate() const { + using namespace utils::errors; + + serif_ensure(m_mesh.GetSequence() == m_mesh_sequence && m_mesh.GetNE() == m_attributes.Size() && (m_space == nullptr || (m_space->GetMesh() == &m_mesh && m_space->GetSequence() == m_space_sequence && m_space->GetVSize() == full_size())), "Rebuild term bindings after discretization updates."); + + for (int element = 0; element < m_mesh.GetNE(); ++element) { + serif_ensure(m_mesh.GetAttribute(element) == m_attributes[element],"Term binding domain attributes changed unexpectedly."); + } + } + + /** + * @brief Check supported coefficient count and finiteness. + * @param values vector of coefficients to validate + */ + void validate_values(const ::mfem::Vector& values) const { + using namespace utils::errors; + + validate(); + serif_ensure(values.Size() == size(), "Term binding size is {}, expected {} supported DOFs.", values.Size(), size()); + + for (int i = 0; i < values.Size(); ++i) { + serif_ensure(std::isfinite(values[i]), "Term binding contains a non-finite coefficient."); + } + } + + /** + * @brief Restrict a full coefficient vector to supported DOFs. + * @param full full coefficient vector + * @return vector of supported coefficients + */ + [[nodiscard]] ::mfem::Vector gather(const ::mfem::Vector& full) const { + using namespace utils::errors; + + validate(); + serif_ensure(full.Size() == full_size(), "Full coefficient vector does not match its space."); + + ::mfem::Vector result(size()); + for (int i = 0; i < size(); ++i) { + result[i] = full[m_reduced_to_full[i]]; + } + + return result; + } + + /** + * @brief Expand supported coefficients with zero outside their support. + * @param reduced vector of supported coefficients + * @return full coefficient vector with zeros outside the support + */ + [[nodiscard]] ::mfem::Vector scatter(const ::mfem::Vector& reduced) const { + using namespace utils::errors; + + validate_values(reduced); + + ::mfem::Vector result(full_size()); + result = 0.0; + for (int i = 0; i < size(); ++i) { + result[m_reduced_to_full[i]] = reduced[i]; + } + + return result; + } + + /** + * @brief Extract element coefficients, including primal orientation transforms. + * @param element element index + * @param reduced vector of supported coefficients + * @param values vector to store extracted element coefficients + */ + void element_values(const int element, const ::mfem::Vector& reduced, ::mfem::Vector& values) const { + using namespace utils::errors; + + validate(); + + (void)supports(element); + + serif_ensure(reduced.Size() == size(), "Element extraction requires a supported coefficient vector."); + + if (m_space == nullptr) { + values.SetSize(1); + values[0] = supports(element) ? reduced[0] : 0.0; + return; + } + + ::mfem::Array dofs; + const auto* transformation = m_space->GetElementVDofs(element, dofs); + values.SetSize(dofs.Size()); + + for (int i = 0; i < dofs.Size(); ++i) { + const int full = dofs[i] < 0 ? -1 - dofs[i] : dofs[i]; + const int index = m_full_to_reduced[full]; + values[i] = index < 0 ? 0.0 : (dofs[i] < 0 ? -reduced[index] : reduced[index]); + } + + if (transformation != nullptr) { + transformation->InvTransformPrimal(values); + } + } + + /** + * @brief Apply dual orientation transforms and add an element residual. + * @param element element index + * @param values vector of element coefficients + * @param full full coefficient vector to which the element residual will be added + */ + void add_element(const int element, ::mfem::Vector& values, ::mfem::Vector& full) const { + using namespace utils::errors; + + validate(); + if (m_space == nullptr) { + serif_ensure(values.Size() == 1 && full.Size() == 1, "Global scalar assembly requires one value."); + if (supports(element)) { full[0] += values[0]; } + return; + } + + // Cast to void to noop that this value wont be used + (void)supports(element); + + ::mfem::Array dofs; + const auto* transformation = m_space->GetElementVDofs(element, dofs); + + serif_ensure(values.Size() == dofs.Size() && full.Size() == full_size(),"Element assembly vector sizes do not match the space."); + + if (transformation != nullptr) { + transformation->TransformDual(values); + } + full.AddElementVector(dofs, values); + } + + private: + ::mfem::FiniteElementSpace* m_space; ///< Borrowed pointer to the finite-element space associated with this dof map, or nullptr for a global scalar. This must outlive this class + ::mfem::Mesh& m_mesh; ///< Borrowed reference to the mesh associated with this dof map. This must outlive this class + long m_space_sequence; ///< Sequence number of the finite-element space at construction, or -1 for a global scalar + long m_mesh_sequence; ///< Sequence number of the mesh at construction + ::mfem::Array m_attributes; ///< Array of attribute indices associated with this dof map + ::mfem::Array m_element_support; ///< Array of 0/1 values indicating whether each element is supported by this dof map + ::mfem::Array m_reduced_to_full; ///< Mapping from reduced to full degrees of freedom + ::mfem::Array m_full_to_reduced; ///< Mapping from full to reduced degrees of freedom + }; +} // namespace serif::discretization::backend::mfem diff --git a/src/include/serif/discretization/backend/mfem/dof_selection.hpp b/src/include/serif/discretization/backend/mfem/dof_selection.hpp new file mode 100644 index 00000000..8e94db07 --- /dev/null +++ b/src/include/serif/discretization/backend/mfem/dof_selection.hpp @@ -0,0 +1,155 @@ +#pragma once + +/** + * @file dof_selection.hpp + * @author Emily M. Boudreaux + * @brief This file contains the @ref DofSelection which allows for the selection of a subset of DOFs from a term binding, with boundary attribute information. + */ + +#include "serif/discretization/backend/mfem/term_binding.hpp" +#include "serif/discretization/backend/mfem/dof_map.hpp" + +namespace serif::discretization::backend::mfem { + + /** + * @brief A selection of supported DOFs from a term binding, with boundary attribute information. + * The purpose here is to allow for the selection of a subset of DOFs from a term binding, + * such as those associated with a particular boundary, while still retaining the ability to + * validate the selection against the underlying mesh and term binding. + */ + class DofSelection final { + public: + /** + * @brief Construct a DOF selection from a term binding and a list of reduced DOFs. + * @param binding A reference to the term binding from which the DOFs are selected. Further, this binding must outlive the DOF selection. + * @param reduced_dofs A list of reduced DOFs to select. This list should be in the reduced space of the term binding. Further, these must be sorted in + * a strictly increasing order and must be unique. This is enforced at construction time. + */ + DofSelection( + const DofMap& binding, + const ::mfem::Array& reduced_dofs, + const bool track_coordinates = false + ) : + m_binding(binding), + m_dofs(reduced_dofs), + m_boundary_attributes(binding.mesh().GetNBE()), + m_nodes_sequence(track_coordinates ? binding.mesh().GetNodesSequence() : -1) { + using namespace utils::errors; + + binding.validate(); + int previous = -1; + + for (const int dof : m_dofs) { + serif_ensure(dof > previous && dof < binding.size(), "Selected DOFs must be sorted, unique, and inside the term binding."); + previous = dof; + } + + for (int boundary = 0; boundary < m_boundary_attributes.Size(); ++boundary) { + m_boundary_attributes[boundary] = binding.mesh().GetBdrAttribute(boundary); + } + } + + /** + * @brief Return the underlying term binding associated with this DOF selection. + * @return A reference to the underlying term binding. + */ + [[nodiscard]] const DofMap& binding() const noexcept { return m_binding; } + + /** + * @brief Return the list of reduced DOFs in this selection. + * @return A reference to the list of reduced DOFs. Note these are in the reduced space of the binding provided at construction time. + */ + [[nodiscard]] const ::mfem::Array& reduced_dofs() const noexcept { return m_dofs; } + + /** + * @brief Return the number of selected DOFs in this selection. + * @return The number of selected DOFs. + */ + [[nodiscard]] int size() const noexcept { return m_dofs.Size(); } + + /** + * @brief Return the list of boundary attributes associated with this selection. + * @return A reference to the list of boundary attributes. These are in the same order as the boundaries in the underlying mesh. + */ + [[nodiscard]] const ::mfem::Array& boundary_attributes() const noexcept { return m_boundary_attributes; } + + /** + * @brief Report whether this selection becomes stale when the fixed mesh coordinates change. + * Point constraints require this tracking even when a caller supplies the selected DOFs explicitly. + * @return True if coordinate updates invalidate the selection. + */ + [[nodiscard]] bool tracks_coordinates() const noexcept { return m_nodes_sequence >= 0; } + + /** + * @brief Validate the DOF selection against the underlying term binding and mesh. + * This checks that the underlying term binding and mesh have not changed since the DOF selection was constructed. + * If they have changed, an exception is thrown. This is important to ensure that the DOF selection remains valid + * and consistent with the underlying data structures. + */ + void validate() const { + using namespace utils::errors; + + m_binding.validate(); + const auto& mesh = m_binding.mesh(); + serif_ensure(m_nodes_sequence < 0 || mesh.GetNodesSequence() == m_nodes_sequence, "Rebuild point selections after fixed-mesh coordinate changes."); + + serif_ensure(mesh.GetNBE() == m_boundary_attributes.Size(), "Rebuild DOF selections after boundary changes."); + + for (int boundary = 0; boundary < m_boundary_attributes.Size(); ++boundary) { + serif_ensure(mesh.GetBdrAttribute(boundary) == m_boundary_attributes[boundary], "Rebuild DOF selections after boundary attributes change."); + } + } + + private: + const DofMap& m_binding; ///< Reference to the underlying term binding. This must outlive the DOF selection. + ::mfem::Array m_dofs; ///< List of reduced DOFs in this selection. These are in the reduced space of the binding provided at construction time. + ::mfem::Array m_boundary_attributes; ///< List of boundary attributes associated with this selection. These are in the same order as the boundaries in the underlying mesh. + long m_nodes_sequence; ///< Coordinate revision, or -1 for selections independent of geometry. + }; + + /** + * @brief Select DOFs associated with a particular boundary from a term binding. + * @tparam BoundaryT The boundary type to select. This must be a valid boundary in the domain schema of the term binding. + * @tparam TermT The term type of the term binding. This must be a valid term type in the domain schema of the term binding. (e.g. something like Density) + * @tparam SchemaT The domain schema type of the term binding. This must be a valid domain schema type that contains the specified boundary. + * @param binding The term binding from which to select DOFs. This must outlive the returned DOF selection. + * @return A DofSelection object representing the selected DOFs associated with the specified boundary. + */ + template + [[nodiscard]] DofSelection select_boundary_dofs( + const TermBinding& binding + ) { + using namespace utils::errors; + + static_assert(SchemaT::template contains_boundary(), "The boundary must be registered in the domain schema."); + binding.validate(); + + ::mfem::Array marker(binding.size()); + marker = 0; + + ::mfem::Array dofs; + ::mfem::Array selected; + + auto& space = binding.space(); + const auto& mesh = binding.mesh(); + + for (int boundary = 0; boundary < mesh.GetNBE(); ++boundary) { + if (!SchemaT::template boundary_id_matches(mesh.GetBdrAttribute(boundary))) { + continue; + } + + space.GetBdrElementVDofs(boundary, dofs); + + for (const int encoded : dofs) { + const auto reduced = binding.reduced_dof(encoded < 0 ? -1 - encoded : encoded); + serif_ensure(reduced.has_value(), "A selected boundary DOF lies outside the term's support."); + marker[*reduced] = 1; + } + } + + ::mfem::FiniteElementSpace::MarkerToList(marker, selected); + serif_ensure(selected.Size() > 0, "The requested boundary has no DOFs on this serial mesh."); + + return DofSelection(binding, selected); + } +} // namespace serif::discretization::backend::mfem diff --git a/src/include/serif/discretization/backend/mfem/mfem_backend.hpp b/src/include/serif/discretization/backend/mfem/mfem_backend.hpp index e69de29b..52fa61c8 100644 --- a/src/include/serif/discretization/backend/mfem/mfem_backend.hpp +++ b/src/include/serif/discretization/backend/mfem/mfem_backend.hpp @@ -0,0 +1,20 @@ +#pragma once + +/** + * @file mfem_backend.hpp + * @author Emily M. Boudreaux + * @brief This is a utility file which allows one import to bring in all of the mfem discretization backend functionality. This file includes + * - dof_map.hpp + * - dof_selection.hpp + * - point_dofs.hpp + * - tensor_basis.hpp + * - term_binding.hpp + */ + +#include "serif/discretization/backend/mfem/basis_evaluation.hpp" +#include "serif/discretization/backend/mfem/coefficient_tracking.hpp" +#include "serif/discretization/backend/mfem/dof_map.hpp" +#include "serif/discretization/backend/mfem/dof_selection.hpp" +#include "serif/discretization/backend/mfem/point_dofs.hpp" +#include "serif/discretization/backend/mfem/tensor_basis.hpp" +#include "serif/discretization/backend/mfem/term_binding.hpp" \ No newline at end of file diff --git a/src/include/serif/discretization/backend/mfem/point_dofs.hpp b/src/include/serif/discretization/backend/mfem/point_dofs.hpp new file mode 100644 index 00000000..83e9ebf7 --- /dev/null +++ b/src/include/serif/discretization/backend/mfem/point_dofs.hpp @@ -0,0 +1,82 @@ +#pragma once + +/** + * @file point_dofs.hpp + * @author Emily M. Boudreaux + * @brief This file contains the implementation of function select_point_dofs. This function is used to select the components + * from a single H1 vertex in the fixed computational mesh. This is used to implement point constraints in the discretization (e.g. fixed center of mass) + */ + +#include +#include + +#include "serif/discretization/backend/mfem/dof_selection.hpp" + +namespace serif::discretization::backend::mfem { + /** + * @brief Select the dofs associated with a single point in the reference domain. This is used to implement point constraints in the discretization (e.g. fixed center of mass). Note that this function only works for H1 terms and will throw an error if used on a non-H1 term. + * @tparam TermT The term type to select the dofs for. This must be a spatial term with an H1 functional space. + * @tparam SchemaT The schema type associated with the term. + * @param binding The term binding to select the dofs from. + * @param position The position of the point in the reference domain. + * @param tolerance The tolerance for matching the point coordinates. + * @return A DofSelection object representing the selected dofs. + */ + template + requires (blocks::SpatialTermTraits::family == functional_space::Family::h1) + [[nodiscard]] DofSelection select_point_dofs( + const TermBinding& binding, + const std::array& position = {}, + const double tolerance = 1.0e-12 + ) { + using namespace utils::errors; + + // Mostly the start of this is just validation down to... + binding.validate(); + serif_ensure(std::isfinite(tolerance) && tolerance >= 0.0, "Point selection requires a finite nonnegative coordinate tolerance."); + for (const double coordinate : position) { + serif_ensure(std::isfinite(coordinate), "Point selection requires finite reference coordinates."); + } + + auto& space = binding.space(); + for (int element = 0; element < space.GetNE(); ++element) { + serif_ensure(dynamic_cast(space.GetFE(element)) != nullptr, "Point values require a nodal H1 basis; positive/Bernstein coefficients are not nodal values."); + } + // down to here. + + // Then we come through and we try to find the vertex that matches the requested position within the numerical tolerance. e.g. we might pass + // a position of (0, 0, 0) here and a tolerance of 1e-12 and this will then match any vertex within a euclidean distance of 1e-12 from (0, 0, 0) + // There are better ways to do this but this works for now. + int selected_vertex = -1; + for (int vertex = 0; vertex < binding.mesh().GetNV(); ++vertex) { + std::array coordinates; + binding.mesh().GetNode(vertex, coordinates.data()); + bool matches = true; + + for (int axis = 0; axis < 3; ++axis) { + matches = matches && std::abs(coordinates[axis] - position[axis]) <= tolerance; + } + if (matches) { + serif_ensure(selected_vertex < 0, "The requested point matches multiple mesh vertices; use a unique material point or a smaller tolerance."); + selected_vertex = vertex; + } + } + + // Note that here we enforce that the selected position is within tolerance of a mesh vertex. A future implementation might consider some more + // complex scheme that involves finding the bounding elements and enforcing some sort of interpolation scheme. + serif_ensure(selected_vertex >= 0, "The requested point is not a mesh vertex; point constraints currently require a unique nodal H1 vertex."); + + // Then we loop through and select the dofs associated with the selected vertex in the reduced dof space of the provided term binding. + // If that position is outside the support of the term binding this will throw. + ::mfem::Array vertex_dofs, selected; + space.GetVertexDofs(selected_vertex, vertex_dofs); + serif_ensure(vertex_dofs.Size() == 1, "A point constraint requires one scalar nodal DOF at its vertex."); + for (int component = 0; component < space.GetVDim(); ++component) { + const auto reduced = binding.reduced_dof(space.DofToVDof(vertex_dofs[0], component)); + serif_ensure(reduced.has_value(), "The selected point lies outside the term's support."); + selected.Append(*reduced); + } + + return DofSelection(binding, selected, true); + } +} // namespace serif::discretization::backend::mfem diff --git a/src/include/serif/discretization/backend/mfem/tensor_basis.hpp b/src/include/serif/discretization/backend/mfem/tensor_basis.hpp new file mode 100644 index 00000000..aa10529c --- /dev/null +++ b/src/include/serif/discretization/backend/mfem/tensor_basis.hpp @@ -0,0 +1,489 @@ +#pragma once + +/** + * @file tensor_basis.hpp + * @author Emily M. Boudreaux + * @brief Implementation of fast tensor evaluation functions. These do not contribute physics but are used to accelerate evaluations. + */ + +#include +#include +#include +#include +#include +#include +#include + +#include "mfem/fem/fe/fe_h1.hpp" +#include "mfem/fem/fe/fe_l2.hpp" +#include "mfem/fem/fe/fe_rt.hpp" +#include "serif/utils/error/errors.hpp" + +namespace serif::discretization::backend::mfem { + /** + * @enum ReferenceOperation + * @brief Options enum to control what is actually evaluated at the reference points. + */ + enum class ReferenceOperation { + value, ///< Evaluate the basis functions at the reference points + gradient, ///< Evaluate the gradients of the basis functions at the reference points + value_and_gradient, ///< Evaluate both the basis functions and their gradients at the reference points + divergence ///< Evaluate the divergence of the basis functions at the reference points + }; + + /** + * @enum ContractionMode + * @brief Options enum to control how the tensor contractions are performed. + */ + enum class ContractionMode { + automatic, ///< Automatically choose the contraction mode based on the number of points and dofs + dense ///< Use a dense contraction mode + }; + + /** + * @struct TensorWorkspace + * @brief A workspace for tensor contractions to avoid repeated allocations and pack memory (TODO: does this actual pack memory in a useful manner when using mfem vectors? Probably something to look into). + */ + struct TensorWorkspace { + ::mfem::Vector coefficients; ///< The coefficients of the basis functions at the reference points + ::mfem::Vector first; ///< The first contraction of the coefficients with the basis functions at the reference points + ::mfem::Vector second; ///< The second contraction of the coefficients with the basis functions at the reference points + }; + + /** + * @class TensorBasis + * @brief A class that represents a tensor basis for a finite element. This class is used to evaluate the basis functions and their derivatives at a set of reference + * points. + */ + class TensorBasis final { + /** + * @struct Product + * @brief A struct that represents a product of one-dimensional basis functions. + */ + struct Product { + std::array<::mfem::DenseMatrix, 3> factors; ///< The one-dimensional basis functions for each axis + ::mfem::Array dofs; ///< The dofs associated with the product + int output; ///< The output index for the product (used for mapping to the correct component of the output vector) + }; + + public: + /// Return no tensor plan unless its factorization agrees with MFEM's canonical basis. + + /** + * @brief Make a tensor basis for a given finite element, integration rule, reference operation, and number of copies. This function returns a shared pointer to a TensorBasis object that can be used to do efficient evaluation of the basis functions and their derivatives at reference points. + * @param element The finite element for which to construct the tensor basis. + * @param rule The integration rule to use for evaluating the basis functions. + * @param operation The reference operation to apply (e.g., value, divergence). + * @param copies The number of copies of the basis functions. + * @returns A shared pointer to a TensorBasis object if the tensor basis can be constructed, or an empty shared pointer if it cannot. + */ + [[nodiscard]] static std::shared_ptr make( + const ::mfem::FiniteElement& element, + const ::mfem::IntegrationRule& rule, + const ReferenceOperation operation, + const int copies + ) { + const auto* scalar = dynamic_cast(&element); + const auto* rt = dynamic_cast(&element); + const bool scalar_hex = dynamic_cast(&element) || dynamic_cast(&element); + + // TODO: This is an awful lot of checks, compress these into some sort of "is this a valid tensor basis" function that can be used in the tests as well. + if (copies <= 0 || (!scalar_hex && rt == nullptr) || rule.GetNPoints() == 0 || (scalar_hex && (scalar == nullptr || element.GetMapType() != ::mfem::FiniteElement::VALUE || operation == ReferenceOperation::divergence)) || (rt && (copies != 1 || (operation != ReferenceOperation::value && operation != ReferenceOperation::divergence)))) { + return {}; + } + + const int q = static_cast(std::lround(std::cbrt(rule.GetNPoints()))); + if (q * q * q != rule.GetNPoints()) { + return {}; + } + + for (int z = 0; z < q; ++z) { + for (int y = 0; y < q; ++y) { + for (int x = 0; x < q; ++x) { + const auto& point = rule.IntPoint(x + q * (y + q * z)); + if (point.x != rule.IntPoint(x).x || point.y != rule.IntPoint(q * y).y || point.z != rule.IntPoint(q * q * z).z) { + return {}; + } + } + } + } + + auto result = std::make_shared(); + + result->m_points = rule.GetNPoints(); + result->m_dofs = element.GetDof() * copies; + + const int samples = operation == ReferenceOperation::gradient ? 3 : operation == ReferenceOperation::value_and_gradient ? 4 : 1; + + result->m_components = rt ? (operation == ReferenceOperation::value ? 3 : 1) : copies * samples; + + const int closed = element.GetOrder() + 1; + const int block = rt ? element.GetDof() / 3 : element.GetDof(); + const auto& native_map = rt ? rt->GetDofMap() : scalar->GetDofMap(); + + if (native_map.Size() != 0 && native_map.Size() != element.GetDof()) { + return {}; + } + + // TODO: I once again make use of a vector since it is the most obvious conceptual datatype. However, the standard backed implementation of + // vector is not actually a vector of bools; rather, it is a bitset. I dont think this matters for what I am doing but is is probably + // worth being aware of this + std::vector seen(element.GetDof(), false); + + for (int dof = 0; dof < element.GetDof(); ++dof) { + const int encoded = native_map.Size() ? native_map[dof] : dof; + const int native = encoded < 0 ? -1 - encoded : encoded; + + if ((!rt && encoded < 0) || native < 0 || native >= element.GetDof() || seen[native]) { + return {}; + } + + seen[native] = true; + } + + for (int component = 0; component < (rt ? 3 : copies); ++component) { + for (int sample = 0; sample < samples; ++sample) { + auto&[factors, dofs, output] = result->m_products.emplace_back(); + + output = rt ? (operation == ReferenceOperation::divergence ? 0 : component) : component * samples + sample; + const int derivative = rt ? (operation == ReferenceOperation::divergence ? component : -1) : operation == ReferenceOperation::gradient ? sample : operation == ReferenceOperation::value_and_gradient ? sample - 1 : -1; + + int tensor_dofs = 1; + for (int axis = 0; axis < 3; ++axis) { + const int count = rt && axis != component ? closed - 1 : closed; + tensor_dofs *= count; + + ::mfem::DenseMatrix& factor = factors[axis]; + factor.SetSize(q, count); + + ::mfem::Vector values(count); + ::mfem::Vector gradients(count); + + const auto& basis = rt ? (axis == component ? rt->GetBasis1D() : rt->GetOpenBasis1D()) : scalar->GetBasis1D(); + for (int point = 0; point < q; ++point) { + const double coordinate = axis == 0 ? rule.IntPoint(point).x : axis == 1 ? rule.IntPoint(q * point).y : rule.IntPoint(q * q * point).z; + + // Note that here this is just two overloads of Eval that mfem::Poly_1D provides, the first one evaluates the gradient and the latter evaluates the value. + if (axis == derivative) { + basis.Eval(coordinate, values, gradients); + } + else { + basis.Eval(coordinate, values); + } + + for (int dof = 0; dof < count; ++dof) { + factor(point, dof) = axis == derivative ? gradients[dof] : values[dof]; + } + } + } + + if (tensor_dofs != block) { + return {}; + } + + dofs.SetSize(block); + for (int lex = 0; lex < block; ++lex) { + const int index = rt ? component * block + lex : lex; + const int encoded = native_map.Size() ? native_map[index] : index; + + dofs[lex] = rt ? encoded : encoded + component * element.GetDof(); + } + } + } + + return result->agrees_with(element, rule, operation) ? result : nullptr; + } + + /** + * @brief Get the number of quadrature points in the Cartesian rule. + * @returns The number of quadrature points in the Cartesian rule. + */ + [[nodiscard]] int points() const noexcept { + return m_points; + } + + /** + * @brief Get the number of degrees of freedom in the tensor basis. + * @returns The number of degrees of freedom in the tensor basis. + */ + [[nodiscard]] int dofs() const noexcept { + return m_dofs; + } + + /** + * @brief Get the number of canonical output components per point in the tensor basis. + * @returns The number of canonical output components per point in the tensor basis. + */ + [[nodiscard]] int components() const noexcept { + return m_components; + } + + /// Evaluate canonical samples, or their coordinate transpose, with caller-owned scratch. + + /** + * @brief Apply the tensor basis to an input vector and store the results in an output vector. This function performs the tensor contraction of the input + * vector with the tensor basis. + * @param input The input vector to be contracted with the tensor basis. + * @param output The output vector to store the result of the contraction. (Owned by the caller) + * @param workspace The workspace to use for the contraction. (Owned by the caller) + * @param transpose Whether to transpose the contraction (default is false). + * + * @pre The input vector must have a size equal to the number of degrees of freedom in the tensor basis if transpose is false, or equal to the number of points times the number of components if transpose is true. + * @pre The output vector must be constructed + * @post The output vector will have a size equal to the number of points times the number of components if transpose is false, or equal to the number of degrees of freedom in the tensor basis if transpose is true. + * @post The output vector will contain the result of the tensor contraction of the input vector with the tensor basis. + */ + void apply(const ::mfem::Vector& input, ::mfem::Vector& output, TensorWorkspace& workspace, const bool transpose = false) const { + using namespace utils::errors; + + serif_ensure(input.Size() == (transpose ? m_points * m_components : m_dofs), "Tensor contraction input has the wrong coefficient or quadrature size."); + + output.SetSize(transpose ? m_dofs : m_points * m_components); + output = 0.0; + + for (const auto& product : m_products) { + auto& coefficients = workspace.coefficients; + coefficients.SetSize(transpose ? m_points : product.dofs.Size()); + + for (int i = 0; i < coefficients.Size(); ++i) { + const int encoded = transpose ? 0 : product.dofs[i]; + coefficients[i] = transpose ? input[i * m_components + product.output] : encoded < 0 ? -input[-1 - encoded] : input[encoded]; + } + + std::array shape; + for (int axis = 0; axis < 3; ++axis) { + shape[axis] = transpose ? product.factors[axis].Height() : product.factors[axis].Width(); + } + + const ::mfem::Vector* source = &coefficients; + for (int step = 0; step < 3; ++step) { + const int axis = transpose ? 2 - step : step; + auto& destination = step % 2 == 0 ? workspace.first : workspace.second; + + contract_axis(product.factors[axis], axis, shape, *source, destination, transpose); + + shape[axis] = transpose ? product.factors[axis].Width() : product.factors[axis].Height(); + source = &destination; + } + + for (int i = 0; i < source->Size(); ++i) { + if (transpose) { + const int encoded = product.dofs[i]; + output[encoded < 0 ? -1 - encoded : encoded] += (encoded < 0 ? -1.0 : 1.0) * (*source)[i]; + } else { + output[i * m_components + product.output] += (*source)[i]; + } + } + } + } + + private: + /** + * @breif apply one factor along one tensor axis, teh other two indices will remain unchanged + * @param factor The one-dimensional basis function to apply along the specified axis. + * @param axis The axis along which to apply the factor (0, 1, or 2). + * @param shape The shape of the input tensor (number of points along each axis). + * @param input The input vector to be contracted with the factor. + * @param output The output vector to store the result of the contraction. + * @param transpose Whether to transpose the factor matrix before applying it. + * + * @pre The input vector must have a size equal to the product of the dimensions specified in the shape array. + * @post The output vector will have a size equal to the product of the dimensions specified in the next array. + * @post The output vector will contain the result of the contraction of the input vector with the factor along the specified axis. + */ + static void contract_axis( + const ::mfem::DenseMatrix& factor, + const int axis, + const std::array& shape, + const ::mfem::Vector& input, + ::mfem::Vector& output, + const bool transpose + ) { + auto next = shape; + next[axis] = transpose ? factor.Width() : factor.Height(); + output.SetSize(next[0] * next[1] * next[2]); + for (int z = 0; z < next[2]; ++z) { + for (int y = 0; y < next[1]; ++y) { + for (int x = 0; x < next[0]; ++x) { + std::array index{x, y, z}; + + const int row = index[axis]; + double value = 0.0; + + for (int column = 0; column < shape[axis]; ++column) { + index[axis] = column; + value += (transpose ? factor(column, row) : factor(row, column)) * input[index[0] + shape[0] * (index[1] + shape[1] * index[2])]; + } + + output[x + next[0] * (y + next[1] * z)] = value; + } + } + } + } + + /// Verify signed ordering and every requested basis entry before enabling the fast path. + + /** + * @brief Verify that the signed ordering and every requested basis entry agrees with MFEM's canonical basis. + * @param element The finite element to verify against. + * @param rule The integration rule to use for evaluating the basis functions. + * @param operation The reference operation to apply (e.g., value, divergence). + * + * @returns True if the tensor basis agrees with MFEM's canonical basis, false otherwise. + */ + [[nodiscard]] bool agrees_with( + const ::mfem::FiniteElement& element, + const ::mfem::IntegrationRule& rule, + const ReferenceOperation operation + ) const { + const bool vector = element.GetRangeType() == ::mfem::FiniteElement::VECTOR; + const int samples = operation == ReferenceOperation::gradient ? 3 : operation == ReferenceOperation::value_and_gradient ? 4 : 1; + + ::mfem::Vector values(element.GetDof()); + ::mfem::DenseMatrix derivatives(element.GetDof(), 3); + + for (int point = 0; point < m_points; ++point) { + if (vector) { + if (operation == ReferenceOperation::divergence) { + element.CalcDivShape(rule.IntPoint(point), values); + } + else { + element.CalcVShape(rule.IntPoint(point), derivatives); + } + } else { + element.CalcShape(rule.IntPoint(point), values); + if (samples > 1) { + element.CalcDShape(rule.IntPoint(point), derivatives); + } + } + + for (const auto&[factors, dofs, output] : m_products) { + const int q = factors[0].Height(), nx = factors[0].Width(), ny = factors[1].Width(); + + for (int dofID = 0; dofID < dofs.Size(); ++dofID) { + const int encoded = dofs[dofID], native = (encoded < 0 ? -1 - encoded : encoded) % element.GetDof(); + + double expected{0.0}; + if (vector) { + expected = operation == ReferenceOperation::divergence ? values[native] : derivatives(native, output); + } + else { + const int derivative = operation == ReferenceOperation::gradient ? output % samples : output % samples - 1; + expected = samples == 1 || derivative < 0 ? values[native] : derivatives(native, derivative); + } + + if (vector && operation == ReferenceOperation::value) { + for (int axis = 0; axis < 3; ++axis) { + if (axis != output && derivatives(native, axis) != 0.0) { + return false; + } + } + } + + expected *= encoded < 0 ? -1.0 : 1.0; + const double actual = factors[0](point % q, dofID % nx) * factors[1]((point / q) % q, (dofID / nx) % ny) * factors[2](point / (q * q), dofID / (nx * ny)); + + // Similar to above in this file, this is a long series of checks, maybe refactor into a function? + if (!std::isfinite(actual) || !std::isfinite(expected) || std::abs(actual - expected) > 64.0 * std::numeric_limits::epsilon() * std::max(1.0, std::abs(expected))) { + return false; + } + } + } + } + + return true; + } + + int m_points{0}; ///< The number of quadrature points in the Cartesian rule. + int m_dofs{0}; ///< The number of degrees of freedom in the tensor basis. + int m_components{0}; ///< The number of canonical output components per point in the tensor basis. + std::vector m_products; ///< The products of one-dimensional basis functions for each axis. (TOOD: consider if this should be a heap allocated vector of it this could be put on the stack somehow) + }; + + /// Share immutable tensor plans within one operator; copied rules are keyed by their contents. + + /** + * @class TensorBasisCache + * @brief A cache for TensorBasis objects to avoid repeated construction of the same tensor basis. + * This class stores a vector of "Entry" structs which each contain a non-owning pointer to some element, + * options, a shared pointer to a TensorBasis, and the points used to construct the TensorBasis. + */ + class TensorBasisCache final { + + /** + * @struct Entry + * @brief A struct that represents an entry in the TensorBasisCache. Each entry contains a non-owning pointer to a finite element, the reference operation, the number of copies, the points used to construct the TensorBasis, and a shared pointer to the TensorBasis itself. + * + * @todo Consider if we can remove the shared pointer from here, shared pointers introduce overhead with reference counting and that may slow down the cache. + */ + struct Entry { + const ::mfem::FiniteElement* element; + ReferenceOperation operation; + int copies; + std::vector> points; + std::shared_ptr basis; + }; + public: + /// Retrieve a verified plan, caching unsupported combinations as well. + + /** + * @brief Retrieve a verified tensor basis, caching unsupported combinations as well. + * @param element The finite element for which to retrieve the tensor basis. + * @param rule The integration rule to use for evaluating the basis functions. + * @param operation The reference operation to apply (e.g., value, divergence). + * @param copies The number of copies of the basis functions. + * @returns A shared pointer to a TensorBasis object if the tensor basis can be constructed, or an empty shared pointer if it cannot. + * + * @todo I think it would make sense, at some point, to switch this to an std::optional return + */ + [[nodiscard]] std::shared_ptr get( + const ::mfem::FiniteElement& element, + const ::mfem::IntegrationRule& rule, + ReferenceOperation operation, + const int copies + ) { + std::vector> points; + + for (int q = 0; q < rule.GetNPoints(); ++q) { + const auto& point = rule.IntPoint(q); + points.push_back({point.x, point.y, point.z, point.weight}); + } + + for (const auto& [el, op, cps, pts, bs] : m_entries) { + if (el == &element && op == operation && cps == copies && pts == points) { + return bs; + } + } + + const auto basis = TensorBasis::make(element, rule, operation, copies); + + m_entries.push_back({ + .element = &element, + .operation = operation, + .copies = copies, + .points = std::move(points), + .basis = basis + }); + + return basis; + } + + /** + * @brief Get the number of entries in the TensorBasisCache. + * @return The number of entries in the cache. + */ + [[nodiscard]] std::size_t size() const noexcept { + return m_entries.size(); + } + + /** + * @brief Clear the TensorBasisCache, removing all entries. + */ + void clear() { + m_entries.clear(); + } + + private: + std::vector m_entries; ///< The entries in the TensorBasisCache, each containing a non-owning pointer to a finite element, the reference operation, the number of copies, the points used to construct the TensorBasis, and a shared pointer to the TensorBasis itself. + }; +} // namespace serif::discretization::backend::mfem diff --git a/src/include/serif/discretization/backend/mfem/term_binding.hpp b/src/include/serif/discretization/backend/mfem/term_binding.hpp new file mode 100644 index 00000000..798197dc --- /dev/null +++ b/src/include/serif/discretization/backend/mfem/term_binding.hpp @@ -0,0 +1,115 @@ +#pragma once +/** + * @file term_binding.hpp + * @author Emily M. Boudreaux + * @brief A runtime binding between some term descriptor and its MFEM finite-element space. This is used to project + * down the compile time term descriptors into runtime reinterpretations which can then be used to extract out things like + * the list of dofs, the element support, and the mesh. + */ + +#include + +#include "serif/discretization/blocks/traits.hpp" +#include "serif/discretization/domain/schema/concepts.hpp" +#include "serif/discretization/backend/mfem/dof_map.hpp" +#include "serif/utils/error/errors.hpp" + +namespace serif::discretization::backend::mfem { + /** + * @brief A runtime binding between a term descriptor and its MFEM finite-element space. This is used to project + * down the compile time term descriptors into runtime reinterpretations which can then be used to extract + * things like the list of dofs, the element support, and the mesh. Note that this is a child of + * DofMap, therefor all term bindings are also dof maps + * @tparam TermT The term descriptor type to bind to. Note that if the FiniteElementSpace is not a global scalar then this must be a spatial term descriptor (@ref SpatialTerm) + * @tparam SchemaT The domain schema type to bind to. This must be a valid schema (@ref IsSchema) + */ + template + requires (blocks::SpatialTerm || blocks::GlobalScalarTerm) + class TermBinding final : public DofMap { + public: + using Term = TermT; + using Schema = SchemaT; + + /** + * @brief Construct a term binding from an MFEM finite-element space. This requires that the term be a valid spatial term (@ref SpatialTerm) + * @param space The finite-element space to bind to. + */ + explicit TermBinding( + ::mfem::FiniteElementSpace& space + ) requires blocks::SpatialTerm : + DofMap(space, element_support(checked_mesh(space))) { + using Traits = blocks::SpatialTermTraits; + using namespace utils::errors; + + static_assert(Traits::family == functional_space::Family::h1 || Traits::family == functional_space::Family::l2 || Traits::family == functional_space::Family::hdiv, "The MFEM backend currently only supports H1, L2, and H(div) spaces."); + static_assert(Traits::rank <= 1, "Tensor terms are not supported by this backend."); + + constexpr bool hdiv = Traits::family == functional_space::Family::hdiv; + static_assert(!hdiv || Traits::rank == 1, "H(div) bindings require rank one."); + + const int vector_dimension = hdiv || Traits::rank == 0 ? 1 : 3; + const int continuity = hdiv ? ::mfem::FiniteElementCollection::NORMAL : (Traits::family == functional_space::Family::h1 ? ::mfem::FiniteElementCollection::CONTINUOUS: ::mfem::FiniteElementCollection::DISCONTINUOUS); + + serif_ensure(space.GetVDim() == vector_dimension,"Term vector dimension disagrees with its scalar-copy or intrinsic-vector representation."); + serif_ensure(space.FEColl()->GetContType() == continuity,"Term continuity disagrees with its descriptor."); + + if constexpr (hdiv) { + serif_ensure(dynamic_cast(space.FEColl()) != nullptr,"H(div) bindings currently require an RT collection."); + } + + for (int element = 0; element < space.GetNE(); ++element) { + const auto* finite_element = space.GetFE(element); + + serif_ensure( + finite_element != nullptr && finite_element->GetOrder() == Traits::order + (hdiv ? 1 : 0) && + finite_element->GetRangeType() == (hdiv ? ::mfem::FiniteElement::VECTOR : ::mfem::FiniteElement::SCALAR) && + finite_element->GetMapType() == (hdiv ? ::mfem::FiniteElement::H_DIV : ::mfem::FiniteElement::VALUE), + "Term elements must match the registered family order, range, and mapping (RT uses FE order p+1)."); + } + } + + /** + * @brief Construct a term binding from an MFEM mesh. This requires that the term be a valid global scalar term (@ref GlobalScalarTerm). E.g. the only way you + * can construct such a term binding is for things like the total mass, angular momentum, etc... + * @param mesh The mesh to bind to, note that the mesh must be three-dimensional and have a space dimension of three. + */ + explicit TermBinding( + ::mfem::Mesh& mesh + ) requires discretization::blocks::GlobalScalarTerm : + DofMap(mesh, element_support(mesh)) {} + + private: + /** + * @brief Check that the mesh associated with a finite-element space is valid for term binding construction. + * @param space The finite-element space to check. + * @return The mesh associated with the finite-element space. + */ + static ::mfem::Mesh& checked_mesh(const ::mfem::FiniteElementSpace& space) { + using namespace utils::errors; + serif_ensure(space.GetMesh() != nullptr && space.FEColl() != nullptr, "Term bindings require an initialized finite-element space."); + return *space.GetMesh(); + } + + /** + * @brief Translate the mesh into an element support array for the particular term binding + * @param mesh the mesh to translate into an element support array + * @return An array where each entry is 1 if the corresponding element is supported by the term binding and 0 otherwise. This array is parallel to the mesh elements ordering array. + */ + static ::mfem::Array element_support(const ::mfem::Mesh& mesh) { + using namespace utils::errors; + constexpr auto domains = Schema::domain_descriptors(); + ::mfem::Array result(mesh.GetNE()); + + auto ID_matches = [](const auto& domain, const int attribute) -> bool { + return domain.ID == attribute; + }; + + for (int element = 0; element < mesh.GetNE(); ++element) { + const int attribute = mesh.GetAttribute(element); + serif_ensure(std::ranges::any_of(domains, [&](const auto& domain) { return ID_matches(domain, attribute); }), "Element attribute does not match the registered domain."); + result[element] = Schema::template domain_id_belongs_to(attribute); + } + return result; + } + }; +} diff --git a/src/lib/serif/discretization/backend/mfem/basis_evaluation.cpp b/src/lib/serif/discretization/backend/mfem/basis_evaluation.cpp index e69de29b..2b469984 100644 --- a/src/lib/serif/discretization/backend/mfem/basis_evaluation.cpp +++ b/src/lib/serif/discretization/backend/mfem/basis_evaluation.cpp @@ -0,0 +1,174 @@ +#include "mfem/fem/fe/fe_base.hpp" +#include "mfem/fem/intrules.hpp" +#include "mfem/linalg/vector.hpp" +#include "mfem/linalg/densemat.hpp" + +#include "serif/discretization/backend/mfem/basis_evaluation.hpp" + +#include "serif/discretization/backend/mfem/tensor_basis.hpp" +#include "serif/utils/error/errors.hpp" + +namespace serif::discretization::backend::mfem { + void BasisEvaluation::sample( + const ::mfem::FiniteElement& element, + const ::mfem::IntegrationRule& rule, + const ReferenceOperation operation, + const int copies, + TensorBasisCache& cache, + const ContractionMode mode + ) { + using namespace utils::errors; + + const bool vector = element.GetRangeType() == ::mfem::FiniteElement::VECTOR; + const bool operation_supported = vector ? (operation == ReferenceOperation::value || operation == ReferenceOperation::divergence) : operation != ReferenceOperation::divergence; + const bool mapping_supported = vector ? copies == 1 && element.GetMapType() == ::mfem::FiniteElement::H_DIV : element.GetMapType() == ::mfem::FiniteElement::VALUE; + + serif_ensure(element.GetDim() == 3 && copies > 0 && operation_supported && mapping_supported, "Canonical basis sampling requires 3D scalar VALUE bases or intrinsic H(div) bases and a compatible reference operation."); + + const int samples = operation == ReferenceOperation::gradient ? 3 : operation == ReferenceOperation::value_and_gradient ? 4 : 1; + m_components = vector ? (operation == ReferenceOperation::value ? 3 : 1) : copies * samples; + m_points = rule.GetNPoints(); + m_width = copies * element.GetDof(); + m_tensor = mode == ContractionMode::automatic ? cache.get(element, rule, operation, copies) : nullptr; + + if (m_tensor) { + m_dense.SetSize(0, 0); + m_mapping.SetSize(m_points * m_components, m_components); + m_mapping = 0.0; + + for (int q = 0; q < m_points; ++q) { + for (int component = 0; component < m_components; ++component) { + m_mapping(q * m_components + component, component) = 1.0; + } + } + return; + } + + m_mapping.SetSize(0, 0); + m_dense.SetSize(m_points * m_components, m_width); + m_dense = 0.0; + + ::mfem::Vector shape(element.GetDof()); + ::mfem::DenseMatrix gradients(element.GetDof(), 3); + + for (int q = 0; q < m_points; ++q) { + if (vector) { + if (operation == ReferenceOperation::divergence) { + element.CalcDivShape(rule.IntPoint(q), shape); + } else { + element.CalcVShape(rule.IntPoint(q), gradients); + } + } else { + element.CalcShape(rule.IntPoint(q), shape); + if (samples > 1) { + element.CalcDShape(rule.IntPoint(q), gradients); + } + } + + for (int copy = 0; copy < copies; ++copy) { + for (int sample = 0; sample < (vector ? m_components : samples); ++sample) { + const int derivative = operation == ReferenceOperation::gradient ? sample : sample - 1; + + for (int dof = 0; dof < element.GetDof(); ++dof) { + double value; + + if (vector) { + value = operation == ReferenceOperation::divergence ? shape[dof] : gradients(dof, sample); + } else { + value = samples == 1 || derivative < 0 ? shape[dof] : gradients(dof, derivative); + } + + m_dense(q * m_components + copy * samples + sample, copy * element.GetDof() + dof) = value; + } + } + } + } + } + + void BasisEvaluation::set_constant(const int points) { + m_tensor.reset(); + m_components = 1; + m_points = points; + m_width = 1; + m_dense.SetSize(points, 1); + m_dense = 1.0; + m_mapping.SetSize(0, 0); + } + + + void BasisEvaluation::reset_from(const BasisEvaluation& reference, const int components) { + m_tensor = reference.m_tensor; + m_points = reference.m_points; + m_width = reference.m_width; + m_components = components; + + if (m_tensor) { + m_mapping.SetSize(m_points * components, m_tensor->components()); + m_dense.SetSize(0, 0); + } else { + m_dense.SetSize(m_points * components, m_width); + m_mapping.SetSize(0, 0); + } + } + + void BasisEvaluation::evaluate(const ::mfem::Vector& coefficients, ::mfem::Vector& values, BasisWorkspace& workspace) const { + using namespace utils::errors; + + serif_ensure( + coefficients.Size() == m_width, + "Basis evaluation coefficient count does not match its element." + ); + values.SetSize(height()); + if (!m_tensor) { + m_dense.Mult(coefficients, values); + return; + } + + m_tensor->apply(coefficients, workspace.reference_values, workspace.tensor); + values = 0.0; + for (int q = 0; q < m_points; ++q) { + for (int output = 0; output < m_components; ++output) { + for (int input = 0; input < m_tensor->components(); ++input) { + values[q * m_components + output] += m_mapping(q * m_components + output, input) * workspace.reference_values[q * m_tensor->components() + input]; + } + } + } + } + + void BasisEvaluation::integrate(const ::mfem::Vector& values, ::mfem::Vector& coefficients, BasisWorkspace& workspace) const { + using namespace utils::errors; + + serif_ensure( + values.Size() == height(), + "Basis transpose quadrature size does not match its rule." + ); + coefficients.SetSize(m_width); + if (!m_tensor) { + m_dense.MultTranspose(values, coefficients); + return; + } + + workspace.reference_values.SetSize(m_points * m_tensor->components()); + workspace.reference_values = 0.0; + for (int q = 0; q < m_points; ++q) { + for (int input = 0; input < m_tensor->components(); ++input) { + for (int output = 0; output < m_components; ++output) { + workspace.reference_values[q * m_tensor->components() + input] += m_mapping(q * m_components + output, input) * values[q * m_components + output]; + } + } + } + m_tensor->apply(workspace.reference_values, coefficients, workspace.tensor, true); + } + + int BasisEvaluation::height() const noexcept { + return m_points * m_components; + } + + int BasisEvaluation::width() const noexcept { + return m_width; + } + + bool BasisEvaluation::uses_tensor() const noexcept { + return static_cast(m_tensor); + } +} \ No newline at end of file