Crouzeix-Raviart FE of degree 3

General Discussion
1

Hello there,
I am trying to implement a new finite element, Crouzeix-Raviart of degree 3 in dimension 2. But I actually don’t speak C++, so I kind of tried to combine P3 and P1nc (Crouzeix-Raviart of degree 1) and changed the basis of P1nc according to the definition.
Element_P3nc.zip (2.6 KB)
Obvouisly it doesn’t work yet. Could anyone help and take a look at what I did wrong and what I should change?
Thank you very much!

Definition of CR-element of degree 3:
Theorem 3.3 in https://arxiv.org/pdf/2105.14981.pdf
So we have the normal Lagrange basis of degree 3, except on the vertices, where we use the edge bubble, 1/2 * (5*(1 - (lambda_Ai)*2)^3 - 3 * (1 - (lambda_Ai)*2)) with lambda being the baricentric coordinate of the node Ai.

2

I think the new finite element

  • plugin Element_P3pnc for new 2d finite element P3pnc (P3 + 2 bulles) nonconforming (continuite of P2 mod)
    and add 2 examples with this new finite element
    examples/plugin/cavityNewtowP3pnc.edp examples/plugin/testFE-P3pnc.edp

is this finite element

for more detail see tiger reach result of string “P3pnc” in the community.

3

Oh, I didn’t see the new update.
Yes, it seems like the same I tried to implement.
Thank you and have a nice day!

4

@frederichecht unfortunately, the P3_pnc does not really fit since in our case, we have 10 nodes instead of 12 (the same on the vertices/edged, but only one inside). Would it possible/easy to just erase the 2 additional nodes from the middle in the skript of P3_pnc?

These are the nodes we need: {R2(0. , 0. ), R2( 1., 0.), R2(0. , 1.),
R2(2 / 3., 1 / 3.), R2(1 / 3., 2 / 3.), R2(0. , 2 / 3.),
R2(0. , 1 / 3.), R2(1 / 3., 0. ), R2(2 / 3., 0. ),
R2(1 / 3., 1 / 3.)}

5

yes, but i have no time do this effectively I have add a bubble to P3

the only work is to recompile the corresponding array (10x10)

double C1[12][12] ={
{-1.2,   6,        1.2,    -8.4,    10.8,     10.8,    -46.8,   -10.8,    20.4,    12,     0,     -84  },
{-1.2,   1.2,      6,      -73.2,   75.6,     73.2,    -63.6,   -73.2,    37.2,    12,     -84,   -84  },
{1.2,    -15.6,    -15.6,  298.8,   -205.2,   -306,    349.2,    198,     -154.8,  -180,   252,   504  },
{1.2,    -1.2,     6,      73.2,    -63.6,    -73.2,   75.6,     37.2,    -73.2,   12,     84,    84   },
{6,      -1.2,     1.2,    10.8,    -46.8,    -8.4,    10.8,     20.4,    -10.8,   12 ,    84,    0    },
{-15.6,  1.2,      -15.6,  -306,    349.2,    298.8,   -205.2,   -154.8,   198,    -180,  -504,   -252 },
{6,      1.2,      -1.2,   -46.8,   10.8,     20.4,    -10.8,    -8.4,     10.8,   12,    -84,    0    },
{1.2,    6,        -1.2,   20.4,    -10.8,    -46.8,   10.8,     10.8,     -8.4,   12,    0,      84   },
{-15.6,  -15.6,    1.2,    97.2,    -54,      97.2,    -54,      46.8,     46.8,   -180,  252,    -252 },
{-9.6,   4.8,      4.8,    189.6,   21.6,     -192,    86.4,     -24,     -81.6,   60,    84,     168  },
{4.8,    -9.6,     4.8,    -192,    86.4,     189.6,   21.6,     -81.6,   -24,     60,   -168,    -84  },
{4.8,    4.8,      -9.6,   2.4,     -108,     2.4,    -108,      105.6,   105.6,   60,   84,      -84  },

}; */

which computed in exemple

testFE-P3pnc.edp (5.7 KB)

after the modification is trivial!
the array is

	  -6   6   6  54  54 -54  18 -54 -42  12
	  -6   6   6  54  54 -54 -42 -54  18  12
	  30 -30 -30 -270 -270 270  90 270  90 -180
	   6  -6   6 -54 -42  54  54  18 -54  12
	   6  -6   6 -54  18  54  54 -42 -54  12
	 -30  30 -30 270  90 -270 -270  90 270 -180
	   6   6  -6  18 -54 -42 -54  54  54  12
	   6   6  -6 -42 -54  18 -54  54  54  12
	 -30 -30  30  90 270  90 270 -270 -270 -180
	   0   0   0   0   0   0   0   0   0  60
6

I have code this finite element in GitHub - FreeFem/FreeFem-sources at develop
see plugin Element_P3nc.cpp add the small test examples/plugin/testFE-P3nc.edp

This is a very first version

7

@frederichecht Perfect!
But I noticed, that some of the basis are not what I have in my case. So instead of l0^3 it would be 1/2 * (5* l0^3 -3*l0) and the same for l1^3 and l2^3. I actually just replaced this to calculate the Matrix C1 and then also in Element_P3nc for l12 (line 147). Does it actually work with just replacing the basis? Do I also need to chance something for D12 (line 229) because of the change of the basis or anything else?

The new Matrix C1 would be:
calculate_C-Matrix.edp (2.4 KB)

And my changes to the Element_P3nc:
Element_P3nc.cpp.zip (3.3 KB)

Does this all make sense programming wise?

8

There was an error in calculating the Matrix C1, so here is the new version:

calculate_C-Matrix.edp (2.4 KB)

And with this the new Element_P3nc:
Element_P3nc.cpp.zip (3.2 KB)

But while testing this with testFE-P3nc.edp (3.8 KB), I get an error for the last part:
err=136418
current line = 150
Exec error : exec assert
– number :1
Exec error : exec assert
– number :1
err code 8 , mpirank 0

9

I sorry I Have no computer so I will correct in week,