General concepts » Solid mechanics

Finite-strain solid mechanics: kinematics, constitutive laws, and Newton solution.

Introduction

The Solid module implements finite-strain (hyperelastic) solid mechanics on top of the variational form language. It is organized exactly like the underlying mathematics:

Layer	Directory	Contents
Kinematics	`Solid/Kinematics/`	Deformation measures and invariants at a point
Constitutive laws	`Solid/Constitutive/`	Strain-energy densities and their derivatives
Integrators	`Solid/Integrators/`	The internal virtual work residual and tangent
Point plumbing	`Solid/Local/`	Constitutive-point inputs (state, fibers, activation)
Postprocessing	`Solid/Fields/`	Stress and strain fields for output
Linear theory	`Solid/Linear/`	Small-strain elasticity integrators

Working examples: BlockGravity.cpp (quasi-static NeoHookean block under ramped gravity), CantileverBeam.cpp, and ActiveContractionPlaneWave.cpp (active fiber contraction driven by an activation wave).

The residual/tangent sign convention and finite-difference consistency checks used by Solid are part of Rodin's shared Numerical contracts.

Kinematics

With displacement , the deformation gradient is , the right Cauchy–Green tensor , and the Green–Lagrange strain . Isotropic response is expressed through the invariants , , ; a fiber direction adds the anisotropic invariant (squared fiber stretch).

These live in Solid/Kinematics/ (KinematicState, Invariants) and are computed at constitutive points — laws never see meshes or quadrature, only local state.

Constitutive laws

A hyperelastic law is a strain-energy density ; stress and tangent follow by differentiation: , . All laws derive from HyperElasticLaw and return, at a constitutive point, the stress and the consistent tangent.

Law	Energy	Regime / notes
`Hooke`	quadratic in	Linear elasticity
`SaintVenantKirchhoff`	quadratic in	Large rotations, small strains; loses ellipticity in strong compression (a property of the model, not a bug)
`NeoHookean`	-based, compressible	General-purpose rubber-like
`MooneyRivlin`	-based	Rubber with second-invariant sensitivity
`HolzapfelOgden`	isotropic matrix + exponential fiber term in	Anisotropic soft tissue, one fiber family
`ActiveFiberLaw`	1D fiber element with internal state	Active (muscle) contraction
`ActiveContraction<Passive, Active>`	passive law + active fiber law	Composition wrapper

Constructing a law and the internal virtual work form (from BlockGravity.cpp):

Solid::NeoHookean law(lambda, mu);
auto ivw = Solid::InternalVirtualWork(law, u);   // u: current displacement GridFunction

The weak form and Newton's method

Static equilibrium in weak form is

and Newton's method solves, at each iterate, . The tangent contains both the material stiffness (from ) and the geometric (initial-stress) stiffness; omitting the latter forfeits quadratic convergence. InternalVirtualWork packages residual and tangent together — the expression ivw(du, v) contributes both to the Newton problem:

// Newton linearization:  K δu = -F_int(u) + F_ext
Problem newton(du, v);
newton = ivw(du, v)
       - Integral(bodyForce, v)
       + DirichletBC(du, zero).on(bottomBC);

SparseLU linearSolver(newton);
NewtonSolver solver(linearSolver);
solver.setMaxIterations(50)
      .setAbsoluteTolerance(1e-10)
      .setRelativeTolerance(1e-8);
solver.solve(u);      // iterates du-solves, accumulating into u

Note the structure: the increment is the trial function; the Dirichlet condition on is homogeneous once satisfies the constraint; loads are ramped in steps (incremental loading) because a full load applied to the undeformed state may lie outside Newton's basin of attraction.

New solid terms should be documented and tested as residual/tangent pairs; the tangent must be the derivative of the residual assembled into the Newton problem.

Constitutive-point inputs

Laws that need more than kinematics (time step, activation, previous internal state, fiber directions) receive them through the constitutive point* mechanism (Solid/Local/ConstitutivePoint.h): the integrators stamp tags — cell index, quadrature-point index, and any user-relevant quantities — onto the point, and a user-supplied input callable maps tags to law inputs:

auto ivw = Solid::InternalVirtualWork(law, u).setInput(activeInput);

Direct constitutive evaluation binds inputs once, then reads the prepared law through const methods:

Solid::KinematicState state(3);
state.setDisplacementGradient(gradient);
Solid::ConstitutivePoint point(state);
Solid::NeoHookean law(lambda, mu);
law.setConstitutivePoint(point);
Math::SpatialMatrix<Real> stress, increment;
law.getFirstPiolaKirchhoffStress(stress);
law.getMaterialTangent(increment, deformationIncrement);
const Real energy = law.getStrainEnergyDensity();

The law owns its kinematic snapshot and keeps intermediate data private. Changing a point or its borrowed state affects evaluation only after the next non-const setConstitutivePoint() call. Binding compares the consumed kinematics and auxiliary inputs exactly and reuses unchanged data; no address or cache identity determines validity. Spatially graded laws sample their coefficient on each binding. Copies have independent prepared state, and assembly uses separate working laws per worker. The same law must not be rebound concurrently without external synchronization; const reads of an already prepared law perform no cache updates.

An output supplied to InternalVirtualWork::setOutput() receives the constitutive point and a const prepared law, from which it can read stresses, energy or active residuals. Both references are valid for the callback; retain a law copy if its prepared response is needed afterwards.

This keeps laws pure functions of local state (unit-testable against finite differences at a single point) and keeps integrators ignorant of constitutive details. Solid/Local/FiberKinematics.h carries preferred directions; fibers are material (reference-configuration) vectors — their pushforward happens inside laws and fields, never in user setup.

Internal variables

Active laws evolve internal state (active extension, cross-bridge stiffness-like variables) by local ODEs. Two facts worth knowing:

The preferred architecture for internal variables in Rodin is to make them first-class fields (grid functions on discontinuous spaces, coupled into the global Newton residual) rather than hiding local solves and Schur-condensed tangents inside the law evaluation.
Quasi-static problems can latch: after activation ceases, a contracted configuration can itself be an equilibrium, so the model never relaxes (no sliding, no decay path). The remedies are a spontaneous-decay regularization parameter, or a genuinely dynamic formulation (inertia + Newmark-type integration) in which the elastic restoring force drives relaxation. If a quasi-static simulation "won't relax", check the formulation regime before suspecting the law.

Postprocessing fields

Solid/Fields/ provides stress/strain functions for output and coupling: FirstPiolaKirchhoffStress ( ), CauchyStress ( ), GreenLagrangeStrain ( ). They are ordinary form-language functions: project them onto an output space and write them with IO::XDMF.

A note on units

The laws are unit-agnostic; consistency is the user's contract. Mixing a Pa-scale modulus with kPa-scale loads or activation is a silent modeling error no assertion can catch — fix a unit system per example and state it in a comment near the parameters.