DirichletBC for higher order polynomial space

Question

In the following code, I define two Dirichlet bc on a unit square domain, one on the left boundary ("x[0]<=DOLFIN_EPS") and one on the "top and bottom" (defined by "on_boundary and x[0]>=DOLFIN_EPS and x[0]<=1-DOLFIN_EPS"). It seems to me that FEniCS doesn't interpret the top-bottom boundary condition correctly. If A=(0,0) is the point in the lower left corner and B=(0.1, 0) is its closest point on the "bottom boundary", then A is contained in the "left" and B in the "top-bottom" region. However, using a quadratic function space, FEniCS doesn't apply the bc to the midpoint of the AB facet -- which is used in the computation -- although this point is contained in the "top-bottom" region. It can be vizualized by projecting the solution to a finer mesh where the (A+B)/2 point is visible and the solution is nonzero there. In this case, it really doesn't matter but I encountered this in a problem where it changed the solution completely..

from dolfin import *
mesh=UnitSquareMesh(10,10)
mesh_ref=refine(mesh)

V = FunctionSpace(mesh, "CG", 2)
V_ref = FunctionSpace(mesh_ref, "CG", 1)

def left(x, on_boundary):
   return on_boundary and x[0]<=DOLFIN_EPS
def bt(x, on_boundary):
   return on_boundary and x[0]>=DOLFIN_EPS and x[0]<=1.0-DOLFIN_EPS

left_bc = DirichletBC(V, Expression('x[1]*(1.0-x[1])') , left)
bt_bc = DirichletBC(V, Constant(0.0), bt)
bcs = [left_bc, bt_bc]

u,v=TrialFunction(V), TestFunction(V)
F = inner(grad(u), grad(v))*dx
a,L=system(F)

sol=Function(V)
solve(a == L, sol, bcs)
plot(sol, title="solution")
# interpolate to finer mesh:
plot(interpolate(sol, V_ref), title="interpolation on refined mesh", interactive=True)

Answer 1

Hi

I don't think there's anything wrong here, except poor interpolation near the corner. All values on the refined mesh (also those on the boundary) are interpolated from the coarse one, which is why you see the kink. Create and apply boundary conditions for the refined mesh as well and you should be fine

u_ref = interpolate(sol, V_ref)
left_bc2 = DirichletBC(V_ref, Expression('x[1]*(1.0-x[1])') , left)
bt_bc2 = DirichletBC(V_ref, Constant(0.0), bt)
bcs2 = [left_bc2, bt_bc2]
[bc.apply(u_ref.vector()) for bc in bcs2]
plot(u_ref, title="interpolation on refined mesh", interactive=True)

In reference to your comment I see now that your bt inside method does not cover all facets on the bottom and top boundaries (misses corner dofs). Try

def up_down(x, on_boundary):
    return on_boundary and near(x[1]*(1-x[1]), 0)

Comments

Mikael, no, unfortunately, it is not about interpolation. I'm sending a longer but probably more ilustrative example where it is clear that the solution is completely wrong -- Navier-Stokes equation on a square with an inflow on the left boundary. Try to comment the second line of the up_down function and uncomment the third line and see the difference. Now the solution is wrong, when you comment/uncomment, then it will be right. The point is that in the first case, the flow runs away from the room by a tiny hole around the left lower and upper corner -- the hole is, in fact, in the middle between two points in the original mesh..

from dolfin import *
mesh =UnitSquareMesh(20,20)
mesh_ref=refine(mesh)
max_time = 10 # maximal time
t=0.0
nu = 0.001 # kinematic viscosity
dt=0.01 # time step

# Function spaces 
V = VectorFunctionSpace(mesh, "CG", 2)
Q = FunctionSpace(mesh, "CG", 1)
mixed=V*Q
V_ref=VectorFunctionSpace(mesh_ref, "CG", 1)

def left(x, on_boundary):
    return on_boundary and x[0]<=DOLFIN_EPS
def up_down(x, on_boundary):
    return on_boundary and x[0]>=DOLFIN_EPS and x[0]<=1.0-DOLFIN_EPS
    #return on_boundary and (x[1] <=DOLFIN_EPS or x[1] >= 1-DOLFIN_EPS)

v_in = Expression(('x[1]*(1.0-x[1])', '0'))  # inflow
inflow = DirichletBC(mixed.sub(0), v_in, left)
noslip = DirichletBC(mixed.sub(0), ('0', '0'), up_down)
bcs = [inflow, noslip]

u_old = Function(V) # old velocity
press = Function(Q) # pressure
sol = Function(mixed) # (u,p)
u_ref=Function(V_ref)

A,B = TrialFunction(mixed), TestFunction(mixed)
u,p = split(A)
v,q = split(B)

# PDE
F = (1/dt)*inner(u - u_old, v)*dx + inner(grad(u)*u_old, v)*dx + nu*inner(grad(u), grad(v))*dx - p*div(v)*dx - q*div(u)*dx
a, L = system(F)

while t < max_time:
    solve(a == L, sol, bcs)
    (u_new, press_new)=sol.split()
    u_old.assign(u_new)
    plot(u_old, title="old velocity", rescale=True)
    u_ref.assign(interpolate(u_old, V_ref))
    plot(u_ref, title="refined velocity")
    t+=dt

---

Yes Mikael, thanks. So, do I understand it right, that
(a) after applying the DirichletBC, FEniCS assembles all mesh facets that are fully contained in the defined region, and
(b) then applies the bc (in case of higher-order polynomial functions), to all "nodes" on every facet?
Peter

---

a) The bcs are applied to matrices and vectors after assembling. All facets take part in the assemble regardless of bcs.

solve(a == L, u, bcs)

is short for

A = assemble(a)
b = assemble(L)
[bc.apply.(A, b) for fc in bcs]
solve(A, u.vector(), b)

b) To all dofs on facets of the boundary. Note there are different kinds of DirichletBCs. See documentation.

---

Thanks, Mikael. So, if I understand it right, if some facet is not completely contained in the bc region (only partially), then no bc is applied to any dofs on this facet... This was probably the core of my question.

---

Like I said, there are different kinds of DirichletBCs. With pointwise you should be able to mark some, but not all dofs on a facet. With default topological this is not possible.

---

I see. Thanks a lot, Mikael!