title: Python: Delaunay triangulation of a cylinder (or any periodic surface)
body: <p>A year ago, I asked about properly triangulating a periodic shape on a plane: the annulus (<a href="https://stackoverflow.com/questions/48372595/getting-a-proper-delaunay-triangulation-of-an-annulus-using-python">Getting a proper Delaunay triangulation of an annulus (using python)</a>).</p>
  
  <p>I now want to expand this to triangulating a cylinder (or in general, any periodic surface). I try a straightforward extension of the 2D code:</p>
  
  <pre><code>from scipy.spatial import Delaunay
  NZ = 14
  NTheta = 14
  
  R = 1 #radius of cylinder 
  L = 3 #length of cylinder 
  
  #define base rectangle (u,v)
  u=np.linspace(0, 2*np.pi, NTheta) #periodic direction
  v=np.linspace(0, L, NZ)
  # u=u[:-1] #leave out one point
  u,v=np.meshgrid(u,v)
  u=u.flatten()
  v=v.flatten()
  
  #evaluate the parameterization at the flattened u and v
  x=R*np.cos(u)
  y=R*np.sin(u)
  z=v
  
  #define 2D points, as input data for the Delaunay triangulation of U
  points2D=np.vstack([u,v]).T
  tri = Delaunay(points2D, incremental=True)#triangulate the rectangle U
  triSimplices = tri.simplices
  
  xyz0 = np.vstack([x,y,z]).T
  </code></pre>
  
  <p>I create a cylinder via parameterisation, and obtain the triangulation via <code>scipy.spatial.Delaunay()</code> of the base domain -- the rectangle. Obviously, this triangulation does not know about the periodicity. I can see this clearly by moving the last row, and plotting:
  <a href="https://i.stack.imgur.com/nPtB9.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/nPtB9.png" alt="enter image description here"></a></p>
  
  <p>To correct this, I try a straightforward extension of the 2D solution-- I add an extra point in 3D, re-triangulate and remove the unwanted simplices. </p>
  
  <pre><code>Tri1 = Delaunay(points2D) #triangulate the rectangle U
  Tri2 = Delaunay(xyz0)
  
  ## we add a central (0,0,L/2) point to xy0 to fill it up with triangles
  last_pt = xyz0.shape[0]
  xy1 = np.vstack((xyz0,(0,0,L/2)))  # add ctr point
  Tri3 = Delaunay(xyz1)
  print(Tri3.points.shape, Tri3.simplices.shape)
  print(Tri1.points.shape, Tri1.simplices.shape)
  print(Tri2.points.shape, Tri2.simplices.shape)
  
  ## remove the simplices that contain the central point
  mask = ~(Tri3.simplices==last_pt).any(axis=1)
  triSimplices = Tri3.simplices[mask,:]
  </code></pre>
  
  <p>However, the extension of the 2D code to 3D seems to have a big problem -- triangulations in 3D give tetrahedra, not triangles! Moreover, it seems to be more sensitive to the choice of the central point. In short, I am stuck. </p>
  
  <p>So, what is the proper way of triangulating such a periodic surface?</p>