About
In Scikit-FEM, the supermeshing feature allows handling contact problems even when two bodies have different or nonconforming meshes.
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. In the previous approach, Nitsche method, which is penalty based method to handle contact between surfaces, is shown. In this article, I will show how to implement Mortar method, which is Lagrangian based method, with scikit-fem.
Algorithm
Mortar Method
Basic Formulation
Since contact I assume is tied contact, the weak form on the contact is
\begin{aligned}
a_c(\mathbf{u}_1, \mathbf{u}_2; \mathbf{v}_\lambda)
&= \int_{\Gamma_c} \mathbf{v}_\lambda \cdot (\mathbf{u}_1 - \mathbf{u}_2) \, \mathrm{d}\Gamma \\[6pt]
&= \int_{\Gamma_c} \mathbf{v}_\lambda \cdot \mathbf{u}_1 \, \mathrm{d}\Gamma
- \int_{\Gamma_c} \mathbf{v}_\lambda \cdot \mathbf{u}_2 \, \mathrm{d}\Gamma
\end{aligned}The discretization and test function are
\mathbf{u}_1 =
\sum_j
\phi_j^{(1)} \mathbf{d}_j^{(1)}
\\
\mathbf{u}_2 =
\sum_k
\phi_k^{(2)} \mathbf{d}_k^{(2)}
\\
\boldsymbol{\lambda}
=
\sum_i
\mu_i \lambda_i
\\
\int_{\Gamma_c}
\mu_i
\left(
\sum_j \phi_j^{(1)} d_j^{(1)}
-
\sum_k \phi_k^{(2)} d_k^{(2)}
\right)
d\Gamma
= 0
where d is displacement, \phi is shape function and \mu shape function of Lagrange multipliers
Each B matrices are introduced with
(\mathbf{B}_1)_{ij}
=
\int_{\Gamma_c}
\mu_i \phi_j^{(1)} d\Gamma
,
\qquad
(\mathbf{B}_2)_{ik}
=
\int_{\Gamma_c}
\mu_i \phi_k^{(2)} d\Gamma
equation between us is
\mathbf{B}_1 \mathbf{u}_1
-
\mathbf{B}_2 \mathbf{u}_2
= 0The global formulation is
\underbrace{
\begin{bmatrix}
\mathbf{K}_1 & \mathbf{0} & \mathbf{B}_1^{\mathrm{T}} \\
\mathbf{0} & \mathbf{K}_2 & \mathbf{B}_2^{\mathrm{T}} \\
\mathbf{B}_1 & \mathbf{B}_2 & \mathbf{0}
\end{bmatrix}
}_{\text{global stiffness matrix } \mathbf{A}_{\mathrm{aug}}}
\begin{bmatrix}
\mathbf{u}_1 \\[4pt]
\mathbf{u}_2 \\[4pt]
\boldsymbol{\lambda}
\end{bmatrix}
=
\begin{bmatrix}
\mathbf{f}_1 \\[4pt]
\mathbf{f}_2 \\[4pt]
\mathbf{0}
\end{bmatrix}The difficulty of this task
For sure this formulation itself is difficult, but what worse is that the stability of solution for this formula because there is a problem that is so called “saddle point”.
To solve this equation, there are several techniques as stabilizers should be introduced, Dual Mortar, Augmented Lagrangian for example, but I will show it later on. Here, I will use implementation with Augmented Lagrangian.
Implementation
"""
Tied Contact (mortar) for two overlapping 3D elastic blocks in scikit-fem.
Features:
- Fine mortar with vector P0 multiplier on the contact intersection (stable baseline)
- Optional coarse mortar (aggregate fine lambda dofs -> coarse dofs)
- Two solvers:
(A) Pure Lagrange multiplier (KKT saddle-point)
(B) Augmented Lagrangian (AL) outer iteration (more robust)
- PyVista VTU output as MultiBlock (top + bottom in one .vtu/.vtm)
Dependencies:
pip install scikit-fem scipy numpy pyvista
Notes:
- This is TIED contact (u_top == u_bot on overlap). Not unilateral contact.
- For AL, gamma needs scaling. Start with gamma ~ E / h where h is contact mesh size.
"""
import numpy as np
from scipy.sparse import bmat, coo_matrix
from scipy.sparse.linalg import minres
import skfem
from skfem import Basis, FacetBasis
from skfem.helpers import dot, ddot, sym_grad
from skfem.supermeshing import intersect, elementwise_quadrature
from skfem import (
ElementTetP1,
ElementVectorH1,
ElementTetP0,
ElementVector,
)
from common import (
assemble_top_pressure_load,
make_default_meshes,
make_elasticity_form,
make_vector_p1_bases,
)
# -----------------------------
# User knobs
# -----------------------------
SOLVER = "AL" # "KKT" or "AL"
USE_COARSE_MORTAR = False
COARSE_GRID_NX = 8 # coarse grouping in x
COARSE_GRID_NY = 8 # coarse grouping in y
GAMMA = None # None -> auto pick ~ E/h ; else set float (e.g. 1e12)
AL_MAXITER = 40
AL_RTOL = 1e-8
TOTAL_FORCE = 100.0 # [N] downward on top z=1 face
OUT_FILE = "tied_contact-mortar.vtu"
# -----------------------------
# Mesh + boundaries
# -----------------------------
mesh_top, mesh_bot = make_default_meshes()
# -----------------------------
# Bulk FE spaces and elasticity
# -----------------------------
elem_u, 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)
K2 = a_elastic.assemble(basis_bot)
n1, n2 = K1.shape[0], K2.shape[0]
# -----------------------------
# Mortar assembly on intersection
# Constraint: u_top = u_bot on overlap (all 3 comps)
# Fine multiplier: vector P0 on contact intersection triangles
# -----------------------------
def assemble_mortar_tied(mesh_top, mesh_bot, elem_u):
# trace contact surfaces to 2D tri meshes and intersect in xy
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_u_top = FacetBasis(mesh_top, elem_u, facets=orig1[t1], quadrature=quad1)
fb_u_bot = FacetBasis(mesh_bot, elem_u, facets=orig2[t2], quadrature=quad2)
elem_lam = ElementVector(ElementTetP0())
fb_lam = FacetBasis(mesh_top, elem_lam, facets=orig1[t1], quadrature=quad1)
@skfem.BilinearForm
def b_dot(u, v, w):
return dot(v, u)
B1 = skfem.asm(b_dot, fb_u_top, fb_lam)
B2 = skfem.asm(b_dot, fb_u_bot, fb_lam)
# return also intersection triangle centers in xy for coarse grouping
# m12 is a MeshTri; its element midpoints can be obtained from p,t:
pts = m12.p.T # (npts, 2)
tri = m12.t.T # (nelems, 3)
centers_xy = pts[tri].mean(axis=1) # (nelems, 2)
# For P0 lambda (vector), dofs correspond to elements, 2D tri elements.
# fb_lam.N == 3 * nelems (vector of P0), but B1/B2 are (nl, n_u) with nl = 3 * nelems.
# We'll use centers per element for grouping and later expand to 3 components.
return B1, B2, centers_xy
B1_f, B2_f, centers_xy = assemble_mortar_tied(mesh_top, mesh_bot, elem_u)
nl_f = B1_f.shape[0] # = 3 * n_contact_tris
ntris = centers_xy.shape[0]
# -----------------------------
# Optional coarse mortar: aggregate lambda dofs
# Build P: (nl_f, nl_c) so that c_c = P^T c_f
# Then Bc = P^T B
# -----------------------------
def make_coarse_P_from_xy_centers(centers_xy, nx=8, ny=8):
x = centers_xy[:, 0]
y = centers_xy[:, 1]
xmin, xmax = x.min(), x.max()
ymin, ymax = y.min(), y.max()
# avoid division by zero if degenerate
dx = (xmax - xmin) if (xmax > xmin) else 1.0
dy = (ymax - ymin) if (ymax > ymin) else 1.0
ix = np.clip(((x - xmin) / dx * nx).astype(int), 0, nx - 1)
iy = np.clip(((y - ymin) / dy * ny).astype(int), 0, ny - 1)
gid = ix + nx * iy # element group id (0..nx*ny-1), some groups may be empty
# compress group ids to only used groups
uniq, inv = np.unique(gid, return_inverse=True)
gid_compact = inv
n_groups = uniq.size
# element-wise areas for weights (optional; here uniform)
w = np.ones_like(gid_compact, dtype=float)
# build P_elem: (ntris, n_groups)
rows = np.arange(ntris, dtype=np.int64)
cols = gid_compact.astype(np.int64)
data = w
P_elem = coo_matrix((data, (rows, cols)), shape=(ntris, n_groups)).tocsr()
# expand to vector components (x,y,z): P = kron(I3, P_elem)
# => shape (3*ntris, 3*n_groups)
# We'll do block-diagonal assembly:
blocks = []
for _ in range(3):
blocks.append(P_elem)
P = bmat([[blocks[0], None, None],
[None, blocks[1], None],
[None, None, blocks[2]]], format="csr")
return P
if USE_COARSE_MORTAR:
P = make_coarse_P_from_xy_centers(centers_xy, nx=COARSE_GRID_NX, ny=COARSE_GRID_NY)
B1 = (P.T @ B1_f).tocsr()
B2 = (P.T @ B2_f).tocsr()
else:
P = None
B1, B2 = B1_f, B2_f
nl = B1.shape[0]
# -----------------------------
# Loads: uniform pressure on top z=1 face so total force = TOTAL_FORCE
# -----------------------------
f1, area, pressure = assemble_top_pressure_load(mesh_top, elem_u, TOTAL_FORCE)
f2 = np.zeros(n2)
# -----------------------------
# Dirichlet: fix bottom z=-0.5 face (ALL components)
# -----------------------------
D_bot = basis_bot.get_dofs("dirichlet").all()
D_global_u = D_bot + n1 # shift to global (u2 block in stacked u=[u1;u2])
# -----------------------------
# Solver A: Pure KKT (Lagrange multiplier saddle system)
# -----------------------------
def solve_kkt(K1, K2, B1, B2, f1, f2, D_global_u, Mlam=None):
F = np.hstack([f1, f2, np.zeros(B1.shape[0])])
K = bmat([
[K1, None, B1.T],
[None, K2, -B2.T],
[B1, -B2, Mlam],
], format="csr")
# Condense ONLY displacement Dirichlet dofs
Kc, Fc, uc, I = skfem.condense(K, F, D=D_global_u)
uc, info = minres(Kc, Fc, rtol=1e-10, maxiter=20000)
if info != 0:
print("[KKT] MINRES info =", info)
x = np.zeros(K.shape[0])
x[I] = uc
return x
# -----------------------------
# Solver B: Augmented Lagrangian (AL) outer loop
# Minimize energy + <lam, c(u)> + (gamma/2)||c(u)||^2
# where c(u)=B1 u1 - B2 u2 (on lambda dofs)
# -----------------------------
def estimate_h_contact(mesh_top):
# crude h estimate from volume mesh spacing (works OK for init_tensor)
# use average edge in z near contact:
# take unique z coordinates and find minimum spacing
z = np.unique(mesh_top.p[2, :])
if z.size >= 2:
hz = np.min(np.diff(np.sort(z)))
else:
hz = 1.0
return hz
def solve_al(K1, K2, B1, B2, f1, f2, D_global_u,
E, gamma=None, max_iter=40, rtol=1e-8):
n1, n2 = K1.shape[0], K2.shape[0]
nl = B1.shape[0]
f = np.hstack([f1, f2])
# pick gamma ~ E / h
if gamma is None:
h = estimate_h_contact(mesh_top)
gamma = 10.0 * E / max(h, 1e-12) # factor 10 as a starter
gamma = float(gamma)
print(f"[AL] gamma = {gamma:.3e}")
# build Kuu + gamma * C^T C
CtC_11 = (B1.T @ B1).tocsr()
CtC_22 = (B2.T @ B2).tocsr()
CtC_12 = (-(B1.T @ B2)).tocsr()
CtC_21 = (-(B2.T @ B1)).tocsr()
Kuu = bmat([
[K1 + gamma * CtC_11, gamma * CtC_12],
[gamma * CtC_21, K2 + gamma * CtC_22 ],
], format="csr")
lam = np.zeros(nl)
u = np.zeros(n1 + n2)
for it in range(max_iter):
# rhs = f - C^T lam, with C=[B1, -B2]
rhs = f - np.hstack([B1.T @ lam, -B2.T @ lam])
Kc, Fc, uc, I = skfem.condense(Kuu, rhs, D=D_global_u)
uc, info = minres(Kc, Fc, rtol=1e-10, maxiter=20000)
if info != 0:
print("[AL] MINRES info =", info)
u[:] = 0.0
u[I] = uc
u1 = u[:n1]
u2 = u[n1:]
c = (B1 @ u1) - (B2 @ u2)
lam = lam + gamma * c
rel = np.linalg.norm(c) / (np.linalg.norm(B1 @ u1) + 1e-30)
print(f"[AL] it={it:02d} ||c||={np.linalg.norm(c):.3e} rel={rel:.3e}")
if rel < rtol:
break
# pack "full" unknown vector compatible with KKT layout: [u1; u2; lam]
x = np.hstack([u[:n1], u[n1:], lam])
return x
# -----------------------------
# Run solve
# -----------------------------
if SOLVER.upper() == "KKT":
# Optional small stabilization on lambda block (usually not needed with P0, but can help numerically)
Mlam = None
x = solve_kkt(K1, K2, B1, B2, f1, f2, D_global_u, Mlam=Mlam)
elif SOLVER.upper() == "AL":
x = solve_al(K1, K2, B1, B2, f1, f2, D_global_u,
E=E, gamma=GAMMA, max_iter=AL_MAXITER, rtol=AL_RTOL)
else:
raise ValueError("SOLVER must be 'KKT' or 'AL'")
# -----------------------------
# Diagnostics (correct slicing)
# -----------------------------
u1 = x[:n1]
u2 = x[n1:n1+n2]
lam_vec = x[n1+n2:n1+n2+nl]
print("---- Diagnostics ----")
print("‖B1 u1 - B2 u2‖ =", np.linalg.norm(B1 @ u1 - B2 @ u2))
print("‖K1 u1 + B1^T λ - f1‖ =", np.linalg.norm(K1 @ u1 + B1.T @ lam_vec - f1))
print("‖K2 u2 - B2^T λ - f2‖ =", np.linalg.norm(K2 @ u2 - B2.T @ lam_vec - f2))
# -----------------------------
# PyVista export (MultiBlock)
# -----------------------------
def skfem_tet_to_pyvista_grid(mesh: "skfem.MeshTet"):
import pyvista as pv
points = mesh.p.T
cells_tet = mesh.t.T.astype(np.int64)
n = cells_tet.shape[0]
cells = np.hstack([np.full((n, 1), 4, dtype=np.int64), cells_tet]).ravel()
celltypes = np.full(n, pv.CellType.TETRA, dtype=np.uint8)
return pv.UnstructuredGrid(cells, celltypes, points)
def export_pyvista_multiblock(mesh_top, mesh_bot, x, n1, out_path="tied_contact.vtm"):
try:
import pyvista as pv
except ImportError:
# Fallback if PyVista is unavailable: export a single combined VTU via meshio.
import meshio
n_top_nodes = mesh_top.p.shape[1]
n_bot_nodes = 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_nodes]))]
u1f = x[:n1].reshape(n_top_nodes, 3)
u2f = x[n1:n1 + 3 * n_bot_nodes].reshape(n_bot_nodes, 3)
u_all = np.vstack([u1f, u2f])
node_tag = np.zeros(n_top_nodes + n_bot_nodes, 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_nodes
bot_nodes_dir = np.unique(mesh_bot.facets[:, mesh_bot.boundaries["dirichlet"]]) + n_top_nodes
top_nodes_force = np.unique(mesh_top.facets[:, mesh_top.boundaries["force"]])
node_tag[np.hstack([top_nodes_contact, bot_nodes_contact])] = 1
node_tag[bot_nodes_dir] = 2
node_tag[top_nodes_force] = 3
alt_path = out_path.rsplit(".", 1)[0] + ".vtu"
meshio.Mesh(
points=points,
cells=cells,
point_data={"u": u_all, "node_tag": node_tag},
).write(alt_path)
print("wrote", alt_path, "(fallback: meshio)")
return
n_top_nodes = mesh_top.p.shape[1]
n_bot_nodes = mesh_bot.p.shape[1]
ndof_u1 = 3 * n_top_nodes
ndof_u2 = 3 * n_bot_nodes
u1 = x[:ndof_u1].reshape(n_top_nodes, 3)
u2 = x[n1:n1+ndof_u2].reshape(n_bot_nodes, 3)
g1 = skfem_tet_to_pyvista_grid(mesh_top)
g2 = skfem_tet_to_pyvista_grid(mesh_bot)
g1.point_data["u"] = u1
g2.point_data["u"] = u2
g1.point_data["part"] = np.zeros(n_top_nodes, dtype=np.int32)
g2.point_data["part"] = np.ones(n_bot_nodes, dtype=np.int32)
# mark some boundary nodes for convenience (optional)
top_nodes_contact = np.unique(mesh_top.facets[:, mesh_top.boundaries["contact"]])
bot_nodes_contact = np.unique(mesh_bot.facets[:, mesh_bot.boundaries["contact"]])
top_nodes_force = np.unique(mesh_top.facets[:, mesh_top.boundaries["force"]])
bot_nodes_dir = np.unique(mesh_bot.facets[:, mesh_bot.boundaries["dirichlet"]])
tag1 = np.zeros(n_top_nodes, dtype=np.int32)
tag2 = np.zeros(n_bot_nodes, dtype=np.int32)
tag1[top_nodes_contact] = 1
tag2[bot_nodes_contact] = 1
tag1[top_nodes_force] = 3
tag2[bot_nodes_dir] = 2
g1.point_data["node_tag"] = tag1
g2.point_data["node_tag"] = tag2
mb = pv.MultiBlock([g1, g2])
mb.save(out_path)
print("✅ wrote", out_path)
export_pyvista_multiblock(mesh_top, mesh_bot, x, n1, OUT_FILE)
How we should regard the result
The results look physically reasonable overall.
- Nitsche
- Total applied load (z): -100.0 N
- Reaction on bottom Dirichlet boundary (z): +99.99999999998681 N
- Force-balance error: 1.3e-11 (essentially zero)
- Mortar
- Total applied load (z): -100.0 N
- Reaction on bottom Dirichlet boundary (z): +99.99998819958344 N
- Force-balance error: 1.18e-05 (relative 1.18e-07)
- Tied constraint residual ||B1 u1 - B2 u2|| = 5.29e-13 (very small)