# Implementation of the Penalty approach in dolfinx - Hertzian contact

Hello everyone,

I am currently trying to transfer the " Hertzian contact with a rigid indenter using a penalty approach" example (Comet-FEniCS Example) into FEniCSx but I have difficulties formulating the penalty-formulation.

So far the code looks like this and fails in the last line:

``````import numpy as np
from dolfinx import fem, mesh
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
import ufl

N=30
domain = mesh.create_unit_cube(MPI.COMM_WORLD, N, N, N//2)

def bottom(x):
return np.isclose(x,0)

def top(x):
return np.isclose(x, N//2)

fdim = domain.topology.dim -1
bottom_facets = mesh.locate_entities_boundary(domain,fdim,bottom)
top_facets = mesh.locate_entities_boundary(domain,fdim,top)

marked_facets = np.hstack([bottom_facets, top_facets])
marked_values = np.hstack([np.full_like(bottom_facets, 1), np.full_like(top_facets, 2)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])

V = fem.VectorFunctionSpace(domain, ("CG",1))
u = fem.Function(V, name="Displacement")
du = ufl.TrialFunction(V)
u_ = ufl.TestFunction(V)

u_zero = np.array((0,)*domain.geometry.dim, dtype=ScalarType)
tdim = domain.topology.dim
fdim = tdim - 1
bc = fem.dirichletbc(u_zero, fem.locate_dofs_topological(V, fdim, facet_tag.find(1)), V)

R = 0.5
d = 0.02

x = ufl.SpatialCoordinate(domain)
circle = fem.Function(V)
circle.interpolate(lambda x: np.stack((x, x, -d+(pow(x,2)+pow(x, 2))/2/R)))

E = fem.Constant(domain, ScalarType(10.))
nu = fem.Constant(domain, ScalarType(0.3))
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
def eps(v):
def sigma(v):
return lmbda*ufl.tr(eps(v))*ufl.Identity(3) + 2.0*mu*eps(v)
def ppos(x):
return (x+np.abs(x))/2.
pen = fem.Constant(domain, ScalarType(1e4))

form = ufl.inner(sigma(u), eps(u_))*ufl.dx + pen*ufl.dot(u_, ppos(u-circle))*ds
``````

The error I get is the following:

``````ufl.log.UFLException: Can't add expressions with different shapes.
``````

Since it is no longer possible to use the Expression class, I interpolated an anonymous in 3D. I am sure there is a proper way of implementing this expression but I have not found something comparable.

This Post comes quite close to my problem, but there was no answer to:

Are there particular ways that penalty methods have to be set up in FEniCSx if the level of penetration is a function of absolute position as opposed to the displacement (as is done in the Hertzian contact example)?

I use FEniCSx 0.6.0 inside a Docker container.

You could simply use the `ufl.SpatialCoordinate` to define your `circle`, i.e.

``````import numpy as np
from dolfinx import fem, mesh
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
import ufl

N = 30
domain = mesh.create_unit_cube(MPI.COMM_WORLD, N, N, N//2)

def bottom(x):
return np.isclose(x, 0)

def top(x):
return np.isclose(x, N//2)

fdim = domain.topology.dim - 1
bottom_facets = mesh.locate_entities_boundary(domain, fdim, bottom)
top_facets = mesh.locate_entities_boundary(domain, fdim, top)

marked_facets = np.hstack([bottom_facets, top_facets])
marked_values = np.hstack(
[np.full_like(bottom_facets, 1), np.full_like(top_facets, 2)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(
domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])

ds = ufl.Measure('ds', domain=domain,

V = fem.VectorFunctionSpace(domain, ("CG", 1))
u = fem.Function(V, name="Displacement")
du = ufl.TrialFunction(V)
u_ = ufl.TestFunction(V)

u_zero = np.array((0,)*domain.geometry.dim, dtype=ScalarType)
tdim = domain.topology.dim
fdim = tdim - 1
bc = fem.dirichletbc(u_zero, fem.locate_dofs_topological(
V, fdim, facet_tag.find(1)), V)

R = 0.5
d = 0.02

x = ufl.SpatialCoordinate(domain)
circle = -d+(pow(x, 2)+pow(x, 2))/2/R

E = fem.Constant(domain, ScalarType(10.))
nu = fem.Constant(domain, ScalarType(0.3))
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)

def eps(v):

def sigma(v):
return lmbda*ufl.tr(eps(v))*ufl.Identity(3) + 2.0*mu*eps(v)

def ppos(x):
return (x+abs(x))/2.

pen = fem.Constant(domain, ScalarType(1e4))

form = ufl.inner(sigma(u), eps(u_))*ufl.dx + pen * \
ufl.dot(u_, ppos(u-circle))*ds
``````

Thanks Dokken!

After setting up the nonlinear problem i now have a TypeError…

This code:

``````form = ufl.inner(sigma(u), eps(u_))*ufl.dx + pen * \
ufl.dot(u_, ppos(u-circle))*ds
J = ufl.derivative(form, u, du)

problem = fem.petsc.NonlinearProblem(form, u, bcs=[bc], J=J)
``````

``````Traceback (most recent call last):
File "/home/fenics/shared/mwe.py", line 62, in <module>
problem = fem.petsc.NonlinearProblem(form, u, bcs=[bc], J=J)
File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/petsc.py", line 683, in __init__
self._L = _create_form(F, form_compiler_options=form_compiler_options,
File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/forms.py", line 176, in form
return _create_form(form)
File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/forms.py", line 171, in _create_form
return _form(form)
File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/forms.py", line 165, in _form
return formcls(ufcx_form, V, coeffs, constants, subdomains, mesh, module.ffi, code)
File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/forms.py", line 51, in __init__
TypeError: __init__(): incompatible constructor arguments.
``````

this is because your definition of `ds` is wrong.
See following working code:

``````import numpy as np
from dolfinx import fem, mesh
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
import ufl

N = 30
domain = mesh.create_unit_cube(MPI.COMM_WORLD, N, N, N//2)

def bottom(x):
return np.isclose(x, 0)

def top(x):
return np.isclose(x, N//2)

fdim = domain.topology.dim - 1
bottom_facets = mesh.locate_entities_boundary(domain, fdim, bottom)
top_facets = mesh.locate_entities_boundary(domain, fdim, top)

marked_facets = np.hstack([bottom_facets, top_facets])
marked_values = np.hstack(
[np.full_like(bottom_facets, 1), np.full_like(top_facets, 2)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(
domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])

ds = ufl.Measure('ds', domain=domain,

V = fem.VectorFunctionSpace(domain, ("CG", 1))
u = fem.Function(V, name="Displacement")
du = ufl.TrialFunction(V)
u_ = ufl.TestFunction(V)

u_zero = np.array((0,)*domain.geometry.dim, dtype=ScalarType)
tdim = domain.topology.dim
fdim = tdim - 1
bc = fem.dirichletbc(u_zero, fem.locate_dofs_topological(
V, fdim, facet_tag.find(1)), V)

R = 0.5
d = 0.02

x = ufl.SpatialCoordinate(domain)
circle = -d+(x**2+x**2)/2/R

E = fem.Constant(domain, ScalarType(10.))
nu = fem.Constant(domain, ScalarType(0.3))
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)

def eps(v):

def sigma(v):
return lmbda*ufl.tr(eps(v))*ufl.Identity(3) + 2.0*mu*eps(v)

def ppos(x):
return (x+abs(x))/2.

pen = fem.Constant(domain, ScalarType(1e4))

form = ufl.inner(sigma(u), eps(u_))*ufl.dx + pen * \
ufl.dot(u_, ppos(u-circle))*ds
J = ufl.derivative(form, u, du)

problem = fem.petsc.NonlinearProblem(form, u, bcs=[bc], J=J)

``````

Could you please share the full working code?

In case someone wants to use the code.
It works now.

``````import numpy as np
import ufl
import math
from dolfinx import fem, mesh, log
from dolfinx.io import XDMFFile
from mpi4py import MPI
from petsc4py.PETSc import ScalarType, Options
from dolfinx import nls

N=30
domain = mesh.create_box(MPI.COMM_WORLD,[[0.0,0.0,0.0], [1.0, 1.0, -1.0]], [N, N, N//2], mesh.CellType.hexahedron)
x = ufl.SpatialCoordinate(domain)

def top(x):
return np.isclose(x, 0)

def symmetry_x(x):
return np.isclose(x, 0)

def symmetry_y(x):
return np.isclose(x, 0)

def bottom(x):
return np.isclose(x, -1.0)

def wall_x(x):
return np.isclose(x, 1.0)

def wall_y(x):
return np.isclose(x, 1.0)

tdim = domain.topology.dim
fdim = tdim - 1

bottom_facets = mesh.locate_entities_boundary(domain,fdim,bottom)
top_facets = mesh.locate_entities_boundary(domain,fdim,top)
symx_facets = mesh.locate_entities_boundary(domain, fdim, symmetry_x)
symy_facets = mesh.locate_entities_boundary(domain, fdim, symmetry_y)
wallx_facets = mesh.locate_entities_boundary(domain, fdim, wall_x)
wally_facets = mesh.locate_entities_boundary(domain, fdim, wall_y)

marked_facets = np.hstack([bottom_facets, top_facets, symx_facets, symy_facets, wallx_facets, wally_facets])
marked_values = np.hstack([np.full_like(bottom_facets, 1), np.full_like(top_facets, 2), np.full_like(symx_facets, 3), \
np.full_like(symy_facets, 4), np.full_like(wallx_facets, 5), np.full_like(wally_facets, 6)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])

V = fem.VectorFunctionSpace(domain, ("CG",1))
V2 = fem.FunctionSpace(domain, ("CG", 1))
V0 = fem.FunctionSpace(domain, ("DG", 0))

u = fem.Function(V, name="Displacement")
du = ufl.TrialFunction(V)
u_ = ufl.TestFunction(V)
gap = fem.Function(V2, name="Gap")
p = fem.Function(V0, name="Contact pressure")

u_zero = np.array((0,)*domain.geometry.dim, dtype=ScalarType)
symx_dof = fem.locate_dofs_topological(V.sub(0), fdim, symx_facets)
symy_dof = fem.locate_dofs_topological(V.sub(1), fdim, symy_facets)

bc1 = fem.dirichletbc(u_zero, fem.locate_dofs_topological(V, fdim, facet_tag.find(1)), V)
bc2 = fem.dirichletbc(ScalarType(0), symx_dof, V.sub(0))
bc3 = fem.dirichletbc(ScalarType(0), symy_dof, V.sub(1))
bc = [bc1, bc2, bc3]

R = 0.5
d = 0.02

circle = -d+(pow(x, 2)+pow(x, 2))/2/R

E = fem.Constant(domain, ScalarType(10.))
nu = fem.Constant(domain, ScalarType(0.3))
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
def eps(v):
def sigma(v):
return lmbda*ufl.tr(eps(v))*ufl.Identity(3) + 2.0*mu*eps(v)
def ppos(x):
return (x+np.abs(x))/2.
pen = fem.Constant(domain, ScalarType(1e5))

form = ufl.inner(sigma(u), eps(u_))*ufl.dx + pen * \
ufl.dot(u_, ppos(u-circle))*ds(2)
J = ufl.derivative(form, u, du)

problem = fem.petsc.NonlinearProblem(form, u, bc, J=J)
solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD, problem)
ksp = solver.krylov_solver
opts = Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "cg"
opts[f"{option_prefix}pc_type"] = "gamg"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()

log.set_log_level(log.LogLevel.INFO)
n, converged = solver.solve(u)
assert(converged)
print(f"Number of interations: {n:d}")

sig = sigma(u)[2, 2]
stress_expr = fem.Expression(sig, V0.element.interpolation_points())
stresses = fem.Function(V0)
stresses.interpolate(stress_expr)
maxstress = max(np.abs(stresses.vector.array))

gap = circle-u
gap_expr = fem.Expression(gap, V2.element.interpolation_points())
gapval = fem.Function(V2)
gapval.interpolate(gap_expr)
maxgap = max(np.abs(gapval.vector.array))

a = math.sqrt(R*d)
F = 4/3.*float(E)/(1-float(nu)**2)*a*d
p0 = 3*F/(2*math.pi*a**2)

with XDMFFile(MPI.COMM_WORLD, "cube_dolfinx.xdmf", "w") as xdmf:
xdmf.write_mesh(domain)
xdmf.write_meshtags(facet_tag)
xdmf.write_function(u, 0.)

print(f'Contactarea:           {a:.3f} mm²')
print(f'Force:                 {F:.3f} N')
print(f'max Pressure:          {p0:.3f} MPa')
print(f'Computed max Pressure: {maxstress:.3} MPa')

``````
2 Likes

is here a `subdomain_id=2` missing?

EDIT: I see it in your form: `form = ufl.inner(sigma(u), eps(u_))*ufl.dx + pen * \ ufl.dot(u_, ppos(u-circle))*ds(2)`

I assume it does the same work as

``````subdomain_id=2
``````

Exactly. Both expressions should work in theory…
I tried both but the verision in the MWE allows to call multiple ds() in the formulation without defining a new measure each time.
At least as I understand it.

1 Like

Ah, yes it is so. Thank you!