Об измельчении триангуляции вблизи границы области Текст научной статьи по специальности «Математика»

ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА

Сер. 10. 2011. Вып. 1

ИНФОРМАТИКА

УДК 519.6

Е. П. Арсентьева

ОБ ИЗМЕЛЬЧЕНИИ ТРИАНГУЛЯЦИИ ВБЛИЗИ ГРАНИЦЫ ОБЛАСТИ

1. Введение. Невырожденная триангуляция области важна для интерполяции и аппроксимации функций, заданных в этой области, а также при решении задач математической физики методом конечных элементов [1-6]. Курантовское приближение для функций и G С'1,“(Г2), a G (0,1], заданных на плоской многоугольной области И, имеет порядок h1+a, где h - максимальная длина ребра триангуляции [1, 2]. Если же функция и или ее производные неограниченно растут при приближении к границе, то для достижения качественной аппроксимации следует измельчить треугольники вблизи границы. Известны методы измельчения сеток (метод Рапперта, см. [3]) и критерий выбора зоны измельчения (критерий Доффлера, см. [4]). В работе [5] рассмотрены алгоритмы триангуляции двумерной области со сгущением вложенных сеток. В [6] построены структурированные сетки в двумерных областях с помощью методов отображений, учитывающих форму ячеек.

Как известно, в случаях, когда углы в треугольниках сетки близки к нулю, погрешность округления сильно увеличивается из-за роста числа обусловленности матрицы, используемой в методе конечных элементов. Цель данной работы - предложить новые варианты измельчения треугольных сеток вблизи границы области с сохранением свойства невырожденности. Здесь предлагаются потенциально бесконечные невырожденные триангуляции, аппроксимирующие границу области с любой наперед заданной точностью. При этом используется метод параметрического задания пограничной полосы и метод дробления четных треугольников. Применяя известные конечно-элементные аппроксимации (например, функции Р. Куранта [1, 2]), с помощью предлагаемых триангуляций можно получить эффективные приближения для функций с сингулярностями на границе области.

2. Описание алгоритма измельчения триангуляции вблизи границы области. Введем предварительные обозначения. Множество натуральных чисел обозначим N, множество целых чисел - Z, а множество вещественных чисел - R1; кроме того, положим No =f {0,1, 2,...}.

Арсентьева Евгения Петровна — аспирант математико-механического факультета Санкт-Петербургского государственного университета. Научный руководитель: доктор физико-математических наук, проф. Ю. К. Демьянович. Количество опубликованных работ: 4. Научные направления: вычислительная математика, математические методы в механике и физике. E-mail: evgeniya. [email protected].

© Е. П. Арсентьева, 2011

На плоскости М2 с декартовой системой координат рассмотрим односвязную область П с гладкой границей дП, задаваемой параметрически векторным уравнением

дП : г = р(в),

где р(в) - 2Ь-периодическая вектор-функция класса С2; в дальнейшем, не нарушая общности, считаем, что в Є [—Б, Б), Б > 0, в - натуральный параметр (т. е. |в| -длина дуги), возрастание которого согласовано с (правой) ориентацией границы дО (как правило, рассматриваем ситуацию лишь в окрестности точки в = 0).

Пусть п(в) - внутренняя единичная нормаль к этой кривой в точке в; таким образом, система ортогональных векторов (р '(в), п(в)) - правая, а поскольку в - натуральный параметр, справедливо соотношение

п(в) = [р '(в)^. (1)

Для любого вектора т введем ортогональный ему вектор , так что |т| = |т^|, а система векторов (ш, ш^) является правой.

Пусть Н(в) - радиус кривизны рассматриваемой кривой в точке в. Введем обозначение є* = ІМ|Д(в)|.

Для точек М(г), отстоящих от кривой дП на расстояние, меньшее числа є*, существуют и единственны числа в Є [—Б, Б) и т Є (—є*,є*) такие, что г = р(в) + тп(в). Множество таких точек обозначим П,

П = {г = р(в) + тп(в) | в Є [—Б, Б), т Є (-є*,є*)},

и будем называть его множеством однозначной представимости.

В дальнейшем радиус-вектор любой точки М будем обозначать гм; для любой упорядоченной пары точек (М, М) на плоскости М2 определим правую полуплоскость Щм,м) для упорядоченной пары точек (М, М) соотношением

^-(м,м) = {А | (гА — гм) • (гN — гм)^ > 0, А Є м2}.

Пусть три точки А, В, С лежат на ориентированной кривой дП; будем говорить, что на ней тройка (А, В, С) положительно ориентирована, если переход от А к В и от В к С происходит в положительном направлении по этой кривой.

Пусть є Є (0,є*); рассмотрим кривые

Хі : г = р(в) + 2-ієп(в), і Є N0,

и положим

П(е) = {г = р(в) + тп(в) | в Є [—Б, Б), т Є (-є, є)}.

Ясно, что область П(е) лежит во множестве П, а кривые Хі разбивают эту область на части

Пі = {г = р(в) + тп(в)|є • 2-і-1 < т < є • 2-і, в Є [-Б, Б)}, і = 0,1, 2,...,

называемые полосами. Заметим, что правая ориентация границы определяет (правую) ориентацию (замкнутых) кривых Хі.

Итак,

П0 = {г = р(в) + тп(в) | є/2 < т < є, в Є [—Б, Б)} - нулевая полоса,

П1 = {г = р(в) + тп(в) | є/4 < т < є/2, в Є [—Б, Б)} - первая полоса

и т. д.

С каждой полосой свяжем триангулированный двусвязный многоугольник, называемый триангулированным слоем (или просто слоем) и обозначаемый Ті, і = 0,1, 2,....

Пусть нулевой триангулированный слой То задан таким образом, что

1) вершины треугольников располагаются лишь на кривых Хо и Хі (у каждого треугольника две вершины лежат на одной из упомянутых кривых, а одна вершина -на другой);

2) число вершин триангуляции на кривой Хо равно числу вершин на кривой Хі; обозначим его Nо;

3) число треугольников триангуляции равно 2Nо, где Nо > 2.

Треугольники, у которых две вершины лежат на кривой Хі, назовем четными треугольниками триангулированного слоя То, а остальные (т. е. треугольники с одной вершиной на Хі ) - нечетными треугольниками этого слоя. Отрезки, соединяющие вершины, лежащие на различных кривых, назовем внутренними, а остальные будем называть внешними отрезками слоя То. Ясно, что

4) число четных треугольников равно числу нечетных треугольников;

5) к четным треугольникам по внутренним отрезкам примыкают нечетные треугольники, а к нечетным по внутренним отрезкам - четные;

6) в рассматриваемом слое из двух треугольников, имеющих общую сторону, один является четным, а другой - нечетным;

7) четные и нечетные треугольники чередуются;

8) внешние отрезки слоя То с вершинами на кривой Хо образуют замкнутую кривую (ее будем называть внутренней многоугольной границей слоя То); внешние отрезки слоя То с вершинами на кривой Хі тоже образуют замкнутую кривую (ее будем называть внешней многоугольной границей рассматриваемого слоя). Обе многоугольные границы имеют правую ориентацию (так же, как кривые Хо и Хі).

Построение триангулированных слоев проведем индуктивно, описав алгоритм построения слоя Ті+і с использованием структуры слоя Т , і = 0,1, 2,.... Заметим сразу, что получаемые в результате триангулированные слои обладают свойствами, аналогичными свойствам 1)-8). В частности, і-й слой составлен из треугольников, вершины которых располагаются лишь на кривых Х^ и Х^+і, причем количества вершин на кривых Х^ и Х^+1 совпадают. Сначала дадим описание алгоритма построения слоя Ті, считая, что слой То со свойствами 1)-8) построен; алгоритм построения слоя Ті+і с использованием структуры слоя Т для і = 1, 2,... может быть описан аналогичным образом.

Пусть ДАВ С - один из четных треугольников слоя То с вершиной В на кривой Хо

и с вершинами А и С на кривой Хі (рисунок). Очевидно, что

єє ГА = 2П(в^) + р(йа), Г в = єп(вв) + р(зв), Г с = 2П(вс) + р(«с), (2)

где в а, вв, вс - соответствующие значения параметра в. Через середину В' отрезка АС проведем прямые, параллельные прямым АВ и ВС до пересечения с кривой Х2; точки пересечения обозначим В'' и В''' соответственно. Очевидны соотношения

гв' = (га + гс)/2, гв = (га + гс)/2+p*(гс - гв),

гв " = (га + гс)/2 + и (гв - га), (3)

где I* и р* - компоненты решений (р*, в*) и ^*,в*) уравнений

(га + ГС)/2 + р*(гс - Гв) = р(в*) + еп(в*)/4,

(га + гс)/2 + г*(гв - га) = р(«*) +еп(в*)/4 (4)

соответственно.

Построим прямолинейные отрезки ЛБ", Б”Б', Б'Б'" и Б'"С. Проделаем аналогичные действия для всех четных треугольников нулевого слоя. Построение триангулированного слоя 71 заканчивается соединением прямолинейными отрезками каждой пары соседних точек, полученных на кривой Л2. Замкнутая ломаная, образованная упомянутыми отрезками, явится внешней многоугольной границей слоя 71; заметим, что внутренней многоугольной границей слоя 71 служит внешняя многоугольная граница слоя

70.

Нетрудно видеть, что свойства 1)-8) выполнены и для слоя 71 (конечно, в них следует заменить Л1 на Л2, Л0 на Л1, 70 на 71 и К0 на К1).

Аналогичным образом по триангулированному слою 71 построим триангулированный слой 72 .

Вообще, по триангулированному слою 7г построим триангулированный слой 71+1, г = 0,1,.... Указанный процесс построения слоя 71+1 будем называть алгоритмом инициализации г + 1 -го слоя (или просто алгоритмом инициализации, если номер появляющегося слоя ясен из контекста).

Четными треугольниками слоя 7г назовем те его треугольники, две вершины которых лежат на кривой Лг+1; остальные треугольники этого слоя (т. е. треугольники с одной вершиной на Л1+1) - нечетными треугольниками слоя 71.

Очевидно, что для слоя 7г выполнены свойства 1)-8) (в них следует заменить Ло на Лг, Л1 на Л1+1, 70 на 7 и К0 на Кг).

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

Предположим, что выполнено следующее условие:

(А) у каждой бинарной пары триангулированного слоя То прямолинейные продолжения боковых сторон пересекаются в правой полуплоскости, определяемой точками (E, C), где E и C - несовпадающие вершины треугольников рассматриваемой бинарной пары на кривой Ai.

Проведем дополнительное построение: продолжим боковые стороны бинарной пары треугольников до пересечения и обозначим буквами 5 (с различными индексами) образованные ими углы. Это построение проделаем для каждой бинарной пары слоя То.

Пусть все упомянутые углы 5, а также углы всех треугольников, составляющих нулевой слой, лежат в интервале (в, п — в), где в - фиксированное число, в G (0,п/4).

Покажем, что при сформулированных далее условиях использование изложенного алгоритма приводит к таким триангулированным слоям, углы треугольников в которых лежат в интервале (в/2, п — в/2).

Рассмотрим бинарную пару (AABC, AEDA) треугольников из слоя То и введем обозначения для соответствующих углов:

a d=f ¿BAC, в = ¿ABC, y = ¿ACB,

ш d=f ¿DEA, ф d=f ¿EDA, ф d=f ¿DAE, 5 d=f ¿(ED,cB).

Отметим, что ¿AB'B" = a, ¿B"B'B"' = в, ¿B"'B'C = Y по построению. Кроме того, рассмотрим углы y ' = ¿B'AB'', в'' = ¿AB''B', a'' d=f ¿B'B''B''', y ''' = ¿B''B'''B', a' d=f ¿B'CB''', в''' = ¿CB'''B', y '' = ¿AB''D''', 5' d=f ¿B''AD''' треугольников AAB''B', AB''B'B''', AB'B'''C, AD'''AB'' из слоя T1, полученных с помощью алгоритма инициализации из четных треугольников AABC, AEDA слоя То.

Для доказательства упомянутой выше невырожденности триангуляции требуются оценки разности углов y и y ', Y'', Y'''; a и a', a''; в и в'', в'''; 5 и 5'; ввиду того, что выкладки весьма громоздки, ограничимся оценкой разности первых двух из них. Используя скалярное произведение векторов, имеем

АС-СВ (гс - гА) ■ (гс ~ г в)

cos7 = ^------==г = ---------¡—¡------------------------(5)

\AC\-\CB\ \rc — га\-\ro — гв \

, АС-АВ" (гс - гА) ■ (гв" - гА)

COS7 = ---г-- ч = --------:-:-------(6)

\Ac\-\AB''\ \rc — га\-\tb » — га\ '

Из (5) и(6) получаем

где

eos y ' — eos y = i— -------------• (R-7 — R-), (7)

\rc — га\

R=,rc % В/= Гд" ГА. (8)

\rc — Гв \ \гв " — га \

3. Триангуляция вблизи прямолинейного участка границы. Пусть рассматриваемый участок границы дП задается параметрическим уравнением

р(в) = Г1 в + Г0, (9)

где в - натуральный параметр, в € (-6, +6), а Г1 и П0 - ортогональные единичные векторы, п0 ± г1, |п0| = 1, |г1| = 1, причем система (г1, п0) - правая.

В этом случае из (2) и (9) получаем

ГА = £П0/2 + Г1вА + Г0, Гв = £Щ + Г1вв + Г0, Г с = еПс>/2 +Г1вс + Г0 (10)

и, следовательно,

гс - Гв = -£П0/2 + Г1(вс - вв), |гс - Гв | = [(е/2)2 + (вс - вв)2]1/2,

так что

_ -еп0/2 + Г1 (вс-«в)

[(е/2)2 + (вс - вв)2]1/2'

Перейдем к отысканию И/ (см. вторую формулу в (8)). Для отыскания разности Гв" - га воспользуемся формулой (3), в которой Ь* отыскивается из уравнения (4); последнее в рассматриваемом случае из-за соотношений (9) и (10) принимает вид

е в А + вс (е вв - вА \

-По н---- ---г 1 + и (^-П0 Н-- ---г 1 j = пв*.

Умножая скалярно на пэ и используя ортогональность векторов пэ и Г1, определяем, что Теперь из (3) имеем г в" — г а = гс~гв, откуда следует, что

гс-гв гс-гв 2 2

Согласно первому из соотношений (8), определяем, что И7 = И, а в силу формулы (7), имеем сов 7' = сов 7; ввиду очевидного соотношения ^,1' € (0, п) приходим к равенству углов 7' = 7.

Аналогично можно получить, что в случае прямолинейного участка границы все углы треугольников слоя 70 равны соответствующим углам треугольников из слоев 7^, построенных с помощью алгоритма инициализации; равными оказываются и углы 6 и 6', являющиеся соответствующими углами четных треугольников бинарной пары. Заметим, что все подобные равенства можно получить аналитическими выкладками, а также непосредственно из геометрических сображений: из равенства углов при пересечении двух параллельных прямых третьей и из подобия треугольников.

Предыдущие рассуждения очевидным образом распространяются и на переход от г-го слоя к г + 1-му: угол, относящийся к г-му слою, оказывается равным соответствующему углу, относящемуся к г + 1-му слою. Выбирая какой-либо угол нулевого слоя и перебирая г = 0,1, 2,..., получаем бесконечную последовательность соответствующих (в данном случае равных) углов; такую последовательность будем называть цепочкой соответствующих углов.

Как видно из п. 2, алгоритм построения триангулированных слоев для криволинейной границы тот же; поэтому понятие цепочка соответствующих углов (теперь, вобще говоря, не равных друг другу) сохраняется и для криволинейной границы. Из полученных в п. 4 результатов следует, что при определенных условиях переход от г-го слоя к г + 1-му ведет к столь малому изменению (уменьшающемуся с ростом г) соответствующих углов, что общая накапливаемая погрешность не выводит углы из промежутка (в/2,п - в/2).

4. Криволинейная граница. Не нарушая общности, будем считать вА = 0, вв,вс € (-п, +п), где п > 0; введем декартову систему координат так, чтобы р(вА) = 0,

p '(sa) = i, [p '{sa)]^ = j, где i и j - направляющие единичные орты осей абсцисс и ординат соответственно.

В дальнейшем (при £ — 0) символы o(£) и о(£) (с различными индексами) используются для обозначения бесконечно малых (векторной и скалярной соответственно) более высокого порядка, чем бесконечно малая £, т. е.

o(£)/£ — 0, о(£)/£ — 0;

символы o(1) и о(1) применяются для обозначения векторной и скалярной бесконечно малых (без указания их порядка).

Предположим, что вектор-функция p(s) достаточно гладкая и для нее справедлива дифференцируемая асимптотика в ^-окрестности нуля, т. е. справедливы представления

p(s) = i • s + r2 ■ s2 + r3 • s3 + ô0(s3), p'{s) = i + 2r2 • s + 3r3 • s2 + ôi(s2),

p"(s) = 2r2 + 6r3 • s +ô2(s), s G (-Г), rj). (11)

Ввиду того, что s - натуральный параметр, справедливо свойство p '(s) ■ p ''(s) = 0; отсюда и из равенств (11) при s =0 имеем i ■ г2 =0, так что можно положить

Г2 = gj,

где g - некоторая константа; итак, предыдущие формулы принимают вид

p(s) = i • s + gj • s2 + r3 • s3 +ü0(s3), p'(s) = i + 2gj ■ s + 3r3 • s2 + ôi(s2),

P"(s) = 2#j + 6r3 • s + ü2(s), s G {—rj, rj). (12)

Замечание 1. В случае прямолинейного участка границы константа g и бесконечно малые равны нулю.

В изучаемом случае из (1) и (2) получаем

г .

ГА = 2-Ь

Гв = p(sB) +e[p'(sB)}±, rc = p(sc) + ^[/,/(sc)]~L- (13)

Пусть выполнено следующее условие:

(Б) переменные sB, sc, s* имеют асимптотические представления

sb = sb,i£ + ов (г), sc = sc,\£ + oc (г), s* = s*,i£ + о* (г),

здесь sb,i, sCji и s*,i - некоторые числа, равномерно ограниченные для всех рассматриваемых углов y, y '.

Лемма. Если выполнено условие (Б), то справедливо представление

t*(e) = —— + t\£ + 0(t*)(e),

где ti = g(sB,i - sc,i)(3sb,i + sc,i)/2.

Далее, проделав математические выкладки, аналогичные выкладкам в прямолинейном случае, и учитывая (12) и (13), получим следующую оценку:

|R' - R| < 8|0*(e)|(1 + 2|sB,i| + o(1)), в которой О* (e) определяется формулой

O*(e) — 2(tie + 0(t*)). (14)

В силу формулы (7), находим

| cos7' — cos7| < 8|O*(e)|(1 + 2|sBji| + o(1)). (15)

Аналогично оцениваются остальные разности между косинусами углов. Из оценок можно получить оценки вида

|7i+i - 7i| < Ce/2i (16)

для цепочки соответствующих углов 70,71,72,..., где 7i - угол, относящийся к г-му слою, и аналогичные оценки для цепочек соответствующих углов ai, f3i, Si, <pi, ^i, wi, г — 0,1,.... Эти оценки позволяют установить следующее утверждение:

Теорема. Если выполнены условия А), Б), то существует число e1 > 0, зависящее от дО. и в, такое, что при e G (0,ei) все углы треугольников полученной триангуляции лежат в интервале (в/2,п — в/2).

Доказательство. Используя (16) в неравенстве

S — i

7 — 70| |7i+i — 7il

i=0

получаем

S—i

7s — 701 < Ce^2 2—i ^ 2Ce,

i=0

откуда и следует доказываемое утверждение. ■

Замечание 2. Для остальных цепочек доказательство аналогично (ввиду стандартных способов их построения достаточно рассмотреть цепочки, порождаемые углами а, в, S, у>, ф, и).

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

Литература

1. Сьярле Ф. Метод конечных элементов для эллиптических задач / пер. с англ. Б. И. Квасова; под ред. Н. Н. Яненко. М.: Мир, 1980. 512 с.

2. Демьянович Ю. К. Симплициальные распространения сеточных функций // Методы вычислений / под ред. З. И. Царьковой. Л.: Изд-во Ленингр. ун-та, 1973. Вып. 8. С. 32—50.

3. Ruppert J. A Delaunay refinement algorithm for quality 2-dimentional mesh generation // J. of Algorithms. 1995. Vol. 18, N 3. P. 548-585.

4. Dorfler W., Rumpf M. An adaptive strategy for elliptic problems including aposteriori controlled boundary approximation algorithm // Mathematics of computation. 1998. Vol. 67, N 224. P. 1361-1382.

5. Ищенко А. В., Киреев И. В. Алгоритм построения двумерных вложенных сеток // Журн. Сибирск. федерал. ун-та. Сер. Математика и физика. 2009. Вып. 2. С. 83-90.

6. Азаренок Б. Н. О построении структурированных сеток в двумерных невыпуклых областях с помощью отображений // Журн. вычисл. математики и матем. физики. 2009. Т. 49, № 5. С. 1-13.

Статья рекомендована к печати проф. Л. А. Петросяном.

Статья принята к печати 14 октября 2010 г.