It seems like you want to apply the Dirichlet BC on a specific point which means that you need to define a point first. Lets say you want to fix a point on a top-right corner in front of the block. You can define the point as a function like:
#FIX the point on the top right corner in front face
def top_right_front_point(x, on_boundary):
tol = DOLFIN_EPS
return (abs(x[0] - 5) < tol) and (abs(x[1] - 0) < tol) and (abs(x[2] - 1) < tol)
and then impose Dirichlet BC (e.g. Zero Displacement) on this point as:
DirichletBC(FS_disp, Constant((0., 0.,0.)), top_right_front_point, method="pointwise")
In addition, you have not applied any load on the cube to see a deformation when solving the linear elasticity problem that should be defined in the linear part of the variational problem. With that being said, I changed your original code a little to fix the bottom face , apply a pressure load on the left face and zero displacement at the point on the top-right corner (On the front face of the cube) as following:
from dolfin import *
mesh = BoxMesh(Point(-5,-10,0),Point(5,0,1),5,5,1)
class Bottom(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], -10.0) and on_boundary
class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], -5.0) and on_boundary
#FIX the point on the top right corner in front face
def top_right_front_point(x, on_boundary):
tol = DOLFIN_EPS
return (abs(x[0] - 5) < tol) and (abs(x[1] - 0) < tol) and (abs(x[2] - 1) < tol)
E = Constant(210e3)
nu = Constant(0.3)
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
def eps(v):
return sym(grad(v))
def sigma(v):
return 2.0*mu*eps(v) + lmbda*tr(eps(v))*Identity(3)
FS_disp = VectorFunctionSpace(mesh, "CG", 1)
facets = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
facets.set_all(0)
Bottom().mark(facets,2)
Left().mark(facets,3)
ds = Measure("ds", subdomain_data=facets)
n = FacetNormal(mesh)
bc_u = [DirichletBC(FS_disp, Constant((0.0,0.0,0.0)), facets,2),DirichletBC(FS_disp, Constant((0., 0.,0.)), top_right_front_point, method="pointwise") ]
DispTrial, DispTest = TrialFunction(FS_disp), TestFunction(FS_disp)
disp = Function(FS_disp, name = "displacement")
E_disp = inner(grad(DispTest),sigma(DispTrial))*dx
l = inner(-Constant(1e4) * n, DispTest) * ds(3) #Apply the pressure load on the left face
solve(E_disp == l, disp, bc_u)
File("Disp.pvd") << disp