// heatB.edp --- 熱方程式を後退Euler法 (陰解法) で解く
// https://m-katsurada.sakura.ne.jp/program/fem/heatB.edp
// 菊地文雄, 有限要素法概説, サイエンス社のPoisson方程式の問題の非定常版
int i,m=10;
real Tmax=1, tau=0.01, t;
// 次の2行の行頭の // を削除すると、m, tau, theta が実行時に入力できる
// cout << "m dt theta: "; cin >> m >> tau >> theta;
// cout << "m=" << m << ", tau=" << tau << ", theta=" << theta << endl;
mesh Th=square(m,m);
plot(Th,wait=true);
fespace Vh(Th,P1);
func f=1; func g1=0; func g2=0;
func u0=sin(pi*x)*sin(pi*y);
Vh u=u0, uold, v;
plot(u,cmm="t=0",wait=1);
problem heat(u,v,init=i)=
  int2d(Th)(u*v)+int2d(Th)(tau*(dx(u)*dx(v)+dy(u)*dy(v)))-int2d(Th)(uold*v)
 -int2d(Th)(tau*f*v)-int1d(Th,2,3)(tau*g2*v)
 +on(1,4,u=g1);
for (i=0;i<Tmax/tau;i++) {
  uold=u;
  t=(i+1)*tau;
  heat;
  plot(u,cmm="t="+t,wait=0); // ps="heat"+i+".ps" として保存することも可能
}
plot(u,wait=1);