// ------ FVFastRecons2D.cpp ------
// S. CLAIN 2013/01
#include "FVFastRecons2D.h"


FVFastRecons2D::FVFastRecons2D()  // constructor
{
_AValue=NULL;
_ADerivative=NULL;
_ptr_s=NULL;_Vertex2DVect=NULL;_Edge2DVect=NULL;_Cell2DVect=NULL;
_list_Point=new vector< FVPoint2D<double> >;
_list_Point->resize(0); 
_ref_point=0.;_degree=0; _nb_row=_nb_col=0.;
_reconstruction_type= REC_NULL;
_reconstruction_purpose= REC_NULL;
}
FVFastRecons2D::~FVFastRecons2D() //destructor
{

if(_AValue){delete _AValue; _AValue=NULL;}
if(_ADerivative){delete _ADerivative; _ADerivative=NULL;}
if(_list_Point){ delete _list_Point;_list_Point=NULL;}
}
FVFastRecons2D::FVFastRecons2D(const FVFastRecons2D & frec) // copy constructor
{ 
_ptr_s=frec._ptr_s;
_Vertex2DVect=frec._Vertex2DVect;
_Edge2DVect=frec._Edge2DVect;
_Cell2DVect=frec._Cell2DVect; 
_ref_point=frec._ref_point;  
_degree=frec._degree;
 _nb_row=frec._nb_row;
 _nb_col=frec._nb_col;
 
_reconstruction_type=frec._reconstruction_type;
_reconstruction_purpose=frec._reconstruction_purpose;
if(frec._AValue)  
    {
    _AValue=new FVVect<double>; 
    _AValue->resize(frec._AValue->size()); 
    (*_AValue)=(*frec._AValue);
    }
     else {_AValue=NULL;}
if(frec._ADerivative)  
    {
    _ADerivative=new  FVVect<FVPoint2D<double> >;
    _ADerivative->resize(frec._ADerivative->size());
    (*_ADerivative)=(*frec._ADerivative);
    }
     else {_ADerivative=NULL;}
if(frec._list_Point) 
    {
    _list_Point=new vector< FVPoint2D<double> >;
    _list_Point->resize(frec._list_Point->size());
    (*_list_Point)=(*frec._list_Point);
    }
     else {_list_Point=NULL;}  
   
}
FVFastRecons2D & FVFastRecons2D::operator =(const  FVFastRecons2D &frec) // assigment operator
{
_ptr_s=frec._ptr_s;
_Vertex2DVect=frec._Vertex2DVect; 
_Edge2DVect=frec._Edge2DVect;
_Cell2DVect=frec._Cell2DVect; 
_ref_point=frec._ref_point; 
_degree=frec._degree; 
 _nb_row=frec._nb_row;
 _nb_col=frec._nb_col;
_reconstruction_type=frec._reconstruction_type;
_reconstruction_purpose=frec._reconstruction_purpose;
if(frec._AValue) 
    {
    _AValue=new FVVect<double>; 
    _AValue->resize(frec._AValue->size());
    (*_AValue)=(*frec._AValue);
    }
     else {_AValue=NULL;}
if(frec._ADerivative)  
    {
    _ADerivative=new  FVVect<FVPoint2D<double> >;
    _ADerivative->resize(frec._ADerivative->size());
    (*_ADerivative)=(*frec._ADerivative);
    }
     else {_ADerivative=NULL;}
if(frec._list_Point)
    {
    _list_Point=new vector< FVPoint2D<double> >;
    _list_Point->resize(frec._list_Point->size());
    (*_list_Point)=(*frec._list_Point);
    }
     else {_list_Point=NULL;}    
return *this;
}
void FVFastRecons2D::clean()
{
if(_AValue){delete _AValue; _AValue=NULL;}
if(_ADerivative){delete _ADerivative; _ADerivative=NULL;}
if(_list_Point){ delete _list_Point;_list_Point=NULL;}
_ptr_s=NULL;_Vertex2DVect=NULL;_Edge2DVect=NULL;_Cell2DVect=NULL;
_list_Point->resize(0);
_ref_point=0.;_degree=0; _nb_row=_nb_col=0.;
_reconstruction_type= REC_NULL;
_reconstruction_purpose=REC_NULL;
}
// setStencil     
void FVFastRecons2D::setStencil(FVStencil *ptr_s)
{
_ptr_s=ptr_s;
_ref_point=0;_degree=0;_Ncoef=0;
if(ptr_s->getReferenceType()==FVVERTEX2D) 
    _ref_point=((FVVertex2D *)(_ptr_s->getReferenceGeometry()))->coord;
if(ptr_s->getReferenceType()==FVEDGE2D)   
    _ref_point=((FVEdge2D *)(_ptr_s->getReferenceGeometry()))->centroid;    
if(ptr_s->getReferenceType()==FVCELL2D)   
    _ref_point=((FVCell2D *)(_ptr_s->getReferenceGeometry()))->centroid;  
}
void FVFastRecons2D::setStencil(FVStencil *ptr_s, size_t degree)
{
_ptr_s=ptr_s;
_ref_point=0;   
_degree=degree; 
_Ncoef=((_degree+2)*(_degree+1))/2-1;
if(ptr_s->getReferenceType()==FVVERTEX2D) 
    _ref_point=((FVVertex2D *)(_ptr_s->getReferenceGeometry()))->coord;
if(ptr_s->getReferenceType()==FVEDGE2D)   
    _ref_point=((FVEdge2D *)(_ptr_s->getReferenceGeometry()))->centroid;    
if(ptr_s->getReferenceType()==FVCELL2D)   
    _ref_point=((FVCell2D *)(_ptr_s->getReferenceGeometry()))->centroid;   
}
FVStencil * FVFastRecons2D::getStencil()
{    
if(!_ptr_s) {cout<<"No stencil store"<<endl;fflush(NULL);return(NULL);}    
return(_ptr_s);    
}
void FVFastRecons2D::doMatrix(FVMesh2D  &m) 
{
if(_reconstruction_type==REC_NULL) return; 
FVVect<double> Vertex2DVect(m.getNbVertex());Vertex2DVect=0.;
FVVect<double> Edge2DVect(m.getNbEdge());Edge2DVect=0.;
FVVect<double> Cell2DVect(m.getNbCell());Cell2DVect=0.; 
FVRecons2D rec(_ptr_s,_degree);
rec.setVectorVertex2D(Vertex2DVect);
rec.setVectorEdge2D(Edge2DVect);
rec.setVectorCell2D(Cell2DVect);
rec.setReconstructionType(_reconstruction_type); 
rec.doMatrix();
FVVertex2D *ptr_v;
FVEdge2D *ptr_e;
FVCell2D *ptr_c;
//
if(_reconstruction_type==REC_CONSERVATIVE) 
    { 
     _nb_col=_ptr_s->getNbGeometry()+1;
     _nb_row=_list_Point->size(); 
    }
if(_reconstruction_type==REC_NON_CONSERVATIVE)  
    {          
     _nb_col=_ptr_s->getNbGeometry();
     _nb_row=_list_Point->size();  
    }       
if(!_AValue)  {_AValue=new FVVect<double>; }
if(!_ADerivative)  {_ADerivative=new FVVect<FVPoint2D<double> >; }
bool doValue=false,doDerivative=false;
doValue=(_reconstruction_purpose==REC_VALUE)||(_reconstruction_purpose==REC_VALUE_AND_DERIVATIVE);
doDerivative=(_reconstruction_purpose==REC_DERIVATIVE)||(_reconstruction_purpose==REC_VALUE_AND_DERIVATIVE);
if(doValue) _AValue->resize(_nb_row*_nb_col);
if(doDerivative) _ADerivative->resize(_nb_row*_nb_col);
// do the reconstruction
for(size_t j=0;j<_ptr_s->getNbGeometry();j++)
    {      
    if(_ptr_s->getType(j)==FVVERTEX2D) 
        {
         ptr_v=(FVVertex2D *) _ptr_s->getGeometry(j);  
        Vertex2DVect[ptr_v->label-1]=1;   
        }
    if(_ptr_s->getType(j)==FVEDGE2D) 
        {
         ptr_e=(FVEdge2D *) _ptr_s->getGeometry(j);   
         Edge2DVect[ptr_e->label-1]=1;   
        }        
    if(_ptr_s->getType(j)==FVCELL2D) 
        {
         ptr_c=(FVCell2D *) _ptr_s->getGeometry(j);   
         Cell2DVect[ptr_c->label-1]=1;   
        } 
    rec.computeCoef();    
    for(size_t i=0;i<_nb_row;i++)
        {
        // compute all the coeff for the diferent point 
        if(doValue) (*_AValue)[i*_nb_col+j]=rec.getValue((*_list_Point)[i]);
        if(doDerivative) (*_ADerivative)[i*_nb_col+j]=rec.getDerivative((*_list_Point)[i]);
        }
    
    if(_ptr_s->getType(j)==FVVERTEX2D) 
        {
         ptr_v=(FVVertex2D *) _ptr_s->getGeometry(j);   
         Vertex2DVect[ptr_v->label-1]=0;   
        }
    if(_ptr_s->getType(j)==FVEDGE2D) 
        {
         ptr_e=(FVEdge2D *) _ptr_s->getGeometry(j);   
         Edge2DVect[ptr_e->label-1]=0;   
        }        
    if(_ptr_s->getType(j)==FVCELL2D) 
        {
         ptr_c=(FVCell2D *) _ptr_s->getGeometry(j);   
         Cell2DVect[ptr_c->label-1]=0;   
        }       
    }
// treat the las colum apart
if(_reconstruction_type==REC_CONSERVATIVE) 
    {      
    if(_ptr_s->getReferenceType()==FVVERTEX2D) 
        {
         ptr_v=(FVVertex2D *) _ptr_s->getReferenceGeometry();   
         Vertex2DVect[ptr_v->label-1]=1;   
        }
    if(_ptr_s->getReferenceType()==FVEDGE2D) 
        {
         ptr_e=(FVEdge2D *) _ptr_s->getReferenceGeometry();   
         Edge2DVect[ptr_e->label-1]=1;   
        }        
    if(_ptr_s->getReferenceType()==FVCELL2D) 
        {
         ptr_c=(FVCell2D *) _ptr_s->getReferenceGeometry();   
         Cell2DVect[ptr_c->label-1]=1;   
        }         
    rec.computeCoef();    
    
    for(size_t i=0;i<_nb_row;i++)
        {
        // compute all the coeff for the diferent point 
        if(doValue) (*_AValue)[i*_nb_col+_nb_col-1]=rec.getValue((*_list_Point)[i]);
        if(doDerivative) (*_ADerivative)[i*_nb_col+_nb_col-1]=rec.getDerivative((*_list_Point)[i]);
        } 
    if(_ptr_s->getReferenceType()==FVVERTEX2D) 
        {
         ptr_v=(FVVertex2D *) _ptr_s->getReferenceGeometry();   
         Vertex2DVect[ptr_v->label-1]=0;   
        }
    if(_ptr_s->getReferenceType()==FVEDGE2D) 
        {
         ptr_e=(FVEdge2D *) _ptr_s->getReferenceGeometry();   
         Edge2DVect[ptr_e->label-1]=0;   
        }        
    if(_ptr_s->getReferenceType()==FVCELL2D) 
        {
         ptr_c=(FVCell2D *) _ptr_s->getReferenceGeometry();   
         Cell2DVect[ptr_c->label-1]=0;   
        }        
    }    
}
double FVFastRecons2D::getValue(size_t i)
{
double val=0,aux=0;
FVVertex2D *ptr_v;
FVEdge2D *ptr_e;
FVCell2D *ptr_c;
if((_reconstruction_purpose!=REC_VALUE)&&(_reconstruction_purpose!=REC_VALUE_AND_DERIVATIVE))
    {
     cout<<"reconstruction for value is not valid"<<endl;
     return(0);
    }
for(size_t j=0;j<_ptr_s->getNbGeometry();j++)
    {
    if(_ptr_s->getType(j)==FVVERTEX2D) 
        {
         ptr_v=(FVVertex2D *) _ptr_s->getGeometry(j);   
         aux= (*_Vertex2DVect)[ptr_v->label-1];   
        }
    if(_ptr_s->getType(j)==FVEDGE2D) 
        {
         ptr_e=(FVEdge2D *) _ptr_s->getGeometry(j);   
         aux= (*_Edge2DVect)[ptr_e->label-1];   
        }        
    if(_ptr_s->getType(j)==FVCELL2D) 
        {
         ptr_c=(FVCell2D *) _ptr_s->getGeometry(j);   
         aux= (*_Cell2DVect)[ptr_c->label-1];   
        }       
    val+=(*_AValue)[i*_nb_col+j]*aux;
    }
 if(_reconstruction_type==REC_CONSERVATIVE)    
    {
    if(_ptr_s->getReferenceType()==FVVERTEX2D) 
        {
         ptr_v=(FVVertex2D *) _ptr_s->getReferenceGeometry();    
         aux= (*_Vertex2DVect)[ptr_v->label-1];   
        }
    if(_ptr_s->getReferenceType()==FVEDGE2D) 
        {
         ptr_e=(FVEdge2D *) _ptr_s->getReferenceGeometry();   
         aux= (*_Edge2DVect)[ptr_e->label-1];   
        }        
    if(_ptr_s->getReferenceType()==FVCELL2D) 
        {
         ptr_c=(FVCell2D *) _ptr_s->getReferenceGeometry(); 
         aux= (*_Cell2DVect)[ptr_c->label-1];   
        }        
    val+=(*_AValue)[i*_nb_col+_nb_col-1]*aux;   
    }

return(val);    
} 
FVPoint2D<double> FVFastRecons2D::getDerivative(size_t i)
{
if((_reconstruction_purpose!=REC_DERIVATIVE)&&(_reconstruction_purpose!=REC_VALUE_AND_DERIVATIVE))
    {
     cout<<"reconstruction for derivative is not valid"<<endl;
     return(0);
    }    
FVPoint2D<double> Dval=0;
double aux=0;
FVVertex2D *ptr_v;
FVEdge2D *ptr_e;
FVCell2D *ptr_c;
for(size_t j=0;j<_ptr_s->getNbGeometry();j++)
    {
    if(_ptr_s->getType(j)==FVVERTEX2D) 
        {
         ptr_v=(FVVertex2D *) _ptr_s->getGeometry(j);   
         aux= (*_Vertex2DVect)[ptr_v->label-1];   
        }
    if(_ptr_s->getType(j)==FVEDGE2D) 
        {
         ptr_e=(FVEdge2D *) _ptr_s->getGeometry(j);   
         aux= (*_Edge2DVect)[ptr_e->label-1];   
        }        
    if(_ptr_s->getType(j)==FVCELL2D) 
        {
         ptr_c=(FVCell2D *) _ptr_s->getGeometry(j);   
         aux= (*_Cell2DVect)[ptr_c->label-1];   
        }   
    Dval+= (*_ADerivative)[i*_nb_col+j]*aux;  
    }
if(_reconstruction_type==REC_CONSERVATIVE)    
    {
    if(_ptr_s->getReferenceType()==FVVERTEX2D) 
        {
         ptr_v=(FVVertex2D *) _ptr_s->getReferenceGeometry();    
         aux= (*_Vertex2DVect)[ptr_v->label-1];   
        }
    if(_ptr_s->getReferenceType()==FVEDGE2D) 
        {
         ptr_e=(FVEdge2D *) _ptr_s->getReferenceGeometry();   
         aux= (*_Edge2DVect)[ptr_e->label-1];   
        }        
    if(_ptr_s->getReferenceType()==FVCELL2D) 
        {
         ptr_c=(FVCell2D *) _ptr_s->getReferenceGeometry(); 
         aux= (*_Cell2DVect)[ptr_c->label-1];   
        }        
    Dval+=(*_ADerivative)[i*_nb_col+_nb_col-1]*aux;   
    }    
return(Dval);
}

double FVFastRecons2D::getCoefValue(size_t stencil_index,size_t k)
{
if(!_AValue)   
    {
     cout<<"reconstruction for value is not valid"<<endl;
     return(0);
    }    
if(stencil_index>_nb_col)
     {
     cout<<"index stencil is not valid"<<endl;
     return(0);
    }  
if(k>_nb_row)
     {
     cout<<"index point is not valid"<<endl;
     return(0);
    }     
return((*_AValue)[stencil_index+k*_nb_col]);    
}
FVPoint2D<double> FVFastRecons2D::getCoefDerivative(size_t stencil_index,size_t k)
{
if(!_ADerivative) 
    {
     cout<<"reconstruction for derivative is not valid"<<endl;
     return((FVPoint2D<double>) 0);
    }
if(stencil_index>_nb_col)
     {
     cout<<"index stencil is not valid"<<endl;
     return(0);
    }  
if(k>_nb_row)
     {
     cout<<"index point is not valid"<<endl;
     return(0);
    }     
return((*_ADerivative)[stencil_index+k*_nb_col]);        
}