# Understanding integral evaluation in FEniCS

**URL:** <https://fenicsproject.discourse.group/t/understanding-integral-evaluation-in-fenics/8823>\
**Category:** mathematics\
**Created:** [July 21, 2022, 11:39am UTC](https://fenicsproject.discourse.group/t/understanding-integral-evaluation-in-fenics/8823 "2022-07-21T11:39:06Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![glimachaves](https://avatars.discourse-cdn.com/v4/letter/g/4af34b/32.png) [@glimachaves](https://fenicsproject.discourse.group/u/glimachaves)\
**Post date:** [July 21, 2022, 11:39am UTC](https://fenicsproject.discourse.group/t/understanding-integral-evaluation-in-fenics/8823/1 "2022-07-21T11:39:06Z")

</div>

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](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/7/7d040eef8f05aad8cfdd0c93c0ebb01feb96a3eb.png)

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:

```auto
#!/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.

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [July 21, 2022, 1:40pm UTC](https://fenicsproject.discourse.group/t/understanding-integral-evaluation-in-fenics/8823/2 "2022-07-21T13:40:38Z")

</div>

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:

```python
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

> [@glimachaves](#):
>
> ` outfile.write(u, 0.)`

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