Hdiv for Stokes, error at boundary


I am trying to implement the Stokes problem for Discontinuous Galerkin (and eventually HDG) H(div) elements. The code is posted below. I do not however get the correct order of convergence. When I plot the norm of the error it is clear that something is wrong with the boundary integrals, but I cannot figure out what. I have attached the plot of the norm of the error.

from ngsolve import *
from ngsolve.internal import *

basic xfem functionality

basic geometry features (for the background mesh)

from netgen.geom2d import SplineGeometry
import numpy as np
import netgen.gui

from numpy import pi

import scipy
import scipy.sparse as sp

geo = SplineGeometry()
geo.AddRectangle( (0, 0), (0.99, 0.99), bcs = (“top”, “left”, “bottom”, “right”))

nu = 1
uex = CoefficientFunction((sin(piy), -cos(pix)))
force = CoefficientFunction((-pisin(pix)cos(piy) + nu * pi2sin(piy), -pisin(piy)cos(pix) - nu * pi2cos(pix)))
pex = CoefficientFunction(cos(pi*x)cos(piy))

l22_list = []

def l2error(refinement = 0):
h0 = 0.1
order = 2
alpha = 10

mesh = Mesh(geo.GenerateMesh(maxh = h0*2**(-refinement), quad_dominated = False))

Wh = HDiv(mesh, order = order, dirichlet="top|left|bottom|right", dgjumps = True BDM = True) #Velocity space
Qh = L2(mesh, order = 0) #Pressure space

Vh = FESpace([Wh,Qh], dgjumps = True)

u,p = Vh.TrialFunction()
v,q = Vh.TestFunction()

jump_u = u-u.Other()
jump_v = v-v.Other()

n = specialcf.normal(2)
h = specialcf.mesh_size

flux_u = 0.5*(grad(u) + grad(u.Other()))*n 
flux_v = 0.5*(grad(v) + grad(v.Other()))*n 

a = BilinearForm(Vh)

#Volume integrals
a += SymbolicBFI(nu*InnerProduct(grad(u), grad(v)))
a += SymbolicBFI(-div(u)*q) #consistent symmetrizing term
a += SymbolicBFI(-div(v)*p) 

#Interior facet integrals
a += SymbolicBFI(-nu*InnerProduct(flux_u,jump_v), VOL, skeleton = True)
a += SymbolicBFI(-nu*InnerProduct(flux_v,jump_u), VOL, skeleton = True) #consistent symmetrizing term
a += SymbolicBFI(h**(-1) * alpha * order**2 * InnerProduct(jump_u,jump_v), VOL, skeleton=True) #consistent stabilizing term

#Boundary facet integrals

a += SymbolicBFI(-nu*InnerProduct(grad(u)*n,v), BND, skeleton = True)
a += SymbolicBFI(-nu*InnerProduct(grad(v)*n,u), BND, skeleton = True) #symmetrizing term
a += SymbolicBFI(h**(-1) * alpha * order**2 * InnerProduct(u,v), BND, skeleton=True) #stabilizing term
a += SymbolicBFI(v*n*p, BND, skeleton=True)
a += SymbolicBFI(u*n*q, BND, skeleton=True) #symmetrizing term 

#Right-hand side

f = LinearForm(Vh)
f += SymbolicLFI(force*v)
f += SymbolicLFI(-nu*InnerProduct(grad(v)*n,uex),BND, skeleton = True) #symmetrizing term
f += SymbolicLFI(h**(-1) * alpha * order**2 * InnerProduct(uex,v), BND, skeleton=True) #stabilizing term
f += SymbolicLFI(uex*n*q, BND, skeleton = True) #symmetrizing term


gfu = GridFunction(Vh)

gfuex = GridFunction(Wh)
gfpex = GridFunction(Qh)

a_inv = a.mat.Inverse(freedofs = Vh.FreeDofs(), inverse="umfpack")

gfu.vec.data = a_inv*f.vec
vel = gfu.components[0]
pres = gfu.components[1]
#Draw(pres, mesh, "p")
#Draw(vel, mesh, "u")
#Draw(div(vel), mesh, "div(u)")
#Draw(Norm(vel), mesh, "|u|")

dv = grad(vel)
dg = grad(gfuex)

#eu = sum((dv-dg)**2)
eu = (vel - gfuex)**2

err = vel-gfuex

Draw(Norm(err), mesh, "|e|")

l22 = Integrate(eu, mesh)
print("h = ", h0*2**(-refinement),",  L2^2 error in velocity: ", l22)

if len(l22_list)>1:
    print("EOC: ", (np.log(sqrt(l22_list[-2]))-(np.log(sqrt(l22_list[-1]))))/np.log(2))

for i in range(2,3):

Best regards,


Hi Alvar,

You use Dirichlet-Flags for the Hdiv-FESpace, but don’t set the Dirichlet-Values. Hence, your discrete solution will have zero values at the dirichlet dofs of you Hdiv space, i.e. u·n is zero on all boundaries. Your Nitsche-type terms are not acting here as u·n and v·n vanishes with the Dirichlet-Flags. Either remove the dirichlet-flags from the space or (preferably) use gfu.components[0].Set(uex) to set the correct Dirichlet values.

· You may want to use “Grad” instead of “grad” (the latter is not the Jacobi matrix for hdiv elements).
· As you are normal-continuous, you may want to add a tangential projection for all jump and boundary terms of u and v.