1. はじめに
.NET、もっと限定すればC#で数値解析できるライブラリが欲しくてIvyFEMを作っています。
- 図面作成
- 図面ファイルの読み書き
- 計算
- 結果表示
の一連の流れをこのライブラリを使って実装することができます。
基本機能は次の通りです。
☑ 単純な2D(ポリゴン)の図面作成
☑ 有限要素（三角形要素）分割
☑ 有限要素行列の作成 （*1）
☑ リニアシステムを解く（LAPACKE, Lis、独自実装）（*1）
☑ サーモグラフィーのような分布図
*1 いま用意しているのは
電磁気学：H面導波管の伝達問題
力学： 弾性体の構造解析
線形弾性体
超弾性体
Saint Venant Kirchhoff
Mooney-Rivlin (非圧縮、微圧縮)
Ogden (非圧縮、微圧縮)
多点拘束(Multipoint Constraint, MPC)(直線)
剛体との接触(直線、円)
2. セットアップ
(1) IvyFEMを次のページからzipをダウンロードします
IvyFEM
https://github.com/ryujimiya/IvyFEM
(2) Visual Studio 2017で、WPFアプリケーションを作成します。プラットフォームターゲットはx64に設定します。
サンプルアプリケーション IvyFEMProtoApp
https://github.com/ryujimiya/IvyFEMProtoApp/
(3) (1)で取得したIvyFEM.dllと関連ライブラリ一式をプロジェクトに追加し、IvyFEM.dllを参照に追加します
プロジェクトに追加したライブラリ一式は、出力ディレクトリに常にコピーするように設定します。
(4) OpenTK.GLControlをNuGetでインストールします
3. MainWindow
3.1. MainWindow.xaml
IvyFEMライブラリで描画する対象、OpenTK.GLControl(名前glControl)を追加します。
MainWindow.xaml
<Grid>
<WindowsFormsHost>
<OpenTK:GLControl x:Name="glControl" Load="glControl_Load" Resize="glControl_Resize" Paint="glControl_Paint" MouseClick="glControl_MouseClick" KeyDown="glControl_KeyDown" KeyUp="glControl_KeyUp" MouseMove="glControl_MouseMove" MouseDown="glControl_MouseDown" MouseWheel="glControl_MouseWheel"/>
</WindowsFormsHost>
</Grid>
GLControlが無いと警告がでますが、Visual Studioの自動修正機能を使って修正できます。修正すると、
MainWindow.xaml
xmlns:OpenTK="clr-namespace:OpenTK;assembly=OpenTK.GLControl"
が追加されると思います。
3.2. MainWindow.xaml.cs
glControlのイベントハンドラなどを実装します。
MainWindow.xaml.cs
/// <summary>
/// MainWindow.xaml の相互作用ロジック
/// </summary>
public partial class MainWindow : Window
{
/// <summary>
/// 問題
/// </summary>
private Problem Problem = null;
/// <summary>
/// カメラ
/// </summary>
internal Camera2D Camera { get; private set; } = new Camera2D();
/// <summary>
/// 場を描画する?
/// </summary>
internal bool IsFieldDraw { get; set; } = false;
/// <summary>
/// 描画アレイ
/// </summary>
internal DrawerArray DrawerArray { get; private set; } = new DrawerArray();
internal FieldDrawerArray FieldDrawerArray { get; private set; } = new FieldDrawerArray();
/// <summary>
/// ウィンドウがロードされた
/// </summary>
/// <param name="sender"></param>
/// <param name="e"></param>
private void Window_Loaded(object sender, RoutedEventArgs e)
{
Problem = new Problem();
}
/// <summary>
/// ウィンドウがアクティブになった
/// </summary>
/// <param name="sender"></param>
/// <param name="e"></param>
private void Window_Activated(object sender, EventArgs e)
{
glControl.MakeCurrent();
}
/// <summary>
/// glControlの起動時に実行される。
/// </summary>
/// <param name="sender"></param>
/// <param name="e"></param>
private void glControl_Load(object sender, EventArgs e)
{
GL.Enable(EnableCap.DepthTest);
}
/// <summary>
/// glControlのサイズ変更時に実行される。
/// </summary>
/// <param name="sender"></param>
/// <param name="e"></param>
private void glControl_Resize(object sender, EventArgs e)
{
glControl_ResizeProc();
}
internal void glControl_ResizeProc()
{
int width = glControl.Size.Width;
int height = glControl.Size.Height;
Camera.WindowAspect = ((double)width / height);
GL.Viewport(0, 0, width, height);
GL.MatrixMode(MatrixMode.Projection);
GL.LoadIdentity();
OpenGLUtils.SetProjectionTransform(Camera);
}
/// <summary>
/// glControlの描画時に実行される。
/// </summary>
/// <param name="sender"></param>
/// <param name="e"></param>
private void glControl_Paint(object sender, System.Windows.Forms.PaintEventArgs e)
{
GL.ClearColor(Color4.White);
GL.Clear(ClearBufferMask.ColorBufferBit | ClearBufferMask.DepthBufferBit);
GL.Enable(EnableCap.PolygonOffsetFill);
GL.PolygonOffset(1.1f, 4.0f);
GL.MatrixMode(MatrixMode.Modelview);
GL.LoadIdentity();
OpenGLUtils.SetModelViewTransform(Camera);
ConstraintDrawerArray.Draw();
if (IsFieldDraw)
{
FieldDrawerArray.Draw();
}
else
{
DrawerArray.Draw();
}
glControl.SwapBuffers();
}
}
DrawerArrayとFieldDrawerArrayがIvyFEMで実装されている描画ライブラリの実体で、DrawerArrayは図面(CAD)描画用、FieldDrawerArrayは計算結果描画用となります。
4. CAD図面の描画
弾性体の曲げの実装の前に、CADの実装例を示します。
CadObject2D（名前cad2D)を生成し、頂点を追加していきます。
次に生成したcad2Dに対応するCadObject2DDrawer(名前drawer)を生成し、最後にDrawerArrayに追加すれば完了です。
あとはMainWindowで実装したglControl_Paint()イベントが発生すればDrawerArray.Draw()が実行されて図面が描画されます
Problem.Cad.cs
public void MakeBluePrint(MainWindow mainWindow)
{
double WaveguideWidth = 1.0;
double InputWGLength = 1.0 * WaveguideWidth;
CadObject2D cad2D = new CadObject2D();
{
IList<OpenTK.Vector2d> pts = new List<OpenTK.Vector2d>();
pts.Add(new OpenTK.Vector2d(0.0, WaveguideWidth)); // 頂点1
pts.Add(new OpenTK.Vector2d(0.0, 0.0)); // 頂点2
pts.Add(new OpenTK.Vector2d(InputWGLength, 0.0)); // 頂点3
pts.Add(new OpenTK.Vector2d(InputWGLength, (-InputWGLength))); // 頂点4
pts.Add(new OpenTK.Vector2d((InputWGLength + WaveguideWidth), (-InputWGLength))); // 頂点5
pts.Add(new OpenTK.Vector2d((InputWGLength + WaveguideWidth), WaveguideWidth)); // 頂点6
var res = cad2D.AddPolygon(pts);
var resCircle = cad2D.AddCircle(new OpenTK.Vector2d(
InputWGLength + 0.5 * WaveguideWidth, 0.5 * WaveguideWidth),
0.25 * WaveguideWidth,
res.AddLId);
}
mainWindow.IsFieldDraw = false;
var drawerArray = mainWindow.DrawerArray;
drawerArray.Clear();
IDrawer drawer = new CadObject2DDrawer(cad2D);
mainWindow.DrawerArray.Add(drawer);
mainWindow.Camera.Fit(drawerArray.GetBoundingBox(mainWindow.Camera.RotMatrix33()));
mainWindow.glControl_ResizeProc();
mainWindow.glControl.Invalidate();
mainWindow.glControl.Update();
}
【描画結果】
5. 弾性体の曲げシミュレーション
Problem.Elastic.cs
public void ElasticProblem(MainWindow mainWindow)
{
CadObject2D cad2D = new CadObject2D();
{
IList<OpenTK.Vector2d> pts = new List<OpenTK.Vector2d>();
pts.Add(new OpenTK.Vector2d(0.0, 0.0));
pts.Add(new OpenTK.Vector2d(5.0, 0.0));
pts.Add(new OpenTK.Vector2d(5.0, 1.0));
pts.Add(new OpenTK.Vector2d(0.0, 1.0));
var res = cad2D.AddPolygon(pts);
}
double eLen = 0.1;
Mesher2D mesher2D = new Mesher2D(cad2D, eLen);
FEWorld world = new FEWorld();
world.Mesh = mesher2D;
uint quantityId;
{
uint dof = 2; // Vector2
uint feOrder = 1;
quantityId = world.AddQuantity(dof, feOrder);
}
world.TriIntegrationPointCount = TriangleIntegrationPointCount.Point3;
{
world.ClearMaterial();
uint maId = 0;
{
var ma = new LinearElasticMaterial();
ma.SetYoungPoisson(10.0, 0.3);
ma.GravityX = 0;
ma.GravityY = 0;
ma.MassDensity = 1.0;
maId = world.AddMaterial(ma);
}
uint lId = 1;
world.SetCadLoopMaterial(lId, maId);
}
uint[] zeroEIds = { 4 };
var zeroFixedCads = world.GetZeroFieldFixedCads(quantityId);
zeroFixedCads.Clear();
foreach (uint eId in zeroEIds)
{
uint dof = 2; // Vector2
var fixedCad = new FieldFixedCad(eId, CadElementType.Edge, FieldValueType.Vector2, dof);
zeroFixedCads.Add(fixedCad);
}
FieldFixedCad fixedCadXY;
{
// FixedDofIndex 0: X 1: Y
var fixedCadDatas = new[]
{
new { CadId = (uint)2, CadElemType = CadElementType.Edge,
FixedDofIndexs = new List<uint> { 0, 1 }, Values = new double[] { 0.0, 0.0 } }
};
IList<FieldFixedCad> fixedCads = world.GetFieldFixedCads(quantityId);
fixedCads.Clear();
foreach (var data in fixedCadDatas)
{
uint dof = 2; // Vector2
var fixedCad = new ConstFieldFixedCad(data.CadId, data.CadElemType,
FieldValueType.Vector2, dof, data.FixedDofIndexs, data.Values);
fixedCads.Add(fixedCad);
}
fixedCadXY = world.GetFieldFixedCads(quantityId)[0];
}
world.MakeElements();
uint valueId = 0;
var fieldDrawerArray = mainWindow.FieldDrawerArray;
{
uint dof = 2; // Vector2
world.ClearFieldValue();
valueId = world.AddFieldValue(FieldValueType.Vector2, FieldDerivativeType.Value,
quantityId, dof, false, FieldShowType.Real);
mainWindow.IsFieldDraw = true;
fieldDrawerArray.Clear();
{
IFieldDrawer faceDrawer = new FaceFieldDrawer(valueId, FieldDerivativeType.Value, false, world);
fieldDrawerArray.Add(faceDrawer);
}
IFieldDrawer edgeDrawer = new EdgeFieldDrawer(valueId, FieldDerivativeType.Value, false, world);
fieldDrawerArray.Add(edgeDrawer);
IFieldDrawer edgeDrawer2 = new EdgeFieldDrawer(valueId, FieldDerivativeType.Value, true, world);
fieldDrawerArray.Add(edgeDrawer2);
mainWindow.Camera.Fit(fieldDrawerArray.GetBoundingBox(mainWindow.Camera.RotMatrix33()));
mainWindow.glControl_ResizeProc();
//mainWindow.glControl.Invalidate();
//mainWindow.glControl.Update();
//WPFUtils.DoEvents();
}
double t = 0;
double dt = 0.05;
for (int iTime = 0; iTime <= 200; iTime++)
{
double[] fixedValueXY = fixedCadXY.GetDoubleValues();
fixedValueXY[0] = 0;
fixedValueXY[1] = Math.Sin(t * 2.0 * Math.PI * 0.1);
var FEM = new Elastic2DFEM(world);
{
//var solver = new IvyFEM.Linear.LapackEquationSolver();
//solver.Method = IvyFEM.Linear.LapackEquationSolverMethod.Dense;
//solver.Method = IvyFEM.Linear.LapackEquationSolverMethod.Band;
//solver.Method = IvyFEM.Linear.LapackEquationSolverMethod.PositiveDefiniteBand;
//FEM.Solver = solver;
}
{
//var solver = new IvyFEM.Linear.LisEquationSolver();
//solver.Method = IvyFEM.Linear.LisEquationSolverMethod.Default;
//FEM.Solver = solver;
}
{
var solver = new IvyFEM.Linear.IvyFEMEquationSolver();
solver.Method = IvyFEM.Linear.IvyFEMEquationSolverMethod.NoPreconCG;
//solver.Method = IvyFEM.Linear.IvyFEMEquationSolverMethod.CG;
//solver.Method = IvyFEM.Linear.IvyFEMEquationSolverMethod.ICCG;
//solver.Method = IvyFEM.Linear.IvyFEMEquationSolverMethod.NoPreconBiCGSTAB;
FEM.Solver = solver;
}
FEM.Solve();
double[] U = FEM.U;
world.UpdateFieldValueValuesFromNodeValues(valueId, FieldDerivativeType.Value, U);
fieldDrawerArray.Update(world);
mainWindow.glControl.Invalidate();
mainWindow.glControl.Update();
WPFUtils.DoEvents();
t += dt;
}
}
【計算結果】
5.1. 図面作成
CadObject2D
5.2. 有限要素法メッシュオブジェクト
Mesher2D
5.3. 有限要素法解析空間作成
FEWorld
5.4. 材質の設定
線形材料オブジェクトLinearElasticMaterialを指定します。
5.5. 境界条件設定
FixedCad
5.6. メッシュ計算
World.MakeElements()
5.7. 場の値を格納する場所を確保
FieldValue
5.8. 場の描画オブジェクト
場の描画オブジェクトFieldDrawerの生成(FaceFieldDrawer, EdgeFieldDrawer, etc)し、場の描画用アレイ(FieldDrawerArray)に追加します。
5.9. 有限要素法オブジェクト
問題に応じた有限要素法オブジェクトを生成します。
今回は弾性体の変位計算用のElastic2DFEM(名前FEM)を使います。
5.10. 連立方程式計算ライブラリの指定
IEquationSolver(FEM.Solverにセットする)
IvyFEM独自の実装のほかにLapackやLisを使用することができます。(問題によって適用できる計算方法とそうでないものがあるので注意が必要です）
5.11. 計算！
FEM.Solve()
5.12. 場の値を更新
World.UpdateFieldValueValuesFromNodeValues() と FieldDrawer.Update()
5.13. glControl再描画
glControlのInvalidate()とUpdate()を実行します。
6. まとめ
IvyFEMを用いて弾性体の曲げのシミュレーションをしました。
ここではすでにIvyFEMに用意されているシミュレーションルーチンを使って計算しましたが、シミュレーションの現場では独自に計算方法を実装したいと思われることもあると思います。その場合は、IvyFEMのコアな部分にも触れる必要があります。IvyFEMの一連のFEM群が参考になると思います。