<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 はじめに
静的な場合の有限要素法問題における変位は以下のような連立方程式をとくことにより得られます。
$KU=F$ 
$HU=R$
ここで、$K$は剛性行列、$U$は変位ベクトルです。
一方、$F$は外力ベクトル(NEUMANN条件)であり、$HU=R$は固定端やローラー条件(DIRICHLET条件)です。
この資料は三脚のメッシュをpythonでインポートし、境界条件を設定し剛性行列$K$と外力ベクトル$F$を作成して変位$U$を計算する過程を説明したものです。
インストールとモジュールのインポート
インストール方法
Linux(Ubuntu)
まずはGetfem++のpythonインターフェースのモジュールをインポートします。Ubuntu系のOSではインストールコマンドは
$ sudo apt-get update; sudo apt-get install python-getfem++
です。さらに、境界条件設定に使うためnumpyをインポートします。numpyはpython-getfem++と依存関係にあるため、自動的にインストールされます。インストールが終わったら、結果出力とGUI表示のため、gmshとipython-notebookをインストールしましょう。
$ sudo apt-get install gmsh python-pip; sudo pip install notebook
このファイルがあるディレクトリに移動
$ ipython notebook
で、Ipython Notebookを起動した後このファイルを開いてください。
Windows
Python 2.7系の32bit版をインストール(Getfem++のWindows版が32bit版しかないので、32bit版でないと動かない)
https
Step1:  パラメータの設定
まずはパラメータを設定します。この例では、このファイルからの相対パスが './mesh/tripod.mesh'であるメッシュファイルを使用します。メッシュは四面体要素です。ヤング係数$E = 1.0\times10^3$、ポアソン比$\nu = 0.3$とします。
Step2:  メッシュの読み込み
メッシュを外部ファイルから読み込みます。メッシュの入力データはGetfem++の独自のフォーマットになっていますが、解読は容易だと思います。
Step3:  次のコマンドで、mに設定したメッシュをgmshのスクリプトでpng画像に打ち出し確認します。
Step4:  以下のようにipythonの機能を使用して、システムのコマンドラインでの操作を行いIpython上に画像を表示してみます。もちろん、gmshで直接posファイルを表示することもできます。
Step5:  これは横から見たメッシュ図です。三脚のメッシュであることがわかります。今回はこの頂点部分に荷重(NEUMANN条件)、下端に固定条件(DIRICHLET条件)を与え変位を計算することにします。
FEMと立体求積法の設定
変数mにセットされたメッシュからそれぞれ変位応答格納用の変数と各節点データ操作用の変数を作成しておきます。GetFEM++の特徴として積分法の選択肢の多さがあります。
Step6:  mfuとmfdにはそれぞれLagrange要素$Q_3$と$Q_1$が入ります。$Q_3$は変位用、$Q_1$はMises応力などのデータ用に準備したものです。
FEM手法として古典的なLagrange要素$P_k$を割り当てます。
| degree               | dimension            | d.o.f. number                        | class        | vectorial                | $\tau$-equivalent        | Polynomial   |
|
Step7:  立体求積法には15積分点・5次のtetrahedronを使用します。
<img src="./image/getfemlistintmethodtetrahedron5.png">
Step8:  領域の定義
次に領域を定義し、境界条件を設定する布石とします。今回は三脚の上端にNEUMANN条件を、下端にDIRICHLET条件を設定します。
大き流れとしては以下の通りです。
節点の座標を取得
節点が領域に属しているかの真偽表の配列を作成
領域に属している節点のみ節点番号を抜き出す
節点番号から要素の面の情報を抜き出す
要素番号と面番号の情報を元に領域に番号をつける
節点の座標を取得
メッシュmの節点の座標を出力するためには、m.pts()を使用します。この関数により行方向がXYZ軸、列方向が各節点を意味している座標の配列が得られます。
Step9:  Gmshで先ほど出力したメッシュを確認してみるとY座標が最大になっている部分に三脚の上端が、Y座標が最小になっている部分に三脚の下端があることがわかります。それぞれの座標をPで確認しておきましょう。
Step10:  節点が領域に属しているかの真偽表の配列を作成
Y座標上の最大値と最小値がわかりましたので、その部分に節点が属しているかの真偽表の配列を作ります。頂部(ctop)は最大値$13$との差が$10^{-6}$の場合に節点が領域に含まれているとみなし、下端は最小値$-10$との差が$10^{-6}$の場合に領域に属しているとみなします。
Step11:  領域に属している節点のみ節点番号を抜き出す
次に領域に属している節点のみ節点番号を抜き出します。意図した操作はNumpyのcompress関数とPythonのrange関数で可能です。ちなみに、m.nbptsはメッシュmの節点数(Number of Points)を意味しています。
Step12:  節点番号から要素の面の情報を抜き出す
次に節点番号から要素の面を格納した配列を作成します。この操作はメッシュmに対してm.faces_from_pid()の関数で操作可能です。この配列には１行目に要素番号、２行目に要素内での面番号を格納しています。
Step13:  要素番号と面番号の情報を元に領域に番号をつける
最後に、面の情報を元にメッシュmに境界領域を作成し付番します。使用する関数はm.set_regionです。
Step14:  外力ベクトルと剛性マトリックスの組み立て
以上で設定した条件を使用して外力ベクトルと剛性行列を設定します。
Step15:  固定条件の設定
DIRICHLET条件(固定端条件)の設定をします。DIRICHLET条件は以下のような式で表される境界条件です。
$HU = R$
固定条件を与える場合は対象領域の節点に
$\left[\begin{array}{ccc}
1\
 & 1\
 &  & 1
\end{array}\right]\left{ \begin{array}{c}
U_{x}\
U_{y}\
U_{z}
\end{array}\right} =\left{ \begin{array}{c}
0\
0\
0
\end{array}\right} $
の関係を与える必要があります。すなわち、対象領域のみに
$H=\left[\begin{array}{ccc}
1\
 & 1\
 &  & 1
\end{array}\right]$、$R=\left{ \begin{array}{c}
0\
0\
0
\end{array}\right} $
になる行列を与えてれる関数が必要です。それが以下のgf.asm_dirichlet関数です。
Step16:  さて、Getfem++では、はじめにご説明した連立方程式の問題を以下のように言い換えます。
$KU = F$ ただし、$HU = R$
$\Downarrow$
$(N^TKN)U^ = N^TF$ ただし、$U = NU^ + U_0$
この$N$は例えば条件が固定端の場合にはその自由度の行と列を削除する効果があります。この$N$と$U_0$は次の関数で計算できます。
Step17:  変位の計算
以上で計算するべき連立方程式の剛性行列と外力ベクトルが定まりましたので、連立方程式を解きます。Getfem++のモジュールで剛性行列からプリコンディショナーを用意し、CG法を利用して変位を計算します。
Step18:  ポスト処理
以上で計算した変位をポスト処理結果として出力します。まずスライスモジュールを使用して、メッシュの外側の部分のみ抜き出します。次にスライスオブジェクトに変位の計算結果と凡例を渡して、posファイルに出力します。これで、'Displacement'と凡例のあるposファイルが'./pos/m_result.pos'というファイルに出力されたことになります。
<ASSISTANT_TASK:>
Python Code:
import getfem as gf
import numpy as np
Explanation: はじめに
静的な場合の有限要素法問題における変位は以下のような連立方程式をとくことにより得られます。
$KU=F$ 
$HU=R$
ここで、$K$は剛性行列、$U$は変位ベクトルです。
一方、$F$は外力ベクトル(NEUMANN条件)であり、$HU=R$は固定端やローラー条件(DIRICHLET条件)です。
この資料は三脚のメッシュをpythonでインポートし、境界条件を設定し剛性行列$K$と外力ベクトル$F$を作成して変位$U$を計算する過程を説明したものです。
インストールとモジュールのインポート
インストール方法
Linux(Ubuntu)
まずはGetfem++のpythonインターフェースのモジュールをインポートします。Ubuntu系のOSではインストールコマンドは
$ sudo apt-get update; sudo apt-get install python-getfem++
です。さらに、境界条件設定に使うためnumpyをインポートします。numpyはpython-getfem++と依存関係にあるため、自動的にインストールされます。インストールが終わったら、結果出力とGUI表示のため、gmshとipython-notebookをインストールしましょう。
$ sudo apt-get install gmsh python-pip; sudo pip install notebook
このファイルがあるディレクトリに移動
$ ipython notebook
で、Ipython Notebookを起動した後このファイルを開いてください。
Windows
Python 2.7系の32bit版をインストール(Getfem++のWindows版が32bit版しかないので、32bit版でないと動かない)
https://www.python.org/ftp/python/2.7.10/python-2.7.10.msi
以下のファイルで、Python2.7系のGetfem++モジュールをインストールします。
http://download.gna.org/getfem/misc/getfem_python-4.1.win32-py2.7.exe
あとは、パスを通したpipコマンドで、numpyとipythonをインストールします。
管理者権限のコマンドライン
pip install notebook
pip install numpy
numpyのインストールはCコンパイラが必要になるので、メッセージをよく読みインストールをしてください。
モジュールのインポート
慣例的にGetfem++はgfとしてインポートし、numpyはnpとしてインポートします。
End of explanation
file_msh = './mesh/tripod.mesh'
E = 1e3
Nu = 0.3
Lambda = E*Nu/((1+Nu)*(1-2*Nu))
Mu = E/(2*(1+Nu))
Explanation: パラメータの設定
まずはパラメータを設定します。この例では、このファイルからの相対パスが './mesh/tripod.mesh'であるメッシュファイルを使用します。メッシュは四面体要素です。ヤング係数$E = 1.0\times10^3$、ポアソン比$\nu = 0.3$とします。
End of explanation
m = gf.Mesh('load',file_msh)
m.set('optimize_structure')
Explanation: メッシュの読み込み
メッシュを外部ファイルから読み込みます。メッシュの入力データはGetfem++の独自のフォーマットになっていますが、解読は容易だと思います。
End of explanation
m.export_to_pos('./pos/m.pos')
Explanation: 次のコマンドで、mに設定したメッシュをgmshのスクリプトでpng画像に打ち出し確認します。
End of explanation
%%writefile gscript
Print "./png/m.png";
Exit;
!gmsh ./pos/m.pos gscript
from IPython.core.display import Image
Image('./png/m.png')
Explanation: 以下のようにipythonの機能を使用して、システムのコマンドラインでの操作を行いIpython上に画像を表示してみます。もちろん、gmshで直接posファイルを表示することもできます。
End of explanation
mfu = gf.MeshFem(m,3) # displacement
mfd = gf.MeshFem(m,1) # data
Explanation: これは横から見たメッシュ図です。三脚のメッシュであることがわかります。今回はこの頂点部分に荷重(NEUMANN条件)、下端に固定条件(DIRICHLET条件)を与え変位を計算することにします。
FEMと立体求積法の設定
変数mにセットされたメッシュからそれぞれ変位応答格納用の変数と各節点データ操作用の変数を作成しておきます。GetFEM++の特徴として積分法の選択肢の多さがあります。
End of explanation
mfu.set_fem(gf.Fem('FEM_PK(3,1)'))
mfd.set_fem(gf.Fem('FEM_PK(3,0)'))
Explanation: mfuとmfdにはそれぞれLagrange要素$Q_3$と$Q_1$が入ります。$Q_3$は変位用、$Q_1$はMises応力などのデータ用に準備したものです。
FEM手法として古典的なLagrange要素$P_k$を割り当てます。
| degree               | dimension            | d.o.f. number                        | class        | vectorial                | $\tau$-equivalent        | Polynomial   |
|:--------------------:|:--------------------:|:------------------------------------:|:------------:|:------------------------:|:------------------------:|:------------:|
| $K,0\leq K\leq255$   | $P,1\leq K\leq255$   | $\dfrac{\left(K+P\right)!}{K!P!}$    | $C^0$        | No$\left(Q=1\right)$     | Yes$\left(M=Id\right)$   | Yes          |
End of explanation
mim = gf.MeshIm(m,gf.Integ('IM_TETRAHEDRON(5)'))
Explanation: 立体求積法には15積分点・5次のtetrahedronを使用します。
<img src="./image/getfemlistintmethodtetrahedron5.png">
End of explanation
P = m.pts()
print P
Explanation: 領域の定義
次に領域を定義し、境界条件を設定する布石とします。今回は三脚の上端にNEUMANN条件を、下端にDIRICHLET条件を設定します。
大き流れとしては以下の通りです。
節点の座標を取得
節点が領域に属しているかの真偽表の配列を作成
領域に属している節点のみ節点番号を抜き出す
節点番号から要素の面の情報を抜き出す
要素番号と面番号の情報を元に領域に番号をつける
節点の座標を取得
メッシュmの節点の座標を出力するためには、m.pts()を使用します。この関数により行方向がXYZ軸、列方向が各節点を意味している座標の配列が得られます。
End of explanation
print 'Ymax = ', P[1,:].max()
print 'Ymin = ', P[1,:].min()
Explanation: Gmshで先ほど出力したメッシュを確認してみるとY座標が最大になっている部分に三脚の上端が、Y座標が最小になっている部分に三脚の下端があることがわかります。それぞれの座標をPで確認しておきましょう。
End of explanation
ctop = (abs(P[1,:] - 13) < 1e-6)
print 'ctop = ', ctop
cbot = (abs(P[1,:] + 10) < 1e-6)
print 'cbot = ', cbot
Explanation: 節点が領域に属しているかの真偽表の配列を作成
Y座標上の最大値と最小値がわかりましたので、その部分に節点が属しているかの真偽表の配列を作ります。頂部(ctop)は最大値$13$との差が$10^{-6}$の場合に節点が領域に含まれているとみなし、下端は最小値$-10$との差が$10^{-6}$の場合に領域に属しているとみなします。
End of explanation
pidtop = np.compress(ctop,range(0,m.nbpts()))
pidbot = np.compress(cbot,range(0,m.nbpts()))
print 'pidtop = ', pidtop
print 'pidbot = ', pidbot
Explanation: 領域に属している節点のみ節点番号を抜き出す
次に領域に属している節点のみ節点番号を抜き出します。意図した操作はNumpyのcompress関数とPythonのrange関数で可能です。ちなみに、m.nbptsはメッシュmの節点数(Number of Points)を意味しています。
End of explanation
ftop = m.faces_from_pid(pidtop)
fbot = m.faces_from_pid(pidbot)
print 'ftop = ', ftop
print 'fbot = ', fbot
Explanation: 節点番号から要素の面の情報を抜き出す
次に節点番号から要素の面を格納した配列を作成します。この操作はメッシュmに対してm.faces_from_pid()の関数で操作可能です。この配列には１行目に要素番号、２行目に要素内での面番号を格納しています。
End of explanation
NEUMANN_BOUNDARY = 1
DIRICHLET_BOUNDARY = 2
m.set_region(NEUMANN_BOUNDARY,ftop)
m.set_region(DIRICHLET_BOUNDARY,fbot)
Explanation: 要素番号と面番号の情報を元に領域に番号をつける
最後に、面の情報を元にメッシュmに境界領域を作成し付番します。使用する関数はm.set_regionです。
End of explanation
nbd = mfd.nbdof()
F = gf.asm_boundary_source(NEUMANN_BOUNDARY, mim, mfu, mfd, np.repeat([[0],[-100],[0]],nbd,1))
K = gf.asm_linear_elasticity(mim, mfu, mfd, np.repeat([Lambda], nbd), np.repeat([Mu], nbd))
sl = gf.Slice(('boundary',), mfu, 1)
sl.export_to_pos('./pos/F.pos',mfu,F,'F')
%%writefile gscript
Print "./png/F.png";
Exit;
!gmsh ./pos/F.pos gscript
from IPython.core.display import Image
Image('./png/F.png')
Explanation: 外力ベクトルと剛性マトリックスの組み立て
以上で設定した条件を使用して外力ベクトルと剛性行列を設定します。
End of explanation
(H,R) = gf.asm_dirichlet(DIRICHLET_BOUNDARY, mim, mfu, mfd, mfd.eval('[[1,0,0],[0,1,0],[0,0,1]]'), mfd.eval('[0,0,0]'))
Explanation: 固定条件の設定
DIRICHLET条件(固定端条件)の設定をします。DIRICHLET条件は以下のような式で表される境界条件です。
$HU = R$
固定条件を与える場合は対象領域の節点に
$\left[\begin{array}{ccc}
1\
 & 1\
 &  & 1
\end{array}\right]\left{ \begin{array}{c}
U_{x}\
U_{y}\
U_{z}
\end{array}\right} =\left{ \begin{array}{c}
0\
0\
0
\end{array}\right} $
の関係を与える必要があります。すなわち、対象領域のみに
$H=\left[\begin{array}{ccc}
1\
 & 1\
 &  & 1
\end{array}\right]$、$R=\left{ \begin{array}{c}
0\
0\
0
\end{array}\right} $
になる行列を与えてれる関数が必要です。それが以下のgf.asm_dirichlet関数です。
End of explanation
(N,U0) = H.dirichlet_nullspace(R)
Nt = gf.Spmat('copy',N)
Nt.transpose()
KK = Nt*K*N
FF = Nt*(F-K*U0)
Explanation: さて、Getfem++では、はじめにご説明した連立方程式の問題を以下のように言い換えます。
$KU = F$ ただし、$HU = R$
$\Downarrow$
$(N^TKN)U^ = N^TF$ ただし、$U = NU^ + U_0$
この$N$は例えば条件が固定端の場合にはその自由度の行と列を削除する効果があります。この$N$と$U_0$は次の関数で計算できます。
End of explanation
P = gf.Precond('ildlt',KK)
UU = gf.linsolve_cg(KK,FF,P)
U = N*UU+U0
Explanation: 変位の計算
以上で計算するべき連立方程式の剛性行列と外力ベクトルが定まりましたので、連立方程式を解きます。Getfem++のモジュールで剛性行列からプリコンディショナーを用意し、CG法を利用して変位を計算します。
End of explanation
sl = gf.Slice(('boundary',), mfu, 1)
sl.export_to_pos('./pos/m_result.pos', mfu, U, 'Displacement')
%%writefile gscript
Print "./png/m_result.png";
Exit;
!gmsh ./pos/m_result.pos gscript
from IPython.core.display import Image
Image('./png/m_result.png')
Explanation: ポスト処理
以上で計算した変位をポスト処理結果として出力します。まずスライスモジュールを使用して、メッシュの外側の部分のみ抜き出します。次にスライスオブジェクトに変位の計算結果と凡例を渡して、posファイルに出力します。これで、'Displacement'と凡例のあるposファイルが'./pos/m_result.pos'というファイルに出力されたことになります。
End of explanation
<END_TASK>