/***************************************************/
/* Tarefa metodos EDP parabolica                   */
/* Dissipacao Calor na Placa                       */
/* Metodo Iterativo: Gauss Seidel				   */
/* Metodo no tempo e espaco: Crank Nicolson        */		   
/*												   */
/* Andreza Lukosiunas							   */
/* MAP-2320 Prof Alexandre M. Roma				   */
/***************************************************/


#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define xMax 4
#define yMax 4
#define tMax 4
#define phi 3.1415927
//#define mi 0.04
//#define mi 0.5
//#define mi 1.0
#define mi 2.5


/* Prototipo das funcoes **************************/

double U(double T, double X, double Y);
double f(double T, double X, double Y);
double ci(double ti, double x, double y);
double cc(double t, double x, double y, int a);
double min(double dx, double dy);
void *mallocSafe (size_t nbytes);
/**************************************************/

int main()
{
	
	FILE *arquivo, *file;
	if ((arquivo = fopen("/Users/andreza/Desktop/euler1.txt", "w")) == NULL)
	{
		printf("\nProblemas na abertura do arquivo da tabela\n");
		exit(-1);
	}
	fprintf(arquivo, "q	h		erro global	e1/e2		log2\n");
	
	if ((file = fopen("/Users/andreza/Desktop/matriz.txt", "w")) == NULL)
	{
		printf("\nProblemas na abertura do arquivo da tabela\n");
		exit(-1);
	}
	
	
	/*********** Variaveis **************/
	double dx, dy, dt, ti, tf, t, xi, xf, yi, yf; 
	double **UE, **u, **u_a, u_old;
	double *x, *y;
	double e1, e2, eAux, r, s, error, tol, diag, eMAX = 0.0;
	int q, N, i, j, k, m, n, cont = 1, gaussIt;
	double A, B, C, D, E;
	
	/*********** Inicializando **********/
	e2 = 0.0; e1 = 0.0; ti = 0.0; tf = 1.0; xi = 0.0; xf = 1.0; yi = 0.0; yf = 1.0;
	t = ti; 
	
	/* Refina Malha */
	for (q=0; q<5; q++) 
	{
		t = ti;
		N = pow(2,q);/* Refina malha em potencia de 2 */
		dx = (xf-xi)/(N*xMax);
		dy = (yf-yi)/(N*yMax);
		dt = min( dx, dy); 
		tol = 0.1*dt;
		
		diag = 1 + mi*dt/(dx*dx) + mi*dt/(dy*dy);
		
		/* Alocar memoria *********************************/
		UE = mallocSafe((N*yMax + 1)*sizeof(double*));
		for (k = 0; k < (N*yMax + 1); k++) 
		{
			UE[k] = mallocSafe((N*xMax + 1)*sizeof(double));
		}
		/***/
		u = mallocSafe((N*yMax + 1)*sizeof(double*));
		for (k = 0; k < (N*yMax + 1); k++) 
		{
			u[k] = mallocSafe((N*xMax+ 1)*sizeof(double));
		}
		/***/
		u_a = mallocSafe((N*yMax + 1)*sizeof(double*));
		for (k = 0; k < (N*yMax + 1); k++) 
		{
			u_a[k] = mallocSafe((N*xMax+ 1)*sizeof(double));
		}
		/***/
		x = mallocSafe((N*xMax + 1)*sizeof(double*));
		y = mallocSafe((N*yMax + 1)*sizeof(double*));
		/**************************************************/
		
		/* Preenche particao x e y ************************/
		for (i = 0; i < (N*xMax + 1); i++) 
		{
			x[i] = xi + i*dx;
		}
		for (j = 0; j < (N*yMax + 1); j++) 
		{
			y[j] = yi + j*dy;
		}
		/**************************************************/
		
		/* Condicoes Iniciais *****************************/
		for (i = 0; i < (N*xMax + 1); i++) 
		{
			for (j = 0; j < (N*yMax + 1); j++) 
			{
				u_a[i][j] = ci(ti, x[i], y[j]);
				u[i][j] = u_a[i][j];
				
			}
		}
		/**************************************************/
		
		
		/* Método de Crank Nicolson com Gauss_Seidel *****/
		
		while (t <= tf) 
		{
			if (t == tf) {
				break;
			}
			
			t = ti + dt*cont;
			//t = t + dt;
			cont++;
			
			if (t > tf) {
				t = tf;
			}
			
			/* Condicoes de contorno */
			for(m=0; m<(N*yMax + 1); m++)
			{ 
				u[0][m] = cc(t, xi, y[m], 1);   
				u[(N*xMax)][m] = cc(t, xf, y[m], 2);   
			}
			for(n=0; n < (N*xMax + 1); n++)
			{ 
				u[n][0] = cc(t, x[n], yi, 3);     
				u[n][(N*yMax)] = cc(t, x[n], yf, 4);      
			}
			
			gaussIt = 0;
			/*Gauss-Seidel Iterations*/
			do
			{ 
				eMAX =0.0; error = 0.0; 				
				for(i=1; i<(N*xMax); i++) 
				{ 
					for(j=1; j<(N*yMax); j++)
					{ 
						u_old = u[i][j];
						
						A = 0.5*dt*(f(t,x[i],y[j])+f(t+dt,x[i],y[j]));
 						B = (0.5*dt*mi*(u_a[i+1][j]-2.0*u_a[i][j]+u_a[i-1][j]))/(dx*dx);
						C = (0.5*dt*mi*(u_a[i][j+1]-2.0*u_a[i][j]+u_a[i][j-1]))/(dy*dy);
						D = (0.5*dt*mi*(u[i+1][j]+u[i-1][j]))/(dx*dx);
						E = (0.5*dt*mi*(u[i][j+1]+u[i][j-1]))/(dy*dy);
		
						u[i][j] = (u_a[i][j] + A + B + C + D + E )/diag;
						
						error = fabs((u[i][j] - U(t, x[i], y[j])));
						
						if (error > eMAX) 
						{
							eMAX = error;
						}
					}
				}
				gaussIt++;
				
			}while(eMAX > tol && gaussIt < 1500);
			
			for(i=0; i<(N*xMax)+1; i++) 
			{ 
				for(j=0; j<(N*yMax)+1; j++)
				{ 
					u_a[i][j] = U(t, x[i], y[j]);
				}
			}
		}
		
		e2 = 0.0; eAux = 0.0;
		/* Solução Exata e Erro Global *******************/
		for (i = 0; i < (N*xMax + 1); i++)
		{
			for (j = 0; j < (N*yMax + 1); j++)
			{
				UE[i][j] = U(t, x[i], y[j]);
				
				/* Erro Global */
				eAux = fabs(u[i][j] - UE[i][j]); /* Norma do Maximo */
				if (eAux > e2) 
				{
					e2 = eAux;
				}
			}
		}
		/**********************************************************/
		printf("eAux: %lf t: %lf q: %d\n", e2, t, q);
		
		r=e1/e2;
		s=log2(r);
		e1=e2;
		
		/* Imprime Tabela Convergencia ******************/
		if (q==0) {
			fprintf(arquivo,"%d	%E	%E\n", q, dx, e2);
		}else{
			fprintf(arquivo,"%d	%E	%E	%E	%E\n", q, dx, e2, r, s);
		}
	}
	return 0;
}

/* Funcoes implementadas **********************************************/

/* Solucao Exata */
double U(double T, double X, double Y)
{
	double A;
	
	//A = sin(X)*cos(Y)*exp(-2*T);
	//A = (1 - Y)*exp(X + T);
	//A = 1 + cos(4*phi*X)*cos(3*phi*Y)*exp(-(phi*phi)*T); 
	//A = exp(-T)*(pow(X, 4) + pow(Y, 4));
	//A = exp(-T)*(cos(X)+sin(Y));
	A = cos(X + Y)*exp(-T);
	
	return A;
}

/* Termo forcante */
double f(double T, double X, double Y)
{
	double C;
	C = (-exp(-T)*cos(X+Y)) - (mi*exp(-T)*(-cos(T+Y)-cos(X+Y)));
	return C;
}

/* Condicao Inicial */
double ci(double ti, double x, double y)
{
	return U(ti, x, y);
	
}

/* Condicao de Contorno */
double cc(double t, double x, double y, int a)
{
	/* CC de xi: 1 e xf: 2 */
	if (a == 1) 
	{
		return U(t, x, y);
	} 
	
	if (a == 2) 
	{
		return U(t, x, y);
	} 
	
	/* CC de yi: 3 e yf: 4 */
	if (a == 3) 
	{
		return U(t, x, y);
	} 
	
	if (a == 4) 
	{
		return U(t, x, y);
	}
	
	return 0;
}

/* Minimo */
double min(double dx, double dy)
{
	if (dx <= dy) 
	{
		return dx;
	} 
	else {
		return dy;
	}
}

/* Alocar memoria */
void *mallocSafe (size_t nbytes)
{
	void *ptr;
	ptr = malloc (nbytes);
	if (ptr == NULL) {
		printf ("Socorro! malloc devolveu NULL!\n");
		exit (EXIT_FAILURE);
	}
	return ptr;
}
/***************************************************************************/