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)
insert pair ((i1,i2,i3,i4), ielem) in hash-table
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.