Performing queries on mesh structure in a clean way

Question

I have the following problem, which can be stated geometrically as follows:

List all pairs of vertex coordinates in a mesh region of measure 2 (a surface consisting of triangles) where these points are on two triangles sharing an edge, and where these resulting vertex points are not on the shared edge.

Best I came around with is the following:

With[{mesh = DiscretizeRegion[Rectangle[], MaxCellMeasure -> 1/9]}, 
 Function[line, 
     Select[MeshPrimitives[mesh, 2], 
       MemberQ[Sort /@ Partition[First@#, 2, 1, 1], line] &] // 
      If[Length@# == 2, Complement[Flatten[#[[All, 1]], 1], line], 
        Nothing] &]@*Sort@*First /@ MeshPrimitives[mesh, 1] // 
  Show[mesh, Graphics[{Red, Line@#}]] &]

Frankly, I'm not at all pleased of all the munging and unintuitiveness on my code, which looks for all lines on the mesh, then finds all pairs of triangles sharing it as an edge (note Sort unintuitiveness there!), and then excludes the vertices on the shared side before spitting out the result. Surely there would be a more intuitive way to do this?

Another way to formulate this problem would be to consider the mesh to be a (well-behaved, planar) graph, and query for all pairs of neighboring vertices on the dual graph, and to extract desired "real" triangle vertices from there. Is making such a query even possible in an efficient manner with Mma built-in graph functionality? (I doubt it.)

Answer 1

With the Mathematica package "IGraphM`" by Szabolcs, you can do the following.

Needs["IGraphM`"];
mesh = DiscretizeRegion[Rectangle[], MaxCellMeasure -> 1/900];
pairs = Complement[Union @@ #, Intersection @@ #] & /@ Partition[
    MeshCells[mesh, 2, "Multicells" -> True][[1, 1,
     Flatten[IGMeshCellAdjacencyMatrix[mesh, 2, 2]["NonzeroPositions"]]
     ]],
    2
    ];
Show[
 mesh,
 Graphics[{Red, 
   Line[Partition[MeshCoordinates[mesh][[Flatten[pairs]]], 2]]}]
 ]

Basically, IGMeshCellAdjacencyMatrix[mesh, 2, 2] computes the adjacency matrix of faces, (this matrix has $A_{ij} = 1$ if $i \neq j$ and if triangle $i$ and $j$ share an edge and $A_{ij} = 0$ otherwise. Try also the routine IGMeshCellAdjacencyGraph. For more explanations on IGMeshCellAdjacencyMatrix see also this post. The implementation of CellAdjacencyMatrix should also be a bit more optimized than those of "IGraphM`".

There are certainly faster ways. For example, mapping the function Complement[Union @@ #, Intersection @@ #] & should be a severe performance bottleneck. But that's something to start from.

A quicker approach is

A = IGMeshCellAdjacencyMatrix[mesh, 0, 2].IGMeshCellAdjacencyMatrix[mesh, 2, 1];
pairs = SparseArray[Clip[
 UpperTriangularize[A.Transpose[A],1], 
 {1, 2}, {0, 0}
 ]]["NonzeroPositions"];

For some MeshRegions (I could not figure out which one but those planar meshes that are generated by DelaunayMesh and VoronoiMesh are among them) have a property called "WingData". For each edge, it's a list of eight integers in roughly this order:

1. index of triangle to the right
2. index of triangle to the left
3. index of first vertex of the edge
4. index of second vertex of the edge
5. index of edge emanating from first vertex and belonging to right triangle?
6. index of edge emanating from first vertex and belonging to left triangle?
7. index of edge emanating from second vertex and belonging to right triangle?
8. index of edge emanating from second vertex and belonging to left triangle?

So, another way to get the pairs of triangles is

R = DelaunayMesh[MeshCoordinates[mesh]];
R["WingData"][[All, 1 ;; 2]]