// define el borde del dominio
border C1(t=-1,1){x=t;y=-1;label=1;}
border C2(t=-1,1){x=1;y=t;label=2;}
border C3(t=-1,1){x=-t;y=1;label=3;}
border C4(t=-1,1){x=-1;y=-t;label=4;}

// declaramos una malla, y definimos un espacio asociado a la malla
mesh Th; 
fespace Vh(Th,P1);

// datos del problema
func f = 2*pi^2*sin(pi*x)*sin(pi*y);

// solución exacta (para poder calcular los errores)
func u = sin(pi*x)*sin(pi*y);
func dxu = pi*cos(pi*x)*sin(pi*y);
func dyu = pi*sin(pi*x)*cos(pi*y);

// declaración de variables
int n=1;
real errL2old,errH1old,errL2,errH1,ordenL2,ordenH1,hold;
real nv,nvold;
Vh uh,v;


// definición del problema
problem poisson(uh,v) = int2d(Th)(dx(uh)*dx(v) + dy(uh)*dy(v)) 
   - int2d(Th)(f*v) + on(1,2,3,4,uh=0); 


// crea un archivo, si ya existe borra el contenido
{ ofstream file("errores.txt");
file << "Hmax" << "   " << "#nodos" << "   " << "err L2" << "   "
   << "err H1" << "   " << "orden L2" << "   " << "orden H1" << endl;
}


// iteración
for(int i=1;i<10;i++){
n=2*n;

Th = buildmesh(C1(n)+C2(n)+C3(n)+C4(n)); // malla no estructurada
// Th = square(n,n,[2*x-1,2*y-1]); // malla estructurada
nv=Th.nv; // numero de vértices

poisson;   // resulve el problema

// calcula los errores
errL2 = sqrt(int2d(Th, qforder=6)((u-uh)^2));
errH1 = errL2 + sqrt(int2d(Th, qforder=6)((dxu-dx(uh))^2+(dyu-dy(uh))^2));

// estima el orden de convergencia
if(i>1){
// orden respecto del número de vértices 
// (en general, es la mitad del orden respecto de hmax)
ordenL2 = -(log(errL2)-log(errL2old))/(log(nv)-log(nvold));
ordenH1 = -(log(errH1)-log(errH1old))/(log(nv)-log(nvold));

// orden respecto de hmax 
// ordenL2 = -(log(errL2)-log(errL2old))/(log(Th.hmax)-log(hold));
// ordenH1 = -(log(errH1)-log(errH1old))/(log(Th.hmax)-log(hold));
}

// guarda los datos actuales
errL2old = errL2;
errH1old = errH1;
hold = Th.hmax;
nvold = Th.nv;


// escribe los datos en el archivo errores.txt
{ ofstream file("errores.txt", append);
if(i>1)
file << Th.hmax << "   " << Th.nv << "   " << errL2 << "   "
   << errH1 << "   " << ordenL2 << "   " << ordenH1 << endl;
else
file << Th.hmax << "   " << Th.nv << "   " << errL2 << "   "
   << errH1 << "   " << 0 << "          " << 0 << endl; 
}
}

// dibuja la última malla y última solución (discreta)
// si la malla es muy fina, puede tardar un poco
plot(Th);
plot(uh);