Разработка и реализация алгоритмов трехмерной триангуляции сложных пространственных областей: итерационные методы (стр. 4 из 6)
Таким образом, алгоритм можно разбить на четыре этапа:
1) построение триангуляции Делоне без ограничений;
2) восстановление "ребер" ограничений ("edge recovery");
3) восстановление поверхностей ограничений ("face recovery");
4) отсечение "лишних" тетраэдров, оказавшихся вне границы заданной области ("trimming").
Рис. 9. Три возможных варианта пересечения тетраэдра плоскостью и соответствующие способы разбиения образующихся частей на тетраэдры: а) плоскость не пересекает вершины; б) плоскость пересекает две вершины; в) плоскость пересекает только одну вершину.
Подробнее остановимся на каждом этапе.
Этап 1. Несмотря на кажущуюся простоту, этот этап имеет несколько важных нюансов. Во-первых, как уже было сказано, все "углы" поверхностей ограничения (включая концы "незакрепленных" "ребер"), должны в обязательном порядке войти в набор узлов, по которым строится триангуляция. Во-вторых, поскольку многогранник, получающийся при использовании описанного выше алгоритма, не обязательно будет выпуклым, вполне возможна ситуация, при которой некоторые участки поверхностей ограничений окажутся вне этого многогранника. Чтобы не допустить этого, можно использовать слой дополнительных узлов, окружающий заданную область. На 4 этапе все тетраэдры, образованные с этими узлами, удаляются наряду с остальными тетраэдрами, оказавшимися "по ту сторону границы".
Этап 2. Восстановление "ребер" производится с помощью дополнительных узлов[4]. Заметим, что если в двумерном случае всегда возможно построить триангуляцию Делоне с ограничениями без использования дополнительных узлов, то в трехмерном случае это, как правило, невозможно.
Дополнительные узлы вставляются в точки пересечений "ребер" ограничений с гранями и ребрами тетраэдров. Заметим, что вариантов пересечения "ребра" с тетраэдром всего шесть (рис. 10). Соответственно, для каждого варианта используется свой шаблон разбиения тетраэдра. При этом встает вопрос согласования триангуляции на гранях соседних тетраэдров, поскольку некоторые случаи (рис. 10. в) допускают различные способы разбиения граней. Это проблема решается установкой следующих двух правил: 1) из всех вариантов всегда выбирается самое короткое ребро; 2) если ребра равны, проводится ребро от узла с самым меньшим глобальным номером (например).
Вспомним, что тетраэдр может пересекаться не одним "ребром", а сразу несколькими. Чтобы избежать ввода дополнительных шаблонов (поскольку, теоретически, "ребер", пересекающих тетраэдр, может быть неограниченно много), достаточно обрабатывать каждое "ребро" отдельно, по очереди. При обработке каждого "ребра" остальные, еще не обработанные, не учитываются, а после обработки "ребро" уже аппроксимировано цепочкой ребер триангуляции и естественным образом входит в сетку, поэтому указанная проблема исчезает.
Этап 3. На предыдущем этапе все "ребра" ограничений были аппроксимированы цепочками ребер сетки. Теперь нужно добиться того, чтобы все поверхности ограничений были представлены множеством смежных граней. Это также осуществляется вставкой дополнительных вершин - в точках пересечения поверхностей ограничения с ребрами сетки (см. рис. 9). Опять встает вопрос об однозначности разбиения граней, но теперь он несколько сложнее. В случае, если поверхность пересекает 3 ребра тетраэдра (рис. 9. а), то возможен такой вариант проведения ребер, при котором нижняя часть элемента - пятигранник (усеченный тетраэдр) - не разбивается на тетраэдры. Разрешить эту ситуацию можно с помощью дополнительного узла, вставляемого внутрь пятигранника. Тогда согласование ребер можно осуществлять по правилам, изложенным выше.
Этап 4. Удаление лишних тетраэдров, лежащих вне границы заданной области, не представляет какой-либо проблемы. Зато можно сказать, что после их удаления граница области будет полностью аппроксимирована гранями и ребрами построенной триангуляции. То же самое относится и к внутренним ребрам и поверхностям ограничений, если таковые были.
Рис. 10. Возможные варианты пересечения "ребра"ограничения с тетраэдром и соответствующие способы разбиения тетраэдра на части: а) вершина - ребро; б) вершина - грань; в) ребро - ребро; г) ребро - грань; д) ребро - вырожденный случай; е) грань - грань.
На рис. 11 представлен результат триангуляции с помощью описанного алгоритма области, являющей собой объединение икосаэдра и додекаэдра[5].
Рис. 11. Триангуляция области, представляющей собой объединение додекаэдра и икосаэдра (триангуляция Делоне с ограничениями)
Качество сеток, построенных данным методом, находится на среднем уровне (тетраэдры у границ могут иметь очень плохую форму), поэтому обычно дополнительно прибегают к одному из методов оптимизации.
В работах Б. Джо [23, 24] предложены другие варианты алгоритма, не использующие дополнительных точек и полностью основанные на локальных трансформациях, аналогичных "трейду".
3.3. Особенности технической реализации алгоритмов на основе критерия Делоне
Общий недостаток всех методов на основе критерия Делоне - крайне высокая чувствительность к точности машинных вычислений. Многие вычислительные процедуры, используемые в этих методах (нахождение центра и радиуса описанной окружности, проверка компланарности и коллинеарности векторов и т.п.) представляют собой плохо обусловленные задачи, а их неизбежное интенсивное использование столь же неизбежно влечет к накоплению ошибок округления, что в итоге может привести к ошибкам в структуре сетки или зацикливанию алгоритма.
Типичная ошибка, допускаемая в реализации этого класса алгоритмов, заключается в использовании приближенных вычислений и операций, например, использование для сравнения не нуля (x>0.0), а заданной точности (x>EPS или x>-EPS). Такой подход вполне оправдан и верен во многих других случаях, но только не в методах на основе критерия Делоне.
К сожалению, простое увеличение точности не дает существенных результатов. Использование числовых типов с удвоенной точностью перестает работать уже на задачах средней сложности (сетки с несколькими тысячами узлов). Частичным решением этой проблемы может быть использование числовых типов с фиксированной запятой. За счет отдельного хранения в 24-байтной структуре целой и десятичной частей числа (в виде целых чисел), на указанных выше плохо обусловленных задачах можно добиться точности вплоть до 10 знака, что является уже более-менее допустимым результатом [11].
Также возможно использование так называемой "точной арифметики", модули для реализации которой уже разработаны для многих популярных прикладных языков программирования (С++, Фортран и др.). В точной арифметике все иррациональные числа представляются как набор предикатов (арифметических и алгебраических операторов) от рациональных/целых чисел и установленных констант (таких, как
, 
и т.п.). То есть, например, в точной арифметике 
- это действительно 
, а не 1,73205080756887729...Недостатком обоих подходов является (весьма существенное) снижение скорости вычислений, так как ни точная арифметика, ни числа с фиксированной запятой пока аппаратно не поддерживаются процессорами современных персональных компьютеров. Поэтому все эти операции приходится реализовывать на уровне подпрограмм, что приводит к дополнительным затратам ресурсов.
Иной путь улучшения ситуации - оптимизация процедур вычисления. Так, например, сотрудники исследовательского центра IBM Павло Кавальканти и Улисс Мелло сумели создать устойчивый алгоритм построения триангуляции Делоне, сведя все геометрические операции к двум предикатам: проверке условия "пустой сферы" и нахождению положения точки относительно заданной плоскости [10].
Впервые идея метода исчерпывания была предложена Рейнальдом Лонером (R. Lohner), а его трехмерный вариант разработал профессор Гонконгского университета С.Х. Ло. Алгоритм Ло вот уже многие годы успешно используется в программном комплексе ANSYS для дискретизации произвольных объемных областей.
Общая идея этого класса методов заключается в последовательном изымании из заданной области фрагментов тетраэдрической формы до тех пор, пока вся область не окажется "исчерпана". Впрочем, этот легко представимый в воображении процесс на самом деле имеет мало общего с практической реализацией. Английское название "advancing front" (что можно перевести как "продвижение фронта"), пожалуй, лучше отражает сущность алгоритмов.
Отправной точкой всех алгоритмов этого класса является некая триангуляция границы заданной области, причем не грубая, а наиболее полно соответствующая требованиям разработчика, поскольку в дальнейшем никаких изменений этой триангуляции не будет. Заметим, что подобное требование к начальным данным может быть как положительной стороной метода, так и отрицательной. Если никаких ограничений на границу области не накладывается, то это приводит к необходимости ее отдельной триангуляции, что само по себе является самостоятельной задачей. С другой стороны, если необходимо обеспечить согласование триангуляции в сложной области, либо же триангуляция поверхности задана изначально, то использование метода исчерпывания будет самым удачным выбором.