Why does this constraint based approach to the Dirichlet condition work on UnitSquare but not UnitCube?

Question

Hi Fenics Q&A,

I thought I had understood how to apply constraints to a PDE but now I am not so sure. As a test problem I have attempted to implement a Dirichlet condition:

ds: u = u_v 

using a Lagrange multplier.

Following a suggestion elsewhere on this forum I constrain the Lagrange Multiplier field for interior points in the space to be 0 since only the value of the field on the boundary enters the weak form.

This seems to work correctly for a UnitSquareMesh. However, solving the weakly constrained problem does not appear to work for the a UnitCubeMesh (setting dim=3 in the script which follows). Any help explaining my muddled thinking would be very much appreciated!

Regards,

Colin

from dolfin import *


# define problem
dim = 2
u_v = Expression('1+0.1*sin(6*x[0])')
Npts = 11

if dim==2:
    mesh = UnitSquareMesh(Npts,Npts)
    f = Expression('10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)') 
elif dim==3:
    mesh = UnitCubeMesh(Npts,Npts,Npts)
    f = Expression('10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)+ pow(x[2] - 0.5, 2)) / 0.02)') 

V = FunctionSpace(mesh,'Lagrange',2)
W = VectorFunctionSpace(mesh,'Lagrange',2,dim=2)


def is_on_boundary(x,on_boundary):
    return on_boundary    

def not_on_boundary(x,on_boundary):
    return not on_boundary



# weak method: define equation
(u,lm) = TrialFunction(W)
(v,v_lm) = TestFunction(W)   

bc = DirichletBC(W.sub(1), Constant(0), not_on_boundary)    
a = inner(grad(u), grad(v))*dx + lm*v*ds+v_lm*u*ds
L = f*v*dx + v_lm*u_v*ds    

U = Function(W)
solve(a == L, U, bc)

ul,lm=U.split()
plot(ul)


# strong method 
u = TrialFunction(V)
v = TestFunction(V)   

bc = DirichletBC(V, u_v, is_on_boundary)    
a = inner(grad(u), grad(v))*dx 
L = f*v*dx     

u = Function(V)
solve(a == L, u, bc)

plot(u)

Answer 1

Ok never mind, I should have read the manual entry for DirichletBC before posting:

Note: when using “pointwise”, the boolean argument on_boundary in SubDomain::inside will always be false.

The ‘check_midpoint’ variable can be used to decide whether or not the midpoint of each facet should be checked when a user-defined SubDomain is used to define the domain of the boundary condition. By default, midpoints are always checked. Note that this variable may be of importance close to corners, in which case it is sometimes important to check the midpoint to avoid including facets “on the diagonal close” to a corner. This variable is also of importance for curved boundaries (like on a sphere or cylinder), in which case it is important not to check the midpoint which will be located in the interior of a domain defined relative to a radius.