Mini element in 3D

NGSolve
1

I konw that the syntax for 2D mini element is:

velocity = VectorH1(mesh,order = 1)
velocity.SetOrder(TRIG,3)   # mini elements
velocity.Update()
pressure = H1(mesh,order = 1)

How about the right syntax for 3D mini element in Ngsolve? Are they same or not? many thanks!

2

yes, except you need

velocity.SetOrder(TET,4)

at least it should work like this, don’t know if someone has tried it out.

3

For 3D NS equations,

I tried

velocity.SetOrder(TET,4)

the codes will break down after several time steps.
I don’t know reasons. I guess this syntax add exact one 4th order polynomial in the interior of tetrahedral element.

However, I found

velocity.SetOrder(TRIG,3)

is stable, can I understand it as adding a bubble function for each face of terahedral elements? Since each face just is a 2D domain.

4

what do you mean with break down ?

Maybe you first do some tests with Stokes before moving to Navier-Stokes,

Joachim

5

Sorry, I made mistakes due to coarse mesh(radius = 0.3; maxh = 0.1) for Cylinder 3D tube domain.

BUT,
When radius = 0.3, maxh = 0.05;
I followed your advice to test Transient Stokes. I find the results for BOTH

velocity.SetOrder(TRIG,3)

&

velocity.SetOrder(TET,4)

are same & good.
SO, does it mean both them are correct?

6 
velocity.SetOrder(TET,4)

is stable,

while

velocity.SetOrder(TRIG,3)

needs conditions on the mesh, but this seems non-trivial to me.

Also in terms of matrix entries (and static condensation), the volume bubble is advantageous.

7

Hello, it appears that the “Set” method of functions in 3D Mini element space is yielding unexpected “nan” values. Please refer to the example below, where the output is:

2031 nan in u_mini
   0 nan in u_p1

The version of ngsolve used is NGSolve-6.2.2305.

from ngsolve import *
from netgen.occ import OCCGeometry, Box
import numpy as np

geo = OCCGeometry(Box((0,0,0), (1, 1, 1)))
mesh = Mesh(geo.GenerateMesh(maxh=0.1))

V1 = VectorH1(mesh, order=1)
V2 = VectorH1(mesh, order=1)
V2.SetOrder(TET, 4)
V2.Update()

u = CF((x, y, z))
u_p1 = GridFunction(V1)
u_mini = GridFunction(V2)
u_p1.Set(u)
u_mini.Set(u)

print(f"{len(np.where(np.isnan(u_mini.vec))[0]):4d} nan in u_mini")
print(f"{len(np.where(np.isnan(u_p1.vec))[0]):4d} nan in u_p1")