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

feat(mapping): The mapping code was all moved to the backend, this module now just brings the MFEM backend in

Browse Source 
The mapping code is fundamentally backend dependent. Previously we had DomainMapper here; however, since it depende on MFEM we have moved it to the backend. The mapping module now just includes parts of the backend so that prior code will not break
This commit is contained in:
Emily Boudreaux
2026-11-09 13:16:36 -05:00
parent 30af52ccea
commit c6d1686c5d
14 changed files with 4 additions and 1293 deletions
Download Patch File Download Diff File Expand all files Collapse all files

95
src/include/serif/mapping/coefficients.hpp
View File

@@ -1,95 +0,0 @@
#pragma once
/**
* @file coefficients.hpp
* @author Emily M. Boudreaux
*
* @brief File containing the PhysicalPositionFunctionCoefficient class, which is a custom MFEM coefficient that evaluates a user-defined function at physical positions determined by a domain mapping.
*
* @todo This should probably be moved into a backend/mfem module
*/
#include <functional>
#include "mfem/fem/gridfunc.hpp"
#include "mfem/fem/coefficient.hpp"
#include "mfem/fem/eltrans.hpp"
#include "mfem/fem/intrules.hpp"
#include "mfem/linalg/vector.hpp"
#include "serif/mapping/domain_mapper.hpp"
#include "serif/mapping/evaluator.hpp"
namespace serif::mapping {
/**
* @class PhysicalPositionFunctionCoefficient
*
* @brief A custom MFEM coefficient that evaluates a user-defined function at physical positions determined by a domain mapping.
*
* @par Example Usage
* Lets say we have some analytic density profile of the form
*
* \f[
* \rho(r) = \rho_0 \exp(-r^2 / R^2)
* \f]
*
* We might want to evaluate this function over some computational domain but in physical space since
* the analytic form is setup to be evaluated in physical space. We cannot simply calculate the
* position of each dof on the mesh and evaluate the function there since the mesh may be deformed. Rather
* if we use this coefficient we can offload the mapping to the mapping module
*
* \code{.cpp}
* // assuming there is already some defined mapper
*
* auto density = [](const mfem::Vector& x) -> double {
* double r = x.Norml2();
* return rho0 * std::exp(-r*r / (R*R));
* };
*
* PhysicalPositionFunctionCoefficient physical_density(mapper, displacement, compactification_coordinate, density);
* // Then there would need to be some loopup code to get the rigth element transformation and integration point
* // for the physical position you want to evaluate the function at.
* double denity = physical_density.Eval(transformation, ip);
*
* \endcode
*/
class PhysicalPositionFunctionCoefficient : public mfem::Coefficient {
public:
using Func = std::function<double(const mfem::Vector& x)>; ///< The type of user-defined function must be a scalar valued function of a 3D vector.
/**
* @brief Construct a PhysicalPositionFunctionCoefficient with a given domain mapper, displacement, compactification coordinate, and user-defined function.
*
* @param mapper The domain mapper that defines the mapping from reference to physical coordinates.
* @param displacement The grid function representing the displacement field.
* @param compactification_coordinate The grid function representing the compactification coordinate.
* @param func The user-defined function to evaluate at physical positions.
*
* @todo The displacement grid function really should not be needed in this constructor since all required information
* is already in the mapper. Same goes for the compactification coordinate. I kept these to maintain compatibility
* with some older code; however, we should consider removing them.
*/
PhysicalPositionFunctionCoefficient(
const DomainMapper& mapper,
const mfem::GridFunction& displacement,
const mfem::GridFunction& compactification_coordinate,
Func func
);
/**
* @brief Evaluate the user-defined function at a given point in the reference domain, mapped to physical space.
*
* @param transformation The element transformation that defines the mapping from reference to physical coordinates.
* @param ip The integration point in the reference domain at which to evaluate the function.
*
* @return The value of the user-defined function evaluated at the physical position corresponding to the given reference point.
*/
double Eval(
mfem::ElementTransformation& transformation,
const mfem::IntegrationPoint& ip
) override;
private:
Func m_f;
GridFunctionMappingEvaluator m_mapping; ///< The mapping evaluator that computes physical positions from reference coordinates.
};
}

82
src/include/serif/mapping/compactification/compactification.hpp
View File

@@ -1,82 +0,0 @@
#pragma once
/**
* @file compactification.hpp
* @author Emily M. Boudreaux
*
* @brief File containing the definition of the ExteriorDomainMap interface, which is used to define the mapping from reference coordinates to physical coordinates in the exterior domain.
*/
#include <string_view>
#include <memory>
#include "serif/mapping/status.hpp"
#include "serif/mapping/compactification/input_result.hpp"
#include "serif/mapping/compactification/variation.hpp"
namespace serif::mapping::compactification {
/**
* @class ExteriorDomainMap
*
* @brief Abstract interface for some mapping from reference coordinates to physical coordinates in the exterior domain.
* Currently the only concrete implementation of this is the kelvin compactification.
*/
class ExteriorDomainMap {
public:
virtual ~ExteriorDomainMap() = default;
/**
* @brief Check if the mapping is affine in the displacement, that is if the mapping is linear in the displacement and does not depend on the displacement.
*
* @return True if the mapping is affine in the displacement, false otherwise.
*/
[[nodiscard]] virtual bool is_affine_in_displacement() const noexcept = 0;
/**
* @brief Evaluate the mapping at a given point in the reference domain.
*
* @param input The input to the mapping, which contains the reference coordinates and the displacement.
* @param result The result of the mapping, which contains the physical coordinates and the Jacobian of the mapping.
*
* @return The status of the mapping evaluation, which indicates if the evaluation was successful or if there was an error.
*
*/
[[nodiscard]] virtual MappingStatus evaluate(
const ExteriorMapInput &input,
ExteriorMapResult &result
) const = 0;
/**
* @brief Evaluate the variation of the mapping at a given point in the reference domain.
*
* @param input The input to the mapping, which contains the reference coordinates and the displacement.
* @param result The result of the mapping, which contains the physical coordinates and the Jacobian of the mapping.
* @param direction The direction of the variation, which is a vector in the reference domain.
* @param variation The variation of the mapping, which contains the variation of the physical coordinates and the variation of the Jacobian.
*
* @return The status of the mapping evaluation, which indicates if the evaluation was successful or if there was an error.
*/
[[nodiscard]] virtual MappingStatus evaluate_variation(
const ExteriorMapInput &input,
const ExteriorMapResult &result,
const ExteriorMapDirection &direction,
ExteriorMapVariation &variation
) const = 0;
/**
* @brief Return the name of the mapping, which is used for debugging and logging.
*
* @return The name of the mapping.
*/
[[nodiscard]] virtual std::string_view get_name() const noexcept = 0;
/**
* @brief Clone the mapping, which creates a new instance of the mapping with the same parameters.
*
* @return A unique pointer to the cloned mapping.
*/
[[nodiscard]] virtual std::unique_ptr<const ExteriorDomainMap> clone() const = 0;
};
}

36
src/include/serif/mapping/compactification/input_result.hpp
View File

@@ -1,36 +0,0 @@
#pragma once
/**
* @file input_result.hpp
* @author Emily M. Boudreaux
*
* @brief File containing the input and result structures for the compactification mapping.
*/
#include "mfem/linalg/vector.hpp"
#include "mfem/linalg/densemat.hpp"
namespace serif::mapping::compactification {
/**
* @struct ExteriorMapInput
*
* @brief Structure representing the input to the exterior domain mapping. This includes the reference position, displaced position, displacement Jacobian, compactification coordinate gradient, and compactification coordinate.
*/
struct ExteriorMapInput {
const mfem::Vector &reference_position; ///< Position in reference space to evaluate the mapping at.
const mfem::Vector &displaced_position; ///< Position in displaced space to evaluate the mapping at (note this is displaced but not yet compactified).
const mfem::DenseMatrix &displacement_jacobian; ///< Jacobian of the displacement mapping at the reference position.
const mfem::Vector &compactification_coordinate_gradient; ///< Gradient of the compactification coordinate at the reference position.
double compactification_coordinate; ///< The compactification coordinate at the reference position. This ranges from 0 to 1, where 0 is the inner surface and 1 is the outer surface of the reference domain.
};
/**
* @struct ExteriorMapResult
*
* @brief Structure representing the result of the exterior domain mapping. This includes the physical position and the Jacobian of the mapping.
*/
struct ExteriorMapResult {
mfem::Vector physical_position; ///< Position in physical space after applying the mapping.
mfem::DenseMatrix mapping_jacobian; ///< Jacobian of the mapping at the reference position, which maps from reference coordinates to physical coordinates.
};
}

197
src/include/serif/mapping/compactification/kelvin.hpp
View File

@@ -1,197 +0,0 @@
#pragma once
/**
* @file kelvin.hpp
* @author Emily M. Boudreaux
*
* @brief File containing the definition of the KelvinCompactification class, which is a concrete implementation of the ExteriorDomainMap interface for the Kelvin compactification.
*/
#include <string_view>
#include "mfem/linalg/vector.hpp"
#include "serif/mapping/status.hpp"
#include "serif/mapping/compactification/compactification.hpp"
#include "serif/utils/user_options/mapping/compactification/kelvin.hpp"
namespace serif::mapping::compactification {
/**
* @class KelvinCompactification
*
* @brief Concrete implementation of the ExteriorDomainMap interface for the Kelvin compactification.
*
* Broadly, for some reference position \f$\mathbf{x_{r}}\f$ on some reference domain \f$\Omega_{r}\f$ the Kelvin compactification maps to a physical position \f$\mathbf{x_{p}}\f$ on some physical domain \f$\Omega_{p}\f$ through the mapping
* \f[
* \mathbf{x_{p}} = \frac{\mathbf{x_{r}}}{1 - \frac{|\mathbf{x_{r}}|}{R_{\infty}}}
* \f]
*
* This results in the inner surface of the domain being mapped to the exterior surface of the physical domain while the outer
* surface of the exterior domain is mapped to infinity. There are some implications of this, namely that we
* loose accuracy in all things as we approach the outer surface of the domain. However, it is nice in so far
* as we can apply true boundary conditions at infinity.
*/
class KelvinCompactification final : public ExteriorDomainMap {
public:
/**
* @brief Construct a KelvinCompactification with the given options.
*
* @param options The options for the Kelvin compactification, which must be valid.
*
* @throws std::invalid_argument If the reference stellar radius is not finite or positive
* @throws std::invalid_argument If the reference infinity radius is not finite or positive
* @throws std::invalid_argument If the reference infinity radius is not greater than the reference stellar radius
* @throws std::invalid_argument If the coordinate tolerance is not finite or positive or is greater than or equal to one
*/
explicit KelvinCompactification(
const utils::user_options::KelvinCompactificationOptions& options
);
/**
* @brief Check if the mapping is affine in the displacement. For the Kelvin compactification, this is always true
*
* @return false
*/
[[nodiscard]] bool is_affine_in_displacement() const noexcept override;
/**
* @brief Evaluate the mapping at a given point in the reference domain.
*
* @param[in] input The input point in the reference domain to evaluate the mapping at.
* @param[out] result The result of the mapping evaluation, which will be filled in by the method.
*
* @return A status indicating the success or failure of the mapping evaluation. Mapping may fail if
* - the size of the input reference position is less than or equal to zero
* - the input displaced position size does not match the input reference position size
* - the input compactification coordinate gradient size does not match the input reference position size
* - the displacement_jacobian height does not match the input reference position size
* - the displacement_jacobian width does not match the input reference position size
* - the input reference, displaced, compactification coordinate gradient or jacobian contains non-finite values
* - There is some issue computing the radial factors (this error will be passed through)
*
* The mapping jacobian comes from trivial application of the product rule
* recall: \f$r_i = scale * x_i\f$ where \f$r\f$ is the physical position and \f$x\f$ is the displaced position.
* then we can differentiate wrt. \f$X_j\f$ holding nothing fixed. Note the capital \f$X\f$ here, this is the mesh coordinate not the displaced position.
* Lets call this jacobian \f$F\f$
* \f[
* F = \frac{\partial r_i}{\partial X_{j}}
* \f]
*
* \f$F\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[
* F = scale * \frac{\partial x_i}{\partial X_j} + x_i * \frac{\partial scale}{\partial X_j}
* \f]
*
* In the implementation you can see the displacement jacobian (\f$\frac{\partial x_i}{\partial X_j}\f$) and the scale derivative (\f$\frac{\partial scale}{\partial X_j}\f$) being applied to compute the mapping jacobian.
*/
[[nodiscard]] MappingStatus evaluate(
const ExteriorMapInput& input,
ExteriorMapResult& result
) const override;
/**
* @brief Evaluate the variation of the mapping at a given point in the reference domain.
*
* @param[in] input The input point in the reference domain to evaluate the mapping at.
* @param[in] result The result of the mapping evaluation, which will be used in the variation computation.
* @param[in] direction The direction in which to evaluate the variation.
* @param[out] variation The result of the variation evaluation, which will be filled in by the method.
*
* @return A status indicating the success or failure of the variation evaluation. The reasons for a failure are effectively the same
* as the reasons for a failure in the evaluate method just extended to cover the additional inputs.
*/
[[nodiscard]] MappingStatus evaluate_variation(
const ExteriorMapInput& input,
const ExteriorMapResult& result,
const ExteriorMapDirection& direction,
ExteriorMapVariation& variation
) const override;
/**
* @brief Get the name of the mapping, which is "KelvinCompactification" for this class.
* @return The name of the mapping.
*/
[[nodiscard]] std::string_view get_name() const noexcept override;
/**
* @brief Clone the mapping, returning a unique pointer to a new instance of the mapping.
* @return A unique pointer to a new instance of the mapping.
*/
[[nodiscard]] std::unique_ptr<const ExteriorDomainMap> clone() const override;
/**
* @brief Get the reference stellar radius for the mapping, which is the radius of the inner surface of the reference domain.
* @return The reference stellar radius for the mapping.
*/
[[nodiscard]] double get_reference_stellar_radius() const noexcept;
/**
* @brief Get the reference infinity radius for the mapping, which is the radius of the outer surface of the reference domain.
* @return The reference infinity radius for the mapping.
*/
[[nodiscard]] double get_reference_infinity_radius() const noexcept;
/**
* @brief Get the coordinate tolerance for the mapping, which is the tolerance for the compactification coordinate.
* @return The coordinate tolerance for the mapping.
*/
[[nodiscard]] double get_coordinate_tolerance() const noexcept;
/**
* @brief Evaluate the mapping jacobian variation factors at some given compactification coordinate.
*
* @param[in] compactification_coordinate The compactification coordinate at which to evaluate the mapping jacobian variation factors. (Note this runs from 0 to 1, where 0 is the inner surface and 1 is the outer surface of the reference domain)
* @param[in] coordinateGradient The gradient of the compactification coordinate with respect to the mesh coordinates, which must be finite and have size 3.
*
* @param[out] scale The scale factor for the mapping jacobian, which will be filled in by the method.
* @param[out] scaleGradient The gradient of the scale factor with respect to the mesh coordinates, which will be filled in by the method.
*
* @return The status of the mapping evaluation.
*/
[[nodiscard]] MappingStatus evaluate_jacobian_variation_factors(
double compactification_coordinate,
const mfem::Vector &coordinateGradient,
double &scale,
mfem::Vector &scaleGradient
) const;
private:
/**
* @struct RadialFactors
*
* @brief A struct to hold the radial factors for the Kelvin compactification mapping.
*/
struct RadialFactors {
double compactification_coordinate; ///< The compactification coordinate, which runs from 0 to 1, where 0 is the inner surface and 1 is the outer surface of the reference domain.
double computational_radius; ///< The computational radius, which is the radius of the reference domain at the given compactification coordinate.
double scale; ///< The scale factor for the mapping jacobian.
double scale_derivative; ///< The derivative of the scale factor with respect to the compactification coordinate.
};
private:
/**
* @brief Compute the radial factors for the Kelvin compactification mapping at a given compactification coordinate.
*
* @param[in] compactification_coordinate The compactification coordinate at which to compute the radial factors. (Note this runs from 0 to 1, where 0 is the inner surface and 1 is the outer surface of the reference domain)
* @param[out] radial_factors The radial factors, which will be filled in by the method.
*
* @return The status of the mapping evaluation.
*
* @note See the implementation file for more details on the actual algorithm used here.
*/
[[nodiscard]] MappingStatus compute_radial_factors(
double compactification_coordinate,
RadialFactors &radial_factors
) const;
private:
utils::user_options::KelvinCompactificationOptions m_options; ///< The options for the Kelvin compactification, which must be valid.
const std::string_view m_name = "KelvinCompactification"; ///< The name of the mapping, which is "KelvinCompactification" for this class.
};
}

35
src/include/serif/mapping/compactification/variation.hpp
View File

@@ -1,35 +0,0 @@
#pragma once
/**
* @file variation.hpp
* @author Emily M. Boudreaux
*
* @brief File containing the structures for the variation of the compactification mapping.
*/
#include "mfem/linalg/vector.hpp"
#include "mfem/linalg/densemat.hpp"
namespace serif::mapping::compactification {
/**
* @struct ExteriorMapDirection
*
* @brief Structure representing the direction of the variation in the exterior domain mapping.
* This includes the variation of the displaced position and the variation of the displacement Jacobian.
*/
struct ExteriorMapDirection {
const mfem::Vector &displaced_position_variation; ///< Variation of the position in displaced space.
const mfem::DenseMatrix &displacement_jacobian_variation; ///< Variation of the Jacobian of the displacement mapping.
};
/**
* @struct ExteriorMapVariation
*
* @brief Structure representing the variation of the exterior domain mapping.
* This includes the variation of the physical position and the variation of the mapping Jacobian.
*/
struct ExteriorMapVariation {
mfem::Vector physical_position_variation; ///< Variation of the position in physical space.
mfem::DenseMatrix mapping_jacobian_variation; ///< Variation of the Jacobian of the mapping.
};
}

92
src/include/serif/mapping/context.hpp
View File

@@ -1,92 +0,0 @@
#pragma once
/**
* @file context.hpp
* @author Emily M. Boudreaux
*
* @brief File containing the mapping context structures for volume and face mappings.
*
* @todo Like a lot of the parts of the mapping module written before I started factoring mfem stuff into
* backend modules, a good amount of this might be movable to some mfem backend
*
* All of these different "context" structures are built to contain sufficient information that the
* caller which gets one of these back can evaluate the mapping and its derivatives properly
* at some point.
*/
#include "mfem/linalg/densemat.hpp"
#include "mfem/linalg/vector.hpp"
namespace serif::mapping {
/**
* @struct VolumeQuadratureContext
*
* @brief A structure representing the context of a volume quadrature point, including the inverse Jacobian of the mapping and the determinant of the mapping.
*
* @todo Add a worked example here
*/
struct VolumeQuadratureContext {
mfem::DenseMatrix J_element_to_physical_inverse; ///< The inverse Jacobian of the mapping from the element to the physical space.
double mapping_determinate; ///< The determinant of the mapping from the element to the physical space,
double physical_integration_weight; ///< The physical integration weight for the quadrature point.
};
/**
* @struct FaceQuadratureContext
*
* @brief A structure representing the context of a face quadrature point, including the normal vector, surface weight, and scale factor for the dot product of the velocity and normal vector.
*
* @todo Add a worked example here
*/
struct FaceQuadratureContext {
mfem::Vector normal; ///< The normal vector at the face quadrature point.
double ds; ///< The surface weight for the face quadrature point.
double v_dot_n_scale; ///< The scale factor for the dot product of the velocity and normal vector at the face quadrature point.
};
/**
* @struct MappingPointContext
*
* @brief A structure representing the context of a mapping point, including the reference position, displaced position, physical position, displacement Jacobian, mapping Jacobian, inverse mapping Jacobian, mapping determinant, and compactification status.
*
* @todo Add a worked example here
*/
struct MappingPointContext {
mfem::Vector reference_position; ///< The reference position of the mapping point.
mfem::Vector displaced_position; ///< The displaced position of the mapping point.
mfem::Vector physical_position; ///< The physical position of the mapping point.
mfem::DenseMatrix displacement_jacobian; ///< The displacement Jacobian of the mapping point.
mfem::DenseMatrix mapping_jacobian; ///< The mapping Jacobian of the mapping point.
mfem::DenseMatrix inverse_mapping_jacobian; ///< The inverse mapping Jacobian of the mapping point.
double mapping_determinate{0.0}; ///< The determinant of the mapping Jacobian at the mapping point.
bool compactified{false}; ///< A flag indicating whether the mapping point is compactified (true) or not (false).
};
/**
* @struct VolumeMappingContext
*
* @brief A structure representing the context of a volume mapping, including the mapping point context and volume quadrature context.
*
* @todo Add a worked example here
*/
struct VolumeMappingContext {
MappingPointContext mapping; ///< The mapping point context for the volume mapping.
VolumeQuadratureContext quadrature; ///< The volume quadrature context for the volume mapping.
};
/**
* @struct FaceMappingContext
*
* @brief A structure representing the context of a face mapping, including the mapping point context and face quadrature context.
*
* @todo Add a worked example here
*/
struct FaceMappingContext {
MappingPointContext mapping; ///< The mapping point context for the face mapping.
FaceQuadratureContext quadrature; ///< The face quadrature context for the face mapping.
mfem::Vector reference_normal; ///< The reference normal vector for the face mapping.
double reference_surface_weight{0.0}; ///< The reference surface weight for the face mapping.
double physical_surface_weight{0.0}; ///< The physical surface weight for the face mapping.
};
}

295
src/include/serif/mapping/data.hpp
View File

@@ -1,295 +0,0 @@
#pragma once
/**
* @file data.hpp
* @author Emily M. Boudreaux
*
* @brief File containing some misc data structures for the mapping module. This was the very first file
* I ported over from MeanField and I think because of that its organization is actually quite wonky.
* For one its int this poorly named "data" file, which could mean any number of things. Regardless,
* these are basically low leval containers for useful stuff and also some scratch workspaces
* for mapping operations.
*
* @todo As said above there is some refactoring to make this cleaner which could happen.
*/
#include "mfem/fem/fe/fe_base.hpp"
#include "mfem/fem/fespace.hpp"
#include "mfem/linalg/vector.hpp"
#include "mfem/linalg/densemat.hpp"
#include "mfem/mesh/mesh.hpp"
#include <cstdint>
#include <map>
#include <tuple>
#include "serif/mapping/compactification/kelvin.hpp"
#include "serif/mapping/context.hpp"
#include "serif/mapping/data.hpp"
namespace serif::mapping {
/**
* @enum FaceElementSide
*
* @brief An enumeration representing the side of a face element, either ELEMENT_1 or ELEMENT_2. There is no inherent ordering
* to these; however, broadly I have been thinking about ELEMENT_1 as the "inner" side and ELEMENT_2 as the "outer" side.
*/
enum class FaceElementSide : uint8_t {
ELEMENT_1, ///< Represents the first side of a face element.
ELEMENT_2 ///< Represents the second side of a face element.
};
/**
* @class ElementDisplacementData
*
* @brief A class representing the displacement data for a finite element, including the finite element, degree of freedom matrix, dimension, and ordering.
*/
class ElementDisplacementData {
public:
/**
* @brief Constructs an ElementDisplacementData object with the given finite element, displacement degrees of freedom, and ordering.
*
* @param element The finite element associated with the displacement data.
* @param displacement_dofs The vector containing the displacement degrees of freedom.
* @param ordering The ordering type for the degrees of freedom (default is byNODES).
*
* @throws std::invalid_argument if the finite element has zero or negative degrees of freedom,
* or if the size of the displacement_dofs vector is not a positive multiple of the element's degree of freedom count.
*
* If ordering is set to byNODES then the dof matrix will be populated such that for row i and column j, with n total dofs, the value is set to
* dof_matrix(i, j) = displacement_dofs(i + j * n). Otherwise dof_matrix(i, j) = displacement_dofs(j + i * dim) for some number of dimensions dim (read as the size of the displacement_dof vector / dof count
* and implicitly truncated to an integer).
*/
ElementDisplacementData(
const mfem::FiniteElement& element,
mfem::Vector& displacement_dofs,
mfem::Ordering::Type ordering = mfem::Ordering::byNODES
);
/**
* @brief Returns the finite element associated with the displacement data.
* @return A constant reference to the finite element.
*/
[[nodiscard]] const mfem::FiniteElement& get_element() const noexcept;
/**
* @brief Returns the degree of freedom matrix for the displacement data.
* @return A constant reference to the degree of freedom matrix.
*/
[[nodiscard]] const mfem::DenseMatrix &get_dof_matrix() const noexcept;
/**
* @brief Returns the dimension of the displacement data.
* @return The dimension of the displacement data.
*/
[[nodiscard]] int get_dimension() const noexcept;
/**
* @brief Returns the degree of freedom count for the displacement data.
* @return The degree of freedom count for the displacement data.
*/
[[nodiscard]] int get_dof_count() const noexcept;
/**
* @brief Returns the ordering type for the degrees of freedom in the displacement data.
* @return The ordering type for the degrees of freedom.
*/
[[nodiscard]] mfem::Ordering::Type get_ordering() const noexcept;
private:
const mfem::FiniteElement *m_element; ///< Pointer to the finite element associated with the displacement data.
mfem::DenseMatrix m_dof_matrix; ///< Matrix representing the degrees of freedom for the displacement data.
int m_dimension; ///< The dimension of the displacement data (number of components).
mfem::Ordering::Type m_ordering; ///< The ordering type for the degrees of freedom (byNODES or byVDIM).
};
/**
* @struct CompactificationPointData
*
* @brief A structure representing the compactification point data, including the compactification coordinate and coordinate gradient.
*/
struct CompactificationPointData {
double compactification_coordinate{0.0}; ///< The compactification coordinate value at the point.
mfem::Vector coordinate_gradient; ///< The gradient of the compactification coordinate at the point.
};
/**
* @brief Creates an ElementDisplacementData object from a finite element and a vector of displacement degrees of freedom.
*
* @param element The finite element associated with the displacement data.
* @param displacement_dofs The vector containing the displacement degrees of freedom.
* @return An ElementDisplacementData object constructed from the provided finite element and displacement degrees of freedom.
*
* @note This is just a convenience function which calls the constructor for ElementDisplacementData with the ordering set to byNODES. If you want to use a different ordering, call the constructor directly.
*/
[[nodiscard]] ElementDisplacementData element_displacement_data_from_element_vdofs(
const mfem::FiniteElement& element,
mfem::Vector &displacement_dofs
);
/**
* @class ElementCompactificationData
*
* @brief A class representing the compactification data for a finite element, including the finite element and degrees of freedom.
*/
class ElementCompactificationData {
public:
/**
* @brief Constructs an ElementCompactificationData object with the given finite element and degrees of freedom.
*
* @param element The finite element associated with the compactification data.
* @param dofs The vector containing the degrees of freedom for the compactification data.
*
* @throws std::invalid_argument if the finite element is not a scalar finite element, does not provide a gradient, or has zero or negative degrees of freedom,
* or if the size of the dofs vector does not match the finite element's degree of freedom count.
*/
ElementCompactificationData(
const mfem::FiniteElement& element,
const mfem::Vector& dofs
);
/**
* @brief Returns the finite element associated with the compactification data.
* @return A constant reference to the finite element.
*/
[[nodiscard]] const mfem::FiniteElement& get_element() const noexcept;
/**
* @brief Returns the degrees of freedom for the compactification data.
* @return A constant reference to the vector of degrees of freedom.
*/
[[nodiscard]] const mfem::Vector &get_dofs() const noexcept;
/**
* @brief Returns the degree of freedom count for the compactification data.
* @return The degree of freedom count for the compactification data.
*/
[[nodiscard]] int get_dof_count() const noexcept;
private:
const mfem::FiniteElement *m_element; ///< Pointer to the finite element associated with the compactification data.
mfem::Vector m_dofs; ///< Vector representing the degrees of freedom for the compactification data.
};
struct ElementMappingData {
const ElementDisplacementData &displacement;
const ElementCompactificationData &compactification;
};
/**
* @class DMWorkspace
*
* @brief A class representing a workspace for displacement mapping operations, including various vectors and matrices used in the mapping process.
* primarily this is just a separate class so we can more easily make this code thread-safe in the future. At construction of
* an instance of this class loads of Heap allocation is done (or really whatever mfem does internally for allocation for vector which dependents on
* the selected device) and then that allocation can be reused.
*/
class DMWorkspace {
public:
/**
* @brief Constructs a DMWorkspace object with the specified dimension.
*
* @param dimension The dimension of the workspace (default is 3).
*
* @throws std::invalid_argument if the specified dimension is not positive and finite.
*/
explicit DMWorkspace (int dimension = 3);
/**
* @brief Sets the dimension of the workspace.
*
* @param dimension The new dimension of the workspace.
*
* @throws std::invalid_argument if the specified dimension is not positive and finite.
*
* @post All internal vectors are sized to dimension
* @post All internal matrices are sized to dimension x dimension
* @post The compactification coordinate is set to 0.0
*/
void set_dimension(int dimension);
/**
* @brief Returns the dimension of the workspace.
* @return The dimension of the workspace.
*/
[[nodiscard]] int get_dimension() const noexcept;
public: // This is a workspace class, these members are kept public so we do not need indirection through getters to access them
mfem::Vector m_shape; ///< Pre-allocated vector for shape function values at a quadrature point.
mfem::DenseMatrix m_reference_dshape; ///< Pre-allocated matrix for reference shape function derivatives at a quadrature point.
mfem::DenseMatrix m_mesh_dshape; ///< Pre-allocated matrix for mesh shape function derivatives at a quadrature point.
mfem::DenseMatrix m_reference_field_jacobian; ///< Pre-allocated matrix for the reference field Jacobian at a quadrature point.
mfem::Vector m_field_value; ///< Pre-allocated vector for field values at a quadrature point.
mfem::DenseMatrix m_field_jacobian; ///< Pre-allocated matrix for field Jacobian at a quadrature point.
mfem::Vector m_compactification_shape; ///< Pre-allocated vector for compactification shape function values at a quadrature point.
mfem::DenseMatrix m_compactification_dshape; ///< Pre-allocated matrix for compactification shape function derivatives at a quadrature point.
CompactificationPointData m_compactification_point; ///< Pre-allocated data for compactification point.
mfem::Vector m_reference_normal; ///< Pre-allocated vector for reference normal at a quadrature point.
mfem::Vector m_mapped_normal; ///< Pre-allocated vector for mapped normal at a quadrature point.
mfem::DenseMatrix m_full_element_jacobian; ///< Pre-allocated matrix for full element Jacobian at a quadrature point.
mfem::Vector m_vector_temp; ///< Pre-allocated temporary vector for intermediate calculations.
mfem::DenseMatrix m_matrix_temp_1; ///< Pre-allocated temporary matrix for intermediate calculations.
mfem::DenseMatrix m_matrix_temp_2; ///< Pre-allocated temporary matrix for intermediate calculations.
compactification::ExteriorMapResult m_exterior_result; ///< Pre-allocated result for exterior mapping at a quadrature point.
compactification::ExteriorMapVariation m_exterior_variation; ///< Pre-allocated variation for exterior mapping at a quadrature point.
private:
int m_dimension; ///< The number of spatial dimensions for the workspace (e.g., 2D or 3D).
};
/**
* @struct DMFixedMeshGeometry
*
* @brief A structure representing the fixed mesh geometry at a specific point, including position, Jacobian, inverse Jacobian, and weight.
*/
struct DMFixedMeshGeometry {
mfem::Vector position; ///< The physical position of the point in the mesh.
mfem::DenseMatrix jacobian; ///< The Jacobian matrix of the mapping from the reference element to the physical element at the point.
mfem::DenseMatrix inverse_jacobian; ///< The inverse of the Jacobian matrix of the mapping from the reference element to the physical element at the point.
double weight{0.0}; ///< The mapping weight associated with the point, typically used for integration purposes.
};
/**
* @struct DMGeometryCacheStatistics
*
* @brief A structure representing statistics for the geometry cache, including hit and miss counts and the number of points cached.
*/
struct DMGeometryCacheStatistics {
std::uint64_t hits{0}; ///< The number of cache hits, indicating how many times a requested geometry was found in the cache.
std::uint64_t misses{0}; ///< The number of cache misses, indicating how many times a requested geometry was not found in the cache and had to be computed.
std::size_t points{0}; ///< The number of points currently stored in the cache.
};
/**
* @struct DMReferenceDisplacementSample
*
* @brief A structure representing a sample of the reference displacement, including the value and gradient of the displacement at a specific point.
*/
struct DMReferenceDisplacementSample {
const mfem::Vector& value; ///< The value of the reference displacement at the point.
const mfem::DenseMatrix& gradient; ///< The gradient of the reference displacement at the point.
};
/**
* @struct DMGeometryCache
*
* @brief A structure representing a cache for geometry data, including the mesh, sequence numbers, and cached points.
*/
struct DMGeometryCache {
bool enabled{true}; ///< A flag indicating whether the geometry cache is enabled (true) or disabled (false).
const mfem::Mesh* mesh{nullptr}; ///< A pointer to the mesh associated with the cached geometry data.
long sequence{-1}; ///< The sequence number of the cached geometry data.
long nodesSequence{-1}; ///< The sequence number of the cached nodes data.
std::map<std::tuple<int, double, double, double>, DMFixedMeshGeometry> points; ///< A map storing the cached geometry points.
uint64_t hits{0}; ///< The number of cache hits.
uint64_t misses{0}; ///< The number of cache misses.
};
}

166
src/include/serif/mapping/domain_mapper.hpp
View File

@@ -1,166 +0,0 @@
#pragma once
#include <cstdint>
#include <memory>
#include "mfem/linalg/vector.hpp"
#include "mfem/linalg/densemat.hpp"
#include "mfem/fem/eltrans.hpp"
#include "mfem/fem/intrules.hpp"
#include "serif/mapping/compactification/kelvin.hpp"
#include "serif/mapping/context.hpp"
#include "serif/mapping/data.hpp"
#include "serif/mapping/variation.hpp"
#include "serif/utils/user_options/mapping/domain_mapper.hpp"
namespace serif::mapping {
class DomainMapper {
public:
DomainMapper(const DomainMapper&) = delete;
DomainMapper& operator=(const DomainMapper&) = delete;
~DomainMapper() = default;
DomainMapper(DomainMapper&&) = default;
DomainMapper& operator=(DomainMapper&&) = default;
DomainMapper(
utils::user_options::DomainMapperOptions options,
std::unique_ptr<const compactification::ExteriorDomainMap> exterior_map
);
std::unique_ptr<DomainMapper> clone() const;
void set_geometry_cache_enabled(bool enabled) const;
[[nodiscard]] bool is_geometry_cache_enabled() const noexcept;
void clear_geometry_cache() const;
[[nodiscard]] DMGeometryCacheStatistics get_geometry_cache_statistics() const;
[[nodiscard]] const DMFixedMeshGeometry* get_fixed_mesh_geometry(mfem::ElementTransformation& transformation, const mfem::IntegrationPoint& ip) const;
[[nodiscard]] MappingStatus evaluate_point(
const ElementMappingData& element_data,
mfem::ElementTransformation &transformation,
const mfem::IntegrationPoint& ip,
DMWorkspace& workspace,
MappingPointContext& context
) const;
[[nodiscard]] MappingStatus evaluate_volume(
const ElementMappingData& element_data,
mfem::ElementTransformation& transformation,
const mfem::IntegrationPoint& ip,
DMWorkspace& workspace,
VolumeMappingContext& context,
const DMReferenceDisplacementSample *displacement_sample = nullptr
) const;
[[nodiscard]] MappingStatus evaluate_face(
const ElementMappingData& element_data,
mfem::FaceElementTransformations& transformation,
FaceElementSide side,
const mfem::IntegrationPoint& ip,
DMWorkspace& workspace,
FaceMappingContext& context
) const;
[[nodiscard]] MappingStatus 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;
[[nodiscard]] MappingStatus 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;
[[nodiscard]] MappingStatus 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;
[[nodiscard]] bool is_compactified_element(const mfem::ElementTransformation& transformation) const noexcept;
[[nodiscard]] int get_dimension() const noexcept;
[[nodiscard]] const compactification::ExteriorDomainMap& get_exterior_map() const noexcept;
private:
[[nodiscard]] MappingStatus 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 = nullptr
) const;
void validate_element_data(
const ElementMappingData& element_data
) const;
void 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;
[[nodiscard]] MappingStatus 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;
[[nodiscard]] MappingStatus 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;
// These are static since they do not use any class members or workspace. Really these could and
// maybe should be free functions...
[[nodiscard]] static mfem::ElementTransformation& select_face_element_transformation(
const mfem::FaceElementTransformations& transformations,
FaceElementSide side
);
[[nodiscard]] static const mfem::IntegrationPoint& select_face_element_integration_point(
const mfem::FaceElementTransformations& transformations,
FaceElementSide side
);
private:
mutable std::unique_ptr<DMGeometryCache> m_geometryCache;
utils::user_options::DomainMapperOptions m_options;
std::unique_ptr<const compactification::ExteriorDomainMap> m_exterior_map;
};
}

93
src/include/serif/mapping/evaluator.hpp
View File

@@ -1,93 +0,0 @@
#pragma once
#include <memory>
#include "mfem/fem/gridfunc.hpp"
#include "mfem/fem/eltrans.hpp"
#include "mfem/fem/intrules.hpp"
#include "mfem/linalg/vector.hpp"
#include "mfem/general/array.hpp"
#include "serif/mapping/domain_mapper.hpp"
#include "serif/mapping/context.hpp"
#include "serif/mapping/status.hpp"
namespace serif::mapping {
// Note this is NOT thread safe.
class GridFunctionMappingEvaluator {
public:
/*
* The evaluator references the supplied grid functions and caches copies of
* their element-local DOFs. You need to make sure you call InvalidateCache() or Refresh() after either
* grid function's values are modified. Changes to the finite-element-space sequence
* should be detected automatically.
*
* This object owns mutable workspace and cache state and *is not* thread-safe.
*/
GridFunctionMappingEvaluator(
const DomainMapper& mapper,
const mfem::GridFunction& displacement,
const mfem::GridFunction& compactification_coordinate
);
void invalidate_cache() noexcept;
void refresh();
[[nodiscard]] MappingStatus evaluate_point(
mfem::ElementTransformation &transformation,
const mfem::IntegrationPoint &integration_point,
MappingPointContext &context
);
[[nodiscard]] MappingStatus evaluate_volume(
mfem::ElementTransformation &transformation,
const mfem::IntegrationPoint &integration_point,
VolumeMappingContext &context
);
[[nodiscard]] MappingStatus evaluate_face(
mfem::FaceElementTransformations &transformation,
FaceElementSide side,
const mfem::IntegrationPoint &integration_point,
FaceMappingContext &context
);
[[nodiscard]] VolumeQuadratureContext get_quadrature_context(
mfem::ElementTransformation &transformation,
const mfem::IntegrationPoint &integration_point
);
[[nodiscard]] FaceQuadratureContext get_face_quadrature_context(
mfem::FaceElementTransformations &transformation,
const mfem::IntegrationPoint &integration_point,
FaceElementSide side = FaceElementSide::ELEMENT_1
);
void get_physical_point(
mfem::ElementTransformation &transformation,
const mfem::IntegrationPoint &integration_point,
mfem::Vector &physical_position
);
private:
void validate_field_bindings() const;
[[nodiscard]] bool invalidate_for_changed_spaces();
void load_element(int element_id);
private:
const DomainMapper &m_mapper;
const mfem::GridFunction &m_displacement;
const mfem::GridFunction &m_compactification_coordinate;
const mfem::FiniteElementSpace *m_displacement_space;
const mfem::FiniteElementSpace *m_compactification_space;
long m_displacement_space_sequence;
long m_compactification_space_sequence;
DMWorkspace m_workspace;
mfem::Array<int> m_displacement_dofs;
mfem::Array<int> m_compactification_dofs;
mfem::Vector m_element_displacement;
mfem::Vector m_element_compactification;
std::unique_ptr<ElementDisplacementData> m_displacement_data;
std::unique_ptr<ElementCompactificationData> m_compactification_data;
int m_cached_element_id{-1};
};
}

4
src/include/serif/mapping/mapping.hpp Normal file
View File

@@ -0,0 +1,4 @@
#pragma once
#include "serif/backend/mfem/mapping/domain_mapper.hpp"
#include "serif/backend/mfem/mapping/transformations.hpp"

58
src/include/serif/mapping/status.hpp
View File

@@ -1,58 +0,0 @@
#pragma once
#include <unordered_map>
#include <cstdint>
#include <utility>
#include "serif/utils/error/codes.hpp"
#include "serif/meta/enums.hpp"
namespace serif::mapping
{
enum class MappingStatus : uint16_t {
VALID,
// Start this at 400 so this is parallel to the mapping error codes. There must be a better way to
// do this but it is not worth the time to sort out now
INVALID_DIMENSION = 400,
NON_FINITE_MAPPING_INPUT,
INVALID_REFERENCE_RADIUS,
AT_COMPACTIFIED_INFINITY,
OUTSIDE_REFERENCE_DOMAIN,
NON_FINITE_MAPPING_RESULT,
NON_POSITIVE_DETERMINATE
};
static auto MappingStatus_enum_string_map = meta::enum_to_enum_string_map<MappingStatus>();
static auto MappingStatus_string_enum_map = meta::enum_to_string_enum_map<MappingStatus>();
inline std::unordered_map<MappingStatus, utils::errors::SERiFErrorCode> make_MappingStatus_SERiFErrorCode_map() {
std::unordered_map<MappingStatus, utils::errors::SERiFErrorCode> ret;
for (const auto& [MS, MS_name] : MappingStatus_enum_string_map) {
if (std::to_underlying(MS) < 400) {
continue; // Skip VALID since it is not in the error code set
}
ret.insert({MS, utils::errors::SERiFErrorCodeNameToTypeMap.at(MS_name)});
}
return ret;
}
inline std::unordered_map<utils::errors::SERiFErrorCode, MappingStatus> make_SERIFErrorCode_MappingStatus_map() {
std::unordered_map<utils::errors::SERiFErrorCode, MappingStatus> ret;
for (const auto& [MS, MS_name] : MappingStatus_enum_string_map) {
if (std::to_underlying(MS) < 400) {
continue; // Skip VALID since it is not in the error code set
}
ret.insert({utils::errors::SERiFErrorCodeNameToTypeMap.at(MS_name), MS});
}
return ret;
}
// We make these parallel lists so to ensure that we can lookup between SERiF error codes
// and equivalent MappingStatus. Note again this is ugly and bad and I am ashamed of this code....
// anyhoodles there is likely a better way but pick your battles amiright?
static auto MappingStatus_SERiFErrorCode_map = make_MappingStatus_SERiFErrorCode_map();
static auto SERiFErrorCode_MappingStatus_map = make_SERIFErrorCode_MappingStatus_map();
}

87
src/include/serif/mapping/transformations.hpp
View File

@@ -1,87 +0,0 @@
#pragma once
#include "serif/mapping/variation.hpp"
#include "serif/mapping/context.hpp"
namespace serif::mapping {
void map_hdiv_flux_to_physical(
const MappingPointContext& context,
const mfem::Vector& reference_flux,
mfem::Vector& physical_flux
);
void map_physical_flux_to_hdiv_reference(
const MappingPointContext& context,
const mfem::Vector& physical_flux,
mfem::Vector& reference_flux
);
void map_reference_gradient_to_physical(
const MappingPointContext& context,
const mfem::Vector& reference_gradient,
mfem::Vector& physical_gradient
);
void map_physical_gradient_to_reference(
const MappingPointContext& context,
const mfem::Vector& physical_gradient,
mfem::Vector& reference_gradient
);
void map_reference_vector_gradient_to_physical(
const MappingPointContext& context,
const mfem::DenseMatrix& reference_vector_gradient,
mfem::DenseMatrix& physical_vector_gradient
);
void map_physical_vector_gradient_to_reference(
const MappingPointContext& context,
const mfem::DenseMatrix& physical_vector_gradient,
mfem::DenseMatrix& reference_vector_gradient
);
[[nodiscard]] double map_hdiv_divergence_to_physical(
const MappingPointContext& context,
double reference_divergence
);
void compute_hdiv_mass_tensor(
const MappingPointContext& context,
mfem::DenseMatrix& mass_tensor
);
void compute_scalar_diffusion_tensor(
const MappingPointContext& context,
mfem::DenseMatrix& diffusion_tensor
);
void map_hcurl_field_to_physical(
const MappingPointContext& context,
const mfem::Vector& reference_field,
mfem::Vector& physical_field
);
void map_physical_field_to_hcurl_reference(
const MappingPointContext& context,
const mfem::Vector& physical_field,
mfem::Vector& reference_field
);
void map_hcurl_curl_to_physical(
const MappingPointContext& context,
const mfem::Vector &reference_curl,
mfem::Vector &physical_curl
);
void map_physical_curl_to_hcurl_reference(
const MappingPointContext& context,
const mfem::Vector& physical_curl,
mfem::Vector& reference_curl
);
void compute_hdiv_mass_tensor_variation(
const MappingPointContext& context,
const MappingPointVariation& variation,
mfem::DenseMatrix& mass_tensor_variation
);
}

29
src/include/serif/mapping/variation.hpp
View File

@@ -1,29 +0,0 @@
#pragma once
#include "mfem/linalg/densemat.hpp"
#include "mfem/linalg/vector.hpp"
namespace serif::mapping {
struct MappingPointVariation {
mfem::Vector displacement_variation;
mfem::Vector physical_position_variation;
mfem::DenseMatrix displacement_jacobian_variation;
mfem::DenseMatrix mapping_jacobian_variation;
mfem::DenseMatrix inverse_mapping_jacobian_variation;
double mapping_determinate_variation{0.0};
};
struct VolumeMappingVariation {
MappingPointVariation mapping;
mfem::DenseMatrix inverse_element_jacobian_variation;
double weight_variation{0.0};
};
struct FaceMappingVariation {
MappingPointVariation mapping;
mfem::Vector physical_normal_variation;
double physical_surface_weight_variation{0.0};
double normal_flux_scale_variation{0.0};
};
}

28
src/include/serif/mapping/volume_cache.hpp
View File

@@ -1,28 +0,0 @@
#pragma once
#include "mfem/linalg/densemat.hpp"
#include "serif/mapping/context.hpp"
#include <vector>
namespace serif::mapping {
class VolumeMappingCache {
public:
void set_size(int point_count, int dimension);
void store(int point, VolumeMappingContext& context);
void load(int point, VolumeMappingContext& context) const;
void load_inverse_jacobian(int point, mfem::DenseMatrix& inverse) const;
[[nodiscard]] int get_point_count() const noexcept;
[[nodiscard]] int get_dimension() const noexcept;
private:
void validate_point(int point) const;
[[nodiscard]] const double* get_point_data(int point) const;
std::vector<double> m_data;
int m_point_count{0};
int m_dimension{0};
int m_point_stride{0};
};
}