Contact between two deformable bodies

kumar · April 10, 2026, 8:18am

Dear all,

I am currently working on a contact problem with friction in FEniCSx. While exploring examples and previous discussions, I noticed that many responses suggest starting with the Hertzian contact example from the FEniCS tutorial:

I have gone through this example, but it models contact using an obstacle equation rather than a physical second body. In my case, I would like to simulate contact between two deformable bodies.

Could anyone suggest resources, examples, or recommended approaches for setting up such a problem in FEniCSx?

Any guidance would be greatly appreciated. Thank you for your time.

kumar · April 10, 2026, 2:03pm

I wanted to add a small follow‑up: unfortunately, I haven’t been able to locate concrete FEniCSx examples that involve two bodies in contact, with or without friction. Most examples—including the Hertzian contact demo—use an obstacle formulation or a rigid indenter, so I’m still unsure about the recommended starting point for a fully contact setup.

If anyone has already explored this, or can point to an example, a minimal research code snippet, or even just confirm the feasibility in FEniCSx, that would really help me understand how to begin structuring such a problem.

Thanks again for your time and any guidance.

nate · April 10, 2026, 5:13pm

Bilateral contact is horrendously difficult. The most advanced implementation I know of with DOLFINx is GitHub - Wells-Group/asimov-contact: Contact models for ASiMoV · GitHub.

dokken · April 10, 2026, 10:41pm

You can also try the third medium contact method, which is for instance covered in: First order finite element formulations for third medium contact | Computational Mechanics | Springer Nature Link
Here is an example implementation:

"""
Author: Jørgen S. Dokken
SPDX-License-Identifier: MIT
Third medium contact example from DOI: 10.1007/s00466-025-02628-y
"""

from mpi4py import MPI
import dolfinx
import dolfinx.fem.petsc
import ufl
import numpy as np

tol = 1.0e-6


def thirdmedium(x):
    return (x[0] >= 0.1 - tol) & (x[2] >= 0.1 - tol) & (x[2] <= 0.6 + tol)


def top(x):
    return np.isclose(x[2], 0.7)


def left(x):
    return np.isclose(x[0], 0.0)


BODY_marker = 1
AIR_marker = 2
mesh = dolfinx.mesh.create_box(
    MPI.COMM_WORLD,
    [[0, 0, 0], [2, 0.7, 0.7]],
    [40, 14, 14],
    ghost_mode=dolfinx.mesh.GhostMode.shared_facet,
)
num_cells_local = (
    mesh.topology.index_map(mesh.topology.dim).size_local
    + mesh.topology.index_map(mesh.topology.dim).num_ghosts
)
markers = np.full(num_cells_local, BODY_marker, dtype=np.int32)
markers[dolfinx.mesh.locate_entities(mesh, mesh.topology.dim, thirdmedium)] = AIR_marker
ct = dolfinx.mesh.meshtags(mesh, mesh.topology.dim, np.arange(num_cells_local), markers)


LEFT_marker = 2
NO_TRACTION_marker = 3
POTENTIAL_CONTACT_marker = 4


mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
num_facets_local = (
    mesh.topology.index_map(mesh.topology.dim - 1).size_local
    + mesh.topology.index_map(mesh.topology.dim - 1).num_ghosts
)
all_body_facets = dolfinx.mesh.compute_incident_entities(
    mesh.topology, ct.find(BODY_marker), mesh.topology.dim, mesh.topology.dim - 1
)
all_air_facets = dolfinx.mesh.compute_incident_entities(
    mesh.topology, ct.find(AIR_marker), mesh.topology.dim, mesh.topology.dim - 1
)
interface_facets = np.intersect1d(all_body_facets, all_air_facets)
all_exterior_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)

facet_markers = np.zeros(num_facets_local, dtype=np.int32)
facet_markers[np.intersect1d(all_exterior_facets, all_body_facets)] = NO_TRACTION_marker
facet_markers[dolfinx.mesh.locate_entities(mesh, mesh.topology.dim - 1, left)] = (
    LEFT_marker
)
facet_markers[interface_facets] = POTENTIAL_CONTACT_marker
f_to_c = mesh.topology.connectivity(mesh.topology.dim - 1, mesh.topology.dim)

ft_indices = np.flatnonzero(facet_markers)

ft = dolfinx.mesh.meshtags(
    mesh, mesh.topology.dim - 1, ft_indices, facet_markers[ft_indices]
)
ft.name = "facet_markers"
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "markers.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_meshtags(ct, mesh.geometry)
    xdmf.write_meshtags(ft, mesh.geometry)


V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1, (mesh.geometry.dim,)))
d = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(mesh.geometry.dim))

third_medium_mesh, medium_map = dolfinx.mesh.create_submesh(
    mesh, mesh.topology.dim, ct.find(AIR_marker)
)[0:2]
P = dolfinx.fem.functionspace(third_medium_mesh, ("Lagrange", 1, (3,)))
W = ufl.MixedFunctionSpace(V, P)
p = dolfinx.fem.Function(P)


u = dolfinx.fem.Function(V)
I = ufl.Identity(len(u))
x = ufl.SpatialCoordinate(mesh)
phi = x + u
F = I + ufl.grad(u)

J = ufl.det(F)
C = F.T * F


# 3. Separate Material Properties
E = 5
nu = 0.3
K = E / (3 * (1 - 2 * nu))
mu = E / (2 * (1 + nu))
K_body = dolfinx.fem.Constant(mesh, K)
mu_body = dolfinx.fem.Constant(mesh, mu)
Psi_body = K_body / 2 * ufl.ln(J) ** 2 + mu_body / 2 * (J ** (-2 / 3) * ufl.tr(C) - 3)

dx = ufl.Measure("dx", domain=mesh, subdomain_data=ct)
dxVol = dx(BODY_marker)
dxThird = dx(AIR_marker)
ds = ufl.Measure("ds", domain=mesh, subdomain_data=ft)

b = dolfinx.fem.Constant(
    mesh, np.zeros(mesh.geometry.dim, dtype=dolfinx.default_scalar_type)
)
t = dolfinx.fem.Constant(
    mesh, np.zeros(mesh.geometry.dim, dtype=dolfinx.default_scalar_type)
)
Pi = (
    Psi_body * dxVol
    - ufl.inner(b, phi) * dxVol
    - ufl.inner(t, phi) * ds((NO_TRACTION_marker))
)


gamma = dolfinx.fem.Constant(mesh, 4.0e-6)
Pi_third = gamma * Psi_body * dxThird
beta_1 = dolfinx.fem.Constant(mesh, 1.0)*gamma
beta_2 = dolfinx.fem.Constant(mesh, 100.0)*gamma
f0_skew = 0.5 * (F[0, 1] - F[1, 0])
f1_skew = 0.5 * (F[0, 2] - F[2, 0])
f2_skew = 0.5 * (F[1, 2] - F[2, 1])
fi_skew = ufl.as_vector((f0_skew, f1_skew, f2_skew))
L_i = np.zeros(3)
for dim in range(mesh.geometry.dim):
    x_i_max = mesh.comm.allreduce(mesh.geometry.x[:, dim].max(), op=MPI.MAX)
    x_i_min = mesh.comm.allreduce(mesh.geometry.x[:, dim].min(), op=MPI.MIN)
    L_i[dim] = x_i_max - x_i_min
d = dolfinx.fem.Constant(mesh, np.max(L_i))

skew_term = fi_skew - 1 / d * p
Pi_fi = (
    beta_1 / 2  * ufl.dot(skew_term, skew_term) + beta_2 / 2 * ufl.inner(ufl.grad(p), ufl.grad(p))
) * dxThird


left_dofs = dolfinx.fem.locate_dofs_topological(
    V, mesh.topology.dim - 1, ft.find(LEFT_marker)
)
bc_left = dolfinx.fem.dirichletbc(
    np.zeros(mesh.geometry.dim, dtype=dolfinx.default_scalar_type), left_dofs, V
)
mesh.topology.create_connectivity(0, mesh.topology.dim)
node = dolfinx.mesh.locate_entities(mesh, 0, lambda x: np.isclose(x[0], 2.0) & np.isclose(x[1], 0.0) & np.isclose(x[2], 0.7))
dofs_point_z = dolfinx.fem.locate_dofs_topological(V.sub(2), 0, node)
dofs_point_y = dolfinx.fem.locate_dofs_topological(V.sub(1), 0, node)
applied_z = dolfinx.fem.Constant(mesh, 0.0)
applied_y = dolfinx.fem.Constant(mesh, 0.0)
bc_point_z = dolfinx.fem.dirichletbc(applied_z, dofs_point_z, V.sub(2))
bc_point_y = dolfinx.fem.dirichletbc(applied_y, dofs_point_y, V.sub(1))
bcs = [bc_left, bc_point_z, bc_point_y]

dup = ufl.TestFunctions(W)
residual = ufl.derivative(Pi + Pi_third + Pi_fi, [u, p], dup)
problem = dolfinx.fem.petsc.NonlinearProblem(
    ufl.extract_blocks(residual),
    [u, p],
    bcs=bcs,
    entity_maps=[medium_map],
    petsc_options_prefix="snes",
    petsc_options={
        "snes_type": "newtonls",
        "snes_monitor": None,
        "snes_converged_reason": None,
        "snes_error_if_not_converged": True,
        "ksp_type": "preonly",
        "pc_type": "lu",
        "pc_factor_mat_solver_type": "mumps",
        "snes_linesearch_type": "bt",
        "snes_atol": 5e-5,
        "snes_rtol": 5e-5,
        "snes_max_it": 100,
    },
)

# Do partial loading
full_disp = np.array([0.0, 0.5, -2])
tm = dolfinx.fem.functionspace(mesh, ("DG", 0))
tm_func = dolfinx.fem.Function(tm, name="medium")
tm_func.x.array[:] = ct.values
with dolfinx.io.VTXWriter(mesh.comm, "solution.bp", [u, tm_func]) as bp:
    loading_steps = 10
    loading = np.linspace(0, 1, loading_steps + 1)[1:]
    for i,load in enumerate(loading):
        applied_z.value = full_disp[2] * load
        applied_y.value = full_disp[1] * load
        problem.solve()
        bp.write(i)

Depending on what kind of elastic model you want for the body, it might be easier to get convergence that it was from the test in Wriggers paper, which I haven’t found a satisfactory linear loading method for (probably need a adaptive loading method).

kumar · April 11, 2026, 6:49am

Thanks for both the responses. This is really helpful to know the first step. I will look into the same.

dokken · April 13, 2026, 4:55am

Tweaking the solver parameters a bit, i got something that seems ok-ish.

"""
Author: Jørgen S. Dokken
SPDX-License-Identifier: MIT
Third medium contact example from DOI: 10.1007/s00466-025-02628-y
"""

from mpi4py import MPI
import dolfinx
import dolfinx.fem.petsc
import ufl
import numpy as np

tol = 1.0e-6


def thirdmedium(x):
    return (x[0] >= 0.1 - tol) & (x[2] >= 0.1 - tol) & (x[2] <= 0.6 + tol)


def top(x):
    return np.isclose(x[2], 0.7)


def left(x):
    return np.isclose(x[0], 0.0)


BODY_marker = 1
AIR_marker = 2
mesh = dolfinx.mesh.create_box(
    MPI.COMM_WORLD,
    [[0, 0, 0], [2, 0.7, 0.7]],
    [40, 14, 14],
    ghost_mode=dolfinx.mesh.GhostMode.shared_facet, #cell_type=dolfinx.mesh.CellType.hexahedron,
)
num_cells_local = (
    mesh.topology.index_map(mesh.topology.dim).size_local
    + mesh.topology.index_map(mesh.topology.dim).num_ghosts
)
markers = np.full(num_cells_local, BODY_marker, dtype=np.int32)
markers[dolfinx.mesh.locate_entities(mesh, mesh.topology.dim, thirdmedium)] = AIR_marker
ct = dolfinx.mesh.meshtags(mesh, mesh.topology.dim, np.arange(num_cells_local), markers)


LEFT_marker = 2
NO_TRACTION_marker = 3
POTENTIAL_CONTACT_marker = 4


mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
num_facets_local = (
    mesh.topology.index_map(mesh.topology.dim - 1).size_local
    + mesh.topology.index_map(mesh.topology.dim - 1).num_ghosts
)
all_body_facets = dolfinx.mesh.compute_incident_entities(
    mesh.topology, ct.find(BODY_marker), mesh.topology.dim, mesh.topology.dim - 1
)
all_air_facets = dolfinx.mesh.compute_incident_entities(
    mesh.topology, ct.find(AIR_marker), mesh.topology.dim, mesh.topology.dim - 1
)
interface_facets = np.intersect1d(all_body_facets, all_air_facets)
all_exterior_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)

facet_markers = np.zeros(num_facets_local, dtype=np.int32)
facet_markers[np.intersect1d(all_exterior_facets, all_body_facets)] = NO_TRACTION_marker
facet_markers[dolfinx.mesh.locate_entities(mesh, mesh.topology.dim - 1, left)] = (
    LEFT_marker
)
facet_markers[interface_facets] = POTENTIAL_CONTACT_marker
f_to_c = mesh.topology.connectivity(mesh.topology.dim - 1, mesh.topology.dim)

ft_indices = np.flatnonzero(facet_markers)

ft = dolfinx.mesh.meshtags(
    mesh, mesh.topology.dim - 1, ft_indices, facet_markers[ft_indices]
)
ft.name = "facet_markers"
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "markers.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_meshtags(ct, mesh.geometry)
    xdmf.write_meshtags(ft, mesh.geometry)


V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1, (mesh.geometry.dim,)))
d = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(mesh.geometry.dim))

third_medium_mesh, medium_map = dolfinx.mesh.create_submesh(
    mesh, mesh.topology.dim, ct.find(AIR_marker)
)[0:2]
P = dolfinx.fem.functionspace(third_medium_mesh, ("Lagrange", 1, (3,)))
W = ufl.MixedFunctionSpace(V, P)
p = dolfinx.fem.Function(P)


u = dolfinx.fem.Function(V)
I = ufl.Identity(len(u))
x = ufl.SpatialCoordinate(mesh)
phi = x + u
F = I + ufl.grad(u)

J = ufl.det(F)
C = F.T * F


# 3. Separate Material Properties
E = 5
nu = 0.3
K = E / (3 * (1 - 2 * nu))
mu = E / (2 * (1 + nu))
K_body = dolfinx.fem.Constant(mesh, K)
mu_body = dolfinx.fem.Constant(mesh, mu)
Psi_body = K_body / 2 * ufl.ln(J) ** 2 + mu_body / 2 * (J ** (-2 / 3) * ufl.tr(C) - 3)

dx = ufl.Measure("dx", domain=mesh, subdomain_data=ct)
dxVol = dx(BODY_marker)
dxThird = dx(AIR_marker)
ds = ufl.Measure("ds", domain=mesh, subdomain_data=ft)

b = dolfinx.fem.Constant(
    mesh, np.zeros(mesh.geometry.dim, dtype=dolfinx.default_scalar_type)
)
t = dolfinx.fem.Constant(
    mesh, np.zeros(mesh.geometry.dim, dtype=dolfinx.default_scalar_type)
)
Pi = (
    Psi_body * dxVol
    - ufl.inner(b, phi) * dxVol
    - ufl.inner(t, phi) * ds((NO_TRACTION_marker))
)


gamma = dolfinx.fem.Constant(mesh, 2.0e-6)
Pi_third = gamma * Psi_body * dxThird
beta_1 = dolfinx.fem.Constant(mesh, 1.0e-1)
beta_2 = dolfinx.fem.Constant(mesh, 1.0e-3)
f0_skew = 0.5 * (F[0, 1] - F[1, 0])
f1_skew = 0.5 * (F[0, 2] - F[2, 0])
f2_skew = 0.5 * (F[1, 2] - F[2, 1])
fi_skew = ufl.as_vector((f0_skew, f1_skew, f2_skew))
L_i = np.zeros(3)
for dim in range(mesh.geometry.dim):
    x_i_max = mesh.comm.allreduce(mesh.geometry.x[:, dim].max(), op=MPI.MAX)
    x_i_min = mesh.comm.allreduce(mesh.geometry.x[:, dim].min(), op=MPI.MIN)
    L_i[dim] = x_i_max - x_i_min
d = dolfinx.fem.Constant(mesh, np.max(L_i))

skew_term = fi_skew - 1 / d * p
Pi_fi = (
    beta_1 / 2  * ufl.dot(skew_term, skew_term) + beta_2 / 2 * ufl.inner(ufl.grad(p), ufl.grad(p))
) * dxThird


left_dofs = dolfinx.fem.locate_dofs_topological(
    V, mesh.topology.dim - 1, ft.find(LEFT_marker)
)
bc_left = dolfinx.fem.dirichletbc(
    np.zeros(mesh.geometry.dim, dtype=dolfinx.default_scalar_type), left_dofs, V
)
mesh.topology.create_connectivity(0, mesh.topology.dim)
node = dolfinx.mesh.locate_entities(mesh, 0, lambda x: np.isclose(x[0], 2.0) & np.isclose(x[1], 0.0) & np.isclose(x[2], 0.7))
dofs_point_z = dolfinx.fem.locate_dofs_topological(V.sub(2), 0, node)
dofs_point_y = dolfinx.fem.locate_dofs_topological(V.sub(1), 0, node)
applied_z = dolfinx.fem.Constant(mesh, 0.0)
applied_y = dolfinx.fem.Constant(mesh, 0.0)
bc_point_z = dolfinx.fem.dirichletbc(applied_z, dofs_point_z, V.sub(2))
bc_point_y = dolfinx.fem.dirichletbc(applied_y, dofs_point_y, V.sub(1))
bcs = [bc_left, bc_point_z, bc_point_y]

dup = ufl.TestFunctions(W)
residual = ufl.derivative(Pi + Pi_third + Pi_fi, [u, p], dup)
problem = dolfinx.fem.petsc.NonlinearProblem(
    ufl.extract_blocks(residual),
    [u, p],
    bcs=bcs,
    entity_maps=[medium_map],
    petsc_options_prefix="snes",
    petsc_options={
        "snes_type": "newtonls",
        "snes_monitor": None,
        "snes_converged_reason": None,
        "snes_error_if_not_converged": False,
        "ksp_type": "preonly",
        "pc_type": "lu",
        "pc_factor_mat_solver_type": "mumps",
        "snes_linesearch_type": "bt",
        "snes_atol": 1e-6,
        "snes_rtol": 1e-6,
        "snes_max_it": 25,
    },
)

# Do partial loading
full_disp = np.array([0.0, 0.5, -2])
tm = dolfinx.fem.functionspace(mesh, ("DG", 0))
tm_func = dolfinx.fem.Function(tm, name="medium")
tm_func.x.array[:] = ct.values
u_prev = u.x.array.copy()
tm_prev = tm_func.x.array.copy()
num_iterations = 0
with dolfinx.io.VTXWriter(mesh.comm, "solution.bp", [u, tm_func]) as bp:
    loading_steps = 33
    load = 1. / loading_steps
    dl = load
    n = 0
    MAX_FAILURE = 5
    NUM_SUCCESSIVE_SOLVES = 0
    last_load = load
    while load <= 1:
        print(f"Loading step: {load:.3e}", flush=True)
        applied_z.value = full_disp[2] * load
        applied_y.value = full_disp[1] * load
        problem.solve()
        reason = problem.solver.getConvergedReason()
        num_iterations += problem.solver.getIterationNumber()
        n += 1
        if reason < 0:
            load = last_load + dl/(n+1)
            u.x.array[:] = u_prev.copy()
            tm_func.x.array[:] = tm_prev.copy()
            NUM_SUCCESSIVE_SOLVES = 0
        else:
            n = 0
            last_load = load
            NUM_SUCCESSIVE_SOLVES += 1
            bp.write(load)
            load += NUM_SUCCESSIVE_SOLVES * dl
            u_prev[:] = u.x.array.copy()
            tm_prev[:] = tm_func.x.array.copy()
        if n > MAX_FAILURE:
            print("Too many failures, aborting.")
            break
print(f"Total number of iterations: {num_iterations}")

I can’t get it to run to the full loading (only 82%), but as the paper doesn’t state the tolerances used or what update scheme I’m not surprised as this is a very non-linear problem.

kumar · April 14, 2026, 5:25pm

Dear Mr Dokken,
Thank you for your suggestion regarding the contact method. It provides a very interesting perspective on handling contact between two bodies, and I find the approach very insightful.

At the moment, I am working with an elasto-plastic formulation based on von Mises material model (Elasto-plastic analysis of a 2D von Mises material — Computational Mechanics Numerical Tours with FEniCSx). I find it challenging to directly integrate the third medium approach in the setting, particularly in terms of consistency between the energy-based formulation and the history-dependent plasticity model.

Additionally, I am using the elasto-plastic setting in my application involving problem of contact between an electrode and a sheet in a spot welding simulation, where both contact and material nonlinearity play important roles.

However, I will continue exploring the third medium method and attempt to integrate it into this formulation. In addition, I would greatly appreciate any recommendations you may have, based on your experience, on alternative contact approaches in FEniCSx for such problems.

Thanks again for your time and considerations.

kumar · April 15, 2026, 9:07pm

Dear Mr. Dokken,

I tried running your third medium contact example, but I could not open the .bp file in ParaView (it crashes once i tried to open the folder), so I switched the part of the code output to XDMF:

with dolfinx.io.XDMFFile(mesh.comm, "full_output.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)

    loading_steps = 33
    load = 1. / loading_steps
    dl = load
    n = 0
    MAX_FAILURE = 5
    NUM_SUCCESSIVE_SOLVES = 0
    last_load = load

    max_load = 0.3
    while load <= max_load:
        print(f"Loading step: {load:.3e}", flush=True)
        applied_z.value = full_disp[2] * load
        applied_y.value = full_disp[1] * load
        problem.solve()
        reason = problem.solver.getConvergedReason()
        num_iterations += problem.solver.getIterationNumber()
        n += 1
        if reason < 0:
            load = last_load + dl/(n+1)
            u.x.array[:] = u_prev.copy()
            tm_func.x.array[:] = tm_prev.copy()
            NUM_SUCCESSIVE_SOLVES = 0
        else:
            n = 0
            last_load = load
            NUM_SUCCESSIVE_SOLVES += 1
            #---------------------------------------------------------
            ### WRITE SOLUTION (time = load)
            xdmf.write_function(u, load)
            xdmf.write_function(tm_func, load)
           #---------------------------------------------------------
            load += NUM_SUCCESSIVE_SOLVES * dl
            u_prev[:] = u.x.array.copy()
            tm_prev[:] = tm_func.x.array.copy()

        if n > MAX_FAILURE:
            print("Too many failures, aborting.")
            break
print(f"Total number of iterations: {num_iterations}")

This works, but I am having issues with visualization of only the domain without third medium. If I first apply Warp By Vector , I see deformation of the full domain including the third medium. However, if I first apply Threshold using the medium field to remove the third medium, I can no longer properly apply Warp By Vector afterwards.

Am I missing the correct ParaView workflow, or is there a recommended way (possibly via code changes) to visualize only the solid body deformation? (Like in your visualised output using an xdmf file). Thanks for your time.

dokken · April 15, 2026, 9:14pm

What version of Paraview are you running, and are you selecting the ADIOS2VTXReader?

You should use the “Extract block” filter first to extract only the f (so it doesn’t say partial)

kumar · April 16, 2026, 7:32am

Dear Mr. Dokken,

Thank you for your response. I found that the issue was due to the old ParaView version (v5.12) when opening the .bp file. With a newer version, the visualization works correctly.

I also have a follow-up question. My problem involves mechanical, thermal, and electrical interactions. While the third medium approach seems suitable for mechanical contact, it may not be appropriate for thermal and electrical fields, as the artificial medium could introduce non-physical conduction effects.

As I am still new to contact modeling in FEniCSx, I would greatly appreciate any guidance you might be willing to share. In particular, I was wondering if the implementation provided by Mr. Nate (asimov-contact/python/demos at main · Wells-Group/asimov-contact · GitHub) would be a reasonable starting point to explore, or if there are other approaches or resources you would recommend, as I have not been able to find much related discussion on this in the forum.

Thank you very much for your time.

dokken · April 16, 2026, 7:46am

Yes, third medium doesn’t seem right if one need thermal and electrical field coupling.
In the repo @nate mentions, Sarah Roggendorf and I implemented contact methods (including thermal coupling) in DOLFINx. However, at the moment I do not have funding or time to keep it up to date with the latest release.

Topic		Replies	Views
Add contact in Fenics variational formulation	2	458	July 6, 2020
Is there a contact model in FEniCS for two imported physical volumes?	9	3925	June 5, 2021
Simulate elastoplastic contact using Fenicsx General dolfinx	3	126	December 2, 2025
Impact of two deformable meshes with velocity	3	851	January 22, 2021
Meshed rigid indenter in contact with a hyperelastic body in FEniCSx variational formulation boundary-conditions	1	51	May 8, 2026

Contact between two deformable bodies

Related topics