DelFEM を使ってみる2010年5月25日 ![]() | |
はじめにDelFEM を使ってみる。 注意 以下、大部分を想像で書いております。正しくない部分があるかも。 環境DelFEM 1.2.0, MinGW。 ファイル
プログラムの使い方ここで作成したプログラムは、マウスで棒をゆらゆら揺らすだけのものである。操作方法はウインドウに書いてあるが、以下の通りである。
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 のオプション。 | |
PENGUINITIS |