Projection of scalar field on DG0 and DG1 on boundary of mesh

Question

I want to calculate a scalar product of the mesh normal and a velocity field on the boundar of a mesh:

v_Expression = Expression(('0.0','-1'), degree=2, domain=mesh)

v = interpolate(v_Expression, SomeVectorSpace)
N = FacetNormal(mesh)

dot(v,N)

I have tried different vectorspaces: DG0 and DG1. In the image you can see the results.

For DG0 (top image): it is working correctly.
For DG1 (bottom image) : i dont understand the results. Some values are dropped or become zero.

Some values become zero and the absolute value is divided by two (should be 1, but is 0.5). Also some non-zero values appear in the domain, when i use a more complex geometry meshed with gmesh.
Is this a bug? Is there a way to handle this?

Here is the working example:

from fenics import *

mesh = RectangleMesh(Point(-1,0), Point(1, 1), 10, 10, "right/left")

dA = Measure("ds", domain=mesh)
dV = Measure("dx", domain=mesh)

N = FacetNormal(mesh)
h = FacetArea(mesh)

pwd = './results/'

v_Expression = Expression(('0.0','-1'), degree=2, domain=mesh)

###################################################################
###################################################################

file_v_dot_n_DG0 = File(pwd + 'v_dot_n_DG0.pvd')

ScalarSpaceDG0 = FunctionSpace(mesh,'DG',0)
VectorSpaceDG0 = VectorFunctionSpace(mesh,'DG',0)


v_DG0 = interpolate(v_Expression, VectorSpaceDG0)
TestFunctionDG0 = TestFunction(ScalarSpaceDG0)
v_dot_n_DG0_array = assemble( dot(v_DG0,N)/h*TestFunctionDG0*dA).array()


v_dot_n_DG0 = Function(ScalarSpaceDG0)
v_dot_n_DG0.vector()[:] = v_dot_n_DG0_array


file_v_dot_n_DG0 << (v_dot_n_DG0)

###################################################################
###################################################################

file_v_dot_n_DG1 = File(pwd + 'v_dot_n_DG1.pvd')

ScalarSpaceDG1 = FunctionSpace(mesh,'DG',1)
VectorSpaceDG1 = VectorFunctionSpace(mesh,'DG',1)

v_DG1 = interpolate(v_Expression, VectorSpaceDG1)


TestFunctionDG1 = TestFunction(ScalarSpaceDG1)
v_dot_n_DG1_array = assemble( dot(v_DG1,N)/h*TestFunctionDG1*dA).array()


v_dot_n_DG1 = Function(ScalarSpaceDG1)
v_dot_n_DG1.vector()[:] = v_dot_n_DG1_array


file_v_dot_n_DG1 << (v_dot_n_DG1)

Second question: In the end i want to project the results on CG1. Is there a easy way to calculate the mean values on the node between two DG elements? Or are there any other ways?

file_v_dot_n_DG0_projection_on_CG1 = File(pwd + 'v_dot_n_DG0_projection_on_CG1.pvd')

ScalarSpaceCG1 = FunctionSpace(mesh,'CG',1)

v_dot_n_DG0_projection_on_CG1 = project(v_dot_n_DG0, ScalarSpaceCG1)

file_v_dot_n_DG0_projection_on_CG1 << (v_dot_n_DG0_projection_on_CG1)

Answer 1

Unfortunately I do not think you understand your problem. You are attempting to project a function which is discontinuous on elements, onto a function defined continuous on each element. Recall the DG space of functions:

$$ V^{DG} := \lbrace v_h \in L_2(\Omega) : \left. v_h \right|_{\kappa} \in \mathcal{P}^p (\kappa), \kappa \in \mathcal{T}_h \rbrace $$

Your code is solving the problem: find $u_h \in V^{DG}$ such that

$$ (u_h, v_h) = \langle b \cdot n, v_h \rangle_{\partial \Omega} \quad \forall v_h \in V^{DG}. $$

If your space is DG0 and you have a piecewise constant approximation, the boundary terms of the DG formulation vanish, and you recover your cell-wise constant result as above.

If you solve for any polynomial approximation $p \geq 1$ then your problem is not well posed and you are trying to project a function which is discontinuous on each element onto the element and you get large numerical errors. Try for example setting the polynomial order of your DG approximation space to 2, 3 and 4.

Consider the following example for projecting $b \cdot n$ onto a CG space. Bear in mind that this problem is also ilposed as $n$ is discontinuous at the vertices of the square mesh.

from dolfin import *

mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "CG", 1, restriction="facet")

vel = Expression(("0.0", "-1.0"), degree=4, domain=mesh)
n = FacetNormal(mesh)

du, v = TrialFunction(V), TestFunction(V)
u = Function(V)
a = du*v*ds
L = dot(vel, n)*v*ds

solve(a == L, u)

plot(u, interactive=True)
File("out.pvd") << u