## We have moved the forum to https://forum.ngsolve.org . This is an archived version of the topics until 05/05/23. All the topics were moved to the new forum and conversations can be continued there. This forum is just kept as legacy to not invalidate old links. If you want to continue a conversation just look for the topic in the new forum.

#### Notice

The forum is in read only mode.

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

2 years 8 months ago #3632
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:
Code:
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
Code:
n=specialcf.normal(2)
is a column vector nor a row vector.

In addition, if one codes like
Code:
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
2 years 8 months ago #3634 by hvwahl
Hi Di Yang,

I think that the problem is that you're using the wrong FESpace. VectorL2 should be the one you need:
Code:
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
The following user(s) said Thank You: Younghigh
2 years 8 months ago #3638 by cwinters
Hi Di Yang,

concerning your question about coefficient functions, the input is taken row-wise. For
Code:
cf = CoefficientFunction((1,2,3,4), dims=(2,2))
you can get out the components by
Code:
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
Code:
cf[0,1](mesh(0.5,0.5))
to get the output 2.

Best,
Christoph
The following user(s) said Thank You: Younghigh
Time to create page: 0.133 seconds