// FVMesh3D.cpp
// Last Modification: May 2014 (R. COSTA)
#include <iostream>
#include <cstdio>
#include <math.h>
#include "FVMesh3D.h"
#include "Gmsh.h"
#include "FVTools.h"
#include "FVSparseV.h"
// defautl constructor
FVMesh3D::FVMesh3D()
{
    _nb_vertex=0;_nb_cell=0;_nb_edge=0,_nb_face=0;_dim=0;_nb_boundary_face=0;
}
// constructor with a xml file
FVMesh3D::FVMesh3D(const char *filename)
{
    FVMesh3D::read(filename);
}
// load the mesh from a xml file
size_t FVMesh3D::read(const char *filename)
{
    _spxml.readXML(filename,_xml); // read the file and put in the XML string
    string key, value,element;
    StringMap attribute; 
    //StringMap::iterator iter;
    size_t code;
	// open  balise FVLIB 
    if (_spxml.openBalise("FVLIB")!=OkOpenBalise)
	{cout<<" No open VFLIB balise found:"<<"  in "<<__FILE__<< " line "<<__LINE__<<endl;return(NoOpenBalise);}    
	// find the FIELD  balise
    code=_spxml.openBalise("MESH");
    if (code!=OkOpenBalise)
	{cout<<" No open MESH balise found:"<<"  in "<<__FILE__<< " line "<<__LINE__<<endl; return(NoOpenBalise);} 
    attribute=_spxml.getAttribute();           
    // read the name of the mesh  
    key=string("dim");
    if(key.empty()) 
	{cout<<" No dim attribute found :"<<"  in "<<__FILE__<< " line "<<__LINE__<<endl; return(NoAttribute);}     
    value= attribute[key];    
    _dim= (unsigned) atoi( value.c_str()); 
    if(_dim!=3)
	{
#ifdef _DEBUGS   
		cout<<" dimension do not match:"<<"  in "<<__FILE__<< " line "<<__LINE__<<endl;
#endif
		return(FVWRONGDIM);
	} 
    key=string("name");
    if(key.empty()) 
	{cout<<" No name attribute found :"<<"  in "<<__FILE__<< " line "<<__LINE__<<endl; return(NoAttribute);} 
    _name= attribute[key];
	// ------------------------ we read the vertices  ------------------
    code=_spxml.openBalise("VERTEX");    
    if (code!=OkOpenBalise)
	{cout<<" No open VERTEX balise found:"<<"  in "<<__FILE__<< " line "<<__LINE__<<endl; return(NoOpenBalise);}  
    attribute=_spxml.getAttribute();     
	// read the number of vertices
	key=string("nbvertex");
	if(key.empty()) 
	{cout<<" No nbvertex attribute found :"<<"  in "<<__FILE__<< " line "<<__LINE__<<endl; return(NoAttribute);}     
	value= attribute[key];
	_nb_vertex=(unsigned) atoi( value.c_str()); 
	// read the vertices data
    _vertex.resize(_nb_vertex);
    _spxml.data();
    {
        size_t beginDATA=_spxml.getPosition();
        size_t lengthDATA=_spxml.getLength();
        char  *thedata,*ptr;
        thedata=(char *) malloc(lengthDATA+1);
        _xml.copy(thedata,lengthDATA,beginDATA);
        thedata[lengthDATA]=0;
        ptr=thedata;
        size_t count=0;
        while(count<_nb_vertex)// read the data and put it in the valarray
		{
            _vertex[count].label= strtod(ptr, &ptr);
            _vertex[count].code= strtod(ptr, &ptr);
            _vertex[count].coord.x= strtod(ptr, &ptr);
            _vertex[count].coord.y= strtod(ptr, &ptr);  
            _vertex[count].coord.z= strtod(ptr, &ptr);          
			count++; 
		}
        free(thedata);
	}      
	// close Balise VERTEX
    code=_spxml.closeBalise("VERTEX");
    if (code!=OkCloseBalise)
	{cout<<" No close VERTEX balise found:"<<"  in "<<__FILE__<< " line "<<__LINE__<<endl; return(NoCloseBalise);}   
	//----------------  we read the edges --------------------------        
	code=_spxml.openBalise("EDGE");    
    if (code!=OkOpenBalise)
	{cout<<" No open EDGE balise found:"<<"  in "<<__FILE__<< " line "<<__LINE__<<endl; return(NoOpenBalise);}    
    attribute=_spxml.getAttribute();     
	// read the number of cell
	key=string("nbedge");
	if(key.empty()) 
	{cout<<" No nbedge attribute found :"<<"  in "<<__FILE__<< " line "<<__LINE__<<endl; return(NoAttribute);}     
	value= attribute[key];
	_nb_edge=(unsigned) atoi( value.c_str());  
	// read the EDGE data
    _edge.resize(_nb_edge);
    _spxml.data();
    {
        size_t beginDATA=_spxml.getPosition();
        size_t lengthDATA=_spxml.getLength();
        char  *thedata,*ptr;
        thedata=(char *) malloc(lengthDATA+1);
        _xml.copy(thedata,lengthDATA,beginDATA);
        thedata[lengthDATA]=0;
        ptr=thedata;
        size_t count=0,_label;
        while(count<_nb_edge)// read the data and put it in the valarray
		{
            _label= strtod(ptr, &ptr);
            _edge[_label-1].label=_label;
            _edge[_label-1].code= strtod(ptr, &ptr); 
            _edge[_label-1].nb_vertex= strtod(ptr, &ptr);   // should be always 2
            _edge[_label-1].firstVertex= &(_vertex[strtod(ptr, &ptr)-1]); 
            _edge[_label-1].secondVertex=&(_vertex[strtod(ptr, &ptr)-1]);     
			count++; 
		}
        free(thedata) ;   
	}           
    // close  Balise   EDGE
    code=_spxml.closeBalise("EDGE");
    if (code!=OkCloseBalise)
	{cout<<" No close EDGE balise found:"<<"  in "<<__FILE__<< " line "<<__LINE__<<endl; return(NoCloseBalise);}         
	//----------------  we read the faces --------------------------        
	code=_spxml.openBalise("FACE");    
    if (code!=OkOpenBalise)
	{cout<<" No open FACE balise found:"<<"  in "<<__FILE__<< " line "<<__LINE__<<endl; return(NoOpenBalise);}    
    attribute=_spxml.getAttribute();     
	// read the number of cell
	key=string("nbface");
	if(key.empty()) 
	{cout<<" No nbface attribute found :"<<"  in "<<__FILE__<< " line "<<__LINE__<<endl; return(NoAttribute);}     
	value= attribute[key];
	_nb_face=(unsigned) atoi( value.c_str());  
	// read the EDGE data
    _face.resize(_nb_face);
    _spxml.data();
    {
        size_t beginDATA=_spxml.getPosition();
        size_t lengthDATA=_spxml.getLength();
        char  *thedata,*ptr;
        thedata=(char *) malloc(lengthDATA+1);
        _xml.copy(thedata,lengthDATA,beginDATA);
        thedata[lengthDATA]=0;
        ptr=thedata;
        size_t count=0,_label;
        while(count<_nb_face)// read the data and put it in the valarray
		{
            _label= strtod(ptr, &ptr);
            _face[_label-1].label=_label;
            _face[_label-1].code= strtod(ptr, &ptr);
            _face[_label-1].nb_edge= strtod(ptr, &ptr);  
            for(size_t i=0;i<_face[_label-1].nb_edge;i++)
				_face[_label-1].edge[i]= &(_edge[strtod(ptr, &ptr)-1]); 
            count++; 
		}
		free(thedata);
	}         
    // close  Balise   FACE
    code=_spxml.closeBalise("FACE");
    if (code!=OkCloseBalise)
	{cout<<" No close FACEbalise found:"<<"  in "<<__FILE__<< " line "<<__LINE__<<endl; return(NoCloseBalise);}  
	//----------------  we read the cells --------------------------
    code=_spxml.openBalise("CELL");    
    if (code!=OkOpenBalise)
	{cout<<" No open CELL balise found:"<<"  in "<<__FILE__<< " line "<<__LINE__<<endl; return(NoOpenBalise);}    
    attribute=_spxml.getAttribute();     
	// read the number of cell
	key=string("nbcell");
	if(key.empty()) 
	{cout<<" No nbcell attribute found :"<<"  in "<<__FILE__<< " line "<<__LINE__<<endl; return(NoAttribute);}     
	value= attribute[key];
	_nb_cell=(unsigned) atoi( value.c_str());  
	// read the CELL data
    _cell.resize(_nb_cell);
    _spxml.data();
    {
        size_t beginDATA=_spxml.getPosition();
        size_t lengthDATA=_spxml.getLength();
        char  *thedata,*ptr;
        thedata=(char *) malloc(lengthDATA+1);
        _xml.copy(thedata,lengthDATA,beginDATA);
        thedata[lengthDATA]=0;
        ptr=thedata;
        size_t count=0,_label;
        while(count<_nb_cell)// read the data and put it in the valarray
		{
            _label= strtod(ptr, &ptr);
            _cell[_label-1].label=_label;
            _cell[_label-1].code= strtod(ptr, &ptr);
            _cell[_label-1].nb_face= strtod(ptr, &ptr);  
            for(size_t i=0;i<_cell[_label-1].nb_face;i++)
				_cell[_label-1].face[i]= &(_face[strtod(ptr, &ptr)-1]);        
            count++; 
		}
		free(thedata);      
	}
    // close Balise CELL
    code=_spxml.closeBalise("CELL");
    if (code!=OkCloseBalise)
	{cout<<" No close CELL balise found:"<<"  in "<<__FILE__<< " line "<<__LINE__<<endl; return(NoCloseBalise);}          
	//---------------------------------------------                
    // close Balise MESH
    code=_spxml.closeBalise("MESH");
    if (code!=OkCloseBalise)
	{cout<<" No close MESH balise found:"<<"  in "<<__FILE__<< " line "<<__LINE__<<endl; return(NoCloseBalise);} 
    // close Balise FVLIB
    code=_spxml.closeBalise("FVLIB");
    if (code!=OkCloseBalise)
	{cout<<" No close FVLIB balise found:"<<"  in "<<__FILE__<< " line "<<__LINE__<<endl; return(NoCloseBalise);}        
    //------- we have store all the data, now we complete the vertices and cells information 
	FVMesh3D::complete_data();
	return(FVOK);
}    
// compuyte the geometrical data associated to the mesh
void FVMesh3D::complete_data()
{
	// initialize the list of cell for each vertex    
	for(size_t i=0;i<_nb_vertex;i++)
    {
		_vertex[i].nb_cell=0;
    }  
	// we have to determine the list of neighbor cell for each vertex further
	// EDGE
	// compute the centroid and length of edge
	for(size_t i=0;i<_nb_edge;i++)
    {
		_edge[i].centroid=(_edge[i].firstVertex->coord+_edge[i].secondVertex->coord)*0.5;
		_edge[i].mass_center=_edge[i].centroid;
		FVPoint3D<double> u;
		u=_edge[i].firstVertex->coord-_edge[i].secondVertex->coord;
		_edge[i].length=Norm(u);
    }
	// We have completly finish with the edge
	// FACE
	// compute geometric stuff for the face
	for(size_t i=0;i<_nb_face;i++)
    {
		// conpute perimeter of the face 
		_face[i].perimeter=0.0;
		for(size_t j=0;j<_face[i].nb_edge;j++)
        {
			_face[i].perimeter+=_face[i].edge[j]->length; 
        }
		// compute the area of the face
		FVPoint3D<double> P,P1,P2,Paux;
		FVPoint3D<double> GPP;
		FVGaussPoint2D G2D;
		// build the list of vertex pointer for the face  
		pos_v=0;
		bool _still_exist;
		for(size_t j=0;j<_face[i].nb_edge;j++)
        {
			_still_exist=false; 
			for(size_t k=0;k<pos_v;k++)
				if(_face[i].edge[j]->firstVertex==_face[i].vertex[k]) _still_exist=true;
			if(!_still_exist)
			{
				_face[i].vertex[pos_v]=_face[i].edge[j]->firstVertex;
				pos_v++;
			}
			_still_exist=false;  
			for(size_t k=0;k<pos_v;k++)
				if(_face[i].edge[j]->secondVertex==_face[i].vertex[k]) _still_exist=true;
			if(!_still_exist)
			{
				_face[i].vertex[pos_v]=_face[i].edge[j]->secondVertex;
				pos_v++;
			}    
        }
		_face[i].nb_vertex=pos_v;         
		// compute the centroid of the face
		_face[i].centroid=0.0;
		for(size_t k=0;k<_face[i].nb_vertex;k++)
			_face[i].centroid+=_face[i].vertex[k]->coord;    
		_face[i].centroid/=_face[i].nb_vertex; 
		// face area and mass center
		_face[i].mass_center=0.0;
		_face[i].area=0.0;
		FVPoint3D<double> u,v,w;
		for(size_t j=0;j<_face[i].nb_edge;j++)
        {
			P1=_face[i].edge[j]->firstVertex->coord;
			P2=_face[i].edge[j]->secondVertex->coord;
			u=P1-_face[i].centroid; 
			v=P2-_face[i].centroid;
			w=CrossProduct(u,v);
			_face[i].LocalArea[j]=Norm(w)*0.5;
			_face[i].area+=_face[i].LocalArea[j];
			Paux=0.0;
			for (size_t k=1;k<=G2D.getNbPoint(5);k++)
            { 
				GPP=G2D.getPoint(5,k);
				P=GPP.x*P1+GPP.y*P2+GPP.z*_face[i].centroid;
				Paux+=G2D.getWeight(5,k)*P;
            }        
			_face[i].mass_center+=Paux*_face[i].LocalArea[j]; 
        }  
		_face[i].mass_center/=_face[i].area;   
    }
    //  left and right cell, normal vector will be determined after the loop on cells
	// end loop on the faces
	size_t pos;
	for(size_t i=0;i<_nb_cell;i++)
    {
		// compute surface of the cell 
		// determine the left and right cell for the face
		_cell[i].surface=0.0;
		for(size_t j=0;j<_cell[i].nb_face;j++)
        {
			_cell[i].surface+=_cell[i].face[j]->area; 
			pos=_cell[i].face[j]->label-1;
			if(!(_face[pos].leftCell))  
				_face[pos].leftCell=&(_cell[i]);
			else
				_face[pos].rightCell=&(_cell[i]); 
        }
		// compute the cell2face vector
		// build the list of the vertex pointer for a cell     
		pos_v=0;
		for(size_t j=0;j<_cell[i].nb_face;j++)
			for(size_t k=0;k<_cell[i].face[j]->nb_edge;k++)
            {
				bool _still_exist;   
				_still_exist=false;  
				for(size_t m=0;m<pos_v;m++)
					if(_cell[i].face[j]->edge[k]->firstVertex==_cell[i].vertex[m]) _still_exist=true;
				if(!_still_exist)
				{
					_cell[i].vertex[pos_v]=_cell[i].face[j]->edge[k]->firstVertex;
					pos_v++;
				}
				_still_exist=false;  
				for(size_t m=0;m<pos_v;m++)
					if(_cell[i].face[j]->edge[k]->secondVertex==_cell[i].vertex[m]) _still_exist=true;
				if(!_still_exist)
				{
					_cell[i].vertex[pos_v]=_cell[i].face[j]->edge[k]->secondVertex;
					pos_v++;
				}  
            }
		_cell[i].nb_vertex=pos_v; 
		// build the list of the cell pointer for a vertex     
		// compute the centroid of the cell
		_cell[i].centroid=0.0;
		for(size_t j=0;j<_cell[i].nb_vertex;j++)
        {
			pos=_cell[i].vertex[j]->label-1;
			_vertex[pos].cell[_vertex[pos].nb_cell]=&(_cell[i]); 
			_vertex[pos].nb_cell++;
			_cell[i].centroid+=_cell[i].vertex[j]->coord;
			if(_vertex[pos].nb_cell>=NB_CELL_PER_VERTEX_3D)
			{
				cout<<"Error: overflow in class FVVertex3D, too many Cells, found "<<_vertex[pos].nb_cell<<" on vertex "<<pos+1<<". Line: "<<__LINE__<<". File: "<<__FILE__<<endl;
				exit(0);
			}
        }
		_cell[i].centroid/=_cell[i].nb_vertex; 
		// compute the volume of the cell
		_cell[i].volume=0.;
		_cell[i].mass_center=0.;    
		FVPoint3D<double> P,P1,P2,Paux;
		FVPoint4D<double> GPPP;
		FVGaussPoint3D G3D; 
		double Vaux;
		FVPoint3D<double> u,v,w;
		for(size_t j=0;j<_cell[i].nb_face;j++)
        {    
			_cell[i].cell2face[j]= _cell[i].face[j]->centroid-_cell[i].centroid;    
			for(size_t k=0;k<_cell[i].face[j]->nb_edge;k++)
            {
				P1=_cell[i].face[j]->edge[k]->firstVertex->coord;
				P2=_cell[i].face[j]->edge[k]->secondVertex->coord;
				u=_cell[i].cell2face[j]; 
				v=P1-_cell[i].centroid; 
				w=P2-_cell[i].centroid; 
				Vaux=fabs(Det(u,v,w))/6; 
				_cell[i].volume+=Vaux;
				Paux=0;
				for(size_t r=1;r<=G3D.getNbPoint(5);r++)
                {              
					GPPP=G3D.getPoint(5,r);
					P=GPPP.x*P1+GPPP.y*P2+GPPP.z*_cell[i].face[j]->centroid+GPPP.t*_cell[i].centroid;
					Paux+=G3D.getWeight(5,r)*P;
                }
				_cell[i].mass_center+=Paux*Vaux;    
            }           
        }
		_cell[i].mass_center/=_cell[i].volume;   
    }
	//  we compute the normal from left to rigth for each sub-triangle   
	_boundary_face.resize(0);
	_nb_boundary_face=0;
	FVPoint3D<double> u,v,w;
	double no;
	for(size_t i=0;i<_nb_face;i++)
    {
		_face[i].normal=0.0;
		double surf=0.0;
		for(size_t j=0;j<_face[i].nb_edge;j++)
        {
			u=_face[i].edge[j]->firstVertex->coord-_face[i].centroid;
			v=_face[i].edge[j]->secondVertex->coord-_face[i].centroid;
			w=CrossProduct(u,v);
			no=Norm(w);
			w/=no;
			_face[i].LocalNormal[j]=w;
			if(!_face[i].leftCell){cout<<" ERROR no left cell for face "<<i<<endl;fflush(NULL);}
			u=_face[i].centroid-_face[i].leftCell->centroid;          
			if(w*u<0) _face[i].LocalNormal[j]*=-1.;
			surf+=no;
			_face[i].normal+=no*_face[i].LocalNormal[j];          
		}     // build the list of boundary face
		if(!(_face[i].rightCell)) {_boundary_face.push_back(&(_face[i]));_nb_boundary_face++;} 
		_face[i].normal/=surf;
    }      
}
// write the mesh into a xml file
size_t FVMesh3D::write(const char *filename)
{
	if((_nb_cell==0) || (_nb_face==0) || (_nb_edge==0) || (_nb_vertex==0))
	{
		cout<<" error in file:"<<filename<<"  in "<<__FILE__<< " line "<<__LINE__<<endl; 
		cout<<" there is no mesh to save"<<endl;
		return(FVERROR);
	}    
	ofstream  out_file;    
	out_file.open(filename);
	if(!out_file.is_open())
	{
		cout<<" error in file:"<<filename<<"  in "<<__FILE__<< " line "<<__LINE__<<endl; 
		cout<<" can not create the file"<<endl;
		return(FVNOFILE);
	}
	out_file<<"<?xml version=\"1.0\" encoding=\"ISO-8859-1\"?>"<<endl;
	out_file<<"<FVLIB>"<<endl; 
	out_file<<"    <MESH dim=\"3\"    name=\""<<_name<<"\">"<<endl;
	out_file<<"         <VERTEX nbvertex=\""<<_nb_vertex<<"\">"<<endl;
	for(size_t i=0;i<_nb_vertex;i++)
    {
		out_file<< setw(FVCHAMPINT)<<_vertex[i].label<< setw(FVCHAMPINT)<< _vertex[i].code;  
		out_file<<scientific << setprecision(FVPRECISION) << setw(FVCHAMP) << _vertex[i].coord.x;
		out_file<<scientific << setprecision(FVPRECISION) << setw(FVCHAMP) << _vertex[i].coord.y;
		out_file<<scientific << setprecision(FVPRECISION) << setw(FVCHAMP) << _vertex[i].coord.z<<endl; 
    }
	out_file<<"         </VERTEX>"<<endl; 
	out_file<<"         <EDGE nbedge=\""<<_nb_edge<<"\">"<<endl;
	for(size_t i=0;i<_nb_edge;i++)
    {
		out_file<< setw(FVCHAMPINT)<<_edge[i].label<< setw(FVCHAMPINT)<< _edge[i].code
		<< setw(FVCHAMPINT)<< _edge[i].nb_vertex; 
		out_file<<setw(FVCHAMPINT) << _edge[i].firstVertex->label<<setw(FVCHAMPINT) << _edge[i].secondVertex->label<<endl; 
    } 
	out_file<<"         </EDGE>"<<endl; 
	out_file<<"         <FACE nbface=\""<<_nb_face<<"\">"<<endl;
	for(size_t i=0;i<_nb_face;i++)
    {
		out_file<< setw(FVCHAMPINT)<<_face[i].label<< setw(FVCHAMPINT)<< _face[i].code
		<< setw(FVCHAMPINT)<< _face[i].nb_edge;
		for(size_t j=0;j<_face[i].nb_edge;j++)        
			out_file<<setw(FVCHAMPINT) << _face[i].edge[j]->label<<setw(FVCHAMPINT); 
		out_file<<endl; 
    } 
	out_file<<"         </FACE>"<<endl; 
	out_file<<"         <CELL nbcell=\""<<_nb_cell<<"\">"<<endl;
	for(size_t i=0;i<_nb_cell;i++)
    {
		out_file<< setw(FVCHAMPINT)<<_cell[i].label<< setw(FVCHAMPINT)<< _cell[i].code
		<< setw(FVCHAMPINT)<< _cell[i].nb_face;
		for(size_t j=0;j<_cell[i].nb_face;j++)        
			out_file<<setw(FVCHAMPINT) << _cell[i].face[j]->label<<setw(FVCHAMPINT); 
		out_file<<endl; 
    } 
	out_file<<"         </CELL>"<<endl;      
	out_file<<"    </MESH>"<<endl;
	out_file<<"</FVLIB>"<<endl;          
	out_file.close();    
	return(FVOK);
}
// convert a Gmsh struct into a FVMesh3D
void FVMesh3D::Gmsh2FVMesh(Gmsh &m) 
{
	if (m.getDim()!=3)
	{
		cout<<" Dimension don't match. The Gmsh mesh should contain only 0 or 1D elements"<<endl;
		return;
	}  
	_nb_vertex=m.getNbNode();
	_name="mesh3D  from gmsh file"; 
	_vertex.resize(_nb_vertex); 
	for(size_t i=0;i<_nb_vertex;i++)
    {
		_vertex[i].label=m.getNode(i)->label;
		_vertex[i].coord.x=m.getNode(i)->coord.x;
		_vertex[i].coord.y=m.getNode(i)->coord.y;
		_vertex[i].coord.z=m.getNode(i)->coord.z;
		_vertex[i].code=0;  // default value
    }     
	_nb_cell=0;_nb_edge=0;_nb_face=0;   
	_cell.resize(0);_face.resize(0);_edge.resize(0);
	_cell.reserve(m.getNbElement());_face.reserve(m.getNbElement()*3);_edge.reserve(m.getNbElement()*5);
	FVCell3D cell;
	FVFace3D face;
	FVEdge3D edge;
	size_t pos_dim3=0;
	bool exist;
	for(size_t i=0;i<m.getNbElement();i++)
    {  
		GMElement *ptr_el=m.getElement(i);  
		if  (ptr_el->type_element==15) 
		{
			pos_dim3++;
			ptr_el->label=_vertex[ptr_el->node[0]-1].label;// label is no longer GMElement label but  give 
			// the label of the vertex in the _vertex table
		}
		// 1D element
		if  (ptr_el->type_element==1)
		{              
			pos_dim3++;
			edge.firstVertex=&(_vertex[ptr_el->node[0]-1]);
			edge.secondVertex=&(_vertex[ptr_el->node[1]-1]);
			edge.nb_vertex=2;
			exist=false;
			for(size_t j=0;j<_edge.size();j++)// look if the edge still exist
				if(isEqual(&_edge[j],&edge)) {exist=true;ptr_el->label=j+1;}
			if(!exist)
			{
				ptr_el->label=edge.label=_edge.size()+1;    // label is no longer GMElement label but  give
				// the label of the edge in the _edge table            
				edge.code=0;  // we see the code later
				_edge.push_back(edge); // it is a real new edge so save it
			}    
		}
		// 2D element         
		if  (ptr_el->type_element==2 ||ptr_el->type_element==3 )  // treat   a triangle or a quadrangle
		{
			pos_dim3++;  
			// first we create the edges
			for(size_t j=0;j<=ptr_el->type_element;j++)
			{
				edge.firstVertex=&(_vertex[ptr_el->node[j]-1]);
				edge.secondVertex=&(_vertex[ptr_el->node[(j+1)%(ptr_el->type_element+1)]-1]);
				edge.nb_vertex=2;
				exist=false;
				for(size_t k=0;k<_edge.size();k++)// look if the edge still exist
					if(isEqual(&_edge[k],&edge))  
					{
						face.edge[j]=&_edge[k];exist=true;                    
					} // if exist we keep the pointer
				if(!exist) // if not exist create it
				{  
					edge.label=_edge.size()+1;             
					edge.code=0;  // we see the code later
					_edge.push_back(edge); // it is a real new edge so save it
					face.edge[j]=&_edge[_edge.size()-1]; // and keep the pointer
				}
			}
			// at that stage we have all the edge 
			// we creat the face if necessary and set the code by dimension reduction
			face.nb_edge=ptr_el->type_element+1;  
			exist=false;
			for(size_t j=0;j<_face.size();j++)// look if the face still exist
				if(isEqual(&_face[j],&face))  {exist=true; ptr_el->label=j+1;}// the face still exist,  
			if(!exist) 
			{
				ptr_el->label=face.label=_face.size()+1;  
				face.code=0;  // we see the code later
				_face.push_back(face); // it is a real new edge so save it
				
			}  
		}    
		//  now we treat the 3D elements         
		// treat  tetrahedron 
		if  (ptr_el->type_element==4) 
		{     
			for(size_t j=0;j<4;j++)  //loop on the face
			{
				FVFace3D face;     // construct a Face of triangle
				for(size_t k=0;k<3;k++)  // loop on the edge
				{
					FVEdge3D edge;  //seek the Edge
					size_t n1,n2,nd[3];
					if(j==0) {nd[0]=1;nd[1]=2;nd[2]=3;}
					if(j==1) {nd[0]=0;nd[1]=2;nd[2]=3;}      
					if(j==2) {nd[0]=0;nd[1]=1;nd[2]=3;}   
					if(j==3) {nd[0]=0;nd[1]=1;nd[2]=2;}                    
					n1=nd[k];n2=nd[(k+1)%3]; 
					edge.firstVertex=&(_vertex[ptr_el->node[n1]-1]);
					edge.secondVertex=&(_vertex[ptr_el->node[n2]-1]);
					edge.nb_vertex=2;             
					exist=false;
					for(size_t m=0;m<_edge.size();m++)// look if the edge still exist
					{
						if(isEqual(&_edge[m],&edge))  
						{ 
							face.edge[k]=&_edge[m];exist=true;                    
						} // if exist we keep the pointer
					}         
					if(!exist) // if not exist create it
					{  
						edge.label=_edge.size()+1;        
						edge.code=0;  // we see the code later
						_edge.push_back(edge); // it is a real new edge so save it
						face.edge[k]=&_edge[_edge.size()-1]; // and keep the pointer
					}                     
				}
				// at that stage, the face has all its edge
				// check if the face still exists
				exist=false;
				face.nb_edge=3;             
				for(size_t m=0;m<_face.size();m++)// look if the face still exist
					if(isEqual(&_face[m],&face))  
					{
						cell.face[j]=&_face[m];exist=true;                    
					} // if exist we keep the pointer
				if(!exist) // if not exist create it
				{  
					face.label=_face.size()+1;        
					face.code=0;  // we see the code later
					_face.push_back(face); // it is a real new face so save it
					cell.face[j]=&_face[_face.size()-1]; // and keep the pointer
				}  
			}    
			// at that stage we have all the face and edge now check the code
			cell.nb_face=4;   
			exist=false;
			for(size_t j=0;j<_cell.size();j++)// look if the cell still exist
			{    
				if(isEqual(&_cell[j],&cell))                 
				{
					exist=true;
					if (ptr_el->code_physical/MINUS_ONE_DIM==3)               
					{_cell[j].setCode2Vertex(ptr_el->code_physical%MINUS_ONE_DIM);}                  
					if (ptr_el->code_physical/MINUS_ONE_DIM==2)               
					{_cell[j].setCode2Edge(ptr_el->code_physical%MINUS_ONE_DIM);}
					if (ptr_el->code_physical/MINUS_ONE_DIM==1)
					{_cell[j].setCode2Face(ptr_el->code_physical%MINUS_ONE_DIM);}  
					if (ptr_el->code_physical/MINUS_ONE_DIM==0)                      
						_cell[j].code=ptr_el->code_physical;
				}
			}      
			if(!exist) 
			{
				//cout<<"the cell (gmsh label)"<<ptr_el->label<< " does not exist"<<endl;  fflush(NULL);
				//cout<<"the edge are "<<cell.edge[0]->label<<" "<<cell.edge[1]->label<<" "<<cell.edge[2]->label<<" "<<endl;
				
				ptr_el->label=cell.label=_cell.size()+1;       
				if (ptr_el->code_physical/MINUS_ONE_DIM==3)               
				{cell.setCode2Vertex(ptr_el->code_physical%MINUS_ONE_DIM);}              
				if (ptr_el->code_physical/MINUS_ONE_DIM==2)               
				{cell.setCode2Edge(ptr_el->code_physical%MINUS_ONE_DIM);}
				if (ptr_el->code_physical/MINUS_ONE_DIM==1)
				{cell.setCode2Face(ptr_el->code_physical%MINUS_ONE_DIM);}   
				if (ptr_el->code_physical/MINUS_ONE_DIM==0)                      
					cell.code=ptr_el->code_physical;
				_cell.push_back(cell); // it is a real new cell so save it
				
			}  
		}
		// treat  hexahedron 
		if(ptr_el->type_element==5) 
		{
			//cout<<"creating hexahedron element"<<endl;fflush(NULL);               
			for(size_t j=0;j<6;j++)  //loop on the face
			{
				FVFace3D face;     // construct a Face of triangle
				for(size_t k=0;k<4;k++)  // loop on the edge
				{
					FVEdge3D edge;  //seek the Edge
					size_t n1,n2,nd[4];
					if(j==0) {nd[0]=0;nd[1]=1;nd[2]=2;nd[3]=3;}
					if(j==1) {nd[0]=4;nd[1]=5;nd[2]=6;nd[3]=7;}
					if(j==2) {nd[0]=0;nd[1]=3;nd[2]=7;nd[3]=4;}
					if(j==3) {nd[0]=1;nd[1]=2;nd[2]=6;nd[3]=5;}                  
					if(j==4) {nd[0]=2;nd[1]=3;nd[2]=7;nd[3]=6;}
					if(j==5) {nd[0]=0;nd[1]=1;nd[2]=5;nd[3]=4;}                  
					n1=nd[k];n2=nd[(k+1)%4];
					edge.firstVertex=&(_vertex[ptr_el->node[n1]-1]);
					edge.secondVertex=&(_vertex[ptr_el->node[n2]-1]);
					edge.nb_vertex=2;             
					exist=false;
					for(size_t m=0;m<_edge.size();m++)// look if the edge still exist
						if(isEqual(&_edge[m],&edge))  
						{
							face.edge[k]=&_edge[m];exist=true;                    
						} // if exist we keep the pointer
					if(!exist) // if not exist create it
					{  
						edge.label=_edge.size()+1;      
						edge.code=0;  // we see the code later
						_edge.push_back(edge); // it is a real new edge so save it
						face.edge[k]=&_edge[_edge.size()-1]; // and keep the pointer
					}
				}
				// at that stage, the face has all its edge
				// check if the face still exists
				exist=false;
				face.nb_edge=4;  
				for(size_t m=0;m<_face.size();m++)// look if the face still exist
					if(isEqual(&_face[m],&face))  
					{
						cell.face[j]=&_face[m];exist=true;                    
					} // if exist we keep the pointer
				if(!exist) // if not exist create it
				{  
					face.label=_face.size()+1;             
					face.code=0;  // we see the code later
					_face.push_back(face); // it is a real new face so save it
					cell.face[j]=&_face[_face.size()-1]; // and keep the pointer
				}  
			}    
			// at that stage we have all the face and edge now chesk the code
			cell.nb_face=6;   
			exist=false;
			for(size_t j=0;j<_cell.size();j++)// look if the cell still exist
			{    
				if(isEqual(&_cell[j],&cell))                 
				{
					exist=true;
					if (ptr_el->code_physical/MINUS_ONE_DIM==3)               
					{_cell[j].setCode2Vertex(ptr_el->code_physical%MINUS_ONE_DIM);}                  
					if (ptr_el->code_physical/MINUS_ONE_DIM==2)               
					{_cell[j].setCode2Edge(ptr_el->code_physical%MINUS_ONE_DIM);}
					if (ptr_el->code_physical/MINUS_ONE_DIM==1)
					{_cell[j].setCode2Face(ptr_el->code_physical%MINUS_ONE_DIM);}  
					if (ptr_el->code_physical/MINUS_ONE_DIM==0)                      
						_cell[j].code=ptr_el->code_physical;
				}
			}      
			if(!exist) 
			{
				ptr_el->label=cell.label=_cell.size()+1;       
				if (ptr_el->code_physical/MINUS_ONE_DIM==3)               
				{cell.setCode2Vertex(ptr_el->code_physical%MINUS_ONE_DIM);}              
				if (ptr_el->code_physical/MINUS_ONE_DIM==2)               
				{cell.setCode2Edge(ptr_el->code_physical%MINUS_ONE_DIM);}
				if (ptr_el->code_physical/MINUS_ONE_DIM==1)
				{cell.setCode2Face(ptr_el->code_physical%MINUS_ONE_DIM);}   
				if (ptr_el->code_physical/MINUS_ONE_DIM==0)                      
					cell.code=ptr_el->code_physical;
				_cell.push_back(cell); // it is a real new cell so save it
			}  
		}    
		// treat  prism 
		if  (ptr_el->type_element==6) 
		{
			for(size_t j=0;j<3;j++)  //loop on the face for the quadrangle
			{
				FVFace3D face;     // construct a Face of triangle
				for(size_t k=0;k<4;k++)  // loop on the edge
				{
					FVEdge3D edge;  //seek the Edge               
					size_t n1,n2,nd[4];
					if(j==0) {nd[0]=0;nd[1]=1;nd[2]=4;nd[3]=3;}
					if(j==1) {nd[0]=0;nd[1]=2;nd[2]=5;nd[3]=3;}
					if(j==2) {nd[0]=1;nd[1]=2;nd[2]=5;nd[3]=4;}                
					n1=nd[k];n2=nd[(k+1)%4];               
					edge.firstVertex=&(_vertex[ptr_el->node[n1]-1]);
					edge.secondVertex=&(_vertex[ptr_el->node[n2]-1]);
					edge.nb_vertex=2;
					exist=false;
					for(size_t m=0;m<_edge.size();m++)// look if the edge still exist
						if(isEqual(&_edge[m],&edge))  
						{
							face.edge[k]=&_edge[m];exist=true;                    
						} // if exist we keep the pointer
					if(!exist) // if not exist create it
					{  
						edge.label=_edge.size()+1;             
						edge.code=0;  // we see the code later
						_edge.push_back(edge); // it is a real new edge so save it
						face.edge[k]=&_edge[_edge.size()-1]; // and keep the pointer
					}
				}
				// at that stage, the face has all its edge
				// check if the face still exists
				exist=false;
				face.nb_edge=4;  
				for(size_t m=0;m<_face.size();m++)// look if the face still exist
					if(isEqual(&_face[m],&face))  
					{
						cell.face[j]=&_face[m];exist=true;                    
					} // if exist we keep the pointer
				if(!exist) // if not exist create it
				{  
					face.label=_face.size()+1;             
					face.code=0;  // we see the code later
					_face.push_back(face); // it is a real new face so save it
					cell.face[j]=&_face[_face.size()-1]; // and keep the pointer
				}  
			}  
			for(size_t j=3;j<5;j++)  //loop on the face for the triangle
			{
				FVFace3D face;     // construct a Face of triangle
				for(size_t k=0;k<3;k++)  // loop on the edge
				{
					FVEdge3D edge;  //seek the Edge
					size_t n1,n2,nd[3];
					if(j==3) {nd[0]=3;nd[1]=4;nd[2]=5;}
					if(j==4) {nd[0]=0;nd[1]=1;nd[2]=2;}            
					n1=nd[k];n2=nd[(k+1)%3];                                      
					edge.firstVertex=&(_vertex[ptr_el->node[n1]-1]);
					edge.secondVertex=&(_vertex[ptr_el->node[n2]-1]);
					edge.nb_vertex=2;
					exist=false;
					for(size_t m=0;m<_edge.size();m++)// look if the edge still exist
						if(isEqual(&_edge[m],&edge))  
						{
							face.edge[k]=&_edge[m];exist=true;                    
						} // if exist we keep the pointer
					if(!exist) // if not exist create it
					{  
						edge.label=_edge.size()+1;             
						edge.code=0;  // we see the code later
						_edge.push_back(edge); // it is a real new edge so save it
						face.edge[k]=&_edge[_edge.size()-1]; // and keep the pointer
					}
				}
				// at that stage, the face has all its edge
				// check if the face still exists
				exist=false;
				face.nb_edge=3;  
				for(size_t m=0;m<_face.size();m++)// look if the face still exist
					if(isEqual(&_face[m],&face))  
					{
						cell.face[j]=&_face[m];exist=true;                    
					} // if exist we keep the pointer
				if(!exist) // if not exist create it
				{  
					face.label=_face.size()+1;             
					face.code=0;  // we see the code later
					_face.push_back(face); // it is a real new face so save it
					cell.face[j]=&_face[_face.size()-1]; // and keep the pointer
				}  
			}
			// at that stage we have all the face and edge now chesk the code
			cell.nb_face=5;   
			exist=false;
			for(size_t j=0;j<_cell.size();j++)// look if the cell still exist
			{    
				if(isEqual(&_cell[j],&cell))                 
				{
					exist=true;
					if (ptr_el->code_physical/MINUS_ONE_DIM==3)               
					{_cell[j].setCode2Vertex(ptr_el->code_physical%MINUS_ONE_DIM);}                  
					if (ptr_el->code_physical/MINUS_ONE_DIM==2)               
					{_cell[j].setCode2Edge(ptr_el->code_physical%MINUS_ONE_DIM);}
					if (ptr_el->code_physical/MINUS_ONE_DIM==1)
					{_cell[j].setCode2Face(ptr_el->code_physical%MINUS_ONE_DIM);}  
					if (ptr_el->code_physical/MINUS_ONE_DIM==0)                      
						_cell[j].code=ptr_el->code_physical;
				}
			}      
			if(!exist) 
			{
				ptr_el->label=cell.label=_cell.size()+1;       
				if (ptr_el->code_physical/MINUS_ONE_DIM==3)               
				{cell.setCode2Vertex(ptr_el->code_physical%MINUS_ONE_DIM);}              
				if (ptr_el->code_physical/MINUS_ONE_DIM==2)               
				{cell.setCode2Edge(ptr_el->code_physical%MINUS_ONE_DIM);}
				if (ptr_el->code_physical/MINUS_ONE_DIM==1)
				{cell.setCode2Face(ptr_el->code_physical%MINUS_ONE_DIM);}   
				if (ptr_el->code_physical/MINUS_ONE_DIM==0)                      
					cell.code=ptr_el->code_physical;
				_cell.push_back(cell); // it is a real new cell so save it
			}                
		}   
		// treat  pyramid 
		if  (ptr_el->type_element==7) 
		{    
			for(size_t j=0;j<1;j++)  //loop on the face for the  quadrangle
			{
				FVFace3D face;     // construct a Face of triangle
				for(size_t k=0;k<4;k++)  // loop on the edge
				{
					FVEdge3D edge;  //seek the Edge               
					size_t n1,n2,nd[4];
					if(j==0) {nd[0]=0;nd[1]=1;nd[2]=2;nd[3]=3;}             
					n1=nd[k];n2=nd[(k+1)%4];               
					edge.firstVertex=&(_vertex[ptr_el->node[n1]-1]);
					edge.secondVertex=&(_vertex[ptr_el->node[n2]-1]);
					edge.nb_vertex=2;
					exist=false;
					for(size_t m=0;m<_edge.size();m++)// look if the edge still exist
						if(isEqual(&_edge[m],&edge))  
						{
							face.edge[k]=&_edge[m];exist=true;                    
						} // if exist we keep the pointer
					if(!exist) // if not exist create it
					{  
						edge.label=_edge.size()+1;             
						edge.code=0;  // we see the code later
						_edge.push_back(edge); // it is a real new edge so save it
						face.edge[k]=&_edge[_edge.size()-1]; // and keep the pointer
					}
				}
				// at that stage, the face has all its edge
				// check if the face still exists
				exist=false;
				face.nb_edge=4;  
				for(size_t m=0;m<_face.size();m++)// look if the face still exist
					if(isEqual(&_face[m],&face))  
					{
						cell.face[j]=&_face[m];exist=true;                    
					} // if exist we keep the pointer
				if(!exist) // if not exist create it
				{  
					face.label=_face.size()+1;             
					face.code=0;  // we see the code later
					_face.push_back(face); // it is a real new face so save it
					cell.face[j]=&_face[_face.size()-1]; // and keep the pointer
				}  
			}  
			for(size_t j=1;j<5;j++)  //loop on the face for the triangle
			{
				FVFace3D face;     // construct a Face of triangle
				for(size_t k=0;k<3;k++)  // loop on the edge
				{
					FVEdge3D edge;  //seek the Edge
					size_t n1,n2,nd[3];
					if(j==1) {nd[0]=4;nd[1]=0;nd[2]=1;}
					if(j==2) {nd[0]=4;nd[1]=1;nd[2]=2;}    
					if(j==3) {nd[0]=4;nd[1]=2;nd[2]=3;}
					if(j==4) {nd[0]=4;nd[1]=3;nd[2]=0;}                    
					n1=nd[k];n2=nd[(k+1)%3];                                      
					edge.firstVertex=&(_vertex[ptr_el->node[n1]-1]);
					edge.secondVertex=&(_vertex[ptr_el->node[n2]-1]);
					edge.nb_vertex=2;
					exist=false;
					for(size_t m=0;m<_edge.size();m++)// look if the edge still exist
						if(isEqual(&_edge[m],&edge))  
						{
							face.edge[k]=&_edge[m];exist=true;                    
						} // if exist we keep the pointer
					if(!exist) // if not exist create it
					{  
						edge.label=_edge.size()+1;             
						edge.code=0;  // we see the code later
						_edge.push_back(edge); // it is a real new edge so save it
						face.edge[k]=&_edge[_edge.size()-1]; // and keep the pointer
					}
				}
				// at that stage, the face has all its edge
				// check if the face still exists
				exist=false;
				face.nb_edge=3;  
				for(size_t m=0;m<_face.size();m++)// look if the face still exist
					if(isEqual(&_face[m],&face))  
					{
						cell.face[j]=&_face[m];exist=true;                    
					} // if exist we keep the pointer
				if(!exist) // if not exist create it
				{  
					face.label=_face.size()+1;             
					face.code=0;  // we see the code later
					_face.push_back(face); // it is a real new face so save it
					cell.face[j]=&_face[_face.size()-1]; // and keep the pointer
				}  
			}               
			// at that stage we have all the face and edge now chesk the code
			cell.nb_face=5;   
			exist=false;
			for(size_t j=0;j<_cell.size();j++)// look if the cell still exist
			{    
				if(isEqual(&_cell[j],&cell))                 
				{
					exist=true;
					if (ptr_el->code_physical/MINUS_ONE_DIM==3)               
					{_cell[j].setCode2Vertex(ptr_el->code_physical%MINUS_ONE_DIM);}                  
					if (ptr_el->code_physical/MINUS_ONE_DIM==2)               
					{_cell[j].setCode2Edge(ptr_el->code_physical%MINUS_ONE_DIM);}
					if (ptr_el->code_physical/MINUS_ONE_DIM==1)
					{_cell[j].setCode2Face(ptr_el->code_physical%MINUS_ONE_DIM);}  
					if (ptr_el->code_physical/MINUS_ONE_DIM==0)                      
						_cell[j].code=ptr_el->code_physical;
				}
			}      
			if(!exist) 
			{
				ptr_el->label=cell.label=_cell.size()+1;       
				if (ptr_el->code_physical/MINUS_ONE_DIM==3)               
				{cell.setCode2Vertex(ptr_el->code_physical%MINUS_ONE_DIM);}              
				if (ptr_el->code_physical/MINUS_ONE_DIM==2)               
				{cell.setCode2Edge(ptr_el->code_physical%MINUS_ONE_DIM);}
				if (ptr_el->code_physical/MINUS_ONE_DIM==1)
				{cell.setCode2Face(ptr_el->code_physical%MINUS_ONE_DIM);}   
				if (ptr_el->code_physical/MINUS_ONE_DIM==0)                      
					cell.code=ptr_el->code_physical;
				_cell.push_back(cell); // it is a real new cell so save it
				
			}   
		}            
	}
	_nb_edge=_edge.size();
	_nb_face=_face.size();
	_nb_cell=_cell.size();  
	// all the edge and cell are created    
	// now we treat the code cascade backward
	for(size_t cont=pos_dim3+1;cont>0;cont--) // we skip of 1 becaus NEGATIVE NUMBER DOES NOT ESXIST with size_t
	{  
		size_t i=cont-1; // astuce to allow deacresing counter   
		GMElement *ptr_el=m.getElement(i); 
		// We treat the 2D element 
		if  (ptr_el->type_element==2 ||ptr_el->type_element==3 )  // treat   a triangle or a quadrangle
        {
			// we set the vertex code 
			if (ptr_el->code_physical/MINUS_ONE_DIM==2)
			{_face[ptr_el->label-1].setCode2Vertex(ptr_el->code_physical%MINUS_ONE_DIM);}  
			// we set the edge code             
			if (ptr_el->code_physical/MINUS_ONE_DIM==1)
			{_face[ptr_el->label-1].setCode2Edge(ptr_el->code_physical%MINUS_ONE_DIM);}  
			// we set the facex code 
			if (ptr_el->code_physical/MINUS_ONE_DIM==0)                      
				_face[ptr_el->label-1].code=ptr_el->code_physical;   
        }
		// we treat the edge
		if(ptr_el->type_element==1)
        { 
			if (ptr_el->code_physical/MINUS_ONE_DIM==1)
			{_edge[ptr_el->label-1].setCode2Vertex(ptr_el->code_physical%MINUS_ONE_DIM);}  // we set the Vertex code 
			if (ptr_el->code_physical/MINUS_ONE_DIM==0)                      
				_edge[ptr_el->label-1].code=ptr_el->code_physical;   // we set the edge code 
        }       
		//we treat the node       
		if(ptr_el->type_element==15) 
		{
            _vertex[ptr_el->label-1].code=ptr_el->code_physical;  // we set the Vertex code 
		}
	}       
	FVMesh3D::complete_data();
}
// load a Gmsh mesh format
void loadMesh(Gmsh &m,vector<A2B> &c2v,vector<A2B> &f2v, vector<A2B> &e2v)
{
	A2B edge2vertex,face2vertex,cell2vertex;    
	for(size_t i=0;i<m.getNbElement();i++)
    {  
		GMElement *ptr_el=m.getElement(i);  
		// Vertex
		if  (ptr_el->type_element==15) 
        {
			; 
        }
		// edge
		if  (ptr_el->type_element==1)
        {     
			edge2vertex.clean();    
			edge2vertex.setLabel(e2v.size());
			edge2vertex.setTypeA(1);    
			edge2vertex.setCode(ptr_el->code_physical);edge2vertex.setLayer(1);
			edge2vertex.addLabel(ptr_el->node[0]);edge2vertex.addLabel(ptr_el->node[1]);  
			e2v.push_back(edge2vertex);
        }    
		// triangle         
		if  (ptr_el->type_element==2)
        {
			face2vertex.clean();
			face2vertex.setLabel(f2v.size());
			face2vertex.setTypeA(2);   
			face2vertex.setLayer(1);face2vertex.setCode(ptr_el->code_physical);
			face2vertex.addLabel(ptr_el->node[0]);face2vertex.addLabel(ptr_el->node[1]);
			face2vertex.addLabel(ptr_el->node[2]);
			f2v.push_back(face2vertex);
        }
		//  quadrilateral    
		if  (ptr_el->type_element==3 )  
        {   
			face2vertex.clean();
			face2vertex.setLabel(f2v.size());
			face2vertex.setTypeA(3);   
			face2vertex.setLayer(1);face2vertex.setCode(ptr_el->code_physical);
			face2vertex.addLabel(ptr_el->node[0]);face2vertex.addLabel(ptr_el->node[1]);
			face2vertex.addLabel(ptr_el->node[2]);face2vertex.addLabel(ptr_el->node[3]);
			f2v.push_back(face2vertex);
        }
		// treat  tetrahedron 
		if  (ptr_el->type_element==4)    
        {              
			cell2vertex.clean();
			cell2vertex.setLabel(c2v.size());
			cell2vertex.setTypeA(4);
			cell2vertex.setLayer(1);cell2vertex.setCode(ptr_el->code_physical);   
			cell2vertex.addLabel(ptr_el->node[0]);cell2vertex.addLabel(ptr_el->node[1]);
			cell2vertex.addLabel(ptr_el->node[2]);cell2vertex.addLabel(ptr_el->node[3]);
			c2v.push_back(cell2vertex);
        }         
		// treat  hexahedron 
		if  (ptr_el->type_element==5) 
        {              
			cell2vertex.clean();
			cell2vertex.setLabel(c2v.size());
			cell2vertex.setTypeA(5);   
			cell2vertex.setLayer(1); cell2vertex.setCode(ptr_el->code_physical);   
			cell2vertex.addLabel(ptr_el->node[0]);cell2vertex.addLabel(ptr_el->node[1]);
			cell2vertex.addLabel(ptr_el->node[2]);cell2vertex.addLabel(ptr_el->node[3]);
			cell2vertex.addLabel(ptr_el->node[4]);cell2vertex.addLabel(ptr_el->node[5]);   
			cell2vertex.addLabel(ptr_el->node[6]);cell2vertex.addLabel(ptr_el->node[7]);     
			c2v.push_back(cell2vertex);
        }         
		// treat  prism 
		if  (ptr_el->type_element==6) 
        {              
			cell2vertex.clean();
			cell2vertex.setLabel(c2v.size());
			cell2vertex.setTypeA(6);   
			cell2vertex.setLayer(1);cell2vertex.setCode(ptr_el->code_physical);   
			cell2vertex.addLabel(ptr_el->node[0]);cell2vertex.addLabel(ptr_el->node[1]);
			cell2vertex.addLabel(ptr_el->node[2]);cell2vertex.addLabel(ptr_el->node[3]);
			cell2vertex.addLabel(ptr_el->node[4]);cell2vertex.addLabel(ptr_el->node[5]);    
			c2v.push_back(cell2vertex);
        }         
		// treat  pyramid 
		if  (ptr_el->type_element==7)  
        {              
			cell2vertex.clean();
			cell2vertex.setLabel(c2v.size());
			cell2vertex.setTypeA(7);    
			cell2vertex.setLayer(1);cell2vertex.setCode(ptr_el->code_physical);   
			cell2vertex.addLabel(ptr_el->node[0]);cell2vertex.addLabel(ptr_el->node[1]);
			cell2vertex.addLabel(ptr_el->node[2]);cell2vertex.addLabel(ptr_el->node[3]);
			cell2vertex.addLabel(ptr_el->node[4]);   
			c2v.push_back(cell2vertex);       
        }         
    }
}
void makeFace(vector<A2B> &c2v,vector<A2B> &f2v, vector<A2B> &c2f,vector<A2B> &f2c)
{
	A2B face2vertex,cell2vertex;
	A2B face2cell,cell2face;    
	for(size_t i=0;i<c2v.size();i++)
    {
		cell2vertex=c2v[i];
		cell2face.clean();
		cell2face.setLabel(i);
		// treat  tetrahedron 
		if(cell2vertex.getTypeA()==4)
        {
			for(size_t j=0;j<3;j++)
            {
				face2vertex.clean();face2vertex.setLabel(f2v.size());  
				face2cell.clean();face2cell.setLabel(f2v.size());
				face2vertex.setTypeA(2);   
				if(cell2vertex.getCode()>=MINUS_ONE_DIM)
                {face2vertex.setLayer(2);face2vertex.setCode(cell2vertex.getCode());}
				face2vertex.addLabel(cell2vertex.labelB[j]);
				face2vertex.addLabel(cell2vertex.labelB[(j+1)%3]);
				face2vertex.addLabel(cell2vertex.labelB[(j+2)%3]);
				face2cell.addLabel(i);
				cell2face.addLabel(f2v.size());
				f2v.push_back(face2vertex); 
				f2c.push_back(face2cell);
            }
        }
		// treat  hexahedron         
		if(cell2vertex.getTypeA()==5)
        {
			face2vertex.clean();face2vertex.setLabel(f2v.size());  
			face2cell.clean();face2cell.setLabel(f2v.size());
			face2vertex.setTypeA(3);   
			if(cell2vertex.getCode()>=MINUS_ONE_DIM)
            {face2vertex.setLayer(2);face2vertex.setCode(cell2vertex.getCode());}
			face2vertex.addLabel(cell2vertex.labelB[0]);
			face2vertex.addLabel(cell2vertex.labelB[1]);
			face2vertex.addLabel(cell2vertex.labelB[2]);
			face2vertex.addLabel(cell2vertex.labelB[3]);        
			face2cell.addLabel(i);
			cell2face.addLabel(f2v.size());
			f2v.push_back(face2vertex); 
			f2c.push_back(face2cell);              
			for(size_t j=0;j<4;j++) //the lateral face
            {
				face2vertex.clean();face2vertex.setLabel(f2v.size());  
				face2cell.clean();face2cell.setLabel(f2v.size());
				face2vertex.setTypeA(3);   
				if(cell2vertex.getCode()>=MINUS_ONE_DIM)
                {face2vertex.setLayer(2);face2vertex.setCode(cell2vertex.getCode());}
				face2vertex.addLabel(cell2vertex.labelB[j]);
				face2vertex.addLabel(cell2vertex.labelB[(j+1)%4]);
				face2vertex.addLabel(cell2vertex.labelB[4+(j+1)%4]);
				face2vertex.addLabel(cell2vertex.labelB[4+j]);        
				face2cell.addLabel(i);
				cell2face.addLabel(f2v.size());
				f2v.push_back(face2vertex); 
				f2c.push_back(face2cell);
            }
			face2vertex.clean();face2vertex.setLabel(f2v.size());  
			face2cell.clean();face2cell.setLabel(f2v.size());
			face2vertex.setTypeA(3);   
			if(cell2vertex.getCode()>=MINUS_ONE_DIM)
            {face2vertex.setLayer(2);face2vertex.setCode(cell2vertex.getCode());}
			face2vertex.addLabel(cell2vertex.labelB[4]);
			face2vertex.addLabel(cell2vertex.labelB[5]);
			face2vertex.addLabel(cell2vertex.labelB[6]);
			face2vertex.addLabel(cell2vertex.labelB[7]);        
			face2cell.addLabel(i);
			cell2face.addLabel(f2v.size());
			f2v.push_back(face2vertex); 
			f2c.push_back(face2cell);    
        }        
        // treat  prism 
		if  (cell2vertex.getTypeA()==6) 
        { 
            // upface
			face2vertex.clean();face2vertex.setLabel(f2v.size());  
			face2cell.clean();face2cell.setLabel(f2v.size());
			face2vertex.setTypeA(2);   
			if(cell2vertex.getCode()>=MINUS_ONE_DIM)
            {face2vertex.setLayer(2);face2vertex.setCode(cell2vertex.getCode());}
			face2vertex.addLabel(cell2vertex.labelB[0]);
			face2vertex.addLabel(cell2vertex.labelB[1]);
			face2vertex.addLabel(cell2vertex.labelB[2]);       
			face2cell.addLabel(i);
			cell2face.addLabel(f2v.size());
			f2v.push_back(face2vertex); 
			f2c.push_back(face2cell);              
			for(size_t j=0;j<3;j++) //the lateral face
            {
				face2vertex.clean();face2vertex.setLabel(f2v.size());  
				face2cell.clean();face2cell.setLabel(f2v.size());
				face2vertex.setTypeA(3);   
				if(cell2vertex.getCode()>=MINUS_ONE_DIM)
                {face2vertex.setLayer(2);face2vertex.setCode(cell2vertex.getCode());}
				face2vertex.addLabel(cell2vertex.labelB[j]);
				face2vertex.addLabel(cell2vertex.labelB[(j+1)%3]);
				face2vertex.addLabel(cell2vertex.labelB[3+(j+1)%3]);
				face2vertex.addLabel(cell2vertex.labelB[3+j]);        
				face2cell.addLabel(i);
				cell2face.addLabel(f2v.size());
				f2v.push_back(face2vertex); 
				f2c.push_back(face2cell);
            }
			//downface
			face2vertex.clean();face2vertex.setLabel(f2v.size());  
			face2cell.clean();face2cell.setLabel(f2v.size());
			face2vertex.setTypeA(2);   
			if(cell2vertex.getCode()>=MINUS_ONE_DIM)
            {face2vertex.setLayer(2);face2vertex.setCode(cell2vertex.getCode());}
			face2vertex.addLabel(cell2vertex.labelB[3]);
			face2vertex.addLabel(cell2vertex.labelB[4]);
			face2vertex.addLabel(cell2vertex.labelB[5]);      
			face2cell.addLabel(i);
			cell2face.addLabel(f2v.size());
			f2v.push_back(face2vertex); 
			f2c.push_back(face2cell);             
        }
		// treat  pyramid 
		if  (cell2vertex.getTypeA()==7)  
        { 
			for(size_t j=0;j<4;j++) //the lateral face
            {
				face2vertex.clean();face2vertex.setLabel(f2v.size());  
				face2cell.clean();face2cell.setLabel(f2v.size());
				face2vertex.setTypeA(2);   
				if(cell2vertex.getCode()>=MINUS_ONE_DIM)
                {face2vertex.setLayer(2);face2vertex.setCode(cell2vertex.getCode());}
				face2vertex.addLabel(cell2vertex.labelB[4]);
				face2vertex.addLabel(cell2vertex.labelB[(j+0)%4]);
				face2vertex.addLabel(cell2vertex.labelB[(j+1)%4]);       
				face2cell.addLabel(i);
				cell2face.addLabel(f2v.size());
				f2v.push_back(face2vertex); 
				f2c.push_back(face2cell);
            }   
            // the down  face
			face2vertex.clean();face2vertex.setLabel(f2v.size());  
			face2cell.clean();face2cell.setLabel(f2v.size());
			face2vertex.setTypeA(3);   
			if(cell2vertex.getCode()>=MINUS_ONE_DIM)
            {face2vertex.setLayer(2);face2vertex.setCode(cell2vertex.getCode());}
			face2vertex.addLabel(cell2vertex.labelB[0]);
			face2vertex.addLabel(cell2vertex.labelB[1]);
			face2vertex.addLabel(cell2vertex.labelB[2]);
			face2vertex.addLabel(cell2vertex.labelB[3]);        
			face2cell.addLabel(i);
			cell2face.addLabel(f2v.size());
			f2v.push_back(face2vertex); 
			f2c.push_back(face2cell);              
        }
		c2f.push_back(cell2face);    
    }
}
void makeEdge(vector<A2B> &f2v,vector<A2B> &e2v, vector<A2B> &f2e,vector<A2B> &e2f)
{
	A2B edge2vertex,face2vertex;
	A2B edge2face,face2edge;    
	for(size_t i=0;i<f2v.size();i++)
    {
		face2vertex=f2v[i];
		face2edge.clean();
		face2edge.setLabel(i);
		// triangle         
		if  (face2vertex.getTypeA()==2)
        {
			for(size_t j=0;j<3;j++)
            {
				edge2vertex.clean();edge2vertex.setLabel(e2v.size());  
				edge2face.clean();edge2face.setLabel(e2v.size());
				edge2vertex.setTypeA(1);   
				if((face2vertex.getLayer()==1)||(face2vertex.getCode()>=MINUS_ONE_DIM))
                {edge2vertex.setLayer(2);edge2vertex.setCode(face2vertex.getCode());}
				if((face2vertex.getLayer()==2)||(face2vertex.getCode()>=MINUS_TWO_DIM))
                {edge2vertex.setLayer(3);edge2vertex.setCode(face2vertex.getCode());}                
				edge2vertex.addLabel(face2vertex.labelB[j]);
				edge2vertex.addLabel(face2vertex.labelB[(j+1)%3]);
				edge2face.addLabel(i);
				face2edge.addLabel(e2v.size());
				e2v.push_back(edge2vertex); 
				e2f.push_back(edge2face);
            }            
        }
		//  quadrilateral    
		if  (face2vertex.getTypeA()==3 )  
        {
			for(size_t j=0;j<4;j++)
            {
				edge2vertex.clean();edge2vertex.setLabel(e2v.size());  
				edge2face.clean();edge2face.setLabel(e2v.size());
				edge2vertex.setTypeA(1);
				if((face2vertex.getLayer()==1)||(face2vertex.getCode()>=MINUS_ONE_DIM))
                {edge2vertex.setLayer(2);edge2vertex.setCode(face2vertex.getCode());}
				if((face2vertex.getLayer()==2)||(face2vertex.getCode()>=MINUS_TWO_DIM))
                {edge2vertex.setLayer(3);edge2vertex.setCode(face2vertex.getCode());}  
				edge2vertex.addLabel(face2vertex.labelB[j]);
				edge2vertex.addLabel(face2vertex.labelB[(j+1)%4]);
				edge2face.addLabel(i);
				face2edge.addLabel(e2v.size());
				e2v.push_back(edge2vertex); 
				e2f.push_back(edge2face);
            }                
        }
		f2e.push_back(face2edge);    
    }        
}
void makeFaceMap(vector<size_t> &fMap,vector<A2B> &f2v,vector<A2B> &c2f,vector<A2B> &f2c)
{
	FVSparseV<size_t> markFaceSum,markFaceProd;    
	A2B face2vertex_ref,face2vertex;   
	size_t sum,prod;
	bool doublon;
	for(size_t i=0;i<f2v.size();i++)
    {
		size_t j;
		fMap[i]=i;
		doublon=false;
		face2vertex_ref=f2v[i];
		sum=face2vertex_ref.getSum();prod=face2vertex_ref.getProd();
		if(!markFaceSum.getValue(sum))   continue; // it is not tagged by the sum
		if(!markFaceProd.getValue(prod)) continue; // it is not tagged by the prod
		// ok here the cell has potentially a doublon find it
		for(j=0;j<i;j++)
        {
			face2vertex=f2v[j];
			if(sum!=face2vertex.getSum()) continue;    
			if(prod!=face2vertex.getProd()) continue; 
			if(!face2vertex.isEqual(face2vertex)) continue;
			// here we know that we have a doublon
			doublon=true;
        }
		if(doublon)
        {
			// we make the mapping + the physical code   
			fMap[i]=j;
			// change the code
			if((face2vertex_ref.getCode()==0) && (face2vertex.getCode()<MINUS_ONE_DIM)) 
				face2vertex_ref.setCode(face2vertex.getCode());
			if((face2vertex_ref.getCode()==0) && (face2vertex.getLayer()==2)
			   && (face2vertex.getCode()<MINUS_TWO_DIM)) 
				face2vertex_ref.setCode(face2vertex.getCode()%MINUS_ONE_DIM);        
			face2vertex.setLabel(0); // forgot this element
			// change the link
			//face2cell=f2c[j];
        }
		else
        {
			markFaceSum.setValue(sum,1);
			markFaceProd.setValue(prod,1);
        }
    }
}
void makeEdgeMap(vector<size_t> &eMap,vector<A2B> &f2e,vector<A2B> &e2f)
{
	//
}
// convert a Gmsh struct into a FVMesh3D   
void FVMesh3D::FastGmsh2FVMesh(Gmsh &m)
{
	if (m.getDim()!=3)
	{
		cout<<" Dimension don't match. The Gmsh mesh should contain only 0 or 1D elements"<<endl;
		return;
	}  
	_nb_vertex=m.getNbNode();
	_name="mesh3D  from gmsh file"; 
	_vertex.resize(_nb_vertex); 
	for(size_t i=0;i<_nb_vertex;i++)
    {
		_vertex[i].label=m.getNode(i)->label;
		_vertex[i].coord.x=m.getNode(i)->coord.x;
		_vertex[i].coord.y=m.getNode(i)->coord.y;
		_vertex[i].coord.z=m.getNode(i)->coord.z;
		_vertex[i].code=0;  // default value
    }     
	_nb_cell=0;_nb_edge=0;_nb_face=0;  
	A2B edge2vertex,face2vertex,cell2vertex;
	A2B edge2face,face2edge;
	A2B face2cell,cell2face;
	vector<A2B> e2v,f2v,c2v;
	vector<A2B> f2e,e2f;
	vector<A2B> c2f,f2c;
	vector<size_t> eMap,fMap;
	loadMesh(m,c2v,f2v,e2v);
	makeFace(c2v,f2v,c2f,f2c);
	makeEdge(f2v,e2v,f2e,e2f);
	makeFaceMap(fMap,f2v,c2f,f2c);
	makeEdgeMap(eMap,f2e,e2f);
	_cell.resize(0);_face.resize(0);_edge.resize(0);
	_cell.reserve(m.getNbElement());_face.reserve(m.getNbElement()*3);_edge.reserve(m.getNbElement()*5);
	//FVMesh3D::complete_data();
}
// make the dual mesh
void FVMesh3D::makeDualMesh(FVMesh3D &pm)
{
	size_t nb_primal_vertex,nb_primal_cell,nb_primal_edge,nb_primal_face,nb_primal_boundary_face,index,count,edges_count,nb_inner_edges_boundary_faces;
	size_t nb_dual_faces=0,nb_dual_boundary_faces=0;
	vector<size_t> edges_on_cell(NB_EDGE_PER_FACE_3D*NB_FACE_PER_CELL_3D);
	bool go_on;
	FVVertex3D v;  
	FVEdge3D e;
	FVFace3D f;
	FVCell3D c;
	FVCell3D *ptr_c;
	FVFace3D *ptr_f,*ptr_ff;
	FVEdge3D *ptr_e,*ptr_ee;
	FVVertex3D *ptr_v;
	// copy the data
	_name="dual"+pm.getName();
	pm._dualMesh=this; 
	_dualMesh=&pm; // the dual mesh of the dual mesh is the primal mesh
	nb_primal_vertex=pm.getNbVertex();
	nb_primal_cell=pm.getNbCell();
	nb_primal_edge=pm.getNbEdge();
	nb_primal_face=pm.getNbFace();
	nb_primal_boundary_face=pm.getNbBoundaryFace();
	// resize vectors
	_vertex.resize(nb_primal_vertex+nb_primal_cell+nb_primal_boundary_face);
	_cell.resize(nb_primal_face);
	// get the number of dual edges
	// one edge per primal vertex in primal cell and in primal boundary face and per primal edge
	size_t nb_dual_edges=nb_primal_edge;
	for(size_t i=0;i<nb_primal_cell;i++)
		nb_dual_edges+=pm.getCell(i)->nb_vertex;
	for(size_t i=0;i<nb_primal_boundary_face;i++)
		nb_dual_edges+=pm.getBoundaryFace(i)->nb_vertex;
	_edge.resize(nb_dual_edges);
	// get the number of dual faces
	// one face per primal edge in primal cell and in primal boundary face
	for(size_t i=0;i<nb_primal_cell;i++)
		nb_dual_faces+=pm.getCell(i)->nb_face+pm.getCell(i)->nb_vertex-2;
	for(size_t i=0;i<nb_primal_boundary_face;i++)
	{
		nb_dual_faces+=pm.getBoundaryFace(i)->nb_edge;
		nb_dual_boundary_faces+=pm.getBoundaryFace(i)->nb_edge;
	}
	_face.resize(nb_dual_faces);
	_boundary_face.resize(nb_dual_boundary_faces);
	//// DUAL VERTICES
	// copy the vertices of the primal mesh into the dual mesh
	// the number of vertices is nb_primal_vertex with the same order of the primal vertices
	for(size_t i=0;i<nb_primal_vertex;i++)
	{
		_vertex[i]=*pm.getVertex(i);
		_vertex[i].nb_cell=0;
	}
	// take the centroid of the cells
	// the centroids of the primal mesh are also vertices of the dual mesh
	index=nb_primal_vertex;
	for(size_t i=0;i<nb_primal_cell;i++)     
	{
		v.label=index+i+1;
		v.coord=pm.getCell(i)->centroid; 
		v.code=pm.getCell(i)->code;
		v.nb_cell=0;
		_vertex[index+i]=v;
	} 
	// take the centroid of the boundary face
	// the centroids of the boundary face of the primal mesh are also vertices of the dual mesh
	index=nb_primal_vertex+nb_primal_cell;
	for(size_t i=0;i<nb_primal_boundary_face;i++)     
	{
		v.label=index+i+1;
		v.coord=pm.getBoundaryFace(i)->centroid;
		v.code=pm.getBoundaryFace(i)->code;
		v.nb_cell=0;   
		_vertex[index+i]=v;
	}  
	// the number of vertices is nb_primal_vertex+nb_primal_cell+nb_primal_boundary_edge using the boundary edge order  
	_nb_vertex=_vertex.size();
	//// DUAL CELLS
	// initialize the dual cells
	for(size_t i=0;i<nb_primal_face;i++)     
	{
		ptr_f=pm.getFace(i);
		_cell[i].label=i+1;
		_cell[i].nb_vertex=0;
		_cell[i].nb_face=0;
		_cell[i].code=ptr_f->code;
		_cell[i].primalFace=ptr_f;
		// associate the vertices to this dual cell and the dual cell to the vertices
		for(size_t j=0;j<ptr_f->nb_vertex;j++)
		{
			_cell[i].vertex[_cell[i].nb_vertex++]=&_vertex[ptr_f->vertex[j]->label-1];
			ptr_f->vertex[j]->cell[ptr_f->vertex[j]->nb_cell++]=&_cell[i];
			if(ptr_f->vertex[j]->nb_cell>=NB_CELL_PER_VERTEX_3D)
			{
				cout<<"Error: overflow in class FVVertex3D, too many Cells, found "<<ptr_f->vertex[j]->nb_cell<<" on vertex "<<ptr_f->vertex[j]->label<<". Line: "<<__LINE__<<". File: "<<__FILE__<<endl;
				exit(0);
			}
		}
		// add the vertex from the centroid of the left cell
		_cell[i].vertex[_cell[i].nb_vertex++]=&_vertex[nb_primal_vertex+ptr_f->leftCell->label-1];
		_vertex[nb_primal_vertex+ptr_f->leftCell->label-1].cell[_vertex[nb_primal_vertex+ptr_f->leftCell->label-1].nb_cell++]=&_cell[i];
		if(_vertex[nb_primal_vertex+ptr_f->leftCell->label-1].nb_cell>=NB_CELL_PER_VERTEX_3D)
		{
			cout<<"Error: overflow in class FVVertex3D, too many Cells, found "<<_vertex[nb_primal_vertex+ptr_f->leftCell->label-1].nb_cell<<" on vertex "<<_vertex[nb_primal_vertex+ptr_f->leftCell->label-1].label<<". Line: "<<__LINE__<<". File: "<<__FILE__<<endl;
			exit(0);
		}
		// add the vertex from the right cell (case exists)
		if(ptr_f->rightCell)
		{
			_cell[i].vertex[_cell[i].nb_vertex++]=&_vertex[nb_primal_vertex+ptr_f->rightCell->label-1];
			_vertex[nb_primal_vertex+ptr_f->rightCell->label-1].cell[_vertex[nb_primal_vertex+ptr_f->rightCell->label-1].nb_cell++]=&_cell[i];
			if(_vertex[nb_primal_vertex+ptr_f->rightCell->label-1].nb_cell>=NB_CELL_PER_VERTEX_3D)
			{
				cout<<"Error: overflow in class FVVertex3D, too many Cells, found "<<_vertex[nb_primal_vertex+ptr_f->rightCell->label-1].nb_cell<<" on vertex "<<_vertex[nb_primal_vertex+ptr_f->rightCell->label-1].label<<". Line: "<<__LINE__<<". File: "<<__FILE__<<endl;
				exit(0);
			}
		}
		// it is missing the vertex from the centroid of boundary face
	}
	// add the vertices from the centroids of the boundary faces to the dual cells
	for(size_t i=0;i<nb_primal_boundary_face;i++)
	{
		ptr_f=pm.getBoundaryFace(i);
		_cell[ptr_f->label-1].vertex[_cell[ptr_f->label-1].nb_vertex++]=&_vertex[nb_primal_vertex+nb_primal_cell+i];
		_vertex[nb_primal_vertex+nb_primal_cell+i].cell[_vertex[nb_primal_vertex+nb_primal_cell+i].nb_cell++]=&_cell[ptr_f->label-1];
		if(_vertex[nb_primal_vertex+nb_primal_cell+i].nb_cell>=NB_CELL_PER_VERTEX_3D)
		{
			cout<<"Error: overflow in class FVVertex3D, too many Cells, found "<<_vertex[nb_primal_vertex+nb_primal_cell+i].nb_cell<<" on vertex "<<_vertex[nb_primal_vertex+nb_primal_cell+i].label<<". Line: "<<__LINE__<<". File: "<<__FILE__<<endl;
			exit(0);
		}
	}
	//// DUAL EDGES
	// the dual edges from the primal edges
	index=0;
	for(size_t i=0;i<nb_primal_edge;i++)
	{
		ptr_e=pm.getEdge(i);
		e.label=index+1;
		e.firstVertex=&_vertex[ptr_e->firstVertex->label-1];
		e.secondVertex=&_vertex[ptr_e->secondVertex->label-1];
		e.code=pm.getEdge(i)->code;
		e.nb_vertex=2;
		_edge[index]=e;
		index++;
	}
	// the dual edges from cell vertices
	for(size_t i=0;i<nb_primal_cell;i++)     
	{
		ptr_c=pm.getCell(i);
		for(size_t j=0;j<ptr_c->nb_vertex;j++)
		{
			e.label=index+1;
			e.firstVertex=&_vertex[nb_primal_vertex+pm.getCell(i)->label-1]; // the centroid
			e.secondVertex=&_vertex[ptr_c->vertex[j]->label-1]; // the vertex on cell
			e.code=pm.getCell(i)->code;
			e.nb_vertex=2;
			_edge[index]=e;
			index++;
		}
	}
	nb_inner_edges_boundary_faces=index;
	// the dual edges from boundary
	for(size_t i=0;i<nb_primal_boundary_face;i++)     
	{
		ptr_f=pm.getBoundaryFace(i);
		for(size_t j=0;j<ptr_f->nb_vertex;j++)
		{
			e.label=index+1;
			e.firstVertex=&_vertex[nb_primal_vertex+nb_primal_cell+i]; // the centroid of the face
			e.secondVertex=&_vertex[ptr_f->vertex[j]->label-1]; //the vertex on edge
			e.code=ptr_f->code;
			e.nb_vertex=2;
			//e.nb_cell=1;   
			_edge[index]=e;
			index++;
		}
	}
	//// DUAL FACES
	// dual faces from primal cells
	edges_count=0;
	index=0;
	count=0;
	for(size_t i=0;i<nb_primal_cell;i++)
	{
		ptr_c=pm.getCell(i);
		// clear the vector of the edges of the cell
		for(size_t l=0;l<edges_count;l++)
			edges_on_cell[l]=0;
		edges_count=0;
		for(size_t j=0;j<ptr_c->nb_face;j++)
		{
			ptr_f=ptr_c->face[j];
			for(size_t k=0;k<ptr_f->nb_edge;k++)
			{
				ptr_e=ptr_f->edge[k];
				// avoid the same edge is use more than once
				go_on=true;
				for(size_t l=0;l<edges_count;l++)
				{
					if(ptr_e->label==edges_on_cell[l])
					{
						go_on=false;
						break;
					}
				}
				// if not, add this edge to the vector
				if(go_on)
				{
					edges_on_cell[edges_count++]=ptr_e->label;
					f.nb_cell=2;
					f.nb_edge=0;
					f.nb_vertex=0;
					// link the dual face to the mother primal cell
					f.dualCell=ptr_c;
					// each dual face has exactly 3 vertices
					f.vertex[f.nb_vertex++]=&_vertex[ptr_e->firstVertex->label-1]; // first vertex of the edge
					f.vertex[f.nb_vertex++]=&_vertex[ptr_e->secondVertex->label-1]; // second vertex of the edge
					f.vertex[f.nb_vertex++]=&_vertex[nb_primal_vertex+i]; // vertex associated to the centroid of the cell
					// each dual face has exactly 3 edges
					f.edge[f.nb_edge++]=&_edge[ptr_e->label-1]; // common edge to the primal cell
					// find the other two edges associated to this face
					for(size_t l=0;l<ptr_c->nb_vertex;l++)
					{
						ptr_v=ptr_c->vertex[l];
						if(ptr_v->label==ptr_e->firstVertex->label || ptr_v->label==ptr_e->secondVertex->label)
							f.edge[f.nb_edge++]=&_edge[nb_primal_edge+count+l];
					}
					if(f.nb_edge!=3)
					{
						cout<<"Error: I could not find all the edges! Line: "<<__LINE__<<". File: "<<__FILE__<<endl;
						exit(0);
					}
					f.label=index+1;
					f.code=ptr_c->code;
					_face[index]=f;
					index++;
					//// associate this face to the dual cells and the dual cells to this face
					for(size_t m=0;m<ptr_c->nb_face;m++)
					{
						ptr_ff=ptr_c->face[m];
						for(size_t n=0;n<ptr_ff->nb_edge;n++)
						{
							ptr_ee=ptr_ff->edge[n];
							if(ptr_ee->label==ptr_e->label)
							{
								_cell[ptr_ff->label-1].face[_cell[ptr_ff->label-1].nb_face++]=&_face[index-1];
								if(_face[index-1].leftCell)
									_face[index-1].rightCell=&_cell[ptr_ff->label-1];
								else
									_face[index-1].leftCell=&_cell[ptr_ff->label-1];
								break;
							}
						}
						// leave the loop if the two cells which share this face were already found
						if(_face[index-1].rightCell!=NULL && _face[index-1].leftCell!=NULL)
							break;
					}
				}
			}
		}
		count+=ptr_c->nb_vertex;
	}
	// dual faces from primal boundary faces
	count=0;
	for(size_t i=0;i<nb_primal_boundary_face;i++)
	{
		ptr_f=pm.getBoundaryFace(i);
		for(size_t j=0;j<ptr_f->nb_edge;j++)
		{
			ptr_e=ptr_f->edge[j];
			f.nb_cell=1;
			f.nb_edge=0;
			f.nb_vertex=0;
			// link the dual face to the mother primal cell
			f.dualCell=ptr_f->leftCell;
			// each dual face has exactly 3 vertices
			f.vertex[f.nb_vertex++]=&_vertex[ptr_e->firstVertex->label-1]; // first vertex of the edge
			f.vertex[f.nb_vertex++]=&_vertex[ptr_e->secondVertex->label-1]; // second vertex of the edge
			f.vertex[f.nb_vertex++]=&_vertex[nb_primal_vertex+nb_primal_cell+i]; // vertex associated to the centroid of the cell
			// each dual face has exactly 3 edges
			f.edge[f.nb_edge++]=&_edge[ptr_e->label-1]; // common edge to the primal cell
			// find the other two edges associated to this face
			for(size_t l=0;l<ptr_f->nb_vertex;l++)
			{
				ptr_v=ptr_f->vertex[l];
				if(ptr_v->label==ptr_e->firstVertex->label || ptr_v->label==ptr_e->secondVertex->label)
					f.edge[f.nb_edge++]=&_edge[nb_inner_edges_boundary_faces+count+l];
			}
			if(f.nb_edge!=3)
			{
				cout<<"Error: I could not find all the edges! Line: "<<__LINE__<<". File: "<<__FILE__<<endl;
				exit(0);
			}
			f.label=index+1;
			f.code=ptr_f->code;
			_face[index++]=f;
			_boundary_face[_nb_boundary_face++]=&_face[index-1];
			//// associate this face to the dual cells and the dual cells to this face
			_face[index-1].rightCell=NULL;
			_face[index-1].leftCell=&_cell[ptr_f->label-1];
			_cell[ptr_f->label-1].face[_cell[ptr_f->label-1].nb_face++]=&_face[index-1];
		}
		count+=ptr_f->nb_edge;
	}
	//// COMPLETE THE DATA
	_nb_vertex=_vertex.size();
	_nb_edge=_edge.size();
	_nb_face=_face.size();
	_nb_cell=_cell.size();
	_dim=3;
	_nb_boundary_face=0;
	for(size_t i=0;i<nb_primal_boundary_face;i++)
	{
		ptr_f=pm.getBoundaryFace(i);
		_nb_boundary_face+=ptr_f->nb_edge;
	}
	//// LINK PRIMAL-DUAL
	// make the link primalFace->dualCell for the faces of the primal mesh
	for(size_t i=0;i<nb_primal_face;i++)
		pm.getFace(i)->dualCell=&_cell[i];
	// the link dual-primal was already performed
	//// COMPLETE WITH THE GEOMETRICAL DATA
	FVPoint3D<double> u,vv,w;
	double no;
	// EDGE
	// compute the centroid and length of edge
	for(size_t i=0;i<_nb_edge;i++)
    {
		_edge[i].centroid=(_edge[i].firstVertex->coord+_edge[i].secondVertex->coord)*0.5;
		_edge[i].mass_center=_edge[i].centroid;
		u=_edge[i].firstVertex->coord-_edge[i].secondVertex->coord;
		_edge[i].length=Norm(u);
    }
	// FACE
	// compute geometric stuff for the face
	for(size_t i=0;i<_nb_face;i++)
    {
		// conpute perimeter of the face 
		_face[i].perimeter=0.0;
		for(size_t j=0;j<_face[i].nb_edge;j++)
			_face[i].perimeter+=_face[i].edge[j]->length; 
		// compute the area of the face
		FVPoint3D<double> P,P1,P2,Paux;
		FVPoint3D<double> GPP;
		FVGaussPoint2D G2D;
		// compute the centroid of the face
		_face[i].centroid=0.0;
		for(size_t k=0;k<_face[i].nb_vertex;k++)
			_face[i].centroid+=_face[i].vertex[k]->coord;    
		_face[i].centroid/=_face[i].nb_vertex; 
		// face area and mass center
		_face[i].mass_center=0.0;
		_face[i].area=0.0;
		for(size_t j=0;j<_face[i].nb_edge;j++)
        {
			P1=_face[i].edge[j]->firstVertex->coord;
			P2=_face[i].edge[j]->secondVertex->coord;
			u=P1-_face[i].centroid; 
			vv=P2-_face[i].centroid;
			w=CrossProduct(u,vv);
			_face[i].LocalArea[j]=Norm(w)*0.5;
			_face[i].area+=_face[i].LocalArea[j];
			Paux=0.0;
			for (size_t k=1;k<=G2D.getNbPoint(5);k++)
            { 
				GPP=G2D.getPoint(5,k);
				P=GPP.x*P1+GPP.y*P2+GPP.z*_face[i].centroid;
				Paux+=G2D.getWeight(5,k)*P;
            }        
			_face[i].mass_center+=Paux*_face[i].LocalArea[j]; 
        }  
		_face[i].mass_center/=_face[i].area; 
		// we compute the normal from left to rigth for each sub-triangle  
		_face[i].normal=0.0;
		double surf=0.0;
		for(size_t j=0;j<_face[i].nb_edge;j++)
        {
			u=_face[i].edge[j]->firstVertex->coord-_face[i].centroid;
			vv=_face[i].edge[j]->secondVertex->coord-_face[i].centroid;
			w=CrossProduct(u,vv);
			no=Norm(w);
			w/=no;
			_face[i].LocalNormal[j]=w;
			if(!_face[i].leftCell){cout<<" ERROR no left cell for face "<<i<<endl;fflush(NULL);}
			u=_face[i].centroid-_face[i].leftCell->centroid;          
			if(w*u<0) _face[i].LocalNormal[j]*=-1.;
			surf+=no;
			_face[i].normal+=no*_face[i].LocalNormal[j];          
		}
		_face[i].normal/=surf;
    } 
    // CELL
	// compute geometric stuff for the cell
	for(size_t i=0;i<_nb_cell;i++)
    {
		// compute surface of the cell 
		_cell[i].surface=0.0;
		for(size_t j=0;j<_cell[i].nb_face;j++)
			_cell[i].surface+=_cell[i].face[j]->area; 
		// compute the centroid of the cell
		_cell[i].centroid=0.0;
		for(size_t j=0;j<_cell[i].nb_vertex;j++)
			_cell[i].centroid+=_cell[i].vertex[j]->coord;
		_cell[i].centroid/=_cell[i].nb_vertex; 
		// compute the volume of the cell
		_cell[i].volume=0.;
		_cell[i].mass_center=0.;    
		FVPoint3D<double> P,P1,P2,Paux;
		FVPoint4D<double> GPPP;
		FVGaussPoint3D G3D; 
		double Vaux;
		FVPoint3D<double> u,v,w;
		for(size_t j=0;j<_cell[i].nb_face;j++)
        {    
			_cell[i].cell2face[j]= _cell[i].face[j]->centroid-_cell[i].centroid;    
			for(size_t k=0;k<_cell[i].face[j]->nb_edge;k++)
            {
				P1=_cell[i].face[j]->edge[k]->firstVertex->coord;
				P2=_cell[i].face[j]->edge[k]->secondVertex->coord;
				u=_cell[i].cell2face[j]; 
				v=P1-_cell[i].centroid; 
				w=P2-_cell[i].centroid; 
				Vaux=fabs(Det(u,v,w))/6; 
				_cell[i].volume+=Vaux;
				Paux=0;
				for(size_t r=1;r<=G3D.getNbPoint(5);r++)
                {              
					GPPP=G3D.getPoint(5,r);
					P=GPPP.x*P1+GPPP.y*P2+GPPP.z*_cell[i].face[j]->centroid+GPPP.t*_cell[i].centroid;
					Paux+=G3D.getWeight(5,r)*P;
                }
				_cell[i].mass_center+=Paux*Vaux;    
            }           
        }
		_cell[i].mass_center/=_cell[i].volume;   
    }	
}
// make the dual mesh for post processing proposes
void FVMesh3D::makePostDualMesh(FVMesh3D &pm)
{
	size_t nb_primal_vertex,nb_primal_cell,nb_primal_edge,nb_primal_face,nb_primal_boundary_face,index,count,edges_count,nb_inner_edges_boundary_faces;
	size_t nb_dual_faces=0,nb_dual_boundary_faces=0;
	vector<size_t> edges_on_cell(NB_EDGE_PER_FACE_3D*NB_FACE_PER_CELL_3D);
	bool go_on;
	FVVertex3D v;   
	FVEdge3D e;
	FVFace3D f;
	FVCell3D c;
	FVCell3D *ptr_c;
	FVFace3D *ptr_f,*ptr_ff;
	FVEdge3D *ptr_e,*ptr_ee;
	FVVertex3D *ptr_v;
	// copy the data
	_name="dual"+pm.getName();
	pm._dualMesh=this; 
	_dualMesh=&pm; // the dual mesh of the dual mesh is the primal mesh
	nb_primal_vertex=pm.getNbVertex();
	nb_primal_cell=pm.getNbCell();
	nb_primal_edge=pm.getNbEdge();
	nb_primal_face=pm.getNbFace();
	nb_primal_boundary_face=pm.getNbBoundaryFace();
	// resize vectors
	_vertex.resize(nb_primal_vertex+nb_primal_cell+nb_primal_boundary_face);
	_cell.resize(nb_primal_face*2); //<------
	// get the number of dual edges
	// one edge per primal vertex in primal cell and in primal boundary face and per primal edge
	size_t nb_dual_edges=nb_primal_edge;
	for(size_t i=0;i<nb_primal_cell;i++)
		nb_dual_edges+=pm.getCell(i)->nb_vertex;
	for(size_t i=0;i<nb_primal_boundary_face;i++)
		nb_dual_edges+=pm.getBoundaryFace(i)->nb_vertex;
	_edge.resize(nb_dual_edges);
	// get the number of dual faces
	// one face per primal edge in primal cell and in primal boundary face
	for(size_t i=0;i<nb_primal_cell;i++)
		nb_dual_faces+=pm.getCell(i)->nb_face+pm.getCell(i)->nb_vertex-2;
	for(size_t i=0;i<nb_primal_boundary_face;i++)
	{
		nb_dual_faces+=pm.getBoundaryFace(i)->nb_edge;
		nb_dual_boundary_faces+=pm.getBoundaryFace(i)->nb_edge;
	}
	nb_dual_faces+=nb_primal_face;
	_face.resize(nb_dual_faces);
	_boundary_face.resize(nb_dual_boundary_faces);
	//// DUAL VERTICES
	// copy the vertices of the primal mesh into the dual mesh
	// the number of vertices is nb_primal_vertex with the same order of the primal vertices
	for(size_t i=0;i<nb_primal_vertex;i++)
	{
		_vertex[i]=*pm.getVertex(i);
		_vertex[i].nb_cell=0;
	}
	// take the centroid of the cells
	// the centroids of the primal mesh are also vertices of the dual mesh
	index=nb_primal_vertex;
	for(size_t i=0;i<nb_primal_cell;i++)     
	{
		v.label=index+i+1;
		v.coord=pm.getCell(i)->centroid; 
		v.code=pm.getCell(i)->code;
		v.nb_cell=0;
		_vertex[index+i]=v;
	} 
	// take the centroid of the boundary face
	// the centroids of the boundary face of the primal mesh are also vertices of the dual mesh
	index=nb_primal_vertex+nb_primal_cell;
	for(size_t i=0;i<nb_primal_boundary_face;i++)     
	{
		v.label=index+i+1;
		v.coord=pm.getBoundaryFace(i)->centroid;
		v.code=pm.getBoundaryFace(i)->code;
		v.nb_cell=0;   
		_vertex[index+i]=v;
	}
	// the number of vertices is nb_primal_vertex+nb_primal_cell+nb_primal_boundary_edge using the boundary edge order  
	_nb_vertex=_vertex.size();
	//// DUAL CELLS
	// initialize the dual cells
	for(size_t i=0;i<nb_primal_face;i++)     
	{
		ptr_f=pm.getFace(i);
		//// first dual cell associated to this face
		_cell[i*2].label=i*2+1;
		_cell[i*2].nb_vertex=0;
		_cell[i*2].nb_face=0;
		_cell[i*2].code=ptr_f->code;
		_cell[i*2].primalFace=ptr_f;
		//// second dual cell associated to this face
		_cell[i*2+1].label=i*2+2;
		_cell[i*2+1].nb_vertex=0;
		_cell[i*2+1].nb_face=0;
		_cell[i*2+1].code=ptr_f->code;
		_cell[i*2+1].primalFace=ptr_f;
		// associate the vertices of the face to the two dual cells
		for(size_t j=0;j<ptr_f->nb_vertex;j++)
		{
			_cell[i*2].vertex[_cell[i*2].nb_vertex++]=&_vertex[ptr_f->vertex[j]->label-1];
			_cell[i*2+1].vertex[_cell[i*2+1].nb_vertex++]=&_vertex[ptr_f->vertex[j]->label-1];
			ptr_f->vertex[j]->cell[ptr_f->vertex[j]->nb_cell++]=&_cell[i];
			ptr_f->vertex[j]->cell[ptr_f->vertex[j]->nb_cell++]=&_cell[i*2+1];
			if(ptr_f->vertex[j]->nb_cell>=NB_CELL_PER_VERTEX_3D)
			{
				cout<<"Error: overflow in class FVVertex3D, too many Cells, found "<<ptr_f->vertex[j]->nb_cell<<" on vertex "<<ptr_f->vertex[j]->label<<". Line: "<<__LINE__<<". File: "<<__FILE__<<endl;
				exit(0);
			}
		}
		// add the vertex from the centroid of the left cell to the first dual cell
		_cell[i*2].vertex[_cell[i*2].nb_vertex++]=&_vertex[nb_primal_vertex+ptr_f->leftCell->label-1];
		_vertex[nb_primal_vertex+ptr_f->leftCell->label-1].cell[_vertex[nb_primal_vertex+ptr_f->leftCell->label-1].nb_cell++]=&_cell[i*2];
		if(_vertex[nb_primal_vertex+ptr_f->leftCell->label-1].nb_cell>=NB_CELL_PER_VERTEX_3D)
		{
			cout<<"Error: overflow in class FVVertex3D, too many Cells, found "<<_vertex[nb_primal_vertex+ptr_f->leftCell->label-1].nb_cell<<" on vertex "<<_vertex[nb_primal_vertex+ptr_f->leftCell->label-1].label<<". Line: "<<__LINE__<<". File: "<<__FILE__<<endl;
			exit(0);
		}
		// add the vertex from the centroid of the right cell to the second dual cell (if exists)
		if(ptr_f->rightCell)
		{
			_cell[i*2+1].vertex[_cell[i*2+1].nb_vertex++]=&_vertex[nb_primal_vertex+ptr_f->rightCell->label-1];
			_vertex[nb_primal_vertex+ptr_f->rightCell->label-1].cell[_vertex[nb_primal_vertex+ptr_f->rightCell->label-1].nb_cell++]=&_cell[i*2+1];
			if(_vertex[nb_primal_vertex+ptr_f->rightCell->label-1].nb_cell>=NB_CELL_PER_VERTEX_3D)
			{
				cout<<"Error: overflow in class FVVertex3D, too many Cells, found "<<_vertex[nb_primal_vertex+ptr_f->rightCell->label-1].nb_cell<<" on vertex "<<_vertex[nb_primal_vertex+ptr_f->rightCell->label-1].label<<". Line: "<<__LINE__<<". File: "<<__FILE__<<endl;
				exit(0);
			}
		}
		// it is missing the vertex from the centroid of boundary face
	}
	// add the vertices from the centroids of the boundary faces to the dual cells
	for(size_t i=0;i<nb_primal_boundary_face;i++)
	{
		ptr_f=pm.getBoundaryFace(i);
		_cell[(ptr_f->label-1)*2+1].vertex[_cell[(ptr_f->label-1)*2+1].nb_vertex++]=&_vertex[nb_primal_vertex+nb_primal_cell+i];
		_vertex[nb_primal_vertex+nb_primal_cell+i].cell[_vertex[nb_primal_vertex+nb_primal_cell+i].nb_cell++]=&_cell[(ptr_f->label-1)*2+1];
		if(_vertex[nb_primal_vertex+nb_primal_cell+i].nb_cell>=NB_CELL_PER_VERTEX_3D)
		{
			cout<<"Error: overflow in class FVVertex3D, too many Cells, found "<<_vertex[nb_primal_vertex+nb_primal_cell+i].nb_cell<<" on vertex "<<_vertex[nb_primal_vertex+nb_primal_cell+i].label<<". Line: "<<__LINE__<<". File: "<<__FILE__<<endl;
			exit(0);
		}
	}
	//// DUAL EDGES
	// the dual edges from the primal edges
	index=0;
	for(size_t i=0;i<nb_primal_edge;i++)
	{
		ptr_e=pm.getEdge(i);
		e.label=index+1;
		e.firstVertex=&_vertex[ptr_e->firstVertex->label-1];
		e.secondVertex=&_vertex[ptr_e->secondVertex->label-1];
		e.code=pm.getEdge(i)->code;
		e.nb_vertex=2;
		_edge[index]=e;
		index++;
	}
	// the dual edges from the cell vertices
	for(size_t i=0;i<nb_primal_cell;i++)     
	{
		ptr_c=pm.getCell(i);
		for(size_t j=0;j<ptr_c->nb_vertex;j++)
		{
			e.label=index+1;
			e.firstVertex=&_vertex[nb_primal_vertex+pm.getCell(i)->label-1]; // the centroid
			e.secondVertex=&_vertex[ptr_c->vertex[j]->label-1]; // the vertex on cell
			e.code=pm.getCell(i)->code;
			e.nb_vertex=2;
			_edge[index]=e;
			index++;
		}
	}
	nb_inner_edges_boundary_faces=index;
	// the dual edges from boundary
	for(size_t i=0;i<nb_primal_boundary_face;i++)     
	{
		ptr_f=pm.getBoundaryFace(i);
		for(size_t j=0;j<ptr_f->nb_vertex;j++)
		{
			e.label=index+1;
			e.firstVertex=&_vertex[nb_primal_vertex+nb_primal_cell+i]; // the centroid of the face
			e.secondVertex=&_vertex[ptr_f->vertex[j]->label-1]; //the vertex on edge
			e.code=ptr_f->code;
			e.nb_vertex=2;  
			_edge[index]=e;
			index++;
		}
	}
	//// DUAL FACES
	// dual faces from primal cells
	edges_count=0;
	index=0;
	count=0;
	for(size_t i=0;i<nb_primal_cell;i++)
	{
		ptr_c=pm.getCell(i);
		// clear the vector of the edges of the cell
		for(size_t l=0;l<edges_count;l++)
			edges_on_cell[l]=0;
		edges_count=0;
		for(size_t j=0;j<ptr_c->nb_face;j++)
		{
			ptr_f=ptr_c->face[j];
			for(size_t k=0;k<ptr_f->nb_edge;k++)
			{
				ptr_e=ptr_f->edge[k];
				// avoid the same edge is use more than once
				go_on=true;
				for(size_t l=0;l<edges_count;l++)
				{
					if(ptr_e->label==edges_on_cell[l])
					{
						go_on=false;
						break;
					}
				}
				// if not, add this edge to the vector and build the face
				if(go_on)
				{
					edges_on_cell[edges_count++]=ptr_e->label;
					f.nb_cell=2;
					f.nb_edge=0;
					f.nb_vertex=0;
					// link this dual face to the mother primal cell
					f.dualCell=ptr_c;
					// each dual face has exactly 3 vertices
					f.vertex[f.nb_vertex++]=&_vertex[ptr_e->firstVertex->label-1]; // first vertex of the edge
					f.vertex[f.nb_vertex++]=&_vertex[ptr_e->secondVertex->label-1]; // second vertex of the edge
					f.vertex[f.nb_vertex++]=&_vertex[nb_primal_vertex+i]; // vertex associated to the centroid of the cell
					// each dual face has exactly 3 edges
					f.edge[f.nb_edge++]=&_edge[ptr_e->label-1]; // common edge to the primal cell
					// find the other two edges associated to this face
					for(size_t l=0;l<ptr_c->nb_vertex;l++)
					{
						ptr_v=ptr_c->vertex[l];
						if(ptr_v->label==ptr_e->firstVertex->label || ptr_v->label==ptr_e->secondVertex->label)
							f.edge[f.nb_edge++]=&_edge[nb_primal_edge+count+l];
					}
					if(f.nb_edge!=3)
					{
						cout<<"Error: I could not find all the edges! Line: "<<__LINE__<<". File: "<<__FILE__<<endl;
						exit(0);
					}
					f.label=index+1;
					f.code=ptr_c->code;
					_face[index]=f;
					index++;
					//// associate this face to the dual cells and the dual cells to this face
					for(size_t m=0;m<ptr_c->nb_face;m++)
					{
						ptr_ff=ptr_c->face[m];
						for(size_t n=0;n<ptr_ff->nb_edge;n++)
						{
							ptr_ee=ptr_ff->edge[n];
							if(ptr_ee->label==ptr_e->label)
							{
								if(ptr_ff->leftCell->label==ptr_c->label)
								{
									_cell[(ptr_ff->label-1)*2].face[_cell[(ptr_ff->label-1)*2].nb_face++]=&_face[index-1];
									if(_face[index-1].leftCell)
										_face[index-1].rightCell=&_cell[(ptr_ff->label-1)*2];
									else
										_face[index-1].leftCell=&_cell[(ptr_ff->label-1)*2];
									break;
								}
								else if(ptr_ff->rightCell->label==ptr_c->label)
								{
									_cell[(ptr_ff->label-1)*2+1].face[_cell[(ptr_ff->label-1)*2+1].nb_face++]=&_face[index-1];
									if(_face[index-1].leftCell)
										_face[index-1].rightCell=&_cell[(ptr_ff->label-1)*2+1];
									else
										_face[index-1].leftCell=&_cell[(ptr_ff->label-1)*2+1];
									break;
								}
								else
								{
									cout<<"Error: I could not find the cell associated to this edge! Line: "<<__LINE__<<". File: "<<__FILE__<<endl;
									exit(0);
								}
							}
						}
						// leave the loop if the two cells which share this face were already found
						if(_face[index-1].rightCell!=NULL && _face[index-1].leftCell!=NULL)
							break;
					}
				}
			}
		}
		count+=ptr_c->nb_vertex;
	}
	// dual faces from primal boundary faces
	count=0;
	for(size_t i=0;i<nb_primal_boundary_face;i++)
	{
		ptr_f=pm.getBoundaryFace(i);
		for(size_t j=0;j<ptr_f->nb_edge;j++)
		{
			ptr_e=ptr_f->edge[j];
			f.nb_cell=1;
			f.nb_edge=0;
			f.nb_vertex=0;
			// link this dual face to the mother primal cell
			f.dualCell=ptr_f->leftCell;
			// each dual face has exactly 3 vertices
			f.vertex[f.nb_vertex++]=&_vertex[ptr_e->firstVertex->label-1]; // first vertex of the edge
			f.vertex[f.nb_vertex++]=&_vertex[ptr_e->secondVertex->label-1]; // second vertex of the edge
			f.vertex[f.nb_vertex++]=&_vertex[nb_primal_vertex+nb_primal_cell+i]; // vertex associated to the centroid of the cell
			// each dual face has exactly 3 edges
			f.edge[f.nb_edge++]=&_edge[ptr_e->label-1]; // common edge to the primal cell
			// find the other two edges associated to this face
			for(size_t l=0;l<ptr_f->nb_vertex;l++)
			{
				ptr_v=ptr_f->vertex[l];
				if(ptr_v->label==ptr_e->firstVertex->label || ptr_v->label==ptr_e->secondVertex->label)
					f.edge[f.nb_edge++]=&_edge[nb_inner_edges_boundary_faces+count+l];
			}
			if(f.nb_edge!=3)
			{
				cout<<"Error: I could not find all the edges! Line: "<<__LINE__<<". File: "<<__FILE__<<endl;
				exit(0);
			}
			f.label=index+1;
			f.code=ptr_f->code;
			_face[index++]=f;
			_boundary_face[_nb_boundary_face++]=&_face[index-1];
			//// associate this face to the dual cell and the dual cell to this face
			_face[index-1].rightCell=NULL;
			_face[index-1].leftCell=&_cell[(ptr_f->label-1)*2+1];
			_cell[(ptr_f->label-1)*2+1].face[_cell[(ptr_f->label-1)*2+1].nb_face++]=&_face[index-1];
		}
		count+=ptr_f->nb_edge;
	}
	// dual faces from primal faces
	for(size_t i=0;i<nb_primal_face;i++)
	{
		ptr_f=pm.getFace(i);
		f.nb_cell=2;
		f.nb_edge=0;
		f.nb_vertex=0;
		// link this edge to the mother primal cell
		if(ptr_f->rightCell)
			f.rightDualCell=ptr_f->rightCell;
		f.leftDualCell=ptr_f->leftCell;
		// contrary to the other dual faces, this faces have the same vertices as primal faces
		for(size_t j=0;j<ptr_f->nb_vertex;j++)
			f.vertex[f.nb_vertex++]=&_vertex[ptr_f->vertex[j]->label-1];
		// contrary to the other dual faces, this faces have the same edges as primal faces
		for(size_t j=0;j<ptr_f->nb_edge;j++)
			f.edge[f.nb_edge++]=&_edge[ptr_f->edge[j]->label-1];
		f.label=index+1;
		f.code=ptr_f->code;
		_face[index++]=f;
		//// associate this face to the dual cell and the dual cell to this face
		_face[index-1].leftCell=&_cell[(ptr_f->label-1)*2];
		_face[index-1].rightCell=&_cell[(ptr_f->label-1)*2+1];
		_cell[(ptr_f->label-1)*2].face[_cell[(ptr_f->label-1)*2].nb_face++]=&_face[index-1];
		_cell[(ptr_f->label-1)*2+1].face[_cell[(ptr_f->label-1)*2+1].nb_face++]=&_face[index-1];
	}
	//// LINK PRIMAL-DUAL
	// in this case, the faces have rightDualCell and leftDualCell
	// the dualCells is set to null
	// make the link primalFace->rightDualCell and primalFace->leftDualFace for the faces of the primal mesh
	for(size_t i=0;i<nb_primal_face;i++)
	{
		pm.getFace(i)->leftDualCell=&_cell[i*2];
		pm.getFace(i)->rightDualCell=&_cell[i*2+1];
	}
	// the link DUAL-PRIMAL was already perfomed
	// COMPLETE THE DATA
	_nb_vertex=_vertex.size();
	_nb_edge=_edge.size();
	_nb_face=_face.size();
	_nb_cell=_cell.size();
	_dim=3;
	_nb_boundary_face=0;
	for(size_t i=0;i<nb_primal_boundary_face;i++)
	{
		ptr_f=pm.getBoundaryFace(i);
		_nb_boundary_face+=ptr_f->nb_edge;
	}
	//// COMPLETE WITH THE GEOMETRICAL DATA
	FVPoint3D<double> u,vv,w;
	double no;
	// EDGE
	// compute the centroid and length of edge
	for(size_t i=0;i<_nb_edge;i++)
    {
		_edge[i].centroid=(_edge[i].firstVertex->coord+_edge[i].secondVertex->coord)*0.5;
		_edge[i].mass_center=_edge[i].centroid;
		u=_edge[i].firstVertex->coord-_edge[i].secondVertex->coord;
		_edge[i].length=Norm(u);
    }
	// FACE
	// compute geometric stuff for the face
	for(size_t i=0;i<_nb_face;i++)
    {
		// conpute perimeter of the face 
		_face[i].perimeter=0.0;
		for(size_t j=0;j<_face[i].nb_edge;j++)
			_face[i].perimeter+=_face[i].edge[j]->length; 
		// compute the area of the face
		FVPoint3D<double> P,P1,P2,Paux;
		FVPoint3D<double> GPP;
		FVGaussPoint2D G2D;
		// compute the centroid of the face
		_face[i].centroid=0.0;
		for(size_t k=0;k<_face[i].nb_vertex;k++)
			_face[i].centroid+=_face[i].vertex[k]->coord;    
		_face[i].centroid/=_face[i].nb_vertex; 
		// face area and mass center
		_face[i].mass_center=0.0;
		_face[i].area=0.0;
		for(size_t j=0;j<_face[i].nb_edge;j++)
        {
			P1=_face[i].edge[j]->firstVertex->coord;
			P2=_face[i].edge[j]->secondVertex->coord;
			u=P1-_face[i].centroid; 
			vv=P2-_face[i].centroid;
			w=CrossProduct(u,vv);
			_face[i].LocalArea[j]=Norm(w)*0.5;
			_face[i].area+=_face[i].LocalArea[j];
			Paux=0.0;
			for (size_t k=1;k<=G2D.getNbPoint(5);k++)
            { 
				GPP=G2D.getPoint(5,k);
				P=GPP.x*P1+GPP.y*P2+GPP.z*_face[i].centroid;
				Paux+=G2D.getWeight(5,k)*P;
            }        
			_face[i].mass_center+=Paux*_face[i].LocalArea[j]; 
        }  
		_face[i].mass_center/=_face[i].area; 
		// we compute the normal from left to rigth for each sub-triangle  
		_face[i].normal=0.0;
		double surf=0.0;
		for(size_t j=0;j<_face[i].nb_edge;j++)
        {
			u=_face[i].edge[j]->firstVertex->coord-_face[i].centroid;
			vv=_face[i].edge[j]->secondVertex->coord-_face[i].centroid;
			w=CrossProduct(u,vv);
			no=Norm(w);
			w/=no;
			_face[i].LocalNormal[j]=w;
			if(!_face[i].leftCell){cout<<" ERROR no left cell for face "<<i<<endl;fflush(NULL);}
			u=_face[i].centroid-_face[i].leftCell->centroid;          
			if(w*u<0) _face[i].LocalNormal[j]*=-1.;
			surf+=no;
			_face[i].normal+=no*_face[i].LocalNormal[j];          
		}
		_face[i].normal/=surf;
    } 
    // CELL
	// compute geometric stuff for the cell
	for(size_t i=0;i<_nb_cell;i++)
    {
		// compute surface of the cell 
		_cell[i].surface=0.0;
		for(size_t j=0;j<_cell[i].nb_face;j++)
			_cell[i].surface+=_cell[i].face[j]->area; 
		// compute the centroid of the cell
		_cell[i].centroid=0.0;
		for(size_t j=0;j<_cell[i].nb_vertex;j++)
			_cell[i].centroid+=_cell[i].vertex[j]->coord;
		_cell[i].centroid/=_cell[i].nb_vertex; 
		// compute the volume of the cell
		_cell[i].volume=0.;
		_cell[i].mass_center=0.;    
		FVPoint3D<double> P,P1,P2,Paux;
		FVPoint4D<double> GPPP;
		FVGaussPoint3D G3D; 
		double Vaux;
		FVPoint3D<double> u,v,w;
		for(size_t j=0;j<_cell[i].nb_face;j++)
        {    
			_cell[i].cell2face[j]= _cell[i].face[j]->centroid-_cell[i].centroid;    
			for(size_t k=0;k<_cell[i].face[j]->nb_edge;k++)
            {
				P1=_cell[i].face[j]->edge[k]->firstVertex->coord;
				P2=_cell[i].face[j]->edge[k]->secondVertex->coord;
				u=_cell[i].cell2face[j]; 
				v=P1-_cell[i].centroid; 
				w=P2-_cell[i].centroid; 
				Vaux=fabs(Det(u,v,w))/6; 
				_cell[i].volume+=Vaux;
				Paux=0;
				for(size_t r=1;r<=G3D.getNbPoint(5);r++)
                {              
					GPPP=G3D.getPoint(5,r);
					P=GPPP.x*P1+GPPP.y*P2+GPPP.z*_cell[i].face[j]->centroid+GPPP.t*_cell[i].centroid;
					Paux+=G3D.getWeight(5,r)*P;
                }
				_cell[i].mass_center+=Paux*Vaux;    
            }           
        }
		_cell[i].mass_center/=_cell[i].volume;   
    }
}
// end of FVMesh3D.cpp