#include "Gmsh.h"
#include <cstring>
#include <cstdlib>


Gmsh::Gmsh()
{
	_if_is_open=false,_of_is_open=false; _dim=0; _nb_node=0;_nb_element=0;
}  

Gmsh::Gmsh(const char *filename)
{
    Gmsh::readMesh(filename);
    
}

Gmsh::~Gmsh()
{
	Gmsh::close();
} 
void Gmsh::show()
{
	cout<<"number of node:"<<_nb_node<<"number of element:"<<_nb_element<<" dim="<<_dim<<endl;    
	for(size_t i=0;i<_nb_node;i++)
    {
		cout<<"node "<<i<<" with coordinate"<< _node[i].coord.x<<"," << _node[i].coord.y
		<<"," << _node[i].coord.z<<endl;
    }
	for(size_t i=0;i<_nb_element;i++)
    {
		cout<<"element "<<_element[i].label<<" of type"<< _element[i].type_element
		<<" and nb node:" << _element[i].nb_node<<" code elementary:"<<_element[i].code_elementary
		<<", code physical:"<<_element[i].code_physical<<endl;
		cout<<"node labels are:";
		for(size_t j=0;j<_element[i].nb_node;j++)
			cout<<" "<<_element[i].node[j];    
		cout<<endl;
    }
}
/*--------- read and write mesh -----------*/
void Gmsh::readMesh(const char *filename)
{
	if(_if_is_open)
	{
		cout<<" error opening file:"<<filename<<"  in "<<__FILE__<< " line "<<__LINE__<<endl; 
		cout<<" an open file exists yet"<<endl;
		exit(0);
	}    
	_if.open(filename);
	if(!_if.is_open())
	{
		cout<<" error opening file:"<<filename<<"  in "<<__FILE__<< " line "<<__LINE__<<endl; 
		cout<<" file does not exist"<<endl;
		exit(0);
	}
	_if_is_open=true;    
	_dim=0;
	size_t filesize;
	_if.seekg( -1, ios_base::end );
	filesize=_if.tellg();
	char *ptr,*file_string;
	file_string=(char *)malloc(filesize+2);
	size_t nb_code,pos;
	_if.seekg(0);
	pos=0;
	while (!_if.eof())// read the file
	{
		file_string[pos]=(char) _if.get();
		pos++;
	}     
	
	// the file is now in the string
	ptr=file_string;
	ptr= strstr( ptr,"$Nodes")+6;
	_nb_node= strtod(ptr, &ptr);
	_node.resize(_nb_node); 
	for(size_t i=0;i<_nb_node;i++)
    {
		_node[i].label  = strtod(ptr, &ptr);
		_node[i].coord.x= strtod(ptr, &ptr);
		_node[i].coord.y= strtod(ptr, &ptr);     
		_node[i].coord.z= strtod(ptr, &ptr);     
    }
	ptr= strstr( ptr,"$Elements")+9;
	_nb_element= strtod(ptr, &ptr);
	_element.resize(_nb_element); 
	for(size_t i=0;i<_nb_element;i++)
    {
		_element[i].label  = strtod(ptr, &ptr);
		_element[i].type_element= strtod(ptr, &ptr);
		nb_code = strtod(ptr, &ptr); 
		if((nb_code < 1))
		{cout<<"invalid number of code in msh file"<<endl;exit(0);}
		if(nb_code ==1)   
		{_element[i].code_physical=0;}
		else
			_element[i].code_physical= strtod(ptr, &ptr);     
		_element[i].code_elementary= strtod(ptr, &ptr);     
		if (nb_code >2) // eliminate the other tag if any
		{
			for(size_t j=0;j<nb_code-2;j++)  pos=strtod(ptr, &ptr);
		}
		switch(_element[i].type_element)
		{
			case 15:   // node
				_element[i].nb_node=1;
				_element[i].dim=0;
				break;
			case 1:  // line
			case 8:
				_element[i].nb_node=2;
				_element[i].dim=1;
				break;
			case 2:  //triangle
			case 9:    
				_element[i].nb_node=3;
				_element[i].dim=2;
				break;
			case 3:  // quadrangle
			case 10:    
				_element[i].nb_node=4;
				_element[i].dim=2;
				break;
			case 4:  //tetrahedron
			case 11:
				_element[i].nb_node=4;
				_element[i].dim=3;
				break;
			case 5:  // hexahedron
			case 12:
				_element[i].nb_node=8;
				_element[i].dim=3;
				break;
			case 6:  // prism
			case 13:
				_element[i].nb_node=6;
				_element[i].dim=3;
				break;
			case 7:  // pyramid
			case 14:
				_element[i].nb_node=5;
				_element[i].dim=3;
				break;    
			default:
				cout<<"error code element in file msh"<<endl; exit(0);
				break;
		}
		if (_dim<  _element[i].dim) _dim= _element[i].dim;  
		for(size_t j=0;j<_element[i].nb_node;j++) _element[i].node[j]= strtod(ptr, &ptr);
    }
	free(file_string);    
	// we have fill all the gmsh structure        
}

void Gmsh::writeMesh(const char *filename)
{
	if(_of_is_open==true)
    { 
		cout<<"a file is still open: close it before writing a mesh"<<endl;
    }
	
	if((_nb_node==0) || (_nb_element==0))
	{
		cout<<" error in file:"<<filename<<"  in "<<__FILE__<< " line "<<__LINE__<<endl; 
		cout<<" there is no mesh to save"<<endl;
		return;
	}       
	_of.open(filename);
	if(!_of.is_open())
	{
		cout<<" error in file:"<<filename<<"  in "<<__FILE__<< " line "<<__LINE__<<endl; 
		cout<<" can not create the file"<<endl;
		return;
	}
	_of_is_open=true;_nb_save=0;           
	_of<<"$MeshFormat"<<endl;
	_of<<"2.2 0 8"<<endl;
	_of<<"$EndMeshFormat"<<endl;
	_of<<"$Nodes"<<endl;
	_of<<_nb_node<<endl;
	for(size_t i=0;i<_nb_node;i++)
    {
		_of<<setw(FVCHAMPINT)<<_node[i].label;
		_of<<setprecision(FVPRECISION) <<setw(FVCHAMP)<<_node[i].coord.x;
		_of<<setprecision(FVPRECISION) <<setw(FVCHAMP)<<_node[i].coord.y;
		_of<<setprecision(FVPRECISION) <<setw(FVCHAMP)<<_node[i].coord.z<<endl;
    }
	_of<<"$EndNodes"<<endl;
	_of<<"$Elements"<<endl;
	_of<<_nb_element<<endl;
	for(size_t i=0;i<_nb_element;i++)
    {
		_of<<_element[i].label<<setw(FVCHAMPINT)<<_element[i].type_element<<setw(FVCHAMPINT)<<2;
		_of<<setw(FVCHAMPINT)<<_element[i].code_physical<<setw(FVCHAMPINT)<<_element[i].code_elementary;
		for(size_t j=0;j<_element[i].nb_node;j++)
            _of<<setw(FVCHAMPINT)<<_element[i].node[j];
		_of<<endl;
	}
	_of<<"$EndElements"<<endl;
}

void Gmsh::close()
{
	if(_if_is_open) {_if.close();_if_is_open=false;}
	if(_of_is_open) {_of.close();_of_is_open=false;} 
} 



void Gmsh::writeVector(FVVect<double> &u, const size_t type, const char *name, double time)
{
	if(!_of_is_open)
    { 
		cout<<"write the mesh file first"<<endl;
    }   
	if(type== VERTEX)
    {
		_of<<"$NodeData"<<endl;   
		_of<<"1"<<endl;
		_of<<"\""<<name<<"\""<<endl;
		_of<<"1"<<endl;
		_of<<time<<endl; 
		_of<<"3"<<endl;
		_of<<_nb_save++<<endl;
		_of<<"1"<<endl;
		_of<<_nb_node<<endl;
		for(size_t i=0;i<_nb_node;i++)
			_of<<setw(FVCHAMPINT)<<i+1<<setw(FVCHAMP)<<setprecision(FVPRECISION)<<u[i]<<endl;
		_of<<"$EndNodeData"<<endl;   
    }
	
	if(type== CELL)    
    {
		_of<<"$ElementData"<<endl; 
		_of<<"1"<<endl;
		_of<<"\""<<name<<"\""<<endl;
		_of<<"1"<<endl;
		_of<<time<<endl;
		_of<<"3"<<endl;
		_of<<_nb_save++<<endl;
		_of<<"1"<<endl;
		_of<<_nb_element<<endl;
		for(size_t i=0;i<_nb_element;i++)
			_of<<setw(FVCHAMPINT)<<i+1<<setw(FVCHAMP)<<setprecision(FVPRECISION)<<u[i]<<endl;
		_of<<"$EndElementData"<<endl;   
    }
	
}

void Gmsh::writeVector(const FVVect<double> &u,const FVVect<double> &v,  const size_t type,const char *name, double time) 
{
	if(!_of_is_open)
    { 
		cout<<"write the mesh file first"<<endl;
    } 
	if(type== VERTEX)
    {
		_of<<"$NodeData"<<endl;   
		_of<<"1"<<endl;
		_of<<"\""<<name<<"\""<<endl;
		_of<<"1"<<endl;
		_of<<time<<endl; 
		_of<<"3"<<endl;
		_of<<_nb_save++<<endl;
		_of<<"3"<<endl;
		_of<<_nb_node<<endl;
		for(size_t i=0;i<_nb_node;i++)
			_of<<setw(FVCHAMPINT)<<i+1<<setw(FVCHAMP)<<u[i]<<setw(FVCHAMP)<<v[i]<<setw(FVCHAMP)<<"0"<<endl;
		_of<<"$EndNodeData"<<endl;   
    }
    
	if(type== CELL)    
    {
		_of<<"$ElementData"<<endl; 
		_of<<"1"<<endl;
		_of<<"\""<<name<<"\""<<endl;
		_of<<"1"<<endl;
		_of<<time<<endl;
		_of<<"3"<<endl;
		_of<<_nb_save++<<endl;
		_of<<"3"<<endl;
		_of<<_nb_element<<endl;
		for(size_t i=0;i<_nb_element;i++)
			_of<<setw(FVCHAMPINT)<<i+1<<setw(FVCHAMP)<<u[i]<<setw(FVCHAMP)<<v[i]<<setw(FVCHAMP)<<"0"<<endl;
		_of<<"$EndElementData"<<endl;   
    }    
}

void Gmsh::writeVector(const FVVect<double> &u,const FVVect<double> &v, const FVVect<double> &w, const size_t type,const char *name, double time)   
{
	if(!_of_is_open)
    { 
		cout<<"write the mesh file first"<<endl; 
    } 
	if(type== VERTEX)
    {
		_of<<"$NodeData"<<endl;   
		_of<<"1"<<endl;
		_of<<"\""<<name<<"\""<<endl;
		_of<<"1"<<endl;
		_of<<time<<endl; 
		_of<<"3"<<endl;
		_of<<_nb_save++<<endl;
		_of<<"3"<<endl;
		_of<<_nb_node<<endl;
		for(size_t i=0;i<_nb_node;i++)
			_of<<setw(FVCHAMPINT)<<i+1<<setw(FVCHAMP)<<u[i]<<setw(FVCHAMP)<<v[i]<<setw(FVCHAMP)<<w[i]<<endl;
		_of<<"$EndNodeData"<<endl;   
    }
    
	if(type== CELL)    
    {
		_of<<"$ElementData"<<endl; 
		_of<<"1"<<endl;
		_of<<"\""<<name<<"\""<<endl;
		_of<<"1"<<endl;
		_of<<time<<endl;
		_of<<"3"<<endl;
		_of<<_nb_save++<<endl;
		_of<<"3"<<endl;
		_of<<_nb_element<<endl;
		for(size_t i=0;i<_nb_element;i++)
			_of<<setw(FVCHAMPINT)<<i+1<<setw(FVCHAMP)<<u[i]<<setw(FVCHAMP)<<v[i]<<setw(FVCHAMP)<<w[i]<<endl;
		_of<<"$EndElementData"<<endl;   
    }    
}


void Gmsh::writeVector(FVVect<FVPoint1D<double> >&u, const size_t type, const char *name, double time)
{
	if(!_of_is_open)
    { 
		cout<<"write the mesh file first"<<endl;
    }   
	if(type== VERTEX)
    {   
		_of<<"$NodeData"<<endl;   
		_of<<"1"<<endl;
		_of<<"\""<<name<<"\""<<endl;
		_of<<"1"<<endl;
		_of<<time<<endl; 
		_of<<"3"<<endl;
		_of<<_nb_save++<<endl;
		_of<<"1"<<endl;
		_of<<_nb_node<<endl;
		for(size_t i=0;i<_nb_node;i++)
			_of<<setw(FVCHAMPINT)<<i+1<<setw(FVCHAMP)<<u[i].x<<endl;
		_of<<"$EndNodeData"<<endl;   
    }
    
	if(type== CELL)    
    {    
		_of<<"$ElementData"<<endl; 
		_of<<"1"<<endl;
		_of<<"\""<<name<<"\""<<endl;
		_of<<"1"<<endl;
		_of<<time<<endl;
		_of<<"3"<<endl;
		_of<<_nb_save++<<endl;
		_of<<"1"<<endl;
		_of<<_nb_element<<endl;
		for(size_t i=0;i<_nb_element;i++)
			_of<<setw(FVCHAMPINT)<<i+1<<setw(FVCHAMP)<<u[i].x<<endl;
		_of<<"$EndElementData"<<endl;   
    }
	
}

void Gmsh::writeVector(FVVect<FVPoint2D<double> >&u, const size_t type,const char *name, double time) 
{
	if(!_of_is_open)
    { 
		cout<<"write the mesh file first"<<endl;
    } 
	if(type== VERTEX)
    {
		_of<<"$NodeData"<<endl;   
		_of<<"1"<<endl;
		_of<<"\""<<name<<"\""<<endl;
		_of<<"1"<<endl;
		_of<<time<<endl; 
		_of<<"3"<<endl;
		_of<<_nb_save++<<endl;
		_of<<"3"<<endl;
		_of<<_nb_node<<endl;
		for(size_t i=0;i<_nb_node;i++)
			_of<<setw(FVCHAMPINT)<<i+1<<setw(FVCHAMP)<<u[i].x<<setw(FVCHAMP)<<u[i].y<<setw(FVCHAMP)<<"0"<<endl;
		_of<<"$EndNodeData"<<endl;   
    }
    
	if(type== CELL)    
    {  
		_of<<"$ElementData"<<endl; 
		_of<<"1"<<endl;
		_of<<"\""<<name<<"\""<<endl;
		_of<<"1"<<endl;
		_of<<time<<endl;
		_of<<"3"<<endl;
		_of<<_nb_save++<<endl;
		_of<<"3"<<endl;
		_of<<_nb_element<<endl;
		for(size_t i=0;i<_nb_element;i++)
			_of<<setw(FVCHAMPINT)<<i+1<<setw(FVCHAMP)<<u[i].x<<setw(FVCHAMP)<<u[i].y<<setw(FVCHAMP)<<"0"<<endl;
		_of<<"$EndElementData"<<endl;   
    }    
}

void Gmsh::writeVector(FVVect<FVPoint3D<double> >&u, const size_t type,const char *name, double time)   
{
	if(!_of_is_open)
    { 
		cout<<"write the mesh file first"<<endl; 
    } 
	if(type== VERTEX)
    {
		_of<<"$NodeData"<<endl;   
		_of<<"1"<<endl;
		_of<<"\""<<name<<"\""<<endl;
		_of<<"1"<<endl;
		_of<<time<<endl; 
		_of<<"3"<<endl;
		_of<<_nb_save++<<endl;
		_of<<"3"<<endl;
		_of<<_nb_node<<endl;
		for(size_t i=0;i<_nb_node;i++)
			_of<<setw(FVCHAMPINT)<<i+1<<setw(FVCHAMP)<<u[i].x<<setw(FVCHAMP)<<u[i].y<<setw(FVCHAMP)<<u[i].z<<endl;
		_of<<"$EndNodeData"<<endl;   
    }
    
	if(type== CELL)    
    {
		_of<<"$ElementData"<<endl; 
		_of<<"1"<<endl;
		_of<<"\""<<name<<"\""<<endl;
		_of<<"1"<<endl;
		_of<<time<<endl;
		_of<<"3"<<endl;
		_of<<_nb_save++<<endl;
		_of<<"3"<<endl;
		_of<<_nb_element<<endl;
		for(size_t i=0;i<_nb_element;i++)
			_of<<setw(FVCHAMPINT)<<i+1<<setw(FVCHAMP)<<u[i].x<<setw(FVCHAMP)<<u[i].y<<setw(FVCHAMP)<<u[i].z<<endl;
		_of<<"$EndElementData"<<endl;   
    }    
}






















































//
//------------- CONVERTISOR  -----------------
//
void Gmsh::FVMesh2Gmsh(FVMesh1D &m) // constructor which convert a FVMesh1D into a gmah struct
{
	_dim=1;   // a one dimension mesh
	FVVertex1D* ptr_v;    
	_nb_node=m.getNbVertex();
	_node.resize(_nb_node); 
	for((ptr_v=m.beginVertex());(ptr_v=m.nextVertex());)
    {
		size_t i=ptr_v->label;
		_node[i-1].label  = i; 
		_node[i-1].coord.x= ptr_v->coord.x;
		_node[i-1].coord.y= 0;     
		_node[i-1].coord.z= 0; 
    }
	FVCell1D *ptr_c;
	_nb_element=m.getNbCell();
	_element.resize(_nb_element);
	for((ptr_c=m.beginCell());(ptr_c=m.nextCell());)
    {
		size_t i=ptr_c->label;
		_element[i-1].label=i;
		_element[i-1].type_element=1; // we use line2
		_element[i-1].code_physical=ptr_c->code;
		_element[i-1].code_elementary=ptr_c->code;
		_element[i-1].dim=_dim;  
		_element[i-1].nb_node=2;
		_element[i-1].node[0]= ptr_c->firstVertex->label;
		_element[i-1].node[1]= ptr_c->secondVertex->label;
    }
    
}



void Gmsh::FVMesh2Gmsh(FVMesh2D &m) // constructor which convert a FVMesh2D into a gmah struct 
{
	_dim=2;   // a two dimension mesh
	FVVertex2D* ptr_v;    
	_nb_node=m.getNbVertex();
	_node.resize(_nb_node); 
	for((ptr_v=m.beginVertex());(ptr_v=m.nextVertex());)
    {
		size_t i=ptr_v->label;
		_node[i-1].label  = i; 
		_node[i-1].coord.x= ptr_v->coord.x;
		_node[i-1].coord.y= ptr_v->coord.y;     
		_node[i-1].coord.z= 0; 
    }
	FVCell2D *ptr_c;
	double order[4],doux;
	size_t noux;
	FVPoint2D<double> Paux;
	_nb_element=m.getNbCell();
	_element.resize(_nb_element);
	for((ptr_c=m.beginCell());(ptr_c=m.nextCell());)
    {
		size_t i=ptr_c->label;
		_element[i-1].label=i;
		switch(ptr_c->nb_vertex)
        {
			case 3:
				_element[i-1].type_element=2; // we use triangle
				break;
			case 4:
				_element[i-1].type_element=3; // we use quadrangle
				break;
			default:
				cout<<"cell does not correspond to a gmsh element. Abort convertion"<<endl;
				_nb_node=0;_nb_element=0;
				return;
        }    
		_element[i-1].code_physical=ptr_c->code;
		_element[i-1].code_elementary=ptr_c->code;
		_element[i-1].dim=_dim;  
		_element[i-1].nb_node=ptr_c->nb_vertex;
		for(size_t j=0;j<ptr_c->nb_vertex;j++)
        {
			_element[i-1].node[j]=ptr_c->vertex[j]->label; 
			Paux.x=_node[_element[i-1].node[j]-1].coord.x-ptr_c->centroid.x;
			Paux.y=_node[_element[i-1].node[j]-1].coord.y-ptr_c->centroid.y;
			Paux/=Norm(Paux);
			order[j]=(1-Paux.x)*0.5;if(Paux.y<0) order[j]*=-1;   
        }  
		// reordering the node (bubble sort) to cast in the gmah framework  following order
		for(size_t j=ptr_c->nb_vertex-1;j>0;j--)
        {
			for(size_t k=0;k<j;k++)
            {
				if(order[k+1]<order[k])
				{
					doux=order[k+1];noux=_element[i-1].node[k+1];
					order[k+1]=order[k];_element[i-1].node[k+1]=_element[i-1].node[k];
					order[k]=doux;_element[i-1].node[k]=noux;
				}
            }
        }   
    } 
	
}



void Gmsh::FVMesh2Gmsh(FVMesh2D &m,size_t level) // constructor which convert a FVMesh2D into a gmah struct 
{
	_dim=2;   // a two dimension mesh
	FVVertex2D *ptr_v;    
	FVEdge2D *ptr_e;
	FVCell2D *ptr_c;
	double order[4],doux;
	size_t noux,shift_vertex,shift_edge;
	FVPoint2D<double> Paux;
	_nb_node=m.getNbVertex();
	_node.resize(_nb_node); 
	for((ptr_v=m.beginVertex());(ptr_v=m.nextVertex());)
    {
		size_t i=ptr_v->label;
		_node[i-1].label=i; 
		_node[i-1].coord.x=ptr_v->coord.x;
		_node[i-1].coord.y=ptr_v->coord.y;     
		_node[i-1].coord.z=0.0; 
    }
	_nb_element=m.getNbCell()+m.getNbVertex()+m.getNbEdge();
	_element.resize(_nb_element);
	//------- add the vertex -------
	shift_vertex=0; 
	for((ptr_v=m.beginVertex());(ptr_v=m.nextVertex());)
    { 
		size_t i=ptr_v->label;shift_vertex++;
		_element[i-1].label=i;
		_element[i-1].code_physical=ptr_v->code;
		_element[i-1].code_elementary=ptr_v->code;
		_element[i-1].dim=_dim;  
		_element[i-1].nb_node=1;
		_element[i-1].type_element=15;
		_element[i-1].node[0]=ptr_v->label;
    }     
	// ------- add the edge --------
	shift_edge=0;
	for((ptr_e=m.beginEdge());(ptr_e=m.nextEdge());)
    {        
		size_t i=ptr_e->label+shift_vertex;shift_edge++;
		_element[i-1].label=i+shift_vertex;
		_element[i-1].code_physical=ptr_e->code;
		_element[i-1].code_elementary=ptr_e->code;
		_element[i-1].dim=_dim;  
		_element[i-1].nb_node=2;
		_element[i-1].type_element=1;
		_element[i-1].node[0]=ptr_e->firstVertex->label;
		_element[i-1].node[1]=ptr_e->secondVertex->label;
    } 
	// -------- add the cell -------
	for((ptr_c=m.beginCell());(ptr_c=m.nextCell());)
    {
		size_t i=ptr_c->label+shift_vertex+shift_edge;
		_element[i-1].label=i+shift_vertex+shift_edge;
		_element[i-1].code_physical=ptr_c->code;
		_element[i-1].code_elementary=ptr_c->code;
		_element[i-1].dim=_dim;  
		_element[i-1].nb_node=ptr_c->nb_vertex;       
		switch(ptr_c->nb_vertex)
        {
			case 3:
				_element[i-1].type_element=2; // we use triangle
				break;
			case 4:
				_element[i-1].type_element=3; // we use quadrangle
				break;
			default:
				cout<<"2D cell "<<ptr_c->label<<" does not correspond to a gmsh element. "<<endl;
				cout<<"nb vertex="<<ptr_c->nb_vertex<<endl;
				_nb_node=0;_nb_element=0;
				return;
        }    
		for(size_t j=0;j<ptr_c->nb_vertex;j++)
        {
			_element[i-1].node[j]=ptr_c->vertex[j]->label; 
			Paux.x=_node[_element[i-1].node[j]-1].coord.x-ptr_c->centroid.x;
			Paux.y=_node[_element[i-1].node[j]-1].coord.y-ptr_c->centroid.y;
			Paux/=Norm(Paux);
			order[j]=(1-Paux.x)*0.5;if(Paux.y<0) order[j]*=-1;   
        }  
		// reordering the node (bubble sort) to cast in the gmah framework  following order
		for(size_t j=ptr_c->nb_vertex-1;j>0;j--)
        {
			for(size_t k=0;k<j;k++)
            {
				if(order[k+1]<order[k])
				{
					doux=order[k+1];noux=_element[i-1].node[k+1];
					order[k+1]=order[k];_element[i-1].node[k+1]=_element[i-1].node[k];
					order[k]=doux;_element[i-1].node[k]=noux;
				}
            }
        }   
    } 
	
}






void Gmsh::FVMesh2Gmsh(FVMesh3D &m) // constructor which convert a FVMesh3D into a gmah struct
{
	_dim=3;   // a three dimension mesh
	FVVertex3D* ptr_v;
	FVEdge3D* ptr_e;
	FVFace3D* ptr_f;
	FVCell3D *ptr_c;
	FVPoint3D<double> Paux,t1,t2;
	double order[5],doux;
	size_t noux;
	size_t it;
	_nb_node=m.getNbVertex();
	_node.resize(_nb_node); 
	//cout<<"find "<<_nb_node<<" nodes"<<endl;fflush(NULL);
	for((ptr_v=m.beginVertex());(ptr_v=m.nextVertex());)
    {
		size_t i=ptr_v->label;
		_node[i-1].label  = i; 
		_node[i-1].coord.x= ptr_v->coord.x;
		_node[i-1].coord.y= ptr_v->coord.y;     
		_node[i-1].coord.z= ptr_v->coord.z;        
    }
	_nb_element=m.getNbCell();
	_element.resize(_nb_element);
	//cout<<"find "<<_nb_element<<" elements"<<endl;fflush(NULL);
	for((ptr_c=m.beginCell());(ptr_c=m.nextCell());)
    {       
		size_t i=ptr_c->label;
		_element[i-1].label=i;
		_element[i-1].code_physical=ptr_c->code;
		_element[i-1].code_elementary=ptr_c->code;
		_element[i-1].dim=_dim;  
		_element[i-1].nb_node=ptr_c->nb_vertex;
		//cout<<"element "<<i<<" with "<<ptr_c->nb_vertex<<" nodes"<<endl;fflush(NULL);
		switch(_element[i-1].nb_node)
        {
			case 4:   
				_element[i-1].type_element=4; // we use tetrahedron
				for(size_t j=0;j<4;j++)  // the order is not important
					_element[i-1].node[j]=ptr_c->vertex[j]->label;
				break;
			case 5:  
				_element[i-1].type_element=7; // pyramid
				it=0;
				for((ptr_f=ptr_c->beginFace());(ptr_f=ptr_c->nextFace());it++)
				{  
					if( ptr_c->face[it]->nb_vertex==4) break; 
				}
				if (it>=ptr_c->nb_face)  {cout<<" not a pyramidal element"<<endl;return;}
				ptr_f=ptr_c->face[it];      
				// ok, I have the 4-vertex-face now I load the 5 nodes
				for(size_t j=0;j<4;j++)  _element[i-1].node[j]=ptr_f->vertex[j]->label;//_element[i-1].node[j]=ptr_c->vertex[j]->label;
				// check the order
				double norm1,norm2,norm3;
				double cos12,cos13;
				norm1=Norm(ptr_f->vertex[1]->coord-ptr_f->vertex[0]->coord);
				norm2=Norm(ptr_f->vertex[2]->coord-ptr_f->vertex[0]->coord);        
				norm3=Norm(ptr_f->vertex[3]->coord-ptr_f->vertex[0]->coord);
				cos12=Dot(ptr_f->vertex[1]->coord-ptr_f->vertex[0]->coord,ptr_f->vertex[2]->coord-ptr_f->vertex[0]->coord)/norm1/norm2;
				cos13=Dot(ptr_f->vertex[1]->coord-ptr_f->vertex[0]->coord,ptr_f->vertex[3]->coord-ptr_f->vertex[0]->coord)/norm1/norm3;
				if(cos12<cos13) {_element[i-1].node[2]=ptr_f->vertex[3]->label;_element[i-1].node[3]=ptr_f->vertex[2]->label;}
				// the "positive" vertex is the top of the pyramid
				for(size_t k=0;k<5;k++)
				{
					bool in_table=false;                
					for(size_t j=0;j<4;j++)
					{
						if(_element[i-1].node[j]==ptr_c->vertex[k]->label){in_table=true;}
					}
					if(!in_table){_element[i-1].node[4]=ptr_c->vertex[k]->label;break;} 
				}        
				break;
			case 6:
				_element[i-1].type_element=6; // prism
				it=0;
				for((ptr_f=ptr_c->beginFace());(ptr_f=ptr_c->nextFace());it++)
				{if(ptr_f->nb_vertex==3) break;}
				for(size_t j=0;j<3;j++)  
				{_element[i-1].node[j]=ptr_f->vertex[j]->label;_element[i-1].node[j+3]=0;}          
				// at that stage the basement face is ordered
				// we see the complementary point
				for(size_t j=0;j<3;j++)
				{
					it=_element[i-1].node[j]; // take a node of the reference base
					for((ptr_f=ptr_c->beginFace());(ptr_f=ptr_c->nextFace());)
						for((ptr_e=ptr_f->beginEdge());(ptr_e=ptr_f->nextEdge());)
						{
							if((ptr_e->firstVertex->label!=it) && (ptr_e->secondVertex->label!=it)) continue;
							if(ptr_e->firstVertex->label==it) noux=ptr_e->secondVertex->label; else noux=ptr_e->firstVertex->label;
							bool still_exist=false;
							for(size_t k=0;k<3;k++)
								if(noux==_element[i-1].node[k]) still_exist=true;
							if(!still_exist) {_element[i-1].node[j+3]=noux;}    
						}
				}
				// check if we have really 6 nodes             
				for(size_t j=0;j<6;j++) 
				{
					if(! _element[i-1].node[j])
					{
						cout<<" Not a prism element: abort convertion"<<endl; return;
					}
				}    
				break;    
			case 8:
				_element[i-1].type_element=5; // cube
				ptr_f=ptr_c->face[0]; // reference face
				for(size_t j=0;j<4;j++)  _element[i-1].node[j]=ptr_f->vertex[j]->label;
				// we have to sort the four first points.
				t1= _node[_element[i-1].node[0]-1].coord-ptr_f->centroid;t1/=Norm(t1);// first tangente
				t2=CrossProduct(t1,ptr_f->LocalNormal[0]);t2/=Norm(t2);// second tangente   
				for(size_t j=0;j<4;j++)
				{
					Paux=_node[_element[i-1].node[j]-1].coord-ptr_f->centroid;
					Paux/=Norm(Paux);
					order[j]=(1-Paux*t1)*0.5;if(Paux*t2<0) order[j]*=-1;  
				}  
				
				for(size_t j=3;j>0;j--)
				{
					for(size_t k=0;k<j;k++)
					{
						if(order[k+1]<order[k])
						{
							doux=order[k+1];noux=_element[i-1].node[k+1];
							order[k+1]=order[k];_element[i-1].node[k+1]=_element[i-1].node[k];
							order[k]=doux;_element[i-1].node[k]=noux;
						}
					}
				}           
				// at that stage the basement face is ordered
				// we see the complementary point
				for(size_t j=0;j<4;j++)
				{
					it=_element[i-1].node[j]; // take a node of the reference base
					for((ptr_f=ptr_c->beginFace());(ptr_f=ptr_c->nextFace());)
						for((ptr_e=ptr_f->beginEdge());(ptr_e=ptr_f->nextEdge());)
						{
							if((ptr_e->firstVertex->label!=it) && (ptr_e->secondVertex->label!=it)) continue;
							if(ptr_e->firstVertex->label==it) noux=ptr_e->secondVertex->label; else noux=ptr_e->firstVertex->label;
							bool still_exist=false;
							for(size_t k=0;k<4;k++)
								if(noux==_element[i-1].node[k]) still_exist=true;
							if(!still_exist) _element[i-1].node[j+4]=noux;    
						}
				}
				break;           
			default:
				cout<<"cell does not correspond to a gmsh element. Abort convertion!"<<endl;
				_nb_node=0;_nb_element=0;
				return;
        } 
    }      
}











void Gmsh::FVMesh2Gmsh(FVMesh3D &m,size_t level) // constructor which convert a FVMesh3D into a gmah struct
{
	if(!level) { Gmsh::FVMesh2Gmsh(m); return;}   
	_dim=3;   // a three dimension mesh
	FVVertex3D* ptr_v;
	FVEdge3D* ptr_e;
	FVFace3D* ptr_f;
	FVCell3D *ptr_c;
	FVPoint3D<double> Paux,t1,t2;
	double order[5],doux;
	size_t noux,shift_vertex,shift_edge,shift_face;
	size_t it;
	_nb_node=m.getNbVertex();
	_node.resize(_nb_node); 
	
	for((ptr_v=m.beginVertex());(ptr_v=m.nextVertex());)
    {
		size_t i=ptr_v->label;
		_node[i-1].label  = i; 
		_node[i-1].coord.x= ptr_v->coord.x;
		_node[i-1].coord.y= ptr_v->coord.y;     
		_node[i-1].coord.z= ptr_v->coord.z;        
    }
	_nb_element=m.getNbCell()+m.getNbVertex()+m.getNbEdge()+m.getNbFace();
	_element.resize(_nb_element);
	//------- add the vertex -------
	shift_vertex=0; 
	for((ptr_v=m.beginVertex());(ptr_v=m.nextVertex());)
    { 
		size_t i=ptr_v->label;shift_vertex++;
		_element[i-1].label=i;
		_element[i-1].code_physical=ptr_v->code;
		_element[i-1].code_elementary=ptr_v->code;
		_element[i-1].dim=_dim;  
		_element[i-1].nb_node=1;
		_element[i-1].type_element=15;
		_element[i-1].node[0]=ptr_v->label;
    } 
	// ------- add the edge --------
	shift_edge=0;
	for((ptr_e=m.beginEdge());(ptr_e=m.nextEdge());)
    {        
		size_t i=ptr_e->label+shift_vertex;shift_edge++;
		_element[i-1].label=i+shift_vertex;
		_element[i-1].code_physical=ptr_e->code;
		_element[i-1].code_elementary=ptr_e->code;
		_element[i-1].dim=_dim;  
		_element[i-1].nb_node=2;
		_element[i-1].type_element=1;
		_element[i-1].node[0]=ptr_e->firstVertex->label;
		_element[i-1].node[1]=ptr_e->secondVertex->label;
    } 
	// -------- add the face -------
	shift_face=0;
	for((ptr_f=m.beginFace());(ptr_f=m.nextFace());)
    {   
		size_t i=ptr_f->label+shift_vertex+shift_edge;shift_face++;
		_element[i-1].label=i+shift_vertex+shift_edge;
		_element[i-1].code_physical=ptr_f->code;
		_element[i-1].code_elementary=ptr_f->code;
		_element[i-1].dim=_dim;  
		_element[i-1].nb_node=ptr_f->nb_vertex;
		switch(_element[i-1].nb_node)
        {
			case 3:
				_element[i-1].type_element=2; // triangle
				for(size_t j=0;j<3;j++)  
					_element[i-1].node[j]=ptr_f->vertex[j]->label;
				break;
			case 4:
				_element[i-1].type_element=3; // triangle
				for(size_t j=0;j<4;j++)  
					_element[i-1].node[j]=ptr_f->vertex[j]->label;
				break;
			default:
				cout<<"face "<<ptr_f->label<<" does not correspond to a gmsh element. Abort convertion"<<endl;
				cout<<"find "<<  _element[i-1].nb_node<<" nodes, "; 
				_nb_node=0;_nb_element=0;
				return;
        }
    } 
	// -------- add the cell ------
	for((ptr_c=m.beginCell());(ptr_c=m.nextCell());)
    {
		size_t i=ptr_c->label+shift_vertex+shift_edge+shift_face;
		_element[i-1].label=i+shift_vertex+shift_edge+shift_face;
		_element[i-1].code_physical=ptr_c->code;
		_element[i-1].code_elementary=ptr_c->code;
		_element[i-1].dim=_dim;  
		_element[i-1].nb_node=ptr_c->nb_vertex;
		//cout<<"element "<<i<<" with "<<ptr_c->nb_vertex<<" nodes"<<endl;fflush(NULL);
		switch(_element[i-1].nb_node)
        {
			case 4:   
				_element[i-1].type_element=4; // we use tetrahedron
				for(size_t j=0;j<4;j++)  // the order is not important
					_element[i-1].node[j]=ptr_c->vertex[j]->label;
				break;
			case 5:  
				_element[i-1].type_element=7; // pyramid
				it=0;
				for((ptr_f=ptr_c->beginFace());(ptr_f=ptr_c->nextFace());it++)
				{  
					if( ptr_c->face[it]->nb_vertex==4) break; 
				}
				if (it>=ptr_c->nb_face)  {cout<<" not a pyramidal element"<<endl;return;}
				ptr_f=ptr_c->face[it];      
				// ok, I have the 4-vertex-face now I load the 5 nodes
				for(size_t j=0;j<4;j++)  _element[i-1].node[j]=ptr_f->vertex[j]->label;//_element[i-1].node[j]=ptr_c->vertex[j]->label;
				// check the order
				double norm1,norm2,norm3;
				double cos12,cos13;
				norm1=Norm(ptr_f->vertex[1]->coord-ptr_f->vertex[0]->coord);
				norm2=Norm(ptr_f->vertex[2]->coord-ptr_f->vertex[0]->coord);        
				norm3=Norm(ptr_f->vertex[3]->coord-ptr_f->vertex[0]->coord);
				cos12=Dot(ptr_f->vertex[1]->coord-ptr_f->vertex[0]->coord,ptr_f->vertex[2]->coord-ptr_f->vertex[0]->coord)/norm1/norm2;
				cos13=Dot(ptr_f->vertex[1]->coord-ptr_f->vertex[0]->coord,ptr_f->vertex[3]->coord-ptr_f->vertex[0]->coord)/norm1/norm3;
				if(cos12<cos13) {_element[i-1].node[2]=ptr_f->vertex[3]->label;_element[i-1].node[3]=ptr_f->vertex[2]->label;}
				// the "positive" vertex is the top of the pyramid
				for(size_t k=0;k<5;k++)
				{
					bool in_table=false;                
					for(size_t j=0;j<4;j++)
					{
						if(_element[i-1].node[j]==ptr_c->vertex[k]->label){in_table=true;}
					}
					if(!in_table){_element[i-1].node[4]=ptr_c->vertex[k]->label;break;} 
				}               
				break;
			case 6:
				_element[i-1].type_element=6; // prism
				it=0;
				for((ptr_f=ptr_c->beginFace());(ptr_f=ptr_c->nextFace());it++)
				{if(ptr_f->nb_vertex==3) break;}
				for(size_t j=0;j<3;j++)  
				{_element[i-1].node[j]=ptr_f->vertex[j]->label;_element[i-1].node[j+3]=0;}          
				// at that stage the basement face is ordered
				// we see the complementary point
				for(size_t j=0;j<3;j++)
				{
					it=_element[i-1].node[j]; // take a node of the reference base
					for((ptr_f=ptr_c->beginFace());(ptr_f=ptr_c->nextFace());)
						for((ptr_e=ptr_f->beginEdge());(ptr_e=ptr_f->nextEdge());)
						{
							if((ptr_e->firstVertex->label!=it) && (ptr_e->secondVertex->label!=it)) continue;
							if(ptr_e->firstVertex->label==it) noux=ptr_e->secondVertex->label; else noux=ptr_e->firstVertex->label;
							bool still_exist=false;
							for(size_t k=0;k<3;k++)
								if(noux==_element[i-1].node[k]) still_exist=true;
							if(!still_exist) {_element[i-1].node[j+3]=noux;}    
						}
				}
				// check if we have really 6 nodes             
				for(size_t j=0;j<6;j++) 
				{
					if(! _element[i-1].node[j])
					{
						cout<<" Not a prism element: abort convertion"<<endl; return;
					}
				}    
				break;    
			case 8:
				_element[i-1].type_element=5; // cube
				ptr_f=ptr_c->face[0]; // reference face
				for(size_t j=0;j<4;j++)  _element[i-1].node[j]=ptr_f->vertex[j]->label;
				// we have to sort the four first points.
				t1= _node[_element[i-1].node[0]-1].coord-ptr_f->centroid;t1/=Norm(t1);// first tangente
				t2=CrossProduct(t1,ptr_f->LocalNormal[0]);t2/=Norm(t2);// second tangente   
				for(size_t j=0;j<4;j++)
				{
					Paux=_node[_element[i-1].node[j]-1].coord-ptr_f->centroid;
					Paux/=Norm(Paux);
					order[j]=(1-Paux*t1)*0.5;if(Paux*t2<0) order[j]*=-1;  
				}  
				
				for(size_t j=3;j>0;j--)
				{
					for(size_t k=0;k<j;k++)
					{
						if(order[k+1]<order[k])
						{
							doux=order[k+1];noux=_element[i-1].node[k+1];
							order[k+1]=order[k];_element[i-1].node[k+1]=_element[i-1].node[k];
							order[k]=doux;_element[i-1].node[k]=noux;
						}
					}
				}           
				// at that stage the basement face is ordered
				// we see the complementary point
				for(size_t j=0;j<4;j++)
				{
					it=_element[i-1].node[j]; // take a node of the reference base
					for((ptr_f=ptr_c->beginFace());(ptr_f=ptr_c->nextFace());)
						for((ptr_e=ptr_f->beginEdge());(ptr_e=ptr_f->nextEdge());)
						{
							if((ptr_e->firstVertex->label!=it) && (ptr_e->secondVertex->label!=it)) continue;
							if(ptr_e->firstVertex->label==it) noux=ptr_e->secondVertex->label; else noux=ptr_e->firstVertex->label;
							bool still_exist=false;
							for(size_t k=0;k<4;k++)
								if(noux==_element[i-1].node[k]) still_exist=true;
							if(!still_exist) _element[i-1].node[j+4]=noux;    
						}
				}
				break;           
			default:
				cout<<"find "<< _element[i-1].nb_node<<" nodes,  ";   
				cout<<"3D cell "<<ptr_c->label<<" does not correspond to a gmsh element. Abort convertion"<<endl;
				_nb_node=0;_nb_element=0;
				return;
        } 
    }
}
void Gmsh::Comsol2Gmsh( Comsol &com) // convert comsol into gmsh
{
	_node.resize(_nb_node=com.getNbNode());
	_element.resize(_nb_element=com.getNbElement());
	//cout<<"nb node"<<_nb_node<<" nb_element"<<_nb_element<<endl;
	_dim=com.getDim();
	_name=com.getName();
	for(size_t i=0;i<_nb_node;i++)  // copy the nodes
    {
		_node[i]=(*com.getNode(i));  
    }
	for(size_t i=0;i<_nb_element;i++)  // copy the elements
    {
		CSElement CSel;
		CSel=(*com.getElement(i)); 
		_element[i].label=CSel.label;
		_element[i].type_element=CSel.type_element;
		_element[i].code_physical=CSel.code_physical; 
		_element[i].code_elementary=CSel.code_elementary;
		_element[i].nb_node=CSel.nb_node;
		_element[i].dim=CSel.dim;
		//cout<<"nb node per element="<<_element[i].nb_node<<endl;
		for(size_t j=0;j<CSel.nb_node;j++)
			_element[i].node[j]=CSel.node[j];
    }    
}