Say I have a simple elliptic PDE: $$ -\nabla\cdot(K\nabla p) = f \;\;\;\text{in}\;\Omega $$
with the appropriate boundary conditions. I solve for $p$ using a FEM (a discontinuous Galerkin method to be exact). Let's denote this as $p^{DG}$.
Now I compute the velocity $\mathbf{u}^{DG} = -K\nabla p^{DG}$. This velocity will be used in the solution of a transport equation. Using this computed velocity is causing oscillations and instability in my solution of the transport equation, and a common technique in the area I'm working in right now is to project the velocity into the $H(div)$ space so that the normal component is continuous.
I've been reading up on this technique and am following the paper "Superconvergence and H(div) projection for discontinuous Galerkin methods." (Riviere and Bastian).
I think I understand conceptually how to do the projection, but am getting caught up on some of the details.
The paper states as follows: Say we have a regular triangulation $\mathcal{T}_{h}$ of the domain $\Omega$. For each $E\in\mathcal{T}_{h}$, we seek a reconstructed velocity/flux $\mathbf{u}^{*} \in BDM_{k-1}(E)$ (Brezzi, Douglas, Marini space) s.t.
$\int_{e}\mathbf{u}^{*}\cdot\mathbf{n}_{e}z = \int_{e}\mathbf{u}^{DG}\cdot\mathbf{n}_{e}z, \;\;\;\;\;\;\;z\in\mathbb{P}_{k-1}(e), \;\forall e\in\partial E $
$\int_{E}\mathbf{u}^{*}\cdot\nabla w = \int_{e}\mathbf{u}^{DG}\cdot\nabla w, \;\;\;\;\;\;w\in\mathbb{P}_{k-2}(E)$
$\int_{E}\mathbf{u}^{*}\cdot\mathbf{S}(\phi) = \int_{E}\mathbf{u}^{DG}\cdot\mathbf{S}(\phi), \;\;\;\phi\in M_{k}(E) = \left\{\phi \in \mathbb{P}_{k}(E): \phi|_{\partial E} = 0\right\} \;\text{and} \;\mathbf{S}(\phi) = (\phi_{y}, -\phi_{x})$
1-3 should uniquely define $\mathbf{u}^{*}$. So on each element, I need to solve a small linear system. I am familiar with $L^{2}$ projections, which I think are much more straightforward. The main thing I am struggling with here is how to generate all the of the necessary equations.
For an example, say I have $k=3$. The $BDM_{2}(E)$ space has dimension 12. So I need to generate a 12x12 system.
Looking at (2), then $z \in \mathbb{P}_{1}(E)$. Since $E$ is triangular, I have exactly three basis functions and thus, I have 3 unique equations.
For (1), I am a bit confused. Since $k=3$, we have $w\in\mathbb{P}_{2}$. For discontinuous Galerkin (and CG), there are typically 6 quadratic polynomials on each element. So does this mean I get 6 equations from (1)? Do I just assign 2 of the polynomials arbitrarily to each edge...?
I understand (3) mostly. I need to generate some cubic bubble functions, but since I don't know how many equations I get from (1), I'm not sure how many bubble functions I need.
Can someone shed some light on this topic? Is the method above a "standard" way of doing it? There is a surprisingly small amount of resources available on this, which is probably the most frustrating part so far and I have worked much with $H(div)$ elements before. If anyone knows of some relevant sources here, I'd also be really grateful.