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

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.

Henrik Schumacher · Answer 1 · 2018-04-25T13:08:24.510

4

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

edited Apr 25 '18 at 13:08

answered Apr 25 '18 at 08:50

Henrik Schumacher

106,770
7
179
309

for pts and M, can I substitute them with what I have for pts and dmesh? Will that change the code? – J.Kwan Apr 25 '18 at 09:03
I replace M by dmesh. Should work fine for any triangle mesh dmesh. – Henrik Schumacher Apr 25 '18 at 09:07
I can also replace pts with the pts that I had written. Is that correct? – J.Kwan Apr 25 '18 at 09:11
Yes, you can... – Henrik Schumacher Apr 25 '18 at 09:13
Ok. I will try that and if I have any more questions, I will input my questions through comment. – J.Kwan Apr 25 '18 at 09:14
I can't get the getVertexOpposingEdges to work. It says "CCompilerDriver`CreateLibrary: A C compiler cannot be found on your system. Please consult the documentation to learn how to set up suitable compilers." It also says "Compile: A library could not be generated from the compiled function." Any idea why I get these errors? – J.Kwan Apr 25 '18 at 09:21
Ahm you have no C-compiler installed. No problem. Just change CompilationTarget -> "C" to CompilationTarget -> "WVM". – Henrik Schumacher Apr 25 '18 at 09:22
Will that solve the problem? If so, I will try that tomorrow. – J.Kwan Apr 25 '18 at 09:53
Yes, that's supposed to solve the problem. That's why I posted it. – Henrik Schumacher Apr 25 '18 at 09:56
2

MeshCells[dmesh, 2, "Multicells" -> True] I think you will like that. – user21 Apr 25 '18 at 13:04
@user21 Do you have an idea how long I have been looking for something like that? (Of course you do...) It would be great if that would be added to the documentation (adjusting the syntax highlighting would also be fine...) – Henrik Schumacher Apr 25 '18 at 13:07
@user21 And of course: Many thanks! That's indeed very helpful! – Henrik Schumacher Apr 25 '18 at 13:15
@HenrikSchumacher, was not my decision not to document that.... But I'd send a suggestion to WRI... – user21 Apr 25 '18 at 13:31

score 3 · Answer 2 · answered Apr 30 '18 at 05:03

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.

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

2 Answers2