Pythonを使用したFEM解析の概要

Binder

このチュートリアルでは、Pythonを使用して有限要素法解析の前処理、求解、および後処理を行う方法を説明します。これは、前処理と求解に GetFEM と呼ばれるPythonインタフェースを持つ有限要素法ライブラリを使用します。vtkファイルを meshio を使ってロードし、前処理と後処理で matplotlib を使って視覚化します。このチュートリアルは PyConJP 2019 talk で使われました。以下のYouTubeでトークを見ることができます。このチュートリアルは、次の official GetFEM page tutorial に基づいています。

[1]:
from IPython.display import YouTubeVideo

YouTubeVideo("6JuB1GiDLQQ", start=512)
[1]:

インストール

Pythonインターフェースを含むGetFEMは aptitude update と aptitude install python3-getfem++ を実行することで端末からインストールできます。

sudo aptitude install python3-getfem++

このチュートリアルには、 requirements.txt の追加パッケージが必要です。これらの環境はすでに Dockerfile で設定されているため、構築する必要はありません。

問題の設定

問題は、Math Worksのホームページに掲載されている 「Poisson’s Equation on Unit Disk」 を参照しています。

\[−\Delta u=1 \ {\rm on}\ \Omega, u=0 \ {\rm on}\ \delta \Omega\]

|pdedemo1\_01|

GetEMの使用方法

有限要素問題を解くためにGetFEMを使用する場合、次の手順を実行します。GetFEMの使用の詳細については、 このページ を参照してください。

  • MesherObject を定義する
  • Mesh を定義する
  • MeshFem を定義する
  • MeshIm を定義する
  • Model を定義して設定する
  • Modelオブジェクトを解く
  • Modelオブジェクトから値を取得する

初期化

GetFEMは次のようにインポートできます(numpyもインポートしなければなりません)。

[2]:
import getfem as gf
import numpy as np

メッシュ生成

We use GetFEM’s MesherObject to create a mesh from the geometric information to be analyzed. This object represent a geometric object to be meshed by the experimental meshing procedure of GetFEM. We can represent a disk by specifying a radius, a 2D center, and using 「ball」 geometry.

[3]:
center = [0.0, 0.0]
radius = 1.0

mo = gf.MesherObject("ball", center, radius)

mo で表される幾何形状に対してGetFEMの実験メッシャを呼び出すことによりメッシュオブジェクト mesh を作成することができます。要素のおおよその直径は h で与えられ、メッシュの次数は K (曲線境界の \(K > 1\) )です。

[4]:
h = 0.1
K = 2
mesh = gf.Mesh("generate", mo, h, K)

境界の選択

境界条件を定義するために、円の外周に境界番号を設定します。

[5]:
outer_faces = mesh.outer_faces()
OUTER_BOUND = 1
mesh.set_region(OUTER_BOUND, outer_faces)

メッシュ描画

We visualize the created mesh to check its quality. We can output mesh objects, but matplotlib can only display triangles. Therefore, we make a Slice object with a slice operation of ("none",), which does not cut the mesh. With a refinement of 1, this serves to convert the mesh to triangles.

[6]:
sl = gf.Slice(("none",), mesh, 1)

export_to_vtkメソッドを使用すると、sliceをVTKファイルにエクスポートできます。

[7]:
sl.export_to_vtk("sl.vtk", "ascii")

VTKファイルはParaviewまたはmayavi 2を使用してレンダリングできます。今回はjupyter notebookに表示するために、meshioを読み込んでmatplotlibで描画します。

[8]:
import pyvista as pv
from pyvirtualdisplay import Display

display = Display(visible=0, size=(1280, 1024))
display.start()
p = pv.Plotter()
m = pv.read("sl.vtk")
p.add_mesh(m, show_edges=True)
pts = m.points
p.show(window_size=[512, 384], cpos="xy")
display.stop()
_images/demo_unit_disk_23_0.png
[8]:
<pyvirtualdisplay.display.Display at 0x7f6be0527c40>

有限要素法及び積分法の定義

We define the finite element and integration methods to use. We create a MeshFem that defines the degree of freedom of the mesh as rank 1 (scalar).

[9]:
mfu = gf.MeshFem(mesh, 1)

Next we set the finite element used. Classical finite element means a continuous Lagrange element. Setting elements_degree to 2 means that we will use quadratic (isoparametric) elements.

[10]:
elements_degree = 2
mfu.set_classical_fem(elements_degree)

最後に定義するのは、積分法 mim です。GetFEMにはデフォルトの積分方法がないため、積分方法を定義するにはこれが必須です。もちろん、積分法の次数は、選択された有限要素法の便利な積分を行うのに十分なものを選択しなければなりません。ここでは、 elements_degree の2乗で十分です。

[11]:
mim = gf.MeshIm(mesh, pow(elements_degree, 2))

モデル定義

The model object in GetFEM gathers the variables of the models (the unknowns), the data and what are called the model bricks. The model bricks are some parts of the model (linear or nonlinear terms) applied on a single variable or linking several variables. A model brick is an object that is supposed to represent a part of a model. It aims to represent some integral terms in a weak formulation of a PDE model. They are used to make the assembly of the (tangent) linear system (see The model object for more details).

\[[K] \left\{ u \right\} = \left\{ F \right\}\]
[12]:
md = gf.Model("real")
md.add_fem_variable("u", mfu)

Poisson方程式

Poisson方程式を定義するために、 Laplacian項とRHSソース項を定義しなければなりません。 add_Laplacian_brick() を用いて、Laplacian項(GetFEMではブリックと呼ばれています)を加えることができます。

[13]:
md.add_Laplacian_brick(mim, "u")
[13]:
0

GetFEMで定数を定義する場合は、 add_fem_data() メソッドを使用します。

[14]:
F = 1.0
md.add_fem_data("F", mfu)

定数値は set_variable() メソッドで設定できます。ここでは、自由度の大きさのベクトル (ndarray) を渡します。

[15]:
md.set_variable("F", np.repeat(F, mfu.nbdof()))

先に定義した定数 F を用いて add_source_term_brick() メソットでRHSという項を定義しました。

[16]:
md.add_source_term_brick(mim, "u", "F")
[16]:
1

最後に、境界にDirichlet条件を設定します。

[17]:
md.add_Dirichlet_condition_with_multipliers(mim, "u", elements_degree - 1, OUTER_BOUND)
[17]:
2

モデルを解く

モデルが正しく定義されると、次のようにして問題を簡単に解くことができます。

[18]:
md.solve()
[18]:
(0, 1)

解のエクスポート/可視化

有限要素の問題が解かれました。 variable メソッドを使うと、 \(u\) という解を得ることができます。

[19]:
U = md.variable("u")

Sliceオブジェクトのメッシュを使用して、計算された u を出力できます。

[20]:
sl.export_to_vtk("u.vtk", "ascii", mfu, U, "U")

display = Display(visible=0, size=(1280, 1024))
display.start()
p = pv.Plotter()
m = pv.read("u.vtk")
p.add_mesh(m, show_edges=True)
pts = m.points
p.show(window_size=[512, 384], cpos="xy")
display.stop()
_images/demo_unit_disk_47_0.png
[20]:
<pyvirtualdisplay.display.Display at 0x7f6b9815f1c0>

理論解

この問題の理論解は、次の式で与えられます。

\[u(x, y) = \dfrac{1-x^2-y^2}{4}\]
[21]:
evalue = mfu.eval("(1-x*x-y*y)/4")

L 2ノルムとH 1ノルムの誤差は、 compute を使用して計算できます。

[22]:
L2error = gf.compute(mfu, U - evalue, "L2 norm", mim)
H1error = gf.compute(mfu, U - evalue, "H1 norm", mim)
print("Error in L2 norm : ", L2error)
print("Error in H1 norm : ", H1error)
Error in L2 norm :  1.965329030567132e-06
Error in H1 norm :  0.00010936971229957788

このように、エラーのサイズは許容範囲内です。

[ ]: