Getting a proper Delaunay triangulation of an annulus (using python)
I am trying to triangulate an annulus using the scipy.spatial.Delaunay() function, but cannot get the desired result. Here is my code:
from scipy.spatial import Delaunay
NTheta = 26
NR = 8
a0 = 1.0
#define base rectangle (r,theta) = (u,v)
u=np.linspace(0, 2*np.pi, NTheta)
v=np.linspace(1*a0, 3*a0, NR)
u,v=np.meshgrid(u,v)
u=u.flatten()
v=v.flatten()
#evaluate the parameterization at the flattened u and v
x=v*np.cos(u)
y=v*np.sin(u)
#define 2D points, as input data for the Delaunay triangulation of U
points2D=np.vstack([u,v]).T
xy0 = np.vstack([x,y]).T
Tri1 = Delaunay(points2D) #triangulate the rectangle U
Tri2 = Delaunay(xy0) #triangulate the annulus
#plt.scatter(x, y)
plt.triplot(x, y, Tri1.simplices, linewidth=0.5)
plt.show()
plt.triplot(x, y, Tri2.simplices, linewidth=0.5)
plt.show()
I get the following:
The triangulation of the annulus itself clearly gives unwanted triangles. The triangulation of the base rectangle seems to give the proper result, until you realise that the annulus is not actually closed, by stretching the annulus (i.e., moving its nodes) a bit.
So, my question is, how do I get the proper triangulation that accounts for the non-trivial topology? Can I remove simplices from the triangulation of the annulus -- for example, based on the length of the bonds -- or somehow stitch the two ends of the base rectangle together? Is there a simple way of doing this?
Answer:
I accepted the answer below but it does not completely solve the question as asked. I still don't know how to tile a periodic surface using scipy.Delaunay (i.e., the qhull routine). However, using a mask as defined below, one can create a new list of triangle simplices, and that should serve for many purposes. However, one cannot use this list with the other methods defined in the scipy.Delaunay class. So, be careful!
qhull works with the convex hull. So it can't work directly with that concave interior. In fig2 it is filling the interior with triangles. That may be more obvious if we add a (0,0) point to xy0.
last_pt = xy0.shape[0]
xy1 = np.vstack((xy0,(0,0))) # add ctr point
Tri3 = Delaunay(xy1)
print(Tri3.points.shape, Tri3.simplices.shape)
plt.triplot(Tri3.points[:,0], Tri3.points[:,1], Tri3.simplices, linewidth=0.5)
plt.show()
Remove the simplices that contain that center point:
mask = ~(Tri3.simplices==last_pt).any(axis=1)
plt.triplot(Tri3.points[:,0], Tri3.points[:,1], Tri3.simplices[mask,:], linewidth=0.5)
plt.show()
To stitch the two ends together, removing a value from u seems to work:
u = u[:-1]
In a FEM model you might leave the center elements in place, but give them the appropriate 'neutral' properties (insulating or whatever works).
