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
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
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... 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"]. Nov 22 '17 at 10:02
• @HenrikSchumacher I only need the connectivity graph. Vertex positions are irrelevant. 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). 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
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
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. 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 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. 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. Nov 30 '17 at 12:13
• Putting your functions to good use :-) Dec 4 '17 at 21:34

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[] ;; ran[]]],
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
pts = RandomReal[1, {10, 2}];
R = VoronoiMesh[pts, MeshCellLabel -> {2 -> "Index"}]
Graph[
VertexLabels -> "Name"
]  Speed test

SeedRandom
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.

Here's another way.

Data from OP:

SeedRandom
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
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!! Feb 23 '18 at 19:23
• 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. 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. 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... Mar 8 '18 at 11:53
• @ChipHurst Hah! Brilliant! The idea to apply Unitize to anything nonnumerical would never have come to my mind. Sep 7 '18 at 21:43

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

BlockRandom[SeedRandom;
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. 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... Nov 22 '17 at 14:19

In version 12.1, the function MeshConnectivityGraph[] is now built-in. Using the examples in Henrik's answer:

Table[Show[MeshConnectivityGraph[mesh, k, VertexLabels -> "Index"]],
{k, 0, 2}] // GraphicsRow Note VoronoiMesh is the dual of the DelaunayMesh, so

Show[mesh, 