Contact boundary condition

Hello Everyone,
I am trying to conduct three point bending test. i want to apply contact boundary condition like the figure show. how to perform it in FEniCS？

I am trying to apply a point boundary condition. But the result is not right. I think that the boundary conditions are different from the experimental ones. Here is my attempt

``````from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
#creat mesh
##---------L=4W------------------------
wt=10.01
lt=4.0*wt
wn=0.35
an=4.9
mesh = Mesh()
editor = MeshEditor()
editor.open(mesh,'triangle', 2, 2)
editor.init_vertices(16)
editor.init_cells(16)

editor.close()
mesh=refine(mesh)
##---------mesh end -------------------
#--------------------------------------
# change num_refine for mesh refinement
#--------------------------------------
def refine_elemnt(mesh):
makers = MeshFunction("bool", mesh, mesh.topology().dim())
makers.set_all(False)
for cell in cells(mesh):
if abs(cell.midpoint().x()-23.0)<3.0:
makers[cell.index()] = True
mesh=refine(mesh,makers)
return mesh
num_refine = 5
for i in range(num_refine):
mesh=refine(mesh)
for i in range(0):
mesh=refine_elemnt(mesh)
#Import mesh file
#mesh = Mesh('mesh.xml')
#mesh = RectangleMesh(Point(0.0, 0.0), Point(46.0, 10.01), 400, 80)
## Output: the infomation of mesh
print(mesh.hmax())
print(mesh.hmin())
print(mesh.num_cells())
print(mesh.num_vertices())
print(mesh)
file_mesh = File("meshwang.pvd")
file_mesh << mesh

# Define function space

W = VectorFunctionSpace(mesh, 'CG', 1)
u, v = TrialFunction(W), TestFunction(W)

# Introduce manually the material parameters
E0 =3.04e3             #N/mm/mm test:3.04e3
nu=0.36                 #  test:0.36
mu = E0/(2*(1.0+nu))     #N/mm/mm  test: E0/(2*(1.0+nu))      hirshikesh: 8.0e3
lmbda = E0*nu/((1.0+nu)*(1.0-2.0*nu))       #N/mm/mm test:E0*nu/((1.0+nu)*(1.0-2.0*nu))  hirshikesh:12.0e3

# Constituive functions
def epsilon(u):
def sigma(u):
return 2.0*mu*epsilon(u)+lmbda*tr(epsilon(u))*Identity(len(u))

# Boundary conditions
#bot = CompiledSubDomain("near(x[1], 0.0) && on_boundary")
def top(x, on_boundary):
tol = DOLFIN_EPS
return (abs(x[1] - wt) < tol) and on_boundary
def pointl(x, on_boundary):
tol = DOLFIN_EPS
return (abs(x[0] - 0.0) < tol) and (abs(x[1] -0.0) < tol)
def pointr(x, on_boundary):
tol = DOLFIN_EPS
return (abs(x[0] - 4.0*wt) < tol) and (abs(x[1] -0.0) < tol)
def pointm(x, on_boundary):
tol = DOLFIN_EPS
return (abs(x[0] - 2.0*wt) <= tol) and (abs(x[1] -10.01) <= tol)
load = Expression("t", t = 0.0, degree=1)
bc_ptl = DirichletBC(W, Constant((0.0,0.0)), pointl,method="pointwise")
bc_ptr = DirichletBC(W.sub(1), Constant(0.0), pointr,method="pointwise")
bc_u = [bc_ptl, bc_ptr,bc_ptm ]

boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
AutoSubDomain(top).mark(boundaries,1)
ds = Measure("ds")(subdomain_data=boundaries)
n = FacetNormal(mesh)
#plot(mesh)
# Apply a point load equal to 100 N downward in Y direction at coordinate (23.0,10.01):

def __init__(self, **kwargs):
super().__init__(degree=kwargs["degree"])

self.point = kwargs['point']
self.value = kwargs['value']

def eval(self, value, x):

if near(x[0], self.point[0]) and near(x[1], self.point[1]):
value[0] = self.value[0]
value[1] = self.value[1]
else:
value[0] = 0
value[1] = 0

def value_shape(self):
return (2,)
Tra = Load(point=(20.02, 10.02), value=(0.0, -100.0), degree=1)
# Variational form
f = Constant((0.0, 0.0))
L =  dot(f, v)*ds(1)

t = 0
forcex = []
forcey = []
x_dofs = W.sub(0).dofmap().dofs()
y_dofs = W.sub(1).dofmap().dofs()
# Staggered scheme
while t<0.35:    # real 0.007  test 0.38 mm
# Compute solution
u = Function(W)
solve(a == L, u, bc_u)

bc_ptm.apply(f_ext_unknown)
Fy = 0
for i in y_dofs:
Fy += f_ext_unknown[i]
forcey.append(-Fy)