【问题标题】:Getting a proper Delaunay triangulation of an annulus (using python)获得环的适当 Delaunay 三角剖分(使用 python)
【发布时间】:2018-06-30 13:17:22
【问题描述】:

我正在尝试使用 scipy.spatial.Delaunay() 函数对环进行三角测量,但无法获得所需的结果。这是我的代码:

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

我得到以下信息:

环本身的三角剖分清楚地给出了不需要的三角形。基本矩形的三角剖分似乎给出了正确的结果,直到您通过稍微拉伸环(即移动其节点)来意识到环实际上并未闭合

所以,我的问题是,我如何获得能够解释非平凡拓扑的正确三角剖分?我可以从环的三角剖分中移除单纯形——例如,基于键的长度——或者以某种方式将基本矩形的两端缝合在一起吗?有没有一种简单的方法可以做到这一点?

答案:

我接受了下面的答案,但它并没有完全解决所提出的问题。我仍然不知道如何使用scipy.Delaunay(即qhull 例程)平铺周期性表面。但是,使用下面定义的掩码,可以创建一个新的三角形单纯形列表,并且可以用于多种用途。但是,不能将此列表与scipy.Delaunay 类中定义的其他方法一起使用。所以,要小心!

【问题讨论】:

    标签: python numpy scipy triangulation delaunay


    【解决方案1】:

    qhull 与凸包一起使用。所以它不能直接与那个凹形的内部一起工作。在图 2 中,它用三角形填充内部。如果我们将 (0,0) 点添加到 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()
    

    删除包含该中心点的单纯形:

    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 中删除一个值似乎可行:

    u = u[:-1]
    

    在 FEM 模型中,您可能会将中心元素留在原处,但要赋予它们适当的“中性”属性(绝缘或任何有效的)。

    【讨论】:

    • 太棒了!你介意详细说明它为什么用三角形填充内部吗?
    • 我认为我无法补充qhull.org 所说的内容。关键是convex hull,在这种情况下是对应于最大v(外环)的点。它填满了整个船体。
    • 明白。另一件事:设置u = u[:-1] 对我不起作用?您具体在哪里实施?
    • 我一开始就这样做了,所以u 不会一直运行到2*pi。你一定做了类似的事情来显示第三个数字中的差距。
    • 我们可以仔细检查,但我认为qhull 不希望这些点在网格上甚至是有序的。顺序可能会影响单个三角形,但凸包应该相同。
    猜你喜欢
    • 2019-04-18
    • 2015-01-07
    • 2021-02-14
    • 2020-05-06
    • 2017-06-18
    • 2021-05-14
    • 2013-05-12
    • 1970-01-01
    • 2019-09-28
    相关资源
    最近更新 更多