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

refactor(discretization): cleaned up discretization and brought MFEM dependent code into backend

Browse Source
This commit is contained in:
Emily Boudreaux
2026-11-09 13:55:27 -05:00
parent 65761a3689
commit 994c818813
47 changed files with 261 additions and 5102 deletions
Download Patch File Download Diff File Expand all files Collapse all files

39
src/include/serif/dimensions/runtime/concepts.hpp
View File

@@ -1,39 +0,0 @@
#pragma once
/**
* @file concepts.hpp
* @author Emily M. Boudreaux
*
* @brief Runtime concept for a type that can be projected to a runtime type.
* This is useful when converting from the compile time type system to something
* that can be inspected at runtime.
*/
#include <string_view>
#include <type_traits>
#include "serif/dimensions/quantities.hpp"
#include "serif/dimensions/runtime/runtime.hpp"
namespace serif::dimensions::runtime {
/**
* @concept RuntimeIdentifiedThermodynamicQuantity
* @breif Concept to limit template arguments such that they are compatible with the expected
* ThermodynamicQuantity interface. Specifically, a valid type must contain a member
* identifier which can be converted to std::string_view and that member, when converted
* to a string_view must be longer than 0 length.
*/
template <typename Quantity>
concept RuntimeIdentifiedThermodynamicQuantity =
ThermodynamicQuantityType<Quantity> &&
requires {{Quantity::identifier} -> std::convertible_to<std::string_view>; } &&
(std::string_view{Quantity::identifier}.size() > 0);
/**
* @brief Extract a compile time identifier from a thermodynamic quantity and project it into a runtime
* string_view
* @tparam Quantity Type which @ref RuntimeIdentifiedThermodynamicQuantity allows
*/
template <RuntimeIdentifiedThermodynamicQuantity Quantity>
inline constexpr ThermodynamicQuantityID thermodynamicQuantityID{std::string_view{Quantity::identifier}};
}

61
src/include/serif/dimensions/runtime/runtime.hpp
View File

@@ -1,61 +0,0 @@
#pragma once
/**
* @file runtime.hpp
* @author Emily M. Boudreaux
*
* @brief runtime usable types representing the compile time dimension type system.
* Note that generally this is not a dimensional analysis library. Though perhaps
* in the future it will grow to something like that. Rather, this system
* is intended to ensure that developers are explicit which what numeric values mean.
* There is more detail on this in the file comment for @ref dimensions.hpp and @ref quantities.hpp
*/
#include <string_view>
#include <span>
#include <cstdint>
namespace serif::dimensions::runtime {
/**
* @class ThermodynamicQuantityID
*
* @brief Runtime resolvable ID for a thermodynamic quantity. Used in conjunction with
* concepts to resolve compile time thermodynamic types to runtime values.
*/
class ThermodynamicQuantityID final {
public:
explicit constexpr ThermodynamicQuantityID(const std::string_view name) noexcept : m_name(name) {};
[[nodiscard]] constexpr std::string_view name() const noexcept {
return m_name;
}
[[nodiscard]] friend constexpr bool operator==(const ThermodynamicQuantityID &, const ThermodynamicQuantityID &) noexcept = default;
private:
std::string_view m_name; // We can use a string view here since the constructor is constexpr and the string view will be valid for the lifetime of the program.
};
/**
* @struct RuntimeQuantityValue
*
* @brief Runtime resolved type for a quantity value. This holds both an ID and a value. This is the runtime
* equivalent to a @ref QuantityValue. Note that where these differ is that QuantityValue
* is built with template types. This allows the compiler to resolve it to a plain numeric type and still do optimizations
* That is generally not the case with RuntimeQuantityValue
*/
struct RuntimeQuantityValue final {
ThermodynamicQuantityID id;
double value{};
};
// TODO: This should be refactored into the eos module rather than the dimensions module
struct RuntimeRelationDescriptor final {
ThermodynamicQuantityID outputQuantity;
std::span<const ThermodynamicQuantityID> inputQuantities;
std::uint64_t partialDerivativeMask;
[[nodiscard]] constexpr bool hasPartialDerivative(const std::size_t inputIndex) const noexcept {
return inputIndex < inputQuantities.size() && (partialDerivativeMask & (std::uint64_t{1} << inputIndex)) != 0;
}
};
}

231
src/include/serif/discretization/backend/mfem/basis_evaluation.hpp
View File

@@ -1,231 +0,0 @@
#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 <memory>
#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 <typename MapT>
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<const TensorBasis> 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 <typename MapT>
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];
}
}
}

95
src/include/serif/discretization/backend/mfem/boundary_trace.hpp
View File

@@ -1,95 +0,0 @@
#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;
}

117
src/include/serif/discretization/backend/mfem/coefficient_tracking.hpp
View File

@@ -1,117 +0,0 @@
#pragma once
/**
* @file coefficient_tracking.hpp
* @brief Snapshot and comparison policy for MFEM coefficient views.
*/
#include <algorithm>
#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 <typename ViewT>
[[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 <typename ViewT>
[[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 <typename ViewT>
[[nodiscard]] static const void* identity(const ViewT& view);
};
} // namespace serif::discretization::backend::mfem
template <typename ViewT>
::mfem::Vector serif::discretization::backend::mfem::CoefficientTracking::capture(const ViewT& view) {
return ::mfem::Vector(view.values);
}
template <typename ViewT>
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 <typename ViewT>
[[nodiscard]] const void* serif::discretization::backend::mfem::CoefficientTracking::identity(const ViewT& view) {
return &view.binding;
}

306
src/include/serif/discretization/backend/mfem/dof_map.hpp
View File

@@ -1,306 +0,0 @@
#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 <cmath>
#include <optional>
#include <stdexcept>
#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<int>& 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<SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(dynamic_cast<::mfem::ParFiniteElementSpace*>(&space) == nullptr,"DofMap does not yet assemble parallel true DOFs.");
#endif
serif_ensure<SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(space.Conforming() && space.GetTrueVSize() == space.GetVSize(), "Term bindings currently require conforming serial spaces.");
serif_ensure<SERiFErrorCode::INVALID_DIMENSION, std::invalid_argument>(space.GetMesh()->Dimension() == 3 && space.GetMesh()->SpaceDimension() == 3,"Term bindings currently require a three-dimensional mesh.");
serif_ensure<SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(element_support.Size() == space.GetNE(), "Element support does not match the space.");
m_full_to_reduced = -1;
::mfem::Array<int> 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<int>& 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<SERiFErrorCode::INVALID_DIMENSION, std::invalid_argument>(mesh.Dimension() == 3 && mesh.SpaceDimension() == 3, "Term bindings currently require a three-dimensional mesh.");
serif_ensure<SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(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<SERiFErrorCode::CACHE_FE_BINDING_ERROR, std::logic_error>(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<int> reduced_dof(const int full_dof) const {
using namespace utils::errors;
serif_ensure<SERiFErrorCode::INVALID_DOF_COUNT, std::out_of_range>(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<int>{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<SERiFErrorCode::INVALID_ELEMENT_ID_ERROR, std::out_of_range>(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<SERiFErrorCode::CACHE_FE_BINDING_ERROR, std::logic_error>(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<SERiFErrorCode::CACHE_FE_BINDING_ERROR>(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<SERiFErrorCode::INVALID_DOF_COUNT, std::invalid_argument>(values.Size() == size(), "Term binding size is {}, expected {} supported DOFs.", values.Size(), size());
for (int i = 0; i < values.Size(); ++i) {
serif_ensure<SERiFErrorCode::NON_FINITE_VALUE_ERROR, std::domain_error>(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<SERiFErrorCode::INVALID_DOF_COUNT, std::invalid_argument>(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<SERiFErrorCode::INVALID_DOF_COUNT, std::invalid_argument>(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<int> 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<SERiFErrorCode::INVALID_DOF_COUNT, std::invalid_argument>(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<int> dofs;
const auto* transformation = m_space->GetElementVDofs(element, dofs);
serif_ensure<SERiFErrorCode::INVALID_DOF_COUNT, std::invalid_argument>(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<int> m_attributes; ///< Array of attribute indices associated with this dof map
::mfem::Array<int> m_element_support; ///< Array of 0/1 values indicating whether each element is supported by this dof map
::mfem::Array<int> m_reduced_to_full; ///< Mapping from reduced to full degrees of freedom
::mfem::Array<int> m_full_to_reduced; ///< Mapping from full to reduced degrees of freedom
};
} // namespace serif::discretization::backend::mfem

155
src/include/serif/discretization/backend/mfem/dof_selection.hpp
View File

@@ -1,155 +0,0 @@
#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<int>& 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<SERiFErrorCode::INVALID_DOF_COUNT, std::invalid_argument>(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<int>& 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<int>& 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<SERiFErrorCode::CACHE_FE_BINDING_ERROR, std::logic_error>(m_nodes_sequence < 0 || mesh.GetNodesSequence() == m_nodes_sequence, "Rebuild point selections after fixed-mesh coordinate changes.");
serif_ensure<SERiFErrorCode::CACHE_FE_BINDING_ERROR, std::logic_error>(mesh.GetNBE() == m_boundary_attributes.Size(), "Rebuild DOF selections after boundary changes.");
for (int boundary = 0; boundary < m_boundary_attributes.Size(); ++boundary) {
serif_ensure<SERiFErrorCode::CACHE_FE_BINDING_ERROR, std::logic_error>(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<int> m_dofs; ///< List of reduced DOFs in this selection. These are in the reduced space of the binding provided at construction time.
::mfem::Array<int> 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 <domain::IsBoundary BoundaryT, blocks::SpatialTerm TermT, domain::schema::IsSchema SchemaT>
[[nodiscard]] DofSelection select_boundary_dofs(
const TermBinding<TermT, SchemaT>& binding
) {
using namespace utils::errors;
static_assert(SchemaT::template contains_boundary<BoundaryT>(), "The boundary must be registered in the domain schema.");
binding.validate();
::mfem::Array<int> marker(binding.size());
marker = 0;
::mfem::Array<int> dofs;
::mfem::Array<int> selected;
auto& space = binding.space();
const auto& mesh = binding.mesh();
for (int boundary = 0; boundary < mesh.GetNBE(); ++boundary) {
if (!SchemaT::template boundary_id_matches<BoundaryT>(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<SERiFErrorCode::INVALID_DOMAIN_ERROR, std::invalid_argument>(reduced.has_value(), "A selected boundary DOF lies outside the term's support.");
marker[*reduced] = 1;
}
}
::mfem::FiniteElementSpace::MarkerToList(marker, selected);
serif_ensure<SERiFErrorCode::INVALID_DOMAIN_ERROR, std::invalid_argument>(selected.Size() > 0, "The requested boundary has no DOFs on this serial mesh.");
return DofSelection(binding, selected);
}
} // namespace serif::discretization::backend::mfem

28
src/include/serif/discretization/backend/mfem/mass_diagonal.hpp
View File

@@ -1,28 +0,0 @@
#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);
}
}

20
src/include/serif/discretization/backend/mfem/mfem_backend.hpp
View File

@@ -1,20 +0,0 @@
#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"

82
src/include/serif/discretization/backend/mfem/point_dofs.hpp
View File

@@ -1,82 +0,0 @@
#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 <array>
#include <cmath>
#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 <blocks::SpatialTerm TermT, domain::schema::IsSchema SchemaT>
requires (blocks::SpatialTermTraits<TermT>::family == functional_space::Family::h1)
[[nodiscard]] DofSelection select_point_dofs(
const TermBinding<TermT, SchemaT>& binding,
const std::array<double, 3>& 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<SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(std::isfinite(tolerance) && tolerance >= 0.0, "Point selection requires a finite nonnegative coordinate tolerance.");
for (const double coordinate : position) {
serif_ensure<SERiFErrorCode::NON_FINITE_VALUE_ERROR, std::invalid_argument>(std::isfinite(coordinate), "Point selection requires finite reference coordinates.");
}
auto& space = binding.space();
for (int element = 0; element < space.GetNE(); ++element) {
serif_ensure<SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(dynamic_cast<const ::mfem::NodalFiniteElement*>(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<double, 3> 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<SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(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<SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(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<int> vertex_dofs, selected;
space.GetVertexDofs(selected_vertex, vertex_dofs);
serif_ensure<SERiFErrorCode::INVALID_DOF_COUNT, std::invalid_argument>(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<SERiFErrorCode::INVALID_DOMAIN_ERROR, std::invalid_argument>(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

489
src/include/serif/discretization/backend/mfem/tensor_basis.hpp
View File

@@ -1,489 +0,0 @@
#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 <algorithm>
#include <array>
#include <cmath>
#include <limits>
#include <memory>
#include <vector>
#include <utility>
#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<int> 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<const TensorBasis> make(
const ::mfem::FiniteElement& element,
const ::mfem::IntegrationRule& rule,
const ReferenceOperation operation,
const int copies
) {
const auto* scalar = dynamic_cast<const ::mfem::NodalTensorFiniteElement*>(&element);
const auto* rt = dynamic_cast<const ::mfem::RT_HexahedronElement*>(&element);
const bool scalar_hex = dynamic_cast<const ::mfem::H1_HexahedronElement*>(&element) || dynamic_cast<const ::mfem::L2_HexahedronElement*>(&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<int>(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<TensorBasis>();
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<bool> since it is the most obvious conceptual datatype. However, the standard backed implementation of
// vector<bool> 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<bool> 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<SERiFErrorCode::INVALID_DOF_COUNT, std::invalid_argument>(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<int, 3> 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<int, 3>& 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<int, 3> 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<double>::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<Product> 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<std::array<double, 4>> points;
std::shared_ptr<const TensorBasis> 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<const TensorBasis> get(
const ::mfem::FiniteElement& element,
const ::mfem::IntegrationRule& rule,
ReferenceOperation operation,
const int copies
) {
std::vector<std::array<double, 4>> 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<Entry> 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

115
src/include/serif/discretization/backend/mfem/term_binding.hpp
View File

@@ -1,115 +0,0 @@
#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 <algorithm>
#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 <typename TermT, domain::schema::IsSchema SchemaT>
requires (blocks::SpatialTerm<TermT> || blocks::GlobalScalarTerm<TermT>)
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<Term> :
DofMap(space, element_support(checked_mesh(space))) {
using Traits = blocks::SpatialTermTraits<Term>;
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<SERiFErrorCode::INVALID_DIMENSION, std::invalid_argument>(space.GetVDim() == vector_dimension,"Term vector dimension disagrees with its scalar-copy or intrinsic-vector representation.");
serif_ensure<SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(space.FEColl()->GetContType() == continuity,"Term continuity disagrees with its descriptor.");
if constexpr (hdiv) {
serif_ensure<SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(dynamic_cast<const ::mfem::RT_FECollection*>(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<SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(
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<Term> :
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<SERiFErrorCode::CACHE_FE_BINDING_ERROR>(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<int> element_support(const ::mfem::Mesh& mesh) {
using namespace utils::errors;
constexpr auto domains = Schema::domain_descriptors();
::mfem::Array<int> 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<SERiFErrorCode::INVALID_DOMAIN_ERROR>(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<typename Term::Support>(attribute);
}
return result;
}
};
}

23
src/include/serif/discretization/blocks/base.hpp
View File

@@ -41,29 +41,6 @@ namespace serif::discretization::blocks {
static constexpr int static_block_size = DYNAMIC_BLOCK_SIZE;
};
/**
* @struct generated_value_block
* @brief Final class for generated value blocks, derived from value_block_base. These are blocks that are generated at runtime
* @todo This may be removable at this point
*/
template <typename GeneratedValue>
struct generated_value_block final : value_block_base {
using GeneratedType = GeneratedValue;
static constexpr int static_block_size = static_cast<int>(GeneratedValue::scalarArity);
};
/**
* @struct generated_residual_block
* @brief Final class for generated residual blocks, derived from residual_block_base. These are blocks that are generated at runtime
* @todo This may be removable at this point
*/
template <typename GeneratedResidual>
struct generated_residual_block final : residual_block_base {
using GeneratedType = GeneratedResidual;
static constexpr int static_block_size = static_cast<int>(GeneratedResidual::scalarArity);
};
/**
* @struct residual_block
* @brief Final class for residual blocks with a fixed index, derived from residual_block_base.

29
src/include/serif/discretization/blocks/concepts.hpp
View File

@@ -87,33 +87,4 @@ namespace serif::discretization::blocks {
*/
template <typename T>
concept GlobalScalarTerm = Term<T> && HasSupport<T> && !HasFunctionalSpace<T> && std::derived_from<std::remove_cvref_t<T>, global_scalar>;
/**
* @brief Concept for a physical spatial term. A physical spatial term is a spatial term (@ref serif::discretization::blocks::SpatialTerm "SpatialTerm") that has a physical quantity associated with it (@ref serif::discretization::blocks::HasQuantity "HasQuantity").
*/
template <typename T>
concept PhysicalSpatialTerm =
SpatialTerm<T> &&
HasQuantity<T>;
/**
* @brief Helper struct to determine if a type is a spatial term (@ref serif::discretization::blocks::SpatialTerm "SpatialTerm").
* @tparam T The type to check.
*/
template <typename T>
struct IsSpatialTerm : std::bool_constant<SpatialTerm<T>> {};
/**
* @brief Helper struct to determine if a type is a global scalar term (@ref serif::discretization::blocks::GlobalScalarTerm "GlobalScalarTerm").
* @tparam T The type to check.
*/
template <typename T>
struct IsGlobalScalarTerm : std::bool_constant<GlobalScalarTerm<T>> {};
/**
* @brief Helper struct to determine if a type is a physical spatial term (@ref serif::discretization::blocks::PhysicalSpatialTerm "PhysicalSpatialTerm").
* @tparam T The type to check.
*/
template <typename T>
struct IsPhysicalSpatialTerm : std::bool_constant<PhysicalSpatialTerm<T>> {};
}

119
src/include/serif/discretization/blocks/fields.hpp
View File

@@ -15,6 +15,7 @@
* if that belongs here instead.
*/
#include <array>
#include <string_view>
#include "serif/dimensions/dimensions.hpp"
@@ -70,27 +71,47 @@ namespace serif::discretization::blocks {
};
/**
* @brief The barotropic_constant field is a field that represents the barotropic constant of the system. It is a global scalar that is associated with the mass normalization term. This is used to enforce the mass normalization condition in the weak form when
* we are in the barotropic case. This field will not generalize as useful to the non-barotropic case. Note that we use compile time reflection to get the name of the field and the terms.
* @brief This field holds the effective potential level on the stellar isobar. The level is solved together with the fixed total mass.
* The barotropic_constant name is retained as an internal coordinate identity. In the density formulation its value is Us,
* the gravitational potential minus the centrifugal potential on the surface; it is not the polytropic EOS constant K.
*/
struct barotropic_constant final : field {
static constexpr std::string_view name = meta::type_name<barotropic_constant>(); ///< Compile time reflected name of the field
/**
* @brief the mass_normalization term represents the mass normalization condition of the system. That is that it enforces that the total mass of the system is equal to the target mass in a weak sense.
* This is both a term and a global scalar.
* @brief The value of this global scalar is the effective surface-potential level Us. Its paired residual enforces the prescribed total mass.
* The surface equation supplies the dependence on Us, while the mass equation integrates the density over the star.
*/
struct mass_normalization final : term<barotropic_constant>, global_scalar {
static constexpr std::string_view name = meta::type_name<mass_normalization>(); ///< Compile time reflected name of the term
using Support = domain::StellarDomains; ///< The support of the mass normalization term is the stellar domains since we only want to enforce the mass normalization condition in the stellar domains.
using Quantity = dimensions::SpecificEnthalpy; ///< As we are in teh barotropic case the mass normalization term is associated with the specific enthalpy. This will nto generalize to the non-barotropic case.
static constexpr std::string_view symbol = "C"; ///< Compile time symbol used for the term
using Quantity = dimensions::GravitationalPotential; ///< Physical units of the solved effective potential level, distinct from the mass residual units.
static constexpr std::string_view symbol = "Us"; ///< Compile time symbol for the effective surface-potential level.
};
static constexpr mass_normalization mass_normalization_term{}; ///< Static constexpr instance of the mass_normalization term. This is used to access the term in the discretization.
using Terms = utils::types::TypeList<mass_normalization>; ///< Type list of terms for compile time reflection and iteration over the terms of the field.
};
/**
* @brief Uniform translation reaction, one global scalar per Cartesian axis, to be paired with BarycenterConstraint
* (for example when an external companion potential acts on the star). It must vanish in a physical equilibrium.
*/
struct frame_acceleration final : field {
static constexpr std::string_view name = "frame acceleration";
template <int Axis>
struct component final : term<frame_acceleration>, global_scalar {
static_assert(Axis >= 0 && Axis < 3, "A Cartesian reaction component requires an axis from 0 through 2.");
static constexpr std::string_view name = meta::type_name<component>();
static constexpr std::string_view stable_id = std::array<std::string_view, 3>{
"field.frame_acceleration.x", "field.frame_acceleration.y", "field.frame_acceleration.z"}[Axis];
static constexpr std::string_view symbol = std::array<std::string_view, 3>{"gamma_x", "gamma_y", "gamma_z"}[Axis];
using Support = domain::StellarDomains;
using Quantity = dimensions::Acceleration;
};
using Terms = utils::types::TypeList<component<0>, component<1>, component<2>>;
};
/**
* @brief The rotational field, represented as the angular velocity term, given that we currently only model rigid rotation this is a global scalar. A future extension which models differential rotation will
* likely need to modify this to some vector field or other more complex representation of the rotation.
@@ -113,26 +134,6 @@ namespace serif::discretization::blocks {
using Terms = utils::types::TypeList<angular_velocity>; ///< Type list of terms for compile time reflection and iteration over the terms of the field.
};
/**
* @brief field used to inject a border row used to enforce a central density constraint
*/
struct central_density final : field {
static constexpr std::string_view name = meta::type_name<central_density>(); ///< Compile time reflected name of the field
/**
* @brief The border term represents the border row used to enforce a central density constraint. This is a global scalar as it is a single value over the entire domain.
*/
struct border final : term<central_density>, global_scalar {
static constexpr std::string_view name = meta::type_name<border>(); ///< Compile time reflected name of the term
using Support = domain::StellarDomains; ///< We only enforce this over the stellar domain as there is only material to "slosh" around in the stellar domain
using Quantity = dimensions::SpecificEnthalpy; ///< As we are in the barotropic case this border term is associate with the specific enthalpy.
static constexpr std::string_view symbol = "ρ_c"; ///< Compile time symbol used for the term
};
static constexpr border border_term{}; ///< Static constexpr instance of the border term. This is used to access the term in the discretization.
using Terms = utils::types::TypeList<border>; ///< Type list of terms for compile time reflection and iteration over the terms of the field.
};
/**
* @brief The density field used to represent the distribution of mass in the system. This is a scalar field that is associated with the mass density of the system. The density field is used to enforce the mass conservation condition in the weak form.
*/
@@ -141,12 +142,12 @@ namespace serif::discretization::blocks {
/**
* @brief The mass density term for the system. Note that the naming here could be a bit confusing as FEM also uses the term "mass" to represent the canonical mass matrix. This is a more standard physics reading of what is going on here.
* We define that this space is a second order scalar L2 space. This allows for the possibility of discontinuous density fields.
* We define that this space is a third order scalar H1 space. This gives density a continuous surface trace and a weak spatial gradient.
*/
struct mass final : term<density> {
static constexpr std::string_view name = meta::type_name<mass>();
using FunctionalSpace = functional_space::DiscreteFunctionalSpace<functional_space::L2, 2, 0>; ///< second order scalar L2 space. Read the template parameters as follows: 1) L2 space, 2) second order, 3) scalar field (rank 0)
using FunctionalSpace = functional_space::DiscreteFunctionalSpace<functional_space::H1, 3, 0>; ///< third order scalar H1 space. Read the template parameters as follows: 1) H1 space, 2) third order, 3) scalar field (rank 0)
using Support = domain::StellarDomains; ///< The support of the mass density term is the stellar domains since we do not allow material to exist in the vacuum domain.
using Quantity = dimensions::MassDensity; ///< The quantity of the mass density term is the mass density.
static constexpr std::string_view symbol = "ρ"; ///< Compile time symbol used for the term. Note that this is a UTF-8 character and may not render correctly in all environments. (TODO: Should we introduce an ASCII fallback?, something like make this a tuple of form <main, fallback>?)
@@ -225,40 +226,19 @@ namespace serif::discretization::blocks {
};
/**
* @brief The gravity field is represented as a mixed system, in many ways this field is the motivation for the entire design of this file. Since it is a
* mixed system we have two separate terms that make it up, the potential and the gradient. These are not derived from one and other, rather they are
* weakly enforced to agree. This buys us an additional order of convergence in the potential which in tern allows us to use a lower resolution mesh.
* @brief The gravity field contains the scalar potential. The primal formulation solves the potential and obtains its gradient by differentiation.
*/
struct gravity final : field {
static constexpr std::string_view name = meta::type_name<gravity>(); ///< Compile time reflected name of the field
/**
* @brief The gravitational gradient term is a 2nd order Hdiv vector field that is used to represent the negative gravitational acceleration of the system.
* Note that when I started developing this code I made the frankly wacky choice to let this be the true gradient of the potential, which is the negative
* acceleration. Broadly this may make sense mathematically; however, it means that developers must be extra careful to remember the sign convention.
*/
struct gradient final : term<gravity> {
static constexpr std::string_view name = meta::type_name<gradient>(); ///< Compile time reflected name of the term
using FunctionalSpace = functional_space::DiscreteFunctionalSpace<functional_space::Hdiv, 2, 1>; ///< 2nd order Hdiv vector field. Read the template parameters as follows: 1) Hdiv space, 2) 2nd order, 3) rank 1 tensor field
using Support = domain::AllDomains; ///< We support gravity on all domains, primarily since the only real boundary condition we have in the general case is that the potential goes to zero at infinity, which is in the vacuum domain.
using Quantity = dimensions::Acceleration; ///< This is a dimension of acceleration, despite the wacky sign convention I adopted.
static constexpr std::string_view symbol = "∇φ"; ///< Compile time symbol used for the term. Note that this is a UTF-8 character and may not render correctly in all environments. (TODO: Should we introduce an ASCII fallback?, something like make this a tuple of form <main, fallback>?)
struct value final : value_block_base {}; ///< The value block is a block that contains the value of the gravitational gradient term. It is a vector field that is associated with the negative gravitational acceleration of the system.
struct residual final : residual_block_base {}; ///< The residual block is a block that contains the residual of the gravitational gradient term. It is a vector field that is associated with the negative gravitational acceleration of the system.
};
/**
* @brief The gravitational potential term is a 2nd order L2 scalar field that is used to represent the gravitational potential of the system. Note that
* this is not exactly the same as the if you were to take the scalar potential of the gravitational gradient term. These are solved for by the code separately and
* weakly enforced to agree. Generally if you need the most accurate gradient you should use the gradient term and if you need the most accurate
* potential you should use this term (this is because differentiating this term loses an order of accuracy in the differentiation).
* @brief The gravitational potential term is a fourth order H1 scalar field that is used to represent the gravitational potential of the system.
* Its continuous trace supports the surface potential condition, and its spatial gradient supplies gravity in the primal formulation.
*/
struct potential final : term<gravity> {
static constexpr std::string_view name = meta::type_name<potential>(); ///< Compile time reflected name of the term
using FunctionalSpace = functional_space::DiscreteFunctionalSpace<functional_space::L2, 2, 0>; ///< 2nd order L2 scalar field. Read the template parameters as follows: 1) L2 space, 2) 2nd order, 3) rank 0 tensor field
using FunctionalSpace = functional_space::DiscreteFunctionalSpace<functional_space::H1, 4, 0>; ///< fourth order H1 scalar field. Read the template parameters as follows: 1) H1 space, 2) fourth order, 3) rank 0 tensor field
using Support = domain::AllDomains; ///< We support gravity on all domains, primarily since the only real boundary condition we have in the general case is that the potential goes to zero at infinity, which is in the vacuum domain.
using Quantity = dimensions::GravitationalPotential; ///< This is a dimension of gravitational potential, which is a scalar field that is associated with the gravitational potential of the system.
static constexpr std::string_view symbol = "φ"; ///< Compile time symbol used for the term. Note that this is a UTF-8 character and may not render correctly in all environments. (TODO: Should we introduce an ASCII fallback?, something like make this a tuple of form <main, fallback>?)
@@ -267,38 +247,7 @@ namespace serif::discretization::blocks {
struct residual final : residual_block_base {}; ///< The residual block is a block that contains the residual of the gravitational potential term. It is a scalar field that is associated with the gravitational potential of the system.
};
static constexpr potential potential_term{}; ///< Static constexpr instance of the potential term. This is used to access the term in the discretization.
static constexpr gradient gradient_term{}; ///< Static constexpr instance of the gradient term. This is used to access the term in the discretization.
using Terms = utils::types::TypeList<potential, gradient>; ///< Type list of terms for compile time reflection and iteration over the terms of the field.
};
/**
* @brief The enthalpy field is a scalar field is a C0 continuous scalar field. Note that this field likely will not survive the transition to the more general
* non-barotropic case. In the barotropic case the gradient of enthalpy is parallel to the gradient of the pressure. This make enthalpy a particularly
* useful quantity as we can represent the 0 pressure surface condition as the 0 enthalpy condition and we can break down mass constraint into
* a game of optimizing the Bernoulli constant.
*/
struct enthalpy final : field {
static constexpr std::string_view name = meta::type_name<enthalpy>(); ///< Compile time reflected name of the field
/**
* @brief the specific enthalpy is the term we specifically work with rather than the total enthalpy of dof. In fact heuristically I would
* say that generally we should always work in specific quantities rather than total quantities as there is not a well defined amount of material
* at some dof.
*/
struct specific final : term<enthalpy> {
static constexpr std::string_view name = meta::type_name<specific>(); ///< Compile time reflected name of the term
using FunctionalSpace = functional_space::DiscreteFunctionalSpace<functional_space::H1, 3, 0>; //< 3 vector H1 space (C0 continuous). Read the template parameters as follows: 1) H1 space, 2) 3rd order, 3) rank 0 tensor field
using Support = domain::StellarDomains; ///< There is only material in the stellar domain therefore we only support the enthalpy field in the stellar domain.
using Quantity = dimensions::SpecificEnthalpy; ///< This has the dimension of specific enthalpy,
static constexpr std::string_view symbol = "h"; ///< Compile time symbol used for the term.
struct value final : value_block_base { }; ///< The value block is a block that contains the value of the specific enthalpy term. It is a scalar field that is associated with the specific enthalpy of the system.
struct residual final : residual_block_base { }; ///< The residual block is a block that contains the residual of the specific enthalpy term. It is a scalar field that is associated with the specific enthalpy of the system.
};
static constexpr specific specific_term{}; ///< Static constexpr instance of the specific term. This is used to access the term in the discretization.
using Terms = utils::types::TypeList<specific>; ///< Type list of terms for compile time reflection and iteration over the terms of the field.
using Terms = utils::types::TypeList<potential>; ///< Type list of terms for compile time reflection and iteration over the terms of the field.
};
}

47
src/include/serif/discretization/blocks/selection.hpp
View File

@@ -1,47 +0,0 @@
#pragma once
#include <limits>
#include <vector>
#include "serif/discretization/blocks/system_layout.hpp"
#include "serif/utils/error/errors.hpp"
namespace serif::discretization::blocks {
/**
* @brief Select named packed coordinates in their requested order, independently for values and residuals.
*
* @tparam Kind The kind of represented values to select, e.g. RepresentedValueKind::values or RepresentedValueKind::residuals.
* @tparam Layout The block layout type, e.g. SystemLayout.
* @tparam Coordinates The coordinate tags to select, e.g. keys::Density,
*
* @param layout The block layout from which to select coordinates.
* @param keys A list of coordinate tags to select. Note that this is not used at runtime, but rather to select the correct overloads of layout.size and layout.offset.
*
* @returns A vector of packed indices corresponding to the requested coordinates in the order requested.
*
* Given some layout (see @ref SystemLayout) and a list of coordinate tags return the corresponding packed indices in the order requested.
* The result is suitable for indexing a contiguous vector of values or residuals.
*/
template <RepresentedValueKind Kind, typename Layout, typename... Coordinates>
[[nodiscard]] std::vector<int> coordinate_indices(
const Layout& layout,
utils::dependencies::KeyList<Coordinates...> /*keys*/
) {
utils::errors::serif_ensure<utils::errors::SERiFErrorCode::INVALID_DOF_COUNT>(layout.template total_size<Kind>() <= static_cast<std::size_t>(std::numeric_limits<int>::max()),"Named coordinate selection exceeds the backend index range; use a layout with representable packed indices.");
std::vector<int> result;
result.reserve((std::size_t{0} + ... + layout.template size<Kind, Coordinates>()));
const auto append = [&]<typename Coordinate>() {
const int offset = static_cast<int>(layout.template offset<Kind, Coordinate>());
const int size = static_cast<int>(layout.template size<Kind, Coordinate>());
for (int index = 0; index < size; ++index) {
result.push_back(offset + index);
}
};
(append.template operator()<Coordinates>(), ...);
return result;
}
}

16
src/include/serif/discretization/blocks/traits.hpp
View File

@@ -11,13 +11,6 @@
#include "serif/discretization/blocks/concepts.hpp"
namespace serif::discretization::blocks {
/**
* @brief Extract the field type associated with a given term type. Will have a type of whatever TermT::field_type is.
* @tparam TermT The term type from which to extract the field type.
*/
template <Term TermT>
using FieldOf = std::remove_cvref_t<TermT>::field_type;
/**
* @brief Extract the functional space type associated with a given term type. Will have a type of whatever TermT::FunctionalSpace is.
* @tparam TermT The term type from which to extract the functional space type.
@@ -45,15 +38,6 @@ namespace serif::discretization::blocks {
template <HasQuantity TermT>
using QuantityOf = std::remove_cvref_t<TermT>::Quantity;
/**
* @brief Extract the terms associated with a given field type. Will have a type of whatever FieldT::Terms is.
* @tparam FieldT The field type from which to extract the terms.
*
* @todo For consistency with other types consider renaming FieldT::Terms to FieldT::terms. This would require updating all field definitions to use the new naming convention. CLion or some other IDE should be able to handle this refactor easily.
*/
template <typename FieldT>
using TermsOf = FieldT::Terms;
/**
* @brief Helper struct representing the traits of some spatial term type, including its functional space, support, family, order, and rank.
* @tparam TermT The spatial term type from which to extract the functional space and support types.

6
src/include/serif/discretization/domain/ids/concepts.hpp
View File

@@ -49,10 +49,4 @@ namespace serif::discretization::domain::ids {
*/
template <typename T>
concept IsBoundaryID = is_boundary_id_v<T>;
/**
* @brief Concept for types that are either domain IDs or boundary IDs in the discretization framework.
*/
template <typename T>
concept IsDomainOrBoundaryID = IsDomainID<T> || IsBoundaryID<T>;
}

34
src/include/serif/discretization/domain/schema/utils.hpp
View File

@@ -1,34 +0,0 @@
#pragma once
/**
* @file utils.hpp
* @author Emily M. Boudreaux
*
* @brief File containing utility functions for working with DomainSchemas.
*/
#include <mfem.hpp>
#include "serif/discretization/domain/concepts.hpp"
#include "serif/discretization/domain/schema/concepts.hpp"
namespace serif::discretization::domain::schema {
/**
* @brief Build the 0/1 marker array MFEM wants when an integrator or an essential
* BC should apply only on part of the mesh.
* @tparam DomainT The domain or set type.
* @tparam SchemaT The schema type.
* @param mesh The MFEM mesh.
* @return An MFEM Array<int> with 0/1 markers.
*/
template <IsDomainOrSet DomainT, IsSchema SchemaT>
[[nodiscard]] mfem::Array<int> make_attribute_marker(const mfem::Mesh &mesh) {
static_assert(SchemaT::template contains_domain<DomainT>(), "Requested marker domain is not completely registered in the supplied DomainSchema.");
mfem::Array<int> marker(mesh.attributes.Max());
for (int domainID = 1; domainID <= marker.Size(); ++domainID) {
marker[domainID - 1] = SchemaT::template domain_id_belongs_to<DomainT>(domainID) ? 1 : 0;
}
return marker;
}
}

7
src/include/serif/discretization/domain/schema/validation/mesh.hpp
View File

@@ -1,7 +0,0 @@
#pragma once
#include "serif/backend/mfem/discretization/schema.hpp"
namespace serif::discretization::domain::schema::validation {
using serif::backend::mfem::discretization::validate_schema;
}

2
src/include/serif/discretization/forms/operands.hpp
View File

@@ -103,7 +103,7 @@ namespace serif::discretization::forms {
* - @ref Operation::NormalTrace (\f$\int_{\partial\Omega} \left(\hat{\mathbf{n}}\cdot\mathbf{v}_i\right)\lambda\,dS\f$)
*
* @note Declaring an operation tag does not itself provide a backend implementation. The current MFEM volume-equation backend supports
* Value and Divergence for its supported field spaces. Gradient, Curl, and NormalTrace require additional backend support;
* Value and Divergence for its supported field spaces, and Gradient for scalar/vector H1 fields. Curl and NormalTrace require additional backend support;
* NormalTrace also requires an appropriate boundary integration path.
*
* In the equation system you may find documentation for the "contributions" system, which uses the term in the test operand to determine

30
src/include/serif/discretization/forms/value.hpp
View File

@@ -69,6 +69,12 @@ namespace serif::discretization::forms {
using type = VectorValue;
};
/** @brief Vector gradients use G[component][derivative] = partial_derivative(u_component, x_derivative). */
template <>
struct OperationValue<1, Operation::Gradient> {
using type = MatrixValue;
};
/**
* @brief The curl of a rank 1 tensor is a rank 1 tensor (a vector of doubles)
*/
@@ -128,7 +134,20 @@ namespace serif::discretization::forms {
* @todo This introduces an assumption of three-dimensionality, we say "this is either a scalar or a tensor of dimension 3". We should probably relax this assumption.
*/
template <typename ValueT>
inline constexpr int component_count = std::same_as<ValueT, double> ? 1 : 3;
struct ValueComponents;
template <>
struct ValueComponents<double> : std::integral_constant<int, 1> {};
// OperandValues uses void to request values instead of derivatives; that branch stores no derivative entries.
template <>
struct ValueComponents<void> : std::integral_constant<int, 0> {};
template <typename Entry, std::size_t Size>
struct ValueComponents<std::array<Entry, Size>> : std::integral_constant<int, Size * ValueComponents<Entry>::value> {};
template <typename ValueT>
inline constexpr int component_count = ValueComponents<std::remove_cv_t<ValueT>>::value;
/**
* @brief Alias for the partial derivative of some output type with respect to some input type. This is a convenience alias that allows you to easily access the partial derivative type without having to explicitly use the OperandResult struct.
@@ -138,7 +157,8 @@ namespace serif::discretization::forms {
* @todo Should probably be consolidated with the other PartialDerivative system... No, should *definitely* be consolidated...
*/
template <typename OutputT, typename InputT>
using PartialDerivative = std::conditional_t<std::same_as<OutputT, double>, InputT, std::conditional_t<std::same_as<InputT, double>, OutputT, MatrixValue>>;
using PartialDerivative = std::conditional_t<std::same_as<OutputT, double>, InputT, std::conditional_t<std::same_as<InputT, double>, OutputT,
std::array<std::array<double, component_count<InputT>>, component_count<OutputT>>>>;
/**
* @brief Read one entry of a scalar/vector/matrix value. This is a convenience function that allows you to easily access the components of a value type without having to explicitly use the OperandResult struct.
@@ -152,7 +172,9 @@ namespace serif::discretization::forms {
if constexpr (std::same_as<std::remove_cv_t<ValueT>, double>) {
return (value);
} else {
return (value[index]);
using Entry = typename std::remove_cv_t<ValueT>::value_type;
constexpr int stride = component_count<Entry>;
return component(value[index / stride], index % stride);
}
}
@@ -170,7 +192,7 @@ namespace serif::discretization::forms {
if constexpr (std::same_as<OutputT, double>) {
return component(value, input);
} else if constexpr (std::same_as<InputT, double>) {
return value[output];
return component(value, output);
} else {
return value[output][input];
}

6
src/include/serif/discretization/functional_space/space.hpp
View File

@@ -151,12 +151,6 @@ namespace serif::discretization::functional_space {
.mapping = Mapping::covariant_piola
};
/**
* @brief Concept to validate a valid space descriptor. A space descriptor is valid if it is of type SpaceDescriptor. This is a very simple concept that just checks that the type is exactly SpaceDescriptor.
*/
template <typename SpaceDescriptorCandidate>
concept IsSpaceDescriptor = std::is_same_v<SpaceDescriptor, SpaceDescriptorCandidate>;
/**
* @brief Compile-time descriptor for a discrete functional space. The difference between a FunctionalSpace and a discrete functional space is that
* a discrete functional space carries some additional information which is needed to actually construct a finite element space (e.g. polynomial order, and rank).

44
src/include/serif/discretization/quadrature/backend/mfem/concepts.hpp
View File

@@ -1,44 +0,0 @@
#pragma once
/**
* @file concepts.hpp
* @author Emily M. Boudreaux
*
* @brief File containing concepts for the quadrature backend module. These concepts are used to constrain template parameters in the quadrature backend module.
*/
#include "mfem/fem/intrules.hpp"
namespace serif::discretization::quadrature::backend::mfem {
/**
* @brief A HasReferenceSetIntegrationRule is an integrator that can be configured with a reference to an integration rule.
*/
template <typename IntegratorT>
concept HasReferenceSetIntegrationRule = requires (IntegratorT &integrator, const ::mfem::IntegrationRule &integration_rule) {
integrator.SetIntegrationRule(integration_rule);
};
/**
* @brief A HasPointerSetIntegrationRule is an integrator that can be configured with a pointer to an integration rule.
*/
template <typename IntegratorT>
concept HasPointerSetIntegrationRule = requires (IntegratorT &integrator, const ::mfem::IntegrationRule &integration_rule) {
integrator.SetIntRule(&integration_rule);
};
/**
* @brief A HasSetIntRule is an integrator that can be configured with an integration rule by value.
*/
template <typename IntegratorT>
concept HasSetIntRule = requires (IntegratorT &integrator, const ::mfem::IntegrationRule &integration_rule) {
integrator.SetIntRule(integration_rule);
};
/**
* @brief A ConfigurableIntegrator is an integrator that can be configured with an integration rule in some way (either by reference, by pointer, or by value).
*/
template <typename IntegratorT>
concept ConfigurableIntegrator = HasReferenceSetIntegrationRule<IntegratorT> ||
HasPointerSetIntegrationRule<IntegratorT> ||
HasSetIntRule<IntegratorT>;
}

32
src/include/serif/discretization/quadrature/backend/mfem/configure.hpp
View File

@@ -1,32 +0,0 @@
#pragma once
/**
* @file configure.hpp
* @author Emily M. Boudreaux
*
* @brief File containing utility functions for configuring MFEM integrators with integration rules.
*/
#include "mfem/fem/intrules.hpp"
#include "serif/discretization/quadrature/backend/mfem/concepts.hpp"
namespace serif::discretization::quadrature::backend::mfem {
/**
* @brief Configure an MFEM integrator with an integration rule. This function will call the appropriate SetIntegrationRule method on the integrator based on its capabilities.
* @tparam IntegratorT The type of the integrator to configure.
* @param integrator The integrator to configure.
* @param integration_rule The integration rule to use for configuration.
*/
template <ConfigurableIntegrator IntegratorT>
void set_integration_rule(IntegratorT &integrator, const ::mfem::IntegrationRule &integration_rule) {
if constexpr (HasReferenceSetIntegrationRule<IntegratorT>) {
integrator.SetIntegrationRule(integration_rule);
} else if constexpr (HasPointerSetIntegrationRule<IntegratorT>) {
integrator.SetIntRule(&integration_rule);
} else {
integrator.SetIntRule(integration_rule);
}
}
}

166
src/include/serif/discretization/quadrature/backend/mfem/mfem_resolver.hpp
View File

@@ -1,166 +0,0 @@
#pragma once
/**
* @file mfem_resolver.hpp
* @author Emily M. Boudreaux
*
* @brief File containing the MFEMResolver class. This class is used to resolve the appropriate MFEM integration rule for a given quadrature context.
*
* @note We may want to consider renaming this file to runtime to stick to convention. The reason I did not initially is that this is more specific than pure runtime
* projection of compile time invariants (which is what most of the other files named runtime.hpp do). Something to think about.
*
* @note You will see ::mfem throughout. This is an important syntax to use here since MFEM's code lives in the mfem namespace; however, we also have code living in an
* mfem namespace. Outside of translation units that include using namespace serif::discretization::quadrature::backend and outside of the
* serif::discretization::quadrature::backend::mfem namespace itself these two things do not conflict. However, in those two cases there is a degeneracy the compiler
* cannot resolve between mfem:: and mfem::. Using ::mfem explicitly tells the compiler that we are looking from the root of the namespace tree not the current namespace.
*/
#include <span>
#include <format>
#include <unordered_map>
#include "mfem/fem/geom.hpp"
#include "mfem/fem/intrules.hpp"
#include "serif/discretization/quadrature/rules.hpp"
#include "serif/discretization/quadrature/backend/order.hpp"
#include "serif/discretization/domain/concepts.hpp"
#include "serif/discretization/quadrature/policy.hpp"
#include "serif/discretization/quadrature/query.hpp"
#include "serif/discretization/quadrature/backend/mfem/concepts.hpp"
#include "serif/discretization/quadrature/backend/mfem/configure.hpp"
#include "serif/discretization/quadrature/backend/mfem/mfem_static_order.hpp"
#include "serif/utils/format/strings.hpp"
namespace serif::discretization::quadrature::backend::mfem {
/**
* @brief The MFEMQuadratureContext struct is used to describe the context in which we are trying to resolve an MFEM integration rule. This includes the geometry type, the order of the geometry weight, any dynamic orders, the domain or set of domains, and the mapping kind.
*
* The idea here is that there is some information we can know at compile time (e.g. the role the particular integration is playing) and some we can only know
* at runtime (e.g. a manual quadrature boost requested by the user) so QuadratureContext contains all of that information which can then be turned into
* a single quadrature order. This is born out of the recognition that there is really not one good choice for quadrature order; rather,
* different systems have different accuracy requirements and should generally use the lowest order that is sufficient for their needs.
*/
struct MFEMQuadratureContext final {
QuadratureRole role = QuadratureRole::discretization; ///< The role of the quadrature rule.
::mfem::Geometry::Type geometry = ::mfem::Geometry::INVALID; ///< The geometry type of the element for which we are trying to resolve an MFEM integration rule.
int geometry_weight_order = 0; ///< The order of the geometry weight.
std::span<const int> dynamic_orders{}; ///< The dynamic orders of the quadrature rule. Dynamics orders come from runtime information.
domain::DomainOrSet domain = domain::AllDomains{}; ///< The domain or set of domains for which we are trying to resolve an MFEM integration rule.
MappingKind mapping = MappingKind::none; ///< The mapping kind of the quadrature rule.
};
/**
* @brief The MFEMRule struct is used to describe the result of resolving an MFEM integration rule. This includes the resolution (which contains the order and any other relevant information) and a pointer to the integration rule itself.
*/
struct MFEMRule final {
Resolution resolution{}; ///< The resolution of the MFEM integration rule, which contains the order and any other relevant information.
const ::mfem::IntegrationRule *integration_rule = nullptr; ///< A pointer to the MFEM integration rule itself.
};
/**
* @brief The MFEMResolver class is used to resolve the appropriate MFEM integration rule for a given quadrature context. It uses a policy to determine the appropriate order and other relevant information for the integration rule.
*/
class MFEMResolver final {
public:
/**
* @brief Constructs an MFEMResolver with the given policy.
*
* @param policy The policy used to determine the appropriate order and other relevant information for the integration rule.
*/
explicit MFEMResolver (Policy policy) noexcept : m_policy(std::move(policy)) {} // Do we actually want to use move rvalue semantics here?
/**
* @brief Resolves the appropriate MFEM integration rule for a given quadrature context and form type.
*
* @tparam FormT The form type for which we are trying to resolve an MFEM integration rule.
*
* @param context The quadrature context describing the situation in which we are trying to resolve an MFEM integration rule.
*
* @return An MFEMRule struct containing the resolution and a pointer to the MFEM integration rule itself.
*/
template <forms::Form FormT>
[[nodiscard]] MFEMRule resolve(const MFEMQuadratureContext& context) const {
if (context.geometry_weight_order < 0) {
utils::errors::serif_error<std::invalid_argument, utils::errors::SERiFErrorCode::INVALID_QUADRATURE_ORDER>("MFEM geometry weight order must be non-negative. Received an order of {}", context.geometry_weight_order);
}
const int dynamic_order = sum_dynamic_orders(context.dynamic_orders);
const int base_order = mfem_form_static_order<FormT> + dynamic_order + context.geometry_weight_order;
const Query query {
.form = forms::form_descriptor<FormT>,
.role = context.role,
.domain = context.domain,
.mapping = context.mapping,
.base_order = base_order
};
const Resolution resolution = m_policy.resolve(query);
const ::mfem::IntegrationRule &integration_rule = ::mfem::IntRules.Get(context.geometry, resolution.order);
return MFEMRule{.resolution = resolution, .integration_rule = &integration_rule};
}
/**
* @brief Configures an MFEM integrator with the appropriate integration rule for a given quadrature context and form type. This function will call the appropriate SetIntegrationRule method on the integrator based on its capabilities.
*
* @tparam FormT The form type for which we are trying to configure the integrator.
* @tparam IntegratorT The type of the integrator to configure.
*
* @param integrator The integrator to configure.
* @param context The quadrature context describing the situation in which we are trying to configure the integrator.
*
* @return The resolution of the MFEM integration rule used to configure the integrator.
*/
template <forms::Form FormT, ConfigurableIntegrator IntegratorT>
Resolution configure(IntegratorT &integrator, const MFEMQuadratureContext& context) const {
const MFEMRule rule = resolve<FormT>(context);
set_integration_rule(integrator, *rule.integration_rule);
return rule.resolution;
}
private:
Policy m_policy; ///< The policy used to determine the appropriate order and other relevant information for the integration rule.
};
}
/**
* @brief Specialization of std::formatter for MFEMQuadratureContext. This allows for easy formatting of MFEMQuadratureContext objects using std::format.
*/
template <>
struct std::formatter<serif::discretization::quadrature::backend::mfem::MFEMQuadratureContext> : std::formatter<std::string> {
template <typename FormatContext>
auto format(const serif::discretization::quadrature::backend::mfem::MFEMQuadratureContext& mfem_context, FormatContext& fc) {
std::unordered_map<::mfem::Geometry::Type, std::string> MFEM_GEOM_TYPE_NAME_MAP {
{::mfem::Geometry::INVALID, "INVALID"},
{::mfem::Geometry::POINT, "POINT"},
{::mfem::Geometry::SEGMENT, "SEGMENT"},
{::mfem::Geometry::SQUARE, "SQUARE"},
{::mfem::Geometry::TETRAHEDRON, "TETRAHEDRON"},
{::mfem::Geometry::CUBE, "CUBE"},
{::mfem::Geometry::PRISM, "PRISM"},
{::mfem::Geometry::PYRAMID, "PYRAMID"},
{::mfem::Geometry::NUM_GEOMETRIES, std::format("NUM_GEOMETRIES={}", static_cast<int>(mfem::Geometry::NUM_GEOMETRIES))}
};
return std::formatter<std::string>::format(
std::format(
"MFEMQuadratureContext(form={}, geom_type={}, geom_weight_order={}, dynamic_orders={}, domains={}, mapping={})",
mfem_context.role,
MFEM_GEOM_TYPE_NAME_MAP.at(mfem_context.geometry),
mfem_context.geometry_weight_order,
serif::utils::misc::format::join(mfem_context.dynamic_orders, ", "),
mfem_context.domain,
mfem_context.mapping
), fc);
}
};

39
src/include/serif/discretization/quadrature/backend/mfem/mfem_static_order.hpp
View File

@@ -1,39 +0,0 @@
#pragma once
/**
* @file mfem_static_order.hpp
* @author Emily M. Boudreaux
*
* @brief File containing the mfem_form_static_order trait. This trait is used to determine the static order of an MFEM integrator for a given form type.
*/
#include "serif/discretization/blocks/concepts.hpp"
#include "serif/discretization/forms/operands.hpp"
#include "serif/discretization/forms/concept.hpp"
#include "serif/utils/types/type_list.hpp"
#include "serif/discretization/quadrature/backend/order.hpp"
namespace serif::discretization::quadrature::backend::mfem {
/**
* @brief A helper struct to compute the static order of an MFEM integrator for a given list of operands.
* @tparam OperandListT The list of operand types.
*/
template <typename OperandListT>
struct MFEMOperandListOrder;
/**
* @brief A specialization of the MFEMOperandListOrder struct for a list of operand types. This struct is used to compute the static order of an MFEM integrator for a given list of operands.
* @tparam OperandTs The operand types that make up the list.
*/
template <forms::IsOperand... OperandTs>
struct MFEMOperandListOrder<utils::types::TypeList<OperandTs...>> final {
static constexpr int value = (resolve_operand_order<OperandTs>() + ... + 0); ///< The static order of the MFEM integrator for the given list of operands. This is computed by summing the static orders of each operand in the list.
};
/**
* @brief A trait to determine the static order of an MFEM integrator for a given form type.
* @tparam FormT The form type to check.
*/
template <forms::Form FormT>
inline constexpr int mfem_form_static_order = MFEMOperandListOrder<typename FormT::operand_types>::value;
}

124
src/include/serif/discretization/quadrature/backend/order.hpp
View File

@@ -1,124 +0,0 @@
#pragma once
/**
* @file order.hpp
* @author Emily M. Boudreaux
*
* @brief File containing the resolve_operand_order function and related utilities. These are used to determine the static order of an integrator for a given operand or form type.
*/
#include <type_traits>
#include "serif/discretization/functional_space/space.hpp"
#include "serif/discretization/blocks/concepts.hpp"
#include "serif/discretization/forms/operands.hpp"
#include "serif/utils/error/errors.hpp"
namespace serif::discretization::quadrature::backend::mfem {
/**
* @brief A helper struct to extract the functional space traits from a DiscreteFunctionalSpace.
* @tparam FunctionalSpaceT The functional space type to extract traits from.
*/
template <typename FunctionalSpaceT>
struct DiscreteFunctionalSpaceTraits;
/**
* @brief A specialization of the DiscreteFunctionalSpaceTraits struct for a DiscreteFunctionalSpace. This struct is used to extract the functional space traits from a DiscreteFunctionalSpace.
*
* @tparam SpaceV The functional space family.
* @tparam OrderV The order of the functional space.
* @tparam RankV The rank of the functional space.
*
* @todo Could we remove this?
*/
template <auto SpaceV, auto OrderV, auto RankV>
struct DiscreteFunctionalSpaceTraits<functional_space::DiscreteFunctionalSpace<SpaceV, OrderV, RankV>> final {
static constexpr auto space = SpaceV; ///< The functional space family.
static constexpr auto order = static_cast<int>(OrderV); ///< The order of the functional space.
static constexpr std::size_t rank = static_cast<std::size_t>(RankV); ///< The rank of the functional space.
};
/**
* @brief A helper alias to extract the functional space traits from a DiscreteFunctionalSpace.
* @tparam FunctionalSpaceT The functional space type to extract traits from.
*/
template <blocks::HasFunctionalSpace TermT>
using FunctionalSpaceTraits = DiscreteFunctionalSpaceTraits<typename TermT::FunctionalSpace>;
/**
* @brief A function to resolve the static order of an MFEM integrator for a given operand type.
*
* @tparam OperandT The operand type to resolve the static order for.
*
* @return The static order of the MFEM integrator for the given operand type.
*
* @todo We already have some of this information in the form_traits file, we should probably extract the order boost from there.
*/
template <forms::IsOperand OperandT>
[[nodiscard]] consteval int resolve_operand_order() noexcept {
using Operand = std::remove_cvref_t<OperandT>;
using Term = Operand::term_type;
using Operation = Operand::operation_type;
if constexpr (!blocks::HasFunctionalSpace<Term>) {
static_assert(std::same_as<Operation, forms::Operation::Value>, "Non-Spatial terms only support the value operation at the moment.");
return 0;
} else {
using Space = FunctionalSpaceTraits<Term>;
constexpr auto family = Space::space.family;
constexpr int order = Space::order;
if constexpr (std::same_as<Operation, forms::Operation::Value>) {
if constexpr (family == functional_space::Family::hdiv) {
return order + 1;
} else {
return order;
}
} else if constexpr (std::same_as<Operation, forms::Operation::Gradient>) {
static_assert(family == functional_space::Family::h1, "Only H1 spaces support the gradient operation at the moment.");
return order > 0 ? order - 1 : 0;
} else if constexpr (std::same_as<Operation, forms::Operation::Divergence>) {
static_assert(Space::rank == 1 && (family == functional_space::Family::hdiv || family == functional_space::Family::h1), "Divergence requires a vector H1 or H(div) term.");
if constexpr (family == functional_space::Family::hdiv) {
return order;
} else {
return order > 0 ? order - 1 : 0;
}
} else if constexpr (std::same_as<Operation, forms::Operation::Curl>) {
static_assert(family == functional_space::Family::hcurl, "Only H(curl) spaces support the curl operation at the moment.");
return order > 0 ? order - 1 : 0;
} else if constexpr (std::same_as<Operation, forms::Operation::NormalTrace>) {
static_assert(family == functional_space::Family::hdiv, "Only H(div) spaces support the normal trace operation at the moment.");
return order;
} else {
static_assert(false, "Unsupported operation type. This indicates that an operation was used which is not currently registered in the operand configuration system. If you are a developer and you are seeing this error, please register your operation in the resolve_operand_order function in discretization/quadrature/detail.hpp. If you are a user and you are seeing this error, please report it to the SERiF developers.");
return 0;
}
}
}
/**
* @brief Sums the dynamic orders in a span of integers.
* @param dynamic_orders The span of dynamic orders to sum.
* @return The sum of the dynamic orders.
*/
[[nodiscard]] inline int sum_dynamic_orders(const std::span<const int> dynamic_orders) {
int total = 0;
for (const int order : dynamic_orders) {
if (order < 0) {
utils::errors::serif_error<std::invalid_argument, utils::errors::SERiFErrorCode::INVALID_QUADRATURE_ORDER>("Dynamic quadrature order cannot be negative. Received {}.", order);
}
total += order;
}
return total;
}
}

130
src/include/serif/discretization/quadrature/rules.hpp
View File

@@ -43,19 +43,6 @@ namespace serif::discretization::quadrature {
kelvin ///< The quadrature rule has a Kelvin mapping (compactified).
};
/**
* @brief An enumeration of the different modes that the quadrature system can operate in.
*
* The idea here is that different modes might require or allow different levels of accuracy. For example, a fast mode might prioritize speed over accuracy,
* while a production mode might prioritize accuracy over speed.
*/
enum class Mode {
fast, ///< The quadrature system is operating in fast mode.
production, ///< The quadrature system is operating in production mode.
reference, ///< The quadrature system is operating in reference mode.
convergence ///< The quadrature system is operating in convergence mode. (note in convergence mode the system may inject additional loops to test convergence of the quadrature rule). Convergence should only ever be used for testing and validation, not in production.
};
/**
* @brief A struct to control the behavior of a quadrature rule for a given form or role. Rule control always
* contains boost (default to 0) and may contain some fixed order (think of this as an override). If the fixed order is empty
@@ -248,16 +235,6 @@ namespace serif::discretization::quadrature {
RuleControl m_fallbackControl; ///< The fallback rule control for the RuleSet. This is used when no specific form or role control is found.
};
/**
* @brief A struct to hold the rule controls for each quadrature role. This is used to group the rule controls together for convenience.
*/
struct RoleControls {
RuleControl discretization; ///< The rule control for the discretization role.
RuleControl preconditioner; ///< The rule control for the preconditioner role.
RuleControl diagnostic; ///< The rule control for the diagnostic role.
RuleControl projection; ///< The rule control for the projection role.
};
/**
* @brief A struct to hold the resolution of a quadrature query. This contains the base order, boost, and final order of the quadrature rule, as well as a flag indicating whether a fixed order was used.
*/
@@ -268,83 +245,10 @@ namespace serif::discretization::quadrature {
bool used_fixed_order; ///< A flag indicating whether a fixed order was used.
};
/**
* @brief A struct to hold the options for a quadrature term. This contains an optional fixed order and an additional boost to the order of the quadrature rule.
*/
struct QuadratureTermOptions {
std::optional<int> fixed_order; ///< An optional fixed order for the quadrature rule. If this is set, it will override the default order determined by the form and role.
int additional_boost = 0; ///< An additional boost to the order of the quadrature rule. This is added to the default order determined by the form and role.
};
/**
* @brief A struct to hold the options for the quadrature manifest. This contains flags for enabling the manifest, including repeated queries, and an optional output file.
*/
struct QuadratureManifestOptions {
bool enabled = false; ///< A flag indicating whether the quadrature manifest is enabled.
bool include_repeated_queries = false; ///< A flag indicating whether to include repeated queries in the quadrature manifest.
std::optional<std::string> output_file; ///< An optional output file for the quadrature manifest. If this is set, the manifest will be written to the specified file.
};
/**
* @brief A struct to hold the options for quadrature validation. This contains flags for requiring explicit base orders and MFEM rules, rejecting negative boosts, and reporting unused overrides.
*/
struct QuadratureValidationOptions {
bool require_explicit_base_order = false; ///< A flag indicating whether explicit base orders are required.
bool require_explicit_mfem_rule = false; ///< A flag indicating whether explicit MFEM rules are required.
bool reject_negative_boosts = true; ///< A flag indicating whether negative boosts should be rejected.
bool report_unused_overrides = true; ///< A flag indicating whether unused overrides should be reported.
};
/**
* @brief A struct to hold the options for quadrature roles. This contains the quadrature term options for each role.
*/
struct QuadratureRoleOptions {
QuadratureTermOptions discretization; ///< The quadrature term options for the discretization role.
QuadratureTermOptions preconditioner; ///< The quadrature term options for the preconditioner role.
QuadratureTermOptions diagnostic; ///< The quadrature term options for the diagnostic role.
QuadratureTermOptions projection; ///< The quadrature term options for the projection role.
};
/**
* @brief A struct to hold the options for the quadrature system. This contains the mode, global boost, fallback fixed order, quadrature term options for each term, quadrature role options, convergence boosts, quadrature manifest options, and quadrature validation options.
*/
struct QuadratureOptions {
Mode mode = Mode::production; ///< The mode of the quadrature system. This determines the level of accuracy and speed of the quadrature system.
int global_boost = 0; ///< The global boost to the order of the quadrature rule. This is added to the default order determined by the form and role. Note that different systems have different accuracy requirements and should generally use the lowest order that is sufficient for their needs.
std::optional<int> fallback_fixed_order; ///< The fallback fixed order for the quadrature rule. This is used if no other order is specified.
QuadratureTermOptions gravity_hdiv_mass; ///< The quadrature term options for the gravity H(div) mass term.
QuadratureTermOptions gravity_divergence; ///< The quadrature term options for the gravity divergence term.
QuadratureTermOptions gravity_source; ///< The quadrature term options for the gravity source term.
QuadratureTermOptions gravity_force; ///< The quadrature term options for the gravity force term.
QuadratureTermOptions gravity_boundary; ///< The quadrature term options for the gravity boundary term.
QuadratureTermOptions centrifugal; ///< The quadrature term options for the centrifugal term.
QuadratureTermOptions density_projection; ///< The quadrature term options for the density projection term.
QuadratureTermOptions eos_closure; ///< The quadrature term options for the EOS closure term.
QuadratureTermOptions hydrostatic_equilibrium; ///< The quadrature term options for the hydrostatic equilibrium term.
QuadratureTermOptions isobaric_surface; ///< The quadrature term options for the isobaric surface term.
QuadratureTermOptions mesh_extension; ///< The quadrature term options for the mesh extension term.
QuadratureTermOptions mass_conservation; ///< The quadrature term options for the mass conservation term.
QuadratureTermOptions mass_normalization; ///< The quadrature term options for the mass normalization term.
QuadratureTermOptions center_of_mass; ///< The quadrature term options for the center of mass term.
QuadratureTermOptions quadrupole; ///< The quadrature term options for the quadrupole term.
QuadratureTermOptions gravitational_energy; ///< The quadrature term options for the gravitational energy term.
QuadratureTermOptions pressure_integral; ///< The quadrature term options for the pressure integral term.
QuadratureTermOptions pressure_force; ///< The quadrature term options for the pressure force term.
QuadratureTermOptions virial; ///< The quadrature term options for the virial term.
QuadratureTermOptions error_norm; ///< The quadrature term options for the error norm term.
QuadratureRoleOptions roles; ///< The quadrature role options for each quadrature role.
std::vector<int> convergence_boosts = {0, 2, 4}; ///< The set of convergence boosts used for a tests (only used when the quadrature system is in convergence mode).
QuadratureManifestOptions manifest; ///< The quadrature manifest options for the quadrature system. This contains flags for enabling the manifest, including repeated queries, and an optional output file.
QuadratureValidationOptions validation; ///< The quadrature validation options for the quadrature system. This contains flags for requiring explicit base orders and MFEM rules, rejecting negative boosts, and reporting unused overrides.
};
}
/**
* @brief Specializations of the std::formatter for the QuadratureRole, MappingKind, Mode, RuleControl, FormRule, and RuleSet types. This allows for easy formatting of these types using std::format.
* @brief Specializations of the std::formatter for the QuadratureRole, MappingKind, RuleControl, FormRule, and RuleSet types. This allows for easy formatting of these types using std::format.
*/
template <>
struct std::formatter<serif::discretization::quadrature::QuadratureRole> : std::formatter<std::string_view> {
@@ -352,6 +256,9 @@ struct std::formatter<serif::discretization::quadrature::QuadratureRole> : std::
auto format(const serif::discretization::quadrature::QuadratureRole& role, FormatContext& ctx) const {
std::string_view roleStr;
switch (role) {
case serif::discretization::quadrature::QuadratureRole::discretization:
roleStr = "discretization";
break;
case serif::discretization::quadrature::QuadratureRole::preconditioner:
roleStr = "preconditioner";
break;
@@ -398,35 +305,6 @@ struct std::formatter<serif::discretization::quadrature::MappingKind> : std::for
}
};
/**
* @brief A specialization of the std::formatter for the Mode type. This allows for easy formatting of the Mode type using std::format.
*/
template <>
struct std::formatter<serif::discretization::quadrature::Mode> : std::formatter<std::string_view> {
template <typename FormatContext>
auto format(const serif::discretization::quadrature::Mode& mode, FormatContext& ctx) const {
std::string_view modeStr;
switch (mode) {
case serif::discretization::quadrature::Mode::fast:
modeStr = "fast";
break;
case serif::discretization::quadrature::Mode::production:
modeStr = "production";
break;
case serif::discretization::quadrature::Mode::reference:
modeStr = "reference";
break;
case serif::discretization::quadrature::Mode::convergence:
modeStr = "convergence";
break;
default:
modeStr = "unknown";
break;
}
return std::formatter<std::string_view>::format(modeStr, ctx);
}
};
/**
* @brief A specialization of the std::formatter for the RuleControl type. This allows for easy formatting of the RuleControl type using std::format.
*/

104
src/include/serif/discretization/reconstruction/surface_power.hpp Normal file
View File

@@ -0,0 +1,104 @@
#pragma once
#include <cmath>
#include <span>
#include <string_view>
#include "serif/discretization/domain/physical_domains.hpp"
#include "serif/utils/error/errors.hpp"
namespace serif::discretization::reconstruction {
/** @brief Name the boundary whose material-side elements use reconstruction. */
struct SurfaceRegion {
std::string_view boundary_name;
};
/** @brief Select material elements touching a typed boundary, resolved through the mesh schema. */
template <typename BoundaryT = domain::StellarSurfaceBoundary>
struct SurfaceAdjacent : SurfaceRegion {
/** @brief Retain the scientific boundary name without exposing a mesh attribute. */
constexpr SurfaceAdjacent() : SurfaceRegion{BoundaryT::name} {}
};
/**
* @brief Request density proportional to inward distance raised to density_exponent near a vacuum surface.
* @note Coefficients remain physical nodal densities. For a Qp element with zero surface trace,
* the reconstructed space is t^n P_(p-1) radially and ordinary Qp tangentially; t runs from
* zero on the surface face to one on the inner face. Deeper elements are unchanged.
* The exponent is n (not n-1) for a vacuum polytrope. A finite-pressure or general EOS
* surface needs a separately validated profile; this request does not infer it from an EOS.
*/
template <typename TermT>
struct SurfacePower {
using Term = TermT;
double density_exponent{1.0};
SurfaceRegion region{SurfaceAdjacent<>{}};
};
/**
* @brief One-dimensional nodal reconstruction, evaluated only while preparing reference operators.
* @note A new form supplies nodal values and inward derivatives, including both endpoints.
* The complete basis must reproduce constants and have compatible traces. Its archive codec
* and admissibility certificate must be added explicitly before it can be used in a saved solve.
*/
class RadialProfile {
public:
/** @brief Release a profile through its immutable shared owner. */
virtual ~RadialProfile() = default;
/** @brief Identify the numerical form in diagnostics and archive validation. */
[[nodiscard]] virtual std::string_view name() const noexcept = 0;
/** @brief Evaluate every radial nodal function and its inward-coordinate derivative. */
virtual void evaluate(std::span<const double> nodes, double inward_coordinate,
std::span<double> values, std::span<double> derivatives) const = 0;
};
/** @brief Complete nodal basis whose homogeneous surface trace spans t^n times degree-(p-1) polynomials. */
class SurfacePowerProfile final : public RadialProfile {
public:
/** @brief Require a finite exponent with a finite density gradient at the surface. */
explicit SurfacePowerProfile(double density_exponent) : m_exponent(density_exponent) {
using namespace utils::errors;
serif_ensure<SERiFErrorCode::INVALID_ARGUMENT>(std::isfinite(m_exponent) && m_exponent >= 1.0,
"SurfacePower requires a finite density exponent >= 1; received {}. Choose the polytropic index of a vacuum model in the supported range.", m_exponent);
}
/** @brief Return the stable numerical-form name. */
[[nodiscard]] std::string_view name() const noexcept override { return "surface_power"; }
/** @brief Return the density's leading surface exponent, rather than the weight's exponent. */
[[nodiscard]] double density_exponent() const noexcept { return m_exponent; }
/** @brief Evaluate the complete lifting rho_s + sum_i psi_i (rho_i-rho_s), including analytic endpoint limits. */
void evaluate(std::span<const double> nodes, double inward_coordinate,
std::span<double> values, std::span<double> derivatives) const override {
using namespace utils::errors;
serif_ensure<SERiFErrorCode::INVALID_ARGUMENT>(nodes.size() >= 2 && values.size() == nodes.size() && derivatives.size() == nodes.size() &&
nodes.front() == 0.0 && nodes.back() == 1.0 && std::isfinite(inward_coordinate) && inward_coordinate >= 0.0 && inward_coordinate <= 1.0,
"SurfacePower evaluation requires closed inward nodes on [0,1] and matching output sizes; received coordinate {} and {} nodes.", inward_coordinate, nodes.size());
for (std::size_t node = 1; node < nodes.size(); ++node) {
serif_ensure<SERiFErrorCode::INVALID_ARGUMENT>(std::isfinite(nodes[node]) && nodes[node] > nodes[node - 1],
"SurfacePower requires strictly increasing radial nodes.");
}
values[0] = 1.0;
derivatives[0] = 0.0;
for (std::size_t node = 1; node < nodes.size(); ++node) {
double cardinal = 1.0;
double cardinal_derivative = 0.0;
for (std::size_t other = 1; other < nodes.size(); ++other) {
if (other == node) { continue; }
const double denominator = nodes[node] - nodes[other];
cardinal_derivative = cardinal_derivative * (inward_coordinate - nodes[other]) / denominator + cardinal / denominator;
cardinal *= (inward_coordinate - nodes[other]) / denominator;
}
const double ratio = inward_coordinate / nodes[node];
const double power = std::pow(ratio, m_exponent);
const double power_derivative = m_exponent / nodes[node] * std::pow(ratio, m_exponent - 1.0);
values[node] = power * cardinal;
derivatives[node] = power_derivative * cardinal + power * cardinal_derivative;
serif_ensure<SERiFErrorCode::INVALID_ARGUMENT>(std::isfinite(values[node]) && std::isfinite(derivatives[node]),
"SurfacePower exponent {} cannot be evaluated finitely on these radial nodes.", m_exponent);
values[0] -= values[node];
derivatives[0] -= derivatives[node];
}
}
private:
double m_exponent;
};
}

13
src/include/serif/eos/evaluation.hpp
View File

@@ -59,4 +59,17 @@ namespace serif::eos {
return eos.partial_derivative(relations::Relation<OutputQ, dimensions::QuantityOfT<InputValues>...>{}, dimensions::WithRespectTo<InputQ>{}, inputValues...);
}
/** @brief Dispatch an analytic second partial without imposing an EOS-specific state. */
template <dimensions::ThermodynamicQuantityType OutputQ, dimensions::ThermodynamicQuantityType FirstInput,
dimensions::ThermodynamicQuantityType SecondInput, models::EOSModel EOS, dimensions::QuantityValueType... InputValues>
requires models::EOSSupportingSecondPartialDerivative<EOS,
relations::Relation<OutputQ, dimensions::QuantityOfT<InputValues>...>, FirstInput, SecondInput>
[[nodiscard]] constexpr dimensions::SecondPartialDerivative<OutputQ, FirstInput, SecondInput>
second_partial_derivative(const EOS& eos, const InputValues... inputValues)
noexcept(noexcept(eos.second_partial_derivative(relations::Relation<OutputQ, dimensions::QuantityOfT<InputValues>...>{},
dimensions::WithRespectTo<FirstInput>{}, dimensions::WithRespectTo<SecondInput>{}, inputValues...))) {
return eos.second_partial_derivative(relations::Relation<OutputQ, dimensions::QuantityOfT<InputValues>...>{},
dimensions::WithRespectTo<FirstInput>{}, dimensions::WithRespectTo<SecondInput>{}, inputValues...);
}
}

23
src/include/serif/eos/models/concepts.hpp
View File

@@ -29,6 +29,11 @@ namespace serif::eos::models {
template <typename EOS, typename RelationT, typename InputQuantity>
concept EOSSupportingPartialDerivative = supports_partial_derivative_v<EOS, RelationT, InputQuantity>;
/** @brief Require a declared EOS relation and its analytic second partial. */
template <typename EOS, typename RelationT, typename FirstInput, typename SecondInput>
concept EOSSupportingSecondPartialDerivative = EOSSupportingRelation<EOS, RelationT> &&
ImplementsSecondPartialDerivative_t<std::remove_cvref_t<EOS>, std::remove_cvref_t<RelationT>, FirstInput, SecondInput>::value;
/**
* @brief A concept to determine if a type is an equation of state model that can be used in the current SERiF seed initialization system. This is satisfied for types that are both valid EOSModels and
* which also support the relation for computing specific enthalpy from density.
@@ -38,22 +43,4 @@ namespace serif::eos::models {
template <typename EOSCandidate>
concept StructureSeedEOS = EOSModel<EOSCandidate> &&
EOSSupportingRelation<EOSCandidate, relations::SpecificEnthalpyFromDensity>;
/**
* @brief A concept to determine if an EOS type can be used as a barotropic closure in the current SERiF system. This is satisfied for types that are
* valid EOS models, support the relation for computing density from specific enthalpy, and support the partial derivative of that relation with respect to specific enthalpy.
*/
template <typename EOSCandidate>
concept BarotropicClosureEOS = EOSModel<EOSCandidate> &&
EOSSupportingRelation<EOSCandidate, relations::DensityFromSpecificEnthalpy> &&
EOSSupportingPartialDerivative<EOSCandidate, relations::DensityFromSpecificEnthalpy, dimensions::SpecificEnthalpy>;
/**
* @brief A concept to determine if an EOS type can be used to compute pressure forces in the current SERiF system. This is satisfied for types that are valid EOS models, support
* the relation for computing pressure from specific enthalpy, and support the partial derivative of that relation with respect to specific enthalpy.
*/
template <typename EOSCandidate>
concept PressureForceEOS = EOSModel<EOSCandidate> &&
EOSSupportingRelation<EOSCandidate, relations::PressureFromSpecificEnthalpy> &&
EOSSupportingPartialDerivative<EOSCandidate, relations::PressureFromSpecificEnthalpy, dimensions::SpecificEnthalpy>;
}

19
src/include/serif/eos/models/polytropic.hpp
View File

@@ -41,17 +41,21 @@ namespace serif::eos::models {
*
* A polytropic EOS supports the following relations:
* - PressureFromDensity
* - DensityFromPressure
* - PressureFromSpecificEnthalpy
* - SpecificEnthalpyFromDensity
* - SpecificEnthalpyFromPressure
* - DensityFromSpecificEnthalpy
* - PressurePotentialFromDensity
*/
using Relations = relations::RelationCatalog<
relations::PressureFromDensity,
relations::DensityFromPressure,
relations::PressureFromSpecificEnthalpy,
relations::SpecificEnthalpyFromDensity,
relations::SpecificEnthalpyFromPressure,
relations::DensityFromSpecificEnthalpy
relations::DensityFromSpecificEnthalpy,
relations::PressurePotentialFromDensity
>;
/**
@@ -96,6 +100,13 @@ namespace serif::eos::models {
*/
[[nodiscard]] dimensions::PressureValue evaluate(relations::PressureFromDensity /*relation*/, dimensions::DensityValue density) const;
/** @brief Evaluate the density at a nonnegative pressure, including the vacuum trace. */
[[nodiscard]] dimensions::DensityValue evaluate(relations::DensityFromPressure /*relation*/, dimensions::PressureValue pressure) const;
/** @brief Evaluate the pressure primitive; for this EOS it equals the specific enthalpy and vanishes in vacuum. */
[[nodiscard]] dimensions::BarotropicPressurePotentialValue evaluate(
relations::PressurePotentialFromDensity /*relation*/, dimensions::DensityValue density) const;
/**
* @brief Evaluate the pressure from the specific enthalpy.
* @param relation The relation to evaluate (note that this is not used in the function but is used to constrain the overload being called).
@@ -160,6 +171,12 @@ namespace serif::eos::models {
dimensions::DensityValue density
) const;
/** @brief Evaluate d²P/dρ²; for n > 1 its vacuum limit is singular and is rejected. */
[[nodiscard]] dimensions::SecondPartialDerivative<dimensions::Pressure, dimensions::MassDensity, dimensions::MassDensity>
second_partial_derivative(relations::PressureFromDensity /*relation*/,
dimensions::WithRespectTo<dimensions::MassDensity> /*first_input*/,
dimensions::WithRespectTo<dimensions::MassDensity> /*second_input*/, dimensions::DensityValue density) const;
private:
double m_polytropic_index; ///< The polytropic index of the equation of state.
double m_polytropic_constant; ///< The polytropic constant of the equation of state.

15
src/include/serif/eos/models/traits.hpp
View File

@@ -102,6 +102,21 @@ namespace serif::eos::models {
template <typename EOSCandidate>
inline constexpr bool is_equation_of_state_model_v = IsEquationOfStateModel_t<std::remove_cvref_t<EOSCandidate>>::value;
/** @brief A second partial is available only for two actual inputs of the relation. */
template <typename EOS, typename RelationT, typename FirstInput, typename SecondInput>
struct ImplementsSecondPartialDerivative_t : std::false_type {};
/** @brief Require an analytic implementation with the exact physical return type. */
template <typename EOS, typename Output, typename... Inputs, typename FirstInput, typename SecondInput>
struct ImplementsSecondPartialDerivative_t<EOS, relations::Relation<Output, Inputs...>, FirstInput, SecondInput> :
std::bool_constant<
(std::same_as<FirstInput, Inputs> || ...) && (std::same_as<SecondInput, Inputs> || ...) &&
requires(const std::remove_cvref_t<EOS>& eos, dimensions::QuantityValue<Inputs>... inputs) {
{ eos.second_partial_derivative(relations::Relation<Output, Inputs...>{},
dimensions::WithRespectTo<FirstInput>{}, dimensions::WithRespectTo<SecondInput>{}, inputs...) }
-> std::same_as<dimensions::SecondPartialDerivative<Output, FirstInput, SecondInput>>;
}> {};
/**
* @brief A variable template to determine if an EOS type supports a given relation.
* @tparam EOS The EOS type to check.

8
src/include/serif/eos/relations/relations.hpp
View File

@@ -12,6 +12,14 @@
namespace serif::eos::relations {
using PressureFromDensity = Relation<dimensions::Pressure, dimensions::MassDensity>; ///< A relation to compute pressure from density.
using DensityFromPressure = Relation<dimensions::MassDensity, dimensions::Pressure>; ///< A relation to compute density at a prescribed pressure.
/**
* @brief Evaluate A(density), whose derivative on the positive barotropic branch is (dP/density)/density.
*
* An additive constant is arbitrary. A finite surface value is required when a consumer
* subtracts A at a prescribed pressure. This mechanical primitive need not be thermodynamic enthalpy.
*/
using PressurePotentialFromDensity = Relation<dimensions::BarotropicPressurePotential, dimensions::MassDensity>;
using PressureFromSpecificEnthalpy = Relation<dimensions::Pressure, dimensions::SpecificEnthalpy>; ///< A relation to compute pressure from specific enthalpy.
using SpecificEnthalpyFromDensity = Relation<dimensions::SpecificEnthalpy, dimensions::MassDensity>; ///< A relation to compute specific enthalpy from density.
using SpecificEnthalpyFromPressure = Relation<dimensions::SpecificEnthalpy, dimensions::Pressure>; ///< A relation to compute specific enthalpy from pressure.

38
src/include/serif/eos/runtime/concepts.hpp
View File

@@ -1,38 +0,0 @@
#pragma once
/**
* @file concepts.hpp
* @author Emily M. Boudreaux
*
* @brief File containing concepts for the eos runtime module. These concepts are used to constrain template parameters in the eos runtime module.
*/
#include "serif/eos/models/concepts.hpp"
#include "serif/eos/runtime/details.hpp"
namespace serif::eos::runtime {
/**
* @brief A concept to determine if a type is a valid runtime EOS model. A type is a valid RuntimeEOSModel if it satisfies the EOSModel concept and if its Relations catalog is supported by the runtime system.
* @tparam Candidate The type to check.
*/
template <typename Candidate>
concept RuntimeEOSModel = eos::models::EOSModel<Candidate> &&
RuntimeCatalogIsSupported<typename std::remove_cvref_t<Candidate>::Relations>::value;
/**
* @brief A helper alias to adapt a runtime EOS model to the runtime catalog dispatch system. This is used to simplify the usage of runtime EOS models in the runtime system.
* @tparam EOSCandidate The runtime EOS model type to adapt.
*/
template <RuntimeEOSModel EOSCandidate>
using RuntimeEOSAdapter = RuntimeCatalogDispatch<EOSCandidate, typename EOSCandidate::Relations>;
/**
* @brief A helper function to get the runtime relation descriptors for a runtime EOS model. This is used to simplify the usage of runtime EOS models in the runtime system.
* @tparam EOSCandidate The runtime EOS model type to get the descriptors for.
* @return A span of runtime relation descriptors.
*/
template <RuntimeEOSModel EOSCandidate>
[[nodiscard]] constexpr std::span<const dimensions::runtime::RuntimeRelationDescriptor> runtimeRelationDescriptors() noexcept {
return RuntimeCatalogStorage<EOSCandidate, typename EOSCandidate::Relations>::descriptors;
}
}

444
src/include/serif/eos/runtime/details.hpp
View File

@@ -1,444 +0,0 @@
#pragma once
/**
* @file details.hpp
* @author Emily M. Boudreaux
*
* @brief File containing implementation details for the eos runtime module. These details are used internally by the runtime system and are not intended to be used directly by users.
*/
#include <array>
#include <type_traits>
#include <tuple>
#include <expected>
#include "serif/dimensions/runtime/concepts.hpp"
#include "serif/eos/exceptions.hpp"
#include "serif/eos/evaluation.hpp"
#include "serif/eos/relations/catalog.hpp"
#include "serif/utils/error/codes.hpp"
namespace serif::eos::runtime {
/**
* @brief A helper trait to determine if a relation has runtime quantity identifiers. This is used to ensure that all relations in a catalog have unique runtime quantity identifiers.
* @tparam RelationT The relation type to check.
*/
template <typename RelationT>
struct HasRuntimeQuantityIdentifiers : std::false_type {};
/**
* @brief A specialization of the HasRuntimeQuantityIdentifiers trait for a relation. This specialization checks if the output and input quantities of the relation have runtime quantity identifiers.
* @tparam OutputT The output quantity type of the relation.
* @tparam InputTs The input quantity types of the relation.
*/
template <typename OutputT, typename... InputTs>
struct HasRuntimeQuantityIdentifiers<relations::Relation<OutputT, InputTs...>>
: std::bool_constant<
dimensions::runtime::RuntimeIdentifiedThermodynamicQuantity<OutputT> &&
(dimensions::runtime::RuntimeIdentifiedThermodynamicQuantity<InputTs> && ...)> {};
/**
* @brief A helper trait to extract the output and input quantity types from a relation. This is used to simplify the usage of relations in the runtime system.
* @tparam RelationT The relation type to extract the quantity types from.
*/
template <typename RelationT>
struct RuntimeRelationQuantities;
/**
* @brief A specialization of the RuntimeRelationQuantities trait for a relation. This specialization extracts the output and input quantity types from the relation and stores them in a tuple.
* @tparam OutputT The output quantity type of the relation.
* @tparam InputTs The input quantity types of the relation.
*/
template <typename OutputT, typename... InputTs>
struct RuntimeRelationQuantities<relations::Relation<OutputT, InputTs...>> {
using Type = std::tuple<OutputT, InputTs...>;
};
/**
* @brief A helper alias to extract the quantity types from a catalog of relations. This is used to simplify the usage of catalogs in the runtime system.
* @tparam RelationTs The relation types in the catalog.
*/
template <typename... RelationTs>
using RuntimeCatalogQuantityTuple = decltype(std::tuple_cat(std::declval<typename RuntimeRelationQuantities<RelationTs>::Type>()...));
/**
* @brief A helper function to determine if two runtime quantity identifiers are compatible. This is used to ensure that all relations in a catalog have unique runtime quantity identifiers.
* @tparam FirstQuantityT The first quantity type to check.
* @tparam SecondQuantityT The second quantity type to check.
* @return True if the quantity identifiers are compatible, false otherwise.
*/
template <typename FirstQuantityT, typename SecondQuantityT>
[[nodiscard]] consteval bool runtimeQuantityIdentifiersAreCompatible() {
if constexpr (std::same_as<FirstQuantityT, SecondQuantityT>) {
return true;
} else {
return dimensions::runtime::thermodynamicQuantityID<FirstQuantityT> !=
dimensions::runtime::thermodynamicQuantityID<SecondQuantityT>;
}
}
/**
* @brief A helper function to determine if the runtime quantity identifiers in a tuple of quantity types are unambiguous. This is used to ensure that all relations in a catalog have unique runtime quantity identifiers.
* @tparam QuantityTupleT The tuple of quantity types to check.
* @tparam First The index of the first quantity type to check.
* @tparam Offsets The indices of the remaining quantity types to check.
* @return True if the quantity identifiers are unambiguous, false otherwise.
*/
template <typename QuantityTupleT, std::size_t First, std::size_t... Offsets>
[[nodiscard]] consteval bool runtimeQuantityIdentifierIsUnambiguous(std::index_sequence<Offsets...>) {
return (
runtimeQuantityIdentifiersAreCompatible<
std::tuple_element_t<First, QuantityTupleT>,
std::tuple_element_t<First + 1 + Offsets, QuantityTupleT>>() && ...
);
}
/**
* @brief A helper function to determine if the runtime quantity identifiers in a tuple of quantity types are unambiguous. This is used to ensure that all relations in a catalog have unique runtime quantity identifiers.
* @tparam QuantityTupleT The tuple of quantity types to check.
* @tparam Indices The indices of the quantity types to check.
* @return True if the quantity identifiers are unambiguous, false otherwise.
*/
template <typename QuantityTupleT, std::size_t... Indices>
[[nodiscard]] consteval bool runtimeQuantityIdentifiersAreUnambiguous(std::index_sequence<Indices...>) {
return (
runtimeQuantityIdentifierIsUnambiguous<QuantityTupleT, Indices>(std::make_index_sequence<std::tuple_size_v<QuantityTupleT> - Indices - 1>{}) && ...
);
}
/**
* @brief A helper trait to determine if the runtime quantity identifiers in a catalog of relations are unambiguous. This is used to ensure that all relations in a catalog have unique runtime quantity identifiers.
* @tparam RelationTs The relation types in the catalog.
*/
template <bool QuantitiesAreIdentified, typename... RelationTs>
struct RuntimeRelationsAreSupported : std::false_type {};
/**
* @brief A specialization of the RuntimeRelationsAreSupported trait for a catalog of relations. This specialization checks if the runtime quantity identifiers in the catalog are unambiguous.
* @tparam RelationTs The relation types in the catalog.
*/
template <typename... RelationTs>
struct RuntimeRelationsAreSupported<true, RelationTs...> :
std::bool_constant<runtimeQuantityIdentifiersAreUnambiguous<RuntimeCatalogQuantityTuple<RelationTs...>>(
std::make_index_sequence<std::tuple_size_v<RuntimeCatalogQuantityTuple<RelationTs...>>>{}
)> {};
/**
* @brief A helper trait to determine if a catalog of relations is supported by the runtime system. This is used to ensure that all relations in a catalog have unique runtime quantity identifiers.
* @tparam CatalogT The catalog type to check.
*/
template <typename CatalogT>
struct RuntimeCatalogIsSupported : std::false_type {};
/**
* @brief A specialization of the RuntimeCatalogIsSupported trait for a catalog of relations. This specialization checks if the relations in the catalog have runtime quantity identifiers and if the runtime quantity identifiers are unambiguous.
* @tparam RelationTs The relation types in the catalog.
*/
template <typename... RelationTs>
struct RuntimeCatalogIsSupported<relations::RelationCatalog<RelationTs...>> :
RuntimeRelationsAreSupported<(HasRuntimeQuantityIdentifiers<RelationTs>::value && ...), RelationTs...> {};
/**
* @brief A helper trait to determine if a catalog of relations is supported by the runtime system. This is used to ensure that all relations in a catalog have unique runtime quantity identifiers.
* @tparam CatalogT The catalog type to check.
*/
template <typename EOS, typename RelationT>
struct RuntimeRelationStorage;
/**
* @brief A specialization of the RuntimeRelationStorage trait for a relation. This specialization extracts the output and input quantity types from the relation and stores them in a tuple.
* @tparam EOS The EOS model type that supports the relation.
* @tparam OutputT The output quantity type of the relation.
* @tparam InputTs The input quantity types of the relation.
*/
template <typename EOS, typename OutputT, typename... InputTs>
struct RuntimeRelationStorage<EOS, relations::Relation<OutputT, InputTs...>> {
using RelationType = relations::Relation<OutputT, InputTs...>; ///< The relation type that this storage is for.
static_assert(sizeof...(InputTs) <= 64, "Runtime EOS relation descriptors support at most 64 inputs. If you need more than this you will need to implement your own custom RuntimeRelationStorage specialization");
/**
* @brief A static array of the runtime quantity identifiers for the input quantities of the relation. This is used to ensure that all relations in a catalog have unique runtime quantity identifiers.
*/
static constexpr std::array<dimensions::runtime::ThermodynamicQuantityID, sizeof...(InputTs)> inputQuantityIDs {
dimensions::runtime::thermodynamicQuantityID<InputTs>...
};
/**
* @brief A helper function to create a bitmask of the partial derivatives that are supported by the EOS model for the relation. This is used to ensure that all relations in a catalog have unique runtime quantity identifiers.
* @tparam Indices The indices of the input quantities to check.
* @return A bitmask of the partial derivatives that are supported by the EOS model for the relation.
*/
template <std::size_t... Indices>
[[nodiscard]] static consteval std::uint64_t makePartialDerivativeMask(std::index_sequence<Indices...>) {
using InputTuple = std::tuple<InputTs...>;
/* This is perhaps overly "clever" and self-indulgent.
*
* Instead of changing that I will provide a brief explanation of what this is doing.
* Effectively we are accumulating a bit mask. std::uint64_t on the left initializes the bit mask to all 0s
* Then we fold over all the indices of the input quantities
*
* As a reminder if you have not run into any of the other folding comments I have left scattered in the code base. C++17 introduced the folding operator,
* what we have here is what is called a right binary fold. Essentially we start with a value then repeatedly logically or it with the result of an iteration
* over the type pack Indices (which is some list of integers from 0 to sizeof...(InputTs) - 1). The | ... | is what indicates this, the first pipe shows we or
* the initial state (a 64 bit integer with all bits set to 0) with the result of the first iteration, then we or that result with the result of the second iteration, and so on until we have iterated over all indices.
* There is persistent confusion I feel with this since we do not have a loop variable as a standard for loop would; however, this is simply the syntax available
* in C++ (as of C++23). The looping is implicit in so far as the compiler implicitly generates the loop over Indices.
*
* We use std::tuple_element_t to get the type of the input quantity at the current index (recall that the fold loop goes over all indices)
* So for example if the tuple is <int, float, int, int> and we are at the second iteration of the folding loop (index 1) then std::tuple_element_t will evaluate to float
* This is contrived though since we are not storing primitives in the tuple but rather types derived from the ThermodynamicQuantity concept.
*
* Note where we evaluate this type. In the template argument for EOSSupportingPartialDerivative. Basically we are asking if, for the given EOS
* and relation type, does it support a partial derivative with respect to the input quantity at the current index. This evaluates at compile
* time to either true or false. There is then a ternary operator that selects a 1 (bit shifted by the current index to fit in the right bitmask location) or a 0 (which we don't need to bitshift)
*
* What that means is when all is said and done for an input type list made of n types there is an n bit integer (padded out to 64 bits) where each bit corresponds
* to whether the EOS supports a partial derivative with respect to the input quantity at that index. The first input quantity corresponds to the least significant bit and the last input quantity corresponds to the most significant bit.
*/
return (std::uint64_t{0} | ... | (models::EOSSupportingPartialDerivative<EOS, RelationType, std::tuple_element_t<Indices, InputTuple>> ? (std::uint64_t{1} << Indices) : std::uint64_t{0}));
}
static constexpr std::uint64_t partialDerivativeMask = makePartialDerivativeMask(std::index_sequence_for<InputTs...>{}); ///< A bitmask of the partial derivatives that are supported by the EOS model for the relation. Each bit corresponds to an input quantity, with the least significant bit corresponding to the first input quantity and the most significant bit corresponding to the last input quantity.
/**
* @brief A runtime relation descriptor for the relation. This is used to ensure that all relations in a catalog have unique runtime quantity identifiers.
*/
static constexpr dimensions::runtime::RuntimeRelationDescriptor descriptor {
.outputQuantity = dimensions::runtime::thermodynamicQuantityID<OutputT>,
.inputQuantities = std::span<const dimensions::runtime::ThermodynamicQuantityID>{inputQuantityIDs},
.partialDerivativeMask = partialDerivativeMask
};
};
/**
* @brief A helper struct to store the runtime relation descriptors for a catalog of relations. This is used to ensure that all relations in a catalog have unique runtime quantity identifiers.
* @tparam EOS The EOS model type that supports the relations in the catalog.
* @tparam CatalogT The catalog type to store the descriptors for.
*/
template <typename EOS, typename CatalogT>
struct RuntimeCatalogStorage;
/**
* @brief A specialization of the RuntimeCatalogStorage struct for a catalog of relations. This specialization stores the runtime relation descriptors for each relation in the catalog.
* @tparam EOS The EOS model type that supports the relations in the catalog.
* @tparam RelationTs The relation types in the catalog.
*/
template <typename EOS, typename... RelationTs>
struct RuntimeCatalogStorage<EOS, relations::RelationCatalog<RelationTs...>> {
static constexpr std::array descriptors {RuntimeRelationStorage<EOS, RelationTs>::descriptor...};
};
/**
* @brief A helper function to create a std::unexpected with a SERiFErrorCode. This is used to simplify the usage of std::expected in the runtime system.
* @param code The SERiFErrorCode to wrap in a std::unexpected.
* @return A std::unexpected containing the given SERiFErrorCode.
*/
[[nodiscard]] inline std::expected<double, utils::errors::SERiFErrorCode> runtimeEvaluationFailure(
const utils::errors::SERiFErrorCode code
) {
return std::unexpected{code};
}
/**
* @brief A helper function to evaluate a relation at runtime. This is used to simplify the usage of std::expected in the runtime system.
* @tparam EOS The EOS model type that supports the relation.
* @tparam OutputT The output quantity type of the relation.
* @tparam InputTs The input quantity types of the relation.
*
* @param eos The EOS model instance to use for evaluation.
* @param relation The relation to evaluate.
* @param inputValues The input values to use for evaluation.
* @return A std::expected containing the result of the evaluation or a SERiFErrorCode if the evaluation failed.
*/
template <typename EOS, typename OutputT, typename... InputTs>
[[nodiscard]] std::expected<double, utils::errors::SERiFErrorCode> evaluateRuntimeRelation(
const EOS& eos,
relations::Relation<OutputT, InputTs...> /*relation*/,
const std::span<const dimensions::runtime::RuntimeQuantityValue> inputValues
) {
/* Once again we find ourself at a use of folding. I recommend you read the comment earlier in this file regarding folding, or find the C++ docs on folding.
*
* The general premise here is that we call the evaluation function for the EOS with the given relation and input values. We fold over all indices using the final ... operator.
* This is an implicit loop generated by the compiler over the so called "parameter pack".
*/
const auto invoke_evaluate = [&]<std::size_t... Indices>(std::index_sequence<Indices...>) {
return eos::evaluate<OutputT>(eos, dimensions::QuantityValue<InputTs>{inputValues[Indices].value}...).value();
};
try {
return invoke_evaluate(std::make_index_sequence<sizeof...(InputTs)>{});
} catch (const EOSEvaluationError& e) { // We may want to reconsider using a try-catch for this, but this is good enough for now. Realistically if we want to not use a try catch we will need to change the signature of the EOS evaluation function to return a std::expected instead of throwing an exception. This is a larger change that I don't want to make right now.
return std::unexpected<utils::errors::SERiFErrorCode>{utils::errors::SERiFErrorCode::EOS_ERROR};
}
}
/**
* @brief A helper function to evaluate a partial derivative of a relation at runtime. This is used to simplify the usage of std::expected in the runtime system.
* @tparam InputQuantityT The input quantity type to take the partial derivative with respect to.
* @tparam EOS The EOS model type that supports the relation.
* @tparam OutputT The output quantity type of the relation.
* @tparam InputTs The input quantity types of the relation.
*
* @param eos The EOS model instance to use for evaluation.
* @param relation The relation to evaluate.
* @param withRespectTo The input quantity type to take the partial derivative with respect to.
* @param inputValues The input values to use for evaluation.
* @return A std::expected containing the result of the evaluation or a SERiFErrorCode if the evaluation failed.
*/
template <typename InputQuantityT, typename EOS, typename OutputT, typename... InputTs>
[[nodiscard]] bool tryRuntimePartialDerivative(
const EOS& eos,
relations::Relation<OutputT, InputTs...> /*relation*/,
const dimensions::runtime::ThermodynamicQuantityID withRespectTo,
const std::span<const dimensions::runtime::RuntimeQuantityValue> inputValues,
std::expected<double, utils::errors::SERiFErrorCode>& result
) {
if (withRespectTo != dimensions::runtime::thermodynamicQuantityID<InputQuantityT>) {
return false;
}
if constexpr (models::EOSSupportingPartialDerivative<EOS, relations::Relation<OutputT, InputTs...>, InputQuantityT>) {
const auto invoke_evaluate_partial_derivative = [&]<std::size_t... Indices>(std::index_sequence<Indices...>) {
return eos::partial_derivative<OutputT, InputQuantityT>(eos, dimensions::QuantityValue<InputTs>{inputValues[Indices].value}...).value();
};
try {
result = invoke_evaluate_partial_derivative(std::index_sequence_for<InputTs...>{});
} catch (const EOSEvaluationError& /*e*/) {
result = std::unexpected{utils::errors::SERiFErrorCode::EOS_ERROR};
}
} else {
result = runtimeEvaluationFailure(utils::errors::SERiFErrorCode::EOS_UNSUPPORTED_DERIVATIVE);
}
return true;
}
/**
* @brief A helper function to evaluate a partial derivative of a relation at runtime. This is used to simplify the usage of std::expected in the runtime system.
* @tparam EOS The EOS model type that supports the relation.
* @tparam OutputT The output quantity type of the relation.
* @tparam InputTs The input quantity types of the relation.
*
* @param eos The EOS model instance to use for evaluation.
* @param relation The relation to evaluate.
* @param withRespectTo The input quantity type to take the partial derivative with respect to.
* @param inputValues The input values to use for evaluation.
* @return A std::expected containing the result of the evaluation or a SERiFErrorCode if the evaluation failed.
*/
template <typename EOS, typename OutputT, typename... InputTs>
[[nodiscard]] std::expected<double, utils::errors::SERiFErrorCode> evaluateRuntimePartialDerivative (
const EOS& eos,
relations::Relation<OutputT, InputTs...> relation,
const dimensions::runtime::ThermodynamicQuantityID withRespectTo,
const std::span<const dimensions::runtime::RuntimeQuantityValue> inputValues
) {
std::expected<double, utils::errors::SERiFErrorCode> result = runtimeEvaluationFailure(
utils::errors::SERiFErrorCode::EOS_UNSUPPORTED_DERIVATIVE
);
[[maybe_unused]] const bool matched = (tryRuntimePartialDerivative<InputTs>(eos, relation, withRespectTo, inputValues, result) || ...); // Here we fold over all input and try to find one that we can evaluate.
return result;
}
/**
* @brief A helper function to determine if a relation matches the given output quantity and input quantities at runtime. This is used to simplify the usage of std::expected in the runtime system.
* @tparam EOS The EOS model type that supports the relation.
* @tparam RelationT The relation type to check.
*
* @param outputQuantityID The output quantity type to check.
* @param inputValues The input values to check.
* @return True if the relation matches the given output quantity and input quantities, false otherwise.
*/
template <typename EOS, typename RelationT>
[[nodiscard]] bool runtimeRelationMatches(
const dimensions::runtime::ThermodynamicQuantityID outputQuantityID,
const std::span<const dimensions::runtime::RuntimeQuantityValue> inputValues
) {
const dimensions::runtime::RuntimeRelationDescriptor& descriptor = RuntimeRelationStorage<EOS, RelationT>::descriptor;
if (descriptor.outputQuantity != outputQuantityID || descriptor.inputQuantities.size() != inputValues.size()) {
return false;
}
for (std::size_t i = 0; i < inputValues.size(); ++i) {
if (descriptor.inputQuantities[i] != inputValues[i].id) {
return false;
}
}
return true;
}
/**
* @brief A helper struct to dispatch runtime evaluation of relations in a catalog. This is used to simplify the usage of std::expected in the runtime system.
* @tparam EOS The EOS model type that supports the relations in the catalog.
* @tparam CatalogT The catalog type to dispatch.
*/
template <typename EOS, typename CatalogT>
struct RuntimeCatalogDispatch;
/**
* @brief A specialization of the RuntimeCatalogDispatch struct for a catalog of relations. This specialization dispatches runtime evaluation of relations in the catalog.
* @tparam EOS The EOS model type that supports the relations in the catalog.
* @tparam RelationTs The relation types in the catalog.
*/
template <typename EOS, typename...RelationTs>
struct RuntimeCatalogDispatch<EOS, relations::RelationCatalog<RelationTs...>> {
/**
* @brief A helper function to evaluate a relation in the catalog at runtime. This is used to simplify the usage of std::expected in the runtime system.
* @param object A pointer to the EOS model instance to use for evaluation.
* @param outputQuantity The output quantity type to evaluate.
* @param inputValues The input values to use for evaluation.
* @return A std::expected containing the result of the evaluation or a SERiFErrorCode if the evaluation failed.
*/
[[nodiscard]] static std::expected<double, utils::errors::SERiFErrorCode> evaluate(
const void *object,
const dimensions::runtime::ThermodynamicQuantityID& outputQuantity,
const std::span<const dimensions::runtime::RuntimeQuantityValue> inputValues
) {
const auto &eos = *static_cast<const EOS*>(object); // Not sure if polymorphic type erasure is the right tool here. It works but we may want to revisit it if it precent compiler optimizations. This does however let us pass any eos model to the runtime view without needing to know the type at compile time (which is therefore helpful when building extension systems in other languages like python). Further, it may make sense to take the performance hit here and use something like a dynamic_cast to ensure that the object is actually of the correct type. This would be a runtime check but it would be a more robust check than just blindly casting to the expected type.
std::expected<double, utils::errors::SERiFErrorCode> result = runtimeEvaluationFailure(
utils::errors::SERiFErrorCode::EOS_UNSUPPORTED_RELATION
);
/* Fold over all relations in the catalog and check if any of them match the requested output quantity and input quantities. If the match is found then
* evaluate the relation and return the result and set matched to true, if no match is found then set matched to false.
*/
[[maybe_unused]] const bool matched = ((runtimeRelationMatches<EOS, RelationTs>(outputQuantity, inputValues) ? (result = evaluateRuntimeRelation(eos, RelationTs{}, inputValues), true) : false) || ...);
return result;
}
/**
* @brief A helper function to evaluate a partial derivative of a relation in the catalog at runtime. This is used to simplify the usage of std::expected in the runtime system.
* @param object A pointer to the EOS model instance to use for evaluation.
* @param outputQuantity The output quantity type to evaluate.
* @param withRespectTo The input quantity type to take the partial derivative with respect to.
* @param inputValues The input values to use for evaluation.
* @return A std::expected containing the result of the evaluation or a SERiFErrorCode if the evaluation failed.
*/
[[nodiscard]] static std::expected<double, utils::errors::SERiFErrorCode> partial_derivative(
const void *object,
const dimensions::runtime::ThermodynamicQuantityID& outputQuantity,
const dimensions::runtime::ThermodynamicQuantityID& withRespectTo,
const std::span<const dimensions::runtime::RuntimeQuantityValue> inputValues
) {
const auto &eos = *static_cast<const EOS*>(object); // Same Comment as above
std::expected<double, utils::errors::SERiFErrorCode> result = runtimeEvaluationFailure(
utils::errors::SERiFErrorCode::EOS_UNSUPPORTED_RELATION
);
[[maybe_unused]] const bool matched = ((runtimeRelationMatches<EOS, RelationTs>(outputQuantity, inputValues) ? (result = evaluateRuntimePartialDerivative(eos, RelationTs{}, withRespectTo, inputValues), true) : false) || ...);
return result;
}
};
}

318
src/include/serif/eos/runtime/views.hpp
View File

@@ -1,318 +0,0 @@
#pragma once
/**
* @file views.hpp
* @author Emily M. Boudreaux
*
* @brief File containing the EOSView class. This class is used to provide a runtime view of an EOS model and its relations.
*/
#include <span>
#include "serif/eos/runtime/concepts.hpp"
namespace serif::eos::runtime {
/**
* @brief A class to provide a runtime view of an EOS model and its relations. This class is used to provide a uniform interface for evaluating relations and partial derivatives at runtime.
*/
class EOSView final {
public:
/**
* @brief Constructs an EOSView from a runtime EOS model. This constructor is used to create a view of the EOS model and its relations at runtime.
* @tparam EOS The type of the EOS model to view.
* @param eos The EOS model instance to view.
*/
template <typename EOS>
requires RuntimeEOSModel<std::remove_cvref_t<EOS>>
explicit EOSView(const EOS& eos) noexcept :
m_object(std::addressof(eos)),
m_relations(runtimeRelationDescriptors<std::remove_cvref_t<EOS>>()),
m_evaluate(&RuntimeEOSAdapter<std::remove_cvref_t<EOS>>::evaluate),
m_partialDerivative(&RuntimeEOSAdapter<std::remove_cvref_t<EOS>>::partial_derivative) {}
/**
* @brief Deleted constructor to prevent constructing an EOSView from an rvalue EOS model. This is to prevent dangling references and ensure that the EOS model instance outlives the EOSView.
* @tparam EOS The type of the EOS model to view.
* @param eos The EOS model instance to view.
*/
template <typename EOS>
requires RuntimeEOSModel<std::remove_cvref_t<EOS>>
explicit EOSView(const EOS&&) = delete;
/**
* @brief Returns a span of the runtime relation descriptors for the EOS model.
* @return A span of the runtime relation descriptors.
*
* @todo Move this to the implementation file to clean things up a bit.
*/
[[nodiscard]] std::span<const dimensions::runtime::RuntimeRelationDescriptor> relations() const noexcept {
return m_relations;
}
/**
* @brief Determines if the EOS model supports a relation with the given output quantity and input quantities at runtime. This is used to check if a relation can be evaluated at runtime.
* @param outputQuantity The output quantity to check.
* @param inputQuantities The input quantities to check.
* @return True if the relation is supported, false otherwise.
*/
[[nodiscard]] bool supports(
const dimensions::runtime::ThermodynamicQuantityID outputQuantity,
const std::span<const dimensions::runtime::ThermodynamicQuantityID> inputQuantities
) const noexcept {
return find_relation(outputQuantity, inputQuantities) != nullptr;
}
/**
* @brief Evaluates a relation with the given output quantity and input quantities at runtime. This is used to evaluate relations at runtime.
* @param outputQuantity The output quantity to evaluate.
* @param inputQuantities The input quantities to evaluate.
* @return A std::expected containing the result of the evaluation or a SERiFErrorCode if the evaluation failed.
*/
[[nodiscard]] std::expected<dimensions::runtime::RuntimeQuantityValue, utils::errors::SERiFErrorCode> try_evaluate(
const dimensions::runtime::ThermodynamicQuantityID outputQuantity,
const std::span<const dimensions::runtime::RuntimeQuantityValue> inputQuantities
) const {
const auto validation = validate_relation_request(outputQuantity, inputQuantities);
if (!validation.has_value()) {
return std::unexpected{validation.error()};
}
auto result = m_evaluate(m_object, outputQuantity, inputQuantities);
if (!result.has_value()) {
return std::unexpected{result.error()};
}
return dimensions::runtime::RuntimeQuantityValue{outputQuantity, *result};
}
/**
* @brief Evaluates a relation with the given output quantity and input quantities at runtime. This is used to evaluate relations at runtime.
* @tparam OutputQuantity The output quantity to evaluate.
* @tparam InputValues The input values to evaluate.
* @return A std::expected containing the result of the evaluation or a SERiFErrorCode if the evaluation failed.
*/
template <dimensions::runtime::RuntimeIdentifiedThermodynamicQuantity OutputQuantity, dimensions::ThermodynamicQuantityValueType... InputValues>
[[nodiscard]] std::expected<dimensions::QuantityValue<OutputQuantity>, utils::errors::SERiFErrorCode> try_evaluate(
const InputValues... inputValues
) const {
constexpr bool inputsHaveRuntimeIdentifiers = (dimensions::runtime::RuntimeIdentifiedThermodynamicQuantity<dimensions::QuantityOfT<InputValues>> && ...);
static_assert(inputsHaveRuntimeIdentifiers, "Every runtime EOS input quantity needs a stable identifier.");
const std::array<dimensions::runtime::RuntimeQuantityValue, sizeof...(InputValues)> runtimeInputs {
dimensions::runtime::RuntimeQuantityValue{
.id = dimensions::runtime::thermodynamicQuantityID<dimensions::QuantityOfT<InputValues>>, .value = inputValues.value()
}...
};
auto result = try_evaluate(
dimensions::runtime::thermodynamicQuantityID<OutputQuantity>, std::span<const dimensions::runtime::RuntimeQuantityValue>{runtimeInputs}
);
if (!result.has_value()) {
return std::unexpected<utils::errors::SERiFErrorCode>{result.error()};
}
return dimensions::QuantityValue<OutputQuantity>{result->value};
}
/**
* @brief Evaluates a partial derivative of a relation with the given output quantity and input quantities at runtime. This is used to evaluate partial derivatives at runtime.
* @param outputQuantity The output quantity to evaluate.
* @param withRespectTo The input quantity with respect to which the partial derivative is taken.
* @param inputValues The input values to evaluate.
* @return A std::expected containing the result of the evaluation or a SERiFErrorCode if the evaluation failed.
*/
[[nodiscard]] std::expected<double, utils::errors::SERiFErrorCode> try_partial_derivative(
const dimensions::runtime::ThermodynamicQuantityID outputQuantity,
const dimensions::runtime::ThermodynamicQuantityID withRespectTo,
const std::span<const dimensions::runtime::RuntimeQuantityValue> inputValues
) const {
const auto validation = validate_relation_request(outputQuantity, inputValues);
if (!validation.has_value()) {
return std::unexpected<utils::errors::SERiFErrorCode>{validation.error()};
}
const dimensions::runtime::RuntimeRelationDescriptor &descriptor = **validation;
bool derivativeAvailable = false;
for (std::size_t index = 0; index < descriptor.inputQuantities.size(); ++index) {
if (descriptor.inputQuantities[index] == withRespectTo) {
derivativeAvailable = descriptor.hasPartialDerivative(index);
break;
}
}
if (!derivativeAvailable) {
return runtime_failure<double>(utils::errors::SERiFErrorCode::EOS_UNSUPPORTED_DERIVATIVE);
}
return m_partialDerivative(m_object, outputQuantity, withRespectTo, inputValues);
}
/**
* @brief Evaluates a partial derivative of a relation with the given output quantity and input quantities at runtime. This is used to evaluate partial derivatives at runtime.
* @tparam OutputQuantity The output quantity to evaluate.
* @tparam InputQuantity The input quantity with respect to which the partial derivative is taken.
* @tparam InputValues The types of the input values.
* @param inputValues The input values to evaluate.
* @return A std::expected containing the result of the evaluation or a SERiFErrorCode if the evaluation failed.
*/
template <dimensions::runtime::RuntimeIdentifiedThermodynamicQuantity OutputQuantity, dimensions::runtime::RuntimeIdentifiedThermodynamicQuantity InputQuantity, dimensions::QuantityValueType... InputValues>
[[nodiscard]] std::expected<dimensions::PartialDerivative<OutputQuantity, InputQuantity>, utils::errors::SERiFErrorCode> try_partial_derivative(
const InputValues... inputValues
) const {
constexpr bool inputHaveRuntimeIdentifiers = (dimensions::runtime::RuntimeIdentifiedThermodynamicQuantity<dimensions::QuantityOfT<InputValues>> && ...);
static_assert(inputHaveRuntimeIdentifiers, "Every runtime EOS input quantity needs a stable identifier.");
const std::array<dimensions::runtime::RuntimeQuantityValue, sizeof...(InputValues)> runtimeInputs{
dimensions::runtime::RuntimeQuantityValue{.id = dimensions::runtime::thermodynamicQuantityID<dimensions::QuantityOfT<InputValues>>, .value = inputValues.value()}...
};
auto result = try_partial_derivative(
dimensions::runtime::thermodynamicQuantityID<OutputQuantity>, // output
dimensions::runtime::thermodynamicQuantityID<InputQuantity>, // input
std::span<const dimensions::runtime::RuntimeQuantityValue>{runtimeInputs} // where to evaluate
);
if (!result.has_value()) {
return std::unexpected<utils::errors::SERiFErrorCode>{result.error()};
}
return dimensions::PartialDerivative<OutputQuantity, InputQuantity>{*result};
}
private: // Type aliases
/**
* @brief A type alias for a function pointer to a runtime evaluation function. This is used to store the evaluation function for the EOS model.
*/
using RuntimeEvaluateFunction = std::expected<double, utils::errors::SERiFErrorCode> (*)(
const void*,
const dimensions::runtime::ThermodynamicQuantityID&,
std::span<const dimensions::runtime::RuntimeQuantityValue>
);
/**
* @brief A type alias for a function pointer to a runtime partial derivative evaluation function. This is used to store the partial derivative evaluation function for the EOS model.
*/
using RuntimePartialDerivativeFunction = std::expected<double, utils::errors::SERiFErrorCode> (*)(
const void*,
const dimensions::runtime::ThermodynamicQuantityID&,
const dimensions::runtime::ThermodynamicQuantityID&,
std::span<const dimensions::runtime::RuntimeQuantityValue>
);
private: // Private methods
/**
* @brief Finds a relation with the given output quantity and input quantities at runtime. This is used to find relations at runtime.
* @param outputQuantity The output quantity to find the relation for.
* @param inputQuantities The input quantities to find the relation for.
* @return A pointer to the runtime relation descriptor if found, nullptr otherwise.
*
* @todo move to implementation file to clean things up a bit.
*/
[[nodiscard]] const dimensions::runtime::RuntimeRelationDescriptor *find_relation(
const dimensions::runtime::ThermodynamicQuantityID outputQuantity,
const std::span<const dimensions::runtime::ThermodynamicQuantityID> inputQuantities
) const noexcept {
for (const auto& descriptor : m_relations) {
if (descriptor.outputQuantity != outputQuantity || descriptor.inputQuantities.size() != inputQuantities.size()) {
continue;
}
bool matches = true;
for (std::size_t index = 0; index < inputQuantities.size(); ++index) {
if (descriptor.inputQuantities[index] != inputQuantities[index]) {
matches = false;
break;
}
}
if (matches) {
return std::addressof(descriptor);
}
}
return nullptr;
}
/**
* @brief Validates a relation request with the given output quantity and input quantities at runtime. This is used to validate relations at runtime.
* @param outputQuantity The output quantity to validate.
* @param inputValues The input values to validate.
* @return A std::expected containing a pointer to the runtime relation descriptor if valid, or a SERiFErrorCode if invalid.
*/
[[nodiscard]] std::expected<const dimensions::runtime::RuntimeRelationDescriptor *, utils::errors::SERiFErrorCode> validate_relation_request(
const dimensions::runtime::ThermodynamicQuantityID outputQuantity,
const std::span<const dimensions::runtime::RuntimeQuantityValue> inputValues
) const {
bool outputAvailable = false;
bool inputCountAvailable = false;
for (const auto& descriptor : m_relations) {
if (descriptor.outputQuantity != outputQuantity) {
continue;
}
outputAvailable = true;
if (descriptor.inputQuantities.size() != inputValues.size()) {
continue;
}
inputCountAvailable = true;
bool matches = true;
for (std::size_t index = 0; index < inputValues.size(); ++index) {
if (descriptor.inputQuantities[index] != inputValues[index].id) {
matches = false;
break;
}
}
if (matches) {
return std::addressof(descriptor);
}
}
if (!outputAvailable) {
return runtime_failure<const dimensions::runtime::RuntimeRelationDescriptor *>(
utils::errors::SERiFErrorCode::EOS_UNSUPPORTED_RELATION
);
}
if (!inputCountAvailable) {
return runtime_failure<const dimensions::runtime::RuntimeRelationDescriptor *>(
utils::errors::SERiFErrorCode::EOS_WRONG_INPUT_COUNT
);
}
return runtime_failure<const dimensions::runtime::RuntimeRelationDescriptor *>(
utils::errors::SERiFErrorCode::EOS_WRONG_INPUT_QUANTITY
);
}
/**
* @brief A helper function to handle runtime failures. This is used to simplify the usage of std::expected in the runtime system.
* @tparam Value The type of the value to return.
* @param code The SERiFErrorCode to return.
* @return A std::unexpected containing the SERiFErrorCode.
*/
template <typename Value>
[[nodiscard]] static std::expected<Value, utils::errors::SERiFErrorCode> runtime_failure(
const utils::errors::SERiFErrorCode code
) {
return std::unexpected{code};
}
private: // Private members
const void *m_object; ///< Type erased pointer to the EOS model instance.
std::span<const dimensions::runtime::RuntimeRelationDescriptor> m_relations; ///< A span of the runtime relation descriptors for the EOS model.
RuntimeEvaluateFunction m_evaluate; ///< A function pointer to the runtime evaluation function for the EOS model.
RuntimePartialDerivativeFunction m_partialDerivative; ///< A function pointer to the runtime partial derivative evaluation function for the EOS model.
};
}

174
src/lib/serif/discretization/backend/mfem/basis_evaluation.cpp
View File

@@ -1,174 +0,0 @@
#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<SERiFErrorCode::INVALID_ARGUMENT, std::invalid_argument>(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<SERiFErrorCode::INVALID_DOF_COUNT, std::invalid_argument>(
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<SERiFErrorCode::INVALID_DOF_COUNT, std::invalid_argument>(
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<bool>(m_tensor);
}
}

28
src/lib/serif/eos/models/polytropic.cpp
View File

@@ -23,7 +23,7 @@ namespace serif::eos::models {
validate_finite<EOSEvaluationError, SERiFErrorCode::EOS_OUTSIDE_DOMAIN>(m_polytropic_index);
validate_greater_than<EOSEvaluationError, SERiFErrorCode::EOS_OUTSIDE_DOMAIN>(m_polytropic_index, 1.0);
validate_finite<EOSEvaluationError, SERiFErrorCode::EOS_OUTSIDE_DOMAIN>(m_polytropic_constant);
validate_greater_than<EOSEvaluationError, SERiFErrorCode::EOS_OUTSIDE_DOMAIN>(m_polytropic_constant, 0.0);
utils::errors::serif_ensure<SERiFErrorCode::EOS_OUTSIDE_DOMAIN, EOSEvaluationError>(m_polytropic_constant > 0.0, "The polytropic constant must be positive. Received {}.", m_polytropic_constant);
}
[[nodiscard]] double Polytrope::polytropic_index() const noexcept {
@@ -52,6 +52,27 @@ namespace serif::eos::models {
};
}
[[nodiscard]] DensityValue Polytrope::evaluate(DensityFromPressure, const PressureValue pressure) const {
validate_nonnegativity<EOSEvaluationError, SERiFErrorCode::EOS_INVALID_THERMODYNAMIC_INPUT>(pressure.value());
if (pressure.value() == 0.0) {
return DensityValue{0.0};
}
return DensityValue{std::pow(pressure.value() / m_polytropic_constant, m_polytropic_index / (m_polytropic_index + 1.0))};
}
[[nodiscard]] SecondPartialDerivative<Pressure, MassDensity, MassDensity> Polytrope::second_partial_derivative(
PressureFromDensity, WithRespectTo<MassDensity>, WithRespectTo<MassDensity>, const DensityValue density
) const {
validate_nonnegativity<EOSEvaluationError, SERiFErrorCode::EOS_INVALID_THERMODYNAMIC_INPUT>(density.value());
utils::errors::serif_ensure<SERiFErrorCode::EOS_INVALID_THERMODYNAMIC_INPUT, EOSEvaluationError>(
density.value() > 0.0 || m_polytropic_index == 1.0,
"The second pressure derivative is singular at zero density for polytropic index {}. Apply the vacuum density trace directly and evaluate bulk Newton coefficients only at positive density.",
m_polytropic_index);
const double exponent = 1.0 / m_polytropic_index;
return SecondPartialDerivative<Pressure, MassDensity, MassDensity>{
m_polytropic_constant * (1.0 + exponent) * exponent * std::pow(density.value(), exponent - 1.0)};
}
[[nodiscard]] SpecificEnthalpyValue Polytrope::evaluate(SpecificEnthalpyFromDensity, const DensityValue density) const {
validate_nonnegativity<EOSEvaluationError, SERiFErrorCode::EOS_INVALID_THERMODYNAMIC_INPUT>(density.value());
if (density.value() == 0.0) {
@@ -65,6 +86,11 @@ namespace serif::eos::models {
};
}
/** @brief Reuse the polytropic enthalpy relation, which is also this EOS's pressure primitive. */
[[nodiscard]] BarotropicPressurePotentialValue Polytrope::evaluate(PressurePotentialFromDensity, const DensityValue density) const {
return BarotropicPressurePotentialValue{evaluate(SpecificEnthalpyFromDensity{}, density).value()};
}
[[nodiscard]] DensityValue Polytrope::evaluate(DensityFromSpecificEnthalpy, SpecificEnthalpyValue specific_enthalpy) const {
validate_finite<EOSEvaluationError, SERiFErrorCode::EOS_INVALID_THERMODYNAMIC_INPUT>(specific_enthalpy.value());

22
src/lib/serif/mapping/coefficients.cpp
View File

@@ -1,22 +0,0 @@
#include "serif/mapping/coefficients.hpp"
#include "serif/utils/error/errors.hpp"
namespace serif::mapping {
using namespace serif::utils::errors;
PhysicalPositionFunctionCoefficient::PhysicalPositionFunctionCoefficient(
const DomainMapper &mapper,
const mfem::GridFunction &displacement,
const mfem::GridFunction &compactification_coordinate,
Func func
) :
m_f(std::move(func)),
m_mapping(mapper, displacement, compactification_coordinate) { }
double PhysicalPositionFunctionCoefficient::Eval(mfem::ElementTransformation &transformation, const mfem::IntegrationPoint &ip) {
transformation.SetIntPoint(&ip);
MappingPointContext context;
serif_ensure<SERiFErrorCode::INVALID_MAPPING>(m_mapping.evaluate_point(transformation, ip, context) == MappingStatus::VALID, "Mapping evaluation failed during PhysicalPositionCoefficient evaluation.");
return m_f(context.physical_position);
}
}

273
src/lib/serif/mapping/compactification/kelvin.cpp
View File

@@ -1,273 +0,0 @@
#include <cmath>
#include "serif/utils/checks/finite.hpp"
#include "serif/mapping/compactification/kelvin.hpp"
#include "serif/utils/user_options/mapping/compactification/kelvin.hpp"
#include "serif/utils/error/errors.hpp"
namespace serif::mapping::compactification {
using namespace serif::utils::user_options;
using namespace serif::utils::errors;
using namespace serif::utils::checks;
KelvinCompactification::KelvinCompactification(const KelvinCompactificationOptions& options) : m_options(options) {
if (!std::isfinite(m_options.r_star_ref)) {
serif_error<std::invalid_argument, SERiFErrorCode::NON_FINITE_MAPPING_INPUT>("Kelvin compactification requires a finite reference stellar radius. Received {} ", m_options.r_star_ref);
}
if (!std::isfinite(m_options.r_inf_ref)) {
serif_error<std::invalid_argument, SERiFErrorCode::NON_FINITE_MAPPING_INPUT>("Kelvin compactification requires a finite reference infinity radius. Received {} ", m_options.r_star_ref);
}
if (m_options.r_star_ref <= 0 || m_options.r_inf_ref < m_options.r_star_ref)
{
serif_error<std::invalid_argument, SERiFErrorCode::IMPROPERLY_ORDERED_REFERENCE_RADII>("Kelvin compactification requires 0 < r_star_ref < r_inf_ref. Got r_star_ref={} and r_inf_ref={}", m_options.r_star_ref, m_options.r_inf_ref);
}
if (!std::isfinite(m_options.coordinate_tolerance) || m_options.coordinate_tolerance < 0.0 || m_options.coordinate_tolerance > 1.0)
{
serif_error<std::invalid_argument, SERiFErrorCode::INVALID_COORDINATE_TOLERANCE>("Kelvin compactification coordinate tolerance must be finite and on [0, 1). Received: {}", m_options.coordinate_tolerance);
}
}
bool KelvinCompactification::is_affine_in_displacement() const noexcept {
return true; // Always true for kelvin compactification
}
MappingStatus KelvinCompactification::evaluate(const ExteriorMapInput& input, ExteriorMapResult& result) const {
const int dimension = input.reference_position.Size();
if (dimension <= 0 || input.displaced_position.Size() != dimension || input.compactification_coordinate_gradient.Size() != dimension) {
return MappingStatus::INVALID_DIMENSION;
}
if (input.displacement_jacobian.Height() != dimension || input.displacement_jacobian.Width() != dimension) {
return MappingStatus::INVALID_DIMENSION;
}
const bool reference_pos_is_finite = is_finite(input.reference_position);
const bool displaced_pos_is_finite = is_finite(input.displaced_position);
const bool coord_grad_is_finite = is_finite(input.compactification_coordinate_gradient);
const bool displacement_jacobian_is_finite = is_finite(input.displacement_jacobian);
if (!(reference_pos_is_finite && displaced_pos_is_finite && coord_grad_is_finite && displacement_jacobian_is_finite)) {
return MappingStatus::NON_FINITE_MAPPING_INPUT;
}
RadialFactors factors{};
// The key radial factors we use are the scale (which stretches the finite computational domain into the infinite physical domain) and the scale derivative (which is used to compute the mapping jacobian)
const MappingStatus factor_status = compute_radial_factors(input.compactification_coordinate, factors);
if (factor_status != MappingStatus::VALID) {
return factor_status;
}
result.physical_position.SetSize(dimension);
result.mapping_jacobian.SetSize(dimension, dimension);
for (int i = 0; i < dimension; ++i) {
// Note how the physical position is just the product of the displaced position and the scale factor.
result.physical_position(i) = factors.scale * input.displaced_position(i);
// The mapping jacobian comes from trivial application of the product rule
// recall: r_i = scale * x_i where r is the physical position and x is the displaced position.
// then we can differentiate wrt. X_j holding nothing fixed. Note the capital X here, this is the mesh coordinate not the displaced position.
// Lets call this jacobian F
// F = \frac{\partial r_i}{\partial X_{j}}
// F then tells us how the physical position changes as we move along mesh coordinates
// Lets then apply this to the function we have for the kelvin compactification
// F = scale * \frac{\partial x_i}{\partial X_j} + x_i * \frac{\partial scale}{\partial X_j}
// Below you can see the displacement jacobian (\frac{\partial x_i}{\partial X_j}) and the scale derivative (\frac{\partial scale}{\partial X_j}) being applied to compute the mapping jacobian.
for (int j = 0; j < dimension; ++j) {
const double scale_gradient = factors.scale_derivative * input.compactification_coordinate_gradient(j);
result.mapping_jacobian(i, j) = factors.scale * input.displacement_jacobian(i, j) + input.displaced_position(i) * scale_gradient;
}
}
if (!is_finite(result.physical_position) || !is_finite(result.mapping_jacobian)) {
return MappingStatus::NON_FINITE_MAPPING_RESULT;
}
const double mapping_determinant = result.mapping_jacobian.Det();
if (!std::isfinite(mapping_determinant)) {
return MappingStatus::NON_FINITE_MAPPING_RESULT;
}
if (mapping_determinant <= 0.0) {
return MappingStatus::NON_POSITIVE_DETERMINATE;
}
return MappingStatus::VALID;
}
MappingStatus KelvinCompactification::evaluate_variation(const ExteriorMapInput& input, const ExteriorMapResult& result, const ExteriorMapDirection& direction, ExteriorMapVariation& variation) const {
const int dimension = input.reference_position.Size();
if (dimension <= 0 || input.displaced_position.Size() != dimension || input.compactification_coordinate_gradient.Size() != dimension) {
return MappingStatus::INVALID_DIMENSION;
}
if (input.displacement_jacobian.Height() != dimension || input.displacement_jacobian.Width() != dimension) {
return MappingStatus::INVALID_DIMENSION;
}
if (result.physical_position.Size() != dimension || result.mapping_jacobian.Height() != dimension || result.mapping_jacobian.Width() != dimension) {
return MappingStatus::INVALID_DIMENSION;
}
if (direction.displaced_position_variation.Size() != dimension || direction.displacement_jacobian_variation.Height() != dimension || direction.displacement_jacobian_variation.Width() != dimension) {
return MappingStatus::INVALID_DIMENSION;
}
const bool reference_position_is_finite = is_finite(input.reference_position);
const bool displaced_position_is_finite = is_finite(input.displaced_position);
const bool compactification_coord_grad_is_finite = is_finite(input.compactification_coordinate_gradient);
const bool displacement_jacobian_is_finite = is_finite(input.displacement_jacobian);
if (!(reference_position_is_finite && displaced_position_is_finite && compactification_coord_grad_is_finite && displacement_jacobian_is_finite)) {
return MappingStatus::NON_FINITE_MAPPING_INPUT;
}
const bool result_physical_position_is_finite = is_finite(result.physical_position);
const bool result_mapping_jacobian_is_finite = is_finite(result.mapping_jacobian);
const bool direction_displaced_position_is_finite = is_finite(direction.displaced_position_variation);
const bool direction_displacement_jacobian_is_finite = is_finite(direction.displacement_jacobian_variation);
if (!(result_physical_position_is_finite && result_mapping_jacobian_is_finite && direction_displaced_position_is_finite && direction_displacement_jacobian_is_finite)) {
return MappingStatus::NON_FINITE_MAPPING_INPUT; // Note we are checking a parameter called result here; however, this is the result from an evaluate call, it is an input to this function.
}
RadialFactors factors{};
const MappingStatus factor_status = compute_radial_factors(input.compactification_coordinate, factors);
if (factor_status != MappingStatus::VALID) {
return factor_status;
}
variation.physical_position_variation.SetSize(dimension);
variation.mapping_jacobian_variation.SetSize(dimension, dimension);
for (int i = 0; i < dimension; ++i) {
variation.physical_position_variation(i) = factors.scale * direction.displaced_position_variation(i);
for (int j = 0; j < dimension; ++j) {
const double scale_gradient = factors.scale_derivative * input.compactification_coordinate_gradient(j);
variation.mapping_jacobian_variation(i, j) = factors.scale * direction.displacement_jacobian_variation(i, j) + direction.displaced_position_variation(i) * scale_gradient;
}
}
if (!is_finite(variation.physical_position_variation) || !is_finite(variation.mapping_jacobian_variation)) {
return MappingStatus::NON_FINITE_MAPPING_RESULT;
}
return MappingStatus::VALID;
}
[[nodiscard]] std::string_view KelvinCompactification::get_name() const noexcept {
return m_name;
}
[[nodiscard]] double KelvinCompactification::get_reference_stellar_radius() const noexcept {
return m_options.r_star_ref;
}
[[nodiscard]] double KelvinCompactification::get_reference_infinity_radius() const noexcept {
return m_options.r_inf_ref;
}
[[nodiscard]] double KelvinCompactification::get_coordinate_tolerance() const noexcept {
return m_options.coordinate_tolerance;
}
[[nodiscard]] MappingStatus KelvinCompactification::evaluate_jacobian_variation_factors(double compactification_coordinate, const mfem::Vector& coordinateGradient, double& scale, mfem::Vector& scaleGradient) const {
if (coordinateGradient.Size() <= 0) {
return MappingStatus::INVALID_DIMENSION;
}
if (!is_finite(coordinateGradient)) {
return MappingStatus::NON_FINITE_MAPPING_INPUT;
}
RadialFactors factors{};
const MappingStatus status = compute_radial_factors(compactification_coordinate, factors);
if (status != MappingStatus::VALID) {
return status;
}
scale = factors.scale;
scaleGradient = coordinateGradient;
scaleGradient *= factors.scale_derivative;
if (!is_finite(scaleGradient)) {
return MappingStatus::NON_FINITE_MAPPING_RESULT;
}
return MappingStatus::VALID;
}
std::unique_ptr<const ExteriorDomainMap> KelvinCompactification::clone() const {
return std::make_unique<const KelvinCompactification>(m_options);
}
MappingStatus KelvinCompactification::compute_radial_factors(double compactification_coordinate, RadialFactors& radial_factors) const {
if (!std::isfinite(compactification_coordinate)) {
return MappingStatus::NON_FINITE_MAPPING_INPUT;
}
const double tolerance = m_options.coordinate_tolerance;
if (compactification_coordinate < -tolerance || compactification_coordinate > 1.0 + tolerance) {
return MappingStatus::OUTSIDE_REFERENCE_DOMAIN;
}
double coordinate = compactification_coordinate;
if (coordinate < 0.0) {
coordinate = 0.0; // Not sure this is the best behavior here...
}
if (coordinate >= 1.0 - tolerance) {
return MappingStatus::AT_COMPACTIFIED_INFINITY;
}
// Here we need to do some transformations from the options defined on the mesh to useful computational coordinates
// r_inf_ref is the computational / reference radius of the infinity surface (the edge of the entire domain) and r_star_ref is the radius of the spherical
// stellar model inscribed within. Therefore radial extent is the computational radial distance between the stellar surface and the infinity surface.
// Note that this is separate from the parameterize compactification coordinate.
const double radial_extent = m_options.r_inf_ref - m_options.r_star_ref;
// This places us at the correct spot in computational space given the current compactification coordinate. Say you have compactification = 0.5,
// an r_star_ref of 2 and a radial extent of 3, this this will place you at 2 + 0.5 * 3 = 3.5 in computational space, which is half way between the stellar surface and the infinity surface.
const double computational_radius = m_options.r_star_ref + coordinate * radial_extent;
if (!std::isfinite(computational_radius) || computational_radius <= 0.0) {
return MappingStatus::INVALID_REFERENCE_RADIUS;
}
// Invert the exterior coordinate so it runs from 0 at the star to 1 at compactified infinity
const double one_minus_coordinate = 1.0 - coordinate;
const double denominator = computational_radius * one_minus_coordinate;
if (!std::isfinite(denominator) || denominator <= 0.0) {
return MappingStatus::NON_FINITE_MAPPING_RESULT;
}
// The scale here is the factor which stretches the finite computational domain into the infinite physical domain. Properties we need this to have
// include that it should go to 1 at the stellar surface and go to infinity at the compactified infinity.
// Mathematically this is scale = |r|/|x| where r is the physical radius and x is the computational radius.
// Put another way, scale is the ratio of the target physical radius for the current compactification coordinate
// to the current mesh radius.
const double scale = m_options.r_star_ref / denominator;
const double scale_derivative = scale * (1.0 / one_minus_coordinate - radial_extent / computational_radius);
if (!std::isfinite(scale) || !std::isfinite(scale_derivative)) {
return MappingStatus::NON_FINITE_MAPPING_RESULT;
}
radial_factors.compactification_coordinate = coordinate;
radial_factors.computational_radius = computational_radius;
radial_factors.scale = scale;
radial_factors.scale_derivative = scale_derivative;
return MappingStatus::VALID;
}
}

153
src/lib/serif/mapping/data.cpp
View File

@@ -1,153 +0,0 @@
#include <stdexcept>
#include "mfem/linalg/vector.hpp"
#include "mfem/fem/fe/fe_base.hpp"
#include "mfem/fem/fespace.hpp"
#include "serif/mapping/data.hpp"
#include "serif/utils/error/errors.hpp"
namespace serif::mapping {
using namespace serif::utils::errors;
using mfem::Ordering;
ElementDisplacementData::ElementDisplacementData(
const mfem::FiniteElement& element,
mfem::Vector& displacement_dofs,
const Ordering::Type ordering
) : m_element(&element), m_dimension(0), m_ordering(ordering) {
const int dof_count = element.GetDof();
if (dof_count <= 0) {
serif_error<std::invalid_argument, SERiFErrorCode::INVALID_DOF_COUNT>("ElementDisplacementData: Displacement elements must have at least one degree of freedom. Got an element with {} dofs", dof_count);
}
if (displacement_dofs.Size() <= 0 || displacement_dofs.Size() % dof_count != 0) {
serif_error<std::invalid_argument, SERiFErrorCode::INVALID_DOF_COUNT>("ElementDisplacementData: The displacement vector size ({}) must be a positive multiple of the element degree of freedom count ({}).", displacement_dofs.Size(), dof_count);
}
m_dimension = displacement_dofs.Size() / dof_count;
m_dof_matrix.SetSize(dof_count, m_dimension);
switch (ordering) {
case Ordering::byNODES: {
for (int componentID = 0; componentID < m_dimension; ++componentID) {
for (int dofID = 0; dofID < dof_count; ++dofID) {
m_dof_matrix(dofID, componentID) = displacement_dofs(dofID + componentID * dof_count);
}
}
break;
}
case Ordering::byVDIM: {
for (int dofID = 0; dofID < dof_count; ++dofID) {
for (int componentID = 0; componentID < m_dimension; ++componentID) {
m_dof_matrix(dofID, componentID) = displacement_dofs(componentID + dofID * m_dimension);
}
}
break;
}
default:
serif_error<std::invalid_argument, SERiFErrorCode::INVALID_DOF_ORDERING>("Unsupported ordering");
}
}
const mfem::FiniteElement& ElementDisplacementData::get_element() const noexcept {
return *m_element;
}
const mfem::DenseMatrix& ElementDisplacementData::get_dof_matrix() const noexcept {
return m_dof_matrix;
}
int ElementDisplacementData::get_dimension() const noexcept {
return m_dimension;
}
int ElementDisplacementData::get_dof_count() const noexcept {
return m_element->GetDof();
}
Ordering::Type ElementDisplacementData::get_ordering() const noexcept {
return m_ordering;
}
ElementDisplacementData element_displacement_data_from_element_vdofs(const mfem::FiniteElement& element, mfem::Vector& displacement_dofs) {
return ElementDisplacementData(element, displacement_dofs, Ordering::byNODES);
}
ElementCompactificationData::ElementCompactificationData(
const mfem::FiniteElement& element,
const mfem::Vector& dofs
) : m_element(&element), m_dofs(dofs) {
if (element.GetRangeType() != mfem::FiniteElement::SCALAR) {
serif_error<std::invalid_argument, SERiFErrorCode::INVALID_COMPACTIFICATION_RANK>("Compactification coordinate requires a scalar finite element.");
}
if (element.GetMapType() != mfem::FiniteElement::VALUE) {
serif_error<std::invalid_argument, SERiFErrorCode::INVALID_COMPACTIFICATION_TYPE>("Compactification coordinate requires a value-mapped scalar finite element");
}
if (element.GetDerivType() != mfem::FiniteElement::GRAD) {
serif_error<std::invalid_argument, SERiFErrorCode::INVALID_COMPACTIFICATION_DERIVATIVE_TYPE>("Compactification coordinate finite element must provide a gradient");
}
if (element.GetDof() <= 0) {
serif_error<std::invalid_argument, SERiFErrorCode::INVALID_DOF_COUNT>("Compactification coordinatge finite element must have a positive non-zero number of degrees of freedom. Recived {}", element.GetDof());
}
if (dofs.Size() != element.GetDof()) {
serif_error<std::invalid_argument, SERiFErrorCode::INVALID_DOF_COUNT>("Compactification coordinate degree of freedom count ({}) does not match its finite element dof count ({})", dofs.Size(), element.GetDof());
}
}
const mfem::FiniteElement& ElementCompactificationData::get_element() const noexcept {
return *m_element;
}
const mfem::Vector& ElementCompactificationData::get_dofs() const noexcept {
return m_dofs;
}
int ElementCompactificationData::get_dof_count() const noexcept {
return m_dofs.Size();
}
DMWorkspace::DMWorkspace(const int dimension) {
set_dimension(dimension);
}
void DMWorkspace::set_dimension(const int dimension) {
if (dimension <= 0) {
serif_error<std::invalid_argument, SERiFErrorCode::INVALID_DIMENSION>("Domain mapping workspace dimension must be positive and finite. Received {}", dimension);
}
m_dimension = dimension;
m_field_value.SetSize(dimension);
m_field_jacobian.SetSize(dimension, dimension);
m_reference_field_jacobian.SetSize(dimension, dimension);
m_compactification_point.compactification_coordinate = 0.0;
m_compactification_point.coordinate_gradient.SetSize(dimension);
m_reference_normal.SetSize(dimension);
m_mapped_normal.SetSize(dimension);
m_full_element_jacobian.SetSize(dimension, dimension);
m_vector_temp.SetSize(dimension);
m_matrix_temp_1.SetSize(dimension, dimension);
m_matrix_temp_2.SetSize(dimension, dimension);
m_exterior_result.physical_position.SetSize(dimension);
m_exterior_result.mapping_jacobian.SetSize(dimension, dimension);
m_exterior_variation.physical_position_variation.SetSize(dimension);
m_exterior_variation.mapping_jacobian_variation.SetSize(dimension, dimension);
}
int DMWorkspace::get_dimension() const noexcept {
return m_dimension;
}
}

730
src/lib/serif/mapping/domain_mapper.cpp
View File

@@ -1,730 +0,0 @@
#include "mfem/linalg/densemat.hpp"
#include "serif/mapping/domain_mapper.hpp"
#include "serif/utils/error/errors.hpp"
#include "serif/utils/checks/finite.hpp"
#include <array>
#include <tuple>
// Anonymous namespace for translation unit-local helper functions and types
namespace {
// TODO: Put something in here using variadic args to validate multiple things satisfy some predicate
}
namespace serif::mapping {
using namespace serif::utils::errors;
using namespace serif::utils::checks;
DomainMapper::DomainMapper(
utils::user_options::DomainMapperOptions options,
std::unique_ptr<const compactification::ExteriorDomainMap> exterior_map
) :
m_geometryCache(std::make_unique<DMGeometryCache>()),
m_options(options),
m_exterior_map(std::move(exterior_map)) {
if (m_options.dimension <= 0) {
serif_error<std::invalid_argument, SERiFErrorCode::INVALID_DIMENSION>("Domain Mapper requires a positive, non-zero, dimension. Received dimension: {}", m_options.dimension);
}
if (m_options.vacuum_element_attribute_ID <= 0) {
serif_error<std::invalid_argument, SERiFErrorCode::INVALID_VACUUM_ELEMENT_ATTRIBUTE_ID>("Domain Mapper requires a positive, non-zero, vacuum element attribute ID. Received attribute ID: {}", m_options.vacuum_element_attribute_ID);
}
if (!m_exterior_map) {
serif_error<std::invalid_argument, SERiFErrorCode::MISSING_EXTERIOR_MAP>("Domain Mapper requires a valid exterior map. One was not provided or a nullptr was provided.");
}
}
std::unique_ptr<DomainMapper> DomainMapper::clone() const {
if (!m_exterior_map) {
serif_error<std::logic_error, SERiFErrorCode::MISSING_EXTERIOR_MAP>("Domain Mapper clone failed because the exterior map is missing or a nullptr.");
}
auto result = std::make_unique<DomainMapper>(m_options, m_exterior_map->clone());
result->set_geometry_cache_enabled(is_geometry_cache_enabled());
return result;
}
void DomainMapper::set_geometry_cache_enabled(const bool enabled) const {
m_geometryCache->enabled = enabled;
}
bool DomainMapper::is_geometry_cache_enabled() const noexcept {
return m_geometryCache->enabled;
}
void DomainMapper::clear_geometry_cache() const {
auto &cache = *m_geometryCache;
cache.points.clear();
cache.mesh = nullptr;
cache.sequence = cache.nodesSequence = -1;
cache.hits = cache.misses = 0;
}
DMGeometryCacheStatistics DomainMapper::get_geometry_cache_statistics() const {
const auto& cache = *m_geometryCache;
return {.hits = cache.hits, .misses = cache.misses, .points = cache.points.size()};
}
const DMFixedMeshGeometry* DomainMapper::get_fixed_mesh_geometry(
mfem::ElementTransformation& transformation,
const mfem::IntegrationPoint& ip
) const {
transformation.SetIntPoint(&ip);
auto &cache = *m_geometryCache;
const auto *mesh = transformation.mesh;
auto should_return_null = [&]() -> bool {
const bool is_cache_disabled = !cache.enabled;
const bool is_mesh_null = mesh == nullptr;
const bool is_not_element = transformation.ElementType == mfem::ElementTransformation::ELEMENT;
const bool is_element_no_lt_0 = transformation.ElementNo < 0;
const bool is_element_no_gt_NE = transformation.ElementNo >= mesh->GetNE();
const bool is_dim_invalid = transformation.GetSpaceDim() != m_options.dimension;
return (is_cache_disabled || is_mesh_null || is_not_element || is_element_no_lt_0 || is_element_no_gt_NE || is_dim_invalid);
};
if (should_return_null()) {
return nullptr;
}
const std::array<double, 3> coordinates{ip.x, m_options.dimension > 1 ? ip.y : 0.0, m_options.dimension > 2 ? ip.z: 0.0};
for (const double coordinate : coordinates) {
if (!std::isfinite(coordinate)) {
return nullptr;
}
}
if (cache.mesh != mesh || cache.sequence != mesh->GetSequence() || cache.nodesSequence != mesh->GetNodesSequence()) {
cache.points.clear();
cache.mesh = mesh;
cache.sequence = mesh->GetSequence();
cache.nodesSequence = mesh->GetNodesSequence();
}
const auto key = std::tuple{transformation.ElementNo, coordinates[0], coordinates[1], coordinates[2]};
const auto iter = cache.points.find(key);
if (iter != cache.points.end()) {
++cache.hits;
return &iter->second;
}
DMFixedMeshGeometry geometry;
transformation.Transform(ip, geometry.position);
geometry.jacobian = transformation.Jacobian();
geometry.inverse_jacobian = transformation.InverseJacobian();
geometry.weight = transformation.Weight();
++cache.misses;
if (!is_finite(geometry.position) || !is_finite(geometry.jacobian) || !is_finite(geometry.inverse_jacobian) || !is_finite(geometry.weight)) {
return nullptr;
}
cache.points.emplace(key, std::move(geometry));
return &cache.points.at(key);
}
MappingStatus DomainMapper::evaluate_point(
const ElementMappingData &element_data,
mfem::ElementTransformation &transformation,
const mfem::IntegrationPoint &ip,
DMWorkspace &workspace,
MappingPointContext &context
) const {
const auto* meshGeometry = get_fixed_mesh_geometry(transformation, ip);
return evaluate_point_impl(element_data, transformation, ip, workspace, context, meshGeometry);
}
MappingStatus DomainMapper::evaluate_volume(
const ElementMappingData &element_data,
mfem::ElementTransformation &transformation,
const mfem::IntegrationPoint &ip,
DMWorkspace &workspace,
VolumeMappingContext &context,
const DMReferenceDisplacementSample *displacement_sample
) const {
const auto *meshGeometry = get_fixed_mesh_geometry(transformation, ip);
const MappingStatus point_status = evaluate_point_impl(
element_data,
transformation,
ip,
workspace,
context.mapping,
meshGeometry,
displacement_sample
);
if (point_status != MappingStatus::VALID) {
return point_status;
}
mfem::Mult(context.mapping.mapping_jacobian,meshGeometry ? meshGeometry->jacobian : transformation.Jacobian(), workspace.m_full_element_jacobian);
context.quadrature.J_element_to_physical_inverse.SetSize(m_options.dimension, m_options.dimension);
mfem::CalcInverse(workspace.m_full_element_jacobian, context.quadrature.J_element_to_physical_inverse);
context.quadrature.mapping_determinate = context.mapping.mapping_determinate;
context.quadrature.physical_integration_weight = ip.weight * (meshGeometry ? meshGeometry->weight : transformation.Weight()) * context.mapping.mapping_determinate;
if (!is_finite(context.quadrature.J_element_to_physical_inverse) || !std::isfinite(context.quadrature.physical_integration_weight)) {
return MappingStatus::NON_FINITE_MAPPING_RESULT;
}
if (context.quadrature.physical_integration_weight <= 0.0) {
return MappingStatus::NON_POSITIVE_DETERMINATE;
}
return MappingStatus::VALID;
}
MappingStatus DomainMapper::evaluate_face(
const ElementMappingData &element_data,
mfem::FaceElementTransformations &transformation,
const FaceElementSide side,
const mfem::IntegrationPoint &ip,
DMWorkspace &workspace,
FaceMappingContext &context
) const {
transformation.SetAllIntPoints(&ip);
mfem::ElementTransformation& element_transformation = select_face_element_transformation(transformation, side);
const mfem::IntegrationPoint& element_integration_point = select_face_element_integration_point(transformation, side);
const MappingStatus point_status = evaluate_point(element_data, element_transformation, element_integration_point, workspace, context.mapping);
if (point_status != MappingStatus::VALID) {
return point_status;
}
// Get the normal vector from the face in reference space
workspace.m_reference_normal.SetSize(m_options.dimension);
mfem::CalcOrtho(transformation.Jacobian(), workspace.m_reference_normal);
if (side == FaceElementSide::ELEMENT_2) {
workspace.m_reference_normal *= -1.0; // Flip the normal for the second element (e.g. the "other" side of the face)
}
const double reference_normal_magnitude = workspace.m_reference_normal.Norml2();
if (!is_finite(reference_normal_magnitude) || reference_normal_magnitude <= 0.0) {
return MappingStatus::NON_FINITE_MAPPING_RESULT;
}
context.reference_normal.SetSize(m_options.dimension);
context.reference_normal = workspace.m_reference_normal;
context.reference_normal /= reference_normal_magnitude;
// Get the normal vector from the face in physical space
context.mapping.inverse_mapping_jacobian.MultTranspose(workspace.m_reference_normal, workspace.m_mapped_normal);
workspace.m_mapped_normal *= context.mapping.mapping_determinate;
const double mapped_normal_magnitude = workspace.m_mapped_normal.Norml2();
if (!is_finite(mapped_normal_magnitude) || mapped_normal_magnitude <= 0.0) {
return MappingStatus::NON_FINITE_MAPPING_RESULT;
}
context.quadrature.normal.SetSize(m_options.dimension);
context.quadrature.normal = workspace.m_mapped_normal;
context.quadrature.normal /= mapped_normal_magnitude;
context.reference_surface_weight = ip.weight * reference_normal_magnitude;
context.physical_surface_weight = ip.weight * mapped_normal_magnitude;
context.quadrature.ds = context.reference_surface_weight;
context.quadrature.v_dot_n_scale = mapped_normal_magnitude / reference_normal_magnitude;
if (!is_finite(context.quadrature.normal) || !is_finite(context.reference_surface_weight) || !is_finite(context.physical_surface_weight) || !is_finite(context.quadrature.v_dot_n_scale)) {
return MappingStatus::NON_FINITE_MAPPING_RESULT;
}
return MappingStatus::VALID;
}
MappingStatus DomainMapper::evaluate_point_variation(
const ElementMappingData &element_data,
const ElementDisplacementData &direction,
mfem::ElementTransformation &transformation,
const mfem::IntegrationPoint &ip,
const MappingPointContext &context,
DMWorkspace &workspace,
MappingPointVariation &variation
) const {
return evaluate_point_variation_impl(element_data, direction, transformation, ip, context, workspace, variation, nullptr);
}
MappingStatus DomainMapper::evaluate_volume_variation(
const ElementMappingData &element_data,
const ElementDisplacementData &direction,
mfem::ElementTransformation &transformation,
const mfem::IntegrationPoint &ip,
const VolumeMappingContext &context,
DMWorkspace &workspace,
VolumeMappingVariation &variation
) const {
mfem::Mult(context.quadrature.J_element_to_physical_inverse, context.mapping.mapping_jacobian, workspace.m_matrix_temp_2);
const MappingStatus point_status = evaluate_point_variation_impl(element_data, direction, transformation, ip, context.mapping, workspace, variation.mapping, &workspace.m_matrix_temp_2);
if (point_status != MappingStatus::VALID) {
return point_status;
}
mfem::Mult(context.quadrature.J_element_to_physical_inverse, variation.mapping.mapping_jacobian_variation, workspace.m_matrix_temp_1);
variation.inverse_element_jacobian_variation.SetSize(m_options.dimension, m_options.dimension);
mfem::Mult(workspace.m_matrix_temp_1, context.mapping.inverse_mapping_jacobian, variation.inverse_element_jacobian_variation);
variation.inverse_element_jacobian_variation *= -1.0; // Account for the negative sign in the derivative of the inverse of a matrix
variation.weight_variation = (context.quadrature.physical_integration_weight / context.mapping.mapping_determinate) * variation.mapping.mapping_determinate_variation;
if (!is_finite(variation.inverse_element_jacobian_variation) || !is_finite(variation.weight_variation)) {
return MappingStatus::NON_FINITE_MAPPING_RESULT;
}
return MappingStatus::VALID;
}
MappingStatus DomainMapper::evaluate_face_variation(
const ElementMappingData &element_data,
const ElementDisplacementData &direction,
mfem::FaceElementTransformations &transformation,
FaceElementSide side,
const mfem::IntegrationPoint &ip,
const FaceMappingContext &context,
DMWorkspace &workspace,
FaceMappingVariation &variation
) const {
transformation.SetAllIntPoints(&ip);
mfem::ElementTransformation &element_transformation = select_face_element_transformation(transformation, side);
const mfem::IntegrationPoint &element_integration_point = select_face_element_integration_point(transformation, side);
const MappingStatus point_status = evaluate_point_variation(element_data, direction, element_transformation, element_integration_point, context.mapping, workspace, variation.mapping);
if (point_status != MappingStatus::VALID) {
return point_status;
}
workspace.m_reference_normal.SetSize(m_options.dimension);
mfem::CalcOrtho(transformation.Jacobian(), workspace.m_reference_normal);
if (side == FaceElementSide::ELEMENT_2) {
workspace.m_reference_normal *= -1.0; // Flip the normal for the second element (e.g. the "other" side of the face)
}
const double reference_normal_magnitude = workspace.m_reference_normal.Norml2();
if (!is_finite(reference_normal_magnitude) || reference_normal_magnitude <= 0) {
return MappingStatus::NON_FINITE_MAPPING_RESULT;
}
context.mapping.inverse_mapping_jacobian.MultTranspose(workspace.m_reference_normal, workspace.m_vector_temp);
workspace.m_mapped_normal = workspace.m_vector_temp;
workspace.m_mapped_normal *= context.mapping.mapping_determinate;
variation.physical_normal_variation.SetSize(m_options.dimension);
variation.mapping.inverse_mapping_jacobian_variation.MultTranspose(workspace.m_reference_normal, variation.physical_normal_variation);
variation.physical_normal_variation *= context.mapping.mapping_determinate;
variation.physical_normal_variation.Add(variation.mapping.mapping_determinate_variation, workspace.m_vector_temp);
const double mapped_normal_magnitude = workspace.m_mapped_normal.Norml2();
if (!is_finite(mapped_normal_magnitude) || mapped_normal_magnitude <= 0) {
return MappingStatus::NON_FINITE_MAPPING_RESULT;
}
const double mapped_normal_magnitude_variation = context.quadrature.normal * variation.physical_normal_variation;
variation.physical_normal_variation.Add(-mapped_normal_magnitude_variation, context.quadrature.normal);
variation.physical_normal_variation /= mapped_normal_magnitude;
variation.physical_surface_weight_variation = ip.weight * mapped_normal_magnitude_variation;
variation.normal_flux_scale_variation = mapped_normal_magnitude_variation / reference_normal_magnitude;
if (!is_finite(variation.physical_normal_variation) || !std::isfinite(variation.physical_surface_weight_variation) || !std::isfinite(variation.normal_flux_scale_variation)) {
return MappingStatus::NON_FINITE_MAPPING_RESULT;
}
return MappingStatus::VALID;
}
bool DomainMapper::is_compactified_element(const mfem::ElementTransformation &transformation) const noexcept {
return transformation.Attribute == m_options.vacuum_element_attribute_ID;
}
int DomainMapper::get_dimension() const noexcept {
return m_options.dimension;
}
const compactification::ExteriorDomainMap& DomainMapper::get_exterior_map() const noexcept {
return *m_exterior_map;
}
MappingStatus DomainMapper::evaluate_point_impl(
const ElementMappingData &element_data,
mfem::ElementTransformation &transformation,
const mfem::IntegrationPoint &ip,
DMWorkspace &workspace,
MappingPointContext &context,
const DMFixedMeshGeometry *geometry,
const DMReferenceDisplacementSample *displacement_sample
) const {
validate_element_data(element_data);
//TODO: Factor into its own function
//TODO: Consider if these should return or throw.
if (workspace.get_dimension() != m_options.dimension) {
serif_error<std::invalid_argument, SERiFErrorCode::INVALID_DIMENSION>("Workspace dimension does not match Domain Mapper dimension. Workspace dimension: {}, Domain Mapper dimension: {}", workspace.get_dimension(), m_options.dimension);
}
if (transformation.GetSpaceDim() != m_options.dimension) {
serif_error<std::invalid_argument, SERiFErrorCode::INVALID_DIMENSION>("Element transformation space dimension does not match Domain Mapper dimension. Element transformation space dimension: {}, Domain Mapper dimension: {}", transformation.GetSpaceDim(), m_options.dimension);
}
if (transformation.GetGeometryType() != element_data.displacement.get_element().GetGeomType()) {
serif_error<std::invalid_argument, SERiFErrorCode::INVALID_ELEMENT_GEOMETRY>("Element transformation geometry type does not match element geometry type. Element transformation geometry type: {}, Element geometry type: {}", transformation.GetGeometryType(), element_data.displacement.get_element().GetGeomType());
}
context.reference_position.SetSize(m_options.dimension);
if (geometry) {
context.reference_position = geometry->position;
} else {
transformation.Transform(ip, context.reference_position);
}
// Get the displacement field value and its Jacobian at the integration point. Note these may be in current workspace to avoid repeated allocations and Multiplies.
if (displacement_sample) {
// TODO: Generally go through and refactor these large checks into their own functions, at the very least it will be easier to read
if (element_data.displacement.get_element().GetMapType() != mfem::FiniteElement::VALUE || displacement_sample->value.Size() != m_options.dimension || displacement_sample->gradient.Height() != m_options.dimension || displacement_sample->gradient.Width() != m_options.dimension) {
return MappingStatus::INVALID_DIMENSION;
}
workspace.m_field_value = displacement_sample->value;
workspace.m_field_jacobian.SetSize(m_options.dimension);
mfem::Mult(displacement_sample->gradient, geometry ? geometry->inverse_jacobian : transformation.InverseJacobian(), workspace.m_field_jacobian);
} else {
evaluate_displacement_field(element_data.displacement, transformation, ip, workspace, workspace.m_field_value, workspace.m_field_jacobian, geometry ? &geometry->inverse_jacobian : nullptr);
}
if (!is_finite(context.reference_position) || !is_finite(workspace.m_field_value) || !is_finite(workspace.m_field_jacobian)) {
return MappingStatus::NON_FINITE_MAPPING_RESULT;
}
context.displaced_position.SetSize(m_options.dimension);
// Get the position of the point in physical space by adding the displacement to the reference position. Note MFEM really dislikes raw arithmetic operators
// so we need to first assign the reference position then use the in place += operator.
context.displaced_position = context.reference_position;
context.displaced_position += workspace.m_field_value;
context.displacement_jacobian.SetSize(m_options.dimension, m_options.dimension);
context.displacement_jacobian = workspace.m_field_jacobian;
// Ensure that the diagonal of the displacement Jacobian is incremented by 1.0 to account for the identity mapping from reference to physical space.
// recall that r = x + d (where d is the workspace.m_field_value and x is context.reference_position) then we can differentiate this
// component wise to find the gradient of the displaced position wrt. the mesh coordinate (reference position). E.g as you move along
// the mesh coordinate how much does the physical coordinate change and in what direction. Lets call this F
// F = \frac{\partial r_i}{\partial x_j} where r is the displaced position and x is the mesh position.
// We then have F = \frac{\partial x_i}{\partial x_j} + \frac{\partial d_i}{x_j} where d is the displacement (recall r = x + d)
// By definition the first term is the identity matrix. The second term we get out of EvaluateField. Thus why we need to add the identity matrix here
for (int i = 0; i < m_options.dimension; ++i)
context.displacement_jacobian(i, i) += 1.0;
context.compactified = is_compactified_element(transformation);
// This branch only runs for vacuum elements
if (context.compactified) {
// There are two things that we need to the mapping. First is a reference coordinate which stroid embeds into the mesh at mesh generation time, this is
// parameterized from 0 - 1 where 0 is the model surface and 1 is the mesh exterior (what will becomes the compactified infinity, note also we never actually evaluate at s=1; rather we define some arbitrary small tolerance to approach s=1). Lets call this s. We also need
// the gradient of s as we move along the mesh coordinates. All of this is stashes within workspace.m_compactification_point.
const MappingStatus coordinate_status = evaluate_compactification_coordinate(element_data.compactification, transformation, ip, workspace, workspace.m_compactification_point, geometry ? &geometry->inverse_jacobian : nullptr);
if (coordinate_status != MappingStatus::VALID) {
return coordinate_status;
}
const compactification::ExteriorMapInput exterior_input{
.reference_position = context.reference_position,
.displaced_position = context.displaced_position,
.displacement_jacobian = context.displacement_jacobian,
.compactification_coordinate_gradient = workspace.m_compactification_point.coordinate_gradient,
.compactification_coordinate = workspace.m_compactification_point.compactification_coordinate,
};
// This apply whatever the exterior map is to generate the new physical exterior coordinate and jacobian between physical and reference space.
// In general we have only implemented a kelvin mapping; however, in future additional mappings may be implemented.
const MappingStatus exterior_status = m_exterior_map->evaluate(exterior_input, workspace.m_exterior_result);
if (exterior_status != MappingStatus::VALID) {
return exterior_status;
}
context.physical_position = workspace.m_exterior_result.physical_position;
context.mapping_jacobian = workspace.m_exterior_result.mapping_jacobian;
} else {
context.physical_position = context.displaced_position;
context.mapping_jacobian = context.displacement_jacobian;
}
// Validation work
if (!is_finite(context.physical_position) || !is_finite(context.mapping_jacobian)) {
return MappingStatus::NON_FINITE_MAPPING_RESULT;
}
context.mapping_determinate = context.mapping_jacobian.Det();
if (!std::isfinite(context.mapping_determinate)) {
return MappingStatus::NON_FINITE_MAPPING_RESULT;
}
if (context.mapping_determinate <= 0.0) {
// This is the most common error we see come out of this function, specifically it is common when we try to deform the mesh too much in one step.
return MappingStatus::NON_POSITIVE_DETERMINATE;
}
context.inverse_mapping_jacobian.SetSize(m_options.dimension, m_options.dimension);
// Finally we get the inverse jacobian.
mfem::CalcInverse(context.mapping_jacobian, context.inverse_mapping_jacobian);
if (!is_finite(context.inverse_mapping_jacobian)) {
return MappingStatus::NON_FINITE_MAPPING_RESULT;
}
return MappingStatus::VALID;
}
void DomainMapper::validate_element_data(const ElementMappingData &element_data) const {
const ElementDisplacementData &displacement = element_data.displacement;
const ElementCompactificationData &compactification = element_data.compactification;
if (displacement.get_dimension() != m_options.dimension) {
serif_error<std::invalid_argument, SERiFErrorCode::INVALID_DIMENSION>("Element displacement dimension does not match Domain Mapper dimension. Element displacement dimension: {}, Domain Mapper dimension: {}", displacement.get_dimension(), m_options.dimension);
}
if (displacement.get_element().GetDim() != m_options.dimension) {
serif_error<std::invalid_argument, SERiFErrorCode::INVALID_DIMENSION>("Element displacement finite element dimension does not match Domain Mapper dimension. Element displacement finite element dimension: {}, Domain Mapper dimension: {}", displacement.get_element().GetDim(), m_options.dimension);
}
if (compactification.get_element().GetDim() != m_options.dimension) {
serif_error<std::invalid_argument, SERiFErrorCode::INVALID_DIMENSION>("Element compactification finite element dimension does not match Domain Mapper dimension. Element compactification finite element dimension: {}, Domain Mapper dimension: {}", compactification.get_element().GetDim(), m_options.dimension);
}
if (displacement.get_element().GetGeomType() != compactification.get_element().GetGeomType()) {
serif_error<std::invalid_argument, SERiFErrorCode::INVALID_ELEMENT_GEOMETRY>("Element displacement finite element geometry type does not match element compactification finite element geometry type. Element displacement finite element geometry type: {}, Element compactification finite element geometry type: {}", displacement.get_element().GetGeomType(), compactification.get_element().GetGeomType());
}
if (compactification.get_element().GetRangeType() != mfem::FiniteElement::SCALAR) {
serif_error<std::invalid_argument, SERiFErrorCode::INVALID_COMPACTIFICATION_RANK>("Element compactification finite element range type is not scalar. Element compactification finite element range type: {}", compactification.get_element().GetRangeType());
}
if (compactification.get_element().GetMapType() != mfem::FiniteElement::VALUE) {
serif_error<std::invalid_argument, SERiFErrorCode::INVALID_COMPACTIFICATION_TYPE>("Element compactification finite element map type is not value. Element compactification finite element map type: {}", compactification.get_element().GetMapType());
}
if (compactification.get_element().GetDerivType() != mfem::FiniteElement::GRAD) {
serif_error<std::invalid_argument, SERiFErrorCode::INVALID_COMPACTIFICATION_DERIVATIVE_TYPE>("Element compactification finite element derivative type is not gradient. Element compactification finite element derivative type: {}", compactification.get_element().GetDerivType());
}
if (compactification.get_dof_count() != compactification.get_element().GetDof()) {
serif_error<std::invalid_argument, SERiFErrorCode::INVALID_DOF_COUNT>("Element compactification finite element DOF count does not match the number of DOFs in the finite element. Element compactification finite element DOF count: {}, Finite element DOF count: {}", compactification.get_dof_count(), compactification.get_element().GetDof());
}
}
// Note that there is nothing preventing this function from being used with some field other than displacement; perhaps we should tighten this?
void DomainMapper::evaluate_displacement_field(
const ElementDisplacementData &field,
mfem::ElementTransformation &transformation,
const mfem::IntegrationPoint &ip,
DMWorkspace &workspace,
mfem::Vector &value,
mfem::DenseMatrix &jacobian,
const mfem::DenseMatrix *inverse_mesh_jacobian
) const {
const mfem::FiniteElement &element = field.get_element();
const mfem::DenseMatrix &dof_matrix = field.get_dof_matrix();
workspace.m_shape.SetSize(element.GetDof());
element.CalcShape(ip, workspace.m_shape);
value.SetSize(m_options.dimension);
dof_matrix.MultTranspose(workspace.m_shape, value);
jacobian.SetSize(m_options.dimension, m_options.dimension);
if (inverse_mesh_jacobian != nullptr || element.GetMapType() == mfem::FiniteElement::VALUE) {
workspace.m_reference_dshape.SetSize(element.GetDof(), m_options.dimension);
element.CalcDShape(ip, workspace.m_reference_dshape);
// Contract DOFs before applying fixed-mesh geometry.
mfem::MultAtB(dof_matrix, workspace.m_reference_dshape, workspace.m_reference_field_jacobian);
const mfem::DenseMatrix &inverseMeshJacobian = inverse_mesh_jacobian != nullptr ? *inverse_mesh_jacobian : transformation.InverseJacobian();
mfem::Mult(workspace.m_reference_field_jacobian, inverseMeshJacobian, jacobian);
} else {
// Retain the original finite-element-specific physical-gradient
// path for mapping types without the ordinary VALUE pullback.
workspace.m_mesh_dshape.SetSize(element.GetDof(), m_options.dimension);
element.CalcPhysDShape(transformation, workspace.m_mesh_dshape);
mfem::MultAtB(dof_matrix, workspace.m_mesh_dshape, jacobian);
}
}
MappingStatus DomainMapper::evaluate_compactification_coordinate(
const ElementCompactificationData &compactification,
mfem::ElementTransformation &transformation,
const mfem::IntegrationPoint &ip,
DMWorkspace &workspace,
CompactificationPointData &point_data,
const mfem::DenseMatrix *inverse_mesh_jacobian
) const {
const mfem::FiniteElement& element = compactification.get_element();
const mfem::Vector& dofs = compactification.get_dofs();
const int dof_count = element.GetDof();
if (workspace.get_dimension() != m_options.dimension || transformation.GetSpaceDim() != m_options.dimension || element.GetDim() != m_options.dimension) {
return MappingStatus::INVALID_DIMENSION;
}
if (dofs.Size() != dof_count) {
return MappingStatus::INVALID_DIMENSION;
}
for (int i = 0; i < dof_count; ++i) {
if (!is_finite(dofs(i))) {
return MappingStatus::NON_FINITE_MAPPING_INPUT;
}
}
workspace.m_compactification_shape.SetSize(dof_count);
workspace.m_compactification_dshape.SetSize(dof_count, m_options.dimension);
element.CalcShape(ip, workspace.m_compactification_shape);
if (inverse_mesh_jacobian != nullptr) {
workspace.m_reference_dshape.SetSize(dof_count, m_options.dimension);
element.CalcDShape(ip, workspace.m_reference_dshape);
mfem::Mult(workspace.m_reference_dshape, *inverse_mesh_jacobian, workspace.m_compactification_dshape);
} else {
element.CalcPhysDShape(transformation, workspace.m_compactification_dshape);
}
point_data.compactification_coordinate = dofs * workspace.m_compactification_shape;
point_data.coordinate_gradient.SetSize(m_options.dimension);
workspace.m_compactification_dshape.MultTranspose(dofs, point_data.coordinate_gradient);
if (!is_finite(point_data.compactification_coordinate)) {
return MappingStatus::NON_FINITE_MAPPING_RESULT;
}
for (int d = 0; d < point_data.coordinate_gradient.Size(); ++d) {
if (!is_finite(point_data.coordinate_gradient(d))) {
return MappingStatus::NON_FINITE_MAPPING_RESULT;
}
}
return MappingStatus::VALID;
}
MappingStatus DomainMapper::evaluate_point_variation_impl(
const ElementMappingData &element_data,
const ElementDisplacementData &direction,
mfem::ElementTransformation &transformation,
const mfem::IntegrationPoint &ip,
const MappingPointContext &context,
DMWorkspace &workspace,
MappingPointVariation &variation,
const mfem::DenseMatrix *inverse_mesh_jacobian
) const {
validate_element_data(element_data);
const ElementMappingData direction_data {
.displacement = direction,
.compactification = element_data.compactification
};
validate_element_data(direction_data);
if (element_data.displacement.get_dof_count() != direction.get_dof_count()) {
serif_error<std::invalid_argument, SERiFErrorCode::INVALID_DOF_COUNT>("Element displacement DOF count does not match direction DOF count. Element displacement DOF count: {}, Direction DOF count: {}", element_data.displacement.get_dof_count(), direction.get_dof_count());
}
if (workspace.get_dimension() != m_options.dimension) {
serif_error<std::invalid_argument, SERiFErrorCode::INVALID_DIMENSION>("Workspace dimension does not match Domain Mapper dimension. Workspace dimension: {}, Domain Mapper dimension: {}", workspace.get_dimension(), m_options.dimension);
}
if (context.compactified != is_compactified_element(transformation)) {
serif_error<std::invalid_argument, SERiFErrorCode::COMPACTIFICATION_MISMATCH>("Mapping context compactification status does not match element transformation compactification status. Mapping context compactified: {}, Element transformation compactified: {}", context.compactified, is_compactified_element(transformation));
}
if (inverse_mesh_jacobian == nullptr) {
if (const auto* geometry = get_fixed_mesh_geometry(transformation, ip)) {
inverse_mesh_jacobian = &geometry->inverse_jacobian;
}
}
evaluate_displacement_field(direction, transformation, ip, workspace, workspace.m_field_value, workspace.m_field_jacobian, inverse_mesh_jacobian);
if (!is_finite(workspace.m_field_jacobian) || !is_finite(workspace.m_field_value)) {
return MappingStatus::NON_FINITE_MAPPING_INPUT;
}
variation.displacement_variation = workspace.m_field_value;
variation.displacement_jacobian_variation = workspace.m_field_jacobian;
if (context.compactified) {
const MappingStatus coordinate_status = evaluate_compactification_coordinate(element_data.compactification, transformation, ip, workspace, workspace.m_compactification_point, inverse_mesh_jacobian);
if (coordinate_status != MappingStatus::VALID) {
return coordinate_status;
}
const compactification::ExteriorMapInput exterior_input{
.reference_position = context.reference_position,
.displaced_position = context.displaced_position,
.displacement_jacobian = context.displacement_jacobian,
.compactification_coordinate_gradient = workspace.m_compactification_point.coordinate_gradient,
.compactification_coordinate = workspace.m_compactification_point.compactification_coordinate,
};
workspace.m_exterior_result.physical_position = context.physical_position;
workspace.m_exterior_result.mapping_jacobian = context.mapping_jacobian;
const compactification::ExteriorMapDirection exterior_direction {
.displaced_position_variation = variation.displacement_variation,
.displacement_jacobian_variation = variation.displacement_jacobian_variation
};
const MappingStatus exterior_status = m_exterior_map->evaluate_variation(exterior_input, workspace.m_exterior_result, exterior_direction, workspace.m_exterior_variation);
if (exterior_status != MappingStatus::VALID) {
return exterior_status;
}
variation.physical_position_variation = workspace.m_exterior_variation.physical_position_variation;
variation.mapping_jacobian_variation = workspace.m_exterior_variation.mapping_jacobian_variation;
} else {
variation.physical_position_variation = variation.displacement_variation;
variation.mapping_jacobian_variation = variation.displacement_jacobian_variation;
}
mfem::Mult(context.inverse_mapping_jacobian, variation.mapping_jacobian_variation, workspace.m_matrix_temp_1);
double trace = 0.0;
for (int dimID = 0; dimID < m_options.dimension; ++dimID) {
trace += workspace.m_matrix_temp_1(dimID, dimID);
}
variation.mapping_determinate_variation = context.mapping_determinate * trace;
variation.inverse_mapping_jacobian_variation.SetSize(m_options.dimension, m_options.dimension);
mfem::Mult(workspace.m_matrix_temp_1, context.inverse_mapping_jacobian, variation.inverse_mapping_jacobian_variation);
variation.inverse_mapping_jacobian_variation *= -1.0;
if (!is_finite(variation.physical_position_variation) || !is_finite(variation.mapping_jacobian_variation) || !std::isfinite(variation.mapping_determinate_variation) || !is_finite(variation.inverse_mapping_jacobian_variation)) {
return MappingStatus::NON_FINITE_MAPPING_RESULT;
}
return MappingStatus::VALID;
}
mfem::ElementTransformation &DomainMapper::select_face_element_transformation(
const mfem::FaceElementTransformations &transformations,
const FaceElementSide side
) {
if (side == FaceElementSide::ELEMENT_1) {
serif_assert<SERiFErrorCode::MISSING_FACE_ELEMENT>(transformations.Elem1 != nullptr, "Face element transformation for ELEMENT_1 is null.");
return *transformations.Elem1;
}
serif_assert<SERiFErrorCode::MISSING_FACE_ELEMENT>(transformations.Elem2 != nullptr, "Face element transformation for ELEMENT_2 is null.");
return *transformations.Elem2;
}
const mfem::IntegrationPoint &DomainMapper::select_face_element_integration_point(
const mfem::FaceElementTransformations &transformations,
const FaceElementSide side
) {
mfem::ElementTransformation& element_transformation = select_face_element_transformation(transformations, side);
return element_transformation.GetIntPoint();
}
}

149
src/lib/serif/mapping/evaluator.cpp
View File

@@ -1,149 +0,0 @@
#include "serif/mapping/evaluator.hpp"
#include "serif/utils/error/errors.hpp"
namespace serif::mapping {
using namespace utils::errors;
GridFunctionMappingEvaluator::GridFunctionMappingEvaluator(
const DomainMapper &mapper,
const mfem::GridFunction &displacement,
const mfem::GridFunction &compactification_coordinate
) :
m_mapper(mapper),
m_displacement(displacement),
m_compactification_coordinate(compactification_coordinate),
m_displacement_space(displacement.FESpace()),
m_compactification_space(compactification_coordinate.FESpace()),
m_displacement_space_sequence(m_displacement_space != nullptr ? m_displacement_space->GetSequence() : -1),
m_compactification_space_sequence(m_compactification_space != nullptr ? m_compactification_space->GetSequence() : -1),
m_workspace(mapper.get_dimension()) {
serif_assert<SERiFErrorCode::MISSING_DISPLACEMENT_FINITE_ELEMENT_SPACE>(m_displacement_space != nullptr, "Displacement grid function is not associated with a finite-element space.");
serif_assert<SERiFErrorCode::MISSING_DISPLACEMENT_FINITE_ELEMENT_SPACE>(m_compactification_space != nullptr, "Compactification-coordinate grid function is not associated with a finite-element space.");
serif_assert<SERiFErrorCode::MISMATCHED_MESHING_MAPS>(m_displacement_space->GetMesh() == m_compactification_space->GetMesh(), "Displacement and compactification-coordinate grid functions are not associated with the same mesh.");
serif_assert<SERiFErrorCode::INVALID_DIMENSION>(m_displacement.VectorDim() == mapper.get_dimension(), "Displacement grid function vector dimension does not match the mapper dimension. ({} vs. {})", m_displacement.VectorDim(), mapper.get_dimension());
serif_assert<SERiFErrorCode::INVALID_COMPACTIFICATION_RANK>(m_compactification_coordinate.VectorDim() == 1, "Compactification-coordinate grid function must be a scalar field. Instead it has vector dimension {}.", m_compactification_coordinate.VectorDim());
serif_assert<SERiFErrorCode::INVALID_DIMENSION>(m_displacement_space->GetMesh()->SpaceDimension() == mapper.get_dimension(), "Displacement finite-element space mesh dimension does not match the mapper dimension. ({} vs. {})", m_displacement_space->GetMesh()->SpaceDimension(), mapper.get_dimension());
}
void GridFunctionMappingEvaluator::invalidate_cache() noexcept {
m_displacement_data.reset();
m_compactification_data.reset();
m_cached_element_id = -1;
}
void GridFunctionMappingEvaluator::validate_field_bindings() const {
serif_assert<SERiFErrorCode::CACHE_FE_BINDING_ERROR>(m_displacement.FESpace() == m_displacement_space, "Displacement grid function is no longer associated with the cached finite-element space.");
serif_assert<SERiFErrorCode::CACHE_FE_BINDING_ERROR>(m_compactification_coordinate.FESpace() == m_compactification_space, "Compactification-coordinate grid function is no longer associated with the cached finite-element space.");
}
bool GridFunctionMappingEvaluator::invalidate_for_changed_spaces() {
const long displacement_sequence = m_displacement_space->GetSequence();
const long compactification_sequence = m_compactification_space->GetSequence();
if (displacement_sequence == m_displacement_space_sequence && compactification_sequence == m_compactification_space_sequence) {
return false;
}
invalidate_cache();
m_displacement_space_sequence = displacement_sequence;
m_compactification_space_sequence = compactification_sequence;
return true;
}
void GridFunctionMappingEvaluator::refresh() {
validate_field_bindings();
if (invalidate_for_changed_spaces()) {
return;
}
const int element_id = m_cached_element_id;
invalidate_cache();
if (element_id >= 0) {
load_element(element_id);
}
}
void GridFunctionMappingEvaluator::load_element(const int element_id) {
validate_field_bindings();
(void)invalidate_for_changed_spaces();
if (element_id == m_cached_element_id) {
return;
}
const mfem::FiniteElementSpace& displacement_space = *m_displacement_space;
const mfem::FiniteElementSpace& compactification_space = *m_compactification_space;
serif_ensure<SERiFErrorCode::INVALID_ELEMENT_ID_ERROR>(element_id >= 0, "Element ID must be non-negative. Instead {} was provided.", element_id);
serif_ensure<SERiFErrorCode::INVALID_ELEMENT_ID_ERROR>(element_id < displacement_space.GetMesh()->GetNE(), "Element ID must be less than the number of elements in the mesh ({}). Instead {} was provided.", displacement_space.GetMesh()->GetNE(), element_id);
mfem::DofTransformation *displacement_transformation = displacement_space.GetElementVDofs(element_id, m_displacement_dofs);
compactification_space.GetElementDofs(element_id, m_compactification_dofs);
m_displacement.GetSubVector(m_displacement_dofs, m_element_displacement);
m_compactification_coordinate.GetSubVector(m_compactification_dofs, m_element_compactification);
if (displacement_transformation != nullptr) {
displacement_transformation->InvTransformPrimal(m_element_displacement);
}
const mfem::FiniteElement &displacement_element = *displacement_space.GetFE(element_id);
const mfem::FiniteElement &compactification_element = *compactification_space.GetFE(element_id);
m_displacement_data = std::make_unique<ElementDisplacementData>(element_displacement_data_from_element_vdofs(displacement_element, m_element_displacement));
m_compactification_data = std::make_unique<ElementCompactificationData>(compactification_element, m_element_compactification);
m_cached_element_id = element_id;
}
MappingStatus GridFunctionMappingEvaluator::evaluate_point(mfem::ElementTransformation &transformation, const mfem::IntegrationPoint &integration_point, MappingPointContext &context) {
load_element(transformation.ElementNo);
const ElementMappingData data{
.displacement = *m_displacement_data,
.compactification = *m_compactification_data
};
return m_mapper.evaluate_point(data, transformation, integration_point, m_workspace, context);
}
MappingStatus GridFunctionMappingEvaluator::evaluate_volume(mfem::ElementTransformation &transformation, const mfem::IntegrationPoint &integration_point, VolumeMappingContext &context) {
load_element(transformation.ElementNo);
const ElementMappingData data{
.displacement = *m_displacement_data,
.compactification = *m_compactification_data
};
return m_mapper.evaluate_volume(data, transformation, integration_point, m_workspace, context);
}
MappingStatus GridFunctionMappingEvaluator::evaluate_face(mfem::FaceElementTransformations &transformations, FaceElementSide side, const mfem::IntegrationPoint &integration_point, FaceMappingContext &context) {
mfem::ElementTransformation *transformation = side == FaceElementSide::ELEMENT_1 ? transformations.Elem1 : transformations.Elem2;
serif_ensure<SERiFErrorCode::ELEMENT_NULLPTR_ERROR>(transformation != nullptr, "Face evaluation requested for a side that does not exist.");
load_element(transformation->ElementNo);
const ElementMappingData data{
.displacement = *m_displacement_data,
.compactification = *m_compactification_data
};
return m_mapper.evaluate_face(data, transformations, side, integration_point, m_workspace, context);
}
VolumeQuadratureContext GridFunctionMappingEvaluator::get_quadrature_context(mfem::ElementTransformation &transformation, const mfem::IntegrationPoint &integration_point) {
VolumeMappingContext context;
serif_ensure<SERiFErrorCode::INVALID_DOMAIN_ERROR>(evaluate_volume(transformation, integration_point, context) == MappingStatus::VALID, "Failed to evaluate the mapping at the given volume quadrature point.");
return context.quadrature;
}
FaceQuadratureContext GridFunctionMappingEvaluator::get_face_quadrature_context(mfem::FaceElementTransformations &transformation, const mfem::IntegrationPoint &integration_point, FaceElementSide side) {
FaceMappingContext context;
serif_ensure<SERiFErrorCode::INVALID_DOMAIN_ERROR>(evaluate_face(transformation, side, integration_point, context) == MappingStatus::VALID, "Failed to evaluate the mapping at the given face quadrature point.");
return context.quadrature;
}
void GridFunctionMappingEvaluator::get_physical_point(mfem::ElementTransformation &transformation, const mfem::IntegrationPoint &integration_point, mfem::Vector &physical_position) {
MappingPointContext context;
serif_ensure<SERiFErrorCode::INVALID_DOMAIN_ERROR>(evaluate_point(transformation, integration_point, context) == MappingStatus::VALID, "Failed to evaluate the mapping at the given point.");
physical_position = context.physical_position;
}
}

136
src/lib/serif/mapping/transformations.cpp
View File

@@ -1,136 +0,0 @@
#include "serif/mapping/transformations.hpp"
#include "serif/utils/error/errors.hpp"
#include "mfem/linalg/vector.hpp"
#include "mfem/linalg/densemat.hpp"
#include "serif/utils/checks/finite.hpp"
namespace serif::mapping {
using namespace serif::utils::errors;
void map_hdiv_flux_to_physical(const MappingPointContext &context, const mfem::Vector &reference_flux, mfem::Vector &physical_flux) {
serif_assert<SERiFErrorCode::INVALID_DIMENSION>(reference_flux.Size() == context.mapping_jacobian.Width(), "Reference flux vector has incompatible dimension with the mapping jacobian ({} <=> {}).", reference_flux.Size(), context.mapping_jacobian.Width());
physical_flux.SetSize(reference_flux.Size());
context.mapping_jacobian.Mult(reference_flux, physical_flux);
physical_flux /= context.mapping_determinate;
}
void map_physical_flux_to_hdiv_reference(const MappingPointContext &context, const mfem::Vector &physical_flux, mfem::Vector &reference_flux) {
serif_assert<SERiFErrorCode::INVALID_DIMENSION>(physical_flux.Size() == context.inverse_mapping_jacobian.Height(), "Physical flux vector has incompatible dimension with the inverse mapping jacobian ({} <=> {}).", physical_flux.Size(), context.inverse_mapping_jacobian.Height());
reference_flux.SetSize(physical_flux.Size());
context.inverse_mapping_jacobian.Mult(physical_flux, reference_flux);
reference_flux /= context.mapping_determinate;
}
void map_reference_gradient_to_physical(const MappingPointContext &context, const mfem::Vector &reference_gradient, mfem::Vector &physical_gradient) {
serif_assert<SERiFErrorCode::INVALID_DIMENSION>(reference_gradient.Size() == context.inverse_mapping_jacobian.Height(), "Reference gradient vector has incompatible dimension with the inverse mapping jacobian ({} <=> {}).", reference_gradient.Size(), context.inverse_mapping_jacobian.Height());
physical_gradient.SetSize(reference_gradient.Size());
context.inverse_mapping_jacobian.Mult(reference_gradient, physical_gradient);
}
void map_physical_gradient_to_reference(const MappingPointContext &context, const mfem::Vector &physical_gradient, mfem::Vector &reference_gradient) {
serif_assert<SERiFErrorCode::INVALID_DIMENSION>(physical_gradient.Size() == context.mapping_jacobian.Height(), "Physical gradient vector has incompatible dimension with the mapping jacobian ({} <=> {}).", physical_gradient.Size(), context.mapping_jacobian.Height());
reference_gradient.SetSize(physical_gradient.Size());
context.mapping_jacobian.MultTranspose(physical_gradient, reference_gradient);
}
void map_reference_vector_gradient_to_physical(const MappingPointContext &context, const mfem::DenseMatrix &reference_vector_gradient, mfem::DenseMatrix &physical_vector_gradient) {
serif_assert<SERiFErrorCode::INVALID_DIMENSION>(reference_vector_gradient.Width() == context.inverse_mapping_jacobian.Height(), "Reference vector gradient matrix has incompatible dimension with the inverse mapping jacobian ({} <=> {}).", reference_vector_gradient.Width(), context.inverse_mapping_jacobian.Height());
physical_vector_gradient.SetSize(reference_vector_gradient.Height(), context.inverse_mapping_jacobian.Width());
mfem::Mult(reference_vector_gradient, context.inverse_mapping_jacobian, physical_vector_gradient);
}
void map_physical_vector_gradient_to_reference(const MappingPointContext &context, const mfem::DenseMatrix &physical_vector_gradient, mfem::DenseMatrix &reference_vector_gradient) {
serif_assert<SERiFErrorCode::INVALID_DIMENSION>(physical_vector_gradient.Width() == context.mapping_jacobian.Height(), "Physical vector gradient matrix has incompatible dimension with the mapping jacobian ({} <=> {}).", physical_vector_gradient.Width(), context.mapping_jacobian.Height());
reference_vector_gradient.SetSize(physical_vector_gradient.Height(), context.mapping_jacobian.Width());
mfem::Mult(physical_vector_gradient, context.mapping_jacobian, reference_vector_gradient);
}
double map_hdiv_divergence_to_physical(const MappingPointContext &context, double reference_divergence) {
return reference_divergence / context.mapping_determinate;
}
void compute_hdiv_mass_tensor(const MappingPointContext &context, mfem::DenseMatrix &mass_tensor) {
const int dimension = context.mapping_jacobian.Height();
serif_assert<SERiFErrorCode::JACOBIAN_SHAPE_ERROR>((context.mapping_jacobian.Width() == dimension) && (context.mapping_jacobian.Height() == dimension), "Mapping jacobian is not square ({} x {}). It is expected to have shape {}x{}", context.mapping_jacobian.Height(), context.mapping_jacobian.Width(), dimension, dimension);
serif_assert<SERiFErrorCode::INVERTED_JACOBIAN>(context.mapping_determinate > 0.0, "Mapping jacobian determinant is non-positive ({}). The mapping may be inverted or degenerate.", context.mapping_determinate);
mass_tensor.SetSize(dimension, dimension);
mfem::MultAtB(context.mapping_jacobian, context.mapping_jacobian, mass_tensor);
mass_tensor *= 1/context.mapping_determinate;
}
void compute_scalar_diffusion_tensor(const MappingPointContext &context, mfem::DenseMatrix &diffusion_tensor) {
const int dimension = context.mapping_jacobian.Height();
serif_assert<SERiFErrorCode::JACOBIAN_SHAPE_ERROR>((context.inverse_mapping_jacobian.Width() == dimension) && (context.inverse_mapping_jacobian.Height() == dimension), "Inverse mapping jacobian is not square ({} x {}). It is expected to have shape {}x{}", context.inverse_mapping_jacobian.Height(), context.inverse_mapping_jacobian.Width(), dimension, dimension);
serif_assert<SERiFErrorCode::INVERTED_JACOBIAN>(context.mapping_determinate > 0.0, "Mapping jacobian determinant is non-positive ({}). The mapping may be inverted or degenerate.", context.mapping_determinate);
diffusion_tensor.SetSize(dimension, dimension);
mfem::MultAtB(context.inverse_mapping_jacobian, context.inverse_mapping_jacobian, diffusion_tensor);
diffusion_tensor *= context.mapping_determinate;
}
void map_hcurl_field_to_physical(const MappingPointContext &context, const mfem::Vector &reference_field, mfem::Vector &physical_field) {
serif_assert<SERiFErrorCode::INVALID_DIMENSION>(reference_field.Size() == context.inverse_mapping_jacobian.Height(), "Reference field vector has incompatible dimension with the inverse mapping jacobian ({} <=> {}).", reference_field.Size(), context.inverse_mapping_jacobian.Height());
physical_field.SetSize(reference_field.Size());
context.inverse_mapping_jacobian.MultTranspose(reference_field, physical_field);
}
void map_physical_field_to_hcurl_reference(const MappingPointContext &context, const mfem::Vector &physical_field, mfem::Vector &reference_field) {
serif_assert<SERiFErrorCode::INVALID_DIMENSION>(physical_field.Size() == context.mapping_jacobian.Height(), "Physical field vector has incompatible dimension with the mapping jacobian ({} <=> {}).", physical_field.Size(), context.mapping_jacobian.Height());
reference_field.SetSize(physical_field.Size());
context.mapping_jacobian.MultTranspose(physical_field, reference_field);
}
void map_hcurl_curl_to_physical(const MappingPointContext &context, const mfem::Vector &reference_curl, mfem::Vector &physical_curl) {
serif_assert<SERiFErrorCode::INVALID_DIMENSION>(reference_curl.Size() == context.mapping_jacobian.Width(), "Reference curl vector has incompatible dimension with the mapping jacobian ({} <=> {}).", reference_curl.Size(), context.mapping_jacobian.Height());
physical_curl.SetSize(reference_curl.Size());
context.mapping_jacobian.Mult(reference_curl, physical_curl);
physical_curl /= context.mapping_determinate;
}
void map_physical_curl_to_hcurl_reference(const MappingPointContext &context, const mfem::Vector &physical_curl, mfem::Vector &reference_curl) {
serif_assert<SERiFErrorCode::INVALID_DIMENSION>(physical_curl.Size() == context.inverse_mapping_jacobian.Width(), "Physical curl vector has incompatible dimension with the inverse mapping jacobian ({} <=> {}).", physical_curl.Size(), context.inverse_mapping_jacobian.Height());
reference_curl.SetSize(physical_curl.Size());
context.inverse_mapping_jacobian.Mult(physical_curl, reference_curl);
reference_curl *= context.mapping_determinate;
}
void compute_hdiv_mass_tensor_variation(const MappingPointContext &context, const MappingPointVariation &variation, mfem::DenseMatrix &mass_tensor_variation) {
const double determinate = context.mapping_determinate;
const double determinate_variation = variation.mapping_determinate_variation;
const int dimension = context.inverse_mapping_jacobian.Width();
mass_tensor_variation.SetSize(dimension, dimension);
serif_assert<SERiFErrorCode::INVERTED_JACOBIAN>(determinate > 0.0, "Mapping jacobian determinant is non-positive ({}). The mapping may be inverted or degenerate.", determinate);
serif_assert<SERiFErrorCode::NON_FINITE_DETERMINATE>(utils::checks::is_finite(determinate), "Mapping jacobian determinant is non-finite ({}). The mapping may be inverted or degenerate.", determinate);
const mfem::DenseMatrix& jacobian = context.mapping_jacobian;
const mfem::DenseMatrix& jacobian_variation = variation.mapping_jacobian_variation;
const double inverse_determinate = 1.0 / determinate;
const double determinate_scale = determinate_variation * inverse_determinate;
for (int row = 0; row < dimension; ++row) {
for (int col = 0; col < dimension; ++col) {
double gram{0.0};
double gram_variation{0.0};
for (int inner = 0; inner < dimension; ++inner) {
gram += jacobian(inner, row) * jacobian(inner, col);
gram_variation += jacobian(inner, row) * jacobian_variation(inner, col) + jacobian_variation(inner, row) * jacobian(inner, col);
}
mass_tensor_variation(row, col) = inverse_determinate * (gram_variation - determinate_scale * gram);
}
}
}
}

112
src/lib/serif/mapping/volume_cache.cpp
View File

@@ -1,112 +0,0 @@
#include "serif/mapping/volume_cache.hpp"
#include "serif/utils/error/errors.hpp"
#include "mfem/linalg/densemat.hpp"
#include "mfem/linalg/vector.hpp"
namespace serif::mapping {
using namespace utils::errors;
namespace {
// Replace all of these with a fast serialization class at some point. Probably in something like a serif/serialize module
void pack_vector(double*& destination, const mfem::Vector& vector) {
const int dim = vector.Size();
std::copy_n(vector.HostRead(), dim, destination);
destination += dim;
}
void pack_matrix(double*& destination, const mfem::DenseMatrix& matrix) {
serif_assert<SERiFErrorCode::MATRIX_SHAPE_ERROR>(matrix.Height() == matrix.Width(), "Matrix must be square to pack");
const int dim = matrix.Height();
std::copy_n(matrix.HostRead(), dim * dim, destination);
destination += dim * dim;
}
void unpack_vector(const double*& source, mfem::Vector& vector, const int dim) {
vector.SetSize(dim);
std::copy_n(source, dim, vector.HostWrite());
source += dim;
}
void unpack_matrix(const double*& source, mfem::DenseMatrix& matrix, const int dim) {
matrix.SetSize(dim, dim);
std::copy_n(source, dim * dim, matrix.HostWrite());
source += dim * dim;
}
}
void VolumeMappingCache::set_size(int point_count, int dimension) {
serif_assert<SERiFErrorCode::INVALID_DIMENSION>(point_count < 0 || dimension < 1 || dimension > 3, "Invalid point count or dimension for volume mapping cache. Point count: {}, dimension: {}", point_count, dimension);
const int stride = 3 * dimension + 4 * dimension * dimension + 4;
m_data.resize(static_cast<std::size_t>(point_count) * stride);
m_point_count = point_count;
m_dimension = dimension;
m_point_stride = stride;
}
void VolumeMappingCache::validate_point(const int point) const {
serif_assert<SERiFErrorCode::OUT_OF_RANGE_ERROR>(point < 0 || point >= m_point_count, "Point index {} is out of range for volume mapping cache with {} points", point, m_point_count);
}
const double *VolumeMappingCache::get_point_data(const int point) const {
validate_point(point);
return m_data.data() + static_cast<std::size_t>(point) * m_point_stride;
}
void VolumeMappingCache::store(int point, VolumeMappingContext &context) {
validate_point(point);
double *data = m_data.data() + static_cast<std::size_t>(point) * m_point_stride;
// What we build here is a FIFO data structure. we pack
// all the raw data from the mapping context into a single contiguous array of doubles.
// This allows us to avoid repeated allocations as the vector will grow to some appropriate size and then stay roughly there.
pack_vector(data, context.mapping.reference_position);
pack_vector(data, context.mapping.displaced_position);
pack_vector(data, context.mapping.physical_position);
pack_matrix(data, context.mapping.displacement_jacobian);
pack_matrix(data, context.mapping.mapping_jacobian);
pack_matrix(data, context.mapping.inverse_mapping_jacobian);
pack_matrix(data, context.quadrature.J_element_to_physical_inverse);
*data++ = context.mapping.mapping_determinate;
*data++ = context.mapping.compactified ? 1.0 : 0.0; // its possible this could cause issues, have not run into but when unpacked I worry that exact equality to some double value could be a problem. If there are wierd issues this may be a place to look,
*data++ = context.quadrature.mapping_determinate;
*data = context.quadrature.physical_integration_weight;
}
void VolumeMappingCache::load(int point, VolumeMappingContext &context) const {
const double *data = get_point_data(point);
// FIFO data structure so need to unpack in the same order as we packed.
// TODO: Redo this with some sort of fast serialization class so we do not need to worry about ordering
unpack_vector(data, context.mapping.reference_position, m_dimension);
unpack_vector(data, context.mapping.displaced_position, m_dimension);
unpack_vector(data, context.mapping.physical_position, m_dimension);
unpack_matrix(data, context.mapping.displacement_jacobian, m_dimension);
unpack_matrix(data, context.mapping.mapping_jacobian, m_dimension);
unpack_matrix(data, context.mapping.inverse_mapping_jacobian, m_dimension);
unpack_matrix(data, context.quadrature.J_element_to_physical_inverse, m_dimension);
context.mapping.mapping_determinate = *data++;
context.mapping.compactified = *data++ != 0.0; // See comment in store()
context.quadrature.mapping_determinate = *data++;
context.quadrature.physical_integration_weight = *data;
}
void VolumeMappingCache::load_inverse_jacobian(int point, mfem::DenseMatrix &inverse) const {
// Again this is unsafe, and really should be replaced. This was a quick cache scheme that solved the direct issue I had; however, it should be replaced
// with something that has at least a modicum of safety built in (we are pretty good about this in other parts of SERiF but its about 8pm and
// Im quite tired....
const double *data = get_point_data(point) + 3 * m_dimension + 3 * m_dimension * m_dimension; // skip to the inverse jacobian, this is just pointer math...Think back to CS-101. Do that have that in the future, are we still learning about memory like this is a PDP-11?
unpack_matrix(data, inverse, m_dimension);
}
int VolumeMappingCache::get_point_count() const noexcept {
return m_point_count;
}
int VolumeMappingCache::get_dimension() const noexcept {
return m_dimension;
}
}