2. プログラム作成から出力まで

サンプル問題の実行

とりあえずサンプルプログラムを実行してみよう

今回は以下の2次元ポアソン方程式を考える

1. ファイル作成
Command Shell上で以下を入力すればファイルを作成できる。

python
f = open("poisson.py", "w")
f.close()
[Ctrl]+[Z]を入力

2. プログラム作成
ホームディレクトリ(Windowsの場合は"c:/FEniCS")でpoisson.pyを開く（テキストエディタで可）。
以下を入力

from dolfin import * 
# mesh
mesh = UnitSquare(30, 30)
V = FunctionSpace(mesh, 'Lagrange', 1)
# boundary condition
u0 = Expression('0')
def u0_boundary(x, on_boundary):
  return on_boundary
bc = DirichletBC(V, u0, u0_boundary)
# Poisson equation
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(1.0)
a = inner(nabla_grad(u), nabla_grad(v))*dx
L = f*v*dx
# solve the equation
u = Function(V)
solve(a == L, u, bc,solver_parameters={'linear_solver': 'cg','preconditioner': 'ilu'})
# plot solution and mesh
plot(u)
plot(mesh)
interactive()
# Dump solution to file in VTK format
file = File("poisson.pvd")
file << u

3. プログラム実行
Command Shell上で

 python poisson.py

と入力すればプログラムが実行できる

プログラム解説

from dolfin import * 

FEniCSプログラムの最初に書く必須項目。

# mesh
mesh = UnitSquare(30, 30)
V = FunctionSpace(mesh, 'Lagrange', 1)

メッシュ作成
mesh = UnitSquare(30, 30)
…(0,1)×(0,1)の正方領域に縦横30ずつの格子点を配置し、その格子点を節点とした三角形メッシュを作成し、そのデータをmeshに格納
V = FunctionSpace(mesh, 'Lagrange', 1)
…関数空間Vをmesh上Lagrange型補間関数(P1要素)による有限次元関数空間で定義

# boundary condition
u0 = Expression('0')
def u0_boundary(x, on_boundary):
  return on_boundary
bc = DirichletBC(V, u0, u0_boundary)

境界条件
u0 = Expression('0')
…境界条件をu=0に設定
def u0_boundary(x, on_boundary):
return on_boundary
…境界情報変数u0_boundaryを定義
bc = DirichletBC(V, u0, u0_boundary)
…境界情報（関数空間、u0）をu0_boundaryに格納

# Poisson equation
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(1.0)
a = inner(nabla_grad(u), nabla_grad(v))*dx
L = f*v*dx

関数の設定
u = TrialFunction(V)
…関数空間V上の未知変数uの設定
v = TestFunction(V)
…関数空間V上の試験関数vの設定
f = Constant(1.0)
…関数fを定数関数f=1で設定

# solve the equation
u = Function(V)
solve(a == L, u, bc,solver_parameters={'linear_solver': 'cg','preconditioner': 'ilu'})
# plot solution and mesh
plot(u)
plot(mesh)
interactive()
# Dump solution to file in VTK format
file = File("poisson.pvd")
file << u