/////////////////////////////////
// Resolution de l'equation    //
// de la chaleur dans un carre //
/////////////////////////////////

//Definition du maillage
mesh Th;

//Nombre de mailles suivant x et y
int Nbnoeuds=10;

Th=square(Nbnoeuds,Nbnoeuds,[x,y]);

//Solution exacte
func u=x*(1.-x)*y*(1.-y);

//Fonction de x et de y
func f=2.*x*(1-x)+2.*y*(1-y);

//Definition de l'espace des elements finis P1 associe
//au maillage Th
fespace Vh(Th,P1);

//uh et vh sont des elements de Vh
Vh uh,vh,erreur;


//Definition du probleme variationnel
problem chaleur(uh,vh,solver=LU)=
	int2d(Th)(dx(uh)*dx(vh)+dy(uh)*dy(vh))
	-int2d(Th)(f*vh)
	+on(1,2,3,4,uh=0)
;




int Nbiter=10;
//Resolution du probleme
chaleur;
erreur=u-uh;

//Calcul de l'erreur L2
real erreurL2=sqrt(int2d(Th)(erreur^2));
real erreurH10=sqrt(int2d(Th)(dx(erreur)^2+dy(erreur)^2));

//On affiche le resultat
string legende="Nombre de noeuds="+Nbnoeuds
		+", Erreur L2 ="+erreurL2
		+", Erreur H10="+erreurH10;
plot(uh,wait=1,cmm=legende);


for (int i=0;i<Nbiter;i++){
	Nbnoeuds=Nbnoeuds+10;
	Th=square(Nbnoeuds,Nbnoeuds,[x,y]);
	chaleur;
	erreur=u-uh;
	erreurL2=sqrt(int2d(Th)(erreur^2));
	erreurH10=sqrt(int2d(Th)(dx(erreur)^2+dy(erreur)^2));
	legende="Nombre de noeuds="+Nbnoeuds
		+", Erreur L2 ="+erreurL2
		+", Erreur H10="+erreurH10;
	plot(uh,wait=1,cmm=legende);
	}