Do-nothing boundary for Stokes behaving weirdly


I am implementing Stokes equation on surfaces, from Lehrenfeld-Lederer-Schöberl-2020. My code is posted below. I am now trying to have an inflow boundary and a outflow/do-nothing boundary. The solution is a constant parabolic velocity-profile in the x-direction, and a linear pressure in the x-direction [attachment=undefined]u_y.png[/attachment] [attachment=undefined]u_y.png[/attachment] [attachment=undefined]u_y.png[/attachment] [attachment=undefined]u_x.png[/attachment]. However, when I solve the problem something strange happens. At the do-nothing boundary I get a sinusoidal profile in the y-direction, and the x-direction is also slightly perturbed. I have attached plots. I suspect this has to do with how I set the boundary-conditions on the do-nothing boundary, but I believe I have followed LLS-2020 virtually to the point.

I am very thankful for any possible help/ideas.

Best regards,



First the definition of the solver:

import re
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 netgen.csg import *
from netgen.meshing import MeshingStep

from numpy import pi

def SurfaceStokes_HdivHDG(mesh, force, g, order_u, order_p, order_g, bcdict, exactsoldict = None, mass=False, refinements = 0, pLagrange = False):

Mesh is a predefined surface mesh with boundary conditions set
force is the right-hand side f in the PDE
g is the divergence of the velocity field

order_u, order_p and order_g are the orders for the velocity u, the pressure p and the geometric mapping.
bcdict is a dicitonary of the form

    bcdict = {'dir': [], 'donothing': []}
where the lists should have the entries of the correct DOFS and the corresponding u-value. 
Any of the entries can be replaced by 'None'.
exactsoldict is also a dictionary of the form

    exactsoldict = {'velocity': uex, 'pressure': pex}
where pex can be set to 'None' if not needed/known. Any of the entries can be replaced by 'None'.

mass is a boolean which determines whether the Stokes problem showuld include a mass matrix, 
necessary for uniqueness if the surface is closed.
refinements is the number of times the mesh should be refined before calculating, necessary for convergence tests.
# Extract eventual solutions
if exactsoldict != None:
    uex = exactsoldict['velocity']
    pex = exactsoldict['pressure']
    uex = None
    pex = None

# Refine the mesh and give it the correct geometric order
for i in range(refinements):

# Defines the mixed Finite element space. Here Wh is the velocity space, What is the facet space defined for HDG,
# Qh is the pressure space and Nh is the real line, needed for Lagrange multiplier to p, for uniqueness. They are 
# put together in the mixed space Vh.
Wh = HDivSurface(mesh, order = order_u)
What = HCurl(mesh, order = order_u, flags={"orderface": 0})
Qh = SurfaceL2(mesh, order=order_p)
if pLagrange:
    Nh = NumberSpace(mesh)
    Vh = FESpace([Wh, What, Qh, Nh])
    Vh = FESpace([Wh, What, Qh])

if bcdict != None:
    dirichlet = ""
    for key in bcdict:
        if key != "outflow":
            dirichlet += "|"+key
    dirichlet = dirichlet[1:]
    mark_boundary_dofs([Vh.FreeDofs(False), Vh.FreeDofs(True)], "inflow|wall", Vh, clear=True)

gamma = 4*order_u*(order_u+1)    

print("Step 1 done")

# Defines the trial- and testfunctions
if pLagrange:
    u, uhat, p, r = Vh.TrialFunction()
    v, vhat, q, s = Vh.TestFunction()
    u, uhat, p = Vh.TrialFunction()
    v, vhat, q = Vh.TestFunction()
uhat = uhat.Trace()
vhat = vhat.Trace()

# Sets up the normal n, the face-tangential vector tau and the tangential co-normal mu
n = specialcf.normal(3)
tau = specialcf.tangential(3)
mu = Cross(n,tau)
# Defines the tangential projection, needed for HDG as the Hdiv space is normal continuous.
tang = lambda v: v-(v*mu)*mu

# Defines the projection matrix proj, the mesh size h and the stability parameter gamma.
nmat = CoefficientFunction(n, dims = (3,1))
proj = Id(3) - nmat*nmat.trans
h = specialcf.mesh_size

# Defines the symmetric gradients epsu and epsv, AKA the displacement matrices
gradu = proj * CoefficientFunction((grad(u),), dims = (3,3)).trans * proj # full gradient
gradv = proj * CoefficientFunction((grad(v),), dims = (3,3)).trans * proj # full gradient

epsu = 0.5 * (gradu + gradu.trans)
epsv = 0.5 * (gradv + gradv.trans)

print("Step 1.1 done")

# Defines the Stokes bilinear form

a = BilinearForm(Vh)
a += SymbolicBFI(2*InnerProduct(epsu,epsv), BND)
a += SymbolicBFI(-2*InnerProduct(epsu*mu, tang(v.Trace()-vhat)), BND, element_boundary = True)
a += SymbolicBFI(-2*InnerProduct(epsv*mu, tang(u.Trace()-uhat)), BND, element_boundary = True) #Symmetrizing
a += SymbolicBFI(gamma/h * InnerProduct(tang(v.Trace()-vhat), tang(u.Trace()-uhat)), BND, element_boundary = True) #Stabilizing

print("Step 1.2 done")

a += SymbolicBFI(- div(u.Trace())*q.Trace(), BND) #Enforces divergence condition
a += SymbolicBFI(- div(v.Trace())*p.Trace(), BND)

print("Step 1.3 done")
if pLagrange:
    a += SymbolicBFI(s.Trace()*p.Trace(), BND) #Uniqueness of pressure through Lagrange multiplier
    a += SymbolicBFI(r.Trace()*q.Trace(), BND) #Symmetrizing

print("Step 1.4 done")

# Adds mass form if necessary
if mass:
    a += SymbolicBFI(u.Trace()*v.Trace(), BND)

print("Step 1.5 done")
print("Step 2 done")
#Defines linear form for right-hand side

f = LinearForm(Vh)
f += SymbolicLFI( proj * force * v.Trace(), BND, bonus_intorder = order_u)    
f += SymbolicLFI(-g*q.Trace(), BND) #Enforces divergence condition
if pex != None and pLagrange: 
    f += SymbolicLFI (s.Trace()*pex, BND) #Ensures the pressure has the correct mean value

print("Step 3 done")

# Homogenize the problem

gfu = GridFunction(Vh)

res = gfu.vec.CreateVector() = f.vec

if bcdict != None:
    dof_in = BitArray(Vh.components[0].FreeDofs())
    mark_boundary_dofs(dof_in, "inflow", Vh.components[0], clear = False)
    u_in = bcdict['inflow']

    mass = BilinearForm(Vh.components[0])
    mass += SymbolicBFI((Vh.components[0].TestFunction().Trace()*mu) * (Vh.components[0].TrialFunction().Trace()*mu), BND, element_boundary = True)

    mass_rhs = LinearForm(Vh.components[0])
    mass_rhs += SymbolicLFI((Vh.components[0].TestFunction().Trace()*mu) *  (u_in * mu), BND, element_boundary = True)
    gfu.components[0] = mass.mat.Inverse(dof_in)*mass_rhs.vec
    #Draw(gfu.components[0], mesh, 'u_in')
    res -= a.mat*gfu.vec

# Inverts the matrix and solves the problem
ainv = a.mat.Inverse(Vh.FreeDofs(), inverse = "umfpack") += ainv * res

uh = gfu.components[0]
ph = gfu.components[2]

# Calculates eventual L2-errors of the solution

if uex != None:
    l2u = sqrt(Integrate((uex-uh)**2, mesh, BND))
    l2u = None
if pex != None:
    l2p = sqrt(Integrate((pex-ph)**2, mesh, BND))
    l2p = None
return {'velocity': uh, 'pressure': ph, 'l2u': l2u, 'l2p': l2p}

def mark_boundary_dofs(bitarrays, boundary, Vh, component=None, clear = False):
# bitarray should be a bitarray for the DOFS of the space Vh
if type(bitarrays) != list:
bitarrays = [bitarrays]

# Define what part of the space you want to mark the dofs for (if applicable), and
# start counting dofs at the correct spot

offset = 0
if component!= None:
    X = Vh.components[component]
    if components>0:
        for i in range(component):
            offset += Vh.components[i].ndof
    X = Vh

# Searches through all the relevant dofs, check if they match the correct type of boundary, and either
# clears or sets the dof.
for el in X.Elements(BBND):
    bnd_name = mesh.GetBBoundaries()[el.index]
    match =, bnd_name)
    if match != None:
        for dofs in el.dofs:
            for ba in bitarrays:
                if clear:

def SetFreeDofs(V):
for el in V.components[0].Elements(BBND):
for dofs in el.dofs:
for el in V.components[1].Elements(BBND):
for dofs in el.dofs:

And now the specific example

import netgen.gui

geo = CSGeometry()
bottom = Plane (Pnt(0,0,0), Vec(0,0, 1) )

surface = SplineSurface(bottom)
pts = [(0,0,0),(0,1,0),(1,1,0),(1,0,0)]
geopts = [surface.AddPoint(*p) for p in pts]
for p1,p2,bc in [(0,1,“inflow”), (1, 2,“wall”),(2,3,“outflow”),(3,0,“wall”)]:

ngmesh = geo.GenerateMesh(maxh=1)
mesh = Mesh(ngmesh)
for i in range(3):

uex = CoefficientFunction(( y*(1 - y), 0, 0 )).Compile(True)
g = CoefficientFunction( 0 ).Compile(True)
pex = CoefficientFunction( 1 - x ).Compile(True)
force = CoefficientFunction(( 1.00000000000000, 0, 0 )).Compile(True)

uin = uex
bcdict = {‘inflow’: uin, ‘wall’: None, ‘outflow’: None}

exactsoldict = {‘velocity’: uex, ‘pressure’: pex}

solutiondict = SurfaceStokes_HdivHDG(mesh, force, g, order_u=2, order_p=1, order_g=2, bcdict=bcdict, exactsoldict = exactsoldict, refinements = 0, mass=False, pLagrange = False)

import netgen.gui
uh = solutiondict[‘velocity’]
ph = solutiondict[‘pressure’]

Draw(uex, mesh, ‘uex’)
Draw(uh, mesh, ‘u’)
Draw(Norm(uh), mesh, ‘|u|’)
Draw(Norm(uh - uex), mesh, ‘|uh - uex|’)
Draw(Norm(uex), mesh, ‘|uex|’)
Draw(pex, mesh, ‘pe [attachment=undefined]p.png[/attachment]x’)
Draw(ph, mesh, ‘p’)
Draw(g - div(uh), mesh, ‘g’)