This page was generated from unit-9.3-highorder/highorder.ipynb.

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')

called 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)

ndof =  489

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)

[(1, 0.018470559015371108), (2, 0.0016028677378615212), (3, 0.00011616801306570434), (4, 5.803815227864814e-06), (5, 4.056875366717925e-07), (6, 1.0483434394610086e-08), (7, 7.453635532436866e-10), (8, 1.2046110241144944e-11), (9, 8.645687137533428e-13), (10, 2.9273665456945194e-13), (11, 5.71196758841869e-13), (12, 1.3652831039644462e-12)]

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