Tensor Product Space and HCurl/HDiv

Main Forum Old-Forum - do not post here
1

I recently started to work with NGSolve.
I would like to work with the tensor product space. However it always crashes when it comes to HCurl and HDiv. Or at least if I try to evaluate a function of a tensor product space with HDiv/HCurl. One small example would be:

[code]from ngsolve.TensorProductTools import *
from ngsolve import *
from netgen.geom2d import unit_square
import netgen.gui
%gui tk

meshx = Mesh(unit_square.GenerateMesh(maxh=0.15))
mesht = Mesh(SegMesh(20,0,1,periodic=True) )
tpmesh = Mesh(MakeTensorProductMesh(mesht,meshx))
Draw(tpmesh)

n=1
m=1

fesx = HDiv(meshx,order=n)
fest = H1(mesht,order=m)
tpfes = TensorProductFESpace([fest,fesx])

u = GridFunction(tpfes)

gtest = x
g = ProlongateCoefficientFunction(gtest,0,tpfes)

u.Set(g)[/code]

It crashes with the “u.Set()”. When I looked into the code it seemed to me, that the “.Set” uses the “.GetFE”-routine which isn’t implemented for the two spaces HDiv and HCurl.
My problem is that I haven’t found the right Set/evaluation routines for those two spaces. Could you give me some advise about the parts of the code I should look into to fix this problem?

2

Hi,

I did a quick debug of your problem:

The tensor-product Set assumes scalar finite elements (in tpintrule.hpp, line 155). It should be rewritten using the Evaluators of the component spaces.

Gerhard may help here.

Don’t worry, also the HCurl and HDiv - spaces have GetFE - functions.

Btw, your tensor-product space is vector valued, so it does not make sense to set a function to the scalar function x.

Joachim

3

Thank you very much for the quick reply!

So it would be best to distinguish between the finite element spaces and use for HCurl the BaseHCurlFiniteElement in line 155 and 156 of tpintrule.hpp? Would you use HDivFiniteElement in the HDiv Case? (However I didn’t find the implementation of “CalcShape” for HDivFiniteElement in hdivfe.cpp. Only for FE_RTTrig0, FE_BDMTrigl and so on.)

Thank you for the remark. I assumed falsely that it would be projected to a vector valued function. (I got a bit confused with the gradx/grady naming of the gradiant of the tensor product space.)

4

Hi jhauser,

It should be possible now, to construct an H1 times HDiv ( and H1 times HCurl) and to set gridfunctions in that space to a given coefficient function. I changed tthe code from your post to give you a first example.

Since one of your spaces is called fest, I expect that you have space-time applications in mind? Because currently we only support operator application for BilinearForms defined in TensorProductFESpace, no system matrix can yet be assembled. ( Since it was intended for DG operator application initially, so mainly a historical reason).

Additionally I’m not sure, if vectorial spaces are treated correctly when applying your BilinearForm.

However, I think most of the above mentioned things an be solved.

Gerhard

5

Sorry, I forgot to give you an example code also. Here it is…
https://ngsolve.org/media/kunena/attachments/780/SetExample.py

Attachment: SetExample.py

6

I updated my ngsolve version yesterday and today. However your example still crashes with the “.Set”. The two lines in tpintrule.hpp (155&156) haven’t changed as well.

Yes, space-time is the goal. But i knew that there will be work before things run smoothly. There isn’t a software yet which is open source for space-time with HCurl/HDiv. (At least as far as I know) Which is why I am trying to get things running with NGSolve.
Is there anything I should keep in mind before I tackle the problem of the BilinearForm?

7

Hi again, sorry for the delay, the changes were merged some minutes ago.
The two lines in tpintrule haven’t change at all, but the function “CalcShape” isn’t called anymore.
The changes are in tpfes.cpp in the functions TransferToTPMesh and TransferToStdMesh.
Should run now.

8

I just merged Gerhard’s fixes, it will be in the nightly release tonight

9

First: Thank you for the fixes!

After some trial and error work the assembling for special differential equations, seems to be working as well.

Now I am am trying to understand the behavior of the interpolation error a bit better. I compare CoefficientFunctions projected on the L2 space (3D) and ProlongateCoefficientFunctions(1D+2D) first projected into the tensor product space and then into the L2. For that I have two questions for understanding the evaluation of coefficient functions in the tp space:

  1. Does the line
CoefficientFunction->Evaluate(MapedIntegrationRule,Matrix);

especially

cfstd->Evaluate(tpmir,result);

give me the evaluation of the CoefficientFunction in the mapped interpolation points? Or rather some coefficients for the basis functions corresponding to the integration rule (and therefore some interpolation)?

  1. I am a bit confused about what’s going on in TransferToTPMesh. I interpret the lines 678 to 687 as something like
    elmatout = A^{-1} A_2^{-T}A^T resultasmat B B^{-1} B_2^{-T}
    which means
    elmatout = A^{-1} A_2^{-T}A^T resultasmat B_2^{-T}
    for
    A = shapesx, A_2 = shapesx2,
    B = shapesy, B_2 = spapesy2.
    Was there transposing intentionally left out or did I misinterpret the CalcInverse?
10

Hi there,

Concerning your questions:
1.) cfstd->Evaluate(tpmir,result);
Probably you refer to line 662 in tpfes.cpp.
Yes, the function evaluates cfstd at the integration points in tpmir.

2.) The lines you mentioned in your post are correct. Your manipulations will only work when the shape matrices quadratic, otherwise they’re not invertible. On the other hand, if they were, use that the mass matrix is symmetric, i.e. you could exchange line 687 also by

elmatout = Trans(elmatx)*elmatrhs*elmaty;

Then your calculations end up in A_2^{-T} resultasmat B_2^{-T}.