/* メッシュを細かくすると厳密解に近づくことを示す。\\
与えられた関数 $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 Th0=square(100,100); //領域$\Omega=]0,1[^2$の三角形分割
fespace Vh0(Th0,P1); // 比較の基準となる要素空間
Vh0 uex; // 厳密解 $u_{ex}$ を入れる

func f = 2*((x-x^2)+(y-y^2));
uex = x*(x-1)*y*(y-1);

for (int i=1; i<10; i++) {
  mesh Th=square(10*i,10*i); //メッシュを細かくする
  fespace Vh(Th,P1); //有限要素空間
  Vh u,v; // $u,v\in V_{h}(\Omega)$
  solve Poisson(u,v) =
      int2d(Th)(  dx(u)*dx(v) + dy(u)*dy(v)) 
    - int2d(Th) ( f*v )  //$-\int_{\Omega}fv dxdy$
    + on(1,2,3,4,u=0)  ;
  //$u$ と $u_{ex}$ を重ねる
  plot(u,uex,cmm="i="+i,wait=1);
}