Получение правильной триангуляции Делоне кольца (с использованием python)

Я пытаюсь триангулировать кольцо с помощью функции 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. Так что будь осторожен!


person ap21    schedule 21.01.2018    source источник


Ответы (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 вы можете оставить центральные элементы на месте, но придать им соответствующие «нейтральные» свойства (изолирующие или любые другие).

введите здесь описание изображения введите здесь описание изображения

person hpaulj    schedule 22.01.2018
comment
Большой! Не могли бы вы уточнить, почему он заполнил интерьер треугольниками? - person ap21; 22.01.2018
comment
Я не думаю, что могу что-то добавить к тому, что говорит qhull.org. Ключ — это convex hull, которые в данном случае являются точками, соответствующими самому большому v, внешнему кольцу. Он заполняет весь этот корпус. - person hpaulj; 22.01.2018
comment
Понял. Другое дело: настройка u = u[:-1] у меня не работает? Где именно вы это реализуете? - person ap21; 23.01.2018
comment
Я делаю это в самом начале, поэтому u не доходит до 2*pi. Вы должны были сделать что-то подобное, чтобы показать разрыв на 3-м рисунке. - person hpaulj; 23.01.2018
comment
Мы можем перепроверить, но я не думаю, что qhull ожидает, что точки будут в сетке или даже упорядочены. Порядок может влиять на отдельные треугольники, но выпуклая оболочка должна быть одинаковой. - person hpaulj; 23.01.2018
comment
Я согласен с этим пунктом выше. И последнее: можно ли применить эту маску постоянно к полному объекту Delaunay triangulation, чтобы мне не нужно было носить ее с собой через программу? Или даже просто Delaunay.vertices? - person ap21; 23.01.2018
comment
Может быть, у вас есть ответ на мой последний вопрос? Можно ли наносить маску постоянно? - person ap21; 26.01.2018