4. H( )-conforming HDG methods for Stokes#
Want to use
-conforming finite elements for velocitySince
, we obtain exact divergence-free velocitiesuse HDG to discretize viscosity term
Since
is normally continuous, HDG is needed only for the tangential componentFacet variable
lives in the tangential space on the facet
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