ElementTetP3 and ElementTetP4

[xida2020]

Thanks to #651 and #741, I create ElementTetP3 and ElementTetP4 for the tetrahedral mesh.
Then I use the two elements to solve Laplace eigenvalue problem, and the results are expected.

The results using P3 elements on the uniform meshes:
Iteration, d.o.f,the first approximation eigenvalue, error
0,125,[30.56421926],[0.95540606]
1,729,[29.62225938],[0.01344617]
2,2197,[29.60945679],[0.00064359]
3,4913,[29.60888247],[6.92623446e-05]
4,9261,[29.60882524],[1.2034608e-05]
5,15625,[29.60881606],[2.85504285e-06]
The results using P4 elements on the uniform meshes:
Iteration, d.o.f,the first approximation eigenvalue, error
0,125,[30.56421926],[0.95540606]
1,729,[29.62225938],[0.01344617]
2,2197,[29.60945679],[0.00064359]
3,4913,[29.60888247],[6.92623447e-05]
4,9261,[29.60882524],[1.2034608e-05]
5,15625,[29.60881606],[2.85504284e-06].

However, when I use the two elements solve the eigenvalue problem on the adaptive meshes, there are something wrong in the codes. I have no idea which parts go wrong.
The results using P3 elements on the adaptive meshes:
iteration, mesh parameter h, d.o.f, the first eigenvalue, error
11,0.34641016151377557,4096,29.610201448711088,0.001388245443013858
10,0.34641016151377557,5743,34.443115113394015,4.834301910125941
9,0.3464101615137755,6409,35.80059383901864,6.191780635750565
8,0.3464101615137755,6961,37.299082439771475,7.690269236503401
7,0.3464101615137755,7627,37.38224725261639,7.773434049348314
6,0.3464101615137755,10477,38.1063522038288,8.497539000560725
5,0.3464101615137755,10933,37.88516277898642,8.276349575718346
4,0.3464101615137755,15706,38.30431903207835,8.695505828810273
3,0.3464101615137755,19801,38.769389504354116,9.160576301086042
2,0.3464101615137755,22972,38.8688762403918,9.260063037123725
1,0.3464101615137755,26455,38.66783041160825,9.059017208340176
0,0.34641016151377535,37546,38.77234628275499,9.163533079486914

The results using P4 elements on the adaptive meshes:
iteration, mesh parameter h, d.o.f, the first eigenvalue, error
5,0.34641016151377557,9261,29.608825237876108,1.2034608033673067e-05
4,0.34641016151377557,13997,39.87218445498515,10.263371251717079
3,0.34641016151377557,17621,40.68697153428219,11.07815833101412
2,0.34641016151377557,21785,39.92128229332504,10.312469090056968
1,0.34641016151377557,26881,38.67292861657085,9.064115413302776
0,0.34641016151377557,45685,39.39686921698523,9.788056013717156

In addition, using P2 elements on the adaptive mesh, the resulting eigenvalues convergence to the exact one.
iteration, mesh parameter h, d.o.f, the first eigenvalue, error
11,0.34641016151377557,1331,29.706286231004054,0.09747302773597966
10,0.3464101615137755,2027,29.64864481627044,0.039831613002366595
9,0.2828427124746191,2871,29.628758231498857,0.019945028230782924
8,0.2828427124746191,3681,29.621735385771252,0.012922182503178448
7,0.2828427124746191,4323,29.61733919810704,0.008525994838965545
6,0.28284271247461906,4641,29.616680199672533,0.007866996404459314
5,0.20000000000000007,7425,29.61482206177459,0.0060088585065152245
4,0.20000000000000007,9897,29.613791385720766,0.004978182452692437
3,0.20000000000000007,15793,29.61119642109418,0.0023832178261073977
2,0.17320508075688779,21461,29.6101262149508,0.0013130116827255733
1,0.14142135623730964,30749,29.609464822981806,0.0006516197137322877
0,0.14142135623730964,34373,29.609343738098502,0.0005305348304283086

[kinnala]

Thanks for the clear example. I need to look more carefully, but I have a guess:

• There is no routine in scikit-fem to orient facets and edges of a tetrahedral mesh.
• If there are more than one DOF per edge, the orientation is important to preserve continuity.
• For P3 and P4 there are more than one DOF per edge and therefore orientation will be required.

Why it works for uniform mesh? My guess is that the uniform refinement algorithm accidentally preserves orientation of the original mesh. This is not by design and not tested in any way!

This is just a guess. I'll look more carefully next week.

[xida2020]

Thanks for your quick reply.

[kinnala]

I confirmed the orientation issue by a simple test.

You can see visually how those sine-like edge bubbles are wrong way around.

I don't have a fix yet but this is obviously something I want to get fixed in the future.

[kinnala]

Here is my snippet:

from skfem import *

m = MeshTet().refined()
import numpy as np
m1 = MeshTet.init_tensor(
    np.linspace(0, 1, 3),
    np.linspace(0, 1, 3),
    np.linspace(0, 1, 3),
)

basis = Basis(m, ElementTetP3())

from skfem.models import laplace
from skfem.models import mass

A = laplace.assemble(basis)
M = mass.assemble(basis)

L, x = solve(A, M)

basis.plot(np.real(x[:, 0]), visuals='vedo', nrefs=3).show()

It uses your ElementTetP3 implementation.

[xida2020]

Thank you @kinnala , and it would be nice to have the possibility of getting the proper orientation of ElementTetP3 or ElementTetP4 in the future.

[xida2020]

May I add the two tetrahedral Lagrangian elements to scikit-fem?

[kinnala]

Yes, but I think we should mark them with some sort experimental feature warning until the orientation issue is solved. I think we cannot get all the usual element tests to pass until this.

[kinnala]

There is example in skfem.autodiff of the warning.