How to efficiently get mesh cell/face connectivity?

Question

If I have 3 lists: points, faces and cells that represents a mesh, where:

• points is a list of x,y,z coordinates, for example

   [-0.05, -0.05,  0.  ],
   [-0.05,  0.05,  0.  ],
   [ 0.05, -0.05,  0.  ],
   ...
  

• faces is a list of list of points, where each row is a list of point indices that represents a face

   [ 1,  4,  8,  7],
   [ 7,  8,  6,  3],
   [ 4,  0,  5,  8],
   ...
  

• cells is a list of list of points, where each row is a list of point indices that represents a cell

   [ 1,  7,  8,  4, 11, 16, 17, 10],
   [ 7,  3,  6,  8, 16, 15, 14, 17],
   [ 4,  8,  5,  0, 10, 17, 12,  9],
   ...
  

For each face, I need to get the two cells ids that are sharing the face (in case of boundary face, there is only one cell that owns that face).

For a small mesh, I can iterate over cells, and for each cell iterate over all faces and check if a face points are in current cell (two for loops and a face-cells map), which is of course awfully slow.

Is there a smarter way to efficiently get face connectivity of my mesh?

Answer 1

I'm answering this question because it's basic meshing knowledge that should be available and I believe neither answer is satisfactory:

• @Francler In an unstructured mesh, there is no theoretical upper bound on the maximum number of elements that contain a given vertex. You could allocate for a practical upper bound (say 100), but do you have the memory for it? What if the mesh has 10M points? 100M points? What if it's a silly disk with a point in the middle starred against the boundary (hundreds of triangles for one point)? And the computational cost is O(V_MN) where V_M is the maximum valence of a point (in the sense of elements containing it).
• @MPIchael's answer only applies to structured grids

The solution is to use hash tables. A very simple hash table is comprised of two arrays:

• the head of size M
• the list of size (4+2+1)xN

and a hash function : $\mathbb N^4$ -> [1,N].

N must be large enough to store all edges in the mesh. You can use Euler's relations to estimate a value. M should be proportional to N and offers a trade-off between memory and speed. If M is larger than N, your algorithm will be fastest but take more memory. M can be much smaller than N, to reduce memory footprint, but it'll slow down the method. The hash function should map numbers as randomly as possible. A common construction is to take any "chaotic" function and apply modulo N. If you ever need to increase speed of this method, this is where to look for improvements.

If you don't need to store the two cells, but only need to know the two cells at a given point (to construct neighbour arrays), you can have list of size (4+1+1)xN.

You'll find plenty of literature on hash tables, but the main idea is:

• You want to know if the face (i,j,k,l) is in the hash table (this is called the key). You compute the hash I = hashfunc(i,j,k,l) and go to head(I).
  • If head(I) = 0, there exists no element in the hash-table that has the same hash as (i,j,k,l).
  • Else, if J = head(I) > 0, there exists at least one element in the hash-table that has the same hash as the key (i,j,k,l). The first such element is found at list(J,:). You can compare list(J,1:4) to (i,j,k,l) to see if they're the same face. If not, K = list(J,7) is the index of the next element in the list with the same hash (if it exists, 0 otherwise).

It is most efficient to sort keys (i,j,k,l) on insertion and to use quick comparison, but you can also implement any-order comparison ; this is because two elements see a given face in the opposite order, and starting at a different first vertex.

Above, I described lookup, but insertion is very similar. On insertion, you store the key (i,j,k,l) but also the element that provided this key. It is stored at either list(J,5) if list(J,5) = 0, or at list(J,6).

The final algorithm is:

for ielem = 1,nelem
    for i = 1,6
        i1,i2,i3,i4 = nodes of i-th face of ielem
        if (i1,i2,i3,i4) is in hash-table at I
            // list(I,5) is my neighbour
            // I can be put at list(I,6)
        else
            insert pair ((i1,i2,i3,i4), ielem) in hash-table
        endif
    endfor
endfor

Its average complexity is O(Nelem Ncol) where Ncol is the average number of collisions in the hash-table and is dictated by M/N. This can be as low as 1.

Another name for hash tables is "maps", so that may be a good keyword to search documentation of your preferred language with for an implementation. Otherwise, it's not very difficult to implement. In C++, you have std::unordered_maps.

Answer 2

Assuming that you are working with cubic cells, you can create a vector called point2cell, such that point2cell[iPoint] gives you the indeces of the cells sharing the point iPoint. This can be done by looping over the cells vector. Take a look at the following pseudo-code:

for iCell = 1:nCell
  for iPointLoc = 1:8
     iPoint = cells[iCell][iPointLoc]
     point2cell[iPoint].insert( iCell )

Next, you can loop over the faces, take one point, take all the cells sharing that point, and looking at the matching between a face and a cell:

for iFace = 1:nFace
  for iPointLoc = 1:4
    iPoint = faces[iFace][iPointLoc]
    for iCellLoc = 1:8
       iCell = point2Cell[iPoint][iCellLoc]
       if( cells[iCell] matches faces[iFace])
         face2cell[iFace].insert( iCell ) 

The first loop is linear in the number of cells while the second is linear in the number of faces.

Answer 3

I think you can achive order N performance if you sort the list of points, list of faces and list of cells in the same manner, i.e., if you sort row-column-layer by midpoints of your three entity types. If that is the case, you can approximate the offset that you will need to find your associated elements. Then you can iterate in your neighbourhood of the approximated offset to find what you need.

Example to clearify: Lets say you have a 2D grid with 9 quadrilaterals. You sort the 16 Points row-column wise, the faces also row-column-wise by their centers and finally the 9 cells row-column wise by their midpoints. If you then want to find the two neighbour cells for a face, you can calculate the offset efficiently. This should take you from Order $O(N^2)$ to $O(N)$.

This of course assumes structured grids.