3 おまけ: gnuplot で可視化

従来、FreeFem++ では、等高線描画やベクトル場の表示などは出来るが、 グラフの鳥瞰図描画などはサポートしていなかった (今は plot(...,dim=3,...); とするだけで出来る)。 そこで gnuplot を使ってグラフを描く方法を紹介する。

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;

``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 で可視化するには、 例えば次のようにする (図 3)。

poisson-disk.g

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

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