Nitsche Method : Fully tied 3D Contact

Introduction

In this article, we take advantage of this capability to solve an optimization problem that fully constrains three-dimensional displacements using the mortar method.

This approach offers a key advantage — it enables the treatment of various combinations of bodies without altering their mesh structures.

Algorithm

Nitche Method

Instead of enforcing the contact constraint through a Lagrange multiplier (which introduces extra unknowns), Nitsche’s method weakly imposes the constraint by modifying the weak form with consistent terms.

The key idea is to:

• – Add penalty-like terms to stabilize the interface,
• – Add consistent adjoint terms to maintain symmetry and consistency.

This results in an augmented bilinear form on the contact surface.

In the general form of force equilibrium equation, we will use below weak form

\int_{\Omega} \boldsymbol{\varepsilon}(\boldsymbol{v}) : \boldsymbol{\sigma}(\boldsymbol{u}) , d\Omega = \int_{\Omega} \boldsymbol{v} \cdot \boldsymbol{f} , d\Omega
+
\int_{\Gamma_t} \boldsymbol{v} \cdot \bar{\boldsymbol{t}} , d\Gamma

And we would obtain linear equation from the weak form

\begin{align}
\mathbf{K}u =F_{N}+F_{body}\\
\begin{bmatrix}
\mathbf{K}_1 & \mathbf{0} \\
\mathbf{0} & \mathbf{K}_2
\end{bmatrix}
\begin{bmatrix}
\mathbf{u}_1 \\[4pt]
\mathbf{u}_2
\end{bmatrix}
=
\begin{bmatrix}
\mathbf{f}_1 \\[4pt]
\mathbf{f}_2
\end{bmatrix}
\end{align}

In the Nitsche method, a penalty-like matrix is added to the stiffness matrix obtained from the above weak formulation by using bilinear form on the contact surface as it is mentioned above,

  \begin{aligned}
  a_\Gamma(u_1, u_2; v_1, v_2)
  =
  \int_{\Gamma}
  \Big[
  \frac{\alpha}{h} (u_1 - u_2) \cdot (v_1 - v_2)\\
  - \frac{1}{2} \big( (\boldsymbol{t}_1(u_1) - \boldsymbol{t}_2(u_2)) \cdot (v_1 - v_2) \big)\\
  - \frac{1}{2} \big( (\boldsymbol{t}_1(v_1) - \boldsymbol{t}_2(v_2)) \cdot (u_1 - u_2) \big)
  \Big] \ ds
  \end{aligned}

\boldsymbol{t}_i(u_i) = \boldsymbol{C} : \boldsymbol{\varepsilon}(u_i)  \boldsymbol{n},
\qquad
\boldsymbol{\varepsilon}(u_i) = \tfrac{1}{2}(\nabla u_i + \nabla u_i^\mathsf{T})\\

where

\boldsymbol{C} is the elasticity tensor,
\boldsymbol{n} is the unit normal vector on the interface ( \Gamma ),
h is the characteristic element size,
and \alpha is the penalty parameter.

From this form, we would get K_{Nitsche}.

(K+K_{Nitsche})u=F_{N}+F_{body}

The meaning of 1st term is that it prevents jump of displacement on the interface. Additionally, the meaning of 2nd and 3rd terms are kind of energy (like compliance).

Implementation with scikit-fem

Code

import numpy as np
import scipy.sparse as sp
from scipy.sparse.linalg import minres, spsolve

import skfem
from skfem import FacetBasis
from skfem.helpers import dot, mul, sym_grad, trace
from skfem.supermeshing import elementwise_quadrature, intersect

from common import (
    assemble_top_pressure_load,
    make_default_meshes,
    make_elasticity_form,
    make_vector_p1_bases,
)

# ----------------------------
# 0) Meshes / bases / material
# ----------------------------
mesh_top, mesh_bot = make_default_meshes()
elem, basis_top, basis_bot = make_vector_p1_bases(mesh_top, mesh_bot)
a_elastic, E, nu, mu, lam = make_elasticity_form(E=210e9, nu=0.30)

K1 = a_elastic.assemble(basis_top).tocsr()
K2 = a_elastic.assemble(basis_bot).tocsr()

# ----------------------------
# 1) Supermesh on contact facets (non-matching)
# ----------------------------
m1t, orig1 = mesh_top.trace("contact", mtype=skfem.MeshTri, project=lambda p: p[[0, 1]])
m2t, orig2 = mesh_bot.trace("contact", mtype=skfem.MeshTri, project=lambda p: p[[0, 1]])

m12, t1, t2 = intersect(m1t, m2t)
quad1 = elementwise_quadrature(m1t, m12, t1)
quad2 = elementwise_quadrature(m2t, m12, t2)

fb_top = FacetBasis(mesh_top, elem, facets=orig1[t1], quadrature=quad1)
fb_bot = FacetBasis(mesh_bot, elem, facets=orig2[t2], quadrature=quad2)

# ----------------------------
# 2) Symmetric Nitsche for tied interface: u_top = u_bot
# ----------------------------
gamma = 20.0
alpha0 = lam + 2.0 * mu


@skfem.BilinearForm
def a_nitsche_tied(u1, u2, v1, v2, w):
    # Important skfem detail:
    # for CompositeBasis(fb_top, fb_bot), default_parameters() comes from
    # bases[0].default_parameters(). Therefore w.n and w.h here are inherited
    # from fb_top only, not from both sides of the interface.
    #
    # In this file that means:
    # - w.n is the top-side normal (here approximately [0, 0, -1])
    # - w.h is the top-side facet length
    #
    # Under this convention, the interface traction contribution should use
    # opposite signs across the interface to match the two outward normals.
    n = w.n

    ju = u1 - u2
    jv = v1 - v2

    eps_u1 = sym_grad(u1)
    eps_u2 = sym_grad(u2)
    eps_v1 = sym_grad(v1)
    eps_v2 = sym_grad(v2)

    t1 = 2.0 * mu * mul(eps_u1, n) + lam * trace(eps_u1) * n
    t2 = 2.0 * mu * mul(eps_u2, n) + lam * trace(eps_u2) * n
    s1 = 2.0 * mu * mul(eps_v1, n) + lam * trace(eps_v1) * n
    s2 = 2.0 * mu * mul(eps_v2, n) + lam * trace(eps_v2) * n

    # TODO:
    # Replace this one-sided penalty length by a two-sided scale derived from
    # both meshes, e.g. h_top/h_bot -> h_eff. At the moment w.h is inherited
    # from fb_top only.
    penalty = gamma * alpha0 / w.h

    return (
        penalty * dot(ju, jv)
        - 0.5 * dot(t1 - t2, jv)
        - 0.5 * dot(s1 - s2, ju)
    )


# NOTE:
# Assembled on the supermesh product basis. Because CompositeBasis forwards
# default parameters from the first basis, the interface normal/mesh size used
# here are tied to fb_top unless passed explicitly another way.
B = skfem.asm(a_nitsche_tied, fb_top * fb_bot).tocsr()
K = sp.bmat([[K1, None], [None, K2]], format="csr") + B


def report_interface_geometry() -> None:
    n_top = np.asarray(fb_top.normals, dtype=float)
    n_bot = np.asarray(fb_bot.normals, dtype=float)
    h_top = np.asarray(fb_top.mesh_parameters(), dtype=float)
    h_bot = np.asarray(fb_bot.mesh_parameters(), dtype=float)

    normals_top_unique = np.unique(np.round(n_top.reshape(3, -1).T, 6), axis=0)
    normals_bot_unique = np.unique(np.round(n_bot.reshape(3, -1).T, 6), axis=0)

    print("interface geometry diagnostics:")
    print("  top normals unique =", normals_top_unique)
    print("  bot normals unique =", normals_bot_unique)
    print("  top h min/max =", float(h_top.min()), float(h_top.max()))
    print("  bot h min/max =", float(h_bot.min()), float(h_bot.max()))


report_interface_geometry()


def stack_vector_field(u_top_nodes: np.ndarray, u_bot_nodes: np.ndarray) -> np.ndarray:
    return np.hstack([u_top_nodes.reshape(-1), u_bot_nodes.reshape(-1)])


def eval_affine_on_nodes(mesh, A: np.ndarray, b: np.ndarray) -> np.ndarray:
    x = np.asarray(mesh.p.T, dtype=float)
    return x @ A.T + b[None, :]


def report_interface_patch_tests() -> None:
    coords_top = np.asarray(mesh_top.p.T, dtype=float)
    coords_bot = np.asarray(mesh_bot.p.T, dtype=float)

    tests: list[tuple[str, np.ndarray]] = []

    # Rigid translations: interface residual should be ~0 because sym_grad = 0
    # and jump = 0.
    for comp, name in enumerate(("tx", "ty", "tz")):
        ut = np.zeros_like(coords_top)
        ub = np.zeros_like(coords_bot)
        ut[:, comp] = 1.0
        ub[:, comp] = 1.0
        tests.append((f"rigid_{name}", stack_vector_field(ut, ub)))

    # Continuous affine fields: jump = 0 on the interface. For a consistent
    # tied-interface formulation, the interface residual B @ u should also be
    # small for these patch-test fields. If these are large while rigid modes
    # are small, that points to a sign/normal/weighting issue in the interface
    # consistency terms rather than a generic assembly failure.
    affine_cases = {
        "affine_exx": (np.array([[1.0, 0.0, 0.0],
                                  [0.0, 0.0, 0.0],
                                  [0.0, 0.0, 0.0]], dtype=float), np.zeros(3)),
        "affine_eyy": (np.array([[0.0, 0.0, 0.0],
                                  [0.0, 1.0, 0.0],
                                  [0.0, 0.0, 0.0]], dtype=float), np.zeros(3)),
        "affine_shear_xy": (np.array([[0.0, 1.0, 0.0],
                                       [1.0, 0.0, 0.0],
                                       [0.0, 0.0, 0.0]], dtype=float), np.zeros(3)),
    }
    for name, (A, b) in affine_cases.items():
        ut = eval_affine_on_nodes(mesh_top, A, b)
        ub = eval_affine_on_nodes(mesh_bot, A, b)
        tests.append((name, stack_vector_field(ut, ub)))

    print("interface patch diagnostics:")
    for name, u_mode in tests:
        iface_res = np.asarray(B @ u_mode, dtype=float)
        print(f"  {name}: ||B u||2 = {np.linalg.norm(iface_res):.6e}")


report_interface_patch_tests()

# ----------------------------
# 3) External load + Dirichlet
# ----------------------------
F_top, area, pressure = assemble_top_pressure_load(mesh_top, elem, total_force=100.0)

F = np.zeros(K.shape[0])
F[:basis_top.N] = F_top

D_bot = basis_bot.get_dofs("dirichlet").all()
D = D_bot + basis_top.N

# Bottom displacement dofs are stacked as [ux0, uy0, uz0, ux1, ...] in the
# second block, so we can recover vector reactions by component modulo 3.
fixed_dofs = np.asarray(D, dtype=int)

Kc, Fc, x0, I = skfem.condense(K, F, D=D)

try:
    u_red = spsolve(Kc, Fc)
except Exception:
    u_red, info = minres(Kc, Fc, rtol=1e-10, maxiter=20000)
    if info != 0:
        raise RuntimeError(f"MINRES did not converge, info={info}")

u = np.zeros(K.shape[0])
u[I] = u_red

res = K @ u - F
r_free = res[I]
reac = np.asarray(res[fixed_dofs], dtype=float)

reac_vec = np.zeros(3, dtype=float)
for comp in range(3):
    reac_vec[comp] = float(np.sum(reac[(fixed_dofs % 3) == comp]))

ext_vec = np.zeros(3, dtype=float)
for comp in range(3):
    ext_vec[comp] = float(np.sum(F[np.arange(K.shape[0]) % 3 == comp]))

balance_vec = ext_vec + reac_vec
print("K:", K.shape, "nnz:", K.nnz)
# Full residual includes reactions on condensed Dirichlet dofs, so it is not a
# clean solve-quality indicator by itself.
print("||res||2 =", np.linalg.norm(res))
print("relative ||res||2 / ||F||2 =", np.linalg.norm(res) / (np.linalg.norm(F) + 1e-30))
print("||r_free||2 =", np.linalg.norm(r_free))
print("relative ||r_free||2 / ||F||2 =", np.linalg.norm(r_free) / (np.linalg.norm(F) + 1e-30))
print("external resultant =", ext_vec)
print("dirichlet reaction resultant =", reac_vec)
print("balance external + reaction =", balance_vec)
print("||balance||2 =", np.linalg.norm(balance_vec))


# ----------------------------
# 4) Export combined VTU
# ----------------------------
def export_combined_vtu(mesh_top, mesh_bot, u, basis_top, file_path="combined.vtu"):
    import meshio

    n_top = mesh_top.p.shape[1]
    n_bot = mesh_bot.p.shape[1]
    points = np.vstack([mesh_top.p.T, mesh_bot.p.T])
    cells = [("tetra", np.vstack([mesh_top.t.T, mesh_bot.t.T + n_top]))]

    u_top = u[:basis_top.N].reshape(n_top, 3)
    u_bot = u[basis_top.N:].reshape(n_bot, 3)
    u_all = np.vstack([u_top, u_bot])

    node_tag = np.zeros(n_top + n_bot, dtype=np.int32)
    top_nodes_contact = np.unique(mesh_top.facets[:, mesh_top.boundaries["contact"]])
    bot_nodes_contact = np.unique(mesh_bot.facets[:, mesh_bot.boundaries["contact"]]) + n_top
    node_tag[np.hstack([top_nodes_contact, bot_nodes_contact])] = 1

    bot_nodes_dir = np.unique(mesh_bot.facets[:, mesh_bot.boundaries["dirichlet"]]) + n_top
    node_tag[bot_nodes_dir] = 2

    top_nodes_force = np.unique(mesh_top.facets[:, mesh_top.boundaries["force"]])
    node_tag[top_nodes_force] = 3

    meshio.Mesh(
        points=points,
        cells=cells,
        point_data={"u": u_all, "node_tag": node_tag},
    ).write(file_path)

    print(f"wrote {file_path}")


export_combined_vtu(mesh_top, mesh_bot, u, basis_top, "tied_interface_nitsche.vtu")

  K: (18000, 18000) nnz: 673907
  ||res||2 = 5.299997797110098
  relative ||res||2 / ||F||2 = 1.0341738192729362
  ||r_free||2 = 1.5306548416229466e-12
  relative ||r_free||2 / ||F||2 = 2.986724191494849e-13
  external resultant = [   0.    0. -100.]
  dirichlet reaction resultant = [-1.38555833e-13  8.69304628e-14  1.00000000e+02]
  balance external + reaction = [-1.38555833e-13  8.69304628e-14 -1.36282097e-11]
  ||balance||2 = 1.3629191226373108e-11