TaitSurfacePressureForceField
\ This component computes the volume enclosed by a surface mesh \ and apply a pressure force following Tait's equation: $P = P_0 - B((V/V_0)^\\gamma - 1)$.\n\ This ForceField can be used to apply :\n\ * a constant pressure (set $B=0$ and use $P_0$)\n\ * an ideal gas pressure (set $\\gamma=1$ and use $B$)\n\ * a pressure from water (set $\\gamma=7$ and use $B$) This component computes the volume enclosed by a surface mesh and apply a pressure force following Tait's equation: $P = P_0 - B((V/V_0)^\gamma - 1)$. This ForceField can be used to apply : * a constant pressure (set $B=0$ and use $P_0$) * an ideal gas pressure (set $\gamma=1$ and use $B$) * a pressure from water (set $\gamma=7$ and use $B$)
The TaitSurfacePressureForceField computes the volume enclosed by a surface mesh and applies pressure forces based on Tait's equation: $P = P_0 - B((V/V_0)^\gamma - 1)$. It supports various types of pressures, including constant, ideal gas, and water.
- module
- Sofa.Component.MechanicalLoad
- namespace
- sofa::component::mechanicalload
- include
- sofa/component/mechanicalload/TaitSurfacePressureForceField.h
- inherits
-
- ForceField
- templates
-
- sofa::defaulttype::Vec3Types
- description
The TaitSurfacePressureForceField in the SOFA framework computes pressure forces on a surface mesh based on Tait's equation, which is given by:
egin{equation}
P = P_0 - B igg(rac{V}{V_0} igg)^eta - 1 \right)
\end{equation}
where:
- $P$ is the pressure,
- $P_0$ is the rest pressure when $V=V_0$ (the reference volume),
- $B$ is the bulk modulus (resistance to uniform compression),
- $\beta$ is a compressibility factor, and
- $V_0$ is the initial volume.
This equation can be used to model various types of pressure forces:
- Constant Pressure: Set $B = 0$, and use only $P_0$. The pressure becomes constant and independent of the current volume.
- Ideal Gas Pressure: Set $\beta=1$ and use $B$. This simplifies Tait's equation to a linear relationship with respect to the volume change.
- Water Pressure: Set $\beta=7$ and use $B$. This models the pressure behavior for water or similar incompressible fluids, where the bulk modulus is large and the compressibility factor is significant.
Governing Equations and Operators:
-
*Pressure Computation:
The current volume $V$ enclosed by the surface mesh is computed from the positions of its vertices. Tait's equation then determines the pressure force $F$ applied at each point on the surface. This computation involves calculating the enclosed volume by summing up the contributions from all triangles in the surface. -
*Force and Stiffness Matrices:
- The
addForcefunction calculates the total force due to the computed pressure:
egin{equation}
F_i = P \cdot n_i / 6,
\end{equation}
where $F_i$ is the force at vertex $i$, $P$ is the current pressure, and $n_i$ is the normal vector to the triangle containing vertex $i$. - The
addDForcefunction calculates the derivative of this force with respect to displacements:
egin{equation}
dF = K \cdot dV,
\end{equation}
where $dF$ is the change in force, $K$ is the stiffness matrix, and $dV$ is the change in volume. - The
addKToMatrixfunction computes the contribution to the global stiffness matrix from the pressure forces. This involves constructing a stiffness matrix that reflects how changes in vertex positions affect the enclosed volume and thus the pressure force.
- The
Constitutive Laws:
- Volume Calculation: For each triangle with vertices $a, b,$ and $c$, the volume increment is computed as:
egin{equation}
\Delta V = \frac{1}{6} (\mathbf{b}- extbf{a}) imes ( extbf{c}- extbf{a}) ullet extbf{x}_i,
\end{equation}
where $x_i$ is the position vector of a vertex. - Stiffness Calculation: The stiffness matrix entry for each triangle, contributing to the global matrix, is computed as:
egin{equation}
K_{ij} = -B\beta \left(\frac{V}{V_0}\right)^{eta-1}.
\end{equation}
Numerical Methods and Discretization Choices:
- Volume Calculation: Uses the divergence theorem to compute the volume of a closed surface from its triangular mesh representation. The normal vectors at each vertex are computed and summed up to determine the enclosed volume.
- Force Application: The pressure force is distributed uniformly over all vertices, with contributions weighted by the areas subtended by the triangles sharing those vertices.
Role in FEM Pipeline:
- This component integrates into the broader SOFA simulation pipeline as follows:
- Computes and applies forces during
addForcecalls, contributing to the nonlinear residual $R$. - Updates internal state during events (e.g., updating injected volume).
- Contributes to stiffness matrices in
addKToMatrix, which are assembled into the global tangent stiffness matrix $K(x)$ for the nonlinear solve step.
- Computes and applies forces during
- The component ensures that force and volume calculations are consistent with variational principles by maintaining the weak form of the governing PDE.
Summary:
The TaitSurfacePressureForceField implements a physically motivated model to apply pressure forces based on Tait's equation, which accounts for various compressible and incompressible fluid behaviors. The component integrates seamlessly into SOFA’s FEM simulation framework by providing contributions to force and stiffness matrices that are consistent with the weak form of continuum mechanics equations.
Data Fields
| Name | Type | Default | Help |
|---|---|---|---|
d_p0 |
Real | |
IN: Rest pressure when V = V0 |
d_B |
Real | |
IN: Bulk modulus (resistance to uniform compression) |
d_gamma |
Real | |
IN: Bulk modulus (resistance to uniform compression) |
d_injectedVolume |
Real | |
IN: Injected (or extracted) volume since the start of the simulation |
d_maxInjectionRate |
Real | |
IN: Maximum injection rate (volume per second) |
d_initialVolume |
Real | |
OUT: Initial volume, as computed from the surface rest position |
d_currentInjectedVolume |
Real | |
OUT: Current injected (or extracted) volume (taking into account maxInjectionRate) |
d_v0 |
Real | |
OUT: Rest volume (as computed from initialVolume + injectedVolume) |
d_currentVolume |
Real | |
OUT: Current volume, as computed from the last surface position |
d_currentPressure |
Real | |
OUT: Current pressure, as computed from the last surface position |
d_currentStiffness |
Real | |
OUT: dP/dV at current volume and pressure |
d_pressureTriangles |
SeqTriangles | |
OUT: list of triangles where a pressure is applied (mesh triangles + tessellated quads) |
d_initialSurfaceArea |
Real | |
OUT: Initial surface area, as computed from the surface rest position |
d_currentSurfaceArea |
Real | |
OUT: Current surface area, as computed from the last surface position |
d_drawForceScale |
Real | |
DEBUG: scale used to render force vectors |
d_drawForceColor |
sofa::type::RGBAColor | |
DEBUG: color used to render force vectors |
d_volumeAfterTC |
Real | |
OUT: Volume after a topology change |
d_surfaceAreaAfterTC |
Real | |
OUT: Surface area after a topology change |
Links
| Name | Type | Help |
|---|---|---|
l_topology |
link to the topology container |
Methods
void
init
()
void
storeResetState
()
void
reset
()
void
handleEvent
(core::objectmodel::Event * event)
void
addForce
(const core::MechanicalParams * mparams, DataVecDeriv & d_f, const DataVecCoord & d_x, const DataVecDeriv & d_v)
void
addDForce
(const core::MechanicalParams * mparams, DataVecDeriv & d_df, const DataVecDeriv & d_dx)
void
addKToMatrix
(const core::MechanicalParams * mparams, const sofa::core::behavior::MultiMatrixAccessor * matrix)
SReal
getPotentialEnergy
(const core::MechanicalParams * , const DataVecCoord & )
void
addKToMatrixT
(const core::MechanicalParams * mparams, MatrixWriter mwriter)
void
buildStiffnessMatrix
(sofa::core::behavior::StiffnessMatrix * matrix)
void
buildDampingMatrix
(core::behavior::DampingMatrix * )
void
draw
(const core::visual::VisualParams * vparams)
void
updateFromTopology
()
virtual
void
computePressureTriangles
()
virtual
void
computeMeshVolumeAndArea
(Real & volume, Real & area, const int & x)
virtual
void
computePressureAndStiffness
(Real & pressure, Real & stiffness, Real currentVolume, Real v0)
void
computeStatistics
(const int & x)
virtual
{
"name": "TaitSurfacePressureForceField",
"namespace": "sofa::component::mechanicalload",
"module": "Sofa.Component.MechanicalLoad",
"include": "sofa/component/mechanicalload/TaitSurfacePressureForceField.h",
"doc": "\\ This component computes the volume enclosed by a surface mesh \\ and apply a pressure force following Tait's equation: $P = P_0 - B((V/V_0)^\\\\gamma - 1)$.\\n\\ This ForceField can be used to apply :\\n\\ * a constant pressure (set $B=0$ and use $P_0$)\\n\\ * an ideal gas pressure (set $\\\\gamma=1$ and use $B$)\\n\\ * a pressure from water (set $\\\\gamma=7$ and use $B$)\n\nThis component computes the volume enclosed by a surface mesh\nand apply a pressure force following Tait's equation: $P = P_0 - B((V/V_0)^\\gamma - 1)$.\nThis ForceField can be used to apply :\n * a constant pressure (set $B=0$ and use $P_0$)\n * an ideal gas pressure (set $\\gamma=1$ and use $B$)\n * a pressure from water (set $\\gamma=7$ and use $B$)",
"inherits": [
"ForceField"
],
"templates": [
"sofa::defaulttype::Vec3Types"
],
"data_fields": [
{
"name": "d_p0",
"type": "Real",
"xmlname": "p0",
"help": "IN: Rest pressure when V = V0"
},
{
"name": "d_B",
"type": "Real",
"xmlname": "B",
"help": "IN: Bulk modulus (resistance to uniform compression)"
},
{
"name": "d_gamma",
"type": "Real",
"xmlname": "gamma",
"help": "IN: Bulk modulus (resistance to uniform compression)"
},
{
"name": "d_injectedVolume",
"type": "Real",
"xmlname": "injectedVolume",
"help": "IN: Injected (or extracted) volume since the start of the simulation"
},
{
"name": "d_maxInjectionRate",
"type": "Real",
"xmlname": "maxInjectionRate",
"help": "IN: Maximum injection rate (volume per second)"
},
{
"name": "d_initialVolume",
"type": "Real",
"xmlname": "initialVolume",
"help": "OUT: Initial volume, as computed from the surface rest position"
},
{
"name": "d_currentInjectedVolume",
"type": "Real",
"xmlname": "currentInjectedVolume",
"help": "OUT: Current injected (or extracted) volume (taking into account maxInjectionRate)"
},
{
"name": "d_v0",
"type": "Real",
"xmlname": "v0",
"help": "OUT: Rest volume (as computed from initialVolume + injectedVolume)"
},
{
"name": "d_currentVolume",
"type": "Real",
"xmlname": "currentVolume",
"help": "OUT: Current volume, as computed from the last surface position"
},
{
"name": "d_currentPressure",
"type": "Real",
"xmlname": "currentPressure",
"help": "OUT: Current pressure, as computed from the last surface position"
},
{
"name": "d_currentStiffness",
"type": "Real",
"xmlname": "currentStiffness",
"help": "OUT: dP/dV at current volume and pressure"
},
{
"name": "d_pressureTriangles",
"type": "SeqTriangles",
"xmlname": "pressureTriangles",
"help": "OUT: list of triangles where a pressure is applied (mesh triangles + tessellated quads)"
},
{
"name": "d_initialSurfaceArea",
"type": "Real",
"xmlname": "initialSurfaceArea",
"help": "OUT: Initial surface area, as computed from the surface rest position"
},
{
"name": "d_currentSurfaceArea",
"type": "Real",
"xmlname": "currentSurfaceArea",
"help": "OUT: Current surface area, as computed from the last surface position"
},
{
"name": "d_drawForceScale",
"type": "Real",
"xmlname": "drawForceScale",
"help": "DEBUG: scale used to render force vectors"
},
{
"name": "d_drawForceColor",
"type": "sofa::type::RGBAColor",
"xmlname": "drawForceColor",
"help": "DEBUG: color used to render force vectors"
},
{
"name": "d_volumeAfterTC",
"type": "Real",
"xmlname": "volumeAfterTC",
"help": "OUT: Volume after a topology change"
},
{
"name": "d_surfaceAreaAfterTC",
"type": "Real",
"xmlname": "surfaceAreaAfterTC",
"help": "OUT: Surface area after a topology change"
}
],
"links": [
{
"name": "l_topology",
"target": "BaseMeshTopology",
"kind": "single",
"xmlname": "topology",
"help": "link to the topology container"
}
],
"methods": [
{
"name": "init",
"return_type": "void",
"params": [],
"is_virtual": false,
"is_pure_virtual": false,
"is_static": false,
"access": "public"
},
{
"name": "storeResetState",
"return_type": "void",
"params": [],
"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": "handleEvent",
"return_type": "void",
"params": [
{
"name": "event",
"type": "core::objectmodel::Event *"
}
],
"is_virtual": false,
"is_pure_virtual": false,
"is_static": false,
"access": "public"
},
{
"name": "addForce",
"return_type": "void",
"params": [
{
"name": "mparams",
"type": "const core::MechanicalParams *"
},
{
"name": "d_f",
"type": "DataVecDeriv &"
},
{
"name": "d_x",
"type": "const DataVecCoord &"
},
{
"name": "d_v",
"type": "const DataVecDeriv &"
}
],
"is_virtual": false,
"is_pure_virtual": false,
"is_static": false,
"access": "public"
},
{
"name": "addDForce",
"return_type": "void",
"params": [
{
"name": "mparams",
"type": "const core::MechanicalParams *"
},
{
"name": "d_df",
"type": "DataVecDeriv &"
},
{
"name": "d_dx",
"type": "const DataVecDeriv &"
}
],
"is_virtual": false,
"is_pure_virtual": false,
"is_static": false,
"access": "public"
},
{
"name": "addKToMatrix",
"return_type": "void",
"params": [
{
"name": "mparams",
"type": "const core::MechanicalParams *"
},
{
"name": "matrix",
"type": "const sofa::core::behavior::MultiMatrixAccessor *"
}
],
"is_virtual": false,
"is_pure_virtual": false,
"is_static": false,
"access": "public"
},
{
"name": "getPotentialEnergy",
"return_type": "SReal",
"params": [
{
"name": "",
"type": "const core::MechanicalParams *"
},
{
"name": "",
"type": "const DataVecCoord &"
}
],
"is_virtual": false,
"is_pure_virtual": false,
"is_static": false,
"access": "public"
},
{
"name": "addKToMatrixT",
"return_type": "void",
"params": [
{
"name": "mparams",
"type": "const core::MechanicalParams *"
},
{
"name": "mwriter",
"type": "MatrixWriter"
}
],
"is_virtual": false,
"is_pure_virtual": false,
"is_static": false,
"access": "public"
},
{
"name": "buildStiffnessMatrix",
"return_type": "void",
"params": [
{
"name": "matrix",
"type": "sofa::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": "updateFromTopology",
"return_type": "void",
"params": [],
"is_virtual": true,
"is_pure_virtual": false,
"is_static": false,
"access": "protected"
},
{
"name": "computePressureTriangles",
"return_type": "void",
"params": [],
"is_virtual": true,
"is_pure_virtual": false,
"is_static": false,
"access": "protected"
},
{
"name": "computeMeshVolumeAndArea",
"return_type": "void",
"params": [
{
"name": "volume",
"type": "Real &"
},
{
"name": "area",
"type": "Real &"
},
{
"name": "x",
"type": "const int &"
}
],
"is_virtual": true,
"is_pure_virtual": false,
"is_static": false,
"access": "protected"
},
{
"name": "computePressureAndStiffness",
"return_type": "void",
"params": [
{
"name": "pressure",
"type": "Real &"
},
{
"name": "stiffness",
"type": "Real &"
},
{
"name": "currentVolume",
"type": "Real"
},
{
"name": "v0",
"type": "Real"
}
],
"is_virtual": false,
"is_pure_virtual": false,
"is_static": false,
"access": "protected"
},
{
"name": "computeStatistics",
"return_type": "void",
"params": [
{
"name": "x",
"type": "const int &"
}
],
"is_virtual": true,
"is_pure_virtual": false,
"is_static": false,
"access": "protected"
}
],
"description": "The `TaitSurfacePressureForceField` is a component in the SOFA framework that computes and applies pressure forces on a surface mesh based on Tait's equation: $P = P_0 - B((V/V_0)^\\gamma - 1)$. It supports different types of pressures, such as constant, ideal gas, or water pressure by adjusting parameters like bulk modulus ($B$), rest pressure ($P_0$), and compressibility factor ($γ$). The component relies on a `BaseMeshTopology` for the mesh topology information. \n\nIt calculates the current volume enclosed by the surface mesh and computes the corresponding pressure based on Tait's equation, applying this force to the mechanical system. Key functionalities include:\n\n- Computing the initial and current volumes and surface areas.\n- Handling events to update injected volume over time.\n- Applying forces and constructing stiffness matrices for integration into the simulation pipeline.\n- Rendering visualization of applied pressures.\n\nThe component provides various data fields such as `p0`, `B`, `gamma`, and `injectedVolume` for user-defined parameters, while also providing computed outputs like `currentPressure`, `currentStiffness`, and volume statistics. It integrates well with other mechanical load components and simulation pipelines in SOFA.",
"maths": "The `TaitSurfacePressureForceField` in the SOFA framework computes pressure forces on a surface mesh based on Tait's equation, which is given by:\n\n\begin{equation}\n P = P_0 - B \bigg(\frac{V}{V_0} \bigg)^\beta - 1 \\right)\n\\end{equation}\n\nwhere:\n- $P$ is the pressure,\n- $P_0$ is the rest pressure when $V=V_0$ (the reference volume),\n- $B$ is the bulk modulus (resistance to uniform compression),\n- $\\beta$ is a compressibility factor, and\n- $V_0$ is the initial volume.\n\nThis equation can be used to model various types of pressure forces:\n- **Constant Pressure**: Set $B = 0$, and use only $P_0$. The pressure becomes constant and independent of the current volume.\n- **Ideal Gas Pressure**: Set $\\beta=1$ and use $B$. This simplifies Tait's equation to a linear relationship with respect to the volume change.\n- **Water Pressure**: Set $\\beta=7$ and use $B$. This models the pressure behavior for water or similar incompressible fluids, where the bulk modulus is large and the compressibility factor is significant.\n\n### Governing Equations and Operators:\n\n1. **Pressure Computation:*\n The current volume $V$ enclosed by the surface mesh is computed from the positions of its vertices. Tait's equation then determines the pressure force $F$ applied at each point on the surface. This computation involves calculating the enclosed volume by summing up the contributions from all triangles in the surface.\n\n2. **Force and Stiffness Matrices:*\n - The `addForce` function calculates the total force due to the computed pressure:\n \begin{equation}\n F_i = P \\cdot n_i / 6,\n \\end{equation}\n where $F_i$ is the force at vertex $i$, $P$ is the current pressure, and $n_i$ is the normal vector to the triangle containing vertex $i$.\n - The `addDForce` function calculates the derivative of this force with respect to displacements:\n \begin{equation}\n dF = K \\cdot dV,\n \\end{equation}\n where $dF$ is the change in force, $K$ is the stiffness matrix, and $dV$ is the change in volume.\n - The `addKToMatrix` function computes the contribution to the global stiffness matrix from the pressure forces. This involves constructing a stiffness matrix that reflects how changes in vertex positions affect the enclosed volume and thus the pressure force.\n\n### Constitutive Laws:\n- **Volume Calculation**: For each triangle with vertices $a, b,$ and $c$, the volume increment is computed as:\n \begin{equation}\n \\Delta V = \\frac{1}{6} (\\mathbf{b}-\textbf{a}) \times (\textbf{c}-\textbf{a}) \bullet \textbf{x}_i,\n \\end{equation}\n where $x_i$ is the position vector of a vertex.\n- **Stiffness Calculation**: The stiffness matrix entry for each triangle, contributing to the global matrix, is computed as:\n \begin{equation}\n K_{ij} = -B\\beta \\left(\\frac{V}{V_0}\\right)^{\beta-1}.\n \\end{equation}\n\n### Numerical Methods and Discretization Choices:\n- **Volume Calculation**: Uses the divergence theorem to compute the volume of a closed surface from its triangular mesh representation. The normal vectors at each vertex are computed and summed up to determine the enclosed volume.\n- **Force Application**: The pressure force is distributed uniformly over all vertices, with contributions weighted by the areas subtended by the triangles sharing those vertices.\n\n### Role in FEM Pipeline:\n- This component integrates into the broader SOFA simulation pipeline as follows:\n - Computes and applies forces during `addForce` calls, contributing to the nonlinear residual $R$.\n - Updates internal state during events (e.g., updating injected volume).\n - Contributes to stiffness matrices in `addKToMatrix`, which are assembled into the global tangent stiffness matrix $K(x)$ for the nonlinear solve step.\n- The component ensures that force and volume calculations are consistent with variational principles by maintaining the weak form of the governing PDE.\n\n### Summary:\nThe `TaitSurfacePressureForceField` implements a physically motivated model to apply pressure forces based on Tait's equation, which accounts for various compressible and incompressible fluid behaviors. The component integrates seamlessly into SOFA’s FEM simulation framework by providing contributions to force and stiffness matrices that are consistent with the weak form of continuum mechanics equations.",
"abstract": "The TaitSurfacePressureForceField computes the volume enclosed by a surface mesh and applies pressure forces based on Tait's equation: $P = P_0 - B((V/V_0)^\\gamma - 1)$. It supports various types of pressures, including constant, ideal gas, and water.",
"sheet": "\n# TaitSurfacePressureForceField\n\n## Overview\nThe `TaitSurfacePressureForceField` is a component in the SOFA framework that computes pressure forces on a surface mesh based on Tait's equation. It handles the computation of the enclosed volume and applies pressure forces according to different types, such as constant, ideal gas, or water pressure.\n\n## Mathematical Model\nThe `TaitSurfacePressureForceField` uses Tait's equation to compute the pressure force:\n\\[\nP = P_0 - B \\left(\\frac{V}{V_0}\\right)^\\gamma - 1 \\\n\\]\nwhere $P$ is the pressure, $P_0$ is the rest pressure when $V=V_0$, $B$ is the bulk modulus (resistance to uniform compression), and $\\gamma$ is a compressibility factor. The volume $V$ enclosed by the surface mesh is computed from the positions of its vertices.\n\n## Parameters and Data\nThe significant Data fields exposed by this component are:\n- `p0`: Rest pressure when $V=V_0$\n- `B`: Bulk modulus (resistance to uniform compression)\n- `gamma`: Compressibility factor\n- `injectedVolume`: Injected or extracted volume since the start of the simulation\n- `maxInjectionRate`: Maximum injection rate (volume per second)\n- `initialVolume`: Initial volume, as computed from the surface rest position\n- `currentInjectedVolume`: Current injected or extracted volume (taking into account maxInjectionRate)\n- `v0`: Rest volume (as computed from initialVolume + injectedVolume)\n- `currentVolume`: Current volume, as computed from the last surface position\n- `currentPressure`: Current pressure, as computed from the last surface position\n- `currentStiffness`: $dP/dV$ at current volume and pressure\n- `pressureTriangles`: List of triangles where a pressure is applied (mesh triangles + tessellated quads)\n- `initialSurfaceArea`: Initial surface area, as computed from the surface rest position\n- `currentSurfaceArea`: Current surface area, as computed from the last surface position\n- `drawForceScale`: Scale used to render force vectors\n- `drawForceColor`: Color used to render force vectors\n- `volumeAfterTC`: Volume after a topology change\n- `surfaceAreaAfterTC`: Surface area after a topology change\n\n## Dependencies and Connections\nThe component requires a link to the `BaseMeshTopology` for mesh topology information.\n\n## Practical Notes\nConfiguration caveats include setting appropriate values for $B$, $P_0$, and $\\gamma$ based on the desired pressure type (constant, ideal gas, or water). The maximum injection rate should be set carefully to avoid numerical instability."
}