Журнал вычислительной математики и математической физики, 2022, T. 62, № 8, стр. 1237-1250

Построение несимплициальных сеток делоне посредством аппроксимации радикальными разбиениями

В. А. Гаранжа 12*, Л. Каменски 2***, Л. Н. Кудрявцева 12**

1 ВЦ им. А.А. Дородницына ФИЦ ИУ РАН
119333 Москва, ул. Вавилова, 40, Россия

2 МФТИ
141701 М.о., Долгопрудный, Институтский пер., 9, Россия

* E-mail: garan@ccas.ru
*** E-mail: l.kamenski@arcor.de
** E-mail: liukudr@yandex.ru

Поступила в редакцию 01.12.2021
После доработки 31.12.2021
Принята к публикации 10.01.2022

Полный текст (PDF)

Аннотация

Рассматривается построение полиэдрального разбиения Делоне как предела последовательности степенных диаграмм (так называемых радикальных разбиений), а двойственная диаграмма Вороного получается как предел последовательности взвешенных разбиений Делоне. С использованием известной концепции подъема точек на параболоид задача сводится к построению пары двойственных выпуклых многогранников как предела последовательности пар общих двойственных выпуклых многогранников, вписанных и описанных вокруг кругового параболоида. При этом последовательность первичных многогранников должна сходиться к описанному многограннику, а последовательность двойственных многогранников сходится к вписанному многограннику. Для построения последовательностей пар взаимно двойственных многогранников нас интересует случай, когда вершины первичных многогранников могут перемещаться или сливаться. Это значит, что для двойственных многогранников не допускается появление новых граней. Данные правила по сути определяют преобразование множества заданных шаров, определяющих степенную диаграмму, во множество шаров Делоне, используя перемещение шара, изменение радиуса и удаление шара как допустимые операции. Хотя строгое обоснование (теоремы существования) для этой задачи пока недоступно, мы предлагаем функционал, измеряющий отклонение общего выпуклого многогранника от многогранника, вписанного в параболоид. Этот функционал является дискретным функционалом Дирихле для функции степени (так называемой степени точки относительно сферы), которая является линейным интерполянтом расстояния двойственных вершин от параболоида по вертикали. Абсолютный минимимум этого функционала достигается в случае, когда степень постоянна, т.е. вписанный многогранник может быть получен из минимизирующего многогранника с помощью параллельного переноса. Функционал Дирихле для двойственной поверхности не является квадратичной, поскольку неизвестными являются вершины первичного многогранника. Следовательно, преобразование множества шаров в шары Делоне в общем случае не является единственным. В данной работе мы сосредоточились на экспериментальном подтверждении работоспособности предложенного подхода и оставили в стороне проблемы качества получающейся сетки. Нулевое значение градиента предложенного функционала определяет многообразие, описывающее эволюцию сфер Делоне. Следовательно, сетки Делоне–Вороного могут быть оптимизированы с использованием этого многообразия в качестве ограничения. Численные примеры иллюстрируют построение многоугольных сеток Делоне в плоских областях. Библ. 19. Фиг. 15.

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

1. ВВЕДЕНИЕ

В своем знаменитом докладе “Sur la sphere vide” на конгрессе по геометрии в Торонто в 1924 г., а затем в своих работах [1], [2], Борис Николаевич Делоне ввел понятие нормального разбиения пространства, заданного дискретным множеством вершин, состоящего из $d$-мерных симплексов с пустыми описанными шарами, пустыми в том смысле, что они не содержат внутри себя никаких вершин. Знаменитая лемма Делоне гласит, что если условие пустого шара выполняется локально для любых двух соседних симплексов, имеющих общую $(d - 1)$-мерную грань, то пустыми являются все описанные шары. Для доказательства этого результата Делоне использовал понятие степени $\tau $ точки ${\mathbf{a}}$ относительно шара $B$ радиуса $R$ с центром ${\mathbf{c}}$:

${{\tau }_{B}}({\mathbf{a}}) = {{\left| {{\mathbf{c}} - {\mathbf{a}}} \right|}^{2}} - {{R}^{2}}.$

В 1937 г. Делоне представил L-разбиение (см. [3], [4]), которое получается путем перемещения, сжатия и раздувания пустого шара, расположенного среди вершин дискретного множества в ${{\mathbb{R}}^{d}}$ (фиг. 1а). Таким образом, определяются все возможные пустые шары с $d$-мерным набором вершин на их поверхности, так называемые L-шары.

Фиг. 1.

L-разбиение: (а) – пустой шар, движущийся через множество точек, (б) – пустые шары с $d$-мерными множествами точек на границе.

Выпуклая оболочка множества вершин, лежащих на L-сфере, представляет собой выпуклый многогранник, названный Делоне L-полиэдром (фиг. 1б). Множество L-полиэдров определяет нормальное разбиение пространства, называемое L-разбиением. Лемма Делоне естественным образом обобщается на это многогранное разбиение. В наши дни L-разбиение называют разбиением Делоне, а L-шары, т.е. пустые шары с $d$-мерным набором вершин на их границе, называют шарами Делоне. Выпуклая оболочка центров сфер Делоне, проходящих через одну из точек дискретного множества, определяет многогранник Вороного для этой точки. Заметим, что как разбиения Делоне, так и разбиения Вороного для заданного набора вершин единственны. Отношения двойственности между разбиениями Делоне и Вороного просты и элегантны, чего нельзя сказать о триангуляциях Делоне. Например, для набора вершин, лежащих на одной сфере, один и тот же центр окружности (вершина Вороного) может быть порожден несколькими симплексами Делоне, тогда как соответствие между многогранником Делоне и соответствующей вершиной Вороного всегда взаимно однозначно.

Каждый многогранник Делоне может быть разбит на симплексы, тем самым задавая нормальное разбиение, которое обычно называется триангуляцией Делоне, при условии, что триангуляции отдельных многогранников Делоне согласованы между собой. Стоит обратить внимание на существенное различие: в разбиении Делоне имеется хотя бы один замкнутый пустой шар, в то время как ребро/грань в триангуляции Делоне имеет хотя бы один открытый пустой шар. Ребра/грани триангуляции, которые не имеют замкнутых пустых шаров – это ребра/грани, добавленные при триангуляции многогранников Делоне. Эти дополнительные ребра не имеют взаимно однозначного соответствия с двойственными гранями Вороного, а дополнительные грани соответственно не имеют взаимно однозначного соответствия с двойственными ребрами Вороного.

Заметим, что построение, основанное на правильной триангуляции разбиения Делоне, исключает плоские вырожденные симплексы, так называемые конверты (slivers). Строго говоря, конверт – это вырожденный симплекс с $(d - m)$-мерным множеством вершин ($m > 0$), для которого можно построить описанную сферу конечного радиуса. Очевидно, что такая сфера не является единственной. Каждый такой симплекс порождается неправильной триангуляцией $(d - m)$-мерной грани разбиения Делоне. Проблема идентификации конвертов усугубляется, если речь идет о приближенных вычислениях, что, например, имеет место в стандартной арифметике с плавающей точкой.

При численном моделировании конверты недопустимы, поскольку они могут ухудшить точность аппроксимации конечных элементов и конечных объемов. Стандартное современное решение для устранения конвертов заключается в локальном нарушении условия пустой сферы Делоне. Если условие Делоне должно быть выполнено, что особенно актуально при использовании методов конечных объемов с использованием ячеек Вороного, выполняющих принцип максимума (см. [5], [6]), возникает печально известная и неприятная особенность современных алгоритмов построения триангуляций Делоне: порождение искусственных конвертов. Например, триангуляция вершин кубической решетки может привести к появлению значительного количества конвертов (фиг. 2), количество и расположение которых непредсказуемо зависит от ошибок округления входных данных, создавая хотя и весьма живописные, но крайне нежелательные узоры даже для самых простых входных данных (фиг. 3). Получение триангуляции Делоне с помощью радикального разбиения сложнее, но позволяет избавиться от искусственных конвертов, не нарушая условия Делоне.

Фиг. 2.

Конверты (slivers) в сетках Делоне генератора TetGen, вызванные округлением во входных данных для простой кубической сетки (а) и более сложного примера неравномерной декартовой сети вершин (б).

Фиг. 3.

“Sliver Art № 3”: с правильно подобранной цветовой картой визуализация выборки конвертов (а) напоминает известную серию Казимира Малевича (б). (Казимир Малевич, “Динамический супрематизм № 38”, общественное достояние, Wikimedia Commons.)

Заметим, однако, что то же обозначение “конверт” иногда используется для симплекса плохой формы, который далек от плоского, но вычисление его центра и радиуса окружности устойчиво. Такие конверты являются законными симплексами Делоне плохой формы и не должны исключаться из триангуляции. Разница между корректным, но почти плоским симплексом, и неправильной триангуляцией грани многогранника Делоне должна быть четко определена. Естественную идентификацию дает степенная диаграмма (радикальное разбиение) для множества шаров, которая восходит к работам Рене Декарта. Радикальное разбиение определяется как множество выпуклых многогранников, построенных пересечениями полупространств, определяемых $(d - 1)$-мерными плоскостями равной степени, ортогональными к отрезкам, соединяющим центры всех пар шаров (подробнее см. разд. 2). Для двух пересекающихся сфер радикальная плоскость всегда проходит через их общее множество. Следовательно, для множества шаров Делоне радикальное разбиение совпадает с разбиением Делоне.

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

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

Фиг. 4.

(а) – Конверт, возникший из-за ошибок округления арифметики с плавающей точкой, создает лишнюю вершину радикального разбиения. (б) – Корректная радикальная ячейка, совпадающая с ячейкой Делоне.

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

Триангуляции Делоне имеют множество применений в самых разных областях (см. [8]). Одним из очень важных свойств для численного моделирования является выполнение принципа максимума для дискретного оператора Лапласа на сетках Делоне (см. [5], [6]). Наиболее значимой для построения расчетных сеток Вороного без их обрезания около границ расчетной области является возможность контролировать расположение шаров Делоне в ключевых областях вычислительной области, в частности, на границах и вблизи них (см. [9]–[11]). Мы предлагаем вычислительную схему, которая потенциально может служить инструментом управления размещения шаров Делоне и обеспечивать построение сеток Делоне–Вороного в обратную сторону, когда набор шаров генерирует набор вершин для построения сеток Делоне.

2. СТЕПЕННАЯ ДИАГРАММА И ПОДЪЕМ ТОЧЕК НА ПАРАБОЛОИД

Идея подъема точек на параболоид берет начало в работах Г.Ф. Вороного (см. [12]), который показал, что триангуляция Делоне в ${{\mathbb{R}}^{d}}$ является проекцией граней выпуклого многогранника $P \in {{\mathbb{R}}^{{d + 1}}}$, вписанного в круговой параболоид $\Pi $. Выпуклое тело $P{\kern 1pt} {\text{*}}$, построенное как пересечение верхних полупространств над касательными плоскостями к $\Pi $ в вершинах $P$, называется женератрисой Вороного. Проекция его граней на ${{\mathbb{R}}^{d}}$ задает диаграмму Вороного. Более общая концепция подъема (см. [13], Edelsbrunner-2) основана на построении пары $P$, $P{\kern 1pt} {\text{*}}$ выпуклых многогранников, удовлетворяющих отношению полярности (см. [14]) относительно параболоида $\Pi $.

Рассмотрим систему шаров $\mathcal{B} = {\kern 1pt} {{B}_{1}},\; \ldots ,\;{{B}_{n}}$, определенных центрами ${{{\mathbf{c}}}_{i}}\;{{\mathbb{R}}^{d}}$ и радиусами ${{R}_{i}} \geqslant 0$, поднятую систему точек ${{\mathcal{E}}_{l}} = \left\{ {{{{\mathbf{p}}}_{1}},\; \ldots ,\;{{{\mathbf{p}}}_{n}}} \right\}$ в ${{\mathbb{R}}^{{d + 1}}}$, где

${\mathbf{p}}_{i}^{{\text{т}}} = \left( {{\mathbf{c}}_{i}^{{\text{т}}},\frac{1}{2}\left( {{{{\left| {{{{\mathbf{c}}}_{i}}} \right|}}^{2}} - R_{i}^{2}} \right)} \right) = \left( {{\mathbf{c}}_{i}^{{\text{т}}},{{h}_{i}}} \right),$
и нижнюю выпуклую оболочку ${{\mathcal{E}}_{l}}$, определяемую выпуклой функцией ${{x}_{{d + 1}}} = {v}({{x}_{1}},\; \ldots ,\;{{x}_{d}})$. Функция, двойственная по Лежандру–Юнгу–Фенхелю (см. [15]) к функции ${v}$, обозначается через ${v}{\kern 1pt} {\text{*}}$. Если быть математически точным, то функцию ${v}({\mathbf{x}})$ следует доопределять $ + \infty $ вне выпуклой оболочки ${\text{conv}}({{{\mathbf{c}}}_{1}},\; \ldots ,\;{{{\mathbf{c}}}_{n}})$ так, чтобы ее эпиграф (надграфик) представлял собой замкнутое множество. Двойственная функция ${v}{\kern 1pt} {\text{*}}({\mathbf{x}})$ определена везде, и ее график содержит неограниченные грани. В дальнейшем мы будем рассматривать ограниченные задачи, в которых неограниченные грани по существу исключены из постановки задачи.

Проекция граней графика ${v}$ определяет взвешенную триангуляцию Делоне $\mathcal{W}$ в ${{\mathbb{R}}^{d}}$ (см. [16]). Вершинами графика ${v}$ являются пары $({{{\mathbf{c}}}_{i}},{{h}_{i}})$, а вершинами взвешенной триангуляции Делоне – ${{{\mathbf{c}}}_{i}}$. ${{T}_{k}}$ обозначает $k$-й взвешенный многогранник Делоне.

Проекция граней графика ${v}{\kern 1pt} {\text{*}}$ определяет радикальное разбиение (степенную диаграмму) $\mathcal{R}$ для системы шаров в ${{\mathbb{R}}^{d}}$ (см. [16]) (фиг. 5). Проекция $k$-й вершины графа ${v}{\kern 1pt} *$ на ${{x}_{{d + 1}}} = 0$ обозначается через ${{{\mathbf{v}}}_{k}}$. Эта точка двойственна к ${{T}_{k}}\mathcal{W}$, а вершина ${{{\mathbf{c}}}_{i}}$ двойственна ячейке ${{D}_{i}}$ из $\mathcal{R}$.

Фиг. 5.

(а) – Взвешенная триангуляция Делоне. (б) – Степенная диаграмма и двойственные многогранники при подъеме на параболоид. Рисунок построен с помощью detri2 Ханга Си.

3. ДВИЖЕНИЕ ШАРОВ КАК ПРЕОБРАЗОВАНИЕ ПАР ДВОЙСТВЕННЫХ МНОГОГРАННИКОВ

Рассмотрим следующую задачу: переместить и масштабировать шары $\mathcal{B}$ так, чтобы все вершины графика ${v}{\kern 1pt} {\text{*}}$ сошлись к поверхности параболоида ${{x}_{{d + 1}}} = \Pi ({\mathbf{x}}) = {{\left( {x_{1}^{2} + \; \cdots \; + x_{d}^{2}} \right)} \mathord{\left/ {\vphantom {{\left( {x_{1}^{2} + \; \cdots \; + x_{d}^{2}} \right)} 2}} \right. \kern-0em} 2}$. Это означает, что проекция графика ${v}{\kern 1pt} {\text{*}}$ в конечном итоге сходится к разбиению Делоне. Заметим, что количество вершин ${{{v}}_{k}}$ может меняться во время движения шаров. Далее удобно использовать обозначение $\left( {{{{\mathbf{x}}}^{{\text{т}}}},{{x}_{{d + 1}}}} \right)$ для произвольной точки в ${{\mathbb{R}}^{{d + 1}}}$, где ${{{\mathbf{x}}}^{{\text{т}}}} = \left\{ {{{x}_{1}},\; \ldots ,\;{{x}_{d}}} \right\}$.

Для множества шаров $\mathcal{B}$ мы строим первичный и двойственный многогранники $P$ и $P$, определяемые выпуклыми кусочно-линейными функциями ${v}({\mathbf{x}})$ и ${v}{\kern 1pt} {\text{*}}({\mathbf{x}})$ соответственно. Согласно соотношению полярности относительно кругового параболоида [14], вершина $\left( {{{{\mathbf{c}}}_{i}},{{h}_{i}}} \right)$ определяет двойственную ей плоскость грани графика ${v}{\kern 1pt} {\text{*}}({\mathbf{x}})$:

${\mathbf{c}}_{i}^{{\text{т}}}x = {{h}_{i}} + {{x}_{{d + 1}}}.$
Вершина $({{{\mathbf{v}}}_{k}},{{z}_{k}})$ графика ${v}{\kern 1pt} {\text{*}}({\mathbf{x}})$ является пересечением по крайней мере $d + 1$ таких плоскостей:
(1)
${\mathbf{v}}_{k}^{{\text{т}}}({{{\mathbf{c}}}_{i}} - {{{\mathbf{c}}}_{p}}) = {{h}_{i}} - {{h}_{p}},$
где $i \ne p$ – индексы всех плоскостей, пересекающихся в $({{{\mathbf{v}}}_{k}},{{z}_{k}})$. Следовательно,
${{z}_{k}} = {\mathbf{v}}_{k}^{{\text{т}}}{{{\mathbf{c}}}_{i}} - {{h}_{i}} = {\mathbf{v}}_{k}^{{\text{т}}}{{{\mathbf{c}}}_{i}} - \frac{1}{2}{\mathbf{c}}_{i}^{2} + \frac{1}{2}R_{i}^{2} = - \frac{1}{2}\left( {{{{\left| {{{{\mathbf{c}}}_{i}} - {{{\mathbf{v}}}_{k}}} \right|}}^{2}} - R_{i}^{2}} \right) + \frac{1}{2}{\mathbf{v}}_{k}^{2}.$
Другими словами,
(2)
${{z}_{k}} - \Pi ({{{\mathbf{v}}}_{k}}) = - \frac{1}{2}{{\tau }_{i}}({{{\mathbf{v}}}_{k}}),$
где
${{\tau }_{i}}({\mathbf{y}}) = {{\left| {{{{\mathbf{c}}}_{i}} - {\mathbf{y}}} \right|}^{2}} - R_{i}^{2}$
является степенью точки ${\mathbf{y}}$ относительно шара ${{B}_{i}}$. Следовательно, вертикальное расстояние вершины $({{{\mathbf{v}}}_{k}},{{z}_{k}})$ от параболоида $\Pi $ полностью определяется величиной степени. Другая интерпретация(1) состоит в том, что для вершины ${{{\mathbf{v}}}_{k}}$, дуальной к взвешенному многограннику Делоне ${{T}_{k}}$, равенство
${{\tau }_{i}}({{{\mathbf{v}}}_{k}}) = {{\tau }_{p}}({{{\mathbf{v}}}_{k}})$
выполняется для всех вершин ${{T}_{k}}$. Это означает, что в нашей ситуации можно опустить индексы и использовать обозначение $\tau ({{{\mathbf{v}}}_{k}})$ для значения степени.

Соотношение (1) подразумевает, что градиент ${v}{\kern 1pt} {\text{*}}(x)$ на $i$-й грани своего графика равен ${{{\mathbf{c}}}_{i}}$. Верно и двойственное утверждение: градиент функции ${v}(x)$ на $k$-й грани графика равен ${{v}_{k}}$. Для любого выпуклого многогранника ${{T}_{k}}$ существует $d$ линейных независимых векторов ${{{\mathbf{c}}}_{i}} - {{{\mathbf{c}}}_{p}}$. В двух измерениях ${{T}_{k}}$ в простейшем случае – треугольник, и это множество просто равно $\left\{ {{{{\mathbf{c}}}_{2}} - {{{\mathbf{c}}}_{1}},\;{{{\mathbf{c}}}_{3}} - {{{\mathbf{c}}}_{1}}} \right\}$ (фиг. 6а). Когда $\tau ({{v}_{k}}) > 0$, можно связать с ${{v}_{k}}$ сферу с радиусом $\sqrt {\tau ({{v}_{k}})} $, которая называется ортосферой. В наших численных экспериментах в разд. 5 мы рисуем искусственные ортосферы с радиусами, определяемыми формулами $\sqrt {\left| {\tau ({{v}_{k}})} \right|} $. Как показано ниже, эти сферы визуализируют отклонение радикального разбиения от разбиения Делоне.

Фиг. 6.

(а) – Взвешенный треугольник Делоне ${{T}_{k}}$, двойственная вершина ${{{\mathbf{v}}}_{k}}$ и ортоокружность. (б) – Радикальная ячейка разбивается на четыре треугольника Делоне; центр ${{{\mathbf{c}}}_{i}}$ аппроксимируется облаком из центров четырех кругов Делоне.

Пусть ${\tilde {v}}{\kern 1pt} {\text{*}}$ – спроецированный вариант функции ${v}{\kern 1pt} {\text{*}}$, построенный следующим образом: рассмотрим множество вершин $({{{\mathbf{v}}}_{j}},{{z}_{j}})$ графа ${v}{\kern 1pt} {\text{*}}$ и спроецируем их на параболоид $\Pi $, задавая

${{\tilde {z}}_{j}} = \frac{1}{2}{{\left| {{{{\mathbf{v}}}_{j}}} \right|}^{2}}\quad {\text{и}}\quad {{{\mathbf{\tilde {v}}}}_{j}} = {{{\mathbf{v}}}_{j}}.$
Как показано выше, эта проекция изменяет вертикальную составляющую ${{z}_{j}}$ на $\tau ({{{\mathbf{v}}}_{j}}){\text{/}}2$.

Вычисление нижней выпуклой оболочки системы точек $\left\{ {({{{{\mathbf{\tilde {v}}}}}_{k}},{{{\tilde {z}}}_{k}})} \right\}$, $k = 1,\;2,\; \ldots ,\;{{n}_{{v}}}$, приводит к графику функции ${\tilde {v}}{\kern 1pt} {\text{*}}$.

4. ФУНКЦИОНАЛ ДИРИХЛЕ ДЛЯ ФУНКЦИИ СТЕПЕНИ

Для оценки меры близости текущего радикального разбиения и разбиения Делоне мы рассматриваем функционал Дирихле для разности $v{\text{*}}$ и ${\tilde {v}}{\kern 1pt} {\text{*}}$:

(3)
$F(X) = \frac{1}{2}\int\limits_\Omega {{{{\left| {\nabla {v}{\kern 1pt} {\text{*}} - \nabla {\tilde {v}}{\kern 1pt} *} \right|}}^{2}}} d{\mathbf{x}}.$
Здесь $X$ – вектор неизвестных, состоящий из ${{{\mathbf{c}}}_{i}}$ и ${{R}_{i}}$, а $\Omega $ – ограниченная область определения функции ${\tilde {v}}{\kern 1pt} *$. Из соотношения (2) следует
$F(X) = \frac{1}{8}\int\limits_\Omega {{{{\left| {\nabla \tau ({\mathbf{x}})} \right|}}^{2}}} d{\mathbf{x}}.$
Пусть $\tau ({\mathbf{x}})$ обозначает кусочно-линейную функцию, которая совпадает с $\tau ({{{\mathbf{v}}}_{k}})$ при ${{{\mathbf{v}}}_{k}}$ и является линейной на каждой ячейке ${{\tilde {T}}_{j}}$ вспомогательного разбиения Делоне. Предполагая, что каждая грань графа ${v}{\kern 1pt} {\text{*}}$ проецируется на параболоид независимо от других граней, $F(X)$ можно переписать как
(4)
${{F}_{I}}(X) = \frac{1}{2}\sum\limits_i {\sum\limits_{{{{\tilde {T}}}_{j}} \in {{D}_{i}}} {{{{\left| {{{{\mathbf{c}}}_{i}} - {{{\mathbf{s}}}_{j}}} \right|}}^{2}}} } \operatorname{vol} {{\tilde {T}}_{j}},$
где ${{{\mathbf{s}}}_{j}}$ – центр описанной сферы симплекса Делоне ${{\tilde {T}}_{j}}$. Приведенное выше равенство является очевидным следствием того, что
${{\left. {\nabla {v}{\kern 1pt} {\text{*}}} \right|}_{{{{D}_{i}}}}} = {{{\mathbf{c}}}_{i}},\quad {{\left. {\nabla {\tilde {v}}{\kern 1pt} {\text{*}}} \right|}_{{{{{\tilde {T}}}_{j}} \in {{D}_{i}}}}} = {{{\mathbf{s}}}_{j}}.$
Чтобы пояснить эту формулу, рассмотрим симплекс Делоне ${{\tilde {T}}_{j}}$ с вершинами ${{{\mathbf{v}}}_{1}},\; \ldots ,\;{{{\mathbf{v}}}_{{d + 1}}}$. Градиент ${{{\mathbf{g}}}_{j}}$ от ${\tilde {v}}{\kern 1pt} {\text{*}}$ определяется следующим образом:
${{({{{\mathbf{v}}}_{m}} - {{{\mathbf{v}}}_{l}})}^{{\text{т}}}}{{g}_{j}} = \frac{1}{2}\left( {{\mathbf{v}}_{m}^{2} - {\mathbf{v}}_{l}^{2}} \right),\quad m,l \leqslant d + 1,$
или
${{\left| {{{g}_{j}} - {{{\mathbf{v}}}_{m}}} \right|}^{2}} = {{\left| {{{g}_{j}} - {{{\mathbf{v}}}_{l}}} \right|}^{2}},$
который как раз и является набором уравнений для центра сферы ${{{\mathbf{s}}}_{j}}$ из ${{\tilde {T}}_{j}}$. Заметим, что в общем случае функционал ${{F}_{I}}(X)$ не совпадает с $F(X)$, так как независимое проецирование граней на параболоид может приводить к вписанному многограннику, который может терять выпуклость на границах раздела проецируемых граней. Однако при сходимости ${{F}_{I}}$ совпадает с $F$.

Равенство ${{F}_{I}}(X) = 0$ означает, что для каждого радикального многогранника ${{D}_{i}}$ его двойственная вершина $*{{{\mathbf{c}}}_{i}}$ совпадает со всеми центрами шаров Делоне триангуляции Делоне его множества вершин, а значит, ${{D}_{i}}$ является многогранником Делоне.

Рассмотрим следующий алгоритм.

• Для $n = 0,\;1,\; \ldots $

• По набору шаров ${{\mathcal{B}}^{n}}$ строятся первичные и двойственные функции ${{{v}}^{n}}$ и ${{{v}}^{n}}{\kern 1pt} *$ с использованием подъема точек на параболоид. Эти функции определяют взвешенную триангуляцию Делоне ${{\mathcal{W}}^{n}}$ и радикальное разбиение ${{\mathcal{R}}^{n}}$ соответственно.

• Строится спроецированная функция ${{{\tilde {v}}}^{n}}{\kern 1pt} *$. Эта функция определяет локальную триангуляцию Делоне ${{\tilde {\mathcal{T}}}^{n}}$ множества вершин каждой ячейки радикального разбиения ${{\mathcal{R}}^{n}}$.

• Выполняется шаг минимизации на основе приближенного градиентного спуска для функционала Дирихле ${{F}_{I}}(X)$, чтобы получить множество шаров ${{\mathcal{B}}^{{n + 1}}}$.

• Процесс повторяется до сходимости.

В этом алгоритме последовательность взвешенных триангуляций Делоне ${{\mathcal{W}}^{n}}$ сходится к диаграмме Вороного $\mathcal{V}$, а последовательность радикальных разбиений ${{\mathcal{R}}^{n}}$ сходится к разбиению Делоне $\mathcal{T}$. Заметим, что предельная триангуляция Делоне $\tilde {\mathcal{T}}$ состоит из симплексов Делоне, а $\mathcal{T}$ – из многогранных ячеек Делоне для того же набора точек. На каждом шаге алгоритма нет необходимости строить согласованную триангуляцию Делоне. При необходимости ее можно построить один раз по готовому разбиению Делоне.

Для упрощения алгоритма, во избежание трудоемкого вычисления точного градиента, величины ${{c}_{i}}$ находятся с помощью локальной минимизации функционала (4), рассматривая его как квадратичную функцию от ${{{\mathbf{c}}}_{i}}$:

${\mathbf{c}}_{i}^{{{\text{new}}}} = {{\left( {\sum\limits_{{{{\tilde {T}}}_{j}} \in {{D}_{i}}} {{{{\mathbf{s}}}_{j}}} \operatorname{vol} {{{\tilde {T}}}_{j}}} \right) } \mathord{\left/ {\vphantom {{\left( {\sum\limits_{{{{\tilde {T}}}_{j}} \in {{D}_{i}}} {{{{\mathbf{s}}}_{j}}} \operatorname{vol} {{{\tilde {T}}}_{j}}} \right) } {\sum\limits_{{{{\tilde {T}}}_{j}} \in {{D}_{i}}} {\operatorname{vol} {{{\tilde {T}}}_{j}}} }}} \right. \kern-0em} {\sum\limits_{{{{\tilde {T}}}_{j}} \in {{D}_{i}}} {\operatorname{vol} {{{\tilde {T}}}_{j}}} }},$
в то время как значение ${{R}_{i}}$ пересчитывается с помощью простой аппроксимации методом наименьших квадратов
(5)
$R_{i}^{{{\text{new}}}} = \frac{1}{{\sqrt M }}{{{{{\left( {\sum\limits_{m = 1}^M {{{{\left| {{\mathbf{c}}_{i}^{{{\text{new}}}} - {{{\mathbf{v}}}_{m}}} \right|}}^{2}}} } \right)}}^{1}}} \mathord{\left/ {\vphantom {{{{{\left( {\sum\limits_{m = 1}^M {{{{\left| {{\mathbf{c}}_{i}^{{{\text{new}}}} - {{{\mathbf{v}}}_{m}}} \right|}}^{2}}} } \right)}}^{1}}} 2}} \right. \kern-0em} 2},$
где ${{{\mathbf{v}}}_{1}},\; \ldots ,\;{{{\mathbf{v}}}_{M}}$ – вершины ${{D}_{i}}$. Заметим, что соотношение (5) эквивалентно равенству
$\sum\limits_{m = 1}^M \tau ({{{\mathbf{v}}}_{m}}) = 0,$
что, в свою очередь, является необходимым условием минимума относительно ${{R}_{i}}$ для локального функционала
$\sum\limits_{m = 1}^M {{{\tau }^{2}}} ({{{\mathbf{v}}}_{m}})$
при условии, что центры зафиксированы. Новые положения и радиусы вычисляются с использованием релаксационного параметра $0 < \theta < 1$:
${\mathbf{c}}_{i}^{{n + 1}} = {\mathbf{c}}_{i}^{n}(1 - \theta ) + {\mathbf{c}}_{i}^{{{\text{new}}}}\theta ,$
$R_{i}^{{n + 1}} = R_{i}^{n}(1 - \theta ) + R_{i}^{{{\text{new}}}}\theta .$
Этот эвристический алгоритм достаточно эффективен для начальных итераций и хорошо работает в двумерном случае, как показано экспериментально. Для трехмерного случая эксперименты пока не проводились, возможно, нужно будет использовать метод градиентного спуска.

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

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

5. ЧИСЛЕННЫЕ ЭКСПЕРИМЕНТЫ

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

На фиг. 7 показаны начальное радикальное разбиение и конечное разбиение Делоне. Заметим, что ограничения удовлетворяются за счет появления больших ячеек вблизи фиксированных кругов, что, в свою очередь, сжимает некоторые свободные круги почти до нуля.

Фиг. 7.

Начальное радикальное разбиение ${{\mathcal{R}}_{0}}$ (а) и конечное разбиение Делоне $\mathcal{T}$ (б).

Обратим внимание, что регулярность слоев ячеек Делоне и Вороного вблизи внутренней границы нарушается из-за появления небольших ребер Делоне. В текущей версии алгоритма отсутствует механизм для их устранения. На фиг. 8 показано, как выглядят радикальные разбиения вместе с порождающими их кругами ${{B}_{i}}$.

Фиг. 8.

Начальное радикальное разбиение ${{\mathcal{R}}_{0}}$ с начальным набором кругов (а) и конечное разбиение Делоне $\mathcal{T}$ с кругами Делоне (б).

Начальная взвешенная триангуляция Делоне ${{\mathcal{W}}_{0}}$ и конечная сетка Вороного показаны на фиг. 9. Вновь можно наблюдать проблемы с качеством ячеек в конечной триангуляции Вороного. Заметим, что два слоя ячеек Вороного около внутренней окружности на самом деле являются слоями четырехугольников, они разбиты на треугольники при отрисовке.

Фиг. 9.

Начальная взвешенная сетка Делоне ${{\mathcal{W}}_{0}}$ (а) и конечная триангуляция Вороного (б).

Увеличенные фрагменты сеток на фиг. 10 позволяют наглядно увидеть, как работает алгоритм.

Фиг. 10.

Фрагменты начального радикального разбиения (а) и соответствующие фрагменты полученного разбиения Делоне (б).

На фиг. 11 показана эволюция фрагмента сетки с добавленными ортокругами, т.е. с кругами, центрированными на двойственных вершинах ${{{\mathbf{v}}}_{k}}$ с $\sqrt {\left| {\tau ({{{\mathbf{v}}}_{k}})} \right|} $ в качестве радиусов. В случае $\tau ({{{\mathbf{v}}}_{k}}) < 0$, красные окружности уже не являются настоящими ортоокружностями; они введены только для оценки отклонения радикальных ячеек от ячеек Делоне.

Фиг. 11.

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

Далее мы используем тот же алгоритм для построения логотипа конференции NUMGRID-2020 (см. [17]). Аббревиатура NG создается с помощью набора фиксированных защитных кругов, включенных в квадратную решетку из кругов. Вокруг букв разбросан квазислучайный набор кругов. Начальное радикальное разбиение и сошедшееся решение (разбиение Делоне) показаны на фиг. 12.

Фиг. 12.

Логотип NG: начальное радикальное разбиение ${{\mathcal{R}}_{0}}$ (а) и конечное разбиение Делоне $\mathcal{T}$ (б).

На фиг. 13. показано, как выглядят радикальные разбиения вместе с порождающими их кругами ${{B}_{i}}$. Обратим внимание, что конечные круги на самом деле являются кругами Делоне.

Фиг. 13.

Логотип NG: начальное радикальное разбиение ${{\mathcal{R}}_{0}}$ с начальным набором кругов (а) и конечное разбиение Делоне $\mathcal{T}$ с кругами Делоне (б).

На фиг. 14 показаны начальная взвешенная триангуляция Делоне и конечная триангуляция Вороного.

Фиг. 14.

Логотип NG: начальная взвешенная сетка Делоне ${{\mathcal{W}}_{0}}$ (а) и конечная триангуляция Вороного (б).

Наконец, на фиг. 15 показаны увеличенные фрагменты начального радикального разбиения и сошедшегося радикального разбиения, которое является разбиением Делоне.

Фиг. 15.

Логотип NG: фрагменты начального радикального разбиения (а) и конечного разбиения Делоне (б).

6. ОБСУЖДЕНИЕ

В двумерном случае мы численно показали, что радикальное разбиение может превращаться в многоугольное разбиение Делоне через эволюцию набора окружностей. Постановка задачи многомерна, поэтому ожидается, что алгоритм применим и в трехмерном случае. Хотя для этой задачи пока нет результата о существовании решения, мы не считаем это критическим недостатком. Как только добавление новых шаров становится допустимой операцией, результат существования становится тривиальным, поскольку спроецированная функция ${\tilde {v}}{\kern 1pt} {\text{*}}$ на каждой итерации является решением задачи. Однако проблема минимального добавления шаров для построения решения остается открытой. Численные эксперименты показывают, что способность локально добавлять новые круги/шары может быть также важна для достижения заданного качества сетки при наличии ограничений.

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

Идея использования значения дискретного функционала Дирихле для произвольного подъема (для меры грубости поднятой поверхности) не нова. В [19] был установлен принцип минимальной грубости для двумерных триангуляций Делоне. На данный момент не ясно, как свойство минимальной грубости может быть связано с представленными результатами.

На практике функционалы качества сетки должны быть оптимизированы с использованием многообразия $\nabla F = 0$ в качестве ограничения. Очевидные требования к качеству связаны с устранением малых ребер и граней Делоне и Вороного, что является предметом текущего исследования. Чтобы построить хорошие сетки Делоне–Вороного, необходимо следовать заданной функции локального размера, а также использовать стратегию, в которой устраняются малые ребра/грани Делоне, удаляются малые шары Делоне и малые ребра/грани Вороного, а также отдается предпочтение многогранной, а не симплициальной сетке Вороного.

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

  1. Delone (Delaunay) B.N. Proc. Inter. Congr. Math. (Toronto 1924) V. 1. P. 695–700. Univ. Toronto Press, 1928.

  2. Delaunay B.N. Sur la sphere vide // Bull. Acad. Sci. URSS. VII. Ser. 1934. V. 6. P. 793–800.

  3. Делоне Б.Н. Геометрия положительных квадратичных форм // Успехи матем. наук. 1937. Т. 3. С. 16–62.

  4. Делоне Б.Н. Геометрия положительных квадратичных форм // Успехи матем. наук. 1938. Т. 4. С. 102–164.

  5. Gärtner K., Kamenski L. Why do we need Voronoi cells and Delaunay meshes? In Garanzha V.A., Kamenski L., and Si H. (eds), Numerical Geometry, Grid Generation and Scientific Computing, Lect. Notes Comput. Sci. Eng. V. 131. P. 45–60. Springer, Cham, 2019.

  6. Гертнер K., Каменски Л. Зачем нужны сетки Вороного–Делоне? Основные свойства метода конечных объемов с использованием ячеек Вороного // Ж. вычисл. матем. и матем. физ. 2019. V. 59. № 12. P. 2007–2023.

  7. Si H. TetGen, a Delaunay-based tetrahedral mesh generator // ACM Transact. on Math. Software. 2015. V. 41. № 2. P. 1–36.

  8. Aurenhammer F., Klein R., Lee D.-T. Voronoi diagrams and delaunay triangulations. World Sci. Publ. Co., Inc., USA, 2013.

  9. Garanzha V.A., Kudryavtseva L.N., Tsvetkova V.O. Structured orthogonal near-boundary Voronoi mesh layers for planar domains. In Garanzha V.A., Kamenski L., Si H. (eds), Numerical Geometry, Grid Generation and Scientific Computing, Lect. Notes Comput. Sci. Eng. V. 131. P. 25–44. Springer, Cham, 2019.

  10. Гаранжа В.А., Кудрявцева Л.Н., Цветкова В.О. Построение гибридных расчетных сеток Вороного. Алгоритмы и нерешенные проблемы // Ж. вычисл. матем. и матем. физ. 2019. V. 59. № 12. P. 2024–2044.

  11. Abdelkader A., Bajaj C., Ebeida M., Mahmoud A., Mitchell S., Owens J., Rushdi A. VoroCrust: Voronoi meshing without clipping // ACM Trans. Graph. 2020. V. 39. № 3. P. 23.

  12. Voronoi G.F. french Nouvelles applications des paramétres continus a la théorie de formes quadratiques // J. Reine Angew. Math. 1908. V. 134. P. 198–287.

  13. Edelsbrunner H., Seidel R. Voronoi diagrams and arrangements // Discrete Comput. Geom. 1986. V. 1. P. 25–44.

  14. Александров А.Д. Выпуклые многогранники. Mосква–Ленингpад: Гостехиздат, 1950.

  15. Fenchel W. On conjugate convex functions // Canad. J. Math. 1949. V. 1. P. 73–77.

  16. Edelsbrunner H. Geometry and topology for mesh generation. Cambridge Monographs on Applied and Computational Mathematics. Cambridge Univer. Press, 2001.

  17. Garanzha V.A., Kamenski L., Si H. (eds) Numerical geometry, grid generation and scientific computing. V. 143 of Lect. Notes Comput. Sci. Eng., Springer, Cham, 2021.

  18. Alexa M. Conforming weighted Delaunay triangulation // ACM Trans. Graph. 2020. V. 39. P. 6.

  19. Rippa S. Minimal roughness property of the Delaunay triangulations // Comput. Aided Geom. Design. 1990. V. 7. P. 489–497.

Дополнительные материалы отсутствуют.