Error in method = pointwise after definition of facets

General
1

I am trying to use method=pointwise in dirichletBC. However, after definition of facets, it is showing error: A SubDomain is required for pointwise search. Below is the MWC

from dolfin import *
mesh = BoxMesh(Point(-5,-10,0),Point(5,0,1),50,50,5)

class Corner(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] ==-5.0 and x[1]==-10.0

class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], -10.0) and on_boundary

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)
Corner().mark(facets, 1)
Bottom().mark(facets,2)

bc_u = [DirichletBC(FS_disp.sub(1), Constant(0.0), facets,2),
        DirichletBC(FS_disp.sub(1), Constant(0.0), facets,1, method = 'pointwise') ]

DispTrial, DispTest = TrialFunction(FS_disp), TestFunction(FS_disp)
disp = Function(FS_disp, name = "displacement")
E_disp = inner(grad(DispTest),sigma(DispTrial))*dx

problem_disp = LinearVariationalProblem(lhs(E_disp), rhs(E_disp), disp, bc_u)
solver_disp = LinearVariationalSolver(problem_disp)

solver_disp.solve()

File("Disp.pvd") << disp`
2

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
3

Can we not mark this point in facets as top_right_front_point().mark(facets,4) and apply in DirichletBC? I want to add this specific point on the facets.

4

Again, if you want to apply Dirichlet BC on a point, this could be done as I suggested above. This is how you can refer to a specific point when you want to apply the BC in a so called “pointwise” manner. However, I am not sure what you mean by adding point on the facet.

5

It doesn’t make sense to use the method pointwise with a facet.

The method pointwise uses the fact that:

  1. In legacy dolfin, a degree of freedom has a corresponding coordinate in physical space.
  2. Method=pointwise checks if the coordinate of a dof matches some geometrical requirement.

The standard method; topological can do two things:

  1. given a facet, if finds all dofs associated with that facet and applies the boundary condition to them,
  2. given a subdomain (some expression that takes in the mesh coordinates and returns a true/false value) it checks if the vertices of the facets satisfies this criterion (it can also check the midpoint of the facet). If all of these satisfies the condition, all relate dofs has boundary conditions applied to them.