Dear all,
I am trying to implement a Hertz contact in axisymmetric coordinates, following the technics described in this tutorial, while translating it into dolfinx v0.10.
Here is my code:
import numpy as np
import ufl
from dolfinx import mesh, fem
from dolfinx.fem.petsc import NonlinearProblem
from mpi4py import MPI
from ufl import tr, sym, grad, Identity, inner
from dolfinx import plot
import pyvista
# --- 1. MESH & MAPPING ---
N_r, N_z = 30, 30
R_domain, H = 1.0, 1.0
# Create simple rectangle
domain = mesh.create_rectangle(MPI.COMM_WORLD, [[0.0, -H], [R_domain, 0.0]], [N_r, N_z], cell_type=mesh.CellType.triangle)
def map_coordinates(x):
y = np.zeros_like(x)
# Jacobian is never zero
y[:, 0] =R_domain*( 0.2 * x[:, 0]/R_domain + 0.8 * x[:, 0]**2/R_domain**2)
y[:, 1] =-H*( -0.2 * x[:, 1]/H + 0.8 * x[:, 1]**2/H**2)
return y
domain.geometry.x[:] = map_coordinates(domain.geometry.x)
# --- 2. FACET TAGGING ---
f_dim = domain.topology.dim - 1
def top(x): return np.isclose(x[1], 0.0)
def bottom(x): return np.isclose(x[1], -H)
def axis(x): return np.isclose(x[0], 0.0)
# Locate and tag
top_facets = mesh.locate_entities_boundary(domain, f_dim, top)
bottom_facets = mesh.locate_entities_boundary(domain, f_dim, bottom)
axis_facets = mesh.locate_entities_boundary(domain, f_dim, axis)
# Combine into meshtags (avoiding duplicates)
indices = np.hstack([top_facets, bottom_facets, axis_facets])
values = np.hstack([np.full(len(top_facets), 1, np.int32),np.full(len(bottom_facets), 2, np.int32),np.full(len(axis_facets), 3, np.int32)])
# Sort indices for meshtags
sorted_idx = np.argsort(indices)
facet_tags = mesh.meshtags(domain, f_dim, indices[sorted_idx], values[sorted_idx])
print(np.unique(facet_tags.values))
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tags)
# --- 3. SPACES & BCs ---
V = fem.functionspace(domain, ("CG", 1, (2,)))
u = fem.Function(V)
du = ufl.TrialFunction(V)
u_ = ufl.TestFunction(V)
# Ur = 0 on axis (ID 3)
bc_axis = fem.dirichletbc(fem.Constant(domain, 0.0), axis_facets, V.sub(0))
# Uz = 0 on bottom (ID 2)
bc_bottom = fem.dirichletbc(fem.Constant(domain, 0.0), bottom_facets, V.sub(1))
bcs = [bc_axis, bc_bottom]
# --- 4. AXISYMMETRIC PHYSICS ---
E, nu = 10.0, 0.3
mu = fem.Constant(domain, E / (2*(1+nu)))
lmbda = fem.Constant(domain, E*nu / ((1+nu)*(1-2*nu)))
# --- Mesh Size Calculation ---
# 1. Get the topological dimension
tdim = domain.topology.dim
# 2. Get the indices of all local cells
num_cells = domain.topology.index_map(tdim).size_local
cells = np.arange(num_cells, dtype=np.int32)
# 3. Call the .h() method on the OBJECT (domain)
h_values = domain.h(tdim, cells)
# 4. Now compute the average
avg_h = np.mean(h_values)
# Scaled Penalty: E / avg_h ensures numerical stability
pen_float = 1e5 * float(E) / avg_h
pen = fem.Constant(domain, pen_float)
print('penalty=',pen_float)
R_ind = 0.5
x_c = ufl.SpatialCoordinate(domain)
# Numerical safety for r=0
r = x_c[0]
r_safe = ufl.conditional(r > 1e-10, r, 1e-10)
d_ind = fem.Constant(domain, 0.05)
obstacle = -d_ind + (r_safe**2)/(2*R_ind)
def eps_axi(v): return sym(grad(v))
def eps_theta(v): return v[0] / r_safe
def div_axi(v): return tr(eps_axi(v)) + eps_theta(v)
def sigma_axi(v): return lmbda * div_axi(v) * Identity(2) + 2 * mu * eps_axi(v)
def sigma_theta(v): return lmbda * div_axi(v) + 2 * mu * eps_theta(v)
# --- 5. VARIATIONAL FORM ---
#Contact
penetration = ufl.max_value(0, -obstacle + u[1])
#Energy
Pi_int = 0.5*(inner(sigma_axi(u), eps_axi(u)) + sigma_theta(u) * eps_theta(u)) * 2.0 * np.pi*r_safe * ufl.dx
Pi_con = 0.5 * pen * penetration**2 * 2*np.pi*r_safe * ds(1)
Pi = Pi_int + Pi_con
F = ufl.derivative(Pi, u, u_)
J = ufl.derivative(F, u, du)
# --- 1. Update the Constant ---
d_target = 0.05
n_steps = 10
d_increments = np.linspace(0, d_target, n_steps + 1)[1:]
#d_increments = [0]
# --- 2. Robust Solver Options ---
petsc_opts = {
"ksp_type": "preonly",
"pc_type": "lu",
"pc_factor_mat_solver_type": "mumps",
"snes_type": "newtonls",
"snes_linesearch_type": "bt",
"snes_max_it": 50,
"snes_atol": 1e-8,
"snes_rtol": 1e-8,
}
problem = NonlinearProblem(F, u, bcs=bcs, J=J, petsc_options_prefix="hertzcontact_", petsc_options=petsc_opts)
print(f"Starting indentation loop (Total target: {d_target})...")
for i, d_val in enumerate(d_increments):
# Update the indentation depth constant
d_ind.value = d_val
print(f"\nStep {i+1}/{n_steps}: Indentation = {d_val:.4f}")
try:
n_iters, converged = problem.solve()
if converged:
print(f" Converged in {n_iters} iterations.")
else:
print(f" Warning: Step {i+1} did not converge!")
break # Stop if we lose convergence
except Exception as e:
print(f" Critical Failure: {e}")
break
There should be an error in the code but I can’t find it.
Besides, I don’t understand the following errors I get:
corrupted double-linked list (not small)
Aborted (core dumped)
And
[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range
[0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
[0]PETSC ERROR: or see FAQ — PETSc 3.24.3 documentation and FAQ — PETSc 3.24.3 documentation
[0]PETSC ERROR: configure using --with-debugging=yes, recompile, link, and run
[0]PETSC ERROR: to get more information on the crash.
[0]PETSC ERROR: Run with -malloc_debug to check if memory corruption is causing the crash.
Abort(59) on node 0 (rank 0 in comm 0): application called MPI_Abort(MPI_COMM_WORLD, 59) - process 0
It seems also that the Spyder kernel “is dying” during each simulation.
The previous errors are gotten from running the code in the terminal.
I am using dolfinx 0.10.
Thank you in advance