Post date : 2024-09-04. ( update : 2024-10-05 )
FreeFEM では弱形式を problem で記述することで自動的に連立一次方程式の行列と右辺のベクトルを生成し計算してくれるが，varf + matrix で明示的に表記したほうが計算も速く，並列ソルバーや領域分割 (hpddm など)も使用できる．
方程式の左辺と右辺の切り分け，境界条件の設定を確認するためにテストプログラムを作成し，確認した．
problem laplaceP(u,v) = int2d(Th)(Grad(u)'*Grad(v)) - int2d(Th)(f*v) + on(1,2,3,4,u=ubc0);
->
varf laplace(u,v) = int2d(Th)(Grad(u)'*Grad(v)) + on(1,2,3,4,u=ubc0);
varf RHS(u,v)=int2d(Th)(f*v) + on(1,2,3,4,u=ubc0);
matrix A = laplace(Vh,Vh); set(A,solver=CG);
real[int] b=RHS(0,Vh);
u[] = A^-1*b;
デフォルト設定，tgv=1e30, tgv=-1, problem の4通りで計算した．
int n=5;
mesh Th=square(n,n);
fespace Vh(Th,P1);
Vh u, v;
plot(Th,wait=1,ps="Th.eps");
// \Delta u = -f
real ubc0=10;
macro Grad(u) [dx(u),dy(u)]// EOM
macro f() (1) //
varf laplace(u,v) = int2d(Th)(Grad(u)'*Grad(v)) + on(1,2,3,4,u=ubc0);
varf RHS(u,v)=int2d(Th)(f*v) + on(1,2,3,4,u=ubc0);
// --- default setting ---
matrix A = laplace(Vh,Vh); set(A,solver=CG);
real[int] b=RHS(0,Vh);
u[] = A^-1*b;
{ofstream f1("A.dat"); f1 << A;}
{ofstream f2("b.dat"); f2 << b;}
{ofstream f3("u.dat"); f3 << u[];}
plot(u,wait=1,value=1,ps="matrix-print.eps");
// --- tgv=1e30 ---
real tgv1=1e30;
matrix A1 = laplace(Vh,Vh,tgv=tgv1); set(A1,solver=CG);
real[int] b1=RHS(0,Vh,tgv=tgv1);
u[] = A1^-1*b1;
{ofstream f1("A-tgv"+tgv1+".dat"); f1 << A1;}
{ofstream f2("b-tgv"+tgv1+".dat"); f2 << b1;}
{ofstream f3("u-tgv"+tgv1+".dat"); f3 << u[];}
plot(u,wait=1,value=1,ps="matrix-print-tgv"+tgv1+".eps");
// --- tgv=-1 ---
real tgv2=-1;
matrix A2 = laplace(Vh,Vh,tgv=tgv2); set(A2,solver=CG);
real[int] b2=RHS(0,Vh,tgv=tgv2);
u[] = A2^-1*b2;
{ofstream f1("A-tgv"+tgv2+".dat"); f1 << A2;}
{ofstream f2("b-tgv"+tgv2+".dat"); f2 << b2;}
{ofstream f3("u-tgv"+tgv2+".dat"); f3 << u[];}
plot(u,wait=1,value=1,ps="matrix-print-tgv"+tgv2+".eps");
// --- problem ---
solve laplaceP(u,v) = int2d(Th)(Grad(u)'*Grad(v)) - int2d(Th)(f*v) + on(1,2,3,4,u=ubc0);
{ofstream f3("u-prob.dat"); f3 << u[];}
plot(u,wait=1,value=1,ps="matrix-print-prob.eps");
-- Square mesh : nb vertices =36 , nb triangles = 50 , nb boundary edges 20
GC: converge in 3 g=2.82422e-28 rho= 1.32993 gamma= 2.1796e-30
GC: converge in 9 g=2.04309e-63 rho= 1.2398 gamma= 9.78178e-29
GC: converge in 3 g=3.77144e-29 rho= 0.833333 gamma= 3.68304e-26
-- Solve :
n=5 -> dx=0.2
A
tgv=1e30:
0 1 2 3 4 5 6 7 8
0:1e+30, -0.5, , , , , -0.5, 0, ,
1: -0.5, 1e+30, -0.5, , , , , -1, 0,
2: , -0.5, 1e30, -0.5, , , , , -1,
3: , , , 1e30, -0.5, , ,
4: , , , -0.5, 1e30, -0.5, ,
5: , , , , , 1e30, ,
6: -0.5, , , , , , , -1,
7: 0, -1, , , , , -1, 4, -1,
8: , 0, -1, , , , , -1, 4,
tgv=-1:
0 1 2 3 4 5 6 7 8
0: , , , , , , , ,
1: , 1, 0, , , , , 0, 0,
2: , 0, 1, 0, , , , , 0,
3: , , , 1, 0, , ,
4: , , , 0, 1, 0, ,
5: , , , , , 1, 0,
6: , , , , , 0, 1, ,
7: , 0, , , , , , 1, 0,
8: , 0, -1, , , , , -1, 4,
b
tgv=1e30
1e+31 1e+31 1e+31 1e+31 1e+31
1e+31 1e+31 0.04 0.04 0.04
0.04 1e+31 1e+31 0.04 0.04
0.04 0.04 1e+31 1e+31 0.04
0.04 0.04 0.04 1e+31 1e+31
0.04 0.04 0.04 0.04 1e+31
1e+31 1e+31 1e+31 1e+31 1e+31
1e+31
tgv=-1
10 10 10 10 10
10 10 0.04 0.04 0.04
0.04 10 10 0.04 0.04
0.04 0.04 10 10 0.04
0.04 0.04 0.04 10 10
0.04 0.04 0.04 0.04 10
10 10 10 10 10
10
u
10 10 10.03333333 10.04666667 10.04666667
10.03333333 10 10 10.04666667 10.06666667
10.06666667 10.04666667 10 10 10.04666667
10.06666667 10.06666667 10.04666667 10 10
10.03333333 10.04666667 10.04666667 10.03333333 10
10 10 10 10 10
10
わかること
tgv
デフォルトは 1e30 (tgv=1e30 と同じ)
tgv=1e30 の場合
境界上の節点では A の対角成分が 1e30， 境界での b は境界値 x 1e30
tgv=-1 の場合
境界上の節点ではAの行成分は対角成分 (=1) のみ，境界での b は境界値
コーナー(例えば (0,j), (i,0))では A の成分がない
[NG] varf RHS(u,v)=int2d(Th)(f*v);
左辺の行列 A は正しいが，右辺 b に境界条件がセットされない
[NG] varf laplace0(u,v)=int2d(Th)(Grad(u)'*Grad(v))-int2d(Th)(f*v)+on(1,2,3,4,u=ubc0);でreal[int] b0=laplace0(0,Vh);
右辺の符号が逆転する (\Delta u = f となってしまう)
real[int] b01=laplace0(0,Vh); b01=-b01; とすると，右辺 b の境界値の符号だけ反転してしまう
さらに，varf BC01(u,v) = on(1,2,3,4,u=ubc0); real[int] rhbc = BC01(0,Vh); b01 = rhbc ? rhbc : b01; を追加すればOK
tgv をデフォルト値から変更する場合，A だけでなく b でも記述が必要
n<5 ではメッシュが作れない
ubc0=0 としても右辺 b の境界部分はきちんと 0 になる
注意：
matrix solver は Au=b を想定しているが、IPOPT は Au+b=0 を想定している