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).
EDIT:
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.