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:

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

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))

def grad(a,b):
    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
mean_gradU = 0.5*(grad(ux, uy) + grad(ux_other, uy_other))
mean_gradV = 0.5*(grad(vx, vy) + grad(vx_other, vy_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


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
a & c \
b & d \
a  & b \
c & d \

Please help me, thank u.

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()

mean_grad_u = 0.5 * (Grad(u) - Grad(u.Other()))
mean_grad_v = 0.5 * (Grad(v) - Grad(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,

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


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


to get the output 2.