はじめに

前回の記事では、FEM解析の理論について説明ました。
本記事は実装編ということで、Grasshopperの環境で、四角形１次要素の平面応力解析を実装していきます。
コードを記述する部分は、C#Scriptコンポーネントを使用していきます。
今回作成したコード群は下記に公開しています。
https://github.com/hrdsh-std/QiitaRepo/tree/main/2-2

完成したGrasshopperのキャンバスは下記。記載のセクションごとに説明していきます。

こんな感じで解析できるようになります！

メッシュの作成

今回実装するのは、平面応力解析であり２次元の解析なので、XY平面に限定して実装していきます。
XY平面上に適当なサイズの四角形を作成し、QuadRemeshでメッシュを作成し、MeshEdgesでメッシュのアウトラインを表示します。

メッシュが作成できました。

RhinoにおけるMeshの構成について少し解説しておきます。
MeshをDeconstructMeshで分解すると、VerticesとFacesに分解されます。VerticesはMeshの構成節点で、FacesはMeshを構成する面を表し、面を構成する節点のインデックスを保持します。この構造を利用して、FEM解析用の節点と要素の定義に利用します。

材料特性・断面の設定

材料特性・断面の情報を設定します。今回は、
材料は鉄、厚さは10mmとして、下記のように設定しました。

境界条件の設定

先ほどDeconstructして得た節点から、X座標が0の点のXY変位を拘束する設定をします。
また、Meshの中の節点のうちどの節点に境界条件を与えるか指定するために、インデックスも取得しておきます。

C＃Scriptコンポーネントに、取得した境界条件を与える節点と、変位の拘束条件のBooleanToggleを入力につないで境界条件をタプルにひとまとめにしておきます。入力を右クリックして、意図したデータ型に設定します。

C#ScriptのRunScript内のコードは下記。
境界条件のID,適用する節点ID、XYの変位拘束を保持します。

Boundary

    private void RunScript(
	List<int> VerticesIDs,
	bool DX,
	bool DY,
	ref object Boundaries)
    {
        var res = new List<(int ID,int VertexID,bool DX,bool DY)>();
        for(int i = 0; i < VerticesIDs.Count ; i++ )
        {
            res.Add((i,VerticesIDs[i],DX,DY));
        }
        Boundaries = res;
    }

荷重の設定

荷重の設定も基本的には境界条件と同じように設定していきます。
片持ち張り先端に荷重値をかける設定としました。X座標が部材長以上の節点を選択して荷重を与えています。
１節点に作用する荷重は作用させる荷重値/節点数としています。（線荷重とする場合は節点の負担幅ごとに荷重を分配すべきですが、簡単のために単純に節点数で割っています。）

こちらも境界条件と同様にタプルで管理します。

Boundary

     private void RunScript(
	List<int> VerticesIDs,
	double FX,
	double FY,
	ref object Loads)
    {
        var res = new List<(int ID,int VertexID,double FX,double FY)>();
        for(int i = 0; i < VerticesIDs.Count ; i++ )
        {
            res.Add((i,VerticesIDs[i],FX,FY));
        }
        Loads = res;
    }

解析実行

下準備が済んだので、ソルバーの部分を作っていきます。
これまで作成した、メッシュ、材料特性、断面、境界条件、荷重　をC#Scriptの入力につなぎます。
出力は変位Dispとしておきます。

コードを書いていく前に、行列計算ライブラリをインポートします。
C#Scriptコンポーネントの左上の段ボール箱の中から、MathNet.NumericsをNuGetでインストールします。

インストールに成功すると下記のコードが追加されます。

#r "nuget: MathNet.Numerics"

インストールしたライブラリを使用するため、下記のコードを追加します。
RhinoのAPIにも Rhino.Geometry.Matrix という同名のクラスが存在しており、MathNetの Matrixとのコンフリクトを避ける必要があります。

using MathNet.Numerics.LinearAlgebra;
using Matrix = MathNet.Numerics.LinearAlgebra.Matrix<double>;

ここからはC#Scriptの中のコードを要素ごとに分割して説明していきます。

Dマトリクスの計算

Dマトリクスを計算する関数を準備します。Dマトリクスは材質から定まる値で、Ｅとνから求めます。
このコードは、RunScriptの外に記述します。式は下記の通り。

$$
\mathbf{D} =
\frac{E}{(1 - \nu^2)}
\begin{bmatrix}
1 - \nu & \nu & 0
\\ \nu & 1 - \nu & 0
\\ 0 & 0 & \frac{1 - \nu}{2}
\end{bmatrix}
$$

calcD

    private Matrix<double> calcD(double E , double nu)
    {
        Matrix<double> D = Matrix<double>.Build.DenseOfArray(new double[,]
        {
            {1, nu, 0},
            {nu, 1, 0},
            {0, 0, (1-nu)/2}
        });
        D = D * E / (1 - Math.Pow(nu, 2));
        return D;
    }

Jマトリクスの計算

ヤコビアンの計算も関数にまとめておきます。
細かな式展開は前回の記事を参照ください。
$$
J =
\begin{bmatrix}
\frac{\partial x}{\partial\xi} & \frac{\partial y}{\partial\xi}\
\\ \frac{\partial x}{\partial\eta} & \frac{\partial y}{\partial\eta}\
\end{bmatrix}
$$

calcJ

    private Matrix<double> calcJ(Matrix<double> x,Matrix<double> y,double xi,double eta)
    {
        Matrix<double> dndxi = Matrix<double>.Build.DenseOfArray(new double[,]
        {
                {(-1+eta)/4,(1-eta)/4,(1+eta)/4,(-1-eta)/4 }
        });
        Matrix<double> dndeta = Matrix<double>.Build.DenseOfArray(new double[,]
        {
                {(-1+xi)/4,(-1-xi)/4,(1+xi)/4,(1-xi)/4 }
        });

        Matrix<double> dxdxi = dndxi * x;
        Matrix<double> dydxi = dndxi * y;
        Matrix<double> dxdeta = dndeta * x;
        Matrix<double> dydeta = dndeta * y;

        Matrix<double> J = Matrix<double>.Build.DenseOfArray(new double[,]
        {
            {dxdxi[0,0],dydxi[0,0]},
            {dxdeta[0, 0],dydeta[0, 0]}
        });
        return J;
    }

Bマトリクスの計算

Bマトリクスの計算も関数にまとめておきます。
式は下記の通り。
$$
\mathbf{B} =
\begin{bmatrix}
\frac{\partial N_1}{\partial x} & 0 &\cdots& \frac{\partial N_4}{\partial x} & 0
\\ 0 & \frac{\partial N_1}{\partial y} &\cdots& 0 & \frac{\partial N_4}{\partial y}
\\ \frac{\partial N_1}{\partial y} & \frac{\partial N_1}{\partial x} &\cdots& \frac{\partial N_4}{\partial y} & \frac{\partial N_4}{\partial x}
\end{bmatrix}
$$

$$
\begin{Bmatrix}
\frac{\partial N_i}{\partial x}
\\ \frac{\partial N_i}{\partial y}
\end{Bmatrix}
=[\mathbf{J}]^{-1}
\begin{Bmatrix}
\frac{\partial N_i}{\partial \xi}
\\ \frac{\partial N_i}{\partial \eta}
\end{Bmatrix}
$$

calcB

    private Matrix<double> calcB(Matrix<double> x,Matrix<double> y,double xi,double eta)
    {
        Matrix<double> dndxi = Matrix<double>.Build.DenseOfArray(new double[,]
        {
                {(-1+eta)/4,(1-eta)/4,(1+eta)/4,(-1-eta)/4 }
        });
        Matrix<double> dndeta = Matrix<double>.Build.DenseOfArray(new double[,]
        {
                {(-1+xi)/4,(-1-xi)/4,(1+xi)/4,(1-xi)/4 }
        });
        
        Matrix<double> J = calcJ(x,y,xi,eta);
        Matrix<double> B = Matrix<double>.Build.Dense(3, 8);

        for (int i = 0; i < 4; i++)
        {
            Matrix<double> Bi = J.Solve(Matrix<double>.Build.DenseOfArray(new double[2, 1] { { dndxi[0, i] }, { dndeta[0, i] } }));
            B[0, 2 * i] = Bi[0, 0];
            B[1, 2 * i + 1] = Bi[1, 0];
            B[2, 2 * i] = Bi[1, 0];
            B[2, 2 * i + 1] = Bi[0, 0];
        }
        return B; 
    }

Keマトリクスの計算

求めたD、Bマトリクスから数値積分を使って、要素剛性マトリクスを計算します。
$$
\mathbf{K}=
t \sum_{i=1}^{2}\sum_{j=1}^{2} w_iw_j\mathbf{B}(\xi_i,\eta_i)^\mathsf{T}\mathbf{D}\mathbf{B}(\xi_i,\eta_i)det|\mathbf{J(\xi_i,\eta_i)}|
$$

calcKe

    private Matrix<double> calcKe(Matrix<double> x,Matrix<double> y,double section ,Matrix<double> D,double[,]gps)
    {
        Matrix<double> Ke = Matrix<double>.Build.DenseOfArray(new double[8, 8]);//８ｘ８のゼロ行列
        for (int i = 0; i < gps.GetLength(0); i++)
        {
            double xi = gps[i, 0];
            double eta = gps[i, 1];
            double wi = gps[i, 2];
            double wj = gps[i, 3];
            Matrix<double> J = calcJ(x,y, xi, eta);
            Matrix<double> B = calcB(x,y, xi, eta);
            Matrix<double> Kei = B.Transpose() * D * B * J.Determinant() * section * wi * wj;
            Ke += Kei;
        }
        return Ke;
    }

積分点GPS(Gauss Points)はScript_Instanceのプロパティとして定義して使用します。

gps

    double[,] gps = new double[,]
    {
        {-1.0 / Math.Sqrt(3.0), -1.0 / Math.Sqrt(3.0), 1.0, 1.0 },
        { 1.0 / Math.Sqrt(3.0), -1.0 / Math.Sqrt(3.0), 1.0, 1.0},
        { 1.0 / Math.Sqrt(3.0), 1.0 / Math.Sqrt(3.0), 1.0, 1.0},
        { -1.0 / Math.Sqrt(3.0), 1.0 / Math.Sqrt(3.0), 1.0, 1.0}//(xi,eta,wi,wj)
    };

これで必要な要素は準備できました。

Kマトリクスの計算

次に、計算した要素剛性マトリクスを組み上げて、全体剛性マトリクスを作成します。
各節点は並進２成分の自由度を持つので、全体剛性マトリクスＫのサイズは、（節点数 × 2）行 × （節点数 × 2）列となります。
Meshの面は　face.A　のようにして、構成節点のインデックスを取得することができるので、Vertices[face.A].XY で構成節点のXY座標を取得しています。
各要素でループを回して、ループの中で各節点のXY座標から要素剛性マトリクスを計算して、全体剛性マトリクスの必要な箇所に流し込む処理を書きます。

RunScriptの中身-1

    private void RunScript(
	Mesh Mesh,
	double E,
	double Nu,
	double Section,
	List<object> Boundaries,
	List<object> Loads,
	ref object Disp)
    {
        Matrix<double> D = calcD(E,Nu);
        var faces = Mesh.Faces;
        var vertices = Mesh.Vertices;

        //剛性マトリクスの計算
        Matrix<double> K = Matrix<double>.Build.Dense(vertices.Count * 2, vertices.Count * 2);
        for(int i = 0; i < faces.Count; i++)
        {
            MeshFace face = faces.GetFace(i);
            Matrix<double> x = Matrix<double>.Build.DenseOfArray(new double[,]
            {
                {vertices[face.A].X},{vertices[face.B].X},{vertices[face.C].X},{vertices[face.D].X}
            });
            Matrix<double> y = Matrix<double>.Build.DenseOfArray(new double[,]
            {
                {vertices[face.A].Y},{vertices[face.B].Y},{vertices[face.C].Y},{vertices[face.D].Y}
            });
            Matrix<double> Ke = calcKe(x,y,Section,D,gps);

            for (int j = 0; j < 4; j++)
            {
                for (int k = 0; k < 4; k++)
                {
                    int nodej = GetFaceVertexIndex(face,j);
                    int nodek = GetFaceVertexIndex(face,k);
                    K[2 * nodej, 2 * nodek] += Ke[2 * j, 2 * k];
                    K[2 * nodej, 2 * nodek + 1] += Ke[2 * j, 2 * k + 1];
                    K[2 * nodej + 1, 2 * nodek] += Ke[2 * j + 1, 2 * k];
                    K[2 * nodej + 1, 2 * nodek + 1] += Ke[2 * j + 1, 2 * k + 1];
                }
            }
        }

構成節点のインデックスをループ処理の中で取得しやすくするために関数を定義しています。

GetFaceVertexIndex

    private int GetFaceVertexIndex(MeshFace face, int i)
    {
        return i switch
        {
            0 => face.A,
            1 => face.B,
            2 => face.C,
            3 => face.IsQuad ? face.D : face.C, // 三角形のときD=C
            _ => throw new ArgumentOutOfRangeException(nameof(i))
        };
    }

荷重マトリクスの作成

入力から受け取った荷重条件をまとめたタプルから荷重マトリクスを組み上げます。

RunScriptの中身-2

        //荷重条件の設定
        Matrix<double> f = Matrix<double>.Build.Dense(vertices.Count * 2, 1);
        for(int i = 0 ; i < Loads.Count ; i++ )
        {
            var load = ((int ID,int VertexID,double FX,double FY)) Loads[i];
            f[load.VertexID * 2, 0] = load.FX;
            f[load.VertexID * 2 + 1, 0] = load.FY;
        }

境界条件の適用

入力から受け取った境界条件をまとめたタプルから全体剛性マトリクスの、対応する箇所の数値を境界条件に応じて書き換えていきます。

RunScriptの中身-3

        //境界条件の設定
        for(int i = 0; i < Boundaries.Count; i++)
        {
            var  boundary = ((int ID,int VertexID,bool DX,bool DY)) Boundaries[i];
            if(boundary.DX)
            {
                for (int j = 0; j < vertices.Count * 2; j++)
                {
                    K[boundary.VertexID * 2, j] = 0.0;
                    K[j, boundary.VertexID * 2] = 0.0;
                }
                K[boundary.VertexID * 2, boundary.VertexID * 2] = 1.0; 
                f[boundary.VertexID * 2,0] = 0.0;    
            }
            if(boundary.DY)
            {
                for (int j = 0; j < vertices.Count * 2; j++)
                {
                    K[boundary.VertexID * 2 + 1, j] = 0.0;
                    K[j, boundary.VertexID * 2 + 1] = 0.0;
                }
                K[boundary.VertexID * 2 + 1, boundary.VertexID * 2 + 1] = 1.0;    
                f[boundary.VertexID * 2 + 1,0] = 0.0;     
            }
        }

剛性方程式を解く！

いよいよあとは、行列計算をするだけ！

:RunScriptの中身-4

        Matrix<double> d = K.Solve(f);
        Disp = d.Enumerate().ToList();

変形の描画

元のメッシュの各節点の座標に変位を足して、Rhino上に図化してみましょう。
こちらもC#Scriptコンポーネントで、実装します。
MathNet.Numericsのインストールを再度行います。
入力は元のMesh・変位Disp・変形倍率DeformationRatioで、出力はDeformedMeshとします。

DeformedMesh

    private void RunScript(
	Mesh Mesh,
	List<double> Disp,
	double DeformationRatio,
	ref object DeformedMesh)
    {
        //変形後のメッシュの描画
        Mesh deformedMesh = Mesh.DuplicateMesh();
        int vertexCount = deformedMesh.Vertices.Count;

        // 各頂点に変位を加える
        for (int i = 0; i < vertexCount; i++)
        {
            Point3d original = deformedMesh.Vertices[i];
            Vector3d vec = new Vector3d(Disp[i * 2],Disp[i * 2 + 1],0.0);

            // 変位を反映
            Point3d moved = new Point3d(original) + vec * DeformationRatio;
            deformedMesh.Vertices.SetVertex(i, moved);
        }

        DeformedMesh = deformedMesh;
    }

完成しました！それっぽい変形が出ています。

検証

まずは、10mm x 10mm 長さ 100mm　の片持ち梁で検証してみます。
平面保持の梁理論によれば、片持ち梁の変形の理論値は下式により求まります。
$$
\delta = PL^3/3EI
$$
これと今回実装した、四角形１次要素の解析結果を比較してみましょう。
メッシュを細かくするほど理論値に近くなっていくことが分かります。
正しく実装できていそうですね。
四角形１次要素では、変形を１次式で表現するので、曲げ変形を表現できず、曲げ変形が卓越する変形では、メッシュが荒いと剛性が硬く表現されてしまいます。

次は長さを短くして、30mmとして解析し理論値と比較してみます。
今度はせん断変形を無視した梁理論の理論値を早い段階で追い越しています。

一般的な商用の解析ソフトだと、メッシュを変えてすぐに結果が返ってくるという体験はできないので、Grasshopperで構造解析をするのは楽しいですね。

余談

全体剛性マトリクスを計算した後、行列計算をする部分を下記のように書き換えると、数学的には同じ計算をしていますが、計算時間がびっくりするくらい変わります。

Matrix<double> d = K.Solve(f);
↓
Matrix<double> d = K.Inverse() * f;

4.0秒で計算できていたものが1.3分かかるようになりました。
違いは一度逆行列を計算しているかどうかですが、この逆行列を求める計算がマトリクスの規模が大きくなるほど、飛躍的に計算コストが増大していきます。
それを解決するために、いくつか方法があって、Solve()のメソッドにはその手法が含まれているというわけです。気になる方は、「LU分解」のキーワードで調べてみてください。

まとめ

本記事では、Grasshopper上でFEM解析をイチから実装しました。
次回は、応力を描画する処理について記事にしたいと思います。

よくわからない箇所、詳しく説明してほしい箇所があれば、コメントください！
なんでもリアクションをいただけると、とてもうれしくて記事を続けるモチベーションになります。

記事の中では、説明のためコードを要素ごとにバラバラに示したので、C#Scriptコンポーネントごとのコードを貼り付けておきます。

Bounadries

// Grasshopper Script Instance
#region Usings
using System;
using System.Linq;
using System.Collections;
using System.Collections.Generic;
using System.Drawing;

using Rhino;
using Rhino.Geometry;

using Grasshopper;
using Grasshopper.Kernel;
using Grasshopper.Kernel.Data;
using Grasshopper.Kernel.Types;
#endregion

public class Script_Instance : GH_ScriptInstance
{
    #region Notes
    /* 
      Members:
        RhinoDoc RhinoDocument
        GH_Document GrasshopperDocument
        IGH_Component Component
        int Iteration

      Methods (Virtual & overridable):
        Print(string text)
        Print(string format, params object[] args)
        Reflect(object obj)
        Reflect(object obj, string method_name)
    */
    #endregion

    private void RunScript(
	List<int> VerticesIDs,
	bool DX,
	bool DY,
	ref object Boundaries)
    {
        var res = new List<(int ID,int VertexID,bool DX,bool DY)>();
        for(int i = 0; i < VerticesIDs.Count ; i++ )
        {
            res.Add((i,VerticesIDs[i],DX,DY));
        }
        Boundaries = res;
    }
}

Loads

// Grasshopper Script Instance
#region Usings
using System;
using System.Linq;
using System.Collections;
using System.Collections.Generic;
using System.Drawing;

using Rhino;
using Rhino.Geometry;

using Grasshopper;
using Grasshopper.Kernel;
using Grasshopper.Kernel.Data;
using Grasshopper.Kernel.Types;
#endregion

public class Script_Instance : GH_ScriptInstance
{
    #region Notes
    /* 
      Members:
        RhinoDoc RhinoDocument
        GH_Document GrasshopperDocument
        IGH_Component Component
        int Iteration

      Methods (Virtual & overridable):
        Print(string text)
        Print(string format, params object[] args)
        Reflect(object obj)
        Reflect(object obj, string method_name)
    */
    #endregion

    private void RunScript(
	List<int> VerticesIDs,
	double FX,
	double FY,
	ref object Loads)
    {
        var res = new List<(int ID,int VertexID,double FX,double FY)>();
        for(int i = 0; i < VerticesIDs.Count ; i++ )
        {
            res.Add((i,VerticesIDs[i],FX,FY));
        }
        Loads = res;
    }
}

Solver

// Grasshopper Script Instance
#region Usings
using System;
#r "nuget: MathNet.Numerics"

using System.Linq;
using System.Collections;
using System.Collections.Generic;
using System.Drawing;

using Rhino;
using Rhino.Geometry;
using Rhino.Geometry.Collections;


using Grasshopper;
using Grasshopper.Kernel;
using Grasshopper.Kernel.Data;
using Grasshopper.Kernel.Types;

using MathNet.Numerics.LinearAlgebra;


#endregion

public class Script_Instance : GH_ScriptInstance
{
    #region Notes
    /* 
      Members:
        RhinoDoc RhinoDocument
        GH_Document GrasshopperDocument
        IGH_Component Component
        int Iteration

      Methods (Virtual & overridable):
        Print(string text)
        Print(string format, params object[] args)
        Reflect(object obj)
        Reflect(object obj, string method_name)
    */
    #endregion

    double[,] gps = new double[,]
    {
        {-1.0 / Math.Sqrt(3.0), -1.0 / Math.Sqrt(3.0), 1.0, 1.0 },
        { 1.0 / Math.Sqrt(3.0), -1.0 / Math.Sqrt(3.0), 1.0, 1.0},
        { 1.0 / Math.Sqrt(3.0), 1.0 / Math.Sqrt(3.0), 1.0, 1.0},
        { -1.0 / Math.Sqrt(3.0), 1.0 / Math.Sqrt(3.0), 1.0, 1.0}//(xi,eta,wi,wj)
    };
    
    private void RunScript(
	Mesh Mesh,
	double E,
	double Nu,
	double Section,
	List<object> Boundaries,
	List<object> Loads,
	ref object Disp)
    {
        Matrix<double> D = calcD(E,Nu);
        var faces = Mesh.Faces;
        var vertices = Mesh.Vertices;

        //剛性マトリクスの計算
        Matrix<double> K = Matrix<double>.Build.Dense(vertices.Count * 2, vertices.Count * 2);
        for(int i = 0; i < faces.Count; i++)
        {
            MeshFace face = faces.GetFace(i);
            Matrix<double> x = Matrix<double>.Build.DenseOfArray(new double[,]
            {
                {vertices[face.A].X},{vertices[face.B].X},{vertices[face.C].X},{vertices[face.D].X}
            });
            Matrix<double> y = Matrix<double>.Build.DenseOfArray(new double[,]
            {
                {vertices[face.A].Y},{vertices[face.B].Y},{vertices[face.C].Y},{vertices[face.D].Y}
            });
            Matrix<double> Ke = calcKe(x,y,Section,D,gps);

            for (int j = 0; j < 4; j++)
            {
                for (int k = 0; k < 4; k++)
                {
                    int nodej = GetFaceVertexIndex(face,j);
                    int nodek = GetFaceVertexIndex(face,k);
                    K[2 * nodej, 2 * nodek] += Ke[2 * j, 2 * k];
                    K[2 * nodej, 2 * nodek + 1] += Ke[2 * j, 2 * k + 1];
                    K[2 * nodej + 1, 2 * nodek] += Ke[2 * j + 1, 2 * k];
                    K[2 * nodej + 1, 2 * nodek + 1] += Ke[2 * j + 1, 2 * k + 1];
                }
            }
        }

        //荷重条件の設定
        Matrix<double> f = Matrix<double>.Build.Dense(vertices.Count * 2, 1);
        for(int i = 0 ; i < Loads.Count ; i++ )
        {
            var load = ((int ID,int VertexID,double FX,double FY)) Loads[i];
            f[load.VertexID * 2, 0] = load.FX;
            f[load.VertexID * 2 + 1, 0] = load.FY;
        }


        //境界条件の設定
        for(int i = 0; i < Boundaries.Count; i++)
        {
            var  boundary = ((int ID,int VertexID,bool DX,bool DY)) Boundaries[i];
            if(boundary.DX)
            {
                for (int j = 0; j < vertices.Count * 2; j++)
                {
                    K[boundary.VertexID * 2, j] = 0.0;
                    K[j, boundary.VertexID * 2] = 0.0;
                }
                K[boundary.VertexID * 2, boundary.VertexID * 2] = 1.0; 
                f[boundary.VertexID * 2,0] = 0.0;    
            }
            if(boundary.DY)
            {
                for (int j = 0; j < vertices.Count * 2; j++)
                {
                    K[boundary.VertexID * 2 + 1, j] = 0.0;
                    K[j, boundary.VertexID * 2 + 1] = 0.0;
                }
                K[boundary.VertexID * 2 + 1, boundary.VertexID * 2 + 1] = 1.0;    
                f[boundary.VertexID * 2 + 1,0] = 0.0;     
            }
        }
        

        //行列計算
        Matrix<double> d = K.Solve(f);
        //Matrix<double> d = K.Inverse() * f;
        Disp = d.Enumerate().ToList();
    }

    //Dマトリクスの算定
    private Matrix<double> calcD(double E , double nu)
    {
        Matrix<double> D = Matrix<double>.Build.DenseOfArray(new double[,]
        {
            {1, nu, 0},
            {nu, 1, 0},
            {0, 0, (1-nu)/2}
        });
        D = D * E / (1 - Math.Pow(nu, 2));
        return D;
    }
    //Jマトリクスの算定
    private Matrix<double> calcJ(Matrix<double> x,Matrix<double> y,double xi,double eta)
    {
        Matrix<double> dndxi = Matrix<double>.Build.DenseOfArray(new double[,]
        {
                {(-1+eta)/4,(1-eta)/4,(1+eta)/4,(-1-eta)/4 }
        });
        Matrix<double> dndeta = Matrix<double>.Build.DenseOfArray(new double[,]
        {
                {(-1+xi)/4,(-1-xi)/4,(1+xi)/4,(1-xi)/4 }
        });

        Matrix<double> dxdxi = dndxi * x;
        Matrix<double> dydxi = dndxi * y;
        Matrix<double> dxdeta = dndeta * x;
        Matrix<double> dydeta = dndeta * y;

        Matrix<double> J = Matrix<double>.Build.DenseOfArray(new double[,]
        {
            {dxdxi[0,0],dydxi[0,0]},
            {dxdeta[0, 0],dydeta[0, 0]}
        });
        return J;
    }
    private Matrix<double> calcB(Matrix<double> x,Matrix<double> y,double xi,double eta)
    {
        Matrix<double> dndxi = Matrix<double>.Build.DenseOfArray(new double[,]
        {
                {(-1+eta)/4,(1-eta)/4,(1+eta)/4,(-1-eta)/4 }
        });
        Matrix<double> dndeta = Matrix<double>.Build.DenseOfArray(new double[,]
        {
                {(-1+xi)/4,(-1-xi)/4,(1+xi)/4,(1-xi)/4 }
        });
        
        Matrix<double> J = calcJ(x,y,xi,eta);
        Matrix<double> invJ = J.Inverse();

        Matrix<double> B = Matrix<double>.Build.Dense(3, 8);

        for (int i = 0; i < 4; i++)
        {
            //負荷が大きいので、逆行列は計算しない
            //Matrix<double> Bi = invJ * Matrix<double>.Build.DenseOfArray(new double[2, 1] { { dndxi[0, i] }, { dndeta[0, i] } });
            Matrix<double> Bi = J.Solve(Matrix<double>.Build.DenseOfArray(new double[2, 1] { { dndxi[0, i] }, { dndeta[0, i] } }));
            B[0, 2 * i] = Bi[0, 0];
            B[1, 2 * i + 1] = Bi[1, 0];
            B[2, 2 * i] = Bi[1, 0];
            B[2, 2 * i + 1] = Bi[0, 0];
            
        }
        return B; 
    }
    private Matrix<double> calcKe(Matrix<double> x,Matrix<double> y,double section ,Matrix<double> D,double[,]gps)
    {
        Matrix<double> Ke = Matrix<double>.Build.DenseOfArray(new double[8, 8]);//８ｘ８のゼロ行列
        for (int i = 0; i < gps.GetLength(0); i++)
        {
            double xi = gps[i, 0];
            double eta = gps[i, 1];
            double wi = gps[i, 2];
            double wj = gps[i, 3];
            Matrix<double> J = calcJ(x,y, xi, eta);
            Matrix<double> B = calcB(x,y, xi, eta);
            Matrix<double> Kei = B.Transpose() * D * B * J.Determinant() * section * wi * wj;
            Ke += Kei;
        }
        return Ke;
    }
    private int GetFaceVertexIndex(MeshFace face, int i)
    {
        return i switch
        {
            0 => face.A,
            1 => face.B,
            2 => face.C,
            3 => face.IsQuad ? face.D : face.C, // 三角形のときD=C
            _ => throw new ArgumentOutOfRangeException(nameof(i))
        };
    }
}