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

[1]:

from ngsolve.fem import CompilePythonModule
from pathlib import Path

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

called mymodule

In file included from /usr/bin/../include/netgen/fem.hpp:75,
                 from /usr/bin/../include/netgen/comp.hpp:10,
                 from /root/src/ngsolve/docs/i-tutorials/unit-9.3-highorder/mymodule.cpp:1:
/usr/bin/../include/netgen/diffop.hpp: In function ‘std::shared_ptr<ngfem::DifferentialOperator> ngfem::CreateFunctionDiffOp(const DOP&, const F&, int)’:
/usr/bin/../include/netgen/diffop.hpp:1127:24: warning: ‘template<class DOP, class F> class ngfem::T_FunctionDiffOp’ is deprecated: guess it never got over the first experimental use [-Wdeprecated-declarations]
 1127 |     return make_shared<T_FunctionDiffOp<DOP, F>> (func, dim);
      |                        ^~~~~~~~~~~~~~~~
/usr/bin/../include/netgen/diffop.hpp:1062:3: note: declared here
 1062 |   T_FunctionDiffOp : public DifferentialOperator
      |   ^~~~~~~~~~~~~~~~

[2]:

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:

[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)
Draw (grad(gfu)[0], mesh);

and solve the standard problem:

[6]:

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

[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.018470561568516793), (2, 0.001602867608186379), (3, 0.00011616801501656098), (4, 5.799927993927836e-06), (5, 4.0563658382219457e-07), (6, 1.0475203105369778e-08), (7, 7.452332085368543e-10), (8, 1.2039759522275085e-11), (9, 8.858089559311819e-13), (10, 3.0963984468444877e-13), (11, 4.670388275962591e-13), (12, 1.3230051795246677e-12)]

[8]:

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

../../_images/i-tutorials_unit-9.3-highorder_highorder_12_0.png

Exercises:

Extend MyHighOrderFESpace by high order quadrilateral elements.

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

[ ]: