// potential2d-v0.edp
//   http://nalab.mind.meiji.ac.jp/~mk/complex2/potential2d-v0.edp
//   2次元非圧縮ポテンシャル流
//     速度ポテンシャル，速度を求め、等ポテンシャル線, 速度場を描く

border Gamma(t=0,2*pi) { x = cos(t); y = sin(t); } // 円盤領域
int m=40;
mesh Th=buildmesh(Gamma(m));
plot(Th, wait=1, ps="Th.eps");
// 次の2行は区分1次多項式を使うという意味
fespace Vh(Th,P1);
Vh phi, v, v1, v2;
// 境界条件の設定
func Vn=x+2*y;// Ωが単位円で, V=(1,2) のとき V・n=x+2y
func Vn2=((x>0&&y>0) || (x<0&&y<0))*(x+2*y); // 右上と左下のみ

// 速度ポテンシャルφを求め、その等高線 (等ポテンシャル線) を描く
solve Laplace(phi,v) =
  int2d(Th)(dx(phi)*dx(v)+dy(phi)*dy(v)) -int1d(Th,Gamma)(Vn2*v);
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);