- Thank you received: 0
Forum Message
We have moved the forum to https://forum.ngsolve.org . This is an archived version of the topics until 05/05/23. All the topics were moved to the new forum and conversations can be continued there. This forum is just kept as legacy to not invalidate old links. If you want to continue a conversation just look for the topic in the new forum.
Notice
The forum is in read only mode.
Sparsity pattern of Hdiv mass matrix
3 years 9 months ago #2931
by THaubold
Sparsity pattern of Hdiv mass matrix was created by THaubold
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
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
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
Code:
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..
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
3 years 9 months ago #2932
by joachim
Replied by joachim on topic Sparsity pattern of Hdiv mass matrix
Hi Tim,
you can compute and process element-matrices directly within Python, without going over the testout file:
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
you can compute and process element-matrices directly within Python, without going over the testout file:
Code:
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
3 years 9 months ago - 3 years 9 months ago #2933
by THaubold
Replied by THaubold on topic Sparsity pattern of Hdiv mass matrix
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
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
Attachments:
Last edit: 3 years 9 months ago by THaubold.
3 years 9 months ago #2934
by THaubold
Replied by THaubold on topic Sparsity pattern of Hdiv mass matrix
Couldn't upload png or jpg, so i added a pdf file.
.
.
3 years 9 months ago #2937
by joachim
Replied by joachim on topic Sparsity pattern of Hdiv mass matrix
By replacing the last group by scaled Legendre, the last block got sparsified, and conditioning improved as well. Thank you for reporting !
For reason of having the RT (and not only RT-like as in the construction with Sabine) we changed the H(div) basis last year. One can certainly optimize also here the number of non-zero diagonals by adjusting the Jacobi indices. Pls let us know about your findings.
For reason of having the RT (and not only RT-like as in the construction with Sabine) we changed the H(div) basis last year. One can certainly optimize also here the number of non-zero diagonals by adjusting the Jacobi indices. Pls let us know about your findings.
Time to create page: 0.153 seconds
