# About the facet integral in the DG methods for solving the Stokes equation

To solve the Stokes equation, we use the DG methods. And the following formula is an important facet integral:

$$\sum_{F\in\mathcal F_h}\oint_F\left{\nabla\bm{u}\right}\bm{n}_F\cdot \left[\bm{v}\right]\mathrm{d}F.$$

If coding in the NGSolve, I try this:

def outer(W,n):
return CoefficientFunction((W[0]*n[0], W[0]*n[1], W[1]*n[0], W[1]*n[1]), dims=(2,2))

return CoefficientFunction((a.Diff(x), b.Diff(x), a.Diff(y), b.Diff(y)), dims=(2,2))

fes = L2(mesh, order=order_p, dim=2, dgjumps=True)
U, V = fes.TnT()
ux, uy = U
vx, vy = V
n = specialcf.normal(2)
ux_other, uy_other = U.Other()
vx_other, vy_other = V.Other()

U_vec = CoefficientFunction((ux, uy))
U_Other = CoefficientFunction((ux_other, uy_other))
V_vec = CoefficientFunction((vx, vy))
V_Other = CoefficientFunction((vx_other, vy_other))

jump_U = U_vec - U_Other
jump_V = V_vec - V_Other

a = BilinearForm(fes)
a += SymbolicBFI( InnerProduct( outer(jump_V, n), mean_gradU ),\
VOL, skeleton=True )
a += SymbolicBFI( InnerProduct( outer(V_vec, n), grad(ux,uy) ),\
BND, skeleton=True )

I cannot sure it is correct because I have no idea that the normal

n=specialcf.normal(2)

is a column vector nor a row vector.

In addition, if one codes like

CoefficientFunction((a,b,c,d), dims=(2,2))

the corresponding matrix is
$$\left( \begin{array}{cc} a & c \ b & d \ \end{array} \right),$$
nor
$$\left( \begin{array}{cc} a  & b \ c & d \ \end{array} \right).$$

Best wishes,
Di Yang

Hi Di Yang,

I think that the problem is that youâ€™re using the wrong FESpace. VectorL2 should be the one you need:

from ngsolve import *
from netgen.geom2d import unit_square

mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))
V = VectorL2(mesh, order=3, dgjumps=True)
u, v = V.TnT()
n = specialcf.normal(mesh.dim)

jump_u = u - u.Other()
jump_v = v - v.Other()

a = BilinearForm(V)
a += InnerProduct(mean_grad_u * n, jump_v) * dx(skeleton=True)
a += InnerProduct(Grad(u) * n, v) * ds(skeleton=True)

Best wishes,
Henry

Hi Di Yang,

concerning your question about coefficient functions, the input is taken row-wise. For

cf = CoefficientFunction((1,2,3,4), dims=(2,2))

you can get out the components by

cf[i,j]

and I think evaluating (or drawing) the components is probably the quickest way to check the alignment.
Assume you have a mesh of the unit square, you could call

cf[0,1](mesh(0.5,0.5))

to get the output 2.

Best,
Christoph