Построение диаграммы Вороного и определение границ области в методе естественных соседей Текст научной статьи по специальности «Математика»

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

Том 13, № 3, 2008

Построение диаграммы Вороного и определение границ области в методе естественных соседей

С.Н. КАРАБЦЕВ, С. В. СТУКОЛОВ Кемеровский государственный университет, Россия e-mail: skarab@kemsu.ru, serg@kemsu.ru

In the present work, we present the algorithms for a fast construction of the Voronoi diagram and Delaunay triangulation. For determination of a boundary of the numerical domain, the alpha-shape algorithm is used. A comparison of the computational speed of the proposed algorithm and its theoretical estimation is carried out. The proposed algorithms are applied to a natural element method for solving fluid dynamics problems with free surfaces.

Введение

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

Большую группу составляют методы, в основе которых лежит дискретизация расчетной области конечными элементами. Наиболее известные представители этой группы — метод конечных элементов (МКЭ) [3] и метод конечных объемов (МКО) [4]. Указанные методы дают возможность проводить расчеты на областях с произвольной геометрией и сложными граничными условиями. Основным недостатком МКЭ и МКО является то, что для каждого временного шага сетка, на которой строится решение, сохраняет свою узловую связность, что в свою очередь может привести к ее вырожденности в задачах с подвижными границами. Для преодоления указанной трудности требуется перестройка сетки на каждом временном шаге. Например, корректное разбиение пространства на треугольные конечные элементы может производиться различными способами и является неединственным [5]. Неединственность триангуляции приводит к

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2008.

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

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

На протяжении последних трех десятилетий разрабатываются новые численные методы, которые аппроксимируют уравнения в частных производных, основываясь только на наборе узлов, без знания дополнительной информации о структуре сетки. Данные методы известны как бессеточные методы (БМ) [9, 10]. Основным преимуществом их является возможность проводить численные расчеты задач гидродинамики на многосвязных областях и в областях с сильной деформацией границ течений. К этой группе относятся такие методы, как метод сглаженных частиц (Smoothed Particle Hydrodynamics, SPH) [11], полунеявный метод движущихся частиц (Moving Particle Semi-implicit, MPS) [12, 13], метод лагранжево-эйлеровых частиц [14]. К недостаткам бессеточных методов можно отнести их сравнительно невысокую точность и трудность введения граничных условий. Дополнительно, в стандартных бессеточных методах для построения интерполяционных функций требуется обеспечить узловую связность. Поэтому стали разрабатываться методы, сочетающие в себе сеточный и бессеточный подходы и использующие преимущества каждой из методологий. Характерные представители этой группы методов — это бессеточный метод конечных элементов (Meshless Finite Element Method, MFEM) [15] и метод естественных соседей (Natural Element Method, NEM) [16, 17]. Методы MFEM и NEM сохраняют преимущества классического метода конечных элементов, такие как простота функции формы в области определения, непрерывность между элементами, простота введения граничных условий. При этом указанные методы имеют все преимущества бессеточных методов в силу того, что функции формы MFEM и NEM зависят только от положения узловых точек.

Идея методов NEM и MFEM заключается в применении интерполяций Сибсона [18] и Лапласа [19, 20], построения осуществляется через координаты естественных соседей. В вычислительной геометрии концепция естественных соседей связана с понятием ячеек Вороного первого и второго порядков — выпуклого многоугольника, определяемого множеством точек на плоскости (в пространстве). Дискретизация области ячейками Вороного обеспечивает быстрый поиск естественных соседей для заданной точки. Нахождение естественных соседей на основе информации о смежных узлах, полученных из диаграммы Вороного, существенно снижает временные затраты, необходимые для построения интерполирующих функций. Скорость выполнения этого шага является критической для NEM и MFEM, так как данные методы, основываясь на методе Галеркина, требуют вычисления функций формы для каждой точки при численном интегрировании по расчетной области задачи. Соответственно, скорость работы метода естественных соседей и бессеточного метода конечных элементов зависит от эффективности алгоритма, реализующего разбиение области ячейками Вороного, что особенно

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

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

1. Постановка задачи

Численные методы, применяемые в механике сплошных сред, требуют дискретизации области течения множеством материальных точек, представляющих некоторый объем жидкости согласно гипотезе о сплошности среды. Пусть на плоскости в пространстве Я2 задано множество материальных точек Р = {р1,р2, ...,рп}. Граница течения задается последовательным набором материальных точек. Метод естественных соседей и бессеточный метод конечных элементов являются лагранжевыми методами, поэтому при решении нестационарных задач на каждом шаге по времени необходимо определять новое положение свободных границ и получать разбиение области ячейками Вороного и треугольниками Делоне.

Основные определения

Диаграмма Вороного — геометрическое разбиение области на многоугольники, обладающие следующим свойством: для любого центра рг системы точек можно указать область пространства, все точки которой ближе к данному центру, чем к любому другому центру системы. Такая область называется многогранником Вороного, или областью Вороного (рис. 1, а). Строгое определение многоугольника Вороного Тг вводится следующим образом:

Тг = {х € Я2 : d (х, хг) < d (х, х^) V? = г}. (1)

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

Рис. 1. Диаграмма Вороного (а); двойственность триангуляции Делоне и диаграммы Вороного (б)

Используя формулу Эйлера для планарного полносвязного графа, можно доказать, что при числе узлов n > 3 диаграмма Вороного содержит не менее 2n — 5 вершин и не менее 3n — 6 ребер [21, 22].

Важным свойством диаграммы Вороного является ее двойственность к триангуляции Делоне. Для того чтобы получить триангуляцию Делоне, необходимо соединить отрезками все пары точек, многоугольники Вороного которых имеют общее ребро (рис. 1, б) [5, 21]. Определение естественных соседей для вычисления интерполяций Сибсона и Лапласа вводится через разбиение Вороного и триангуляцию Делоне: для ячейки Вороного Ti, xi £ P натуральные соседи для xi £ P есть вершины треугольников Делоне, инцидентных к xi.

Для построения диаграммы Вороного существует множество методов. Самый простой и известный — метод пересечения серединных перпендикуляров (полуплоскостей), идея которого заключается в поочередном построении многоугольников, входящих в диаграмму [21]. Для двух точек pi и pj множество точек, более близких к pi, чем к Pj, есть полуплоскость H (pi ,pj), определяемая прямой, перпендикулярной к отрезку pipj, и делящая его пополам, и содержащая точку pi. Множество точек, более близких к pi, чем к любой другой точке из множества, получается в результате пересечения

n — 1 полуплоскостей: T (i) = П H (pi ,Pj). Временная сложность алгоритма составляет

i=j

O (n2 log n).

Следующий класс методов основывается на инкрементальном алгоритме, асимптотическая оценка которого составляет O (n^/n) [23]. Суть метода состоит в следующем: строится диаграмма Вороного для некоторого начального числа точек, затем в нее поочередно вставляются оставшиеся из множества точки и перестраивается лишь часть диаграммы.

С развитием вычислительной геометрии был разработан более сложный, но превосходящий оценку O (n log n) алгоритм. Этот алгоритм основан на принципе "разделяй и властвуй" ("devide and conqueror") [21]. Множество точек P делится на два приблизительно равных подмножества Pi и P2. Рекурсивно строится диаграмма Вороного для множеств P1 и P2, а затем с помощью "разделяющей цепи" восстанавливается диаграмма для исходного множества P. Несмотря на всю привлекательность данного метода, он вызывает большие трудности при численной реализации.

В 1987 году Стивен Фортун [24] на основе метода "sweep line" (метод заметающей плоскости) предложил алгоритм для вычисления диаграммы Вороного, временная сложность которого составляет O (n log n). Данный метод быстрый по скорости работы и относительно простой в реализации.

2. Идея метода "sweep line"

Группа методов "sweep line" нашла свое применение во многих практических задачах вычислительной геометрии. Все они сводятся к задаче построения графа определенного вида. Идея методов "sweep line" заключается в перемещении горизонтальной линии "sweep line" l ("линии зачистки") сверху вниз по плоскости. Во время движения такой линии формируется информация о графе, который необходимо вычислить. Вид найденного графа меняется при достижении "линией зачистки" специальных точек на плоскости — событий. При построении диаграммы Вороного необходимо сохранять точки пересечения диаграммы с линией "sweep line". Трудность заключается в том, что вычис-

ленная на текущий момент часть диаграммы зависит от точек pi Е P, которые могут находиться как выше, так и ниже прямой l. При достижении "линией зачистки" самой верхней вершины ячейки Вороного невозможно определить узел pi Е P, вокруг которого строится данная ячейка. Для решения указанных проблем С. Фортун предложил модифицированный вариант метода "sweep line" построения диаграммы Вороного на плоскости [24]. В данной работе описывается несколько измененный вид этого алгоритма [22]. Вместо того чтобы искать диаграмму Вороного через пересечение с линией l, строится только та часть диаграммы, на которую точки, находящиеся ниже l, не оказывают влияния.

Пусть задано множество точек на плоскости P = (pi,p2, ...,pn}. Для преодоления проблемы взаимного влияния точек, расположенных по разные стороны от прямой l, необходимо разделить полуплоскость, лежащую над "линией зачистки", на две области. Одна область содержит точки, расположенные над "sweep line", находящиеся ближе к некоторой точке pi Е P, чем к самой "sweep line", другая — точки, расположенные ближе к "sweep line", чем к некоторой точке из заданного множества P. Множество точек q, равноудаленных от "sweep line" и самой близкой точки pi Е P, называется "береговой линией" ("beach line"). По определению, множество точек, равноудаленных от заданной точки pi Е P, лежащей выше горизонтальной линии, и самой линии, образует параболу. Таким образом, "береговая линия" представляет собой последовательность параболических дуг. Уравнение параболы, относительно фокуса F (xf,Vf) и директрисы y = L, записывается в следующем виде:

Из (2) видно, что парабола, ветви которой направлены вверх, сужается по мере приближения точки p к "линии зачистки" и постепенно вырождается в вертикально направленный луч, выходящий из точки p.

Точка пересечения двух парабол, называемая контрольной точкой (breakpoint), равноудалена от точек pi,pj Е P и "береговой линии". Следовательно, она лежит на ребре ячейки Вороного, проходящего между точками pi и pj Е P. Таким образом, алгоритм Фортуна состоит в моделировании роста "береговой линии" во время движения "линии зачистки" сверху вниз и отслеживании точного пути, пройденного контрольной точкой, перемещающейся вдоль ребра ячейки Вороного.

При движении "линии зачистки" вниз параболы, формирующие "береговую линию", изменяются непрерывно. Как и во всех алгоритмах, использующих идею метода "sweep line", необходимо проводить дискретное моделирование движения "линии зачистки" по точкам плоскости, в которых возникают события, изменяющие структуру "береговой линии". Эти события могут быть двух видов:

событие точки (site events): "линия зачистки" достигает новую точку pi Е P .В результате новая дуга вставляется в "береговую линию", появляется две новых контрольных точки;

событие круга (circle events, vertex events — событие вершины): длина дуги параболы стремится к нулю вследствие сближения двух контрольных точек. В результате дуга исчезает с "береговой линии", а в точке столкновения образуется вершина ячейки Вороного.

Реализация метода "sweep line" включает в себя обработку описанных выше типов событий.

(2)

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

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

В работе [22] доказано следующее утверждение: новая парабола может появиться на "береговой линии" только после возникновения и обработки "события точки". Следствием этого является тот факт, что "береговая линия" состоит из 2к — 1 параболических дуг, здесь к — число точек, находящихся выше "линии зачистки", для которых еще не построена замкнутая ячейка. Основным преимуществом "событий точки" является то, что они известны до построения самой диаграммы. Для определения этих событий необходимо отсортировать по убыванию значения ординаты массива точек рг £ Р.

Событие круга. В отличие от "события точки" "событие круга" не известно априори: оно генерируется динамически по ходу работы алгоритма. Данный тип события

Рис. 2. Возникновение "события точки": а — прямая I не достигла точки рг; б — прямая I достигла точки рг, дуга "береговой линии" разделилась на две части бесконечно малой дугой; в — прямая I пересекла точку рг, дуга на "береговой линии" расширяется

а

в

Рис. 3. Отслеживание ребра диаграммы контрольными точками

Рис. 4. Возникновение "события круга": a — "sweep line" не достигла нижней точки окружности; б — "событие круга", образование вершины q диаграммы; в — отслеживание нового ребра диаграммы

связан с удалением дуги из "береговой линии" при уменьшении длины дуги до точки вследствие движения пары контрольных точек навстречу друг другу. Каждое "событие вершины" формируется объектами, являющимися соседними на "береговой линии". События генерируются тройкой точек p, pj, pk, параболические дуги которых последовательно расположены на "береговой линии" слева направо (рис. 4).

Пусть а' — дуга, которая в следующий момент времени будет удалена, а а и а" — соседние с ней дуги, не являющиеся частями одной параболы. Эти три дуги определяются тремя различными точками pi,pj,pk. В момент времени, когда дуга а' исчезает с "береговой линии", все три параболы проходят через общую точку q, т. е. точка q равноудалена от прямой l и трех заданных точек. Тогда существует окружность с центром в точке q, проходящая через pi, pj, pk, причем нижняя точка окружности касается прямой l. В данной окружности не существует ни одной точки из P: такая точка была бы ближе к q, чем q к l, тем самым противореча тому факту, что q находится на "береговой линии". Отсюда следует, что q есть вершина ячейки Вороного. В момент удаления дуги с "береговой линии" две контрольные точки сталкиваются, соединяя два новых ребра диаграммы, так как контрольная точка движется по ребру ячейки Вороного. В работе [22] приведено доказательство, что исчезновение дуги из "береговой линии" может происходить только при возникновении "события круга".

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

3. Реализация алгоритма

Структура "береговой линии" меняется в зависимости от последовательности поступающих событий: при возникновении "события точки" новое ребро диаграммы начинает расти, при возникновении "события круга" два растущих ребра сталкиваются, образуя

Рис. 5. Представление ячейки Вороного вокруг точки pi двусвязным списком ребер

вершину ячейки Вороного. Для реализации алгоритма необходима структура данных, способная хранить части диаграммы, вычисленные на текущий момент времени. Как правило, для большинства алгоритмов метода "sweep line" используются список с приоритетом для хранения очереди событий и некоторая динамическая структура для описания состояния "береговой линии" — двухсвязный список, или двоичное дерево поиска.

Строящаяся диаграмма Вороного хранится в двусвязном списке ребер DCEL (Double-Connected Edge List) [22] (рис. 5). Каждый элемент списка имеет указатель на узел p Е P (табл. 1), указатель на начальную и конечную вершины ребра (табл. 2), указатель на ребро, являющееся соседним для данного ребра, но направленным в противоположную сторону (табл. 3). Ребра ячейки, построенной вокруг узла p Е P, ориентированы против хода часовой стрелки, что дает дополнительные преимущества при реализации численных методов, подобных МКЭ. Для получения данных о ячейке Воро-

Таблица 1. Массив материальных точек

Точка Координаты точки Произвольное ребро ячейки

p1 (5, 7) e2.3

p2 (3, 5) e4.1

p3 (7, 8) e2.1

Таблица 2. Вершины диаграммы Вороного

Вершина Координаты вершины

v1 (8, 4.5)

v2 (5.5, 9)

v3 (1, 9)

Таблица 3. Структура DCEL

Ребро Начало Конец Близнец Инцидентный узел Следующее ребро Предыдущее ребро

e1.2 v1 v2 e2.1 p1 e2.3 e3.1

e2.3 v2 v3 e3.2 p1 e3.1 e1.2

e3.1 v3 v1 e3.1 p1 e1.2 e2.3

e2.1 v2 v1 e1.2 p3 e1.5 e6.2

e1.3 v1 v3 e3.1 p2 e3.7 e4.1

Рис. 6. Представление "береговой линии" двоичным деревом поиска

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

"Береговая линия" представляется сбалансированным двоичным деревом поиска T (рис. 6). Его листья соответствуют дугам на "береговой линии": крайний левый лист представляет крайнюю левую дугу, следующий лист представляет вторую крайнюю левую дугу, и т. д. Для реализации алгоритма нет необходимости хранить параболы в явном виде, так как каждый лист дерева содержит указатель на узел pi £ P, представляющий соответствующую дугу, которую можно вычислить по текущему положению "линии зачистки". Внутренние узлы дерева T соответствуют контрольным точкам на "береговой линии" и представляются в виде упорядоченной пары < pi,pj >, где pi обозначает параболу, располагающуюся слева от точки, а pj — справа. При использовании такого представления "береговой линии" поиск дуги, располагающейся непосредственно над новой точкой, занимает O (log n) операций. Для этого необходимо сравнивать значение координаты по оси абсцисс "события точки" со значением абсциссы контрольной точки, которая может быть вычислена из пары < pi,pj > и текущей позиции "sweep line".

В узлах дерева хранятся указатели на структуры, используемые во время "зачистки". Каждый лист T, представляющий дугу а, хранит один указатель на узел в очереди событий. Этот узел соответствует "событию круга", при возникновении которого дуга а будет удалена. Дополнительно каждый внутренний узел содержит указатель на ребро диаграммы, которое отслеживается

Очередь событий Q представляется очередью с приоритетом, где элементы очереди сортируются по убыванию значения ординаты. Очередь хранит два типа событий: "события точки", которые известны заранее, до начала работы алгоритма, и "события круга", которые вычисляются по мере движения линии вниз. Для "события точки" хранятся координаты самих точек. Для "события круга" хранятся координаты самой нижней точки окружности, а также указатель на лист дерева T, представляющий дугу, которая должна исчезнуть.

Следующий псевдокод описывает алгоритм основной программы для построения диаграммы Вороного методом "sweep line".

1. Инициализация очереди событий Q значениями координат точек pi £ P.

2. [ЦИКЛ ПОКА] очередь не пуста [ДЕЛАТЬ].

3. Получить событие из очереди с максимальным значением координаты у.

4. [ЕСЛИ] событие qi является "событием точки', [ТО] обработать_событие_точки (qi).

5. [ИНАЧЕ] обработать_событие_круга (qi).

6. Удалить событие qi из очереди Q.

7. [КОНЕЦ ЦИКЛА ПОКА].

4. Алгоритм определения границы области

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

Эффективным методом определения границ области является метод а-8Ьаре [25]. Пусть задано множество точек Р = {'1,'2, ■■■,'Рп} на плоскости. Задача нахождения границы множества подразумевает распознавание его формы. Данная задача не имеет единственного решения. Авторы работы [21] упоминают, что Лагу1э был первым, кто затронул проблему вычисления выпуклой оболочки множества точек Р — границы наименьшей выпуклой области, которая охватывает Р. Математические определения формы были введены позднее Х. Эдельсбруннером [25]. Он предложил естественное обобщение для выпуклых оболочек, так называемые а-формы. Для двухмерного пространства а-формы обусловливаются окружностями Делоне, полученными из триангуляции Делоне, а в трехмерном — сферами Делоне. По своей сути а-формы представляют целое семейство форм в диапазоне от точки до выпуклой оболочки; а-форма параметризуется действительным числом а £ [0, то). При а = то а-форма — выпуклая оболочка множества Р, при а = 0 а-форма состоит только из множества точек Р. С уменьшением значения а от бесконечности до нуля а-форма уменьшается за счет образования полостей внутри рассматриваемого множества.

Для 0 < а < то а-кругом (а-Ьа11) называется открытый (геометрическое множество точек, не включающее границу) круг радиуса а. а-Круг Ь называется пустым, если Ь П Р = 0. Для любого подмножества точек Т С Р в пространстве Д^, где |Т| = к +1 < d +1, выпуклая оболочка Ду = еопуТ называется к-симплексом. Для 0 < к < 2 к-симплекс называется а-открытым (а-ехроэе^, если существует пустой а-шар Ь, где Т = дЬ П Р. На рис. 7, а изображен а-открытый симплекс, на рис. 7, б симплекс не является а-открытым, так как окружность не является пустой.

Фиксированным значением параметра а определяется множество Рк,а а-открытых к-симплексов, где 0 < к < 2. Тогда определение а-формы вводится следующим образом: а-формой множества Р являются многоугольники из множества симплексов Л2 а, ребра из и вершины из Р0,а [25]. По определению, граница дРа а-формы множества

а

Рис. 7. Симплекс а-открытый (а); симплекс, не являющийся а-открытым (б)

точек P состоит из всех а-открытых k-симплексов. Чтобы найти границу области, необходимо построить все возможные k-симплексы и определить среди них а-открытые. На практике для вычисления границы используется триангуляция Делоне. Связь между а-формой и триангуляцией Делоне определяется следующим образом: если k-симплекс Ду множества P является а-открытым, то A у принадлежит триангуляции Делоне. Соответственно, граница множества P принадлежит множеству триангуляции Делоне. Для определения симплексов, принадлежащих а-форме, используется понятие а-теста: описанная вокруг симплекса Ду окружность радиуса оу с центром в точке ^у должна быть пустой и оу < а.

В данной работе для нахождения границ области используется алгоритм, предложенный Эдельсбруннером [25]. В соответствии с этим алгоритмом, для заданного набора точек P и 0 < а < ж необходимо вычислить а-комплекс Ca (P) — подмножество триангуляции Делоне, обладающее одним из следующих свойств: k-симплекс Ду удовлетворяет условию а-теста или k-симплекс является частью поверхности (k + 1)-симплекса большей размерности из Ca (S). Тогда для нахождения границ области по заданному набору точек необходимо произвести:

— построение триангуляции Делоне;

— проверку каждого треугольника на принадлежность Ca (S).

Чтобы определить, является ли ребро треугольника граничным или внутренним, необходимо проверить треугольники, содержащие данное ребро, на принадлежность а-комплексу Ca (P). Если один из двух треугольников принадлежит а-комплексу, а второй не принадлежит, или существует только один треугольник из а-комплекса, то ребро является граничным. Если оба ребра из Ca (P), то ребро внутреннее. Во всех остальных случаях считается, что симплекс не принадлежит области.

5. Результаты

Алгоритм метода "sweep line" реализован на языке Fortran с применением динамических структур типа списков и деревьев. Тестирование метода проводилось на областях с различной геометрией по произвольно или специально заданному набору точек на плоскости. На рис. 8, a показано разбиение области ячейками Вороного, в которой положение частиц получено из решения нестационарной задачи о движении солитона по дну с расположенным на нем уступом в момент опрокидывания гребня волны [1]. Как видно из рисунка, граничные ячейки имеют бесконечные ребра на плоскости. Для решения данной проблемы область заключается в прямоугольник большего размера, и ребра отсекаются по границе этого прямоугольника.

На рис. 9 дан график зависимости времени построения диаграммы Вороного от количества точек. Тестирование проводилось на персональном компьютере Intel Pentium 3Ghz, 512Mb RAM. Пунктирная кривая соответствует серии численных расчетов, а сплошная кривая — график функции вида y = an1'4 log2 (n), где константа a =1/ (21 •lO6). Из хорошего качественного совпадения двух кривых можно сделать вывод, что число операций в реализованном алгоритме составляет O (n1'4 log(n)), что несколько больше числа операций в теоретической оценке вида O (n log (n)). Данное несовпадение можно объяснить двумя причинами. Во-первых, в теоретической оценке не учитывается число операций, необходимых для отсечения ребер бесконечной длины, во-вторых, в теоретическую оценку не входят специфические операции, подготавливающие данные для дальнейшего их использования в бессеточных методах NEM и MFEM.

а

б

в

г

Рис. 8. Диаграмма Вороного (а); первоначальная триангуляция Делоне (б); триангуляция Делоне при a = 0.25 (в); триангуляция Делоне при a = 0.0596 (г)

Алгоритм метода a-shape реализован на языке Fortran. В табл. 4 представлено время определения границы в зависимости от количества точек в области.

Оценку числа операций для алгоритма a-shape дать сложно, все зависит от расположения точек на плоскости. В общем случае отстоящая на расстояние от "основного"

Таблица 4. Время определения границы области

Количество точек 503 1986 4319 6048 8180 10000 12624

Время, с 9.4в-4 3.74в-3 8.14в-3 1.15в-2 1.53в-2 1.84в-2 2.34в-2

Рис. 9. График зависимости времени построения диаграммы Вороного от количества точек в области

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

На рис. 8, б—г показан процесс определения границы области по заданному набору точек. На первом шаге из диаграммы Вороного (рис. 8, а) восстанавливается первоначальная триангуляция Делоне (рис. 8, б). Затем область очищается от элементов, не удовлетворяющих критерию а-теста. Для определения геометрии границы заданной области необходимо подобрать значение параметра а, дающего "приемлемые" результаты. Универсальных методов выбора значения а нет, существуют лишь некоторые рекомендации, полученные из опыта применения алгоритма и основанные на соотношениях расстояний между самыми близкими и самыми удаленными точками области. На рис. 8, в приведена область после работы алгоритма а-эЬаре со значением а = 0.25, на рис. 8, г а = 0.0596. Как видно из последнего рисунка, граница области восстановилась так, что гребень волны соприкасается с ее подошвой, чего не было в распределении точек, взятом из решения нестационарной задачи о движении солитона. Этот нежелательный эффект возникает всегда, если расстояние между двумя точками области меньше, чем значение параметра а. Существуют некоторые модификации алгоритма а-эЬаре, основанные на плотности распределения точек в области [26]. Реализация и тестирование данных алгоритмов показывают, что они способны справиться с подобными трудностями, но их область применения ограничивается вычислительной геометрией и распознаванием образов. Для решения задач механики жидкости они не подходят.

Заключение

Развитие современных бессеточных численных методов для моделирования течений вязкой и идеальной несжимаемых жидкостей является привлекательной задачей: такие методы способны проводить расчеты на многосвязных областях и областях с большими деформациями. Несмотря на отсутствие расчетной сетки, в бессеточных методах для построения интерполирующих функций необходимо обеспечить узловую связность, с помощью которой определяется область влияния расчетных узлов друг на друга. Как правило, связность для каждого метода строится по определенному критерию, например, в методе SPH определяется, какие узлы попадают в описанную окружность радиуса r с центром в заданном узле. Если не применяются специальные алгоритмы, то методы построения узловой связности неэффективны и занимают достаточно много времени. С увеличением числа расчетных узлов рост временных затрат происходит по нелинейному закону. Помимо этого, такие задачи, как определение динамической нагрузки, создаваемой жидкостью на твердые стенки, требует дополнительно знания упорядоченной последовательности расчетных узлов для вычисления интегралов по границам. Предложенный в данной работе подход на основе построения диаграммы Вороного позволяет преодолевать все указанные трудности для большинства бессеточных методов. На основании структуры данных, полученной при построении диаграммы Вороного, по заданному множеству узлов строится триангуляция Делоне, определяются границы расчетной области и естественные соседи, выстраивается упорядоченный массив граничных узлов. Использование диаграммы Вороного позволяет существенно сократить область поиска расчетных узлов для большинства бессеточных методов, таких как SPH, MPS, MFEM, NEM.

Реализованные алгоритмы построения диаграммы Вороного, триангуляции Делоне и определения границы области включены в комплекс программ для расчета течений идеальной и вязкой несжимаемой жидкости бессеточным методом конечных элементов и методом естественных соседей. Апробация описанного в данной статье подхода проводилась в работах [17, 27, 28]. Применение в методах NEM и MFEM алгоритма "sweep line" для определения структуры данных и a-shape для построения границы области позволило значительно сократить временные затраты поиска естественных соседей и определения положения свободной границы области на каждом временном шаге.

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

[1] Afanasiev K.E., Stukolov S.V. Simulation of problems with free surface by a boundary element method // Comput. Technology. 2003. Vol. 8. Part 3. Special Issue. P. 3-33.

[2] Бреббиа К., Телес Ж., Вроувел Л. Методы граничных элементов: пер. с англ. М.: Мир, 1987. 524 с.

[3] Зенкевич О., Морган К. Конечные элементы и аппроксимация. М.: Мир, 1986. 317 с.

[4] Ильин В.П. Методы конечных разностей и конечных объемов для эллиптических уравнений. Новосибирск: Ин-т математики СО РАН, 2000. 345 с.

[5] Скворцов А.В. Триангуляция Делоне и ее применение. Томск: Изд-во Томск. у-та, 2002. 128 с.

[6] Численное моделирование течений жидкости с поверхностными волнами / Г.С. Хаким-зянов, Ю.И. Шокин, В.Б. Барахнин, Н.Ю. Шокина. Новосибирск: Изд-во СО РАН, 2001. 394 с.

[7] ЛисЕйкин В.Д., Молородов Ю.И., ХАкимзянов Г.С. Об интерактивном комплексе программ построения двумерных структурных сеток // Вычисл. технологии. 2000. Т. 5, № 1. С. 70-84.

[8] Методы римановой геометрии в задачах построения разностных сеток / Ю.И. Шокин,

B.Д. Лисейкин, А.С. Лебедев, Н.Т. Данаев, И.А. Китаева. Новосибирск: Наука, 2005. 256 с.

[9] Li S., Liu W.K. Meshfree and particle methods and their applications // Appl. Mech. Rev. 2002. Vol. 55. P. 1-34.

[10] Афанасьев К.Е. Численное моделирование задач гидродинамики со свободными границами современными численными методами // Матер. III Междунар. научной летней школы "Гидродинамика больших скоростей и численное моделирование". Кемерово: ИНТ, 2006. P. 53-65.

[11] MoNAGHAN J. Smoothed particle hydrodynamics // Ann. Rev. Astron and Astrophysics. 1992. Vol. 30. P. 543-574.

[12] Koshizuka S., Tamako H., Oka Y. A particle method for incompressible viscous flow with fluid fragmentation // Comput. Fluid Dynamics J. 1995. Vol. 4. P. 29-46.

[13] Афанасьев К.Е., Макарчук Р.С., Ильясов А.Е., Попов А.Ю. Численное моделирование течений жидкости со свободными границами методами SPH и MPS // Вычисл. технологии. 2006. Т. 11. Спецвып.: Избранные докл. семинара по числ. методам и информ. технологиям Кемеровского гос. ун-та. С. 26-45.

[14] Франк А.М. Дискретные модели несжимаемой жидкости. М.: Физматлит, 2001. 208 с.

[15] Facundo P. The meshless finite element method applied to a lagrangian particle formulation of fluid flows: In Partial Fulfillment of the Requirements for the Degree Doctor of Philosophy. Instituto de Desarrollo tecnologico para la industria quimica (INTEC) universidad nacional del litoral noviembre. 2003. 157 p.

[16] Sukumar N., Moran B., Belytschko T. The natural element method in solid mechanics // Int. J. Num. Methods Eng. 1998. Vol. 43, N 5. P. 839-887.

[17] Афанасьев К.Е., Каравцев С.Н., Рейн Т.С., Стуколов С.В. Метод естественных соседей на основе интерполяции Сибсона // Вестн. Томск. у-та. 2006. № 19. С. 210-219.

[18] Sibson R. A vector identity for the Dirichlet Tessellation // Math. Proc. Cambridge Phil. Soc. 1980. Vol. 87, N 1. P. 151-155.

[19] Беликов В.В., Иванов В.Д., Конторович В.К., Корытник С.А., Семенов А.Ю. Несибсоновская интерполяция — новый метод интерполяции значений функции на произвольной системе точек // Журн. вычисл. математики и мат. физики. 1997. Т. 37, № 1.

C. 11-17.

[20] Belikov V.V., Semenov A.Y. Non-sibsonian interpolation on arbitrary system of points in euclidean space and adaptative generating isolines algorithm // Appl. Numer. Mathematics. 2000. Vol. 32. P. 371-387.

[21] Препарата Ф., Шеймос М. Вычислительная геометрия: введение. М.: Мир, 1989. 450 с.

[22] De Berg M., Van Kreveld M. Computational Geometry. Algorithms and Applications. Second, Revised Edition. Berlin: Springer-Verlag, 2000.

[23] Green P., Sibson R. Computing dirichlet tessellations in the plane // Comp. J. 1977. Vol. 21, N 22. P. 168-173.

[24] Fortune S.J. A sweepline algorithm for Voronoi diagrams // J. Algorithmica. 1987. N 2. P. 153-174.

[25] Edelsbrunner H., Maoke E.P. Three-dimensional alpha shapes. // ACM Trans. Graph. 1994. Vol. 13, N 1. P. 43-72.

[26] Teiohmann M., Capps M. Surface reconstruction with anisotropic density-scaled alpha shapes // Proceedings of the IEEE Visualization Conference. 1998.

[27] КАРАБЦЕВ С.Н., Рейн Т.С. Применение метода естественных соседей к решению задач механики жидкости со свободными поверхностями // Матер. III Междунар. научной летней школы "Гидродинамика больших скоростей и численное моделирование". Кемерово: ИНТ, 2006. С. 393-401.

[28] КАРАБЦЕВ С.Н., Стуколов С.В. Эффективный алгоритм генерации конечноэлемент-ной сетки для метода естественных соседей // Матер. III Междунар. научной летней школы "Гидродинамика больших скоростей и численное моделирование". Кемерово: ИНТ, 2006. С. 401-409.

Поступила в редакцию 15 июня 2007 г., в переработанном виде — 5 марта 2008 г.