Stress recovery from a discrete displacement field#

The FE solve produces a nodal displacement vector \(\mathbf{u}\) of shape \((N,)\) indexed by the DOF map. Stress is not stored in the result file — it’s a derived quantity recovered on demand from \(\mathbf{u}\) plus the model’s element kernels and material table. This page documents how femorph-solver does that recovery, why the recovery is per-element rather than purely nodal, what the multi-element averaging step does and doesn’t preserve, and the failure modes the resulting stress field exhibits when a user reads StaticResult.stress.

The implementation lives in femorph_solver.result._stress_recovery and the public femorph_solver.recover.compute_nodal_stress / femorph_solver.recover.compute_nodal_strain helpers expose it as a free-function path.

Per-element strain#

For an isoparametric element with shape functions \(N_i(\xi)\) and nodal displacements \(\mathbf{u}_e\), the strain at a parametric coordinate \(\xi\) is

\[\boldsymbol{\varepsilon}(\xi) \;=\; \mathbf{B}(\xi)\, \mathbf{u}_e,\]

with \(\mathbf{B}(\xi) = \partial \mathbf{N}/\partial \mathbf{x}\) the 6 × (n_dof_per_node × n_nodes) strain-displacement matrix. The pull-back from parametric to physical coordinates uses the Jacobian \(\mathbf{J}(\xi)\): \(\partial \mathbf{N}/\partial \mathbf{x} = \mathbf{J}^{-\top}\, \partial \mathbf{N}/\partial \boldsymbol{\xi}\).

Strain at the element nodes (rather than at Gauss points) is the canonical recovery target, for two reasons:

  • Engineering quantities of interest — peak von Mises, principal stresses, fatigue scalars — are evaluated at the nodes a user identifies on the mesh.

  • Multi-element averaging at shared nodes (next section) requires the same parametric coordinate to be reachable from every element that touches the node, and “the element’s own nodes” is the natural choice.

The element kernel’s strain_batch() method evaluates \(\mathbf{B}(\xi_\text{node})\, \mathbf{u}_e\) at every node of the element and returns a \((n_\text{nodes}, 6)\) Voigt-format strain array.

Constitutive contraction:

\[\boldsymbol{\sigma}(\xi) \;=\; \mathbf{D}\, \boldsymbol{\varepsilon}(\xi) \;-\; \mathbf{D}\, \boldsymbol{\varepsilon}_\text{th}(\xi),\]

with \(\mathbf{D}\) the linear constitutive matrix for the element’s material id and \(\boldsymbol{\varepsilon}_\text{th}\) the optional thermal-expansion strain (zero for non-thermal problems).

Multi-element averaging#

A node shared by \(k\) elements has \(k\) per-element stress samples — one per touching element — that may disagree because the underlying FE basis is \(C^0\) (continuous in displacement) but not \(C^1\) (continuous in strain). The discrepancy at the node is the first-order error indicator the FE result carries about itself; the averaging step tries to recover the smooth field underneath.

femorph-solver uses the simple-average smoothing — every element that touches a node contributes its stress sample with weight one:

\[\bar{\boldsymbol{\sigma}}_\text{node} \;=\; \frac{1}{k} \sum_{e \in N(\text{node})} \boldsymbol{\sigma}_e(\text{node}),\]

where \(N(\text{node})\) is the set of elements containing the node. This is the same convention MAPDL uses for PowerGraphics plotting, NASTRAN’s MATRIX recovery, and Abaqus’s default fieldOutput.

What averaging preserves:

  • Smoothness across element boundaries when the underlying field is smooth — as the mesh refines, the per-element samples agree at every node and the average converges to the analytic stress.

  • Per-component independence — averaging happens componentwise in Voigt notation, so principal stresses and von Mises are computed on the averaged tensor at each node afterwards.

What averaging does *not* preserve:

  • The exact-equilibrium property of the per-element stress. After averaging, \(\nabla \cdot \bar{\boldsymbol{\sigma}} \neq \mathbf{0}\) in general — even on a problem where the per-element stress is exactly equilibrating. This is the classical “stress-recovery error” that drives every Z²-style estimator.

  • Discontinuity at material interfaces. When two elements on opposite sides of a material boundary share a node, the averaging mixes their stresses — which is wrong physically because the actual stress state at a bonded interface is jump-discontinuous in the strain components transverse to the interface. Use StaticResult.elastic_strain_per_element for the un-averaged dictionary if you need to inspect the per-element values around such interfaces.

The averaging is therefore a smoothing — it suppresses the \(O(h^p)\) error indicator into a \(O(h^{p+1})\) super-convergent estimate of the smooth field, at the cost of hiding the indicator that told you whether the mesh was good enough.

Z²-style super-convergent estimate#

The Zienkiewicz–Zhu (ZZ, often Z²) super-convergent patch recovery (Zienkiewicz & Zhu 1992, Int. J. Numer. Methods Eng. 33 (7), 1331–1364) is the gold-standard for nodal stress recovery on linear-elastic problems. It fits a higher-order polynomial through the Gauss-point stress samples of the elements adjacent to a given node (the patch) and evaluates the polynomial at the node — the result converges one order faster than naive averaging on smooth solutions.

femorph-solver’s current shipping recovery uses the simple average rather than full ZZ. The rationale:

  • The simple average is exactly the convention every commercial verification suite (NAFEMS LE-1 through LE-10, MAPDL VM corpus, NASTRAN VG corpus) compares against, so stress agreement to the published reference is interpreted consistently.

  • For most engineering meshes the difference between simple- average and ZZ recovery is below the engineering judgment threshold (a few percent at the peak).

  • ZZ recovery requires the patch geometry of every node, which means another pass over the connectivity that the current C++ kernel doesn’t accelerate.

A ZZ upgrade is on the roadmap for the case where it actually matters: post-processing for fatigue / fracture analyses where the stress gradient drives the answer.

Shell-fibre stress recovery#

Shell elements (QUAD4_SHELL / SHELL281) carry membrane plus bending strain at the mid-surface, not 6-component Voigt stress. The recovery path is two-step:

  1. Recover the mid-surface in-plane membrane strain \(\boldsymbol{\varepsilon}^\text{mem}\) and the bending curvature \(\boldsymbol{\kappa}\) from the element kernel.

  2. Evaluate the through-thickness stress at the requested fibre position \(z \in [-t/2, +t/2]\) (top fibre, mid fibre, bottom fibre, or any user-supplied \(z\)):

    \[\boldsymbol{\sigma}(z) \;=\; \mathbf{D}\, \bigl( \boldsymbol{\varepsilon}^\text{mem} + z\, \boldsymbol{\kappa} \bigr).\]

The two surfaces (top and bottom) typically carry equal-and- opposite bending contributions; the peak stress on a bending- dominated plate sits at one of the surfaces, not at the mid-surface. The StaticResult.stress() API forwards a fibre keyword to the QUAD4_SHELL kernel for this; see the element page for the surface-position convention.

Per-element strain access (no averaging)#

For diagnostics — comparing per-element stress against the average to see how rough the mesh is, or inspecting the stress discontinuity at a material interface — StaticResult.elastic_strain_per_element() returns the raw {elem_num: (n_nodes_in_elem, 6)} dict that StaticResult.elastic_strain() averages. Same data, no smoothing.

This is the canonical input for any user-side stress estimator (Z² patch recovery, energy-norm error, p-refinement decisions). The library provides the per-element samples; the user keeps control of the smoothing choice.

Failure modes#

Three patterns to recognise in a recovered stress field:

  1. Sharp peaks at concentrated-load nodes. A node where Model.apply_force put a point load carries the singular reaction of the FE approximation. The stress there is not a real engineering stress; the local field needs Saint-Venant smoothing or a distributed-load representation. See Saint-Venant.

  2. Re-entrant-corner singularities. Stress at a re-entrant corner of a 2D plane mesh diverges as \(r^{-\lambda}\) where \(\lambda\) depends on the corner angle (Williams 1952). No FE basis can represent the singularity exactly, so the recovered stress at the corner node is the FE basis’s best polynomial fit to an unbounded field — meaningless as an engineering quantity, but a sharp diagnostic that the mesh needs h-refinement (or a fracture-mechanics treatment if the stress gradient is what matters).

  3. Material-interface discontinuity smoothed away. As noted above, averaging across a bonded material interface mixes stresses across the boundary. When that boundary is load-bearing, switch to per-element strain inspection.

References#

  • R. D. Cook, M. E. Plesha, M. S. Witt, Concepts and Applications of Finite Element Analysis, 4th ed., Wiley 2002, §6.14 — averaged nodal stress recovery.

  • O. C. Zienkiewicz & J. Z. Zhu, “The superconvergent patch recovery and a posteriori error estimates. Part 1: The recovery technique,” International Journal for Numerical Methods in Engineering 33 (7), 1992, 1331–1364 — the Z² patch-recovery algorithm.

  • O. C. Zienkiewicz & J. Z. Zhu, “The superconvergent patch recovery and a posteriori error estimates. Part 2: Error estimates and adaptivity,” International Journal for Numerical Methods in Engineering 33 (7), 1992, 1365–1382.

  • M. L. Williams, “Stress singularities resulting from various boundary conditions in angular corners of plates in extension,” Journal of Applied Mechanics 19, 1952, 526–528 — re-entrant-corner singularity exponents.

See also#