// FVMesh2D.cpp
// May 2014
#include <iostream>
#include <cstdio>
#include <math.h>
#include "FVMesh2D.h"
#include "Gmsh.h"
// default constructor
FVMesh2D::FVMesh2D()
{
    _nb_vertex=0;_nb_cell=0;_nb_edge=0;_dim=0;_dualMesh=NULL;
}  
// constructor with a xml file
FVMesh2D::FVMesh2D(const char *filename)
{
    FVMesh2D::read(filename);
}
// load th mesh from a xml file
size_t FVMesh2D::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 dim 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!=2)
	{
#ifdef _DEBUGS
		cout<<" dimension do not match:"<<"  in "<<__FILE__<< " line "<<__LINE__<<endl;
#endif 
		return(FVWRONGDIM);
	} 
    // read the mesh       
    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[lengthDATA+1],*ptr;
        char  *thedata,*ptr;
        thedata=(char *)malloc(sizeof(char)*(lengthDATA+1));
        _xml.copy(thedata,lengthDATA,beginDATA);
        thedata[lengthDATA]=0;
		//for(size_t i=0;i<lengthDATA;i++) cout<<thedata[i]; cout<<endl;
        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);   
			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[lengthDATA+1],*ptr;
        char  *thedata,*ptr;
        thedata=(char *)malloc(sizeof(char)*(lengthDATA+1));
        _xml.copy(thedata,lengthDATA,beginDATA);
        thedata[lengthDATA]=0;
		//for(size_t i=0;i<lengthDATA;i++) cout<<thedata[i]; cout<<endl;
        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 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[lengthDATA+1],*ptr;
        char  *thedata,*ptr;
        thedata=(char *)malloc(sizeof(char)*(lengthDATA+1));
        _xml.copy(thedata,lengthDATA,beginDATA);
        thedata[lengthDATA]=0;
		//for(size_t i=0;i<lengthDATA;i++) cout<<thedata[i]; cout<<endl;
        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_edge= strtod(ptr, &ptr);  
            for(size_t i=0;i<_cell[_label-1].nb_edge;i++)
				_cell[_label-1].edge[i]= &(_edge[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    
	FVMesh2D::complete_data();
	return(FVOK);     
}
// compute the geometrical data of the mesh
void FVMesh2D::complete_data()
{
	// initialize the vertices and compute the centroid and length
	for(size_t i=0;i<_nb_edge;i++)
    {
		_edge[i].leftCell=NULL;_edge[i].rightCell=NULL;
		_edge[i].nb_cell=2;
		_edge[i].centroid=(_edge[i].firstVertex->coord+_edge[i].secondVertex->coord)*0.5;
		_edge[i].mass_center=_edge[i].centroid;
		FVPoint2D<double> u;
		u=_edge[i].firstVertex->coord-_edge[i].secondVertex->coord;
		_edge[i].length=Norm(u);
    }
	// contruct the list of cell for each vertex    
	for(size_t i=0;i<_nb_vertex;i++)
    {
		_vertex[i].nb_cell=0;
    }
	for(size_t i=0;i<_nb_cell;i++)
    {
		_cell[i].perimeter=0.;
		_cell[i].area=0.;
		_cell[i].centroid=0.;
		_cell[i].onBoundary=false;
		// conpute centroid  and perimeter for the cell 
		// determine the left and right cell for the vertices
		for(size_t j=0;j<_cell[i].nb_edge;j++)
        { 
			_cell[i].onBoundary|=(!(_cell[i].edge[j]->rightCell));// set true if on boundary    
			_cell[i].perimeter+=_cell[i].edge[j]->length; 
			
			if(!(_edge[_cell[i].edge[j]->label-1].leftCell)  )  
				_edge[_cell[i].edge[j]->label-1].leftCell=&(_cell[i]);
			else
				_edge[_cell[i].edge[j]->label-1].rightCell=&(_cell[i]);  
        }
		// update the list of the vertices->cell pointer   
		pos_v=0;
		for(size_t j=0;j<_cell[i].nb_edge;j++)
        {
			bool _still_exist;   
			_still_exist=false; 
			for(size_t k=0;k<pos_v;k++)
				if(_cell[i].edge[j]->firstVertex==_cell[i].vertex[k])  _still_exist=true;
			if(!_still_exist) {_cell[i].vertex[pos_v]=_cell[i].edge[j]->firstVertex;pos_v++;}
			_still_exist=false;  
			for(size_t k=0;k<pos_v;k++)
				if(_cell[i].edge[j]->secondVertex==_cell[i].vertex[k])  _still_exist=true;
			if(!_still_exist) {_cell[i].vertex[pos_v]=_cell[i].edge[j]->secondVertex;pos_v++;}    
        }
		_cell[i].nb_vertex=pos_v;   
		if(pos_v>NB_CELL_PER_VERTEX_2D)
			cout<<"Alert: number of cells per vertex exceeding. Memory overlap! Line: "<<__LINE__<<". File: "<<__FILE__<<endl; 
		for(size_t j=0;j<_cell[i].nb_vertex;j++)
        {
			size_t pos; 
			pos=_cell[i].vertex[j]->label-1;
			_vertex[pos].cell[_vertex[pos].nb_cell]=&(_cell[i]); 
			_vertex[pos].nb_cell++;
        }   
        // compute the centroid of the cell  
		for(size_t k=0;k<_cell[i].nb_vertex;k++)
			_cell[i].centroid+=_cell[i].vertex[k]->coord; 
		_cell[i].centroid/=_cell[i].nb_vertex;      
		// compute the area
		FVPoint2D<double> P,P1,P2,auxP;
		FVPoint3D<double> GPP;
		FVGaussPoint2D G2D; 
		double Saux;
		_cell[i].mass_center=0.;
		for(size_t j=0;j<_cell[i].nb_edge;j++)
        {
			_cell[i].cell2edge[j]= _cell[i].edge[j]->centroid-_cell[i].centroid;      
			FVPoint2D<double> u,v; 
			P1=_cell[i].edge[j]->firstVertex->coord;
			P2=_cell[i].edge[j]->secondVertex->coord;   
			u=P1-_cell[i].centroid; 
			v=P2-_cell[i].centroid;
			Saux=fabs(Det(u,v))*0.5;
			_cell[i].area+=Saux;
			auxP=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*_cell[i].centroid;
				auxP+=G2D.getWeight(5,k)*P;
            }    
			_cell[i].mass_center+=auxP*Saux;
        }
		_cell[i].mass_center/=_cell[i].area;           
    }
	//  we compute the normal from left to rigth  
	_nb_boundary_edge=0;
	for(size_t i=0;i<_nb_edge;i++) 
    {
		FVPoint2D<double> u,v;  
		double no;
		u=_edge[i].firstVertex->coord-_edge[i].secondVertex->coord;  
		no=Norm(u);
		_edge[i].normal.x=u.y;_edge[i].normal.y=-u.x;
		_edge[i].normal/=no;
		if(!_edge[i].leftCell) 
        {
			cout<<"Critical ERROR: no left cell!!!"<<endl;
			if(!_edge[i].rightCell) cout<<"and righcell"<<endl;fflush(NULL);   
        }
		v=_edge[i].centroid-_edge[i].leftCell->centroid;
		if(_edge[i].normal*v<0) _edge[i].normal*=-1.;  
		if(! (_edge[i].rightCell)) {_boundary_edge.push_back(&(_edge[i]));_nb_boundary_edge++;} 
    }
}
// write the mesh to a xml file
size_t FVMesh2D::write(const char *filename)
{
	if((_nb_cell==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=\"2\"    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<<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<<"         <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_edge;           
		for(size_t j=0;j<_cell[i].nb_edge;j++)        
			out_file<<setw(FVCHAMPINT) << _cell[i].edge[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 FVMesh2D
void FVMesh2D::Gmsh2FVMesh(Gmsh &m)
{   
	if (m.getDim()!=2)
	{
		cout<<" Dimension don't match. The Gmsh mesh should contain only 0D or 1D elements"<<endl;
		return;
	}  
	_nb_vertex=m.getNbNode();
	_name="mesh2D  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].code=0;  // default value
    }   
	_nb_cell=0;_nb_edge=0;   
	_cell.resize(0);_edge.resize(0);
	_cell.reserve(m.getNbElement());_edge.reserve(m.getNbElement()*3);
	FVCell2D cell;FVEdge2D edge;
	size_t pos_dim2=0;
	bool exist;
	//cout<<"step1"<<endl;
	
	for(size_t i=0;i<m.getNbElement();i++)
    {  
		GMElement *ptr_el=m.getElement(i);  
		if  (ptr_el->type_element==GMSH_NODE) 
		{
			pos_dim2++;
			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
		} 
		if  (ptr_el->type_element==GMSH_EDGE)
		{    
			pos_dim2++;
			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
			}   
		}
		if(ptr_el->type_element==GMSH_TRI ||ptr_el->type_element==GMSH_QUAD )  // treat   a triangle or a quadrangle
		{
			for(size_t j=0;j<=ptr_el->type_element;j++)
			{
				FVEdge2D edge;
				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))  
					{
						cell.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
					cell.edge[j]=&_edge[_edge.size()-1]; // and keep the pointer
				}
			}
			// at that stage we have all the edge now chesk the code
			cell.nb_edge=ptr_el->type_element+1;  
			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==2)               
					{_cell[j].setCode2Vertex(ptr_el->code_physical%MINUS_ONE_DIM);}
					if (ptr_el->code_physical/MINUS_ONE_DIM==1)
					{_cell[j].setCode2Edge(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;    // label is no longer GMElement label but  give
				// the label of the cell in the _cell table   
				if (ptr_el->code_physical/MINUS_ONE_DIM==2)               
					cell.setCode2Vertex(ptr_el->code_physical%MINUS_ONE_DIM);
				if (ptr_el->code_physical/MINUS_ONE_DIM==1)
					cell.setCode2Edge(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 edge so save it
				
			}  
		}      
	}   
	_nb_edge=_edge.size();    
	_nb_cell=_cell.size();   
	// all the edge and cell are created    
	// now we treat the code cascade backward
	for(size_t cont=pos_dim2;cont>0;cont--) // we skip of 1 because NEGATIVE NUMBER DOES NOT EXIST with size_t
	{  
		size_t i=cont-1; // astuce to allow deacresing counter   
		GMElement *ptr_el=m.getElement(i); 
		// 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 
		}
        
	}
	//cout<<"step3"<<endl;        
	// complete the data
	/*
	 cout<<"I have made "<< _nb_edge<<"edges  and "<<_nb_cell<<" cells"<<endl;
	 
	 for(size_t i=0;i<_nb_edge;i++)
	 {
	 cout<<" edge no "<< _edge[i].label<<" with "<<_edge[i].nb_vertex<<" vertices "
	 <<_edge[i].firstVertex->label<<" and "<<_edge[i].secondVertex->label<<endl,
	 cout<<"the code is "<< _edge[i].code<<endl;
	 }   
	 cout<<"======================"<<endl;
	 for(size_t i=0;i<_nb_cell;i++)
	 {
	 cout<<" cell no "<< _cell[i].label<<" with "<<_cell[i].nb_edge<<" edges "<<", the code is "<< _cell[i].code<<endl;
	 cout<<" list od the edges"<<endl;
	 for(size_t j=0;j<_cell[i].nb_edge;j++)
	 cout<<" edge no "<<_cell[i].edge[j]->label<<", ";
	 cout<<endl;
	 } 
	 exit(0);
	 */
	FVMesh2D::complete_data();     
}
// make the dual mesh
void FVMesh2D::makeDualMesh(FVMesh2D &pm)
{
	size_t nb_primal_vertex,nb_primal_cell,nb_primal_edge,nb_primal_boundary_edge,index;
	FVVertex2D v;   
	FVEdge2D e;
	FVCell2D c;
	FVCell2D *ptr_c;
	FVEdge2D *ptr_e;  
	FVVertex2D *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_boundary_edge=pm.getNbBoundaryEdge();
	//vector<size_t> PrimalCell2DualEdge(nb_primal_cell);
	// resize vectors
	_vertex.resize(nb_primal_vertex+nb_primal_cell+nb_primal_boundary_edge);
	_cell.resize(nb_primal_edge);
	// get the number of dual edges
	size_t nb_dual_edges=0;
	for(size_t i=0;i<nb_primal_cell;i++)
		nb_dual_edges+=pm.getCell(i)->nb_vertex;
	nb_dual_edges+=2*nb_primal_boundary_edge;
	_edge.resize(nb_dual_edges);
	//// 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 edge
	// the centroids of the boundary edges 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_edge;i++)     
    {
		ptr_e=pm.getBoundaryEdge(i);
		v.label=index+i+1;
		v.coord=ptr_e->centroid;
		v.code=ptr_e->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_edge;i++)     
    {
		_cell[i].label=i+1;
		_cell[i].nb_vertex=0;
		_cell[i].nb_edge=0;
		_cell[i].code=pm.getEdge(i)->code;
		_cell[i].primalEdge=pm.getEdge(i);
    }
	//// DUAL EDGES
	// the dual edges from cell
	index=0;
	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.leftCell=NULL;
			e.rightCell=NULL;
			e.nb_vertex=2;
			e.nb_cell=2;
			_edge[index]=e;
			index++;
        }
	}
	// the dual edges from boundary
	for(size_t i=0;i<nb_primal_boundary_edge;i++)     
    {
		ptr_e=pm.getBoundaryEdge(i);
		// first edge
		e.label=index+1;
		e.firstVertex=&_vertex[nb_primal_vertex+nb_primal_cell+i]; // the centroid of the edge
		e.secondVertex=&_vertex[ptr_e->firstVertex->label-1]; //the vertex on edge
		e.code=ptr_e->code;
		e.nb_vertex=2;
		e.nb_cell=1;   
		_edge[index]=e;
		index++;
		// second edge
		e.label=index+1;
		e.firstVertex=&_vertex[nb_primal_vertex+nb_primal_cell+i];// the centroid
		e.secondVertex=&_vertex[ptr_e->secondVertex->label-1];//the vertex on cell
		e.code=ptr_e->code;
		e.nb_vertex=2;
		e.nb_cell=1; 
		_edge[index]=e;
		index++;
    }
	//%%%% talvez possa ser optimizada esta parte
	//// COMBINE CELL AND EDGE
	// find the edges of each dual cell and the two dual cells associated to each dual edge
	size_t NoEdge=0; 
	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++)
        {
			NoEdge++; // for each primal cell*vertex we have this index in this order   
			for(size_t k=0;k<ptr_c->nb_edge;k++)
            {
				ptr_e=ptr_c->edge[k];
				if((ptr_e->firstVertex->label==ptr_c->vertex[j]->label) || (ptr_e->secondVertex->label==ptr_c->vertex[j]->label))
				{
					_cell[ptr_e->label-1].edge[_cell[ptr_e->label-1].nb_edge++]=&_edge[NoEdge-1];  
					if(_edge[NoEdge-1].leftCell)
						_edge[NoEdge-1].rightCell=&_cell[ptr_e->label-1];
					else
						_edge[NoEdge-1].leftCell=&_cell[ptr_e->label-1]; 
				}
            }
        }
    }
	// the dual edges from boundary
	for(size_t i=0;i<nb_primal_boundary_edge;i++)     
    {
		ptr_e=pm.getBoundaryEdge(i);
		// first edge
		NoEdge++;  
		_cell[ptr_e->label-1].edge[_cell[ptr_e->label-1].nb_edge++]=&_edge[NoEdge-1];  
		_edge[NoEdge-1].leftCell=&_cell[ptr_e->label-1]; 
		_edge[NoEdge-1].rightCell=NULL;     
		// second edge
		NoEdge++;
		_cell[ptr_e->label-1].edge[_cell[ptr_e->label-1].nb_edge++]=&_edge[NoEdge-1];  
		_edge[NoEdge-1].leftCell=&_cell[ptr_e->label-1]; 
		_edge[NoEdge-1].rightCell=NULL;      
    }
	// make the link edge->dualCell for the edges of the primal mesh
	for(size_t i=0;i<nb_primal_edge;i++)
		pm.getEdge(i)->dualCell=&_cell[i];
	// make the link edge->dualCell for the edges of the dual mesh
	FVEdge2D *ptr_eLeft_prim,*ptr_eRight_prim;
	for(size_t i=0;i<_edge.size();i++)
    {   
		if(_edge[i].rightCell)
        {
			ptr_eLeft_prim=pm.getEdge(_edge[i].leftCell->label-1);
			ptr_eRight_prim=pm.getEdge(_edge[i].rightCell->label-1); 
			if(ptr_eLeft_prim->leftCell==ptr_eRight_prim->leftCell)
				_edge[i].dualCell=ptr_eLeft_prim->leftCell;
			else if(ptr_eRight_prim->rightCell && (ptr_eLeft_prim->leftCell==ptr_eRight_prim->rightCell))
                _edge[i].dualCell=ptr_eLeft_prim->leftCell;        
			else if(ptr_eLeft_prim->rightCell && (ptr_eLeft_prim->rightCell==ptr_eRight_prim->leftCell))
                _edge[i].dualCell=ptr_eLeft_prim->rightCell;
			else if(ptr_eLeft_prim->rightCell && (ptr_eLeft_prim->rightCell==ptr_eRight_prim->rightCell))
                _edge[i].dualCell=ptr_eLeft_prim->rightCell;
			else
            {
				cout<<"Alert: no dual edge found! Line: "<<__LINE__<<". File: "<<__FILE__<<endl;fflush(NULL);
				cout<<"cell left left ="  <<ptr_eLeft_prim->leftCell<<endl;   
				cout<<"cell left right =" <<ptr_eLeft_prim->rightCell<<endl;  
				cout<<"cell right left =" <<ptr_eRight_prim->leftCell<<endl;  
				cout<<"cell right right ="<<ptr_eRight_prim->rightCell<<endl;               
            }  
        }
		else
			_edge[i].dualCell=pm.getEdge(_edge[i].leftCell->label-1)->leftCell;    
    }
	/*    
	 for(size_t i=0;i<_edge.size();i++)
	 {
	 cout<<"edge "<<_edge[i].label<<" with first vertex:"<<_edge[i].firstVertex->label
	 <<", second vertex:"<<_edge[i].secondVertex->label<<endl;fflush(NULL);
     cout<<"  -- left cell:"<< _edge[i].leftCell->label;
     if(_edge[i].rightCell)  
	 {cout<<", right cell:"<< _edge[i].rightCell->label<<endl;fflush(NULL);}
	 else 
	 {cout<<", right cell is NULL"<<endl;fflush(NULL);}
	 }  
	 */    
	/*
	 for(size_t i=0;i<_cell.size();i++)
	 {
     cout<<"cell "<<_cell[i].label<<" with "<<_cell[i].nb_edge<<" edges"<<endl;   
     for(size_t j=0;j<_cell[i].nb_edge;j++)
	 cout<<" -- "<<_cell[i].edge[j]->label;
     cout<<endl;fflush(NULL);
	 }
	 */
	//  complete the dual cell seting the vertex->cell and cell->vertex    
	bool exist_vertex,exist_cell;
	for(size_t i=0;i<_cell.size();i++)     
    {    
		ptr_c=&_cell[i];
		for(size_t j=0;j<ptr_c->nb_edge;j++)
		{
			// first vertex
			ptr_v=ptr_c->edge[j]->firstVertex;
			// vexter->cell
			exist_vertex=false;
			for(size_t k=0;k<ptr_c->nb_vertex;k++)
			{
				if(ptr_v==ptr_c->vertex[k])
				{
					exist_vertex=true;
					break;
				}
			}
			if(!exist_vertex)
				ptr_c->vertex[ptr_c->nb_vertex++]=ptr_v;
			// cell->vertex
			exist_cell=false;
			for(size_t k=0;k<ptr_v->nb_cell;k++)
			{
				if(ptr_c==ptr_v->cell[k])
				{
					exist_cell=true;
					break;
				}
			}
			if(!exist_cell)
				ptr_v->cell[ptr_v->nb_cell++]=ptr_c;   
			// second vertex
			ptr_v=ptr_c->edge[j]->secondVertex;      
			// vertex->cell
			exist_vertex=false;
			for(size_t k=0;k<ptr_c->nb_vertex;k++)
			{
				if(ptr_v==ptr_c->vertex[k])
				{
					exist_vertex=true;
					break;
				}
			}
			if(!exist_vertex)
				ptr_c->vertex[ ptr_c->nb_vertex++]=ptr_v;
			// cell->vertex
			exist_cell=false;
			for(size_t k=0;k<ptr_v->nb_cell;k++)
			{
				if(ptr_c==ptr_v->cell[k])
				{
					exist_cell=true;
					break;
				}
			}
			if(!exist_cell)
				ptr_v->cell[ptr_v->nb_cell++]=ptr_c;
		}
		
    }
	/*    
	 for(size_t i=0;i<_cell.size();i++)
	 {
     cout<<"cell "<<_cell[i].label<<" with "<<_cell[i].nb_vertex<<" vertices"<<endl;   
     for(size_t j=0;j<_cell[i].nb_vertex;j++)
	 cout<<" -- "<<_cell[i].vertex[j]->label;
     cout<<endl;fflush(NULL);
	 }   
	 for(size_t i=0;i<_vertex.size();i++)
	 {
     cout<<"vertex "<<_vertex[i].label<<" with "<<_vertex[i].nb_cell<<" cells"<<endl;   
     for(size_t j=0;j<_vertex[i].nb_cell;j++)
	 cout<<" -- "<<_vertex[i].cell[j]->label;
     cout<<endl;fflush(NULL);
	 }   
	 */
	// COMPLETE THE DATA    
	_nb_edge=_edge.size();
	_nb_cell=_cell.size();
	for(size_t i=0;i<_nb_edge;i++)
    {
		_edge[i].centroid=(_edge[i].firstVertex->coord+_edge[i].secondVertex->coord)*0.5;
		_edge[i].length=Norm(_edge[i].firstVertex->coord-_edge[i].secondVertex->coord);
    }
	for(size_t i=0;i<_nb_cell;i++)
    {
		_cell[i].perimeter=0.0;
		_cell[i].area=0.0;
		_cell[i].centroid=0.0;
		_cell[i].onBoundary=false;
		for(size_t j=0;j<_cell[i].nb_edge;j++)
        { 
			_cell[i].onBoundary|=(!(_cell[i].edge[j]->rightCell)); // set true if on boundary  
			_cell[i].perimeter+=_cell[i].edge[j]->length; 
        }    
		for(size_t k=0;k<_cell[i].nb_vertex;k++) 
            _cell[i].centroid+=_cell[i].vertex[k]->coord; 
		_cell[i].centroid/=_cell[i].nb_vertex; 
		// compute mass center and the surface
		FVPoint2D<double> P,P1,P2,auxP;
		FVPoint3D<double> GPP;
		FVGaussPoint2D G2D; 
		double Saux;
		_cell[i].mass_center=0.0;
		for(size_t j=0;j<_cell[i].nb_edge;j++)
        {
			_cell[i].cell2edge[j]=_cell[i].edge[j]->centroid-_cell[i].centroid;      
			FVPoint2D<double> u,v; 
			P1=_cell[i].edge[j]->firstVertex->coord;
			P2=_cell[i].edge[j]->secondVertex->coord;   
			u=P1-_cell[i].centroid; 
			v=P2-_cell[i].centroid;
			Saux=fabs(Det(u,v))*0.5;
			_cell[i].area+=Saux;
			auxP=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*_cell[i].centroid;
				auxP+=G2D.getWeight(5,k)*P;
            }    
			_cell[i].mass_center+=auxP*Saux;
        }
		_cell[i].mass_center/=_cell[i].area;         
    }
	_nb_boundary_edge=0;
	_boundary_edge.resize(nb_primal_boundary_edge*2);
	for(size_t i=0;i<_nb_edge;i++) 
    {
		FVPoint2D<double> u,v;  
		double no;
		u=_edge[i].firstVertex->coord-_edge[i].secondVertex->coord;
		no=Norm(u);
		_edge[i].normal.x=u.y;
		_edge[i].normal.y=-u.x;
		_edge[i].normal/=no;
		v=_edge[i].centroid-_edge[i].leftCell->centroid;
		if(_edge[i].normal*v<0)
			_edge[i].normal*=-1.0; 
		if(!(_edge[i].rightCell))
		{
			_boundary_edge[_nb_boundary_edge]=&(_edge[i]);
			_nb_boundary_edge++;
		} 
	}     
}
// remove a cell from the mesh
void FVMesh2D::removeCell(size_t index)
{
	FVCell2D *ptr_c,*ptr_cc;
	FVEdge2D* ptr_e;
	FVVertex2D* ptr_v;
	size_t label;
	ptr_c=&_cell[index];
	//// find the vertices to be removed
	for(size_t i=0;i<ptr_c->nb_vertex;i++)
	{
		ptr_v=ptr_c->vertex[i];
		label=ptr_v->label;
		// there is more than one cell sharing this vertex, we cannot eliminate it
		if(ptr_v->nb_cell>1)
		{
			for(size_t j=0;j<ptr_v->nb_cell;j++)
			{
				if(ptr_v->cell[j]->label==ptr_c->label)
				{
					for(size_t k=j;k<ptr_v->nb_cell-1;k++)
						ptr_v->cell[k]=&_cell[ptr_v->cell[k+1]->label-1];
					ptr_v->cell[ptr_v->nb_cell-1]=NULL;
					ptr_v->nb_cell--;
					break;
				}
			}
		}
		// there is only one cell with this vertex, we can eliminate it
		else
		{
			//// reallocate the pointers before delete
			// for edges
			for(size_t k=0;k<_nb_edge;k++)
			{
				ptr_e=&_edge[k];
				if(ptr_e->firstVertex->label>ptr_v->label && ptr_e->firstVertex->label!=1)
					ptr_e->firstVertex=&_vertex[ptr_e->firstVertex->label-2];
				if(ptr_e->secondVertex->label>ptr_v->label && ptr_e->secondVertex->label!=1)
					ptr_e->secondVertex=&_vertex[ptr_e->secondVertex->label-2];
			}
			// for cells
			for(size_t k=0;k<_nb_cell;k++)
			{
				ptr_cc=&_cell[k];
				for(size_t l=0;l<ptr_cc->nb_vertex;l++)
				{
					if(ptr_cc->vertex[l]->label>label)
						ptr_cc->vertex[l]=&_vertex[ptr_cc->vertex[l]->label-2];
				}
			}
			//// delete the vertex
			_vertex.erase(_vertex.begin()+label-1);
			_nb_vertex--;
			for(size_t j=label-1;j<_nb_vertex;j++)
				_vertex[j].label-=1;
		}
	}
	//// find the edges to be removed
	for(size_t i=0;i<ptr_c->nb_edge;i++)
	{
		ptr_e=ptr_c->edge[i];
		label=ptr_e->label;
		// the edge only has one cell, we can eliminate it
		if(!ptr_e->rightCell)
		{
			//// reallocate the pointers before delete
			for(size_t k=0;k<_nb_cell;k++)
			{
				ptr_cc=&_cell[k];
				for(size_t l=0;l<ptr_cc->nb_edge;l++)
				{
					if(ptr_cc->edge[l]->label!=1 && ptr_cc->edge[l]->label>label)
						ptr_cc->edge[l]=&_edge[ptr_cc->edge[l]->label-2];
				}
			}
			// delete
			_edge.erase(_edge.begin()+label-1);
			_nb_edge--;
			for(size_t j=label-1;j<_nb_edge;j++)
				_edge[j].label-=1;
		}
		// the edge has two cells, we cannot eliminate it
		else if(ptr_e->rightCell->label==ptr_c->label)
			ptr_e->rightCell=NULL;
		else if(ptr_e->leftCell->label==ptr_c->label)
		{
			ptr_e->leftCell=ptr_e->rightCell;
			ptr_e->rightCell=NULL;
		}
	}
	//// reallocate the pointers before delete
	for(size_t i=0;i<_nb_vertex;i++)
	{
		ptr_v=&_vertex[i];
		for(size_t j=0;j<ptr_v->nb_cell;j++)
		{
			if(ptr_v->cell[j]->label-1>index)
				ptr_v->cell[j]=&_cell[ptr_v->cell[j]->label-2];
		}
	}
	for(size_t i=0;i<_nb_edge;i++)
	{
		ptr_e=&_edge[i];
		if(ptr_e->leftCell->label-1>index)
			ptr_e->leftCell=&_cell[ptr_e->leftCell->label-2];
		if(ptr_e->rightCell && ptr_e->rightCell->label-1>index)
			ptr_e->rightCell=&_cell[ptr_e->rightCell->label-2];
	}
	//// remove the cell
	_cell.erase(_cell.begin()+index);
	_nb_cell--;
	for(size_t j=index;j<_nb_cell;j++)
		_cell[j].label-=1;
}
// end of FVMesh.cpp