Meshing options to generate number of the sides of and element (tetgen-triangle)

Question

I wrote a finite element code in fortran 90.

This code is really fast, except the meshing process.

I used triangle and tetgen for meshing in 2D and 3D, respectively, so this process is fast, of course.

For example, for the unit square [0,1]x[0,1] in 2D I have a file with the coordinates of its nodes (for example, a mesh with 5 nodes):

1   0.0 0.0  # coordinates of node 1
2   1.0 0.0  # coordinates of node 2
3   1.0 1.0  # coordinates of node 3
4   0.0 1.0  # coordinates of node 4
5   0.5 0.5  # coordinates of node 5

called coordinate.dat, which have 4 elements (triangles) with nodes called element.dat

1   1 5 4  # vertices of triangle 1
2   1 2 5  # vertices of triangle 2
3   2 3 5  # vertices of triangle 3
4   5 2 4  # vertices of triangle 4

I also have a file where each row i is the number of its initial an final node, called edge.dat:

1   1 2  # initial and final node of edge 1
2   2 3  # initial and final node of edge 2
3   3 4  # initial and final node of edge 3
4   4 1  # initial and final node of edge 4
5   1 5  # initial and final node of edge 5
6   5 2  # initial and final node of edge 6
7   2 5  # initial and final node of edge 7
8   5 4  # initial and final node of edge 8

With this files, I need to generate the following information:

(1) Given an element (triangle or tetrahedron), I need to know the number of its sides (edges and faces, respectively). For example, I need to generate the following structure or file, called struct1.dat:

1   5 8 4  # triangle 1 has the edges number 5, 8 and 4
2   1 6 5  # triangle 2 has the edges number 1, 6 and 5
3   6 2 7  # triangle 2 has the edges number 6, 2 and 7
4   7 3 8  # triangle 4 has the edges number 7, 3 and 8

(2) Furthermore, given a side (edge or face) I need to know the element numbers of the 2 elements (or only one if the side is on the boundary) shared by that side. For example, I need to generate the following structure (or file) called struct2.dat:

1   2 0  # edge number 1 is only on element 2
2   3 0  # edge number 2 is only on element 3
3   4 0  # edge number 3 is only on element 4
4   1 0  # edge number 4 is only on element 1
5   1 2  # edge number 5 is sharing by elements 1 and 2
6   3 2  # edge number 6 is sharing by elements 3 and 2
7   4 3  # edge number 7 is sharing by elements 4 and 3
8   1 4  # edge number 8 is sharing by elements 1 and 4

For both of these structures, struct1.dat and struct2.dat, my code is very slow because I used a brute force approach with a lot of loops..

I am looking for an algorithm (a paper, or better: a subroutine in fortran available for download) optimized for this? I want to continue using triangle and tetgen, but I am willing to listen to other options.

Answer 1

I have implemented a method that works pretty well for reconstructing edge information in surfacic meshes (or facet information in volumetric meshes), but I am working in C++. I'll try to explain that in an "abstract algorithmic language" that could (hopefully) be translated into Fortran:

Suppose that I have a triangle mesh of nt triangles.

Step 1: I first construct an array of 3*nt "edge records", each "edge record" has the following fields: (i,j,t) and encodes an edge (i,j) seen from a given triangle (t). I also need to ensure that (i < j).

(note: in Fortran, this could be replaced with three arrays of size 3*nt, one array I[], one array J[] and one array T[], but maybe Fortran90 has "structures"/"records")

To construct the structure:

E : array of 3*nt "edge records"
e = 1
For t = 1 to nt       
  get i,j,k the vertices indices of triangle t        
  E[e].i = min(i,j); E[e].j = max(i,j); E[e].t = t; e = e + 1
  E[e].i = min(j,k); E[e].j = max(j,k); E[e].t = t; e = e + 1
  E[e].i = min(k,l); E[e].j = max(k,l); E[e].t = t; e = e + 1

Step 2: Sort the records by lexicographic sort, i.e. given two records E[e1] and E[e2], E[e1] is supposed to be before E[e2] if E[e1].i < E[e2].i OR (E[e1].i = E[e1].i and E[e1].j < E[e2].j)

(in Fortran, I am unsure of how this can be achieved, you will need to find a generic sorting routine, or maybe will need to implement it on your own if you have three separate I[], J[], T[] arrays. If Fortran90 has "structure/records", maybe the sort() function of the C runtime can be called from Fortran). In the worst case, it is possible to save the array E in a file, call the system sort program on it and read it back.

Step 3: Once the "edges record" array is sorted, the 'same edge' seen from two different triangles will be represented by two records that are contiguous. Then you can generate unique edge indices by iterating on it, and incrementing the current edge index each time a different edge is encountered. Edges that have two triangles incident to them will be represented by a sequence of two "edge records" (i,j,t1), (i,j,t2). Then it is easy to generate the arrays that give for each triangle the indexes of its three edges, and for each edge the index of the (up to two) triangles adjacent to it. It also allows to detect invalid meshes, with more than two triangles incident to the same edge.

There is also an alternative method, that gives for each vertex the list of triangles incident to it. It is done by "chaining" the "corners" of the triangles (using an array of size 3*nt).

C++ implementation in my "geogram" programming library: http://alice.loria.fr/software/geogram/doc/html/namespaceGEO.html#af10149ff831242b94bf3c92ac233111a

Download: https://gforge.inria.fr/frs/?group_id=5833