CC BY
34Моделирование деформаций упругих объектов в виртуальных окружениях в реальном времени, используя метод конечных элементов и предвычисленные функции Грина
Никитина Л.Д., Никитин H.H., Фролов П.В., Геббельс Г., Гёбель М. Фраунгоферовский Институт Медиакоммупикаций, Санкт Августин, Германия
Клименко C.B. Институт Физики Высоких Энергий, Протвино, Россия
Нильсон Г.
Аризонский Университет, Тэмпе, Аризона, США
Аннотация
Моделирование упругих деформаций объектов чрезвычайно важно для приложений виртуального окружения, в которых требуется реалистично представлять поведение трехмерных объектов. В случае взаимодействия пользователя с объектом особенно желательно вычислять эти деформации так быстро, как возможно. В данной работе мы представляем новый подход для моделирования упругих объектов в виртуальных окружениях в условиях реального времени. Наш подход использует метод конечных элементов и предвычисленные функции Грина. Деформации интерактивно визуализируются в виртуальных окружениях обратной проекции с полным погружением, таких как CyberStage и в системах с частичным погружением, таких как Responsive Workbench. Взаимодействуя с объектом посредством виртуального луча, пользователь может интерактивно прикладывать силы к объекту, производя его деформации.
Ключевые слови: упругие деформации, виртуальные окружения, метод конечных элементов, функции Грина.
1 Введение
В большинстве приложений моделирования, таких как инженерные, медицинские и т.д. имеется большая необходимость в подходах, позволяющих моделирование материальных свойств объекта и его поведения под действием приложенных к нему сил. В медицинских приложениях требуется моделирование мягких тканей; в этой области основоположниками являются Бро-Нильсен и Котин [1], которые использовали метод конечных элементов для моделирования упругих деформаций. Технические приложения нуждаются, в основном, в моделировании упругих деформаций объектов, компоненты которых состоят из изотроп-
ных материалов, в условиях применимости закона Гуна [2]. Обзорная работа [3] описывает большую часть предшествующих работ по моделированию деформируемых объектов. Недавними достижениями являются применение метода граничных элементов и формулы Шермана-Моррисона-Вудбари для моделирования линейно-упругих однородных объектов [4], и использование итеративных динамических схем [5, 6] для моделирования нелинейной упругости. Для нелинейно-упругих объектов достаточно высокой детализации до сих пор не удается моделировать деформации в реальном времени. Для упрощенной, линейной формулировки теории упругости существует возможность решать системы уравнений достаточно большого размера в реальном времени синхронно с графической обработкой, или, для ускорения метода, часть этих вычислений может быть выполнена предварительно. Также в зависимости от внутреннего строения объекта можно использовать метод конечных элементов или метод граничных элементов. В данной работе мы коротко обсудим эти различные подходы и затем опишем методы и техники, использованные в нашей системе, которые концентрируются на применении линейной теории упругости, метода конечных элементов и предвычислении полного базиса решений (функций Грина). Система разработана в программной среде Аванго и используется для моделирования в реальном времени упругих деформаций трехмерных объектов в виртуальном окружении, позволяя интерактивное, определенное пользователем приложение силы.
2 Системы виртуального окружения
Предназначение виртуальных окружений заключается в том, чтобы обеспечить отдельных пользователей, или группы инженеров, дизайнеров, медиков виртуальным рабочим пространством, в котором они могут наблюдать, исследовать и создавать в реальном времени необходимые им виртуальные данные. Большинство таких виртуальных окружений имеют сходную аппаратурную конфигурацию. Центральное место занимает графический обработчик: SGI суперкомпьютер или кластер Linux PC, который обрабатывает данные и отображает их на экран, а также устройство слежения, которое измеряет положение и ориентацию головы пользователя. Данные, произведенные устройством слежения, читаются графическим обработчиком, чтобы определить перспективно правильное изображение для любой точки зрения пользователя. Аппаратурная конфигурация, используемая для интерактивной визуализации деформируемых объектов в подходе данной работы включает устройства ввода для взаимодействия и стереоскопические системы обратной проекции в качестве устройств вывода. Далее мы описываем наиболее часто используемые проекционные системы Responsive Workbench и CyberStage, а также систему разработки виртуальных окружений Аванго, используемую для визуализации.
Рис.1. Двусторонний Responsive Workbench.
Responsive Workbench™ - это стереоскопическая проекционная система высокого разрешения, разработанная в отделе виртуальных окружений Немецкого Национального Исследовательского Центра Информационных Технологий [7]. Responsive Workbench использует CRT-проектор, который отображает стереоизображения через зеркало на плоскость экрана. Для каждого глаза вычисляется соответствующее изображение, графический обработчик быстро чередует эти
два изображения на экране. Пользователи носят жидкокристаллические шаттр-очки, у которых левое стекло становится непрозрачным, в то время как изображение для правого глаза отображается на экране, и наоборот. Тем самым производится стереоскопический эффект. Электромагнитный датчик фирмы Polhemus, у которого измеряются все 6 степеней свободы, т.е. положение и ориентация, прикреплен к шаттр-очкам. Это позволяет системе вычислить перспективно правильное изображение для любого положения пользователя. Пользователи взаимодействуют напрямую с 3-мерными виртуальными объектами, находящимися в поле зрения, используя устройства ввода, у которых аналогичным образом измеряются положение и ориентация. Типичными областями приложений Responsive Workbench являются инженерная и медицинская визуализация.
Новый двусторонний Responsive Workbench имеет горизонтальный и вертикальный экраны, которые гладко примыкают друг к другу (рис.1). Дополнительный вертикальный экран значительно увеличивает поле зрения и позволяет наблюдать виртуальные объекты на уровне глаз пользователя, что ранее не было возможно. Используя современные персональные компьютеры и графические рабочие станции, двусторонний Responsive Workbench производит стерео-изображение с разрешением 1280x1024x96Гц. Благодаря расширенному полю зрения двусторонний Responsive Workbench используется в более широком диапазоне приложений, чем односторонний.
CyberStage™ - это проекционная система объемного изображения и звука, в которой достигается эффект полного погружения [8, 9]. Она имеет размеры небольшой комнаты (ЗмхЗмх2.4м) и включает 4-стороннюю проекцию стерео-изображения и 8-канальную проекцию пространственного звука, контролируемые положением головы пользователя. Суперкомпьютер SGI Опух 2 с четырьмя графическими картами Infinite Reality 2 и двенадцатью процессорами MIPS R12000 производит четыре изображения с разрешением 1280x1024x120Гц. Каждая карта в пиковом режиме генерирует 11 млн. треугольников в секунду.
Для создания иллюзии присутствия в виртуальных пространствах, система CyberStage предоставляет различные интерфейсы и метафоры взаимодействия, визуально и акустически откликаясь на действия пользователя. Эти интерфейсы используются для навигации в виртуальных пространствах и для взаимодействия с виртуальными объектами. Программное обеспечение для CyberStage и Responsive Workbench - это система разработки виртуальных окружений Аванго [10].
Аванго™ - среда программирования, предназначенная для создания распределенных взаимодействую-
щих приложений виртуального окружения. Аванго использует язык программирования С++ для определения двух категорий объектов - узлы и датчики. Узлы составляют объектно-ориентированный граф сцены, который осуществляет представление изображаемой геометрической модели. Датчики обеспечивают связь с реальным миром и используются для ввода данных с внешних устройств в приложение. Все объекты Аванго - полевые контейнеры, представляющие информацию о своем состоянии в виде совокупности полей. Объекты Аванго обладают универсальным потоковым интерфейсом, который позволяет записывать в поток сами объекты и информацию об их состоянии и впоследствии реконструировать эти данные из потока. Аванго использует язык сценариев Scheme [11] для объединения объектов в граф сцены, а также вводит понятие графа потоков данных, концептуально ортогонального графу сцены и служащего для определения взаимодействия между объектами и для ввода в сцену реальных данных, необходимых для моделирования интерактивного поведения. Scheme - универсальный язык программирования, происходящий от Алгола и Лиспа. Это язык высокого уровня, который умеет оперировать со структурными данными, такими как строки, списки и векторы. Все объекты Аванго могут быть созданы и управляемы посредством Scheme. Аванго основан на графической системе (rix OpenGL Performer, которая обеспечивает наиболее качественное исполнение приложения, и удовлетворяет специальным требованиям, возникающим при разработке приложений для виртуального окружения. Расширенные графические возможности, такие как стирание невидимых граней, переключение степени детализации и взаимодействие с графическим оборудованием, целиком выполняются системой Performer. Если аппаратное обеспечение позволяет, Performer использует множественные процессоры и множественные графические конвейеры. Мы используем реализацию Аванго на графических станциях SGI и Linux-PC.
Наша цель - представить деформируемые модели в виртуальных окружениях. В современном математическом моделировании это обычно делается с помощью дискретных методов, таких как мет,од конечных элементов. При этом тело разбивается на конечное множество элементов, и определяются условия физического равновесия для каждого из элементов. Эти условия записываются в виде очень большой системы уравнений с десятками тысяч неизвестных. Такая система обычно является вырожденной, но после наложения соответствующих граничных условий становится хорошо определенной и имеет однозначное решение. Для нахождения решения необходимо чрезвычайно интенсивное вычисление, требующее большого объема оперативной памяти и времени исполнения. Для приложений вирту-
альной реальности желательно оптимизировать алгоритм решения, в особенности по отношению к его скорости исполнения. В следующих разделах мы приводим необходимые детали о методах, используемых для моделирования деформируемых объектов в реальном времени.
3 Математические основы метода
Теория упругости [2] характеризует деформации упругих объектов с помощью смещений й(х), преобразующих точки объекта следующим образом: х —► х + й(х). Тензор деформации определяется частными производными: е^ = | + + где предполагается суммирование по повторяющимся индексам. Тензор деформации характеризует изменение расстояний в деформированном объекте, так что длина малого элемента ёх преобразуется следующим образом: \сЩ2 -Л \с1х\2 + 2щйхгйх^. В случае малых смещений можно сохранить в тензоре деформации только линейные члены, записывая просто еу = | + , результирующее приближение называется теорией малых деформаций. Удобно поместить компоненты тензора деформации в столбец е = (ехх,еуу,егг,2еху,2ехг,2еуг)Т, и записать его определение через смещения как е = Аи, где Д - символическая матрица дифференцирования:
€хх еУУ
f'ZZ 2 €ху 2 €xz
дх 0 0
0 JL ву 0
0 0 в dz
д_ JL 0
ву дх
JL 0 JL
dz дх
0 JL JL
dz ву
(1)
Распределение сил в упругом теле характеризуется следующим образом. Рассмотрим воображаемый срез тела и малый элемент с1Б на нем (этот вектор направлен вдоль нормали к срезу, модуль \clSl равен площади элемента). Сила df, действующая сквозь этот элемент между частями тела, разделенными срезом, определяется как = а^йБу Коэффициент иц называется тензором напряжения. Он симметричен: оц = также как тензор деформации, и может быть записан аналогично: а = ((тхх,ауу,агг,аху,ахг,ауг)т.
Тензоры напряжения и деформации связаны законом Гука: <7 = О«-, где матрица не зависит от деформаций и характеризует свойства материала. Это линейное соотношение справедливо для малых деформаций при дополнительном предположении, что неде-формированное тело обладает нулевым тензором напряжения, т.е. материал не является преднапряжен-
ным. Для однородного изотропного материала матрица Ю представляется в виде линейной комбинации I) = А£>1 + где 1?1,2 - постоянные 6x6 матрицы:
Di
1 1 1 1 1 1 1 1 1
D2
2 2 2
1 1 1
А, ГПа ц, ГПа
керамика сталь резина 124.12 119.42 1.12 124.12 79.23 0.047
S 2 1
S 4
3
2 T Si
Рис.2. Тетраэдральный элемент.
где V - объем тетраэдра и Sk - площади его граней, направленные по внешним нормалям (рис.2). Явные формулы, выражающие V, Sk в терминах координат узлов хк:
(пустые места заполнены нулями) и А, ц - постоянные Ламе, например
Si : $2 : S3 :
S 4 :
aijk
№ i{d L
(^123> ^123> ^12з)>
234> "234> "234 h 3411^341! ^341 )i
V ■■
1111 Xi X2 X3 X4 VI У2 УЯ УА Z\ Z2 Z3 Z4
1 1 1 1 1 1 1 1 1
У г Уз У к II Zi Zj ZL > aijk - Xj, Xj xk
Zj Xi Xj Xfe Уг У, Ук
(ср. с [13, 1]). Частные производные записываются сле-
4
дующим образом: j^- = —3V S3ku\ и 1 к=1
Энергия деформированного тела задается интегралом Е = /с1У(^етсг + II), где первое слагаемое представляет чисто упругую энергию, а второе ассоциируется с потенциальным полем внешних сил / = — VII (например, поле тяжести). Условия равновесия, соот-
дац л
ветствующие минимуму энергии, имеют вид = В терминах неизвестных смещений и это дает систему линейных дифференциальных уравнений в частных производных 2-го порядка.
При подстановке в правую часть этой системы сил вида /¿(ж) = 5{к6(х — х,о), где ё^к - символ Кронекера и ё{) функция Дирака, решение системы для данных а?о,к: щ(х;хо,к) = хо) называется функцией Грина системы [12]. Эта функция описывает, как воздействие силы, приложенной локально в точке хо, распространяется в точку х. Если эта функция известна, решение системы для произвольных правых частей /(ж) можно найти как щ{х) = [d3xoGiк(x,xo)fk(xo)■
Метод конечных элементов [13] вводит дискретизацию тела, используя, например, разбиение тела на тетраэдры, и представляет непрерывное поле смещений й(х) с помощью линейной интерполяции величин, заданных в узлах:
4
= ^2мк(х)йк, Мк(х) = 1 - 3^7Зк(х-Хк), к=1
^xx
еУУ ^zz = bV) Si S2 S3 S4
^xz Si S2 S3 S4
Ц-. А
где йк - столбцы (ukx,uky,ukz)T и Sk, Sk обозначают 3x3 блоки
Sk
skx 0 0 sky skx 0
0 sky 0 sk = Skz 0 Skx
0 0 Skz 0 Skz Sky
В этом соотношении Д - дискретизованный аналог матрицы дифференцирования из уравнения (1), ассоциированный с тетраэдрическим разбиением.
Замечание: для проверки трансляционной инвариантности системы, мы подставляем в это соотношение постоянное смещение йк = с и немедленно получаем е = О в силу тождества: Si + S2 + S3 + S4 =0 (геометрическая теорема "о еже"). Аналогичным образом можно показать, что тензор деформации обращается в нуль при инфинитезимальных вращениях, т.е. при линейных преобразованиях вида йк = aft х хк, где а - угол и п - ось вращения.
Упругая энергия для одного тетраэдра записывается как /•;,„ = Ijit1 Кi,iи. Кi,i = ATDAV, после суммирования по всем тетраэдрам и добавлении потенциальной энергии это дает Е = |итКи + U, где К = YLtei transfer(iftet) называется матрицей жесткости. Функция "transfer" [1] выполняет переход от локальных индексов (1234), помечающих вершины одного тетраэдра, к глобальным индексам, помечающим узлы модели. Матрица К симметрична и положительно
S
3
полуопродолотта: Кт = К, итКи > 0, она имеет нулевые направления К и = 0, соответствующие описанным выше степеням свободы твердого тела (трансляциям и иттфиттитозималытым вращениям). Чтобы исключить эти степени свободы, должны быть зафиксированы граничные условия, для чего достаточно положить и = 0 в п > 3 узлах, но лежащих па одной прямой. Устранение строк и столбцов, соответствующих этим узлам, из матрицы К делает ее положительно определенной: итКи >0, и ф 0, что соответствует строгому минимуму энергии. Условие минимума имеет вид системы линейных уравнений К и = /, которая должна быть ротона: и = КФункция Грина в этих обозначениях равна обратной матрице А"-1, в которой (у)-эломопт описывает, как воздействие единичной силы, приложенной к _7~тому узлу, распространяется в г-тый узел.
В более общем виде можно рассмотреть границу тола Г = дП и область границы Г^ с Г, где заданы ненулевые смещения: и = Ъ. Ставится задача о нахождении смещения и па Г\Г(, (и, если требуется, во внутренних областях П). Такие граничные задачи сводятся к задачам с заданной силой, используя следующее разбиение матрицы А":
троппых материалов. Полученные интегральные уравнения являются сингулярными, и для их численного решения требуются специальные регуляризации. В методе граничных элементов необходимо вводить разбиение только поверхности объекта, практически можно использовать ту же триангуляцию, что и для графической обработки. В методе коночных элементов разбиение необходимо вводить также во внутренних областях объекта, даже в том случае, если по требуется проводить визуализацию этих областей. Применение метода граничных элементов ограничено объектами, внутренность которых однородна и изотропна. Для объектов со сложной внутренней структурой необходимо использовать метод коночных элементов. Метод коночных элементов порождает большие и разреженные линейные системы, в то время как эквивалентные им системы в методе граничных элементов имеют моиыпий размер и являются плотными. Оба метода могут использоваться в режиме предварительных вычислений. В области совмостой применимости они производят один и тот же результат - функцию Грина па поверхности. Этот результат можно затем использовать для ускорения вычислений в реальном времени. В дальнейшем мы описываем такую оптимизированную схему для метода коночных элементов.
А В
ВТ С
Л
ГД° /б ~ сила, действующая па Г;,, в то время как для других узлов предполагается отсутствие внешней силы. Отсюда мы получаем следующие две системы: Вти + СЬ = Д, определение силы /(,, которую можно использовать, например, для генерации тактильных ощущений |14|; Аи = —ВЬ, система па и, непосредственно используемая при визуализации деформаций объекта. Дальнейшее упрощение |1| возможно, если изображается только грйницЕ (не внутренне области) тола. В последней системе участвуют смещения и как па границе, так и для внутренних узлов, однако, когда эта система ротона: и = —А^ВЬ, необходимо сохранить в матрице А^1В только тот блок, который соответствует граничным узлам, т.о. функцию Грина, описывающую распространение воздействия из граничных точек в граничные точки.
Замечание: к аналогичной формулировке приводит использование метода граничных элементов |4|. Этот метод использует аналитическую переформулировку исходных дифференциальных уравнений в частных производных к интегральному виду, который включает только поверхностные переменные (смещения и распределения сил по поверхности). Такая переформулировка возможна только для случая однородных изо-
4 Оптимизированная схема вычислений
Хранение матрицы жесткости.
N V е г 1 ~ 103..104
Рис.3. Хранение матрицы жесткости.
Матрица К обычно имеет очень большой размер: КсМхКсМ, \(1оГ ЗХуопр. Куойя ~ 103..104, однако она является весьма разреженной, содержащей обычно более 99% пулой. Эффективное храпение этой матрицы возможно в силу следующего наблюдения: элемент Ку с г ф j можно представить как отрезок, соединяющий г^'-тыо узлы, и если этот отрезок по совпадает с ребром одного из тетраэдров в разбиении, соответствующий элемент обращается в пуль. В результате этого,
недиагональные элементы Кц можно хранить на ребрах разбиения, в то время как диагональные элементы К и могут храниться в узлах. Такой способ хранения требует объема памяти Д/д для матрицы К, а также дополнительного объема Мд для хранения данных о разбиении (координаты узлов, индексирующие массивы) и необходимых рабочих массивов:
Mx=(6Nverts+9Nedges)-sizeof(double) ~ 600Kb ... 6МЬ
M0=18(Nverts)-sizeof(double)+(2Nedges+10Ntets)
•sizeof(long int) ~ 440Kb ... 4.4Mb.
Произведение g = К и, необходимое в дальнейших вычислениях, можно найти, используя следующий алгоритм:
для каждого узла i { f)i = Кцщ }
для каждого ребра ij { §i+ = KijUj; gj+ = КТ.щ }
Замечание: В узлах хранятся симметричные 3x3 матрицы К и] ребра рассматриваются как ориентированные отрезки, в которых хранится матрица Кц, если ij совпадает с направлением ребра, и транспонированная матрица К?, если ij противоположно этому направлению.
Решение линейной системы Ки = / может находиться, используя либо явные, либо итеративные методы [14]. Явные методы, такие как метод Гаусса, LU-разложение и др. требуют значительного времени для вычислений и большого объема оперативной памяти. В частности, один шаг метода Гаусса в терминах нашей модели хранения данных можно представить как соединение между собой всех узлов, связанных с г-тым узлом, и устранение /-того узла. На промежуточных стадиях этого процесса почти все оставшиеся узлы будут связаны между собой, приводя чрезвычайно большому потреблению памяти.
Среди итеративных методов наиболее часто используется метод сопряженных градиентов (СГ). Этот метод, по существу, является последовательной минимизацией энергии вдоль направлений {pi}, обладающих свойством if-ортогональности: pjKpj = 0, i ф j. Теоретически, для квадратичных энергий, через конечное число шагов Nsteps=Ndof эта последовательность исчерпывает пространство решений, и алгоритм останавливается на точном решении. На практике Nsteps« Ndof дает достаточную точность. Обычно используется предварительно обусловленная (preconditioned) версия метода сопряженных градиентов (ПСГ), обладающая лучшей сходимостью [14].
Сравнение разных схем вычисления. Описанные методы могут применяться в реальном времени для вычисления решения, одновременно с графической обработкой, или часть этих вычислений может быть
выполнена предварительно. Обе схемы обсуждались в [1] и мы также опробовали обе. Схема вычислений в реальном времени больше подходит для малых моделей. В частности, Nverts ~ 200 могут быть обработаны в реальном времени (20 кадров в секунду), Nverts ~ 3000 требуют порядка 1 сек релаксационного времени для полной визуальной сходимости итераций, в то время как для Nverts~104 это время составляет около 1 мин (для машин типа SGI/Onyx2 300MHz MIPS R12000). Метод сходится быстрее для тел, размеры которых во всех трех измерениях приблизительно равны (причины этого обсуждаются в приложении), и для материалов, обладающих Л < ¡л, по свойствам сходных с керамикой.
В работе [5] использовалась итеративная динамическая схема для моделирования деформаций объектов в рамках линейной и нелинейной теории упругости. В этой схеме вычисляются упругие силы, действующие на узлы модели, и дальнейшая эволюция системы просчитывается в реальном времени с помощью стандартных методов интегрирования по законам Ньютона, моделируя должным образом распределения масс и сил трения. При этом нелинейные члены в упругой силе, отвечающие квадратичному члену в тензоре деформаций (см. выше), значительно увеличивают объем вычислений. В работе [5] для объекта, содержащего 1,400 узлов, была достигнута скорость вычислений 45 итераций в секунду для линейной модели и 8 итераций в секунду для нелинейной модели, на компьютере типа PC Pentium III 500 MHz. В то же время, по нашим данным, для сходимости итераций до относительных погрешностей е ~ 10 :i (достаточных для визуализации) требуется от десяти итераций до нескольких тысяч итераций, в зависимости от плавности взаимодействия, формы и материальных свойств объекта. Это приводит к релаксационным эффектам той же длительности, как и полученные нами, которые внешне сходны с физическими процессами затухания колебаний, но протекают медленнее, чем в реальности, из-за недостаточной мощности вычислителя.
Другое свойство, которое обычно рассматривается как преимущество вычислений в реальном времени, это возможность выполнять интерактивные изменения системной матрицы К, включая те, которые вызваны вариацией граничных условий и изменением топологии модели (разрезы). Однако, В НбДЙВНИХ работах [4, 15] эти возможности были также реализованы для предвы-численных моделей, используя ускоренный алгоритм обращения матрицы К при ее локальных изменениях.
В нашей работе мы производим предварительное обращение матрицы К и представление решений в виде и = К^1 /, в результате получаются следующие характеристики: объект из 10,000 узлов деформируется при скорости графики 85 кадров в секунду на Linux PC
Athlon 1.3 GHz. Метод является квазистациоттарттым, т.е. па каждый кадр находится точное состояние равновесия объекта, отвечающее данному взаимодейстию. Основная проблема этого подхода состоит в том, что матрица К^1 больше пе является разреженной и требует довольно много памяти для храпения. К счастью, пе всю матрицу необходимо хранить, а только те строки, которые отвечают видимым узлам, и только те столбцы, которые отвечают узлам, в которых может быть приложена внешняя сила:
*
* *
Is
Для большинства приложений достаточно выделить из К-1 блок (K-^s размера (3Nsvert) х (3Nsvert), соответствующий поверхностным узлам. Для этой цели мы решаем систему К и = / с /.; = ё^, где к ограничено па поверхностные узлы. Решение вычисляется используя метод ПСГ. Заметим, что внутренние узлы дают вклад в элементы К^1 во время предварительной стадии вычислений, и затем могут быть опущены. Полученные решения (столбцы Кпишутся в файл, для i £ поверхностным узлам. Далее, используя симметрию {K-^s и уменьшая точность матричных элементов до float (двойная точность необходима для устойчивости предварительного вычисления, в то время как одинарная точность достаточна для визуализации), мы уменьшаем объем памяти, необходимой для храттия (К-^з до приемлемой величины:
Ms=4.5(Nsverts)2-sizeof(float) ~ 6МЬ ... 450МЬ,
для числа поверхностных узлов в интервале Nsverts ~ 600 ... 5000.
Вычисление произведения из = (К'~1)sfs, которое необходимо выполнять в реальном времени, в общем случае все еще является интенсивным. Дальнейшая оптимизация возможна, если использовать дополнительные предположения о характере силовых распределений. В частности, если сила локализована в нескольких узлах поверхности, можно использовать ускоренный алгоритм умножения |1|. Для силы, локализованной в единственном узле, этот алгоритм сводится к извлечению соответствующего столбца из (КВ другом случае сила может быть распределена в широких областях поверхности, по вводится небольшое число типов силовых распределений, которые представляют различные типы взаимодействия (растяжение, сдвиг, кручение и т.п.). В этом случае имеется небольшое число линейных комбинаций столбцов которые можно вычислить предварительно, составить из них
новую матрицу А, сохранить ее, и затем комбинировать ее столбцы в реальном времени, используя короткое вычисление (в простейшем случае сводящееся к извлечению необходимого столбца из матрицы А). Еще один пример, появляющийся в инженерных приложениях, это деформация упругого тела, определяемая по нескольким контрольным элементам, прикрепленным к границе. Мы рассмотрим эту модель более детально.
Упругие деформации с помощью контрольных элементов
Рассмотрим граничную задачу Аи = — ВЪ, где смещения Ъ па участке границы задаются аффинными преобразованиями Ъ = Т + МЬо. Эти преобразования описываются 12 параметрами, а именно, 3-мя переносами Т и 3x3 матрицей общих линейных преобразований М
o f f - l i n e
N v e r t ~1 0 0 0 0 3 s e c
f o r ( i = 0 ; i < d o f ; i + + )
£=10 t ~1 m i n / d o f
Nsvert~5400
c r e a t e g e o m
b u i l d s y s
ш
s e t b c
o n - l i n e ( A v a n g o )
Mat1 у_у Mat2
f i e l d H a s C h a n g e d : f i n d r,.
6)
t i .
J
хз Е ck3v'm
7
Рис.4, а) схема вычислений; б) базис решения (¡у-файл); в) оптимизация вычислений, выполняемых в реальном времени.
Рис.5, интерактивно деформируемая модель трубки в демонстрационном приложении Аванго.
(включающей, в частности, 3 вращения), определяя 12-мерттое пространство па контрольный элемент. Пусть к = 1..12 - базис в этом пространстве, соответствующий единице, помещаемой последовательно в ячейки Т и М, в то время как остальные ячейки заполнены пулями. Пусть - соответствующие решения системы, т.е. = —А^ВЬ^'К Тогда, в силу линейности системы, форма упругого тела для данного взаимодействия Ь = Y^k будет задаваться выражением и = fkU^ ■ Фактически, нам необходимо вычислить решения для элементарных взаимодействий, соответствующих базисным векторам
,в предварительном режиме, и затем комбинировать эти решения, чтобы найти форму
любого заданного взаимодействия, используя короткое вычисление в реальном времени. Объем памяти па контрольный элемент, требуемый для храпения решений мал:
Mu=36Nsverts-sizeof(float) ~ 100Kb ... 700Kb.
Данная схема вычислений показана па рис.4а. Для модели с числом узлов Nverts~10,000 инициализация системы занимает порядка 3 сек, затем система решается в цикле для различных базисных векторов
с относительной погрешностью е. = 10_6, требуя около 1 мин па решение и около 20 мин па предвычислепие системы с двумя контрольными элементами. Время вычисления зависит от формы модели и упругих свойств материала. Полученные решения для поверхностных узлов Nsvorts~5400 пишутся в файл, который в сжатом виде занимает 1Mb. Мы используем графический формат Open Inventor, чтобы иметь возможность контролировать решения визуально (рис.46) и чтобы упростить перенос данных в виртуальные окружения. Визуализация выполняется в реальном времени посредством программного модуля, разработанного с помощью системы виртуального окружения Аванго.
Реализация в Аванго: Взаимодействие с виртуальными
объектами в крупномасштабных системах виртуального окружения выполняется с помощью спе-
циального электромагнитного устройства (трехмерной указки), визуально представленной в виртуальном мире, например, как зеленый луч. Аванго предоставляет удобный интерфейс, называемый MatrixDragger, который па нижнем уровне определяет пересечение луча с объектом и копирует матрицу устройства взаимодействия в матрицу объекта, что позволяет с легкостью перемещать виртуальные объекты. В нашем приложении MatrixDragger прикреплены к контрольным элементам и их матрицы присоединены к нашему программному модулю посредством полевых связей. Всякий раз, когда матрицы изменяются, активируется метод fieldHas Changed, который извлекает коэффициенты с/. из элементов матриц. Вычисление суммы Ylk ckU^ откладывается до вызова метода evaluate, который активируется только один раз па кадр, если какое-либо поле изменилось. Различные методы могут использоваться для ускорения этого вычисления. В частности, если используется только одно устройство взаимодействия, тогда только один контрольный элемент может переместиться за один акт взаимодействия, при этом изменяются только соответствующие ему коэффициенты Ск, и только соответствующая часть суммы должна быть перевычислепа, см рис.4в. Дополнительные возможности для ускорения предоставляют использование арифметики указателей, помещение ДАННЫХ в область памяти быстрого доступа (cache) и параллелизация вычислений. Описанная модель интерактивно деформируется в виртуальном окружении при скорости графики 20 кадров в секунду, при этом все вычисления, необходимые для деформации, требуют около 30% времени одного кадра, для компьютеров типа SGI/Onyx2 300MHz MIPS R12000. Характеристики исполнения па Linux PC Athlon 1.3MHz лучше: 85 кадров в секунду.
Заключение
В данной работе мы описали подход, применимый для моделирования упругих деформаций объектов в реальном времени. Этот подход использует метод конечных элементов для представления уравнений теории упругости. Большая часть интенсивных вычислений выполняется заранее. Полученные данные, сохраненные в файле, позволяют ускорить в тысячи раз оставшуюся часть вычислений, которые выполняются в реальном времени. Это дает возможность моделировать упругие деформации обьектов высокой сложности в виртуальных окружениях.
Список литературы
|1| М. Bro-Niolson, S. Cotin, Real-Time Volumetrie Deformable Models for Surgery Simulation using Finite Elements and Condensation, Eurographics'96, Vol.15 (1996) No.3, C57-C66.
|2| Л.Д.Ландау, Е.М.Лифитиц, Теория упругости, Теоретическая физика, Т.7, Москва, Наука, 1970.
|3| Gibson S.F., Mirtich В., A survey of deformable models in computer graphics, Technical Report TR-97-18, Mitsubishi Electric Research Laboratories, Cambridge, MA, November 1997.
|4| James D., Pai D., Artdefo Accurate Real Time Deformable Objects, Computer Graphics, vol.33, pp.65-72, 1999.
|5| Picinbono G., Delingette H., Nicholas H.-A., NonLinear Anisotropic Elasticity for Real-Time Surgery Simulation, INRIA Research Report 4028, October 2000.
|6| Y.Zhuang, J.F.Canny, Real-time global deformations, In Algorithmic and Computational Robotics, pp.97-107, Natick MA, 2001. A.K.Peters.
|7| Krüger W., Bohn C., Fröhlich В., Schiith H., Strauss W., Wesche G., The Responsive Workbench: A Virtual Work Environment, IEEE Computer, pp. 1215, 1994.
|8| Eckel G., Göbel M., Hasenbrink F., Heiden W., Lechner U., Tramberend H., Wesche G., Wind J. Benches and Caves. In: Bullinger H..J., Riedel О. (eds.) Proc. 1st Int. Immersive Projection Technology Workshop. Springer-Verlag, London, 1997.
|9| Cruz-Neira C. Surround-Screen Projection-Based Virtual Reality: The Design and Implementation of the CAVE. Computer Graphics Proc., Annual Conference Series, 1993, pp.135-142.
|10| Tramberend H., Avocado: A Distributed Virtual Reality Framework, Proc. of the IEEE Virtual Reality. 1999.
1111 R. Kent Dybvig. The Scheme programming language: ANSI Scheme. P T R Prentice-Hall, Englewood Cliffs, N.J 07632, USA, Edition 2, 1996.
1121 В.С.Владимиров, Уравнения математической физики, Издание 5, Москва, Наука, 1988.
|13| О. С. Zienkiewcz, R. L. Taylor, The Finite Element Method, Vol.1, Edition 5, Butterworth-Heinemann, Oxford 2000.
|14| A. O. Frank, I. A. Twombly, J. D. Smith, Finite Element Methods for real-time Haptic Feedback of Soft-Tissue Models in Virtual Reality Simulators, Proceedings of VR2001.
|15| U. Meier et al, Real-Time Simulation of Minimally-Invasive Surgery with Cutting Based on Boundary Element Methods, Lecture Notes in Computer Science (2001) V.2208, p.1263.
|16| Hirota, G., Fisher, S. and Lin, M., Simulation of Non-penetrating Elastic Bodies Using Distance Fields. UNC Technical Report TR00-018, 2000,
http://www.cs.unc.edu/~sfisher/Research/research.html
Приложение: За пределами теории малых деформаций
линейндя упругпсть
нелинейння упругость
Рис.6. Различные типы упругих объектов.
Имеются два типа нелинейностей, которые могут усложнить моделирование: материальные и геометрические. Материальные нелинейности, вызванные нарушением линейного соотношения между тензорами напряжения и деформации (закон Гука), обычно, по крайней мере в инженерных приложениях, возникают только в экстремальных условиях сильных напряжений, когда материал начинает терять свои упругие свойства. Геометрические нелинейности проявляются при больших смещениях точек тела друг относительно друга и связаны с квадратичным членом в тензоре деформаций, см. раздел 3. Для тел, размеры которых примерно одинаковы во всех трех измерениях и которые мы будем называть ЗОЕ-телами, см. рис.6, геометрические нелинейности обычно пренебрежимо малы, поскольку большие относительные смещения для этого типа тел могут возникать только под действием очень больших напряжений. В противоположность атому, при деформациях длинных стержней (IDE-тел) и цилиндрическом изгиб« тонких пластин (20Е-тел) большие смещения возможны без создания большого напряжения. Существенная часть теории упругости посвящена чтим специальным случаям.
Большинство разработанных в настоящее время численных методов, пригодных для моделирования деформаций упругих объектов в реальном времени, основаны на линейной теории малых деформаций. Итеративные методы обладают худшей сходимостью в случае IDE или 2DE тел, поскольку для таких тел существуют направления, вдоль которых энергия изменяется намного медленнее, чем вдоль других. Именно эти направления, отвечающие большим смещениям без создания большого напряжения, являются предпочтительными при деформации тел, в то время как шаг итерации, определяемый другими направлениями, оказывается малым, что замедляет сходимость метода. Также, для IDE или 2DE объектов решения линейной теории часто проявляют нефизические свойства. Например, если один конец длинного стержня зафиксирован, в то время как другой вращается вокруг своей оси, выполняя один полный оборот, граничное условие Ь возвращается к недеформированному состоянию, при этом решение теории малых деформаций и = —А~^ИЬ также непрерывно возвращается в недеформированное состояние,
3DE
2DE
сложным образом проникая сквозь себя. Эта проблема вызвана опущенным квадратичным членом в определении тензора деформаций, который не является малым для данных деформаций. Этот член также отвечает за полную вращательную инвариантность в теории, так что энергия линейной теории оказывается инвариантной только относительно инфинитезимальных, но не относительно полных вращений. Неправильное поведение линейной системы по отношению ко вращениям объекта является индикатором проблем, связанных с нелинейными эффектами
И-
Ниже мы приводим несколько подходов, которые способны учесть геометрически нелинейные эффекты для моделирования упругих деформаций. Мы подразделяем подходы на те, которые должны целиком использоваться в реальном времени, те, которые основаны на предварительных вычислениях, и смешанные.
Методы, работающие в реальном времени:
1а) использование тензора деформаций с квадратичным членом, как было сделано в [5];
16) введение системы локальных координат [2], в которых квадратичные члены в тензоре деформаций малы и ими можно пренебречь.
Методы, использующие предварительные вычисления (применимые для не очень больших углов вращения, например, при более точном моделировании ЗВЕ-тел):
2а) в схеме, использующей контрольные элементы, выделение средней матрицы вращения всех контрольных элементов и перенесение ее на вращение модели как целого. Это уменьшит углы вращения контрольных элементов по отношению к среднему положению, значительно уменьшая нелинейные члены (пропорциональные квадратам углов).
26) решение нелинейной системы, используя теорию возмущений по нелинейному члену и предвычисляя достаточное количество коэффициентов разложения.
Смешанные методы: За) подразбиение IDE и 2DE тел на множество 3DE частей, для которых применима теория малых деформаций и поведение которых может быть описано в терминах предвычисленных функций Грина, состыковка частей вдоль границ, выражение энергии системы в терминах меньшего набора граничных координат и ее минимизация в реальном времени, используя подходящий нелинейный минимизатор [16].
36) рассматривая нелинейное поведение длинных стержней и тонких пластин, можно использовать методы, основанные на частично аналитическом разрешении дифференциальных уравнений в частных производных, проведенном в нелинейной теории упругости [2].
