# 9.3 High Order Finite ElementsΒΆ

[1]:

from ngsolve.fem import CompilePythonModule
from pathlib import Path

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

called mymodule

[2]:

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



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

[3]:

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


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

[4]:

print ("ndof = ", fes.ndof)

ndof =  489

[5]:

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

Draw (gfu)


and solve the standard problem:

[6]:

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

[7]:

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.01847055901538545), (2, 0.0016028677378614217), (3, 0.00011616801306597612), (4, 5.80381522616505e-06), (5, 4.0568753743986986e-07), (6, 1.0483438418049003e-08), (7, 7.453675123934696e-10), (8, 1.2049863405103921e-11), (9, 8.874250841087325e-13), (10, 4.3306659616668924e-13), (11, 5.126904771937532e-13), (12, 1.231766820985857e-12)]

[8]:

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


Exercises:

Extend MyHighOrderFESpace by high order quadrilateral elements.

[ ]: