/* メッシュを細かくすると厳密解に近づくことを示す。\\
与えられた関数 $f(x,y)=1$に対して次を考える。$$
-\Delta u = f\quad \textrm{in }\Omega; \qquad u=0 \quad \textrm{on }\partial\Omega$$
を考える。
*/

mesh Th=square(100,100); //メッシュ
fespace Vh(Th,P1); //有限要素空間
Vh u,v; // $u,v\in V_{h}(\Omega)$

func f=x*y;
varf l(u,v) = int2d(Th)(f*v)+on(1,2,3,4,u=0);
Vh F; F[] = l(0,Vh); // $F=\int_{\Omega}f\varphi_i$
real cpu=clock();
varf a1(u,v,solver=LU) = int2d(Th)( dx(u)*dx(v) + dy(u)*dy(v))
              + on(1,2,3,4,u=0) ;
matrix A1=a1(Vh,Vh); // LU分解
u[]=A1^-1*F[]; // $u=A^{-1}F$
cout << "LU: CPU time = " << clock()-cpu << endl;
plot(u);
cpu=clock();
varf a2(u,v,solver=CG) = int2d(Th)( dx(u)*dx(v) + dy(u)*dy(v))
              + on(1,2,3,4,u=0) ;
matrix A2=a2(Vh,Vh); // Conjugate Gradient
u[]=A2^-1*F[]; // $u=A^{-1}F$
cout << "CG: CPU time = " << clock()-cpu << endl;
plot(u);
cpu=clock();
varf a3(u,v,solver=Crout) = int2d(Th)( dx(u)*dx(v) + dy(u)*dy(v))
              + on(1,2,3,4,u=0) ;
matrix A3=a3(Vh,Vh); //Crout分解
u[]=A3^-1*F[]; // $u=A^{-1}F$
cout << "Crout: CPU time = " << clock()-cpu << endl;
plot(u);
cpu=clock();
varf a4(u,v,solver=Cholesky) = int2d(Th)( dx(u)*dx(v) + dy(u)*dy(v))
              + on(1,2,3,4,u=0) ;
matrix A4=a4(Vh,Vh); //Cholesky分解
u[]=A4^-1*F[]; // $u=A^{-1}F$
cout << "Cholesky: CPU time = " << clock()-cpu << endl;
plot(u);
cpu=clock();
varf a5(u,v,solver=GMRES) = int2d(Th)( dx(u)*dx(v) + dy(u)*dy(v))
              + on(1,2,3,4,u=0) ;
matrix A5=a5(Vh,Vh); //GMRES
u[]=A5^-1*F[]; // $u=A^{-1}F$
cout << "GMRES: CPU time = " << clock()-cpu << endl;
plot(u);
cpu=clock();
varf a6(u,v,solver=UMFPACK) = int2d(Th)( dx(u)*dx(v) + dy(u)*dy(v))
              + on(1,2,3,4,u=0) ;
matrix A6=a6(Vh,Vh); //UMPACK
u[]=A6^-1*F[]; // $u=A^{-1}F$
cout << "UMPACK: CPU time = " << clock()-cpu << endl;
plot(u);