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

Main Forum Old-Forum - do not post here
1

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

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
[tex]\left(
\begin{array}{cc}
a & c \
b & d \
\end{array}
\right),[/tex]
nor
[tex]\left(
\begin{array}{cc}
a  & b \
c & d \
\end{array}
\right).[/tex]

Please help me, thank u.

Best wishes,
Di Yang

2

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,
Henry

3

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