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.

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.