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

\(\newcommand{\opdiv}{\operatorname{div}}\)

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

  • Want to use \(H(\opdiv)\)-conforming finite elements for velocity

  • Since \(Q_h = \opdiv V_h\), we obtain exact divergence-free velocities

  • use HDG to discretize viscosity term

  • Since \(u_h \in {\mathcal BDM_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-\widehat{v})_t - \int_{\partial T} \frac{\partial v}{\partial n} (u-\widehat{u})_t + \int_{\partial T} \frac{\alpha p^2}{h} (u-\widehat{u})_t (v-\widehat{v})_t \]
\[ b(u, q) = \int \opdiv \, 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