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?
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)