Data Fields

Methods

{
  "name": "SarcomereContractionForceField",
  "namespace": "sofa::component",
  "module": "SofaTestPlugin",
  "include": "forcefield/sarcomere/SarcomereContractionForceField.h",
  "doc": "Applies functional contraction (stress is computed from a function of strain) force along fiber directions in tetrahedral meshes",
  "inherits": [
    "BaseTetraForceField"
  ],
  "templates": [
    "sofa::defaulttype::Vec3Types"
  ],
  "data_fields": [
    {
      "name": "d_depol",
      "type": "VecReal",
      "xmlname": "depol",
      "help": "Depolarisation time at each node (from loader)"
    },
    {
      "name": "d_apd",
      "type": "VecReal",
      "xmlname": "apd",
      "help": "Action potential duration at each node (from loader)"
    },
    {
      "name": "d_sigma0",
      "type": "Real",
      "xmlname": "sigma0",
      "help": "Parameter sigma0"
    },
    {
      "name": "d_k0",
      "type": "Real",
      "xmlname": "k0",
      "help": "Parameter k0"
    },
    {
      "name": "d_hp",
      "type": "Real",
      "xmlname": "heart_period",
      "help": ""
    },
    {
      "name": "d_katp",
      "type": "Real",
      "xmlname": "katp",
      "help": ""
    },
    {
      "name": "d_krs",
      "type": "Real",
      "xmlname": "krs",
      "help": ""
    },
    {
      "name": "d_alpha",
      "type": "Real",
      "xmlname": "alpha",
      "help": "Parameter alpha"
    },
    {
      "name": "d_n0",
      "type": "Real",
      "xmlname": "n0",
      "help": "Parameter n0"
    },
    {
      "name": "d_Es",
      "type": "Real",
      "xmlname": "Es",
      "help": "Parameter Es"
    },
    {
      "name": "d_mu",
      "type": "Real",
      "xmlname": "mu",
      "help": "Parameter mu"
    },
    {
      "name": "d_ini_cond_v",
      "type": "VecDeriv2",
      "xmlname": "ini_cond_v",
      "help": "Initial condition of the contraction ODE at each gauss point"
    }
  ],
  "links": [],
  "methods": [
    {
      "name": "init",
      "return_type": "void",
      "params": [],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "public"
    },
    {
      "name": "bwdInit",
      "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": "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": "step1_compute_strain_fiber",
      "return_type": "void",
      "params": [
        {
          "name": "x",
          "type": "const DataVecCoord &"
        }
      ],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "public"
    },
    {
      "name": "step2_solve_active_stress",
      "return_type": "void",
      "params": [],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "public"
    },
    {
      "name": "step3_add_force",
      "return_type": "void",
      "params": [
        {
          "name": "f",
          "type": "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": "updateMatrixData",
      "return_type": "void",
      "params": [],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "public"
    },
    {
      "name": "createODEInformation",
      "return_type": "void",
      "params": [
        {
          "name": "",
          "type": "int"
        },
        {
          "name": "t",
          "type": "HuxleyODEInformation &"
        }
      ],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "public"
    },
    {
      "name": "averageValueAtTetrahedronCenter",
      "return_type": "Real",
      "params": [
        {
          "name": "tetra",
          "type": "const Tetrahedron &"
        },
        {
          "name": "valueAtNode",
          "type": "const VecReal &"
        }
      ],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "public"
    },
    {
      "name": "get_electrical_input",
      "return_type": "Real",
      "params": [
        {
          "name": "info",
          "type": "const HuxleyODEInformation &"
        }
      ],
      "is_virtual": false,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "public"
    }
  ],
  "description": "The `SarcomereContractionForceField` is a specialized force field designed for simulating muscle contraction in tetrahedral meshes within the SOFA framework. It applies functional contraction forces along fiber directions, where stress is computed as a function of strain. The component inherits from `BaseTetraForceField`, indicating it operates on tetrahedron-based topologies.\n\n### Role and Purpose\n- **Simulation**: Simulates muscle contraction dynamics in 3D models using tetrahedral meshes.\n- **Mechanics**: Incorporates electromechanical coupling by solving a system of ODEs to determine the active stress due to sarcomere contractions. The contraction force is added along fiber directions based on computed strain and stress values.\n\n### Interactions with Other Components\n- **`BaseTetraForceField`**: Inherits methods for handling tetrahedral topologies, including initialization of tetrahedron information and edge data.\n- **Topology Data**: Requires access to the mesh topology (`l_topology`) to compute strain at each node and integrate forces into the system.\n- **Mechanical State**: Uses mechanical state information from `core::behavior::MechanicalState` for position and velocity data, essential for force computation.\n\n### Practical Usage Guidance\n- **Input Parameters**:\n  - `d_depol`: Depolarization time at each node (from loader).\n  - `d_apd`: Action potential duration at each node (from loader).\n  - Various material parameters (`sigma0`, `k0`, `heart_period`, `katp`, `krs`, `Es`, `mu`, `alpha`, and `n0`) to configure the contraction behavior.\n\n- **Data Fields**:\n  - `d_ini_cond_v`: Initial condition of the contraction ODE at each Gauss point, essential for setting up the electromechanical coupling.\n  - `m_ODEInfo`: Information related to electromechanics per tetrahedron, including strain (`e1D`), stiffness (`kc`), and contractile stress (`tauc`).\n\n### Methods Overview\n- **Initialization**:\n  - `init()`: Initializes the force field by setting up ODE information for each tetrahedron based on input parameters.\n- **Force Application**:\n  - `addForce()`: Computes strain, solves active stress via an ODE system, and adds contraction forces to the model.\n  - `step1_compute_strain_fiber()` calculates Green-Lagrange strain along fiber directions.\n  - `step2_solve_active_stress()` resolves the coupled electromechanical system for each tetrahedron.\n  - `step3_add_force()` integrates computed active stresses into the mechanical forces.\n- **Matrix Updates**:\n  - `updateMatrixData()`: Updates matrix data used in force calculations, ensuring accurate integration of contraction forces.",
  "maths": "<h2>Governing Equations and Operators</h2>\n\n<p>The <code>SarcomereContractionForceField</code> component implements a specialized force field for simulating muscle contraction in tetrahedral meshes within the SOFA framework. It applies functional contraction forces along fiber directions where stress is computed as a function of strain. The mathematical and physical principles underlying this implementation are outlined below.</p>\n\n<h3>Contraction Stress Calculation</h3>\n\n<p>The component computes active stress ($\\sigma_c$) due to muscle contraction by solving an ODE system for each tetrahedron in the mesh. The ODE describes the evolution of fiber strain ($e_{1D}$), stiffness ($k_c$), and contractile stress ($\\tau_c$). The following steps outline the mathematical formulation:</p>\n\n<ol>\n<li><strong>Strain Calculation</strong>: Given the deformation gradient $F$, the Green-Lagrange strain tensor $E$ is computed as:\n  \n\begin{equation}\n    E = \\frac{1}{2} (C - I) = \\frac{1}{2} (F^T F - I)\n\tag{1}\\label{eq:GLstrain}\n\\end{equation}\n\n<p>where $I$ is the identity tensor and $C = F^T F$. The strain along the fiber direction ($e_{1D}$) is then projected from the Green-Lagrange strain tensor:</p>\n  \n\begin{equation}\n    e_{1D} = \\mathbf{f}_d^T E \\mathbf{f}_d\n\tag{2}\\label{eq:e1D}\n\\end{equation}\n\n<p>where $\\mathbf{f}_d$ is the fiber direction vector.</p>\n</li>\n<li><strong>Solving the ODE for Active Stress</strong>: The active stress ($\\sigma_c$) is computed based on the evolution of the strain and stiffness parameters. For each tetrahedron, an ODE system determines the values of $k_c$, $\\tau_c$, and $e_{1D}$. Specifically, the new value for $e_{c}$ (fiber strain) at time step $n+1$ is given by:\n  \n\begin{equation}\n    e_{c}^{n+1} = e_{c}^n \\left(1 - \\Delta t \\frac{E_s}{\\mu}\\right) + \\Delta t \\frac{E_s}{\\mu} (e_{1D} - \\tau_c)\n\tag{3}\\label{eq:ec}\n\\end{equation}\n\n<p>where $E_s$ is the elasticity modulus of the elastic component and $\\mu$ is the sarcomere viscosity. The active stress ($\\sigma_c$) is then computed as:</p>\n  \n\begin{equation}\n    \\sigma_c = E_s (e_{1D} - e_{c})\n\tag{4}\\label{eq:sigma_c}\n\\end{equation}\n</li>\n<li><strong>Force Integration</strong>: The computed active stress is then integrated into the mechanical forces acting on each tetrahedron. This involves calculating the contraction force for each node based on the strain and deformation gradient:</p>\n  \n\begin{equation}\n    \\mathbf{f}_{\\text{contraction}} = -V_{rest} F^{-1} (0.25 \\mathbf{T}_c \\sigma_c) S_j\n\tag{5}\\label{eq:force}\n\\end{equation}\n\n<p>where $F^{-1}$ is the inverse of the deformation gradient, $S_j$ are shape functions, $V_{rest}$ is the rest volume, and $\\mathbf{T}_c$ is a contraction tensor that describes the fiber orientation.</p>\n</li>\n</ol>\n\n<h3>Constitutive Laws and Kinematic Models</h3>\n\n<p>The <code>SarcomereContractionForceField</code> incorporates several constitutive laws and kinematic models:</p>\n\n<ul>\n<li><strong>Fiber Direction Projection**: The Green-Lagrange strain tensor is projected onto the fiber direction to determine $e_{1D}$ (Equation \\ref{eq:e1D}).</li>\n<li><strong>Active Stress Computation**: The active stress is determined from the ODE system describing the evolution of sarcomere stiffness and contractile stress (Equations \\ref{eq:ec} and \\ref{eq:sigma_c}).</li>\n<li><strong>Force Integration**: The force due to contraction is computed by integrating the active stress along fiber directions using the deformation gradient and shape functions (Equation \\ref{eq:force}).</li>\n</ul>\n\n<h3>Role in the Global FEM Pipeline</h3>\n\n<p>The <code>SarcomereContractionForceField</code> plays a key role in several stages of the FEM simulation pipeline:</p>\n\n<ul>\n<li><strong>Initialization**: During initialization, it sets up ODE information for each tetrahedron based on input parameters (e.g., depolarization times and action potential durations).</li>\n<li><strong>Force Computation**: It computes strain along fiber directions and solves the active stress through an ODE system.</li>\n<li><strong>Force Integration**: The computed forces are integrated into the mechanical model, affecting nodal displacements in the subsequent time integration steps.</li>\n</ul>\n\n<h3>Numerical Methods and Discretization Choices</h3>\n\n<p>The component employs several numerical methods for discretization:</p>\n\n<ul>\n<li><strong>Implicit Time Integration**: The active stress is computed using an implicit ODE solver, ensuring stability in the simulation.</li>\n<li><strong>Tetrahedral Elements**: Tetrahedral elements are used to represent the mesh topology and facilitate strain calculations along fiber directions.</li>\n<li><strong>Gauss Points**: Strain and stress computations are performed at Gauss points within each tetrahedron to ensure accuracy.</li>\n</ul>\n\n<h3>Fitting into Variational/Lagrangian Mechanics Framework</h3>\n\n<p>The <code>SarcomereContractionForceField</code> fits into the broader variational and Lagrangian mechanics framework by:</p>\n\n<ul>\n<li><strong>Variational Formulation**: The strain energy is derived from a variational principle, where forces are computed as gradients of potential energy.</li>\n<li><strong>Lagrangian Mechanics**: The ODE system for active stress is formulated in terms of Lagrangian coordinates to capture the evolution of strain and stiffness parameters over time.</li>\n</ul>\n\n<p>This component effectively bridges analytical mechanics with practical numerical simulation by providing a physically consistent and numerically stable model for simulating muscle contraction dynamics within tetrahedral meshes.</p>",
  "abstract": "The `SarcomereContractionForceField` applies functional contraction forces along fiber directions in tetrahedral meshes by computing stress from strain using an ODE system for each tetrahedron.",
  "sheet": "\n# SarcomereContractionForceField\n\n## Overview\n\nThe `SarcomereContractionForceField` is a specialized force field designed to simulate muscle contraction dynamics within tetrahedral meshes in the SOFA framework. It applies functional contraction forces along fiber directions, where stress is computed as a function of strain using an ODE system for each tetrahedron.\n\n## Mathematical Model\n\nThe component computes active stress ($\\sigma_c$) due to muscle contraction by solving an ODE system for each tetrahedron in the mesh. The following steps outline the mathematical formulation:\n\n1. **Strain Calculation**: Given the deformation gradient $F$, the Green-Lagrange strain tensor $E$ is computed as:\n   \n   \\[ E = \\frac{1}{2} (C - I) = \\frac{1}{2} (F^T F - I) \\]\n\n   where $I$ is the identity tensor and $C = F^T F$. The strain along the fiber direction ($e_{1D}$) is then projected from the Green-Lagrange strain tensor:\n   \n   \\[ e_{1D} = \\mathbf{f}_d^T E \\mathbf{f}_d \\]\n\n2. **Solving the ODE for Active Stress**: The active stress ($\\sigma_c$) is computed based on the evolution of strain and stiffness parameters. For each tetrahedron, an ODE system determines the values of $k_c$, $\\tau_c$, and $e_{1D}$. Specifically, the new value for $e_{c}$ (fiber strain) at time step $n+1$ is given by:\n   \n   \\[ e_{c}^{n+1} = e_{c}^n \\left(1 - \\Delta t \\frac{E_s}{\\mu}\\right) + \\Delta t \\frac{E_s}{\\mu} (e_{1D} - \\tau_c) \\]\n\n   where $E_s$ is the elasticity modulus of the elastic component and $\\mu$ is the sarcomere viscosity. The active stress ($\\sigma_c$) is then computed as:\n   \n   \\[ \\sigma_c = E_s (e_{1D} - e_{c}) \\]\n\n3. **Force Integration**: The computed active stress is integrated into the mechanical forces acting on each tetrahedron, given by:\n   \n   \\[ \\mathbf{f}_{\\text{contraction}} = -V_{rest} F^{-1} (0.25 \\mathbf{T}_c \\sigma_c) S_j \\]\n\n   where $F^{-1}$ is the inverse of the deformation gradient, $S_j$ are shape functions, $V_{rest}$ is the rest volume, and $\\mathbf{T}_c$ is a contraction tensor that describes the fiber orientation.\n\n## Parameters and Data\n\nThe significant data fields exposed by the component include:\n\n- `d_depol`: Depolarization time at each node (from loader).\n- `d_apd`: Action potential duration at each node (from loader).\n- `d_sigma0`, `d_k0`, `d_hp`, `d_katp`, `d_krs`, `d_alpha`, `d_n0`, `d_Es`, and `d_mu`: Material parameters to configure the contraction behavior.\n- `d_ini_cond_v`: Initial condition of the contraction ODE at each Gauss point, essential for setting up the electromechanical coupling.\n\n## Practical Notes\n\nThe component requires careful configuration of input parameters such as depolarization times and action potential durations. Numerical stability is ensured through implicit time integration methods used in solving the ODE system."
}