# Penalty Contact with circular shaped boundary Condition

Hello everyone,

i am currently working on a contact model for a multibody contact of a rigid circular shaped indenter and an elastic domain. The indenter should be realized as a boundary condition because of a considerable difference in stiffness (steel compared to rubber).

This topic refers to:

However, I have difficulties creating a penalty term which use both x and y displacements of the rigid indenter as an additional energy term in the formulation. I tried lots of different ways to implement a circle equation but nothing worked. My idea was to use something like this:

``````def in_circle(x):
return np.less_equal( ((x[0]-center_x)**2 + (x[1]-center_y)**2), ((D/2)**2))
``````

The MWE should realize the sketch below, where a circular boundary should cause a static deformation realized with linear elasticity and the penalty approach:

The code for the dolfinx MWE + mesh generation using gmsh:

``````import gmsh
import numpy as np
from dolfinx import fem
from dolfinx.io import XDMFFile, gmshio
from mpi4py import MPI
from petsc4py.PETSc import ScalarType, Options
from ufl import dx, inner, derivative, tr, Identity, sym, grad, TrialFunction, TestFunction, Measure

gmsh.initialize()
lc=1e-1

# generate points for inner contour
# generate inner contour
# generate points for outer contour
# generate outer contour
# generate curveloops
# generate surface
gmsh.model.geo.synchronize()
# generate physical groups for boundary conditions
inside, outside, surface = 1, 2, 3
gmsh.model.addPhysicalGroup(1, [1, 2, 3, 4], inside, name="Inner Contour")
gmsh.model.addPhysicalGroup(1, [5, 6, 7, 8], outside, name="Outer Contour")
# finish gmsh
gmsh.model.mesh.generate(2)
gmsh.write("rectangle.msh")
# create a DOLFINx mesh
dim = 2
gmsh_model_rank = 0
domain_com = MPI.COMM_SELF
domain, cell_markers, facet_markers = gmshio.model_to_mesh(gmsh.model, domain_com, gmsh_model_rank, dim)
domain.name = "Stator"
cell_markers.name = f"{domain.name}_cells"
facet_markers.name = f"{domain.name}_facets"

with XDMFFile(domain.comm, f"rectangle_mesh_{MPI.COMM_SELF.rank}.xdmf", "w") as file:
file.write_mesh(domain)
file.write_meshtags(cell_markers)
domain.topology.create_connectivity(domain.topology.dim - 1, domain.topology.dim)
file.write_meshtags(facet_markers)

gmsh.finalize()

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

# define inner boundry

# define function space
V = fem.VectorFunctionSpace(domain, ("CG",1))

# definition of circular shape
########## Insert 'circle' here ##########

# test and trial functions
u = fem.Function(V, name="Displacement")
du = TrialFunction(V)
u_ = TestFunction(V)

# define outer contour as homogenous boundary condition
u_zero = np.array((0,)*domain.geometry.dim, dtype=ScalarType)
bc = fem.dirichletbc(u_zero, fem.locate_dofs_topological(V, fdim, facet_markers.find(outside)), V)

# constituive relation
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*tr(eps(v))*Identity(2) + 2.0*mu*eps(v)
def ppos(x):
return (x+np.abs(x))/2.

pen = fem.Constant(domain, ScalarType(1e4))
form = inner(sigma(u), eps(u_))*dx
#form = inner(sigma(u), eps(u_))*dx +pen*dot(u_, ppos(u-circle))*ds(1)
J = derivative(form, u, du)

problem = fem.petsc.NonlinearProblem(form, u, bcs=[bc], J=J)
from dolfinx import nls
solver = nls.petsc.NewtonSolver(domain.comm, problem)

with XDMFFile(domain.comm, "rec_output.xdmf", "w") as file:
file.write_mesh(domain)
file.write_function(u)

``````

I use FEniCSx 0.6.0 inside a Docker container.