// potential2d.edp (version 0)
//   2次元非圧縮ポテンシャル流
//     速度ポテンシャル、速度を求め
//     等ポテンシャル線, 速度場を描く

border Gamma(t=0,2*pi) { x = cos(t); y = sin(t); } // 円盤領域
int m=40;
mesh Th=buildmesh(Gamma(m));
fespace Vh(Th,P1);
Vh phi, v, u;
func Vn=x+2*y; // Ωが単位円で, V=(1,2) のとき V・n=x+2y

// 速度ポテンシャルφを求め、その等高線 (等ポテンシャル線) を描く
solve Laplace(phi,v) =
  int2d(Th)(dx(phi)*dx(v)+dy(phi)*dy(v))
  -int1d(Th,Gamma)(Vn*v);
plot(phi,ps="contourpotential.eps",wait=1);

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