/*
* 2D HEAT EQUATION SOLVER
*
* Copyright (c) 2008 Tony Saad <saadtony@gmail.com>
* (Written by Tony Saad <saadtony@gmail.com>)
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
*/
/*--------------------------------------------------------*/
/* INCLUDE HEADERS
/*--------------------------------------------------------*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <ctime>
/*--------------------------------------------------------*/
/* STRUCTURES
/*--------------------------------------------------------*/
struct ProblemInfo {
 int numberOfNodesX;
 int numberOfNodesY;
 double endSimulationTime;
 double thermalDiffusivity;
 double CFL;
};
/*--------------------------------------------------------*/
/* FUNCTION PROTOTYPES
/*--------------------------------------------------------*/
bool processUserInput(ProblemInfo &theProblemInfo);
void setupBoundaryConditions(double** theArray, int arraySizeX, int arraySizeY);
void initializeInteriorNodes(double** theArray, int arraySizeX, int arraySizeY);
void initializeArray(double** theArray, int arraySizeX, int arraySizeY);
int* makeIntArray(int arraySize);
int** make2DIntArray(int arraySizeX, int arraySizeY);
double* makeDoubleArray(int arraySize);
double** make2DDoubleArray(int arraySizeX, int arraySizeY);
void saveToFile(ProblemInfo theProblemInfo, double** theData,char* theFileName);

/*--------------------------------------------------------*/
/* MAIN PROGRAM
/*--------------------------------------------------------*/
int main(int argc, char* argv[]) {

 /***********************************************************/
 /* VARIABLES
 /***********************************************************/
 int nNodesX,
 nNodesY;

 int ix = 0,
 jy =0;

 double endSimulationTime = 1.0,
 CFL = 0.25,
 deltaT,
 deltaX,
 deltaY,
 deltaXSquared,
 deltaYSquared;

 clock_t startCPUTime = 0,
 endCPUTime = 0,
 totalCPUTime = 0;

 ProblemInfo probInfo;

 /* process user input */
 if (!processUserInput(probInfo)) return 0;

 startCPUTime = clock();

 nNodesX = probInfo.numberOfNodesX;
 nNodesY = probInfo.numberOfNodesY;
 deltaX =(double) 1.0/probInfo.numberOfNodesX;
 deltaY =(double) 1.0/probInfo.numberOfNodesY;
 deltaXSquared = deltaX*deltaX;
 deltaYSquared = deltaY*deltaY;

 CFL = probInfo.CFL;
 deltaT = CFL*(deltaXSquared*deltaYSquared)/((deltaXSquared+deltaYSquared)*probInfo.thermalDiffusivity);

 /* memory allocation */
 double** uOld = make2DDoubleArray(nNodesX,nNodesY);
 initializeArray(uOld,nNodesX,nNodesY);

 double** uNew = make2DDoubleArray(nNodesX,nNodesY);
 setupBoundaryConditions(uNew,nNodesX,nNodesY);

 /***********************************************************/
 /* COMPUTATION
 /***********************************************************/
 int iter = 0;
 while(iter++*deltaT <= probInfo.endSimulationTime) {

 if (iter % 100 == 0) {
 printf("Iteration: %i - Current Time = %lf s\n", iter-1, (iter-1)*deltaT);
 }

 /* calculate the new values of the dependant variable */
 for (ix = 1; ix < (nNodesX-1); ix++) {
 for (jy = 1; jy < (nNodesY-1) ;jy++){
 uNew[ix][jy] = uOld[ix][jy] + CFL*(uOld[ix+1][jy] + uOld[ix-1][jy] - 4*uOld[ix][jy] + uOld[ix][jy-1] + uOld[ix][jy+1]);
 }
 }

 /* set uOld = uNew */
 for (ix = 1; ix < (nNodesX-1); ix++) {
 for (jy = 1; jy < (nNodesY-1) ;jy++){
 uOld[ix][jy] = uNew[ix][jy];
 }
 }
 }

 endCPUTime = clock();

 /* compute the total wall time */
 totalCPUTime= endCPUTime-startCPUTime;

 /* compute the maximum wall time */
 printf("Wall clock time %lf s\n",(double) totalCPUTime/CLOCKS_PER_SEC);

 saveToFile(probInfo,uNew,"C:\\2D Heat Equation.dat");

 /* works only in Windows */
 system("pause");
 return 0;
}

/*--------------------------------------------------------*/
/* FUNCTIONS
/*--------------------------------------------------------*/
bool processUserInput(ProblemInfo &theProblemInfo) {
 printf("Enter the number of nodes in the x direction:");
 fflush(stdout);
 scanf("%d",&theProblemInfo.numberOfNodesX); 

 printf("Enter the number of nodes in the y direction:");
 fflush(stdout);
 scanf("%d",&theProblemInfo.numberOfNodesY); 

 if (theProblemInfo.numberOfNodesX == 0 || theProblemInfo.numberOfNodesY==0) return 0;

 printf("Enter thermalDiffusivity:");
 fflush(stdout);
 scanf("%lf",&theProblemInfo.thermalDiffusivity); 

 printf("Enter end simulation time:");
 fflush(stdout);
 scanf("%lf",&theProblemInfo.endSimulationTime); 

 printf("Enter the Courant-Friedrichs-Lewy (CFL) number \n (For stability, the CFL number should be less than 0.25 \n based on a different number of nodes in the X and Y directions):");
 fflush(stdout);
 scanf("%lf",&theProblemInfo.CFL);

 return 1;
}

void setupBoundaryConditions(double** theArray,int arraySizeX, int arraySizeY) {
 /* set boundary conditions for ix, jy = 0 and ix, jy = n-1 */
 /* Here, we are using Dirichlet conditions */
 int ix, jy;
 double leftBC = 10,
 rightBC = 10,
 topBC = 20,
 bottomBC = 20;

 /* setup the bottom and top BCs, jy = 0 and jy = n-1 or arraySizeY - 1 */
 for (ix = 0; ix < arraySizeX; ix++) {
 theArray[ix][0] = bottomBC; //bottom BC
 theArray[ix][arraySizeY-1] = topBC; //top BC
 } 

 /* setup the left and right BCs, ix = 0 and ix = arraySizeX - 1 */
 for (jy = 0; jy < arraySizeY; jy++) {
 theArray[0][jy] = leftBC; //left BC
 theArray[arraySizeX-1][jy] = rightBC; //right BC
 }

 /* set the values at the corner nodes as averages of both sides*/
 // bottom left
 theArray[0][0] = (leftBC + bottomBC)*0.5;
 // top left
 theArray[0][arraySizeY-1] = 0.5*(topBC + leftBC);
 // top right
 theArray[arraySizeX-1][arraySizeY-1] = 0.5*(topBC + rightBC);
 // bottom right
 theArray[arraySizeX-1][0] = 0.5*(bottomBC + rightBC);
}

void initializeInteriorNodes(double** theArray, int arraySizeX, int arraySizeY) {
 int ix, jy;
 for (ix = 1; ix <= arraySizeX - 2; ix++) {
 for (jy = 1; jy <= arraySizeY - 2; jy++) {
 theArray[ix][jy] = 0.0;
 }
 }
}

void initializeArray(double** theArray, int arraySizeX, int arraySizeY) {
 setupBoundaryConditions(theArray,arraySizeX,arraySizeY);
 initializeInteriorNodes(theArray,arraySizeX,arraySizeY);
}

double** make2DDoubleArray(int arraySizeX, int arraySizeY) {
 double** theArray;
 theArray = (double**) malloc(arraySizeX*sizeof(double*));
 for (int ix = 0; ix < arraySizeX; ix++) {
 theArray[ix] =(double*) malloc(arraySizeY*sizeof(double));
 }
 return theArray;
}

double* makeDoubleArray(int arraySize) {
 double* theArray;
 theArray = (double*) malloc(arraySize*sizeof(double));
 return theArray;
}

int* makeIntArray(int arraySize) {
 int* theArray;
 theArray = (int*) malloc(arraySize*sizeof(int));
 return theArray;
}

int** make2DIntArray(int arraySizeX, int arraySizeY) {
 int** theArray;
 theArray = (int**) malloc(arraySizeX*sizeof(int*));
 for (int ix = 0; ix < arraySizeX; ix++) {
 theArray[ix] =(int*) malloc(arraySizeY*sizeof(int));
 }
 return theArray;
}

void saveToFile(ProblemInfo theProblemInfo, double** theData, char* theFileName) {
 int nNodesX = theProblemInfo.numberOfNodesX;
 int nNodesY = theProblemInfo.numberOfNodesY;
 double deltaX =(double) 1.0/theProblemInfo.numberOfNodesX;
 double deltaY =(double) 1.0/theProblemInfo.numberOfNodesY;
 FILE* theDataFile;
 theDataFile = fopen(theFileName,"w");
 for (int ix = 0; ix < nNodesX; ix++) {
 for (int jy = 0; jy < nNodesY ;jy++){
 fprintf(theDataFile,"%lf %lf %lf\n",ix*deltaX,jy*deltaY,theData[ix][jy]);
 }
 fprintf(theDataFile,"\n");
 }
 fclose(theDataFile);
}