TetrahedronHyperelasticityFEMForceField
Generic Hyperelastic Tetrahedral finite elements. Compute Finite Element forces based on tetrahedral elements.
The `TetrahedronHyperelasticityFEMForceField` computes internal forces for hyperelastic materials using tetrahedral finite elements, supporting various material models such as Neo-Hookean, Ogden, and StVenant-Kirchhoff.
- module
- Sofa.Component.SolidMechanics.FEM.HyperElastic
- namespace
- sofa::component::solidmechanics::fem::hyperelastic
- include
- sofa/component/solidmechanics/fem/hyperelastic/TetrahedronHyperelasticityFEMForceField.h
- inherits
-
- ForceField
- templates
-
- sofa::defaulttype::Vec3Types
- description
Governing Equations and Operators
The TetrahedronHyperelasticityFEMForceField component computes the internal forces in a hyperelastic material using tetrahedral finite elements. It contributes to the global stiffness matrix (K), tangent stiffness matrix, and internal force vector (fint) of the system.
The primary governing equation for this component is derived from the principle of virtual work:
\[ \delta W = \int_\Omega \sigma : \varepsilon(\delta u) d\Omega = 0 \]where \(\sigma\) is the Cauchy stress tensor, \(\varepsilon(u)\) is the strain measure derived from displacement field \(u\), and \(\delta u\) represents a virtual displacement.
The internal force vector fint and stiffness matrix K are computed by discretizing the domain into tetrahedral elements. For each tetrahedron, the deformation gradient \(F = \nabla u + I\), right Cauchy-Green deformation tensor \(C = F^T F\), and strain energy density function \(W(C)\) are used to derive the stress:
\[ \sigma = 2 \frac{\partial W}{\partial C} : F^{-1} \]The internal force contribution for a tetrahedron \(T\) is given by:
\[ f^e_{int}(u) = -J_e \sum_{i=0}^{3} N_i(u) \cdot P_i(C, F^{-1}) \]where \(N_i\) are the shape functions, and \(P_i = 2 J_e \frac{\partial W}{\partial C} : B^T F^{-1}\). Here, \(B\) is a strain-displacement matrix.
Constitutive and Kinematic Laws
The constitutive law depends on the selected hyperelastic material model. The component supports various models such as:
- Neo-Hookean: \[ W_{NH} = \frac{1}{2}\mu (I_1 - 3) + K(I_3^{\alpha} - 1) \]
- Ogden: \[ W_{O} = \sum_i^m \frac{2}{\alpha_i}\mu_i (\lambda_i^{\alpha_i}-1) + K(I_3 - 1) \]
- StVenant–Kirchhoff: \[ W_{SK} = \frac{1}{2}\lambda (I_1 - 3) + G(I_4 - 1) \]
where \(I_1, I_3\) are the first and third invariants of C.
where \(\lambda_i\) are the principal stretches, and \(K\), \(\mu_i\), and \(\alpha_i\) are material parameters.
where \(G\), and \(\lambda\) are Lame's constants.
Role in the FEM Pipeline
- Assembly Phase: Assembles global stiffness matrix and internal forces for each tetrahedron element.
- Time Integration: Provides fint, which is used in implicit time integration schemes such as backward Euler.
- Nonlinear Resolution: Computes the nonlinear residual \(R(x_{n+1}) = M \frac{x_{n+1} - x_n}{\Delta t} - f^e_{int}(x_{n+1}) - f_{ext}\).
- Linear Resolution: Computes the tangent stiffness matrix \(K(x) = \frac{1}{\Delta t} M - \Delta t K(x)\), which is required for Newton–Raphson iterations.
Numerical Methods and Discretization Choices
- FEM Discretization: Uses tetrahedral elements with linear interpolation functions. Deformation gradient and strain energy calculations are performed for each element based on the nodal displacements.
- Numerical Integration: Gauss quadrature is employed within each tetrahedron to integrate internal forces and stiffness contributions.
Integration into Variational/Lagrangian Mechanics Framework
The TetrahedronHyperelasticityFEMForceField component fits into the broader variational mechanics framework by discretizing the principle of virtual work. It derives from Lagrangian mechanics, where the kinetic and potential energies are defined through the strain energy density function for hyperelastic materials:
The Lagrangian \(L\) is then used to derive the equations of motion:
\[ M \ddot{u} + f_{int}(u) = f_{ext} \]This system is solved using implicit time integration and nonlinear solvers, ensuring stability and consistency in simulations involving hyperelastic materials.
Data Fields
| Name | Type | Default | Help |
|---|---|---|---|
d_stiffnessMatrixRegularizationWeight |
bool | |
Regularization of the Stiffness Matrix (between true or false) |
d_materialName |
sofa::helper::OptionsGroup | |
the name of the material to be used. Possible options are: 'ArrudaBoyce', 'Costa', 'MooneyRivlin', 'NeoHookean', 'Ogden', 'StVenantKirchhoff', 'VerondaWestman', 'StableNeoHookean' |
d_parameterSet |
SetParameterArray | |
The global parameters specifying the material |
d_anisotropySet |
SetAnisotropyDirectionArray | |
The global directions of anisotropy of the material: vector containing anisotropic directions. The vector size is 0 if the material is isotropic, 1 if it is transversely isotropic and 2 for orthotropic materials |
Links
| Name | Type | Help |
|---|---|---|
l_topology |
link to the topology container |
Methods
void
setMaterialName
(int materialName)
void
setparameter
(const SetParameterArray & param)
void
setdirection
(const SetAnisotropyDirectionArray & direction)
void
createTetrahedronRestInformation
(Index , TetrahedronRestInformation & t, const Tetrahedron & , const int & , const int & )
void
init
()
void
addForce
(const core::MechanicalParams * mparams, DataVecDeriv & d_f, const DataVecCoord & d_x, const DataVecDeriv & d_v)
void
addDForce
(const core::MechanicalParams * mparams, DataVecDeriv & d_df, const DataVecDeriv & d_dx)
SReal
getPotentialEnergy
(const core::MechanicalParams * , const DataVecCoord & )
void
addKToMatrix
(sofa::linearalgebra::BaseMatrix * mat, SReal k, unsigned int & offset)
void
buildStiffnessMatrix
(core::behavior::StiffnessMatrix * matrix)
void
buildDampingMatrix
(core::behavior::DampingMatrix * )
void
draw
(const core::visual::VisualParams * vparams)
void
computeBBox
(const core::ExecParams * params, bool onlyVisible)
Mat<3, 3, SReal>
getPhi
(int tetrahedronIndex)
void
testDerivatives
()
void
updateTangentMatrix
()
void
instantiateMaterial
()
{
"name": "TetrahedronHyperelasticityFEMForceField",
"namespace": "sofa::component::solidmechanics::fem::hyperelastic",
"module": "Sofa.Component.SolidMechanics.FEM.HyperElastic",
"include": "sofa/component/solidmechanics/fem/hyperelastic/TetrahedronHyperelasticityFEMForceField.h",
"doc": "Generic Hyperelastic Tetrahedral finite elements.\n\nCompute Finite Element forces based on tetrahedral elements.",
"inherits": [
"ForceField"
],
"templates": [
"sofa::defaulttype::Vec3Types"
],
"data_fields": [
{
"name": "d_stiffnessMatrixRegularizationWeight",
"type": "bool",
"xmlname": "matrixRegularization",
"help": "Regularization of the Stiffness Matrix (between true or false)"
},
{
"name": "d_materialName",
"type": "sofa::helper::OptionsGroup",
"xmlname": "materialName",
"help": "the name of the material to be used. Possible options are: 'ArrudaBoyce', 'Costa', 'MooneyRivlin', 'NeoHookean', 'Ogden', 'StVenantKirchhoff', 'VerondaWestman', 'StableNeoHookean'"
},
{
"name": "d_parameterSet",
"type": "SetParameterArray",
"xmlname": "ParameterSet",
"help": "The global parameters specifying the material"
},
{
"name": "d_anisotropySet",
"type": "SetAnisotropyDirectionArray",
"xmlname": "AnisotropyDirections",
"help": "The global directions of anisotropy of the material: vector containing anisotropic directions. The vector size is 0 if the material is isotropic, 1 if it is transversely isotropic and 2 for orthotropic materials"
}
],
"links": [
{
"name": "l_topology",
"target": "BaseMeshTopology",
"kind": "single",
"xmlname": "topology",
"help": "link to the topology container"
}
],
"methods": [
{
"name": "setMaterialName",
"return_type": "void",
"params": [
{
"name": "materialName",
"type": "int"
}
],
"is_virtual": false,
"is_pure_virtual": false,
"is_static": false,
"access": "public"
},
{
"name": "setparameter",
"return_type": "void",
"params": [
{
"name": "param",
"type": "const SetParameterArray &"
}
],
"is_virtual": false,
"is_pure_virtual": false,
"is_static": false,
"access": "public"
},
{
"name": "setdirection",
"return_type": "void",
"params": [
{
"name": "direction",
"type": "const SetAnisotropyDirectionArray &"
}
],
"is_virtual": false,
"is_pure_virtual": false,
"is_static": false,
"access": "public"
},
{
"name": "createTetrahedronRestInformation",
"return_type": "void",
"params": [
{
"name": "",
"type": "Index"
},
{
"name": "t",
"type": "TetrahedronRestInformation &"
},
{
"name": "",
"type": "const Tetrahedron &"
},
{
"name": "",
"type": "const int &"
},
{
"name": "",
"type": "const int &"
}
],
"is_virtual": false,
"is_pure_virtual": false,
"is_static": false,
"access": "public"
},
{
"name": "init",
"return_type": "void",
"params": [],
"is_virtual": false,
"is_pure_virtual": false,
"is_static": false,
"access": "public"
},
{
"name": "addForce",
"return_type": "void",
"params": [
{
"name": "mparams",
"type": "const core::MechanicalParams *"
},
{
"name": "d_f",
"type": "DataVecDeriv &"
},
{
"name": "d_x",
"type": "const DataVecCoord &"
},
{
"name": "d_v",
"type": "const DataVecDeriv &"
}
],
"is_virtual": false,
"is_pure_virtual": false,
"is_static": false,
"access": "public"
},
{
"name": "addDForce",
"return_type": "void",
"params": [
{
"name": "mparams",
"type": "const core::MechanicalParams *"
},
{
"name": "d_df",
"type": "DataVecDeriv &"
},
{
"name": "d_dx",
"type": "const DataVecDeriv &"
}
],
"is_virtual": false,
"is_pure_virtual": false,
"is_static": false,
"access": "public"
},
{
"name": "getPotentialEnergy",
"return_type": "SReal",
"params": [
{
"name": "",
"type": "const core::MechanicalParams *"
},
{
"name": "",
"type": "const DataVecCoord &"
}
],
"is_virtual": false,
"is_pure_virtual": false,
"is_static": false,
"access": "public"
},
{
"name": "addKToMatrix",
"return_type": "void",
"params": [
{
"name": "mat",
"type": "sofa::linearalgebra::BaseMatrix *"
},
{
"name": "k",
"type": "SReal"
},
{
"name": "offset",
"type": "unsigned int &"
}
],
"is_virtual": false,
"is_pure_virtual": false,
"is_static": false,
"access": "public"
},
{
"name": "buildStiffnessMatrix",
"return_type": "void",
"params": [
{
"name": "matrix",
"type": "core::behavior::StiffnessMatrix *"
}
],
"is_virtual": false,
"is_pure_virtual": false,
"is_static": false,
"access": "public"
},
{
"name": "buildDampingMatrix",
"return_type": "void",
"params": [
{
"name": "",
"type": "core::behavior::DampingMatrix *"
}
],
"is_virtual": false,
"is_pure_virtual": false,
"is_static": false,
"access": "public"
},
{
"name": "draw",
"return_type": "void",
"params": [
{
"name": "vparams",
"type": "const core::visual::VisualParams *"
}
],
"is_virtual": false,
"is_pure_virtual": false,
"is_static": false,
"access": "public"
},
{
"name": "computeBBox",
"return_type": "void",
"params": [
{
"name": "params",
"type": "const core::ExecParams *"
},
{
"name": "onlyVisible",
"type": "bool"
}
],
"is_virtual": false,
"is_pure_virtual": false,
"is_static": false,
"access": "public"
},
{
"name": "getPhi",
"return_type": "Mat<3, 3, SReal>",
"params": [
{
"name": "tetrahedronIndex",
"type": "int"
}
],
"is_virtual": false,
"is_pure_virtual": false,
"is_static": false,
"access": "public"
},
{
"name": "testDerivatives",
"return_type": "void",
"params": [],
"is_virtual": false,
"is_pure_virtual": false,
"is_static": false,
"access": "protected"
},
{
"name": "updateTangentMatrix",
"return_type": "void",
"params": [],
"is_virtual": false,
"is_pure_virtual": false,
"is_static": false,
"access": "protected"
},
{
"name": "instantiateMaterial",
"return_type": "void",
"params": [],
"is_virtual": false,
"is_pure_virtual": false,
"is_static": false,
"access": "protected"
}
],
"description": "The `TetrahedronHyperelasticityFEMForceField` is a force field component in the SOFA framework designed for simulating hyperelastic materials using tetrahedral finite elements. It inherits from `core::behavior::ForceField`, enabling it to compute and apply forces on nodes based on their position and deformation. This component supports several material models, including ArrudaBoyce, NeoHookean, MooneyRivlin, Ogden, StVenantKirchhoff, Costa, VerondaWestman, and StableNeoHookean, which can be specified via the `d_materialName` data field.\n\nThe component interacts with a topology container to obtain mesh information (tetrahedra and edges) through the `l_topology` link. This allows it to define internal data structures for each tetrahedron (`TetrahedronRestInformation`) and edge (`EdgeInformation`). The force computation is based on deformation gradients, strain energy calculations, and elasticity tensor applications specific to the chosen hyperelastic material model.\n\nThe component provides methods such as `addForce`, `addDForce`, and `buildStiffnessMatrix` for integrating with the mechanical solver. It also offers options for regularization of the stiffness matrix (`d_stiffnessMatrixRegularizationWeight`) and setting global parameters (`d_parameterSet`). Anisotropic material properties can be configured using `d_anisotropySet`. The `init` method initializes internal data structures based on topology information and material settings.\n\nThis component is useful for simulating soft tissue mechanics, where hyperelastic models are commonly employed to capture nonlinear deformations. It is part of the solid mechanics simulation workflow in SOFA and can be used alongside other components for comprehensive biomechanical modeling.",
"maths": "<h2>Governing Equations and Operators</h2>\n\n<p>The <code>TetrahedronHyperelasticityFEMForceField</code> component computes the internal forces in a hyperelastic material using tetrahedral finite elements. It contributes to the global stiffness matrix (K), tangent stiffness matrix, and internal force vector (f<sub>int</sub>) of the system.</p>\n\n<p>The primary governing equation for this component is derived from the principle of virtual work:</p>\n\n\\[ \\delta W = \\int_\\Omega \\sigma : \\varepsilon(\\delta u) d\\Omega = 0 \\]\n\n<p>where <span class=\"math-inline\">\\(\\sigma\\)</span> is the Cauchy stress tensor, <span class=\"math-inline\">\\(\\varepsilon(u)\\)</span> is the strain measure derived from displacement field <span class=\"math-inline\">\\(u\\)</span>, and <span class=\"math-inline\">\\(\\delta u\\)</span> represents a virtual displacement.</p>\n\n<p>The internal force vector f<sub>int</sub> and stiffness matrix K are computed by discretizing the domain into tetrahedral elements. For each tetrahedron, the deformation gradient <span class=\"math-inline\">\\(F = \\nabla u + I\\)</span>, right Cauchy-Green deformation tensor <span class=\"math-inline\">\\(C = F^T F\\)</span>, and strain energy density function <span class=\"math-inline\">\\(W(C)\\)</span> are used to derive the stress:</p>\n\n\\[ \\sigma = 2 \\frac{\\partial W}{\\partial C} : F^{-1} \\]\n\n<p>The internal force contribution for a tetrahedron <span class=\"math-inline\">\\(T\\)</span> is given by:</p>\n\n\\[\n f^e_{int}(u) = -J_e \\sum_{i=0}^{3} N_i(u) \\cdot P_i(C, F^{-1})\n\\]\n\n<p>where <span class=\"math-inline\">\\(N_i\\)</span> are the shape functions, and <span class=\"math-inline\">\\(P_i = 2 J_e \\frac{\\partial W}{\\partial C} : B^T F^{-1}\\)</span>. Here, <span class=\"math-inline\">\\(B\\)</span> is a strain-displacement matrix.</p>\n\n<h3>Constitutive and Kinematic Laws</h3>\n\n<p>The constitutive law depends on the selected hyperelastic material model. The component supports various models such as:</p>\n<ul>\n<li><strong>Neo-Hookean:</strong></li>\n\\[\n W_{NH} = \\frac{1}{2}\\mu (I_1 - 3) + K(I_3^{\\alpha} - 1)\n\\]\n\n<p>where <span class=\"math-inline\">\\(I_1, I_3\\)</span> are the first and third invariants of C.</p>\n<li><strong>Ogden:</strong></li>\n\\[\n W_{O} = \\sum_i^m \\frac{2}{\\alpha_i}\\mu_i (\\lambda_i^{\\alpha_i}-1) + K(I_3 - 1)\n\\]\n\n<p>where <span class=\"math-inline\">\\(\\lambda_i\\)</span> are the principal stretches, and <span class=\"math-inline\">\\(K\\)</span>, <span class=\"math-inline\">\\(\\mu_i\\)</span>, and <span class=\"math-inline\">\\(\\alpha_i\\)</span> are material parameters.</p>\n<li><strong>StVenant–Kirchhoff:</strong></li>\n\\[\n W_{SK} = \\frac{1}{2}\\lambda (I_1 - 3) + G(I_4 - 1)\n\\]\n\n<p>where <span class=\"math-inline\">\\(G\\)</span>, and <span class=\"math-inline\">\\(\\lambda\\)</span> are Lame's constants.</p>\n</ul>\n\n<h2>Role in the FEM Pipeline</h2>\n\n<ul>\n<li><strong>Assembly Phase:</strong> Assembles global stiffness matrix and internal forces for each tetrahedron element.</li>\n<li><strong>Time Integration:</strong> Provides f<sub>int</sub>, which is used in implicit time integration schemes such as backward Euler.</li>\n<li><strong>Nonlinear Resolution:</strong> Computes the nonlinear residual <span class=\"math-inline\">\\(R(x_{n+1}) = M \\frac{x_{n+1} - x_n}{\\Delta t} - f^e_{int}(x_{n+1}) - f_{ext}\\)</span>.</li>\n<li><strong>Linear Resolution:</strong> Computes the tangent stiffness matrix <span class=\"math-inline\">\\(K(x) = \\frac{1}{\\Delta t} M - \\Delta t K(x)\\)</span>, which is required for Newton–Raphson iterations.</li>\n</ul>\n\n<h3>Numerical Methods and Discretization Choices</h3>\n\n<ul>\n<li><strong>FEM Discretization:</strong> Uses tetrahedral elements with linear interpolation functions. Deformation gradient and strain energy calculations are performed for each element based on the nodal displacements.</li>\n<li><strong>Numerical Integration:</strong> Gauss quadrature is employed within each tetrahedron to integrate internal forces and stiffness contributions.</li>\n</ul>\n\n<h3>Integration into Variational/Lagrangian Mechanics Framework</h3>\n\n<p>The <code>TetrahedronHyperelasticityFEMForceField</code> component fits into the broader variational mechanics framework by discretizing the principle of virtual work. It derives from Lagrangian mechanics, where the kinetic and potential energies are defined through the strain energy density function for hyperelastic materials:</p>\n\n\\[\n T = \\frac{1}{2} \\int_\\Omega \\rho \\dot{u}^T \\dot{u} d\\Omega \\\\ V(u) = \\int_\\Omega W(C) d\\Omega\n\\]\n\n<p>The Lagrangian <span class=\"math-inline\">\\(L\\)</span> is then used to derive the equations of motion:</p>\n\n\\[\n M \\ddot{u} + f_{int}(u) = f_{ext}\n\\]\n\n<p>This system is solved using implicit time integration and nonlinear solvers, ensuring stability and consistency in simulations involving hyperelastic materials.</p>",
"abstract": "The `TetrahedronHyperelasticityFEMForceField` computes internal forces for hyperelastic materials using tetrahedral finite elements, supporting various material models such as Neo-Hookean, Ogden, and StVenant-Kirchhoff.",
"sheet": "# TetrahedronHyperelasticityFEMForceField\n\n## Overview\nThe `TetrahedronHyperelasticityFEMForceField` is a force field component in the SOFA framework designed for simulating hyperelastic materials using tetrahedral finite elements. It inherits from `core::behavior::ForceField`, enabling it to compute and apply forces on nodes based on their position and deformation. This component supports several material models, including ArrudaBoyce, NeoHookean, MooneyRivlin, Ogden, StVenantKirchhoff, Costa, VerondaWestman, and StableNeoHookean.\n\n## Mathematical Model\nThe governing equation for this component is derived from the principle of virtual work:\n\n\begin{equation}\n \\delta W = \\int_\\Omega \\sigma : \\varepsilon(\\delta u) d\\Omega = 0\n\\end{equation}\n\nwhere $\\sigma$ is the Cauchy stress tensor, $\\varepsilon(u)$ is the strain measure derived from displacement field $u$, and $\\delta u$ represents a virtual displacement. The internal force vector $f_{int}$ and stiffness matrix $K$ are computed by discretizing the domain into tetrahedral elements.\n\nFor each tetrahedron, the deformation gradient $F = \\nabla u + I$, right Cauchy-Green deformation tensor $C = F^T F$, and strain energy density function $W(C)$ are used to derive the stress:\n\n\begin{equation}\n \\sigma = 2 \\frac{\\partial W}{\\partial C} : F^{-1}\n\\end{equation}\n\nThe internal force contribution for a tetrahedron $T$ is given by:\n\n\begin{equation}\n f^e_{int}(u) = -J_e \\sum_{i=0}^{3} N_i(u) \\cdot P_i(C, F^{-1})\n\\end{equation}\n\nwhere $N_i$ are the shape functions, and $P_i = 2 J_e \\frac{\\partial W}{\\partial C} : B^T F^{-1}$. Here, $B$ is a strain-displacement matrix.\n\n### Constitutive Laws\nThe constitutive law depends on the selected hyperelastic material model. The component supports various models such as:\n\n- **Neo-Hookean:**\n \begin{equation}\n W_{NH} = \\frac{1}{2}\\mu (I_1 - 3) + K(I_3^\\alpha - 1)\n \\end{equation}\n where $I_1, I_3$ are the first and third invariants of $C$.\n\n- **Ogden:**\n \begin{equation}\n W_{O} = \\sum_i^m \\frac{2}{\\alpha_i}\\mu_i (\\lambda_i^{\\alpha_i}-1) + K(I_3 - 1)\n \\end{equation}\n where $\\lambda_i$ are the principal stretches, and $K$, $\\mu_i$, and $\\alpha_i$ are material parameters.\n\n- **StVenant–Kirchhoff:**\n \begin{equation}\n W_{SK} = \\frac{1}{2}\\lambda (I_1 - 3) + G(I_4 - 1)\n \\end{equation}\n where $G$, and $\\lambda$ are Lame's constants.\n\n## Parameters and Data\nThe significant data fields exposed by the component include:\n- `d_stiffnessMatrixRegularizationWeight`: Regularization of the Stiffness Matrix (between true or false).\n- `d_materialName`: The name of the material to be used. Possible options are: 'ArrudaBoyce', 'Costa', 'MooneyRivlin', 'NeoHookean', 'Ogden', 'StVenantKirchhoff', 'VerondaWestman', 'StableNeoHookean'.\n- `d_parameterSet`: The global parameters specifying the material.\n- `d_anisotropySet`: The global directions of anisotropy of the material: vector containing anisotropic directions. The vector size is 0 if the material is isotropic, 1 if it is transversely isotropic and 2 for orthotropic materials.\n\n## Dependencies and Connections\nThe component requires a topology container (`l_topology`) to obtain mesh information (tetrahedra and edges). This allows it to define internal data structures for each tetrahedron (`TetrahedronRestInformation`) and edge (`EdgeInformation`).\n\n## Practical Notes\n- The stiffness matrix regularization can be enabled or disabled via the `d_stiffnessMatrixRegularizationWeight` parameter.\n- Proper material parameters must be set to ensure numerical stability and accurate simulation results."
}