25
$\begingroup$

In addition to the accepted answer, see also the answer by Chip Hurst. 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

I need this Graph:

enter image description here

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?

Improve this question
$\endgroup$
6
  • $\begingroup$ You really need only the connectivity graph, not the vertex positions? I think I can do something about it... $\endgroup$ Commented Nov 22, 2017 at 10:00
  • $\begingroup$ 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"]. $\endgroup$ Commented Nov 22, 2017 at 10:02
  • $\begingroup$ @HenrikSchumacher I only need the connectivity graph. Vertex positions are irrelevant. $\endgroup$ Commented Nov 22, 2017 at 10:08
  • $\begingroup$ @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). $\endgroup$ Commented Nov 22, 2017 at 10:09
  • $\begingroup$ Related: mathematica.stackexchange.com/questions/90109/… mathematica.stackexchange.com/questions/52056/… $\endgroup$ Commented Nov 25, 2017 at 5:33

6 Answers 6

Reset to default
20
$\begingroup$

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 = Developer`ToPackedArray; ClearAll[getCellCellAdjacencyList]; getCellCellAdjacencyList[R_MeshRegion, d_] := 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"}], Graph[getCellCellAdjacencyList[R, d], VertexLabels -> "Name"] }, {d, 0, 2}], ImageSize -> Large]

enter image description here

And some timings for comparison:

SeedRandom[123] pts = RandomReal[1, {10000, 2}]; R = VoronoiMesh[pts]; // RepeatedTiming getCellCellAdjacencyList[R, 0]; // RepeatedTiming getCellCellAdjacencyList[R, 1]; // RepeatedTiming getCellCellAdjacencyList[R, 2]; // 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]]], {} ]; CellAdjacencyMatrix[R_MeshRegion, 0, 0] := If[MeshCellCount[R, 1] > 0, With[{A = CellAdjacencyMatrix[R, 0, 1]}, With[{B = A.Transpose[A]}, SparseArray[B - DiagonalMatrix[Diagonal[B]]] ] ], {} ]; CellAdjacencyMatrix[R_MeshRegion, d1_, d2_] := If[(MeshCellCount[R, d1] > 0) && (MeshCellCount[R, d2] > 0), With[{B = CellAdjacencyMatrix[R, d1, 0].CellAdjacencyMatrix[R, 0, d2]}, SparseArray[ If[d1 == d2, UnitStep[B - DiagonalMatrix[Diagonal[B]] - d1], UnitStep[B - (Min[d1, d2] + 1)] ] ] ], {} ]; CellAdjacencyLists[R_MeshRegion, d1_, d2_] := If[(MeshCellCount[R, d1] > 0) && (MeshCellCount[R, d2] > 0), Module[{i1, i2, data}, data = If[d1 == d2, UpperTriangularize[CellAdjacencyMatrix[R, d1, d2], 1]["NonzeroPositions"], CellAdjacencyMatrix[R, d1, d2]["NonzeroPositions"] ]; If[Length[data] > 0, {i1, i2} = Transpose[data]; Transpose[ { Transpose[{ConstantArray[d1, {Length[i1]}], i1}], Transpose[{ConstantArray[d2, {Length[i2]}], i2}] } ], {} ] ], {} ]; CellAdjacencyGraph[R_MeshRegion, d1_, d2_] := Graph[ Join[MeshCellIndex[R, d1], MeshCellIndex[R, d2]], UndirectedEdge @@@ CellAdjacencyLists[R, d1, 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]

enter image description here

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}

Improve this answer
$\endgroup$
7
  • $\begingroup$ Thank you for this! I haven't forgotten, and I will accept an answer in due time. $\endgroup$ Commented Nov 28, 2017 at 17:40
  • $\begingroup$ Would you mind if I incorporated some of this into IGraph/M (with credits of course)? github.com/szhorvat/IGraphM $\endgroup$ Commented Nov 28, 2017 at 17:43
  • $\begingroup$ @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. $\endgroup$ Commented Nov 28, 2017 at 20:18
  • 1
    $\begingroup$ 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. $\endgroup$ Commented Nov 30, 2017 at 12:13
  • 1
    $\begingroup$ Putting your functions to good use :-) $\endgroup$ Commented Dec 4, 2017 at 21:34
13
$\begingroup$

I need three compiled helper functions:

getEdgesFromPolygons = Compile[{{f, _Integer, 1}}, Table[ { Min[Compile`GetElement[f, i], Compile`GetElement[f, Mod[i + 1, Length[f], 1]]], Max[Compile`GetElement[f, i], Compile`GetElement[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 = Developer`ToPackedArray; getPolygonPolygonAdjacencyList[R_MeshRegion] := Module[{pts, polygons, edgesfrompolygons, edges, edgelookupcontainer, polyranges, polygonsneighedges, edgepolygonadjacencymatrix, acc}, 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]}]; polygonsneighedges = takeSortedThread[extractIntegerFromSparseMatrix[ edgelookupcontainer["NonzeroValues"], edgelookupcontainer["RowPointers"], Flatten@edgelookupcontainer["ColumnIndices"], edgelookupcontainer["Background"], edgesfrompolygons[[All, 1]], edgesfrompolygons[[All, 2]]], polyranges]; edgepolygonadjacencymatrix = Transpose@With[{ n = Length[edges], m = Length[polygons], data = ToPack[Flatten[polygonsneighedges]] }, SparseArray @@ {Automatic, {m, n}, 0, {1, {acc, Transpose[{data}]}, ConstantArray[1, Length[data]]}} ]; Select[(edgepolygonadjacencymatrix["AdjacencyLists"]), Length[#] == 2 &] ]

Testing with OP's example:

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

enter image description here

enter image description here

Speed test

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

{0.625, Null}

{0.086, Null}

Edit

Slight improvement by merging Sort into takeThread (takeThread replaced by takeSortedThread).

Slight improvement by replacing Extract with extractIntegerFromSparseMatrix.

Improve this answer
$\endgroup$
11
$\begingroup$

Here's another way.

Data from OP:

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

Get the adjacency matrix:

conn = mesh["ConnectivityMatrix"[2, 1]]; adj = conn.Transpose[conn];

Find the cell centroids for visualization purposes:

centers = PropertyValue[{mesh, 2}, MeshCellCentroid]; g = AdjacencyGraph[adj, PlotTheme -> "Scientific", VertexCoordinates -> centers]; Show[mesh, g]

enter image description here

Using the same profiling code as Henrik, we have

SeedRandom[123] pts = RandomReal[1, {10000, 2}]; R = VoronoiMesh[pts]; // RepeatedTiming getCellCellAdjacencyList[R, 2]; // RepeatedTiming RepeatedTiming[ conn = R["ConnectivityMatrix"[2, 1]]; conn . Transpose[conn]; ]

{0.632, Null}

{0.042, Null}

{0.012, Null}

Improve this answer
$\endgroup$
14
  • 2
    $\begingroup$ So this was already built in ... it's really a pity that all this stuff is undocumented!! $\endgroup$ Commented Feb 23, 2018 at 19:23
  • 1
    $\begingroup$ 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. $\endgroup$ Commented Feb 23, 2018 at 19:29
  • 1
    $\begingroup$ 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. $\endgroup$ Commented Feb 24, 2018 at 12:47
  • 1
    $\begingroup$ @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... $\endgroup$ Commented Mar 8, 2018 at 11:53
  • 1
    $\begingroup$ @ChipHurst Hah! Brilliant! The idea to apply Unitize to anything nonnumerical would never have come to my mind. $\endgroup$ Commented Sep 7, 2018 at 21:43
10
$\begingroup$

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"]], Union[Sort /@ Flatten[MapIndexed[Thread[UndirectedEdge[#2[[1]], #1]] &, vm["FaceFaceConnectivity"]]]], PlotTheme -> "ClassicDiagram", VertexCoordinates -> Map[Mean, vm["FaceCoordinates"]]]]

mesh and its face-face connectivity graph

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

Improve this answer
$\endgroup$
2
  • $\begingroup$ vm["EdgeEdgeConnectivity"] and vm["VertexVertexConnectivity"] might also be of interest. $\endgroup$ Commented Nov 22, 2017 at 14:12
  • $\begingroup$ 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... $\endgroup$ Commented Nov 22, 2017 at 14:19
4
$\begingroup$

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

mesh connectivities

Improve this answer
$\endgroup$
1
$\begingroup$

Note VoronoiMesh is the dual of the DelaunayMesh, so

Show[mesh, AdjacencyGraph[DelaunayMesh[pts]["AdjacencyMatrix"], VertexCoordinates -> MeshCoordinates[DelaunayMesh[pts]], EdgeStyle -> Red]]

mesh and graph

But I don't know how to get the right vertex label..

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.