2
votes

A year ago, I asked about properly triangulating a periodic shape on a plane: the annulus (Getting a proper Delaunay triangulation of an annulus (using python)).

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:

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

I create a cylinder via parameterisation, and obtain the triangulation via scipy.spatial.Delaunay() 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: enter image description here

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.

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,:]

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.

So, what is the proper way of triangulating such a periodic surface?

What is your expected output an open or closed cylinder? Are you aware of the matplotlib.tri.Triangulation method that implements Delaunay triangulation for matplotlib ? Here an example of usage. - jlandercy
@jlandercy I am looking for a closed cylinder. Yes, I am aware of matplotlib's method, but I can't see how that will do anything different unless I specifically ask it to do so. - ap21
What is wrong with the initial triangulation? - lrineau
@lrineau It doesn't reflect the fact that the surface is periodic. The triangulation is a discrete representation of a continuous surface, so it should reflect any periodicities as well. How do I make it do so? - ap21