Физическое моделирование воды.
Автор: Петр С. Попов
Эффект водной глади является желанным гостем в трехмерной графике. В краткой заметке я попытаюсь рассказать про физическое моделирование. Побудительными мотивами являются следующие обстоятельства. Во-первых, довольно абстрактная mIRC лекция Yann L., посвященная моделированию воды. Безусловно очень компетентного и знающего человека. Однако, в лекции допущены грубые неточности. Его реализация (которая, вопреки заявленному, не имеет к уравнению Навье-Стокса ни малейшего отношения) не является оптимальной с точки зрения сеточных методов решения дифференциальных уравнений. В короткий промежуток времени мне на глаза попалось и другое воплощение. Тоже весьма некачественное.
В заметке будет рассказано о самом общеупотребительном сеточном методе решения уравнения колебаний:
(*) d2U / dt2 = d2U / dx2 + d2U / dy2
Здесь t - время, x,y - координаты. Решая это уравнение, мы получим очень реалистичную поверхность воды. Кроме того, в статье будет приведен очень быстрый (и правильный с точки зрения математики) расчет нормалей без вычисления обратного квадратного корня и векторных произведений.
Не надо пугаться страшных непонятных буковок в уравнении (*). Дальнейшее изложение не требует понимания того, что именно там написано.
Итак, некоторые определения: вершина (нормаль и координата), и массив, содержащий числа. Эти самые числа будут интерпретированы как высота поверхности воды в регулярной сетке. Нам надо завести два таких массива, и мы будем переключаться между ними по надобности.
struct vertex { float coo[3]; float nor[3]; }; struct field { float U[128][128]; }; vertex vertices[128][128]; field A,B; field *p=&A,*n=&B;
А вот и основной шаг программы.
void time_step() { int i,j,i1,j1; i1=rand( )%110; j1=rand( )%110; /*1*/ if( ( rand( )&127)==0) for( i=-3;i<4;i++) { for( j=-3;j<4;j++) { float v=6.0f-i*i-j*j; if( v<0.0f)v=0.0f; n->U[i+i1+3][j+j1+3]-=v*0.004f; } } for( i=1;i<127;i++) { for( j=1;j<127;j++) { /*2*/ vertices[i][j].coo[2]=n->U[i][j]; vertices[i][j].nor[0]=n->U[i-1][j]-n->U[i+1][j]; vertices[i][j].nor[1]=n->U[i][j-1]-n->U[i][j+1]; /*3*/ #define vis 0.005f float laplas=( n->U[i-1][j]+ n->U[i+1][j]+ n->U[i][j+1]+ n->U[i][j-1])*0.25f-n->U[i][j]; /*4*/ p->U[i][j]=( ( 2.0f-vis)*n->U[i][j]-p->U[i][j]*( 1.0f-vis)+laplas); } } /*5*/ for( i=1;i<127;i++) { glBegin( GL_TRIANGLE_STRIP); for( j=1;j<127;j++) { glNormal3fv( vertices[i][j].nor); glVertex3fv( vertices[i][j].coo); glNormal3fv( vertices[i+1][j].nor); glVertex3fv( vertices[i+1][j].coo); } glEnd( ); } /*6*/ field *sw=p;p=n;n=sw; }
Сначала /*1*/ мы случайным образом возмущаем водную гладь (c вероятностью 1/127 упадет капля в произвольную точку).
Потом /*2*/ мы обновляем информацию о высоте для OpenGL и считаем нормаль (обратите внимание, как именно мы ее считаем!) Нормализация будет поручена ядру OpenGL.
Потом /*3*/ мы считаем аналог правой части уравнения (*). Так называемый оператор Лапласа. Этакий симметричный крест, причем сумма коэффициентов равна 0. Заметим, что в массиве p мы храним высоту воды на прошлом кадре, а в массиве n - текущую высоту воды.
Потом /*4*/ мы считаем новое значение высоты поверхности воды. Заметим, что мы используем значения на двух временных слоях, значение оператора Лапласа на текущем слое и некоторое число vis. Это - вязкость (вот тут оказывается, что решаем мы не совсем уравнение (*), а чуть другое уравнение). С ней волны затухают. А без нее постоянные падения капель разболтают схему. Замечу также, что все вычисления собраны в один проход. Это хорошо.
Пункт /*5*/ посвящен отрисовке.
В конце /*6*/ мы переключаем временные слои. Только что вычисленные значения становятся "новыми", а предыдущий массив "стареет".
Вот и все. Любопытный читатель может спросить, какое отношение уравнение (*) имеет к данной схеме. Отвечу. Пусть U(x,y) ~ U[ i ][ j ] Частная производная d2 U / dy2 может быть в некотором смысле приближена разностью значений по горизонтали:
d2U / dy2 ~ U[ i ][ j+1 ] - 2* U[ i ][ j ]+ U[ i ][ j-1 ]
Аналогично, частная производная по x может быть приближена "вертикальной" разностью. Аналогично мы можем представить и производную по времени - как разность трех последовательных временных слоев. Мы ищем "самый новый" слой, зная два предыдущих. Шаг /*4*/ получается в некотором роде автоматически. Пока без вязкости. Роль вязкости видна, если положить ее коэффициент равным единице. Тогда формула /*4*/ превращается в чистое усреднение - любое возмущение водной глади будет размываться. А так получается что-то среднее.
Вопросы и Ответы.
1.) Почему для рендеринга OpenGL используется glBegin():glEnd()? - Удобство, краткость. Безусловно, необходимо разделить статические и динамические данные, первые надо поместить в высокоприоритетную память, вторые обновлять таким образом, чтобы согласовать работу CPU и GPU. Использование расширения VBO будет весьма полезным.
2.) В некоторых воплощениях хранят значение на данном слое и "скорость" - разницу между значением на предыдущем слое. Кто прав? - Вопрос исключительно удобства. Можно делать и так и так.
3.) Какой смысл имеет оператор Лапласа? Я видел воплощения, где используется не крест, а шаблон по всем 9 точкам. - Данная величина имеет смысл силы, с которой соседние точки влияют на данную. Крестовой шаблон является с математической точки зрения оптимальным.
4.) Дальнейшие оптимизации возможны? - Да, конечно, хотя основной код неплох. Хороший пример - вычисление текстуры, хранящей водную поверхность, на GPU. На любом ускорителе можно быстро вычислять такую текстуру: http://developer.nvidia.com/docs/IO/1156/ATT/Water.zip . На ускорителях от NVidia (или с помощью грядущих в OpenGL UberBuffers) возможно использовать эту текстуру как массив, содержащий высоту в вершинах. Таким образом, можно перенести алгоритм на GPU. Также возможно считать и физику деформируемых поверхностей: http://developer.nvidia.com/docs/IO/4544/ATT/cg_physics.zip .
5.) В примере волны отражаются от стенок. Можно ли сделать периодическую картину?- Да, заменив циклы от 1 до 127 на циклы от 0 до 128. Точки, выходящие за границу массива, циклически переставляются (-1 переходит 127, 128 переходит в 0). Более того, можно сделать острова в воде, от которых волны отражаются. Для этого надо указать статическую маску, в которой значения водной поверхности не обновляются.
6.) Какова скорость расчета без рендеринга? - Примерно 4000 проходов в секунду для массива 128x128 (P4 2.4, оптимизующий комплятор IntelC) с выполнением пункта /*2*/ или примерно 8000 проходов без оного.
7.) Мне интересна реализация грамотного отражения и преломления в воде. А что мне подсунули! - Возможно, я напишу Cg пример с преломлениями и отражениями. Не надо ожидать многого от реализации менее, чем в 200 строк.
8.) Почему такой большой объем разных дополнительных библиотек? - GLUT является кросспалтформенным стандартом де-факто, без загрузки jpg текстуы (с помощью IJL в данном примере) обойтись сложно, а GLEW содержит все современные OpenGL расширения и поддерживает высокие версии OpenGL (в отличие от стандартных хидеров).
Исходный код: 20030524.rar
24 мая 2003 (Обновление: 15 июня 2009)