4.3 例: 有限要素法で解いた Poisson 方程式の境界値問題の解の可視化

(工事中: 説明が舌足らずですが、とりあえず実例を。なお、 付録 B.2 に差分法の解の可視化例の紹介をしている。 そちらは格子点 $(x_i,y_j)$ での値 $u(x_i,y_j)$ からグラフを描く話で、 本来はそれを詳しく解説すべきものだったかも…)

単位円盤領域 $\Omega:=\{(x,y)\in\R^2; x^2+y^2<1\}$ における Poisson 方程式の境界値問題

$$-\Laplacian u(x,y)=f(x,y)$$   $$\mbox{in $\Omega$}$$$$,\quad u(x,y)=0$$   $$\mbox{on $\Gamma:=\rd\Omega$}$$

は、 次のような FreeFEM++ プログラムで解ける (ここでは $f(x,y)=x y$ としている)。

FreeFEM++ プログラム poisson-disk.edp

// 境界の定義 (単位円), いわゆる正の向き
border Gamma(t=0,2*pi) { x=cos(t); y=sin(t); }
// 三角形要素分割を生成 (境界を50に分割) 
mesh Th = buildmesh(Gamma(50));
plot(Th,ps="mesh.eps");
plot(Th,wait=1);
// 有限要素空間は P1 (区分的１次多項式) 要素
fespace Vh(Th,P1);
Vh u,v;
// Poisson 方程式 -△u=f の右辺
func f = x*y;
// 現在時刻をメモ
real start = clock();
// 問題を解く
solve Poisson(u,v)
 = int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v))-int2d(Th)(f*v)
   +on(Gamma,u=0);
// 可視化
plot(u,ps="poisson-disk.eps");
plot(u);
//
{
  ofstream ug("u-g.txt");
  for (int i=0; i<Th.nt; i++) {
    for (int j=0; j<3; j++) {
      ug << Th[i][j].x << " " << Th[i][j].y << " " << u[][Vh(i,j)]<<endl;
    }
    ug << Th[i][0].x << " " << Th[i][0].y << " " << u[][Vh(i,0)]<<"\n\n\n";
   }
}
// 計算時間を表示
cout << " CPU time= " << clock() - start << endl;

図 9: 円盤の三角形分割とPoisson 方程式の解の等高線表示 (FreeFEM++ による)

FreeFEM++ では、等高線描画やベクトル場の表示などは出来るが、 グラフの鳥瞰図描画などはサポートしていないようである (最近は出来るようになっているみたいなので、 ここに書いてあることは「古い」かも知れないけれど、 参考のため、しばらく残します)。 そこで gnuplot を使ってグラフを描いてみる。

“u-g.txt” は次のようなデータである (3次元空間における三角形が並んでいる)。

0.735953 0.529892 0.00584987
0.876307 0.481754 6.33978e-033
0.809017 0.587785 7.53444e-033
0.735953 0.529892 0.00584987


0.512882 0.46938 0.0104192
0.577982 0.366799 0.00923669
0.656968 0.448345 0.00912418
0.512882 0.46938 0.0104192

(中略)
......
0.656968 0.448345 0.00912418
0.447782 0.737755 0.00678195
0.425779 0.904827 4.83201e-033
0.317582 0.777124 0.00599877
0.447782 0.737755 0.00678195


0.187381 0.982287 3.00659e-033
0.0627905 0.998027 8.0753e-034
0.149059 0.864343 0.00242012
0.187381 0.982287 3.00659e-033

このデータを gnuplot で可視化するには、 例えば次のようにする (図 10)。

poisson-disk.g

set hidden3d
set palette rgbformulae 33,13,10
splot "u-g.txt" with lines palette

図 10: Poisson 方程式の解のグラフ表示 (gnuplot による)