diff --git a/src/include/serif/mapping/coefficients.hpp b/src/include/serif/mapping/coefficients.hpp index 4020496d..204154a0 100644 --- a/src/include/serif/mapping/coefficients.hpp +++ b/src/include/serif/mapping/coefficients.hpp @@ -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 #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; + using Func = std::function; ///< 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. }; } diff --git a/src/include/serif/mapping/context.hpp b/src/include/serif/mapping/context.hpp index 087a3afc..6d8f3f3e 100644 --- a/src/include/serif/mapping/context.hpp +++ b/src/include/serif/mapping/context.hpp @@ -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. }; } \ No newline at end of file diff --git a/src/include/serif/mapping/data.hpp b/src/include/serif/mapping/data.hpp index 29a49269..8c378f37 100644 --- a/src/include/serif/mapping/data.hpp +++ b/src/include/serif/mapping/data.hpp @@ -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, 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. 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, 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. }; }