Boundary conditions are not preserved under projections onto FE-spaces

Question

I have troubles in understanding how fenics applies Dirichlet boundary conditions.

In particular, I have a function u on a domain D and I know the Dirichlet values g of u at the boundary. Say I project u onto a FE-space V to get uh, then I get a difference between uh and uh after having imposed the Dirichlet data.

|uh - bc.apply(uh)| ~ 1e-7

I understand, that an L2-projection does not preserve function values at the nodes but

• the error should get smaller with finer discretizations
• I expect a zero error for linear u and piecewise linear finite elements

which is not the case.

What am I doing wrong? Maybe it's just a conceptual misunderstanding.

Below is my dolfin code to reproduce my statements

import numpy as np
from dolfin import *

mesh = UnitSquareMesh(64, 64)
V = VectorFunctionSpace(mesh, "CG", 1)
u = Expression(('x[1]','0'))

ufun = project(u, V)
ufunbc = project(u, V)

# definition of the boundary
def boundary(x, on_boundary): 
    return (x[1] > 1.0 - DOLFIN_EPS 
        or x[0] > 1.0 - DOLFIN_EPS 
        or x[1] < DOLFIN_EPS 
        or x[0] < DOLFIN_EPS)

# take the function values as dirichlet data
diridata = u 
bc = DirichletBC(V, diridata, boundary)

# apply boundary conditions
bc.apply(ufunbc.vector())

plot(ufun-ufunbc)

print np.linalg.norm(ufun.vector().array() - ufunbc.vector().array())

print np.allclose(ufun.vector().array(), ufunbc.vector().array())

Answer 1

Try using a tighter tolerance in the projection:

solve(inner(TestFunction(V), TrialFunction(V)) * dx == inner(TestFunction(V), u) * dx, ufun,
  solver_parameters = {"linear_solver":"cg", "preconditioner":"sor", "krylov_solver":{"absolute_tolerance":1.0e-16, "relative_tolerance":1.0e-16}})

or a direct solver:

ufun = project(u, V, solver_type = "lu")

Comments

That solved it. Thanks a lot!