Plane strain / plane stress wrappers

PlaneStrain and PlaneStress adapt a 3D material model to a 2D plane analysis:

material = PlaneStrain(NeoHooke(1.0e5, 0.3))   # out-of-plane stretch prescribed
material = PlaneStress(NeoHooke(1.0e5, 0.3))   # out-of-plane stress solved to zero

Both wrappers leave the constitutive evaluation to the wrapped 3D model and return only its in-plane stress and tangent. A small-strain model stays geometrically linear under both wrappers and remains its preassembled stiffness.

The strain measure of a wrapped model determines which out-of-plane quantity is prescribed or is solved for:

For a usage example, see the 2D plate with a hole.

Any 3D material built on material_response can be wrapped natively. Hooking into the wrappers describes the optional stress-only hook and lists bundled exceptions.

The following sections describe the algorithms and maths behind the finite-strain path and the small-strain reduction. The code is provided in src/PlaneStrainStress.jl.

Notation

A 2D problem lives in the x₁–x₂ plane. The 3D embedding adds an x₃-direction that is out-of-plane. Vectors and tensors follow the conventions of Tensors.jl. In-plane quantities use the symbols F (deformation gradient), P (first Piola–Kirchhoff stress), and σ (Cauchy stress); their 3D counterparts carry an overbar: F̄, P̄, etc.

Embedding the 2D deformation gradient

Both wrappers start by embedding the 2D F into a 3D F̄ (using embed_F_2D_to_3D(F, F33)). The choice of F̄₃₃ distinguishes plane strain from plane stress.

Plane strain with F̄₃₃ = 1

The constraint F̄₃₃ = 1 in

\[\bar F \;=\; \begin{pmatrix} F_{11} & F_{12} & 0 \\ F_{21} & F_{22} & 0 \\ 0 & 0 & 1 \end{pmatrix}\]

enforces ε₃₃ = 0 in the reference configuration.

Plane stress with unknown F̄₃₃

Plane stress enforces σ̄₃₃ = 0. Due to the absence of out-of-plane shear (F̄₁₃ = F̄₂₃ = F̄₃₁ = F̄₃₂ = 0) in 2D, this constraint is equivalent to P̄₃₃ = 0 (since P̄ = J σ̄ F̄⁻ᵀ and thus P̄₃₃ = J σ̄₃₃ (F̄⁻ᵀ)₃₃ = J σ̄₃₃ / F̄₃₃). The 3D constitutive update becomes an implicit equation in the scalar F̄₃₃ by defining a scalar residual function g as

\[g(F̄_{33}) \;:=\; \bar P_{33}(F_{11}, F_{12}, F_{21}, F_{22}, F̄_{33}) \;=\; 0 \text{.}\]

A local Newton iteration solves g(F̄₃₃) = 0 at each quadrature point per assembly call, using the previously converged value (see State handling) as the initial guess.

PlaneStrain: PK1, tangent, state

PlaneStrain returns the in-plane PK1 via P = extract_P_2D(P̄). Since F̄₃₃ = 1 is prescribed, the in-plane consistent tangent is the in-plane block of the 3D tangent that material_response returns, so ∂P/∂F[i,j,k,l] = ∂P̄/∂F̄[i,j,k,l] for i,j,k,l ∈ {1,2}.

If a model instead opts into tangent_source= compute_PK1_3D, it calculates the same tangent by differentiating that hook with respect to the 2D F, which propagates F̄₃₃ = 1 through the differentiation:

dP2D_dF2D = Tensors.gradient(F2D_ -> extract_P_2D(Tensor{2,3}(compute_PK1_3D(model, embed_F_2D_to_3D(F2D_, 1.0), dt, state))), F2D)

PlaneStress: local Newton & condensation

PlaneStress wraps a Newton solver around the 3D material. The loop runs inside the wrapper's material_response, i.e., at every quadrature point of every assembly call. The local Newton controls are configurable via keyword arguments tol and maxiter in PlaneStress(model; tol=1e-10, maxiter=20). The core algorithm for this local solve works as follows:

F̄₃₃  ←  state.out33_previous      # initial guess from previous load step

repeat up to maxiter:
    P̄, ∂P̄/∂F̄  ←  material_response(model, embed_F_2D_to_3D(F, F̄₃₃), …)
    if |P̄₃₃| ≤ tol · max(‖P̄‖, 1): break
    F̄₃₃  ←  F̄₃₃  −  P̄₃₃ / ∂P̄₃₃/∂F̄₃₃
state.out33_current ← F̄₃₃

Each iteration therefore requires only one evaluation of material_response, whose tangent contains the Newton derivative ∂P̄₃₃/∂F̄₃₃. A model using the compute_PK1_3D hook instead needs two constitutive evaluations per correction. The full 3D tangent is formed once at the converged F̄₃₃.

If the local Newton solve does not converge within maxiter, or if ∂P̄₃₃/∂F̄₃₃ drops below the threshold sqrt(eps(Float64)) ≈ 1.5e-8, the wrapper throws a PlaneStressConvergenceError (a LocalAssemblyFailure subtype) containing a reason symbol (:small_newton_derivative or :local_newton_nonconvergence). Outer load step controllers can catch this directly or call try_stiffness_matrix, which returns a named tuple with converged = false and the error object.

At the converged F̄₃₃ the wrapper extracts the in-plane PK1 stress and determines the consistent 2D tangent from the full 3D tangent ∂P̄/∂F̄ by static condensation (Schur complement):

\[\frac{\partial P_{ij}}{\partial F_{kl}} \;=\; \frac{\partial \bar P_{ij}}{\partial \bar F_{kl}} \;-\; \frac{\partial \bar P_{ij}}{\partial \bar F_{33}} \frac{\partial \bar P_{33}}{\partial \bar F_{kl}} \;\Big/\; \frac{\partial \bar P_{33}}{\partial \bar F_{33}} \quad \forall\, i,j,k,l \in \{1,2\} \text{,}\]

which is implemented by condense_tangent_2D in src/PlaneStrainStress.jl. Note that if the out-of-plane stiffness component in this full 3D tangent drops below the numerical threshold, this condensation step will independently throw a PlaneStressConvergenceError with the reason :singular_condensation.

State handling

PlaneStress requires an extended state representation: a PlaneStressStateWrapper holds the wrapped 3D state alongside the current and previous out-of-plane component, F̄₃₃ for a finite-strain model and ε̄₃₃ for a small-strain one. The finite-strain solve warm-starts from the previous value; the small-strain solve starts from ε̄₃₃ = 0 and only records the converged value, which compute_stresses then warm-starts from. The 3D state is written only after the local Newton loop converges, and revert_state! rolls back both it and the out-of-plane component when a load step is rejected.

For PlaneStrain, create_state returns the wrapped model's own state and material_response returns the trial state, as for any unwrapped model. A linear PlaneStress material is preassembled and its cells don't hold a quadrature point state, so no PlaneStressStateWrapper is needed.

The small-strain reduction

When the wrapped model declares SmallStrain kinematics, the wrappers drive it with ε and reduce the returned σ and ∂σ/∂ε, keeping the wrapped analysis geometrically linear as in 3D.

The 2D strain is embedded by embed_eps_2D_to_3D(ε, ε33), analogous to embed_F_2D_to_3D:

\[\bar\varepsilon \;=\; \begin{pmatrix} \varepsilon_{11} & \varepsilon_{12} & 0 \\ \varepsilon_{12} & \varepsilon_{22} & 0 \\ 0 & 0 & \bar\varepsilon_{33} \end{pmatrix}\]

The out-of-plane shear components are held at zero, matching F̄₁₃ = F̄₂₃ = F̄₃₁ = F̄₃₂ = 0 on the finite-strain path.

PlaneStrain: restriction

Plane strain prescribes ε̄₃₃ = 0, so the in-plane tangent is the plain restriction of the 3D tangent:

\[D^{2D}_{ijkl} \;=\; D_{ijkl}, \qquad i,j,k,l \in \{1,2\}\]

The restriction stays exact for a nonlinear algorithmic tangent, e.g. for PlaneStrain(J2Plasticity(...)). The recovered σ̄₃₃ is the out-of-plane reaction stress and is dropped from the 2D output.

PlaneStress: scalar Newton on ε̄₃₃

Plane stress makes ε̄₃₃ the unknown and enforces σ̄₃₃ = 0:

ε̄₃₃ ← 0
repeat up to maxiter:
    σ̄, D̄  ←  material_response(model, embed_eps_2D_to_3D(ε, ε̄₃₃), …)
    if |σ̄₃₃| ≤ tol · max(‖σ̄‖, 1): break
    ε̄₃₃  ←  ε̄₃₃ − σ̄₃₃ / D̄₃₃₃₃

The in-plane tangent follows by the same static condensation as the finite-strain path, applied to D̄ instead of ∂P̄/∂F̄. As the residual is linear in ε̄₃₃ for a linear model, the loop converges in one step and the condensed tangent is constant. is_linear is forwarded and the cell is preassembled.

Hooking into the wrappers

Every wrappable model declares kinematics, which selects the reduction path. This includes models that assemble their own element.

A FiniteStrain model may instead offer a stress-only hook for the wrappers to differentiate:

tangent_source(::MyModel) = compute_PK1_3D         # opt in
compute_PK1_3D(material, F̄, dt, state) -> P̄        # 3D first Piola-Kirchhoff stress
update_state_from_3D!(state, material, F̄, dt)      # trial state at the converged F̄

This is purely a performance choice; it pays off for model's whose 3D tangent is expensive to form. Note that an explicit tangent_source declaration is required because compute_PK1_3D always has a generic fallback. Furthermore, an AbstractHyperelastic subtype needs no own methods: the declaration and a compute_PK1_3D derived from Ψ are already provided by the supertype.

compute_PK1_3D is differentiated with ForwardDiff, so the deformation gradient carries ForwardDiff.Dual entries while the committed history stays Float64. A model whose state object or scratch cache mixes both types fails with a MethodError, which try_stiffness_matrix does not catch because it is not a LocalAssemblyFailure. The Developer guide explains how to write state and cache types that work with both.

Bundled model notes

Check the Stable models table for a full list of bundled 3D materials that support the wrapper API. Ogden is a bundled hyperelastic that opts out because AD through its spectral powers is non-finite at repeated principal stretches, where its analytic tangent stays finite.

VEVP_Zhao2021_AT does not support the generic wrapper AD path: its analytic-tangent implementation is Float64/matrix-based. It provides a specialized PlaneStrain(VEVP_Zhao2021_AT(...)) assembly path, but PlaneStress is not implemented for VEVP_Zhao2021_AT. Use VEVP_Zhao2021_AD for plane stress analyses.

For implementing a new 3D model, see the Developer guide.