«Технология построения неструктурированных сеток и монотонная дискретизация уравнения диффузии ...»

Учреждение Российской академии наук

Институт вычислительной математики РАН

На правах рукописи

Данилов Александр Анатольевич

Технология построения

неструктурированных сеток

и монотонная дискретизация

уравнения диффузии

05.13.18 – Математическое моделирование,

численные методы и комплексы программ

ДИССЕРТАЦИЯ

на соискание учёной степени

кандидата физико-математических наук

Научный руководитель д. ф.-м. н.

Василевский Юрий Викторович Москва – 2010 Содержание Введение................................... 4 Обзор методов построения сеток и комплексов программ.... Глава 1. Построение двумерных треугольных сеток....... 1.1. Обозначения............................. 1.2. Алгоритм продвигаемого фронта................. 1.3. Влияние вычислительных погрешностей............. 1.4. Конечность работы алгоритма................... 1.5. Скорость работы алгоритма продвигаемого фронта....... 1.6. Улучшение качества полученной сетки.............. 1.7. Результаты экспериментов..................... 1.8. Выводы к первой главе....................... Глава 2. Построение поверхностных сеток............. 2.1. Представление поверхности..................... 2.2. Взаимодействие с геометрическим ядром САПР......... 2.3. Алгоритм продвигаемого фронта................. 2.4. Результаты экспериментов..................... 2.5. Выводы к второй главе....................... Глава 3. Построение трёхмерных сеток............... 3.1. Алгоритм продвигаемого фронта................. 3.2. Устойчивый метод на основе тетраэдризации Делоне...... 3.3. Улучшение качества полученной сетки.............. 3.4. Результаты экспериментов..................... 3.5. Выводы к третьей главе....................... Глава 4. Сетки с многогранными ячейками и монотонная дис­ кретизация уравнения диффузии.................. 4.1. Уравнение диффузии........................ 4.2. Монотонная нелинейная схема на основе метода конечных объ­ ёмов.................................. 4.3. Схема дискретизации и анализ её монотонности......... 4.4. Результаты экспериментов..................... 4.5. Моделирование стационарного распределения интерферона в лимфатическом узле........................ 4.6. Выводы к четвёртой главе..................... Заключение.................................. Литература.................................. Введение При решении современных научных и инженерных задач широко при­ меняется математическое моделирование. Существуют готовые системы ав­ томатизации инженерных расчётов (CAE – Computer-aided engineering) для анализа и моделирования физических процессов на ЭВМ. Как правило, это большие закрытые коммерческие проекты, разрабатывающиеся на протяже­ нии десятков лет группами в несколько сотен специалистов. С точки зрения пользователя системы CAE предоставляют удобный графический интерфейс для задания параметров и условий задачи, проведения расчётов, отображе­ ния и анализа результатов. Алгоритмы и методы, используемые для решения конкретной задачи, зачастую скрыты от пользователя.

Как правило, процесс моделирования физических процессов с использо­ ванием ЭВМ сводится к трём основным операциям: построение расчётных сеток, дискретизация дифференциального уравнения, решение системы ал­ гебраических уравнений. В диссертационной работе мы детально рассмотрим проблемы построения сеток. Различные методы дискретизаций представлены в работах [54, 55, 58]. Методы решения алгебраических уравнений рассматри­ ваются в работах [53, 56, 60].

Большинство современных комплексов CAE состоит из набора библиотек для графического интерфейса, визуализации, построения расчётных сеток, дискретизации дифференциальных уравнений и решения систем линейных и нелинейных уравнений. Разработку большинства низкоуровневых библиотек можно вынести за пределы всей системы CAE, что позволяет значительно упростить внедрение новых методов и алгоритмов для решения отдельных подзадач.

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

Автор диссертационной работы входит в коллектив исследователей, раз­ рабатывающих пакет библиотек Ani3D с открытым исходным кодом – про­ граммный код библиотеки доступен для просмотра, изучения и изменения.

В этот пакет входят библиотеки для построения и адаптации неструктуриро­ ванных сеток, для построения дискретизаций на неструктурированных сет­ ках, для решения систем линейных и нелинейных алгебраических уравнений.

Также развивается и “двумерная” версия пакета – Ani2D. Программный код обеих библиотек распространяется свободно.

Автор работы является основным разработчиком библиотек для постро­ ения неструктурированных сеток. Конечной целью “трёхмерной” библиоте­ ки является создание технологии надёжного построения гибридных сеток с многогранными ячейками для произвольных областей. В рамках диссерта­ ционной работы автором была разработана технология надёжного построе­ ния неструктурированных тетраэдральных сеток, которая уже включена в пакет Ani3D. Также автором диссертационной работы была разработана но­ вая монотонная нелинейная схема дискретизации уравнения диффузии на неструктурированных сетках с многогранными ячейками.

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

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

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

Научная новизна. В работе предложена технология автоматического надёжного построения тетраэдральных сеток для сложных областей на осно­ ве метода продвигаемого фронта. Проведён анализ влияния вычислительных погрешностей на алгоритмы при их реализации на ЭВМ, представлено теоре­ тическое обоснование конечности работы предложенных алгоритмов. Предло­ жена и исследована монотонная нелинейная схема на основе метода конечных объёмов для уравнения диффузии на неструктурированных сетках с много­ гранными ячейками.

Практическая значимость. Практическая значимость диссертацион­ ной работы заключается в создании и поддержке комплекса программ для автоматического построения тетраэдральных сеток и для приближённого ре­ шения задач диффузии на многогранных сетках. Программы для построения поверхностных треугольных и тетраэдральных сеток включены в библиоте­ ку программ Ani3D и находятся в свободном доступе. Решена практическая задача о стационарном распределении интерферона в лимфоузле.

На защиту выносятся следующие основные результаты:

1. Разработана технология надёжного построения неструктурированных тетраэдральных сеток, реализованы разные способы задания области, в том числе с помощью САПР.

2. Предложена и исследована новая монотонная нелинейная схема дискре­ тизации уравнения диффузии на основе метода конечных объёмов на сетках с многогранными ячейками.

3. На основе предложенных методов разработана модель стационарного распределения интерферона в лимфоузле с учётом геометрических осо­ бенностей лимфоузла.

Апробация работы. Результаты диссертационной работы докладыва­ лись автором и обсуждались на научных семинарах Института вычислитель­ ной математики РАН, Института математического моделирования РАН, Вы­ числительного центра РАН, Механико-математического факультета МГУ им.

М. В. Ломоносова, в Институте прикладной математики и информационных технологий, г. Павия (Италия), в Лос-Аламосской национальной лаборато­ рии, г. Лос-Аламос (США) и на следующих научных конференциях: конфе­ ренции “Тихоновские чтения” (МГУ, Москва, ноябрь 2007, октябрь 2009); кон­ ференции “Ломоносов” (МГУ, Москва, апрель 2008, апрель 2010); конферен­ ции “Ломоносовские чтения” (МГУ, Москва, апрель 2009, апрель 2010); конфе­ ренция молодых учёных “Технологии высокопроизводительных вычислений и компьютерного моделирования” (СПбГУ ИТМО, С.-Петербург, апрель 2009);

конференция “Лобачевские чтения” (КГУ, Казань, ноябрь 2009); международ­ ные конференции “NUMGRID-2006” и “NUMGRID-2008” (ВЦ РАН, Москва, июнь 2006, июнь 2008); международная конференция “SIAM Geosciences 2009” (Лейпциг, Германия, июнь 2009); международный научный семинар “Computational Mathematics and Applications” (Технологический университет Тампе­ ре, Тампере, Финляндия, сентябрь 2009).

Публикации. Основные материалы диссертации опубликованы в 6 пе­ чатных работах, из них 3 статьи в рецензируемых журналах, входящих в перечень ВАК, [8, 9, 49].

Личный вклад автора. В совместной работе [9] вклад автора заклю­ чался в разработке нелинейной схемы дискретизации уравнения диффузии в трёхмерном пространстве, в программной реализации метода и в постановке численных экспериментов.

В совместной работе [43] вклад автора заключался в разработке методов построения треугольных и тетраэдральных неструктурированных сеток.

В совместной работе [44] вклад автора заключался в разработке методов взаимодействия с САПР при построении сеток, в программной реализации алгоритмов и проведении численных экспериментов.

Структура и объём диссертации. Диссертационная работа состоит из введения, обзора методов построения сеток и комплексов программ, четы­ рёх глав, заключения и списка литературы из 73 наименований. Диссертаци­ онная работа содержит 32 рисунка и 15 таблиц. Общий объём диссертацион­ ной работы – 148 страниц.

Благодарности Автор диссертационной работы выражает глубокую признательность на­ учному руководителю Ю. В. Василевскому за продолжительную поддержку, ценные советы и плодотворное обсуждение вопросов. Автор благодарен Р. Гаримелле из Лос-Аламосской национальной лаборатории за идею использова­ ния библиотеки CGM в качестве интерфейса между генератором сеток и гео­ метрическим ядром САПР и за возможность использования при программ­ ной реализации новой монотонной схемы на основе метода конечных объё­ мов библиотеки MSTK для работы с многогранными сетками. Автор также выражает благодарность К. Д. Никитину, И. В. Капырину, К. Н. Липникову, Д. А. Святскому и многим другим за помощь в обсуждении идей и методов.

Особую благодарность автор выражает Г. А. Бочарову за постановку задачи о стационарном распределении интерферона в лимфоузле, за подбор необхо­ димых параметров и за помощь в интерпретации полученных результатов.

Работа над диссертацией была частично поддержана грантами РФФИ 08-01-00159-а, 09-01-00115-а, программой Президиума РАН “Фундаменталь­ ные науки – медицине”, федеральной целевой программой “Научные и научно­ педагогические кадры инновационной России”, а также грантом Upstream Research Center, ExxonMobil corp.

Обзор методов построения сеток и комплексов В диссертационной работе будем рассматривать только конформные сет­ ки, любые два элемента которых либо не имеют общих точек, либо имеют ровно одну общую вершину, либо одно общее ребро, либо одну общую грань.

Важным примером неконформных сеток могут служить сетки типа восьме­ ричное дерево (см. Рис. 1),тем не менее их можно рассматривать как кон­ формные сетки с многогранными ячейками.

На плоскости можно выделить два класса сеток – треугольные и четы­ рёхугольные. В трёхмерном пространстве выделим три класса сеток – тет­ раэдральные, треугольные призматические, и гексаэдральные. При фиксиро­ ванном количестве вершин сетки выполняются следующие оценки на ко­ личество элементов,, и граней,, в сетке. Для неструктурированных тет­ раэдральных сеток, используемых в приложениях, 5.5, 11.

Для треугольных призматических сеток 2, 5, для гексаэд­ С точки зрения количества элементов и граней гексаэдральные сетки являются наиболее удобными, так как получающиеся при дискретизации на них матрицы будут иметь меньше ненулевых значений. С точки зрения авто­ матического построения сеток для сложных областей, гексаэдральные сетки являются самыми неудобными. Существующие на сегодняшний день методы построения сеток, такие как метод отображений, метод продвигаемого фрон­ та, методы заметания (sweeping) и замощения (paving) области, методы на ос­ нове восьмеричных деревьев (octree), применимы либо для достаточно узкого класса областей, либо требуют ручного вмешательства в процесс построения сетки.

Призматические сетки широко применяются в геофизических приложе­ ниях в задачах со слоистыми областями. В таких областях призматическую сетку можно представить как тензорное произведение неструктурированной треугольной сетки в плоскости и одномерной сетки вдоль. В Инсти­ туте вычислительной математики РАН В. Н. Чугуновым разработан генера­ тор призматических сеток [63]. Для построения начальной треугольной сет­ ки используется разработанный автором диссертационной работы генератор треугольных сеток из библиотеки программ Ani2D. В будущем планируется включение этого призматического генератора в состав Ani3D.

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

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

Параллельно с разработкой методов построения гибридных сеток необ­ ходима разработка схем дискретизации на новых сетках. В диссертационной работе будет предложена новая монотонная схема дискретизации уравнения диффузии на основе метода конечных объёмов на сетках с многогранными ячейками.

Проблемам построения сеток уделяется большое внимание, постоянно выходят научные публикации, книги, сборники [37], проводятся ежегодные конференции. Наиболее известные международные конференции, посвящён­ ные проблемам построения сеток, проходят в США – International Meshing Roundtable, и в России – “Численная геометрия, построение сеток и высоко­ производительные вычисления” (NUMGRID), результаты работ которой пуб­ ликуются в специальных выпусках Журнала вычислительной математики и математической физики. Сделаем краткий обзор основных методов и наибо­ лее известных комплексов программ для построения треугольных и тетра­ эдральных сеток. Краткий обзор имеющихся схем дискретизации уравнения диффузии будет сделан в Главе 4.

Построение треугольных сеток на плоскости Существует два основных метода построения треугольной сетки на плос­ кости: 1) метод продвигаемого фронта [27, 30]; 2) метод построения триангуля­ ции Делоне [2, 59]. Идея первого метода будет раскрыта в Главе 1, подробно­ сти использования второго метода хорошо представлены в [59]. Отметим, что оба метода могут быть использованы для построения конформной треуголь­ ной сетки внутри заданного многоугольника на заданном наборе вершин без добавления новых вершин. На практике сетки строятся с добавлением но­ вых вершин внутри области, и качество и свойства сетки зависят от способа выбора положения этих новых вершин [2, 27].

В ходе реализации и тестирования метода продвигаемого фронта в соста­ ве пакетов программ Ani2D и Ani3D автором был разработан набор эвристи­ ческих алгоритмов для выбора положения новых вершин с целью увеличения надёжности и качества полученных сеток при сохранении простоты самого алгоритма. В Главе 1 проведён анализ влияния вычислительных погрешно­ стей, доказана конечность работы алгоритма, и сделаны оценки на время работы алгоритма.

Отметим основные особенности комплекса программ для построения се­ ток, разработанного автором диссертационной работы и включённого в паке­ ты Ani2D и Ani3D. Граница задаётся либо дискретно, либо аналитически с помощью параметризующих функций, предоставленных пользователем; есть возможность контролировать желаемый размер треугольников во всей об­ ласти, размер треугольников может выбираться автоматически на основе за­ данной дискретизации границы и заданной скорости увеличения размера тре­ угольников при продвижении внутрь области. При использовании одинарной точности для хранения координат в памяти ЭВМ для корректной работы ге­ нератора достаточна двойная точность вычислений внутри генератора. Про­ грамма собирается в виде библиотеки, приведены примеры её использования.

Исходный код открыт и свободно распространяется [65, 66].

Существует множество программных реализаций для построения тре­ угольных сеток. Отметим один из наиболее известных надёжных генераторов треугольных сеток, генератор сеток Triangle, разработанный Дж. Р. Шевчуком [33]. Этот генератор сеток основан на методе построения триангуляции Де­ лоне, и позволяет строить конформные триангуляции Делоне, триангуляции Делоне с ограничениями, и треугольные сетки с высоким качеством [36]. Гра­ ница области задаётся дискретно, нет возможности полностью контролиро­ вать желаемый размер элементов внутри области. Генератор может использо­ вать арифметику с произвольной точностью для исключения ошибок округ­ ления при вычислениях [34], например, при вычислении знака определителя 3 3. Код программы может быть скомпилирован как готовый исполняемый файл, или собран в виде библиотеки для дальнейшего подключения к другим программам. Исходный код генератора открыт и свободно распространяется [64].

Построение тетраэдральных сеток в пространстве Обычно процесс построения тетраэдральной сетки делится на два этапа:

1) построение поверхностной триангуляции на границе области; 2) построение тетраэдральной сетки внутри области.

В разных задачах геометрические формы моделей имеют разную слож­ ность. Например, в случае моделирования физических процессов в однород­ ной среде можно использовать область простой формы в виде шара или куба.

В задачах обтекания используются более сложные модели, например, пересе­ чение или объединение нескольких простых фигур. В инженерных задачах широко применяются системы автоматизации проектных работ (САПР) для создания геометрических моделей изделий. Не исключено и применение трёх­ мерного сканирования для получения компьютерной модели реально суще­ ствующего объекта, а также использование картографических данных для создания моделей реального рельефа, например, в геофизических задачах.

Универсального подхода не существует: в разных задачах удобен тот или иной способ задания геометрии модели.

В пакете программ Ani3D предусмотрен богатый выбор способов постро­ ения треугольных поверхностных сеток [49], включающий библиотеку для улучшения качества имеющихся поверхностных триангуляций, разработан­ ную совместно с А. В. Плёнкиным и А. В. Вершининым [43], и библиотеку для задания области на основе подходов конструктивной блочной геометрии, раз­ работанную К. Д. Никитиным [57]. В диссертационной работе будет рассмот­ рен только способ построения треугольных поверхностных сеток для обла­ стей, заданных аналитически, или с помощью САПР [44]. Подробнее об этих методах будет рассказано в Главе 2.

При построении объёмной тетраэдральной сетки важным свойством ме­ тода является возможность сохранения заданного следа сетки на границе. В отличие от двумерных треугольных сеток, в трёхмерном случае существуют такие многогранные области, для которых невозможно построить конформ­ ную тетраэдральную сетку с заданным следом на границе без добавления новых вершин. Пример такой области приведён в [32], к нему мы вернёмся в разделе 3.1.2. Б. Джо в своей работе [17] показал, что для выпуклых мно­ гогранников добавление одной вершины позволяет построить конформную согласованную с границей сетку. П. Л. Джордж показал в [14] существова­ ние конформной согласованной сетки для произвольного многогранника при добавлении конечного числа вершин. Тем не менее задача построения сетки, согласованной с заданной границей, остаётся сложной практической задачей.

При построении тетраэдральных сеток могут использоваться методы, ос­ нованные на тех же идеях, что и в двумерном случае. Классический метод построения тетраэдризации Делоне строит сетку для выпуклой оболочки то­ чек и не гарантирует сохранение границы области. Предложено несколько методов для сохранения границы: локальные модификации сетки [4], измель­ чение сетки [10], построение тетраэдризации Делоне с ограничениями (CDT – Constrained Delaunay triangulations) [35]. Алгоритм продвигаемого фронта, с другой стороны, не испытывает проблем с сохранением заданной границы, так как с неё и начинает [16]. Однако, этот алгоритм может столкнуться с такой конфигурацией фронта, продвинуть которую он не в состоянии. В ли­ тературе встречаются и гибридные методы, основанные на комбинировании идей методов продвигаемого фронта и тетраэдризации Делоне [40].

В диссертационной работе предложена комбинация двух методов. Основ­ ной используемый метод – метод продвигаемого фронта, с его помощью стро­ ится бльшая часть сетки. В качестве дополнительного метода используется упрощённая версия метода, предложенного П. Л. Джорджем в [15]. Подробнее об этих двух методах будет рассказано в Главе 3.

Разработанный и внедрённый в пакет Ani3D комплекс программ позво­ ляет автоматически надёжно конформно строить тетраэдральные сетки внут­ ри заданных многогранных областей с заданным следом сетки на границе.

Как и в двумерном случае, есть возможность контролировать желаемый раз­ мер ячеек сетки во всей области. В случае отсутствия у пользователя ин­ формации о желаемом размере, размер может определяться автоматически на основании заданной дискретной триангуляции границы области и задан­ ной пользователем скорости увеличения размера элементов при продвиже­ нии внутрь области. На первом этапе при построении большей части сетки методом продвигаемого фронта гарантируется корректная конечная работа генератора при неточных вычислениях чисел с плавающей точкой. Использо­ вание неточной арифметики уменьшает часть построенной области по срав­ нению с точной арифметикой. Тем не менее на практике больше 99% области разбивается с помощью алгоритма продвигаемого фронта. Устойчивый метод на основе тетраэдризации Делоне второго этапа может построить тетраэдры, которые из-за накопления вычислительных погрешностей могут быть вырож­ денными или вывернутыми. На третьем этапе качество ячеек сетки улучша­ ется за счёт локальных операций модификации сетки. Для этого использует­ ся библиотека aniMBA из пакета Ani3D, разработанная К. Н. Липниковым и Ю. В. Василевским. Исходный код библиотек пакета Ani3D открыт и свобод­ но распространяется.

Перейдём к краткому обзору нескольких наиболее известных комплек­ сов программ для автоматического построения неструктурированных тетра­ эдральных сеток. Рассмотрим следующие генераторы: TetGen [70], TetMesh (GHS3D) [71], NETGEN [68], CUBIT [67].

Пакет TetGen (основной разработчик – Х. Си) строит тетраэдральные сетки для произвольных многогранных областей, для этого используется ме­ тод построения и измельчения тетраэдризации Делоне. Граница области за­ даётся дискретно. При некоторых существенных ограничениях на поверх­ ностную триангуляцию гарантируется возможность построения сетки, согла­ сованной с заданной границей [12]. Желаемый размер ячеек сетки может полностью контролироваться пользователем. Программа может быть ском­ пилирована в виде исполняемого файла генератора или в виде библиотеки подключена к другим программам. Исходный код открыт и свободно распро­ страняется.

Коммерческий пакет TetMesh, ранее известный как GHS3D, основным автором которого является П. Л. Джордж, включён в ряд CAE систем, та­ ких как ANSYS и MSC.Nastran. Используемый метод основан на построении тетраэдризации Делоне. Граница области задаётся дискретно, гарантируется построение сетки, согласованной с заданной границей. Шаг сетки выбирается автоматически, также есть возможность его контролировать. Для построения начальных поверхностных триангуляций используются отдельные комплексы программ. Исходный код пакета закрыт. Предоставляются разовые лицензии для ознакомления.

Следующий пакет, NETGEN (автор – И. Шоберл), представляет собой готовый комплекс программ для автоматического построения треугольных поверхностных сеток и тетраэдральных сеток. Граница области задаётся с помощью методов конструктивной блочной геометрии, с помощью сплайнов в формате STEP и IGES, или дискретно в формате STL. Построение сет­ ки сводится к построению поверхностной сетки с помощью алгоритма про­ двигаемого фронта или метода триангуляции Делоне, а затем строится тет­ раэдральная сетка внутри области. Согласно документации бльшая часть тетраэдральной сетки строится с помощью быстрого алгоритма построения тетраэдризации Делоне, а остальная часть достраивается с помощью локаль­ ных эвристик. Автору диссертационной работы неизвестно, имеется ли воз­ можность задать фиксированную поверхностную сетку и с помощью пакета получить сетку, согласованную на границе с заданной триангуляцией, также неизвестно, имеется ли теоретическое обоснование надёжности используемых методов. Пользователь может ограничивать желаемый размер элементов гло­ бально или локально вблизи точечных источников или внутри прямоуголь­ ных областей, что даёт частичный контроль над желаемым размером элемен­ тов сетки. Имеется возможность собрать код программы в виде библиотеки.

Исходный код открыт и свободно распространяется.

Последний пакет из нашего обзора – пакет CUBIT, разрабатываемый группой разработчиков из Национальной лаборатории Сандии (США), явля­ ется готовым коммерческим продуктом для построения разных типов сеток в различных областях. Граница области может быть задана как с помощью конструктивной блочной геометрии, так и с помощью сплайнов и фасеток.

Проект активно развивает направление автоматического построения гексаэд­ ральных сеток с помощью методов отображения и заметания области. Для построения тетраэдральных сеток используется генератор TetMesh, который нами уже был рассмотрен. Библиотека для подключения к своим программам распространяется отдельно под именем CAMAL, и также является коммер­ ческим продуктом. Исходный код пакета закрыт. Предоставляются лицензии для некоммерческого использования в университетах на весь пакет CUBIT за исключением генератора тетраэдральных сеток TetMesh.

Отметим, что во всех рассмотренных примерах в качестве основного ме­ тода построения сетки используется метод построения тетраэдризации Де­ лоне. По сравнению с рассмотренными комплексами программ библиотека для построения сеток из пакета Ani3D обладает как недостатками, так и пре­ имуществами. К недостаткам можно отнести относительно медленную ско­ рость работы по сравнению с коммерческими пакетами, отсутствие полно­ ценного графического интерфейса. К преимуществам можно отнести широ­ кие возможности для задания границы области и построения треугольных поверхностных сеток, построение тетраэдральных сеток, согласованных с за­ данной граничной триангуляцией, возможность контроля желаемого шага сетки в пространстве и возможность автоматического выбора шага сетки.

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

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

Будут рассмотрены алгоритмы улучшения построенной сетки, а также про­ ведены некоторые эксперименты.

1.1. Обозначения Упорядоченную пару точек на плоскости будем называть направленным отрезком, и обозначать. Направленный отрезок можно представить как вектор с компонентами (, ) = (, ). Будем использовать сле­ дующее определение векторного произведения двух векторов:

Для трёх точек,, будем использовать следующую сокращённую запись:

[] =. Заметим, что площадь треугольника в таком слу­ чае легко выражается по формуле () = полуплоскостью относительно направленного отрезка будем называть геометрическое место точек, для которых [] > 0. Аналогично, отри­ цательной полуплоскостью будем называть геометрическое место точек, для которых [] < 0.

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

Многоугольником на плоскости назовём ограниченную открытую часть плоскости, граница которой состоит из одной или нескольких непересекаю­ щихся ломаных без самопересечений. Простым многоугольником назовём многоугольник, ограниченный одной замкнутой ломаной. На Рис. 1.1 показа­ ны несколько примеров многоугольников. Вершинами многоугольника будем называть вершины его ломаных, а сторонами многоугольника – отрезки ло­ маных.

Стороны многоугольника можно разбить на два типа: внешние и внут­ ренние. Пусть – сторона многоугольника, – середина, ( ) – -окрестность точки. Рассмотрим = ( ). Если при любом > часть многоугольника лежит по обе стороны от отрезка, то будем назы­ вать внутренней стороной или разрезом. Если при некотором > 0 часть многоугольника лежит только по одну сторону от отрезка, то бу­ дем называть внешней стороной многоугольника. В последнем случае на отрезке можно ввести направление так, чтобы оказалась в положи­ Рис. 1.1. Примеры многоугольников: (а) строго выпуклый, (б ) нестрого выпуклый, (в) невыпуклый, (г) многосвязный, (д ) многокомпонентный, (е) с внутренними разрезами.

тельной полуплоскости относительно, будем называть это направление направлением обхода многоугольника против часовой стрелки. На Рис. 1. показаны направления обхода внешних сторон нескольких многоугольников.

Заметим, что совокупность всех внешних сторон образует одну или несколько замкнутых ломаных.

Фронтом на плоскости будем называть набор направленных отрезков без самопересечений. Каждому многоугольнику можно поставить в соот­ ветствие фронт (). Для этого введём на внешних сторонах многоуголь­ ника направление обхода против часовой стрелки, а для каждой внутренней стороны поставим в соответствие пару направленных отрезков и.

Совокупность всех этих направленных отрезков и образует фронт (). Бу­ дем называть фронт замкнутым, если существует такой многоугольник, что = ().

Рис. 1.2. Обход многоугольника против часовой стрелки: (а) многосвязный многоугольник, (б ) многокомпонентный многоугольник, (в) многоугольник с внутренними разрезами.

Утверждение 1.1.1. Площадь многоугольника может быть вычислена с помощью соответствующего замкнутого фронта = () = по следующей формуле:

где – начало координат.

Доказательство. Введём векторное поле u = (, ), тогда divu = 2. И по формуле Грина где n – внешняя единичная нормаль к. Рассмотрим в один из отрезков,. Векторное поле на этом отрезке изменяется линейно, и поэтому инте­ грал по нему можно заменить на произведение значения подынтегрального выражения в центре отрезка на длину отрезка.

Суммируя по всем отрезкам фронта, получаем (1.1).

Периметром многоугольника назовём сумму длин отрезков соответству­ ющего фронта, иными словами это сумма длин внешних сторон плюс удвоен­ ная сумма длин внутренних сторон.

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

1.2. Алгоритм продвигаемого фронта Будем решать задачу построения для заданной многоугольной области конформной треугольной сетки, согласованной с границей. В этом разделе мы рассмотрим алгоритм продвигаемого фронта.

Текущим фронтом назовём некоторый замкнутый фронт, который бу­ дет являться границей многоугольной области, в которой нужно построить треугольную сетку. В начальный момент времени положим 0 =, 0 = ( 0 ), 0 =. Далее, на -ом шаге будем строить новый треугольник, добавлять его в треугольную сетку, и вычитать из области, в которой сет­ ( +1 ). Когда +1 станет пустым множеством, или что то же самое, +1 =, можно утверждать, что сетка = +1 полностью покроет.

Треугольник будем строить так, что одна из его сторон будет од­ ним из отрезков текущего фронта. Тогда после вычитания из многоуголь­ ной области нового треугольника текущий фронт изменится незначительно:

как минимум один отрезок из фронта будет удалён и не более чем два новых будут добавлены.

Пусть – некоторый направленный отрезок текущего фронта. Чтобы лежал внутри, достаточно, чтобы [] > 0 и не пере­ секал. При этом допускается, чтобы вершина была вершиной фронта, и чтобы стороны треугольника совпадали с отрезками. Будем называть треугольник, удовлетворяющий этим условиям, подходящим. Такой выбор гарантирует, что после добавления его в сетку, она сохранит свой­ ство конформности. Также при этом будет гарантироваться согласованность конечной сетки с начальной областью.

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

Так как множество конечно, то мы можем просто перебирать все тре­ угольники, до тех пор, пока очередной из них не окажется подходящим. В разделе 1.4.1 будет доказано, что хотя бы один подходящий треугольник при таком выборе найдётся.

Однако сложность такого алгоритма очень большая и на практике он бу­ дет неприменим. Рассмотрим некоторую выпуклую открытую окрестность, содержащую точки,,. Выберем из текущего фронта только те отрезки, которые имеют хотя бы одну общую точку с, получившийся фронт может оказаться незамкнутым. Назовём его локальным фронтом, и обозначим.

Любая подобласть пересекается с текущим фронтом тогда и только тогда, когда она пересекается с локальным фронтом.

Построим множество =. В силу выпуклости любой треуголь­ ник, будет лежать в. Для проверки пересечения этого треугольника с текущим фронтом достаточно проверить пересечение с ло­ кальным фронтом. В качестве области можно выбрать круг с центром в и радиуса >. Локальный перебор будет быстрее, но уже не будет гарантиро­ вать существование подходящего треугольника, поэтому надо предусмотреть возможность вернуться к полному перебору, если локальный перебор будет закончен неудачно.

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

С учётом всех предыдущих замечаний запишем алгоритм построения сетки, который назовём Алгоритм 1.1. Будем использовать следующие до­ полнительные обозначения. На каждом шаге вершины фронта обозначим как { }, а вершины построенной сетки как { }. Для некоторого > и точки через () обозначим открытую -окрестность точки. Тогда Алгоритм 1.1 запишется следующим образом:

Алгоритм 1.1 Алгоритм продвигаемого фронта в двумерном случае 1: Положить =, = ().

Выбрать с минимальной длиной.

Определить желаемую длину стороны элемента.

Определить множество кандидатов = { } ().

11:

12:

13:

14:

15:

16:

17:

18:

19:

20:

В следующих разделах мы изучим устойчивость этого алгоритма к вы­ числительным погрешностям, его конечность и его сложность.

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

1.3.1. Проверка пересечения треугольника с отрезком Одной из основных операций, чувствительных к вычислительным по­ грешностям, является операция проверки пересечения треугольника с отрез­ ком. Неправильное определение пересечения из-за вычислительных погреш­ ностей может привести к неправильной работе всего алгоритма. Отметим, что существует только две ситуации получения неправильного результата:

1. Непересекающиеся треугольник и отрезок неверно восприняты как пе­ ресекающиеся;

2. Пересекающиеся треугольник и отрезок неверно восприняты как непе­ ресекающиеся.

С точки зрения Алгоритма 1.1 ошибки первого типа не являются критич­ ными, пересечение треугольника и отрезка означает, что новый треугольник не пройдёт проверку пересечения с текущим фронтом, и соответствующая вершина-кандидат будет забракована. Ошибки второго типа в свою очередь могут привести к неверному прохождению проверки пересечения треугольни­ ка с текущим фронтом, что может привести к самопересекающемуся фронту или к неконформной сетке. Таким образом, ошибки первого типа лишь увели­ чивают перебор вершин-кандидатов, незначительно увеличивая общее время работы программы, а ошибки второго типа могут привести к неправильному результату работы программы.

Разработаем алгоритм проверки пересечения пары треугольник-отрезок, исключающий возникновение ошибок второго типа, и допускающий возникно­ вение ошибок первого типа. В общем случае два плоских выпуклых объекта не пересекаются тогда и только тогда, когда существует прямая, отделяю­ щая их друг от друга. В условиях нашего Алгоритма 1.1 также допускают­ ся общие вершины и общие рёбра, например, отрезок и треугольник могут иметь общую вершину, или отрезок может быть стороной треугольника. В этих случаях конформность сетки нарушаться не будет, поэтому их нужно рассматривать как непересекающиеся.

Пусть на входе алгоритма у нас есть треугольник и отрезок.

Будем предполагать, что [] > 0, так как только такие треугольники будут тестироваться в Алгоритме 1.1. Сначала определим, сколько общих вершин у треугольника и отрезка.

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

Эти условия можно записать в более удобной форме:

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

Если выполнено хотя бы одно из этих условий, то пересечения нет, и наоборот, если нет пересечения, то выполнено хотя бы одно из условий.

В последнем случае, когда обе вершины отрезка совпали с двумя верши­ нами треугольника, ребро является стороной треугольника, и по нашей договорённости треугольник и отрезок считаются непересекающимися.

Предположим, что у нас есть функция (), неточно вычисляющая знак выражения []. Но при этом, если () = 0, то () = sign([]), и () = 0 в спорных ситуациях. Тогда проверку на пе­ ресечение можно выполнять с помощью Алгоритма 1.2.

Алгоритм 1.2 Проверка пересечения треугольника с отрезком 1: Найти общие вершины у треугольника и отрезка если общих вершин нет, тогда если одна общая вершина, тогда 10:

11:

12:

если две общих вершины, тогда 13:

14:

15:

В остальных случаях считаем, что треугольник и отрезок пересекаются.

16:

Отметим, что допустимая неточность вычисления () может при­ вести только к ошибкам первого типа, и Алгоритм 1.2 исключает ошибки второго типа.

1.3.2. Вычисление знака определителя Операция определения знака выражения [] является основной опе­ рацией, используемой как в Алгоритме 1.1, так и в Алгоритме 1.2. Предложим способ неточного определения знака () со следующим свойством: если () = 0, то () = sign([]).

Выражение [] сводится к определителю 2 2:

При вычислении элементов матрицы мы уже вносим некоторую относитель­ ную погрешность. В дополнение к этому в выражения и опять вносится относительная погрешность, и в конце в вносится ещё одна погрешность.

Суммарную абсолютную погрешность можно оценить сверху как (|| + || + || + || + 2|| + 2||), где зависит от машинной точности. С учётом этой допустимой погрешности запишем алгоритм для определения sign([]).

Алгоритм 1.3 Вычисление () Оценить погрешность = (|| + || + || + || + 2|| + 2||).

Вычислить определитель матрицы =.

Иначе вернуть 0.

В Алгоритме 1.3 величина вычисляется с абсолютной погрешностью, не превосходящей. Поэтому если () = 0, то можно быть уверенным, что знак вычислен правильно, и () = sign([]).

На практике, если для хранения координат в памяти ЭВМ использует­ ся тип вещественных чисел с одинарной точностью, а вычисления величин,,,, в Алгоритме 1.3 производятся с двойной точностью, то машинной точности хватит для точного вычисления как всех промежуточных величин, так и последней величины. В этом случае можно положить = 0, и Алго­ ритм 1.3 будет вычислять sign([]) точно.

1.4. Конечность работы алгоритма Проведём анализ Алгоритма 1.1, и покажем конечность его работы, то есть конечность числа операций. В предложенном алгоритме два явных вло­ женных цикла и один неявный за счёт перезапуска перебора в строке 18.

Цикл перебора в строке 10 всегда конечен в силу конечности множества.

Необходимо показать, что после перезапуска перебора подходящий треуголь­ ник всегда будет найден, а также показать, что внешний цикл в строке будет конечным.

1.4.1. Существование подходящего треугольника Докажем, что поиск вершины-кандидата в строке 10 Алгоритма 1. всегда найдёт хотя бы один подходящий треугольник после перезапуска пе­ ребора в строке 18. Сначала докажем вспомогательную лемму.

Лемма 1.4.1. Пусть – многоугольник с вершинами 1 2.... 1 2 – сторона многоугольника, и угол при 1. Тогда существует такая диагональ 1, целиком лежащая внутри, что угол 2 1 <.

Доказательство. Для определённости рассмотрим случай, когда 1 2 – внеш­ няя сторона, и направление 1 2 совпадает в направлением обхода про­ тив часовой стрелки. Остальные случаи рассматриваются аналогично. Пусть + – положительная полуплоскость относительно 1 2, + = { : [1 2 ] > 0}. Рассмотрим множество вершин многоугольника, лежащих в этой по­ луплоскости, = { } +. Покажем, что не пустое. На самом деле, предположим, что оно пустое. Выпустим из середины 1 2 перпендикуляр в сторону + до пересечения с границей. Многоугольник ограничен, поэто­ му пересечение всегда существует. Пусть – ребро, которое мы пересекли.

Если =, то +, +, и, следовательно, + =, и точка пересечения с перпендикуляром не лежит в + – противоречие.

Пусть вершина = arg max 2 1. Выпустим луч 1 до пересече­ ния с границей в точке. Так как угол при вершине 1, а 2 1 <, то = 1. Пусть лежит на стороне. Хотя бы одна из вершин, лежит в, пусть. Если =, то диагональ 1 лежит внутри и угол 2 1 <, утверждение леммы верно. Иначе, – часть границы многоугольника, а 1 целиком лежит внутри (см. Рис. 1.3, а). Рассмотрим множество вершин, попавших в треугольник 1, = { }1.

не пустое, так как.

Пусть вершина = arg max 2 1, если максимальное значение угла достигается на нескольких вершинах, то выберем из них ту, которая ближе к 1. Покажем, что удовлетворяет условию леммы. Действительно, предположим, что диагональ 1 пересекается со стороной многоугольника (см. Рис. 1.3, б ). Заметим, что не пересекает и не пересекает 1, при этом точка лежит в 1. Тогда одна из точек, пусть, также лежит ближе к 1 чем, что противоречит выбору.

С помощью этой леммы можно доказать следующую теорему.

Теорема 1.1. Пусть – многоугольник с вершинами 1 2.... 1 2 – сторона многоугольника. Тогда существует такая третья вершина, что треугольник 1 2 целиком лежит внутри.

Доказательство. Если один из углов при 1 или 2, то, используя Лемму 1.4.1, мы можем так мысленно разрезать диагональю многоугольник на две части, что угол при соответствующей вершине станет меньше. Пусть соседние вершины при 1 и 2 – и соответственно. Рассмотрим область, ограниченную отрезком 1 2 и лучами 1, 2, = { : [1 ] 0, [1 2 ] > 0, [2 ] 0}. Пусть = { }. Так как стороны 1 и 2 не пересекаются, то хотя бы одна из вершин, попадёт в.

Заметим, что величина [1 2 ] с точностью до множителя и знака равна расстоянию от точки до прямой 1 2. Выберем ближайшую к прямой вершину из, = arg min [1 2 ]. Покажем, что треугольник целиком лежит внутри. Предположим, что это не так, и 1 2 пересекает сторону многоугольника (см. Рис. 1.4). Пусть для определённости лежит ближе к прямой 1 2. Но тогда [1 2 ] < [1 2 ]. При этом сторона не пересекает ни 1, ни 1 2, ни 2, то есть, что противоречит выбору.

Наличие повторного выполнения перебора в строке 18 Алгоритма 1. делает перебор в строке 10 полным. Из Теоремы 1.1 с учётом предыдущих замечаний о влиянии вычислительных погрешностей и возможности точного вычисления в Алгоритме 1.3 следует, что после повторного перебора подхо­ дящий треугольник всегда будет найден.

1.4.2. Конечность числа построенных треугольников Покажем, что внешний цикл в строке 2 Алгоритма 1.1 будет конечным, то есть количество построенных треугольников конечно. Для этого мы сде­ лаем оценку на количество треугольников через количество вершин, а ко­ личество вершин ограничим сверху за счёт конечности многоугольника и ограничения снизу на расстояние между вершинами { }.

Пусть – минимальное из попарных расстояний между точками из { }, а – минимальное из попарных расстояний между точками из { }, 0 =. При добавлении нового треугольника в возможны две ситуации: { } и =. В первом случае минимальное расстояние между вершинами сетки не уменьшится больше чем до минимального рас­ стояния между вершинами фронта, и минимальное расстояние между вер­ шинами фронта не уменьшится: +1 min(, ), +1. Рассмотрим случай =. Так как не пересекает, то расстояние от до { } не меньше расстояния от до, которое по построению не меньше, вы­ бранного в строке 6. Расстояние от до { } не меньше. Соответственно, +1 min(, ), +1 min(, ). Предложим некоторый эвристический способ выбора параметров, и. Обозначим = || – длина отрезка, а – высота в из вершины.

Предположим, что у нас есть некоторая скалярная функция (, ) :

R+, и мы хотим чтобы размер элементов в сетке менялся в соответ­ ствии с этой зависимостью. В этом случае будем говорить, что размер тре­ угольников задаётся пользовательской функцией (, ). Потребуем, чтобы функция (, ) была отделена от нуля на : (, ) > min > 0. Возьмём = max( 9, (, )). Такой выбор гарантирует существование точки и огра­ ничивает высоту треугольника снизу: 81 4 = 18.

Если соответствующая функция (, ) не задана, то будем определять желаемую длину треугольника на основе длины его основания и некоторого параметра, отвечающего за скорость увеличения размера треугольников: =, где 1 – параметр, который задаётся пользователем. Будем говорить, что шаг сетки выбирается автоматически. Так как 1, то точка всегда Возьмём параметр = 1, а параметр вычислим по следующей фор­ муле:

где > 0, 0 0, – некоторые параметры. Отметим, что при > 0 величи­ на >, а при < 0 величина > 0. Параметр можно подобрать так, чтобы при = 0 было выполнено = 0. Ограничение > на практике позволяет избегать появления резких перепадов размеров соседних отрезков в новом фронте. А ограничение > 0 позволяет ограничить снизу величину При реализации алгоритма на ЭВМ в пакете Ani2D параметры выбира­ ются одним из двух способов в зависимости от того, задаётся ли размер тре­ угольников пользователем, или используется автоматический выбор размера.

Примеры использования разных способов будут представлены в разделе 1.7.2.

Рассмотрим оба случая отдельно.

Размер треугольников задаётся пользователем. Положим = 2, 0 = 0, = 0. Тогда = 2. По построению > min, следовательно, > 2 min. Тогда Пусть теперь размер треугольников выбирается автоматически. Возьмём для любого.

Мы только что показали, что существует такое > 0, что > для лю­ бого. С помощью этого условия мы можем сделать оценку на максимальное количество вершин. Все вершины сетки { } лежат в замыкании и попар­ ное расстояние между ними больше. Пусть – количество вершин в { }.

Покроем каждую вершину кругом с радиусом 2, тогда круги не будут пересе­ каться, и покроют площадь равную 4. Расширим исходный многоуголь­ ник на, = { : dist(, ) < 2 }, см. Рис. 1.5. Тогда полностью по­ кроет круги. Оценим площадь. У каждой стороны многоугольника длиной появляется прямоугольник, который даёт прирост площади не больше чем 2. Таким образом прирост площади за счёт прямоугольников составит 2 (), где () – периметр многоугольника. Круговые сектора у вершин с внут­ ренним углом дают прирост площади не больше чем ( ) 8. Причём эта оценка верна и для внутренних углов >, в этом случае прирост становится отрицательным, но он полностью покрывается наложением двух прямоуголь­ ников, которые мы учли ранее. Из формулы на сумму углов многоугольника можно оценить вклад всех круговых секторов как 4 (), где () – коли­ чество компонент связности в. Имеем ( ) ()+ 2 ()+ 4 (). От­ Заметим, что аналогичная оценка верна и на количество вершин в фрон­ Треугольную сетку можно рассматривать как планарный граф [61]. С помо­ щью формулы Эйлера можно доказать следующее утверждение.

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

max(0, 2 5), max(2, 3 6).

Доказательство. Будем рассматривать плоскую сетку как планарный граф, для которого верна формула Эйлера:

где – количество вершин, – количество рёбер, – количество граней, вклю­ чая внешнюю, а – количество компонент связности графа. Тогда количество ячеек в сетке = 1. Пусть 3, тогда у внешней грани не меньше трёх рёбер. Остальные грани являются ячейками сетки, поэтому в них тоже не меньше трёх рёбер в каждой. Каждое ребро лежит ровно в двух гранях, по­ этого неравенства предыдущее, получаем 2 + 1, то есть 3 6, а 2 2 4. Из последнего получаем, что = 1 2 5, утверждение доказано. Отметим, что эти оценки неулучшаемы.

Из этого утверждения следует, что количество треугольников в огра­ ничено:

Эта оценка даёт ограничение на количество итераций во внешнем цикле Ал­ горитма 1.1, что в свою очередь доказывает конечность числа операций ал­ горитма.

1.5. Скорость работы алгоритма продвигаемого фронта Проведём краткий анализ скорости работы предложенного Алгоритма 1.1.

Мы уже получили оценку сверху на количество итераций внешнего цикла в строке 2 в (1.4). Обозначим через ( ) и ( ) количество отрезков в и, соответственно. Выпишем все нетривиальные операции, которые используются в Алгоритме 1.1:

1. Выбор отрезка с минимальной длиной из (строка 3);

2. Построение локального фронта (строка 7);

3. Построение списка кандидатов (строка 8);

4. Проверка пересечения треугольника с локальным фронтом (строка 11);

5. Обновление фронта (строка 13);

Операция 4 выполняется за ( ) операций проверки пересечения треуголь­ ника с отрезком (см. Алгоритм 1.2). Операция 3 может быть сделана за 2 ( ) операций проверки принадлежности вершин из окрестности ().

Сложность операций 1, 2, и 5 зависит от используемых структур дан­ ных для хранения фронта. Будем использовать для хранения фронта кучу с минимальным по длине элементом в корне и дерево поиска для быстрого по­ иска локального фронта. В этом случае операция 1 становится тривиальной, сложность операции 2 будет в среднем пропорциональна log( ( ))+ ( ).

Операция 5 состоит из трёх операций добавления или удаления отрезка из фронта, сложность каждой из них в среднем пропорциональна log( ( )).

Более подробно устройство дерева поиска обсудим в следующем разделе.

1.5.1. Быстрый поиск локального фронта В этом разделе мы рассмотрим структуру данных, позволяющую быстро находить подмножество отрезков на плоскости, пересекающихся с заданной окрестностью. Предлагаемая структура основана на классическом четверич­ ном дереве поиска. Зафиксируем некоторый ограничивающий квадрат для фронта, то есть такой квадрат, который полностью покрывает все отрезки фронта, пусть его сторона равна. Будем строить четверичное дерево, у ко­ торого у каждой вершины либо нет потомков, либо их ровно четыре. Каждой вершине дерева будем ставить в соответствие подквадрат – некоторую часть ограничивающего квадрата. А именно, корень дерева будет соответствовать полному квадрату. Для каждой вершины с четырьмя потомками разделим соответствующий ей квадрат на четыре одинаковых квадрата и поставим их в соответствие потомкам. Уровнем вершины в графе назовём количество рё­ бер в пути от этой вершины до корня дерева, так у корня дерева уровень будет равен 0, а у его непосредственных потомков – 1.

Для каждого отрезка из найдём его середину и длину. По длине отрезка найдём такое наименьшее число 0, что 2. Число бу­ дем называть уровнем отрезка. К каждой вершине дерева припишем список отрезков из, у которых уровень совпадает с уровнем вершины в графе, и середина отрезка попадает в соответствующий подквадрат. По построению < 2, то есть меньше стороны квадратов на уровне в дереве поис­ ка. Это означает, что этот отрезок фронта может пересекать только текущий квадрат и 8 его непосредственных соседей. В дереве будем хранить только те листья, в которых есть хотя бы один отрезок. Пусть у нас также имеется оценка на минимальное попарное расстояние между вершинами фронта –.

Тогда длина каждого отрезка не меньше, и максимальный уровень в дереве будет не больше log2. Количество вершин в фронте можно грубо оценить по (1.3), а количество отрезков в фронте – используя Утверждение 1.4.1.

Пусть теперь мы хотим найти локальный фронт, то есть те отрезки, кото­ рые пересекаются с окрестностью (). Для этого мы делаем обход дерева, рассматривая только подквадраты, имеющие общие точки с (), и их непо­ средственных соседей. В каждой такой вершине дерева проверяются припи­ санные ей отрезки. Рассмотрим множество квадратов, соответствующих вер­ шинам на -ом уровне дерева, которые участвуют в переборе. Эти квадраты полностью покрываются кругом, если его радиус увеличить на толщину двой­ ного слоя квадратов, то есть кругом радиуса + 2 22. В этом круге мы можем оценить количество вершин фронта по формуле, аналогичной (1.3), и из результатов Утверждения 1.4.1 получить оценку на количество отрезков:

На практике длина фронта, пересекающего круг радиуса + пропорциональна его радиусу, а из ограничения на минимальную длину от­ резков на уровне можно сделать вывод, что в среднем на уровне будет рассмотрено количество отрезков пропорциональное величине 2 +1. В сум­ ме по всем уровням количество рассмотренных рёбер составит величину по­ рядка 2 +, где – максимальный уровень в дереве, и зависит от.

Сложность такого поиска можно оценить в среднем как величину пропорци­ + log2. Количество операций, необходимые для обхода дерева, не превысят по порядку количество протестированных отрезков, так что их тоже можно включить в эту оценку.

Отметим, что в худшем случае сложность поиска локального фронта может составить величину, пропорциональную общему количеству отрезков во всём фронте, причём независимо от. Пример конфигурации фронта, отображающей возможный из “худших” случаев, приведён на Рис. 1.6.

Операция добавления отрезка в фронт по сложности пропорциональна глубине дерева. Операция поиска и удаления отрезков из фронта, сводится к поиску нужной вершины в дереве за время пропорциональное глубине дереве, а также поиску и удалению в обычном связном списке отрезков, соответству­ ющих найденной вершине. Обычно длина списка небольшая на большинстве конфигураций фронтов, однако в худшем случае может достигать величины пропорциональной количеству отрезков во всём фронте, см. Рис. 1.6.

1.5.2. Анализ сложности алгоритма Подытожим имеющиеся результаты. Пусть на входе задан многоуголь­ ник с площадью, периметром, состоящий из компонент связности, и полностью покрываемый квадратом со стороной. Если пользователем зада­ на функция желаемого размера треугольников (, ), то пусть она отделена от нуля: (, ) > min > 0.

В разделе 1.4.2 было показано, что существует параметр, с помощью которого можно оценить максимальное количество треугольников в области (1.4). Оценим сверху теперь сложность одной итерации внешнего цикла в Алгоритме 1.1. Нас будет интересовать, как растёт сложность алгоритма при уменьшении шага сетки, и, соответственно, уменьшении.

Предположим, что подходящий треугольник был найден в локальном пе­ реборе. На практике выбирается в пределах от до 2, и длины отрезков фронта в этой окрестности сопоставимы с. Оценка в среднем на количество отрезков в локальном фронте ( ) получается не зависящей от. Однако в худшем случае может составить величину пропорциональную ( ), кото­ рая в свою очередь в худшем случае достигает 2 + +. Поиск локального Далее, в переборе для каждого отрезка локального фронта требуется проверить его пересечение с локальным фронтом. Сложность этой операции пропорциональна ( )2, что в среднем является величиной, не зависящей от. Но в худшем случае сложность может достигать величины, пропорцио­ нальной ( )2.

В случае рестарта в строке 18 проводится полный перебор, который со­ ответствует худшей оценке для локального перебора.

) log2. В худшем случае оно пропорционально ( 2 + средней скорости работы будет экспериментально подтверждена в разделе 1.7.1.

Покажем пример набора областей, для которых сложность алгоритма значительно хуже средней оценки. Рассмотрим область в виде расчёски, пред­ ставленную на Рис. 1.6. Будем увеличивать в ней количество зубьев,, при этом площадь будет оставаться примерно одной и той же, периметр будет пропорционален. Параметр будет уменьшаться обратно пропорциональ­ но. Количество отрезков в начальном фронте будет 2 + 1. На каждом шаге локальная окрестность будет почти целиком накрывать всю область, поэтому сложность одной итерации внешнего цикла будет пропорциональ­ на 2. Всего будет построен 2 1 треугольник. Общая сложность работы составит величину пропорциональную 3, что значительно хуже оценки в среднем порядка log.

Вспомним конструктивное доказательство Теоремы 1.1, с помощью ко­ торого можно сразу находить подходящий треугольник и даже не проводить проверку пересечения его с фронтом. Одна итерация внешнего цикла в таком случае сводится к поиску ближайшей к отрезку вершины в положительной полуплоскости. За счёт использования деревьев поиска сложность этой опе­ рации в среднем будет пропорциональна log, и в худшем случае, где – количество вершин в текущем фронте. Если не накладывать ограничений на исходный многоугольник, то нам также может понадобится операция, описанная в Лемме 1.4.1, сложность которой в худшем случае пропорциональ­ на. Общая сложность нового алгоритма будет пропорциональна в среднем log, и 2 в худшем случае. Отметим, что минусом этой модификации ал­ горитма является невозможность контролировать размер элементов, а также низкое качество получаемых сеток.

1.5.3. Быстрая модификация алгоритма В этом разделе мы модифицируем Алгоритм 1.1, и покажем, что при определённых условиях можно гарантировать работу алгоритма без полного перебора, сделаем оценки сложности получившегося алгоритма.

Введём некоторые дополнительные обозначения, которые мы будем ис­ пользовать в этом разделе. Пусть – замкнутый фронт для многоугольника. Обозначим через min и max – минимальное и максимальное значение длин отрезков фронта. Множество вершин фронта обозначим как { }. Те­ перь рассмотрим все такие пары вершин (, ), для которых отрезок лежит внутри. Минимальное расстояние между вершинами среди всех та­ ких пар обозначим как int.

Заметим, что если отрезок пересекался с фронтом, или лежал за пределами, то при дальнейшем продвижении фронта с помощью Алгорит­ ма 1.1 этот отрезок никогда не станет отрезком фронта.

Зафиксируем некоторое > 0. Предположим, что в каждый момент времени min, max < 2, int. Покажем, что в этом случае мы можем продвигать фронт, поддерживая это свойство. В частности, в получившейся сетке стороны всех треугольников также будут лежать в пределах от до 2.

Лемма 1.5.1. Пусть задано некоторое значение > 0, и многоугольник с фронтом с вершинами { }, min, max < 2, int, отрезок – минимальный по длине отрезок фронта. Тогда можно так вы­ брать точку в положительной полуплоскости относительно, что < 2, < 2, не пересекает, и после продвижения фронта будет выполнено min, max < 2, int.

Доказательство. В силу громоздкости доказательства, приведём только ос­ новную его идею, оставляя некоторые утверждения без строгого доказатель­ ства. Из условия int следует, что внутри открытых -окрестностей вер­ шин и нет других вершин из фронта.

Сначала рассмотрим случай || 3. В этом случае можно постро­ ить такую точку в положительной полуплоскости относительно, что = =. И углы при основании треугольника будут не меньше.

Если || > 3, то построим с углами при основании. Тогда Можно показать, что при ограничениях min, max < 2, int либо пересекается с отрезком, выходящим из или, либо не пе­ ресекается с. Рассмотрим первый случай, пусть для определённости тре­ угольник пересекается с отрезком. Можно показать, что внутри нет вершин из, и не пересекает. в этом случае треугольник, у которого как минимум две стороны из, удовлетворяет условиям леммы.

Во втором случае проверим, что отрезки, лежащие в, не меньше. Если это так, то положим = и условие леммы будет выполнено, при этом у только одна сторона лежит в.

Иначе, выберем ближайшую к такую точку. Можно показать, что внутри нет вершин фронта, и не пересекается никаким отрез­ ком из фронта в силу ограничения на максимальную длину отрезков фронта.

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

Отметим, что у построенного нами треугольника либо все углы не мень­ ше, либо хотя бы две его стороны лежали в.

Модифицируем Алгоритм 1.1 с учётом Леммы 1.5.1. В строке 4 выберем в зависимости от длины ||. В строке 6 возьмём = 2, = 0, =. В строке 9 заменим проверку на проверку из Леммы 1.5.1: если { :

(), } =, то добавить к. А строчку 18 удалим из алгоритма.

Теорема 1.2. Пусть – многоугольник с фронтом, для которого min > и int > 1 max. Тогда модифицированный алгоритм построит сетку для за конечное время.

Доказательство. Возьмём = min(min, int ). Тогда в начальном фронте min, max < 2, int. Лемма 1.5.1 гарантирует существование под­ ходящего треугольника в локальном переборе, и после продвижения фронта условие будет сохраняться. На каждом шаге добавляется треугольник со сто­ ронами не меньше и меньше 2.

Покажем, что количество таких треугольников внутри ограничено.

Так как в этот раз у нас нет ограничения на минимальное попарное расстоя­ ние между вершинами сетки, то рассуждения аналогичные рассуждениям из раздела 1.4.2 здесь не работают. В доказательстве Леммы 1.5.1 мы сделали том случае, когда две его стороны уже были в фронте. А значит при добавле­ нии такого треугольника в сетку количество отрезков в фронте уменьшилось хотя бы на один. Каждый раз, когда количество отрезков в фронте увеличи­ вается на один, мы добавляем треугольник с углами не меньше. Мысленно покрасим такие треугольники в синий цвет, а треугольники с углами меньше в красный. Тогда количество красных треугольников,, не может превы­ шать количество синих,, больше чем на количество отрезков в начальном Оценим площадь одного синего треугольника. Все его углы не меньше 6, следовательно, они не больше 3. То есть синус любого его угла не меньше тогда площадь треугольника будет не меньше 1 2.

Оценим вычислительную сложность нового алгоритма. По сравнению с первоначальной версией, в этом алгоритме используется только локальный перебор, а длины отрезков в фронте ограничены сверху. Поэтому все отрезки фронта будут лежать на двух последних уровнях четверичного дерева поиска.

Из формулы (1.5) при =, = 2 получим, что ( ) + 1 < 2 < 0. То есть ( ) ограничено константной, не зависящей ни от, ни от. Таким образом выполнение одной итерации внешнего цикла занимает 1 + 2 log2 операций. А значит даже в худшем случае сложность всего Мы показали, что новая модификация Алгоритма 1.1 работает за вре­ мя пропорциональное log, где – количество треугольников. При этом длины сторон построенных треугольников лежат в пределах от до 2, где минимальная длина стороны или внутренней диагонали исходного много­ площадь многоугольника, – периметр многоугольника. К сожалению, мы не можем сделать никаких дополнительных оценок на наихудшее качество треугольников в полученной сетке.

Новый алгоритм применим только к определённому классу областей и не позволяет контролировать желаемый размер элементов при построении.

Тем не менее эта модификация имеет важный практический интерес в силу хорошей оценки на время работы даже в худшем случае. Возможность рабо­ ты в локальном фронте может быть использована, например, для построения параллельной версии алгоритма, или для использования алгоритма на кри­ волинейных поверхностях, о чём подробнее будет рассказано в Главе 2.

1.6. Улучшение качества полученной сетки В этом разделе мы обсудим способы улучшения качества треугольной сетки. Будем использовать только операцию сдвига внутренних вершин, за­ фиксировав при этом связи между вершинами, и внешние вершины, лежащие на границе области.

Для измерения качества треугольника будем использовать сле­ дующую величину [28]:

где – ориентированная площадь треугольника, = 1 [],, и – длины сторон треугольника, = 4 3. Тогда для невырожденных треугольников будет выполнено 0 < 1. Если треугольник вырожден, то = 0, а если он вывернут, то есть [] < 0, то < 0.

Пусть – множество всех вершин в сетке, – вершины лежащие на границе, неподвижные вершины, а – внутренние вершины. Для каждой внутренней вершины определим множество её соседей ( ). Количе­ ство соседей обозначим как (( )).

Простейший способ улучшения сетки заключается в сдвиге каждой внут­ ренней вершины в центр масс её соседей:

В приведённой формуле все вершины имеют одинаковый вес. После этого все вершины одновременно сдвигаются в новые посчитанные значения.

Надо следить, чтобы сетка оставалась конформной, и треугольники не выворачивались. Для этого можно проверить, что в новой сетке для каж­ дого треугольника было выполнено [] > 0. Если какой-то из треугольников вывернулся, то нужно вернуться на шаг назад, к предыдуще­ му конформному состоянию. Далее можно попробовать опять провести сдвиг вершин, но на этот раз сдвигать вершины не сразу в центр масс, а только на какую-то часть от этого расстояния: new = mass + (1 ), 0 1.

Можно каждый раз уменьшать.

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

1.6.1. Псевдоминимизация функционала качества сетки Рассмотрим следующий функционал:

где функция для треугольника зависит от его качества :

С помощью параметра можно регулировать насколько больше будет “штраф” за плохое качество треугольников. Эта функция растёт с уменьшением каче­ ства треугольников, и имеет разрыв при = 0. Немного подправим функцию так, чтобы не было разрыва [11, 47]. Для этого заменим в формуле на где некоторый параметр.

Теперь функция * всегда положительная, и функция ( ) = не будет иметь разрывов. Это позволяет использовать её в том числе и для сеток с вывернутыми треугольниками. Будем так двигать внутренние вершины, чтобы минимизировать функционал, при этом будем постепенно уменьшать значение параметра. Большие значения для нужны в первую очередь для распутывания сетки с помощью функционала, в то время как при почти нулевом функционал используется для улучшения качества треугольников.

Задача минимизации функционала является довольно сложной задачей.

Поэтому мы не будем его точно минимизировать, а будем использовать бо­ лее простой подход. Каждую внутреннюю вершину будем сдвигать отдельно, зафиксировав при этом все остальные. В этом случае ( ) = (, ), где (, ) – координаты подвижной вершины. Будем сдвигать её в направлении, противоположном градиенту функционала. Значение градиента в точке мож­ но достаточно просто вычислить аналитически.

Выпишем градиент (, ) для треугольника с координатами (, ), (1, 1 ), и (2, 2 ). В этом случае имеем:

= 22 +22 +22 21 22 21 2 +2 2 +21 +22 21 22 21 2, Для плохих и вырожденных треугольников значения * будут находить­ ся вблизи нуля, и значения * будут большие, а значит и модуль градиента функционала тоже может быть большим. Если сдвинуть вершину на большое расстояние, то могут появиться новые вырожденные треугольники и итера­ ционный процесс не будет сходиться. Чтобы избежать этого и улучшить схо­ димость, предлагается уменьшить влияние больших значений модуля функ­ ционала.

Введём функцию, уменьшающую эффект от сдвига на большие расстоя­ ния ( > 0 ) за счёт логарифмической зависимости:

Алгоритм, используемый для псевдоминимизации функционала представ­ лен в Алгоритме 1.4.

Алгоритм 1.4 Псевдоминимизации функционала качества сетки 1: Выбрать начальное значение.

Сдвинуть на расстояние в сторону, противоположную направ­ В программной реализации этого алгоритма в пакете Ani2D используют­ ся следующие параметры. = 8, = 100 – число итераций. Начальное значение = 2 /100, где – средняя длина ребра во всей области. На каждой итерации значение линейно уменьшается до значения = 2 /1000. Парамет­ ры для вычисления = · (| * |) используются следующие: = 0.252 /, 0 = 1/.

Результаты работы этого метода улучшения сетки будут показаны в раз­ деле 1.7.3.

1.7. Результаты экспериментов Приведём ряд приложений основного Алгоритма 1.1, предложенного в этой главе, и реализованного в библиотеке aniAFT из пакетов программ Ani2D и Ani3D. В первой части мы экспериментально оценим среднюю сложность работы алгоритма, во второй части приведём результаты использования раз­ ных методов задания желаемого шага сетки, в третьей части покажем резуль­ тат работы предложенного в разделе 1.6 метода улучшения сетки для области с нерегулярным фронтом.

1.7.1. Скорость работы алгоритма продвигаемого фронта Экспериментально измерим скорость работы алгоритма продвигаемого фронта на плоскости. В качестве области возьмём единичный квадрат, и бу­ дем строить на нём неструктурированные квазиравномерные сетки с шагом. Уменьшая шаг сетки, будем следить за количеством треугольников в сетке,, и временем построения сетки,. Результаты экспериментов пред­ ставлены в Таблице 1.1.

Таблица 1.1. Скорость работы алгоритма продвигаемого фронта.

Из результатов расчётов видно, что время работы пропорционально вели­ чине log. В таблицу также включены столбцы со скоростью построения (треугольники в секунду), и оценкой коэффициента перед log.

1.7.2. Различные способы выбора шага сетки В этом разделе мы проследим, как выбор локального шага сетки, или параметра, при построении нового треугольника в Алгоритме 1.1 влияет на сетку и на качество сетки.

Рассмотрим область в виде единичного квадрата с вырезанным в центре кругом с радиусом 0.1. Реализованные в пакетах Ani2D и Ani3D алгоритмы позволяют задавать желаемый размер элементов сетки с помощью скалярной функции. С помощью этой же функции можно построить дискретную грани­ цу области, заданной аналитически. Подробнее об этом будет рассказано в разделе 2.1 следующей главы.

При построении сеток использовалось два варианта скалярных функ­ ций, отвечающих за желаемый размер треугольников. Первая уменьшалась вблизи двух полуокружностей, напоминающих по форме график функции sin(). Вторая функция была тождественно равна константе на всей обла­ сти. Полученные с их помощью сетки представлены на Рис. 1.7, а и 1.7, б соответственно.

Далее протестируем автоматический выбор шага при заданной дискрет­ ной границе области. В качестве дискретной границы возьмём границу, полу­ ченную при использовании дискретизации границы с постоянной функцией желаемого шага сетки. Протестируем работу алгоритма с автоматическим увеличением шага при разных значениях = 1, 1.05, 1.25, 20. Полученные сетки представлены на Рис 1.7. Выбор большой скорости увеличения шага сетки приводит на практике к очень грубым сеткам с минимальным количе­ ством новых точек внутри области (см. Рис. 1.7, е).

В Таблице 1.2 собрана информация о количестве треугольников для по­ строенных сеток, и о минимальном качестве треугольников в получившихся сетках. Качество треугольников вычисляется по формуле (1.6). Отметим, что во всех рассмотренных примерах, кроме последнего, наихудшее качество не меньше 0.5.

Рис. 1.7. Разные способы выбора шага сетки: (а) заданная пользователем нетривиальная функция, (б ) заданная пользователем константная функция; (в)–(е) автоматический вы­ бор с = 1, 1.05, 1.25, 20, соответственно.

Таблица 1.2. Количество треугольников и наихудшее качество треугольников при разных способах выбора шага сетки.

В дополнение к рассмотренным сеткам с равномерным следом на границе проверим алгоритм для неравномерной граничной дискретизации. В качестве такой дискретизации выберем след триангуляции, полученной с помощью непостоянной функции желаемого размера сетки (см. Рис. 1.7, а). Граничная дискретизация имеет непостоянный, но плавно меняющийся шаг. Результаты автоматического выбора шага с параметрами = 1.05 и = 1.25 показаны на Рис. 1.8. Количество треугольников в получившихся сетках составило и 550 соответственно. Наихудшее качество – 0.280 и 0.645 соответственно.

Рис. 1.8. Автоматический выбор шага сетки при неравномерной дискретизации границы:

(а) = 1, 1.05, (б ) = 1, 1.25.

В следующем разделе мы рассмотрим использование автоматического выбора шага для начального фронта с сильно неравномерным шагом на гра­ нице.

1.7.3. Сложный начальный фронт Рассмотрим пример достаточно сложного и нерегулярного начального фронта. На Рис. 1.9, а показан весь начальный фронт, а на Рис. 1.9, б увеличена малая подобласть всей области. Начальный фронт предоставлен И. В. Капыриным и представляет собой границы Теченского каскада водоё­ мов и окружающих каналов. Отметим, что шаг дискретизации на границе меняется неравномерно, и вообще говоря нужна дополнительная предобра­ ботка фронта для получения высококачественной сетки. Тем не менее попро­ буем построить треугольную сетку с помощью автоматического выбора шага сетки. Возьмём = 1.25.

Рис. 1.9. Дискретная границы области. (а) вся область целиком, (б ) увеличенная часть.

Начальный фронт состоит из 2743 вершин и 5187 направленных отрез­ ков. Напомним, что внутренние разрезы порождают пару разнонаправленных отрезков. После построения сетки с помощью алгоритма продвигаемого фрон­ та была получена сетка из 14807 вершин и 29279 треугольников. Часть сетки, соответствующая Рис. 1.9, б показана на Рис. 1.10, а, на Рис. 1.10, б показана сетка после применения алгоритма улучшения качества сетки, описанного в разделе 1.6.

Рис. 1.10. Треугольная сетка. (а) после алгоритма продвигаемого фронта, (б ) после алго­ ритма улучшения качества сетки.

В Таблице 1.3 приведено распределение качества треугольников в сетке до и после улучшения сетки.

Таблица 1.3. Распределение качества треугольников до (Д) и после (П) улучшения тре­ угольной сетки.

Из таблицы следует, что качество наихудшего элемента увеличилось на порядок, и количество треугольников, для которых < 0.1, теперь равно пяти, при общем количестве треугольников около 30 тысяч. Отметим, что ка­ чество наихудших треугольников в этом примере ограничено неравномерным начальным фронтом.

1.8. Выводы к первой главе Разработан, реализован и протестирован алгоритм построения конформ­ ных согласованных с заданной границей сеток для произвольных многокомпо­ нентных многосвязных многоугольных областей. Доказана конечность числа операций и продемонстрирована устойчивость алгоритма к вычислительным погрешностям при реализации алгоритма на ЭВМ. Получены оценки в сред­ нем и в худшем случае на количество операций алгоритма продвигаемого фронта. Предложены две модификации алгоритма с лучшими оценками на сложность работы, но обладающие недостатками при их практическом ис­ пользовании. Разработан, реализован и протестирован алгоритм улучшения качества сетки на основе псевдоминимизации функционала качества сетки.

Экспериментально подтверждена оценка на сложность алгоритма продвигае­ мого фронта в среднем. Представлены некоторые примеры работы алгорит­ мов, реализованных автором в библиотеке aniAFT из пакета Ani3D.

Построение поверхностных сеток Технология построения треугольных поверхностных сеток, то есть се­ ток на границе трёхмерных областей, необходима для дальнейшего построе­ ния неструктурированных тетраэдральных сеток. Также поверхностные тре­ угольные сетки могут использоваться при решении квазидвумерных задач на криволинейных поверхностях. Примерами таких задач могут служить зада­ чи деформации тонкостенных конструкций, или решение уравнений мелкой воды на поверхности сферы или геоида.

В этой главе мы покажем, что двумерный алгоритм продвигаемого фрон­ та может быть расширен на случай криволинейных поверхностей, обсудим способы построения дискретных криволинейных границ, способы взаимодей­ ствия с системами САПР для получения необходимой информации о слож­ ных областях. Будет проведён анализ конечности и сложности алгоритмов, представлено несколько примеров работы реализованных алгоритмов.

Поверхностью назовём двумерное многообразие в трёхмерном простран­ стве. Простой поверхностью назовём такую поверхность, которая является образом гомеоморфного отображения единичного квадрата.

Кривой на поверхности или в пространстве назовём одномерное многооб­ разие на поверхности или в пространстве, соответственно. Простой кривой назовём кривую, являющуюся образом гомеорфизма единичного отрезка.

Будем по возможности использовать обозначения из раздела 1.1 и для поверхностей. В качестве векторного произведения двух направленных от­ резков на поверхности будем использовать следующую величину: [] =, n, где n – внешняя единичная нормаль к поверхности в точке, а (·, ·) – скалярное умножение. Отметим, что для плоскости новое обозначе­ ние в точности совпадёт с определением из Главы 1. Также можно отметить, что такая величина соответствует значению [ ] в касательной плоскости, где и – проекции и соответственно.

2.1. Представление поверхности Сначала рассмотрим простую поверхность, которую можно представить, как двумерное многообразие – образ гомеоморфного отображения единично­ го квадрата. То есть должна быть определена такая параметризующая функ­ ция, которая переводит точку из квадрата в параметрическом пространстве в точку на поверхности.

Построение такого гомеоморфизма для произвольной поверхности явля­ ется трудной задачей. На практике можно построить гомеоморфизм из неко­ торой ограниченной области в параметрическом пространстве R2.

Разбиение границы при гомеоморфизме порождает разбиение границы простой поверхности, которое можно рассматривать как дискретную границу поверхности. Используя эту границу как фронт, можно применить алгоритм продвигаемого фронта, локально рассматривая поверхность как плоскость.

Сложная поверхность, как правило, может быть разрезана на несколько простых поверхностей. Зафиксировав дискретизацию на разрезах, можно для каждой простой поверхности построить согласованную с этой дискретизацией триангуляцию, тогда общая треугольная сетка для всей сложной поверхности будет конформной.

Составим алгоритм действий на основе описанной выше идеи. Вся по­ верхность разбивается на несколько простых поверхностей, будем назы­ вать их криволинейными гранями. Каждая криволинейная грань параметри­ зуется гладкой функцией.

Под гладкостью функции будем подразумевать существование и непрерыв­ ность производных по, : 1 ( ).

Границу криволинейной грани разобьём на несколько простых кривых,, которые будем называть криволинейными рёбрами. Введём параметриза­ цию на криволинейных рёбрах.

Нам также понадобится отображения, из параметрического про­ странства кривой в параметрическое пространство поверхности:

Отметим, что для произвольных, вычисление этих отображений мо­ жет оказаться сложной задачей. Поэтому на практике такой подход приме­ ним для достаточно простых областей.

При реализации этого метода на ЭВМ в пакете программ Ani3D был использован другой подход. Пусть у нас задана параметризация криволиней­ ных рёбер в параметрическом пространстве поверхности отдельно для каж­ дой криволинейной грани. Пусть грань параметризована функцией (, ), ребро можно рассматривать как кривую в параметрическом пространстве (, ). Например, если мы используем отображение (, ()), тогда в R3 кривая будет параметризована следующим образом:

Такие параметризации кривой могут быть разными для разных криво­ линейных граней. На практике можно допускать незначительное отклонение в пределах допустимой погрешности, что упрощает построение параметриза­ ций и расширяет класс областей, для которых такой подход применим.

Для построения поверхностных сеток на криволинейных гранях сначала строится дискретизация криволинейных рёбер. Рассмотрим параметризован­ ную простую кривую:

Будем строить на ней точки:

Значения параметров [0, 1 ] выбираются с помощью метода бисекции в соответствии с желаемым расстоянием между точками.

После построения дискретизации границы криволинейной грани, мы мо­ жем применить аналог алгоритма продвигаемого фронта для построения тре­ угольной сетки. Так как соседние грани имеют одну и ту же дискретизацию их общего ребра, и построенные на них треугольные сетки согласованы с этой дискретизацией, то составленная из них общая сетка будет конформной.

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

На первом этапе построим дискретизацию для всех криволинейных рё­ бер. Для этого на кривой линии расставляются точки, и кривая приближается ломаной. Длина отрезков регулируется желаемым размером элементов сетки.

На втором этапе для каждой криволинейной грани составляется дискрет­ ная граница из имеющихся рёбер. Для точек на рёбрах восстанавливаются значения параметров (, ) в параметрическом пространстве грани. С помо­ щью алгоритма продвигаемого фронта строится поверхностная триангуля­ ция, согласованная с дискретной границей. Запишем Алгоритм 2.1.

Алгоритм 2.1 Построение поверхностной сетки 1: для всех криволинейных рёбер: начало цикла Выбрать одну из параметризаций простой кривой.

Построить дискретизацию кривой с помощью метода бисекций.

для всех криволинейных граней: начало цикла Составить дискретизацию границы из дискретизаций рёбер.

Применить метод продвигаемого фронта для построения сетки.

Объединить все построенные сетки в одну общую сетку.

Согласованность поверхностных сеток граней с дискретизациями рёбер гарантирует сохранение конформности общей сетки.

2.2. Взаимодействие с геометрическим ядром САПР Информация о границе области может быть получена с помощью гео­ метрического ядра САПР [42]. Большинство из существующих САПР пред­ лагают интерфейс для взаимодействия со своим внутренним геометрическим ядром, позволяя получать информацию о топологии и геометрии области. У каждой системы САПР свой интерфейс для взаимодействия, и общих стан­ дартов в этой области пока нет. Некоторые САПР с открытым исходным кодом хорошо документированы. Одним из примеров такой САПР может слу­ жить открытая система Open CASCADE Technology [69]. Однако, в большин­ стве случаев закрытость исходного кода и отсутствие документации по ин­ терфейсам в свободном доступе усложняет разработку интерфейсов с САПР.

Перспективным направлением видится использование промежуточных библиотек, предоставляющих общий унифицированный интерфейс для раз­ работчика и поддерживающих взаимодействие с разными САПР. Во время разработки коммерческого пакета CUBIT Tool Suite, включающего в себя в том числе и собственное геометрическое ядро ACIS, разработчики стали до­ бавлять новые интерфейсы к другим САПР. Для этого была создана специ­ альная прослойка, предоставляющая общий интерфейс для взаимодействия с САПР. Позднее эта часть кода была вынесена из коммерческого проекта в виде отдельной открытой библиотеки Common Geometry Module (CGM) [72].

На тот момент были разработаны интерфейсы к геометрическим ядрам си­ стем ACIS и Pro/ENGINEER. Открытость исходного кода позволила новым разработчикам использовать и совершенствовать эту библиотеку. В 2006 го­ ду автор диссертационной работы приступил с созданию интерфейса между геометрическим ядром Open CASCADE и библиотекой CGM. Базовая вер­ сия этого интерфейса позволила использовать CGM вместе с САПР Open CASCADE в пакете Ani3D. В настоящее время новой группой разработчи­ ков ведётся дальнейшее развитие проекта CGM под новым именем, CGMA [73]. В этой версии основной упор делается на улучшение взаимодействия с геометрическим ядром Open CASCADE.

CGMA предлагает универсальный интерфейс для общения с геометри­ ческими ядрами САПР. При этом вся информация делится на две части:

Топологическая информация о модели: составные части модели и их топологические отношения между собой.

Геометрическая информация: координаты, размеры, параметры, и па­ раметризующие функции для кривых и поверхностей.

В большинстве ядер САПР используется граничное представление моделей (boundary representation, B-Rep или BREP). Такой метод представления ис­ пользуется и в библиотеке CGMA. Модель состоит из частей, образующих древовидную иерархию. Рассмотрим их от самых простых к сложным.

1. Точка, для которой заданы её координаты в пространстве (,, ).

2. Ребро – часть гладкой параметризованной кривой, ограниченная точ­ ками. Для ребра определена некоторая параметризация (,, ), а также ограничивающие её точки и значения параметра в этих точках.

3. Петля – связный и замкнутый набор рёбер. Не несёт никакой геометри­ ческой информации и является чисто топологическим объектом.

4. Грань – часть гладкой параметризованной поверхности, ограниченная петлёй. Рёбра петли должны лежать на поверхности грани. Грань зада­ на параметризацией поверхности (, ) (,, ).

5. Оболочка – связный и замкнутый набор граней. Также, как и петля, является чисто топологическим объектом.

6. Тело – часть пространства, ограниченная оболочкой.

7. Набор тел, представляющих модель в целом.

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

Особое внимание стоит уделить криволинейным поверхностям. Геомет­ рические ядра САПР часто создают грани с периодичной параметризаци­ ей, а также с параметризацией, имеющей особые точки. Например, Open CASCADE задаёт боковую поверхность цилиндра одной гранью с периодич­ ной параметризацией, а боковую поверхность конуса – гранью с особой точ­ кой в вершине конуса. Более того, поверхность шара задаётся одной гранью с периодичной параметризацией и двумя особыми точками в полюсах, сама грань при этом ограничивается лишь одним ребром, соединяющим два полю­ са. Алгоритм продвигаемого фронта должен быть соответствующим образом модифицирован, чтобы допускать такие особенности. Примеры работы алго­ ритма с геометрическими моделями САПР будут приведены в разделе 2.4.2.

2.3. Алгоритм продвигаемого фронта Алгоритм 1.1 продвигаемого фронта естественно расширяется на случай криволинейных поверхностей. В этом разделе мы лишь отметим основные особенности поверхностного алгоритма продвигаемого фронта.