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
omega, ct, ft = dolfinx.io.gmshio.model_to_mesh(gmsh.model, mesh_comm, gmsh_model_rank, gdim=2)
gmsh.finalize()
Variational formulation#
We will use a formulation of this problem based on [KS24] and [DPEKS24]. 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
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]
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\).
This means that we have to create a map from each facet in \(\Omega\) to
the corresponding facet in \(\Gamma\).
facet_imap = omega.topology.index_map(ft.dim)
num_facets = facet_imap.size_local + facet_imap.num_ghosts
omega_to_gamma = np.full(num_facets, -1)
omega_to_gamma[gamma_to_omega] = np.arange(len(gamma_to_omega))
entity_maps = {gamma: omega_to_gamma}
# -
Next, we define the integration measures
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 = dolfinx.fem.form(ufl.extract_blocks(F), entity_maps=entity_maps)
Similarly, we can write the Jacobian
du, dpsi = ufl.TrialFunctions(W)
jac = ufl.derivative(F, u, du) + ufl.derivative(F, psi, dpsi)
J = dolfinx.fem.form(ufl.extract_blocks(jac), entity_maps=entity_maps)
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]
To solve this problem, we will use PETSc. We use the following wrapper to solve the problem
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
solver = NewtonSolver(
residual,
J,
[u, psi],
bcs=bcs,
max_iterations=25,
petsc_options={
"ksp_type": "preonly",
"pc_type": "lu",
"pc_factor_mat_solver_type": "mumps",
"mat_mumps_icntl_14": 120,
"ksp_error_if_not_converged": True,
},
error_on_nonconvergence=True,
)
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
u_prev = dolfinx.fem.Function(V)
diff = dolfinx.fem.Function(V)
for it in range(max_iterations):
print(f"{it=}/{max_iterations} {normed_diff:.2e}")
# Solve the first iterations inaccurately
solver_tol = 100*tol if it < 3 else tol
converged = solver.solve(solver_tol, 1)
diff.x.array[:] = u.x.array - u_prev.x.array
diff.x.petsc_vec.normBegin(2)
normed_diff = diff.x.petsc_vec.normEnd(2)
if normed_diff <= tol and it >=3:
print(f"Converged at {it=} with increment norm {normed_diff:.2e}<{tol:.2e}")
break
u_prev.x.array[:] = u.x.array
psi_k.x.array[:] = psi.x.array
if not converged:
print(f"Solver did not convert at {it=}, exiting with {converged=}")
break
it=0/25 0.00e+00
Iteration 1: Correction norm 19.75603520918023
Iteration 2: Correction norm 0.0020541893014397712
Iteration 3: Correction norm 0.0032824903576515814
Iteration 4: Correction norm 0.010567366459077002
Iteration 5: Correction norm 0.7826663660541742
Iteration 6: Correction norm 0.28180190008178135
Iteration 7: Correction norm 0.0003624695477901472
it=1/25 1.91e+01
Iteration 1: Correction norm 0.0029443433068172785
Iteration 2: Correction norm 0.0001222423005396523
it=2/25 2.88e-03
Iteration 1: Correction norm 6.454429854631308e-05
it=3/25 6.45e-05
Iteration 1: Correction norm 1.1308753015128175e-05
Iteration 2: Correction norm 2.0675218317521407e-05
Iteration 3: Correction norm 1.198927556775523e-07
it=4/25 3.19e-05
Iteration 1: Correction norm 2.6758383762843034e-05
Iteration 2: Correction norm 6.107902701558581e-06
it=5/25 2.07e-05
Iteration 1: Correction norm 1.4332322165873925e-05
Iteration 2: Correction norm 2.3843265682231737e-06
it=6/25 1.19e-05
Iteration 1: Correction norm 9.47839641376783e-06
Converged at it=6 with increment norm 9.48e-06<1.00e-05
if it == max_iterations - 1:
print(f"Did not converge within {max_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)
References#
J.S. Dokken, Farrell P.E., B. Keith, and T.M. Surowiec. The latent variable proximal point algorithm for problems with pointwise constraints. 2024.
Brendan Keith and Thomas M. Surowiec. Proximal galerkin: a structure-preserving finite element method for pointwise bound constraints. 2024. URL: https://arxiv.org/abs/2307.12444, arXiv:2307.12444.