Third Medium ContactĀ¶
InĀ [1]:
import matplotlib.pyplot as plt
import torch
from matplotlib.colors import ListedColormap
from torchfem import Planar
from torchfem.elements import linear_to_quadratic
from torchfem.materials import HyperelasticPlaneStrain
from torchfem.mesh import rect_quad
torch.set_default_dtype(torch.float64)
Material modelĀ¶
Compressible neo-Hookean energy with LamƩ parameters $\mu$ and $\lambda$ from $E$ and $\nu$
$$ \psi = \frac{\mu}{2}(\mathrm{tr}\,\mathbf{C} - 3) - \mu \ln J + \frac{\lambda}{2} \ln^2 J. $$
Since $\ln J = \frac{1}{2} \ln \det \mathbf{C}$ is blind to the sign of $J$, $\psi$ returns NaN for inverted elements, so the solver cuts back instead of accepting them.
InĀ [2]:
# Material parameters
E = 1000.0
nu = 0.4
lbd = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))
mu = E / (2.0 * (1.0 + nu))
# Energy function
def psi(F, params):
"""Compressible neo-Hookean energy that becomes NaN once an element inverts."""
mu = params[..., 0]
lam = params[..., 1]
C = F.transpose(-1, -2) @ F
lnJ = 0.5 * torch.logdet(C)
trC = C.diagonal(dim1=-2, dim2=-1).sum(-1)
W = 0.5 * mu * (trC - 3) - mu * lnJ + 0.5 * lam * lnJ**2
return torch.where(torch.linalg.det(F) > 0, 1.0, torch.nan) * W
# Bulk material: params = [mu, lbd]
material = HyperelasticPlaneStrain(psi, params=[mu, lbd])
GeometryĀ¶
A C-shaped cantilever of length $L$ and height $H$ with arms of unit thickness is clamped on the left and pulled down at its top tip. The gap and two element columns beyond its mouth form the third medium.
InĀ [3]:
# Dimensions
L = 10.0
H = 5.0
# Create mesh
nodes, elements = rect_quad(33, 16, 32 / 3, H)
nodes, elements = linear_to_quadratic(nodes, elements)
# Create model
cantilever = Planar(nodes, elements, material)
# Create sets
top_tip = torch.isclose(nodes[:, 0], torch.tensor(L)) & (nodes[:, 1] == H)
bottom_tip = torch.isclose(nodes[:, 0], torch.tensor(L)) & (nodes[:, 1] == 1.0)
ec = cantilever.nodes[elements].mean(dim=1)
cutout = (ec[:, 0] > 1.0) & (ec[:, 1] > 1.0) & (ec[:, 1] < 4.0) | (ec[:, 0] > L)
# Load at top tip
cantilever.constraints[top_tip, 1] = True
cantilever.displacements[top_tip, 1] = -8.0
# Constrained displacement at left end
left = nodes[:, 0] == 0.0
cantilever.constraints[left, :] = True
# Show part in light blue and third medium in white
cmap = ListedColormap(["lightblue", "white"])
cantilever.plot(element_property=cutout, cmap=cmap)
Third mediumĀ¶
The gap and the space beyond the mouth are filled with the bulk material scaled by $\gamma$. This third medium transmits forces only once it is compressed to almost zero volume, but being so soft, it distorts excessively before. Quadratic elements therefore carry the regularization energy (Bluhm et al. 2021)
$$ \frac{1}{2} k_r \nabla\nabla\mathbf{u} \mathbin{\vdots} \nabla\nabla\mathbf{u} \quad \text{with} \quad k_r = \alpha L^2 (\lambda + 2\mu), $$
scaled by the P-wave modulus of the bulk, for which $\alpha \in [10^{-7}, 10^{-5}]$ is a reasonable range (Frederiksen et al. 2025). Larger values help the Newton-Raphson iterations but transmit forces before contact.
InĀ [4]:
# Stiffness scaling of the third medium
gamma = 5e-6
# Regularization of the third medium
alpha = 1e-6
k_r = alpha * L**2 * (lbd + 2.0 * mu)
# Scale the bulk material in the third medium
cantilever.material.params[cutout] = gamma * torch.tensor([mu, lbd])
SolveĀ¶
InĀ [5]:
increments = torch.linspace(0.0, 1.0, 151)
# Solve
u, f, Ļ, F, Ī± = cantilever.solve(
increments=increments,
hessian_modulus=k_r * cutout,
return_intermediate=True,
verbose=True,
)
āāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāā
model Planar Ā· 480 elem Ā· 3,070 dof Ā· float64
machine x86_64 Ā· 2 threads Ā· 16 GB RAM
solver direct Ā· scipy Ā· cpu
newton rtol 1e-08 Ā· atol 1e-06 Ā· <=10 it Ā· finite strain Ā· regularized
āāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāā
Increment Load factor Steps Newton Krylov Residual Wall time
1 0.006667 1 7 - 2.05e-12 0.18 s
2 0.01333 1 2 - 2.22e-12 0.07 s
3 0.02 1 2 - 2.31e-12 0.07 s
4 0.02667 1 2 - 2.07e-12 0.07 s
5 0.03333 1 2 - 2.12e-12 0.07 s
6 0.04 1 2 - 2.00e-12 0.07 s
7 0.04667 1 2 - 2.05e-12 0.07 s
8 0.05333 1 2 - 1.98e-12 0.07 s
9 0.06 1 2 - 2.06e-12 0.07 s
10 0.06667 1 2 - 2.01e-12 0.07 s
11 0.07333 1 2 - 2.00e-12 0.07 s
12 0.08 1 2 - 2.05e-12 0.07 s
13 0.08667 1 2 - 2.29e-12 0.07 s
14 0.09333 1 2 - 2.37e-12 0.07 s
15 0.1 1 2 - 2.55e-12 0.07 s
16 0.1067 1 2 - 2.91e-12 0.07 s
17 0.1133 1 2 - 3.24e-12 0.07 s
18 0.12 1 2 - 3.66e-12 0.07 s
19 0.1267 1 2 - 4.24e-12 0.07 s
20 0.1333 1 2 - 4.80e-12 0.07 s
21 0.14 1 2 - 5.42e-12 0.07 s
22 0.1467 1 2 - 6.27e-12 0.07 s
23 0.1533 1 2 - 6.98e-12 0.07 s
24 0.16 1 2 - 8.02e-12 0.07 s
25 0.1667 1 2 - 8.99e-12 0.07 s
26 0.1733 1 2 - 1.01e-11 0.07 s
27 0.18 1 2 - 1.13e-11 0.07 s
28 0.1867 1 2 - 1.26e-11 0.07 s
29 0.1933 1 2 - 1.41e-11 0.07 s
30 0.2 1 2 - 1.58e-11 0.07 s
31 0.2067 1 2 - 1.75e-11 0.07 s
32 0.2133 1 2 - 1.95e-11 0.07 s
33 0.22 1 2 - 2.17e-11 0.07 s
34 0.2267 1 2 - 2.43e-11 0.07 s
35 0.2333 1 2 - 2.72e-11 0.07 s
36 0.24 1 2 - 3.08e-11 0.07 s
37 0.2467 1 2 - 3.51e-11 0.07 s
38 0.2533 1 2 - 4.06e-11 0.07 s
39 0.26 1 2 - 4.77e-11 0.07 s
40 0.2667 1 2 - 5.74e-11 0.07 s
41 0.2733 1 2 - 7.06e-11 0.07 s
42 0.28 1 2 - 8.93e-11 0.06 s
43 0.2867 1 2 - 1.16e-10 0.07 s
44 0.2933 1 2 - 1.55e-10 0.07 s
45 0.3 1 2 - 2.11e-10 0.07 s
46 0.3067 1 2 - 2.96e-10 0.07 s
47 0.3133 1 2 - 4.21e-10 0.07 s
48 0.32 1 2 - 6.11e-10 0.07 s
49 0.3267 1 2 - 8.96e-10 0.07 s
50 0.3333 1 2 - 1.33e-09 0.07 s
51 0.34 1 2 - 1.98e-09 0.07 s
52 0.3467 1 2 - 2.97e-09 0.07 s
53 0.3533 1 2 - 4.44e-09 0.07 s
54 0.36 1 2 - 6.60e-09 0.07 s
55 0.3667 1 2 - 9.73e-09 0.11 s
56 0.3733 1 2 - 1.42e-08 0.07 s
57 0.38 1 2 - 2.03e-08 0.06 s
58 0.3867 1 2 - 2.84e-08 0.07 s
59 0.3933 1 2 - 3.90e-08 0.07 s
60 0.4 1 2 - 5.22e-08 0.06 s
61 0.4067 1 2 - 6.81e-08 0.07 s
62 0.4133 1 2 - 8.65e-08 0.07 s
63 0.42 1 2 - 1.07e-07 0.07 s
64 0.4267 1 2 - 1.29e-07 0.11 s
65 0.4333 1 2 - 1.51e-07 0.07 s
66 0.44 1 2 - 1.72e-07 0.07 s
67 0.4467 1 2 - 1.91e-07 0.07 s
68 0.4533 1 2 - 2.05e-07 0.06 s
69 0.46 1 2 - 2.15e-07 0.07 s
70 0.4667 1 2 - 2.21e-07 0.06 s
71 0.4733 1 2 - 2.21e-07 0.07 s
72 0.48 1 2 - 2.18e-07 0.06 s
73 0.4867 1 2 - 2.12e-07 0.08 s
74 0.4933 1 2 - 2.05e-07 0.07 s
75 0.5 1 2 - 1.98e-07 0.07 s
76 0.5067 1 2 - 1.91e-07 0.07 s
77 0.5133 1 2 - 1.85e-07 0.07 s
78 0.52 1 2 - 1.81e-07 0.07 s
79 0.5267 1 2 - 1.79e-07 0.06 s
80 0.5333 1 2 - 1.80e-07 0.08 s
81 0.54 1 2 - 1.82e-07 0.07 s
82 0.5467 1 2 - 1.88e-07 0.07 s
83 0.5533 1 2 - 1.97e-07 0.07 s
84 0.56 1 2 - 2.08e-07 0.07 s
85 0.5667 1 2 - 2.23e-07 0.07 s
86 0.5733 1 2 - 2.42e-07 0.07 s
87 0.58 1 2 - 2.65e-07 0.07 s
88 0.5867 1 2 - 2.94e-07 0.07 s
89 0.5933 1 2 - 3.30e-07 0.07 s
90 0.6 1 2 - 3.76e-07 0.07 s
91 0.6067 1 2 - 4.36e-07 0.07 s
92 0.6133 1 2 - 5.14e-07 0.07 s
93 0.62 1 2 - 6.16e-07 0.07 s
94 0.6267 1 2 - 7.51e-07 0.07 s
95 0.6333 1 2 - 9.27e-07 0.07 s
96 0.64 1 3 - 4.28e-12 0.09 s
97 0.6467 1 3 - 5.00e-12 0.09 s
98 0.6533 1 3 - 3.34e-12 0.09 s
99 0.66 1 3 - 6.54e-12 0.09 s
100 0.6667 1 3 - 3.85e-12 0.09 s
101 0.6733 1 3 - 9.97e-12 0.09 s
102 0.68 1 3 - 1.32e-11 0.10 s
103 0.6867 1 3 - 3.17e-11 0.09 s
104 0.6933 1 3 - 7.58e-09 0.09 s
105 0.7 1 6 - 1.29e-10 0.16 s
106 0.7067 3 9 - 5.87e-11 0.31 s ā1 ā2
107 0.7133 2 6 - 3.35e-07 0.18 s ā2
108 0.72 2 9 - 2.43e-07 0.26 s ā2
109 0.7267 4 12 - 1.11e-07 0.34 s ā1 ā3
110 0.7333 5 14 - 1.80e-09 0.41 s ā1 ā4
111 0.74 3 14 - 4.13e-08 0.38 s ā3
112 0.7467 5 14 - 1.21e-09 0.41 s ā1 ā4
113 0.7533 3 14 - 3.47e-08 0.38 s ā3
114 0.76 5 14 - 9.14e-10 0.42 s ā1 ā4
115 0.7667 3 13 - 2.71e-08 0.41 s ā3
116 0.7733 5 14 - 6.03e-10 0.42 s ā1 ā4
117 0.78 3 13 - 1.87e-08 0.36 s ā3
118 0.7867 5 14 - 3.27e-10 0.41 s ā1 ā4
119 0.7933 3 12 - 1.12e-08 0.36 s ā3
120 0.8 2 10 - 7.08e-11 0.27 s ā2
121 0.8067 2 8 - 6.37e-08 0.22 s ā2
122 0.8133 4 12 - 1.06e-08 0.36 s ā1 ā3
123 0.82 2 8 - 6.64e-07 0.22 s ā2
124 0.8267 2 8 - 5.79e-08 0.22 s ā2
125 0.8333 4 14 - 8.87e-11 0.39 s ā1 ā3
126 0.84 6 18 - 5.00e-07 0.55 s ā6
127 0.8467 4 17 - 6.04e-07 0.54 s ā4
128 0.8533 5 23 - 7.16e-08 0.62 s ā1 ā4
129 0.86 4 13 - 1.96e-07 0.43 s ā4
130 0.8667 3 12 - 1.56e-10 0.34 s ā3
131 0.8733 2 7 - 4.03e-07 0.20 s ā2
132 0.88 2 9 - 2.33e-10 0.24 s ā2
133 0.8867 4 18 - 3.38e-10 0.51 s ā1 ā3
134 0.8933 9 19 - 2.61e-07 0.62 s ā9
135 0.9 5 16 - 2.22e-10 0.50 s ā5
136 0.9067 3 10 - 3.72e-10 0.29 s ā3
137 0.9133 3 11 - 4.18e-10 0.32 s ā3
138 0.92 2 9 - 2.19e-09 0.26 s ā2
139 0.9267 2 11 - 5.68e-08 0.35 s ā2
140 0.9333 4 27 - 7.93e-10 0.76 s ā1 ā3
141 0.94 8 34 - 2.78e-07 0.95 s ā8
142 0.9467 5 15 - 1.96e-07 0.45 s ā5
143 0.9533 3 10 - 8.36e-10 0.30 s ā3
144 0.96 3 10 - 6.94e-10 0.31 s ā3
145 0.9667 2 7 - 2.11e-08 0.20 s ā2
146 0.9733 2 8 - 7.22e-10 0.22 s ā2
147 0.98 4 19 - 1.32e-09 0.52 s ā1 ā3
148 0.9867 7 16 - 3.18e-08 0.51 s ā7
149 0.9933 4 12 - 9.17e-10 0.36 s ā4
150 1 3 9 - 4.74e-09 0.27 s ā3 āāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāā converged Ā· 150 increments Ā· 820 Newton iterations Ā· 24.87 s
ResultsĀ¶
The reference is a node-to-segment contact solution of the same problem (Faltus et al. 2026, Fig. 1b).
InĀ [6]:
# Plot
cantilever.plot(u[-1], element_property=cutout, cmap=cmap)
InĀ [7]:
reference = torch.tensor(
[
[0.0, 0.0],
[3.16, 0.0],
[5.02, 2.24],
[6.19, 3.90],
[7.66, 6.18],
[8.26, 6.92],
[8.98, 7.51],
[9.66, 7.83],
]
)
plt.plot(reference[:, 0], reference[:, 1], color="black", label="Reference")
plt.plot(-u[:, top_tip, 1], -u[:, bottom_tip, 1], color="deeppink", label="torch-fem")
plt.xlabel("Top Tip Displacement")
plt.ylabel("Bottom Tip Displacement")
plt.grid()
plt.show()
