This page was generated from unit-9.2-C++Assemble/cppassembling.ipynb.

9.2 Implement our own system assembling

In this tutorial we

• write an integrators for \(\int_T f v dx\) and \(\int_T \nabla u \nabla v dx\): myIntegrator.hpp myIntegrator.cpp

• put together element matrices to the global system matrix: myAssembling.cpp

from ngsolve import *
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)

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.2-C++Assemble/myassemblemodule.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
      |   ^~~~~~~~~~~~~~~~

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

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

use our own integrators for element matrix calculation:

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

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

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

use our own matrix assembling function:

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

with TaskManager(pajetrace=10**8):
    # using our integrator
    mymatrix = m.MyAssembleMatrix(fes, m.MyLaplace(CF(1)), parallel=True)

    # using NGSolve built-in symbolic integrator
    # integrator = BilinearForm(grad(u)*grad(v)*dx).integrators[0]
    # mymatrix = m.MyAssembleMatrix(fes, integrator, parallel=True)

# print ("my matrix = ", mymat)

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

We assemble matrix
parallel assembling

from ngsolve.ngstd import Timers
for t in Timers():
    if t["name"].startswith("MyAssemble"):
        print (t["name"],"   ", t["time"])

MyAssemble - calc element matrix     0.06846664528213037
MyAssemble - assemble matrix     0.025100226550213087
MyAssemble - buildmatrixgraph     0.01897116181189274