3.3.1 サンプル・プログラム

potential2d-v2.edp というサンプル・プログラムを用意した。これは単位円盤領域 $\Omega$ の境界 $\rd\Omega$ を

	$\displaystyle \Gamma_1:(x,y)=(\cos\theta,\sin\theta)$ ( $\theta\in[0,\pi/6]$ ) $\displaystyle ,$
	$\displaystyle \Gamma_2:(x,y)=(\cos\theta,\sin\theta)$ ( $\theta\in[\pi/6,\pi]$ ) $\displaystyle ,$
	$\displaystyle \Gamma_3:(x,y)=(\cos\theta,\sin\theta)$ ( $\theta\in[\pi,4\pi/3]$ ) $\displaystyle ,$
	$\displaystyle \Gamma_3:(x,y)=(\cos\theta,\sin\theta)$ ( $\theta\in[4\pi/3,2\pi]$ )

$\displaystyle v_n= \left\{ \begin{array}{ll} -2 & \text{(on $\Gamma_1$)}\\ ... ...n $\Gamma_3$)}\\ 0 & \text{(on $\Gamma_2\cup\Gamma_4$)} \end{array} \right.$

[

l]potential2d-v2.edp

// potential2d-v2.edp
//   https://m-katsurada.sakura.ne.jp/ana/potential2d-v2.edp
//   2次元非圧縮ポテンシャル流
//     速度ポテンシャル，速度を求め、等ポテンシャル線, 速度場を描く

border Gamma1(t=0,pi/6)      { x = cos(t); y = sin(t); }
border Gamma2(t=pi/6,pi)     { x = cos(t); y = sin(t); }
border Gamma3(t=pi,4*pi/3)   { x = cos(t); y = sin(t); }
border Gamma4(t=4*pi/3,2*pi) { x = cos(t); y = sin(t); }
int m=20;
mesh Th=buildmesh(Gamma1(m)+Gamma2(m)+Gamma3(m)+Gamma4(m));
plot(Th, wait=1, ps="Th.eps");
// 次の2行は区分1次多項式を使うという意味
fespace Vh(Th,P1);
Vh phi, psi, v, v1, v2;
// 境界条件の設定
func Vn1=-2.0;                  // Gamma1 では勢いよく注水
func Vn3=1.0;                   // Gamma3 からゆるく排水
func Vn=(x>0)*(-2.0)+(x<0)*1.0; // xの符号でGamma1, Gamma3 を区別してまとめる

// 速度ポテンシャルφを求める
solve Laplace(phi,v) =
  int2d(Th)(dx(phi)*dx(v)+dy(phi)*dy(v))
    -int1d(Th,Gamma1)(Vn1*v)-int1d(Th,Gamma3)(Vn3*v);
//    -int1d(Th,Gamma1,Gamma3)(Vn*v);
// 平均を0にする (解の一意性がないため、時々生じる大きなズレを消す
real mean=int2d(Th)(phi)/int2d(Th)(1);
phi=phi-mean;
// φの等高線 (等ポテンシャル線) を描く
plot(phi,ps="contourpotential.eps",wait=1);

// ベクトル場 (v1,v2)=∇φ を描く (ちょっと雑なやり方)
v1=dx(phi); v2=dy(phi);
plot([v1,v2],ps="vectorfield.eps",wait=1);

// 等ポテンシャル線とベクトル場を同時に描く
plot([v1,v2],phi,ps="both.eps", wait=1);

// 流れ関数を求める
func angle=(atan2(y,x)>=0)*atan2(y,x)+(atan2(y,x)<0)*(atan2(y,x)+2*pi);
solve LaplaceD(psi,v) =
  int2d(Th)(dx(psi)*dx(v)+dy(psi)*dy(v))
    +on(Gamma1,psi=-2*angle(x,y))
    +on(Gamma2,psi=-pi/3)
    +on(Gamma3,psi=angle(x,y)-4*pi/3)
    +on(Gamma4,psi=0);
plot(psi,ps="streamline.eps",wait=1);
plot([v1,v2],phi,psi,ps="all.eps", wait=1);