Back

TetrahedronFEMForceField

sofa::component::solidmechanics::fem::elastic::TetrahedronFEMForceField
BaseLinearElasticityFEMForceField
Doc (from source)

Tetrahedral finite elements. Compute Finite Element forces based on tetrahedral elements. Corotational methods are based on a rotation from world-space to material-space.

Abstract (AI generated)

TetrahedronFEMForceField computes finite element forces using tetrahedral elements for deformable object simulation in SOFA, supporting small and large displacement formulations with corotational methods.

Metadata
module
Sofa.Component.SolidMechanics.FEM.Elastic
namespace
sofa::component::solidmechanics::fem::elastic
include
sofa/component/solidmechanics/fem/elastic/TetrahedronFEMForceField.h
inherits
  • BaseLinearElasticityFEMForceField
templates
  • sofa::defaulttype::Vec3Types
description

Mathematical and Physical Description of TetrahedronFEMForceField Component

Overview

The TetrahedronFEMForceField is a computational model used in solid mechanics simulations, specifically within the Simulation Open-Framework Architecture (SOFA). It employs Finite Element Method (FEM) to compute forces on tetrahedral elements of deformable objects. This component supports both small and large displacements methods, as well as various corotational approaches for handling rotations between world space and material space.

Material Model

The underlying assumption is that the solid object behaves according to a linear elastic model. For an isotropic linear elastic material with Young's modulus $E$ and Poisson’s ratio $
u$, the material stiffness matrix $K$ in terms of Lamé parameters ($ ilde{oldsymbol{C}} = ext{diag}( ilde{oldsymbol{C}}{11}, ilde{oldsymbol{C}}{22}, ilde{oldsymbol{C}}{33}, ilde{oldsymbol{C}}{44})$) is given by:

$$ \tilde{\boldsymbol{C}} = \begin{pmatrix} \lambda + 2\mu & \lambda & \lambda & 0 \\ \lambda & \lambda + 2\mu & \lambda & 0 \\ \lambda & \lambda & \lambda + 2\mu & 0 \\ 0 & 0 & 0 & \mu \end{pmatrix} $$

where $\lambda = E \nu / ((1+\nu)(1-2\nu))$ and $\mu = E/(2(1+
u))$. Here, $\tilde{\boldsymbol{C}}$ is a diagonal matrix representing the material properties in Voigt notation.

Displacement Methods

Small Displacements (Linear)

The small displacements method assumes that deformations are small enough for linear approximations. The strain-displacement relationship is given by:

$$ \boldsymbol{\varepsilon} = \frac{1}{2}(\boldsymbol{\nabla u} + (\boldsymbol{\nabla u})^T) $$

where $\boldsymbol{u}$ is the displacement vector and $\boldsymbol{\varepsilon}$ represents the strain tensor.

Large Displacements (Corotational Methods)

The large displacements method can handle significant deformations by rotating between world space and material space. The rotation matrix $R$ aligns the deformed configuration with the initial one, and the strains are computed in this rotated frame:

$$ \boldsymbol{F} = R J $$

where $J$ is the deformation gradient in the corotational frame.

Strain-Displacement Matrix ($\mathbf{B}$)

The strain-displacement matrix $B$ links nodal displacements to strains within each tetrahedron. For an element with four vertices, the displacement vector $\boldsymbol{d} \in \mathbb{R}^{12}$ (with 3 components per vertex) and the strain vector $\boldsymbol{\varepsilon} \in \mathbb{R}^6$ are related by:

$$ \boldsymbol{\varepsilon} = B \cdot \boldsymbol{d} $$

where $B \in \mathbb{R}^{6 \times 12}$ is the strain-displacement matrix specific to each tetrahedral element.

Stiffness Matrix ($\mathbf{K}$)

The stiffness matrix $K$ for an individual element in corotational methods can be expressed as:

$$ K = R B^T \tilde{C} B R^T $$

where $R$ is the rotation matrix, $B$ is the strain-displacement matrix, and $\tilde{C}$ is the material stiffness matrix. For small displacements, the rotation matrix $R$ can be omitted.

Force Computation

The internal forces $\boldsymbol{f}_{int}$ are computed using the stiffness matrices and nodal displacements:

$$ \boldsymbol{f}_{int} = K (\boldsymbol{u}_d - \boldsymbol{u}_0) $$

where $\boldsymbol{u}_d$ is the current displacement vector, and $\boldsymbol{u}_0$ represents the initial positions.

Plasticity Model

The component supports plasticity models to simulate irreversible deformation. For instance, a creep factor $\eta$ can be defined as:

$$ \dot{p} = \eta (|\varepsilon_{\text{eff}}| - \sigma_y) $$

where $p$ is the accumulated plastic strain, $\eta$ is the creep rate, and $\sigma_y$ is the yield threshold.

Visualization (Von Mises Stress)

The component can compute Von Mises stress for visualization purposes. The von Mises equivalent stress $\sigma_{VM}$ is given by:

$$ \sigma_{VM} = \sqrt{\frac{3}{2}\boldsymbol{s}:\boldsymbol{s}} $$

where $s$ is the deviatoric stress tensor.

Summary

The TetrahedronFEMForceField component provides a versatile framework for simulating linear elastic materials using tetrahedral elements. It supports both small and large deformation regimes, with corotational methods to handle rotations accurately. This makes it suitable for applications ranging from biomedical simulations to real-time virtual material handling.

Data Fields
NameTypeDefaultHelp
data DataTypes
d_updateStiffnessMatrix bool
d_assembling bool
d_gatherPt sofa::helper::OptionsGroup number of dof accumulated per threads during the gather operation (Only use in GPU version)
d_gatherBsize sofa::helper::OptionsGroup number of dof accumulated per threads during the gather operation (Only use in GPU version)
d_drawHeterogeneousTetra bool Draw Heterogeneous Tetra in different color
d_computeVonMisesStress int compute and display von Mises stress: 0: no computations, 1: using corotational strain, 2: using full Green strain. Set listening=1
d_showStressAlpha float Alpha for vonMises visualisation
d_showVonMisesStressPerNode bool draw points showing vonMises stress interpolated in nodes
d_showVonMisesStressPerNodeColorMap bool draw elements showing vonMises stress interpolated in nodes
d_showVonMisesStressPerElement bool draw triangles showing vonMises stress interpolated in elements
d_updateStiffness bool update structures (precomputed in init) using stiffness parameters in each iteration (set listening=1)
Methods
int getRestVolume ()
void getRotation (Mat33 & R, int nodeIdx)
void getRotations (int & vecR)
void getRotations (linearalgebra::BaseMatrix * rotations, int offset)
const int & getRotations ()
void setComputeGlobalMatrix (bool val)
const Transformation & getActualTetraRotation (int index)
const Transformation & getInitialTetraRotation (int index)
const MaterialStiffness & getMaterialStiffness (TetrahedronID tetraId)
const StrainDisplacement & getStrainDisplacement (TetrahedronID tetraId)
const int & getRotatedInitialElements (TetrahedronID tetraId)
void setMethod (int methodName)
void setUpdateStiffnessMatrix (bool val)
void reset ()
void init ()
void reinit ()
void addForce (const core::MechanicalParams * mparams, int & d_f, const int & d_x, const int & d_v)
void addDForce (const core::MechanicalParams * mparams, int & d_df, const int & d_dx)
SReal getPotentialEnergy (const core::MechanicalParams * , const int & x)
void addKToMatrix (sofa::linearalgebra::BaseMatrix * m, SReal kFactor, 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)
void getElementStiffnessMatrix (int * stiffness, int nodeIdx)
void computeStrainDisplacement (StrainDisplacement & J, int a, int b, int c, int d)
int peudo_determinant_for_coef (const int & M)
void computeStiffnessMatrix (StiffnessMatrix & S, StiffnessMatrix & SR, const MaterialStiffness & K, const StrainDisplacement & J, const Transformation & Rot)
void computeMaterialStiffness (int i, int & a, int & b, int & c, int & d) virtual
void computeForce (Displacement & F, const Displacement & Depl, VoigtTensor & plasticStrain, const MaterialStiffness & K, const StrainDisplacement & J)
void computeForce (Displacement & F, const Displacement & Depl, const MaterialStiffness & K, const StrainDisplacement & J, SReal fact)
void initSmall (int i, int & a, int & b, int & c, int & d)
void accumulateForceSmall (int & f, const int & p, int elementIt, int elementIndex)
void applyStiffnessSmall (int & f, const int & x, int i, int a, int b, int c, int d, SReal fact)
void initLarge (int i, int & a, int & b, int & c, int & d)
void computeRotationLarge (Transformation & r, const int & p, const int & a, const int & b, const int & c)
void accumulateForceLarge (int & f, const int & p, int elementIt, int elementIndex)
void initPolar (int i, int & a, int & b, int & c, int & d)
void accumulateForcePolar (int & f, const int & p, int elementIt, int elementIndex)
void initSVD (int i, int & a, int & b, int & c, int & d)
void accumulateForceSVD (int & f, const int & p, int elementIt, int elementIndex)
void applyStiffnessCorotational (int & f, const int & x, int i, int a, int b, int c, int d, SReal fact)
void handleTopologyChange ()
void computeVonMisesStress ()
bool isComputeVonMisesStressMethodSet ()
void computeMinMaxFromYoungsModulus ()
void drawTrianglesFromTetrahedra (const core::visual::VisualParams * vparams, bool showVonMisesStressPerElement, bool drawVonMisesStress, const int & x, const int & youngModulus, bool heterogeneous, int minVM, int maxVM, int vM) virtual
void drawTrianglesFromRangeOfTetrahedra (const int & range, const core::visual::VisualParams * vparams, bool showVonMisesStressPerElement, bool drawVonMisesStress, bool showWireFrame, const int & x, const int & youngModulus, bool heterogeneous, int minVM, int maxVM, int vM) virtual
void handleEvent (core::objectmodel::Event * event) 
{
  "name": "TetrahedronFEMForceField",
  "namespace": "sofa::component::solidmechanics::fem::elastic",
  "module": "Sofa.Component.SolidMechanics.FEM.Elastic",
  "include": "sofa/component/solidmechanics/fem/elastic/TetrahedronFEMForceField.h",
  "doc": "Tetrahedral finite elements.\n\nCompute Finite Element forces based on tetrahedral elements.\n  Corotational methods are based on a rotation from world-space to material-space.",
  "inherits": [
    "BaseLinearElasticityFEMForceField"
  ],
  "templates": [
    "sofa::defaulttype::Vec3Types"
  ],
  "data_fields": [
    {
      "name": "data",
      "type": "DataTypes"
    },
    {
      "name": "d_updateStiffnessMatrix",
      "type": "bool",
      "xmlname": "updateStiffnessMatrix",
      "help": ""
    },
    {
      "name": "d_assembling",
      "type": "bool",
      "xmlname": "computeGlobalMatrix",
      "help": ""
    },
    {
      "name": "d_gatherPt",
      "type": "sofa::helper::OptionsGroup",
      "xmlname": "gatherPt",
      "help": "number of dof accumulated per threads during the gather operation (Only use in GPU version)"
    },
    {
      "name": "d_gatherBsize",
      "type": "sofa::helper::OptionsGroup",
      "xmlname": "gatherBsize",
      "help": "number of dof accumulated per threads during the gather operation (Only use in GPU version)"
    },
    {
      "name": "d_drawHeterogeneousTetra",
      "type": "bool",
      "xmlname": "drawHeterogeneousTetra",
      "help": "Draw Heterogeneous Tetra in different color"
    },
    {
      "name": "d_computeVonMisesStress",
      "type": "int",
      "xmlname": "computeVonMisesStress",
      "help": "compute and display von Mises stress: 0: no computations, 1: using corotational strain, 2: using full Green strain. Set listening=1"
    },
    {
      "name": "d_showStressAlpha",
      "type": "float",
      "xmlname": "showStressAlpha",
      "help": "Alpha for vonMises visualisation"
    },
    {
      "name": "d_showVonMisesStressPerNode",
      "type": "bool",
      "xmlname": "showVonMisesStressPerNode",
      "help": "draw points showing vonMises stress interpolated in nodes"
    },
    {
      "name": "d_showVonMisesStressPerNodeColorMap",
      "type": "bool",
      "xmlname": "showVonMisesStressPerNodeColorMap",
      "help": "draw elements showing vonMises stress interpolated in nodes"
    },
    {
      "name": "d_showVonMisesStressPerElement",
      "type": "bool",
      "xmlname": "showVonMisesStressPerElement",
      "help": "draw triangles showing vonMises stress interpolated in elements"
    },
    {
      "name": "d_updateStiffness",
      "type": "bool",
      "xmlname": "updateStiffness",
      "help": "update structures (precomputed in init) using stiffness parameters in each iteration (set listening=1)"
    }
  ],
  "links": [],
  "methods": [
    {
      "name": "getRestVolume",
      "return_type": "int",
      "params": [],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "public"
    },
    {
      "name": "getRotation",
      "return_type": "void",
      "params": [
        {
          "name": "R",
          "type": "Mat33 &"
        },
        {
          "name": "nodeIdx",
          "type": "int"
        }
      ],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "public"
    },
    {
      "name": "getRotations",
      "return_type": "void",
      "params": [
        {
          "name": "vecR",
          "type": "int &"
        }
      ],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "public"
    },
    {
      "name": "getRotations",
      "return_type": "void",
      "params": [
        {
          "name": "rotations",
          "type": "linearalgebra::BaseMatrix *"
        },
        {
          "name": "offset",
          "type": "int"
        }
      ],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "public"
    },
    {
      "name": "getRotations",
      "return_type": "const int &",
      "params": [],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "public"
    },
    {
      "name": "setComputeGlobalMatrix",
      "return_type": "void",
      "params": [
        {
          "name": "val",
          "type": "bool"
        }
      ],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "public"
    },
    {
      "name": "getActualTetraRotation",
      "return_type": "const Transformation &",
      "params": [
        {
          "name": "index",
          "type": "int"
        }
      ],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "public"
    },
    {
      "name": "getInitialTetraRotation",
      "return_type": "const Transformation &",
      "params": [
        {
          "name": "index",
          "type": "int"
        }
      ],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "public"
    },
    {
      "name": "getMaterialStiffness",
      "return_type": "const MaterialStiffness &",
      "params": [
        {
          "name": "tetraId",
          "type": "TetrahedronID"
        }
      ],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "public"
    },
    {
      "name": "getStrainDisplacement",
      "return_type": "const StrainDisplacement &",
      "params": [
        {
          "name": "tetraId",
          "type": "TetrahedronID"
        }
      ],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "public"
    },
    {
      "name": "getRotatedInitialElements",
      "return_type": "const int &",
      "params": [
        {
          "name": "tetraId",
          "type": "TetrahedronID"
        }
      ],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "public"
    },
    {
      "name": "setMethod",
      "return_type": "void",
      "params": [
        {
          "name": "methodName",
          "type": "int"
        }
      ],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "public"
    },
    {
      "name": "setUpdateStiffnessMatrix",
      "return_type": "void",
      "params": [
        {
          "name": "val",
          "type": "bool"
        }
      ],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "public"
    },
    {
      "name": "reset",
      "return_type": "void",
      "params": [],
      "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": "reinit",
      "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": "int &"
        },
        {
          "name": "d_x",
          "type": "const int &"
        },
        {
          "name": "d_v",
          "type": "const int &"
        }
      ],
      "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": "int &"
        },
        {
          "name": "d_dx",
          "type": "const int &"
        }
      ],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "public"
    },
    {
      "name": "getPotentialEnergy",
      "return_type": "SReal",
      "params": [
        {
          "name": "",
          "type": "const core::MechanicalParams *"
        },
        {
          "name": "x",
          "type": "const int &"
        }
      ],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "public"
    },
    {
      "name": "addKToMatrix",
      "return_type": "void",
      "params": [
        {
          "name": "m",
          "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": "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": "getElementStiffnessMatrix",
      "return_type": "void",
      "params": [
        {
          "name": "stiffness",
          "type": "int *"
        },
        {
          "name": "nodeIdx",
          "type": "int"
        }
      ],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "public"
    },
    {
      "name": "computeStrainDisplacement",
      "return_type": "void",
      "params": [
        {
          "name": "J",
          "type": "StrainDisplacement &"
        },
        {
          "name": "a",
          "type": "int"
        },
        {
          "name": "b",
          "type": "int"
        },
        {
          "name": "c",
          "type": "int"
        },
        {
          "name": "d",
          "type": "int"
        }
      ],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "protected"
    },
    {
      "name": "peudo_determinant_for_coef",
      "return_type": "int",
      "params": [
        {
          "name": "M",
          "type": "const int &"
        }
      ],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "protected"
    },
    {
      "name": "computeStiffnessMatrix",
      "return_type": "void",
      "params": [
        {
          "name": "S",
          "type": "StiffnessMatrix &"
        },
        {
          "name": "SR",
          "type": "StiffnessMatrix &"
        },
        {
          "name": "K",
          "type": "const MaterialStiffness &"
        },
        {
          "name": "J",
          "type": "const StrainDisplacement &"
        },
        {
          "name": "Rot",
          "type": "const Transformation &"
        }
      ],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "protected"
    },
    {
      "name": "computeMaterialStiffness",
      "return_type": "void",
      "params": [
        {
          "name": "i",
          "type": "int"
        },
        {
          "name": "a",
          "type": "int &"
        },
        {
          "name": "b",
          "type": "int &"
        },
        {
          "name": "c",
          "type": "int &"
        },
        {
          "name": "d",
          "type": "int &"
        }
      ],
      "is_virtual": true,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "protected"
    },
    {
      "name": "computeForce",
      "return_type": "void",
      "params": [
        {
          "name": "F",
          "type": "Displacement &"
        },
        {
          "name": "Depl",
          "type": "const Displacement &"
        },
        {
          "name": "plasticStrain",
          "type": "VoigtTensor &"
        },
        {
          "name": "K",
          "type": "const MaterialStiffness &"
        },
        {
          "name": "J",
          "type": "const StrainDisplacement &"
        }
      ],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "protected"
    },
    {
      "name": "computeForce",
      "return_type": "void",
      "params": [
        {
          "name": "F",
          "type": "Displacement &"
        },
        {
          "name": "Depl",
          "type": "const Displacement &"
        },
        {
          "name": "K",
          "type": "const MaterialStiffness &"
        },
        {
          "name": "J",
          "type": "const StrainDisplacement &"
        },
        {
          "name": "fact",
          "type": "SReal"
        }
      ],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "protected"
    },
    {
      "name": "initSmall",
      "return_type": "void",
      "params": [
        {
          "name": "i",
          "type": "int"
        },
        {
          "name": "a",
          "type": "int &"
        },
        {
          "name": "b",
          "type": "int &"
        },
        {
          "name": "c",
          "type": "int &"
        },
        {
          "name": "d",
          "type": "int &"
        }
      ],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "protected"
    },
    {
      "name": "accumulateForceSmall",
      "return_type": "void",
      "params": [
        {
          "name": "f",
          "type": "int &"
        },
        {
          "name": "p",
          "type": "const int &"
        },
        {
          "name": "elementIt",
          "type": "int"
        },
        {
          "name": "elementIndex",
          "type": "int"
        }
      ],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "protected"
    },
    {
      "name": "applyStiffnessSmall",
      "return_type": "void",
      "params": [
        {
          "name": "f",
          "type": "int &"
        },
        {
          "name": "x",
          "type": "const int &"
        },
        {
          "name": "i",
          "type": "int"
        },
        {
          "name": "a",
          "type": "int"
        },
        {
          "name": "b",
          "type": "int"
        },
        {
          "name": "c",
          "type": "int"
        },
        {
          "name": "d",
          "type": "int"
        },
        {
          "name": "fact",
          "type": "SReal"
        }
      ],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "protected"
    },
    {
      "name": "initLarge",
      "return_type": "void",
      "params": [
        {
          "name": "i",
          "type": "int"
        },
        {
          "name": "a",
          "type": "int &"
        },
        {
          "name": "b",
          "type": "int &"
        },
        {
          "name": "c",
          "type": "int &"
        },
        {
          "name": "d",
          "type": "int &"
        }
      ],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "protected"
    },
    {
      "name": "computeRotationLarge",
      "return_type": "void",
      "params": [
        {
          "name": "r",
          "type": "Transformation &"
        },
        {
          "name": "p",
          "type": "const int &"
        },
        {
          "name": "a",
          "type": "const int &"
        },
        {
          "name": "b",
          "type": "const int &"
        },
        {
          "name": "c",
          "type": "const int &"
        }
      ],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "protected"
    },
    {
      "name": "accumulateForceLarge",
      "return_type": "void",
      "params": [
        {
          "name": "f",
          "type": "int &"
        },
        {
          "name": "p",
          "type": "const int &"
        },
        {
          "name": "elementIt",
          "type": "int"
        },
        {
          "name": "elementIndex",
          "type": "int"
        }
      ],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "protected"
    },
    {
      "name": "initPolar",
      "return_type": "void",
      "params": [
        {
          "name": "i",
          "type": "int"
        },
        {
          "name": "a",
          "type": "int &"
        },
        {
          "name": "b",
          "type": "int &"
        },
        {
          "name": "c",
          "type": "int &"
        },
        {
          "name": "d",
          "type": "int &"
        }
      ],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "protected"
    },
    {
      "name": "accumulateForcePolar",
      "return_type": "void",
      "params": [
        {
          "name": "f",
          "type": "int &"
        },
        {
          "name": "p",
          "type": "const int &"
        },
        {
          "name": "elementIt",
          "type": "int"
        },
        {
          "name": "elementIndex",
          "type": "int"
        }
      ],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "protected"
    },
    {
      "name": "initSVD",
      "return_type": "void",
      "params": [
        {
          "name": "i",
          "type": "int"
        },
        {
          "name": "a",
          "type": "int &"
        },
        {
          "name": "b",
          "type": "int &"
        },
        {
          "name": "c",
          "type": "int &"
        },
        {
          "name": "d",
          "type": "int &"
        }
      ],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "protected"
    },
    {
      "name": "accumulateForceSVD",
      "return_type": "void",
      "params": [
        {
          "name": "f",
          "type": "int &"
        },
        {
          "name": "p",
          "type": "const int &"
        },
        {
          "name": "elementIt",
          "type": "int"
        },
        {
          "name": "elementIndex",
          "type": "int"
        }
      ],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "protected"
    },
    {
      "name": "applyStiffnessCorotational",
      "return_type": "void",
      "params": [
        {
          "name": "f",
          "type": "int &"
        },
        {
          "name": "x",
          "type": "const int &"
        },
        {
          "name": "i",
          "type": "int"
        },
        {
          "name": "a",
          "type": "int"
        },
        {
          "name": "b",
          "type": "int"
        },
        {
          "name": "c",
          "type": "int"
        },
        {
          "name": "d",
          "type": "int"
        },
        {
          "name": "fact",
          "type": "SReal"
        }
      ],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "protected"
    },
    {
      "name": "handleTopologyChange",
      "return_type": "void",
      "params": [],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "protected"
    },
    {
      "name": "computeVonMisesStress",
      "return_type": "void",
      "params": [],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "protected"
    },
    {
      "name": "isComputeVonMisesStressMethodSet",
      "return_type": "bool",
      "params": [],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "protected"
    },
    {
      "name": "computeMinMaxFromYoungsModulus",
      "return_type": "void",
      "params": [],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "protected"
    },
    {
      "name": "drawTrianglesFromTetrahedra",
      "return_type": "void",
      "params": [
        {
          "name": "vparams",
          "type": "const core::visual::VisualParams *"
        },
        {
          "name": "showVonMisesStressPerElement",
          "type": "bool"
        },
        {
          "name": "drawVonMisesStress",
          "type": "bool"
        },
        {
          "name": "x",
          "type": "const int &"
        },
        {
          "name": "youngModulus",
          "type": "const int &"
        },
        {
          "name": "heterogeneous",
          "type": "bool"
        },
        {
          "name": "minVM",
          "type": "int"
        },
        {
          "name": "maxVM",
          "type": "int"
        },
        {
          "name": "vM",
          "type": "int"
        }
      ],
      "is_virtual": true,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "protected"
    },
    {
      "name": "drawTrianglesFromRangeOfTetrahedra",
      "return_type": "void",
      "params": [
        {
          "name": "range",
          "type": "const int &"
        },
        {
          "name": "vparams",
          "type": "const core::visual::VisualParams *"
        },
        {
          "name": "showVonMisesStressPerElement",
          "type": "bool"
        },
        {
          "name": "drawVonMisesStress",
          "type": "bool"
        },
        {
          "name": "showWireFrame",
          "type": "bool"
        },
        {
          "name": "x",
          "type": "const int &"
        },
        {
          "name": "youngModulus",
          "type": "const int &"
        },
        {
          "name": "heterogeneous",
          "type": "bool"
        },
        {
          "name": "minVM",
          "type": "int"
        },
        {
          "name": "maxVM",
          "type": "int"
        },
        {
          "name": "vM",
          "type": "int"
        }
      ],
      "is_virtual": true,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "protected"
    },
    {
      "name": "handleEvent",
      "return_type": "void",
      "params": [
        {
          "name": "event",
          "type": "core::objectmodel::Event *"
        }
      ],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "protected"
    }
  ],
  "title": "Tetrahedron FEM Force Field",
  "category": [
    "Solid Mechanics"
  ],
  "briefDescription": "Computes finite element forces based on tetrahedral elements for solid mechanics simulations.",
  "description": "This component calculates finite element forces using a tetrahedral mesh for solid object simulation. It implements both small and large displacement formulations, including corotational methods that use rotations to transform between world space and material space.\n\nKey Features:\n- Supports 'small', 'large', 'polar' and 'svd' displacement methods\n- Can compute Von Mises stress for visualization\n- Allows setting local stiffness factors per element \n- Provides accessors for retrieving rotation matrices, stiffness matrices etc. of individual tetrahedra\n- Implements full system matrix assembly support\n- Supports plasticity models like those in \"Interactive Virtual Materials\" paper\n\nThe component computes forces based on a specified material model (linear elasticity by default) and applies them to the deformable object. It can be used for simulating soft tissues, virtual materials etc. in real-time or offline physics simulations.",
  "parameters": [
    {
      "name": "method",
      "type": "string",
      "defaultValue": "'small'",
      "description": "Displacement method: 'small', 'large' (QR), 'polar' or 'svd'"
    },
    {
      "name": "localStiffnessFactor",
      "type": "vector<Real>",
      "defaultValue": "",
      "description": "Allows specifying different stiffness per element"
    },
    {
      "name": "youngModulus",
      "type": "Real or vector<Real>",
      "defaultValue": "1000.0",
      "description": "Young's modulus of the material (or per element)"
    },
    {
      "name": "poissonRatio",
      "type": "Real or vector<Real>",
      "defaultValue": "0.3",
      "description": "Poisson ratio of the material (or per element)"
    },
    {
      "name": "updateStiffnessMatrix",
      "type": "bool",
      "defaultValue": "false",
      "description": "Update stiffness matrices each iteration"
    },
    {
      "name": "assembling",
      "type": "bool",
      "defaultValue": "false",
      "description": "Enable full system matrix assembly support"
    },
    {
      "name": "plasticMaxThreshold",
      "type": "Real",
      "defaultValue": "0.01",
      "description": "Plasticity max threshold (2-norm of strain)"
    },
    {
      "name": "plasticYieldThreshold",
      "type": "Real",
      "defaultValue": "0.03",
      "description": "Plastic yield threshold"
    },
    {
      "name": "plasticCreep",
      "type": "Real",
      "defaultValue": "1e-4",
      "description": "Plastic creep factor [0,1]"
    }
  ],
  "inputs": [
    {
      "name": "position",
      "description": "Current positions of the object's points",
      "type": "vector<Vec3>"
    },
    {
      "name": "velocity",
      "description": "Current velocities of the object's points",
      "type": "vector<Vec3>"
    }
  ],
  "outputs": [
    {
      "name": "force",
      "description": "Calculated forces applied to each point",
      "type": "vector<Vec3>"
    }
  ],
  "maths": "## Mathematical and Physical Description of TetrahedronFEMForceField Component\n\n### Overview\nThe **TetrahedronFEMForceField** is a computational model used in solid mechanics simulations, specifically within the Simulation Open-Framework Architecture (SOFA). It employs Finite Element Method (FEM) to compute forces on tetrahedral elements of deformable objects. This component supports both small and large displacements methods, as well as various corotational approaches for handling rotations between world space and material space.\n\n### Material Model\nThe underlying assumption is that the solid object behaves according to a linear elastic model. For an isotropic linear elastic material with Young's modulus $E$ and Poisson’s ratio $\nu$, the material stiffness matrix $K$ in terms of Lamé parameters ($\tilde{\boldsymbol{C}} = \text{diag}(\tilde{\boldsymbol{C}}_{11}, \tilde{\boldsymbol{C}}_{22}, \tilde{\boldsymbol{C}}_{33}, \tilde{\boldsymbol{C}}_{44})$) is given by:\n\n\\[\n    \\tilde{\\boldsymbol{C}} =\n        \\begin{pmatrix}\n            \\lambda + 2\\mu &   \\lambda         &   \\lambda          &   0                \\\\\n            \\lambda         & \\lambda + 2\\mu   &   \\lambda          &   0                \\\\\n            \\lambda         &   \\lambda         & \\lambda + 2\\mu    &   0                \\\\\n            0               &   0               &   0                &     \\mu\n        \\end{pmatrix}\n\\]\nwhere $\\lambda = E \\nu / ((1+\\nu)(1-2\\nu))$ and $\\mu = E/(2(1+\nu))$. Here, $\\tilde{\\boldsymbol{C}}$ is a diagonal matrix representing the material properties in Voigt notation.\n\n### Displacement Methods\n#### Small Displacements (Linear)\nThe small displacements method assumes that deformations are small enough for linear approximations. The strain-displacement relationship is given by:\n\n\\[\n    \\boldsymbol{\\varepsilon} = \\frac{1}{2}(\\boldsymbol{\\nabla u} + (\\boldsymbol{\\nabla u})^T)\n\\]\nwhere $\\boldsymbol{u}$ is the displacement vector and $\\boldsymbol{\\varepsilon}$ represents the strain tensor.\n\n#### Large Displacements (Corotational Methods)\nThe large displacements method can handle significant deformations by rotating between world space and material space. The rotation matrix $R$ aligns the deformed configuration with the initial one, and the strains are computed in this rotated frame:\n\n\\[\n    \\boldsymbol{F} = R J \n\\]\nwhere $J$ is the deformation gradient in the corotational frame.\n\n### Strain-Displacement Matrix ($\\mathbf{B}$)\nThe strain-displacement matrix $B$ links nodal displacements to strains within each tetrahedron. For an element with four vertices, the displacement vector $\\boldsymbol{d} \\in \\mathbb{R}^{12}$ (with 3 components per vertex) and the strain vector $\\boldsymbol{\\varepsilon} \\in \\mathbb{R}^6$ are related by:\n\n\\[\n    \\boldsymbol{\\varepsilon} = B \\cdot \\boldsymbol{d}\n\\]\nwhere $B \\in \\mathbb{R}^{6 \\times 12}$ is the strain-displacement matrix specific to each tetrahedral element.\n\n### Stiffness Matrix ($\\mathbf{K}$)\nThe stiffness matrix $K$ for an individual element in corotational methods can be expressed as:\n\n\\[\n    K = R B^T \\tilde{C} B R^T \n\\]\nwhere $R$ is the rotation matrix, $B$ is the strain-displacement matrix, and $\\tilde{C}$ is the material stiffness matrix. For small displacements, the rotation matrix $R$ can be omitted.\n\n### Force Computation\nThe internal forces $\\boldsymbol{f}_{int}$ are computed using the stiffness matrices and nodal displacements:\n\n\\[\n    \\boldsymbol{f}_{int} = K (\\boldsymbol{u}_d - \\boldsymbol{u}_0)\n\\]\nwhere $\\boldsymbol{u}_d$ is the current displacement vector, and $\\boldsymbol{u}_0$ represents the initial positions.\n\n### Plasticity Model\nThe component supports plasticity models to simulate irreversible deformation. For instance, a creep factor $\\eta$ can be defined as:\n\n\\[\n    \\dot{p} = \\eta (|\\varepsilon_{\\text{eff}}| - \\sigma_y)\n\\]\nwhere $p$ is the accumulated plastic strain, $\\eta$ is the creep rate, and $\\sigma_y$ is the yield threshold.\n\n### Visualization (Von Mises Stress)\nThe component can compute Von Mises stress for visualization purposes. The von Mises equivalent stress $\\sigma_{VM}$ is given by:\n\n\\[\n    \\sigma_{VM} = \\sqrt{\\frac{3}{2}\\boldsymbol{s}:\\boldsymbol{s}}\n\\]\nwhere $s$ is the deviatoric stress tensor.\n\n### Summary\nThe **TetrahedronFEMForceField** component provides a versatile framework for simulating linear elastic materials using tetrahedral elements. It supports both small and large deformation regimes, with corotational methods to handle rotations accurately. This makes it suitable for applications ranging from biomedical simulations to real-time virtual material handling.",
  "abstract": "TetrahedronFEMForceField computes finite element forces using tetrahedral elements for deformable object simulation in SOFA, supporting small and large displacement formulations with corotational methods.",
  "sheet": "# TetrahedronFEMForceField\n\n## Overview\nThe **TetrahedronFEMForceField** component calculates finite element forces on tetrahedral elements of deformable objects. It supports both small and large displacement formulations, including corotational methods that handle rotations between world space and material space.\n\n## Mathematical Model\n### Material Model\nFor an isotropic linear elastic material with Young's modulus $E$ and Poisson’s ratio $\nu$, the material stiffness matrix $K$ in terms of Lamé parameters is given by:\n\n\\[\n    \\tilde{\\boldsymbol{C}} =\n        \\begin{pmatrix}\n            \\lambda + 2\\mu &   \\lambda         &   \\lambda          &   0                \\\\\n            \\lambda         & \\lambda + 2\\mu   &   \\lambda          &   0                \\\\\n            \\lambda         &   \\lambda         & \\lambda + 2\\mu    &   0                \\\\\n            0               &   0               &   0                &     \\mu\n        \\end{pmatrix}\n\\]\nwhere $\\lambda = E \\nu / ((1+\\nu)(1-2\\nu))$ and $\\mu = E/(2(1+\\nu))$. Here, $\\tilde{\\boldsymbol{C}}$ is a diagonal matrix representing the material properties in Voigt notation.\n\n### Displacement Methods\n#### Small Displacements (Linear)\nThe small displacements method assumes that deformations are small enough for linear approximations. The strain-displacement relationship is given by:\n\n\\[\n    \\boldsymbol{\\varepsilon} = \\frac{1}{2}(\\boldsymbol{\\nabla u} + (\\boldsymbol{\\nabla u})^T)\n\\]\nwhere $\\boldsymbol{u}$ is the displacement vector and $\\boldsymbol{\\varepsilon}$ represents the strain tensor.\n\n#### Large Displacements (Corotational Methods)\nThe large displacements method can handle significant deformations by rotating between world space and material space. The rotation matrix $R$ aligns the deformed configuration with the initial one, and the strains are computed in this rotated frame:\n\n\\[\n    \\boldsymbol{F} = R J \n\\]\nwhere $J$ is the deformation gradient in the corotational frame.\n\n### Strain-Displacement Matrix ($B$)\nThe strain-displacement matrix $B$ links nodal displacements to strains within each tetrahedron. For an element with four vertices, the displacement vector $d \\in \\mathbb{R}^{12}$ (with 3 components per vertex) and the strain vector $\\boldsymbol{\\varepsilon} \\in \\mathbb{R}^6$ are related by:\n\n\\[\n    \\boldsymbol{\\varepsilon} = B \\cdot d\n\\]\nwhere $B \\in \\mathbb{R}^{6 \\times 12}$ is the strain-displacement matrix specific to each tetrahedral element.\n\n### Stiffness Matrix ($K$)\nThe stiffness matrix $K$ for an individual element in corotational methods can be expressed as:\n\n\\[\n    K = R B^T \\tilde{C} B R^T \n\\]\nwhere $R$ is the rotation matrix, $B$ is the strain-displacement matrix, and $\\tilde{C}$ is the material stiffness matrix. For small displacements, the rotation matrix $R$ can be omitted.\n\n### Force Computation\nThe internal forces $f_{int}$ are computed using the stiffness matrices and nodal displacements:\n\n\\[\n    f_{int} = K (u_d - u_0)\n\\]\nwhere $u_d$ is the current displacement vector, and $u_0$ represents the initial positions.\n\n### Visualization (Von Mises Stress)\nThe component can compute Von Mises stress for visualization purposes. The von Mises equivalent stress $\\sigma_{VM}$ is given by:\n\n\\[\n    \\sigma_{VM} = \\sqrt{\\frac{3}{2}s:s}\n\\]\nwhere $s$ is the deviatoric stress tensor.\n\n## Parameters and Data\n- **updateStiffnessMatrix**: Boolean to update stiffness matrix (default: false).\n- **computeGlobalMatrix**: Boolean to compute global matrix (default: false).\n- **gatherPt**: Number of degrees of freedom accumulated per thread during gather operation (GPU version only).\n- **gatherBsize**: Number of degrees of freedom accumulated per thread during gather operation (GPU version only).\n- **drawHeterogeneousTetra**: Boolean to draw heterogeneous tetrahedra in different colors.\n- **computeVonMisesStress**: Integer to compute and display von Mises stress: 0 - no computations, 1 - using corotational strain, 2 - using full Green strain (default: 0).\n- **showStressAlpha**: Alpha for von Mises visualisation.\n- **showVonMisesStressPerNode**: Boolean to draw points showing von Mises stress interpolated in nodes.\n- **showVonMisesStressPerNodeColorMap**: Boolean to draw elements showing von Mises stress interpolated in nodes.\n- **showVonMisesStressPerElement**: Boolean to draw triangles showing von Mises stress interpolated in elements.\n- **updateStiffness**: Boolean to update structures using stiffness parameters (default: false)."
}