Search code examples
pythonmatplotlibscipytriangulationdelaunay

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


Solution

  • 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