H(\operatorname{div})-conforming HDG methods for Stokes

4. H( $div$ )-conforming HDG methods for Stokes#

Want to use $H (div)$ -conforming finite elements for velocity
Since $Q_{h} = div V_{h}$ , we obtain exact divergence-free velocities
use HDG to discretize viscosity term
Since $u_{h} \in B D M_{k}$ is normally continuous, HDG is needed only for the tangential component
Facet variable $\hat{u}$ lives in the tangential space on the facet

a (u, \hat{u}; v, \hat{v}) = \sum_{T} \int_{T} \nabla u \nabla v - \int_{\partial T} \frac{\partial u}{\partial n} (v - \hat{v})_{t} - \int_{\partial T} \frac{\partial v}{\partial n} (u - \hat{u})_{t} + \int_{\partial T} \frac{α p^{2}}{h} (u - \hat{u})_{t} (v - \hat{v})_{t}

b (u, q) = \int div u q

from netgen.geom2d import SplineGeometry
from ngsolve import *
from ngsolve.webgui import Draw

geo = SplineGeometry()
geo.AddRectangle( (0, 0), (2, 0.41), bcs = ("wall", "outlet", "wall", "inlet"))
geo.AddCircle ( (0.2, 0.2), r=0.05, leftdomain=0, rightdomain=1, bc="cyl")
mesh = Mesh( geo.GenerateMesh(maxh=0.08))
mesh.Curve(5);
Draw (mesh)

BaseWebGuiScene

order=4
VT = HDiv(mesh, order=order, dirichlet="wall|cyl|inlet")
VF = TangentialFacetFESpace(mesh, order=order, dirichlet="wall|cyl|inlet")
Q = L2(mesh, order=order-1)
V = VT*VF*Q

u, uhat, p = V.TrialFunction()
v, vhat, q = V.TestFunction()

n = specialcf.normal(mesh.dim)
h = specialcf.mesh_size
dS = dx(element_boundary=True)

nu = 1e-3

def tang(vec):
    return vec - (vec*n)*n

# Thesis Christoph Lehrenfeld, page 71
a = BilinearForm (V)
a += nu*InnerProduct(Grad(u), Grad(v)) * dx
a += nu*InnerProduct(Grad(u)*n, tang(vhat-v)) * dS
a += nu*InnerProduct(Grad(v)*n, tang(uhat-u)) * dS
a += nu*4*order*order/h * InnerProduct(tang(vhat-v), tang(uhat-u)) * dS
a += div(u)*q*dx + div(v)*p*dx
a.Assemble();

f = LinearForm(V).Assemble()

invstokes = a.mat.Inverse(V.FreeDofs())

gfu = GridFunction(V)
uin = CF( (1.5*4*y*(0.41-y)/(0.41*0.41),0))
gfu.components[0].Set(uin, definedon=mesh.Boundaries("inlet"))

res = f.vec - a.mat*gfu.vec
gfu.vec.data += invstokes * res

Draw (gfu.components[0], mesh, vectors={"grid_size":30})  # velocity
Draw (gfu.components[2], mesh); # pressure