HJA24
HJA24

Reputation: 422

Create triangular mesh of pentagon or higher

I have the coordinates of the corners of a pentagon

import math
import numpy as np

n_angles = 5
r = 1

angles = np.linspace(0, 2 * math.pi, n_angles, endpoint=False)
x = (r * np.cos(angles))
y = (r * np.sin(angles))
corners = np.array(list(zip(x, y)))

From these points I want to create a triangular mesh using scipy.spatial.Delaunay.

import matplotlib.pyplot as plt
from scipy.spatial import Delaunay

tri = Delaunay(points=corners)

plt.triplot(corners[:,0], corners[:,1], tri.simplices)
plt.plot(corners[:,0], corners[:,1], 'o')
for corner in range(len(corners)):
    plt.annotate(text=f'{corner + 1}', xy=(corners[:,0][corner] + 0.05, corners[:,1][corner]))

plt.axis('equal')    
plt.show()

The result looks as follows:

enter image description here

Why are the simplices (1, 3), (1, 4), (2, 4) missing and how can I add them?\

I think I need to specify incremental=True in Delaunay() and use scipy.spatial.Delaunay.add_points but this results in an error:

QhullError: QH6239 Qhull precision error: initial Delaunay input sites are cocircular or cospherical.  Use option 'Qz' for the Delaunay triangulation or Voronoi diagram of cocircular/cospherical points; it adds a point "at infinity".  Alternatively use option 'QJ' to joggle the input

Please advice

Upvotes: 2

Views: 361

Answers (1)

mozway
mozway

Reputation: 260835

scipy.spatial.Delaunay computes a triangular tessellation, meaning non-overlapping triangles.

You rather want all combinations of vertices.

Use the combination of "corners" indices as a source of triangles for triplot.

from itertools import combinations

plt.triplot(corners[:,0], corners[:,1], list(combinations((range(len(corners))), r=3)))

Full code:

import matplotlib.pyplot as plt
from itertools import combinations

plt.triplot(corners[:,0], corners[:,1], list(combinations((range(len(corners))), r=3)))

plt.plot(corners[:,0], corners[:,1], 'o')
for corner in range(len(corners)):
    plt.annotate(text=f'{corner + 1}', xy=(corners[:,0][corner] + 0.05, corners[:,1][corner]))

plt.axis('equal')    
plt.show()

Alternatively, skip the triplot and use a LineCollection from a combination of pairs of points:

import matplotlib.pyplot as plt
from matplotlib import collections  as mc

ax = plt.subplot()
plt.plot(corners[:,0], corners[:,1], 'o')
lines = mc.LineCollection(combinations(corners, r=2))
ax.add_collection(lines)

for corner in range(len(corners)):
    plt.annotate(text=f'{corner + 1}', xy=(corners[:,0][corner] + 0.05, corners[:,1][corner]))

plt.axis('equal')    
plt.show()

Output:

enter image description here

Upvotes: 3

Related Questions