VegaFEM免费试用

最新推荐文章于 2024-05-08 09:48:02 发布

iteye_10018

最新推荐文章于 2024-05-08 09:48:02 发布

阅读量282

点赞数

文章标签： python

何为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的年轻教授，他主页提供的图形学的相关课程也非常值得好好学习。

官网：http://run.usc.edu/vega/

编译

环境：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 777generateVeg.sh

sudo chmod 777VegGenerator.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

运行结果：

iteye_10018

关注

0
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
VegaFEM免费试用

何为FEM 还记得ipod touch5发布的那个忍者游戏么：http://v.youku.com/v_show/id_XNDU5Njc3ODQw.html 这个人物身体动作的模拟用的就是FEM。 FEM是有限元分析（Finite element method）的缩写，最早是应用在工程领域，今年随着计算机性能的不断发展，FEM在图形和游戏领域都有着很好的发展。 ...
复制链接

扫一扫