scipy Python:圆柱(或任何周期性曲面)的Delaunay三角剖分

lx0bsm1f  于 2023-10-20  发布在  Python
关注(0)|答案(1)|浏览(96)

一年前,我问过如何正确地在平面上三角形化一个周期性的形状:环(Getting a proper Delaunay triangulation of an annulus (using python))。
现在我想把它扩展到对一个圆柱体(或者一般来说,任何周期性的表面)进行三角测量。我尝试了一个简单的2D代码扩展:

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

我通过参数化创建了一个圆柱体,并通过scipy.spatial.Delaunay()获得了基本域的三角剖分-矩形。显然,这种三角测量不知道周期性。我可以通过移动最后一行并绘制图来清楚地看到这一点:

为了纠正这一点,我尝试了一个简单的扩展二维解决方案-我在三维中添加一个额外的点,重新三角化并删除不需要的单形。

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

然而,将二维代码扩展到三维似乎有一个大问题--三维中的三角剖分给予四面体,而不是三角形!而且,它似乎对中心点的选择更敏感。简而言之,我被卡住了。
那么,对这样一个周期性曲面进行三角剖分的正确方法是什么呢?

mm9b1k5b

mm9b1k5b1#

下面的简短脚本创建顶点,并使用scipy创建三角形。然后重新分配最后一列

# Create a 2d map of the vertices
num_tri_high = 10
num_tri_in_circum = 9
xx = np.linspace(0, 10, num_tri_in_circum + 1)
yy = np.linspace(0, 10, num_tri_high)
xxx,yyy = np.meshgrid(xx,yy)

# Triangulate vertices
receiver_tri = Delaunay(np.transpose(np.vstack((xxx.flatten(), yyy.flatten()))))
triangles = receiver_tri.simplices

# Set last column equal to 1st (which effectively closes the gap)
tot = (num_tri_in_circum + 1) * num_tri_high
for i in range(num_tri_high+1):
    triangles[triangles == (num_tri_in_circum+1) * i + num_tri_in_circum] = (num_tri_in_circum+1) * i

# Wrap into a cylinder
radius = 2
xp = radius*np.sin(xxx.flatten()/np.max(xxx.flatten())*np.pi*2)
yp = radius*np.cos(xxx.flatten()/np.max(xxx.flatten())*np.pi*2)
zp = yyy.flatten()
tow_scene = mlab.triangular_mesh(xp.flatten(),yp.flatten(),zp.flatten(), triangles, representation='wireframe', color=(0,0,0), opacity=1)
mlab.show()

相关问题