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

mesh Th=square(10,10); //領域$]0,1[^2$の三角形分割
fespace Vh(Th,P1); //１次要素空間
Vh u, v; // $u,v\in V_h$
func f=2*((x-x^2)+(y-y^2)); // $f(x,y)=2((x-x^2)+(y-y^2))$
varf a(u,v) = int2d(Th)( dx(u)*dx(v) + dy(u)*dy(v))
              + on(1,2,3,4,u=0) ;
/* 剛性行列
$A_{ij}=\int_{\Omega}\nabla \varphi_{j}\cdot\nabla\varphi_{i}dxdy
\qquad \textrm{with }\varphi_{i}=0\quad\textrm{on }\partial\Omega$
*/
matrix A=a(Vh,Vh); // 剛性行列
varf l(u,v) = int2d(Th)(f*v)+on(1,2,3,4,u=0);
Vh F; F[] = l(0,Vh);
u = A^-1*F; // $u=A^{-1}F$
 plot(u,wait=1);