2

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: Triangulation of base rectangle Triangulation of base annulus

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. enter image description here

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!

ap21
  • 2,372
  • 3
  • 16
  • 32

1 Answers1

2

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).

enter image description here enter image description here

hpaulj
  • 221,503
  • 14
  • 230
  • 353
  • Great! Would you mind elaborating on why it was filling the interior with triangles? – ap21 Jan 22 '18 at 18:09
  • I don't think I can add to what http://www.qhull.org/ says. The key is the `convex hull`, which in this case are the points corresponding to the largest `v`, the outer ring. It fills this whole hull. – hpaulj Jan 22 '18 at 18:17
  • Understood. A different thing: setting `u = u[:-1]` does not work for me? Where exactly are you implementing this? – ap21 Jan 22 '18 at 23:55
  • I do that right at the start, so `u` doesn't run all the way to `2*pi`. You must have done something like that to show the gap in the 3rd figure. – hpaulj Jan 23 '18 at 00:01
  • We can double check, but I don't think `qhull` expects the points to be on a grid or even ordered. The order may affect the individual triangles, but the convex hull should be the same. – hpaulj Jan 23 '18 at 00:52
  • I agree with that above point. Last thing -- is it possible to apply that mask permanently to the full `Delaunay triangulation` object, so that I don't need to carry it about through the program? Or even just to `Delaunay.vertices`? – ap21 Jan 23 '18 at 03:09
  • Would you happen to have a response to my last question? Can I apply the mask permanently? – ap21 Jan 25 '18 at 22:59