Define the inflow and outflow on interior boundaries

Hello everyone, I’m trying to implement an upwind DG scheme, but this requires defining the inflow and outflow on the interior boundaries, and I don’t know how to represent them. This has been bothering me for a long time.

Define mesh Th = square(128, 128, [2*x-1, 2*y-1]); Assume that the velocity \mathbf{u}_{h}^n has already been obtained. uh=[u1old,u2old] with fespace Vh(Th, [P2,P2]); We have \sigma_{h}^{n+1},w_{h}\in W_{h} with fespace Wh(Th, P2dc);

Now I want to reproduce an upwind DG scheme in W_{h}.
\left( \frac{\sigma_{h}^{n+1}-\sigma_{h}^n}{\delta t},w_{h} \right) + \left( \sum_{e}\left(\hat{\sigma}_{h}^{n+1}\mathbf{u}_{h}^n,[[w_{h}]]_{e}\right) - (\sigma_{h}^{n+1}\mathbf{u}_{h}^n,\nabla w_{h})\right) = 0.

For any function w continuous on each K\in \sum, we define
[[w]] = (w^{-} - w^{+})\cdot \mathbf{n}_e
to be the jump of the function w on an edge e=\partial K^{-}\cap \partial K^{+}, with w^{\pm} denoting the trace of w from K^{\pm} onto e, and \mathbf{n}_{e} denoting the normal vector on e pointing towards K^{+}.

The \hat{\sigma}_{h}^{n+1} is the numerical flux, which is a single-valued function defined at the cell interfaces and in general depending on the values of the numerical solution \sigma_{h}^{n+1} from both sides of the interface. Here, the numerical flux on an edge e can be chosen according to the upwind principle, namely,

\hat{\sigma}_{h}^{n+1} = (\sigma_{h}^{n+1})^{-}\quad \text{on }e_{+}^n,\qquad \hat{\sigma}_{h}^{n+1} = (\sigma_{h}^{n+1})^{+}\quad \text{on }e_{-}^n

where e_{-}^n and e_{+}^n denote the numerical inflow and outflow parts of an edge e at time t=t_{n}, defined by
e_{-}^n = \{\mathbf{x}\in e:(\mathbf{u}_{h}^n\cdot \mathbf{n}_{e})(\mathbf{x})\leq 0\},\quad\quad e_{+}^n = \{\mathbf{x}\in e:(\mathbf{u}_{h}^n\cdot \mathbf{n}_{e})(\mathbf{x})> 0\}.

I want to know how to implement the crucial feature in FreeFem++.
\sum_{e}\left(\hat{\sigma}_{h}^{n+1}\mathbf{u}_{h}^n,[[w_{h}]]_{e}\right)