This question is mainly an inquiry about the usefulness of Discontinuous Galerkin (DG) for the time-independent transport equation of the form $$\sigma u+\beta\cdot\nabla u =f,\;\;\;\text{on }\Omega\subset\mathbb{R}^2$$ when the advection field $\beta$ is a non-constant function of $x$ (or $u$). Boundary conditions are $$ u=g\;\;\;\text{ on the 'inflow' boundary }\;\;\Gamma_-:=\{x\in\partial\Omega:\beta(x)\cdot\mathbf{n}(x)<0\}. $$ Now let $\mathcal{T}_h$ be a triangulation of $\Omega$ and $V_h=\{v\in L^2(\Omega):v_{|K}\in\mathbb{P}_d(K),\forall K\in\mathcal{T}_h\}$ be the usual DG space. Then the local formulation on each triangle $K\in\mathcal{T}_h$ is given by $$\sigma\int_K u_hv +\int_K (\beta\cdot\nabla u_h)v-\int_K fv=\int_{\partial K^-} (\beta\cdot\mathbf{n}_K)(u_h-\xi_h)v,$$ where $\partial K^-$ is in the inflow boundary and $$\xi_h=\begin{cases}0 &\text{ on }\;\;\partial K^-\cap \Gamma_-,\\ \beta\cdot{n}_K(x)\lim_{\epsilon\to 0^+}u_h(x-\epsilon\beta)&\text{ on }\;\;\partial K^-\setminus\Gamma_-. \end{cases} $$
In the case of when $\beta$ is constant, Lesaint and Raviart [1] showed that there is an enumeration of the triangles following the characteristic direction $\beta$ which renders the matrix of the discrete problem block triangular. I understand the interest of DG over CG (Continuous Galerkin) here, where it is easier to invert the system matrix, thus more or less justifying the extra d.o.f (degrees of freedom).
However, I don't think this advantage remains when $\beta$ becomes non-constant and has loops (if $\beta$ is an electric field for instance). Finding a proper enumeration becomes harder and there is more global coupling. So why should I use DG with its extra d.o.f? Am I missing some fix that makes DG advantageous again?
