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

**URL:** <https://community.freefem.org/t/order-of-variables-in-matrix-generated-with-varf-with-a-p1-p1-p0-fespace/3489>\
**Category:** General Discussion\
**Created:** [September 17, 2024, 2:51pm UTC](https://community.freefem.org/t/order-of-variables-in-matrix-generated-with-varf-with-a-p1-p1-p0-fespace/3489 "2024-09-17T14:51:22Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![djahid](https://avatars.discourse-cdn.com/v4/letter/d/bbe5ce/32.png) [@djahid](https://community.freefem.org/u/djahid)\
**Post date:** [September 17, 2024, 2:51pm UTC](https://community.freefem.org/t/order-of-variables-in-matrix-generated-with-varf-with-a-p1-p1-p0-fespace/3489/1 "2024-09-17T14:51:22Z")

</div>

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

```auto
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:

```auto
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:

```auto
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

```auto
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.

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [September 20, 2024, 11:17am UTC](https://community.freefem.org/t/order-of-variables-in-matrix-generated-with-varf-with-a-p1-p1-p0-fespace/3489/2 "2024-09-20T11:17:18Z")

</div>

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

> **[FreeFEM-documentation.pdf](https://doc.freefem.org/pdf/FreeFEM-documentation.pdf)**
>
> 33.81 MB

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

```auto
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…

---

<div class="post-metadata">

**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\
**Post date:** [September 20, 2024, 1:47pm UTC](https://community.freefem.org/t/order-of-variables-in-matrix-generated-with-varf-with-a-p1-p1-p0-fespace/3489/3 "2024-09-20T13:47:04Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![djahid](https://avatars.discourse-cdn.com/v4/letter/d/bbe5ce/32.png) [@djahid](https://community.freefem.org/u/djahid)\
**Post date:** [September 20, 2024, 3:11pm UTC](https://community.freefem.org/t/order-of-variables-in-matrix-generated-with-varf-with-a-p1-p1-p0-fespace/3489/4 "2024-09-20T15:11:51Z")

</div>

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:

```auto
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?

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [September 20, 2024, 4:53pm UTC](https://community.freefem.org/t/order-of-variables-in-matrix-generated-with-varf-with-a-p1-p1-p0-fespace/3489/5 "2024-09-20T16:53:36Z")

</div>

I think you can do as

```auto
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

```auto
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);

```

---

<div class="post-metadata">

**Author:** ![pi.3.141](https://avatars.discourse-cdn.com/v4/letter/p/ecccb3/32.png) [@pi.3.141](https://community.freefem.org/u/pi.3.141)\
**Post date:** [May 29, 2025, 2:28pm UTC](https://community.freefem.org/t/order-of-variables-in-matrix-generated-with-varf-with-a-p1-p1-p0-fespace/3489/6 "2025-05-29T14:28:15Z")

</div>

Hello,  
I have a closely related question. To make it more precise, assume we have

fespace UP(Th,[P1b,P1]);  
fespace XY(Th,[P1,P1]);  
varf UPxXY([u,p],[xh,yh]) = int(Th)(…); // for example …=u_xh+u_yh+p_xh+p_yh);  
matrix M=UPxXY(UP,XY);

One can take UP instead of XY, or vice-versa. Motivated by problem specifics, I want to modify M, say its element M(i,j). For this, I need to know the dof corresponding to ith row and to jth column.

Wondering how to get this information directly or a posteriori. Definitely this depends on the algorithm the developers have used for computing M. I hope anyone can help.

Thank you!