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

mesh Th=square(10,10); //メッシュ
fespace Vh(Th,P1); //有限要素空間
Vh u,v; // $u,v\in V_{h}(\Omega)$
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); // 剛性行列
func f=1;
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$
u[]=A^-1*F[]; // $u=A^{-1}F$
plot(u);