Skip to content

Commit 93bbb13

Browse files
Q-MinhQ-Minh
committed
Document hyper elastic potential
1 parent 46f21e0 commit 93bbb13

File tree

1 file changed

+124
-37
lines changed

1 file changed

+124
-37
lines changed

‎source/pbat/fem/HyperElasticPotential.h‎

Lines changed: 124 additions & 37 deletions
Original file line numberDiff line numberDiff line change
@@ -3,9 +3,9 @@
33
* @author Quoc-Minh Ton-That ([email protected])
44
* @brief Hyper elastic potential energy
55
* @date 2025-02-11
6-
*
6+
*
77
* @copyright Copyright (c) 2025
8-
*
8+
*
99
*/
1010

1111
#ifndef PBA_FEM_HYPER_ELASTIC_POTENTIAL_H
@@ -32,35 +32,64 @@ namespace pbat {
3232
namespace fem {
3333

3434
/**
35-
* @brief
36-
* @tparam TMesh
37-
* @tparam THyperElasticEnergy
35+
* @brief Total hyper elastic potential \f$ U(\mathbf{x}) = \int_\Omega \Psi(\mathbf{F}) d\Omega \f$
36+
*
37+
* where \f$ \mathbf{F} \f$ is the deformation gradient and \f$ \Psi \f$ is the hyper elastic energy
38+
* density.
39+
*
40+
* HyperElasticPotential's depends on a \f$ |Q| \f$-point quadrature specified by element
41+
* indices \f$ e_g \f$, quadrature weights \f$ w_g \f$, and shape function gradients \f$ \nabla
42+
* \phi_g \f$ at quadrature points. This allows users to try out different quadrature rules
43+
* seamlessly.
44+
*
45+
* @tparam TMesh Type satisfying concept CMesh
46+
* @tparam THyperElasticEnergy Type satisfying concept CHyperElasticEnergy
3847
*/
3948
template <CMesh TMesh, physics::CHyperElasticEnergy THyperElasticEnergy>
4049
struct HyperElasticPotential
4150
{
4251
public:
43-
using SelfType = HyperElasticPotential<TMesh, THyperElasticEnergy>;
44-
using MeshType = TMesh;
45-
using ElementType = typename TMesh::ElementType;
46-
using ElasticEnergyType = THyperElasticEnergy;
52+
using SelfType = HyperElasticPotential<TMesh, THyperElasticEnergy>; ///< Self type
53+
using MeshType = TMesh; ///< Mesh type
54+
using ElementType = typename TMesh::ElementType; ///< FEM element type
55+
using ElasticEnergyType = THyperElasticEnergy; ///< Hyper elastic energy density type
4756
static_assert(
4857
MeshType::kDims == ElasticEnergyType::kDims,
4958
"Embedding dimensions of mesh must match dimensionality of hyper elastic energy.");
5059

51-
static auto constexpr kDims = THyperElasticEnergy::kDims;
52-
static int constexpr kOrder = ElementType::kOrder - 1;
60+
static auto constexpr kDims = THyperElasticEnergy::kDims; ///< Number of spatial dimensions
5361

5462
SelfType& operator=(SelfType const&) = delete;
5563

64+
/**
65+
* @brief Construct a new Hyper Elastic Potential object
66+
*
67+
* @param mesh FEM mesh
68+
* @param eg \f$ |Q| \f$ array of element indices at quadrature points
69+
* @param wg \f$ |Q| \f$ array of quadrature weights
70+
* @param GNeg Shape function gradients at quadrature points. See ShapeFunctionGradients().
71+
* @param Y Young's modulus
72+
* @param nu Poisson's ratio
73+
*/
5674
HyperElasticPotential(
5775
MeshType const& mesh,
5876
Eigen::Ref<IndexVectorX const> const& eg,
5977
Eigen::Ref<VectorX const> const& wg,
6078
Eigen::Ref<MatrixX const> const& GNeg,
6179
Scalar Y,
6280
Scalar nu);
63-
81+
/**
82+
* @brief Construct a new Hyper Elastic Potential object
83+
*
84+
* @tparam TDerivedY Eigen dense expression type
85+
* @tparam TDerivednu Eigen dense expression type
86+
* @param mesh FEM mesh
87+
* @param eg \f$ |Q| \f$ array of element indices at quadrature points
88+
* @param wg \f$ |Q| \f$ array of quadrature weights
89+
* @param GNeg Shape function gradients at quadrature points. See ShapeFunctionGradients().
90+
* @param Y \f$ |Q| \f$ Young's moduli
91+
* @param nu \f$ |Q| \f$ Poisson's ratios
92+
*/
6493
template <class TDerivedY, class TDerivednu>
6594
HyperElasticPotential(
6695
MeshType const& mesh,
@@ -69,7 +98,20 @@ struct HyperElasticPotential
6998
Eigen::Ref<MatrixX const> const& GNeg,
7099
Eigen::DenseBase<TDerivedY> const& Y,
71100
Eigen::DenseBase<TDerivednu> const& nu);
72-
101+
/**
102+
* @brief Construct a new Hyper Elastic Potential object
103+
*
104+
* Eagerly computes element elasticity (and its derivatives)
105+
*
106+
* @tparam TDerived Eigen matrix expression type
107+
* @param mesh FEM mesh
108+
* @param eg \f$ |Q| \f$ array of element indices at quadrature points
109+
* @param wg \f$ |Q| \f$ array of quadrature weights
110+
* @param GNeg Shape function gradients at quadrature points. See ShapeFunctionGradients().
111+
* @param x \f$ d \times n \f$ matrix of deformed nodal positions
112+
* @param Y Young's modulus
113+
* @param nu Poisson's ratio
114+
*/
73115
template <class TDerived>
74116
HyperElasticPotential(
75117
MeshType const& mesh,
@@ -79,7 +121,22 @@ struct HyperElasticPotential
79121
Eigen::MatrixBase<TDerived> const& x,
80122
Scalar Y,
81123
Scalar nu);
82-
124+
/**
125+
* @brief Construct a new Hyper Elastic Potential object
126+
*
127+
* Eagerly computes element elasticity (and its derivatives)
128+
*
129+
* @tparam TDerivedx Eigen matrix expression type
130+
* @tparam TDerivedY Eigen dense expression type
131+
* @tparam TDerivednu Eigen dense expression type
132+
* @param mesh FEM mesh
133+
* @param eg \f$ |Q| \f$ array of element indices at quadrature points
134+
* @param wg \f$ |Q| \f$ array of quadrature weights
135+
* @param GNeg Shape function gradients at quadrature points. See ShapeFunctionGradients().
136+
* @param x \f$ d \times n \f$ matrix of deformed nodal positions
137+
* @param Y \f$ |Q| \f$ Young's moduli
138+
* @param nu \f$ |Q| \f$ Poisson's ratios
139+
*/
83140
template <class TDerivedx, class TDerivedY, class TDerivednu>
84141
HyperElasticPotential(
85142
MeshType const& mesh,
@@ -89,9 +146,22 @@ struct HyperElasticPotential
89146
Eigen::MatrixBase<TDerivedx> const& x,
90147
Eigen::DenseBase<TDerivedY> const& Y,
91148
Eigen::DenseBase<TDerivednu> const& nu);
92-
149+
/**
150+
* @brief Precomputes the sparsity pattern of the hessian matrix
151+
*
152+
* Enables parallel sparse hessian assembly in all future operations.
153+
*/
93154
void PrecomputeHessianSparsity();
94-
155+
/**
156+
* @brief Computes the element elasticity and its derivatives at the given shape
157+
*
158+
* @tparam TDerived Eigen matrix expression type
159+
* @param x \f$ d \times n \f$ matrix of deformed nodal positions
160+
* @param bWithGradient Compute gradient
161+
* @param bWithHessian Compute hessian
162+
* @param bUseSpdProjection Project per quadrature point hessians to nearest symmetric positive
163+
* definite (SPD) matrix
164+
*/
95165
template <class TDerived>
96166
void ComputeElementElasticity(
97167
Eigen::MatrixBase<TDerived> const& x,
@@ -103,53 +173,70 @@ struct HyperElasticPotential
103173
* @brief Applies the hessian matrix of this potential as a linear operator on x, adding result
104174
* to y.
105175
*
106-
* @tparam TDerivedIn
107-
* @tparam TDerivedOut
108-
* @param x
109-
* @param y
176+
* @tparam TDerivedIn Input matrix type
177+
* @tparam TDerivedOut Output matrix type
178+
* @param x Input matrix
179+
* @param y Output matrix
110180
*/
111181
template <class TDerivedIn, class TDerivedOut>
112182
void Apply(Eigen::MatrixBase<TDerivedIn> const& x, Eigen::DenseBase<TDerivedOut>& y) const;
113183

114184
/**
115185
* @brief Transforms this matrix-free hessian matrix representation into sparse compressed
116-
* format.
117-
* @return
186+
* column format.
187+
* @return CSCMatrix Sparse compressed column matrix representation of the hessian operator
118188
*/
119189
CSCMatrix ToMatrix() const;
120190

121191
/**
122-
* @brief Transforms this element-wise gradient representation into the global gradient.
123-
* @return
192+
* @brief Transforms this per quadrature point gradient representation into the global gradient.
193+
* @return VectorX Global gradient
124194
*/
125195
VectorX ToVector() const;
126196

127197
/**
128-
* @brief Computes the elastic potential
129-
* @return
198+
* @brief Computes the total elastic potential
199+
* @return Scalar Total elastic potential
130200
*/
131201
Scalar Eval() const;
132202

203+
/**
204+
* @brief Number of input dimensions
205+
*
206+
* Effectively the number of nodes in the system
207+
*
208+
* @return Index
209+
*/
133210
Index InputDimensions() const;
211+
/**
212+
* @brief Number of output dimensions
213+
*
214+
* Effectively the number of nodes in the system
215+
*
216+
* @return Index
217+
*/
134218
Index OutputDimensions() const;
135-
219+
/**
220+
* @brief Checks the validity of the held data
221+
*/
136222
void CheckValidState() const;
137223

138224
MeshType const& mesh; ///< The finite element mesh
139225
Eigen::Ref<IndexVectorX const>
140226
eg; ///< Maps quadrature point index g to its corresponding element e
141-
Eigen::Ref<VectorX const> wg; ///< |#quad.pts.|x# array of quadrature points on the mesh
227+
Eigen::Ref<VectorX const> wg; ///< Vector of quadrature weights \f$ w \in \mathbb{R}^{|Q|} \f$
142228
Eigen::Ref<MatrixX const>
143-
GNeg; ///< |ElementType::kNodes| x |MeshType::kDims * # element quadrature points *
144-
///< #elements| element shape function gradients
145-
146-
VectorX mug; ///< |#quad.pts.|x1 1st Lame coefficient
147-
VectorX lambdag; ///< |#quad.pts.|x1 2nd Lame coefficient
148-
MatrixX Hg; ///< |(ElementType::kNodes*kDims)| x |#quad.pts. *
149-
///< (ElementType::kNodes*kDims)| element hessian matrices at quadrature points
150-
MatrixX Gg; ///< |ElementType::kNodes*kDims| x |#quad.pts.| element gradient vectors at
229+
GNeg; ///< `|ElementType::kNodes| x |MeshType::kDims * # element quadrature points *
230+
///< # elements|` shape function gradients at quadrature points
231+
232+
VectorX mug; ///< 1st Lame coefficients \f$ \mu \in \mathbb{R}^{|Q|} \f$ at quadrature points
233+
VectorX lambdag; ///< 2nd Lame coefficients \f$ \lambda \in \mathbb{R}^{|Q|} \f$ at quadrature
234+
///< points
235+
MatrixX Hg; ///< `|(ElementType::kNodes*kDims)| x |# quad.pts. *
236+
///< ElementType::kNodes*kDims|` element hessian matrices at quadrature points
237+
MatrixX Gg; ///< `|ElementType::kNodes*kDims| x |#quad.pts.|` element gradient vectors at
151238
///< quadrature points
152-
VectorX Ug; ///< |#quad.pts.| x 1 element elastic potentials at quadrature points
239+
VectorX Ug; ///< `|# quad.pts.|` array of elastic potentials at quadrature points
153240
math::linalg::SparsityPattern GH; ///< Directed adjacency graph of hessian
154241
};
155242

0 commit comments

Comments
 (0)