3 おまけ: gnuplot で可視化

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

poisson-g.edp

// poisson-g.edp

// 境界の定義 (単位円), いわゆる正の向き
border Gamma(t=0,2*pi) { x=cos(t); y=sin(t); }
// 三角形要素分割を生成 (境界を50に分割) 
mesh Th = buildmesh(Gamma(50));
plot(Th,ps="Th.eps",wait=true);
// 有限要素空間は 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");
// 計算結果をファイルに出力
{
  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;

入手と実行

curl -O https://m-katsurada.sakura.ne.jp/program/fem/poisson-g.edp
FreeFem++ poisson-g.edp

実行が終わると “u-g.txt” というファイルが作成される。 “u-g.txt” は次のようなデータである (3次元空間における三角形が並んでいる)。

u-g.txt

0.84068 0.327399 0.00414118
0.794983 0.411712 0.00532531
0.751583 0.286674 0.00623707
0.84068 0.327399 0.00414118


0.870933 0.153872 0.00249565
0.84068 0.327399 0.00414118
0.751583 0.286674 0.00623707
0.870933 0.153872 0.00249565

(中略)
......
0.1092 0.735914 0.00288613
0.140979 0.864349 0.00225845
0.0271079 0.868281 0.000433488
0.1092 0.735914 0.00288613


-0.0867634 0.872212 -0.00146171
-0.203275 0.85648 -0.00323151
-0.110736 0.746398 -0.00290419
-0.0867634 0.872212 -0.00146171

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

poisson-g.gp

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

この内容を “gnuplot> ” という gnuplot のプロンプロトに入力すれば良いが、 自動的にやるには、以下の2つのコマンドをターミナルで実行しても良い。

curl -O https://m-katsurada.sakura.ne.jp/program/fem/poisson-g.gp
gnuplot poisson-g.gp

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