Interpretation Stiffness matrix P1edge/P1

[Loic (Balazi)]

Hello all,

I need some help to understand a result. I have the following mesh.

I define the FE space P1edge on the Coarse mesh (mesh done only with the red edges) TH.
I define the FE space P1 on the Fine mesh (all edges/vertices) Th. The boundary edges are labelled 3.

Then I compute the stiffness matrix using these two spaces.

fespace VH(TH,P1edge);
fespace Vh(Th,P1);
varf vtest(u,v) = int1d(Th, 3)(u*v);
matrix  = vtest(VH, Vh);

1.179078516   0   0 1.179078516   0   0
0   0 0.389630518 1.610369482   0   0
0   0 1.52782145 0.47217855   0   0
0   0 1.16256891   0 1.178643716   0
1.610369482 0.389630518   0   0   0   0
0   0   0   0   0   0
0   0   0   0 1.592555475 0.4074445248
0.47217855 1.52782145   0   0   0   0
0   0   0   0 0.4074445248 1.592555475
0 1.16256891   0   0   0 1.178643716

I don’t understand why for example the second term in the first line (the term A(1,2) is nul and same for the lines 4 and 10), since each DoF of Vh located on the boundary edges should interact with the two DoF of VH, as is in the second line, I suppose.

Could someone, explain me (detailed me) the calculation that it is done ?
I can interpret the case in which I use P0edge but not the case in which I use P1edge.

Thank you in advance for your help,

Best regards,

Loïc,

[frederichecht (Frédéric Hecht)]

Without the numbering of mesh it is impossible to given a answer.

And I do not see what is mesh TH
the matrix is 10x 6 => only 3 edge in the mesh TH => MeshS but int1d(Th, 3)(u*v);
do not work with this mesh.

Please send the 2 mesh Th and TH.

[Loic (Balazi)]

Hello,

Sorry I have done a mistake.

My coarse mesh TH is :

My fine mesh Th is:

I use the FE space P1Edge on the coarse mesh TH (2 basis per edge, thus 6 basis/6 column in the matrix).

I use the FE space P1 on the fine mesh Th (thus 10 DoF/10 lines in the matrix).

Then I compute the stiffness matrix between these two spaces,

fespace VH(TH,P1edge);
fespace Vh(Th,P1);
varf vtest(u,v) = int1d(Th)(u*v/lenEdge*2);
matrix  = vtest(VH, Vh);

and I obtain the following stiffness matrix:

1.179078516   -0.1790785162    -0.1790785162   1.179078516   0   0
0   0 0.389630518 1.610369482   0   0
0   0 1.52782145 0.47217855   0   0
0   0 1.16256891 -0.1625689098 1.178643716 -0.1786437161
1.610369482 0.389630518   0   0   0   0
0   0   0   0   0   0
  0   0   0   0 1.592555475 0.4074445248
0.47217855 1.52782145   0   0   0   0
 0   0   0   0 0.4074445248 1.592555475
 -0.1625689098 1.16256891   0   0 -0.1786437161 1.178643716

I want to multiply v in Vh by a linear polynomial that is defined on the edge of the coarse mesh TH.
I have some difficulties to interpret the results.

Could you explain me (detailed me) the calculation that it is done ?
For a given Edge E in TH, is all the DoF of Vh multiplied by the same polynomial ?

Thank you in advance for your help,

Best regards,

Loïc,

[frederichecht (Frédéric Hecht)]

load "msh3"
load "Element_Pkedge"
mesh Th = trunc(square(3,3),x-y<0);
mesh TH = trunc(square(1,1),x-y<0);
plot(Th,TH,wait=1);
int[int] l4=[4];
//mesh TH=Th;//extract(Th,label=l4);
fespace VH(TH,P1edge);
fespace Vh(Th,P1);
varf vtest(u,v) = int1d(Th, 4)(u*v);
matrix H = vtest(VH, Vh);
cout << H << endl; 	

le matrix :

  -- Square mesh : nb vertices  =16 ,  nb triangles = 18 ,  nb boundary edges 12
  -- Square mesh : nb vertices  =4 ,  nb triangles = 2 ,  nb boundary edges 4
#  HashMatrix Matrix (COO) 0x7fcb34caf800
#    n       m        nnz     half     fortran   state  
10 6 8 0 0 0 0 
         0         2 -0.028929219009093901788
         0         3 0.19559588567576058349
         1         2 0.07044162180172902632
         1         3 0.26289171153160428851
         3         2 0.26289171153160428851
         3         3 0.07044162180172902632
         6         2 0.19559588567576058349
         6         3 -0.028929219009093894849

times: compile 0.006945s, execution 0.00024s,  mpirank:0

[Loic (Balazi)]

Hello @frederichecht ,

Thank you for you reply,

What did you illustrate through this example ?

Best regards,

Loïc,

[frederichecht (Frédéric Hecht)]

I do that to understand and see the they are a bug.

So a build 2 mesh Th and TH the same caracteristic but not the same numbering
.I a see I make at erreur in test I

varf vtest(u,v) = int1d(Th)(u*v);

So we compute integral on boundary edge of Th ( small boundary edge of length 1/3 or \sqrt(2)/3)

And we use a quadrature formulas qf3pE (with 3 points exact of P5 ) on each edges …

so in my example 4 vertex (0,1,3,6) on the big edge (the 2 dof are 2,3 and the edge number is 1)

the P1edge are orthogonal polynom of edge such that le sum is 1.
So the interpolation point are (1 ± \sqrt{1/3} )/2

[Loic (Balazi)]

Hello @frederichecht,

Thank you for your message. I understand better the calculation that is done.

One last thing:

1. You use the quadrature formulae qfpE (3 points of interpolations) but you give only two of them (1±√1/3)/2. I think the last one is 0.5, isn’t it ?

2. “The P1edge are orthogonal polynomial of edge such that le sum is 1”: it seems that for a given edge this couple is not unique. Do you know their explicit expression ?

Best regards,

Loïc,

[frederichecht (Frédéric Hecht)]

No it is qfpE2 (2 points pf interpolation)
then the it is Lagrange polynom corresponding to this 2 points (1±√1/3)/2,

[Loic (Balazi)]

Thank you for your message,

and about their expressions ? (I would like to do some calculation by hand to check my implementation).

Best regards,

Loïc,

[frederichecht (Frédéric Hecht)]

You put you 6 points on vertical edge
this point a the gauss point on each sub edge => exact integration for P3 polynom.

on each edge of size 1/3 4 function
first classical P1 on edge l2=x, l1=1-x on 0,1 big edge
second = f1, f2 function on gauss point (1±√1/3)/2. small edge

It is you job to do the computation but freefem++ make no error!

[Loic (Balazi)]

Thank you for your explanations,

These elements of explanation are enough for me to make some calculations,

Best regards,

Loïc,