PythonにはNumpy/Matplotlib/Scipy/Getfem++などのパッケージがあり、それらを組み合わせることにより容易に有限要素法のスクリプトを組むことができます。それぞれのパッケージに関しては各Wikipediaの内容を確認してください。
Numpy (https://ja.wikipedia.org/wiki/NumPy)
Matplotlib(https://ja.wikipedia.org/wiki/Matplotlib)
Scipy (https://ja.wikipedia.org/wiki/SciPy)
Getfem++ (https://ja.wikipedia.org/wiki/Getfem%2B%2B)
それぞれのパッケージはPythonでは以下のように呼び出します。
In [ ]:
import numpy as np
これで、これらのパッケージはそれぞれ、np/sp/gf/pltとしてインポートされました。
各パッケージの機能は
In [ ]:
np.
In [ ]:
np?
本節ではNumpyのチュートリアル
https://docs.scipy.org/doc/numpy-dev/user/quickstart.html
の内容を解説します。
numpyを使用する利点としては多次元配列の操作が容易になることがあげられます。 numpyの配列のクラスはndarrayと呼ばれています。ndarrayの主な属性には主に以下のものがあります。
配列の次元
配列の形状。n行m列の行列については、 形状 （M、N）となります 。
配列の要素数の合計。
配列の要素のタイプ。numpy.int32、numpy.int16、およびnumpy.float64などがあります。
配列の各要素のサイズ（バイト単位）。
配列の実際の要素を含むバッファー。 インデックスで配列の要素にアクセスするため、この属性を使用する必要はありません。
In [ ]:
b = np.array([6, 7, 8])
A = np.array([[1, 2, 0],[2, 5, 6],[0, 6, 9]])
# bを出力します
print "b ="
print b
# Aを出力します
print "A ="
print A
np.arrayについて詳細を確認したい場合は次のようにコマンドを入力します。
In [ ]:
np.array?
このコマンドを入力するとコマンドの呼び出し方法・目的・パラメータの説明が出てきます。 詳細説明から分かるようにarrayには先ほど指定したobject以外にも引数があります。この引数のなかで"引数="となっているものはデフォルト引数と呼ばれ指定されない場合は説明に書かれた値が使用されます。
デフォルト引数を変更したい場合には以下のように指定します。
In [ ]:
A = np.array([[1, 2, 0],[2, 5, 6],[0, 6, 9]], dtype=complex)
print A
dtypeは詳細説明にもあるとおりarrayのデータ型を指定する引数です。これを"dtype=complex"としているため今回は複素数の行列が作成されています。
ipoython notebookで2次元の図を出力する場合にはmatplotlibが使われます。
https://ja.wikipedia.org/wiki/Matplotlib
以下のようにモジュールを呼び出します。２行目の"%matplotlib inline"はmatplotlibの図をiPython Notebook内で描写したグラフを展開するためのものです。
In [ ]:
from matplotlib import pylab as plt
%matplotlib inline
In [ ]:
a = np.array([1.0, 2.0])
b = np.array([3.0, 4.0])
plt.plot(a, b)
In [ ]:
np.linalg?
例えば、先ほど作成した行列$A$の逆行列を計算したい場合はinv関数を使用します。
In [ ]:
np.linalg.inv?
In [ ]:
AA = np.linalg.inv(A)
print AA
Aに対して行列積(np.dot)を行うと以下のように単位行列が作成されます。これにより正しく逆行列が計算されることがわかります。
In [ ]:
np.dot?
In [ ]:
np.dot(A,AA)
また、逆行列だけでなく固有値解析の実行も可能です。
In [ ]:
a, V = np.linalg.eig(A)
# 固有値を出力します
print "固有値="
print a
# 固有値ベクトル
print "固有ベクトル="
print V
plt.plot(V[:,0])
plt.plot(V[:,1])
ScipyはNumPyを基盤にした科学計算ライブラリです。本節ではScipyのチュートリアル
http://docs.scipy.org/doc/scipy/reference/tutorial/
の内容のうち、疎行列による線形代数関係のライブラリの説明をします。
In [ ]:
import numpy as np
import scipy as sp
from scipy import sparse
In [ ]:
sparse?
Scipyで疎行列を操作する際にはNumpyで使用していたarrayオブジェクトとは別のオブジェクトを使用します。疎行列に使用されているオブジェクトは複数あり、相互に変換可能ですが今回はcoo_matrixオブジェクトを紹介します。
In [ ]:
# A sparse matrix in COOrdinate format
sparse.coo_matrix?
マニュアルにある通り疎行列に代入するための値とそのインデックスを準備します。
In [ ]:
# 行のインデックス
I = np.array([0,0,1,1,1,2,2])
# 列のインデックス
J = np.array([0,1,0,1,2,1,2])
# 値
V = np.array([1,2,2,5,6,6,9])
A = sparse.coo_matrix((V,(I,J)),shape=(3,3))
print A
以上のように各インデックスに値の入った疎行列が作成されます。このようにオブジェクトを作成した際には必ずなんのオブジェクトであるかを確認しておきましょう。
In [ ]:
A?
また、これらの疎行列はファイルへの保存と読み込みが可能です。その際にはScipyのioモジュールを使用します。
In [ ]:
from scipy import io
In [ ]:
io?
ここではMatrix Market形式(http://math.nist.gov/MatrixMarket/formats.html) で書き出しを行ってみます。
In [ ]:
io.mmwrite("A.mtx", A)
ファイルの中身を出力してみるとMatrix Market形式で正しく出力されていることがわかります。Matrix Market形式は国際的に使われている疎行列のファイル形式です。
In [ ]:
%cat A.mtx
また、既存のファイルを読み込んで新しく行列を作成することも可能です。
In [ ]:
B = io.mmread("A.mtx")
以下のようにしてファイルから行列が読み込めていることを確認できます。
In [ ]:
print B
以上では、Pythonによる密行列の線形代数までを紹介してきました。しかしながら、有限要素法計算の際にはメモリ節約のため疎行列計算が行われます。その際にはscipyのsparseパッケージとその下のlinalgパッケージを使用します。
In [ ]:
import scipy as sp
from scipy.sparse import linalg
In [ ]:
linalg?
密行列と同じ操作を疎行列の線形代数パッケージで実行します。
In [ ]:
A = A.tocsc()
AA = linalg.inv(A)
print AA
密行列の場合は行列積が"*"になりますので注意してください。先ほどと同じように単位行列が逆行列との積で得られます。
In [ ]:
print A*AA
また、固有値解析も少し複雑ですが以下のように計算できます。
In [ ]:
X = sp.rand(A.shape[0], 3)
linalg.lobpcg(A,X)
密行列の場合と同じ計算結果になることが確認できます。 それぞれの関数の説明は以下で確認してください。
In [ ]:
sp.rand?
In [ ]:
linalg.lobpcg?
In [ ]:
import getfem as gf
import numpy as np
from scipy import io
本節ではメッシュ作成とそれを離散化した行列を作成するパッケージGetfem++について説明します。必要となるオブジェクトは以下の図に示すものになります。
有限要素で何を使用するかを指定する。(ラグランジュ要素Pk、Qkなど)
メッシュの指定(point:節点とconvex:要素などから作られる)
積分方法を指定する。
メッシュの要素に対してFemを設定したオブジェクト、自由度数などの情報も設定する。
メッシュの要素に対して積分方法を指定したオブジェクト
離散化した行列組み立てのためのオブジェクト。MeshFemとMeshImを基に作成する。
モデル全体の情報を設定するオブジェクト。今回は未使用。
マニュアルからわかるように、規則的なメッシュを作成したり、Gmshのメッシュファイルなどからインポートすることも可能です。ただ、今回はメッシュオブジェクトに慣れることを目的としますので、まずは１次元の空のMeshオブジェクトを作成します。
In [ ]:
gf.Mesh?
In [ ]:
m = gf.Mesh('empty', 1)
ご覧のようにこの時点ではメッシュは空です。このメッシュに節点と要素を追加していきます。
In [ ]:
print m
要素を追加していきます。(add_convexメソッド)
In [ ]:
m.add_convex?
GeoTransが必要であることが分かりましたので、GeoTransのマニュアルも確認します。
In [ ]:
gf.GeoTrans?
座標 x = 0.0からx = 2.0の間に平行６面体要素を追加します。
In [ ]:
m.add_convex?
In [ ]:
m.add_convex(gf.GeoTrans('GT_QK(1,1)'),[[0,2]])
In [ ]:
m.add_convex(gf.GeoTrans('GT_QK(1,1)'),[[2,4]])
メッシュmに２つの要素が追加されました。
In [ ]:
print m
ちなみにこのMeshオブジェクトはVTKファイルにExportできます。 1次元のメッシュでは寂しいので、２次元の方眼メッシュを作成しましょう。
In [ ]:
m2 = gf.Mesh('cartesian', np.array([0.0, 1.0, 2.0]), np.array([0.0, 1.0, 2.0]))
In [ ]:
m2.export_to_vtk?
In [ ]:
m2.export_to_vtk('m2.vtk','ascii')
これをParaviewで確認すると4要素の方眼メッシュが作成されていることが分かります。
次に、先ほどの１次元のメッシュを使用して剛性マトリックスと質量マトリックスを作成してみます。
In [ ]:
gf.asm_linear_elasticity?
In [ ]:
gf.MeshFem?
変位用のMeshFemとデータ用のMeshFemを作成します。MeshFemはMesh上にFem要素を割り当てるためのオブジェクトです。
In [ ]:
mfu = gf.MeshFem(m,1) # 変位用MeshFem
mfd = gf.MeshFem(m,1) # データ用MeshFem
今回は１次元のLagrange要素$Q_1$を使用します。
In [ ]:
gf.Fem?
In [ ]:
fem = gf.Fem('FEM_QK(1,1)')
作成したFem要素をMeshFemに設定します。
In [ ]:
mfu.set_fem(fem)
mfd.set_fem(fem)
次は、積分要素を設定します。
In [ ]:
gf.Integ?
今回は１次のガウス積分を使用します。
In [ ]:
im = gf.Integ('IM_GAUSS1D(1)')
In [ ]:
mim = gf.MeshIm(m, im)
In [ ]:
gf.asm_linear_elasticity?
In [ ]:
gf.asm_mass_matrix?
In [ ]:
# ラメの弾性定数
Lambda = np.array([1.0, 1.0, 1.0])
Mu = np.array([1.0, 1.0, 1.0])
# 剛性行列の作成
K = gf.asm_linear_elasticity(mim, mfu, mfd, Lambda, Mu)
print K.full() # 剛性行列の出力
# 質量行列の作成
M = gf.asm_mass_matrix(mim, mfu)
print M.full() # 質量行列の出力
以上により作成した剛性マトリックスと質量マトリックスはSpmatというGetfem++独自の行列オブジェクトで、先ほど紹介したMatrixMarketのフォーマットに出力するメソッドも備えています。
また、自分で微分方程式を定義して計算することも可能です。その際にはgf.asm_volumicを使用します。
In [ ]:
gf.asm_volumic?
線形等方弾性問題の微分方程式は次のように表せます。
$\nabla\cdot\left[\left\{ \nabla\boldsymbol{U}+\left(\nabla\boldsymbol{U}\right)^{T}\right\}\mu + tr\left(\nabla\boldsymbol{U}\right)\boldsymbol{I}\lambda\right]=\boldsymbol{F}$
Getfem++ではテンソルを使用した表現で以下のように微分方程式を定義できます。
In [ ]:
K = gf.asm_volumic("lambda=data$1(#2);"+ \
"mu=data$2(#2);"+ \
"t=comp(vGrad(#1).vGrad(#1).Base(#2));"+ \
"M(#1,#1)+= sym(t(:,i,j,:,i,j,k).mu(k)+ t(:,j,i,:,i,j,k).mu(k)+ t(:,i,i,:,j,j,k).lambda(k))",
mim, mfu, mfd, np.array([1.0, 1.0, 1.0]), np.array([1.0, 1.0, 1.0]))
当然ながら結果は同じになります。以下のようにすると"小さい"行列であれば先ほどのNumpyのarrayオブジェクトが出力されます。
In [ ]:
K.full()
もちろん、構造用のNavierの式の離散化だけではありません。laplacianやhelmholtzなどの関数があります。 つまり、１つのメッシュでさまざまな微分方程式の離散化が可能ということです。
In [ ]:
gf.asm_laplacian?
In [ ]:
gf.asm_helmholtz?
さて、この得られた行列はScipyの疎行列というわけではなく、SpmatというGetfem++(Gmm++)独自の形式です。
In [ ]:
K?
Matrix-Marketフォーマットで出力してみます。
In [ ]:
K.save('mm', "K.mtx")
M.save('mm', "M.mtx")
出力されたファイルはノートブックと同じディレクトリに保存されています。
In [ ]:
K = io.mmread("K.mtx")
M = io.mmread("M.mtx")
In [ ]:
print "K = "
print K
print "M = "
print M
このようにScipyの疎行列に容易に変換できます。そのため、Scipyを使用した静解析や固有値解析を行うことにより有限要素法の解析が可能となります。 さて１次元メッシュで次のような階がnp.arrayオブジェクトで得られたとします。
In [ ]:
outou = np.array([1.0, 2.0, 3.0])
print outou
これをベクトル情報付きのVTKファイルとして出力します。その際にはGetfem++のSliceオブジェクトを使用します。
In [ ]:
gf.Slice?
スライスオブジェクトを作成する際にはメッシュをスライスすることで断面などの表示か可能です。 しかし、今回は特に断面を見る必要もないので何もしません。
In [ ]:
sl = gf.Slice(('none',), m, 1)
スライスオブジェクトにはVTKにMeshFemオブジェクトと応答値であるarrayオブジェクトを引数としてVTKに出力する機能があります。
In [ ]:
sl.export_to_vtk?
In [ ]:
sl.export_to_vtk('m_results.vtk', 'ascii', mfu, outou, 'jiyuunanamae')
In [ ]:
"""
Demonstrates computation of gradient with matplotlib.tri.CubicTriInterpolator.
"""
from matplotlib.tri import Triangulation, UniformTriRefiner,\
CubicTriInterpolator
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import numpy as np
import math
%matplotlib inline
#-----------------------------------------------------------------------------
# Electrical potential of a dipole
#-----------------------------------------------------------------------------
def dipole_potential(x, y):
""" The electric dipole potential V """
r_sq = x**2 + y**2
theta = np.arctan2(y, x)
z = np.cos(theta)/r_sq
return (np.max(z) - z) / (np.max(z) - np.min(z))
#-----------------------------------------------------------------------------
# Creating a Triangulation
#-----------------------------------------------------------------------------
# First create the x and y coordinates of the points.
n_angles = 30
n_radii = 10
min_radius = 0.2
radii = np.linspace(min_radius, 0.95, n_radii)
angles = np.linspace(0, 2*math.pi, n_angles, endpoint=False)
angles = np.repeat(angles[..., np.newaxis], n_radii, axis=1)
angles[:, 1::2] += math.pi/n_angles
x = (radii*np.cos(angles)).flatten()
y = (radii*np.sin(angles)).flatten()
V = dipole_potential(x, y)
# Create the Triangulation; no triangles specified so Delaunay triangulation
# created.
triang = Triangulation(x, y)
# Mask off unwanted triangles.
xmid = x[triang.triangles].mean(axis=1)
ymid = y[triang.triangles].mean(axis=1)
mask = np.where(xmid*xmid + ymid*ymid < min_radius*min_radius, 1, 0)
triang.set_mask(mask)
#-----------------------------------------------------------------------------
# Refine data - interpolates the electrical potential V
#-----------------------------------------------------------------------------
refiner = UniformTriRefiner(triang)
tri_refi, z_test_refi = refiner.refine_field(V, subdiv=3)
#-----------------------------------------------------------------------------
# Computes the electrical field (Ex, Ey) as gradient of electrical potential
#-----------------------------------------------------------------------------
tci = CubicTriInterpolator(triang, -V)
# Gradient requested here at the mesh nodes but could be anywhere else:
(Ex, Ey) = tci.gradient(triang.x, triang.y)
E_norm = np.sqrt(Ex**2 + Ey**2)
#-----------------------------------------------------------------------------
# Plot the triangulation, the potential iso-contours and the vector field
#-----------------------------------------------------------------------------
plt.figure()
plt.gca().set_aspect('equal')
plt.triplot(triang, color='0.8')
levels = np.arange(0., 1., 0.01)
cmap = cm.get_cmap(name='hot', lut=None)
plt.tricontour(tri_refi, z_test_refi, levels=levels, cmap=cmap,
linewidths=[2.0, 1.0, 1.0, 1.0])
# Plots direction of the electrical vector field
plt.quiver(triang.x, triang.y, Ex/E_norm, Ey/E_norm,
units='xy', scale=10., zorder=3, color='blue',
width=0.007, headwidth=3., headlength=4.)
plt.title('Gradient plot: an electrical dipole')
plt.show()
In [ ]:
from tempfile import NamedTemporaryFile
from IPython.display import HTML
VIDEO_TAG = """<video controls>
<source src="data:video/x-m4v;base64,{0}" type="video/mp4">
Your browser does not support the video tag.
</video>"""
def anim_to_html(anim):
if not hasattr(anim, '_encoded_video'):
with NamedTemporaryFile(suffix='.mp4') as f:
anim.save(f.name, fps=20, extra_args=['-vcodec', 'libx264'])
video = open(f.name, "rb").read()
anim._encoded_video = video.encode("base64")
return VIDEO_TAG.format(anim._encoded_video)
def display_animation(anim):
plt.close(anim._fig)
return HTML(anim_to_html(anim))
In [ ]:
import pylab as plt
import numpy as np
from matplotlib import animation
# First set up the figure, the axis, and the plot element we want to animate
fig = plt.figure()
ax = plt.axes(xlim=(0, 2), ylim=(-2, 2))
line, = ax.plot([], [], lw=2)
# initialization function: plot the background of each frame
def init():
line.set_data([], [])
return line,
# animation function. This is called sequentially
def animate(i):
x = np.linspace(0, 2, 1000)
y = np.sin(2 * np.pi * (x - 0.01 * i))
line.set_data(x, y)
return line,
# call the animator. blit=True means only re-draw the parts that have changed.
anim = animation.FuncAnimation(fig, animate, init_func=init,
frames=100, interval=20, blit=True)
# call our new function to display the animation
display_animation(anim)
In [ ]: