# How to obtain the cell-adjacency graph of a mesh?

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"}]

I need this Graph:

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?

• You really need only the connectivity graph, not the vertex positions? I think I can do something about it... – Henrik Schumacher Nov 22 '17 at 10:00
• If you have the initial points, wouldn't it be easier to use DelaunayMesh and MeshPrimitives? As in pi = PositionIndex[MeshPrimitives[DelaunayMesh[pts], 0][[;; , 1]]][[;; , 1]], Graph[Values[pi], MeshPrimitives[DelaunayMesh[pts], 1][[;; , 1]] /. pi, VertexLabels -> "Name"]. – aardvark2012 Nov 22 '17 at 10:02
• @HenrikSchumacher I only need the connectivity graph. Vertex positions are irrelevant. – Szabolcs Nov 22 '17 at 10:08
• @aardvark2012 The VoronoiMesh was just an example. I am looking for a method that works for any mesh. I used a Voronoi mesh as an example because, unlike most other meshes one might encounter in Mathematica, its cells are not triangles (some methods may work only for triangles). – Szabolcs Nov 22 '17 at 10:09
• – hftf Nov 25 '17 at 5:33

I think, I found a general and even faster way, but I haven't tested it for $1$- or $3$-dimensional MeshRegions.

The following function computes the cell-vertex-adjacency matrix A. Two $d$-dimensional cells ($d>0$) are adjacent if they share at least $d$ common points. We can find these pairs by looking for entries $\geq d$ in A.Transpose[A].

ToPack = DeveloperToPackedArray;
Module[{pts, cells, A, lens, n, m, nn},
pts = MeshCoordinates[R];
cells = ToPack[MeshCells[R, d][[All, 1]]];
lens = Length /@ cells;
n = Length[pts];
m = Length[cells];
nn = Total[lens];
A = SparseArray @@ {Automatic, {m, n}, 0, {1, {
ToPack[Join[{0}, Accumulate[lens]]],
ArrayReshape[Flatten[Sort /@ cells], {nn, 1}]
},
ConstantArray[1, nn]}};
SparseArray[
UnitStep[UpperTriangularize[A.Transpose[A], 1] -  d]
]["NonzeroPositions"]
]

A special treatment is necessary for 0-dimensional cells; it's just the edges that we need.

getCellCellAdjacencyList[R_MeshRegion, 0] := ToPack[MeshCells[R, 1][[All, 1]]]

Here are some examples:

SeedRandom[123]
pts = RandomReal[1, {10, 2}];
R = VoronoiMesh[pts];

GraphicsGrid[Table[
{VoronoiMesh[pts, MeshCellLabel -> {d -> "Index"}],
}, {d, 0, 2}], ImageSize -> Large]

And some timings for comparison:

SeedRandom[123]
pts = RandomReal[1, {10000, 2}];
R = VoronoiMesh[pts]; // RepeatedTiming

{0.636, Null}

{0.015, Null}

{0.031, Null}

{0.041, Null}

Edit

It's now rather straight-forward to write methods for the various adjacency matrices, lists, and graphs, even for cells of different dimensions (see below).

Edit 2

As Chip Hurst pointed out, the adjacency matrix of a MeshRegion R for distinct dimensions d1, d2 can be found as pattern SparseArray under R["ConnectivityMatrix"[d1,d2]]. (Its "RowPointers" and "ColumnIndices" must have been computed immediately when the MeshRegion was built.)

Many applications of adjacency matrices, in particular in finite elements, need 1 instead of Pattern as nonzero entries. Even computing vertex rings in a graph by using MatrixPowers of the adjacency matrix is considerably faster with (real) numeric matrices. A remedy could be the function SparseArrayFromPatternArray below. As Chip Hurst has pointed out, we can turn a pattern array into a numerical one with Unitize. I updated my old code to utilize this observation, leading to a tremendous performance boost. Somewhat surprisingly, even the old implementation of CellAdjacencyMatrix[R, 1, 2] tends to be faster than R["ConnectivityMatrix"[1,2]], so that I decided to use the new approach only for the case when either d1 or d2 is equal to 0.

CellAdjacencyMatrix[R_MeshRegion, d_, 0] := If[MeshCellCount[R, d] > 0,
Unitize[R["ConnectivityMatrix"[d, 0]]],
{}
];

CellAdjacencyMatrix[R_MeshRegion, 0, d_] := If[MeshCellCount[R, d] > 0,
Unitize[R["ConnectivityMatrix"[0, d]]],
{}
];

If[MeshCellCount[R, 1] > 0,
With[{B = A.Transpose[A]},
SparseArray[B - DiagonalMatrix[Diagonal[B]]]
]
],
{}
];

If[(MeshCellCount[R, d1] > 0) && (MeshCellCount[R, d2] > 0),
SparseArray[
If[d1 == d2,
UnitStep[B - DiagonalMatrix[Diagonal[B]] - d1],
UnitStep[B - (Min[d1, d2] + 1)]
]
]
],
{}
];

If[(MeshCellCount[R, d1] > 0) && (MeshCellCount[R, d2] > 0),
Module[{i1, i2, data},
data = If[d1 == d2,
];
If[Length[data] > 0,
{i1, i2} = Transpose[data];
Transpose[
{
Transpose[{ConstantArray[d1, {Length[i1]}], i1}],
Transpose[{ConstantArray[d2, {Length[i2]}], i2}]
}
],
{}
]
],
{}
];

Join[MeshCellIndex[R, d1], MeshCellIndex[R, d2]],
VertexLabels -> "Name"
];

Note that CellAdjacencyLists and CellAdjacencyGraph use labels that are compatible with those obtained from MeshCellIndex. Applied to Szabolcs's example MeshRegion, theses graphs look as follows:

GraphicsGrid[
Table[CellAdjacencyGraph[R, d1, d2], {d1, 0, 2}, {d2, 0, 2}],
ImageSize -> Full]

As for comparing the performance of these new implementations to getCellCellAdjacencyList:

{
getCellCellAdjacencyList[R, 0]; // RepeatedTiming // First,
getCellCellAdjacencyList[R, 1]; // RepeatedTiming // First,
getCellCellAdjacencyList[R, 2]; // RepeatedTiming // First
}
{
CellAdjacencyLists[R, 0, 0]; // RepeatedTiming // First,
CellAdjacencyLists[R, 1, 1]; // RepeatedTiming // First,
CellAdjacencyLists[R, 2, 2]; // RepeatedTiming // First
}

{0.015, 0.030, 0.037}

{0.0068, 0.011, 0.0066}

• Thank you for this! I haven't forgotten, and I will accept an answer in due time. – Szabolcs Nov 28 '17 at 17:40
• Would you mind if I incorporated some of this into IGraph/M (with credits of course)? github.com/szhorvat/IGraphM – Szabolcs Nov 28 '17 at 17:43
• @Szabolcs No, I would not. To the contrary: I would feel honored! I would also be thankful for any improvements you might find as I use these things for my own work on surface meshes. – Henrik Schumacher Nov 28 '17 at 20:18
• Just a note that while Transpose[Transpose[edges][[{2, 1}]] is faster than Reverse[edges, {2}], it fails on {}. I think that given that computing the edge list in that function takes two orders of magnitude longer for big meshes, the slowdown due to Reverse is acceptable. – Szabolcs Nov 30 '17 at 12:13
• @Szabolcs Thank you! Good point. I will add a check against empty edge lists. You are right that MeshCells[R, 1] takes pretty long. In an ideal word, the edge list would be cached within the MeshRegion for it is used for quite a number of things. Actually, that's what I do in my own data type for simplicial surfaces. If you subtract the time for getting the edges, then Reverse slows down building the sparse matrix by a factor of two. This is why I prefer a check against emptiness of the edge list. As long as MeshRegions don't cache the edge list, this is a matter of taste, though. – Henrik Schumacher Nov 30 '17 at 12:40

I need three compiled helper functions:

getEdgesFromPolygons = Compile[{{f, _Integer, 1}},
Table[
{
Min[CompileGetElement[f, i], CompileGetElement[f, Mod[i + 1, Length[f], 1]]],
Max[CompileGetElement[f, i], CompileGetElement[f, Mod[i + 1, Length[f], 1]]]
},
{i, 1, Length[f]}
],
RuntimeAttributes -> {Listable},
Parallelization -> True,
CompilationTarget -> "C"
];
takeSortedThread = Compile[{{data, _Integer, 1}, {ran, _Integer, 1}},
Sort[Part[data, ran[[1]] ;; ran[[2]]]],
RuntimeAttributes -> {Listable},
Parallelization -> True,
CompilationTarget -> "C"
];
extractIntegerFromSparseMatrix = Compile[
{{vals, _Integer, 1}, {rp, _Integer, 1}, {ci, _Integer,
1}, {background, _Integer},
{i, _Integer}, {j, _Integer}},
Block[{k},
k = rp[[i]] + 1;
While[k < rp[[i + 1]] + 1 && ci[[k]] != j, ++k];
If[k == rp[[i + 1]] + 1, background, vals[[k]]]
],
RuntimeAttributes -> {Listable},
Parallelization -> True,
CompilationTarget -> "C"
];

The following functions takes a MeshRegion and finds all pairs of neighboring two-dimensional MeshCells. First, it generates a list of all edges (with sorted indices) and creates a lookup table for edges in form of a SparseArray. With the lookup table, we can find the indices of all edges bounding a given polygon, so that we can build the SparseArray edgepolygonadjacencymatrix, whose "AdjacencyLists" is what we are looking for. The method should have linear complexity.

ToPack = DeveloperToPackedArray;
Module[{pts, polygons, edgesfrompolygons, edges, edgelookupcontainer,
pts = MeshCoordinates[R];
polygons = ToPack[MeshCells[R, 2][[All, 1]]];
edgesfrompolygons =  ToPack[Flatten[getEdgesFromPolygons[polygons], 1]];
edges = DeleteDuplicates[edgesfrompolygons];
edgelookupcontainer =
SparseArray[
Rule[Join[edges, Transpose[Transpose[edges][[{2, 1}]]]],
Join[Range[1, Length[edges]], Range[1, Length[edges]]]], {Length[
pts], Length[pts]}];
acc = Join[{0}, Accumulate[ToPack[Length /@ polygons]]];
polyranges = Transpose[{Most[acc] + 1, Rest[acc]}];
edgelookupcontainer["NonzeroValues"],
edgelookupcontainer["RowPointers"],
Flatten@edgelookupcontainer["ColumnIndices"],
edgelookupcontainer["Background"],
edgesfrompolygons[[All, 1]],
edgesfrompolygons[[All, 2]]],
polyranges];
n = Length[edges], m = Length[polygons],
data = ToPack[Flatten[polygonsneighedges]]
},
SparseArray @@ {Automatic, {m, n},
0, {1, {acc, Transpose[{data}]}, ConstantArray[1, Length[data]]}}
];
]

Testing with OP's example:

SeedRandom[123]
pts = RandomReal[1, {10, 2}];
R = VoronoiMesh[pts, MeshCellLabel -> {2 -> "Index"}]
Graph[
VertexLabels -> "Name"
]

Speed test

SeedRandom[123]
pts = RandomReal[1, {10000, 2}];
R = VoronoiMesh[pts,
MeshCellLabel -> {2 -> "Index"}]; // RepeatedTiming

{0.625, Null}

{0.086, Null}

Edit

Slight improvement by replacing Extract with extractIntegerFromSparseMatrix.

Using a few undocumented properties of MeshRegion[] objects, we have the following:

BlockRandom[SeedRandom[123];
vm = VoronoiMesh[RandomReal[1, {10, 2}]]];

Show[vm,
Graph[Range[vm["FaceCount"]],
vm["FaceFaceConnectivity"]]]],
PlotTheme -> "ClassicDiagram",
VertexCoordinates -> Map[Mean, vm["FaceCoordinates"]]]]

I'm not sure why the labels from this version are not consistent with the labels from MeshCellLabel, tho.

• vm["EdgeEdgeConnectivity"] and vm["VertexVertexConnectivity"] might also be of interest. – J. M. is away Nov 22 '17 at 14:12
• Ah, "FaceFaceConnectivity" works for MeshRegions in the plane! I tried these properties several times for two-dimensional MeshRegions in $\mathbb{R}^3$ and it did not work. Seemingly, there is still a lot to be done internally... – Henrik Schumacher Nov 22 '17 at 14:19

Here's another way.

Data from OP:

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

conn = mesh["ConnectivityMatrix"[2, 1]];

Find the cell centroids for visualization purposes:

centers = PropertyValue[{mesh, 2}, MeshCellCentroid];

Show[mesh, g]

Using the same profiling code as Henrik, we have

SeedRandom[123]
pts = RandomReal[1, {10000, 2}];
R = VoronoiMesh[pts]; // RepeatedTiming

RepeatedTiming[
conn = R["ConnectivityMatrix"[2, 1]];
conn . Transpose[conn];
]

{0.632, Null}

{0.042, Null}

{0.012, Null}

• So this was already built in ... it's really a pity that all this stuff is undocumented!! – Szabolcs Feb 23 '18 at 19:23
• I imagine it's partly because this is not a good design for documented functionality. "ConnectivityMatrix"[2, 1] is quite strange looking. – Chip Hurst Feb 23 '18 at 19:26
• Just write it as mesh@"ConnectivityMatrix"[2, 1] and it looks fine. It's exactly the same notation that J/Link uses for method calls. The @ replaces the . that most object-oriented languages would use. – Szabolcs Feb 23 '18 at 19:29
• Oh my, that's so embarassing. Why aren't these things documented? @Szabolcs Of course I don't mind if you make this the accepted answer. – Henrik Schumacher Feb 24 '18 at 12:47
• @ChipHurst Do you know of an efficient way to turn R["ConnectivityMatrix"[2, 1]] into a numerical array with 1 instead of Pattern as nonzero values? That would be very useful for many applications... – Henrik Schumacher Mar 8 '18 at 11:53

Note VoronoiMesh is the dual of the DelaunayMesh, so

Show[mesh,