P1ISOP2-P1

Main Forum Old-Forum - do not post here
1

Dear All,

I’m trying to implement the P1ISOP2 element for the Stokes problem on a square.
To do this I need to work with two different mesh, in particular, I need the test function to be on coarser mesh than the trial function.
I generate both meshes using the MakeStructured2DMesh function with the Quad flag set to False.
Now when I try to build this the bilinear form B = div(u)qdx I get the error that both trial and test space must be defined on the same space.

Do you have any suggestions on how to proceed?

In my mind, a solution could have been to create a custom mixed finite element integrator, as shown in GitHub - NGSolve/mylittlengsolve but I it seems I’m not able to build mylittlengsolve.

Thank you very much for your help!
Kind Regards,
Umberto Zerbinati

Attachment: Example1.py

2

Yes you definitely need to go to the cpp level to do this. A very rough idea of how this could be done:
Create a wrapper finite element space that is defined on the fine grid and does the mapping from fine to coarse grid and evaluation of basis functions on the coarse grid.
Then you can create a mixed bilinearform with

a = BilinearForm(fesFine, testspace=wrapperFesCoarse)

and add the forms. Note that the spaces are both defined on the same mesh (this is expected by the integrators) and the coarse mapping is done in the test space.

What is your problem with building mylittlengsolve?

You do not need much to create a c++ extension, a short introduction is given here:
https://docu.ngsolve.org/v6.2.2007/mylittlengs/1_Basic/pythonExport.html

Best
Christopher

3

Dear Christopher,

thank you very much for the suggestion!
I attach a log with the error I obtain from the make command

Kind regards,
Umberto

https://ngsolve.org/media/kunena/attachments/1202/error.txt

Attachment: error.txt

4

I can only see one line:
[ 9%] Building CXX object CMakeFiles/myngspy.dir/2_myHOFEM/myHOFESpace.cpp.o

5

Sorry, my bad-
This time should be the full file.
Thanks,
Umberto

Attachment: error_2021-06-21-2.txt

6

Ah, yes this was due to changes in recent NGSolve. It is fixed now.
Best
Christopher

7

Works like a charm !
Thank you very much,

8

Dear Christopher,
I’m working on the C++ side, and I was wondering is there a way to get the coordinates of the vertex of the element inside the function CalcShape and CalcDShape.
Kind regards,
Umberto

9

Hi Umberto,

if you have a MappedIntegrationPoint, you can get the ElementTransformation from it, and from the ElementTransformation you can get the Mesh and the ElementId.

A short example is attached.

It uses the brand-new feature of including C++ code in Python scripts.

Joachim

Attachment: cppextension_2021-06-26-2.py

10

Hi Umberto,

we can now do it by refining the mesh, and keeping the pressure in the range of the prolongation operator.

I added the method to the Stokes-tutorial, available in the latest nightly docu:

https://docu.ngsolve.org/nightly/i-tutorials/unit-2.6-stokes/stokes.html#(P1-iso-P2)-P1

Joachim

11

Dear Joachim,

thank you very much for the example on the prolongation operator and for the example on the class MappedIntegrationPoint !!

Kind regards,
Umberto

12

Dear Joachim,
is there a way to implement the prolongation operator also for the custom elements that I’ve implemented using mylittlengsolve ?
Kind regards,
Umberto

13

The finite element space has a shared_ptr prol member, that has to be set in the space constructor to your custom prolongation. To see how such prolongations are implemented you can have a look at the ngsolve/multigrid/prolongation.*pp files.

Best Christopher

14

I manage to implement the prolongation operator on my P1+P0 elements thank you very much for all the help!