ParallelTetrahedronFEMForceField
Parallel implementation of a linear elastic material using tetrahedral finite elements.. Parallel implementation of TetrahedronFEMForceField This implementation is the most efficient when: 1) the number of tetrahedron is large (> 1000) The following methods are executed in parallel: - addDForce - addKToMatrix
The `ParallelTetrahedronFEMForceField` simulates linear elastic materials using tetrahedral finite elements with parallel processing, particularly efficient for simulations involving more than 1000 tetrahedra.
- module
- MultiThreading
- namespace
- multithreading::component::forcefield::solidmechanics::fem::elastic
- include
- MultiThreading/component/solidmechanics/fem/elastic/ParallelTetrahedronFEMForceField.h
- inherits
-
- TaskSchedulerUser
- TetrahedronFEMForceField
- templates
-
- sofa::defaulttype::Vec3Types
The ParallelTetrahedronFEMForceField component in the SOFA framework is designed to simulate linear elastic materials using tetrahedral finite elements. This implementation leverages parallel processing for efficient computation, particularly when dealing with a large number of tetrahedra (>1000). Below is a detailed mathematical and physical description of this component.
Governing Equations and Operators
This component implements the following key operators:
- Internal Force (wzxhzdk:15)
- Stiffness Matrix (wzxhzdk:16 or wzxhzdk:17, where wzxhzdk:18 is the Jacobian of the deformation gradient)
- Residual (wzxhzdk:19)
Constitutive and Kinematic Laws Involved
The component supports two formulations:
1. Small Deformation (Linear Elasticity)
2. Corotational Formulation for Large Deformations
Strain Measures
- For small deformations, the strain tensor klzzwxh:0024 is approximated linearly.
- For large deformations, Green-Lagrange or Almansi strains are used in the corotational framework to account for rigid body rotations.
Stress Tensors and Hyperelastic Potentials
The material behavior is assumed linear elastic with a Hookean constitutive law:
$f_{int}$where wzxhzdk:8 is the fourth-order elasticity tensor, and wzxhzdk:9 is the strain measure.
For small deformations:
$K$where wzxhzdk:10 is the tangent stiffness matrix computed from the constitutive law and strain-displacement relations, and wzxhzdk:11 is the Jacobian of the deformation gradient.
For corotational deformations:
$N_i(X)$where wzxhzdk:12 incorporates the rotation transformation matrix wzxhzdk:13 to decouple rigid body rotations from deformation.
Role in the Global FEM Pipeline
Assembly Phase
- Mass and Force Construction: The internal forces are assembled into the global force vector $f_{int}$, which is part of the semi-discrete dynamical system equation:
wzxhzdk:3
- Element Operators: The stiffness matrix contributions for each tetrahedron are computed and assembled into the global stiffness matrix $K$.
Time Integration Phase
The component supports implicit time integration schemes, typically Backward Euler or Newmark-beta methods. The nonlinear residual is constructed as:
wzxhzdk:4
Nonlinear Resolution Phase
- Newton-Raphson Iteration: The nonlinear residual is solved iteratively using Newton's method. Each iteration involves linearization:
wzxhzdk:5
where wzxhzdk:14 is the Jacobian matrix given by
wzxhzdk:6
Numerical Methods and Discretization Choices
- Spatial Discretization: The domain is discretized into tetrahedral elements. Each node carries 3 degrees of freedom (DOFs) in 3D space.
- Isoparametric Mapping: Shape functions $N_i(X)$ are used to interpolate displacements within each element:
wzxhzdk:7
- Numerical Integration (Gauss Quadrature): Gauss quadrature is used for numerical integration over the reference tetrahedron.
Parallelization Strategy
The component parallelizes the computation of internal forces and stiffness matrix contributions using a task scheduler. The following methods are executed in parallel:
- addDForce
- addKToMatrix
Variational / Lagrangian Mechanics Framework Fit
This component fits into the broader variational mechanics framework by discretizing the weak (variational) form of the governing PDEs using FEM. It ensures that the linearization and assembly steps preserve physical consistency, numerical stability, and extensibility within the SOFA simulation pipeline.
Methods
void
init
()
void
addForce
(const sofa::core::MechanicalParams * mparams, DataVecDeriv & d_f, const DataVecCoord & d_x, const DataVecDeriv & d_v)
void
addDForce
(const sofa::core::MechanicalParams * mparams, DataVecDeriv & d_df, const DataVecDeriv & d_dx)
void
addKToMatrix
(sofa::linearalgebra::BaseMatrix * mat, SReal kFactor, unsigned int & offset)
void
addDForceGeneric
(VecDeriv & df, const VecDeriv & dx, Real kFactor, const int & indexedElements, Function f)
void
addDForceSmall
(VecDeriv & df, const VecDeriv & dx, Real kFactor, const int & indexedElements)
void
addDForceCorotational
(VecDeriv & df, const VecDeriv & dx, Real kFactor, const int & indexedElements)
void
drawTrianglesFromTetrahedra
(const sofa::core::visual::VisualParams * vparams, bool showVonMisesStressPerElement, bool drawVonMisesStress, const VecCoord & x, const VecReal & youngModulus, bool heterogeneous, Real minVM, Real maxVM, int vM)
{
"name": "ParallelTetrahedronFEMForceField",
"namespace": "multithreading::component::forcefield::solidmechanics::fem::elastic",
"module": "MultiThreading",
"include": "MultiThreading/component/solidmechanics/fem/elastic/ParallelTetrahedronFEMForceField.h",
"doc": "Parallel implementation of a linear elastic material using tetrahedral finite elements..\n\nParallel implementation of TetrahedronFEMForceField\nThis implementation is the most efficient when:\n1) the number of tetrahedron is large (> 1000)\nThe following methods are executed in parallel:\n- addDForce\n- addKToMatrix",
"inherits": [
"TaskSchedulerUser",
"TetrahedronFEMForceField"
],
"templates": [
"sofa::defaulttype::Vec3Types"
],
"data_fields": [],
"links": [],
"methods": [
{
"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 sofa::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 sofa::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": "addKToMatrix",
"return_type": "void",
"params": [
{
"name": "mat",
"type": "sofa::linearalgebra::BaseMatrix *"
},
{
"name": "kFactor",
"type": "SReal"
},
{
"name": "offset",
"type": "unsigned int &"
}
],
"is_virtual": false,
"is_pure_virtual": false,
"is_static": false,
"access": "public"
},
{
"name": "addDForceGeneric",
"return_type": "void",
"params": [
{
"name": "df",
"type": "VecDeriv &"
},
{
"name": "dx",
"type": "const VecDeriv &"
},
{
"name": "kFactor",
"type": "Real"
},
{
"name": "indexedElements",
"type": "const int &"
},
{
"name": "f",
"type": "Function"
}
],
"is_virtual": false,
"is_pure_virtual": false,
"is_static": false,
"access": "protected"
},
{
"name": "addDForceSmall",
"return_type": "void",
"params": [
{
"name": "df",
"type": "VecDeriv &"
},
{
"name": "dx",
"type": "const VecDeriv &"
},
{
"name": "kFactor",
"type": "Real"
},
{
"name": "indexedElements",
"type": "const int &"
}
],
"is_virtual": false,
"is_pure_virtual": false,
"is_static": false,
"access": "protected"
},
{
"name": "addDForceCorotational",
"return_type": "void",
"params": [
{
"name": "df",
"type": "VecDeriv &"
},
{
"name": "dx",
"type": "const VecDeriv &"
},
{
"name": "kFactor",
"type": "Real"
},
{
"name": "indexedElements",
"type": "const int &"
}
],
"is_virtual": false,
"is_pure_virtual": false,
"is_static": false,
"access": "protected"
},
{
"name": "drawTrianglesFromTetrahedra",
"return_type": "void",
"params": [
{
"name": "vparams",
"type": "const sofa::core::visual::VisualParams *"
},
{
"name": "showVonMisesStressPerElement",
"type": "bool"
},
{
"name": "drawVonMisesStress",
"type": "bool"
},
{
"name": "x",
"type": "const VecCoord &"
},
{
"name": "youngModulus",
"type": "const VecReal &"
},
{
"name": "heterogeneous",
"type": "bool"
},
{
"name": "minVM",
"type": "Real"
},
{
"name": "maxVM",
"type": "Real"
},
{
"name": "vM",
"type": "int"
}
],
"is_virtual": false,
"is_pure_virtual": false,
"is_static": false,
"access": "protected"
}
],
"component": {
"name": "ParallelTetrahedronFEMForceField",
"description": "The `ParallelTetrahedronFEMForceField` is a parallel implementation of the finite element method (FEM) for simulating linear elastic materials using tetrahedral elements. It extends the functionality of the standard TetrahedronFEMForceField by optimizing computational performance through parallel execution, making it particularly efficient for simulations involving large numbers of tetrahedral elements.",
"namespace": "multithreading::component::forcefield::solidmechanics::fem::elastic",
"inherits_from": [
{
"name": "sofa::component::solidmechanics::fem::elastic::TetrahedronFEMForceField",
"description": "Base class implementing the finite element method with tetrahedral elements."
},
{
"name": "sofa::simulation::TaskSchedulerUser",
"description": "Provides functionality for managing task scheduling in parallel environments."
}
],
"template_parameters": [
{
"type": "DataTypes",
"description": "Specifies the type of coordinates and derivatives used, typically instantiated with `Vec3Types` for 3D simulations."
}
],
"methods": [
{
"name": "init",
"description": "Initializes the component. It calls the base class initialization method and sets up any necessary parallel execution configurations.",
"return_type": "void",
"parameters": [],
"overridden": true,
"implementation_file": "./MultiThreading/component/solidmechanics/fem/elastic/ParallelTetrahedronFEMForceField.inl"
},
{
"name": "addForce",
"description": "Calculates the forces acting on each vertex due to linear elastic material behavior. This method is overridden from the base class and optimized for parallel execution.",
"return_type": "void",
"parameters": [
{
"type": "const sofa::core::MechanicalParams*",
"name": "mparams"
},
{
"type": "DataVecDeriv&",
"name": "d_f"
},
{
"type": "const DataVecCoord&",
"name": "d_x"
},
{
"type": "const DataVecDeriv&",
"name": "d_v"
}
],
"overridden": true,
"implementation_file": "./MultiThreading/component/solidmechanics/fem/elastic/ParallelTetrahedronFEMForceField.inl"
},
{
"name": "addDForce",
"description": "Calculates the derivative of forces with respect to displacements. This method is overridden and optimized for parallel execution using a generic function template `addDForceGeneric`.",
"return_type": "void",
"parameters": [
{
"type": "const sofa::core::MechanicalParams*",
"name": "mparams"
},
{
"type": "DataVecDeriv&",
"name": "d_df"
},
{
"type": "const DataVecDeriv&",
"name": "d_dx"
}
],
"overridden": true,
"implementation_file": "./MultiThreading/component/solidmechanics/fem/elastic/ParallelTetrahedronFEMForceField.inl"
},
{
"name": "addKToMatrix",
"description": "Computes and adds the stiffness matrix to the global matrix. This method is overridden from the base class and optimized for parallel execution.",
"return_type": "void",
"parameters": [
{
"type": "sofa::linearalgebra::BaseMatrix*",
"name": "mat"
},
{
"type": "SReal",
"name": "kFactor"
},
{
"type": "unsigned int&",
"name": "offset"
}
],
"overridden": true,
"implementation_file": "./MultiThreading/component/solidmechanics/fem/elastic/ParallelTetrahedronFEMForceField.inl"
},
{
"name": "addDForceGeneric",
"description": "A generic function template used to calculate the derivative of forces with respect to displacements. It is specialized for different implementations (`addDForceSmall` and `addDForceCorotational`).",
"return_type": "void",
"parameters": [
{
"type": "VecDeriv&",
"name": "df"
},
{
"type": "const VecDeriv&",
"name": "dx"
},
{
"type": "Real",
"name": "kFactor"
},
{
"type": "const VecElement&",
"name": "indexedElements"
},
{
"type": "Function",
"name": "f"
}
],
"implementation_file": "./MultiThreading/component/solidmechanics/fem/elastic/ParallelTetrahedronFEMForceField.inl"
},
{
"name": "drawTrianglesFromTetrahedra",
"description": "Draws tetrahedral elements as triangles for visualization purposes. This method is overridden to optimize the drawing process in a parallel context.",
"return_type": "void",
"parameters": [
{
"type": "const sofa::core::visual::VisualParams*",
"name": "vparams"
},
{
"type": "bool",
"name": "showVonMisesStressPerElement"
},
{
"type": "bool",
"name": "drawVonMisesStress"
},
{
"type": "const VecCoord&",
"name": "x"
},
{
"type": "const VecReal&",
"name": "youngModulus"
},
{
"type": "bool",
"name": "heterogeneous"
},
{
"type": "Real",
"name": "minVM"
},
{
"type": "Real",
"name": "maxVM"
},
{
"type": "sofa::helper::ReadAccessor< sofa::Data<sofa::type::vector<Real>> >",
"name": "vM"
}
],
"overridden": true,
"implementation_file": "./MultiThreading/component/solidmechanics/fem/elastic/ParallelTetrahedronFEMForceField.inl"
}
],
"members": [
{
"name": "m_threadLocal_df",
"type": "std::map<std::thread::id, VecDeriv>",
"description": "Stores thread-local data for derivative calculations to avoid race conditions in parallel execution."
}
],
"registration_information": {
"registered_name": "ParallelTetrahedronFEMForceField",
"factory_registration_code": "registerParallelTetrahedronFEMForceField(sofa::core::ObjectFactory* factory)",
"description": "A parallel implementation of a linear elastic material using tetrahedral finite elements.",
"registration_file": "./MultiThreading/component/solidmechanics/fem/elastic/ParallelTetrahedronFEMForceField.cpp"
},
"notes": {
"optimization_context": "The ParallelTetrahedronFEMForceField is optimized for simulations with a large number of tetrahedral elements (> 1000), where parallel execution significantly improves performance.",
"parallel_implementation_registration": "Parallel implementations are registered in the `multithreading::ParallelImplementationsRegistry`, linking this class to its base implementation (`TetrahedronFEMForceField`).",
"template_instantiation": "The template is instantiated with `Vec3Types` for 3D simulations, making it suitable for a wide range of mechanical simulation scenarios."
}
},
"maths": "The **ParallelTetrahedronFEMForceField** component in the SOFA framework is designed to simulate linear elastic materials using tetrahedral finite elements. This implementation leverages parallel processing for efficient computation, particularly when dealing with a large number of tetrahedra (>1000). Below is a detailed mathematical and physical description of this component.\n\n### Governing Equations and Operators\n\nThis component implements the following key operators:\n\n- **Internal Force** ($f_{int}$)\n- **Stiffness Matrix** ($K$ or $\\mathbf{J}\\mathbf{K}\\mathbf{J}^T$, where $\\mathbf{J}$ is the Jacobian of the deformation gradient)\n- **Residual** ($R(x)$)\n\n### Constitutive and Kinematic Laws Involved\n\nThe component supports two formulations:\n1. **Small Deformation (Linear Elasticity)**\n2. **Corotational Formulation for Large Deformations**\n\n#### Strain Measures\n- For small deformations, the strain tensor $\\varepsilon$ is approximated linearly.\n- For large deformations, Green-Lagrange or Almansi strains are used in the corotational framework to account for rigid body rotations.\n\n#### Stress Tensors and Hyperelastic Potentials\nThe material behavior is assumed linear elastic with a Hookean constitutive law:\n\\[\n \\sigma = C : \\varepsilon,\n\\]\nwhere $C$ is the fourth-order elasticity tensor, and $\\varepsilon$ is the strain measure.\n\nFor small deformations:\n\\[\n K_{ij}^{tetrahedron} = -J^{-T}\\cdot D\\cdot J^{-1},\n\\]\nwhere $D$ is the tangent stiffness matrix computed from the constitutive law and strain-displacement relations, and $J$ is the Jacobian of the deformation gradient.\n\nFor corotational deformations:\n\\[\n K_{ij}^{tetrahedron}(R) = -J^{-T}\\cdot D(R)\\cdot J^{-1},\n\\]\nwhere $D(R)$ incorporates the rotation transformation matrix $R$ to decouple rigid body rotations from deformation.\n\n### Role in the Global FEM Pipeline\n\n#### Assembly Phase\n- **Mass and Force Construction**: The internal forces are assembled into the global force vector $f_{int}$, which is part of the semi-discrete dynamical system equation:\n\\[\n M \\ddot{x}(t) + f_{int}(x(t)) = f_{ext}(t).\n\\]\n- **Element Operators**: The stiffness matrix contributions for each tetrahedron are computed and assembled into the global stiffness matrix $K$.\n\n#### Time Integration Phase\nThe component supports implicit time integration schemes, typically Backward Euler or Newmark-beta methods. The nonlinear residual is constructed as:\n\\[\n R(x_{n+1}) = M \\frac{x_{n+1} - x_n}{\\Delta t} - \\Delta t f_{int}(x_{n+1}) - \\Delta t f_{ext}.\n\\]\n\n#### Nonlinear Resolution Phase\n- **Newton-Raphson Iteration**: The nonlinear residual is solved iteratively using Newton's method. Each iteration involves linearization:\n\\[\n R(x_k + \\delta x) \\approx R(x_k) + J(x_k) \\delta x,\n\\]\nwhere $J(x)$ is the Jacobian matrix given by\n\\[\n J(x) = \\frac{1}{\\Delta t} M - \\Delta t K(x).\n\\]\n\n### Numerical Methods and Discretization Choices\n- **Spatial Discretization**: The domain is discretized into tetrahedral elements. Each node carries 3 degrees of freedom (DOFs) in 3D space.\n- **Isoparametric Mapping**: Shape functions $N_i(X)$ are used to interpolate displacements within each element:\n\\[\n u(X, t) \\approx \\sum_i N_i(X) u_i(t).\n\\]\n- **Numerical Integration (Gauss Quadrature)**: Gauss quadrature is used for numerical integration over the reference tetrahedron.\n\n#### Parallelization Strategy\nThe component parallelizes the computation of internal forces and stiffness matrix contributions using a task scheduler. The following methods are executed in parallel:\n- `addDForce`\n- `addKToMatrix`\n\n### Variational / Lagrangian Mechanics Framework Fit\nThis component fits into the broader variational mechanics framework by discretizing the weak (variational) form of the governing PDEs using FEM. It ensures that the linearization and assembly steps preserve physical consistency, numerical stability, and extensibility within the SOFA simulation pipeline.",
"abstract": "The `ParallelTetrahedronFEMForceField` simulates linear elastic materials using tetrahedral finite elements with parallel processing, particularly efficient for simulations involving more than 1000 tetrahedra.",
"sheet": "# ParallelTetrahedronFEMForceField\n\n## Overview\nThe `ParallelTetrahedronFEMForceField` component simulates linear elastic materials using tetrahedral finite elements with parallel processing. It is particularly efficient for simulations involving more than 1000 tetrahedra, as it leverages parallel execution of internal force and stiffness matrix computations.\n\n## Mathematical Model\nThe component implements the following key operators:\n- **Internal Force** ($f_{int}$)\n- **Stiffness Matrix** ($K$ or $\boldsymbol{J}\boldsymbol{K}\boldsymbol{J}^T$, where $\boldsymbol{J}$ is the Jacobian of the deformation gradient)\n\nThe material behavior follows a Hookean constitutive law:\n\\[\n \\sigma = C : \\varepsilon,\n\\]\nwhere $C$ is the fourth-order elasticity tensor, and $\\varepsilon$ is the strain measure.\n\nFor small deformations:\n\\[\n K_{ij}^{tetrahedron} = -J^{-T}\\cdot D\\cdot J^{-1},\n\\]\nwhere $D$ is the tangent stiffness matrix computed from the constitutive law and strain-displacement relations, and $J$ is the Jacobian of the deformation gradient.\n\nFor corotational deformations:\n\\[\n K_{ij}^{tetrahedron}(R) = -J^{-T}\\cdot D(R)\\cdot J^{-1},\n\\]\nwhere $D(R)$ incorporates the rotation transformation matrix $R$ to decouple rigid body rotations from deformation.\n\n## Parameters and Data\nThe component does not expose any significant data fields. It inherits parameters and functionality from its parent classes `TaskSchedulerUser` and `TetrahedronFEMForceField`.\n\n## Dependencies and Connections\nThis component typically requires a mechanical state (degrees of freedom, positions, velocities) and exchanges data with solvers for force application and stiffness matrix assembly. It fits into the scene graph as part of the solid mechanics simulation pipeline."
}