B. sample1.edp, sample2.edp, sample3.edp

以下に紹介するすべてのプログラムは、単位円版領域において

$$\Laplacian\phi=0$$   in $\Omega$$$,\quad \frac{\rd\phi}{\rd n}=V_n$$   on $\rd\Omega$

を解くプログラムである。境界値 $V_n$ も同じである (すぐ後で説明する)。 コーディングの仕方だけが違っている。

単位円盤領域の境界である単位円周上で、 $-\pi/8\le\theta\le\pi/8$, $\pi/8\le\theta\le\pi/2$, $\pi/2\le\theta\le 3\pi/2$, $3\pi/2\le\theta\le 15\pi/8$ の範囲をそれぞれ $\Gamma_1$, $\Gamma_2$, $\Gamma_3$, $\Gamma_4$ とする。 すなわち

$$\Gamma_1 =\left\{(\cos\theta,\sin\theta)\relmiddle\vert -\pi/8\le\theta\le\pi/8\right\},$$

$$\Gamma_2 =\left\{(\cos\theta,\sin\theta)\relmiddle\vert\pi/8\le\theta\le\pi/2\right\},$$

$$\Gamma_3 =\left\{(\cos\theta,\sin\theta)\relmiddle\vert\pi/2\le\theta\le 3\pi/2\right\},$$

$$\Gamma_4 =\left\{(\cos\theta,\sin\theta)\relmiddle\vert 3\pi/2\le\theta\le 15\pi/8\right\}.$$

$V_n$ は次式で定める:

$$V_n= \left\{ \begin{array}{ll} 4\cos(4\theta) & \text{(on $\Ga... ... & \text{(on $\Gamma_3$)}\\ 0 & \text{(on $\Gamma_4$)}. \end{array} \right.$$

($\Gamma_3$ ではゆっくり入ってくるが、 狭い $\Gamma_1$ では勢いよく出る。$\Gamma_2$, $\Gamma_4$ では出入りがない。)

念のため $\dsp\int_{\rd\Omega}V_n \;\D\sigma=0$ の確認

$$\int_{\rd\Omega}V_n\;\D\sigma =\int_{\Gamma_1}V_n\;\D\sigma +\int_{\Gamma_2}V_n\;\D\sigma +\int... ...mma_4}V_n\;\D\sigma =\int_{\Gamma_1}V_n\;\D\sigma +\int_{\Gamma_3}V_n\;\D\sigma$$

$$=4\int_{-\pi/8}^{\pi/8}\cos(4\theta)\D\theta +\int_{\pi/2}^{3\pi/... ...i/2}^{\pi/2}\cos\theta'\D\theta' +\int_{\pi/2}^{3\pi/2}\cos\theta\;\D\theta =0.$$

デカルト座標で表すと

$$\cos \theta=x,$$

$$\cos(4\theta)=\cos^4\theta-6\cos^2\theta\sin^2\theta+\sin^4\theta =x^4-6x^2y^2+y^4.$$

あるいは

$$\cos(4\theta)=\MyRe\left[e^{4i\theta}\right] =\MyRe\left[(\cos\theta+i\sin\theta)^4\right] =x^4-6x^2y^2+y^4.$$

(繰り返し) sample$\;x$.edp ($x=1,2,3$) はどれもこの問題を解くプログラムである。

もしも、プログラム中で $V_n$ がうまく定義できれば、 sample0.edp と同様に

 solve Laplace(phi,v) =
    int2d(Th)(dx(phi)*dx(v)+dy(phi)*dy(v)) -int1d(Th,Gamma)(Vn*v);

というコードが使える。

しかし $V_n$ を$1$行の関数として定義するのはあまり簡単でない。 一応出来なくはなくて、sample2.edp ではそうしているが、 ここでは sample1.edp のようにすることを勧める。

要点は

$$\int_{\rd\Omega}V_n v\;\D\sigma =\sum_{j=1}^4 \int_{\Gamma_j}V_n v\;\D\sigma =\int_{\Gamma_1}V_n v\;\D\sigma+\int_{\Gamma_3}V_n v\;\D\sigma$$

と分解することである ($\Gamma_2$, $\Gamma_4$ では $V_n=0$ であるから、 $\dsp\int_{\Gamma_j}V_n v\;\D\sigma$ ($j=2,4$) は必要ない)。

$\Gamma_1$, $\Gamma_3$ における $V_n$ を、 それぞれ Vn1, Vn3 と定義すれば

 solve Laplace(phi,v) =
    int2d(Th)(dx(phi)*dx(v)+dy(phi)*dy(v))
   -int1d(Th,Gamma1)(Vn1*v)-int1d(Th,Gamma3)(Vn3*v);

とできる。こうするためには、 もちろんプログラム中で $\Gamma_1$, $\Gamma_2$, $\Gamma_3$, $\Gamma_4$ を用意する必要があるが、それは難しくない (境界を4つに分けた正方形領域のプログラム poisson-kikuchi.edp と同様のことをすれば良い)。 以上の考察から、次のプログラムが得られる。

sample1.edp

// sample1.edp
//   https://m-katsurada.sakura.ne.jp/complex2/sample1.edp
//   2次元非圧縮ポテンシャル流
//     速度ポテンシャル，速度を求め、等ポテンシャル線, 速度場を描く
//   sample1.edp と同じ問題を解く。境界値の指定法が異なる。

border Gamma1(t=-pi/8,pi/8)     { x = cos(t); y = sin(t); }// 勢いよく出る
border Gamma2(t=pi/8,pi/2)      { x = cos(t); y = sin(t); }
border Gamma3(t=pi/2,3*pi/2)    { x = cos(t); y = sin(t); }// ゆっくり入る
border Gamma4(t=3*pi/2,15*pi/8) { x = cos(t); y = sin(t); }
int m=4;
mesh Th=buildmesh(Gamma1(2*m)+Gamma2(3*m)+Gamma3(8*m)+Gamma4(3*m));
plot(Th, wait=1, ps="Th.eps");
// 次の2行は区分1次多項式を使うという意味
fespace Vh(Th,P1);
Vh phi;
// 境界条件の設定
// Gamma1 上で 4 cos(4θ)=4(x^4-6x^2y^2+y^4), Gamma3 上で cos θ=x
func Vn1 = 4*(x^4-6*x^2*y^2+y^4);
func Vn3 = x;
// 整合条件 ∫ V_n ds=0 のチェック
cout << "\int Vn ds=" << int1d(Th,Gamma1)(Vn1)+int1d(Th,Gamma3)(Vn3) << endl;

// 有限要素法で連立1次方程式 A u=f を導く
varf Laplace(phi,v) =
  int2d(Th)(dx(phi)*dx(v)+dy(phi)*dy(v))
  -int1d(Th,Gamma1)(Vn1*v)-int1d(Th,Gamma3)(Vn3*v);

matrix A = Laplace(Vh, Vh);
set(A,solver=CG);
real[int] f = Laplace(0, Vh); f = -f;
// fの定数成分(丸め誤差によるゴミ)を除く (Image A に射影する)
real[int] one(Vh.ndof);
one = 1.0;
one = one / sqrt(one' * one);
f=f-one'*f*one;
cout << "Does f belong to image space? (1, f)="<< (one' * f) << endl;

// 有限要素解を求める
phi[]=A^-1*f;

real meanphi=int2d(Th)(phi)/Th.area; // Th.area は int2d(Th)(1.0)
cout << "mean phi=" << meanphi << endl;
phi=phi-meanphi;
plot(phi,ps="contourpotential.eps",wait=1);

// ベクトル場 (v1,v2)=∇φ を描く (ちょっと雑なやり方)
Vh 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);

図 3: sample1.edp の結果 ( $\theta \in [\pi /2,3\pi /2]$ でゆっくり入り、 $\theta \in [-\pi /8,\pi /8]$で勢い良く出る)

以下、参考までに、$V_n$ を1命令(長いので2行になっている)で書いたプログラム sample2.edp を紹介しておく。 いわゆる偏角の主値 (値が $(-\pi,\pi]$ に属する) を返す関数 atan2(y,x) を使うが、 やや不自然な感じがするのは否めない (角度の範囲が $[0,2\pi)$ でないことにも注意が必要である)。

強引に1行関数を作る

func Vn=(atan2(y,x)>=-pi/8 && atan2(y,x)<=pi/8) * 4.0 * (x^4-6*x*x*y*y+y^4)
       + (x<=0) * x;

sample2.edp

// sample2.edp
//   https://m-katsurada.sakura.ne.jp/complex2/sample2.edp
//   2次元非圧縮ポテンシャル流
//     速度ポテンシャル，速度を求め、等ポテンシャル線, 速度場を描く
//   sample1.edp と同じ問題を解く。境界値の指定法が異なる。

border Gamma(t=0,2*pi) { x = cos(t); y = sin(t); } // 円盤領域
int m=64;
mesh Th=buildmesh(Gamma(m));
plot(Th, wait=1, ps="Th.eps");
// 次の2行は区分1次多項式を使うという意味
fespace Vh(Th,P1);
Vh phi;
// 境界条件の設定 (atan2(y,x) は [-pi,pi) の範囲の値を返す)
// Gamma1 上で 4 cos(4θ)=4(x^4-6x^2y^2+y^4), Gamma3 上で cos θ=x
func Vn=(atan2(y,x)>=-pi/8 && atan2(y,x)<=pi/8) * 4.0 * (x^4-6*x*x*y*y+y^4)
       + (x<=0) * x;
// 整合条件 ∫ V_n ds=0 のチェック
cout << "\int Vn ds=" << int1d(Th,Gamma)(Vn) << endl;

// 有限要素法で連立1次方程式 A u=f を導く
varf Laplace(phi,v) =
  int2d(Th)(dx(phi)*dx(v)+dy(phi)*dy(v)) -int1d(Th,Gamma)(Vn*v);

matrix A = Laplace(Vh, Vh);
set(A,solver=CG);
real[int] f = Laplace(0, Vh); f = -f;
// fの定数成分(丸め誤差によるゴミ)を除く (Image A に射影する)
real[int] one(Vh.ndof);
one = 1.0;
one = one / sqrt(one' * one);
f=f-one'*f*one;
cout << "Does f belong to image space? (1, f)="<< (one' * f) << endl;

// 有限要素解を求める
phi[]=A^-1*f;

real meanphi=int2d(Th)(phi)/Th.area; // Th.area は int2d(Th)(1.0)
cout << "mean phi=" << meanphi << endl;
phi=phi-meanphi;
plot(phi,ps="contourpotential.eps",wait=1);

// ベクトル場 (v1,v2)=∇φ を描く (ちょっと雑なやり方)
Vh 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);

図 4: sample2.edp の結果 ( $\theta \in [\pi /2,3\pi /2]$でゆっくり入り、 $\theta \in [-\pi /8,\pi /8]$ で勢い良く出る)

普通のプログラミング言語であれば、 $V_n$ 程度の場合分けを含む関数は、if 文を使って書くのだろう。 そのように関数 Vn(x,y) を作ってみよう。 しかしこの Vn() を弱形式にそのまま使うことは出来ないので (FreeFEM の文法の仕様上の制限)、 fespace の変数 Vn2 を用意して Vn2=Vn(x,y) と代入して、 Vn2 を弱形式中で用いる (領域内部での値を計算しているが、それはめちゃくちゃな値で、 後の計算で利用もしていない。プログラムは読みやすい)。

sample3.edp

// sample3.edp
//   https://m-katsurada.sakura.ne.jp/complex2/sample3.edp
//   2次元非圧縮ポテンシャル流
//     速度ポテンシャル，速度を求め、等ポテンシャル線, 速度場を描く
//   sample1.edp と同じ問題を解く。境界値の指定法が異なる。

border Gamma(t=0,2*pi) { x = cos(t); y = sin(t); } // 円盤領域
int m=64;
mesh Th=buildmesh(Gamma(m));
plot(Th, wait=1, ps="Th.eps");
// 次の2行は区分1次多項式を使うという意味
fespace Vh(Th,P1);
Vh phi;

// 境界条件の設定
func real Vn(real x, real y)
{
  real theta=atan2(y,x);
  if (theta >= - pi / 8 && theta <= pi / 8)
    return 4.0 * (x^4 - 6*x^2*y^2 + y^4);
  else if (x <= 0)
    return x;
  else
    return 0.0;
}

// 境界値 (弱形式の中で func Vn() は直接指定できないので値を求めておく)
Vh Vn2;
Vn2 = Vn(x,y);
cout << "int Vn ds=" << int1d(Th,Gamma)(Vn2) << endl;
// 有限要素法で連立1次方程式 A u=f を導く
varf Laplace(phi,v) =
  int2d(Th)(dx(phi)*dx(v)+dy(phi)*dy(v)) -int1d(Th,Gamma)(Vn2*v);

matrix A = Laplace(Vh, Vh);
set(A,solver=CG);
real[int] f = Laplace(0, Vh); f = -f;
// fの定数成分(丸め誤差によるゴミ)を除く (Image A に射影する)
real[int] one(Vh.ndof);
one = 1.0;
one = one / sqrt(one' * one);
f=f-one'*f*one;
cout << "Does f belong to image space? (1, f)="<< (one' * f) << endl;

// 有限要素解を求める
phi[]=A^-1*f;

real meanphi=int2d(Th)(phi)/Th.area; // Th.area は int2d(Th)(1.0)
cout << "mean phi=" << meanphi << endl;
phi=phi-meanphi;
plot(phi,ps="contourpotential.eps",wait=1);

// ベクトル場 (v1,v2)=∇φ を描く (ちょっと雑なやり方)
Vh 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);