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

Вычислительные технологии

Том 5, № 2, 2000

АВТОМАТИЗИРОВАННАЯ ТРИАНГУЛЯЦИЯ МНОГОСВЯЗНЫХ ОБЛАСТЕЙ СО СГУЩЕНИЕМ И РАЗРЕЖЕНИЕМ УЗЛОВ *

Ю.В. НЕМИРОВСКИй Институт теоретической и прикладной механики СО РАН

Новосибирск, Россия e-mail: [email protected]

С. Ф. ПЯТАЕВ Институт вычислительного моделирования СО РАН

Красноярск, Россия e-mail: [email protected]

Algorithms of fragmentation of arbitrary piecewise smooth closed boundary contours of multiply connected domain and its triangulation are developed. In the course of fragmentation of the boundary of domain and its triangulation a function of steps is used, which adjusts the sizes of one-dimensional and triangular finite elements according to their position in the domain. Any positive function can appear as the function of steps. The process of triangulation is a consecutive filling of domain with triangular elements. The process of filling starts from the boundary that is preliminarily fragmented into one-dimensional finite elements. Samples of grids for some domains are presented.

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

С целью построения полностью автоматизированного процесса триангуляции произвольных двумерных многосвязных областей в данной работе разработан алгоритм раз-

*Доклад на Международной школе-семинаре “Информационные технологии в задачах математического моделирования”, Новосибирск, 13-18 сентября 1998 г.

© Ю.В. Немировский, С.Ф. Пятаев, 2000.

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

В разделе 2 на основе схемы, предложенной в [4, 5] и представляющей собой последовательное заполнение области треугольными элементами, построен процесс триангуляции области. Начало процесса заполнения идет от границы, предварительно разбитой на одномерные конечные элементы. По мере построения треугольных элементов граница еще не триангулированной области (следуя [4], будем называть ее текущей границей сетки ТГС) представляет собой ряд непрерывных замкнутых ломаных кривых с возможными самопересечениями.

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

h(x,y) = ho +

hi h

o

і=1 і+(A

Ni / y X Ni

+f B

где h0 — основной шаг сетки, n — количество сгущений или разряжений сетки, hi — шаг сетки в центре i-го сгущения, Ai,Bi — размеры области сгущения, Ni — показатель степени, характеризующий градиент сгущения,

Xi = (x — xi) cos ai + (y — yi) sin ai,

Vi = — (x — Xi) sin ai + (y — yi) cos ai,

(хг, у г) — координаты центра г-го сгущения, аг — угол поворота осей г-го сгущения.

1. Разбиение границы многосвязной области

Пусть граница многосвязной области образована N кусочно-гладкими замкнутыми контурами, заданными в некоторой декартовой системе координат Оху в параметрическом виде. Рассмотрим некоторый кусочно-гладкий контур Г (индекс контура опущен), образованный Ь гладкими кривыми 7п, п = 1, ..., Ь, параметрические уравнения которых есть

(x(t), у(г)) = х(^) = хпСО, *_ < * < *+, (1)

где *_, *+ — границы изменения параметра * для 7п. Из условия непрерывности и замкнутости контура Г следует:

Хп(*+) = Хп+1(*-+1), п =1,Ь - 1, Х1 (*-) = Хь(*+).

Параметризация (1) должна быть такой, чтобы для внутренних контуров Г направление обхода при возрастании параметра * было по часовой стрелке, для внешнего контура — против часовой стрелки.

Разбиение Г осуществляется последовательно, начиная с первой гладкой кривой 71: Х = Х1(*). Первый узел на Г есть у1 = Х1(*-). Предположим, что I — 1 первых кривых контура Г уже разбито, последний построенный узел на этих кривых есть упг-1 = Хг_1(*+11) = Хг(*_)

и разбита часть кривой 7г с последним узлом ущ_1+к = Хг (*к), где *к — значение параметра * для последнего узла, *к € [*_,*+).

1.1. Обозначим через вг (*к ,*) длину части кривой 7г, соответствующей значениям *к, *:

«г(*к,*) = / 1 1 ^ € [*к^

•Ч

где Х(*) — производная по * от Х(*).

Новый узел уп1_1+к+1 строится следующим образом: вычисляется значение функции шагов ^(х,у) в последнем построенном узле уп1_1+к и ищется решение *к+1 уравнения

вг(*к,*) = Мупг_1 +к), * € [*к,*+]. (2)

Если решение этого уравнения существует (обратное рассматривается в 1.2), то оно единственно в силу | Хг (*) |> 0. По найденному значению ^+1 вычисляется Хг (?/г+1) и ищется решение *к+1 уравнения

вг (*к,*) = 2[^(упг_1+к) + МХ(4+1))] • (3)

Предположим, что это уравнение также имеет решение (обратное рассматривается в 1.3). Рассмотрим неравенство

— в(4,*к+1) л ^ ^ ^

1 ---ст—1—^ \< е ^ =| Ущ-1+к — Х(*к+1) ^ (4)

вг ( к, *к+1)

левая часть которого есть относительная разность между длиной дуги и длиной d соответствующего ей отрезка и которое характеризует отклонение дуги от отрезка прямой. Значение параметра е задается пользователем (например, е = 0.01). Если неравенство (4) выполнено, то точку Хг(*к+1) объявляем новым узлом уп1_1+к+1 и переходим к построению следующего узла. При невыполнении неравенства (4) правую часть уравнения (3) последовательно уменьшаем на некоторую величину (например, на одну десятую от правой

части) до тех пор, пока (4) не будет выполнено. Эта ситуация возникает, когда длина d од-

номерного элемента, вычисленная в соответствии с функцией шагов, достаточно “велика” для приемлемой аппроксимации этим элементом соответствующей ему дуги, и поэтому производится последовательное уменьшение этой длины до требуемой неравенством (4). Затем переходим к построению следующего узла ущ_1+к+2.

Окончание процедуры построения новых узлов на /-й кривой контура Г связано с отсутствием решения уравнения (2) и описано в 1.2.

1.2. Уравнение (2) не имеет решения, т.е. вг(*к,*+) < ^(упг-1 +к). Это означает, что последний узел уп;_1+к “близок” к Хг(*+) и построение нового узла через Л,(х, у) невозможно. Обозначим через 8г длину остатка 7^, через dk длину последнего построенного элемента:

81 в1 (*к, *1 ), ^ 1 упг_1+к_1 — Ущ_1+к \ .

Если невязка $г удовлетворяет неравенству

$г < е^к,

(5)

где е1 0.1, то последний узел уп1_1+к смещаем в крайнюю точку 7г, т. е. ’ущ_1+к = Хг(*+) и

переходим к построению узлов на следующей кривой 7г+1. При невыполнении (5) рассма-

триваем неравенство

$г < 0•5dk. (6)

Если (6) имеет место, то число узлов на остается прежним, невязка $г разносится по

всем построенным на дугам пропорционально их длинам. Если обозначить через в

длину пройденной при построении узлов части кривой 7г

к_ 1

в = 5] вМ+1> вМ+1 = вг(^*т )> *0 = *_> (7)

г=0

то длины новых дуг в*г^+1 определяются через длины старых дуг в^ ^+1 по формулам

в*;+1 = в!;г+1(1 + 81/в1), г = 0,к - 1. (8)

Разрешая последовательно уравнения

вгК'■**+ ,) = в*‘+1. *0г = *_, г = 01-1, (9)

определяем новые значения *пг, полагая *кг = *+, и вычисляем новые узлы на 7г:

уп1_ + = Хг(**г). г = 0. к - 1. упг_1+к = Хг(*+). (10)

после чего переходим к построению узлов на следующей кривой 7г+1.

Если невязка $г не удовлетворяет условию (6), то рассматриваем новое неравенство

8г < 4. (11)

Если это неравенство выполнено, т. е. длина неотработанной части 7г меньше длины последнего построенного одномерного конечного элемента, но больше его полудлины в силу невыполнения (6), то число узлов на 7г увеличивается на единицу, вводится новая невязка $г:

к_1

8г = вг - в , 5г = ^ в^,г+1 + вк_1,к,

г=0

где вг — длина 7г, и переходим к (8) - (10) с новыми $г, 5г и с дополнением еще одной дуги, длина которой вк к+1 равна длине последней построенной дуги вк_1 к. При этом необходимо в (8) — (10) увеличить значение к на единицу.

В случае, когда (11) не выполнено, рассматриваем уравнение

в г(*к,*) = 4. * € [*к.*+]. (12)

решение которого существует в силу

в г(*к,*+) = $г > dk, в г(*к,*к) = 0 < dk.

Найденное из (12) значение *гк+1 проверяется на выполнение неравенства (4) и далее по алгоритму с той лишь разницей, что в случае невыполнения (4) последовательно уменьшается правая часть не (3), а (12).

1.3. Уравнение (3) не имеет решения. После нахождения ?к+1 из уравнения (2) переходим к неравенству (4), в котором *к+1 заменяется на ?/г+1. В случае его выполнения точку Хг(?к+1) объявляем новым узлом упг_1+к+1 и переходим к построению нового узла на 7г, в противном случае последовательно уменьшаем правую часть уравнения (2) до тех пор, пока не будет выполнено (4), после чего переходим к построению следующего узла.

2. Триангуляция области

2.1. Найдем звено гт1п ТГС, имеющее минимальную длину /т1п и узлы Жт1п,Жт1п- Обо-

— +

значим через гт1п, гт1п соответственно предыдущее и последующее по отношению к 2т1п звенья- Выберем из 2—1п , 2+1п звено ^т, образующее с 2т1п наименьший угол вт1п (вт1п = ш1п(в1,в2); углы в и в2 измеряются против часовой стрелки от 2т1п к 2—1п и от 2+1п к 2т1п соответственно). Обозначим узлы выбранной пары звеньев (это либо г—1п,гт1п либо гт1п,2+1п) через ХГп,Хтт ХТ- Если вт1п < 80° (в противном случае переходим к процедуре 2.4), то:

2.2. Проводим проверку на попадание узлов ТГС в треугольник Д(гт1п, гт1п)- Если таких узлов нет, то переходим к процедуре 2.3, в противном случае из всех попавших в Д(гтт, 2”1п) узлов выбираем узел у*, наиболее близкий к 2т1п, и переходим к процедуре 2.122.3. Рассмотрим круг радиуса ^ | ж”1™ — Ж”11™ | с центром жс, жс = (ж”1™ + Жт1п)/2- Если во внешнюю по отношению к треугольнику Д (гт1п, 2”1п) половину круга узлы ТГС не попадают, то переходим к процедуре 2.13- В противном случае из всех попавших узлов выбираем наиболее близкий к отрезку [ж”1™, ж”1™] узел жт и четырехугольник □(гт1п, 2”1п, Хт) разбиваем на два треугольника так, чтобы минимальный угол получающихся треугольников был максимален- Оба полученных треугольника объявляем элементами, из ТГС удаляем 2т1п, гт1п, определяем связность области, добавляем два новых звена [ж”1™, жт], [жт, ж”1™] и переходим к процедуре 2.1-

2.4. Построим точку ж*:

х* = +кп = 2(хт1п+^ (13)

где п — нормаль к 2т1п, направленная внутрь области:

(ж^п — Х”1п) X п = 1т1пез, ез = (0, 0, 1),

1 3 ^3

кср = 3^ ' М^О, О = жт 1п, 3 = 1, 2; у3 = Жт1п + 2^ 1т1пп-

г=1

На базе звена гт1п построим прямоугольник П, одна из сторон которого есть гт1п, а другая направлена по нормали и ее длина равна 2к, где к — высота в Д(гт1п , ж*), опущенная на 2т1п- Установим с помощью контрольной области П критерии близости строящегося нового узла ж* к ранее построенным узлам и звеньям- Если ж* в определенном смысле близок к узлам или звеньям, то от построения нового узла отказываемся и из близких узлов для построения нового элемента выбираем наилучший-Определим два множества М0 и М1 :

М0 — множество номеров узлов ТГС, попавших в П, за исключением номеров узлов

12 ж 1. их2. :

т1п т1п

М1 — множество номеров звеньев ТГС, пересекающих дП, за исключением номера минимального звена 2т1п:

./1^1 — {п * П дП = 0, = 2т1п}

Введем две дополнительные точки 21 и 22:

+ (—1)*Д/ Т, г — 1, 2,

" (жт1п Жт1п), Д/ Т(к(жт1п) /т1п)-

т1п т1п 1т1п 4

Определим в равнобедренном треугольнике Д(ж*, 21, 22) угол С при основании [21, 22] и положим у — 21, если г”1п — 2—1п, и 2^ — 22, если г”1п — 2++1п (см- 2.1)-Если М0 — 0 или М1 — 0, то переходим к процедуре 2.5-

Определим угол а1 при вершине ж* в треугольнике Д (21, 22, ж*)- Если а1 > 30° и вт1п — С > 20°, то точку ж* объявляем новым узлом и переходим к процедуре 2.14- Если а1 < 30°, то точку ж* переопределяем так, чтобы новый угол а1 — 30°, ж* — ж”1п + 122 — ж”1п | ■ tg75° ■ п. При этом если вт1п — 75° < 20°, то переходим к процедуре 2.2, иначе — к процедуре 2.14-

2.5. Введем целый параметр переключения /ЖД и положим /ЖД — 0- Если М0 — 0 (иначе переходим к процедуре 2.9), то из М0 выбираем номер т, при котором соответствующий узел жт наиболее близок к 2т1п- Для этого определяем расстояния от точек Ж, г £ М0 до 2т1п и выбираем

/т — ш1п —> т. (14)

г€Мо

Рассмотрим в треугольнике Д(жт, 21, 22) угол а при вершине жт (у определены в 2.4)-Если а > 30° (иначе переходим к процедуре 2.8), то проводим проверку на пересечение отрезка [жт,жт1п] со звеньями ТГС, имеющими номера из М1 без номеров звеньев, примыкающих к узлам жт,Жт1п - Обозначим для удобства жт через у*-

2.6. Если пересечений нет, то переходим либо к процедуре 2.12 (при /ЖД — 0), либо к процедуре 2.14 (при /ЖД — 1)- В противном случае переходим к процедуре 2.72.7. С использованием узлов пересекающего звена 2р строим ориентированные против часовой стрелки треугольники Д(2т1п,ж&),ж& — узлы звена 2р, к £ {к1,к2}- При построении треугольников необходимо следить за их существованием- Из этих треугольников (если существуют оба) выбираем Д(2т1п,ж^4), к £ {к1, к2}, у которого минимальный угол больше (в дальнейшем минимальный угол в каком-либо треугольнике Д(2,ж) будем обозначать через а(2,ж))- Обозначив выбранный узел ж^ через у*, проводим проверку на пересечение обеих боковых сторон Д (2т1п, у*) со всеми звеньями ТГС за исключением минимального и звеньев, примыкающих к минимальному углу и к узлу у*, и переходим к процедуре 2.62.8. Ближайший узел жт достаточно далек от 2т1п (так как а < 30°), поэтому ж* смещаем к 2т1п так, чтобы расстояние от ее нового положения (обозначим эту точку через

у*) до 2т1п было равно /т/2; у* — ж”1п + ^/тп; /т определено в (14)- В треугольнике

Д(21, у2, у*) определяем угол С при вершине 2^; у определено в 2.4- Если вт1п — С < 20°, то переходим к процедуре 2.2, иначе проводим проверку на пересечение одной из боковых сторон Д(2т1п,у*) со звеньями, имеющими номера из М1, и переходим к процедуре 2.6, положив /ЖД — 1-

2.9. Если вт1п — С < 20°, то переходим к процедуре 2.2, в противном случае из всех звеньев, пересекающих дП, выберем 2т1 и 2т2, которые наиболее “близки” к 2т1п- Для этого рассмотрим все точки пересечения у1 звеньев 2^, к^ £ М1, с боковой стороной Г1 прямоугольника П, проходящей через ж”1п, и аналогичные точки у2 для Г2- Тогда номера т1 и т2 определяются как числа из М1, доставляющие минимум выражению |у^ — Ж”^,

т- е-

тт у — ж^п! —> т, з — 1 Ар.

кг€М1

Если Г1 и Г2 пересекают различные “ближайшие звенья”, то — 2- Если одну из боковых сторон Г звенья не пересекают или Г1 и Г2 пересекаются одним и тем же звеном, то — 1 и рассматриваем только номер т1- Рассмотрим ориентированные Д(2т.,Ж*), г — 1,Жр; Ж* определено в (13)- Если таких ориентированных треугольников не существует (что возможно только в случае, когда “ближайшее” пересекающее Г1 и Г2 звено единственно и проходит через Г1 и Г2 между 2т1п и ж*), то звено 2т1 обозначаем через и переходим к процедуре 2.7- При существовании Д(2т. ,ж*) рассмотрим углы а^ при вершине ж* в этих треугольниках, г — 1, Жр- Если а^ < 90°, то переходим к процедуре 2.14, в противном случае выберем а^1 — шах^=у^р а^-

2.10. Построим ориентированные треугольники Д(2т1п,жт. ), где Жт. — узлы звена

г1 г1

- Список значений параметра з может быть одним из следующих: если существуют оба треугольника, то з £ {1, 2}, если существует только один треугольник, то з — 1 при использовании узла ж^. и 3 — 2 для второго узла звена 2тг1 - Выберем из Д (2т1п, ж^. ) тот треугольник, у которого минимальный угол больше, а(2т1п,жЛ. ) > а(2т1п,Жтг ) для всех 3 из списка значений этого индекса, и переходим к процедуре 2.112.11. Если а(2т1п,жЛ. ) < а(2т., Ж*), то переходим к процедуре 2.14, в противном

г1 _ г

случае боковые стороны Д(2т1п,жЛг ) проверяем на пересечение со всеми звеньями ТГС

за исключением звеньев, примыкающих к узлу Ж^. и к звену 2т1п- Кроме того, проводим

г1

проверку на попадание узлов ТГС в этот треугольник за исключением узла Ж^. и узлов

г1

Ж”1п, Ж”1п - Если пересечений и узлов внутри треугольника нет, то узел Ж^ переобозначим как у* и перейдем к процедуре 2.12-

Если пересечения имеются или треугольник содержит хотя бы один узел ТГС, то из списка значений параметра з (если он не исчерпан) выбираем новое значение з2 и переходим в начало процедуры 2.11, переобозначив предварительно з2 на з\- Если список параметра 3 исчерпан, то рассматриваем второе пересекающее звено 2т.2 (при условии, что список параметра т^ не исчерпан, г — 1, Ар), и если а^2 > 90°, переходим к процедуре 2.10, переобозначив предварительно г2 на г1- В противном случае (либо а^2 < 90°, либо исчерпан список параметра т^) переходим к процедуре 2.142.12. Треугольник Д(2т1п,у*) объявляем элементом, из ТГС удаляем 2т1п, определяем число связности области, добавляем два новых звена [Жт1п,у*], [у*,Жт1п] и переходим к процедуре 2.12.13. Треугольник Д (2т1п, 2”1п) объявляем элементом, из ТГС удаляем 2т1п, 2”1п, добавляем одно звено [ЖТ\жГ] и переходим к процедуре 2.1.

2.14. Треугольник Д(2т1п,Ж*) объявляем новым элементом, из ТГС удаляем 2т1п, добавляем два новых звена [жт1п,ж*], [ж*,жт1п] и переходим к процедуре 2.1-

Заключение

В качестве примеров работы алгоритмов на представленных ниже рисунках показаны сетки для различных областей-

На рис- 1 приведено кольцо с двумя круговыми вырезами- Сетка имеет 864 элемента и 483 узла- Для детального рассмотрения области сгущения около окружности радиуса 0-05 из сетки была вырезана часть, изображение которой представлено на рис- 2-

Рис. 1.

0.582 0.710

Рис. 2.

Рис. 3.

Рис. 4.

На рис- 3 и 4 показаны разбиения одной области с различными параметрами функции шагов, отвечающими за степень сгущения- Центром области сгущения является центр эллипса- Сетка на рис- 3 имеет 880 элементов и 489 узлов, на рис- 4 — соответственно 1768

и 954.

На рис. 5 представлена сетка с 1219 элементами и 702 узлами. Областью является шестерня с N = 20 зубьями. Центры областей сгущения для зубьев расположены в вершинах зубьев. Для детального рассмотрения сетки на зубе был вырезан один из зубьев, и увеличенная сетка показана на рис. 6.

Рис. 7.

Рис. 8.

Рис. 9.

На рис. 8 приведена сетка для прямоугольника с клиновидным и треугольным вырезами, имеющая 1348 элементов и 726 узлов. Приведенная сетка имеет две области сгущения; центром первой области является центр треугольного выреза, центром второй — вершина клиновидного выреза. На рис. 9, 9 приведены увеличенные вырезы сетки в районе вершины треугольного выреза и в районе вершины клиновидного выреза соответственно.

Список литературы

[1] КАМЕЛЬ Х. А., ЭЙЗЕНШТЕЙН Г. К. Автоматическое построение сетки в двух- и трехмерных составных областях. В “Расчет упругих конструкций с использованием ЭВМ”. Ленинградский гос. ун-т, 1974, 21-35.

[2] Квитка А. Л., Ворошко П.П., Бобрицкая С. Д. Напряженно-деформированное состояние тел вращения. Наук. думка, Киев, 1977.

[3] УмАнский С. Э. Алгоритм и программа триангуляции двумерной области произвольной формы. Проблемы прочности, 6, 1978, 83-87.

[4] МильковА Н. И. Особенности дискретизации области при решении задач концентрации напряжений методом конечных элементов. Машиноведение, 2, 1979, 67-71.

[5] Сакало В. И., Шкурин А. А. Универсальная программа триангуляции двумерной области произвольной формы со сгущениями сетки. Проблемы прочности, 1, 1985, 106-108.

Поступила в редакцию 16 марта 1999 г., в переработанном виде 8 июня 1999 г.