How can I use Q2-P1dc element pair in a quadrilateral mesh to solve Stokes?

Main Forum Old-Forum - do not post here
1

Hi, all

In NGSolve, it provides Pk element w.r.t. a triangulation, or Qk element w.r.t. a quadrilateral mesh. To solve Stokes equation, we try to use Q2-P1dc in a quadrilateral mesh. Is there a method in NGSolve to set P1dc element space for pressure in a quadrilateral mesh?

Best,
Di Yang

2

Hi Di Yang,

You can hack the Q1 L2 space by removing the bilinear parts. Below is a working example.

Best,
Guosheng

Attachment: Q2P1.py

3

Thank u very much. By the way, does NGSolve provide high order BDFM element directly? If not, how can I design it just in NGSolve instead of using other platforms?

The BDFM element is a reduced element from BDM element. Because I have no idea how BDM is designed in the core of NGSolve developing, how to obtain BDFM elements in quadrilateral meshes is still an open issue for me.

4

I think you can hack a vectorQk+1 space to get a DG version of the BDFMk space on quads (no hdiv-conformity), just removing unwanted high order DOFs.

You can take a look at this thesis for the construction of basis function in ngsolve
http://www.numa.uni-linz.ac.at/Teaching/PhD/Finished/zaglmayr

All the finite elements are in the source file

5

How is the Qk element L2basis ranked? Is there a law about its ranking?

6

Thank u for my coding. But I am still confused how to hack a vector Q{k+1} space to get a DG version of the BDFMk space on quads. The scalar Qk+1 space is easy to implement but for a vector space, what is the structure of its basis?

(Q1)^2=span{(1,0), (x,0), (y,0), (xy,0), (0,1), (0,x), (0,y), (0,xy)},
please tell me how to hack it in the codes to get DG version of BDFM1 space. Note that dim(BDFM1)=4 by the way.

7

First, you shall figure out the basis functions for Qk by yourself. They are classical tensor products of 1D Legendre basis. VectorL2 is simply two copies of Qk on quads.

You can do so by visualize the shape functions of the basis on a unit square without going into details of the c++ code, e.g.,

mesh = MakeStructured2DMesh(nx=1,ny=1)
gfu = GridFunction(VectorL2(mesh, order=2))

Draw(gfu)
for i in range(gfu.vec.size):
    gfu.vec[i] = 1
    if i > 1:
         gfu.vec[i-1] = 0
     Redraw()
     input("next?")
8

If I just simply run the following code:

import netgen.gui
from ngsolve import *
import ngsolve.meshes as ngm
mesh = ngm.MakeStructured2DMesh(nx=1,ny=1)
gfu = GridFunction(VectorL2(mesh, order=2))
Draw(gfu)
gfu.vec[0] = 1
Redraw()

there is nothing shown for me in GUI, which confuses me. Can you tell me some reason?

9

Thank u again. Using the VtkOutput and viewing in ParaView can realize your idea. Nice suggestion.

10

If you run

python3 test.py

and want to see the pop up gui, put the following line after Redraw().

input("next?")

I think i-tutorials is the good place for you to get familiar with ngsolve.