In addition to the accepted answer, [see also the answer by Chip Hurst.](https://mathematica.stackexchange.com/a/166491/12) This functionality is built in, but not documented.

---

Given an arbitrary mesh region, how can I efficiently obtain the graph describing the adjacency structure of mesh cells?

For example, given the following mesh,

```
SeedRandom[123]
pts = RandomReal[1, {10, 2}];
mesh = VoronoiMesh[pts, MeshCellLabel -> {2 -> "Index"}]

```

[![enter image description here](https://i.stack.imgur.com/Im87X.png)](https://i.stack.imgur.com/Im87X.png)

I need this `Graph`:

[![enter image description here](https://i.stack.imgur.com/SadF7.png)](https://i.stack.imgur.com/SadF7.png)

It tells me that e.g. cells 4 and 5 are neighbours.

This example was for 2D cells, but the problem can be stated generally for cells of any dimension:

* two edges (1D cells) are adjacent if they have a common point
* two faces (2D cells) are adjacent if they have a common edge
* two 3D cells are adjacent if they have a common face
* ...

I am looking both for specialized methods to obtain the face-adjacency graph of a 2D or 3D mesh, and general methods to obtain the $k$-dimensional cell adjacency graph of a $d > k$ dimensional mesh.

---

Here's a naïve method for face-adjacency. It is too slow for practical use due to its $O(n^2)$ complexity.

```
SimpleGraph[
 RelationGraph[
  Length@Intersection[First@MeshCells[mesh, #1], First@MeshCells[mesh, #2]] >= 2 &,
  MeshCellIndex[mesh, 2]
  ],
 VertexLabels -> "Name"
 ]

```

This method exploits the fact that if two polygons (faces) are adjacent, then they will have (at least) two points in common. For example, `Polygon[{14, 8, 7, 11}]` and `Polygon[{11, 7, 2, 4, 13}]` have points `7` and `11` in common. It generalizes to higher dimensions as well: two 3D cells are adjacent if they have at least 3 points in common.

However, it is quite slow because `RelationGraph` will test each pair of cells for adjacency.

```
SeedRandom[123]
pts = RandomReal[1, {500, 2}];
mesh = VoronoiMesh[pts];

RelationGraph[
   Length@Intersection[First@MeshCells[mesh, #1], 
       First@MeshCells[mesh, #2]] >= 2 &,
   MeshCellIndex[mesh, 2]
   ]; // AbsoluteTiming

(* {2.36978, Null} *)

```

While this method can be sped up by a constant factor with minor tweaks, this won't fix the core problem: quadratic complexity.

```
cells = MeshCells[mesh, 2][[All, 1]];        
RelationGraph[Length@Intersection[#1, #2] >= 2 &, cells]; // AbsoluteTiming
(* {0.815857, Null} *)

```

1 second for only 500 cells is still too slow. Can we do significantly better?