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

docs(mapping): added doxygen comments

Browse Source
This commit is contained in:
Emily Boudreaux
2026-10-12 13:58:37 -04:00
parent 546dc42d3e
commit c440134059
3 changed files with 348 additions and 73 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -1,5 +1,14 @@
#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"
@@ -12,10 +21,54 @@
#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)>;
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,
@@ -23,12 +76,20 @@ namespace serif::mapping {
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;
GridFunctionMappingEvaluator m_mapping; ///< The mapping evaluator that computes physical positions from reference coordinates.
};
}

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

@@ -1,46 +1,92 @@
#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"
// TODO: There is likely some major readability improvements that could be gained by naming these things
// less confusing names
// TODO [Emily]: Put more comments in this file describing what the various terms are
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;
double mapping_determinate;
double physical_integration_weight;
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;
double ds;
double v_dot_n_scale;
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;
mfem::Vector displaced_position;
mfem::Vector physical_position;
mfem::DenseMatrix displacement_jacobian;
mfem::DenseMatrix mapping_jacobian;
mfem::DenseMatrix inverse_mapping_jacobian;
double mapping_determinate{0.0};
bool compactified{false};
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;
VolumeQuadratureContext quadrature;
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;
FaceQuadratureContext quadrature;
mfem::Vector reference_normal;
double reference_surface_weight{0.0};
double physical_surface_weight{0.0};
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.
};
}

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

@@ -1,4 +1,18 @@
#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"
@@ -15,50 +29,148 @@
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,
ELEMENT_2
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:
ElementDisplacementData(const mfem::FiniteElement& element, mfem::Vector& displacement_dofs, mfem::Ordering::Type ordering = mfem::Ordering::byNODES);
/**
* @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;
mfem::DenseMatrix m_dof_matrix;
int m_dimension;
mfem::Ordering::Type m_ordering;
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).
};
// Maybe this should be moved into the compactification module
/**
* @struct CompactificationPointData
*
* @brief A structure representing the compactification point data, including the compactification coordinate and coordinate gradient.
*/
struct CompactificationPointData {
double compactification_coordinate{0.0};
mfem::Vector coordinate_gradient;
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:
ElementCompactificationData(const mfem::FiniteElement& element, const mfem::Vector& dofs);
/**
* @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;
mfem::Vector m_dofs;
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 {
@@ -66,62 +178,118 @@ namespace serif::mapping {
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;
mfem::DenseMatrix m_reference_dshape;
mfem::DenseMatrix m_mesh_dshape;
mfem::DenseMatrix m_reference_field_jacobian;
mfem::Vector m_field_value;
mfem::DenseMatrix m_field_jacobian;
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;
mfem::DenseMatrix m_compactification_dshape;
CompactificationPointData m_compactification_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;
mfem::Vector m_mapped_normal;
mfem::DenseMatrix m_full_element_jacobian;
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;
mfem::DenseMatrix m_matrix_temp_1;
mfem::DenseMatrix m_matrix_temp_2;
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;
compactification::ExteriorMapVariation m_exterior_variation;
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;
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;
mfem::DenseMatrix jacobian;
mfem::DenseMatrix inverse_jacobian;
double weight{0.0};
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};
std::uint64_t misses{0};
std::size_t points{0};
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;
const mfem::DenseMatrix& gradient;
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};
const mfem::Mesh* mesh{nullptr};
long sequence{-1}; // We explicitly use the non explicltiy sized versions here since that what MFEM uses fort these members. Alternative we could use a int64_t and verify that on the compiled system that is the same as a long
long nodesSequence{-1};
std::map<std::tuple<int, double, double, double>, DMFixedMeshGeometry> points; // Currently this is an regular map since to use an unordered map we would need to add a std::hash specialization for tuple<int, double, double, double>. Some small future optimizations may be available if we implement a fast hashing algorithm for that type.
uint64_t hits{0};
uint64_t misses{0};
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.
};
}