Different results surface integration depending on elemtype (with working example)

Question

Dear Fenics community,

I have a question regarding the integration over mesh boundaries, that seems to result in different solutions depending on the finite element type I choose. The context of this question is another question I asked recently here, in which I want to compute the wall shear stress. A key step in this computation is integration over the boundary domain and depending on the element type I choose I get different answers: for instance if I choose first order continuous Galerkin elements my computations are a factor 2 too high.

I can illustrate this by a working example at the bottom of this post, where I integrate a constant vector function (1., 1., 1.) over the boundary domain and divide by the FacetArea. I would imagine that irrespective of the element type this would again result in the vector (1., 1., 1.), but this only seems to work for zeroth-order discontinous Galerkin, whereas, for instance first-order continuous Galerkin results in a constant vector (2., 2., 2.). In addition, I observe some discrepancies on the edges of the mesh. Can anybody explain to me the reason why different element types result in different solutions? I do not mind having to multiply a solution by a certain factor to get the correct answer (for instance the first-order Galerkin by a factor 0.5), but I like to know why I should do so, to avoid making incorrect assumptions. Thanks in advance

#!/usr/bin/env python
# -*- coding: utf-8 -*-

"""
Minimal working example surface integration question
"""

# Import statements
from fenics import *

# Create a 3D cube mesh
meshsize = 10
mesh = UnitCubeMesh(meshsize, meshsize, meshsize)

# Define function spaces
V0 = VectorFunctionSpace(mesh, 'DG', 0)  # Discontinous Galerkin const.
V1 = VectorFunctionSpace(mesh, 'DG', 1)  # Discontinous Galerkin 1-order
V2 = VectorFunctionSpace(mesh, 'DG', 2)  # Discontinous Galerkin 2-order
VC1 = VectorFunctionSpace(mesh, 'CG', 1)  # Continous Galerkin 2-order
VC2 = VectorFunctionSpace(mesh, 'CG', 2)  # Continous Galerkin 2-order

# Define the test functions
v0 = TestFunction(V0)
v1 = TestFunction(V1)
v2 = TestFunction(V2)
vc1 = TestFunction(VC1)
vc2 = TestFunction(VC2)

# Define functions to be integrated
u_0 = Function(V0)
u_0 = Constant((1., 1., 1.))
u_00 = Function(V0)
u_1 = Function(V1)
u_1 = Constant((1., 1., 1.))
u_11 = Function(V1)
u_2 = Function(V2)
u_2 = Constant((1., 1., 1.))
u_22 = Function(V2)
u_c1 = Function(VC1)
u_c1 = Constant((1., 1., 1.))
u_c11 = Function(VC1)
u_c2 = Function(VC2)
u_c2 = Constant((1., 1., 1.))
u_c22 = Function(VC2)


# Integrate the functions elementwise
assemble((1 / FacetArea(mesh)) * inner(v0, u_0) * ds, tensor=u_00.vector())
assemble((1 / FacetArea(mesh)) * inner(v1, u_1) * ds, tensor=u_11.vector())
assemble((1 / FacetArea(mesh)) * inner(v2, u_2) * ds, tensor=u_22.vector())
assemble((1 / FacetArea(mesh)) * inner(vc1, u_c1) * ds, tensor=u_c11.vector())
assemble((1 / FacetArea(mesh)) * inner(vc2, u_c2) * ds, tensor=u_c22.vector())
plot(u_00, title='DG 0 vec')
# plot(u_11, title='DG 1 vec')
# plot(u_22, title='DG 2 vec')
plot(u_c11, title='CG 1 vec')
# plot(u_c22, title='CG 2 vec')

# Only one component
plot(u_00[0], title='DG 0 x-component')
# plot(u_11[0], title='DG 1 x-component')
# plot(u_22[0], title='DG 2 x-component')
plot(u_c11[0], title='CG 1 x-component')
# plot(u_c22[0], title='CG 2 x-component')
interactive()

Answer 1

Hi, you seem to be after a function in V that represents f=Constant((1, 1, 1)) in the space. In general this function is not represented by a vector you assemble above. First of all the problem: Find u in V such that inner(u, v)*ds == inner(f, v)*ds (v being test functions from V), is illposed. However, for V a DG0 space assemble(1/FacetArea(mesh)*inner(f, v)*ds can be viewed as applying the pseudoinverse of the ds-mass matrix, i.e. assemble(inner(TrialFunction(V), v)*ds) to vector assemble(inner(f, v)*ds). Note that the mass matrix is diagonal with nonzero entries being the facet areas. In summary, for the other elements you are not projecting correctly.

Comments

Dear Mirok,

Thanks for your answer! This clarifies why I do not get the anticipated result for function spaces other than DG0 (and probably why I do not get the right wall shear stress calculations using this method, as this was my original problem).

I have one remaining question: I agree with you that I am searching for a vector u in V such that (u, v)ds == inner(f, v)ds. However, you state that this problem is ill posed. How would you propose to obtain this, u, vector? By, for instance, setting a Dirichlet condition for everything but the boundary and then call a solve(A,b,u_)? Is it possible do set such a Dirichlet condition, as it isn't specified on the boundary?

Sjeng

---

Ok, I've worked out that setting a Dirichlet BC = Constant((0., 0., 0.)) for all internal nodes is probably the way to go. However, I am unable to do this. Using the example given in https://fenicsproject.org/qa/830/interior-boundary-conditions I end up with a FacetFunction (facet_domains) that is 0 on the interior and 1 on the exterior.

However, DirichletBC(V, u_in, facet_domains, 0) applies the boundary condition Constant((0., 0., 0.)) to all vertices that belong to any element that has one or more nodes on the interior. Is there another way to do this? (See working code that I borrowed and modified from the above mentioned post below)

#!/usr/bin/env python
# -*- coding: utf-8 -*-

"""
Test to see if we can prescribe Dirichlet BC's on the internal domain
"""

# Import statements
from dolfin import *

def inside(x, on_boundary):
    print(on_boundary)
    # return x[0] > DOLFIN_EPS and abs(x[0] - 1.0) > DOLFIN_EPS \
    #     and x[1] > DOLFIN_EPS and abs(x[1] - 1.0) > DOLFIN_EPS
    return(True)


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

V = FunctionSpace(mesh, "CG", 1)

# initialize mesh connectivity so that Facet.exterior() works
mesh.init()

# define the interior of the domain by looking at each facets
facet_domains = FacetFunction('size_t', mesh)
facet_domains.set_all(0)
for f in facets(mesh):
    if any(ff.exterior() for ff in facets(f)):
        facet_domains[f] = 1

u_in = Constant(-0.0)

bc_in = DirichletBC(V, u_in, facet_domains, 0)
print(dir(bc_in))
print(bc_in.markers())
print(len(bc_in.get_boundary_values()))
print(sum(facet_domains.array()))
print(dir(facet_domains))

u = Function(V)
u.vector()[:] = 0.5

bc_in.apply(u.vector())
plot(facet_domains)
plot(u, axes=True)
interactive()

---

Hi, consider the following

from dolfin import *

mesh = UnitSquareMesh(32, 32)
V = VectorFunctionSpace(mesh, 'CG', 1)
u = TrialFunction(V)
v = TestFunction(V)
f = Constant((1, 1))

a = inner(u, v)*ds
L = inner(f, v)*ds
A = assemble(a, keep_diagonal=True)
b = assemble(L)
A.ident_zeros()

uh = Function(V)
solve(A, uh.vector(), b)

plot(uh, interactive=True)

The important piece here is ident_zeros which makes the matrix A invertible.

---

Thanks! I did not know of the existence of this function, but it seems to work for me.