#include <fstream>
#include <iostream>
#include <stdlib.h>
#include <string.h>
#include "petscksp.h"
using namespace std;

#define PI 3.14159265358979323846264338327950288419716939937510

static char help[] = "[A][x]=[b] posprocesado.\n\n";

int archivoGmshPotEle();
int archivoGmshBxyEle();
int enrutar(const char * );

Vec Pv_solution;
int N_N, N_E;
double **pNode;
const char *ruta;
char *ptrRutaArchivo;

int main(int argc, char** args)
{
	PetscViewer fd;
	char file[PETSC_MAX_PATH_LEN], direccion[PETSC_MAX_PATH_LEN];
	PetscTruth flg;
	ifstream fileSumary, fileNodes;
	double ignore;

	PetscInitialize(&argc,&args,(char *)0,help);
	PetscOptionsGetString(PETSC_NULL,"-r",direccion,PETSC_MAX_PATH_LEN-1,&flg); //Ruta
	ruta = direccion;
	PetscOptionsGetString(PETSC_NULL,"-f",file,PETSC_MAX_PATH_LEN-1,&flg); //Vector x
	PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
	VecLoad(fd,PETSC_NULL,&Pv_solution);
	PetscViewerDestroy(fd);

	int size=strlen(ruta);
	ptrRutaArchivo=new char[size+99];
	enrutar("/temporal/g2f.ene");

	fileSumary.open(ptrRutaArchivo, ios::in);	// input summary file
		if(!fileSumary.is_open())
			{ cout << "P_sistemaEc: Problema al leer archivo temporal "<<ptrRutaArchivo<<endl;exit(0); }
	fileSumary >> N_E; //elementos
	fileSumary >> N_N; // nodos

	enrutar("/temporal/g2f.node");
	fileNodes.open(ptrRutaArchivo, ios::in);	// input nodes file
		if(!fileNodes.is_open())
			{ cout << "femsii_pos: Problema al leer archivo temporal g2f.node\n";exit(0); }
	pNode = new double*[N_N]; for(int i=0;i<N_N;i++) pNode[i] = new double[2];
	for( int i=0;i<=N_N;i++ )
		{ fileNodes >> ignore >> pNode[i][0] >> pNode[i][1]; }
		// caution! first pair of node values are [1][0] and [1][1]

	fileNodes.close();
	fileSumary.close();

	if(archivoGmshPotEle()) cout << "Archivos de potencial creados.\n";
	if(archivoGmshBxyEle()) cout << "Archivos de induccion creados.\n";

	VecDestroy(Pv_solution);
	PetscFinalize();
	return 0;
}

/*! ENRUTAR() coloca en ptrRutaArchivo la suma de *ruta y *archivo.
 */
int enrutar(const char *archivo )
{
	strcpy(ptrRutaArchivo,ruta);
	strncat(ptrRutaArchivo,archivo,strlen(archivo));
	cout << "Archivo enrrutado: " << ptrRutaArchivo << endl;
	return 1;
}

/*! ARCHIVO-GMSH-POT-ELE() crea en la carpeta "ruta" un archivo de resultados GMSH con
 * los valores de potencial magnetico resultantes por elemento. Utilizando el formato VT.
 * Es necesario utilizar el formato VT para que GSMH pueda hacer el rotacional.
 */
int archivoGmshPotEle()
{
	if( ruta == NULL ) { cout << "P_sistemaEc: No hay ruta para el archivo GmshPot.txt\n"; exit(0); }
	int size=strlen(ruta);
	char *ptrRuta1=new char[size+15];
	char *ptrRuta2=new char[size+15];

	int n1, n2, n3, ignorar;

	PetscScalar Ps_valor1,Ps_valor2,Ps_valor3;
	PetscInt *posicion1,*posicion2,*posicion3 ;
	posicion1 = new PetscInt[1];
	posicion2 = new PetscInt[1];
	posicion3 = new PetscInt[1];

	ifstream fileElements;
	ofstream fileGmshRe,fileGmshIm;

	enrutar("/temporal/g2f.ele");
	fileElements.open(ptrRutaArchivo, ios::in);	// input elements file
		if(!fileElements.is_open()){cout<<"P_sistemaEc: Error leyendo archivo g2f.ele\n";exit(0);}

	// Parte Real a mshPotReal.ele
	strcpy(ptrRuta1,ruta); strncat(ptrRuta1,"/mshPotReal.veg",15);
	strcpy(ptrRuta2,ruta); strncat(ptrRuta2,"/mshPotImag.veg",15);
	fileGmshRe.open(ptrRuta1, ios::out);
	fileGmshIm.open(ptrRuta2, ios::out);
	if(!fileGmshRe.is_open()){cout<<"P_sistemaEc: Error creando archivo de salida Gmsh\n";exit(0);}
	if(!fileGmshIm.is_open()){cout<<"P_sistemaEc: Error creando archivo de salida Gmsh\n";exit(0);}
	fileGmshRe << "View \"Pot Real\" {" << endl;
	fileGmshIm << "View \"Pot Imaginario\" {" << endl;

	for(int i=0;i<N_E;i++){
		fileElements >> ignorar >> n1 >> n2 >> n3 >> ignorar;

		posicion1[0] = n1;
		posicion2[0] = n2;
		posicion3[0] = n3;

		VecGetValues(Pv_solution,1,posicion1,&Ps_valor1);
		VecGetValues(Pv_solution,1,posicion2,&Ps_valor2);
		VecGetValues(Pv_solution,1,posicion3,&Ps_valor3);

		fileGmshRe << "VT(" << pNode[n1][0] << "," << pNode[n1][1] << ",0,";
		fileGmshRe << pNode[n2][0] << "," << pNode[n2][1] << ",0,";
		fileGmshRe << pNode[n3][0] << "," << pNode[n3][1] << ",0)";
		fileGmshRe << "{0,0," << PetscRealPart(Ps_valor1);
		fileGmshRe << ",0,0," << PetscRealPart(Ps_valor2);
		fileGmshRe << ",0,0," << PetscRealPart(Ps_valor3) << "};\n";

		fileGmshIm << "VT(" << pNode[n1][0] << "," << pNode[n1][1] << ",0,";
		fileGmshIm << pNode[n2][0] << "," << pNode[n2][1] << ",0,";
		fileGmshIm << pNode[n3][0] << "," << pNode[n3][1] << ",0)";
		fileGmshIm << "{0,0," << PetscImaginaryPart(Ps_valor1);
		fileGmshIm << ",0,0," << PetscImaginaryPart(Ps_valor2);
		fileGmshIm << ",0,0," << PetscImaginaryPart(Ps_valor3) << "};\n";

	}
	fileGmshRe << "};";
	fileGmshRe.close();

	fileGmshIm << "};";
	fileGmshIm.close();

	fileElements.close();
	cout << "Se han creado los archivos de potencial magnetico.\n";
	return 1;
}

/*! ARCHIVO-GMSH-Bxy-ELE() crea en la carpeta "ruta" un archivo de resultados GMSH con
 * los valores de induccion magnetica resultantes por elemento. Utilizando el formato VT.
 * Es necesario utilizar el formato VT para que GSMH pueda realizar operaciones vectoriales.
 */
int archivoGmshBxyEle()
{
	if( ruta == NULL ) { cout << "P_sistemaEc: No hay ruta para el archivo GmshPot.txt\n"; exit(0); }
	int size=strlen(ruta);
	char *ptrRuta1=new char[size+15];
	char *ptrRuta2=new char[size+15];
	int n1, n2, n3, ignorar;
	double bi,bj,bk,ci,cj,ck,delta,xi,yi,xj,yj,xk,yk, Bx, By, centroX, centroY;
	ifstream fileElements;
	ofstream fileGmshRe,fileGmshIm;

	PetscScalar Ps_valor1,Ps_valor2,Ps_valor3;
	PetscInt *posicion1,*posicion2,*posicion3 ;
	posicion1 = new PetscInt[1];
	posicion2 = new PetscInt[1];
	posicion3 = new PetscInt[1];

	enrutar("/temporal/g2f.ele");
	fileElements.open(ptrRutaArchivo, ios::in);	// input elements file
		if(!fileElements.is_open()){cout<<"P_sistemaEc: Error leyendo archivo g2f.ele\n";exit(0);}

	// Parte Real a mshBxyReal.ele
	strcpy(ptrRuta1,ruta); strncat(ptrRuta1,"/mshBxyReal.veg",15);
	strcpy(ptrRuta2,ruta); strncat(ptrRuta2,"/mshBxyImag.veg",15);
	fileGmshRe.open(ptrRuta1, ios::out);
	fileGmshIm.open(ptrRuta2, ios::out);
	if(!fileGmshRe.is_open()){cout<<"P_sistemaEc: Error creando archivo de salida Gmsh\n";exit(0);}
	if(!fileGmshIm.is_open()){cout<<"P_sistemaEc: Error creando archivo de salida Gmsh\n";exit(0);}
	fileGmshRe << "View \"Bxy Real\" {" << endl;
	fileGmshIm << "View \"Bxy Imaginario\" {" << endl;

	for(int i=0;i<N_E;i++){
		fileElements >> ignorar >> n1 >> n2 >> n3 >> ignorar;

		xi=pNode[n1][0]; yi=pNode[n1][1]; // Nodo 1 coordenada x e y
		xj=pNode[n2][0]; yj=pNode[n2][1]; // Nodo 2 coordenada x e y
		xk=pNode[n3][0]; yk=pNode[n3][1]; // Nodo 3 coordenada x e y
		bi=yj-yk; bj=yk-yi; bk=yi-yj; ci=xk-xj; cj=xi-xk; ck=xj-xi;
		delta=0.5*(xi*yj+yi*xk+xj*yk-xk*yj-xj*yi-yk*xi); // NO SE UTILIZA VALOR ABSOLUTO PARA SER COHERENTE CON EL SIGNO DE LA DIRECCION DE MALLADO

		posicion1[0] = n1;
		posicion2[0] = n2;
		posicion3[0] = n3;

		VecGetValues(Pv_solution,1,posicion1,&Ps_valor1);
		VecGetValues(Pv_solution,1,posicion2,&Ps_valor2);
		VecGetValues(Pv_solution,1,posicion3,&Ps_valor3);

		centroX = (xi+xj+xk)/3;
		centroY = (yi+yj+yk)/3;

		// Bx = (dA/dy) ; By = -(dA/dx) // Parte real
		Bx=(ci*PetscRealPart(Ps_valor1)+cj*PetscRealPart(Ps_valor2)+ck*PetscRealPart(Ps_valor3))/(2*delta);
		By=-1*(bi*PetscRealPart(Ps_valor1)+bj*PetscRealPart(Ps_valor2)+bk*PetscRealPart(Ps_valor3))/(2*delta);

		fileGmshRe << "VP(" << centroX << "," << centroY << ",0)";
		fileGmshRe << "{" << Bx << "," << By << ",0};\n";

		// Bx = -(dA/dy) ; By = -(dA/dx) // Parte imaginaria
		Bx=(ci*PetscImaginaryPart(Ps_valor1)+cj*PetscImaginaryPart(Ps_valor2)+ck*PetscImaginaryPart(Ps_valor3))/(2*delta);
		By=-1*(bi*PetscImaginaryPart(Ps_valor1)+bj*PetscImaginaryPart(Ps_valor2)+bk*PetscImaginaryPart(Ps_valor3))/(2*delta);


		fileGmshIm << "VP(" << centroX << "," << centroY << ",0)";
		fileGmshIm << "{" << Bx << "," << By << ",0};\n";

	}

	fileGmshRe << "};";
	fileGmshRe.close();

	fileGmshIm << "};";
	fileGmshIm.close();

	fileElements.close();
	cout << "Se han creado los archivos de induccion magnetica.\n";
	return 1;
}

/*!
\file femsii_pos.cpp
\brief
Este programa acepta (obligatorios) dos parametros de entrada:
-f Archivo binario que contiene el vector solucion generado por solve.bin.
-r Ruta donde seran escritos los archivos generados por pos.bin.

DESCRIPCION DE main(int argc, char **args): Asume los parametros de entrada. Busca los archivos temporados para rellenar las variables N_N (numero de nodos), N_E (numero de elementos) y pNode (vector doble para pares de nodos). Luego lanza las funciones archivoGmshPotEle() y archivoGmshBxyEle().

DESCRIPCION DE archivoGmshPotEle(): Crea los archivos de potencial vector magnetico para Gmsh: mshPotReal.veg y mshPotImag.veg

DESCRIPCION DE archivoGmshBxyEle(): Realiza la operacion rotacional de forma numerica (para cada elemento) y crea los archivos de flujo magnetico para Gmsh: mshBxyReal.veg y mshBxyImag.veg
*/