.NET向け有限要素法CAEライブラリIvyFEMでカスタマイズ方程式を実装する

1. はじめに

.NET（C#）で有限要素法の計算がしたいということで開発を始めたライブラリIvyFEM.dllですが、スタンダードな定式化については対応していますが、それ以外は自分で実装できるような仕組みを持っています（それほど厳密なフレームワークではありませんが）。
ここでは、Poissonの方程式を例にとって実装手順を記します。
Poissonの方程式は
α∇^2(φ) = -F
です。

2. 実装するなクラス

Material：
　媒質のパラメータを保持するクラス
FEM：
　支配方程式の組み立て、連立方程式を解く処理を実装するクラスです。

以下各ファイルの先頭にusing IvyFEMしているものとします。namespaceはアプリのnamespaceで構いません。

3. Material

Materialクラスの派生クラスSampleFEMMaterialを実装します。αとFを材料定数として持っています。

SampleFEMMaterial.cs

    class SampleFEMMaterial : Material
    {
        public double Alpha { get => Values[0]; set => Values[0] = value; }
        public double F { get => Values[1]; set => Values[1] = value; }

        public SampleFEMMaterial() : base()
        {
            int len = 2;
            Values = new double[len];
            for (int i = 0; i < len; i++)
            {
                Values[i] = 0.0;
            }
        }
    }

4. FEM

Poissonの方程式の弱形式から、Ax = bのマトリクスAとベクトルbを求めます。そしてAx = bを解くという手順になります。
FEMはSolve()の実装が必須になっています。
そこから呼ばれるCalcAB()で有限要素マトリクスの構築しています。カスタマイズした定式化を実装したい場合は、CalcAB()を主に実装することになります。

SampleFEM.cs

    class SampleFEM : FEM
    {
        // output
        public double[] U { get; private set; } = null;

        public SampleFEM(FEWorld world)
        {
            World = world;
        }

        private void CalcAB(IvyFEM.Linear.DoubleSparseMatrix A, double[] B)
        {
            uint quantityId = 0;
            IList<uint> feIds = World.GetTriangleFEIds(quantityId);
            foreach (uint feId in feIds)
            {
                TriangleFE triFE = World.GetTriangleFE(quantityId, feId);
                uint elemNodeCnt = triFE.NodeCount;
                int[] nodes = new int[elemNodeCnt];
                for (int iNode = 0; iNode < elemNodeCnt; iNode++)
                {
                    int coId = triFE.NodeCoordIds[iNode];
                    int nodeId = World.Coord2Node(quantityId, coId);
                    nodes[iNode] = nodeId;
                }

                Material ma0 = World.GetMaterial(triFE.MaterialId);
                System.Diagnostics.Debug.Assert(ma0 is SampleFEMMaterial);
                var ma = ma0 as SampleFEMMaterial;
                double k = ma.Alpha;
                double f = ma.F;

                double[,] sNN = triFE.CalcSNN();
                double[,][,] sNuNv = triFE.CalcSNuNv();
                double[,] sNxNx = sNuNv[0, 0];
                double[,] sNyNx = sNuNv[1, 0];
                double[,] sNxNy = sNuNv[0, 1];
                double[,] sNyNy = sNuNv[1, 1];
                double[] sN = triFE.CalcSN();
                for (int row = 0; row < elemNodeCnt; row++)
                {
                    int rowNodeId = nodes[row];
                    if (rowNodeId == -1)
                    {
                        continue;
                    }
                    for (int col = 0; col < elemNodeCnt; col++)
                    {
                        int colNodeId = nodes[col];
                        if (colNodeId == -1)
                        {
                            continue;
                        }
                        double a = k * (sNxNx[row, col] + sNyNy[row, col]);

                        A[rowNodeId, colNodeId] += a;
                    }
                }
                for (int row = 0; row < elemNodeCnt; row++)
                {
                    int rowNodeId = nodes[row];
                    if (rowNodeId == -1)
                    {
                        continue;
                    }
                    B[rowNodeId] += f * sN[row];
                }
            }
        }

        public override void Solve()
        {
            uint quantityId = 0;
            int nodeCnt = (int)World.GetNodeCount(quantityId);
            var A = new IvyFEM.Linear.DoubleSparseMatrix(nodeCnt, nodeCnt);
            var B = new double[nodeCnt];
            CalcAB(A, B);

            DoubleSetFixedCadsCondtion(A, B, new int[] { nodeCnt }, new int[] { 1 });

            double[] X;
            Solver.DoubleSolve(out X, A, B);
            U = X;
        }
    }

5. カスタマイズした方程式を計算する

以上でカスタマイズ方程式のクラスはそろいました。
これらを使って有限要素法シミュレーションをしてみましょう。

Problem.SampleFEM.cs

        public void SampleFEMProblem(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(1.0, 0.0));
                pts.Add(new OpenTK.Vector2d(1.0, 1.0));
                pts.Add(new OpenTK.Vector2d(0.0, 1.0));
                var res = cad2D.AddPolygon(pts);
                uint addLId = res.AddLId;
                System.Diagnostics.Debug.Assert(addLId == 1);

                // ソース
                var res2 = cad2D.AddCircle(new OpenTK.Vector2d(0.5, 0.5), 0.1, addLId);
                System.Diagnostics.Debug.Assert(res2.AddLId == 2);
            }

            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();
            WPFUtils.DoEvents();

            double eLen = 0.05;
            Mesher2D mesher2D = new Mesher2D(cad2D, eLen);

            FEWorld world = new FEWorld();
            world.Mesh = mesher2D;
            uint quantityId;
            {
                uint dof = 1; // スカラー
                uint feOrder = 1;
                quantityId = world.AddQuantity(dof, feOrder);
            }

            {
                world.ClearMaterial();
                SampleFEMMaterial ma1 = new SampleFEMMaterial
                {
                    Alpha = 1.0,
                    F = 0.0
                };
                SampleFEMMaterial ma2 = new SampleFEMMaterial
                {
                    Alpha = 1.0,
                    F = 100.0
                };
                uint maId1 = world.AddMaterial(ma1);
                uint maId2 = world.AddMaterial(ma2);

                uint lId1 = 1;
                world.SetCadLoopMaterial(lId1, maId1);

                uint lId2 = 2;
                world.SetCadLoopMaterial(lId2, maId2);

                uint lId3 = 3;
                world.SetCadLoopMaterial(lId3, maId1);
            }

            uint[] zeroEIds = { 1, 2, 3, 4, 9, 10, 11, 12 };
            var zeroFixedCads = world.GetZeroFieldFixedCads(quantityId);
            zeroFixedCads.Clear();
            foreach (uint eId in zeroEIds)
            {
                uint dof = 1; // スカラー
                var fixedCad = new FieldFixedCad(eId, CadElementType.Edge, FieldValueType.Scalar, dof);
                zeroFixedCads.Add(fixedCad);
            }

            world.MakeElements();

            uint valueId = 0;
            var fieldDrawerArray = mainWindow.FieldDrawerArray;
            {
                uint dof = 1; // スカラー
                world.ClearFieldValue();
                valueId = world.AddFieldValue(FieldValueType.Scalar, FieldDerivativeType.Value,
                    quantityId, dof, false, FieldShowType.Real);
                mainWindow.IsFieldDraw = true;
                fieldDrawerArray.Clear();
                IFieldDrawer faceDrawer = new FaceFieldDrawer(valueId, FieldDerivativeType.Value, true, world,
                    valueId, FieldDerivativeType.Value);
                fieldDrawerArray.Add(faceDrawer);
                IFieldDrawer edgeDrawer = new EdgeFieldDrawer(
                    valueId, FieldDerivativeType.Value, true, false, world);
                fieldDrawerArray.Add(edgeDrawer);
                mainWindow.Camera.Fit(fieldDrawerArray.GetBoundingBox(mainWindow.Camera.RotMatrix33()));
                mainWindow.glControl_ResizeProc();
                //mainWindow.glControl.Invalidate();
                //mainWindow.glControl.Update();
                //WPFUtils.DoEvents();
            }

            {
                var FEM = new SampleFEM(world);
                {
                    //var solver = new IvyFEM.Linear.LapackEquationSolver();
                    //solver.Method = IvyFEM.Linear.LapackEquationSolverMethod.Dense;
                    //solver.IsOrderingToBandMatrix = true;
                    //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();
            }
        }

【計算結果】

6. まとめ

IvyFEMで提供されていない支配方程式を実装する方法をPoissonの方程式を例にとって説明しました。