获取正确的环形Delaunay三角剖分(使用python) [英] Getting a proper Delaunay triangulation of an annulus (using python)

问题描述

我正在尝试使用scipy.spatial.Delaunay()函数对环形空间进行三角剖分，但无法获得所需的结果.这是我的代码:

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?

我接受了以下答案，但并不能完全解决所问的问题.我仍然不知道如何使用scipy.Delaunay(即qhull例程)平铺周期性曲面.但是，使用以下定义的遮罩，可以创建一个新的三角单纯形列表，并且应该用于许多目的.但是，一个列表不能与scipy.Delaunay类中定义的其他方法一起使用.所以，要小心！

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适用于凸包.因此，它不能直接用于那种凹入的内部.在图2中，它用三角形填充内部.如果在xy0上添加(0,0)点，则可能会更明显.

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

要将两端缝合在一起，请从u中删除一个值似乎可行:

To stitch the two ends together, removing a value from u seems to work:

u = u[:-1]

在FEM模型中，您可以将中心元素保留在适当的位置，但为它们提供适当的中性"属性(绝缘或其他有效方法).

In a FEM model you might leave the center elements in place, but give them the appropriate 'neutral' properties (insulating or whatever works).

这篇关于获取正确的环形Delaunay三角剖分(使用python)的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！