DelFEM を使ってみる

2010年5月25日

はじめに

DelFEM を使ってみる。

注意 以下、大部分を想像で書いております。正しくない部分があるかも。

環境

DelFEM 1.2.0, MinGW。

ファイル

プログラムの使い方

ここで作成したプログラムは、マウスで棒をゆらゆら揺らすだけのものである。操作方法はウインドウに書いてあるが、以下の通りである。

  • 左ボタンドラッグ: 回転
  • 中央ボタンドラッグ: 移動
  • 右ボタンドラッグ: 拡大・縮小
  • Shift + 左ボタンドラッグ: 棒を揺らす
  • p: 透視射影と正射影の切り替え
  • w: ワイヤーフレーム表示の切り替え
  • h: ヘルプの表示の切り替え
  • q: 終了

Shift + 左ボタンドラッグは、上下が Y 方向正負、左右が Z 方向正負に棒が動くようになっているので、棒の先端を手前に見ながらやると、マウスと棒の動く方向が一致して揺らしやすい。

プログラムの内容

DelFEM を使うためには、ベースを OpenGL で作る必要があるので、それは自分でなんとかする。DelFEM は C++ のライブラリなので、プログラムは C++ で記述する。ここでは、サンプルの solid3d を参考にする。DelFEM にはカメラ管理クラスが存在するが、ベースに既存のプログラムを流用したため、ここではカメラを自前で管理している。

DelFEM に関しては、以下のファイルをインクルードしている。

#include <delfem/cad_obj2d.h>
#include <delfem/mesh3d.h>
#include <delfem/field.h>
#include <delfem/field_world.h>
#include <delfem/drawer_field.h>
#include <delfem/drawer_field_face.h>
#include <delfem/drawer_field_edge.h>
#include <delfem/eqnsys_solid.h>

初期化は init() で行っている。上のほうの数行は OpenGL の操作についての準備で、以下が DelFEM に関するコードである。

	/* create 2D mesh */
	Cad::CCadObj2D cadObj2D;
	vector<Com::CVector2D> v;
	v.push_back(Com::CVector2D(-0.5, 0.));
	v.push_back(Com::CVector2D(0.5, 0.));
	v.push_back(Com::CVector2D(0.5, 0.1));
	v.push_back(Com::CVector2D(-0.5, 0.1));
	cadObj2D.AddPolygon(v);
	Msh::CMesher2D mesh2D(cadObj2D, 0.05);

	/* extrude 2D mesh */
	Msh::CMesh3D_Extrude mesh3D;
	mesh3D.Extrude(mesh2D, 0.1, 0.05);

	/* add mesh to world */
	world.Clear();
	int meshID = world.AddMesh(mesh3D);

	/* set solid */
	solid.SetDomain_Field(meshID, world);
	solid.SetYoungPoisson(YOUNG_MODULUS, POISSON_RATIO);
	solid.SetTimeIntegrationParameter(DTIME, 0.8);

	solid.UnSetStationary();
	solid.UnSetGeometricalNonLinear();
	solid.UnSetSaveStiffMat();

	solid.SetGravitation(0., -1., 0.);

	/* B. C. */
	const Fem::Field::CIDConvEAMshCad &conv = world.GetIDConverter(meshID);
	bcID = solid.AddFixElemAry(conv.GetIdEA_fromCad(4, Cad::EDGE, 2), world);

	/* drawer */
	setDrawer(0);

処理をひとつずつ見ていこう。まず、棒のメッシュを作る。

	/* create 2D mesh */
	Cad::CCadObj2D cadObj2D;
	vector<Com::CVector2D> v;
	v.push_back(Com::CVector2D(-0.5, 0.));
	v.push_back(Com::CVector2D(0.5, 0.));
	v.push_back(Com::CVector2D(0.5, 0.1));
	v.push_back(Com::CVector2D(-0.5, 0.1));
	cadObj2D.AddPolygon(v);
	Msh::CMesher2D mesh2D(cadObj2D, 0.05);

2 次元の長方形を作り、メッシュを作成している。棒の先端の面ではなくて、側面を作っている。CMesher2D のコンストラクタの第 2 引数は、メッシュサイズである。

	/* extrude 2D mesh */
	Msh::CMesh3D_Extrude mesh3D;
	mesh3D.Extrude(mesh2D, 0.1, 0.05);

2 次元のメッシュを伸ばして、3 次元の棒のメッシュを作る。Extrude の第 2 引数は伸ばす長さ、第 3 引数はメッシュサイズ。

	/* add mesh to world */
	world.Clear();
	int meshID = world.AddMesh(mesh3D);

world にメッシュを追加する。ここで world は Fem::Field::CFieldWorld のオブジェクトである。

	/* set solid */
	solid.SetDomain_Field(meshID, world);
	solid.SetYoungPoisson(YOUNG_MODULUS, POISSON_RATIO);
	solid.SetTimeIntegrationParameter(DTIME, 0.8);

方程式 (Fem::Eqn::CEqn_Solid3D_Linear のオブジェクト) に関する設定。領域の設定、物性の設定、時間積分の設定をしている。最後の行の DTIME は計算の時間刻み、0.8 というのは Newmark-β法のγの値。時間刻みだけ設定してもよいのだが、減衰を入れるために適当にγをいじる。

	solid.UnSetStationary();
	solid.UnSetGeometricalNonLinear();
	solid.UnSetSaveStiffMat();

SetStationary() で静的計算、UnSetStationary() で動く。SetGeometricalNonLinear() にすると幾何学的非線形が考慮され、要素が膨張したりしなくなるが、計算が重いので、ここでは UnSetGeometricalNonLinear() とする。SetSaveStiffMat() とすると、剛性行列を保存する、らしい。

	solid.SetGravitation(0., -1., 0.);

重力。Y 下向きに適当な値を入れる。

	/* B. C. */
	const Fem::Field::CIDConvEAMshCad &conv = world.GetIDConverter(meshID);
	bcID = solid.AddFixElemAry(conv.GetIdEA_fromCad(4, Cad::EDGE, 2), world);

拘束条件の追加。棒の左端を固定している。

	/* drawer */
	setDrawer(0);

描画に関する Fem::Field::View::CDrawerArrayField のオブジェクトの設定をする。setDrawer() の中身は以下のようになっている。

void setDrawer(int wireframe)
{
	drawerArray.Clear();

	int fieldID = solid.GetIdField_Disp();

	if(!wireframe){
		drawerArray.PushBack(
		  new Fem::Field::View::CDrawerFace(fieldID, false, world)
		);
	}

	drawerArray.PushBack(
	  new Fem::Field::View::CDrawerEdge(fieldID, false, world)
	);
}

フィールドのフェイスとエッジを描画するものとして追加している (ワイヤーフレーム表示のときは、エッジのみ追加)。

これで準備は OK。描画に関しては、display() に処理がある。

	glEnable(GL_POLYGON_OFFSET_FILL);
	glPolygonOffset(1., 1.);

	/* solve */
	currentTime += DTIME;
	world.FieldValueExec(currentTime);
	solid.Solve(world);
	drawerArray.Update(world);

	/* draw model */
	glEnable(GL_DEPTH_TEST);
	drawerArray.Draw();
	glDisable(GL_DEPTH_TEST);

はじめの OpenGL の設定は、エッジの表示が変にならないようにするものらしい。次の処理は、時間を時間刻み分進め、方程式を解き、更新。最後にデプステストを ON にして描画している。

最後に、マウスで棒を揺らす処理は、glutMotionFunc() に設定している motion() で行っている。 

	switch(mouseButton){
	case GLUT_LEFT_BUTTON:
		if(modifiers == GLUT_ACTIVE_SHIFT){
			/* swing */
			Fem::Field::CField &bc = world.GetField(bcID);
			positionY += -(y - mouseStartY)*BC_MOVE_FACTOR;
			bc.SetValue(positionY, 1, Fem::Field::VALUE, world, true);
			positionZ += -(x - mouseStartX)*BC_MOVE_FACTOR;
			bc.SetValue(positionZ, 2, Fem::Field::VALUE, world, true);
		}else{
			...

左ボタンでドラッグされていて、Shift キーが押されているとき (modifiers の値は mouse() で glutGetModifiers() から得ている)、棒の左端を変位させる。world.GetField() で固定境界条件のフィールド bc を得て、bc.SetValue() で変位を設定している。SetValue() の第 1 引数が変位 (というか "現在位置" というべきか)、第 2 引数は値を設定する方向で、1 が Y 方向、2 が Z 方向である。

動作速度の調整

タイマーを使っていないので、動作の速さは環境に依存する。速すぎたり遅すぎたりする場合は、DTIME (時間刻み幅) で調整できる。

DOS 窓を開かないようにする

MinGW でふつうにコンパイルすると、プログラムの実行時に DOS 窓が開く。開かないようにしたい場合は、Makefile の LDFLAGS に次のように書く。

LDFLAGS = -Wl,--subsystem,windows

-Wl は、リンカにオプションを渡すための gcc のオプション。