/*
* 2D HEAT EQUATION SOLVER
*
* Copyright © 2008 Tony Saad <remineur@gmail.com>
* (Written by Tony Saad <remineur@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);
}