essential boundary conditions fail when geometric exceeds topological dimension

This code solves the Neumann problem on the unit square using a mixed method. Note that the Neumann problems requires essential boundary conditions (called DirichletBC in FEniCS) for the mixed formulation. First we use a two element mesh of the square setting the geometric dimension and the topological dimension equal to 2. Then we compute the same thing again, but with geometric dimension 3 (viewing the square as a manifold embedded in R^3). The exact solution belongs to the finite element space, so the result should be exact in both cases, but it is completely incorrect when the geometric dimension is greater than the topological dimension. I posted this example earlier if the FEniCS Q&A here.

from dolfin import *
# solve the Neumann problem on a two element mesh on unit square
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])
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.))

Comments (5)