3 point bending test using phase-field method

variational formulation
1

Hi, I’m currently trying to solve the 3 point bending problem shown in the figure below. I think that I have a problem with the boundary conditions. Could you please help me to fix it?
Capture d’écran_2022-06-28_10-58-48

In my mesh, I defined a small segment containing 2 elements as a physical curve for each point. These are my boundary conditions:

#Converting the mesh
mesh_file     = path_file + base_file
os.system("gmsh -2 "+mesh_file+".geo -format msh2")
os.system("dolfin-convert "+mesh_file+".msh "+mesh_file+".xml")
mesh          = Mesh(mesh_file +'.xml')

Regions_file  = path_file + base_file + '_physical_region.xml'
Regions_zone  = MeshFunction("size_t", mesh, Regions_file)

Boundary_file = path_file + base_file + '_facet_region.xml'
Boundary_zone = MeshFunction('size_t', mesh, Boundary_file)

# Defining FEM Space
V   = FunctionSpace(mesh, 'CG', 1)
W   = VectorFunctionSpace(mesh, 'CG', 1)
WW  = FunctionSpace(mesh, 'DG', 0)

# Trial and test function
u, v = TrialFunction(W), TestFunction(W)  #displacement
p, q = TrialFunction(V), TestFunction(V)  #phase-field

# id zones
Bot_id_1 = 100
Bot_id_2 = 200
Top_id = 300

bc_load1 = DirichletBC(W.sub(1),  load, Boundary_zone, Top_id)   # pull definition
bc_load2 = DirichletBC(W.sub(0), Constant(0.0), Boundary_zone, Top_id)   # pull definition
bc_fix1  = DirichletBC(W.sub(1), Constant(0.0), Boundary_zone, Bot_id_1)   # fix definition
bc_fix2  = DirichletBC(W.sub(1), Constant(0.0), Boundary_zone, Bot_id_2)   # fix definition

bc_u = [bc_fix1, bc_fix2, bc_load1, bc_load2]

ds = Measure("ds", domain=mesh, subdomain_data=Boundary_zone, subdomain_id=Top_id)
n = FacetNormal(mesh)