Three point bending test

1

Hello Everyone,
I am trying to conduct three point bending test for phase field modeling. However, I failed to create the boundary condition and point load properly. Would anyone show me properly how to implement BCs. Here is my attempt :

import matplotlib.pyplot as plt
from fenics import *
import mshr


# Load mesh and define functional spaces

Nn = 80
dom1 = mshr.Rectangle(Point(-1.0, -0.50), Point(1.0, 0.50))
dom2 = mshr.Rectangle(Point(0.0-1e-03, -0.5), Point(0.0+0.0+1e-03, 0.0))
dom = dom1 - dom2
mesh = mshr.generate_mesh(dom, Nn)
#plot(mesh)
#plt.show()


dim  = mesh.topology().dim()
V    = FunctionSpace(mesh, 'Lagrange', 1)
W    = VectorFunctionSpace(mesh, 'Lagrange', 1, dim)


class Left(SubDomain):
    def inside(self,x,on_boundary):
        return near(x[0]+1.0, DOLFIN_EPS) and near(x[1]+0.50, DOLFIN_EPS) 

class Right(SubDomain):
    def inside(self,x,on_boundary):
        return near(x[0]-1.0, DOLFIN_EPS) and near(x[1]+0.50, DOLFIN_EPS) 

class Point(SubDomain):
    def inside(self, x , on_boundary):
        return near(x[0], DOLFIN_EPS) and near(x[1]-0.50, DOLFIN_EPS)


def Crack(x):
    return abs(x[0]) < 1e-03 and x[1] <= 0.0
bc_phi = [DirichletBC(V, Constant(1.0), Crack)]


subd = MeshFunction("size_t", mesh, 0)

ptbl = Left()
ptbl.mark(subd, 1)
ptbr = Right()
ptbr.mark(subd, 2)
ptbp = Point()
ptbp.mark(subd, 3)

ud = Expression("-1*t", t=0.0, degree=1)
bcleft= DirichletBC(W, Constant((0.0,0.0)), ptbl, method='pointwise')
bcright= DirichletBC(W, Constant((0.0,0.0)), ptbr, method='pointwise')
bcPoint= DirichletBC(W.sub(1), ud, ptbp)
bc_u = [bcleft, bcright, bcPoint]