メッシュを読み込む

2021年6月20日

はじめに

FEniCS でメッシュを読み込む。メッシュの作成には gmsh を用いる。

環境

• FEniCS 2019.2.0
• gmsh 4.4.1
• Ubuntu 20.04 LTS

問題

1 m x 1 m の領域を考え、左端を 300 K、右端を 500 K とした場合の温度分布を求める。

メッシュを用意する

gmsh でメッシュを作成する。まず、形状を用意する。

square.geo

Point(1) = {0.5, 0, 0, 1.0};
Recursive Delete {
  Point{1}; 
}
Point(1) = {0, 0, 0, 1.0};
Point(2) = {1, 0, 0, 1.0};
Point(3) = {1, 1, 0, 1.0};
Point(4) = {0, 1, 0, 1.0};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Curve Loop(1) = {4, 1, 2, 3};
Physical Curve("left") = {4};
Physical Curve("right") = {2};
Physical Curve("top") = {3};
Physical Curve("bottom") = {1};
Plane Surface(1) = {1};
Physical Surface("fluid") = {1};

コマンドでメッシュを作成する。この際、フォーマットを指定する。

$ gmsh -2 -clmin 0.05 -clmax 0.05 square.geo -format msh2

square.msh ができるので、これを FEniCS 用に変換する。

$ dolfin-convert square.msh square.xml

これにより、square.xml、square_facet_region.xml、square_physical_region.xml が作成される。

コード

mesh.py

from fenics import *

mesh = Mesh("square.xml")
markers = MeshFunction("size_t", mesh, 'square_physical_region.xml')
boundaries = MeshFunction('size_t', mesh, 'square_facet_region.xml')
V = FunctionSpace(mesh, "P", 1)

bc_left = DirichletBC(V, Constant(300.), boundaries, 1)
bc_right = DirichletBC(V, Constant(500.), boundaries, 2)
bc = [bc_left, bc_right]

u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0.)

a = dot(grad(u), grad(v))*dx
L = f*v*dx

T = Function(V, name="T")
solve(a == L, T, bc)

vtkfile = File("result.pvd")
vtkfile << T

実行

$ python3 mesh.py

コードの解説

mesh = Mesh("square.xml")
markers = MeshFunction("size_t", mesh, 'square_physical_region.xml')
boundaries = MeshFunction('size_t', mesh, 'square_facet_region.xml')
V = FunctionSpace(mesh, "P", 1)

メッシュを読み込んでいる。boundaries が境界になる。

境界条件の設定。ここで、DirichletBC の最後の引数は境界の番号である。番号は gmsh のメッシュファイル (この場合は square.msh) で確認できる。

square.msh

$MeshFormat
2.2 0 8
$EndMeshFormat
$PhysicalNames
5
1 1 "left"
1 2 "right"
1 3 "top"
1 4 "bottom"
2 5 "fluid"
...

この場合、1 は "left"、2 は "right" である。