Mesh Vertex numeration in FEniCS

Hello everyone,
I have some doubts about the vertex numeration of each element of a mesh in FEniCS.

I want to reproduce the FreeFem++ code:

 // Current metric
 matrix M = [[(Th[k][1].x-Th[k][0].x)/sqrt(3.),
                (2*Th[k][2].x-Th[k][0].x-Th[k][1].x)/3],
                [(Th[k][1].y-Th[k][0].y)/sqrt(3.),
                 (2*Th[k][2].y-Th[k][0].y-Th[k][1].y)/3]
                 ];

using the FEniCS code:

 Th_cells = np.array(Th.cells())
 cell_k = Th_cells[k]

 M[0, 0] = (Th_coord[cell_k[1]][0] - Th_coord[cell_k[0]][0]) / sqrt(3.)
 M[0, 1] = (2 * Th_coord[cell_k[2]][0] - Th_coord[cell_k[0]][0] -
                   Th_coord[cell_k[1]][0]) / 3
 M[1, 0] = (Th_coord[cell_k[1]][1] - Th_coord[cell_k[0]][1]) / sqrt(3.)
 M[1, 1] = (2 * Th_coord[cell_k[2]][1] - Th_coord[cell_k[0]][1] -
                   Th_coord[cell_k[1]][1]) / 3

Which computes a local metric… but I’ve found that the vertices in each element k are ordered in a different in the two softwares… I cannot find a way to map this ordering… can anyone help me? Thank you !

It is unclear to me what Th and Th_coord is in your FEniCS (or is it FEniCSx) code.

Please make a minimal reproducible example, say on a 1x1 unit square?

1 Like

Sorry, here a reproducible example, on a simple mesh:

from fenics import *
import numpy as np

Th = RectangleMesh(Point(0.0, 0.0), Point(1, 1), 10, 10)


nOfElements = Th.num_cells()

Th_coord = np.array(Th.coordinates())
Th_cells = np.array(Th.cells())

for k in range(nOfElements):

    cell_k = Th_cells[k]

    M = np.zeros((2, 2))

    M[0, 0] = (Th_coord[cell_k[1]][0] - Th_coord[cell_k[0]][0]) / sqrt(3.)
    M[0, 1] = (2 * Th_coord[cell_k[2]][0] - Th_coord[cell_k[0]][0] -
               Th_coord[cell_k[1]][0]) / 3
    M[1, 0] = (Th_coord[cell_k[1]][1] - Th_coord[cell_k[0]][1]) / sqrt(3.)
    M[1, 1] = (2 * Th_coord[cell_k[2]][1] - Th_coord[cell_k[0]][1] -
               Th_coord[cell_k[1]][1]) / 3

    print(M)

    # Some operations on M

I don’t get how can I do the link with FreeFem, since it seems that the vertices are numerated differently… thank you in advance