Turning MeshFunctions into Functions

Question

Hi everyone. In the context of AMR, I want to compute an error estimator of the kind

\[
\eta_K = ||f - g||^2_{L^2(K)}
\]

for each cell \(K\).

I have a CellFunction<std::valarray<double>> that contains values of \(f\) and a VertexFunction<std::valarray<double>> that contains values of \(g\).
I read of the trick of integrating with a DG test function, and this is what I would like to do. To this end, I thought I could interpolate the CellFunction into a DG0 space and the VertexFunction into a CG1 space. Is this possible?

Suppose I manage to do this and eventually obtain the vector with the integrals I want, how do I later match each vector entry to the corresponding cell index?
I will later need to put everything into another CellFunction.

Is there a better way to do this?

Answer 1

If you already have obtained f and g as you described, then the simplest solution is to loop over cells and compute each integral with a 3 point quadrature rule.

See Python example below.

eta = CellFunction("double", mesh)
for cell in cells(mesh):
    a, b, c = cell.entities(0).astype(int)
    eta[cell] += (f[cell] - (1./2)*(g[a] + g[b]))**2
    eta[cell] += (f[cell] - (1./2)*(g[a] + g[c]))**2
    eta[cell] += (f[cell] - (1./2)*(g[b] + g[c]))**2
    eta[cell] *= (cell.volume()/3.)

Comments

Thank you for your hint, I was also considering this opportunity. Can you give me a reference/bibliography for this formula and the associated quadrature error? My final aim is to take this to three dimensions and do the same on a tetrahedron.

---

The formula is exact for quadratic polynomials.

See e.g. G. Dahlquist and Å. Björck, Numerical methods in scientific computing, volume 1, SIAM, 2008. Theorem 5.4.3, page 597.

---

Great, I found the formula and it seems to solve the two-dimensional problem. However, my aim is to apply the technique in 3D. Are you aware of an analogous formula for the three-dimensional case?

---

For example, consider the 3D simplex with vertices (0,0,0), (0,0,1), (0,1,0) and (1,0,0).
Four quadrature points are needed, so choose symmetrically placed points with barycentric coordinates $(t, t, t, 1- 3t)$ and its permutations. The weights must be equal with $w =1/24$. This is a quadrature rule that is exact for linear polynomials, for any choice of $t\in[0,1/3]$.

The rule is exact for quadratic polynomials only for a particular choice of $t$. It's easy to see that $t$ must be be a solution of

$$ 3t^2 + (1-3t)^2 - \frac{2}{5} = 0 $$

The only root in $[0, 1/3]$ is $t = \frac{1}{4}(1- \frac{1}{\sqrt{5}})$.

The dolfin implementation is then

eta = CellFunction("double", mesh)
t = (1-5**-.5)/4
s = 1 - 3 * t
for cell in cells(mesh):
    a, b, c, d = cell.entities(0).astype(int)
    eta[cell] += (f[cell] - t*(g[a] + g[b] + g[c]) - s*g[d])**2
    eta[cell] += (f[cell] - t*(g[a] + g[b] + g[d]) - s*g[c])**2
    eta[cell] += (f[cell] - t*(g[a] + g[c] + g[d]) - s*g[b])**2
    eta[cell] += (f[cell] - t*(g[b] + g[c] + g[d]) - s*g[a])**2
    eta[cell] *= (1./4) * cell.volume()

The weight here is scaled with the determinant of the Jacobian.

---

Thanks for the detailed explanation! I have been traveling in the past few days, I will get back to you here very soon :)