H(div)-HDG on surfaces

Hi,

I am trying to implement H(div)-HDG on surfaces for Stokes equation, following Lehrenfeld-Schöberl-Lederer, 2020. My code is down below. For some reason it never finishes inverting the bilinear form a, or at least not for the approx 10 mins I have let it load. I have implemented it for Vector Laplace and it works fine and quickly, so I don't really get what is the problem. I have tried to remove the stuff from the example code which deals with condensation to try to remove possible error factors. Anyone have an idea of why it doesn't work, or has anyone had a similar problem?

Best regards

Alvar


Code:

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


geo = CSGeometry()
r = 1
sphere = Sphere(Pnt(0,0,0),r).bc("dir")
geo.Add(sphere)


order_u = 2
order_p = 1

mesh = Mesh(geo.GenerateMesh(maxh=0.2))
mesh.Curve(2)
Draw(mesh)


refinements = 1
for i in range(refinements):
mesh.Refine()

print("Done with part 1")

Draw(mesh)

print("Done with part 2")

n = specialcf.normal(3)
tau = specialcf.tangential(3)
mu = Cross(n,tau)


print("Done with part 3")

Wh = HDivSurface(mesh, order = order_u)
Wh = Compress(Wh)
What = HCurl(mesh, order = order_u, orderface = 0)
What = Compress(What)
Qh = SurfaceL2(mesh, order=order_p)
Nh = NumberSpace(mesh)
Vh = FESpace([Wh, What, Qh, Nh])





print("Done with part 4")




nmat = CoefficientFunction(n, dims = (3,1))
proj = Id(3) - nmat*nmat.trans

h = specialcf.mesh_size

gamma = 4*order_u*(order_u+1)



force = CoefficientFunction(((-56 * z ** 2 + 36 * z) * x ** 8 / 2 + 18 * x ** 7 * y ** 2 + (-2 - 168 * z ** 4 + 108 * z ** 3 + (-168 * y ** 2 + 228) * z ** 2 + (108 * y ** 2 - 140) * z) * x ** 6 / 2 + (108 * y ** 4 + 108 * y ** 2 * z ** 2 - 140 * y ** 2 + 2) * x ** 5 / 2 + (-168 * z ** 6 + 108 * z ** 5 + (-336 * y ** 2 + 456) * z ** 4 + (216 * y ** 2 - 280) * z ** 3 + (-168 * y ** 4 + 456 * y ** 2 - 342) * z ** 2 + (108 * y ** 4 - 280 * y ** 2 + 198) * z - 4 * y ** 2 + 6) * x ** 4 / 2 + (108 * y ** 2 * z ** 4 + (216 * y ** 4 - 280 * y ** 2 + 4) * z ** 2 + 108 * y ** 6 - 280 * y ** 4 + 198 * y ** 2 - 6) * x ** 3 / 2 + (-56 * z ** 8 + 36 * z ** 7 + (-168 * y ** 2 + 228) * z ** 6 + (108 * y ** 2 - 140) * z ** 5 + (-168 * y ** 4 + 456 * y ** 2 - 345) * z ** 4 + (108 * y ** 4 - 280 * y ** 2 + 202) * z ** 3 + (-56 * y ** 6 + 228 * y ** 4 - 347 * y ** 2 + 208) * z ** 2 + (36 * y ** 6 - 140 * y ** 4 + 202 * y ** 2 - 113) * z - 2 * y ** 4 + 6 * y ** 2 - 6) * x ** 2 / 2 + (36 * y ** 2 * z ** 6 + (108 * y ** 4 - 140 * y ** 2 + 2) * z ** 4 + (108 * y ** 6 - 280 * y ** 4 + 198 * y ** 2 - 6) * z ** 2 + 36 * y ** 8 - 140 * y ** 6 + 196 * y ** 4 - 108 * y ** 2 + 4) * x / 2 + 1 - 0.5e1 / 0.2e1 * z ** 6 + 2 * z ** 5 + (-10 * y ** 2 + 17) * z ** 4 / 2 + (8 * y ** 2 - 11) * z ** 3 / 2 + (-5 * y ** 4 + 17 * y ** 2 - 24) * z ** 2 / 2 + (4 * y ** 4 - 11 * y ** 2 + 12) * z / 2,-28 * y * ((z ** 2 - 0.9e1 / 0.14e2 * z) * x ** 7 - 0.9e1 / 0.14e2 * x ** 6 * y ** 2 + (0.1e1 / 0.28e2 + 3 * z ** 4 - 0.27e2 / 0.14e2 * z ** 3 + (-0.57e2 / 0.14e2 + 3 * y ** 2) * z ** 2 + (0.5e1 / 0.2e1 - 0.27e2 / 0.14e2 * y ** 2) * z) * x ** 5 + (-0.27e2 / 0.14e2 * y ** 4 - 0.5e1 / 0.28e2 + 0.5e1 / 0.2e1 * y ** 2 - 0.27e2 / 0.14e2 * y ** 2 * z ** 2) * x ** 4 + (3 * z ** 6 - 0.27e2 / 0.14e2 * z ** 5 + (-0.57e2 / 0.7e1 + 6 * y ** 2) * z ** 4 + (5 - 0.27e2 / 0.7e1 * y ** 2) * z ** 3 + (0.337e3 / 0.56e2 - 0.57e2 / 0.7e1 * y ** 2 + 3 * y ** 4) * z ** 2 + (-0.97e2 / 0.28e2 + 5 * y ** 2 - 0.27e2 / 0.14e2 * y ** 4) * z + y ** 2 / 14 - 0.3e1 / 0.28e2) * x ** 3 + (-0.27e2 / 0.14e2 * y ** 2 * z ** 4 + (5 * y ** 2 - 0.5e1 / 0.14e2 - 0.27e2 / 0.7e1 * y ** 4) * z ** 2 - 0.27e2 / 0.14e2 * y ** 6 - 0.107e3 / 0.28e2 * y ** 2 + 0.1e1 / 0.2e1 + 5 * y ** 4) * x ** 2 + (z ** 8 - 0.9e1 / 0.14e2 * z ** 7 + (-0.57e2 / 0.14e2 + 3 * y ** 2) * z ** 6 + (0.5e1 / 0.2e1 - 0.27e2 / 0.14e2 * y ** 2) * z ** 5 + (0.335e3 / 0.56e2 - 0.57e2 / 0.7e1 * y ** 2 + 3 * y ** 4) * z ** 4 + (-0.97e2 / 0.28e2 + 5 * y ** 2 - 0.27e2 / 0.14e2 * y ** 4) * z ** 3 + (-0.191e3 / 0.56e2 - 0.57e2 / 0.14e2 * y ** 4 + 0.337e3 / 0.56e2 * y ** 2 + y ** 6) * z ** 2 + (0.51e2 / 0.28e2 + 0.5e1 / 0.2e1 * y ** 4 - 0.97e2 / 0.28e2 * y ** 2 - 0.9e1 / 0.14e2 * y ** 6) * z - 0.3e1 / 0.28e2 * y ** 2 + 0.3e1 / 0.28e2 + y ** 4 / 28) * x - 0.9e1 / 0.14e2 * y ** 2 * z ** 6 + (-0.27e2 / 0.14e2 * y ** 4 - 0.5e1 / 0.28e2 + 0.5e1 / 0.2e1 * y ** 2) * z ** 4 + (-0.27e2 / 0.14e2 * y ** 6 - 0.107e3 / 0.28e2 * y ** 2 + 0.1e1 / 0.2e1 + 5 * y ** 4) * z ** 2 - 0.15e2 / 0.28e2 - 0.9e1 / 0.14e2 * y ** 8 + 0.65e2 / 0.28e2 * y ** 2 + 0.5e1 / 0.2e1 * y ** 6 - 0.51e2 / 0.14e2 * y ** 4),-28 * z ** 9 * x + 18 * z ** 8 * x + (-168 * x ** 3 - 168 * x * y ** 2 + 36 * y ** 2 + 228 * x) * z ** 7 / 2 + (108 * x ** 3 + 108 * x * y ** 2 - 140 * x) * z ** 6 / 2 + (-168 * x ** 5 + (-336 * y ** 2 + 456) * x ** 3 + 108 * x ** 2 * y ** 2 + (-168 * y ** 4 + 456 * y ** 2 - 345) * x + 108 * y ** 4 - 140 * y ** 2 + 2) * z ** 5 / 2 + (108 * x ** 5 + (216 * y ** 2 - 280) * x ** 3 + (108 * y ** 4 - 280 * y ** 2 + 198) * x) * z ** 4 / 2 + (-56 * x ** 7 + (-168 * y ** 2 + 228) * x ** 5 + 108 * x ** 4 * y ** 2 + (-168 * y ** 4 + 456 * y ** 2 - 357) * x ** 3 + (216 * y ** 4 - 280 * y ** 2 + 4) * x ** 2 + (-56 * y ** 6 + 228 * y ** 4 - 357 * y ** 2 + 219) * x + 108 * y ** 6 - 280 * y ** 4 + 198 * y ** 2 - 6) * z ** 3 / 2 + 18 * x * (x ** 6 + (3 * y ** 2 - 0.35e2 / 0.9e1) * x ** 4 + (3 * y ** 4 - 0.70e2 / 0.9e1 * y ** 2 + 0.101e3 / 0.18e2) * x ** 2 + y ** 6 - 0.35e2 / 0.9e1 * y ** 4 + 0.101e3 / 0.18e2 * y ** 2 - 0.113e3 / 0.36e2) * z ** 2 + (36 * x ** 6 * y ** 2 - 12 * x ** 5 + (108 * y ** 4 - 140 * y ** 2 + 2) * x ** 4 + (-24 * y ** 2 + 34) * x ** 3 + (108 * y ** 6 - 280 * y ** 4 + 198 * y ** 2 - 6) * x ** 2 + (-12 * y ** 4 + 34 * y ** 2 - 36) * x + 36 * y ** 8 - 140 * y ** 6 + 196 * y ** 4 - 108 * y ** 2 + 4) * z / 2 + 2 * x ** 5 + (8 * y ** 2 - 11) * x ** 3 / 2 + (4 * y ** 4 - 11 * y ** 2 + 14) * x / 2)).Compile(True)
uex = proj * CoefficientFunction((-z*z,y,x))
pex = None


u, uhat, p, r = Vh.TrialFunction()
v, vhat, q, s = Vh.TestFunction()

uhat = uhat.Trace()
vhat = vhat.Trace()

tang = lambda v: v-(v*mu)*mu


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("Done with part 5")


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

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

a += SymbolicBFI(s.Trace()*p.Trace(), BND)
a += SymbolicBFI(r.Trace()*q.Trace(), BND)


print("Done with part 6")

f = LinearForm(Vh)
f += SymbolicLFI( proj * force * v.Trace(), BND, bonus_intorder = order_u)
if pex != None:
f += SymbolicLFI (s.Trace()*pex, BND)


print("Done with part 7")

a.Assemble()

print("Done with part 8")
f.Assemble()
print("Done with part 9")

gfu = GridFunction(Vh)

def SetFreeDofs(V):
for el in V.components[0].Elements(BBND):
for dofs in el.dofs:
V.FreeDofs().Clear(dofs)
V.FreeDofs(True).Clear(dofs)
for el in V.components[1].Elements(BBND):
for dofs in el.dofs:
V.FreeDofs().Clear(dofs+V.components[0].ndof)
V.FreeDofs(True).Clear(dofs+V.components[0].ndof)
print("Done with part 9.1")

SetFreeDofs(Vh)

print("Done with part 9.2")
#ainv = a.mat.Inverse(Vh.FreeDofs())

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

print("Done with part 10")
gfu.vec.data = ainv * f.vec

vel = gfu.components[0]

print("Done with part 11")
Draw(uex, mesh, 'uex')
Draw(vel, mesh, 'u')
Draw(Norm(uex), mesh, '|uex|')
Draw(Norm(vel), mesh, '|u|')
Draw(Norm(vel-uex), mesh, '|u-uex|')
Draw(force, mesh, "force")


print("Done with part 12")
print(sqrt(Integrate((uex - vel)**2, mesh, BND)))
print(sqrt(Integrate(uex**2, mesh, BND)))
print(sqrt(Integrate(vel**2, mesh, BND)))

Solved it, it had to do with my meshing of the sphere

Yes, I just wanted to mention this. Easy fix is to move the "from ngsolve import *" below the "from netgen import..." and use

mesh = Mesh(geo.GenerateMesh(perfstepsend=MeshingStep.MESHSURFACE,optsteps2d=3,maxh=0.3))

Best, Philip