How to obtain the neighboring lists of triangles for a given triangulation [duplicate]

Question

I have asked this question before in my previous post, but I'm still unable to generate the coding. I'm new to Mathematica and I am currently still having some troubles generating the code for finding the neighboring lists of triangles for a given triangulation. Here is what I have so far.

(* Number of data points *)
numpts = 5;
(* Test data set *)
pts = {{0, 0}, {1, 0}, {0, 1}, {1, 1}, {0.3, 0.4}};
dmesh = DelaunayMesh[pts];
(* Extract the triangle list from the Delaunay triangulation *)
tris = MeshCells[dmesh, 2];
numtris = Length[tris];
(* This demonstrates how to access the pts from the tris list *)
Print["Number of triangles numtris = ", numtris];
Do[
  Print["Tri ", i,
   " v1=", pts[[tris[[i, 1, 1]]]],
   " v2=", pts[[tris[[i, 1, 2]]]],
   " v3 = ", pts[[tris[[i, 1, 3]]]]],
  {i, 1, numtris}
  ];
(* Data structure to hold neighbor data *)
nghbrs = Table[{0, 0, 0}, {i, 1, numtris}];

What I'm planning to do for the neighboring list is the following -For each vertex of each triangle, i, define the vertices of edge across from the current vertex, v2 and v3. For each triangle except for the current triangle, look for the triangle that shares the same vertices, v2 and v3, as across from the current vertex. If it is found, load this triangle as neighbor for the current triangle. Otherwise, set neighbor to -1.

I want to put what I'm planning to do for the neighboring list into a do loop but I am having some trouble with that. Can anybody help me by showing me what to input for the do loop? I would appreciate if you can help me with this. Thank you.

Answer 1

I think the following will return the result in desired form without any loop constructs. Strictly speaking this does not answer your question but for performance reasons, you don't want to use a loop.

First, we need a simple helper function that eats a triangle and spits out the opposing oriented edges for each triangle vertex in mathematical positive orientation:

getVertexOpposingEdges = Compile[{{t, _Integer, 1}},
   {{t[[2]], t[[3]]}, {t[[3]], t[[1]]}, {t[[1]], t[[2]]}},
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   RuntimeOptions -> "Speed"
   ];

Now we generate a random delaunay mesh and extract the index triples of the triangles without the Polygon-wrappers. These wrappers actually slow down the processing of the triangles.

pts = RandomReal[{-1, 1}, {25, 2}];
dmesh = DelaunayMesh[pts];
n = MeshCellCount[dmesh, 0]
triangles = MeshCells[dmesh, 2, "Multicells" -> True][[1, 1]];

Next, we let getVertexOpposingEdges run on all triangles; since getVertexOpposingEdges is Listable, it threads automatically over the list triangles.

vertexopposingedges = Flatten[getVertexOpposingEdges[triangles], 1];

Now this is the trick: We generate a sparse matrix triangleLeftToEdge such that triangleLeftToEdge[[i,j]] returns -1 if {i,j} is not an oriented edge of dmesh or if {i,j} is an oriented edge but does not have a triangle to its left. Otherwise triangleLeftToEdge[[i,j]] is the index of the unique triangle to the left of the oriented edge {i,j}.

triangleLeftToEdge = SparseArray[
   Rule[
     vertexopposingedges,
     Join @@ Transpose[{Range[Length[triangles]]}[[ConstantArray[1, 3]]]]
   ],
   {n, n}, -1
   ];

This matrix serves us as a lookup table. We can look up the index of the right neigboring triangle of an oriented edge {i,j} with triangleLeftToEdge[[j,i]]. Let's do it all at once for all oriented edges and pack the result into a Length[triangles] by 3 matrix neighbortriangles; the latter is (hopefully) the desired result:

reversededges = Transpose[Reverse[Transpose[vertexopposingedges]]];
neighbortriangles = Partition[Extract[triangleLeftToEdge, reversededges], 3];

You can check it with

Manipulate[
 HighlightMesh[
  dmesh, {
   Style[{{2, i}}, ColorData[97][4]],
   Labeled[Style[Thread[{0, triangles[[i]]}], ColorData[97][2]], 
    "Index"],
   Labeled[
    Style[Thread[{2, DeleteCases[neighbortriangles[[i]], -1]}], 
     ColorData[97][3]], "Index"]
   }
  ],
 {i, 1, Length[triangles], 1}]

Answer 2

There is a simple way to do this if you use ElementMesh:

Needs["NDSolve`FEM`"]
mesh = ToElementMesh[pts];
mesh["ElementConnectivity"]

{{{2, 3, 0}, {8, 5, 1}, {1, 9, 0}, {18, 7, 5}, {9, 2, 4}, {7, 22, 
   17}, {6, 11, 4}, {18, 2, 0}, {5, 11, 3}, {12, 16, 0}, {9, 7, 
   12}, {11, 14, 10}, {28, 38, 16}, {12, 17, 15}, {14, 29, 28}, {0, 
   10, 13}, {6, 31, 14}, {8, 39, 4}, {21, 39, 20}, {0, 26, 19}, {23, 
   27, 19}, {6, 39, 27}, {21, 26, 25}, {0, 25, 26}, {24, 37, 23}, {20,
    24, 23}, {31, 22, 21}, {13, 15, 30}, {36, 15, 33}, {36, 38, 
   28}, {33, 17, 27}, {0, 35, 37}, {37, 29, 31}, {0, 38, 35}, {36, 32,
    34}, {30, 29, 35}, {32, 33, 25}, {34, 13, 30}, {18, 19, 22}}}

The first element is connected to elements 2 and 3. Zeros indicate a boundary. Negative numbers (not present here) indicate an internal boundary.

This is documented on the ElementMesh ref page in the scope section. There is also a companion method called "ContinuousElementConnectivity" which returns a different list structure for mixed type elements.