Order of variables in matrix generated with varf with a [P1, P1, P0] fespace

General Discussion
1

Hello,
I’m trying to compare the matrix generated from the following code with another code I have that solves the same problem

mesh Th=square(2, 2);
fespace Xh(Th,P1);
fespace Mh(Th,P0);
fespace XhxXhxMh(Th,[P1,P1,P0]);

varf vfStokes ([u1,u2,p],[v1,v2,q]) =
    int2d(Th)(
             alpha*( u1*v1 + u2*v2)
            + nu * ( dx(u1)*dx(v1) + dy(u1)*dy(v1)
            +        dx(u2)*dx(v2) + dy(u2)*dy(v2) )
            + p*q*(1e-4)
            - p*dx(v1) - p*dy(v2)
            - dx(u1)*q - dy(u2)*q
           );
varf b([u1,u2,p],[v1,v2,q]) = int2d(Th)(  u1*v1+u2*v2+p*q*0.) ;
matrix A= vfStokes(XhxXhxMh, XhxXhxMh, solver="SPARSESOLVER", factorize=0); 
matrix B= b(XhxXhxMh, XhxXhxMh, solver=CG, eps=1e-20);

Here the A matrix is:

0: 1 0 0.25 -0.5 0 0 0 0 -0.5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1: 0 1 -0.25 0 -0.5 0 0 0 0 -0.5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
2: 0.25 -0.25 1.25e-05 -0.25 0 0 0 0 0 0.25 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
3: -0.5 0 -0.25 2 0 0.25 0 0 0 0 0 -0.5 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 
4: 0 -0.5 0 0 2 -0.25 -0.25 0 0 0 0 0 -0.5 0 -1 0 0 0 0 0 0 0 0 0 0 0
5: 0 0 0 0.25 -0.25 1.25e-05 0 0 0 0 0 -0.25 0 0 0.25 0 0 0 0 0 0 0 0 0 0 0
6: 0 0 0 0 -0.25 0 1.25e-05 0 0.25 0 0 0 0 -0.25 0.25 0 0 0 0 0 0 0 0 0 0 0
7: 0 0 0 0 0 0 0 1.25e-05 0.25 -0.25 0 0 0 -0.25 0 0 0.25 0 0 0 0 0 0 0 0 0
8: -0.5 0 0 0 0 0 0.25 0.25 2 0 0 0 0 -1 0 -0.5 0 0 0 0 0 0 0 0 0 0
9: 0 -0.5 0.25 0 0 0 0 -0.25 0 2 0 0 0 0 -1 0 -0.5 0 0 0 0 0 0 0 0 0
10: 0 0 0 0 0 0 0 0 0 0 1.25e-05 0 -0.25 0.25 0 0 0 0 0 -0.25 0.25 0 0 0 0 0
11: 0 0 0 -0.5 0 -0.25 0 0 0 0 0 1 0 0 0 0 0 0 0 -0.5 0 0 0 0 0 0
12: 0 0 0 0 -0.5 0 0 0 0 0 -0.25 0 1 0 0 0 0 0 0 0 -0.5 0 0 0 0 0
13: 0 0 0 -1 0 0 -0.25 -0.25 -1 0 0.25 0 0 4 0 0 0 0.25 0 -1 0 -1 0 0 0 0
14: 0 0 0 0 -1 0.25 0.25 0 0 -1 0 0 0 0 4 0 0 -0.25 -0.25 0 -1 0 -1 0 0 0
15: 0 0 0 0 0 0 0 0 -0.5 0 0 0 0 0 0 1 0 0 0.25 0 0 -0.5 0 0 0 0
16: 0 0 0 0 0 0 0 0.25 0 -0.5 0 0 0 0 0 0 1 0 0 0 0 0 -0.5 0 0 0
17: 0 0 0 0 0 0 0 0 0 0 0 0 0 0.25 -0.25 0 0 1.25e-05 0 -0.25 0 0 0.25 0 0 0
18: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.25 0.25 0 0 1.25e-05 0 0 -0.25 0.25 0 0 0
19: 0 0 0 0 0 0 0 0 0 0 -0.25 -0.5 0 -1 0 0 0 -0.25 0 2 0 0 0 -0.5 0 0 
20: 0 0 0 0 0 0 0 0 0 0 0.25 0 -0.5 0 -1 0 0 0 0 0 2 0 0 0 -0.5 -0.25
21: 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 -0.5 0 0 -0.25 0 0 2 0 -0.5 0 0.25
22: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 -0.5 0.25 0.25 0 0 0 2 0 -0.5 0
23: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.5 0 -0.5 0 1 0 -0.25
24: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.5 0 -0.5 0 1 0.25
25: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.25 0.25 0 -0.25 0.25 1.25e-05

I’m trying to understand the order of the u1, u2, p variables in this matrix. Because I had thought it would be ordered like [u1_1, u1_2, u1_3, …, u1_NV, u2_1, u2_2, u2_3, …, u2_NV, p_1, p_2, …, p_NT] but that’s not the case.

here is the matrix when I generated it using this order:

1 -0.5 0 -0.5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.25 0 0 0 0 0 0 0 
-0.5 2 -0.5 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.25 0 0.25 0 0 0 0 0 
0 -0.5 1 0 0 -0.5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.25 0 0 0 0 0 
-0.5 0 0 2 -1 0 -0.5 0 0 0 0 0 0 0 0 0 0 0 0 0.25 0 0 0.25 0 0 0 
0 -1 0 -1 4 -1 0 -1 0 0 0 0 0 0 0 0 0 0 0 -0.25 0 0.25 -0.25 0 0.25 0 
0 0 -0.5 0 -1 2 0 0 -0.5 0 0 0 0 0 0 0 0 0 0 0 0 -0.25 0 0 -0.25 0 
0 0 0 -0.5 0 0 1 -0.5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.25 0 0 
0 0 0 0 -1 0 -0.5 2 -0.5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.25 0 0.25 
0 0 0 0 0 -0.5 0 -0.5 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.25 
0 0 0 0 0 0 0 0 0 1 -0.5 0 -0.5 0 0 0 0 0 0 0.25 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 -0.5 2 -0.5 0 -1 0 0 0 0 0.25 0 0 0.25 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 -0.5 1 0 0 -0.5 0 0 0 0 0 0.25 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 -0.5 0 0 2 -1 0 -0.5 0 0 0 -0.25 0 0 0 0.25 0 0 
0 0 0 0 0 0 0 0 0 0 -1 0 -1 4 -1 0 -1 0 -0.25 0 0 -0.25 0.25 0 0 0.25 
0 0 0 0 0 0 0 0 0 0 0 -0.5 0 -1 2 0 0 -0.5 0 0 -0.25 0 0 0 0.25 0 
0 0 0 0 0 0 0 0 0 0 0 0 -0.5 0 0 1 -0.5 0 0 0 0 0 0 -0.25 0 0 
0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 -0.5 2 -0.5 0 0 0 0 -0.25 0 0 -0.25 
0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.5 0 -0.5 1 0 0 0 0 0 0 -0.25 0 
0.25 -0.25 0 0 0 0 0 0 0 0 0.25 0 0 -0.25 0 0 0 0 0.0001 0 0 0 0 0 0 0 
0 0 0 0.25 -0.25 0 0 0 0 0.25 0 0 -0.25 0 0 0 0 0 0 0.0001 0 0 0 0 0 0 
0 0.25 -0.25 0 0 0 0 0 0 0 0 0.25 0 0 -0.25 0 0 0 0 0 0.0001 0 0 0 0 0 
0 0 0 0 0.25 -0.25 0 0 0 0 0.25 0 0 -0.25 0 0 0 0 0 0 0 0.0001 0 0 0 0 
0 0 0 0.25 -0.25 0 0 0 0 0 0 0 0 0.25 0 0 -0.25 0 0 0 0 0 0.0001 0 0 0 
0 0 0 0 0 0 0.25 -0.25 0 0 0 0 0.25 0 0 -0.25 0 0 0 0 0 0 0 0.0001 0 0
0 0 0 0 0.25 -0.25 0 0 0 0 0 0 0 0 0.25 0 0 -0.25 0 0 0 0 0 0 0.0001 0
0 0 0 0 0 0 0 0.25 -0.25 0 0 0 0 0.25 0 0 -0.25 0 0 0 0 0 0 0 0 0.0001

which is basically a matrix of the form

A      0    Bx.T
0      A    By.T
Bx     By   eps*I

So my question is where can I find this variable order information of the matrix A or the fespace XhxXhxMh? so I can compare the matrices

Thank you very much for your help and have a good day.

2

Hello,
I think there is no official rule concerning the ordering of dof.
Nevertheless you can get it a posteriori.
In the FreeFem documentation

p.116-120 you can find ways to get the mesh info,
like Th[k][i] to get the vertex number (for i=0,1,2) in triangle k
and Th[k][i].x Th[k][i].y its coordinates
It induces a first vertex (i=0) associated to triangle k,
a second vertex (i=1) and a third vertex (i=2).

For your fespace XhxXhxMh(Th,[P1,P1,P0])
corresponding to x,y components of velocity in P1, pressure in P0
you have the command XhxXhxMh(k,i) for k the triangle number, i=0,1…6
it gives the 7 dof numbers around the triangle k. They are listed as
x component of velocity for first vertex
x component of velocity for second vertex
x component of velocity for third vertex
y component of velocity for first vertex
y component of velocity for second vertex
y component of velocity for third vertex
pressure component for the triangle

For example, to show these values you can write

for (int k=0;k<Th.nt;k++){//loop on triangles
 for (int i=0;i<7;i++){//loop on dof around triangle
  int dof=XhxXhxMh(k,i);
  cout << "k="<< k << " i=" << i << " dof=" << dof << endl;
 }
}

With this you have the needed information. Then you should be able to reorder the indices as you wish…

3

For vectorial finite elements, the degrees of freedom are always interleaved.

4

Unfortunately that still doesn’t really give me the information necessary to understand the interleaving of the u1, u2 and p in the A matrix. I know how to get the dofs for the triangle. My problem is I don’t really know what the order of the variables is going to be in the matrix. for example in the matrix generated with a 2x2 square (9 vertices, 8 triangles) i get the following distribution in the A matrix:

u1, u2, p, u1, u2, p, p, p, u1, u2, p, u1, u2, u1, u2, u1, u2, p, p, u1, u2, u1, u2, u1, u2, p

it does work as simple interleaving when it’s only [P1, P1] giving a pattern of [u1, u2, u1, u2 …] but when there is an element space (ex [P1, P1, P0]) where the P0 is of different shape from the P1 it doesn’t give any recognizable pattern as seen above. Unless this specific information can be accessed from the user side?

5

I think you can do as

int[int] dofindex(XhxXhxMh.ndof);
for (int k=0;k<Th.nt;k++){//loop on triangles
 for (int i=0;i<3;i++){//loop on vertices around triangle to get x component of velocity
  int dof=XhxXhxMh(k,i);
  dofindex(dof)=Th[k][i];//vertex index
 }
 for (int i=3;i<6;i++){//loop on vertices around triangle to get y component of the velocity
  int dof=XhxXhxMh(k,i);
  dofindex(dof)=Xh.ndof+Th[k][i-3];//vertex index + Xh.ndof
 }
  int dof=XhxXhxMh(k,6);// dof of pressure in the triangle
  dofindex(dof)=2*Xh.ndof+k;//triangle index +2*Xh.ndof
}

matrix Aordered=A(dofindex^-1,dofindex^-1);

In this way the new matrix Aordered should correspond to the order
vx[0],…vx[n-1], vy[0], … vy[n-1], p[0], … ,p[m-1]
where n =Th.nv and m=Th.nt
the order correspond to the order of vertices and of triangles respectively.

Alternatively, to avoid the inversion ^-1 you can do

int[int] dofindex(XhxXhxMh.ndof);
for (int k=0;k<Th.nt;k++){//loop on triangles
 for (int i=0;i<3;i++){//loop on vertices around triangle
  dofindex(Th[k][i])=XhxXhxMh(k,i);
 }
 for (int i=3;i<6;i++){//loop on vertices around triangle
  dofindex(Xh.ndof+Th[k][i-3])=XhxXhxMh(k,i);
 }
  dofindex(2*Xh.ndof+k)=XhxXhxMh(k,6);
}

matrix Aordered=A(dofindex,dofindex);