2. The Forward Finite Element Method

Strong vs weak form, virtual work, minimum potential energy, shape functions, isoparametric mapping, Gauss quadrature, B-matrix, assembly, KU=F, convergence, shear locking, Euler-Bernoulli vs Timoshenko, Kirchhoff vs Reissner-Mindlin, section strains

Contents
1. From Strong Form to Weak Form 2. Virtual Work and Minimum Potential Energy 3. Discretization: Shape Functions 4. Isoparametric Mapping and Gauss Quadrature 5. The B-Matrix and Element Stiffness 6. Assembly, Boundary Conditions, and Solving KU = F 7. Convergence and the Patch Test 8. Beam Kinematics: Euler-Bernoulli vs Timoshenko 9. Plate Kinematics: Kirchhoff vs Reissner-Mindlin 10. Shear Locking: Mechanism and Cures 11. Section Strains: The Bridge to iFEM Interactive: One Element to Global System Interactive: Shear-Locking Laboratory Interactive: From Strong Form to KU = F (Walkthrough) Flashcards

Every inverse method in this guide is built on top of a forward solver's machinery. In the observation-model framing that organizes the whole course — $y = \mathcal{H}(u, p, f, X_0) + b + \eta$ — this module builds the operator that evaluates the physics: given geometry, material, supports, and loads, produce the displacement field and, from it, anything a sensor could read. Module 3 formalizes what happens when you try to run this map backwards; Module 4 shows how iFEM reuses almost every piece of it, verbatim, with the knowns and unknowns swapped.

Core Problem
Given a body $\Omega$ with known elastic properties, prescribed supports, and prescribed loads: compute the displacement field $\mathbf{u}(\mathbf{x})$ — in practice, a finite vector of nodal displacements $\mathbf{U}$ solving $\mathbf{K}\mathbf{U} = \mathbf{F}$ — and understand precisely which parts of the machinery (shape functions, $\mathbf{B}$-operators, section strains) survive unchanged when the problem is inverted.

1. From Strong Form to Weak Form

Start with the simplest structure that has all the ingredients: a 1D elastic bar of length $L$, stiffness $EA$, fixed at $x=0$, with distributed axial load $f(x)$ and an end force $P$ at $x=L$. The strong form states equilibrium pointwise:

$$ (EA\,u')' + f = 0 \;\;\text{on } (0,L), \qquad u(0) = 0, \qquad EA\,u'(L) = P $$

The displacement condition $u(0)=0$ is an essential boundary condition; the force condition $EA\,u'(L)=P$ is a natural one — the names will earn themselves in a moment. The strong form demands that $u$ be twice differentiable and that equilibrium hold at every point. That is a heavy demand, and the standard $C^0$ Lagrange piecewise polynomials — the cheapest functions we can compute with — cannot meet it. (Piecewise polynomials of higher interelement continuity do exist; Section 3 builds $C^1$ Hermite cubics and Section 9 names the $C^1$ triangles. They cost more, which is exactly why the weak form is worth deriving.)

The Four-Step Derivation

  1. Multiply the residual $(EA\,u')' + f$ by an arbitrary test function $w$ with $w(0)=0$ (zero wherever displacement is prescribed).
  2. Integrate the product over the domain: $\int_0^L \left[(EA\,u')' + f\right] w \, dx = 0$.
  3. Integrate by parts once, transferring one derivative from $u$ onto $w$: $\int_0^L (EA\,u')' w \, dx = \left[EA\,u'\,w\right]_0^L - \int_0^L EA\,u'w' \, dx$.
  4. Apply the boundary data. At $x=0$, $w=0$ kills the boundary term. At $x=L$, $EA\,u'(L) = P$, so the term becomes $P\,w(L)$. The traction condition has entered the equation naturally, through the boundary integral. The displacement condition never appears in the equation at all — it must be built into the space of admissible $u$, essentially.
Weak Form — 1D Bar
Find $u$ with $u(0)=0$ such that $$\int_0^L EA\, u' w'\, dx = \int_0^L f\, w\, dx + P\, w(L) \quad \forall\, w,\; w(0)=0$$ Only first derivatives of $u$ and $w$ appear, and equilibrium is enforced in a weighted-average sense rather than pointwise.

The same construction generalizes directly to 3D linear elasticity. The strong form is

$$ -\nabla\cdot\boldsymbol{\sigma} = \mathbf{b} \;\text{in }\Omega, \qquad \mathbf{u}=\bar{\mathbf{u}} \;\text{on }\Gamma_D, \qquad \boldsymbol{\sigma}\mathbf{n}=\bar{\mathbf{t}} \;\text{on }\Gamma_N $$

and the weak form is the virtual-work integral: find $\mathbf{u}$ satisfying $\mathbf{u}=\bar{\mathbf{u}}$ on $\Gamma_D$ such that for all admissible $\mathbf{w}$ (zero on $\Gamma_D$),

$$ \int_\Omega \boldsymbol{\varepsilon}(\mathbf{w})^{\mathsf T} \mathbf{D}\, \boldsymbol{\varepsilon}(\mathbf{u})\, d\Omega = \int_\Omega \mathbf{w}^{\mathsf T}\mathbf{b}\, d\Omega + \int_{\Gamma_N} \mathbf{w}^{\mathsf T}\bar{\mathbf{t}}\, d\Gamma $$

The Three Payoffs

  1. The continuity demand drops. From $C^1$/$C^2$ down to piecewise $C^0$: $u, w \in H^1$ — first derivatives square-integrable is enough. This is exactly what makes piecewise-polynomial approximation legal (Section 3).
  2. Neumann data is absorbed automatically. Tractions enter through the boundary integral; there is no flux matching to enforce at element interfaces.
  3. Symmetry and coercivity. For elasticity the bilinear form $a(\mathbf{u},\mathbf{w}) = \int \boldsymbol{\varepsilon}(\mathbf{w})^{\mathsf T}\mathbf{D}\boldsymbol{\varepsilon}(\mathbf{u})\,d\Omega$ is symmetric and coercive — this is what will make $\mathbf{K}$ symmetric positive definite after constraints (Section 6).
Why weak and strong are equivalent (when the solution is smooth)

Any strong solution satisfies the weak form — the derivation above is just multiplication and integration. The converse needs one extra hypothesis: if a weak solution happens to be $C^2$, reverse the integration by parts to get $-\int_0^L \left[(EA\,u')' + f\right] w\, dx + \left[EA\,u'(L) - P\right] w(L) = 0$ for all admissible $w$ (mind the relative sign — the domain and boundary terms come out with opposite signs). By the fundamental lemma of the calculus of variations, an integrand that integrates to zero against every test function must vanish pointwise — so the PDE and the natural boundary condition are both recovered. The weak form loses nothing; it only relaxes where the derivatives are allowed to live. This is Hughes ch. 1 and Zienkiewicz–Taylor–Zhu ch. 3 material.

Strong form −∇·σ = b in Ω u = ū on ΓD, σn = t̄ on ΓN pointwise equilibrium, C² demanded Weak form ∫ ε(w):Dε(u) = ∫ w·b + ∫ w·t̄ C⁰ suffices, natural BCs absorbed Discrete weak form a(uₖ, wₖ) = ℓ(wₖ) best approximation in energy norm KU = F sparse, symmetric, pos. definite after BCs multiply by w, integrate by parts restrict to finite basis uₖ = ∑ Nₖ dₖ (Galerkin) integrate element-wise, assemble GUARDIAN CONCEPT: well-posedness coercivity completeness + compatibility conditioning
The four-stage pipeline of this module. Each stage has a guardian concept that can fail: an ill-posed continuous problem, a non-coercive form, an incomplete or incompatible element, or a badly conditioned system.

2. Virtual Work and Minimum Potential Energy

The weak form has two classical variational faces. They produce the same equations for linear elasticity, but they are not the same statement — and the difference matters for everything nonlinear that comes later.

Principle of Virtual Work (PVW)
A stress field is in equilibrium if and only if, for every kinematically admissible virtual displacement $\delta\mathbf{u}$, internal virtual work equals external virtual work: $$\int_\Omega \delta\boldsymbol{\varepsilon}^{\mathsf T}\boldsymbol{\sigma}\, d\Omega = \int_\Omega \delta\mathbf{u}^{\mathsf T}\mathbf{b}\, d\Omega + \int_{\Gamma_N} \delta\mathbf{u}^{\mathsf T}\bar{\mathbf{t}}\, d\Gamma$$
Key insight
PVW is a statement about equilibrium, not material law — no constitutive assumption appears in it. It holds for nonlinear and inelastic materials just as well, which is why it is the more general starting point for FE formulations (Bathe ch. 4). Minimum potential energy, below, is the narrower statement that needs an energy to exist.

Aside: Minimum Potential Energy

When a strain energy density exists — linear elasticity gives $U = \tfrac{1}{2}\boldsymbol{\varepsilon}^{\mathsf T}\mathbf{D}\boldsymbol{\varepsilon}$ — and the loading is conservative (a load potential exists), define the total potential energy:

$$ \Pi(\mathbf{u}) = \tfrac{1}{2}\int_\Omega \boldsymbol{\varepsilon}^{\mathsf T}\mathbf{D}\,\boldsymbol{\varepsilon}\, d\Omega - \int_\Omega \mathbf{u}^{\mathsf T}\mathbf{b}\, d\Omega - \int_{\Gamma_N} \mathbf{u}^{\mathsf T}\bar{\mathbf{t}}\, d\Gamma, \qquad \delta\Pi = 0 $$

The first variation $\delta\Pi = 0$ recovers exactly the weak form. For linear elasticity with positive-definite $\mathbf{D}$ and enough essential constraints to remove the rigid-body modes, the second variation is positive definite (coercivity), so the stationary point is the unique minimizer over kinematically admissible fields. Both qualifiers earn their keep: without sufficient essential BCs the form is only semicoercive (Section 6 — $\mathbf{K}$ is singular until you constrain it), and for a general hyperelastic energy, which need not be convex, $\delta\Pi = 0$ delivers stationarity alone — buckling and limit points are precisely the case of several stationary points, and stability has to be established separately. The Ritz method minimizes $\Pi$ over a finite-dimensional span of trial functions. FEM is Ritz with locally supported piecewise polynomials — and local support is precisely what makes the resulting equations sparse. Treat MPE as the special case it is: it presupposes conservative loading and a hyperelastic material, while PVW does not.

Galerkin Orthogonality and the Best-Approximation Property

Because $a(\cdot,\cdot)$ is symmetric and the discrete space $\mathcal{V}_h$ is a subspace of the admissible space, subtracting the discrete weak form from the continuous one gives

$$ a(\mathbf{u}-\mathbf{u}_h,\, \mathbf{w}_h) = 0 \quad \forall\, \mathbf{w}_h \in \mathcal{V}_h $$

The error is orthogonal to the whole discrete space in the energy inner product: $\mathbf{u}_h$ is the projection of the exact solution onto $\mathcal{V}_h$, hence the best approximation in the energy norm that the mesh can express.

Why it matters — "too stiff" is an energy statement
A conforming displacement model with consistent loads underestimates the strain energy of the exact solution: the mesh is "too stiff" in the energy sense. Do not read this pointwise — local displacements can overshoot; the theorem constrains the energy norm, and it requires conforming elements, consistent loads, and homogeneous essential boundary conditions. (Prescribe a nonzero displacement instead and the discrete solution is an affine, not linear, perturbation: Galerkin orthogonality no longer delivers the bound and the inequality can reverse.) The same property explains why shear locking (Section 10) is so dramatic: the locked answer is not a solver failure but the honest, energy-optimal solution within a badly over-constrained subspace.

3. Discretization: Shape Functions

Partition $\Omega$ into elements and approximate the displacement field by nodal interpolation:

$$ \mathbf{u}_h(\mathbf{x}) = \sum_a N_a(\mathbf{x})\, \mathbf{d}_a, \qquad N_a(\mathbf{x}_b)=\delta_{ab}, \qquad \sum_a N_a = 1 $$

Four properties make a set of shape functions usable:

PropertyStatementWhat it buys
Kronecker delta$N_a(\mathbf{x}_b) = \delta_{ab}$$\mathbf{d}_a$ are physical nodal displacements — BCs and sensors attach to real quantities
Partition of unity$\sum_a N_a = 1$rigid translation represented exactly
Completenessexact to polynomial degree $p$linear completeness ⇒ constant strain representable — the ingredient the patch test checks (Section 7)
Compact support$N_a \ne 0$ only on elements touching node $a$sparse $\mathbf{K}$

The two workhorses of this course: the 1D two-node bar on $\xi \in [-1,1]$,

$$ N_1 = \tfrac{1-\xi}{2}, \qquad N_2 = \tfrac{1+\xi}{2} $$

and the 4-node bilinear quadrilateral,

$$ N_a(\xi,\eta) = \tfrac{1}{4}\left(1+\xi_a\xi\right)\left(1+\eta_a\eta\right), \quad \xi_a,\eta_a = \pm 1 $$

How Much Continuity Do You Need?

For a $2m$-th order problem the weak form contains $m$-th derivatives, so trial functions need $C^{m-1}$ interelement continuity and $H^m$ integrability. Second-order elasticity: $m=1$, so $C^0$ suffices — hat functions are legal. Fourth-order Euler-Bernoulli beams and Kirchhoff plates: $m=2$, so $C^1$ is required. In 1D that is achievable with the two-node Hermite cubic element (4 DOFs: $w_1, \theta_1, w_2, \theta_2$; shape functions $H_1 = 1 - 3s^2 + 2s^3$ and companions). In 2D, conforming $C^1$ triangles need high-order polynomials — a historical pain point that motivated the $C^0$ shear-deformable formulations of Sections 8–9, which are exactly the formulations iFEM later builds on.

Key insight
Within one element the interpolation is written in matrix form, $\mathbf{u}_h = \mathbf{N}\,\mathbf{d}^e$ — for a plane quad, $\mathbf{N}$ is the $2\times 8$ shape-function matrix. This matrix is the single object that everything downstream differentiates: the $\mathbf{B}$-matrix of Section 5 is nothing but derivatives of $\mathbf{N}$, and the iFEM strain operators of Module 4 differentiate the same $\mathbf{N}$.

4. Isoparametric Mapping and Gauss Quadrature

Real meshes are not unions of squares. The isoparametric idea — use the same shape functions for the geometry as for the displacement — is what lets arbitrarily shaped (and curved) elements exist at all (Zienkiewicz–Taylor–Zhu ch. 6; the idea entered the literature with Irons and the Swansea school in the 1960s):

$$ \mathbf{x}(\boldsymbol{\xi}) = \sum_a N_a(\boldsymbol{\xi})\,\mathbf{x}_a, \qquad \mathbf{J} = \frac{\partial \mathbf{x}}{\partial \boldsymbol{\xi}}, \qquad \nabla_{\mathbf{x}} N_a = \mathbf{J}^{-\mathsf T}\, \nabla_{\boldsymbol{\xi}} N_a $$

The Jacobian $\mathbf{J}$ ($2\times 2$ in 2D) converts parent-domain derivatives into physical derivatives, and the volume element transforms as $d\Omega = \det\mathbf{J}\, d\xi\, d\eta$. Note the trap this kills: shape-function derivatives $\nabla_\xi N_a$ live in the parent domain — strains need physical derivatives, and only $\mathbf{J}^{-\mathsf T}$ gets you there.

Why it matters — distorted and inverted elements
The map is valid only where $\det\mathbf{J} > 0$ everywhere in the element: interior angles below 180°, no collapsed or inverted nodes. And validity is not a binary: a $\det\mathbf{J}$ near zero relative to the element's own scale destroys accuracy before it destroys invertibility. The qualifier is load-bearing, because $\det\mathbf{J}$ is dimensional — it scales with element area, so a small but perfectly shaped element has a small $\det\mathbf{J}$ and nothing wrong with it. The scale-invariant measure is the scaled Jacobian, $\det\mathbf{J}$ divided by the lengths of the two mapped basis vectors — which is exactly the signed sine of the angle between them: $1$ when the mapped axes are orthogonal, $0$ when they fall parallel, negative when the element is inverted. Read it for what it is, a measure of angular distortion. It is blind to aspect-ratio collapse: a $100{:}1$ sliver scores a flawless $1.000$, exactly like a square, because its mapped axes are still perpendicular. Aspect ratio is a separate quality check, and the lab below will happily print $1.000$ for the thin strip that Section 10 then gets badly wrong. Drag a node past the opposite diagonal in the lab and watch $\det\mathbf{J}$ go negative; watch the scaled Jacobian degrade well before that, because that distortion is angular.

Gauss-Legendre Quadrature

The element integrals are evaluated numerically. $n$ Gauss-Legendre points integrate 1D polynomials exactly through degree $2n-1$:

$n$Points $\xi_i$Weights $w_i$Exact through degree
1$0$$2$1
2$\pm 1/\sqrt{3}$$1,\; 1$3
3$0,\; \pm\sqrt{3/5}$$8/9,\; 5/9$5

Quads and hexes use tensor products of the 1D rule:

$$ \int_{\Omega_e} f\, d\Omega = \int_{-1}^{1}\!\!\int_{-1}^{1} f(\xi,\eta)\, \det\mathbf{J}\, d\xi\, d\eta \;\approx\; \sum_{i=1}^{n}\sum_{j=1}^{n} w_i w_j\, f(\xi_i,\eta_j)\, \det\mathbf{J}(\xi_i,\eta_j) $$

"Full integration" for the bilinear quad means the $2\times 2$ rule. The same quadrature loop also builds consistent load vectors and mass matrices — one loop, three products.

Key insight — "full integration" is a convention, not a guarantee
For a distorted element the stiffness integrand $\mathbf{B}^{\mathsf T}\mathbf{D}\mathbf{B}\,\det\mathbf{J}$ is a rational function of $(\xi,\eta)$, because $\mathbf{B}$ contains $\mathbf{J}^{-1}$. A Gauss-Legendre rule carries no exactness guarantee for a rational integrand — its degree-$2n-1$ guarantee is a statement about polynomials, and nothing replaces it here. "Full" means "exact in the affine (undistorted) limit" — nothing more. Claiming that $2\times 2$ computes the isoparametric stiffness exactly is a classic error.

Reduced integration (one order lower) is cheaper and softens the element — genuinely helpful against locking (Section 10) — but it risks rank deficiency: the element stiffness can be written as a sum of rank-limited contributions, so $\operatorname{rank}(\mathbf{k}^e) \le n_{gp} \times (\text{rows of } \mathbf{B})$. A 1-point plane quad has rank $\le 3$ against the $8 - 3 = 5$ independent deformation modes it must resist — leaving 2 spurious zero-energy (hourglass) modes. You will count them yourself in the lab below.

Top: the isoparametric pipeline — parent square with the $2\times2$ Gauss points at $\pm 1/\sqrt{3}$, mapped through $\mathbf{x}(\boldsymbol{\xi}) = \sum N_a \mathbf{x}_a$ to a distorted physical quad with its $\xi$/$\eta$ isolines and mapped Gauss points. Bottom: each $\mathbf{k}^e$ is scatter-added into the sparse global $\mathbf{K}$ via the connectivity array; shared nodes produce overlapping (darker) blocks.

5. The B-Matrix and Element Stiffness

For plane problems, strains are collected in Voigt notation: $\boldsymbol{\varepsilon} = [\varepsilon_{xx},\, \varepsilon_{yy},\, \gamma_{xy}]^{\mathsf T}$ with engineering shear $\gamma_{xy} = 2\varepsilon_{12}$ — a classic sign-up-for-confusion point carried over from Module 1: mix the tensor and engineering conventions once and every shear term downstream is off by a factor of 2.

The strain-displacement matrix $\mathbf{B}$ maps nodal DOFs to strains, assembled from one $3\times 2$ block per node:

$$ \boldsymbol{\varepsilon} = \mathbf{B}\,\mathbf{d}^{e}, \qquad \mathbf{B}_a = \begin{bmatrix} \partial N_a/\partial x & 0 \\ 0 & \partial N_a/\partial y \\ \partial N_a/\partial y & \partial N_a/\partial x \end{bmatrix} $$

The constitutive matrix for isotropic plane stress:

$$ \mathbf{D}_{\text{plane stress}} = \frac{E}{1-\nu^2}\begin{bmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & (1-\nu)/2 \end{bmatrix} $$

Plane strain swaps the factors: $E/(1-\nu^2) \to E/((1+\nu)(1-2\nu))$, the diagonal $1 \to 1-\nu$, $\nu$ stays, and the shear entry becomes $(1-2\nu)/2$. (Different reduced matrices — never interchange them silently.)

The Element Stiffness Loop

$$ \mathbf{k}^{e} = \int_{\Omega_e} \mathbf{B}^{\mathsf T}\mathbf{D}\,\mathbf{B}\; t\; d\Omega \;\approx\; \sum_{gp} \mathbf{B}^{\mathsf T}\mathbf{D}\,\mathbf{B}\; t\, \det\mathbf{J}\; w_{gp} $$

with $t$ the thickness. In code, this is a loop over Gauss points; at each point:

  1. Evaluate $\nabla_\xi N_a$, build $\mathbf{J}$, invert it.
  2. Form physical derivatives $\nabla_x N_a = \mathbf{J}^{-\mathsf T} \nabla_\xi N_a$ and stack them into $\mathbf{B}$.
  3. Accumulate $\mathbf{B}^{\mathsf T}\mathbf{D}\mathbf{B}\, t \det\mathbf{J}\, w_i w_j$.

Consistent load vectors come from the same loop — and carry the same thickness: $\mathbf{f}^e = \int \mathbf{N}^{\mathsf T}\mathbf{b}\; t\; d\Omega + \int \mathbf{N}^{\mathsf T}\bar{\mathbf{t}}\; t\; d\Gamma$, with $\mathbf{b}$ a force per unit volume and $\bar{\mathbf{t}}$ a force per unit area (the traction $\bar{\mathbf{t}}$ is unrelated to the thickness $t$, despite the letter clash). Keep $t$ in $\mathbf{k}^e$ but drop it here and $\mathbf{K}\mathbf{U}=\mathbf{F}$ is off by a factor of $t$ for any plate that is not unit-thick. Distributed loads become nodal forces through the shape functions, not by ad-hoc lumping.

Why it matters — the eigenvalue test
Three properties of $\mathbf{k}^e$ to check whenever an element implementation misbehaves: it is symmetric; it is positive semi-definite; and with full integration its null space is exactly the rigid-body modes — 3 in 2D (two translations + one rotation), 6 in 3D. Counting the zero eigenvalues of $\mathbf{k}^e$ is the standard element-verification tool (Cook's tradition; presented in Bathe ch. 4). More zeros than rigid modes ⇒ spurious mechanisms; fewer ⇒ a constraint you did not intend. The lab after Section 6 runs this test live.

Everything downstream — including the iFEM machinery of later modules — is built from exactly these $\mathbf{B}$-type operators. iFEM replaces the $\mathbf{D}$-weighted energy with least-squares strain matching, but it differentiates the same shape functions.

6. Assembly, Boundary Conditions, and Solving KU = F

Assembly is scatter-add. A connectivity (LM/location) array maps element-local DOF indices to global ones; each entry of $\mathbf{k}^e$ is added into $\mathbf{K}$ at the mapped row and column. Conceptually,

$$ \mathbf{K} = \sum_e \mathbf{L}_e^{\mathsf T}\, \mathbf{k}^{e}\, \mathbf{L}_e, \qquad \mathbf{F} = \sum_e \mathbf{L}_e^{\mathsf T}\, \mathbf{f}^{e} $$

with Boolean gather matrices $\mathbf{L}_e$ — but this is implemented as loops, never as matrix products. $\mathbf{K}$ is sparse: nonzeros only couple DOFs that share an element. Bandwidth depends on node numbering; storage is CSR or skyline.

Key insight — K is singular until you constrain it
Before boundary conditions, $\mathbf{K}$ has a null space: the rigid-body modes — 3 in 2D, 6 in 3D. "Matrix is singular" is the single most common solver failure, and it almost always means the model can still translate or rotate freely. Related trap: quoting a penalty factor without noting its trade-off (below) — too small leaks constraint error, too large wrecks conditioning.

Three Ways to Impose Essential BCs

MethodMechanismTrade-off
Elimination / partitioningsplit DOFs into free ($f$) and constrained ($c$); solve $\mathbf{K}_{ff}\mathbf{U}_f = \mathbf{F}_f - \mathbf{K}_{fc}\bar{\mathbf{U}}_c$exact, keeps SPD; needs index bookkeeping
Penaltyadd large $\beta$ on the diagonal, $\beta\bar{u}$ to the RHStrivially simple; $\beta$ trades constraint error (too small) against conditioning damage (too large)
Lagrange multipliersappend constraint equations with multiplier unknownsexact, gives reactions directly; produces an indefinite saddle-point system
$$ \begin{bmatrix} \mathbf{K}_{ff} & \mathbf{K}_{fc} \\ \mathbf{K}_{cf} & \mathbf{K}_{cc} \end{bmatrix} \begin{bmatrix} \mathbf{U}_f \\ \bar{\mathbf{U}}_c \end{bmatrix} = \begin{bmatrix} \mathbf{F}_f \\ \mathbf{F}_c + \mathbf{R} \end{bmatrix} \;\Rightarrow\; \mathbf{K}_{ff}\mathbf{U}_f = \mathbf{F}_f - \mathbf{K}_{fc}\bar{\mathbf{U}}_c $$

Solving — and Why You Never Form K⁻¹

Sparse Cholesky (direct) is the default below roughly $10^6$ DOFs; preconditioned conjugate gradients (with algebraic multigrid or incomplete-Cholesky preconditioning) beyond that. For second-order elasticity the condition number grows like $\operatorname{cond}(\mathbf{K}) = O(h^{-2})$ under uniform refinement — the general law for a $2m$-th order problem is $O(h^{-2m})$, so the fourth-order beam and plate discretizations of Sections 8–9 grow like $O(h^{-4})$ — and unpreconditioned CG iteration counts grow with it as the mesh refines. The entire weak-form→$\mathbf{KU=F}$ pipeline of this module is directly executable in a few dozen lines of Python in a modern open-source stack such as DOLFINx (Baratta et al. 2023).

Why it matters — "solve the system", never invert
Write $\mathbf{U} = \mathbf{K}^{-1}\mathbf{F}$ on a whiteboard if you must; never in code. Before constraints the inverse does not exist (rigid modes); after constraints, the explicit inverse of a sparse matrix is dense — forming it throws away the sparsity that assembly worked to create. Factorize or iterate. Every "inverse" method in this course also solves systems; none of them forms an inverse.

Post-Processing and Sanity Checks

Reactions follow from the constrained rows: $\mathbf{R} = \mathbf{K}_{cf}\mathbf{U}_f + \mathbf{K}_{cc}\bar{\mathbf{U}}_c - \mathbf{F}_c$. Stresses are evaluated at the integration points — where $\mathbf{B}$ and $\mathbf{D}$ are already formed — then extrapolated and averaged to nodes for plotting. Accuracy is a separate question: the superconvergent sampling sites are the Gauss points of the reduced rule (the Barlow points of Section 7), which for the bilinear quad this module uses throughout is the element centroid, not the $2\times 2$ points. Two checks worth automating in every solver you ever write: the sum of reactions balances the applied loads, and the energy identity $\tfrac{1}{2}\mathbf{U}^{\mathsf T}\mathbf{K}\mathbf{U} = \tfrac{1}{2}\mathbf{U}^{\mathsf T}\mathbf{F}$ holds. Mind the scope of that second one: it assumes homogeneous essential conditions and configuration-independent loads. With prescribed nonzero displacements the reactions do work, and the identity becomes $\tfrac{1}{2}\mathbf{U}^{\mathsf T}\mathbf{K}\mathbf{U} = \tfrac{1}{2}\mathbf{U}^{\mathsf T}\mathbf{F} + \tfrac{1}{2}\bar{\mathbf{U}}_c^{\mathsf T}\mathbf{R}_c$ — half the reaction work, not all of it: $\mathbf{K}\mathbf{U} = \mathbf{F} + \mathbf{R}$ with $\mathbf{R}$ supported on the constrained DOFs, and the $\tfrac12$ carries through to both sides.

From Solution to Outputs: the Response Operator

The forward solve defines a map from data to displacement, $(\,p, f\,) \mapsto \mathbf{U}$. But no experiment measures $\mathbf{U}$ everywhere. What an instrument reads is an output: a strain gauge reads rows of $\mathbf{B}$ evaluated at its location, a displacement probe reads a few entries of $\mathbf{U}$, a load cell reads reactions. Predicted measurements are an output operator applied to the solution, $\mathbf{y}_{\text{pred}} = \mathbf{Q}\,\mathbf{U} + \mathbf{q}_0$, where the constant offset $\mathbf{q}_0$ vanishes for strain and displacement outputs and is in general nonzero for reaction-type outputs: from the reaction formula above $\mathbf{q}_0 = -\mathbf{F}_c$, which makes a load cell an affine, not linear, function of $\mathbf{U}$ — and which collapses to zero in the common case where no applied load acts on a constrained DOF (as in the lab below, whose loads sit only on free DOFs). This is the bridge to everything that follows: the observation model $y = \mathcal{H}(u, p, f, X_0) + b + \eta$ wraps exactly this forward-solve-plus-output-operator, and Module 3 asks what happens when $y$ is known and one of the arguments is not.

Interactive Tool — One Element to Global System

Drag the blue nodes to distort elements: the picture and the rank/stability panel track the drag live, and $\mathbf{J}$, $\det\mathbf{J}$ and $\mathbf{B}$ are recomputed in the readout below when you release. Switch quadrature, clamp or release the left edge, and read the persistent panel: rank, zero-energy modes, conditioning, and a validation check are recomputed live. The presets walk through the failure modes.

7. Convergence and the Patch Test

When does refining the mesh actually converge to the exact solution? For conforming displacement elements the classical conditions are completeness — the element represents all rigid-body modes and all constant-strain states exactly — plus compatibility — $C^{m-1}$ interelement continuity for a $2m$-th order problem.

A Priori Rates

For smooth solutions and complete degree-$p$ elements on a $2m$-th order problem:

$$ \|\mathbf{u}-\mathbf{u}_h\|_{E(\Omega)} \le C\, h^{p+1-m}\, |\mathbf{u}|_{H^{p+1}(\Omega)}, \qquad \|\mathbf{u}-\mathbf{u}_h\|_{L^2} = O\!\left(h^{\,p+1-m+\min\{m,\; p+1-m\}}\right) $$

The second estimate is the duality (Aubin–Nitsche) bound, and its exponent is written out because the two conditions behind the familiar version are routinely dropped. Duality needs a fully regular dual problem — a smooth or convex domain — and it buys $\min\{m,\, p+1-m\}$ extra orders over the energy rate, not $m$. Only when the element clears $p \ge 2m-1$ does that exponent collapse to the $O(h^{p+1})$ everyone quotes.

The table specializes all of this to second-order elasticity, $m = 1$ — the setting of Sections 1–6 and of every element in the lab above, where $p \ge 1$ clears the bar automatically.

Quantity ($m = 1$)RateWhy
Energy-norm error$O(h^p)$best-approximation + interpolation estimate
Displacement $L^2$ error$O(h^{p+1})$one extra order via duality, given a regular dual problem
Strains / stresses$O(h^p)$differentiation loses one order

Carry those numbers into the $C^1$ theories of Sections 8–9 and each rate is wrong by one order in $h$, not two: Euler-Bernoulli beams and Kirchhoff plates are fourth-order ($m = 2$), so $O(h^{p+1-m})$ takes the energy norm from $O(h^p)$ down to $O(h^{p-1})$ — $O(h^2)$ for the Hermite cubics of Section 3 — and the stress-like quantity, now the curvature, follows it to $O(h^{p-1})$. The $L^2$ displacement rate is the one row that can come through the change of $m$ untouched — but it is a conditional survival, not a free one: it needs $p \ge 2m-1$ on top of the regular dual problem. Hermite cubics ($p=3$, $m=2$) clear that bar exactly, so their $L^2$ rate is still $O(h^{p+1}) = O(h^4)$; an element that falls short of it gets only $2(p+1-m)$, below $p+1$. Keep the two counts apart: changing $m$ costs one order per row, whereas within the beam theory the curvature — a second derivative — trails the displacement rate by two orders, against the single order that separates strain from displacement in elasticity.

Key insight — Barlow points are optimal sampling sites, not magic
Under smooth-solution and regular-mesh hypotheses, strains sampled at the Gauss points of the reduced rule — the Barlow points — can be one order more accurate than elsewhere in the element. Two qualifications, both routinely dropped. First, mind which rule: for the bilinear quad the reduced rule is the single centroid, so the $2\times 2$ integration points are not themselves the Barlow points. That also means the Barlow property is not what explains a code's habit of evaluating stresses at its own integration points — it does that because $\mathbf{B}$, $\mathbf{D}$ and any history variables are already formed there; recovery and smoothing schemes are where Barlow points actually earn their keep. Second, superconvergence is an optimal-sampling result conditional on those hypotheses, not an unconditional property — and the common claim that stresses are "exact at Gauss points" is false.

Singular problems break the rates. At a reentrant corner (the L-shaped domain is the canonical case) the solution behaves like $r^\lambda$ with $\lambda < 1$, and the energy-norm rate degrades to $O(h^{\min(p,\lambda)})$ under uniform refinement regardless of $p$ (the $L^2$ rate degrades too, by a different exponent — always name the norm before quoting a degraded rate). This is the motivation for adaptive $h$-refinement and graded meshes: polynomial degree cannot buy back what a singularity takes away.

The Patch Test

Irons's legacy, with the standard exposition in Zienkiewicz–Taylor–Zhu ch. 8: assemble a small arbitrary patch of distorted elements, prescribe boundary displacements from an arbitrary linear field, solve. The interior must reproduce the linear field and its constant strain to machine precision.

Why it matters — the exact logical status
First, the scope: the guarantee below attaches to the generalized (weak) form of the test, not to the classical constant-strain version described just above. The constant-strain test on its own is neither necessary nor sufficient — it is a practical screen, and elements are known that pass it without converging. With that said: passing a properly formulated patch test is necessary for convergence of nonconforming elements, and — together with stability (no spurious zero-energy modes after assembly) — practically sufficient. Presenting it as unconditionally sufficient is a classic error: a patch-test-passing element with spurious mechanisms still fails. For practitioners: the patch test is the first automated regression test any new element implementation must pass, and the eigenvalue test of Section 5 is the second.

8. Beam Kinematics: Euler-Bernoulli vs Timoshenko

The next three sections deal with structures that are thin in one or two directions — beams and plates. iFEM is built on the shear-deformable members of this family, so the kinematic assumptions deserve care. (Full element derivations — forward and inverse — live in Module 5; here we need the kinematics and the section strains they produce.)

Euler-Bernoulli: the Rotation Is Enslaved

Euler-Bernoulli (EB) assumes plane cross-sections remain plane and perpendicular to the deformed axis. The rotation is then enslaved to the deflection, $\theta = dw/dx$, and the axial displacement field is $u(x,z) = u_0(x) - z\, dw/dx$. Section strains: membrane $e = u_0'$ and curvature $\kappa = d^2w/dx^2$; the transverse shear strain is identically zero.

$$ \text{EB:}\quad \theta = \frac{dw}{dx}, \qquad \kappa = \frac{d^2 w}{dx^2}, \qquad EI\,\frac{d^4 w}{dx^4} = q $$

Fourth order ⇒ $C^1$ elements — the Hermite cubics of Section 3, which have the pleasant property that nodal values are exact for point-loaded prismatic EB beams. (Watch the sign convention: $\kappa = +w''$ here; parts of the beam and plate literatures use $-w''$, and Section 11 shows how an inconsistent sign silently flips iFEM's strain recovery.)

Timoshenko: the Rotation Is Freed

Timoshenko keeps cross-sections plane but not necessarily normal: $\theta(x)$ becomes an independent field, and the gap between the slope and the rotation is the shear strain. Section strains: $e = u_0'$, $\kappa = d\theta/dx$, $\gamma = dw/dx - \theta$. The section constitutive law:

$$ \text{Timoshenko:}\quad \kappa = \frac{d\theta}{dx}, \quad \gamma = \frac{dw}{dx} - \theta, \qquad N = EA\,e, \quad M = EI\,\kappa, \quad Q = kGA\,\gamma $$
Why it matters — k is not a universal constant
The shear correction factor $k$ exists because the theory assumes uniform shear strain through the thickness while the true distribution is nonuniform — parabolic for the elementary rectangular section — and vanishes on traction-free faces. $k = 5/6$ for rectangular sections, $\approx 0.9$ for solid circles — it depends on the cross-section shape and on the matching convention. Those $5/6$ and $0.9$ are static, energy-matched values; matching the thickness-shear cut-off frequency instead gives $\pi^2/12 \approx 0.82$ for the same rectangular section, in beams as in plates (Section 9 shows the same split between Reissner's and Mindlin's plate constants). Quoting "$k = 5/6$" as a universal constant is a trap. The correction originates in Timoshenko's 1921 Philosophical Magazine paper, which added a shear correction to the differential equation for transverse vibrations of prismatic bars.

When Does the Difference Matter?

For a tip-loaded cantilever, $w(L) = \dfrac{PL^3}{3EI} + \dfrac{PL}{kGA}$, so the ratio of the shear contribution to the bending contribution is

$$ \frac{w_{\text{shear}}}{w_{\text{bend}}}\bigg|_{\text{tip-loaded cantilever}} = \frac{3EI}{kGA\,L^2} $$

which for a rectangular section equals $(E/4kG)\,(t/L)^2$ — about 0.8% at $t/L = 1/10$ and 3% at $t/L = 1/5$ in metals ($E/G \approx 2.6$, $k = 5/6$), growing quadratically for stubby beams and for shear-compliant materials. (Do not confuse this ratio with the shear share of the total tip deflection, $3EI/(kGA\,L^2 + 3EI)$, which is what the shear-locking lab below reports: 0.77% and 3.0% at the same two aspect ratios. The two agree to within rounding here and diverge for stubbier beams.) Be precise about that last clause, because the obvious reading of it is wrong. The ratio is fixed by $E/kG$ and $(t/L)^2$, not by absolute stiffness: scale $E$ and $G$ down together — a softer structure — and it does not move at all, and for an isotropic material $E/G = 2(1+\nu)$ cannot exceed $3$ whatever the modulus ($2.6$ for the metals just quoted, at $\nu = 0.3$). What makes shear compliance matter is stubbiness, quadratically, and a $G$ that is low relative to the bending modulus — laminates and sandwich sections with compliant cores, where no Poisson's ratio bounds the ratio.

FE implications: Timoshenko needs only $C^0$ interpolation — but naively interpolated low-order elements shear-lock (Section 10); the anisoparametric Tessler–Dong family is the clean fix. And a counting exercise that pays off in Module 4: EB has 2 section strains $(e, \kappa)$; planar Timoshenko has 3 $(e, \kappa, \gamma)$; a full 3D Timoshenko frame member has 6 — stretching, torsion, two bendings, two transverse shears — exactly the set the beam iFEM of Gherlone et al. (2012) reconstructs.

The same bent beam under both kinematic hypotheses. Left (Euler-Bernoulli): the cross-section stays perpendicular to the deformed midline — $\theta = dw/dx$, $\gamma = 0$. Right (Timoshenko): the section stays plane but rotates less than the midline slope; the shaded wedge is $\gamma = dw/dx - \theta$, carried by $Q = kGA\gamma$. EB: 2 section strains $(e, \kappa)$ — $C^1$ elements. Timoshenko: 3 section strains $(e, \kappa, \gamma)$ — $C^0$ elements, locking risk.

9. Plate Kinematics: Kirchhoff vs Reissner-Mindlin

The plate story parallels the beam story one dimension up — with the same fork in the road and higher stakes, because $C^1$ continuity in 2D is genuinely hard.

Kirchhoff: Thin-Plate Theory

Material normals remain straight, unstretched, and normal to the deformed midsurface. The only unknown is $w(x,y)$; curvatures are second derivatives of $w$; the governing equation is biharmonic:

$$ D\,\nabla^4 w = q, \qquad D = \frac{E\,t^3}{12\,(1-\nu^2)} $$

Fourth order ⇒ $C^1$ continuity of $w$ across element edges — notoriously hard for triangles: conforming $C^1$ triangles require quintic polynomials in the classical single-polynomial (Argyris) construction — piecewise-cubic macroelements such as Hsieh–Clough–Tocher achieve it at lower degree. That difficulty historically drove nonconforming and discrete-Kirchhoff elements, and ultimately the $C^0$ alternative below.

Reissner-Mindlin: First-Order Shear Deformation Theory

Normals remain straight but not necessarily normal. The kinematics are $\mathbf{u} = (u_0 + z\theta_x,\; v_0 + z\theta_y,\; w)$ with independent rotations $\theta_x, \theta_y$ — a second-order system, so $C^0$ suffices.

Convention Declared (read before wiring anything)
This guide uses a plane-labelled convention: $\theta_x$ is the rotation of the material normal in the x–z plane (and $\theta_y$ in the y–z plane), entering the displacement field as $\mathbf{u} = \mathbf{u}_0 + z\boldsymbol{\theta}$, so that $$\boldsymbol{\varepsilon}(z) = \mathbf{e} + z\boldsymbol{\kappa}, \qquad \boldsymbol{\gamma} = \nabla w + \boldsymbol{\theta}$$ That lettering is not the iFEM sources' lettering, and the difference is worth writing down once. Tessler & Spangler name their rotations after the axes they turn about: in NASA/TM-2003-212445 (§2.1) $\theta_x$ and $\theta_y$ are “the rotations of the normal about the negative $x$ and positive $y$ axes”, so it is their $\theta_y$ that multiplies $z$ in the $x$-direction displacement — $u_x = u_0 + z\theta_y$, $u_y = v_0 + z\theta_x$. The translation is a pure swap, $\theta_x^{\text{here}} = \theta_y^{\text{TS}}$ and $\theta_y^{\text{here}} = \theta_x^{\text{TS}}$; with that swap the two equations above are theirs. Copy their symbols without it and every curvature component lands in the wrong slot. Two further conventions are common and both self-consistent: $\mathbf{u} = \mathbf{u}_0 - z\boldsymbol{\theta}$ with $\boldsymbol{\gamma} = \nabla w - \boldsymbol{\theta}$ (giving $\boldsymbol{\varepsilon} = \mathbf{e} - z\boldsymbol{\kappa}$), and the strict right-hand-rule convention of rotations about the positive axes, which swaps the two components and flips the sign of one of them ($u_y = v_0 - z\theta_x$; this is the form used by the later iQS4 element). Mixing any two of them silently corrupts the shear $\mathbf{B}$-matrix and every strain-recovery relation downstream — a plate reconstructed with the wrong curvature sign bends the wrong way. Declare one convention per project and convert at the boundary. Modules 4–6 use the one declared here.

The deformation state is captured by eight section strain measures: 3 membrane, 3 bending, 2 transverse shear:

$$ \mathbf{e} = \begin{bmatrix} \partial u_0/\partial x \\ \partial v_0/\partial y \\ \partial u_0/\partial y + \partial v_0/\partial x \end{bmatrix}, \qquad \boldsymbol{\kappa} = \begin{bmatrix} \partial\theta_x/\partial x \\ \partial\theta_y/\partial y \\ \partial\theta_x/\partial y + \partial\theta_y/\partial x \end{bmatrix}, \qquad \boldsymbol{\gamma} = \begin{bmatrix} \partial w/\partial x + \theta_x \\ \partial w/\partial y + \theta_y \end{bmatrix} $$

In the thin limit the physics enforces the Kirchhoff constraint $\boldsymbol{\gamma} \to \mathbf{0}$, and low-order $C^0$ elements shear-lock (Section 10).

Key insight — Reissner and Mindlin are not the same theory
Reissner (1945) derived a static, stress-based sixth-order theory yielding $k = 5/6$. Mindlin (1951) derived the displacement-based, dynamic theory and chose $k = \pi^2/12$ to match the thickness-shear frequency. The theories agree in engineering practice, but the correction factors differ because the matching criteria differ — one more reason $k$ is not a universal constant.

Well-behaved $C^0$ Mindlin elements exist: MITC4, the discrete shear gap family, and the anisoparametric 3-node MIN3 of Tessler & Hughes (1985) — quadratic $w$ constrained by continuous edge shear conditions, linear rotations. The anisoparametric treatment traces back to Tessler & Hughes (1983) on the four-node Mindlin quadrilateral. MIN3 is worth remembering by name: it is precisely the element that reappears as iMIN3 in the iFEM modules.

10. Shear Locking: Mechanism and Cures

Now that the shear-deformable kinematics exist, we can explain the pathology that haunts them. Take the simplest possible shear-deformable element: the 2-node Timoshenko beam element with linear $w$ and linear $\theta$ (DOFs $w_1, \theta_1, w_2, \theta_2$). Its element strains:

$$ \kappa = \frac{d\theta}{dx} \;\;(\text{constant}), \qquad \gamma = \frac{dw}{dx} - \theta \;\;(\text{linear in } x), \qquad U = \tfrac{1}{2}\int_0^{L} \left( EI\,\kappa^2 + kGA\,\gamma^2 \right) dx $$

$\gamma$ is linear because $\theta$ is linear while $w'$ is constant. Now demand pure bending: $\kappa = \text{const}$ with $\gamma = 0$ everywhere. But $\gamma \equiv 0$ with linear $\theta$ forces $w$ to be quadratic — unavailable in the linear space. The element cannot bend without storing parasitic shear energy, and being the energy-optimal answer (Section 2), it splits the difference: it bends less and pays shear.

The Scaling That Kills

Work the energies for a bending state with rotation amplitude $\bar\theta$: the parasitic linear-in-$x$ shear contributes $\propto kGA\,L_e\,\bar\theta^2/12$ against bending energy $\propto EI\,\bar\theta^2/L_e$, giving the ratio

$$ \frac{U_{\text{parasitic shear}}}{U_{\text{bending}}} \;\sim\; \frac{kGA\,L_e^2}{12\,EI} \;\propto\; \left(\frac{L_e}{t}\right)^{2} \quad \text{(rectangular section, } I = At^2/12\text{)} $$

The spurious constraint stiffens without bound as the structure thins: at fixed mesh, $w_{FE}/w_{\text{exact}} \to 0$ like $(t/L)^2$. With $k = 5/6$ and $G = E/2.6$, the ratio $kGA L_e^2/(12EI) = (kG/E)(L_e/t)^2$ is already $\approx 32$ at $L_e/t = 10$ and $\approx 3200$ at $L_e/t = 100$.

Key insight — locking is not a solver problem
The locked linear system is nonsingular and solves cleanly — the answer is simply wrong. (Thin elements can be ill-conditioned for unrelated reasons — extreme aspect ratios, badly scaled rotational vs translational DOFs — so a condition number is not a locking diagnostic in either direction.) Babuška & Suri's formal definition: locking is loss of uniform convergence with respect to the small parameter $t$; robustness means error bounds uniform in $t$. Diagnosing locking as ill-conditioning sends you to the wrong fix (better solver) instead of the right one (better discrete space).
Why it matters — the model limit vs the discrete limit
"Euler-Bernoulli = Timoshenko with $\gamma = 0$" is true for the continuum theories: the constraint limit is well-behaved. But imposing $\gamma_h = 0$ on low-order interpolants is precisely the mechanism of locking — the discrete space cannot satisfy the constraint without freezing. Distinguish the model limit from the discrete limit, or the sentence above turns from a theorem into a bug.

The Cures, Each With Its Mechanism

  1. Reduced integration of the shear term. 1-point quadrature samples $\gamma$ only at the midpoint — exactly where the linear parasitic part vanishes. Equivalent to an assumed-constant-shear element. Selective reduced integration generalizes this: full on bending, reduced on shear. Caveat, and it is not optional: uniform reduced integration can introduce spurious zero-energy (hourglass) modes — recall $\operatorname{rank}(\mathbf{k}^e) \le n_{gp}\times(\text{strain rows})$ from Section 4. Preset 3 of the lab shows them; the strip preset shows an entire mesh going singular.
  2. Assumed natural strain / MITC: re-interpolate $\gamma$ from tying points chosen so the parasitic part cannot be represented (Bathe's MITC family).
  3. Discrete shear gap: construct the shear field from nodal gap values rather than from $w' - \theta$ directly.
  4. Interdependent (anisoparametric) interpolation — Tessler & Dong's hierarchy: take $w$ one degree higher than $\theta$ (quadratic $w$, linear $\theta$), then condense the internal DOF via an element-level constant-shear constraint. The result is a locking-free 2-node element with no spurious modes. This idea becomes MIN3 for plates (Tessler & Hughes 1985, after the 1983 quadrilateral treatment) and later iMIN3 in iFEM.
  5. Mixed formulations (Hellinger–Reissner, Hu–Washizu): promote the offending stress or strain to a genuinely independent field with its own variation, so the displacement interpolation no longer has to produce it.
  6. $\bar{B}$ / selective projection: stay in the displacement formulation, but replace the offending strain component by its projection onto a lower-order space before it ever reaches the energy. (Reduced integration of the shear term is the quadrature-level version of the same idea.)
  7. Enhanced assumed strain (EAS): keep the displacement field and add element-internal enhanced strain modes that supply the missing deformation, then statically condense them out at element level.

Shear locking has two siblings: volumetric locking (incompressibility, $\nu \to \tfrac12$) and membrane locking (curved shells). Bombarde et al. (2024) give the modern comparative map of all three pathologies and their remedies, and Le & Bui (2023) show that locking-free FSDT plate implementation is still an active research topic — a half-century after the problem was identified, the element technology this module teaches (and that iFEM inherits) continues to be refined rather than replaced.

Single 2-node element, pure end moments exact: w quadratic, γ = 0 FE: w linear — the quadratic is unavailable γₖ(x) = w′ − θₖ ≠ 0 γₖ = 0 at midpoint γₖ ≠ 0 1-point shear quadrature samples γ at the midpoint, where the parasitic linear part vanishes Energy split (log scale) ≈32× L/t = 10 ≈3200× L/t = 100 bending energy (=1) parasitic shear ∝ (L/t)²
Why the linear Timoshenko element locks. The exact pure-bending deflection is quadratic with $\gamma = 0$ everywhere; the linear space cannot represent it, so the sections (red) rotate against a chord that stays straight, generating parasitic shear $\gamma_h$ (orange wedges) that is linear in $x$ and vanishes only at the midpoint. Right: the parasitic shear energy dwarfs the bending energy and grows as $(L/t)^2$ (values for $k=5/6$, $G=E/2.6$).
Interactive Tool — Shear-Locking Laboratory: Timoshenko Cantilever

A tip-loaded cantilever meshed with 2-node Timoshenko elements. Sweep the slenderness $L/t$, refine the mesh, and switch the shear quadrature. Full integration collapses like $(t/L)^2$; 1-point (reduced) shear stays honest. Overlay Euler-Bernoulli to see what it misses on stubby beams.

 

11. Section Strains: The Bridge to iFEM

This module chose Timoshenko and Reissner-Mindlin deliberately: iFEM is built on them. The FSDT plate/shell displacement field is fully parametrized by $(u_0, v_0, w, \theta_x, \theta_y)$, and its deformation state is captured by the 8 section strain measures of Section 9 — membrane $\mathbf{e}(\mathbf{u})$ (3), bending curvatures $\boldsymbol{\kappa}(\mathbf{u})$ (3), transverse shears $\boldsymbol{\gamma}(\mathbf{u})$ (2). These are the strain measures of Tessler & Spangler's (2005) least-squares variational principle — the same eight, up to the rotation relabelling declared in Section 9 — which reconstructs displacements by matching these analytic measures to experimentally measured ones.

The Measurement Link

With strain rosettes on the top and bottom surfaces $z = \pm t/2$, the in-plane surface strains under the through-thickness distribution $\boldsymbol{\varepsilon}(z) = \mathbf{e} + z\boldsymbol{\kappa}$ are:

$$ \boldsymbol{\varepsilon}^{\pm} = \mathbf{e} \pm \frac{t}{2}\,\boldsymbol{\kappa} \quad\Rightarrow\quad \mathbf{e} = \frac{\boldsymbol{\varepsilon}^{+}+\boldsymbol{\varepsilon}^{-}}{2}, \qquad \boldsymbol{\kappa} = \frac{\boldsymbol{\varepsilon}^{+}-\boldsymbol{\varepsilon}^{-}}{t} $$

Semi-sum and semi-difference: a back-to-back sensor pair splits into a membrane part and a bending part. That is the bridge this module exists to build — the forward machinery hands the inverse problem a quantity it can actually measure. The general version (all three membrane strains, all three curvatures, rosette transformations, and what happens when you can only instrument one surface) is derived in Module 5, which is its canonical home.

Why it matters — fix the sign convention before wiring sensors
The recovery formulas above follow directly from the convention declared in Section 9 ($\mathbf{u} = \mathbf{u}_0 + z\boldsymbol{\theta} \Rightarrow \boldsymbol{\varepsilon}(z) = \mathbf{e} + z\boldsymbol{\kappa}$), which is the convention adopted here and in Modules 4–6 — not the lettering of the iFEM sources, which Section 9 converts at the boundary (iQS4 writes $u_x = u + z\theta_y$, $u_y = v - z\theta_x$; Kefal et al. 2016, eqs. 4a–b). What does survive that conversion is precisely the structure these formulas rest on: the same sources write the in-plane strain as $\mathbf{e} + z\boldsymbol{\kappa}$ (eq. 6a) and recover the curvature as the same semi-difference of surface strains (eq. 7b). The rotation labels are negotiable; the semi-sum/semi-difference is not. Adopt the opposite rotation sign ($\mathbf{u} = \mathbf{u}_0 - z\boldsymbol{\theta}$) and the through-thickness distribution becomes $\mathbf{e} - z\boldsymbol{\kappa}$, so the curvature recovery flips to $\boldsymbol{\kappa} = (\boldsymbol{\varepsilon}^{-}-\boldsymbol{\varepsilon}^{+})/t$. Same physics, opposite sign — and a curvature recovered with the wrong sign reconstructs a plate bending the wrong way. Two habits prevent this: write the convention next to the equation, and check the recovered sign against a known load case before trusting any sensor data. One more thing to pin down before wiring: whether your $t$ is the total thickness (as here, with surfaces at $z=\pm t/2$) or the half-thickness — that choice is a factor of two on every curvature. Modules 4–6 follow the iFEM literature and place the surfaces at $z = \pm h$, so their total thickness is $2h$ and the curvature recovery reads $(\boldsymbol{\varepsilon}^{+}-\boldsymbol{\varepsilon}^{-})/2h$; same plate, same physics, different letter.

And a structural observation with a big consequence: transverse shear strains cannot be measured on the surfaces — they vanish there in reality, and FSDT models them as constant through the thickness. iFEM handles this by down-weighting the $\boldsymbol{\gamma}$ terms in its least-squares functional. This is a structural reason why FSDT — not 3D elasticity, not Kirchhoff — is the right kinematic frame for shape sensing: it exposes exactly the strain measures that surface sensors can and cannot see. For beams and frames, the analogous set is the 6 Timoshenko section strains used by Gherlone et al. (2012).

Forward and Inverse, Side by Side

Forward FEM computes $\mathbf{d}$ from loads by minimizing $\Pi = \tfrac12\mathbf{d}^{\mathsf T}\mathbf{K}\mathbf{d} - \mathbf{d}^{\mathsf T}\mathbf{F}$, built from $\mathbf{B}$-operators. iFEM computes $\mathbf{d}$ from measured strains by minimizing a least-squares functional over the same interpolation:

$$ \Phi(\mathbf{d}) = \sum_{i} \left[\, w_e \left\| \mathbf{e}(\mathbf{d}) - \mathbf{e}^{\varepsilon}_i \right\|^2 + w_\kappa\, t^{2} \left\| \boldsymbol{\kappa}(\mathbf{d}) - \boldsymbol{\kappa}^{\varepsilon}_i \right\|^2 + w_\gamma \left\| \boldsymbol{\gamma}(\mathbf{d}) - \boldsymbol{\gamma}^{\varepsilon}_i \right\|^2 \,\right] $$

The $t^2$ on the curvature term is part of the definition, not decoration: entries of $\mathbf{e}$ and $\boldsymbol{\gamma}$ are dimensionless strains while entries of $\boldsymbol{\kappa}$ have units of 1/length, so without it the three residuals are not commensurable and $w_e, w_\kappa, w_\gamma$ stop being comparable numbers. Module 5 writes the same factor as $(2h)^2$ in the half-thickness lettering of the inverse-element literature (e.g. the iQS4 formulation of Kefal et al. 2016, eq. 8c), where it is likewise mandatory.

Forward FEMiFEM (Module 4)
Minimizes$\Pi(\mathbf{d}) = \tfrac12\mathbf{d}^{\mathsf T}\mathbf{K}\mathbf{d} - \mathbf{d}^{\mathsf T}\mathbf{F}$$\Phi(\mathbf{d})$: weighted strain mismatch
Needs loads?yes — $\mathbf{F}$ drives everythingno
Needs material stiffness?yes — $\mathbf{D}$ inside $\mathbf{K}$no — only relative weights $w_e, w_\kappa, w_\gamma$
Normal equationssparse SPD $\mathbf{K}\mathbf{U} = \mathbf{F}$sparse symmetric PSD by construction ($\mathbf{A}^{\mathsf T}\mathbf{W}\mathbf{A}$) — SPD only once rigid modes are constrained and the sensor set gives the residual operator full column rank
Element continuity$C^0$ (FSDT)$C^0$ — interfaces naturally with discrete sensors
Key insight
Everything in Sections 3–6 — shape functions, isoparametric mapping, quadrature, assembly, essential constraints — is reused verbatim by iFEM. What changes is only the functional being minimized: energy with loads becomes least-squares with measurements. If you can implement the forward pipeline of this module, you are one functional away from a shape-sensing code.
z +t/2 −t/2 0 membrane e (constant) bending z·κ (zero at midplane) ε⁺ (top sensor) ε⁻ (bottom sensor) total ε(z) = e + z·κ e = (ε⁺ + ε⁻)/2     κ = (ε⁺ − ε⁻)/t the two surfaces are exactly where iFEM strain sensors sit; transverse shear γ is invisible to them + =
Through-thickness strain decomposition. A constant membrane part plus a linear bending part sum to the total $\varepsilon(z) = e + z\kappa$. Sensors at $z = \pm t/2$ read $\varepsilon^\pm$; semi-sum and semi-difference recover $e$ and $\kappa$ — the section strains iFEM matches.

Flashcards

References

← 1. Continuum Mechanics Primer 3. Inverse Problems & Regularization →