/*  heat.c

    Chris Nevison
	June 24, 1997

	This file contains functions for setting up and solving the
	heat equation on a rectangular surface, using an 
	iterative method.
	The boundary values of the rectangle are held constant.

*/

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "heat.h"

int CreateDblMat(int numrows, int numcols, double initval,
                 double * * * mat)
/* pre:  *mat is uninitialized or NULL                  */
/* post: *mat is a numrows x numcols matrix of doubles  */
/*       with value initval.                            */
/*       if there is insufficient memory, return 0      */
/*       otherwise return 1                             */
{
  int r, c;
  if(*mat = (double **) malloc(numrows * sizeof(double *))){
    for(r = 0; r < numrows; r++){
      if(!((*mat)[r] = 
          (double *) malloc(numcols * sizeof(double))))
        return 0;
    }
    for(r = 0; r < numrows; r++)
      for(c = 0; c < numcols; c++)
        (*mat)[r][c] = initval;
    return 1;
  }
  else{
    return 0;
  }
}

void UpdateMat(double * * matin, int numrows, int numcols,
               double * * matout, double * maxdif)
/* pre:  matin and matout are initialized numrows x numcols */
/*       arrays of doubles.                                 */
/* post: matout contains values such that the boundary      */
/*       values are the same as the boundary values for     */
/*       matin and interior values are the average of the   */
/*       four (or eight) nearest neighbors.  maxdif is the  */
/*       maximum difference between corresponding values    */
/*       of the two arrays.                                 */
{
  int r, c;
  double dif;
  *maxdif = 0.0;
  for(r = 1; r < numrows-1; r++){
    for(c = 1; c < numcols-1; c++){
      matout[r][c] = (matin[r-1][c] + matin[r+1][c] +
                      matin[r][c-1] + matin[r][c+1] +
                      matin[r-1][c-1] + matin[r-1][c+1] +
                      matin[r+1][c-1] + matin[r+1][c+1]
                      ) / 8;
      dif = fabs(matout[r][c] - matin[r][c]);
      if(dif > *maxdif)
        *maxdif = dif;
    }
  }
}

void SolveHeat(double * row0, double * rowmax,
               double * col0, double * colmax, 
               int numrows, int numcols,
               double tol, int maxiter,
               double * * * mat)
/* pre:  row0 and rowmax each have numcols entries,          */
/*       col0 and colmax each have numrows entries,          */
/*       row0[0] = col0[0]; row0[numcols-1] = colmax[0];     */
/*       rowmax[0] = col0[numrows-1];                        */
/*       rowmax[numcols-1] = colmax[[numrows-1];             */
/*       tol > 0; maxiter > 0;                               */
/*       *mat is uninitialized or NULL.                      */
/* post: *mat is a numrows x numcols matrix whose entries    */
/*       of the heat equation with boundary values fixed at  */
/*       the values given by row0, rowmax, col0, and colmax  */
/*       computed using the ?????????? iterative method      */
/*       with accuracy given by tol unless maxiter iterations*/
/*       are reached.                                        */
{
  int r, c, iter;
  double avg, maxdif;
  double * * mat2;

  avg = 0.0;
  for(r = 0; r < numrows; r++)
    avg += col0[r] + colmax[r];
  for(c = 0; c < numcols; c++)
    avg += row0[c] + rowmax[c];
  avg /= (double)(2 * numrows + 2 * numcols - 4);

  if(!(CreateDblMat(numrows, numcols, avg, mat))){
    printf("Not enough memory to create matrix in SolveHeat\n\n");
    exit(-1);
  }
  if(!(CreateDblMat(numrows, numcols, avg, &mat2))){
    printf("Not enough memory to create matrix 2 in SolveHeat\n\n");
    exit(-1);
  }

  for(c = 0; c < numcols; c++){
    (*mat)[0][c] = row0[c];
    (*mat)[numrows-1][c] = rowmax[c];
    mat2[0][c] = row0[c];
    mat2[numrows-1][c] = rowmax[c];
  }
  for(r = 0; r < numrows; r++){
    (*mat)[r][0] = col0[r];
    (*mat)[r][numcols-1] = colmax[r];
    mat2[r][0] = col0[r];
    mat2[r][numcols-1] = colmax[r];
  }
 
  maxdif = tol + 1.0;
  iter = 0;
  while((maxdif > tol) && (iter < maxiter)){
    UpdateMat(*mat, numrows, numcols, mat2, &maxdif);
    UpdateMat(mat2, numrows, numcols, *mat, &maxdif);
    iter += 2;
    printf(".");
  }
  printf("\n");
  if(iter == maxiter)
    printf("number of iterations at limit on SolveHeat!\n");
}

void FindMinMax(int numrows, int numcols, double * * mat,
                double * minval, double * maxval)
/* pre:  mat is a numrows x numcols array of doubles         */
/* post: minval and maxval are the minimum and maximum       */
/*       values, respectively, in mat.                       */
{
  int r, c;
  *minval = mat[0][0];
  *maxval = mat[0][0];

  for(r = 0; r < numrows; r++)
    for(c = 0; c < numcols; c++){
      if(mat[r][c] < *minval)
        *minval = mat[r][c];
      if(mat[r][c] > *maxval)
        *maxval = mat[r][c];
  }
}

void ReadHeatBndry(FILE * infile,
                   int * numrows, int * numcols,
                   double * * row0, double * * rowmax,
                   double * * col0, double * * colmax)
/* pre:  infile is open and ready for reading, formatted to  */
/*       contain the number of rows and number of columns    */
/*       for the rectangular region (two ints), followed by  */
/*       the initial values for the first and last rows and  */
/*       the first and last columns of the regions           */
/*       (2 * numcols + 2 * numrows double values)           */
/* post: The number of rows and number of columns have been  */
/*       assigned to the parameters numrows and numcols;     */
/*       row0 and rowmax are each allocated space for        */
/*       numcols double values and col0 and colmax are       */
/*       each allocated space for numrows double values;     */
/*       the first and last row values have been assigned to */
/*       parameters row0 and rowmax, and the first and last  */
/*       column values have been assigned to the parameters  */
/*       col0 and colmax.                                    */
/*       row0[0] = col0[0]; row0[numcols-1] = colmax[0];     */
/*       rowmax[0] = col0[numrows-1];                        */
/*       rowmax[numcols-1] = colmax[numrows-1]               */
{
  int k;
  fscanf(infile, "%d %d", numrows, numcols);

  *row0   = (double *) malloc(*numcols * sizeof(double));
  *rowmax = (double *) malloc(*numcols * sizeof(double));
  *col0   = (double *) malloc(*numrows * sizeof(double));
  *colmax = (double *) malloc(*numrows * sizeof(double));

  for(k = 0; k < *numcols; k++)
    fscanf(infile, "%lf", &((*row0)[k]));
  for(k = 0; k < *numcols; k++)
    fscanf(infile, "%lf", &((*rowmax)[k]));
  for(k = 0; k < *numrows; k++)
    fscanf(infile, "%lf", &((*col0)[k]));
  for(k = 0; k < *numrows; k++)
    fscanf(infile, "%lf", &((*colmax)[k]));
}