サンプルプログラム解説

potential2d-v0.edp

   1 // potential2d-v0.edp
   2 //   http://nalab.mind.meiji.ac.jp/~mk/complex2/potential2d-v0.edp
   3 //   2次元非圧縮ポテンシャル流
   4 //     速度ポテンシャル，速度を求め、等ポテンシャル線, 速度場を描く
   5 
   6 border Gamma(t=0,2*pi) { x = cos(t); y = sin(t); } // 円盤領域
   7 int m=40;
   8 mesh Th=buildmesh(Gamma(m));
   9 plot(Th, wait=1, ps="Th.eps");
  10 // 次の2行は区分1次多項式を使うという意味
  11 fespace Vh(Th,P1);
  12 Vh phi, v, v1, v2;
  13 // 境界条件の設定
  14 func Vn=x+2*y; // Ωが単位円で, V=(1,2) のとき V・n=x+2y
  15 
  16 // 速度ポテンシャルφを求め、その等高線 (等ポテンシャル線) を描く
  17 solve Laplace(phi,v) =
  18   int2d(Th)(dx(phi)*dx(v)+dy(phi)*dy(v)) -int1d(Th,Gamma)(Vn*v);
  19 plot(phi,ps="contourpotential.eps",wait=1);
  20 
  21 // ベクトル場 (v1,v2)=∇φ を描く (ちょっと雑なやり方)
  22 v1=dx(phi); v2=dy(phi);
  23 plot([v1,v2],ps="vectorfield.eps",wait=1);
  24 
  25 // 等ポテンシャル線とベクトル場を同時に描く
  26 plot([v1,v2],phi,ps="both.eps", wait=1);

• 6行目で領域$\Omega$の境界$\Gamma$を指定している (それで $\Omega$ が決まる)。
• 8行目で、$\Gamma$を$m$分割して、$\Gamma$の囲む範囲を三角形分割して、 それを mesh 型の変数 Th に代入している。
• 11行目、有限要素空間 $V_h$ を区分的1次関数の空間と定義している。 この辺についてはこの講義では説明を省略する (知りたい人は有限要素法のテキストを読んで下さい)。
• 12行目、$\phi$ と試験関数 $v$, 流速ベクトル場 $\bm{v}$ の成分 $v_1,v_2$ を、 Vh の要素とする。
• 14行目、境界条件の設定をしている。 ここで $\bm{v}\cdot\bm{n}=x+2y$ とする理由は上の注意(1)で説明した。
• 17〜18行目、弱形式の定義。ここを修正する必要が生じる可能性は低い。