GrasshopperでFEM解析をイチから実装する　その2の2 平面応力要素(四角形1次要素) 実装編

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

完成したGrasshopperのキャンバスは下記。記載のセクションごとに説明していきます。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4059820/d87ea080-5d62-40c9-8e1f-119c9aed71d4.png)
こんな感じで解析できるようになります！
![平面応力要素　kannsei.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4059820/cb194e35-221d-4787-97a4-8945e974dc4f.gif)

## メッシュの作成
今回実装するのは、平面応力解析であり２次元の解析なので、XY平面に限定して実装していきます。
XY平面上に適当なサイズの四角形を作成し、QuadRemeshでメッシュを作成し、MeshEdgesでメッシュのアウトラインを表示します。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4059820/5b44a098-9550-49c7-b97e-99ffd781f846.png)

メッシュが作成できました。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4059820/82709a43-eab9-4a2b-8b59-dca802447a88.png)

RhinoにおけるMeshの構成について少し解説しておきます。
MeshをDeconstructMeshで分解すると、VerticesとFacesに分解されます。VerticesはMeshの構成節点で、FacesはMeshを構成する面を表し、面を構成する節点のインデックスを保持します。この構造を利用して、FEM解析用の節点と要素の定義に利用します。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4059820/4e1eeba8-76bf-435b-b088-d42f17338cff.png)

## 材料特性・断面の設定
材料特性・断面の情報を設定します。今回は、
材料は鉄、厚さは10mmとして、下記のように設定しました。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4059820/13a9222a-bffc-433b-943f-8882c9a8b757.png)

## 境界条件の設定

先ほどDeconstructして得た節点から、X座標が0の点のXY変位を拘束する設定をします。
また、Meshの中の節点のうちどの節点に境界条件を与えるか指定するために、インデックスも取得しておきます。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4059820/ed43469c-1d37-4c59-8ee5-463ec44d2e36.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4059820/a787d543-8fcc-4cbe-af2d-3a01ddff4a7e.png)



C＃Scriptコンポーネントに、取得した境界条件を与える節点と、変位の拘束条件のBooleanToggleを入力につないで境界条件をタプルにひとまとめにしておきます。入力を右クリックして、意図したデータ型に設定します。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4059820/eebf897d-aa26-4ca8-8501-2c79ee9b3411.png)


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

```C#: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座標が部材長以上の節点を選択して荷重を与えています。
１節点に作用する荷重は作用させる荷重値/節点数としています。（線荷重とする場合は節点の負担幅ごとに荷重を分配すべきですが、簡単のために単純に節点数で割っています。）
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4059820/c27d81ae-8f24-4a94-86d6-e136fc48c508.png)

こちらも境界条件と同様にタプルで管理します。
```C#: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、応力Stressとしておきます。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4059820/1fd0fe75-40fd-437c-8dcc-523f28287bc1.png)


コードを書いていく前に、行列計算ライブラリをインポートします。
C#Scriptコンポーネントの左上の段ボール箱の中から、MathNet.NumericsをNuGetでインストールします。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4059820/c37a92b1-c492-4c94-a01b-904cc7486378.png)

インストールに成功すると下記のコードが追加されます。
```
#r "nuget: MathNet.Numerics"
```

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

```C#
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}
$$

```C#: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}
$$
~~~C#: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}
$$

~~~C#: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)}|
$$
~~~C#: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のプロパティとして定義して使用します。
~~~C#: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座標から要素剛性マトリクスを計算して、全体剛性マトリクスの必要な箇所に流し込む処理を書きます。

~~~C#: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];
                }
            }
        }
~~~
構成節点のインデックスをループ処理の中で取得しやすくするために関数を定義しています。
~~~C#: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))
        };
    }
~~~

### 荷重マトリクスの作成

入力から受け取った荷重条件をまとめたタプルから荷重マトリクスを組み上げます。
~~~C#: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;
        }
~~~

### 境界条件の適用

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

~~~C#: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;     
            }
        }
~~~

### 剛性方程式を解く！
いよいよあとは、行列計算をするだけ！
~~~C#::RunScriptの中身-4
        Matrix<double> d = K.Solve(f);
        Disp = d.Enumerate().ToList();
~~~

### 応力の計算
応力は算出した、変位から下記により求まります。
$$
\sigma =  D\epsilon = DBu
$$
これで求まる応力は積分点での応力なので、節点位置に外挿します。
また、各要素の節点位置で求まる応力は異なるので、滑らかな応力を描画するために、節点位置での応力の平均をとる処理を行います。
~~~C#
Stress = calc_stress(vertices,faces,d,D);
~~~
~~~C#
public DataTree<double> calc_stress(MeshVertexList vertices,MeshFaceList faces,Matrix<double> d,Matrix<double> D)
{
    Matrix<double> stress = Matrix<double>.Build.Dense(vertices.Count, 4);
    DataTree<double> res = new DataTree<double>();
    List<int> node_count = Enumerable.Repeat(0, vertices.Count).ToList();

    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> d_elem = Matrix<double>.Build.Dense(4 * 2, 1);
        for(int j = 0; j <4;j++)
        {
            int nodej = GetFaceVertexIndex(face,j);
            d_elem[j *2 ,0] = d[nodej * 2 ,0];
            d_elem[2 * j + 1, 0] = d[2 * nodej + 1, 0]; 
        }

        
        Matrix<double> stress_gps = Matrix<double>.Build.Dense(4, 3);
        for (int j = 0; j < gps.GetLength(0); j++)
        {
            double xi = gps[j, 0];
            double eta = gps[j, 1];
            Matrix<double> B = Matrix<double>.Build.Dense(3, 8);
            B = calcB(x,y,xi,eta);
            Matrix<double> stress_gp = D * B * d_elem;//積分点での応力を算定3x1matrix
            for (int k = 0; k < 3; k++)
            {
                stress_gps[j, k] += stress_gp[k, 0];
            }
        }

        //積分点応力から節点応力を計算
        Matrix<double> N = Matrix<double>.Build.DenseOfArray(new double[4, 4]{
            { 2 + Math.Sqrt(3.0),-1,2 - Math.Sqrt(3.0),-1},
            { -1,2 + Math.Sqrt(3.0),-1,2 - Math.Sqrt(3.0)},
            {2 - Math.Sqrt(3.0),-1,2 + Math.Sqrt(3.0),-1},
            { -1,2 - Math.Sqrt(3.0),-1,2 + Math.Sqrt(3.0)}
        });

        Matrix<double> stress_node = N / 2 * stress_gps; //4x3matrix要素内の節点

        for(int j = 0; j < 4 ; j++)
        {
            int nodej = GetFaceVertexIndex(face,j);
            double sigxx = stress_node[j, 0];
            double sigyy = stress_node[j, 1];
            double tauxy = stress_node[j, 2];
            double mises = Math.Sqrt(Math.Pow(sigxx, 2) - sigxx * sigyy + Math.Pow(sigyy, 2) + 3 * Math.Pow(tauxy, 2));

            node_count[nodej] += 1;
            stress[nodej, 0] += sigxx;
            stress[nodej, 1] += sigyy;
            stress[nodej, 2] += tauxy;
            stress[nodej, 3] += mises;
        }
    }
    for (int i = 0; i < vertices.Count; i++)
    {
        for(int j =0 ; j < 4; j++)
        {
            stress[i, j] /= node_count[i];
        }
    }
    for (int i = 0; i < stress.ColumnCount; i++)
    {
        GH_Path path = new GH_Path(i);
        for (int j = 0; j < vertices.Count; j++)
        {
            res.Add(stress[j, i], path);
        }
    }

    return res;
}
~~~

### 変形と応力の描画
計算した変形と応力を、Rhino上に図化してみましょう。
こちらもC#Scriptコンポーネントで、実装します。
MathNet.Numericsのインストールを再度行います。
入力は元のMesh・変位Disp・応力・描画する応力の種類・変形倍率DeformationRatioで、出力はDeformedMesh・凡例のColour・Tagとします。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4059820/51611f4a-f271-4aed-8ef1-d80a9ea48064.png)


~~~C#:DeformedMesh
    private void RunScript(
	Mesh Mesh,
	List<Vector3d> Disp,
	DataTree<double> Stress,
	int StressType,
	double DeformationRatio,
	ref object DeformedMesh,
	ref object Colour,
	ref object Tags)
    {
        //変形後のメッシュの描画
        Mesh deformedMesh = Mesh.DuplicateMesh();
        int vertexCount = deformedMesh.Vertices.Count;
        GH_Path path = new GH_Path(StressType);
        List<double> selectedStress = Stress.Branch(path);

        List<Color> nodesColor = GetColors(selectedStress);

        // 各頂点に変位を加える
        for (int i = 0; i < vertexCount; i++)
        {
            Point3d original = deformedMesh.Vertices[i];

            // 変位を反映
            Point3d moved = new Point3d(original) + Disp[i] * DeformationRatio;
            deformedMesh.VertexColors.SetColor(i,nodesColor[i]);
            deformedMesh.Vertices.SetVertex(i, moved);
        }

        //出力するタグの処理
        List<double> tagList = new List<double>();
        List<Color> tagColor = new List<Color>();
        int COLORLEVEL = 10;
        for (int i = 0; i < COLORLEVEL; i++)
        {
            tagList.Add(selectedStress.Min() + (selectedStress.Max() - selectedStress.Min()) / (COLORLEVEL - 1) * i );
            tagColor.Add(GetJetColor((double)i / (COLORLEVEL - 1)));
        }
        tagList.Reverse();
        tagColor.Reverse();

        Colour = tagColor;
        Tags = tagList;

        DeformedMesh = deformedMesh;
    }

    public static List<Color> GetColors(List<double> value)
    {
        List<Color> colors = new List<Color>();
        double minVal = value.Min();
        double maxVal = value.Max();
        foreach (double val in value)
        {
            double ratio = (val - minVal) / (maxVal - minVal);
            colors.Add(GetJetColor(ratio));
        }
        return colors;
    }

    public static Color GetJetColor(double value)
    {
        value = Math.Max(0.0, Math.Min(1.0, value)); // 0-1にクランプ
        double r = Math.Max(0, Math.Min(1, 1.5 - Math.Abs(4 * value - 3)));
        double g = Math.Max(0, Math.Min(1, 1.5 - Math.Abs(4 * value - 2)));
        double b = Math.Max(0, Math.Min(1, 1.5 - Math.Abs(4 * value - 1)));
        Color color = Color.FromArgb((int)(r * 255), (int)(g * 255), (int)(b * 255));
        return color;
    }
~~~
完成しました！それっぽい変形と応力が出ています。
![平面応力要素　henkeiouryoku.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4059820/45355eee-4794-44a0-b35c-91234ff4de64.gif)


## 検証
10mm x 10mm 長さ 100mm　の片持ち梁で検証してみます。
平面保持の梁理論によれば、片持ち梁の変形と最大応力の理論値は下式により求まります。
$$
\begin{gather}
\delta = PL^3/3EI \\\\
\sigma_{max} = M/Z = PL / (bh^2/6)
\end{gather}
$$
これと今回実装した、四角形１次要素の解析結果を比較してみましょう。
メッシュを細かくするほど理論値に近くなっていくことが分かります。
正しく実装できていそうですね。
四角形１次要素では、変形を１次式で表現するので、曲げ変形を表現できず、曲げ変形が卓越する変形では、メッシュが荒いと剛性が硬く表現されてしまいます。
![平面応力要素　test-1.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4059820/9062e718-bd48-4dce-a8a6-e5f1b46fc1d9.gif)

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

### 余談
全体剛性マトリクスを計算した後、行列計算をする部分を下記のように書き換えると、数学的には同じ計算をしていますが、計算時間がびっくりするくらい変わります。
~~~C#:
Matrix<double> d = K.Solve(f);
↓
Matrix<double> d = K.Inverse() * f;
~~~
4.0秒で計算できていたものが1.3分かかるようになりました。
違いは一度逆行列を計算しているかどうかですが、この逆行列を求める計算がマトリクスの規模が大きくなるほど、飛躍的に計算コストが増大していきます。
それを解決するために、いくつか方法があって、Solve()のメソッドにはその手法が含まれているというわけです。気になる方は、「LU分解」のキーワードで調べてみてください。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4059820/dcad1682-be01-44c6-8145-cd8f6a31d43e.png)

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/4059820/bb04ed40-4027-4538-a677-e461e6b74a0c.png)

## まとめ

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

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












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


~~~C#: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;
    }
}
~~~

~~~C#: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;
    }
}
~~~
~~~C#: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;
using Matrix = MathNet.Numerics.LinearAlgebra.Matrix<double>;



#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,
	ref object Stress)
    {
        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;
        List<Vector3d> dispVec = new List<Vector3d>();
        for(int i = 0; i < vertices.Count; i++)
        {
            double x = d[i * 2 ,0];
            double y = d[i * 2 + 1,0];
            double z = 0;
            dispVec.Add(new Vector3d(x,y,z));
        }
        Disp = dispVec;
        Stress = calc_stress(vertices,faces,d,D);
    }

    //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))
        };
    }

    public DataTree<double> calc_stress(MeshVertexList vertices,MeshFaceList faces,Matrix<double> d,Matrix<double> D)
    {
        Matrix<double> stress = Matrix<double>.Build.Dense(vertices.Count, 4);
        DataTree<double> res = new DataTree<double>();
        List<int> node_count = Enumerable.Repeat(0, vertices.Count).ToList();

        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> d_elem = Matrix<double>.Build.Dense(4 * 2, 1);
            for(int j = 0; j <4;j++)
            {
                int nodej = GetFaceVertexIndex(face,j);
                d_elem[j *2 ,0] = d[nodej * 2 ,0];
                d_elem[2 * j + 1, 0] = d[2 * nodej + 1, 0]; 
            }

            
            Matrix<double> stress_gps = Matrix<double>.Build.Dense(4, 3);
            for (int j = 0; j < gps.GetLength(0); j++)
            {
                double xi = gps[j, 0];
                double eta = gps[j, 1];
                Matrix<double> B = Matrix<double>.Build.Dense(3, 8);
                B = calcB(x,y,xi,eta);
                Matrix<double> stress_gp = D * B * d_elem;//積分点での応力を算定3x1matrix
                for (int k = 0; k < 3; k++)
                {
                    stress_gps[j, k] += stress_gp[k, 0];
                }
            }

            //積分点応力から節点応力を計算
            Matrix<double> N = Matrix<double>.Build.DenseOfArray(new double[4, 4]{
                { 2 + Math.Sqrt(3.0),-1,2 - Math.Sqrt(3.0),-1},
                { -1,2 + Math.Sqrt(3.0),-1,2 - Math.Sqrt(3.0)},
                {2 - Math.Sqrt(3.0),-1,2 + Math.Sqrt(3.0),-1},
                { -1,2 - Math.Sqrt(3.0),-1,2 + Math.Sqrt(3.0)}
            });

            Matrix<double> stress_node = N / 2 * stress_gps; //4x3matrix要素内の節点

            for(int j = 0; j < 4 ; j++)
            {
                int nodej = GetFaceVertexIndex(face,j);
                double sigxx = stress_node[j, 0];
                double sigyy = stress_node[j, 1];
                double tauxy = stress_node[j, 2];
                double mises = Math.Sqrt(Math.Pow(sigxx, 2) - sigxx * sigyy + Math.Pow(sigyy, 2) + 3 * Math.Pow(tauxy, 2));

                node_count[nodej] += 1;
                stress[nodej, 0] += sigxx;
                stress[nodej, 1] += sigyy;
                stress[nodej, 2] += tauxy;
                stress[nodej, 3] += mises;
            }
        }
        for (int i = 0; i < vertices.Count; i++)
        {
            for(int j =0 ; j < 4; j++)
            {
                stress[i, j] /= node_count[i];
            }
        }
        for (int i = 0; i < stress.ColumnCount; i++)
        {
            GH_Path path = new GH_Path(i);
            for (int j = 0; j < vertices.Count; j++)
            {
                res.Add(stress[j, i], path);
            }
        }

        return res;
    }
}
~~~


~~~C#:Render
// 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(
	Mesh Mesh,
	List<Vector3d> Disp,
	DataTree<double> Stress,
	int StressType,
	double DeformationRatio,
	ref object DeformedMesh,
	ref object Colour,
	ref object Tags)
    {
        //変形後のメッシュの描画
        Mesh deformedMesh = Mesh.DuplicateMesh();
        int vertexCount = deformedMesh.Vertices.Count;
        GH_Path path = new GH_Path(StressType);
        List<double> selectedStress = Stress.Branch(path);

        List<Color> nodesColor = GetColors(selectedStress);

        // 各頂点に変位を加える
        for (int i = 0; i < vertexCount; i++)
        {
            Point3d original = deformedMesh.Vertices[i];

            // 変位を反映
            Point3d moved = new Point3d(original) + Disp[i] * DeformationRatio;
            deformedMesh.VertexColors.SetColor(i,nodesColor[i]);
            deformedMesh.Vertices.SetVertex(i, moved);
        }

        //出力するタグの処理
        List<double> tagList = new List<double>();
        List<Color> tagColor = new List<Color>();
        int COLORLEVEL = 10;
        for (int i = 0; i < COLORLEVEL; i++)
        {
            tagList.Add(selectedStress.Min() + (selectedStress.Max() - selectedStress.Min()) / (COLORLEVEL - 1) * i );
            tagColor.Add(GetJetColor((double)i / (COLORLEVEL - 1)));
        }
        tagList.Reverse();
        tagColor.Reverse();

        Colour = tagColor;
        Tags = tagList;

        DeformedMesh = deformedMesh;
    }

    public static List<Color> GetColors(List<double> value)
    {
        List<Color> colors = new List<Color>();
        double minVal = value.Min();
        double maxVal = value.Max();
        foreach (double val in value)
        {
            double ratio = (val - minVal) / (maxVal - minVal);
            colors.Add(GetJetColor(ratio));
        }
        return colors;
    }

    public static Color GetJetColor(double value)
    {
        value = Math.Max(0.0, Math.Min(1.0, value)); // 0-1にクランプ
        double r = Math.Max(0, Math.Min(1, 1.5 - Math.Abs(4 * value - 3)));
        double g = Math.Max(0, Math.Min(1, 1.5 - Math.Abs(4 * value - 2)));
        double b = Math.Max(0, Math.Min(1, 1.5 - Math.Abs(4 * value - 1)));
        Color color = Color.FromArgb((int)(r * 255), (int)(g * 255), (int)(b * 255));
        return color;
    }

}

~~~