Mixed FEM - numeration

Dear FreeFEM users,

For a mixed element, such as P2/P1 element for velocity (in space Vh) and pressure (in space Ph) pair. If I loop for (i=0; i<Vh.ndof; i++), how can I know which DOF includes pressure and which one does not ?

If I do varf rhs([u,v,p],[uh,vh,ph]) =…, and later do
real[int] b = rhs(0, Rh), with Rh=[Vh, Vh, Ph]
How are variables u,v, and p numerated in b ?

Alternatively, if I solve for
w=A^-1 * b
how can I give w to [u, v, p], what’s the order?

Thanks in advance!


fespace Wh(Th, [P2, P1]);
Wh [u, v] = [1, 2];
cout << u[] << endl;

You’ll see that the degrees of freedom are interleaved.

Thank you very much! This is interesting, so there is no particular order.
why not u first then v, or [u1, v1, u2, v2, u3, v3…] ?
The interleaved arrangement can guarantee a good matrix structure?

Yes, much better. For example, take a [P1, P1] fespace, instead of having a 2-by-2 block matrix [[A, B]; [C, D]], you end up with a sparse matrix A where all unknowns are 2-by-2 dense small blocks. This is coherent, for example, with PETSc MATBAIJ format.