Sparsity pattern of Hdiv mass matrix

Hi all,

at the moment i try to familiarize myself with the implementation of the HDiv basis functions.
I took a look at the sparsity pattern of the mass matrix by

``````mesh = Mesh(unit_square.GenerateMesh(maxh=0.2))
fes = HDiv(mesh, order=25, dirichlet="top|bottom|right|left")
u = fes.TrialFunction()
v = fes.TestFunction()
SetTestoutFile ("test.out")
a = BilinearForm(fes,printelmat=True)
a += u*v*dx
#And so on..``````

and then importing it into matlab.

It looks for the most parts as expected, but there are some dense blocks at both ends of the matrix (far away from machine precision, i.e. 10^-6). Where do these blocks come from?

Best
Tim

Hi Tim,

you can compute and process element-matrices directly within Python, without going over the testout file:

``````from netgen.geom2d import unit_square
from ngsolve import *

mesh = Mesh(unit_square.GenerateMesh(maxh=2))

fes = HDiv(mesh, order=5)
u,v = fes.TnT()

a = BilinearForm(fes)
a += u*v*dx

igt = a.integrators[0]

ei = ElementId(0)
el = fes.GetFE(ei)
trafo = mesh.GetTrafo(ei)

mat = igt.CalcElementMatrix(el, trafo)
print (mat)

import scipy.linalg
ev,evec = scipy.linalg.eigh(a=mat) ``````

You have to look into the H(div) basis functions, in the file fem/hdivhofe_impl.hpp, starting line 89:

We are using something like Dubiner, multiplied by RT0 functions. This gives us exactly the RT space.

We could use a construction like the Dubiner, but with different Jacobi weights to improve non-zero entries. If it improves conditioning (after diagonal scaling) it is another good argument to change.

The last block should better be a scaled Legendre, maybe this is what you have observed.

Joachim

Hi Joachim,

thanks for the tipp.

I forgot the attachment in the previous post. It’s not just one block, but some dense columns and rows. Is that the scaled Legendre? Looks kinda weird.

(There are no such blocks for example in Beuchler, Pillwein Zaglmayer (2010))

Tim

Attachment: SparsityPattern.pdf