Implement our own system assembling

23. Implement our own system assembling#

In this tutorial we

  • write integrators for \(\int_T f v dx\) and \(\int_T \nabla u \nabla v dx\):

  • put together element matrices to the global system matrix:

The implementation is available on github: TUWien-ASC/NGS-myassembling

from ngsolve import *
from ngsolve.webgui import Draw
import myassembling
dir (myassembling)

---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
Cell In[1], line 3
      1 from ngsolve import *
      2 from ngsolve.webgui import Draw
----> 3 import myassembling
      4 dir (myassembling)

ImportError: dlopen(/Users/joachim/Library/Python/3.14/lib/python/site-packages/myassembling.so, 0x0002): Symbol not found: __ZN4ngla12SparseMatrixIdddEC1EiiRKN6ngcore5TableIimEES6_b
  Referenced from: <347BE42F-B0A0-38A4-BFCA-A4ABA2406FDB> /Users/joachim/Library/Python/3.14/lib/python/site-packages/myassembling.so
  Expected in:     <1A2542D6-D60B-3898-A3E8-FB9EA99A299B> /Applications/Netgen.app/Contents/MacOS/libngla.dylib

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

23.1. use our own integrators for element matrix calculation:#

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

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

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[3], line 2
      1 a = BilinearForm(fes)
----> 2 a += myassembling.MyLaplace(1)
      3 a.Assemble()
      4 
      5 f = LinearForm(fes)

NameError: name 'myassembling' is not defined

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

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[4], line 2
      1 gfu = GridFunction(fes)
----> 2 gfu.vec.data = a.mat.Inverse(fes.FreeDofs()) * f.vec
      3 Draw (gfu);

TypeError: matrix not ready - assemble bilinearform first

23.2. we can inspect the element matrix:#

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

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[5], line 1
----> 1 mylap = myassembling.MyLaplace(1)
      2 ei = ElementId(VOL,17)
      3 mylap.CalcElementMatrix(fes.GetFE(ei), mesh.GetTrafo(ei))

NameError: name 'myassembling' is not defined

23.3. use our own matrix assembling function:#

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

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[6], line 4
      1 for l in range(4): mesh.Refine()
      2 fes.Update()
      3 gfu.Update()
----> 4 f.Assemble();

NameError: name 'f' is not defined

mymatrix = myassembling.MyAssembleMatrix(fes, myassembling.MyLaplace(CF(1)))
# mymatrix = myassembling.MyAssembleMatrix(fes, (grad(u)*grad(v)*dx)[0].MakeBFI())

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

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[7], line 1
----> 1 mymatrix = myassembling.MyAssembleMatrix(fes, myassembling.MyLaplace(CF(1)))
      2 # mymatrix = myassembling.MyAssembleMatrix(fes, (grad(u)*grad(v)*dx)[0].MakeBFI())
      3 
      4 if fes.ndof < 100000:

NameError: name 'myassembling' is not defined

23.4. 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\)

  • Have a look into thread-parallel assembling at NGSolve-itutorials unit 9.2

Contents