9.3 High Order Finite Elements

• Finite elements implement the basis functions: myHOElement.hpp myHOElement.cpp

• Finite element spaces implement the enumeration of degrees of freedom, and creation of elements: myHOFESpace.hpp myHOFESpace.cpp

from ngsolve.fem import CompilePythonModule
from pathlib import Path

m = CompilePythonModule(Path('mymodule.cpp'), init_function_name='mymodule')

from netgen.occ import unit_square
from ngsolve import *
from ngsolve.webgui import Draw

mesh = Mesh(unit_square.GenerateMesh(maxh=0.2, quad_dominated=False))

We can now create an instance of our own finite element space:

fes = m.MyHighOrderFESpace(mesh, order=4, dirichlet="left|bottom|top")

and use it within NGSolve such as the builtin finite element spaces:

print ("ndof = ", fes.ndof)

gfu = GridFunction(fes)
gfu.Set(x*x*y*y)

Draw (gfu)
Draw (grad(gfu)[0], mesh);

and solve the standard problem:

u,v = fes.TnT()
a = BilinearForm(grad(u)*grad(v)*dx).Assemble()
f = LinearForm(10*v*dx).Assemble()
gfu.vec.data = a.mat.Inverse(fes.FreeDofs())*f.vec
Draw (gfu, order=3);

errlist = []
for p in range(1,13):
    fes = m.MyHighOrderFESpace(mesh, order=p)
    func = sin(pi*x)*sin(pi*y)
    gfu = GridFunction(fes)
    gfu.Set(func)
    err = sqrt(Integrate( (func-gfu)**2, mesh, order=5+2*p))
    errlist.append((p,err))
print (errlist)

import matplotlib.pyplot as plt
n,err = zip(*errlist)
plt.yscale('log')
plt.plot(n,err);

Exercises:

Extend MyHighOrderFESpace by high order quadrilateral elements.

http://www.numa.uni-linz.ac.at/Teaching/PhD/Finished/zaglmayr-diss.pdf, page 68 ff