Penalty, Lagrange multipliers, and Nitsche

One constraint, three ways to enforce it

The model above is two rectangular blocks of plane-strain Q4 elements meeting at \(x = W/2\). The two meshes line up geometrically, but the nodes on the seam are duplicated: node \(a\) belongs to the left block, node \(b\) to the right block, and nothing in the assembled stiffness matrix connects them. Solve as-is and you get two independent bodies. What we want is the displacement to be continuous across the seam,

\[ [\![\mathbf u]\!] \;=\; \mathbf u_A - \mathbf u_B \;=\; \mathbf 0 \qquad\text{on }\Gamma, \]

which in discrete form is a set of linear rows acting on the degrees of freedom, \(\mathbf C\mathbf u = \mathbf g\), with \(\mathbf g = \mathbf 0\) here. That is the same algebraic object as a fixed support, a rigid tie, a periodic boundary condition, or a mortar continuity condition. How you enforce those rows is the entire subject.

This matters well beyond two glued rectangles. Non-matching meshes between parts of an assembly, embedded and cut-cell methods where the boundary does not follow element faces, contact between bodies, and weakly imposed Dirichlet data are all the same question wearing different clothes.

The penalty method

Append a quadratic term to the potential energy that charges for any opening of the seam. For a tied node pair the term is \(\tfrac12\alpha\,\lVert\mathbf u_A - \mathbf u_B\rVert^2\), and its second derivative is a small \(4\times4\) block scattered into the global matrix:

\[ \mathbf K \;\leftarrow\; \mathbf K + \alpha\begin{bmatrix}\mathbf I & -\mathbf I\\ -\mathbf I & \mathbf I\end{bmatrix},\qquad\text{globally}\qquad (\mathbf K + \alpha\,\mathbf C^{\mathsf T}\mathbf C)\,\mathbf u = \mathbf f + \alpha\,\mathbf C^{\mathsf T}\mathbf g. \]

Nothing about the unknown vector changes, and if \(\mathbf K\) was positive definite it stays positive definite — you can keep using Cholesky or conjugate gradients. That is why the penalty method is everywhere.

The cost is that \(\alpha\) is a lie you tell the model. The seam is not tied; it is glued with a very stiff spring, and the residual opening is \(O(1/\alpha)\). Drag θ in the figure and watch the jump fall by exactly one decade per decade of θ, with the error following. So push \(\alpha\) up. But the condition number rises one decade per decade too, and once \(\kappa\) approaches \(1/\varepsilon_{\text{mach}} \approx 10^{16}\) the arithmetic stops delivering the digits the formulation promises. Panel A shows the resulting valley: the error falls, bottoms out around \(\theta \sim 10^{6}\), and then climbs back up. Past the bottom you are paying in round-off for a constraint you cannot actually resolve.

Lagrange multipliers

Instead of approximating the constraint, add unknowns that enforce it. Make the potential stationary with respect to both \(\mathbf u\) and a new vector \(\boldsymbol\lambda\), one multiplier per constraint row:

\[ \Pi(\mathbf u,\boldsymbol\lambda) = \tfrac12\mathbf u^{\mathsf T}\mathbf K\mathbf u - \mathbf u^{\mathsf T}\mathbf f + \boldsymbol\lambda^{\mathsf T}(\mathbf C\mathbf u - \mathbf g) \;\Longrightarrow\; \begin{bmatrix}\mathbf K & \mathbf C^{\mathsf T}\\ \mathbf C & \mathbf 0\end{bmatrix}\begin{bmatrix}\mathbf u\\ \boldsymbol\lambda\end{bmatrix} = \begin{bmatrix}\mathbf f\\ \mathbf g\end{bmatrix}. \]

The bottom block row is the constraint, so it holds exactly — the reported jump is at the level of machine round-off, and it stays there no matter how the mesh or the materials change. There is no parameter to choose. The multipliers are not bookkeeping either: \(\boldsymbol\lambda\) is the reaction that enforces the tie, so dividing it by the tributary length of each node recovers the interface traction. In the traction panel it lands on top of the \(\sigma_{xx}\) recovered independently from the volume elements.

What you give up is the character of the matrix. The zero block on the diagonal makes the system a saddle point problem: symmetric but indefinite. With \(m\) constraint rows the augmented matrix has exactly \(m\) negative eigenvalues, which the live table reports. A Cholesky factorisation will fail on it, and conjugate gradients is not applicable — you need LDLᵀ with pivoting, MINRES, or a block-preconditioned Krylov method. The system is also larger, by two unknowns per tied node here.

One more caveat that this example is too well-behaved to show: the multiplier space cannot be chosen freely. If λ is richer than the trace of the displacement space can support, the discrete inf–sup (LBB) condition fails and the constraint over-determines the problem, producing locking or a singular system. With matching meshes and one multiplier per node pair, as here, the pairing is stable. On genuinely non-matching meshes it is the central difficulty, and it is what mortar methods exist to solve.

Nitsche's method

Nitsche's method keeps the penalty term but repairs what makes penalty wrong. Derive the weak form honestly: integrating by parts over two subdomains leaves a flux term on the seam, \(\int_\Gamma [\![\mathbf v]\!]\cdot\langle\boldsymbol\sigma(\mathbf u)\mathbf n\rangle\), which the penalty method simply discards. Put it back, add its symmetric twin so the matrix stays symmetric, and keep a stabilising penalty:

\[ a_N(\mathbf u,\mathbf v) = \underbrace{\int_\Omega \boldsymbol\sigma(\mathbf u):\boldsymbol\varepsilon(\mathbf v)\,d\Omega}_{\text{volume}} \;-\; \underbrace{\int_\Gamma [\![\mathbf v]\!]\cdot\langle\boldsymbol\sigma(\mathbf u)\mathbf n\rangle\,d\Gamma}_{\text{consistency}} \;-\; \underbrace{\int_\Gamma [\![\mathbf u]\!]\cdot\langle\boldsymbol\sigma(\mathbf v)\mathbf n\rangle\,d\Gamma}_{\text{symmetry}} \;+\; \underbrace{\gamma\int_\Gamma [\![\mathbf u]\!]\cdot[\![\mathbf v]\!]\,d\Gamma}_{\text{stabilisation}} \]

with \(\langle\cdot\rangle\) the average across the seam and \(\gamma = \beta\,E_{\max}/h\). The name for what the extra terms buy is variational consistency: substitute the exact solution and the identity holds, because the exact solution has zero jump (killing the last two terms) and continuous traction (the first seam term balancing the volume integration by parts). The penalty method fails this test for every finite \(\alpha\).

That is not an abstraction, and the figure will show it to you. Set the stiffness ratio to 1 and the load to uniaxial stretch. Now the exact solution — uniform strain in each block — lies inside the Q4 space, so a consistent method must reproduce it exactly. Nitsche returns it to fourteen digits at every β it remains stable at, with a condition number around \(10^4\); penalty at the same parameter is wrong in the second decimal place, and would need \(\kappa\sim10^{16}\) to catch up — by which point round-off has already taken the digits back. Add a stiffness contrast or switch to the shear panel and neither method is exact any more, but Nitsche stays roughly an order of magnitude more accurate at equal θ, and it gets there on a better-conditioned matrix. That is the comparison in panel C.

Nitsche's price is β. The two consistency terms are not sign-definite, so they can eat the coercivity of the volume term; the penalty term is what holds the form coercive, and it only does so above a critical \(\beta_{\text{crit}}\) determined by a discrete trace inequality. Below it the matrix acquires negative eigenvalues. A direct solver may still hand you a plausible-looking answer — the live table will show it doing exactly that — but the formulation has lost its stability guarantee and Cholesky and CG are off the table. The useful part is that \(\beta_{\text{crit}}\) depends on the element type and polynomial degree, not on \(h\): refine the mesh in the controls and watch the shaded band stay put. That is why a fixed, modest β such as the 64 used in the source scripts is safe across a mesh study, whereas a penalty α must be re-argued every time.

Which one to reach for

In practice the augmented Lagrangian method sits between the first two: iterate the penalty solve while updating a multiplier estimate, and the constraint converges to machine precision without ever needing \(\alpha\) large enough to wreck the conditioning. It is the standard choice in contact mechanics, where the constraint set is an inequality that changes from iteration to iteration and the exactness of Lagrange is wanted without a growing indefinite system.

What is actually being computed

Everything on this page is solved in your browser when you move a control — there is no precomputed data. Each update assembles the plane-strain Q4 stiffness for both blocks, applies the interface treatment, eliminates the prescribed degrees of freedom, and factorises the result with LU and partial pivoting. Alongside it the same mesh is solved with the seam nodes merged, which is the conforming discretisation the tie is trying to reproduce and the reference every error is measured against. The condition number comes from Hager's 1-norm estimator applied to the LU factors, and the eigenvalue signs from the diagonal of an LDLᵀ factorisation, which is why the table can state whether a matrix is positive definite rather than guessing.

The formulation follows a set of MATLAB research scripts one-for-one — PFEM_two_block_penalty.m, PFEM_two_block_lagrange.m, PFEM_two_block_nitsche.m and the assemble_interface_* routines they call. The mesh generator duplicates the seam nodes exactly the way meshTwoBlocksQ4.m does, the Nitsche terms use the same \(\mathbf n=[1,0]\) convention and the same \(\gamma=\beta E_{\max}/h\) scaling, and the geometry is a rectangle split in half with the left edge held and the right edge driven. Two things are deliberately different: the aspect ratio is landscape so the seam reads clearly on screen, and the default load is a stretch rather than a compression so the seam opens toward you instead of interpenetrating.