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

Вместе с описанием новой функциональности в этой статье отдельное внимание уделяется численным методам и подходам к поиску точек срединных поверхностей.

Определение и область применения

Срединная поверхность — это геометрическое место точек, равноудаленных от двух указанных совокупностей граней. В нашем случае на совокупности граней наложены ограничения: они должны быть связны и иметь гладкость G1. В этом случае мы будем ожидать связную срединную поверхность между ними.

Чаще всего срединные поверхности используются для построения расчетных моделей в виде оболочек в CAE-системах. Также они применяются для разделения сложных моделей на части или для построения симметричных моделей.

В первом разделе рассмотрим схему построения сплайновых срединных поверхностей.

Общий алгоритм

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

Идея метода заключается в триангуляции исходной поверхности и отображении этой триангуляции на срединную поверхность. Топология сетки при ее отображении сохраняется, что позволяет сразу получить первичное сеточное представление срединной поверхности и далее использовать ее для построения сплайновой аппроксимации.

На первом этапе алгоритма выполняется построение сетки на базовом наборе граней.

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

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

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

Завершает алгоритм адаптивная аппроксимация сетки сплайновой поверхностью. Полученная сплайновая поверхность возвращается пользователю как результат построения срединной поверхности.

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

Поиск точки срединной поверхности

Пусть дана точка А на базовой поверхности s1 и нормаль в этой точке n1. Ответную поверхность обозначим s2. Нужно найти сферу, которая касается поверхности s1 в точке А и поверхности s2 в некоторой точке В. Центр этой сферы будет являться искомой точкой срединной поверхности. Схематично постановка задачи изображена на следующем рисунке. На нем и далее поверхности и сферы для простоты восприятия будут показаны в   плоском сечении, но будем помнить, что речь идет о трехмерной задаче.

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

Для решения оптимизационной задачи нужно определиться с тремя составляющими:

• параметризация, которая определит пространство решений;

• целевой функционал, минимум которого будет отвечать искомому решению;

• метод минимизации — алгоритм, по которому будет производиться поиск минимума в выбранном пространстве решений.

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

Первый подход

Поскольку поставленная перед нами задача геометрическая и никакой вероятностной или случайной природы не содержит, статистические методы минимизации здесь для практики не подойдут. Нужно использовать всю имеющую геометрическую информацию. Например, очевидно, что центр сферы будет лежать на прямой n1 — нормали в точке А. Тогда в качестве искомого параметра можно выбрать r — радиус сферы, которая касается s1 в точке А и центр которой лежит на луче n1.

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

Определимся с целевым функционалом. Пусть дан некоторый радиус ri и соответствующая ему точка Ci(ri) на n1. Для того чтобы сфера с центром в Ci касалась s2, нужно, чтобы расстояние от Ci до s2 равнялось радиусу ri или, что то же самое, расстоянию от Ci до А. Вычислим проекцию точки Ci на s2, и обозначим ее как Bi. Функционал тогда можно взять как разность расстояний АСi и СiВ.

Поскольку это одномерный функционал, минимизировать его можно, по большому счету, любым одномерным методом, например методом золотого сечения, но лучше, конечно, выбрать один из градиентных методов, чтобы в ходе итераций использовать больше доступной информации о геометрии модели. Но проблема в том, что для расчета шага минимизации производную в данном случае придется оценивать численно, поскольку функционал содержит операцию проецирования и, кроме того, само наличие операции проецирования в расчете функционала вносит в его расчет некоторую, пусть и незначительную, погрешность. Это приводит к тому, что данный метод не обладает достаточной устойчивостью и требует довольно большого количества итераций (в данном случае 50-100 и больше), что существенно замедляет поиск. Поэтому для практики выбор этого подхода нельзя считать правильным.

Второй подход

Во втором подходе будем искать положение точки B на поверхности s2. Возьмем в качестве параметров координаты uv точки В. Минимизация будет двумерная.

Через любую точку B на s2 можно построить сферу так, что она будет касаться точки А и ее центр будет лежать на n1. Но в каком случае эта сфера будет касаться s2? Это будет тогда, когда нормаль в точке В к поверхности s2 будет совпадать с нормалью в этой точке к поверхности сферы. На рисунке выше эти нормали обозначены n2 и n3 соответственно. Целевой функционал можно сформулировать как квадрат разности векторов n3 и n2 (на рисунке вектор разности обозначен f). Несмотря на то что эта постановка двумерная, в отличие от предыдущей одномерной постановки, она оказывается более выигрышной. Операции вычисления нормали к поверхности и ее производные в геометрическом ядре c3d рассчитываются аналитически. Это означает, что значение функционала и, что важно, его производные будут рассчитаны точно. Целевой функционал, в свою очередь, является сумой квадратов разностей, а для таких, выпуклых, функционалов существует метод минимизации второго порядка Гаусса-Ньютона, для расчета шага минимизации которого требуются только первые частные производные функции, которая стоит под квадратом в функционале. На практике этому методу для поиска касающейся сферы требуется всего 3-4 итерации, за которые он понижает функционал до 10^-12 и ниже. Но у метода Гаусса-Ньютона есть одно важное ограничение. Он гарантирует сходимость только в окрестности решения, и это означает, что ему нужны хорошие начальные приближения.

Множественность решений

На этом моменте вернемся к постановке задачи. Мы не упомянули важный факт: решений может быть несколько и довольно произвольное количество. Сфера может касаться поверхности s2 во многих точках, в том числе и на ее продолжении.

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

• Соответствующая ей касающаяся сфера будет иметь наименьший радиус.

• Сфера будет иметь внутреннее касание поверхности s2, что означает, что радиус сферы будет меньше радиуса кривизны s2 в точке B с учетом знака. Пример внутреннего касания на рисунке выше показан в точке B1, а внешнего — в точке B2.

• Сфера касалась саму поверхность, а не ее продолжение.

• Имела тупой угол между векторами СА и СВ.

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

Поиск начальных приближений

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

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

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

Примеры построений

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

Построение срединной поверхности между группами граней, сверху — трехмерный вид, снизу показано поперечное сечение.

Построение замкнутой срединной поверхности.

Вместо выводов

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

Илья Патрушев

Ведущий математик-программист
C3D Labs