#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include "mpi.h"
#include "header.h"

#ifndef FUNCTIONS_C
#define FUNCTIONS_C

/*
 * compilation of the used functions to solve the laplace equation.
 * variables are named consistently with the input parameters, thus
 * further no explanations are given considering these.
 * */

/*
 * this function performs a check whether
 * the directory 'char *in' is available or not.
 * */
void check_file(char *in, char *tags){
	FILE *file;
	if(0==strcmp(tags, "r")){
		if(NULL==(file=fopen(in, "r"))){
			printf("failed to open %s for reading.\n", in);
			abort();
		}
		fclose(file);
	}
	else if(0==strcmp(tags, "w")){
		if(NULL==(file=fopen(in, "w"))){
			printf("failed to open %s for writing.\n", in);
			abort();
		}
		fclose(file);
	}
	else if(0==strcmp(tags, "rw")){
		if(NULL==(file=fopen(in, "r"))){
			if(NULL==(file=fopen(in, "w"))){
				printf("failed to open %s for writing.\n", in);
				abort();
			}
		}
		else{
			printf("%s already exists.\n", in);
			abort();
		}
	}
	else if(0==strcmp(tags, "a")){
		if(NULL==(file=fopen(in, "a"))){
			printf("failed to open %s for appending\n", in);
			abort();
		}
		fclose(file);
	}
}

int map(int ix, int *x, int *y, int X, int Y){
	int i=0,j=0,k=-1;
	while(i<X && ix!=k){
		j=0;
		while(j<Y && ix!=k){
			k=j*X+i;
			j++;
		}
		i++;
	}
	if(k==ix){
		*(x)=i-1;
		*(y)=j-1;
		return 0;
	}
	else{
		printf("int map(%d,...,X=%d, Y=%d) failed to find coordinates\n", ix, X, Y);
		return 1;
	}
}
/*
 * check if (ix, iy) is set to a potential
 * */
int check_potential(int ix, int iy, int NPHI, pot *phi0){
	int n;
	for(n=0; n<NPHI; n++){
		if(ix>=phi0[n].pos[0] && ix<=phi0[n].pos[1] && iy>=phi0[n].pos[2] && iy<=phi0[n].pos[3])
			return 1;
	}
	return 0;
}

/*
 * check if (ix, iy) is inside of the grounded area
 */
int check_ground(int ix, int iy, int NGROUND, int **ground){
	int n;
	for(n=0; n<NGROUND; n++){
		if(ix>ground[n][0] && ix <ground[n][1] && iy>ground[n][2] && iy<ground[n][3])
			return 1;
	}
	return 0;
}


/*
 * average over grounded lattice sites. separate grounded areas are averaged independently.
 *
 * NOTE: only the boundary of the grounded sites is computed void jacobi() thus only these lattice sites
 * may contribute to the average. the interior of the grounded areas is set to this average.
 * */
void ground_it(int NGROUND, int **ground, int X, double *phi){
	int j, x, y, N=0;
	double buf;
	for(j=0; j<NGROUND; j++){
		for(N=0, buf=0., x=ground[j][0], y=ground[j][2]; x<=ground[j][1]; x++, N++)
			buf+=phi[y*X+x];

		for(x=ground[j][0], y=ground[j][3]; x<=ground[j][1]; x++, N++)
			buf+=phi[y*X+x];

		for(y=ground[j][2]+1, x=ground[j][0]; y<=ground[j][3]-1; y++, N++)
			buf+=phi[y*X+x];

		for(y=ground[j][2]+1, x=ground[j][1]; y<=ground[j][3]-1; y++, N++)
			buf+=phi[y*X+x];

		buf/=(double)N;
		for(x=ground[j][0]; x<=ground[j][1]; x++){
			for(y=ground[j][2]; y<=ground[j][3]; y++){
				phi[y*X+x]=buf;
			}
		}
	}
}

/*
 * Implementation of the jacobi algorithm to solve a 2d-laplace equation.
 * char *dif is the output directory for the precision and iteration-number.
 * char *logf is the output directory for the total required number of iterations,
 * the lattice volume of the solution, the total time for the solution and its actual
 * precision.
 * */
void jacobi(int X, int Y, double NP, double *phi, int NPHI,
		pot *phi0, int NGROUND, int **ground, char *dif, char *logf){
	double *dum, pres=1.+NP, buf=0.;
	int x, y, i, ix;
	double t1, t2;
	FILE *file;
	/*MPI_Wtime() returns the current time of the process in seconds with double-precision. meaning
	 * that this function allows for a much more precise estimate for the required time than any function
	 * defined in time.h where the highest precision is in seconds. For example consider a lattice of the size 8x8,
	 * the solution of this lattice is computed in the order of 10e-4 seconds.
	 * */
	t1 = MPI_Wtime();
	/*
	 * double *dum is used to store the new solution of the laplace equation.
	 * */
	dum=malloc(X*Y*sizeof(double));
	for(ix=0; ix<X*Y; ix++)
		dum[ix]=phi[ix];

	check_file(dif, "w");
	file=fopen(dif, "w");
	i=0;
	/*
	 * loop over either the precision (pres) if NP<1.
	 * or the iteration number (i) if NP>1.
	 * */
	while((i<(int)NP && NP>1.) || (pres>NP && NP<1.)){
		for(x=1; x<X-1; x++){/*loop over the volume of the lattice, excluding the boundary*/
			for(y=1; y<Y-1; y++){
				ix=y*X+x;
				/*
				 * check if (x,y) is inside of the grounded area or if it is
				 * predefined by a potential.
				 * */
				if(check_potential(x, y, NPHI, phi0)==0 && check_ground(x, y, NGROUND, ground)==0){
					dum[ix]=(phi[(y-1)*X+x]+phi[(y+1)*X+x]+phi[y*X+(x-1)]+phi[y*X+(x+1)])/4.;
					/*dum[ix]=(1.-r)*phi[ix]+r*0.25*(phi[(y-1)*X+x]+phi[(y+1)*X+x]+phi[y*X+(x-1)]+phi[y*X+(x+1)]);*/
				}
			}
		}/*end loop over lattice-volume
		average over the grounded area
		*/
		ground_it(NGROUND, ground, X, dum);
		buf=0.;/*buffer for the precision*/
		for(x=1; x<X-1; x++){
			for(y=1; y<Y-1; y++){
				ix=y*X+x;
				if(fabs(phi[ix]-dum[ix])>buf)
					buf=fabs(phi[ix]-dum[ix]);/*find the maximum absolute change from the previous step to this one*/
				phi[ix]=dum[ix];/*copy the new solution*/
			}
		}
		fprintf(file, "%d %.6e\n", i, buf);/*output for the current precision and iteration-number*/
		pres=buf;
		i++;
	}
	fclose(file);
	/*
	 * logfile output: lattice extent, method, volume, iteration-number, precision and total time.
	 * */
	check_file(logf, "a");
	file=fopen(logf, "a");
	fprintf(file, "###################################################\n");
	fprintf(file, "#X = %d\t Y = %d\t jacobi\n", X, Y);
	fprintf(file, "#\tX*Y\tN\tP\tR\tTime\n");
	t2 = MPI_Wtime();
	fprintf(file, "%8d %8d %.2e %.4e\n", X*Y, i-1, pres, t2-t1);
	fclose(file);
	free(dum);/*free memory required for the temporary array*/
}

/**
 * Implementation of the red-black algorithm for the solution of the 2d-laplace equation.
 * Since it is largely the same as the jacobi-algorithm only changes to the algorithm are commented.
 *
 * NOTE: grounded areas are not computed with this method.
 */
void red_black(int X, int Y, double NP, double *phi, int NPHI,
		pot *phi0, char *dif, char *logf){
	double *dum, pres=1.+NP, buf=0.;
	int x, y, i,ix, eo;
	double t1, t2;
	FILE *file;
	t1 = MPI_Wtime();
	dum=malloc(X*Y*sizeof(double));
	for(ix=0; ix<X*Y; ix++)
		dum[ix]=phi[ix];

	check_file(dif, "w");
	file=fopen(dif, "w");
	i=0;
	while((i<(int)NP && NP>1.) || (pres>NP && NP<1.)){
		/*loop over even and odd lattice sites, which are ordered like:
		 * eoeoeo
		 * oeoeoe
		 * eoeoeo
		 * oeoeoe
		 * */
		for(eo=0; eo<2; eo++){
			for(x=1; x<X-1; x++){
				for(y=1; y<Y-1; y++){
					if((x+y)%2==eo && check_potential(x, y, NPHI, phi0)==0){/*compute either even or odd lattice sites*/
						dum[y*X+x]=(phi[(y-1)*X+x]+phi[(y+1)*X+x]+phi[y*X+(x-1)]+phi[y*X+(x+1)])/4.;
					}
				}
			}
			/*partial update of either even or odd lattice sites.*/
			for(x=1, buf=0.; x<X-1; x++){
				for(y=1; y<Y-1; y++){
					if((x+y)%2==eo){
						ix=y*X+x;
						if(fabs(phi[ix]-dum[ix])>buf)
							buf=fabs(phi[ix]-dum[ix]);

						phi[ix]=dum[ix];
					}
				}
			}
		}
		fprintf(file, "%d %.6e\n", i, buf);
		pres=buf;
		i++;
	}
	fclose(file);
	check_file(logf, "a");
	file=fopen(logf, "a");
	fprintf(file, "###################################################\n");
	fprintf(file, "#X = %d\t Y = %d\t red-black\n", X, Y);
	fprintf(file, "#\tX*Y\tN\tP\tR\tTime\n");
	t2 = MPI_Wtime();
	fprintf(file, "%8d %8d %.2e %.4e\n", X*Y, i-1, pres, t2-t1);
	fclose(file);
	free(dum);
}

/**
 * Implementation of the gauss-seidel algorithm for the solution of the 2d-laplace equation.
 * Since it is largely the same as the jacobi-algorithm only changes to the algorithm are commented.
 *
 * NOTE: grounded areas are not computed with this method.
 */
void gauss_seidel(int X, int Y, double NP, double *phi, int NPHI,
		pot *phi0, char *dif, char *logf){
	double  pres=1.+NP, buf=0., buf2=0.;
	int x, y, i, ix;
	double t1, t2;
	FILE *file;
	t1 = MPI_Wtime();
	check_file(dif, "w");
	file=fopen(dif, "w");
	i=0;
	while((i<(int)NP && NP>1.) || (pres>NP && NP<1.)){
		for(y=1, buf2=0.; y<Y-1; y++){
			for(x=1; x<X-1; x++){
				ix=y*X+x;
				if(check_potential(x, y, NPHI, phi0)==0){
					/*the buffer is only needed to compute the current precision*/
					buf=(phi[(y-1)*X+x]+phi[(y+1)*X+x]+phi[y*X+(x-1)]+phi[y*X+(x+1)])/4.;
					if(fabs(phi[ix]-buf)>buf2){
						buf2=fabs(phi[ix]-buf);
					}
					phi[ix]=buf;/*update of the new step within the same loop as the actual computation.*/
				}
			}
		}
		pres=buf2;
		fprintf(file, "%d %.6e\n", i, pres);
		i++;
	}
	fclose(file);
	check_file(logf, "a");
	file=fopen(logf, "a");
	fprintf(file, "###################################################\n");
	fprintf(file, "#X = %d\t Y = %d\t gauss-seidel\n", X, Y);
	fprintf(file, "#\tX*Y\tN\tP\tR\tTime\n");
	t2 = MPI_Wtime();
	fprintf(file, "%8d %8d %.2e %.4e\n", X*Y, i-1, pres, t2-t1);
	fclose(file);
}


#endif