#include "P_sistemaEc.h"
#include <fstream>
#include <iostream>
#include <iomanip>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "petscksp.h"

#define PI 3.14159265358979323846264338327950288419716939937510

using namespace std;

// PETSC global variables
PetscErrorCode ierr;
PetscScalar **Ps_local; // local stiffness matrix
Mat Pm_global; // global stiffness matrix
Vec Pv_rightside; // rigth hand "B" of the ecuation: Ax=B
Vec Pv_solution; // solution
static char help[] = "Bibliotecas basadas en PETSC para elementos finitos";


//! Constructor de Sistema_ecuaciones .
/*! El CONSTRUCTOR crea las cadenas de caracteres (ptrMesh, ptrConducti, etc) con las direcciones
* de los diferentes archivos de trabajo. Ademas se llama a la funcion getOrder() de
* de la clase G2f para crear archivos ordenados con los datos del archivo GMSH.
* Por ultimo se llama a la funcion regiones, que crear variables de tipo map globales con
* valores de permeabilidad y conductividad.
*/
Sistema_ecuaciones::Sistema_ecuaciones(const char *entrada) // CONSTRUCTOR
{
	if( entrada == NULL ){
		cout << "P_sistemaEc: No hay mallado de entrada\n"; exit(0);
	}else{ ruta=entrada; }

// creating folder. Be Careful!!!! check if "mkdir" can work on other operating systems.
 	if(!mkdir(enrrutar("/temporal"))) //create folder for temporal files: g2f and femsii
 		{ cout << "P_sistemaEc: Problema al crear directorio temporal\n";exit(0); }

	// get order to gmsh files by create a instance of G2f class
	G2f mallado( enrrutar("/mesh.msh"), enrrutar("/temporal") );
	mallado.getOrder();

	// introduce information on LTS maps permeabilidad conductividad.
	regiones();

	// start PETSC machine
	PetscInitialize(PETSC_NULL,PETSC_NULL,(char *)0,help);
	
	// get summary from g2f.ene file: count of nodes, elements, edges
	ifstream fileSumary;
	fileSumary.open(enrrutar("/temporal/g2f.ene"), ios::in);
		if(!fileSumary.is_open())
			{ cout <<ruta<< " P_sistemaEc: Problema al leer archivo "<< enrrutar("/temporal/g2f.ene") <<endl;exit(0); }
	fileSumary >> N_E; //elementos
	fileSumary >> N_N; // nodos
	fileSumary >> N_ED; // edges elementos
	fileSumary.close();

	// get air permeability and freq of the system
	ifstream fileExtra;
	fileExtra.open(enrrutar("/extra.txt"), ios::in);
		if(!fileExtra.is_open())
			{ cout << "P_sistemaEc: Problema al leer archivo"<< enrrutar("/extra.txt") <<endl;exit(0); }
	fileExtra >> VACIO;
	fileExtra >> frec;
	fileExtra.close();
}


//! Destructor de Sistema_ecuaciones .
/*! El destructor de Sistema_ecuaciones se encarga de eliminar de la memoria la matriz
* principal de rigidez y el vector de terminos independientes.
* Ademas lanza la orden de finalizar PetscFinalize().
*/
Sistema_ecuaciones::~Sistema_ecuaciones() // DESTRUCTOR
{
	VecDestroy(Pv_rightside);
	MatDestroy(Pm_global);
	PetscFinalize();
}


/*! MKDIR() crea la carpeta que viene en la cadena *ptrTemporal y
 * notifica si lo consigue.
 */
int Sistema_ecuaciones::mkdir( const char *ptrTemporal)
{
	char *ptrComando;
	int answer;
	int size=strlen(ptrTemporal);
	ptrComando=new char[size+6];
	strcpy(ptrComando,"mkdir ");
	strncat(ptrComando,ptrTemporal,size);
	answer=system(ptrComando);
	if(answer==0){
		cout << "P_sistemaEc: Carpeta \"temporal\" creada.\n";
		return 1;}
	if(answer==256){
		cout << "P_sistemaEc: La carpeta \"temporal\" existe.\n";
		return 1;}
	if(answer!=0 && answer!=256){
		cout << "P_sistemaEc: Error al crear la carpera \"temporal\"\n.";
		return 0;}
}


/*! ENRUTAR() es una funcion auxiliar que coloca en ptrRutaArchivo la suma
* de las cadenas de caracteres ruta y *archivo.
*/
const char* Sistema_ecuaciones::enrrutar(const char *archivo )
{
	string cadena=ruta+archivo;
	char *ptrRutaArchivo = new char [cadena.size()+1];
	strcpy (ptrRutaArchivo, cadena.c_str());
	cout << "Archivo enrrutado: " << ptrRutaArchivo << endl;
	return ptrRutaArchivo;
}


/*! La funcion GET-STIFFNESS() utiliza los archivos GMSH ordenados por getOrder() para
 * calcular la matriz global del sistema, tanto parte real como parte imaginaria. Para ello
 * utiliza un buble que llama a la funcion local_matrix() que genera la matriz local de
 * cada elemento, que luego es ensamblada en la matriz global.
 * La funcion local_matrix() necesitara las coordenadas de los nodos que forman el
 * elemento, asi como la region a la que pertenese el elemento.
 */
int Sistema_ecuaciones::getStiffness( int verb )
{

	int element, n1, n2, n3, region, // element, nodes numbers, region of element
	ignore; // inputs to ignore
	ifstream fileElements, fileNodes;


	// nodes to memory --> pNode
	fileNodes.open(enrrutar("/temporal/g2f.node"), ios::in);	// input nodes file
		if(!fileNodes.is_open())
			{ cout << "P_sistemaEc: Problema al leer archivo temporal g2f.node\n";exit(0); }
	//memory for nodes
	pNode = new double*[N_N+1]; 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]; }
	fileNodes.close();


	 // -- PETSC -- //

	ierr = MatCreateSeqAIJ(PETSC_COMM_WORLD,N_N,N_N,10,0,&Pm_global);CHKERRQ(ierr);

	// memory for local stiffness matrix
	Ps_local=new PetscScalar*[3]; for(int i=0;i<3;i++) Ps_local[i] = new PetscScalar[3];
	if( Ps_local==NULL ) { cout << "P_sistemaEc: Problema asignando memoria"; exit(0); }

	// we use a "for" bucle to work with every element
	fileElements.open(enrrutar("/temporal/g2f.ele"), ios::in);	// input elements file
		if(!fileElements.is_open())
			{ cout<<"P_sistemaEc: Problema al leer archivo temporal g2f.ele\n";exit(0); }
	cout << "P_sistemaEc: Creando la matriz de rigidez.\n";
	for( int i=0;i<N_E;i++ ){
		fileElements >> element >> n1 >> n2 >> n3 >> region;

		/*send node coord to local stiffness matrix maker.
		pNode index and the index from node file are the same
		works like: local_matrix( xi, yi,  xj, yj, xk, yk, 'matriz de destino' ) */
		if ( !(
		local_matrix(pNode[n1][0],pNode[n1][1],pNode[n2][0],pNode[n2][1],pNode[n3][0],pNode[n3][1],region)
		))  { cout << "P_sistemaEc: Problema al crear la matriz local del elemento"<<i; exit(0); }

		// if verb==1 activate verbose mode, all local matrix to cout
		if( verb == 1 ) // BE CAREFUL !!!
		{cout<<endl<<" Matriz del elemento "<< element <<endl;
			cout<<"Permeabilidad: "<<permeabilidad[region]<<endl;
			cout<<Ps_local[0][0]<<' '<<Ps_local[0][1]<<' '<<Ps_local[0][2]<<endl;
			cout<<Ps_local[1][0]<<' '<<Ps_local[1][1]<<' '<<Ps_local[1][2]<<endl;
			cout<<Ps_local[2][0]<<' '<<Ps_local[2][1]<<' '<<Ps_local[2][2]<<endl <<endl;
		}

		/* to ensamble PETSC matrix starts in [0] but first gmsh node starts in "1"
		so we will cut "1" from node index */
//		n1=n1-1; n2=n2-1; n3=n3-1;
		// ENSAMBLE NOW to global stiffnesss matrix
		ierr=MatSetValue(Pm_global,n1,n1,Ps_local[0][0],ADD_VALUES);
		CHKERRQ(ierr);
		ierr=MatSetValue(Pm_global,n1,n2,Ps_local[0][1],ADD_VALUES);
		CHKERRQ(ierr);
		ierr=MatSetValue(Pm_global,n1,n3,Ps_local[0][2],ADD_VALUES);
		CHKERRQ(ierr);
		ierr=MatSetValue(Pm_global,n2,n1,Ps_local[1][0],ADD_VALUES);
		CHKERRQ(ierr);
		ierr=MatSetValue(Pm_global,n2,n2,Ps_local[1][1],ADD_VALUES);
		CHKERRQ(ierr);
		ierr=MatSetValue(Pm_global,n2,n3,Ps_local[1][2],ADD_VALUES);
		CHKERRQ(ierr);
		ierr=MatSetValue(Pm_global,n3,n1,Ps_local[2][0],ADD_VALUES);
		CHKERRQ(ierr);
		ierr=MatSetValue(Pm_global,n3,n2,Ps_local[2][1],ADD_VALUES);
		CHKERRQ(ierr);
		ierr=MatSetValue(Pm_global,n3,n3,Ps_local[2][2],ADD_VALUES);
		CHKERRQ(ierr);

	}
	fileElements.close();

	// free memory Ps_local
	delete[] Ps_local;


	MatAssemblyBegin(Pm_global,MAT_FINAL_ASSEMBLY);
	MatAssemblyEnd(Pm_global,MAT_FINAL_ASSEMBLY);
//	ierr = MatView(Pm_global,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);

	return 1;
}

  ///////////////////////////
 // POST MATRIX FUNCTIONS //
///////////////////////////

/*! La funcion LOCAL_MATRIX() es mayormente el trozo de algoritmo que cambia al cambiar
 * el campo de aplicacion del programa (magnetismo,distribucion termica, mecanica).
 * En esta funcion esta reflejada en gran medida la funcion fisica que representa
 * la base del programa Femsii.
 * Se encarga de calcular la matriz local de cada elemento, aceptando la posicion
 * de cada uno de los tres puntos, y propiedades aplicables a cada elemento/region. En
 * nuestro caso la permeabilidad magnetica.
 */
int Sistema_ecuaciones::local_matrix(double xi,double yi,double xj,double yj,double xk,double yk,int region )
{
	double W,bi,bj,bk,ci,cj,ck,per,con,delta,cte;

	// calculating values "b" and "c" for the actual element
	bi=yj-yk; bj=yk-yi; bk=yi-yj; ci=xk-xj; cj=xi-xk; ck=xj-xi;

	// getting the area by determinat
	delta=fabs(0.5*(xi*yj+yi*xk+xj*yk-xk*yj-xj*yi-yk*xi));

	// permeabily can't be less than 1
	if(permeabilidad[region]<1){per=1;cout<<"Error en Permeabilidad";exit(0);}else{per=permeabilidad[region];}
	// conductivity can't be less than 0
	if(conductividad[region]<0){con=0;cout<<"Error en Conductividad";exit(0);}else{con=conductividad[region];}
	
	// Stiffness matrix
	// REAL PART
	cte=4.0*delta*VACIO*per;
	PetscRealPart( Ps_local[0][0] )= (bi*bi+ci*ci)/cte;
	PetscRealPart( Ps_local[0][1] )= (bi*bj+ci*cj)/cte;
	PetscRealPart( Ps_local[0][2] )= (bi*bk+ci*ck)/cte;
	PetscRealPart( Ps_local[1][0] )= (bi*bj+ci*cj)/cte;
	PetscRealPart( Ps_local[1][1] )= (bj*bj+cj*cj)/cte;
	PetscRealPart( Ps_local[1][2] )= (bj*bk+cj*ck)/cte;
	PetscRealPart( Ps_local[2][0] )= (bi*bk+ci*ck)/cte;
	PetscRealPart( Ps_local[2][1] )= (bj*bk+cj*ck)/cte;
	PetscRealPart( Ps_local[2][2] )= (bk*bk+ck*ck)/cte;

	// IMAGINARY PART
	if(frec!=0.0){
		W=2*PI*frec;
		cte=W*con*delta/12.0;
		PetscImaginaryPart( Ps_local[0][0] )= 2.0*cte;
		PetscImaginaryPart( Ps_local[0][1] )= cte;
		PetscImaginaryPart( Ps_local[0][2] )= cte;
		PetscImaginaryPart( Ps_local[1][0] )= cte;
		PetscImaginaryPart( Ps_local[1][1] )= 2.0*cte;
		PetscImaginaryPart( Ps_local[1][2] )= cte;
		PetscImaginaryPart( Ps_local[2][0] )= cte;
		PetscImaginaryPart( Ps_local[2][1] )= cte;
		PetscImaginaryPart( Ps_local[2][2] )= 2.0*cte;
	}
	return 1;
}


/*! REGIONES() Se trata de que esta funcion llene las variables de
 * tipo map "permeabilidad" y "conductividad" con los valores de conductividad y
 * permeabilidad en cada region que da el usuario.
  * Necesario para la funcion local_matrix() .
 */
int Sistema_ecuaciones::regiones()
{
	int N_R, region;
	double valor;
	char verifi;

	// PERMEABILITY
	ifstream filePermeabi;
	filePermeabi.open(enrrutar("/permeabilidad.reg"), ios::in);	// input summary file
	filePermeabi >> N_R;
	for(int i=0;i<N_R;i++){
		filePermeabi >> region >> valor;
		permeabilidad[region]=valor;
	}
	filePermeabi >> verifi;
	if (verifi!='F'){cout<<"P_sistemaEc: Archivo "<<enrrutar("/permeabilidad.reg")<<" corrupto\n"; exit(0);}
	filePermeabi.close();

	// CONDUCTIVITY
	ifstream fileConducti;
	fileConducti.open(enrrutar("/conductividad.reg"), ios::in);	// input summary file
	fileConducti >> N_R;
	for(int i=0;i<N_R;i++){
		fileConducti >> region >> valor;
		conductividad[region]=valor;
	}
	fileConducti >> verifi;
	if (verifi!='F'){cout<<"P_sistemaEc: Archivo "<<enrrutar("/conductividad.reg")<<" corrupto\n"; exit(0);}
	fileConducti.close();
	return 1;
}

/*! PRINT_GLOBAL() muestra en la salida estandar la matriz de rigidez global.
 */
void Sistema_ecuaciones::print_global()
{
	MatView(Pm_global,PETSC_VIEWER_STDOUT_SELF);
}


/*! RIGHT-SIDE() genera Pv_rightside a partir del archivo de densidades de corriente.
 * En el archivo de densidades de corriente se nos proporcionara una lista de
 * regiones y densidades correspondientes. Mediante el archivo g2f.ele se
 * obtienen los nodos asociados a cada region.
 * Tener en cuenta que el archivo de densidades de corriente incluye
 * siempre parte real y parte imaginaria para cada region.
 */
int Sistema_ecuaciones::rightSide()
{
	char verifi;
	int N_C,region,ignorar, *n1, *n2, *n3, *reg; // number of conditions, working node
	double *area, real, imag, xi, yi, xj, yj, xk, yk;
	ifstream fileDensidadCorri, fileElements;
	fileDensidadCorri.open(enrrutar("/densidadCorri.reg"), ios::in);
		if(!fileDensidadCorri.is_open()){cout<<"P_sistemaEc: Error leyendo archivo de condiciones\n";exit(0);}
	fileDensidadCorri >> N_C;
	cout << "P_sistemaEc: Numero de condiciones de contorno: " << N_C << endl;

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

	// -- PETSC -- //
	PetscScalar Ps_valor, cero;
	VecCreate(PETSC_COMM_WORLD,&Pv_rightside);
	PetscObjectSetName((PetscObject) Pv_rightside, "Solution");
	VecSetSizes(Pv_rightside,PETSC_DECIDE,N_N);
	VecSetFromOptions(Pv_rightside);

	n1 = new int[N_E];
	n2 = new int[N_E];
	n3 = new int[N_E];
	reg = new int[N_E];
	area = new double[N_E];
	PetscRealPart(cero)=0.0;
	PetscImaginaryPart(cero)=0.0;

	for(int j=0; j<N_E; j++){
		fileElements >> ignorar >> n1[j] >> n2[j] >> n3[j] >> reg[j];
		xi=pNode[n1[j]][0]; yi=pNode[n1[j]][1]; // Node 1 coord x e y
		xj=pNode[n2[j]][0]; yj=pNode[n2[j]][1]; // Node 2 coord x e y
		xk=pNode[n3[j]][0]; yk=pNode[n3[j]][1]; // Node 3 coord x e y
		area[j]=fabs(0.5*(xi*yj+yi*xk+xj*yk-xk*yj-xj*yi-yk*xi));
	}

	for(int i=0; i<N_C; i++){

		fileDensidadCorri >> region;
		fileDensidadCorri >> real;
		fileDensidadCorri >> imag;

		for(int j=0; j<N_E; j++){
			if( region==reg[j] ){
				PetscRealPart(Ps_valor)=real*area[j]/3;
				PetscImaginaryPart(Ps_valor)=imag*area[j]/3;
				ierr=VecSetValue(Pv_rightside,n1[j],Ps_valor,ADD_VALUES);CHKERRQ(ierr);
				ierr=VecSetValue(Pv_rightside,n2[j],Ps_valor,ADD_VALUES);CHKERRQ(ierr);
				ierr=VecSetValue(Pv_rightside,n3[j],Ps_valor,ADD_VALUES);CHKERRQ(ierr);
			}
		}
	}

	fileDensidadCorri >> verifi;
	if (verifi!='F'){cout<<"P_sistemaEc: Archivo "<<enrrutar("/densidadCorri.reg")<<" corrupto\n";
		fileDensidadCorri.close();
		exit(0);
	}

	ierr=VecAssemblyBegin(Pv_rightside);CHKERRQ(ierr);
	ierr=VecAssemblyEnd(Pv_rightside);CHKERRQ(ierr);

	fileDensidadCorri.close();
	fileElements.close();
	cout << "P_sistemaEc: Densidades de corriente aplicadas.\n";
	
	// Liberar la memoria de n1 n2 n3
	delete n1;
	delete n2;
	delete n3;
	delete reg;

	return 1;
}


/*! DIRICHLET-PTO() aplica dirichlet en Pm_global a patir de valores en puntos.
* ESTA funcion da problemas de rentidimiento. Por ahora se omite el uso de
* Dirichlet sobre una linea.
* CUIDADO: funcion prototipo con vistas a ser implementada en un futuro
* proyecto. No implemetnada en Femsii 1.0
*/
int Sistema_ecuaciones::dirichletPto()
{
	ifstream fileDirichlet;
	int N_C,node; // number of conditions, working node
	PetscScalar Ps_valor,Ps_zero,Ps_one;

	PetscRealPart(Ps_zero)=0.0; PetscImaginaryPart(Ps_zero)=0.0;
	PetscRealPart(Ps_one)=1.0; PetscImaginaryPart(Ps_one)=1.0;
	
	fileDirichlet.open(enrrutar("/dirichlet.pto"), ios::in);	// input summary file
		if(!fileDirichlet.is_open())
			{cout<<"P_sistemaEc: Error leyendo archivo de dirichlet\n";exit(0);}
	fileDirichlet >> N_C;

	for(int i=0;i<N_C;i++){
		fileDirichlet >> node;
		fileDirichlet >> Ps_valor;
		ierr=MatZeroRows(Pm_global,1,&node,Ps_one);CHKERRQ(ierr);
		VecSetValue(Pv_rightside,node,Ps_valor,INSERT_VALUES);
	}
	VecAssemblyBegin(Pv_rightside);
	VecAssemblyEnd(Pv_rightside);
	MatAssemblyBegin(Pm_global,MAT_FINAL_ASSEMBLY);
	MatAssemblyEnd(Pm_global,MAT_FINAL_ASSEMBLY);
	fileDirichlet.close();
	cout << "P_sistemaEc: Condiciones de Dirichlet aplicadas.\n\n";

	ierr = MatView(Pm_global,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);
	return 1;
}


/*! DIRICHLET-LIN() aplica dirichlet en Pm_global a patir de valores en lineas.
* Esta funcion debe realizar un barrido sobre el archivo dirichlet.lin para
* encontrar encontrar las regiones donde existe de antemano un valor dado
* del potencial vector magnetico.
* Para ello busca los nodos de las regiones en el archivo g2f.edge, donde se
* encuentran los elementos de tipo linea. Una vez encontrados los nodos
* asociados fuerza al sistema de ecuaciones (Pm_global y Pv_righside) para
* asignar el valor adecuado al nodo correspondiente.
*/
int Sistema_ecuaciones::dirichletLin()
{
	int N_C,region,ignorar, *n1, *n2, *reg, node;
	ifstream fileDirichlet, fileEdge;
	PetscScalar Ps_valor;
	
	fileDirichlet.open(enrrutar("/dirichlet.lin"), ios::in);	// input summary file
		if(!fileDirichlet.is_open())
			{cout<<"P_sistemaEc: Error leyendo archivo de dirichlet\n";exit(0);}
	fileDirichlet >> N_C;

	fileEdge.open(enrrutar("/temporal/g2f.edge"), ios::in);	// input edge file
		if(!fileEdge.is_open())
			{cout<<"P_sistemaEc: Error leyendo archivo /temporal/g2f.edge\n";exit(0);}

	n1 = new int[N_ED];
	n2 = new int[N_ED];
	reg = new int[N_ED];

	for(int j=0; j<N_ED; j++){
		fileEdge >> ignorar >> n1[j] >> n2[j] >> reg[j];
	}


	for(int i=0; i<N_C; i++){
		fileDirichlet >> region;
		fileDirichlet >> Ps_valor;
		for(int j=0; j<N_ED; j++){
			if(region==reg[j]){
				node = n1[j];
				ierr=MatZeroRows(Pm_global,1,&node,1.0);CHKERRQ(ierr);
				ierr=VecSetValue(Pv_rightside,node,Ps_valor,INSERT_VALUES);CHKERRQ(ierr);
				node = n2[j];
				ierr=MatZeroRows(Pm_global,1,&node,1.0);CHKERRQ(ierr);
				ierr=VecSetValue(Pv_rightside,node,Ps_valor,INSERT_VALUES);CHKERRQ(ierr);
			}
		}
	}
	
	VecAssemblyBegin(Pv_rightside);
	VecAssemblyEnd(Pv_rightside);
	MatAssemblyBegin(Pm_global,MAT_FINAL_ASSEMBLY);
	MatAssemblyEnd(Pm_global,MAT_FINAL_ASSEMBLY);
	
	fileDirichlet.close();
	fileEdge.close();
	
	cout << "P_sistemaEc: Condiciones de Dirichlet aplicadas.\n\n";
	return 1;
}


/*! PARES-HOMOLOGOS() Se trata de leer un archivo en el que hay un listado de regiones homologas.
 * Con ello se debe crear un puntero doble con los nodos relacionados. El procedimiento
 * es leer el archivo que relaciona las regiones y buscar en el archivo g2f.edge
 * los nodos que perteneces a esa region.
 * Cuidado, funcion prototipo EN DESARROLLO, no incluido en Femsii 1.0
*/
int Sistema_ecuaciones::paresHomologos()
{
	int ignorar, NE, NR,**ptrNodoEdge, **ptrRegHomo, **ptrNodoHomo;
	ifstream fileEdge,fileSumary, fileKelvin;
	ofstream fileOut;

	fileKelvin.open(enrrutar("/kelvin.reg"), ios::in);		// input homologian regions file
		if(!fileKelvin.is_open())
			{cout<<"P_sistemaEc: Error leyendo archivo 'kelvin.reg'\n";exit(0);}

	fileSumary.open(enrrutar("/temporal/g2f.ene"), ios::in);	// input summary file
		if(!fileSumary.is_open())
			{cout<<"P_sistemaEc: Error leyendo archivo '/temporal/g2f.ene'\n";exit(0);}

	fileEdge.open(enrrutar("/temporal/g2f.edge"), ios::in);	// input edge file
		if(!fileEdge.is_open())
			{cout<<"P_sistemaEc: Error leyendo archivo '/temporal/g2f.edge'\n";exit(0);}

	fileOut.open(enrrutar("/temporal/homo.node"), ios::out);	// input edge file
		if(!fileOut.is_open())
			{cout<<"P_sistemaEc: Error al abrir el archivo '/temporal/homo.node'\n";exit(0);}

	cout << "P_sistemaEc: Creando archivo de nodos homologos \"/temporal/homo.node\" .\n";

	// all kelvin.reg file to ptrNodoHomo pointer
	fileKelvin >> NR;
	ptrRegHomo = new int*[NR]; for(int i=0;i<NR;i++) ptrRegHomo[i] = new int[2];
	for( int i=0;i<NR;i++ ) fileKelvin >> ptrRegHomo[i][0] >> ptrRegHomo[i][1];
	fileKelvin.close();

	// get node number from g2f.ene and use it to put the domention of ptrNodoHomo and ptrNodoEdge
	fileSumary >> ignorar >> ignorar >> NE; // nos vale la tercera lectura, la de los edge
	ptrNodoHomo = new int*[NE/8+1]; for(int i=0;i<NE/8+1;i++) ptrNodoHomo[i] = new int[2];
	ptrNodoEdge = new int*[NE+1]; for(int i=0;i<NE+1;i++) ptrNodoEdge[i] = new int[4];
	fileSumary.close();

	// all g2f.edge pNodeEdge pointer
	for( int i=0;i<NE;i++ )
		fileEdge >>ptrNodoEdge[i][0]>>ptrNodoEdge[i][1]>>ptrNodoEdge[i][2]>>ptrNodoEdge[i][3];
	fileEdge.close();

	for( int i=0;i<NR;i++ ){
		regNodos(ptrNodoEdge, ptrRegHomo[i][0], 0); // 0 define use list listaRegNodos0
		regNodos(ptrNodoEdge, ptrRegHomo[i][1], 1); // 1 define use list listaRegNodos1
	}

	// Solo para imprimir los resultados
	fileOut << listaRegNodos0.size() << endl;
	list<int>::iterator it0 = listaRegNodos0.begin();
	list<int>::iterator it1 = listaRegNodos1.begin();
	for( it0 = listaRegNodos0.begin(); it0 != listaRegNodos0.end(); it0++){
		fileOut << *it0 << ' ' << *it1 << endl; // it1 se incrementa?? !!!!! CUIDADO mirar
		it1++;
	}
	fileOut << 'F';

	fileOut.close();
	cout << "P_sistemaEc: Lista de nodos homologos creada.\n";
	return 1;
}


/*! REG-NODOS() es llamado varias veces por paresHomologos() busca los nodos de
 * que pertenezcan a una region y los coloca en listaRegNodos0 en caso de que la
 * llamada sea para los nodos de la region origen, y listaRegNodos1 si los nodos
 * son de la region homologa.
 * Cuidado, funcion prototipo EN DESARROLLO, no incluido en Femsii 1.0
 */
int Sistema_ecuaciones::regNodos(int **ptrNodoEdge, int region, int homo)
{ //CUIDADO, COMPROBAR SI CON EL CAMBIO A ELEMENTOS EMPEZANDO EN CERO HAY QUE RETOCAR
	int salir=0, i=0;
	if( homo==0 ){
		while( salir!=1 ){
			if( ptrNodoEdge[i][3]==region){
				i++; // jump to the second node of each region
				while( ptrNodoEdge[i][3]==region ){
					listaRegNodos0.push_back( ptrNodoEdge[i][1] );
					i++;
				}
				listaRegNodos0.push_back( ptrNodoEdge[i-1][2] );
				salir=1;
			}
			i++; // i want to insert en the list every node but the first
		}
	}
	if( homo==1 ){
		while( salir!=1 ){
			if( ptrNodoEdge[i][3]==region){
				i++;
				while( ptrNodoEdge[i][3]==region ){
					listaRegNodos1.push_back( ptrNodoEdge[i][1] );
					i++;
				}
				listaRegNodos1.push_back( ptrNodoEdge[i-1][2] );
				salir=1;
			}
			i++;
		}
	}
	if( homo!=0 && homo!=1 ){ cout << "Error en la funcion Sistema_ecuaciones::regNodos\n"; exit(0);}

	return 1;
}


/*! Funcion para implementar las condiciones kelvin en Femsii atendiado a la
* salida de paresHomologos().
* Cuidado, funcion prototipo EN DESARROLLO, no incluido en Femsii 1.0
*/
// es muy importante que solo aparezca una vez cada nodo en homo.node
int Sistema_ecuaciones::kelvin()
{
	int N_H;// numero de homologos
//	int n1, n2, *ptrFila;
	PetscInt n1, n2, pNcols;
	const PetscInt *pCols;
	ifstream fileHomo;
	PetscScalar *positivo, *negativo, *cero;
	const PetscScalar *pVals;

	positivo = new PetscScalar[1]; PetscRealPart(positivo[0])=1.0; PetscImaginaryPart(positivo[0])=1.0;
	negativo = new PetscScalar[1]; PetscRealPart(negativo[0])=-1.0; PetscImaginaryPart(negativo[0])=-1.0;
	cero = new PetscScalar[1]; PetscRealPart(cero[0])=0.0; PetscImaginaryPart(cero[0])=0.0;
	
	fileHomo.open(enrrutar("/temporal/homo.node"), ios::in);	// input edge file
	if(!fileHomo.is_open())
		{cout<<"P_sistemaEc: Error al abrir el archivo '/temporal/homo.node'\n";exit(0);}

	cout << "hola\n";
//	for(int i=0; i<N_N; i++) ptrFila[i]=i;

	fileHomo >> N_H;
	for(int i=0; i<N_H; i++){
		fileHomo >> n1 >> n2;

		MatGetRow(Pm_global,n2,&pNcols,&pCols,&pVals);
		MatSetValues(Pm_global,1,&n1,pNcols,pCols,pVals,ADD_VALUES);
//		MatRestoreRow(Pm_global,n2,&pNcols,&pCols,&pVals);

		MatAssemblyBegin(Pm_global,MAT_FINAL_ASSEMBLY);
		MatAssemblyEnd(Pm_global,MAT_FINAL_ASSEMBLY);

		ierr=MatZeroRows(Pm_global,1,&n2,positivo[0]);CHKERRQ(ierr);

		ierr=MatSetValues(Pm_global,1,&n2,1,&n1,&negativo[0],ADD_VALUES);CHKERRQ(ierr);
		MatAssemblyBegin(Pm_global,MAT_FINAL_ASSEMBLY);
		MatAssemblyEnd(Pm_global,MAT_FINAL_ASSEMBLY);
		
		ierr=VecSetValue(Pv_rightside,n2,cero[0],INSERT_VALUES);CHKERRQ(ierr);

		VecAssemblyBegin(Pv_rightside);
		VecAssemblyEnd(Pv_rightside);
	}


	cout << "P_sistemaEc: Condiciones kelvin aplicadas.\n";
	fileHomo.close();
	return 1;
}


/*! aBINARIO() crea un archivo binario especifico PETSC con el contenido de una
 * matriz, respetando la estructura y el tipo. Asi mismo tambien se crea un binario
 * con el contenido del vactor de terminos independientes.
 * Estos archivos se crean para ser leidos posteriormente por solve.bin .
 */
int Sistema_ecuaciones::aBinario()
{
	PetscViewer viewerM, viewerV;
	ierr=PetscViewerBinaryOpen(PETSC_COMM_WORLD,enrrutar("/temporal/femsii_b"),FILE_MODE_WRITE,&viewerV);CHKERRQ(ierr);
	ierr=PetscViewerBinaryOpen(PETSC_COMM_WORLD,enrrutar("/temporal/femsii_a"),FILE_MODE_WRITE,&viewerM);CHKERRQ(ierr);
	ierr=VecView(Pv_rightside,viewerV);CHKERRQ(ierr);
	ierr=MatView(Pm_global,viewerM);CHKERRQ(ierr);
	ierr=PetscViewerDestroy(viewerV);CHKERRQ(ierr);
	ierr=PetscViewerDestroy(viewerM);CHKERRQ(ierr);
	ierr=PetscPrintf(PETSC_COMM_WORLD,"Matriz de rigidez y terminos independiente a un binario temporal\n");CHKERRQ(ierr);
	return 1;
}


void _init()
{
	printf("Inside _init()\n"); 
} 
void _fini()
{
	printf("Inside _fini()\n"); 
}
