Coupling PDEs of multiple dimensions#
In the previous section, we considered problems where we solved non-coupled physics on two parts of a domain. However, the problems we actually want to considered are problems that are coupled across domains.
In this section we will cover how to couple PDEs formulated in the domain (2D) with those living on a subset of facets (1D). This can naturally be extended to 3D.
In this section, we will show how to solve the Signorini problem
where \(\mathbf{u}\) is the displacement, \(C\) the stiffness tensor, \(\epsilon\) the symmetric strain tensor and \(\mathbf{f}\) the body force.
In this tutorial we will consider a half circle, where we apply a displacement on the top boundary, and let the curved boundary be a potential contact boundary. We define a rigid surface as a plane at \(y=-h\), where \(h\in \mathbb{R}\)
As seen in Meshes from external sources, we can use GMSH to create such a geometry.
# We inspect the generated mesh and markers
mesh_data = dolfinx.io.gmsh.model_to_mesh(gmsh.model, mesh_comm, gmsh_model_rank, gdim=2)
omega = mesh_data.mesh
ct = mesh_data.cell_tags
ft = mesh_data.facet_tags
gmsh.finalize()
2025-10-23 15:13:12.470 ( 0.956s) [ 7F5270175140]vtkXOpenGLRenderWindow.:1458 WARN| bad X server connection. DISPLAY=
Variational formulation#
We will use a formulation of this problem based on [KS24] and [DFK+25]. We phrase this problem as a minimization problem, where we seek to find the displacement \(\mathbf{u}\) that minimizes the functional
where $\( \mathcal{K} = \{ \mathbf{u}\in V_{\mathbf{u}_D} \vert \mathbf{u}\cdot \hat{\mathbf{n}} \leq g \} \)$
With this re-formulation, we can write a mixed finite element method, where we use two variables, the displacement \(\mathbf{u}\in V(\Omega)\) and a latent variable \(\psi \in Q(\Gamma)\)/
Co-dimension 1 problem
We note that since \(\Gamma\) is the whole curved boundary, we need to solve a mixed dimensional finite element problem.
We choose \(V\) to be a P-th order vector Lagrange space, while \(Q\) is a Pth order scalar Lagrange field.
We start by defining the sub-mesh using the read in mesh markers using dolfinx.mesh.create_submesh().
tdim = omega.topology.dim
fdim = tdim - 1
gdim = omega.geometry.dim
gamma, gamma_to_omega = dolfinx.mesh.create_submesh(omega, fdim, ft.find(potential_contact_marker))[0:2]
dolfinx.mesh.create_submesh() also returns a mapping from the sub-mesh to the parent mesh, as a bi-directional
dolfinx.mesh.EntityMap.
Next, we define the function spaces, and combine them in a block structure
using ufl.MixedFunctionSpace
V = dolfinx.fem.functionspace(omega, ("Lagrange", 1, (omega.geometry.dim,)))
Q = dolfinx.fem.functionspace(gamma, ("Lagrange", 1))
W = ufl.MixedFunctionSpace(V, Q)
We can write the variational formulation as Given \(\alpha_k\), \(\psi_{k-1}\), solve:
Check for convergence.
Update latent variable \(\psi^{k-1}\), \(\alpha_k\).
Notes on the variational form
It is highly non-linear, as it contains \(e^{\psi}\).
We can recognize the traditional part of the Signorini problem as the terms multiplied by \(\alpha_k\).
One can update \(\alpha_k\) at each iteration, but it is not a requirement for convergence, and it is often sufficient to use a constant value.
We have that for the solved problem \(e^{\psi} = \mathbf{u} \cdot \mathbf{n} - g\), which we are guaranteed that this satisfies the Signorini condition.
As discussed in the previous chapter, we need to choose a mesh to integrate over.
As we would like to exploit the definition of the n=ufl.FacetNormal(\Omega) in our
variational problem, we choose the integration domain to be \(\Omega\).
dx = ufl.Measure("dx", domain=omega)
ds = ufl.Measure("ds", domain=omega, subdomain_data=ft, subdomain_id=potential_contact_marker)
Note
Integration over \(\Gamma\)
Note that we have restricted the ds integration measure to the boundary \(\Gamma\),
where \(Q\) is defined.
Next, we define some problem specific parameters
E = dolfinx.fem.Constant(omega, dolfinx.default_scalar_type(2e4))
nu = dolfinx.fem.Constant(omega, 0.4)
alpha = dolfinx.fem.Constant(omega, dolfinx.default_scalar_type(0.1))
As we define the rigid surface as a plane at \(y=-h\), we can define the gap between the undeformed geometry with coordinates (x, y) and the surface as
h = 0.13
x, y = ufl.SpatialCoordinate(omega)
g = y + dolfinx.fem.Constant(omega, dolfinx.default_scalar_type(h))
uD = 0.72
Similarly, we have that \(\hat n = (0, -1)\)
n = ufl.FacetNormal(omega)
n_g = dolfinx.fem.Constant(omega, np.zeros(gdim, dtype=dolfinx.default_scalar_type))
n_g.value[-1] = -1
f = dolfinx.fem.Constant(omega, np.zeros(gdim, dtype=dolfinx.default_scalar_type))
We can write the residual of the variational form
v, w = ufl.TestFunctions(W)
u = dolfinx.fem.Function(V, name="Displacement")
psi = dolfinx.fem.Function(Q, name="Latent_variable")
psi_k = dolfinx.fem.Function(Q, name="Previous_Latent_variable")
mu = E / (2.0 * (1.0 + nu))
lmbda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))
def epsilon(w):
return ufl.sym(ufl.grad(w))
def sigma(w, mu, lmbda):
ew = epsilon(w)
gdim = ew.ufl_shape[0]
return 2.0 * mu * epsilon(w) + lmbda * ufl.tr(ufl.grad(w)) * ufl.Identity(gdim)
F = alpha * ufl.inner(sigma(u, mu, lmbda), epsilon(v)) * dx
F -= alpha * ufl.inner(f, v) * dx
F += -ufl.inner(psi - psi_k, ufl.dot(v, n)) * ds
F += ufl.inner(ufl.dot(u, n_g), w) * ds
F += ufl.inner(ufl.exp(psi), w) * ds - ufl.inner(g, w) * ds
residual = ufl.extract_blocks(F)
Similarly, we can write the Jacobian
du, dpsi = ufl.TrialFunctions(W)
jac = ufl.derivative(F, u, du) + ufl.derivative(F, psi, dpsi)
J = ufl.extract_blocks(jac)
Jacobian for mixed function spaces
Note that we differentiate with respect to the function in the respective “sub space”, \(V\) and \(Q\), but use trial functions form \(W\). This is to be able to extract blocks when creating the form for the Jacobian.
We define the displacement on the top boundary as we have done in previous tutorials.
However, as we are using a ufl.MixedFunctionSpace, we can define the boundary condition
with the appropriate sub-space without a mapping.
def disp_func(x):
values = np.zeros((gdim, x.shape[1]), dtype=dolfinx.default_scalar_type())
values[1] = -uD
return values
u_bc = dolfinx.fem.Function(V)
u_bc.interpolate(disp_func)
bc = dolfinx.fem.dirichletbc(u_bc, dolfinx.fem.locate_dofs_topological(V, fdim, ft.find(displacement_marker)))
bcs = [bc]
# We want to consider the Von-Mises stresses in post-processing, and
# use DOLFINx Expression to interpolate the stresses into an appropriate
# function space.
V_DG = dolfinx.fem.functionspace(omega, ("DG", 1, (omega.geometry.dim,)))
stress_space, stress_to_disp = V_DG.sub(0).collapse()
von_mises = dolfinx.fem.Function(stress_space, name="von_Mises")
u_dg = dolfinx.fem.Function(V_DG, name="u")
s = sigma(u, mu, lmbda) - 1.0 / 3 * ufl.tr(sigma(u, mu, lmbda)) * ufl.Identity(len(u))
von_Mises = ufl.sqrt(3.0 / 2 * ufl.inner(s, s))
stress_expr = dolfinx.fem.Expression(von_Mises, stress_space.element.interpolation_points)
# We can now set up the solver and solve the problem
entity_maps = [gamma_to_omega]
solver = dolfinx.fem.petsc.NonlinearProblem(
residual,
u=[u, psi],
J=J,
bcs=bcs,
entity_maps=entity_maps,
petsc_options={
"snes_monitor": None,
"ksp_type": "preonly",
"pc_type": "lu",
"pc_factor_mat_solver_type": "mumps",
"mat_mumps_icntl_14": 120,
"ksp_error_if_not_converged": True,
"snes_error_if_not_converged": True,
},
petsc_options_prefix="signorini_",
)
# Note that all memory is assigned outside the for-lopp.
# In this problem, we measure the norm of the change in the primal space,
# rather than the for the mixed function.
max_iterations = 25
normed_diff = 0
tol = 1e-5
solver.solve()
iterations = solver.solver.getIterationNumber()
print(f"Converged in {iterations} iterations")
# We compute the von-Mises stresses for the final solution
von_mises.interpolate(stress_expr)
# Additionally, we interpolate the displacement onto a DG function space
# for compatible visualization in Pyvista.
u_dg.interpolate(u)
0 SNES Function norm 2.049913151577e+04
1 SNES Function norm 6.659105853306e-02
2 SNES Function norm 2.447919804152e-02
3 SNES Function norm 1.116150597457e-02
4 SNES Function norm 8.281115673245e-03
5 SNES Function norm 2.347664043339e-03
6 SNES Function norm 2.687040995684e-04
7 SNES Function norm 6.078804962746e-05
Converged in 7 iterations
References#
Jørgen S. Dokken, Patrick E. Farrell, Brendan Keith, Ioannis P.A. Papadopoulos, and Thomas M. Surowiec. The latent variable proximal point algorithm for variational problems with inequality constraints. Computer Methods in Applied Mechanics and Engineering, 445:118181, 2025. doi:10.1016/j.cma.2025.118181.
Brendan Keith and Thomas M. Surowiec. Proximal galerkin: a structure-preserving finite element method for pointwise bound constraints. Foundations of Computational Mathematics, Nov 2024. doi:10.1007/s10208-024-09681-8.