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

Строить делоневскую сетку самостоятельно — задача коварнее, чем кажется: расчётная область часто имеет форму квадрата или кольца, граница которых порождает вырожденные конфигурации, и на них наивные вычисления с плавающей точкой дают сбой (см. ниже). Поэтому низкоуровневое построение поручено зрелой библиотеке 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)$$\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} \leq 0$$ выполнено, точка $p_0$ лежит вне описанной окружности или на её границе; если нет — внутри, и пару соседних треугольников необходимо перестроить. При обходе вершин по часовой стрелке знак определителя меняется на противоположный, поэтому ориентацию треугольника контролирует предикат (4.1)$$\operatorname{orient}(p_1,p_2,p_3) = (x_2 - x_1)(y_3 - y_1) - (y_2 - y_1)(x_3 - x_1)$$.

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

Предикаты (4.1)$$\operatorname{orient}(p_1,p_2,p_3) = (x_2 - x_1)(y_3 - y_1) - (y_2 - y_1)(x_3 - x_1)$$ и (4.2)$$\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} \leq 0$$ — это знаки определителей. На входах общего положения их значение заметно отличается от нуля, и обычная арифметика двойной точности даёт правильный знак. Но расчётная область в форме квадрата порождает два неприятных случая:

• коллинеарность — много граничных точек лежат строго на одной вертикали или горизонтали, и предикат (4.1)$$\operatorname{orient}(p_1,p_2,p_3) = (x_2 - x_1)(y_3 - y_1) - (y_2 - y_1)(x_3 - x_1)$$ обращается в точный ноль;
• концикличность — четыре угла любой ячейки регулярной сетки лежат на одной окружности, и определитель (4.2)$$\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} \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)$$\operatorname{orient}(p_1,p_2,p_3) = (x_2 - x_1)(y_3 - y_1) - (y_2 - y_1)(x_3 - x_1)$$ шагает от соседа к соседу в сторону точки.
• Найденный треугольник делится вставляемой точкой на три (точка на ребре делит два смежных треугольника на четыре), после чего условие Делоне восстанавливается каскадом flip-операций, расходящимся от новой вершины: если по предикату (4.2)$$\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} \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 вершины) и вставлены как рёбра-ограничения, поэтому сетка точно повторяет форму области, а внутренние точки порождают равномерную делоневскую сетку.