まずは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を起動した後このファイルを開いてください。
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としてインポートします。
In [1]:
import getfem as gf
import numpy as np
In [2]:
file_msh = './mesh/tripod.mesh'
E = 1e3
Nu = 0.3
Lambda = E*Nu/((1+Nu)*(1-2*Nu))
Mu = E/(2*(1+Nu))
In [3]:
m = gf.Mesh('load',file_msh)
m.set('optimize_structure')
次のコマンドで、mに設定したメッシュをgmshのスクリプトでpng画像に打ち出し確認します。
In [4]:
m.export_to_pos('./pos/m.pos')
以下のようにipythonの機能を使用して、システムのコマンドラインでの操作を行いIpython上に画像を表示してみます。もちろん、gmshで直接posファイルを表示することもできます。
In [5]:
%%writefile gscript
Print "./png/m.png";
Exit;
In [6]:
!gmsh ./pos/m.pos gscript
from IPython.core.display import Image
Image('./png/m.png')
Out[6]:
これは横から見たメッシュ図です。三脚のメッシュであることがわかります。今回はこの頂点部分に荷重(NEUMANN条件)、下端に固定条件(DIRICHLET条件)を与え変位を計算することにします。
In [7]:
mfu = gf.MeshFem(m,3) # displacement
mfd = gf.MeshFem(m,1) # data
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 |
In [8]:
mfu.set_fem(gf.Fem('FEM_PK(3,1)'))
mfd.set_fem(gf.Fem('FEM_PK(3,0)'))
立体求積法には15積分点・5次のtetrahedronを使用します。
In [9]:
mim = gf.MeshIm(m,gf.Integ('IM_TETRAHEDRON(5)'))
In [10]:
P = m.pts()
print P
Gmshで先ほど出力したメッシュを確認してみるとY座標が最大になっている部分に三脚の上端が、Y座標が最小になっている部分に三脚の下端があることがわかります。それぞれの座標をPで確認しておきましょう。
In [11]:
print 'Ymax = ', P[1,:].max()
print 'Ymin = ', P[1,:].min()
In [12]:
ctop = (abs(P[1,:] - 13) < 1e-6)
print 'ctop = ', ctop
cbot = (abs(P[1,:] + 10) < 1e-6)
print 'cbot = ', cbot
In [13]:
pidtop = np.compress(ctop,range(0,m.nbpts()))
pidbot = np.compress(cbot,range(0,m.nbpts()))
print 'pidtop = ', pidtop
print 'pidbot = ', pidbot
In [14]:
ftop = m.faces_from_pid(pidtop)
fbot = m.faces_from_pid(pidbot)
print 'ftop = ', ftop
print 'fbot = ', fbot
In [15]:
NEUMANN_BOUNDARY = 1
DIRICHLET_BOUNDARY = 2
m.set_region(NEUMANN_BOUNDARY,ftop)
m.set_region(DIRICHLET_BOUNDARY,fbot)
In [16]:
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')
In [17]:
%%writefile gscript
Print "./png/F.png";
Exit;
In [18]:
!gmsh ./pos/F.pos gscript
from IPython.core.display import Image
Image('./png/F.png')
Out[18]:
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関数です。
In [19]:
(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]'))
さて、Getfem++では、はじめにご説明した連立方程式の問題を以下のように言い換えます。
$KU = F$ ただし、$HU = R$
$\Downarrow$
$(N^TKN)U^* = N^TF$ ただし、$U = NU^* + U_0$
この$N$は例えば条件が固定端の場合にはその自由度の行と列を削除する効果があります。この$N$と$U_0$は次の関数で計算できます。
In [20]:
(N,U0) = H.dirichlet_nullspace(R)
Nt = gf.Spmat('copy',N)
Nt.transpose()
KK = Nt*K*N
FF = Nt*(F-K*U0)
In [21]:
P = gf.Precond('ildlt',KK)
UU = gf.linsolve_cg(KK,FF,P)
U = N*UU+U0
In [22]:
sl = gf.Slice(('boundary',), mfu, 1)
sl.export_to_pos('./pos/m_result.pos', mfu, U, 'Displacement')
In [23]:
%%writefile gscript
Print "./png/m_result.png";
Exit;
In [24]:
!gmsh ./pos/m_result.pos gscript
from IPython.core.display import Image
Image('./png/m_result.png')
Out[24]:
In [ ]: