CC BY
25
УДК 519.6242, 681.3
Алгоритм построения двумерных вложенных сеток
Антон В.Ищенко*
Институт математики, Сибирский федеральный университет, Свободный 79, Красноярск, 660041
Россия
Игорь В.Киреев^
Институт вычислительного моделирования СО РАН, Академгородок 50, Красноярск, 660036,
Россия
Получена 10.11.2008, окончательный вариант 15.12.2008, принята к печати 10.12.2008 В работе предлагается алгоритм построения последовательности неравномерных вложенных сеток для двумерных многосвязных областей, содержащий в себе достоинства метода шаблонов и алгоритма бисекции.
Ключевые слова: сетка, триангуляция, конечные элементы.
В последние годы особое значение среди проекционно-сеточных методов приобретают многосеточные алгоритмы конечных элементов [1] как наиболее экономичные. Но для их реализации необходимо построить не только триангуляцию расчётной области, но и последовательность вложенных сеток для неё. В данной работе предлагается алгоритм разбиения двумерных многосвязных областей на треугольники, результатом работы которого является последовательность неравномерных вложенных друг в друга сеток.
Под сеткой [2] будем понимать некоторое топологическое множество точек (вершин, узлов), связанных между собой ребрами - отрезками прямых (а в некоторых случаях и кривых) линий таким образом, что исходная область разбивается на элементы определенной формы. При этом в качестве элементов сетки, если речь идет о геометрически сложных плоских областях, обычно используются треугольники, а сам процесс построения сетки из треугольников обычно называется дискретизацией или триангуляцией.
Все методы триангуляции по принципу построения можно разбить на две большие группы: итерационные [3] и прямые [4]. Предлагаемый в работе алгоритм триангуляции содержит в себе основные достоинства итерационного метода граничной коррекции, прямых методов шаблонов [5] и бисекции [6].
Методы граничной коррекции являются самыми быстрыми из итерационных методов. Построение сеток в этих методах осуществляется в два этапа. На первом этапе производится триангуляция некой простой суперобласти, полностью включающей в себя заданную область. Как правило, эта суперобласть представляет собой прямоугольник, триангуляция которого осуществляется методом на основе одного из шаблонов, что подразумевает разбиение исходной области на подобласти заданного вида. Для каждого вида области используется свой шаблон, то есть принцип размещения узлов и установки связей между
* e-mail: Al357@yandex.ru t e-mail: kiv@icm.krasn.ru © Siberian Federal University. All rights reserved
ними. Каждый шаблон применим только к областям заданного вида. Благодаря такой узкой специализации сетки, построенные на шаблонах, часто бывают высокого качества. На втором этапе все узлы полученной сетки, лежащие вблизи границы заданной области, проецируются на её границу, а узлы, лежащие вне заданной области, — удаляются. Чтобы компенсировать неизбежные геометрические искажения элементов сетки вблизи границ, часто дополнительно проводят еще один этап — этап оптимизации сетки, что в итоге позволяет получить достаточно хорошие результаты.
В основе предлагаемого алгоритма триангуляции лежит единственность представления прямоугольного равнобедренного треугольника в виде объединения двух треугольников, которые получены из исходного разбиением на два равных прямоугольных равнобедренных треугольника медианой, как изображено на рис. 1. Эту процедуру разбиения коротко называют бисекцией.
Как и в методе граничной коррекции, исходную многосвязную двумерную область О помещаем в прямоугольник Д, который изначально разбит на равнобедренные прямоугольные треугольники. Т. е. в методе шаблонов в качестве базового элемента выбираем равнобедренный прямоугольный треугольник, шаблоном для которого является представление в виде объединения двух равных равнобедренных прямоугольных треугольников, полученных бисекцией первоначального, как на рис. 1. После триангулирования объемлющего прямоугольника Д при помощи граничной коррекции получаем сетку для области О.
Триангуляция области О состоит из трёх структур - узлов, сторон (катеты, гипотенузы) и треугольников.
Всякий узел Ык, к = 1,..., Жр сетки характеризуются двумя декартовыми координатами хь и ук и своим порядковым номером. Естественно, что в рассматриваемом случае всякий узел треугольной сетки должен быть только вершиной одного или нескольких треугольников, число которых не более 8.
Всякое ребро сетки (сторона треугольника) определяется узлами, его ограничивающими. Рёбра, соответствующие гипотенузе, выделяются особо, поскольку бисекция в предлагаемом алгоритме делит пополам в первую очередь именно гипотенузу.
Каждый треугольник сетки имеет свой номер и полностью описывается тремя вершинами, одна из которых, соответствующая вершине прямого угла, выделяется особо. Через Жд будем обозначать число всех треугольников сетки объемлющего прямоугольника Д, а через ¡к - длину гипотенузы треугольника Дк, которую можно интерпретировать как его линейную меру; к = 1,... ,Жд. Все треугольники до проведения коррекции являются прямоугольными и равнобедренными.
Одним из основных понятий в многосеточных методах является понятие уровня вложенности треугольника. Все треугольники начального разбиения прямоугольника Д имеют
нулевой уровень вложенности; треугольники, полученные путем разбиения треугольников нулевого уровня на два, имеют первый уровень вложенности и т.д. Линейная мера треугольника при его бисекции уменьшается в раз.
Введём на прямоугольнике К функцию шагов к(М). Для любой точки М € К значение к(М) всегда больше нуля и определяет линейную меру треугольника сетки, которому принадлежит точка М, следующим образом: бисекция треугольников сетки продолжается до тех пор, пока каждый элемент триангуляции Дк не будет удовлетворять соотношению
¡к < тах к(М). (1)
м едк
На практике, вместо этого критерия обычно выполняют проверку неравенства ¡к ^ к(М) лишь в нескольких точках М треугольника Дк, и если хотя бы в одной из них последнее неравенство выполнено, то Дк больше не делят (к = 1, 2,..., Жд). Если известны ктах и кт\п
прямоугольника
— максимальное и минимальное значения к(М) на всей П соответственно, то разница между уровнями вложенности треугольников по порядку совпадает с числом 2log2(kmax/kmin).
Проиллюстрируем на конкретном примере процесс дробления сетки. Пусть в объемлющем прямоугольнике задана триангуляция, изображённая на рис. 2, а, и для одного из её треугольников условие (1) не выполнено. На рис.3 изображена последовательность преоб-
разований, порождённых бисекцией такого треугольника. Треугольники, подлежащие разбиению на следующем шаге, выделены штриховкой.
Сначала в разбиваемом треугольнике проводится высота. Получаем точку на гипотенузе, которую пока ещё нельзя интерпретировать как узел триангуляции, поскольку узлы должны находиться только в вершинах треугольника. Если рассматриваемая гипотенуза является гипотенузой и смежного треугольника, то в последнем проводится высота и множество всех узлов пополняется только что построенной точкой. В этом случае из списка
треугольников сетки удаляется два модифицированных треугольника и вносятся четыре новых, так что общее число треугольников сетки увеличивается на два. При этом старая сетка вкладывается в новую и узлы исходной триангуляция служат узлами модифицированной. На этом очередной этап дробления сетки заканчивается, и переходим к анализу другого треугольника.
Если же построенная середина гипотенузы есть точка катета смежного треугольника, а не его гипотенузы (такое имеет место только в случае, когда смежный треугольник имеет более низкий уровень вложенности), то этот смежный прямоугольный треугольник разбивается сначала высотой на два, а затем изложенная, процедура повторяется со вновь полученным треугольником того же уровня вложенности, что и исходный. Но теперь возникает точка — середина гипотенузы смежного треугольника и процесс повторяется с момента распознавания "катет/гипотенуза". Эти преобразования трансформируют сетку на рис. 2, а в триангуляцию, изображённую на рис. 2, б.
Описанный процесс является рекурсивным и конечным, поскольку мы либо выходим на границу прямоугольника Д, либо сталкиваемся с гипотенузой двух смежных треугольников одного и того же уровня вложенности (рис. 3). По существу доказано следующее утвержде-
Теорема 1. Процесс построения вложенных сеток объемлющего прямоугольника Д, основанный на бисекции какого-либо треугольника сетки, является конечным, а возникающие треугольники наибольшего уровня вложенности будут сгруппированы по 4, образуя квадрат одного из двух типов, изображённых на рис. 4-
Рис. 4. Структуры, образуемые треугольниками наибольшего уровня вложенности
Структурированность построенной триангуляции позволяет по заданной функции шагов к(Ы) получить априорную верхнюю оценку числа возникающих в процессе разбиения треугольников. Имеет место
Теорема 2. Если стратегия триангуляции объемлющего прямоугольника Д такова, что для каждого треугольника Д^ вместо неравенства (1) выполняется равенство
= тах Н(Ы),
м еАк
(2)
то число возникающих треугольников Жд удовлетворяет неравенству
к
Доказательство. Пусть по заданной функции шагов Л(х, у) построено разбиение прямоугольника Я на Жд равнобедренных прямоугольных треугольников Дк, так что Я =
k=N△
и Дк. Обозначим через Як = /4 площадь треугольника Дк из этой триангуляции;
к=1
здесь /к — длина его гипотенузы, удовлетворяющая равенству (2). Но тогда, в силу известных теорем классического математического анализа,
Г Г 4 ¿хс!у ^^ [! 4 ¿хйу Г Г 4 ¿хд/у N△0 4 N△'
. . -ЩГУ)=^1. ЩХ~У) ^^ ИГ = 2^ЬкЦ = 2^ =
У к = 1 дк к= 1 дк к к=1 к к = 1
т.е. выполнено неравенство (3). □
Верхнюю оценку (3) можно использовать для прогнозирования вычислительных затрат. Для функции шагов вида
Л(х,у) = Но + (Л1 - Ло) (1 + (^)Р + (^)')" ,
зависящей от 9 параметров [7], было проведено исследование точности неравенства (3). Результаты численных экспериментов при значениях параметров
Л0 е [1,1.5], Л1 =0.05, х0 = у0 = 2.5, а = Ь = 4, р = д = 2, г = -1
представлены на рис. 5. Вычисление интеграла в правой части неравенства (3) проводилось по квадратурным формулам Гаусса 32-го порядка точности [8].
Жд
0.9 1.0 1.1 1.2 1.3 1.4 1.5 Рис. 5. Реальное {-} и прогнозируемое {— • —} количество Жд треугольников сетки
Обозначим через дП границу области П, содержащейся в объемлющем прямоугольнике Я. После построения по заданной функции шагов Л(М) триангуляции Я, все треугольники по отношению к П разбиваются на три типа: внутренние, внешние и граничные.
В каждом граничном треугольнике определяем стороны, пересекающие дП. Вершины Вк таких сторон будем называть приграничными и вычисляем координаты точек пересечения Вк = N к ^ 8) этих сторон с дП. Для каждого приграничного узла В к вычисляем расстояние от него до д П как минимальное из длин отрезков В к Б]]. Возникает упорядочение приграничных узлов. Начиная с наименее удалённого от границы приграничного узла, последовательно перемещаем приграничные узлы Вк в ближайшую к Вк точку из множества {Б^,..., } С дП.
Одним из отрицательных последствий граничной коррекции узлов сетки является искажение некоторых треугольников сетки, подчас весьма значительное. Поэтому целесообразно оптимизировать узлы сетки, перемещая каждый внутренний узел из П в центр тяжести соответствующего ему —¡-угольника: для каждого внутреннего узла сетки вычисляем его новое положение М(п)(х((П), у(п)) следующим образом [7]: тк тк
(п) 1 С"-1) (") 1 С"-1) С0) С0) !л\
х к = — 1^х к ^ Ук = ^ х к =х k, у к = Ук. (4)
т к - т к -]=1 ]=1
Здесь 4 < тк < 8 - количество узлов сетки, связанных с данным узлом Мк, а п - номер итерации.
Заметим, что все элементы матрицы перехода итерационного процесса (4) неотрицательны, а поскольку область П является связной, то матрица перехода этого процесса будет и неразложимой. Кроме этого, в этом процессе (4) изменяются координаты только внутренних узлов сетки, стало быть, сумма элементов строки матрицы перехода, содержащей ссылки на граничные узлы, будет меньше 1. Поэтому из спектральных свойств неразложимой матрицы с неотрицательными элементами [9, §18.9] следует
Теорема 3. Итерационный процесс (4) оптимизации внутренних узлов сетки сходится.
Рис. 6-7 иллюстрируют процесс граничной коррекции и оптимизации.
Рис. 6. Триангулгяция прямоугольника К, содержащего круг с отверстием: а — результат работы алгоритма на прямоугольнике К; б — результат граничной коррекции сетки прямоугольника К
Обратный процесс — процесс перехода к предыдущему уровню вложенности (вплоть до нулевого) осуществляется при помощи укрупнения треугольников самого высокого на данном этапе уровня вложенности проходит в несколько этапов и основан на теореме 1.
Рис. 7. Триангуляция Округа с отверстием: а — сетка после граничной коррекции; б — сетка после оптимизации
Рис. 8. Процесс объединения: а — исходное состояние, б — шаг 1 и в — шаг 2
В каждом квадрате, образуемом четверкой треугольников самого высокого уровня вложенности, выбрасывается та его диагональ, которая построена последней. Если перед началом процесса объединения все треугольники были прямоугольными, то линейные размеры треугольников на каждом этапе увеличиваются в а/2 раз. Работа алгоритма объединения продемонстрирована на рис. 8, 9.
Основные результаты этой работы были доложены на VII Всероссийской конференции молодых учёных по математическому моделированию и информационным технологиям [10].
Работа выполнена при финансовой поддержке РФФИ: грант 08-01-00621-а "Вычислительные технологии для расчета течений несжимаемой жидкости.
Список литературы
[1] В.В.Шайдуров, Многосеточные методы конечных элементов, М., Наука, 1989.
[2] P.J.Frey, P.L.Geotge, Mesh Generation. Application to Finite Elements, Oxford, Paris, Hermes Science Publishing, 2000.
[3] М.П.Галанин, И.А.Щеглов, Разработка и реализация алгоритмов трехмерной триангуляции сложных пространственных областей: итерационные методы, М., ИПМ им. М.В.Келдыша РАН, 2006, Препринт №9, 32с.
Рис. 9. Процесс объединения: а — наивысший уровень вложенности треугольников; б —
сетка после одного шага алгоритма; в —сетка после двух шагов алгоритма
[4] М.П.Галанин, И.А.Щеглов, Разработка и реализация алгоритмов трехмерной триангуляции сложных пространственных областей: прямые методы, М., ИПМ им. М.В.Келдыша РАН, Препринт №10, 2006, 32с.
[5] S.J.Owen, A Survey of Unstructured Mesh Generation Technology, Proceedings of 7th International Meshing Roundtable, Dearborn, MI., 1998, 239-269.
[6] M.C.Rivara, M.M.Vemere, Cost analysis of the longest-side (triangle bisection) refinement algorithm for triangulations, Engineering with Computers, 12(1996), №3-4, 224-234.
[7] С.Ф.Пятаев, Ю.В. Немировский, Автоматизированная триангуляция многосвязных областей со сгущением и разрежением узлов, Вычислительные технологии, 5(2000), №2, 82-91.
[8] В.И.Крылов, Л.В.Шульгина, Справочная книга по численному интегрированию, М., Наука, 1966.
[9] В.В.Воеводин, Ю.А.Кузнецов, Матрицы и вычисления, М., Наука, 1984.
[10] А.В.Ищенко, Алгоритм построения треугольных вложенных сеток для двумерных многосвязных областей, Материалы VII Всерос. конф. молодых учёных по мат. моделированию и информационным технологиям, Красноярск, ИВТ СО РАН, 2006, 21.
The Algorithm for Generation of Two-Dimensional Embedded Grids
Anton V.Ishenko Igor' V.Kireev
The paper proposes an algorithm for constructing a sequence of embedded grids for two-dimensional areas.
Keywords: embedded grid, two-dimensional areas.
