/****************************************************************************
   * DESCRIPTION:  
   *   Serial HEAT2D Example - C Version
   *   This example is based on a simplified 
   *   two-dimensional heat equation domain decomposition.  The initial 
   *   temperature is computed to be high in the middle of the domain and 
   *   zero at the boundaries.  The boundaries are held at zero throughout 
   *   the simulation.  During the time-stepping, an array containing two 
   *   domains is used; these domains alternate between old data and new data.
   *
   *  The physical region, and the boundary conditions, are suggested
      by this diagram;
  
                     u = 0
               +------------------+
               |                  |
      u = 100  |                  | u = 100
               |                  |
               +------------------+
                     u = 100
  
  Interrior point :
    u[Central] = (1/4) * ( u[North] + u[South] + u[East] + u[West] )
  
   ****************************************************************************/

  #include <stdio.h> 

  #include <stdlib.h>
  #include <math.h>

  #define NN 50
  #define MM 50  

  
  #define ITER_PRINT 100

  #define _MAX_ITER 400 

  #define PRINT_DATA 1

  #define _EPSILON 0.01

  float update(int nx,int ny, float *u, float *unew);

  void inidat(int nx, int ny, float *u, float *unew); 
  void prtdat(int nx, int ny, float *u,const char *fnam);
  
  
  
  
  int main(int argc, char *argv[])
  {

      int N=NN,M=MM;

      float EPSILON=_EPSILON;

      int MAX_ITER= _MAX_ITER;

      if(argc !=4)
      {
          fprintf(stderr,"usage %s N EPSILON MAX_ITER
   ", argv[0]);
          fprintf(stderr,"\t\twhere N is GRID size, EPSILON is  Tolerance, MAX_ITER is max iteration
  "); 
          fprintf(stderr,"\t\texample N = 100, EPSILON = 0.1, MAX_ITER=1000
  "); 
          return -1;
      } 

      N = M = atoi(argv[1]);
      EPSILON = atof(argv[2]);
      MAX_ITER = atoi (argv[3]);

  
      float diff=1.0;
  
      float *u     = (float *)malloc(N * M * sizeof(float));
      float *unew  = (float *)malloc(N * M * sizeof(float));

      if(u==0 || unew ==0)
      {
          perror("Can't allocated data
  ");
          return -1;
      }

      printf ( "
  " );
      printf ( "HEATED_PLATE
  " );
      printf ( "  Serial version
  " );
      printf ( "  A program to solve for the steady state temperature distribution
  " );
      printf ( "  over a rectangular plate.
  " );
      printf ( "
  " );
      printf ( "  Spatial grid of %d by %d points.
  ", M, N );
      printf ( "  The iteration will end until tolerance <= %f
  
  ",EPSILON);

  
      /* Initialize grid and create input file */
      printf("Initializing grid
  ");

      inidat(N, M,u,unew);
  
      prtdat(N, M,u, "initial.dat");

      printf("Start computing
  
  ");

  
      int iter=0;
  
      /* 
       *   iterate until the  new solution unew differs from the old solution u
       *     by no more than EPSILON.
       *     */

      while(diff> EPSILON && iter <MAX_ITER) {
  
          diff= update(N, M, u, unew);

          if(iter%ITER_PRINT==0)

              printf("Iteration %d, diff = %e
   ", iter,diff);

  
          iter++;
      }

      if(diff>EPSILON)
         printf("***Not converged MAX_ITERATION %d reached
  ***", MAX_ITER);

      prtdat(N, M, u, "final.dat");

      free(u);
      free(unew);
  }
  
  
  
  /****************************************************************************
   *  subroutine update
   ****************************************************************************/

  float update(int nx,int ny, float *u, float *unew)

  {
      int ix, iy;

      float diff=0.0;

  
      for (ix = 1; ix < nx-1; ix++) {
          for (iy = 1; iy < ny-1; iy++) {

   
              unew[ix*ny+iy] = (u[(ix+1)*ny+iy] +  u[(ix-1)*ny+iy] + u[ix*ny+iy+1] +  u[ix*ny+iy-1] )*0.25;
              diff = fmax( diff, fabs(unew[ix*ny+iy] - u[ix*ny+iy]));

          }
  
      }

      /**
       * COPY OLD DATA
       */

      for (ix = 1; ix < nx-1; ix++) {
          for (iy = 1; iy < ny-1; iy++) {
              u[ix*ny+iy] = unew[ix*ny+iy]; 
          }
      }   

      return diff;

  }
  
  /*****************************************************************************
   *  Initialize Data
   *****************************************************************************/
  void inidat(int nx, int ny, float *u, float *unew) 
  {
      int ix, iy;
  
      /*
       *Set boundary data and interrior values
       * */

      // interior points
      for (ix = 1; ix < nx-1; ix++) 
          for (iy = 1; iy < ny-1; iy++) { 

              u[ix*ny+iy]=0.0; 

          }

      //boundary left

      for (ix = 1; ix < nx-1; ix++){ 
          u[ix*ny]=100.0; 
  
      }
  
      //boundary right
      for (ix = 1; ix < nx-1; ix++){ 
          u[ix*ny+ (ny-1)]=100.0; 
  
      }
  
      //boundary down

      for (iy = 0; iy < ny; iy++){ 

          u[(nx-1)*(ny)+iy]=100.0; 
  
      }
  
      //boundary top

      for (iy = 0; iy < ny; iy++){ 

          u[iy]=100.0; 

  
      }

  }

  /**************************************************************************
   * Print Data to files
   **************************************************************************/

  void prtdat(int nx, int ny, float *u,const char *fname)

  {

      int ix, iy;
      FILE *fp;
  
      if(ITER_PRINT==0)return;

      fp = fopen(fname, "w");

      // fprintf ( fp, "%d
  ", M );
      // fprintf ( fp, "%d
  ", N );

  
      for (ix = 0 ; ix < nx; ix++) {
          for (iy =0; iy < ny; iy++) {

              fprintf(fp, "%6.2f ", u[ix*ny+iy]);

          }

          fputc ( '
  ', fp);

      }

      printf ("  Data written to the output file %s
  ", fname);

      fclose(fp);
  }