Back

LCPConstraintSolver

sofa::component::constraint::lagrangian::solver::LCPConstraintSolver
ConstraintSolverImpl
Doc (from source)

A Constraint Solver using the Linear Complementarity Problem formulation to solve BaseConstraint based components.

Abstract (AI generated)

The `LCPConstraintSolver` uses Linear Complementarity Problem (LCP) formulation to solve constraints in SOFA simulations, handling mechanical interactions including frictional contacts.

Metadata
module
Sofa.Component.Constraint.Lagrangian.Solver
namespace
sofa::component::constraint::lagrangian::solver
include
sofa/component/constraint/lagrangian/solver/LCPConstraintSolver.h
inherits
  • ConstraintSolverImpl
description

Governing Equations and Operators

The LCPConstraintSolver in the SOFA framework is designed to solve constraint equations formulated as a Linear Complementarity Problem (LCP). The LCP formulation is essential for handling frictional contacts and constraints in deformable bodies. Key operators and governing equations include:

  • Mechanical Constraints: Constraints are represented by the equation: \[ \mathbf{d} = \mathbf{W}\boldsymbol{\lambda} - \mathbf{f}, \] where \(\mathbf{d}\) is the constraint violation, \(\mathbf{W}\) is the compliance matrix (inverse of the stiffness), and \(\boldsymbol{\lambda}\) are the Lagrange multipliers representing contact forces. The vector \(\mathbf{f}\) represents external force contributions.
  • LCP Formulation: The LCP is formulated as: \[ \begin{cases} \boldsymbol{d} = \mathbf{W}\boldsymbol{\lambda} - \mathbf{f}, \\ \boldsymbol{\lambda} \geq 0, \quad \boldsymbol{d} \geq 0, \\ \boldsymbol{\lambda}^T \cdot \boldsymbol{d} = 0. \end{cases}\] The constraints ensure that contact forces and violations are non-negative and complementary.

Constitutive or Kinematic Laws Involved

  • Friction: The frictional behavior is modeled using a Coulomb friction law, where the friction coefficient \(\mu\) governs tangential forces. In the presence of friction, each contact can be represented by a set of equations involving normal and tangential components.
  • Compliance Matrix: The compliance matrix \(\mathbf{W}\) is derived from the constitutive law governing material behavior. For linear elasticity, it typically corresponds to an inverse stiffness matrix. In nonlinear cases, such as hyperelastic materials, it may be more complex.
  • Damping Models: The solver can incorporate regularization terms to stabilize the solution process and avoid numerical instabilities. These terms often take the form of adding a small multiple of the identity matrix to \(\mathbf{W}\).

Role in the Global FEM Pipeline

  • Assembly Phase: During assembly, the solver gathers information on constraints and builds the compliance matrix \(\mathbf{W}\). This involves calling methods like addComplianceInConstraintSpace(), which sum up contributions from various constraint correction components.
  • Time Integration: The solver operates within an implicit time integration scheme, typically using the Backward Euler method. It computes contact forces and updates positions/velocities based on these forces.
  • Nonlinear Solution: The nonlinear equations are solved iteratively using methods like Gauss-Seidel or variants thereof. For frictional contacts, a modified version of the Gauss-Seidel algorithm (NLCP) is used to account for tangential forces and stick-slip transitions.

Numerical Methods or Discretization Choices Encoded

  • Iterative Solvers: The LCPConstraintSolver uses iterative solvers, primarily Gauss-Seidel methods. These are suitable for large systems due to their low memory footprint and simplicity.
  • Built vs Unbuilt LCP: Depending on the configuration, the full system matrix (LCP) can be fully built or partially constructed (unbuilt), with only necessary components being computed at each iteration. This choice allows for performance optimization in different scenarios.
  • Initial Guess and Regularization: The solver supports options to use previous solutions as initial guesses (d_initial_guess) and add regularization terms to the compliance matrix \(\mathbf{W}\). These choices help in stabilizing the iterative solution process.

Integration into Variational / Lagrangian Mechanics Framework

The LCPConstraintSolver fits into the broader variational mechanics framework by solving the discretized weak form of constraint equations. These constraints are derived from the strong PDEs governing continuum mechanics, which are then transformed into a system of algebraic equations via FEM. The solver ensures that these constraints are satisfied in each time step, thus preserving physical consistency and stability.

The LCP formulation is particularly useful for handling unilateral contact conditions and frictional interactions, where the Lagrange multipliers represent forces that can only be active under certain conditions (i.e., when objects come into contact). This approach maintains a variational structure, ensuring that the solution remains consistent with the underlying physical principles.

In summary, LCPConstraintSolver is a critical component for enforcing constraints in deformable body simulations within SOFA, integrating seamlessly into both the mechanical modeling and numerical discretization processes.

Data Fields
NameTypeDefaultHelp
d_displayDebug bool Display debug information.
d_initial_guess bool activate LCP results history to improve its resolution performances.
d_build_lcp bool LCP is not fully built to increase performance in some case.
d_tol SReal residual error threshold for termination of the Gauss-Seidel algorithm
d_maxIt int maximal number of iterations of the Gauss-Seidel algorithm
d_regularizationTerm SReal Add regularization factor times the identity matrix to the compliance W when solving constraints
d_mu SReal Friction coefficient
d_minW SReal If not zero, constraints whose self-compliance (i.e. the corresponding value on the diagonal of W) is smaller than this threshold will be ignored
d_maxF SReal If not zero, constraints whose response force becomes larger than this threshold will be ignored
d_computeConstraintForces bool enable the storage of the constraintForces.
d_showCellWidth SReal Distance between each constraint cells
d_showTranslation type::Vec3 Position of the first cell
Methods
bool prepareStates (const core::ConstraintParams * , MultiVecId res1, MultiVecId res2) virtual
bool buildSystem (const core::ConstraintParams * , MultiVecId res1, MultiVecId res2) virtual
bool solveSystem (const core::ConstraintParams * , MultiVecId res1, MultiVecId res2) virtual
bool applyCorrection (const core::ConstraintParams * , MultiVecId res1, MultiVecId res2) virtual
ConstraintProblem * getConstraintProblem () virtual
void lockConstraintProblem (sofa::core::objectmodel::BaseObject * from, ConstraintProblem * p1, ConstraintProblem * p2) virtual
{
  "name": "LCPConstraintSolver",
  "namespace": "sofa::component::constraint::lagrangian::solver",
  "module": "Sofa.Component.Constraint.Lagrangian.Solver",
  "include": "sofa/component/constraint/lagrangian/solver/LCPConstraintSolver.h",
  "doc": "A Constraint Solver using the Linear Complementarity Problem formulation to solve BaseConstraint based components.",
  "inherits": [
    "ConstraintSolverImpl"
  ],
  "templates": [],
  "data_fields": [
    {
      "name": "d_displayDebug",
      "type": "bool",
      "xmlname": "displayDebug",
      "help": "Display debug information."
    },
    {
      "name": "d_initial_guess",
      "type": "bool",
      "xmlname": "initial_guess",
      "help": "activate LCP results history to improve its resolution performances."
    },
    {
      "name": "d_build_lcp",
      "type": "bool",
      "xmlname": "build_lcp",
      "help": "LCP is not fully built to increase performance in some case."
    },
    {
      "name": "d_tol",
      "type": "SReal",
      "xmlname": "tolerance",
      "help": "residual error threshold for termination of the Gauss-Seidel algorithm"
    },
    {
      "name": "d_maxIt",
      "type": "int",
      "xmlname": "maxIt",
      "help": "maximal number of iterations of the Gauss-Seidel algorithm"
    },
    {
      "name": "d_regularizationTerm",
      "type": "SReal",
      "xmlname": "regularizationTerm",
      "help": "Add regularization factor times the identity matrix to the compliance W when solving constraints"
    },
    {
      "name": "d_mu",
      "type": "SReal",
      "xmlname": "mu",
      "help": "Friction coefficient"
    },
    {
      "name": "d_minW",
      "type": "SReal",
      "xmlname": "minW",
      "help": "If not zero, constraints whose self-compliance (i.e. the corresponding value on the diagonal of W) is smaller than this threshold will be ignored"
    },
    {
      "name": "d_maxF",
      "type": "SReal",
      "xmlname": "maxF",
      "help": "If not zero, constraints whose response force becomes larger than this threshold will be ignored"
    },
    {
      "name": "d_computeConstraintForces",
      "type": "bool",
      "xmlname": "computeConstraintForces",
      "help": "enable the storage of the constraintForces."
    },
    {
      "name": "d_showCellWidth",
      "type": "SReal",
      "xmlname": "showCellWidth",
      "help": "Distance between each constraint cells"
    },
    {
      "name": "d_showTranslation",
      "type": "type::Vec3",
      "xmlname": "showTranslation",
      "help": "Position of the first cell"
    }
  ],
  "links": [],
  "methods": [
    {
      "name": "prepareStates",
      "return_type": "bool",
      "params": [
        {
          "name": "",
          "type": "const core::ConstraintParams *"
        },
        {
          "name": "res1",
          "type": "MultiVecId"
        },
        {
          "name": "res2",
          "type": "MultiVecId"
        }
      ],
      "is_virtual": true,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "public"
    },
    {
      "name": "buildSystem",
      "return_type": "bool",
      "params": [
        {
          "name": "",
          "type": "const core::ConstraintParams *"
        },
        {
          "name": "res1",
          "type": "MultiVecId"
        },
        {
          "name": "res2",
          "type": "MultiVecId"
        }
      ],
      "is_virtual": true,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "public"
    },
    {
      "name": "solveSystem",
      "return_type": "bool",
      "params": [
        {
          "name": "",
          "type": "const core::ConstraintParams *"
        },
        {
          "name": "res1",
          "type": "MultiVecId"
        },
        {
          "name": "res2",
          "type": "MultiVecId"
        }
      ],
      "is_virtual": true,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "public"
    },
    {
      "name": "applyCorrection",
      "return_type": "bool",
      "params": [
        {
          "name": "",
          "type": "const core::ConstraintParams *"
        },
        {
          "name": "res1",
          "type": "MultiVecId"
        },
        {
          "name": "res2",
          "type": "MultiVecId"
        }
      ],
      "is_virtual": true,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "public"
    },
    {
      "name": "getConstraintProblem",
      "return_type": "ConstraintProblem *",
      "params": [],
      "is_virtual": true,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "public"
    },
    {
      "name": "lockConstraintProblem",
      "return_type": "void",
      "params": [
        {
          "name": "from",
          "type": "sofa::core::objectmodel::BaseObject *"
        },
        {
          "name": "p1",
          "type": "ConstraintProblem *"
        },
        {
          "name": "p2",
          "type": "ConstraintProblem *"
        }
      ],
      "is_virtual": true,
      "is_pure_virtual": false,
      "is_static": false,
      "access": "public"
    }
  ],
  "description": "The `LCPConstraintSolver` is a constraint solver in the SOFA framework that uses the Linear Complementarity Problem (LCP) formulation to solve constraints based on `BaseConstraint` components. It inherits from `ConstraintSolverImpl`, indicating its role in solving mechanical systems with constraints, particularly those involving contact or frictional forces. The component interacts with other components such as constraint correction and compliance calculation mechanisms through virtual methods like `buildSystem()`, `solveSystem()`, `applyCorrection()`, etc., allowing for flexible integration within the SOFA scene graph architecture.\n\nKey functionalities include:\n- **Build System**: Assembles the system of constraints, incorporating friction coefficients (`mu`), regularization terms to stabilize the solution process, and compliance matrices (`W`) that describe how objects interact under forces.\n- **Solve System**: Solves the LCP using iterative methods like Gauss-Seidel. It handles both built and unbuilt cases, where the full matrix is constructed or only parts of it are used for performance reasons.\n- **Apply Correction**: Applies the computed constraint forces back to the system, updating positions and velocities accordingly.\n\nThe component also supports options such as initial guess history (`initial_guess`), tolerance thresholds for convergence (`tolerance`, `maxIt`), and optional frictional behavior via a coefficient (`mu`). It allows users to configure how constraints are handled, including ignoring very stiff or weak constraints based on compliance values (`minW`) and force magnitudes (`maxF`).\n\nData fields like `displayDebug`, `computeConstraintForces`, and visualization parameters allow for debugging and monitoring the solver's behavior.",
  "maths": "<h3>Governing Equations and Operators</h3>\n\n<p>The <code>LCPConstraintSolver</code> in the SOFA framework is designed to solve constraint equations formulated as a Linear Complementarity Problem (LCP). The LCP formulation is essential for handling frictional contacts and constraints in deformable bodies. Key operators and governing equations include:</p>\n<ul>\n<li><strong>Mechanical Constraints:</strong> Constraints are represented by the equation:\n\\[ \\mathbf{d} = \\mathbf{W}\\boldsymbol{\\lambda} - \\mathbf{f}, \\]\nwhere <span class=\"math inline\">\\(\\mathbf{d}\\)</span> is the constraint violation, <span class=\"math inline\">\\(\\mathbf{W}\\)</span> is the compliance matrix (inverse of the stiffness), and <span class=\"math inline\">\\(\\boldsymbol{\\lambda}\\)</span> are the Lagrange multipliers representing contact forces. The vector <span class=\"math inline\">\\(\\mathbf{f}\\)</span> represents external force contributions.</li>\n<li><strong>LCP Formulation:</strong> The LCP is formulated as:\n\\[ \\begin{cases}\n    \\boldsymbol{d} = \\mathbf{W}\\boldsymbol{\\lambda} - \\mathbf{f}, \\\\\n    \\boldsymbol{\\lambda} \\geq 0, \\quad \\boldsymbol{d} \\geq 0, \\\\\n    \\boldsymbol{\\lambda}^T \\cdot \\boldsymbol{d} = 0.\n   \\end{cases}\\]\nThe constraints ensure that contact forces and violations are non-negative and complementary.</li>\n</ul>\n<h3>Constitutive or Kinematic Laws Involved</h3>\n<ul>\n<li><strong>Friction:</strong> The frictional behavior is modeled using a Coulomb friction law, where the friction coefficient <span class=\"math inline\">\\(\\mu\\)</span> governs tangential forces. In the presence of friction, each contact can be represented by a set of equations involving normal and tangential components.</li>\n<li><strong>Compliance Matrix:</strong> The compliance matrix <span class=\"math inline\">\\(\\mathbf{W}\\)</span> is derived from the constitutive law governing material behavior. For linear elasticity, it typically corresponds to an inverse stiffness matrix. In nonlinear cases, such as hyperelastic materials, it may be more complex.</li>\n<li><strong>Damping Models:</strong> The solver can incorporate regularization terms to stabilize the solution process and avoid numerical instabilities. These terms often take the form of adding a small multiple of the identity matrix to <span class=\"math inline\">\\(\\mathbf{W}\\)</span>.</li>\n</ul>\n<h3>Role in the Global FEM Pipeline</h3>\n<ul>\n<li><strong>Assembly Phase:</strong> During assembly, the solver gathers information on constraints and builds the compliance matrix <span class=\"math inline\">\\(\\mathbf{W}\\)</span>. This involves calling methods like <code>addComplianceInConstraintSpace()</code>, which sum up contributions from various constraint correction components.</li>\n<li><strong>Time Integration:</strong> The solver operates within an implicit time integration scheme, typically using the Backward Euler method. It computes contact forces and updates positions/velocities based on these forces.</li>\n<li><strong>Nonlinear Solution:</strong> The nonlinear equations are solved iteratively using methods like Gauss-Seidel or variants thereof. For frictional contacts, a modified version of the Gauss-Seidel algorithm (NLCP) is used to account for tangential forces and stick-slip transitions.</li>\n</ul>\n<h3>Numerical Methods or Discretization Choices Encoded</h3>\n<ul>\n<li><strong>Iterative Solvers:</strong> The <code>LCPConstraintSolver</code> uses iterative solvers, primarily Gauss-Seidel methods. These are suitable for large systems due to their low memory footprint and simplicity.</li>\n<li><strong>Built vs Unbuilt LCP:</strong> Depending on the configuration, the full system matrix (LCP) can be fully built or partially constructed (unbuilt), with only necessary components being computed at each iteration. This choice allows for performance optimization in different scenarios.</li>\n<li><strong>Initial Guess and Regularization:</strong> The solver supports options to use previous solutions as initial guesses (<code>d_initial_guess</code>) and add regularization terms to the compliance matrix <span class=\"math inline\">\\(\\mathbf{W}\\)</span>. These choices help in stabilizing the iterative solution process.</li>\n</ul>\n<h3>Integration into Variational / Lagrangian Mechanics Framework</h3>\n<p>The <code>LCPConstraintSolver</code> fits into the broader variational mechanics framework by solving the discretized weak form of constraint equations. These constraints are derived from the strong PDEs governing continuum mechanics, which are then transformed into a system of algebraic equations via FEM. The solver ensures that these constraints are satisfied in each time step, thus preserving physical consistency and stability.</p>\n<p>The LCP formulation is particularly useful for handling unilateral contact conditions and frictional interactions, where the Lagrange multipliers represent forces that can only be active under certain conditions (i.e., when objects come into contact). This approach maintains a variational structure, ensuring that the solution remains consistent with the underlying physical principles.</p>\n<p>In summary, <code>LCPConstraintSolver</code> is a critical component for enforcing constraints in deformable body simulations within SOFA, integrating seamlessly into both the mechanical modeling and numerical discretization processes.</p>",
  "abstract": "The `LCPConstraintSolver` uses Linear Complementarity Problem (LCP) formulation to solve constraints in SOFA simulations, handling mechanical interactions including frictional contacts.",
  "sheet": "# LCPConstraintSolver\n\n## Overview\nThe `LCPConstraintSolver` is a constraint solver that uses the Linear Complementarity Problem (LCP) formulation to handle constraints based on `BaseConstraint` components. It integrates into the SOFA scene graph architecture through methods like `buildSystem()`, `solveSystem()`, and `applyCorrection()`.\n\n## Mathematical Model\nThe LCP formulation is used to solve constraint equations:\n\\[ \\mathbf{d} = \\mathbf{W}\\boldsymbol{\\lambda} - \\mathbf{f}, \\]\nwhere \\(\\mathbf{d}\\) is the constraint violation, \\(\\mathbf{W}\\) is the compliance matrix (inverse of stiffness), and \\(\\boldsymbol{\\lambda}\\) are Lagrange multipliers representing contact forces. The LCP constraints ensure non-negative forces and violations:\n\\[ \\begin{cases}\n    \\boldsymbol{d} = \\mathbf{W}\\boldsymbol{\\lambda} - \\mathbf{f}, \\\\\n    \\boldsymbol{\\lambda} \\geq 0, \\quad \\boldsymbol{d} \\geq 0, \\\\\n    \\boldsymbol{\\lambda}^T \\cdot \\boldsymbol{d} = 0.\n   \\end{cases}\\]\nThe frictional behavior is modeled using a Coulomb friction law with coefficient \\(\\mu\\).\n\n## Parameters and Data\n- **displayDebug**: Display debug information (`bool`, default: `false`)\n- **initial_guess**: Use previous solutions as initial guesses to improve performance (`bool`, default: `false`)\n- **build_lcp**: Build the full LCP system for increased performance in some cases (`bool`, default: `true`)\n- **tolerance**: Residual error threshold for termination of the Gauss-Seidel algorithm (`SReal`, default: 1e-6)\n- **maxIt**: Maximal number of iterations of the Gauss-Seidel algorithm (`int`, default: 100)\n- **regularizationTerm**: Add regularization factor to stabilize the solution process (`SReal`, default: 0.0)\n- **mu**: Friction coefficient (`SReal`, default: 0.5)\n- **minW**: Threshold for ignoring constraints with very small self-compliance (`SReal`, default: 0.0)\n- **maxF**: Threshold for ignoring constraints with large response forces (`SReal`, default: 0.0)\n- **computeConstraintForces**: Enable storage of constraint forces (`bool`, default: `false`)\n- **showCellWidth**: Distance between each constraint cell for visualization (`SReal`, default: 1.0)\n- **showTranslation**: Position of the first cell for visualization (`type::Vec3`, default: (0, 0, 0))\n\n## Dependencies and Connections\nThe `LCPConstraintSolver` typically requires components that define constraints (e.g., contact or frictional forces) and interacts with other solvers and constraint correction mechanisms in the SOFA scene graph.\n\n## Practical Notes\n- The solver uses iterative methods like Gauss-Seidel for solving LCPs, which can be sensitive to initial conditions and regularization terms.\n- Setting `initial_guess` to true can improve performance but may require careful tuning of other parameters.\n- The friction coefficient (`mu`) should be set appropriately based on the material properties being simulated."
}