inhomogeneous Neumann boundary condition for mixed Poisson on manifold

Question

This code computes the solution of Poisson equation with Neumann boundary conditions on the hemisphere using the mixed formulation. (The needed mesh file can be downloaded from here.) The code works great if the Neumann condition is homogeneous, but not if it is inhomogeneous. What am I doing wrong?

"""
Douglas N. Arnold  2014-03-03

 This program solves the mixed Poisson equation

   u = -grad p,  -div u + c0 = -f,  int p = 0

 with Neumann (essential) boundary conditions on a hemispherical surface.
 The constant c0 (equal to the mean value of f) is the Lagrange multiplier
 to enforce the side condition that the mean value of p vanishes.
"""

from dolfin import *
from numpy import log

# forcing function
f = Expression("exp(-10.*(pow(x[0], 2) + pow(x[1]-1/sqrt(2.), 2) + pow(x[2]-1/sqrt(2.), 2)))")

mesh = Mesh("hemisphere_mesh.xml")
global_normal = Expression(("x[0]", "x[1]", "x[2]"))
mesh.init_cell_orientations(global_normal)

# finite element mesh and spaces
V0 = FunctionSpace(mesh, "RT", 3)
V1 = FunctionSpace(mesh, "DG", 2)
R = FunctionSpace(mesh, "Real", 0)
V = MixedFunctionSpace([V0, V1, R])
# essential boundary condition
def bdry(x, on_boundary):
    return on_boundary
#   homogeneous BC works
#bc = DirichletBC(V.sub(0), Constant((0., 0., 0.0)), bdry)
# but inhomogeneous does not (examine computed vector field at boundary)
bc = DirichletBC(V.sub(0), Constant((0., 0., 0.1)), bdry)

# trial and test functions
(u, p, c0) = TrialFunctions(V)
(v, q, c1) = TestFunctions(V)

# bilinear form, right hand side
a = (dot(u, v) - p*div(v) - div(u)*q + c0*q + p*c1) * dx
L = -f*q * dx

# assemble, solve
up = Function(V)
solve(a == L, up, bc)
(u, p, c0) = up.split()

up = project(u, VectorFunctionSpace(mesh, "CG", 1))
# plot the scalar variable
plot(up, interactive=True)

Answer 1

You can check that solution fulfills BC:

# assemble, solve
up = Function(V)
solve(a == L, up, bc)
(u, p, c0) = up.split()

x  = up.vector()
x1 = x.copy()
bc.apply(x1)
x1 -= x
print x1.norm('l2'), x1.norm('linf')

(One would rather like to check this using error norm on boundary but I was not able to write compilable form for this as there are issues with RTs, facet normals and manifolds.)

Funny plots are produced either by L2 projection to CG1 (would anything else work without oscillations?) or interpolation of vertex values (if you omit project step).

Comments

Thanks Jan, but I think the problem is worse than you indicate. You seem to imply that it is a problem with projecting the solution for plotting, but as far as I can tell, if the Neumann boundary condition is inhomogeneous, FEniCS computes the solution incorrectly. Here is a simple test case on the hemisphere of the unit sphere with z>0. I take the exact solution to be p(x,y,z)=z-1/2, which has average value zero on the hemisphere, and solve the Neumann problem in mixed form. I get an error of 85%, while if I solve a homogeneous Neumann problem or a homogeneous Dirichlet problem on the same mesh I get errors less than 1%. So I think this may be a bug in FEniCS.

"""
Douglas N. Arnold  2014-03-03

 This program solves the mixed Poisson equation

   u = -grad p,  -div u + c0 = -f  on hemisphere,  u.n = g on equator, int p = 0

 with Neumann (essential) boundary conditions on a hemispherical surface.
 The constant c0 (equal to the mean value of f) is the Lagrange multiplier
 to enforce the side condition that the mean value of p vanishes.
"""

from dolfin import *
from numpy import log
# exact solution p(x, y, z) = z - 1/2 (which has mean value 0 on hemisphere)
pex = Expression("x[2]-.5")
# on equator exact solution for u = -grad p 
ubc = Constant((0., 0., -1.))
# forcing function f = - Lap p
f = Expression("2*x[2]")
mesh = Mesh("hemisphere_mesh.xml.gz")
global_normal = Expression(("x[0]", "x[1]", "x[2]"))
mesh.init_cell_orientations(global_normal)

# finite element mesh and spaces
V0 = FunctionSpace(mesh, "RT", 3)
V1 = FunctionSpace(mesh, "DG", 2)
R = FunctionSpace(mesh, "Real", 0)
V = MixedFunctionSpace([V0, V1, R])
# essential boundary condition
def bdry(x, on_boundary):
    return on_boundary
bc = DirichletBC(V.sub(0), ubc, bdry)

(u, p, c0) = TrialFunctions(V)
(v, q, c1) = TestFunctions(V)
a = (dot(u, v) - p*div(v) - div(u)*q + c0*q + p*c1) * dx
L = -f*q * dx
up = Function(V)
solve(a == L, up, bc)
(u, p, c0) = up.split()
p0 = project(pex, V1)
nrm = norm(p0)
err = errornorm(p0, p)
perct = err/nrm*100
print "p: L2 norm  {},  L2 error {} = {:.3f}%".format(nrm, err, perct)

---

I thought of an even simpler example makes it very clear that there is an error in the solution. This code computes the solution to a Neumann problem with exact solution p(x, y) = x-.5 on the unit square using RT2 * DG1 elements, and, of course, it gets the solution exactly correct, since the exact solution belongs to the FE space. I only use a mesh with 2 elements to make the point.

The second part of the code solves the identical problem, but now the unit square is viewed as a manifold embedded in 3 space. It is essentially identical except for the addition of a third coordinate equal to zero for the mesh vertices. But this code gives the wrong answer. Here is the output:

Unit square in 2D.  Computed value for p(0.3, 0.7) [should be -0.2]: -0.2000000000
Square embedded in 3D.  Computed value for p(0.3, 0.7, 0.0) [should be -0.2]: 0.0066666667

And here is the code:

from dolfin import *
from numpy import log

# compute the solution of Neumann problem on unit square in 2D

pex = Expression("x[0] - .5")
uex = Constant((-1., 0.))
f = Constant(0.)
mesh = Mesh()
me = MeshEditor()
me.open(mesh, 2, 2)
me.init_vertices(4)
me.init_cells(2)
me.add_vertex(0, 0., 0.)
me.add_vertex(1, 1., 0.)
me.add_vertex(2, 0., 1.)
me.add_vertex(3, 1., 1.)
me.add_cell(0, 0, 1, 3)
me.add_cell(1, 0, 2, 3)
me.close()
mesh.order()
V0 = FunctionSpace(mesh, "RT", 2)
V1 = FunctionSpace(mesh, "DG", 1)
R = FunctionSpace(mesh, "Real", 0)
V = MixedFunctionSpace([V0, V1, R])
# essential boundary condition
def bdry(x, on_boundary):
    return on_boundary
bc = DirichletBC(V.sub(0), uex, bdry)
(u, p, c0) = TrialFunctions(V)
(v, q, c1) = TestFunctions(V)
a = (dot(u, v) - p*div(v) - div(u)*q + c0*q + p*c1) * dx
L = -f*q * dx
up = Function(V)
solve(a == L, up, bc)
(u, p, c0) = up.split()
print "Unit square in 2D.  Computed value for p(0.3, 0.7) [should be -0.2]: {:.10f}".format(p(.3, .7))

# Redo the same computation, but view the square as embedded in R^3

pex = Expression("x[0] - .5")
uex = Constant((-1., 0., 0.))
f = Constant(0.)
mesh = Mesh()
me = MeshEditor()
me.open(mesh, 2, 3)
me.init_vertices(4)
me.init_cells(2)
me.add_vertex(0, 0., 0., 0.)
me.add_vertex(1, 1., 0., 0.)
me.add_vertex(2, 0., 1., 0.)
me.add_vertex(3, 1., 1., 0.)
me.add_cell(0, 0, 1, 3)
me.add_cell(1, 0, 2, 3)
me.close()
mesh.order()
global_normal = Expression(("0.", "0.", "1."))
mesh.init_cell_orientations(global_normal)
V0 = FunctionSpace(mesh, "RT", 2)
V1 = FunctionSpace(mesh, "DG", 1)
R = FunctionSpace(mesh, "Real", 0)
V = MixedFunctionSpace([V0, V1, R])
def bdry(x, on_boundary):
    return on_boundary
bc = DirichletBC(V.sub(0), uex, bdry)
(u, p, c0) = TrialFunctions(V)
(v, q, c1) = TestFunctions(V)
a = (dot(u, v) - p*div(v) - div(u)*q + c0*q + p*c1) * dx
L = -f*q * dx
up = Function(V)
solve(a == L, up, bc)
(u, p, c0) = up.split()
print "Square embedded in 3D.  Computed value for p(0.3, 0.7, 0.0) [should be -0.2]: {:.10f}".format(p(.3, .7, 0.))

---

Yes. This is distinct from issue 14 and should be reported. Note that changing representation makes no difference. Test in dolfin/unit/fem/python/manifolds.py can give a clue where the bug is not.

---

I believe this issue is covered by issue 10 which should be reopened.

Answer 2

Nothing. This is/was a rather awkward bug in the application of boundary conditions for spaces that require cell orientations. Should be fixed now (see https://bitbucket.org/fenics-project/ffc/issue/21/essential-boundary-conditions-fail-when)

Comments

Why doesn't the last Doug's example work in parallel? This is strange

Traceback (most recent call last):
  File "test2.py", line 36, in <module>
    solve(a == L, up, bc)
...
*** Error:   Unable to access vector of degrees of freedom.
*** Reason:  Cannot access a non-const vector from a subfunction.
*** Where:   This error was encountered inside Function.cpp.

---

I don't know, I can't immediately reproduce and I don't see what it has to do with the original issue. Please report a separate issue.

---

Ok, sorry. Nevertheless, used initialization of mesh is not correct in parallel which causes the error.