Understanding integral evaluation in FEniCS

mathematics
1

Hi all,

I need to evaluate the surface integral of a given function, but I’ve been struggling to understand how exactly FEniCS computes an integral.

More specifically, I’m computing the integral of a scalar function depicted in the figure below.

image
image771×751 60 KB

This is a 2D square domain, and the function is in a DG1 space and is equal to 1 inside a central square of 1 x 1, and zero everywhere else. Naturally, I would expect the integral of this function over the domain to yield 1, the value of the area in which it evaluates to 1. However, the value I’m getting is 2.33, as can be obtained using the following MWE:

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

import dolfin as dlf
import numpy as np

mesh = dlf.RectangleMesh(dlf.Point(-2.5, -2.5, 0), dlf.Point(2.5, 2.5, 0), 10, 10, "left/right")
V = dlf.FunctionSpace(mesh, "DG", 1)
u = dlf.Function(V, name="u")

r_c = 0.5
u_exp = dlf.Expression(
    "fabs(x[0]) <= rc & fabs(x[1]) <= rc ? 1 : 0.",
    rc=r_c,
    element=V.ufl_element(),
)
u.assign(dlf.project(u_exp, V))

dx = dlf.Measure("dx", domain=mesh)
integral = dlf.assemble(u*dx)
print("Integral of u = %.2f" % integral)

with dlf.XDMFFile("outfile.xdmf") as outfile:
    outfile.write(u, 0.)

What am I missing here?

Thanks in advance for your help.

2

If you use a DG 0 function, you will get exactly what you want to achieve, as then each cell has only one value.
However, by using a DG-1 function, you get into some subtle issues when using dlf.Expression.
I would suggest using ufl.conditional and dolfin.SpatialCoordiante as it has the desired effect:

import dolfin as dlf
from ufl import conditional
import numpy as np

mesh = dlf.RectangleMesh(dlf.Point(-2.5, -2.5, 0), dlf.Point(2.5, 2.5, 0), 10, 10, "left/right")
V = dlf.FunctionSpace(mesh, "DG", 1)
u = dlf.Function(V, name="u")

r_c = 0.5
x = dlf.SpatialCoordinate(mesh)
u_exp = conditional(abs(x[0])<=r_c, 1, 0) * conditional(abs(x[1])<=r_c, 1, 0)

u = dlf.project(u_exp, V)

dx = dlf.Measure("dx", domain=mesh)
integral = dlf.assemble(u*dx)
print("Integral of u = %.2f" % integral)

with dlf.XDMFFile("outfile.xdmf") as outfile:
    outfile.write_checkpoint(u, "u", 0.)

Note that I’ve also changed

to use write_checkpoint as write will interpolate the solution to a CG-1 space.

1 Like