Getting adjacent cells map for an unstructured polyhedral mesh

Question

I am doing a little project on solving the heat equation using finite-volume method on a solid cube, I converted the polyhedral mesh of the cube to an OpenFOAM mesh.

I have a Python code where I parse the points, faces, owner, neighbour and boundary files of the OF case mesh, I've managed to create a hash list (let's call it CellFacesMap) to map the cell number to the list of its owned faces and faces neighbours to it. Now, I am trying to generate a hash list (let's call it CellAdjacentsMap that maps each cell with all its adjacent cells, in order to apply the discretized heat equation to each cell (and also for generating the volume weighting factors).

The only way I could think of now to do so is to visit each cell in CellFacesMap and for the owned faces I get the cells neighbor to this face (adjacent cells) and for faces neighbor the same cell I get the cells that owns theses faces (also adjacent cells), meaning I visit CellFacesMap over and over for each cell to check its faces.

But this method is costly (time wise, it took 3.5 minutes to generate a map for 4000+ cell mesh). So is there any faster way to generate CellAdjacentsMap?

Answer 1

I'm assuming you are starting with a list of cell definitions, say, as a list of vertices defining the cell and a "type" defining the topology of each cell. As part of the topology definition for each cell, you can easily get the faces for that cell defined, say, as a list of vertices.

The first step in creating the CellAdjacentsMap that you need is to first create a map that I'll call FaceCellMap that lists the cells (either one or two) that each face belongs to. You do this by iterating over all cells, adding each cell's faces to the map, one by one.

There are some implementation details in forming this FaceCellMap. The first is what data structure do you use? If you have a general HashMap capability, the key can be face definition and the value can be a length-two array of cell IDs (one of the entries may be zero). The second issue is how do you compare faces? You can use a list of vertices but you have implement the comparison function in such a way that it is insensitive to the "starting" vertex and whether you traverse the vertex list in clockwise or counterclockwise order.

But after you have created FaceCellMap, creating CellAdjacentsMap is straightforward. You iterate through all the cells. For each cell, you iterate through all its faces. For each face, you check the FaceCellMap for that face to see if it has an adjacent cell. If so you add it to the list of adjacencies for the current cell.

Answer 2

Using a hash map adds a log(n) complexity to all accesses (then traversing the whole mesh will cost n log(n) in general), so clearly it is not the best solution.

Now your question is how you can efficiently store and traverse adjacency information, i.e. for each cell the list of adjacent cells. First thing you can do for mapping a cell to its list of neighbors: you can use an array instead of a hash, this will remove the log(n) cost.

Then, using an array, this means you need to store an array of lists of different lengths (if your mesh has arbitrary polyhedra), which is in general not very efficient: each list has its own memory overhead to be added to the memory used by the elements it contains. A more efficient way of storing this type of information (i.e. an array of lists of different lengths) is the "compressed row storage" (called CRS, CSR or Yale format) representation [1], used to store sparse matrices.

Edit If adjacency information is not present in your input, then you will need to recreate it. Instead of using a hash or associative map as suggested in the other answer, an alternative is to store all facets in an array then sort this array. After sorting, the pairs of facets that correspond to adjacent cells are contiguous in the array. I observed that it is significantly more efficient (at least in C++, I did not test that in Python). See also my answer to this question [2]. To adapt my answer to your polyhedral grid, facet records will have a facet index and as a key, three vertex indices (the lowest vertex id in the facet and the two ones connected to it with an edge, in increasing order). By sorting facet records with the lexicographic order, you will retrieve the pair of adjacent facets contiguously in the array.

The algorithm:

For each cell c
   For each facet f of cell c
       Find v1 the smallest vertex index in facet f
       Find v2 and v3, the vertices before and after v1 in facet f
       If v3 < v2, swap v2 and v3
       Append a new record [v1, v2, v3, c] to the array

Sort the record array with the lexicographic order on (v1,v2,v3)

In the record array, find (by scanning the array) the contiguous 
pairs of record (v1,v2,v3,c1);(v1,v2,v3,c2) with the same (v1,v2,v3). 
For each such pair, c1 and c2 are adjacent.

If input data has facet ids, then records are simply (f,c) couples, and they are sorted by f (instead of (v1,v2,v3)).

[1] https://en.wikipedia.org/wiki/Sparse_matrix

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