Триангуляция

В двумерном случае симплексом служит треугольник, поэтому для заданного набора точек требуется построить такое семейство треугольников, которое целиком покрывает область и не имеет внутренних пересечений. Дополнительно мы требуем, чтобы итоговая сетка удовлетворяла условию Делоне: внутри описанной окружности каждого треугольника не должно быть других вершин триангуляции. Для невырожденного набора точек без заданных рёбер такая триангуляция максимизирует минимальный угол. Поэтому условие Делоне помогает уменьшить число плохо вытянутых треугольников, а устойчивость последующих конечно-элементных расчётов дополнительно зависит от качества входной геометрии, плотности точек и корректной обработки ограничений.

Строить делоневскую сетку самостоятельно — задача коварнее, чем кажется: расчётная область часто имеет форму квадрата или кольца, граница которых порождает вырожденные конфигурации, и на них наивные вычисления с плавающей точкой дают сбой (см. ниже). Поэтому низкоуровневое построение поручено зрелой библиотеке CGAL с точными геометрическими предикатами, а поверх неё реализован конвейер подготовки данных и проверок. Собственная реализация по классической схеме «разделяй и властвуй», с которой проект начинался, перенесена в приложение «Триангуляция методом „разделяй и властвуй“».

Геометрические предикаты

Почти все геометрические решения в алгоритме сводятся к двум знаковым предикатам.

Положение точки относительно ребра. Для направленного ребра с концами $p_1$ $=(x_1,y_1)$, $p_2$ $=(x_2,y_2)$ и точки $p_3$ $=(x_3,y_3)$ знак выражения

определяет, лежит ли $p_3$ слева, справа или строго на прямой $p_1 p_2$. Этот же предикат используется при построении выпуклой оболочки и при контроле ориентации треугольников.

Проверка условия Делоне. Для треугольника с вершинами $p_1, p_2, p_3$, перечисленными против часовой стрелки, и проверяемой точки $p_0$ $=(x_0,y_0)$ условие удобно записать через знак определителя

Если условие (4.2)$\displaystyle \det \begin{pmatrix} x_1 & y_1 & x_1^2 + y_1^2 & 1 \\ x_2 & y_2 & x_2^2 + y_2^2 & 1 \\ x_3 & y_3 & x_3^2 + y_3^2 & 1 \\ x_0 & y_0 & x_0^2 + y_0^2 & 1 \end{pmatrix}$ $\displaystyle \leq 0$ выполнено, точка $p_0$ лежит вне описанной окружности или на её границе; если нет — внутри, и пару соседних треугольников необходимо перестроить. При обходе вершин по часовой стрелке знак определителя меняется на противоположный, поэтому ориентацию треугольника контролирует предикат (4.1)$\displaystyle \operatorname{orient}(p_1,p_2,p_3)$ $\displaystyle = (x_2 - x_1)(y_3 - y_1)$ $\displaystyle - (y_2 - y_1)(x_3 - x_1)$.

Устойчивость на вырожденных данных

Предикаты (4.1)$\displaystyle \operatorname{orient}(p_1,p_2,p_3)$ $\displaystyle = (x_2 - x_1)(y_3 - y_1)$ $\displaystyle - (y_2 - y_1)(x_3 - x_1)$ и (4.2)$\displaystyle \det \begin{pmatrix} x_1 & y_1 & x_1^2 + y_1^2 & 1 \\ x_2 & y_2 & x_2^2 + y_2^2 & 1 \\ x_3 & y_3 & x_3^2 + y_3^2 & 1 \\ x_0 & y_0 & x_0^2 + y_0^2 & 1 \end{pmatrix}$ $\displaystyle \leq 0$ — это знаки определителей. На входах общего положения их значение заметно отличается от нуля, и обычная арифметика двойной точности даёт правильный знак. Но расчётная область в форме квадрата порождает два неприятных случая:

• коллинеарность — много граничных точек лежат строго на одной вертикали или горизонтали, и предикат (4.1)$\displaystyle \operatorname{orient}(p_1,p_2,p_3)$ $\displaystyle = (x_2 - x_1)(y_3 - y_1)$ $\displaystyle - (y_2 - y_1)(x_3 - x_1)$ обращается в точный ноль;
• концикличность — четыре угла любой ячейки регулярной сетки лежат на одной окружности, и определитель (4.2)$\displaystyle \det \begin{pmatrix} x_1 & y_1 & x_1^2 + y_1^2 & 1 \\ x_2 & y_2 & x_2^2 + y_2^2 & 1 \\ x_3 & y_3 & x_3^2 + y_3^2 & 1 \\ x_0 & y_0 & x_0^2 + y_0^2 & 1 \end{pmatrix}$ $\displaystyle \leq 0$ тоже точно равен нулю.

Беда в том, что при вычислении такого «нулевого» определителя в арифметике с плавающей точкой накапливается ошибка округления, и вместо нуля получается маленькое число произвольного знака — «шум». Разные части алгоритма, получив для одной и той же тройки точек противоречивые ответы, начинают строить несогласованную сетку: вдоль границы возникают плоские (нулевой площади) треугольники, нарушается условие Делоне.

Решение — считать знак предиката точно. Используемая библиотека вычисляет определители с помощью отфильтрованной арифметики: быстрый предварительный расчёт с гарантированной оценкой погрешности, а в сомнительных (близких к нулю) случаях — точный пересчёт. В результате знак никогда не бывает ошибочным, а истинный ноль распознаётся как ноль. Важно, что точность нужна только предикатам: алгоритм не строит новых точек — вершины сетки в точности совпадают с исходными, — поэтому точная арифметика построений не требуется, и координаты остаются обычными числами двойной точности. Такой компромисс почти не замедляет работу.

Триангуляция с ограничениями

На практике помимо самих точек нужно учитывать заданные пользователем рёбра — границы расчётной области и контуры отверстий. Такие рёбра должны обязательно присутствовать в итоговой сетке и не разрушаться локальными перестройками. Эта задача известна как триангуляция Делоне с ограничениями: вблизи ребра-ограничения условие Делоне может законно нарушаться, зато сетка гарантированно повторяет заданные контуры. Сама операция вставки ребра описана ниже, в разделе о реализации, а во всех подробностях — в приложении.

Построение сетки расчётной области

Полный конвейер построения сетки для области с границей и отверстиями состоит из шести шагов.

1. Контуры — внешний и контуры отверстий — разбиваются на отрезки длиной не больше целевого шага $\sqrt{S/n}$, где $S$ — площадь области, а $n$ — число внутренних точек; так плотность граничных точек согласуется с плотностью внутренних и снижается риск появления вытянутых треугольников вдоль границы.
2. Внутри области (вне отверстий) методом отбраковки генерируются $n$ случайных точек, при необходимости — с заданным минимальным расстоянием между ними.
3. В пустую триангуляцию по одной вставляются внутренние точки, затем граничные точки всех контуров.
4. Каждый отрезок каждого контура принудительно вставляется в сетку как ребро-ограничение.
5. Полученная триангуляция — она пока покрывает всю выпуклую оболочку — проверяется на корректность, при любом нарушении построение прерывается с ошибкой.
6. Треугольники, центр тяжести которых лежит вне области (внутри отверстия или за внешним контуром), отбрасываются — остаётся сетка собственно расчётной области.

Контроль на пятом шаге служит страховкой и включает четыре проверки:

• полнота вершин — каждая исходная точка присутствует в сетке как вершина, и в сетке нет посторонних вершин;
• число треугольников — оно должно совпадать со следующим из формулы Эйлера значением $T$ $= 2N$ $- 2$ $- h$, где $N$ — число вершин, а $h$ — число вершин на выпуклой оболочке;
• непересекаемость — никакие два треугольника не пересекаются по внутренностям рёбер;
• условие Делоне — проверяется для пар треугольников, общая сторона которых не является ребром-ограничением; около ограничений условие Делоне может законно нарушаться.

Хранение триангуляции

Готовая триангуляция хранится не просто списком треугольников, а структурой, где каждому ребру сопоставляются смежные с ним треугольники. По общему ребру всегда можно быстро найти пару треугольников, потенциально нарушающих условие Делоне, — это и используется при flip-операциях. По этой же структуре рёбра, принадлежащие единственному треугольнику, дают границу сетки: до вырезания отверстий это выпуклая оболочка. На выходе получается набор треугольных симплексов, по которым на дальнейших этапах строятся локальные матрицы жёсткости и собирается глобальная система метода конечных элементов.

Реализация: библиотека CGAL

CGAL (Computational Geometry Algorithms Library) — открытая библиотека вычислительной геометрии на C++ и де-факто стандарт в этой области. Из неё используются классы Delaunay_triangulation_2 (обычная триангуляция Делоне) и Constrained_Delaunay_triangulation_2 (триангуляция с ограничениями) с ядром Exact_predicates_inexact_constructions_kernel. Поэтапный конвейер, описанный выше, вставка ограничений и проверки пятого шага реализованы поверх неё.

Сетку CGAL строит инкрементально — точки вставляются по одной в уже готовую триангуляцию.

• Вместо вспомогательного «супертреугольника», который пришлось бы подбирать по размеру и удалять в конце, используется бесконечная вершина: выпуклая оболочка текущей сетки замыкается фиктивными треугольниками, третьей вершиной которых служит формальная «точка на бесконечности».
• Очередная точка сначала локализуется «прогулкой» по сетке: от стартового треугольника алгоритм по предикату ориентации (4.1)$\displaystyle \operatorname{orient}(p_1,p_2,p_3)$ $\displaystyle = (x_2 - x_1)(y_3 - y_1)$ $\displaystyle - (y_2 - y_1)(x_3 - x_1)$ шагает от соседа к соседу в сторону точки.
• Найденный треугольник делится вставляемой точкой на три (точка на ребре делит два смежных треугольника на четыре), после чего условие Делоне восстанавливается каскадом flip-операций, расходящимся от новой вершины: если по предикату (4.2)$\displaystyle \det \begin{pmatrix} x_1 & y_1 & x_1^2 + y_1^2 & 1 \\ x_2 & y_2 & x_2^2 + y_2^2 & 1 \\ x_3 & y_3 & x_3^2 + y_3^2 & 1 \\ x_0 & y_0 & x_0^2 + y_0^2 & 1 \end{pmatrix}$ $\displaystyle \leq 0$ четвёртая вершина соседнего треугольника попадает в описанную окружность, общая диагональ четырёхугольника заменяется на противоположную (рисунок Замена диагонали).

Ребро-ограничение вставляется в готовую сетку отдельной операцией: все пересечённые отрезком треугольники удаляются, образовавшаяся пустота делится отрезком на две цепочки вершин, и каждая цепочка заново триангулируется; при последующих перестройках flip через ребро-ограничение запрещён. Поскольку к моменту вставки оба конца каждого отрезка уже являются вершинами сетки, а контуры не пересекаются, новые (штейнеровы) точки не появляются. Подробный разбор этой операции — в приложении.

Выбранное ядро сочетает точные предикаты с неточными построениями. Предикаты вычисляются по схеме фильтрации: сначала интервальная арифметика с гарантированной оценкой погрешности и лишь тогда, когда интервал накрывает ноль, — точный пересчёт в многоразрядной арифметике. Именно это обеспечивает корректную обработку коллинеарных и конциклических конфигураций квадрата. Построения же остаются обычными числами двойной точности: новых точек алгоритм не создаёт, поэтому каждая вершина итоговой сетки без потерь сопоставляется исходной точке. На случайном порядке вставки — а наш генератор случайных точек его и даёт — ожидаемая трудоёмкость построения составляет $O(n \log n)$.

Примеры

Все сетки ниже — реальные результаты описанного конвейера (утилита triangulate вычислительного ядра), отрисованные без какой-либо ретуши.

На рисунке «Квадрат» показана триангуляция единичного квадрата при $n$ $= 200$ внутренних точках: после разбиения границы с шагом $\sqrt{S/n}$ сетка содержит 260 вершин и 458 треугольников. Граница квадрата — это десятки строго коллинеарных точек, тот самый вырожденный вход, на котором плавающая арифметика без точных предикатов даёт сбой; здесь же вдоль всей границы сетка корректна, плоских треугольников нет.

Рисунок «Квадрат с отверстием» демонстрирует работу рёбер-ограничений: в квадрате вырезано квадратное отверстие, оба контура вставлены как ограничения ($n$ $= 300$; 396 вершин, 696 треугольников). Сетка точно повторяет контур отверстия и не заходит внутрь — треугольники, чей центр тяжести попал в отверстие, отброшены на финальном шаге конвейера.

Наконец, на рисунке «Кольцо» показана триангуляция кольца (внешний радиус 1, внутренний 0.4) при $n$ $= 1000$: итоговая сетка содержит 1192 вершины и 2192 треугольника. Оба контура аппроксимированы многоугольниками (128 и 64 вершины) и вставлены как рёбра-ограничения, поэтому сетка повторяет заданную полигональную границу области, а внутренние точки порождают близкую к равномерной делоневскую сетку.