I'm studying an error estimator for the equation $\nabla\cdot(\beta u) + cu = f$ and it contains the following term $$||f - cU_h - \Pi(f-c U_h) ||_T$$ where :

  • $\Pi$ is the local orthogonal $L^2$ projection onto the classical polynomial DG finite element space $\{v \in L^2: v_{|K} \in \mathbb{P}^{r}(K) \}$
  • $U_h$ is the finite element solution, which I have on each cell $T$

I'm implementing this in deal.II. Morally what I want to achieve is the computation of that integral over each cell using quadratures, so

$$\sum_{q} \Bigl(f(x_q) - c U_h(x_q) - \Pi(f-cU_h)(x_q) \Bigr) w_q \det(B_T)$$

I recall that $U_h$ is a vector. Luckily, there's already a routine that allows me to compute $U_h(x_q)$, that is named get_function_values()

However, the true issue is the term $\Pi(f-cU_h)(x_q)$, i.e. the projection of that difference, evaluated at the quadrature points.

What is the practical way people follow to compute that kind of expression? I wasn't able to find anything useful on the deal.II documentation (probably I didn't look in the good places).


I tried to follow the second suggestion given by prof. Bangerth. By definition of $L^2$ projection: $$(\Pi(f - cU_h), \phi_i) = (f-cU_h, \phi_i)$$ for every test function $\phi_i$ in my dg space. This means that on every cell I can define and solve $$M_{ij} x_j = (f - c U_h, \phi_i):=F_i$$ where $M_{ij}=(\phi_j,\phi_i)$ is the mass matrix corresponding to a generic cell.

The r.h.s of the equation above is computed by quadrature using the values of $U_h$ at quadrature points (given by the get_function_values() routine), while of course $f$ is a known continuous function evaluated as usual.

I was implementing this by using a MeshWorker framework as done in step-74. That means that the linear system above was put inside the cell_worker lambda expression so that, while looping over all the cells, I compute that projection $x$ for every cell and use that to compute the actual error indicator on each cell.


1 Answer 1


You need to compute the projection $\Pi(f-cU_h)$ as a first step. This is a finite element field, so you can create a global field (in deal.II lingo: Create the corresponding finite element and a DoFHandler object) and compute it by global projection by inverting a mass matrix. Alternatively, because the field is discontinuous, you can also forgo the global field and just compute the projection locally by inverting a local mass matrix. In either case, you end up with a finite element field on the current cell associated with the discontinuous element, which you can then evaluate at quadrature points like any other finite element field.

  • $\begingroup$ Thanks Prof. Bangerth. I edited my question with what I did in these two days: I didn't create any new DoFHandler object to do this, I used the one I was using in the assembly routine. Does it sound okay to you? @WolfgangBangerth $\endgroup$
    – FEGirl
    Commented Aug 22, 2021 at 16:34
  • $\begingroup$ @bobinthebox, this answer is great, but I would like to add something. Depending on the polynomial degree of the solution and $f$, you have to use a quadrature rule with higher degree to avoid introducing truncation errors (from the quadrature). Most people use $2p+1$-st degree quadrature rules. But I found that (experientially) even for a simple nice function like $\sin(\pi x)\cos(\pi x)$ over the unit square, $2p+1$ is not enough. So I use $2p+3$-rd degree or higher quadrature rules. $\endgroup$ Commented Aug 22, 2021 at 17:30
  • $\begingroup$ @bobinthebox It seems reasonable what you describe. $\endgroup$ Commented Aug 22, 2021 at 19:30
  • $\begingroup$ Thanks a lot for your check. My estimator is almost completed. I'm trying to implement the estimator proposed at equation (4.7) of this paper (sci-hub.st/10.1007/978-3-319-01818-8_8). Unfortunately, this involves a reconstruction $\hat{u}$ which is defined elementwise. I've never done such a "reconstruction" and, if I am not mistaken, there's no tutorial which does a similar thing, right? @WolfgangBangerth $\endgroup$
    – FEGirl
    Commented Aug 22, 2021 at 22:17
  • $\begingroup$ There is no tutorial that shows exactly this, but you should take a look at the HDG and WG tutorials (step-51 and 62) to see similar techniques at play. $\endgroup$ Commented Aug 23, 2021 at 2:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.