«Технология построения адаптируемых многогранных сеток и численное решение эллиптических уравнений 2-го порядка в трехмерных областях и на поверхностях ...»

Федеральное государственное бюджетное учреждение наук

и

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

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

Чернышенко Алексей Юрьевич

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

адаптируемых многогранных сеток

и численное решение эллиптических

уравнений 2-го порядка

в трехмерных областях и на поверхностях

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

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

ДИССЕРТАЦИЯ

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

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

Василевский Юрий Викторович Москва – Содержание Введение................................... Обзор используемой терминологии.................. Глава 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. Выводы к второй главе....................... Глава 3. Метод приближенного решения эллиптических урав­ нений 2-го порядка на поверхностях............... 3.1. Предварительные сведения и обозначения............ 3.2. Расширение уравнения на поверхности.............. 3.3. Численные методы......................... 3.4. Численные эксперименты...................... 3.5. Выводы к третьей главе...................... Заключение................................. Литература................................. Введение Процесс решения задач математической физики с использованием ЭВМ можно разделить на три основных этапа: построение расчетной сетки, дис­ кретизация дифференциальных или интегральных уравнений и решение си­ стемы алгебраических уравнений. В настоящей работе центральную часть занимает проблема построения качественных расчетных сеток. Кроме этого, будут рассмотрены дискретизации эллиптических уравнений в трехмерных областях и на поверхностях. Различные методы дискретизаций представлены в работах [4, 5, 7]. Методы решения алгебраических уравнений рассматрива­ ются в работах [6, 8, 9, 82].

В трехмерном пространстве среди прочих выделяются три класса сеток:

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

Так, например, сетки должны приближать границу области с достаточ­ ным порядком точности. В настоящее время при расчетах приемлемым яв­ ляется второй порядок точности. В консервативных дискретизациях, попу­ лярных у инженеров, степени свободы находятся в центрах ячеек, поэтому генераторы сеток должны стремиться минимизировать число ячеек для за­ данной плотности распределения узлов сетки. При фиксированном количе­ стве вершин сетки известны следующие оценки количества ячеек и граней в сетках разных типов. Для неструктурированных тетраэдраль­ ных сеток 5.5, 11, для треугольных призматических сеток точки зрения, гексаэдральные сетки являются наиболее выгодными. Полу­ чающиеся при дискретизации на гексаэдральных сетках шаблоны являются самыми компактными, а, значит, соответствующие матрицы имеют меньшее количество ненулевых элементов. Такие матрицы являются более экономич­ ными с точки зрения требуемого объема памяти и вычислительной сложности решения системы уравнений. С другой стороны, с точки зрения автоматиче­ ского построения сеток для сложных областей, гексаэдральные сетки явля­ ются самыми трудными и неудобными. Существующие на сегодняшний день методы построения гексаэдральных сеток чаще всего применимы для узко­ го класса областей (блочно-структурированные сетки). Кроме этого остается открытым вопрос качественной аппроксимации границы подобными сетками.

Однако, можно отметить некоторые коммерческие генераторы, позволяющие строить качественные гексаэдральные сетки в сложных областях: генератор HEXPRESS компании Numeca [15], а также генератор Hexotic [20], являю­ щийся частью комплекса MeshGems [21].

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

Другим недостатком является то, что для анизотропных областей получае­ мые неструктурированные сетки могут иметь ячейки с очень большими дву­ гранными углами, что значительно ухудшает качество дискретизации. Одна­ ко, с помощью тетраэдральных сеток можно строить расчётные сетки в сколь угодно сложных областях [2, 39]. Существует ряд комплексов программ, поз­ воляющих строить тетраэдральные сетки: открытые Ani3D [17], TetGen [24], NETGEN [23], а также коммерческие TetMesh [25], CUBIT [18] и другие.

В некоторых задачах призматические сетки являются компромиссом меж­ ду тетраэдральными и гексаэдральными сетками. Они не так дороги с вычис­ лительной точки зрения и могут быть эффективными в анизотропных об­ ластях. Они широко применяются в геофизических приложениях в задачах со слоистыми областями. В таких областях призматическую сетку можно представить как тензорное произведение неструктурированной треугольной сетки в плоскости и одномерной сетки вдоль оси. В Институте вы­ числительной математики РАН генератор треугольных призматических се­ ток разработан В.Н.Чугуновым [14]. Известны также коммерческие пакеты FEFLOW [19], MODFLOW-USG [22], а также разрабатываемый расчетный комплекс GeRa [1], позволяющие строить такие сетки.

Перспективным направлением в развитии сеточных генераторов являет­ ся создание надёжной технологии построения гибридных сеток для сложных областей, состоящих преимущественно из гексаэдров, а также из пригранич­ ных многогранных ячеек. В настоящее время можно выделить несколько па­ кетов программ, которые позволяют строить сетки с преимущественно гек­ саэдральными ячейками. Закрытый программный пакет ЛОГОС.ПреПост, разрабатываемый в ВНИИЭФ г.Сарова, позволяет строить сетки с шести­ гранными и многогранными ячейками. Кроме этого, известен коммерческий пакет HEXPRESS/Hybrid компании Numeca [15], позволяющий строить по­ добные сетки.

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

Geomigration of Radionuclides - совместный проект ИВМ РАН и ИБРАЭ РАН в рамках проекта "Прорыв" ГК Росатом Сколотая ячейка представляет собой часть кубической ячейки, получен­ ную в результате ее среза поверхностной сеткой. Скалывание кубической ячейки может осуществляться различными способами. Так, например, Бре­ тоннет и др. [34] применяют алгоритм полигональных сечений (Polygon clipping algorithm) [87], который используется в коммерческом генераторе HEXPRESS [15]. В настоящей работе в алгоритме скалывания используются поверхност­ ные триангуляции. Существует ряд методов построения поверхностной три­ ангуляции, использующих кубические сетки.

Широко известен метод марширующих кубов [71] (Marching cubes, MC).

Однако, при всей своей популярности, известны также проблемы этого мето­ да. Во-первых, в силу неоднозначности выбора триангуляции в кубической ячейке, алгоритм может порождать топологически несвязные сетки. Суще­ ствует ряд работ, посвященных решению этой проблемы, [36, 53, 64, 73] и др.

Во-вторых, классический метод марширующих кубов не может воспроизво­ дить резкие искривления и особенности границы. В своей работе [60] Коб­ бельт и др. предлагают расширенный метод марширующих кубов (Extended marching cubes), который отслеживает сложную геометрию границы, исполь­ зуя дополнительную информацию о нормалях к поверхности. В-третьих, ал­ горитм MC может быть применён только к равномерным кубическим сеткам.

В случае, если размеры двух соседних ячеек отличаются, получаемая триан­ гуляция становится топологически несвязной и может содержать дырки и пустоты. В работах [54, 56, 70, 91] предложены методы для восстановления топологической связности такой триангуляции на сетках типа восьмеричное дерево. Однако, данные алгоритмы также имеют существенные недостатки и сложности реализации, описанные в работе [55]. Там же предложен эле­ гантный алгоритм кубических марширующих квадратов (Cubical marching squares, CMS), который лишен всех описанных недостатков. Данный алго­ ритм позволяет строить конформную триангуляцию поверхности для сеток типа восьмеричное дерево.

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

Во многих прикладных задачах, таких как построение сеток для числен­ ных моделей подземной гидродинамики или построение трехмерных сеток по данным МРТ или КТ, область разбита на непересекающиеся подобласти с различными физическими свойствами, которые граничат между собой про­ извольным образом. Построенная сетка должна правильно отражать нали­ чие различных подобластей. Для построения подобных сеток в настоящей работе предложена модификация метода марширующих кубов для областей с нескольким материалами (Multiple material marching cubes, M3 C [93]) с более точным приближением границы. Эта модификация обобщает преимущества методов CMS и M3 C, а, значит, может быть использована для построения сколотых ячеек на сетках типа восьмеричное дерево в областях с нескольки­ ми материалами. В результате получается многогранная сетка типа восьме­ ричное дерево со сколотыми ячейками, которая состоит преимущественно из гексаэдров.

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

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

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

Классический метод конечных объёмов со степенями свободы в центрах ячеек и линейной двухточечной аппроксимацией диффузионного потока удо­ влетворяет ДПМ [31, 86]. Однако, главные направления тензора диффузии должны быть ортогональны граням сетки, в противном случае метод не име­ ет даже первого порядка точности для анизотропных диффузионных задач или на неструктурированных сетках. Тем не менее именно этот метод являет­ ся наиболее распространённым в моделировании течений в пористых средах в силу своей технологической простоты и монотонности. Многоточечная ап­ проксимации потока (MPFA – Multipoint flux approximation) имеет второй по­ рядок точности благодаря использованию большего числа точек в шаблоне, однако является условно устойчивой и условно монотонной. Ограничения на устойчивость и монотонность конечно-объёмных схем MPFA рассмотрены в [26, 59, 75].

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

Этот подход развивался в работах [40, 66, 67, 74, 95] (см. также ссылки в этих работах), в которых доказывается неотрицательность решения на про­ извольных многогранных сетках и тензорах общего вида. Некоторые из по­ добных схем вводят дополнительные неизвестные на вершинах, ребрах или гранях сетки, значения в которые интерполируются из концентраций в со­ седних ячейках. Выбор схемы интерполяции оказывает большое влияние на точность нелинейной схемы [65, 95]. Определённый метод интерполяции мо­ жет оказаться эффективным для одних задач и неэффективным для других.

В работах [66, 67] предлагаются и исследуются схемы для уравнений диф­ фузии и конвекции-диффузии, не требующие интерполяции решения в узлы сетки.

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

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

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

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

Помимо дифференциальных уравнений в трехмерных областях, во мно­ гих прикладных задачах возникают уравнения в частных производных на поверхностях. Такие уравнения возникают в математических моделях раз­ личных природных явлений: диффузия вдоль межкристаллических границ [72], перенос поверхностно активных веществ на интерфейсах многофазных течений [81], липидные взаимодействия в биомембранах [47], а также во мно­ гих инженерных и биомедицинских приложениях: визуализация векторного поля [43], синтез текстур [89], построение развертки мозга [88], моделирова­ ние жидкости в легких [52] и многих других. Таким образом, в настоящее время большой интерес представляет развитие методов численного решения уравнений на поверхностях.

Можно выделить несколько основных подходов к численному решению уравнений на поверхностях. Один из них требует построения явных триангу­ ляций поверхности и дискретизации уравнений на получаемых сетках. Разви­ тие численных методов, основанных на поверхностных триангуляциях, нача­ лось с работы [45]. В этом классе методов, поверхность аппроксимируется се­ мейством последовательных регулярных триангуляций. Предполагается, что все вершины триангуляций лежат на поверхности. В работе [46] метод из [45] был скомбинирован с Лагранжевым методом отслеживания поверхности и был обобщен для уравнений на движущихся поверхностях. Однако, исполь­ зование поверхностных триангуляций имеет ряд существенных недостатков, в частности: требуется явная параметризация поверхности, получаются доста­ точно сложные дискретизации, если поверхность эволюционирует, то триангу­ ляции требуют перестроения, и некоторые другие. Отметим, что для постро­ ения поверхностных триангуляций, вообще говоря, может быть использован метод, описанный в первой главе настоящей работы. Однако качество полу­ чаемых треугольников не является приемлемым для дискретизаций на таких сетках. В статье [78] был предложен подход для решения уравнений заданных на поверхностях, который не требует параметризации и явной триангуляции поверхности. Метод основан на следах функций из внешнего конечно-эле­ ментного пространства на дискретной поверхности, которая может являться нулем кусочно-линейной непрерывной функции уровня и произвольно накла­ дываться на внешнюю сетку. Данный метод получил развитие в [80] и для эволюционирующих поверхностей в [77, 79]. Отметим, однако, что матрицы алгебраических систем, возникающие в данном подходе могут быть плохо обусловленными (см. [76]) и реализация метода требует численного интегри­ рования вдоль дискретной поверхности. Другой метод, избегающий исполь­ зования триангуляций поверхности, был предложен в [32]. Он предполагает продолжение дифференциального уравнения на поверхности во множество из R3 положительной лебеговой меры. В результате получается формулиров­ ка уравнения в пространстве большей размерности, однако уравнение может быть решено на сетке, не привязанной к поверхности. В случае движущейся поверхности, этот метод позволяет избежать Лагранжева описания эволюции поверхности, используя Эйлеров подход [94].

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

Более подробно плюсы и минусы этого метода описаны в работах [41, 51].

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

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

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

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

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

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

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

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

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

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

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

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

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

1. Предложен алгоритм и разработана технология надежного построения многогранных сеток типа восьмеричное дерево со сколотыми ячейками для сложных областей с несколькими материалами.

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

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

Апробация работы. Результаты диссертационной работы докладыва­ лись автором и обсуждались на научных семинарах Института вычислитель­ ной математики РАН, Института прикладной математики РАН им. М. В. Кел­ дыша, Вычислительного центра РАН им. А. А. Дородницына, Института про­ блем безопасного развития атомной энергетики РАН и на следующих науч­ ных конференциях: конференция молодых ученых “Технологии высокопроиз­ водительных вычислений и компьютерного моделирования” (СПбГУ ИТМО, С.-Петербург, апрель 2009); конференция “Лобачевские чтения” (КГУ, Ка­ зань, ноябрь 2009); конференции “Тихоновские чтения” (МГУ, октябрь 2012);

конференции “Актуальные проблемы прикладной математики и механики” (Абрау-Дюрсо, сентябрь 2012); международная конференция “NUMGRID-2012” (ВЦ РАН, Москва, июнь 2012); международная конференция “CRC-NAA-2013” (Ростов-на-Дону, июнь 2013); международная конференция “Mathematical modeling of natural disasters and technical hazards” (Сьон, Швейцария, август 2013).

Публикации. Основные материалы диссертации опубликованы в 6 пе­ чатных работах, из них 2 статьи в рецензируемых журналах, входящих в перечень ВАК [13, 38], и 4 - в сборниках тезисов конференций [10–12, 37].

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

Структура и объем диссертации. Диссертационная работа состоит из введения, обзора используемой терминологии, трех глав, заключения и списка литературы из 97 наименований. Диссертационная работа содержит 31 рисунок и 11 таблиц. Общий объем диссертационной работы – 125 страниц.

Благодарности Автор диссертационной работы выражает огромную благодарность и глубокую признательность научному руководителю Ю. В. Василевскому за продолжительную поддержку, ценные советы и плодотворное обсуждение во­ просов. Автор выражает глубокую признательность М. А. Ольшанскому за совместную работу, помощь и ценные советы. Автор также выражает бла­ годарность А. А. Данилову, К. Д. Никитину, И. В. Капырину, К. М. Терехову и многим другим за активную помощь в обсуждении идей и методов, используе­ мых в диссертационной работе. Автор благодарен своим родителям и друзьям за поддержку.

Работы над диссертацией была частично поддержана грантами РФФИ 11-01-00971, 12-01-33084, 12-01-31223, федеральной целевой программой “На­ учные и научно-педагогические кадры инновационной России”, грантом Upstream Research Center of ExxonMobil corp. и проектом “Прорыв” ГК “Росатом”.

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

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

Восьмеричное дерево (или восьмидерево) представляет собой древовидную структуру данных, элементы которой удовлетворяют следующим свойствам:

Каждый элемент имеет либо ровно восемь потомков, либо не имеет по­ томков. Элемент, не имеющий потомков, называется листом, или ли­ стовым элементом;

Элемент не может быть потомком самого себя;

Каждый элемент, кроме исходного, называемого корневым, является потомком только одного элемента, называемого родителем;

Корневой элемент не имеет родителя, всегда существует и единственен.

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

Рис. 1. Восьмеричное дерево (а). Пример сетки типа восьмеричное дерево (б ).

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

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

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

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

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

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

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

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

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

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

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

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

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

1.1. Построение сколотых ячеек В этом разделе опишем алгоритм кубических марширующих квадратов (CMS) и процесс создания сколотых ячеек на его основе.

1.1.1. Алгоритм кубических марширующих квадратов (CMS) Входными данными для алгоритма является кубическая сетка, верши­ нам которой приписан знак характеристической функции области. На реб­ рах сетки отмечены точные точки пересечения с границей области, а также нормали к границе в этих точках. При этом на каждом ребре должно быть от­ мечено не более одной точки пересечения. В случае, если на ребре возникает более одной точки пересечения, то все ячейки, разделяющие это ребро, сле­ дует измельчить. Как показано в [60], такой тип данных может быть получен из большинства используемых типов входных данных.

Алгоритм CMS основывается на следующих идеях:

марширующие кубы могут быть развернуты в марширующие квадраты;

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

Как показано на Рис. 1.2, б ), куб может быть развернут в шесть гра­ ней. Для каждой грани строится кусочно-линейный контур пересечения с поверхностью границы области с помощью алгоритма марширующих квад­ ратов (MS, Marching squares). Здесь и далее контуром будем называть набор ломаных. Элементами контура будем называть отрезки ломаных. Если грани свернуть обратно в куб, то контуры будут образовывать замкнутые компонен­ ты, которые триангулируются каким-либо способом. Стоит отметить, что в случае равномерной кубической сетки, полученная триангуляция будет по­ хожа на триангуляцию, построенную классическим методом марширующих кубов.

Рассмотрим процесс построения кусочно-линейного контура поверхности на листовой грани. Каждая грань имеет 4 вершины, которым приписан один из двух знаков характеристической функции. Таким образом, всего существу­ ет 24 = 16 различных вариантов распределения знаков, которые с учетом симметрии и поворотов могут быть сведены к 5 вариантам, изображенным на Рис. 1.2, а). В случаях 4-5 возникает неоднозначность выбора пары то­ чек, которые необходимо соединить контурами. Эту неоднозначность можно решать различными способами. Классический подход предполагает заранее выбрать приоритет связности одного из знаков характеристической функции.

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

Помимо этого, нормали в алгоритме CMS можно использовать для то­ го, чтобы отслеживать особенности поверхности. Рассмотрим грань с двумя точками пересечения на ребрах. Если угол между проекциями нормалей на плоскость грани в этих точках меньше определенного порога, то контуром на данной грани будет отрезок, соединяющий эти точки. В противном случае считается, что поверхность имеет особенность. Тогда из точек на ребрах про­ Рис. 1.2. Марширующие квадраты (а). Развертка куба (б ). Рисунки из [55].

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

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

После того, как на всех гранях куба построены контуры, происходит объединение всех контуров. Затем из них выделяются компоненты, которые образуют циклы. Для этого достаточно сделать обход по ломаным, начав с произвольной вершины и двигаясь к соседней. После этого на всех замкнутых компонентах необходимо построить триангуляцию. Для этого можно исполь­ зовать метод “разделяй и властвуй” (divide-and-conquer triangulation) [83].

Суть метода заключается в следующем. Рассмотрим замкнутую последова­ тельность вершин, соединенных отрезками. Она рекурсивно делится на две Рис. 1.3. Использование нормалей: (а) разрешение неоднозначности, (б ) отслеживание искривлений границы.

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

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

Обычно последовательность вершин может быть разделена нескольки­ ми способами. В таком случае, можно выбирать наилучшую разделяющую прямую по некоторому критерию. Так, например, в [83] предлагается рас­ сматривать отношение минимального расстояния от вершин последователь­ ности до разделяющей плоскости к длине разделяющего отрезка (отрезка разделяющей прямой, ограниченного соответствующими двумя точками по­ следовательности). Тогда выбирается разделяющая прямая, для которой это отношение является наибольшим.

Триангуляция “разделяй и властвуй” может быть описана в виде алго­ ритма 1.1. Для простоты в данном алгоритме будем выбирать первую под­ ходящую разделяющую прямую.

Алгоритм 1.1 Триангуляция “разделяй и властвуй” Присвоить количество вершин в Создать треугольник из вершин Вычислить среднюю нормаль n Положить - диагональ, соединяющая и вершины 10:

11:

12:

если and лежат в разных полупространствах относи­ 13:

14:

15:

16:

17:

18:

Наличие нормалей в точках на ребрах позволяет также отслеживать осо­ бенности поверхности внутри ячейки. Для поиска особенности поверхности используется идея, предложенная в [60]. Обозначим точки на ребрах за s, единичные нормали в этих точках за n. Оценим открытый угол конуса, на­ тянутого на n, по следующей формуле:

Если меньше некоторой пороговой величины, то считается, что поверх­ ность имеет особенность. Пусть n0 и n1 - нормали, образующие наибольший угол и пусть n* = n0 n1. Тогда оценим максимальное отклонение нормалей n от плоскости, натянутой на n0 и n1. Другими словами, если величина больше, чем некоторая пороговая величина, то считается, что поверхность имеет особенную точку в ячейке. На практике авторы из [60] используют следующие значения: = 0.9, = 0.7. Искомая точка p находится как решение системы:

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

В случае, если имеется особенная точка p внутри ячейки, то триангуля­ ция строится методом “веер треугольников” (triangle fan), с центром в точке p. Если точка p не попадает внутрь ячейки, то считаем, что ячейка не имеет особенной точки.

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

Метод поиска особенной точки также позволяет разрешать неоднознач­ ность триангуляции кубической ячейки. Такие неоднозначности могут воз­ никать, когда невозможно определить, нужно ли соединять две замкнутые Рис. 1.4. Неоднозначность триангуляции. Желтые точки - на ребрах, красные - на гра­ нях. (а) Замкнутые компоненты контура, (б ) триангуляция объединения контуров, (в) триангуляция каждого контура независимо.

компоненты контура или нет, см. Рис. 1.4(а). Разрешить такие ситуации в [55] предлагают способом, аналогичным для случая неоднозначности на гра­ ни. Для этого, в каждой из компонент находится особенная точка описанным выше способом. Затем в них строится конусообразная поверхность с центром в этой точке. Если полученные поверхности пересекаются, значит, компонен­ ты должны быть связными. Тогда триангуляция, соединяющая эти компо­ ненты, топологически будет цилиндрической поверхностью, см. Рис. 1.4(б ).

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

Для произвольной кубической ячейки здесь и далее будем обозначать множества ее граней, ребер и вершин через, и, соответственно.

Тогда метод CMS для ячейки может быть представлен в виде следующего алгоритма 1.2:

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

11:

12:

Триангуляция() 13:

14:

В строке 10 имеется в виду неоднозначность, описанная ранее, см. Рис. 1.4.

В строке 13 может быть использована триангуляция “разделяй и властвуй” либо “веер треугольников”, если имеется особенная точка в ячейке.

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

Затем на каждой из граней ячейки строится соответствующий контур. Объ­ единив полученные контуры в замкнутую ломаную, строится ее триангуля­ ция (Рис. 1.5, в). Эта триангуляция делит ячейку на части, причем те части, которые не принадлежат области, удаляются. В результате получается ско­ лотая ячейка (Рис. 1.5, г).

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

Однако, плоские двугранные углы отслеживаются точно.

1.2. Случай составной области Традиционный метод марширующих кубов, а также упомянутые ранее модификации, могут быть применены только в однородных областях. На практике существует множество задач, в которых область разбита на непере­ секающиеся подобласти с различными физическими свойствами. Подобласти могут располагаться произвольным образом, и, в таком случае, одна куби­ ческая ячейка может принадлежать нескольким подобластям. Если приме­ нить алгоритм MC или CMS для каждой подобласти итерационно, то полу­ ченные треугольные интерфейсы будут образовывать пустоты, и сетка по­ лучится топологически несвязной. Дотреппе и др. в своей работе [44] пред­ лагают модификацию метода марширующих тетраэдров на случай области с несколькими материалами. В работе [97] предложен алгоритм построения тет­ раэдральных и гексаэдральных сеток для составных областей, использующий метод дуальных контуров (Dual Contouring [57]). В работе [93] Ву и Салли­ ван предложили модификацию метода марширующих кубов для областей с несколькими материалами. Этот алгоритм позволяет создать топологически связную конформную триангуляцию на равномерной кубической сетке для составных областей. Он использует двумерный метод марширующих квад­ ратов для областей с несколькими материалами (Multiple material marching squares, M3 S). Далее опишем методы марширующих квадратов и кубов для областей с несколькими материалами, а затем опишем предлагаемую моди­ фикацию этих алгоритмов.

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

Квадратная ячейка имеет четыре вершины, каждая из которых имеет определенный индекс подобласти. Таким образом, всего имеется 44 = вариантов распределения индексов подобластей и построения кусочно-линей­ ного контура на этой грани. Однако, они могут быть сгруппированы в следу­ ющие случаи:

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

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

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

Обозначим через - квадратную грань, имеющую ребра и вершины. Метод M3 S может быть описан следующим алгоритмом 1.3:

Алгоритм 1.3 Марширующих квадратов для областей с несколькими мате­ риалами (M3 S) 1: Вход: Грань если имеет различные индексы на концах тогда Присвоить число точек в 10:

Построить интерфейсную точку p на 11:

12:

13:

Определить количество различных индексов в вершинах 14:

15:

Разрешить неоднозначность и выбрать пары точек для соединения 16:

17:

Соединить 2 пары точек так, чтобы вершины с одинаковым индек­ 18:

сом остались связанными 19:

20:

21:

22:

23:

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

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

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

В зависимости от количества интерфейсных точек на гранях, можно опи­ сать алгоритм замыкания контуров:

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

Если куб имеет ровно две интерфейсных точки на гранях (ровно одну он иметь не может, см. [92]), то некоторые компоненты формируют неза­ мкнутые ломаные с концами в этих интерфейсных точках. Соединив их, сделаем все незамкнутые компоненты замкнутыми. Для триангуляции таких компонент также используется метод “разделяй и властвуй”.

Если имеется более двух интерфейсных точек на гранях, то создается точка в центре ячейки (аналогично, будем называть ее интерфейсной точкой ячейки). Это точка соединяется отрезками со всеми интерфейс­ ными точками на гранях, после чего все компоненты становятся связ­ ными. Каждая компонента, содержащая интерфейсную точку ячейки, триангулируется с помощью классического метода “веер треугольников” с центром в этой точке. Если компонента не содержит интерфейсную точку ячейки, то она триангулируется методом “разделяй и властвуй”.

Метод 3 для ячейки может быть описан следующим алгоритмом 1.4:

Алгоритм 1.4 Марширующих кубов для областей с несколькими материа­ лами 1: Вход: ячейка если создана интерфейсная точка p на грани тогда Присвоить число точек в Соединить две интерфейсные точки на гранях 10:

11:

Построить интерфейсную точку p в центре 12:

13:

14:

15:

16:

Выделить замкнутые компоненты контуров 17:

18:

19:

Построить триангуляцию “веер треугольников” для 20:

21:

Построить триангуляцию “разделяй и властвуй” для 22:

23:

24:

25:

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

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

1.2.3. Модификация метода марширующих кубов для областей с несколькими материалами Будем считать, что у нас имеется область, заданная функцией, для каж­ дой точки возвращающей индекс подобласти, в которой она находится. Кро­ ме этого будем считать, что построена сетка типа восьмеричное дерево, ко­ торую мы можем модифицировать. Итак, вместо точек в серединах ребер, интерфейсные точки на ребрах вычисляются с помощью метода бисекции по заданной функции. Как и в методе CMS, на каждом ребре должно быть не более одной интерфейсной точки. Для этого необходимо измельчить все ячей­ ки, содержащие ребра с более чем одной интерфейсной точкой. В противном случае часть области будет отслежена сеткой недостаточно точно. В каче­ стве интерфейсной точки на грани вместо центральной точки используется более точное положение. Для ее вычисления используется идея, аналогичная поиску особенной точки области внутри ячейки. Обозначим интерфейсные точки на ребрах за s, проекции нормалей в этих точках на плоскость грани за n. Искомую точку p будем находить как точку пересечения касательных элементов точек s. Другими словами, p удовлетворяет системе (1.4).

В результате решения системы может получиться точка, не принадле­ жащая соответствующей грани. Это может быть вызвано неточностью вход­ ных данных или наличием особенностей на поверхностях интерфейсов.Кроме этого, это может произойти в случае, если одно ребро ячейки пересекает бо­ лее чем один интерфейс. Как было отмечено ранее, в таком случае следует иерархически разбить ячейки, содержащие такое ребро. Если интерфейсную точку на грани сразу найти не удалось, то можно сделать несколько уровней иерархического измельчения грани, и попробовать найти ее в одной из по­ лученных дочерних граней. Стоить заметить, что этот процесс измельчения можно применить только для поиска интерфейсной точки на грани, но при этом итоговую сетку оставить без изменений, если измельчение сетки неже­ лательно. Для простоты можно воспользоваться альтернативным способом выбора интерфейсной точки на грани, используя, например, центр масс то­ чек s. Однако, при неточном выборе интерфейсной точки на грани, на сетке могут присутствовать визуальные артефакты, например, появление “зубчи­ ков”, см. Рис. 1.7.

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

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

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

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

Когда некоторые интерфейсные точки на ребрах попадают в вершины грани, интерфейсные точки на грани также могут совпасть с вершиной грани, либо попасть на ее ребро. Действительно, для примера рассмотрим грань, в которой две соседние вершины имеют индекс подобласти 1, две другие индексы 2 и 3 соответственно, см. Рис. 1.9. Если интерфейсная точка на ребре совпадет с, а интерфейсная точка на ребре совпадет с, то, очевидно, интерфейсная точка на грани будет лежать на ребре.

Причем, если интерфейсная точка на ребре совпадет, скажем, с вершиной Рис. 1.9. Попадание интерфейсной точки на грани на ребро 1.9, а и в вершину 1.9, б.

, то тогда интерфейсная точка на грани также совпадет с вершиной, см.

Рис. 1.9 (б ).

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

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

Далее, в зависимости от расположения интерфейсной точки на второй грани, рассмотрим несколько вариантов:

Если у нет интерфейсной точки, то на грани не будут построены дополнительные отрезки.

Пусть p - интерфейсная точка на грани и p. Так как на одном ребре может быть не более одной интерфейсной точки, то p = p, а, значит, общее количество интерфейсных точек на гранях ячейки умень­ шится на один. Покажем, что это не приведет к изменению алгоритма замыкания контуров. Отличие может возникнуть, если новое количе­ ство интерфейсных точек на гранях стало равно одному. Однако, это эквивалентно ситуации, когда их число равно нулю, то есть все компо­ ненты контуров являются замкнутыми, см., например, Рис. 1.10(а).

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

Если концы ребра имеют одинаковый индекс, то на нем нет интер­ фейсной точки. Если нужно соединить точки p и p, то этот отрезок будет новым для грани, см. Рис. 1.10(в). Этот отрезок нужно учиты­ вать со стороны соседней по грани ячейки. Как было сказано ранее, это не приводит к изменениям алгоритма для этой ячейки, однако нуж­ но учесть наличие новой точки на ребре с помощью постпроцессинга, описанного ранее.

Рис. 1.10. Попадание на ребро интерфейсной точки на грани Случай, когда интерфейсная точка на грани попадает в вершину, ана­ логичен случаю попадания на ребро, при этом новые точки на ребрах не появляются.

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

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

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

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

Таким образом, метод CMS является частным случаем этой модифика­ ции метода M3 C. Соответственно, полученная модификация имеет преиму­ щества обоих методов:

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

Модифицированный метод будем называть методом кубических марши­ рующих квадратов для областей с несколькими материалами (Multiple material cubical marching squares, MMCMS). Опишем общий алгоритм построения сет­ ки типа восьмеричное дерево со сколотыми ячейками для составной области Будем считать, что входными данными является функция, которая для любой точки возвращает индекс подобласти, которой она принадлежит, а также нормаль к внутренним или внешним границам, если точка находится вблизи границы. Такую функцию можно восстановить из различных типов входных данных. Например, если область задана поверхностной триангуляци­ ей, то функцию можно получить с помощью алгоритма отслеживания лучей (ray-tracing). Стоит отметить, что вершины этой триангуляции, вообще гово­ ря, не переносятся на итоговую сетку. Если точка находится вне области, ей присваивается фиктивный индекс, обозначающий внешнюю подобласть.

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

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

Построение сетки типа восьмеричное дерево со сколотыми ячейками для области может быть представлено в виде следующего алгоритма 1.5:

Алгоритм 1.5 Построения сетки типа восьмеричное дерево со сколотыми ячейками 1: Построить сетку типа восьмеричное дерево для Присвоить вершинам сетки индексы подобластей если Число различных индексов в больше 1 тогда Разделить ячейку на части с помощью триангуляции Удалить ячейки с фиктивным индексом 10:

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

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

Напомним, что сетку типа восьмеричное дерево можно рассматривать как конформную многогранную сетку. Сетка со сколотыми ячейками также принадлежит классу конформных многогранных сеток, если область не раз­ Рис. 1.11. Кубическая ячейка, поделенная на две части триангуляцией деляется на подобласти. В случае, если область разделена на подобласти, некоторые соседние ячейки могут иметь более одной общей грани. Действи­ тельно, когда кубическая ячейка разбивается на части методом MMCMS, граница между этими частями представляет собой набор из нескольких тре­ угольников, см. Рис. 1.11. Поэтому сетку со сколотыми ячейками для состав­ ных областей можно отнести к слабо-конформным многогранным сеткам.

1.4. Анализ алгоритма Проведем анализ Алгоритма 1.5, покажем конечность его работы и оце­ ним вычислительную сложность.

1.4.1. Конечность работы алгоритма Процесс построения сетки типа восьмеричное дерево описан в различной литературе, см., например, [69, 84]. Он конечен при ограниченной глубине из­ мельчения дерева. Процесс дополнительного измельчения сетки также огра­ ничен заданной глубиной измельчения дерева.

Рассмотрим цикл по ячейкам в строке 3. Для каждой ячейки в алгоритме MMCMS в цикле по граням применяется модифицированный алгоритм M3 S.

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

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

Далее, для каждого замкнутого контура строится триангуляция. Очевидно, что метод триангуляции “веер треугольников” конечен, так как требует лишь соединения центра триангуляции с остальными не соседними вершинами. Ме­ тод триангуляции “разделяй и властвуй” требует поиска подходящей разде­ ляющей плоскости. Как отмечено в [83], для произвольных последователь­ ностей вершин существуют примеры, в которых разделяющая плоскость не будет найдена данным алгоритмом. Однако, в [50] автор указывает, что такие случаи могут возникать лишь в более сложных последовательностях вершин, которые не возникают в методе M3 C. В любом случае докажем, что если не удалось построить триангуляцию методом “разделяй и властвуй”, то ее можно построить альтернативным способом.

Утверждение 1.4.1. Пусть имеется ячейка сетки типа восьмеричное де­ рево, причем сетка удовлетворяет условиям, накладываемым алгоритмом MMCMS. Тогда процесс построения триангуляции для замкнутого конту­ ра, полученного алгоритмом, конечен.

Доказательство.

Сначала рассмотрим случай, когда внутри ячейки не создается интер­ фейсная точка. Предположим, что не удалось построить триангуляцию мето­ дом “разделяй и властвуй”. Тогда она строится методом “веер треугольников”, где в качестве центра триангуляции берется особенная точка ячейки, либо центр масс интерфейсных точек. Если же ячейка имеет внутри интерфейс­ ную точку, то все контуры, проходящие через нее, триангулируются этим же способом. Осталось рассмотреть случай, когда ячейка имеет интерфейс­ ную точку, но контур не проходит через нее. Докажем, что если контур не проходит через интерфейсную точку ячейки, то он не проходит также ни че­ рез какую интерфейсную точку на грани. Предположим, что интерфейсная точка грани принадлежит последовательности и рассмотрим две ее соседние точки. Ни одна из них не является интерфейсной точкой ячейки, так как контур не проходит через нее. В то же время, ни одна из них не является другой интерфейсной точкой грани, так как такие точки не могут быть со­ единены в данном случае. Следовательно, оба ее соседа - это интерфейсные точки на ребрах, причем на этой же грани. Этим точкам приписаны ровно по два индекса подобластей, которые они разделяют, причем эти пары ин­ дексов не совпадают. Но тогда эти две точки не могут принадлежать одному контуру, так как контуру соответствует только одна пара индексов. Таким образом, контур не проходит через интерфейсные точки на гранях, и, сле­ довательно, не проходит через грани, на которых есть такая интерфейсная точка. Так как имеется интерфейсная точка внутри ячейки, то согласно алго­ ритму замыкания контуров, число интерфейсных точек на гранях не меньше трех. Следовательно, число граней, на которых нет интерфейсных точек, не больше трех.

Рассмотрим случай, когда таких граней ровно три, и контур проходит через них. Очевидно, что они имеют одну общую вершину, иначе контур был бы незамкнутый. Для наглядности обозначим последовательно вершины за 1, 2,...,, 1, 2,.., 1,...,, причем 1, 1, 1 - интерфейсные точки на ребрах, см. Рис. 1.12. Рассмотрим первую вершину, не совпадающую с Рис. 1.12. Замкнутый контур (а), его триангуляция (б ).

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

Если таких граней две, то рассмотрим их общее ребро. Так как контур замкнутый, то он пересекает это ребро в двух точках (как минимум одна из них совпадает с вершиной ребра), которые разделяют контур на две части.

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

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

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

1.4.2. Вычислительная сложность алгоритма Проведем краткий анализ вычислительной сложности работы Алгорит­ ма 1.5. Будем считать, что область задается функцией : R3 Z, возвраща­ ющую индекс подобласти, которой принадлежит точка (,, ). В зависимо­ сти от типа входных данных, вызов этой функции может быть относительно дорогой, поэтому будем учитывать ее сложность.

Максимальная сложность построения сетки типа восьмеричное дерево порядка ( ln ), где - число вершин дерева.

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

1. Поиск интерфейсной точки на ребре и нормали к ней;

2. Поиск интерфейсной точки на грани и в ячейке;

3. Построение триангуляции замкнутого контура.

Рассмотрим ребро ячейки, концы которого имеют разные индексы. Для нахождения интерфейсной точки на ребре необходимо сделать шагов ме­ тода бисекции, где зависит от желаемой точности сетки. На каждом шаге один раз вызывается функция, значит общая сложность равна ·. В гексаэдральной сетке типа восьмеричное дерево имеется следующая оценка на количество ребер: 3. Предположим, что на каждом ребре может находиться интерфейсная точка, тогда максимальная сложность поиска всех интерфейсных точек пропорциональна · · 3 = ( · ). Для вычисления нормали, в зависимости от типа входных данных, используется один из способов, описанных в [60]. Так, например, если область задается поверхностной триангуляцией, то используется нормаль к ближайшему тре­ угольнику.

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

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

Следовательно, число не превосходит 48, хотя существует более точная оценка.

Заметим, что применение алгоритма MMCMS является локальным, и его сложность для одной ячейки не зависит от общего числа ячеек сетки. При построении сколотой ячейки используются только соседние ячейки. Таким образом, общая сложность алгоритма пропорциональна 1 + 2 ln, причем 1, 2 линейно зависят от. Причем сложность построения сколотых ячеек является линейной.

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

1.5.1. Улучшение качества полученной сетки После построения сетки, некоторые сколотые ячейки могут иметь очень маленький объем. Это отрицательно сказывается при дискретизациях на та­ ких сетках. Поэтому предлагается алгоритм объединения маленьких ячеек с соседними по граням. Кроме критерия маленького объема может быть ка­ кой-либо другой, например, толщина ячейки. Рассмотрим ячейку, которую мы хотим объединить с соседней. Среди соседних ячеек выбираются те, у которых такой же индекс подобласти с данной ячейкой. Далее выбирается ячейка, с которой их общая грань имеет наибольшую площадь. Эта ячейка объединяется с рассматриваемой простым удалением их общей грани. При этом ребра удаленной грани остаются в ячейке. Стоит отметить, что, вообще говоря, у объединяемых ячеек могут быть разные родительские ячейки. Это следует учитывать при динамическом перестроении сетки.

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

которых единичные нормали n1, n2 удовлетворяют условию (n1, n2 ) 0.99.

Количество граней на границе равно 1208 и 708, соответственно.

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

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

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

Как было отмечено ранее, при динамическом перестроении сетки необхо­ димо следить за объединенными ячейками, см. раздел 1.5.1. Если объединен­ ные ячейки имели разных родителей, то когда один из них перестраивается, то также должен быть перестроен второй, после чего можно снова объеди­ нять полученные новые ячейки.

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

Рассмотрим область, состоящую из почти горизонтальных слоев. Вве­ дем оператор проектирования P : R3 R2 на плоскость = 0. Обозначим двумерную область = P() и ее границу. Будем считать, что область представлена набором интерфейсов и слоев. Каждый слой находится между интерфейсами и +1. Кроме этого, в области могут быть разломы, представляющие собой наклонные поверхности. Будем счи­ тать, что каждый интерфейс задан некоторой функцией, которая для любой точки (, ) возвращает -координату соответствующей точки на интерфейсе:

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

санным ранее способом.

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

Обозначим их количество за. Остальные интерфейсы маркируем как пе­ ресекающиеся, они будут использованы позже. Затем строится сетка типа восьмеричное дерево со сколотыми ячейками в области, состоящей из ровных горизонтальных слоев, разделенных интерфейсами. Вообще го­ воря, сетка может сгущаться к различным местам области. Однако, для силь­ но искривленных слоев алгоритм будет надежным только в случае, если измельчение сетки типа восьмеричное дерево в вертикальном направлении происходит полностью. Другими словами, сетка должна быть представима в виде тензорного произведения двумерной сетки типа четвертичное дере­ во в плоскости и одномерной сетки вдоль оси. Такое ограничение обусловлено тем, что, в противном случае, в результате кусочно-линейного отображения этой сетки вдоль оси, некоторые ячейки могут оказаться вывернутыми, см. Рис. 1.14. Область является составной, где каждый из сло­ ев является подобластью. Каждый интерфейс принадлежит плос­ кости = для некоторого числа. Боковой границей области является цилиндрическая поверхность [1, ]. Для аппроксимации боковой гра­ ницы используются сколотые ячейки. Полученную сетку назовем прообразом сетки.

После этого мы можем построить семейство кусочно-линейных отобра­ жений для -координат вершин прообраза сетки, переводящее ровные слои в искривленные. Определим это отображение следующим образом. Рас­ смотрим точку из прообраза (,, ), тогда [, + 1), где определяется положением интерфейса. Тогда линейно по на отрезке образом, точка (,, ) получит новую координату = (,, ). В резуль­ тате такого отображения, все интерфейсы прообраза сетки попадут на интерфейсы.

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

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

Описанный выше метод может быть сформулирован в виде алгоритма 1.6 с использованием введенных ранее обозначений.

Алгоритм 1.6 Построения сетки в слоистой области 1: Вход: Интерфейсы, разломы, боковая граница если не пересекается с другими интерфейсами тогда Создать сетку типа восьмеричное дерево из 1 горизонтальных слоев 10:

Сколоть ячейки, пересекающие [1, ] 11:

12:

13:

14:

15:

16:

17:

18:

19:

20:

21:

22:

23:

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

Рис 1.15.

1.7. Примеры сеток На основе Алгоритма 1.5 был программно реализован сеточный генера­ тор, написанный на языке C++. В этом разделе представлены некоторые при­ меры сеток, построенных с помощью этого генератора. В данных примерах не использовался алгоритм объединения граней, лежащих в одной плоско­ сти, поэтому сетка на поверхности всегда является треугольной. Кроме этого была создана версия генератора на платформе INMOST [1]. Преимущества этой платформы заключаются в том, что она предоставляет возможности для параллельного построения сеток, а также дискретизаций и решения алгебра­ ических систем. Сетки, построенные этой версией генератора, будут также продемонстрированы.

1.7.1. Однородная область В качестве в однородной области рассмотрим модель Стэнфордского кро­ лика [16], см. Рис. 1.16. Входными данными является поверхностная триангу­ ляция области, с помощью которой была восстановлена характеристическая функция области. В полученной сетке 49997 ячеек, из них 14742 сколотых, 68244 вершины и 191059 граней. Сетка сгущается к границе области.

1.7.2. Составная область Больший интерес представляют составные области, ведь построение се­ ток в таких областях является ключевым результатом работы.

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

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

Таблица 1.1. Параметры сетки для области из пересечения примитивов и время построения Уровень Время, сек # ячеек # cколотых ячеек # граней # вершин разделе 1.4.2: = 10, = 1. В таблице 1.1 представлены параметры по­ лученных сеток, а также время их построения. В данном примере основное время работы алгоритма занимает построение сколотых ячеек. Можно заме­ тить квазилинейную зависимость между числом сколотых ячеек и временем работы алгоритма. На Рис. 1.17 изображена сетка на первом уровне измель­ чения, а также некоторые ее срезы. Разные подобласти изображены разным цветом.

Далее, рассмотрим область, состоящую из шести геологических слоев.

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

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

Таблица 1.2. Параметры сетки для области из 6 слоев.

Уровень Время, сек # ячеек # cколотых ячеек # граней # вершин Рис. 1.18. Сетка для 6 слоев (а), ее срез (б ). Увеличения сетки в окрестности интерфейсной точки на грани (в),(г).

1.7.3. Адаптирующиеся сетки Для модели из шести слоев можно построить отображенную сетку, опи­ санную в разделе 1.6. Построим набор динамически изменяющихся сеток та­ кого вида.

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

1.7.4. Приложения сеток Генератор сеток типа восьмеричное дерево со сколотыми ячейками для слоистых областей используется в расчетном комплексе GeRa (Geomigration of Radionuclides - совместный проект ИВМ и ИБРАЭ РАН в рамках проекта “Прорыв” ГК Росатом). В этом комплексе он является частью технологиче­ ской цепочки расчетов. Кроме этих сеток в данном комплексе используются треугольные призматические сетки. На Рис. 1.20 изображено окно графиче­ ского интерфейса комплекса GeRa. В нем можно видеть контур границы об­ ласти, а также отмеченную на ней скважину. Справа изображена поверхность сетки типа восьмеричное дерево со сколотыми ячейками. Сетка сгущается к границе области, а также к скважине.

Рис. 1.19. Срезы динамически перестраивающейся сетки Рис. 1.20. Окно графического интерфейса комплекса GeRa.