Understanding integral evaluation in FEniCS

mathematics
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