Problem formulation

SOFAx v2 targets nonlinear constrained mechanics formulated from Newton's second law. The system is subject to exact constraints that define kinematic relations and exchange power with the mechanical system.

This formulation treats constraints as physical interactions rather than algebraic equations. After discretization and linearization, a coupled primal–dual structure emerges naturally, which guides both the numerical solver design and potential AI integration strategies.

Why this formulation?

The choice of exact constraint enforcement via Lagrange multipliers, rather than penalty methods, provides several advantages:

Physical accuracy: Constraints represent ideal mechanical interactions that exchange power without storage or dissipation, which is naturally captured by Lagrange multipliers
Numerical robustness: Exact enforcement avoids the conditioning issues and parameter tuning required by penalty methods
Solver structure: The coupled primal–dual system enables geometric decompositions and Schur complement strategies that respect the underlying physics
Differentiability: The formulation is fully differentiable, enabling inverse problems and end-to-end learning

As we will see in Constraints & Physical Interactions, this formulation leads to a fundamental geometric decomposition of the configuration space that guides both preconditioning strategies and AI integration points.

Dynamic formulation

The problem is derived from Newton's second law. Consider a mechanical system with generalized coordinates \(q \in \mathbb{R}^n\) and velocities \(v = \dot{q} \in \mathbb{R}^n\). The dynamic behavior is modeled as:

\[ M(q)\,\dot{v} = P(t) - F(q, v) + G(q)^\top \lambda \]

where:

\(M(q) : \mathbb{R}^n \to \mathbb{R}^{n \times n}\) is the inertia matrix
\(v = \dot{q}\) is the velocity vector
\(F(q, v)\) represents internal forces applied to the system depending on the current state
\(P(t)\) gathers known external forces
\(G(q)^\top\) is the matrix containing constraint directions
\(\lambda \in \mathbb{R}^{n_c}\) is the vector of Lagrange multipliers containing constraint force intensities

The constraints are enforced through:

\[ \Phi(q) = 0 \]

where \(\Phi(q) \in \mathbb{R}^{n_c}\) are the constraint functions and \(G(q) = \partial \Phi / \partial q\) is the constraint Jacobian.

Key insight: The constraint forces \(G(q)^\top \lambda\) appear naturally in Newton's second law, reflecting that constraints are physical interactions that exchange power with the system. The Lagrange multipliers \(\lambda\) are not merely algebraic variables but represent the intensity of constraint reactions.

Discrete formulation

After discretization (typically via finite elements) and choice of coordinates, we seek the displacement field \(u(t)\) and the associated multipliers \(\lambda(t)\) such that:

\[ \begin{cases} M(u)\,\ddot u \;+\; B(u)\,\dot u \;+\; f_{\mathrm{int}}(u) \;+\; G(u)^{\top}\,\lambda \;=\; f_{\mathrm{ext}}(t) \\[6pt] \Phi(u) = 0 \end{cases} \]

where \(G(u) = \partial \Phi / \partial u\) is the constraint Jacobian.

This corresponds to constrained elastodynamics with exact constraint enforcement at the continuous level, where constraints are understood as physical interactions that exchange power with the system.

The discrete formulation preserves the coupled structure: the dynamics equation depends on the multipliers \(\lambda\), while the constraint equation \(\Phi(u) = 0\) must be satisfied exactly. This coupling is fundamental and cannot be eliminated without losing the physical interpretation of constraints.

Notation

\(u(t)\): displacement degrees of freedom
\(\dot u(t), \ddot u(t)\): velocity and acceleration
\(M(u)\): mass operator (may depend on configuration for geometric nonlinearity)
\(B(u)\): damping operator
\(f_{\mathrm{int}}(u) = -\partial W(u)/\partial u\): internal forces derived from the strain energy \(W\)
\(f_{\mathrm{ext}}(t)\): external forces
\(\Phi(u) \in \mathbb{R}^{n_c}\): constraint functions (holonomic or non-holonomic)
\(G(u) = \partial \Phi / \partial u\): constraint Jacobian
\(\lambda \in \mathbb{R}^{n_c}\): Lagrange multipliers (constraint force intensities)

Note: In practice, the unknowns solved for at each time step are typically accelerations \(a = \ddot u\) and multipliers \(\lambda\), with velocities and displacements updated via the time integration scheme. See Time Integration for details.

Coupled primal–dual structure

After linearization (as occurs in Newton iterations), the system takes the form:

\[ \begin{bmatrix} A & G^\top \\ G & 0 \end{bmatrix} \begin{bmatrix} \delta x \\ \delta \lambda \end{bmatrix} = - \begin{bmatrix} r_u \\ r_c \end{bmatrix} \]

where:

\(x\) denotes primal unknowns (typically accelerations \(a\))
\(A\) is the tangent operator gathering inertia, damping, and stiffness effects
\(G = \partial \Phi / \partial x\) is the constraint Jacobian
\(r_u\) is the physics residual
\(r_c = \Phi(u)\) is the constraint violation

This block structure is the discrete representation of a primal–dual coupling arising from Newton's second law with constraints, not merely an algebraic convenience. The geometric interpretation of this structure is discussed in Constraints & Physical Interactions.

Implementation principles

The formulation leads to several key implementation principles:

Residual-based problem: Each time step (or equilibrium solve) reduces to finding a root of a nonlinear function. The solver evaluates residuals and applies their Jacobians via automatic differentiation, never assembling global matrices. See Residual Operator.
Exact constraints: Constraints remain exact (no penalty-only formulation) and enter the solve through the coupled variables \((x,\lambda)\). This preserves physical accuracy and enables geometric decompositions.
Projective boundary conditions: Boundary conditions are enforced through projective operators applied to primal variables only. Multipliers remain unconstrained, consistent with the Lagrange multiplier formulation.
Differentiability: The formulation is compatible with implicit time integration, Newton–Krylov solvers, and full differentiation of the simulation pipeline. This enables inverse problems and end-to-end learning, as discussed in AI for Simulation.
Matrix-free operations: All linear algebra is performed via operator actions (JVP/VJP), enabling scalability to large systems and seamless integration with automatic differentiation frameworks like JAX.

Connection to solver structure

This problem formulation directly guides the solver design:

Time integration (Time Integration): Each step solves for \((a_{n+1}, \lambda_{n+1})\) via a nonlinear root-finding problem
Newton–Krylov solver (Newton–Krylov Solver): The coupled system is solved using Newton iterations, with each linearized system solved via Krylov methods (matrix-free)
Preconditioning (Preconditioning): The primal–dual structure enables block preconditioners and Schur complement strategies that respect the geometric decomposition
AI integration (AI for Simulation): The operator-based, matrix-free structure provides natural integration points for learned components