提要

植物的运动模拟是图形学中的一个方向，今天就来讨论下怎么模拟出一个在风中荷叶。

植物力学模型

由于植物的模型在风中会有大变形的情况存在，这里采用的是有限元模拟的方法来实现。

关于基于有限元的物理模拟，可以参考：

有限单元法(The Finite Element Method)

专注网格剖分 - TetGen,NETGEN,Steller

在基于有限元物理模拟中，速度是最大的缺点，特别在四面体特别多的情况下，现在的配置比较不错的电脑（i7+GTX680）加上GPU计算的代码，可以接受的情况是小于10000个四面体的情况，但是对于荷花这样比较精细和复杂的模型就很难做到小于10000个的四面体，这时候就需要用到插值技术，就是在精细模型的外面套一层粗糙的四面体，预处理出精细模型到外层四面体的插值矩阵，有限元的模拟用粗糙的四面体来进行，在模拟的同时，将四面体的位置插值回精细的模型，这样就可以在保证运行速度的情况下得到模拟复杂模型的相对精确的结果。下面就是VegaFEM中自带的乌龟模拟的例子：

植物的话，首先要在建模软件中建立植物的模型，和包围体，然后用细分软件对包围体进行四面体划分，之后就可以加载到模拟程序中了。

关于插值矩阵的计算可以参考VegaFEM自带的手册，VegaFEM里面也自带了一个计算矩阵的工具generateInterpolationMatrix，使用的方法是

./generateInterpolationMatrix mesh.veg finemesh.obj inter_mat.txt

生成的插值矩阵就在inter_mat.txt里了。 

风场定义

可以用一个解析式来定义一个风场，

其中x为空间的一个坐标 (x,y,z) , v表示风在 x 处的向量。

这里的风的大小是一个和时间相关的函数，这里就简单地用sin函数来表示。

在配置文件里定义两个变量，一个是double型的strength表示风的强度，还有一个direction变量表示风的方向。

受力计算

风吹在植物模型上，不能够简单地将力离散到每个节点上，需要根据Tet的位置来决定顶点力的大小，按照常识来看，迎着风的三角面肯定受力，背着风的三角面不受力。具体说来就是植物表面三角面的法向和风向的夹角决定了顶点受力的大小。

在计算风力的时候，需要计算风的方向向量和三角面法向的点乘，作为一个计算的系数。最终三角面上风的大小就是：

F = windStrength * sin(simulationTime) * sin(theta) * windDirection

分配到顶点上的力就是1/3 * F.当然，这也是近似计算，因为三角面的大小不同，实际中受力应该是不同的。

整个模拟的流程相关代码如下

double forceStrength = m_configs.value(QStringLiteral("strength")).toDouble();
auto directionList = m_configs.value(QStringLiteral("direction")).split(" ");
Vector3d windDirection(directionList[0].toDouble(), directionList[1].toDouble(), directionList[2].toDouble());
windDirection.Normalize();

//Calculate forces on every vertices
double sinTime = sin(6 * m_simuTime) + 1;
double sinTime2 = sin(2 * m_simuTime) + 1;
double sinTime3 = sin(m_simuTime) + 1;

m_integrator->ClearExternalForce();
for(int i = 0; i < surfaceNormal.size(); i++)
	{
		Vector3d tmpNormal(surfaceNormal[i].x(), surfaceNormal[i].y(), surfaceNormal[i].z());
		tmpNormal.Normalize();
		double dotResult = tmpNormal.Dot(windDirection);
		if(dotResult < 0)
		{
			//Add windforce in face to vertices
			Vector3d windDrag = dotResult * forceStrength * sinTime * windDirection;
			Vector3d windDrag2 = dotResult * forceStrength * sinTime2 * windDirection;
			Vector3d windDrag3 = dotResult * forceStrength * sinTime3 * windDirection;
			Vector3d averageForce = (windDrag + windDrag2 + windDrag3) * 0.33;
			externalForce[surfaceFace[i]] += averageForce;
			externalForce[surfaceFace[i-1]] += averageForce;
			externalForce[surfaceFace[i-2]] += averageForce;
		}
	}		
m_integrator->AddExternalForce(externalForce);
	
m_integrator->BackwardEulerStep(m_timeStep);
m_flagDumpingImg = true;
m_tetMesh->GetVertices() = m_integrator->GetNodalPosition();
DumpNode(m_tetMesh->GetVertices(),frame_index++);

上面计算了三个sinTime，是为了增加风的随机性，让植物看上去更自然一些。

运行效果如下

参考

GPU-Generated Procedural Wind Animations for Trees - http://http.developer.nvidia.com/GPUGems3/gpugems3_ch06.html

Physically Based Real-Time Translucency for Leaves