/*IMPLEMENTAÇÃO DE UM METODO PARA RESOLVER UMA EDP PARABÓLICA - Equação do calor numa placa*
* Método: Euler Implícito no tempo e diferenças centradas no espaço
* Daniela Passos Maia Moura   nº USP. 5916438
*
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define MMAX 560
#define NMAX 560
#define min(a,b)  (((a) < (b)) ? (a) : (b))

/*Declaração das funções*/
double ue(double x,double y, double t);  /* solução exata, escolhida por mim*/
double ci(double x,double y, double t); /* função Condição Inicial*/
double cc(double x,double y, double t); /* função Condição de Contorno*/
double fe(double x, double y, double t); /*  função termo Forçante*/
void GS(double u[NMAX][MMAX], double u_a[NMAX][MMAX],double hx, double hy, double dt,double f[NMAX][MMAX], int NX, int NY, double t, double x[NMAX], double y[MMAX]); /* resolve o sistema linear pelo método *de Gauss-Seidel.*/
double norma_s(double u[NMAX][MMAX], int NX, int NY, double x[NMAX], double y[MMAX], double t);
double norma2(double u[NMAX][MMAX], int NX, int NY, double x[NMAX], double y[MMAX],double hx, double hy, double t, double A);

/*Programa principal*/
int main()
{

  /*Discretização do Domínio*/
  int NX=8, NY=8;
  int n,i,j;

  double u[NMAX][MMAX],u_a[NMAX][MMAX];      /*matriz atual e matrix no tempo anterior respectivamente*/
  double hx, hy;
  double a1 = 0.0, b1 = 1.0, a2 = 0.0, b2 =1.0; 
  double x[NMAX], y[MMAX]; 
  double f[NMAX][MMAX];          /* vetor termo forçante*/

  /* Arquivos de saída */
  FILE *saida;      /*Arquivo auxiliar para testes*/
  FILE *tabela;     /* Tabela de conergência */
  tabela = fopen("Tabela de Convergencia.txt","w");
  saida = fopen ("Testes.txt","w");
  if (tabela == NULL)
    {
     printf ("Nao foi possivel abrir o arquivo!\n");
     exit (1);
     }
  if (saida == NULL)
   {
    printf ("Nao foi possivel abrir o arquivo!\n");
    exit (1);
    }
    
    /*Cabeçalho do arquivo*/
   fprintf(tabela, "\nDaniela Passos Maia Moura\nMAP2320 - Métodos Numéricos em Equações diferenciais\n\nTABELA DE CONVERGENCIA\n\n");
   fprintf(tabela,"\n           Erro_sup              razão               Erro2                razao \n\n");


   for(n=0; n<5; n++) /* Controla o refinamento da malha.*/
    { 
        
      /*Malha espacial - localização das incógnitas nos nós da malha*/
      hx = (b1 - a1)/NX;
      hy = (b2 - a2)/NY; 

      for(j = 0; j <= NX; j++)
       {
        x[j] = a1 + j*hx;}
       
         
      for(i = 0; i <= NY; i++)
        {
         y[i] = a2+ i*hy;}
       
       
      /*Malha temporal*/
      double to = 0.0 ,tf = 0.5;  /* tempo inicial e final*/
      double t, dt;               /* variável dt: subdivisões do tempo */

      dt = min(hx*hx,hy*hy);
      //dt = min(hx,hy);
      t=to;

      /*Condições iniciais para t=to : 
      *Atribui a u[i][j] e a u_a[i][j] os valores da função ci em cada ponto x[i],y[j] da malha */  
      for(i = 0; i <= NY; i++)
      {
        for(j = 0; j <= NX; j++)
        {
          u_a[i][j] = ci(x[i],y[j],to); 
          u[i][j] = ci(x[i],y[j],to);
          
        }
      }

     /*Laço no tempo (Euler implícito)*/
     while(t<tf)
     { 
       dt=min(dt,tf-t); /* atualizando o passo no tempo*/               
       fprintf(saida,"\nt=%f ; dt= %f ",t ,dt);

       /* Montagem do vetor termo forçante f(x,y,t) a partir da função fe*/
       
       for(i = 0; i <= NY; i++){
         for(j = 0; j <= NX; j++){
           f[i][j] = fe(x[i],y[j],t+dt); 
          
          }
       }

       /*Condições de contorno de Dirichlet:
       /*Atribui valores da função cc à fronteira da malha espacial*/

       /*Na horizontal*/
       for(j = 0; j <= NX; j++)
        {
          u[0][j] = cc(x[0],y[j],t+dt); 
          u[NY][j] = cc(x[NY],y[j],t+dt); 
       
        }
       
       /*Na vertical*/
       for(i = 0; i <= NY; i++)
        {  
         u[i][0] = cc(x[i],y[0],t+dt); 
         u[i][NX] = cc(x[i],y[NX],t+dt); 
         
        }

        /*Resolução do sistema linear resultante da discretização pelo metodo de G-S*/
        GS(u, u_a, hx, hy, dt, f, NX, NY, t,x,y); 

        /*atualizando o passo no tempo*/
        t=t+dt; 
        
        /*Atualizando as matrizes*/
        for(i=0; i<=NY; i++)
         {
          for (j=0; j<=NX; j++)
            {
              u_a[i][j]= u[i][j]; 
             }
         }

    } // fim do while - euler implícito;
 
     /*Cálculo dos erros  (na norma infinito=erro_sup e na norma 2) */ 
             
     double A;            /*Área do domínio, necessário para o cálculo da norma2*/
     A= (b1-a1)*(b2-a2);
     
     double e_sup=0.0, e2=0.0;
     e2= norma2(u, NX, NY,x,y,hx,hy,t,A);
     e_sup=norma_s(u, NX, NY,x,y,t);
  
     /*Calculo da razão entre os erros das etapas n-1 e n */  
     double r_sup, r2;      /* razões entre os erros em cada etapa*/
     double e_sup_a, e2_a;    /* variáveis auxiliares para o cálculo das razões*/

     if (n>0)
     { r_sup = e_sup_a/e_sup;
       r2 = e2_a/e2;

       e_sup_a = e_sup;     /* variável auxiliar recebe o valor do erro na etapa n*/
       e2_a = e2;
      }

     else 
     { e_sup_a = e_sup;     /* variável auxiliar recebe o valor do erro na etapa n*/
       e2_a = e2;

       r_sup = e_sup_a/e_sup;
       r2 = e2_a/e2;
     }

      /*Impressão dos erros na tela e no arquivo Tabela*/ 
      printf("\n\ne1=%f\ne2=%f",e_sup, e2);
      fprintf(tabela,"n=%d:      %f            %f             %f         %f\n",n,e_sup,r_sup,e2,r2);
      
     /*Atualização do refinamento da malha*/
     NX=2*NX; 
     NY=2*NY;
    }                      /*fim do laço que refina a malha*/

    //system("PAUSE");

    return 0;
}


/*Funções auxiliares*/

double ue(double x,double y, double t)  /* solução exata para a validação*/
 { double u;
    u=(exp(-t))*cos(x+y);
   
   return u;
   }

double ci(double x, double y, double t) /* função Condição Inicial*/
  {
   double f;
   f=ue(x,y,t); 
   return f ;
  }

double cc(double x,double y,double t) /* função Condição de Contorno */
  {
   double f;
   f = ue(x,y,t); 
   return f;
       }

double fe(double x, double y, double t) /*função termo forçante*/
  {  
    double f,mi=0.5;
    f= -(exp(-t))*cos(x+y) + mi*2.0*((exp(-t))*cos(x+y));
    return f;
   }

double norma_s(double u[NMAX][MMAX], int NX, int NY, double x[NMAX], double y[MMAX], double t) /*Função norma sup*/
  {
   int i,j;
   double e = 0.0;          
   double e_aux=0.0; /* variável auxiliar*/


   for(i = 0; i <= NY; i++)
   {
      for(j = 0; j <= NX; j++)
       {
        e_aux = fabs((u[i][j]-ue(x[i],y[j],t))); 
        if(e_aux > e)  e=e_aux;
      }
    }
   return e;
   }

   double norma2 (double u[NMAX][MMAX], int NX, int NY, double x[NMAX], double y[MMAX], double hx, double hy , double t, double A)   /*Função norma 2*/
   { 
    int i,j;
    double e=0.0, soma=0.0;
   
  
    for(i = 0; i <= NY; i++){
       for(j = 0; j <= NX; j++){
          
          soma = soma + (pow((u[i][j] - ue(x[i],y[j],t)),2))*hx*hy;
          }
       }
     e = sqrt(soma);
     e = e/A;
     
     return e;
    }


    /* Esta função resolve um sistema linear pelo método de Gauss Seidel */
    void GS(double u[NMAX][MMAX],double u_a[NMAX][MMAX], double hx, double hy, double dt, double f[NMAX][MMAX], int NX, int NY, double t, double x[NMAX], double y[MMAX])
    {
     int itmax=1420, k=0, i,j;  /* numero maximo de iterações, contadores*/
     double mi=0.5,ax,ay,d,e=0.0;
     double aux=0.0, var_aux, var=1.0;  /* variação relativa para o criterio de parada e uma auxiliar para o calculo;*/ 
    
     ax=mi*dt/(hx*hx);
     ay=mi*dt/(hy*hy);
     d=1.0+(2.0*mi*dt*((1.0/(hx*hx)) + (1.0/(hy*hy))));
  
     //e= dt*dt*dt;   // tolerância maxima para o G-S quando dt=min{hx,hy}
     e=dt*dt;        // tolerância maxima para o G-S quando dt=min{hx²,hy²}

     while(var>=e && k<=itmax) /* laço que resolve o sistema*/
      { k=k+1;
        var=0.0;
   
        for(i=1; i<NY; i++)
        { for (j=1; j<NX; j++)
          { aux=u[i][j];
            u[i][j] = (1.0/d)*(u_a[i][j]+ dt*f[i][j] + ax*(u[i-1][j]+u[i+1][j]) + ay*(u[i][j-1]+u[i][j+1]));
              
                      
	    /* calculo da variação entre as iterações*/
            var_aux = fabs((u[i][j]-aux));
		          
	    /* calculando o erro entre a solução do GS e da solução exata no tempo t*/
	    //var_aux = fabs(u[i][j]-(ue(x[i],y[j],t+dt)));
         
            /* calculando a variação máxima entre as coordenadas*/
             if(var_aux >= var)                   var=var_aux;
               
           }
         }   /* fim do laço for*/
         
      }  /* fim do while*/

      printf("\ntotal de iteracoes %d, %f\n",k,var);
      
    }