Hi,

I am working on the stationary Stokes problem. My advisor told me that the stability of Taylor-hood elements for the Stokes problem when using triangular elements requires that each element only has maximum one side on the boundary. Hence, we thought to use a large mesh that satisfies this, and then refine it. However, the error of the solution diverges with more refinements. Any idea’s? The code is found below.

Best regards,

Alvar

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), (1, 1), bcs = (“top”, “left”, “bottom”, “right”))

#uex = CoefficientFunction((2*pi*sin(pi*x)**2*cos(pi*y) sin(piy), -2*pi

*sin(pi*x)

*cos(pi*x)

*sin(pi*y)

**2))**

#force = CoefficientFunction((12*pi3

#force = CoefficientFunction((12*pi

*sin(pi*x)

**2**3

*sin(pi*y)*cos(pi*y) - pi*sin(pi*x)*cos(pi*y) - 4*pi*sin(pi*y)

*cos(pi*x)

**2**3

*cos(pi*y),-12*pi*sin(pi*x)

*sin(pi*y)

**2**3

*cos(pi*x) + 4*pi*sin(pi*x)

*cos(pi*x)

*cos(pi*y)**2 - pi

*sin(pi*y)

*cos(pi*x)))

#pex = CoefficientFunction(cos(pi*x)

*cos(pi*y))

nu = 1

#uex = CoefficientFunction((sin(pi*y), -cos(pi*x)))

#force = CoefficientFunction((-pi*sin(pi*x)*cos(pi*y) + nu * pi**2 sin(piy), -pisin(piy)cos(pix) - nu * pi**2

*cos(pi*x)))

#pex = CoefficientFunction(cos(pi*x)

*cos(pi*y))

#uex = CoefficientFunction((sin(y) + cos(y), 0))

#force = CoefficientFunction((sin(y) + cos(y) + cos(x),0))

#pex = sin(x)

uex = CoefficientFunction((y*(1-y), 0))

force = CoefficientFunction((3,0))

pex = CoefficientFunction(x-1/2)

l22_list = []

l2p_list = []

def my_abs(mesh, v, h0):

a = 0

```
x = np.arange(0,1+h0, h0)
for i in range(len(x)):
for j in range(len(x)):
p1 = np.array([x[i]])
p2 = np.array([x[j]])
mip = mesh(p1, p2)
val = abs(v(mip))
a = max(a, val)
return a
```

def l2error(refinement = 0, structured = False, quad = False):

```
order = 2
alpha = 10
h0 = 0.3
if structured:
mesh = Mesh(geo.GenerateMesh(maxh = h0, quad_dominated = quad))
for i in range(refinement):
mesh.Refine()
#Draw(mesh)
else:
mesh = Mesh(geo.GenerateMesh(maxh = h0*2**(-refinement), quad_dominated = quad))
#Wh = HDiv(mesh, order = order+1, dgjumps = True, BDM = True) #Velocity space
Wh = VectorH1(mesh, order = order, dgjumps = True) #Velocity space
Qh = L2(mesh, order = order-1) #Pressure space
Nh = NumberSpace(mesh) #Lagrange multiplier
Vh = FESpace([Wh,Qh,Nh], dgjumps = True)
u,p,r = Vh.TrialFunction()
v,q,s = Vh.TestFunction()
n = specialcf.normal(2)
h = specialcf.mesh_size
tang = lambda v:v #If no tangential projections are wanted
#tang = lambda v: v-(v*n)*n #If tangential projections are wanted
jump_u = tang(u-u.Other())
jump_v = tang(v-v.Other())
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)
a += SymbolicBFI(p*s)
a += SymbolicBFI(q*r)
#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,tang(v)), BND, skeleton = True)
a += SymbolicBFI(-nu*InnerProduct(Grad(v)*n,tang(u)), BND, skeleton = True) #symmetrizing term
a += SymbolicBFI(h**(-1) * alpha * order**2 * InnerProduct(tang(u),tang(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
f += SymbolicLFI(pex*s)
a.Assemble()
f.Assemble()
gfu = GridFunction(Vh)
gfuex = GridFunction(Wh)
gfuex.Set(uex)
gfpex = GridFunction(Qh)
gfpex.Set(pex)
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)
divu = div(vel)
#eu = sum((dv-dg)**2)
eu = (vel - gfuex)**2
ep = (pres - gfpex)**2
err = vel-gfuex
Draw(Norm(err), mesh, "|e|")
l22 = Integrate(eu, mesh)
l2p = Integrate(ep, mesh)
l22_list.append(l22)
l2p_list.append(l2p)
print("\nRefinement level ", refinement)
print("L2 error in velocity: ", sqrt(l22))
print("L2 error in pressure: ", sqrt(l2p))
print("div(u)-norm", Integrate(divu**2, mesh))
if len(l22_list)>1:
print("EOCu: ", (np.log(sqrt(l22_list[-2]))-(np.log(sqrt(l22_list[-1]))))/np.log(2))
print("EOCp: ", (np.log(sqrt(l2p_list[-2]))-(np.log(sqrt(l2p_list[-1]))))/np.log(2))
```

for i in range(3,4):

l2error(i, structured = True, quad = False)