How to define Dirichlet Boundary Conditions with MeshFunction?

Question

Using the demo_stokes-taylor-hood.py demo I wrote a script to
run a 2D cavity case. First I defined the Dirichlet conditions as follow:

from dolfin import *

mesh = UnitSquare(50,50)

def BCtop(x, on_boundary):
    return x[0] > DOLFIN_EPS  and x[0] < 1. - DOLFIN_EPS and x[1] > 1. - DOLFIN_EPS  and  on_boundary
def BCslip(x, on_boundary):
    return      (x[0] < DOLFIN_EPS and on_boundary)\
            or  (x[0] > 1. - DOLFIN_EPS and on_boundary)\
            or  (x[0] > DOLFIN_EPS  and x[0] < 1. - DOLFIN_EPS and x[1] < DOLFIN_EPS  and  on_boundary)

# Define function spaces
V = VectorFunctionSpace(mesh, "CG", 2)
Q = FunctionSpace(mesh, "CG", 1)
W = V * Q

# Noslip boundary condition for velocity
noslip = Constant((0, 0))
bc0 = DirichletBC(W.sub(0), noslip, BCslip)

# Top flow boundary condition for velocity
inflow = Expression(("sin(x[0]*pi)", "0.0"))
bc1 = DirichletBC(W.sub(0), inflow, BCtop)


# Collect boundary conditions
bcs = [bc0, bc1]

# Define variational problem
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
f = Constant((0, 0))
a = (inner(grad(u), grad(v)) - div(v)*p + q*div(u))*dx
L = inner(f, v)*dx


# Compute solution
w = Function(W)

A = assemble(a)
b = assemble(L)
for bc in bcs: bc.apply(A,b)
solve(A,w.vector(),b)


(u, p) = w.split()
plot(u[0])
plot(u[1])
plot(u)
plot(p)
interactive()

This is working fine. Then I tried using a MeshFunction to define the boundaries but it did not work. Is it not allowed to proceed this way ?

from dolfin import *

mesh = UnitSquare(50,50)

boundary_parts = MeshFunction("uint", mesh, mesh.topology().dim()-1)

class top_boundary(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] > DOLFIN_EPS  and x[0] < 1. - DOLFIN_EPS and x[1] > 1. - DOLFIN_EPS  and  on_boundary
class noslip_boundary(SubDomain):
    def inside(self, x, on_boundary):
        return      (x[0] < DOLFIN_EPS and on_boundary)\
                or  (x[0] > 1. - DOLFIN_EPS and on_boundary)\
                or  (x[0] > DOLFIN_EPS  and x[0] < 1. - DOLFIN_EPS and x[1] < DOLFIN_EPS  and  on_boundary)

noslip_bnd = noslip_boundary()
noslip_bnd.mark(boundary_parts, 0)
top_bnd = top_boundary()
top_bnd.mark(boundary_parts, 1)

# Define function spaces
V = VectorFunctionSpace(mesh, "CG", 2)
Q = FunctionSpace(mesh, "CG", 1)
W = V * Q

# Noslip boundary condition for velocity
noslip = Constant((0, 0))
bc0 = DirichletBC(W.sub(0), noslip, boundary_parts, 0)

# Inflow boundary condition for velocity
inflow = Expression(("sin(x[0]*pi)", "0.0"))
bc1 = DirichletBC(W.sub(0), inflow, boundary_parts, 1)



# Collect boundary conditions
bcs = [bc0, bc1]

# Define variational problem
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
f = Constant((0, 0))
a = (inner(grad(u), grad(v)) - div(v)*p + q*div(u))*dx
L = inner(f, v)*dx


# Compute solution
w = Function(W)

A = assemble(a)
b = assemble(L)
for bc in bcs: bc.apply(A,b)
solve(A,w.vector(),b)


(u, p) = w.split()
plot(u[0])
plot(u[1])
plot(u)
plot(p)
interactive()

Comments

To maximise the chances of an answer, you should post simple but complete code.

---

You need to tell us what is the actual problem because your code passes ok with mesh = UnitSquare(4, 4).

Answer 1

You need to initialize boundary_parts to say 2

boundary_parts = MeshFunction("uint", mesh, mesh.topology().dim()-1)
boundary_parts.set_all(2)

# now we can set other boundary_parts values by SubDomain.mark() method

because

• in DOLFIN 1.0.0 uninitialized MeshFunction takes random values so DirichletBC can randomly apply also to some interior facets

• in newer DOLFIN MeshFunction is automatically initialized to 0 so bc0 is applied to all interior facets as a side effect

Comments

Thank you very much for your answer!!

Answer 2

I made it work after substituting the lines

bc0 = DirichletBC(W.sub(0), noslip, boundary_parts, 0)
...
bc1 = DirichletBC(W.sub(0), inflow, boundary_parts, 1)

by

bc0 = DirichletBC(W.sub(0), noslip, noslip_bnd)
...
bc1 = DirichletBC(W.sub(0), inflow, top_bnd)