ParallelHexahedronFEMForceField
Parallel implementation of a linear elastic material using hexahedral finite elements. Parallel implementation of HexahedronFEMForceField This implementation is the most efficient when: 1) the number of hexahedron is large (> 1000) 2) the global system matrix is not assembled. It is usually the case with a CGLinearSolver templated with GraphScattered types. 3) the method is 'large'. If the method is 'polar' or 'small', addForce is executed sequentially, but addDForce in parallel. The following methods are executed in parallel: - addForce for method 'large'. - addDForce The method addKToMatrix is not executed in parallel. This method is called with an assembled system, usually with a direct solver or a CGLinearSolver templated with types different from GraphScattered. In this case, the most time-consumming step is to invert the matrix. This is where efforts should be put to accelerate the simulation.
The `ParallelHexahedronFEMForceField` component implements a parallel version of linear elastic forces using hexahedral finite elements, optimized for large numbers (>1000) of hexahedra with un-assembled global system matrices.
- module
- MultiThreading
- namespace
- multithreading::component::forcefield::solidmechanics::fem::elastic
- include
- MultiThreading/component/solidmechanics/fem/elastic/ParallelHexahedronFEMForceField.h
- inherits
-
- HexahedronFEMForceField
- TaskSchedulerUser
- templates
-
- sofa::defaulttype::Vec3Types
The ParallelHexahedronFEMForceField component in the SOFA framework implements a parallel version of a linear elastic material using hexahedral finite elements. This component is specifically designed for efficient simulation when dealing with large numbers (>1000) of hexahedral elements and scenarios where the global system matrix is not assembled, typically with iterative solvers like CGLinearSolver templated with GraphScattered types.
Governing Equations and Operators
This component contributes to the internal force vector $\mathbf{f}_{int}$ through its implementation of linear elasticity. It computes the contribution from each hexahedral element to the global residual, which is part of the nonlinear solve process in FEM simulation:
egin{align}
R(x) &= M \frac{x_{n+1} - x_n}{\Delta t} + f_{int}(x_{n+1}) - f_{ext}(t)
\end{align}
Constitutive Laws and Kinematic Relations
The constitutive law implemented is linear elasticity, which can be described by the stress-strain relationship:
egin{equation}
\boldsymbol{\sigma} = \mathcal{C}:\boldsymbol{\varepsilon}
\end{equation}
where $\mathcal{C}$ is the fourth-order elasticity tensor and $\boldsymbol{\varepsilon}$ is the strain tensor.
For hexahedral elements, the displacement field within an element can be expressed as a polynomial interpolation of nodal displacements using shape functions $N_i$:
egin{equation}
u(X,t) \approx \sum_{i=1}^8 N_i(X) u_i(t)
\end{equation}
Numerical Methods and Discretization Choices
Element Assembly
The component performs element assembly in parallel, where each hexahedron contributes to the internal force vector $f_{int}$ through:
egin{align}
f^e_{int} &= \int_{\Omega_e} B_e^T \boldsymbol{\sigma}_e d\Omega
\end{align}
where $B_e$ is the strain-displacement matrix for element $e$. This step includes computing the rotation matrix to handle large deformations and applying the constitutive law.
Parallelization Strategy
- addForce: For method 'large', this function adds force contributions in parallel. It computes displacement, deformation, and stiffness matrices within each hexahedron before accumulating forces into nodal degrees of freedom (DOFs).
- addDForce: This function calculates the derivative of internal forces with respect to displacements, $\mathbf{K} = \frac{\partial f_{int}}{\partial x}$, in parallel by computing local stiffness matrices and applying them to nodal displacements.
Constraints and Mappings
The component does not directly handle constraints or mappings but integrates with the broader SOFA framework where such operations are performed. It ensures consistency within its domain of linear elasticity for hexahedral elements.
Role in Global FEM Pipeline
- Assembly: This component participates in the assembly phase by computing internal forces and stiffness contributions from hexahedra.
- Nonlinear Solve: The computed internal forces contribute to the residual equation, which is solved iteratively using methods like Newton-Raphson.
- Time Integration: While it does not directly implement time integration, its force computations are essential for the implicit integration schemes used in SOFA simulations.
Variational and Lagrangian Mechanics Framework
The component's role fits into a broader framework where the governing equations originate from variational principles. The weak form of the elasticity problem is discretized using FEM with hexahedral elements, ensuring consistency between analytical mechanics and numerical simulation.
Methods
void
init
()
void
addForce
(const sofa::core::MechanicalParams * mparams, DataVecDeriv & f, const DataVecCoord & x, const DataVecDeriv & v)
void
addDForce
(const sofa::core::MechanicalParams * mparams, DataVecDeriv & df, const DataVecDeriv & dx)
void
computeTaskForceLarge
(RDataRefVecCoord & p, int elementId, const Element & elem, const VecElementStiffness & elementStiffnesses, SReal & OutPotentialEnery, sofa::type::Vec<8, Deriv> & OutF)
{
"name": "ParallelHexahedronFEMForceField",
"namespace": "multithreading::component::forcefield::solidmechanics::fem::elastic",
"module": "MultiThreading",
"include": "MultiThreading/component/solidmechanics/fem/elastic/ParallelHexahedronFEMForceField.h",
"doc": "Parallel implementation of a linear elastic material using hexahedral finite elements.\n\nParallel implementation of HexahedronFEMForceField\nThis implementation is the most efficient when:\n1) the number of hexahedron is large (> 1000)\n2) the global system matrix is not assembled. It is usually the case with a CGLinearSolver templated with GraphScattered types.\n3) the method is 'large'. If the method is 'polar' or 'small', addForce is executed sequentially, but addDForce in parallel.\nThe following methods are executed in parallel:\n- addForce for method 'large'.\n- addDForce\nThe method addKToMatrix is not executed in parallel. This method is called with an assembled system, usually with\na direct solver or a CGLinearSolver templated with types different from GraphScattered. In this case, the most\ntime-consumming step is to invert the matrix. This is where efforts should be put to accelerate the simulation.",
"inherits": [
"HexahedronFEMForceField",
"TaskSchedulerUser"
],
"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": "f",
"type": "DataVecDeriv &"
},
{
"name": "x",
"type": "const DataVecCoord &"
},
{
"name": "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": "df",
"type": "DataVecDeriv &"
},
{
"name": "dx",
"type": "const DataVecDeriv &"
}
],
"is_virtual": false,
"is_pure_virtual": false,
"is_static": false,
"access": "public"
},
{
"name": "computeTaskForceLarge",
"return_type": "void",
"params": [
{
"name": "p",
"type": "RDataRefVecCoord &"
},
{
"name": "elementId",
"type": "int"
},
{
"name": "elem",
"type": "const Element &"
},
{
"name": "elementStiffnesses",
"type": "const VecElementStiffness &"
},
{
"name": "OutPotentialEnery",
"type": "SReal &"
},
{
"name": "OutF",
"type": "sofa::type::Vec<8, Deriv> &"
}
],
"is_virtual": false,
"is_pure_virtual": false,
"is_static": false,
"access": "protected"
}
],
"component": {
"name": "ParallelHexahedronFEMForceField",
"namespace": "multithreading::component::forcefield::solidmechanics::fem::elastic",
"description": "The `ParallelHexahedronFEMForceField` component is a parallel implementation of the Hexahedron Finite Element Method (FEM) force field in the SOFA (Simulation Open-Framework Architecture) framework. This implementation is designed to improve computational efficiency for simulations involving large numbers of hexahedral elements, typically greater than 1000. It leverages multithreading techniques to parallelize certain methods and processes, thereby reducing overall computation time.",
"inherited_from": [
"sofa::component::solidmechanics::fem::elastic::HexahedronFEMForceField<DataTypes>",
"sofa::simulation::TaskSchedulerUser"
],
"public_methods": [
{
"name": "addForce",
"description": "Adds force to the simulation. This method is executed in parallel when using the 'large' method, but sequentially for other methods like 'polar' or 'small'. It updates the force field based on mechanical parameters and current state of coordinates and velocities.",
"parameters": [
{
"name": "mparams",
"description": "Mechanical parameters used for the force computation."
},
{
"name": "f",
"description": "Data container for forces to be updated."
},
{
"name": "x",
"description": "Current positions/coordinates of elements in the simulation."
},
{
"name": "v",
"description": "Velocities of elements in the simulation."
}
]
},
{
"name": "addDForce",
"description": "Adds a derivative force to the simulation. This method is executed in parallel, consisting of two main steps: first, forces are computed inside each element; second, these forces are accumulated at vertices based on adjacent hexahedra.",
"parameters": [
{
"name": "mparams",
"description": "Mechanical parameters used for the force derivative computation."
},
{
"name": "df",
"description": "Data container for force derivatives to be updated."
},
{
"name": "dx",
"description": "Derivatives of positions/coordinates in the simulation."
}
]
}
],
"protected_methods": [
{
"name": "computeTaskForceLarge",
"description": "Helper method used to compute forces within each element, adapted to be thread-safe for parallel execution.",
"parameters": [
{
"name": "p",
"description": "Reference to positions/coordinates of elements."
},
{
"name": "elementId",
"description": "Identifier of the element being processed."
},
{
"name": "elem",
"description": "Hexahedral element object."
},
{
"name": "elementStiffnesses",
"description": "Vector containing stiffness matrices for each element."
},
{
"name": "outPotentialEnergy",
"description": "Output parameter to store potential energy."
},
{
"name": "outF",
"description": "Output vector of forces computed within the element."
}
]
}
],
"private_members": [
{
"name": "updateStiffnessMatrices",
"type": "bool",
"description": "Boolean flag used as a cache to avoid repeatedly accessing the 'd_updateStiffnessMatrix' value during execution."
},
{
"name": "m_vertexIdInAdjacentHexahedra",
"type": "sofa::type::vector<HexaAroundVerticesIndex>",
"description": "Vector storing indices of vertices within adjacent hexahedra, aiding in force accumulation at vertices."
},
{
"name": "m_elementsDf",
"type": "sofa::type::vector<sofa::type::Vec<8, Deriv> >",
"description": "Vector containing force derivatives for each element, stored as a class member to avoid reallocation during execution."
},
{
"name": "m_around",
"type": "sofa::type::vector<HexahedraAroundVertex>",
"description": "Cache storing lists of hexahedra around vertices, used in parallel processing of forces at vertices."
}
]
},
"maths": "The `ParallelHexahedronFEMForceField` component in the SOFA framework implements a parallel version of a linear elastic material using hexahedral finite elements. This component is specifically designed for efficient simulation when dealing with large numbers (>1000) of hexahedral elements and scenarios where the global system matrix is not assembled, typically with iterative solvers like `CGLinearSolver` templated with `GraphScattered` types.\n\n### Governing Equations and Operators\n\nThis component contributes to the internal force vector $\\mathbf{f}_{int}$ through its implementation of linear elasticity. It computes the contribution from each hexahedral element to the global residual, which is part of the nonlinear solve process in FEM simulation:\n\n\begin{align*}\nR(x) &= M \\frac{x_{n+1} - x_n}{\\Delta t} + f_{int}(x_{n+1}) - f_{ext}(t)\n\\end{align*}\n\n### Constitutive Laws and Kinematic Relations\n\nThe constitutive law implemented is linear elasticity, which can be described by the stress-strain relationship:\n\begin{equation*}\n\\boldsymbol{\\sigma} = \\mathcal{C}:\\boldsymbol{\\varepsilon}\n\\end{equation*}\nwhere $\\mathcal{C}$ is the fourth-order elasticity tensor and $\\boldsymbol{\\varepsilon}$ is the strain tensor.\n\nFor hexahedral elements, the displacement field within an element can be expressed as a polynomial interpolation of nodal displacements using shape functions $N_i$:\n\begin{equation*}\nu(X,t) \\approx \\sum_{i=1}^8 N_i(X) u_i(t)\n\\end{equation*}\n\n### Numerical Methods and Discretization Choices\n\n#### Element Assembly\n\nThe component performs element assembly in parallel, where each hexahedron contributes to the internal force vector $f_{int}$ through:\n\begin{align*}\nf^e_{int} &= \\int_{\\Omega_e} B_e^T \\boldsymbol{\\sigma}_e d\\Omega\n\\end{align*}\nwhere $B_e$ is the strain-displacement matrix for element $e$. This step includes computing the rotation matrix to handle large deformations and applying the constitutive law.\n\n#### Parallelization Strategy\n\n- **addForce**: For method 'large', this function adds force contributions in parallel. It computes displacement, deformation, and stiffness matrices within each hexahedron before accumulating forces into nodal degrees of freedom (DOFs).\n- **addDForce**: This function calculates the derivative of internal forces with respect to displacements, $\\mathbf{K} = \\frac{\\partial f_{int}}{\\partial x}$, in parallel by computing local stiffness matrices and applying them to nodal displacements.\n\n#### Constraints and Mappings\n\nThe component does not directly handle constraints or mappings but integrates with the broader SOFA framework where such operations are performed. It ensures consistency within its domain of linear elasticity for hexahedral elements.\n\n### Role in Global FEM Pipeline\n\n- **Assembly**: This component participates in the assembly phase by computing internal forces and stiffness contributions from hexahedra.\n- **Nonlinear Solve**: The computed internal forces contribute to the residual equation, which is solved iteratively using methods like Newton-Raphson.\n- **Time Integration**: While it does not directly implement time integration, its force computations are essential for the implicit integration schemes used in SOFA simulations.\n\n### Variational and Lagrangian Mechanics Framework\n\nThe component's role fits into a broader framework where the governing equations originate from variational principles. The weak form of the elasticity problem is discretized using FEM with hexahedral elements, ensuring consistency between analytical mechanics and numerical simulation.",
"abstract": "The `ParallelHexahedronFEMForceField` component implements a parallel version of linear elastic forces using hexahedral finite elements, optimized for large numbers (>1000) of hexahedra with un-assembled global system matrices.",
"sheet": "# ParallelHexahedronFEMForceField\n\n## Overview\nThe `ParallelHexahedronFEMForceField` component is part of the SOFA framework's multithreading module, designed for efficient parallel computation of linear elastic forces using hexahedral finite elements. It inherits from `HexahedronFEMForceField` and `TaskSchedulerUser`, indicating its role in handling large numbers (>1000) of hexahedra with un-assembled global system matrices.\n\n## Mathematical Model\nThe component implements linear elasticity, where the stress-strain relationship is given by:\n\\[ \\boldsymbol{\\sigma} = \\mathcal{C}:\\boldsymbol{\\varepsilon} \\]\nwhere $\\mathcal{C}$ is the fourth-order elasticity tensor and $\\boldsymbol{\\varepsilon}$ is the strain tensor. For hexahedral elements, the displacement field within an element can be expressed as:\n\\[ u(X,t) \\approx \\sum_{i=1}^8 N_i(X) u_i(t) \\]\nThe internal force contribution from each hexahedron to the global residual is computed through:\n\\begin{align*}\nf^e_{int} &= \\int_{\\Omega_e} B_e^T \\boldsymbol{\\sigma}_e d\\Omega\n\\end{align*}\nwhere $B_e$ is the strain-displacement matrix for element $e$. The component computes displacement, deformation, and stiffness matrices within each hexahedron before accumulating forces into nodal degrees of freedom (DOFs).\n\n## Parameters and Data\nThe `ParallelHexahedronFEMForceField` does not expose any significant data fields beyond those inherited from its parent classes.\n\n## Dependencies and Connections\nThis component typically requires a large number of hexahedral elements (>1000) and is optimized for scenarios where the global system matrix is not assembled, such as with iterative solvers like `CGLinearSolver` templated with `GraphScattered` types. It integrates with the broader SOFA framework's mechanical pipeline to compute internal forces and stiffness contributions.\n\n## Practical Notes\nThe component executes `addForce` in parallel for method 'large' and `addDForce` in parallel, while `addKToMatrix` remains sequential. The most time-consuming step is typically matrix inversion with direct solvers or non-GraphScattered types, where efforts should be focused to accelerate the simulation."
}