何为FEM
还记得ipod touch5发布的那个忍者游戏么:http://v.youku.com/v_show/id_XNDU5Njc3ODQw.html
这个人物身体动作的模拟用的就是FEM。
FEM是有限元分析(Finite element method)的缩写,最早是应用在工程领域,今年随着计算机性能的不断发展,FEM在图形和游戏领域都有着很好的发展。
整个物理模拟的过程都围绕着牛顿第二定律进行:
也就是 F = ma。
何为VegaFEM
VegaFEM是一个用于模拟三维变形物体的高效稳定的C/ C + +物理库。它被可用于计算大变形,包括几何和材料非线性模型,也可以有效地模拟线性系统。Vega是开放源代码和自由的。它的3-clause BSD许可证下发布的,这意味着它可以用来无论是在学术研究和商业应用。
Vega的作者是Jernej Barbič,一个MIT的年轻教授,他主页提供的图形学的相关课程也非常值得好好学习。
编译
环境:Ubuntu 12.04
VegaFEM最近发布了2.0版本,新版本增加了布料模拟,model reduction等新特性,建议直接下在最新版本。
下载好源码之后,解压,终端cd进入目录,直接make,过1分钟左右就编译好了。
源码中包含了一些例子,可以用刚才编译好的 interactiveDeformableSimulator 直接模拟。这个Simulator能够从配置文件中读取初始化参数,加载模型,加载单元并计算,最终用GLUT显示和交互。
我们可以来看看 sample 中的配置文件,了解每项的含义。
beam3_tet.config
*volumetricMeshFilename
../../models/beam3/beam3_tet.veg
*solver
#implicitBackwardEuler
implicitNewmark
#centralDifferences
#symplecticEuler
#Euler
*deformableObjectMethod
#StVK
InvertibleFEM
#CLFEM
#LinearFEM
*invertibleMaterial
StVK
#neoHookean
*renderingMeshFilename
../../models/beam3/beam3_tet.obj
*timestep
0.0005
*dampingMassCoef
0.0
#1.0
#10.0
*dampingStiffnessCoef
0.01
#0.001
*deformableObjectCompliance
1.0
*baseFrequency
1.0
*lightingConfigFilename
beam3_tet.lighting
*massMatrixFilename
../../models/beam3/beam3_tet.mass
*cameraRadius
2
*cameraLongitude
65.8
*cameraLattitude
29
*focusPositionX
0.0
*focusPositionY
0.5
*focusPositionZ
0.0
*syncTimestepWithGraphics
0
*fixedVerticesFilename
../../models/beam3/beam3.bou
*numInternalForceThreads
4
*inversionThreshold
0.1
一项项解释,
1.分割好的单元模型文件;
2.求解器
3.力学模型
4.可逆材料模型
5.用于放在opengl渲染的mesh模型
6.时间步长
7.质量阻尼系数
8.刚度阻尼系数
7.和8都是用来计算瑞丽阻尼。
9.变形对象合规性( 不懂)
10.基础频率
11.灯光参数
后面几个是Opengl的参数。
fixedVerticesFilename :固定的顶点。
使用Vega
如果想用Vega写一些自己的程序的话就不能用自带的sumulator 了,局限性实在太大。
下面按照文档中给出的教程一步步来写点代码。
写代码之前首先要生成一个单元网格,要生成单元网格当然需要一个mesh,要创建mesh当然要用blender了。
注:网格剖分的知识请参考:专注网格剖分
打开blinder直接添加一个圆柱体好了~
File -> Export 导出一个 pillar.stl 模型,注意勾选ASCII选项,不然tetgen不认。
接下来要用 tetgen 来生成四面体网格,由于tetgen生成的文件并不是Vega所需要的 *.veg 文件,所有需要写个脚本来转换一下。
在tetgen目录下面
创建一个 generateVeg.sh
./tetgen $1
python VegGenerator.py $1
再创建一个 VegGenerator.py
import sys
import string
tmpSplit = sys.argv[1].split('.')
modelName = tmpSplit[0]
print "Generator " + tmpSplit[0] +".veg ..."
node_file = open(tmpSplit[0] + ".1.node")
ele_file = open(tmpSplit[0] + ".1.ele")
output_file = open(tmpSplit[0] +".veg","w");
node_lines = node_file.readlines()
ele_lines = ele_file.readlines()
output_file.write("*VERTICES\n")
for line in node_lines:
if line.startswith('*') or line.startswith('#') or line.startswith('E'):
output_file.write(line);
continue
line = line.strip()
split = line.split()
output_file.write(split[0])
output_file.write(' ')
output_file.write(split[1])
output_file.write(' ')
output_file.write(split[2])
output_file.write(' ')
output_file.write(split[3])
#print '%s strip=%s' % (line,line.split())
output_file.write('\n')
print line
node_file.close();
output_file.write("*ELEMENTS\n")
output_file.write("TET\n")
for line in ele_lines:
if line.startswith('*') or line.startswith('#') or line.startswith('E'):
output_file.write(line);
continue
line = line.strip()
split = line.split()
for i in split:
output_file.write(i)
output_file.write(' ')
output_file.write('\n')
print line
ele_file.close();
output_file.write("*MATERIAL BODY\nENU, 1000, 10000000, 0.45\n\n*REGION\nallElements, BODY")
output_file.close();
sudo chmod 777 generateVeg.sh
sudo chmod 777 VegGenerator.py
将刚才制作好的模型 pillar.stl 拷贝到tetgen目录,终端执行:
./generateVeg.sh pillar.stl
接下来,当前 目录下就会生成相应的文件:
pillar.veg就是我们需要的。
下面就需要写代码了。
这个工程创建会有一点麻烦。
创建一个文件夹,里面创建一些文件和文件夹:
将vega 的所有头文件拷贝到 inc 文件夹中, 所有 *.a文件拷贝到lib中,main.cpp 内容如下:
#include "volumetricMeshLoader.h"
#include "corotationalLinearFEM.h"
#include "corotationalLinearFEMForceModel.h"
#include "generateMassMatrix.h"
#include "implicitBackwardEulerSparse.h"
#include <iostream>
using namespace std;
int main()
{
char inputFilename[96] = "pillar.veg";
VolumetricMesh * volumetricMesh = VolumetricMeshLoader::load(inputFilename);
if (volumetricMesh == NULL)
cout<<"Error: failed to load mesh."<<endl;
else
{
cout<<"Success. Number of vertices:"<<volumetricMesh->getNumVertices()<<endl;
cout<<" Number of Elements:"<<volumetricMesh->getNumElements()<<endl;
}
//printf("Success. Number of vertices: %d . Number of elements: %d .\n",
//volumetricMesh->getNumVertices(), volumetricMesh->getNumElements());
TetMesh *tetMesh;
if(volumetricMesh->getElementType()==VolumetricMesh::TET)
{
tetMesh = (TetMesh*) volumetricMesh;
cout<<"A tet mesh."<<endl;
}
else
{
cout<<"Error,not a tet maehs."<<endl;
exit(1);
}
CorotationalLinearFEM * deformableModel = new CorotationalLinearFEM(tetMesh);
cout<<"good"<<endl;
// create the class to connect the deformable model to the integrator
ForceModel * forceModel = new CorotationalLinearFEMForceModel(deformableModel);
int r = 3 * tetMesh->getNumVertices(); // total number of DOFs
double timestep = 0.0333; // the timestep
SparseMatrix * massMatrix;
GenerateMassMatrix::computeMassMatrix(tetMesh, &massMatrix, true);
// This option only affects PARDISO and SPOOLES solvers, where it is best
// to keep it at 0, which implies a symmetric, non-PD solve.
// With CG, this option is ignored.
int positiveDefiniteSolver = 0;
// constraining vertices 4, 10, 14 (constrained DOFs are specified 0-indexed):
int numConstrainedDOFs = 9;
int constrainedDOFs[9] = { 12, 13, 14, 30, 31, 32, 42, 43, 44 };
// (tangential) Rayleigh damping
double dampingMassCoef = 0.0; // "underwater"-like damping (here turned off)
double dampingStiffnessCoef = 0.01; // (primarily) high-frequency damping
// initialize the integrator
ImplicitBackwardEulerSparse * implicitBackwardEulerSparse = new
ImplicitBackwardEulerSparse(r, timestep, massMatrix, forceModel,
positiveDefiniteSolver, numConstrainedDOFs, constrainedDOFs,
dampingMassCoef, dampingStiffnessCoef);
// alocate buffer for external forces
double * f = (double*) malloc (sizeof(double) * r);
int numTimesteps = 10;
for(int i=0; i<numTimesteps; i++)
{
// important: must always clear forces, as they remain in effect unless changed
implicitBackwardEulerSparse->SetExternalForcesToZero();
if (i == 0) // set some force at the first timestep
{
for(int j=0; j<r; j++)
f[j] = 0; // clear to 0
f[37] = -500; // apply force of -500 N to vertex 12, in y-direction, 3*12+1 = 37
implicitBackwardEulerSparse->SetExternalForces(f);
}
implicitBackwardEulerSparse->DoTimestep();
}
// alocate buffer to read the resulting displacements
double * u = (double*) malloc (sizeof(double) * r);
implicitBackwardEulerSparse->GetqState(u);
for(int i=0;i<10;i++)
cout<<*(u+i)<<endl;
return 1;
}
代码就是加载.veg文件,然后进行模拟。几个类说明一下。
InetgratorBase:存储单元位移,速度,加速度,提供了修改和计算的方法,包含了瑞丽阻尼系数。
IntegratorBaseSparse: 继承自IntegratorBase,当M,D,K为稀疏矩阵的时候使用。
implicitBackwardEulerSparse:继承自IntegratorBaseSparse,实现了隐式欧拉后向计算。
ForceMode:获取内力和刚度矩阵。
更多的说明请参考文档。
接下来写CMakeLists.txt
PROJECT(MAIN)
CMAKE_MINIMUM_REQUIRED(VERSION 2.8)
set(SRC_LIST main.cpp)
INCLUDE_DIRECTORIES(${MAIN_SOURCE_DIR}/inc)
LINK_DIRECTORIES(${MAIN_SOURCE_DIR}/lib)
ADD_EXECUTABLE(main ${SRC_LIST})
TARGET_LINK_LIBRARIES(main sceneObject integrator elasticForceModel forceModel matrix sparseMatrix loadList insertRows lighting configFile volumetricMesh getopts camera graph isotropicHyperelasticFEM minivector stvk corotationalLinearFEM polarDecomposition massSpringSystem objMesh sparseSolver)
SET(SRC_LIST main.cpp)
终端进入build文件夹
cmake..
make
运行结果: