Node numbering of a refined mesh (truncated with split = 2) for P2 elements post-processing

General Discussion
1

I’m doing a post-processing work which outputs the node values of a 2D P2 element using VTK. Since Freefem++ currently doesn’t support outputing the P2 elements (VTK order = 2) directly, one possible way is to first refine the mesh (truncated with split = 2), then output the refined mesh with P1 elements (VTK order = 1).

The mesh refinement is easy to implement. Assuming we have an original 2D mesh called Th, and a refined mesh ThR is established by using the command:

mesh ThR = Th;
ThR = trunc(ThR, true, split=2);

But, it is a bit difficult to establish the node numbering mapping between the original Th with P2 elements and the refined ThR with P1 elements.

The node numbering for the original P2 elements can be accessed as follows:

mesh Th = square(2, 2);

fespace Wh(Th, P2);

for (int k = 0; k < Th.nt; k++) {
    for (int i = 0; i < 6; i++) {
        cout << Wh(k, i) << " ";
    }
    cout << endl;
}

However, the node numbering rule for the refined mesh ThR with P1 elements is not very clear. The node numbering for the common nodes of Th & ThR are the same, but the rule for the newly added nodes is not clear.

Here is an example.

Fig. 1 gives the P1 node numbering of the mesh Th.

1
1832×838 11.2 KB

Fig. 2 gives the P1 node numbering of the refined mesh ThR.

2
2832×838 21.6 KB

Fig. 3 gives the P2 node numbering of the mesh Th.

3
3832×838 30.7 KB

Of course, we can find the mapping relation by comparing each nodes’ coordinates, and determine whether they’re the same one by Euclidean distance like

sqrt((xi - xj)^2 + (yi - yj)^2 + (zi - zj)^2) < epsilon

Is there any more efficient way to address this problem, since the above method maybe time cost expensive when there are millions of nodes?

2

I don’t have a simple method to solve your problem, but some search trough the neighbouring vertices of a vertex i of ThR (which is not a vertex of Th) should give the two that belong to Th (then the vertex i is the middle of these two) . Note that a vertex of ThR of index i is a vertex of Th if and only if i<Th.nv.
By the way, how do you produce these nice images including vertex indices ?

3

Thanks for your reply.

I guess an interpolation operator = maybe can solve the problem indirectly. Using the interpolation, the nodes with same coordinates maybe have same values.

// Original P2 solution
fespace Vh(Th, P2);
Vh u;

// Refined P1 solution
fespace Rh(ThR, P1);
Rh uR = u;

In this way, we don’t focus on the mapping between the P2 node indices and the refined P1 node indices, but just transferring nodal values of the P2 nodes to the refined P1 nodes.

Regarding the figures with vertex indices, Fig.1 & Fig. 2 are produced directly by Paraview, and Fig. 3 is produced by hand.

4

Your idea can indeed be converted to a real mapping:

mesh Th = square(2, 2);
mesh ThR = trunc(Th, true, split=2);

fespace Vh(Th,P2);
fespace Rh(ThR,P1);

matrix convert=interpolate(Vh,Rh);
cout << convert;

Then convert(i,j) is 1 (0 otherwise)
if (i,j) is a couple with i the Vh label and j the Rh label, of the same vertex.

About Fig 1 & 2, do you save the meshes as .vtu ? How do you plot the label numbers in paraview?