#include <iostream>
#include <math.h>
#include "pmmc.h"
#include "Domain.h"
//#include "PointList.h"
//#include "Array.h"

using namespace std;

int main(int argc, char **argv)
{
	//.......................................................................
	//  printf("Radius = %s \n,"RADIUS);  
	int Nx,Ny,Nz,N;
	int i,j,k,p,q,r,n;
	int ReadSize;  // Number of entries in the input files that are read
	int nspheres;
	double Lx,Ly,Lz;
	//.......................................................................

	//.......................................................................
	string PhiName,PressName,VelName;
	ifstream readnames("Filenames.in");
	readnames >> PhiName;
	readnames >> PressName;
	readnames >> VelName;
	readnames >> ReadSize;
	//.......................................................................
	// Reading the domain information file
	//.......................................................................
	ifstream domain("Domain.in");
	domain >> Nx;
	domain >> Ny;
	domain >> Nz;
	domain >> nspheres;
	domain >> Lx;
	domain >> Ly;
	domain >> Lz;
	//.......................................................................
	N = Nx*Ny*Nz;
	printf("Domain size = %i x %i x %i \n",Nx,Ny,Nz);
	printf("Domain length = %f x %f x %f \n",Lx,Ly,Lz);
	 
	cout << "PhiName = " << PressName << endl;
	cout << "ReadSize = " << ReadSize << endl;
	
	//.......................................................................
	DoubleArray SignDist(Nx,Ny,Nz);
	DoubleArray SignDistCopy(Nx+2,Ny+2,Nz+2);
	DoubleArray Phase(Nx,Ny,Nz);
	DoubleArray Phase_x(Nx,Ny,Nz);
	DoubleArray Phase_y(Nx,Ny,Nz);
	DoubleArray Phase_z(Nx,Ny,Nz);
	DoubleArray Sx(Nx,Ny,Nz);
	DoubleArray Sy(Nx,Ny,Nz);
	DoubleArray Sz(Nx,Ny,Nz);
	DoubleArray Vel_x(Nx,Ny,Nz);
	DoubleArray Vel_y(Nx,Ny,Nz);
	DoubleArray Vel_z(Nx,Ny,Nz);
	DoubleArray Press(Nx,Ny,Nz);
	DoubleArray GaussCurvature(Nx,Ny,Nz);
	DoubleArray MeanCurvature(Nx,Ny,Nz);
	//.......................................................................
	
	// Read in sphere pack
	printf("nspheres =%i \n",nspheres);
	//.......................................................................
	double *cx,*cy,*cz,*rad;
	cx = new double[nspheres];
	cy = new double[nspheres];
	cz = new double[nspheres];
	rad = new double[nspheres];
	//.......................................................................
	int *ReadInfo;
	double *ReadPhi, *ReadPress, *ReadVel;
	ReadInfo = new int[3*ReadSize];
	ReadPhi = new double[ReadSize];
	ReadPress = new double[ReadSize];
	ReadVel = new double[3*ReadSize];
	//.......................................................................
	printf("Reading the sphere packing \n");
	ReadSpherePacking(nspheres,cx,cy,cz,rad);

	SignedDistance(SignDistCopy.data,nspheres,cx,cy,cz,rad,Lx,Ly,Lz,Nx+2,Ny+2,Nz+2,
				   0,0,0,1,1,1);
	
	for (k=0; k<Nz; k++){
		for (j=0; j<Ny; j++){
			for (i=0; i<Nx; i++){
				SignDist(i,j,k) =  SignDistCopy(i+1,j+1,k+1)-1.0;
			}
		}
	}
	
	//.......................................................................
	// Read all of the output files from Lattice Boltzmann Simulation
	//.......................................................................
	printf("Reading data from the Lattice Boltzmann output files \n");
	int infovalue;
	double value;
	ifstream INFO("Info",ios::binary);
	for (n=0; n< ReadSize; n++){
		// Read 3D indices from the info file 
		INFO.read((char*) &infovalue, sizeof(infovalue));
		ReadInfo[3*n] = infovalue;
		INFO.read((char*) &infovalue, sizeof(infovalue));
		ReadInfo[3*n+1] = infovalue;
		INFO.read((char*) &infovalue, sizeof(infovalue));
		ReadInfo[3*n+2] = infovalue-4; // account for pressure boundary conditions 
	}
	INFO.close();
	ifstream PHI(PhiName.c_str(),ios::binary);
	for (n=0; n< ReadSize; n++){
		PHI.read((char*) &value, sizeof(value));
		ReadPhi[n] = value;
	}
	PHI.close();
	ifstream PRESSURE(PressName.c_str(),ios::binary);
	for (n=0; n< ReadSize; n++){
		PRESSURE.read((char*) &value, sizeof(value));
		ReadPress[n] = value;
	}
	PRESSURE.close();
	
#ifdef ReadVelocityFromFile
	ifstream VELOCITY(VelName.c_str(),ios::binary);
	for (n=0; n< ReadSize; n++){
		VELOCITY.read((char*) &value, sizeof(value));
		ReadVel[3*n] = value;
		VELOCITY.read((char*) &value, sizeof(value));
		ReadVel[3*n+1] = value;
		VELOCITY.read((char*) &value, sizeof(value));
		ReadVel[3*n+2] = value;
	}
	VELOCITY.close();
#else
	value = 0.f;
	for (n=0; n< ReadSize; n++){
		ReadVel[3*n] = value;
		ReadVel[3*n+1] = value;
		ReadVel[3*n+2] = value;
	}
#endif
	//.......................................................................
	
	//.......................................................................
	// Reconstruct the global meshes using the info file
	//.......................................................................
	for (int idx=0; idx< ReadSize; idx++){
		// 3D indices
		i = ReadInfo[3*idx];
		j = ReadInfo[3*idx+1];
		k = ReadInfo[3*idx+2];
	//	printf("%i, %i, %i \n",i,j,k);
		// Update the values 
		if (k>0){
			Phase(i,j,k) = ReadPhi[idx];
			Press(i,j,k) = ReadPress[idx];
			Vel_x(i,j,k) = ReadVel[3*idx];
			Vel_y(i,j,k) = ReadVel[3*idx+1];
			Vel_z(i,j,k) = ReadVel[3*idx+2];
		}
	}	
	//.......................................................................
	
	//.......................................................................
	double fluid_isovalue = 0.0;
	double solid_isovalue = 0.0;
	//.......................................................................

/*	//.......................................................................
	double *cx,*cy,*cz,*rad;
	cx = new double[nspheres];
	cy = new double[nspheres];
	cz = new double[nspheres];
	rad = new double[nspheres];
	//...............................
	printf("Reading the sphere packing \n");
	ReadSpherePacking(nspheres,cx,cy,cz,rad);
	//.......................................................................
	//.......................................................................
	// Compute the signed distance function for the sphere packing
	SignedDistance(SignDist.data,nspheres,cx,cy,cz,rad,Lx,Ly,Lz,Nx,Ny,Nz,0,0,0,1,1,1);
*/	//.......................................................................

	/* ****************************************************************
	 VARIABLES FOR THE PMMC ALGORITHM
	 ****************************************************************** */
	//...........................................................................
	// Averaging variables
	//...........................................................................
	double awn,ans,aws,lwns,nwp_volume;
	double sw,vol_n,vol_w,paw,pan;
	double efawns,Jwn;
	double As;
	double dEs,dAwn,dAns;			 // Global surface energy (calculated by rank=0)
	double awn_global,ans_global,aws_global,lwns_global,nwp_volume_global;	
	double As_global;
	//	bool add=1;			// Set to false if any corners contain nw-phase ( F > fluid_isovalue)
	
	int n_nw_pts=0,n_ns_pts=0,n_ws_pts=0,n_nws_pts=0, map=0;
	int n_nw_tris=0, n_ns_tris=0, n_ws_tris=0, n_nws_seg=0;
	
	double s,s1,s2,s3;		// Triangle sides (lengths)
	Point A,B,C,P;
	//	double area;
	int cube[8][3] = {{0,0,0},{1,0,0},{0,1,0},{1,1,0},{0,0,1},{1,0,1},{0,1,1},{1,1,1}};  // cube corners
	//	int count_in=0,count_out=0;
	//	int nodx,nody,nodz;
	// initialize lists for vertices for surfaces, common line
	DTMutableList<Point> nw_pts(20);
	DTMutableList<Point> ns_pts(20);
	DTMutableList<Point> ws_pts(20);
	DTMutableList<Point> nws_pts(20);
	// initialize triangle lists for surfaces
	IntArray nw_tris(3,20);
	IntArray ns_tris(3,20);
	IntArray ws_tris(3,20);
	// initialize list for line segments
	IntArray nws_seg(2,20);
	DTMutableList<Point> tmp(20);
	
	// Initialize arrays for local solid surface
	DTMutableList<Point> local_sol_pts(20);
	int n_local_sol_pts = 0;
	IntArray local_sol_tris(3,18);
	int n_local_sol_tris;
	DoubleArray values(20);
	DTMutableList<Point> local_nws_pts(20);
	int n_local_nws_pts;
	
	DoubleArray CubeValues(2,2,2);
	DoubleArray ContactAngle(20);
	DoubleArray Curvature(20);
	DoubleArray InterfaceSpeed(20);
	DoubleArray NormalVector(60);
	DoubleArray van(3);
	DoubleArray vaw(3);
	DoubleArray vawn(3);
	DoubleArray Gwn(6);
	DoubleArray Gns(6);
	DoubleArray Gws(6);
	
	double iVol = 1.0/Nx/Ny/Nz;
	
	int c;
	//...........................................................................
	int ncubes = (Nx-2)*(Ny-2)*(Nz-2);	// Exclude the "upper" halo
	IntArray cubeList(3,ncubes);
	pmmc_CubeListFromMesh(cubeList, ncubes, Nx, Ny, Nz);
	//...........................................................................
	double Cx,Cy,Cz;
	double dist1,dist2;
	// Extra copies of phase indicator needed to compute time derivatives on CPU
	DoubleArray Phase_tminus(Nx,Ny,Nz);
	DoubleArray Phase_tplus(Nx,Ny,Nz);
	DoubleArray dPdt(Nx,Ny,Nz);
	for (k=0; k<Nz; k++){
		for (j=0; j<Ny; j++){
			for (i=0; i<Nx; i++){
				Phase_tminus(i,j,k) = Phase(i,j,k);
				Phase_tplus(i,j,k) = Phase(i,j,k);
				dPdt(i,j,k) = 0.0;
			}
		}
	}
	
	printf("Computing TCAT Averages \n");
	printf("--------------------------------------------------------------------------------------\n");
	printf("sw pw pn vw[x, y, z] vn[x, y, z] ");			// Volume averages
	printf("awn ans aws Jwn lwns efawns ");	// Interface and common curve averages
	printf("Gwn [xx, yy, zz, xy, xz, yz] ");				// Orientation tensors
	printf("Gws [xx, yy, zz, xy, xz, yz] ");
	printf("Gns [xx, yy, zz, xy, xz, yz] \n");
	printf("--------------------------------------------------------------------------------------\n");
	
	awn = aws = ans = lwns = 0.0;
	nwp_volume = 0.0;
	As = 0.0;
	Jwn = 0.0;
	efawns = 0.0;
	// Compute phase averages
	pan = paw = 0.0;
	vaw(0) = vaw(1) = vaw(2) = 0.0;
	van(0) = van(1) = van(2) = 0.0;
	vawn(0) = vawn(1) = vawn(2) = 0.0;
	Gwn(0) = Gwn(1) = Gwn(2) = 0.0;
	Gwn(3) = Gwn(4) = Gwn(5) = 0.0;
	Gws(0) = Gws(1) = Gws(2) = 0.0;
	Gws(3) = Gws(4) = Gws(5) = 0.0;
	Gns(0) = Gns(1) = Gns(2) = 0.0;
	Gns(3) = Gns(4) = Gns(5) = 0.0;
	vol_w = vol_n =0.0;
	
	// Calculate the time derivative of the phase indicator field
	for (int n=0; n<Nx*Ny*Nz; n++)	dPdt(n) = 0.5*(Phase_tplus(n) - Phase_tminus(n));
	
	pmmc_MeshGradient(Phase,Phase_x,Phase_y,Phase_z,Nx,Ny,Nz);
	pmmc_MeshGradient(SignDist,Sx,Sy,Sz,Nx,Ny,Nz);
	pmmc_MeshCurvature(Phase, MeanCurvature, GaussCurvature, Nx, Ny, Nz);
	
	/// Compute volume averages
	int total = 0;
	// Compute phase averages
	nwp_volume = 0.0;
	pan = paw = 0.0;
	vaw(0) = vaw(1) = vaw(2) = 0.0;
	vol_w = vol_n =0.0;
	for (k=0; k<Nz; k++){
		for (j=0; j<Ny; j++){
			for (i=0; i<Nx; i++){
				if ( SignDist(i,j,k) > 0 ){
					total++;
					// 1-D index for this cube corner
					n = i + j*Nx + k*Nx*Ny;

					// Compute the non-wetting phase volume contribution
					if ( Phase(i,j,k) > 0.0 )
						nwp_volume += 1.0;

					// volume averages over the non-wetting phase
					if ( Phase(i,j,k) > 0.9999 ){
						// volume the excludes the interfacial region
						vol_n += 1.0;
						// pressure
						pan += Press(i,j,k);
						// velocity
						van(0) += Vel_x(i,j,k);
						van(1) += Vel_y(i,j,k);
						van(2) += Vel_z(i,j,k);
					}

					// volume averages over the wetting phase
					if ( Phase(i,j,k) < -0.9999 ){
						// volume the excludes the interfacial region
						vol_w += 1.0;
						// pressure
						paw += Press(i,j,k);
						// velocity
						vaw(0) += Vel_x(i,j,k);
						vaw(1) += Vel_y(i,j,k);
						vaw(2) += Vel_z(i,j,k);
						
					}
				}
			}
		}
	}
	
	// End of the loop to set the values
	awn = aws = ans = lwns = 0.0;
	As = 0.0;
	Gwn(0) = Gwn(1) = Gwn(2) = Gwn(3) = Gwn(4) = Gwn(5) = 0.0;
	Gws(0) = Gws(1) = Gws(2) = Gws(3) = Gws(4) = Gws(5) = 0.0;
	Gns(0) = Gns(1) = Gns(2) = Gns(3) = Gns(4) = Gns(5) = 0.0;
	
	FILE *WN_TRIS;
	WN_TRIS = fopen("wn-tris.out","w");
	
	FILE *NS_TRIS;
	NS_TRIS = fopen("ns-tris.out","w");
	
	FILE *WS_TRIS;
	WS_TRIS = fopen("ws-tris.out","w");
	
	FILE *WNS_PTS;
	WNS_PTS = fopen("wns-pts.out","w");
		
	for (c=0;c<ncubes;c++){
		// Get cube from the list
		i = cubeList(0,c);
		j = cubeList(1,c);
		k = cubeList(2,c);
		
		// Run PMMC
		n_local_sol_tris = 0;
		n_local_sol_pts = 0;
		n_local_nws_pts = 0;
		n_nw_pts=0,n_ns_pts=0,n_ws_pts=0,n_nws_pts=0, map=0;
		n_nw_tris=0, n_ns_tris=0, n_ws_tris=0, n_nws_seg=0;

		// Construct the interfaces and common curve
		pmmc_ConstructLocalCube(SignDist, Phase, solid_isovalue, fluid_isovalue,
				nw_pts, nw_tris, values, ns_pts, ns_tris, ws_pts, ws_tris,
				local_nws_pts, nws_pts, nws_seg, local_sol_pts, local_sol_tris,
				n_local_sol_tris, n_local_sol_pts, n_nw_pts, n_nw_tris,
				n_ws_pts, n_ws_tris, n_ns_tris, n_ns_pts, n_local_nws_pts, n_nws_pts, n_nws_seg,
				i, j, k, Nx, Ny, Nz);
		
		// Compute the velocity of the wn interface
		pmmc_InterfaceSpeed(dPdt, Phase_x, Phase_y, Phase_z, CubeValues, nw_pts, nw_tris,
							NormalVector, InterfaceSpeed, vawn, i, j, k, n_nw_pts, n_nw_tris);
		
//		pmmc_InterfaceSpeed(dPdt, Phase_x, Phase_y, Phase_z, CubeValues, nw_pts, nw_tris,
//							NormalVector, InterfaceSpeed, vawn, i, j, k, n_nw_pts, n_nw_tris);
		
		
		// Compute the average contact angle
		efawns += pmmc_CubeContactAngle(CubeValues,ContactAngle,Phase_x,Phase_y,Phase_z,Sx,Sy,Sz,
										local_nws_pts,i,j,k,n_local_nws_pts);

		// Compute the curvature of the wn interface
		Jwn += pmmc_CubeSurfaceInterpValue(CubeValues, MeanCurvature, nw_pts, nw_tris,
											Curvature, i, j, k, n_nw_pts, n_nw_tris);
		
		// Compute the surface orientation and the interfacial area
		awn += pmmc_CubeSurfaceOrientation(Gwn,nw_pts,nw_tris,n_nw_tris);
		ans += pmmc_CubeSurfaceOrientation(Gns,ns_pts,ns_tris,n_ns_tris);
		aws += pmmc_CubeSurfaceOrientation(Gws,ws_pts,ws_tris,n_ws_tris);
		
		//*******************************************************************
		// Compute the Interfacial Areas, Common Line length
		As += pmmc_CubeSurfaceArea(local_sol_pts,local_sol_tris,n_local_sol_tris);
		lwns +=  pmmc_CubeCurveLength(local_nws_pts,n_local_nws_pts);
		
		//.......................................................................................
		// Write the triangle lists to text file
		for (r=0;r<n_nw_tris;r++){
			A = nw_pts(nw_tris(0,r));
			B = nw_pts(nw_tris(1,r));
			C = nw_pts(nw_tris(2,r));
			fprintf(WN_TRIS,"%f %f %f %f %f %f %f %f %f \n",A.x,A.y,A.z,B.x,B.y,B.z,C.x,C.y,C.z);
		}		
		for (r=0;r<n_ws_tris;r++){
			A = ws_pts(ws_tris(0,r));
			B = ws_pts(ws_tris(1,r));
			C = ws_pts(ws_tris(2,r));
			fprintf(WS_TRIS,"%f %f %f %f %f %f %f %f %f \n",A.x,A.y,A.z,B.x,B.y,B.z,C.x,C.y,C.z);
		}
		for (r=0;r<n_ns_tris;r++){
			A = ns_pts(ns_tris(0,r));
			B = ns_pts(ns_tris(1,r));
			C = ns_pts(ns_tris(2,r));
			fprintf(NS_TRIS,"%f %f %f %f %f %f %f %f %f \n",A.x,A.y,A.z,B.x,B.y,B.z,C.x,C.y,C.z);
		}
		for (p=0; p < n_nws_pts; p++){
			P = nws_pts(p);
			fprintf(WNS_PTS,"%f %f %f \n",P.x, P.y, P.z);
		}
		
	}
	fclose(WN_TRIS);
	fclose(NS_TRIS);
	fclose(WS_TRIS);
	fclose(WNS_PTS);

	printf("paw = %f \n",paw);
	printf("vol_w = %f \n",vol_w);
	printf("efawns = %f \n",efawns);
	printf("lwns = %f \n",lwns);
	printf("efawns = %f \n",efawns/lwns);

	paw /= vol_w;
	pan /= vol_n;
	for (i=0; i<3; i++)	vaw(i) /= vol_w;
	for (i=0; i<3; i++)	van(i) /= vol_n;
	
	Jwn /= awn;
	efawns /= lwns;
	for (i=0; i<3; i++)	vawn(i) /= awn;
	for (i=0; i<6; i++)	Gwn(i) /= awn;
	for (i=0; i<6; i++)	Gns(i) /= ans;
	for (i=0; i<6; i++)	Gws(i) /= aws;
	
	awn = awn*iVol;
	aws = aws*iVol;
	ans = ans*iVol;
	lwns = lwns*iVol;
	sw = 1.0 - nwp_volume/total;

	printf("%.5g %.5g %.5g ",sw,paw,pan);					// saturation and pressure
	printf("%.5g %.5g %.5g ",vaw(0),vaw(1),vaw(2));			// average velocity of w phase
	printf("%.5g %.5g %.5g ",van(0),van(1),van(2));			// average velocity of n phase
	printf("%.5g %.5g %.5g ",awn,ans,aws);					// interfacial areas
	printf("%.5g ",Jwn);									// curvature of wn interface
	printf("%.5g ",lwns);									// common curve length
	printf("%.5g ",efawns);									// average contact angle
	printf("%.5g %.5g %.5g %.5g %.5g %.5g ",
			Gwn(0),Gwn(1),Gwn(2),Gwn(3),Gwn(4),Gwn(5));		// orientation of wn interface
	printf("%.5g %.5g %.5g %.5g %.5g %.5g ",
			Gns(0),Gns(1),Gns(2),Gns(3),Gns(4),Gns(5));		// orientation of ns interface	
	printf("%.5g %.5g %.5g %.5g %.5g %.5g \n",
			Gws(0),Gws(1),Gws(2),Gws(3),Gws(4),Gws(5));		// orientation of ws interface
	
	FILE *PHASE;
	PHASE = fopen("Phase.out","wb");
	fwrite(Phase.data,8,N,PHASE);
	fclose(PHASE);
	
	FILE *SOLID;
	SOLID = fopen("Distance.out","wb");
	fwrite(SignDist.data,8,N,SOLID);
	fclose(SOLID);

}