normal vectors $n^+$ and $n^-$ point in opposite directions. A characteristic property of the DG method is that computed functions are allowed to be discontinuous across $F$ (hence the “D” in DG). For example, for a function $v$, be it a scalar or a vector, its value on $F$ from $K^+$ is $v^+$, while its value on $F$ from $K^-$ is $v^-$, and these $v^+$ and $v^-$ are not necessarily equal. By contrast, the hybrid variables in the HDG method [*are*]{} single-valued on $F$.[]{data-label="fig:interEdge"}](./traceDem.pdf)

On every face, we define *mean* (*average*) *values* ${\left\{ \cdot \right\}}$ and *jumps* ${\llbracket \cdot \rrbracket}$ as $$\left\{
    \begin{aligned}
      {\left\{ {\mathbf{v}}\right\}}_F  & = \dfrac{1}{2} ({\mathbf{v}}^+ + {\mathbf{v}}^-),
      \\
      {\left\{ v \right\}}_F &  = \dfrac{1}{2}(v^+ + v^-),
      \\
      {\llbracket {\mathbf{n}}\times{\mathbf{v}}\rrbracket}_F &  = {\mathbf{n}}^+\times{\mathbf{v}}^+ +
      {\mathbf{n}}^-\times{\mathbf{v}}^-, \\
       {\llbracket {\mathbf{n}}\cdot{\mathbf{v}}\rrbracket}_F &  = {\mathbf{n}}^+\cdot{\mathbf{v}}^+ +
      {\mathbf{n}}^-\cdot{\mathbf{v}}^-,
      \\
      {\llbracket v{\mathbf{t}}\rrbracket}_F &  = v^+ {\mathbf{t}}^+ + v^- {\mathbf{t}}^-,
    \end{aligned}
  \right.$$ where ${\mathbf{n}}^{\pm}$ denotes the outward unit norm vector to $K^\pm$ and ${\mathbf{t}}^{\pm}$ denotes the unit tangent vectors to the boundaries $\partial    K^{\pm}$ such that ${\mathbf{t}}^+\times{\mathbf{n}}^+=1$ and ${\mathbf{t}}^-\times{\mathbf{n}}^-=1$. For the boundary edges, either on $\partial\Omega$ or on $\partial\Omega_S$, these expressions are modified as $$\left\{
    \begin{aligned}
      {\left\{ {\mathbf{v}}\right\}}_F &  = {\mathbf{v}}^+,
      \\
      {\left\{ v 