# 9.2 Implement our own system assembling¶

In this tutorial we

[1]:

from ngsolve import *
import pyngcore   # for timers

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

from ngsolve.fem import CompilePythonModule
from pathlib import Path

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

[1]:

['MyAssembleMatrix',
'MyLaplace',
'MySource',
'__doc__',
'__name__',
'__package__',
'__spec__']

[2]:

mesh = Mesh(unit_square.GenerateMesh(maxh=0.2))
fes = H1(mesh, order=1, dirichlet=".*")
u, v = fes.TnT()


## use our own integrators for element matrix calculation:¶

[3]:

a = BilinearForm(fes)
a += m.MyLaplace(1)
a.Assemble()

f = LinearForm(fes)
f += m.MySource(x)
f.Assemble();

[4]:

gfu = GridFunction(fes)
gfu.vec.data = a.mat.Inverse(fes.FreeDofs()) * f.vec
Draw (gfu);


## we can inspect the element matrix:¶

[5]:

mylap = m.MyLaplace(1)
ei = ElementId(VOL,17)
mylap.CalcElementMatrix(fes.GetFE(ei), mesh.GetTrafo(ei))

[5]:

 0.665986 -0.426433 -0.239553
-0.426433 0.648429 -0.221997
-0.239553 -0.221997 0.46155


## use our own matrix assembling function:¶

[6]:

for l in range(4): mesh.Refine()
fes.Update()
gfu.Update()
f.Assemble();

[7]:

pyngcore.ResetTimers()
print ("num elements =", mesh.ne, ", ndof =", fes.ndof)
# using our integrator
mymatrix = m.MyAssembleMatrix(fes, m.MyLaplace(CF(1)), parallel=False)

# using NGSolve built-in symbolic integrator

# print ("my matrix = ", mymat)

if fes.ndof < 100000:
gfu.vec.data = mymatrix.Inverse(fes.FreeDofs()) * f.vec
Draw (gfu);

num elements = 14336 , ndof = 7329
We assemble matrix
sequential assembling


## Compare timings for¶

• MyLaplace and grad(u)*grad(v)*dx

• sequential and parallel

• for different polynomial orders

[8]:

for t in pyngcore.Timers():
if t["name"].startswith("MyAssemble"):
print (t["name"],"   ", t["time"])

MyAssemble - calc element matrix     0.009151094954780491
MyAssemble - assemble matrix     0.01567438885565707
MyAssemble - buildmatrixgraph     0.0041391719826308114


## Exercise:¶

• Implement MyAssembleVector for building the right hand side vector for $$\int_\Omega f v \, dx$$

• Implement MyNeumannIntegrator for evaluating $$\int_{\partial \Omega} g v \, ds$$

[ ]: