/*====================================================================================
! Purpose(s): 
!
! Contain(s): 
!
! Author(s):  Javier Lopez Lara (jav.lopez@unican.es)
!             Gabriel Barajas   (barajasg@unican.es)
!====================================================================================*/

#include "CArrayMasks.h"
#include <assert.h>
#include <stdio.h>



static double *Buffer=NULL;
static int BufferSize;

TOffsets Offsets;

//----------------------------------------------------------------------------------------
//----------------------------------------------------------------------------------------
//----------------------------------------------------------------------------------------

// Functions that calculate derivatives of masks.
// ---------------------------------------------

double R_Dx(const TRMask &AM, int AValidNodes)
{
  return AM.Mesh->DUDXr(AM.l(), AM.c(), AM.r(), AValidNodes);
}

//----------------------------------------------------------------------------------------

double R_Dx(const TCMask &AM, int AValidNodes) 
{
  return AM.Mesh->DVDXtr(AM.c(), AM.r(), AValidNodes);
};

//----------------------------------------------------------------------------------------

double T_Dx(const TTMask &AM, int AValidNodes)
{
  return AM.Mesh->DVDXc(AM.l(), AM.l(), AM.c(), AM.c(), AM.r(), AM.r(), AValidNodes);
}

//----------------------------------------------------------------------------------------

double T_Dx(const TTRMask &AM, int AValidNodes)
{
  return AM.Mesh->Dx_T(AM, AValidNodes);
}

//----------------------------------------------------------------------------------------

double C_Dx(const TCMask &AM, int AValidNodes)
{
  return AM.Mesh->DVDXc(AM.l(), AM.l(), AM.c(), AM.c(), AM.r(), AM.r(), AValidNodes);
}

//----------------------------------------------------------------------------------------

double C_Dx(const TCMask &AM, double AAlfa, const TVector &AVector, int AValidNodes)
{
  return AM.Mesh->C_Dx(AM, AAlfa, AVector, AValidNodes);
}

//----------------------------------------------------------------------------------------

double C_Dx(const TLMask &AM, int AValidNodes)
{
  return AM.Mesh->Dx_C(AM, AValidNodes);
}

//----------------------------------------------------------------------------------------

TIntCMask CM_Dx(const TRMask &AM, int AValidNodes)
{
  return AM.Mesh->CM_Dx(AM, AValidNodes);
}

//----------------------------------------------------------------------------------------

TIntBLMask BLM_Dx(const TTMask &AM, int AValidNodes)
{
  return TIntBLMask(AM.Mesh->DVDXtr(AM.c(), AM.r(), AValidNodes));
}

//----------------------------------------------------------------------------------------

TIntBLMask BLM_Dx(const TRMask &AM, int AValidNodes)
{
  return TIntBLMask(AM.Mesh->DUDXtr(AM.tl(), AM.l(), AM.t(), AM.c(), AM.tr(), AM.r(), AValidNodes));
}

//----------------------------------------------------------------------------------------

TIntLMask LM_Dx(const TCMask &AM, int AValidNodes)
{
  int provValidNodes=(((AValidNodes & T_B)!=0)?C_B:0) | (((AValidNodes & TR_B)!=0)?R_B:0);
  return TIntLMask(AM.Mesh->DVDXtr(AM.t(), AM.tr(), provValidNodes));
}

//----------------------------------------------------------------------------------------

TIntTMask TM_Dx(const TTRMask &AM, int AValidNodes)
{
  return (AM.Mesh->TM_Dx(AM, AValidNodes));
}

//----------------------------------------------------------------------------------------

TIntRMask RM_Dy(const TTRMask &AM, int AValidNodes)
{
  return (AM.Mesh->RM_Dy(AM, AValidNodes));
}

//----------------------------------------------------------------------------------------

double T_Dy(const TTMask &AM, int AValidNodes)
{
  return AM.Mesh->DVDYt(AM.b(), AM.c(), AM.t(), AValidNodes);
}

//----------------------------------------------------------------------------------------

double T_Dy(const TCMask &AM, int AValidNodes)
{
  return AM.Mesh->DUDYtr(AM.c(), AM.t(), AValidNodes);
}

//----------------------------------------------------------------------------------------

double R_Dy(const TRMask &AM, int AValidNodes)
{
  return AM.Mesh->DUDYc(AM.b(), AM.b(), AM.c(), AM.c(), AM.t(), AM.t(), AValidNodes);
}

//----------------------------------------------------------------------------------------

double R_Dy(const TTRMask &AM, int AValidNodes)
{
  return AM.Mesh->Dy_R(AM, AValidNodes);
}

//----------------------------------------------------------------------------------------

double C_Dy(const TCMask &AM, int AValidNodes)
{
  return AM.Mesh->DUDYc(AM.b(), AM.b(), AM.c(), AM.c(), AM.t(), AM.t(), AValidNodes);
}

//----------------------------------------------------------------------------------------

double C_Dy(const TCMask &AM, double AAlfa, const TVector &AVector, int AValidNodes)
{
  return AM.Mesh->C_Dy(AM, AAlfa, AVector, AValidNodes);
}

//----------------------------------------------------------------------------------------

double B_Dy(const TBMask &AM, double AAlfa, const TVector &AVector, int AValidNodes)
{
  return AM.Mesh->Dy_C(AM, AValidNodes);
}

//----------------------------------------------------------------------------------------

double C_Dy(const TBMask &AM, int AValidNodes)
{
  return AM.Mesh->Dy_C(AM, AValidNodes);
}

//----------------------------------------------------------------------------------------

TIntCMask  CM_Dy(const TTMask &AM, int AValidNodes)
{
  return AM.Mesh->CM_Dy(AM, AValidNodes);
}

//----------------------------------------------------------------------------------------

TIntBLMask BLM_Dy(const TRMask &AM, int AValidNodes)
{
  return TIntBLMask(AM.Mesh->DUDYtr(AM.c(), AM.t(), AValidNodes));
}

//----------------------------------------------------------------------------------------

TIntBLMask BLM_Dy(const TTMask &AM, int AValidNodes)
{
  return TIntBLMask(AM.Mesh->DVDYtr(AM.tr(), AM.t(), AM.r(), AM.c(), AM.br(), AM.b(), AValidNodes));
}

//----------------------------------------------------------------------------------------

TIntBMask  BM_Dy(const TCMask &AM, int AValidNodes)
{
  int provValidNodes=(((AValidNodes & R_B)!=0)?C_B:0) | (((AValidNodes & TR_B)!=0)?T_B:0);
  return TIntBMask(AM.Mesh->DUDYtr(AM.r(), AM.tr(), provValidNodes));
}


//----------------------------------------------------------------------------------------

double C_Div(const TLMask &AMx, const TBMask &AMy, int AValidNodes) 
{
  return C_Dx(AMx, AValidNodes)+C_Dy(AMy, AValidNodes);
}


//----------------------------------------------------------------------------------------

double C_Div(const TCMask &AMx, const TCMask &AMy, int AValidNodes) 
{
  return C_Dx(AMx, AValidNodes)+C_Dy(AMy, AValidNodes);
};

//----------------------------------------------------------------------------------------

//----------------------------------------------------------------------------------------
//----------------------------------------------------------------------------------------
//----------------------------------------------------------------------------------------


TMask::TMask(TMesh &AMesh, double *AArray)
// Desc: Creates a mask of doubles associated to AArray and adds it to the list of masks of AMesh.
{
  Array=AArray;
  AMesh.AddMask(this);
  Mesh=&AMesh;
};

//----------------------------------------------------------------------------------------

TMask::TMask(TMesh &AMesh)
// Desc: Creates a mask of doubles associated to an empty buffer and adds it to the list of masks of AMesh.
{
  Array=AMesh.NewBuffer();
  AMesh.AddMask(this);
  Mesh=&AMesh;
}

//----------------------------------------------------------------------------------------

inline TIntermedMask TMask::Mul(const TMask &AOp2)
{
	TIntermedMask Res(tr()*AOp2.tr());

  return Res;
};

//----------------------------------------------------------------------------------------

inline TIntermedMask TMask::Div(const TMask &AOp2)
{
	TIntermedMask Res(tr()/AOp2.tr());

  return Res;
};

//----------------------------------------------------------------------------------------

inline TIntermedMask TMask::Add(const TMask &AOp2)
{
	TIntermedMask Res(tr()+AOp2.tr());

  return Res;
};

//----------------------------------------------------------------------------------------


inline TIntermedMask TMask::Sub(const TMask &AOp2)
{
	TIntermedMask Res(tr()-AOp2.tr());

  return Res;
};

//----------------------------------------------------------------------------------------


inline TIntermedMask TMask::Mul(const TIntermedMask &AOp2)
{
	TIntermedMask Res(tr()*AOp2.tr);

  return Res;
};

//----------------------------------------------------------------------------------------

inline TIntermedMask TMask::Div(const TIntermedMask &AOp2)
{
	TIntermedMask Res(tr()/AOp2.tr);

  return Res;
};

//----------------------------------------------------------------------------------------

inline TIntermedMask TMask::Add(const TIntermedMask &AOp2)
{
	TIntermedMask Res(tr()+AOp2.tr);

  return Res;
};

//----------------------------------------------------------------------------------------


inline TIntermedMask TMask::Sub(const TIntermedMask &AOp2)
{
	TIntermedMask Res(tr()-AOp2.tr);

  return Res;
};

//----------------------------------------------------------------------------------------

inline void TMask::Assign(const TMask &AOp2)
{
  *(C+OF.tr)=AOp2.tr();
};

//----------------------------------------------------------------------------------------

inline void TMask::Assign(const TIntermedMask &AOp2)
{
  *(C+OF.tr)=AOp2.tr;
};

 //----------------------------------------------------------------------------------------

inline TIntermedMask TMask::Mul(const double AOp2)
{
	TIntermedMask Res(tr()*AOp2);

  return Res;
};

//----------------------------------------------------------------------------------------

inline TIntermedMask TMask::Div(const double AOp2)
{
	TIntermedMask Res(tr()/AOp2);

  return Res;
};

//----------------------------------------------------------------------------------------

inline TIntermedMask TMask::Add(const double AOp2)
{
	TIntermedMask Res(tr()+AOp2);

  return Res;
};

//----------------------------------------------------------------------------------------


inline TIntermedMask TMask::Sub(const double AOp2)
{
	TIntermedMask Res(tr()-AOp2);

  return Res;
};

//----------------------------------------------------------------------------------------

inline void TMask::Assign(const double AOp2)
{
  *(C+OF.tr)=AOp2;
};

//----------------------------------------------------------------------------------------


//----------------------------------------------------------------------------------------
//----------------------------------------------------------------------------------------
//----------------------------------------------------------------------------------------

// The compact notation of the following code is intentional.

TIntBLMask TTRMask::operator+(const TBLMask &AM)    {return TIntBLMask(*C+AM.tr());};
TIntBLMask TTRMask::operator-(const TBLMask &AM)    {return TIntBLMask(*C-AM.tr());};
TIntBLMask TTRMask::operator*(const TBLMask &AM)    {return TIntBLMask(*C*AM.tr());};
TIntBLMask TTRMask::operator/(const TBLMask &AM)    {return TIntBLMask(*C/AM.tr());};

TIntBLMask TBLMask::operator+(const TTRMask &AM)    {return TIntBLMask(*(C+OF.tr)+AM.c());};
TIntBLMask TBLMask::operator-(const TTRMask &AM)    {return TIntBLMask(*(C+OF.tr)-AM.c());};
TIntBLMask TBLMask::operator*(const TTRMask &AM)    {return TIntBLMask(*(C+OF.tr)*AM.c());};
TIntBLMask TBLMask::operator/(const TTRMask &AM)    {return TIntBLMask(*(C+OF.tr)/AM.c());};

TIntLMask TLMask::operator+(const TIntLMask &AM)    {return Add(AM);};
TIntLMask TLMask::operator-(const TIntLMask &AM)    {return Sub(AM);};
TIntLMask TLMask::operator*(const TIntLMask &AM)    {return Mul(AM);};
TIntLMask TLMask::operator/(const TIntLMask &AM)    {return Div(AM);};
void      TLMask::operator=(const TIntLMask &AM)    {Assign(AM);};

TIntLMask TLMask::operator+(const TLMask &AM)       {return Add(AM);};
TIntLMask TLMask::operator-(const TLMask &AM)       {return Sub(AM);};
TIntLMask TLMask::operator*(const TLMask &AM)       {return Mul(AM);};
TIntLMask TLMask::operator/(const TLMask &AM)       {return Div(AM);};
void      TLMask::operator=(const TLMask &AM)       {Assign(AM);};

TIntLMask TLMask::operator+(const double AM)        {return TIntLMask(*(C+OF.tr)+AM);};
TIntLMask TLMask::operator-(const double AM)        {return TIntLMask(*(C+OF.tr)-AM);};
TIntLMask TLMask::operator*(const double AM)        {return TIntLMask(*(C+OF.tr)*AM);};
TIntLMask TLMask::operator/(const double AM)        {return TIntLMask(*(C+OF.tr)/AM);};
void      TLMask::operator=(const double AM)        {*(C+OF.tr)=AM;};

TIntRMask TRMask::operator+(const TIntRMask &AM) {return Add(AM);};
TIntRMask TRMask::operator-(const TIntRMask &AM) {return Sub(AM);};
TIntRMask TRMask::operator*(const TIntRMask &AM) {return Mul(AM);};
TIntRMask TRMask::operator/(const TIntRMask &AM) {return Div(AM);};
void      TRMask::operator=(const TIntRMask &AM) {Assign(AM);};


TIntRMask TRMask::operator+(const TRMask &AM)    {return Add(AM);};
TIntRMask TRMask::operator-(const TRMask &AM)    {return Sub(AM);};
TIntRMask TRMask::operator*(const TRMask &AM)    {return Mul(AM);};
TIntRMask TRMask::operator/(const TRMask &AM)    {return Div(AM);};
void      TRMask::operator=(const TRMask &AM)    {Assign(AM);};
double    TRMask::ITR(void) {return Mesh->ITR(this);};


TIntRMask TRMask::operator+(const double AM)        {return TIntRMask(*(C+OF.tr)+AM);};
TIntRMask TRMask::operator-(const double AM)        {return TIntRMask(*(C+OF.tr)-AM);};
TIntRMask TRMask::operator*(const double AM)        {return TIntRMask(*(C+OF.tr)*AM);};
TIntRMask TRMask::operator/(const double AM)        {return TIntRMask(*(C+OF.tr)/AM);};
void      TRMask::operator=(const double AM)        {*(C+OF.tr)=AM;};

TIntTMask TTMask::operator+(const TIntTMask &AM) {return Add(AM);};
TIntTMask TTMask::operator-(const TIntTMask &AM) {return Sub(AM);};
TIntTMask TTMask::operator*(const TIntTMask &AM) {return Mul(AM);};
TIntTMask TTMask::operator/(const TIntTMask &AM) {return Div(AM);};
void      TTMask::operator=(const TIntTMask &AM) {Assign(AM);};

TIntTMask TTMask::operator+(const TTMask &AM)    {return Add(AM);};
TIntTMask TTMask::operator-(const TTMask &AM)    {return Sub(AM);};
TIntTMask TTMask::operator*(const TTMask &AM)    {return Mul(AM);};
TIntTMask TTMask::operator/(const TTMask &AM)    {return Div(AM);};
void      TTMask::operator=(const TTMask &AM)    {Assign(AM);};
double    TTMask::ITR(void){return Mesh->ITR(this);};

TIntTMask TTMask::operator+(const double AM)        {return TIntTMask(*(C+OF.tr)+AM);};
TIntTMask TTMask::operator-(const double AM)        {return TIntTMask(*(C+OF.tr)-AM);};
TIntTMask TTMask::operator*(const double AM)        {return TIntTMask(*(C+OF.tr)*AM);};
TIntTMask TTMask::operator/(const double AM)        {return TIntTMask(*(C+OF.tr)/AM);};
void      TTMask::operator=(const double AM)        {*(C+OF.tr)=AM;};

TIntBMask TBMask::operator+(const TIntBMask &AM) {return Add(AM);};
TIntBMask TBMask::operator-(const TIntBMask &AM) {return Sub(AM);};
TIntBMask TBMask::operator*(const TIntBMask &AM) {return Mul(AM);};
TIntBMask TBMask::operator/(const TIntBMask &AM) {return Div(AM);};
void      TBMask::operator=(const TIntBMask &AM) {Assign(AM);};

TIntBMask TBMask::operator+(const TBMask &AM)    {return Add(AM);};
TIntBMask TBMask::operator-(const TBMask &AM)    {return Sub(AM);};
TIntBMask TBMask::operator*(const TBMask &AM)    {return Mul(AM);};
TIntBMask TBMask::operator/(const TBMask &AM)    {return Div(AM);};
void      TBMask::operator=(const TBMask &AM)    {Assign(AM);};

TIntBMask TBMask::operator+(const double AM)        {return TIntBMask(*(C+OF.tr)+AM);};
TIntBMask TBMask::operator-(const double AM)        {return TIntBMask(*(C+OF.tr)-AM);};
TIntBMask TBMask::operator*(const double AM)        {return TIntBMask(*(C+OF.tr)*AM);};
TIntBMask TBMask::operator/(const double AM)        {return TIntBMask(*(C+OF.tr)/AM);};
void      TBMask::operator=(const double AM)        {*(C+OF.tr)=AM;};

TIntCMask TCMask::operator+(const TIntCMask &AM) {return Add(AM);};
TIntCMask TCMask::operator-(const TIntCMask &AM) {return Sub(AM);};
TIntCMask TCMask::operator*(const TIntCMask &AM) {return Mul(AM);};
TIntCMask TCMask::operator/(const TIntCMask &AM) {return Div(AM);};
void      TCMask::operator=(const TIntCMask &AM) {Assign(AM);};

TIntCMask TCMask::operator+(const TCMask &AM)    {return Add(AM);};
TIntCMask TCMask::operator-(const TCMask &AM)    {return Sub(AM);};
TIntCMask TCMask::operator*(const TCMask &AM)    {return Mul(AM);};
TIntCMask TCMask::operator/(const TCMask &AM)    {return Div(AM);};
void      TCMask::operator=(const TCMask &AM)    {Assign(AM);};
double    TCMask::IR(void)  {return Mesh->IR(this);};
double    TCMask::IT(void)  {return Mesh->IT(this);};
double    TCMask::ITR(void) {return Mesh->ITR(this);};
double    TCMask::ITR(int AValidNodes) {return Mesh->ITR(this, AValidNodes);};

TIntLMask TCMask::MIL(int AValidNodes) {return Mesh->MIL(this, AValidNodes);};
TIntBMask TCMask::MIB(int AValidNodes) {return Mesh->MIB(this, AValidNodes);};

TIntRMask TCMask::MCopyToR(void) {return TIntRMask(*(C+OF.tr));};
TIntTMask TCMask::MCopyToT(void) {return TIntTMask(*(C+OF.tr));};

TIntCMask TCMask::operator+(const double AM)        {return TIntCMask(*(C+OF.tr)+AM);};
TIntCMask TCMask::operator-(const double AM)        {return TIntCMask(*(C+OF.tr)-AM);};
TIntCMask TCMask::operator*(const double AM)        {return TIntCMask(*(C+OF.tr)*AM);};
TIntCMask TCMask::operator/(const double AM)        {return TIntCMask(*(C+OF.tr)/AM);};
void      TCMask::operator=(const double AM)        {*(C+OF.tr)=AM;};


TIntTRMask TTRMask::operator+(const TIntTRMask &AM) {return Add(AM);};
TIntTRMask TTRMask::operator-(const TIntTRMask &AM) {return Sub(AM);};
TIntTRMask TTRMask::operator*(const TIntTRMask &AM) {return Mul(AM);};
TIntTRMask TTRMask::operator/(const TIntTRMask &AM) {return Div(AM);};
void       TTRMask::operator=(const TIntTRMask &AM) {Assign(AM);};

TIntTRMask TTRMask::operator+(const TTRMask &AM)    {return Add(AM);};
TIntTRMask TTRMask::operator-(const TTRMask &AM)    {return Sub(AM);};
TIntTRMask TTRMask::operator*(const TTRMask &AM)    {return Mul(AM);};
TIntTRMask TTRMask::operator/(const TTRMask &AM)    {return Div(AM);};
void       TTRMask::operator=(const TTRMask &AM)    {Assign(AM);};


TIntTRMask TTRMask::operator+(const double AM)        {return TIntTRMask(*(C+OF.tr)+AM);};
TIntTRMask TTRMask::operator-(const double AM)        {return TIntTRMask(*(C+OF.tr)-AM);};
TIntTRMask TTRMask::operator*(const double AM)        {return TIntTRMask(*(C+OF.tr)*AM);};
TIntTRMask TTRMask::operator/(const double AM)        {return TIntTRMask(*(C+OF.tr)/AM);};
void       TTRMask::operator=(const double AM)        {*(C+OF.tr)=AM;};


TIntBLMask TBLMask::operator+(const TIntBLMask &AM) {return Add(AM);};
TIntBLMask TBLMask::operator-(const TIntBLMask &AM) {return Sub(AM);};
TIntBLMask TBLMask::operator*(const TIntBLMask &AM) {return Mul(AM);};
TIntBLMask TBLMask::operator/(const TIntBLMask &AM) {return Div(AM);};
void       TBLMask::operator=(const TIntBLMask &AM) {Assign(AM);};

TIntBLMask TBLMask::operator+(const TBLMask &AM)    {return Add(AM);};
TIntBLMask TBLMask::operator-(const TBLMask &AM)    {return Sub(AM);};
TIntBLMask TBLMask::operator*(const TBLMask &AM)    {return Mul(AM);};
TIntBLMask TBLMask::operator/(const TBLMask &AM)    {return Div(AM);};
void       TBLMask::operator=(const TBLMask &AM)    {Assign(AM);};

TIntBLMask TBLMask::operator+(const double AM)        {return TIntBLMask(*(C+OF.tr)+AM);};
TIntBLMask TBLMask::operator-(const double AM)        {return TIntBLMask(*(C+OF.tr)-AM);};
TIntBLMask TBLMask::operator*(const double AM)        {return TIntBLMask(*(C+OF.tr)*AM);};
TIntBLMask TBLMask::operator/(const double AM)        {return TIntBLMask(*(C+OF.tr)/AM);};
void       TBLMask::operator=(const double AM)        {*(C+OF.tr)=AM;};



//----------------------------------------------------------------------------------------
//----------------------------------------------------------------------------------------
//----------------------------------------------------------------------------------------

TIntMask::TIntMask(TMesh &AMesh, int *AArray)
// Desc: Creates a mask of ints associated to AArray and adds it to the list of masks of AMesh.
{
  Array=AArray;
  AMesh.AddMask(this);
}


//----------------------------------------------------------------------------------------

TYMask::TYMask(TMesh &AMesh, double *AArray)
// Desc: Creates a Y-mask of doubles associated to AArray and adds it to the list of masks of AMesh.
{
  Array=AArray;
  Height=AMesh.Height;
  AMesh.AddMask(this);
}

//----------------------------------------------------------------------------------------

inline void TYMask::CenterMask(int AJ)
// Sets the position of the Y-Mask at AJ.
{
  C=Array+AJ-1;
  if (AJ<Height) Ot=1;  else Ot=0;
  if (AJ>1)      Ob=-1; else Ob=0;
}

//----------------------------------------------------------------------------------------

TXMask::TXMask(TMesh &AMesh, double *AArray)
// Desc: Creates a X-mask of doubles associated to AArray and adds it to the list of masks of AMesh.
{
  Array=AArray;
  AMesh.AddMask(this);
}

//----------------------------------------------------------------------------------------

inline void TXMask::CenterMask(int AI)
{
  C=Array+AI-1;
};

//----------------------------------------------------------------------------------------
//----------------------------------------------------------------------------------------
//----------------------------------------------------------------------------------------



//----------------------------------------------------------------------------------------
//----------------------------------------------------------------------------------------
//----------------------------------------------------------------------------------------


TMesh::TMesh(int AWidth, int AHeight, double *AXCoefs, double *AYCoefs, double *AXC, double *AYC, double *ADelX, double *ADelY)
// Desc: Creates the mesh object of dimensions AWidth and AHeight. AXCoefs & AYCoefs are two
//       arrays of coefficients needed to calculate derivatives and interpolations. At this
//       moment they are calculated in a Fortran rutine. AX & AY are pointer to the mesh division arrays.
{
  Count=0;
  CountInt=0;
  CountYTCB=0;
  CountXLCR=0;
  BufferCount=0;

  Width=AWidth;
  Height=AHeight;

  IMin=1;
  IMax=Width;
  JMin=1;
  JMax=Height;
  
  XCoefs=AXCoefs;
  YCoefs=AYCoefs;

  FXC=AXC;
  FYC=AYC;
  FDelX=ADelX;
  FDelY=ADelY;

  IncludeMetricInfo=(XCoefs!=NULL && YCoefs!=NULL && FXC!=NULL && FYC!=NULL && FDelX!=NULL && FDelY!=NULL);
  
  // Allocates buffer array.
  if (Buffer==NULL)
  {
  	// It has to have enough space to store MAXNUMCALCMASKS calculated masks.
  	BufferSize=Width*(Height+(MAXNUMCALCMASKS-1)*3);
  	Buffer=(double *)malloc(BufferSize*sizeof(double));
  }
}

//----------------------------------------------------------------------------------------

TMesh::~TMesh(void)
{

}

//----------------------------------------------------------------------------------------

void TMesh::CenterMasks(int AI, int AJ)
// Desc: Sets the position of all masks in the mesh. It also updates the coefficients for
//       computing derivatives and interpolations. 0<AI<=Width, 0<AJ<=Height
{
  I=AI;
  J=AJ;

  int i=I-1;
  int j=J-1;

  IJ=j*Width+i;
  
  // Computes Offset variable. This variable stores the values we have to add to the center value
  // to obtain the surrounding ones.
  // In other words, if centra value is pointed by pointer C, topleft value will be: *(C+Offset.tl).
  // WARNING, the order of the following lines is important.
  if (AJ<Height) Offsets.t=Width;
            else Offsets.t=0;

  if (AJ>1)  Offsets.b=-Width;
        else Offsets.b=0;
        
  if (AI<Width)
  {
  	Offsets.r=1;
  	Offsets.tr=Offsets.t+1;
  	Offsets.br=Offsets.b+1;
  }
  else
  {
  	Offsets.r=0;
  	Offsets.tr=Offsets.t;
  	Offsets.br=Offsets.b;
  }


  if (AI>1)
  {
    Offsets.bl=Offsets.b-1;
    Offsets.l=-1;
    Offsets.tl=Offsets.t-1;
  }
  else
  {
    Offsets.bl=Offsets.b;
    Offsets.l=0;
    Offsets.tl=Offsets.t;
  } //if



  // Sets the position of all the masks.
  int n;
  for (n=0;n<Count;n++)      Masks[n]->CenterMask(IJ);
  for (n=0;n<CountInt;n++)   MasksInt[n]->CenterMask(IJ);
  for (n=0;n<CountYTCB;n++)  MasksYTCB[n]->CenterMask(AJ);
  for (n=0;n<CountXLCR;n++)  MasksXLCR[n]->CenterMask(AI);


  // Computes coefficients for interpolating and derivating magnitudes along the mesh.
  if (IncludeMetricInfo)
  {
    
    YCoef2=YCoefs[j*TOT_COEF+2];
    YCoef3=YCoefs[j*TOT_COEF+3];
    YCoef4=YCoefs[j*TOT_COEF+4];
    YCoef5=YCoefs[j*TOT_COEF+5];

    c_dudyc_b= YCoefs[j*TOT_COEF+1];
    c_dudyc_t= YCoefs[j*TOT_COEF];
    
    
    c_dvdxc_l=    XCoefs[i*TOT_COEF+1];
    c_dvdxc_r=    XCoefs[i*TOT_COEF];

    double XCoef2=XCoefs[i*TOT_COEF+2];
    double XCoef3=XCoefs[i*TOT_COEF+3];
    double XCoef4=XCoefs[i*TOT_COEF+4];
    double XCoef5=XCoefs[i*TOT_COEF+5];

    c_ftr_t= XCoef3*YCoef2;
    c_ftr_c= XCoef3*YCoef3;
    c_ftr_tr=XCoef2*YCoef2;
    c_ftr_r= XCoef2*YCoef3;


    c_dudxtr_t =c_ftr_t* XCoef4;
    c_dudxtr_tr=c_ftr_tr*XCoef5;
    c_dudxtr_c =c_ftr_c *XCoef4;
    c_dudxtr_r =c_ftr_r* XCoef5;

    c_dvdytr_t =c_ftr_t* YCoef4;
    c_dvdytr_tr=c_ftr_tr*YCoef4;
    c_dvdytr_c =c_ftr_c* YCoef5;
    c_dvdytr_r =c_ftr_r* YCoef5;

    c_dudx=XCoef4;
    c_dvdy=YCoef4;

    c_rdudx=XCoef5;
    c_tdvdy=YCoef5;


    XC=FXC[i];
    YC=FYC[j];
    DelX=FDelX[i];
    DelY=FDelY[j];
  } //if


}

//----------------------------------------------------------------------------------------

void TMesh::GoRight(void)
// Desc: Moves all masks in the mesh one step to the right. It also updates the coefficients for
//       computing derivatives and interpolations. This function doesn't check if we go beyond the
//       right limit.
{
  int n;


  for (n=0;n<Count;n++)     Masks[n]->GoRight();
  for (n=0;n<CountXLCR;n++) MasksXLCR[n]->GoRight();
  for (n=0;n<CountInt;n++)  MasksInt[n]->GoRight();

  if (I==1)
  {
  	Offsets.tl=Offsets.t-1;
  	Offsets.l=-1;
  	Offsets.bl=Offsets.b-1;
  }
  if (I==Width-1)
  {
  	Offsets.tr=Offsets.t;
  	Offsets.r=0;
  	Offsets.br=Offsets.b;
  } //if

  I++;
  IJ++;


  
  // Computes coefficients for interpolating and derivating magnitudes along the mesh.
  if (IncludeMetricInfo)
  {
    int i=I-1;
    int j=J-1;

    c_dvdxc_l=    XCoefs[i*TOT_COEF+1];
    c_dvdxc_r=    XCoefs[i*TOT_COEF];

    double XCoef2=XCoefs[i*TOT_COEF+2];
    double XCoef3=XCoefs[i*TOT_COEF+3];
    double XCoef4=XCoefs[i*TOT_COEF+4];
    double XCoef5=XCoefs[i*TOT_COEF+5];

    c_ftr_t= XCoef3*YCoef2;
    c_ftr_c= XCoef3*YCoef3;
    c_ftr_tr=XCoef2*YCoef2;
    c_ftr_r= XCoef2*YCoef3;


    c_dudxtr_t =c_ftr_t *XCoef4;
    c_dudxtr_tr=c_ftr_tr*XCoef5;
    c_dudxtr_c =c_ftr_c *XCoef4;
    c_dudxtr_r =c_ftr_r *XCoef5;

    c_dvdytr_t =c_ftr_t *YCoef4;
    c_dvdytr_tr=c_ftr_tr*YCoef4;
    c_dvdytr_c =c_ftr_c *YCoef5;
    c_dvdytr_r =c_ftr_r *YCoef5;

    c_dudx=XCoef4;
    c_dvdy=YCoef4;

    c_rdudx=XCoef5;
    c_tdvdy=YCoef5;


    XC=FXC[i];
    DelX=FDelX[i];
  } //if

}

//----------------------------------------------------------------------------------------

void TMesh::First(void)
// Desc: Set all masks at the first element of the grid.
{
  CenterMasks(IMin, JMin);
  
}

//----------------------------------------------------------------------------------------

bool TMesh::Next(void)
// Desc: Moves all masks to the next element of the grid. When all elements have been visited, this
//       function returns false.
{
  if (I<IMax)
  {
  	GoRight();
  }
  else
  {
    I=IMin;
    J++;
    if (J==JMax+1) return false;
              else CenterMasks(I, J);
  } //if
  return true;

}

//----------------------------------------------------------------------------------------

void TMesh::CenterMasks2(int AI, int AJ)
// Desc: Sets the position of all masks in the mesh. It also updates the coefficients for
//       computing derivatives and interpolations. 0<AI<=Width, 0<AJ<=Height
{
  I=AI;
  J=AJ;
  IJ=(AJ-1)*Width+AI-1;
  
  // Computes Offset variable. This variable stores the values we have to add to the center value
  // to obtain the surrounding ones.
  // In other words, if center value is pointed by pointer C, topleft value will be: *(C+Offset.tl).
  if (AJ<Height && AJ>0) Offsets.t=Width;
                    else Offsets.t=0;

  if (AJ>1)  Offsets.b=-Width;
        else Offsets.b=0;
        
  if (AI<Width && AI>0) Offsets.r=1;
                   else Offsets.r=0;

  if (AI>1) Offsets.l=-1;   
       else Offsets.l=0;


  Offsets.tl=Offsets.t+Offsets.l;
  Offsets.bl=Offsets.b+Offsets.l;
  
	Offsets.tr=Offsets.t+Offsets.r;
	Offsets.br=Offsets.b+Offsets.r;


  int provI, provJ, provIJ;
  if (AI>0) provI=AI; else provI=1;
  if (AJ>0) provJ=AJ; else provJ=1;
  provIJ=(provJ-1)*Width+provI-1;

  // Sets the position of all the masks.
  int n;
  for (n=0;n<Count;n++)      Masks[n]->CenterMask(provIJ);
  for (n=0;n<CountInt;n++)   MasksInt[n]->CenterMask(provIJ);
  for (n=0;n<CountYTCB;n++)  MasksYTCB[n]->CenterMask(provJ);
  for (n=0;n<CountXLCR;n++)  MasksXLCR[n]->CenterMask(provI);


  // Computes coefficients for interpolating and derivating magnitudes along the mesh.
  if (IncludeMetricInfo)
  {
    int i=provI-1;
    int j=provJ-1;

    YCoef2=YCoefs[j*TOT_COEF+2];
    YCoef3=YCoefs[j*TOT_COEF+3];
    YCoef4=YCoefs[j*TOT_COEF+4];
    YCoef5=YCoefs[j*TOT_COEF+5];

    c_dudyc_b= YCoefs[j*TOT_COEF+1];
    c_dudyc_t= YCoefs[j*TOT_COEF];
    
    
    c_dvdxc_l=    XCoefs[i*TOT_COEF+1];
    c_dvdxc_r=    XCoefs[i*TOT_COEF];

    double XCoef2=XCoefs[i*TOT_COEF+2];
    double XCoef3=XCoefs[i*TOT_COEF+3];
    double XCoef4=XCoefs[i*TOT_COEF+4];
    double XCoef5=XCoefs[i*TOT_COEF+5];

    c_ftr_t= XCoef3*YCoef2;
    c_ftr_c= XCoef3*YCoef3;
    c_ftr_tr=XCoef2*YCoef2;
    c_ftr_r= XCoef2*YCoef3;


    c_dudxtr_t =c_ftr_t* XCoef4;
    c_dudxtr_tr=c_ftr_tr*XCoef5;
    c_dudxtr_c =c_ftr_c *XCoef4;
    c_dudxtr_r =c_ftr_r* XCoef5;

    c_dvdytr_t =c_ftr_t* YCoef4;
    c_dvdytr_tr=c_ftr_tr*YCoef4;
    c_dvdytr_c =c_ftr_c* YCoef5;
    c_dvdytr_r =c_ftr_r* YCoef5;
    
    c_dudx=XCoef4;
    c_dvdy=YCoef4;
    
    c_rdudx=XCoef5;
    c_tdvdy=YCoef5;


    XC=FXC[i];
    YC=FYC[j];
    DelX=FDelX[i];
    DelY=FDelY[j];
  } //if

}

//----------------------------------------------------------------------------------------

void TMesh::GoRight2(void)
// Desc: Moves all masks in the mesh one step to the right. It also updates the coefficients for
//       computing derivatives and interpolations. This function doesn't check if we go beyond the
//       right limit.
{
  int n;


  for (n=0;n<Count;n++)     Masks[n]->GoRight();
  //for (n=0;n<CountXLCR;n++) MasksXLCR[n]->GoRight();
  for (n=0;n<CountInt;n++)  MasksInt[n]->GoRight();

  if (I==1)
  {
  	Offsets.tl=Offsets.t-1;
  	Offsets.l=-1;
  	Offsets.bl=Offsets.b-1;
  }
  if (I==Width-1)
  {
  	Offsets.tr=Offsets.t;
  	Offsets.r=0;
  	Offsets.br=Offsets.b;
  } //if

  I++;
  IJ++;
  
  // Computes coefficients for interpolating and derivating magnitudes along the mesh.
  if (IncludeMetricInfo)
  {
    int i, j;
    if (I>0) i=I-1; else i=0;
    if (J>0) j=J-1; else j=0;

    c_dvdxc_l=    XCoefs[i*TOT_COEF+1];
    c_dvdxc_r=    XCoefs[i*TOT_COEF];

    double XCoef2=XCoefs[i*TOT_COEF+2];
    double XCoef3=XCoefs[i*TOT_COEF+3];
    double XCoef4=XCoefs[i*TOT_COEF+4];
    double XCoef5=XCoefs[i*TOT_COEF+5];

    c_ftr_t= XCoef3*YCoef2;
    c_ftr_c= XCoef3*YCoef3;
    c_ftr_tr=XCoef2*YCoef2;
    c_ftr_r= XCoef2*YCoef3;


    c_dudxtr_t =c_ftr_t *XCoef4;
    c_dudxtr_tr=c_ftr_tr*XCoef5;
    c_dudxtr_c =c_ftr_c *XCoef4;
    c_dudxtr_r =c_ftr_r *XCoef5;

    c_dvdytr_t =c_ftr_t *YCoef4;
    c_dvdytr_tr=c_ftr_tr*YCoef4;
    c_dvdytr_c =c_ftr_c *YCoef5;
    c_dvdytr_r =c_ftr_r *YCoef5;
    
    c_dudx=XCoef4;
    c_dvdy=YCoef4;
    
    c_rdudx=XCoef5;
    c_tdvdy=YCoef5;


    XC=FXC[i];
    DelX=FDelX[i];

  }

}

//----------------------------------------------------------------------------------------

void TMesh::First2(void)
// Desc: Set all masks at the first element of the grid.
{
  BufferCount=0;

  CenterMasks2(0, 0);

  if (CountXLCR>0) printf("XLCRMask not allowed with first2()");
}

//----------------------------------------------------------------------------------------

bool TMesh::Next2(void)
{
  if (I>0 && I<Width)
  {
  	GoRight2();
  }
  else if (I==0)
  {
    I++;
    CenterMasks2(I, J);
  }
  else
  {
    I=0;
    J++;
    if (J==Height+1) return false;
                else CenterMasks2(I, J);
  } //if
  return true;

}

//----------------------------------------------------------------------------------------

double *TMesh::NewBuffer(void)
{
  assert(BufferCount<MAXNUMCALCMASKS);

  double *Res=Buffer+3*Width*BufferCount;

  BufferCount++;


  return Res;             
}

//----------------------------------------------------------------------------------------

double TMesh::ITR(const TCMask *AM, int AValidNodes) 
{
  double Num=(VN(C_B)? (AM->c()* c_ftr_c): 0.0)+
             (VN(T_B)? (AM->t()* c_ftr_t): 0.0)+
             (VN(R_B)? (AM->r()* c_ftr_r): 0.0)+
             (VN(TR_B)?(AM->tr()*c_ftr_tr):0.0);

  double Den=(VN(C_B)? c_ftr_c: 0.0)+
             (VN(T_B)? c_ftr_t: 0.0)+
             (VN(R_B)? c_ftr_r: 0.0)+
             (VN(TR_B)?c_ftr_tr:0.0);

  if (Den>0) return (Num/Den);
        else return 0.0;
}

//----------------------------------------------------------------------------------------

double TMesh::C_Dx(const TCMask &AM, double AAlfa, const TVector &AVector, int AValidNodes)
{
  double AL, AR;
  
  // This coefficients take account of how much importance have the derivatives at one side of the cell and at
  // the opposite one. If AAlfa=0.0, both derivatives have the same importance, and the derivative at the center 
  // is just the mean of both. If AAlfa is 1.0, the one at the side where the flux comes is set at the center.
  if (AVector.x>=0) AL=AAlfa+1.0; else AL=1.0-AAlfa;
  AR=2.0-AL;


  double Res=0.5*(AR*(VN(R_B)?((AM.r()-AM.c())*c_dvdxc_r):0.0)+
                  AL*(VN(L_B)?((AM.c()-AM.l())*c_dvdxc_l):0.0));


  return Res;
}

//----------------------------------------------------------------------------------------

double TMesh::C_Dy(const TCMask &AM, double AAlfa, const TVector &AVector, int AValidNodes)
{
  double AB, AT;
  
  // This coefficients take account of how much importance have the derivatives at one side of the cell and at
  // the opposite one. If AAlfa=0.0, both derivatives have the same importance, and the derivative at the center 
  // is just the mean of both. If AAlfa is 1.0, the one at the side where the flux comes is set at the center.
  if (AVector.y>=0) AB=AAlfa+1.0; else AB=1.0-AAlfa;
  AT=2.0-AB;


  double Res=0.5*(AT*(VN(T_B)?((AM.t()-AM.c())*c_dudyc_t):0.0)+
                  AB*(VN(B_B)?((AM.c()-AM.b())*c_dudyc_b):0.0));


  return Res;
}

//----------------------------------------------------------------------------------------

double TMesh::Dx_C(const TRMask &AM, int AValidNodes)
{

  return (VN(LC_B)?((AM.c()-AM.l())*c_dudx):0.0);

}

//----------------------------------------------------------------------------------------

double TMesh::Dx_C(const TLMask &AM, int AValidNodes)
{

  return (VN(CR_B)?((AM.r()-AM.c())*c_dudx):0.0);

}

//----------------------------------------------------------------------------------------

TIntCMask TMesh::CM_Dx(const TRMask &AM, int AValidNodes)
{
  return (VN(TTR_B)?((AM.tr()-AM.t())*c_rdudx):0.0);
}

//----------------------------------------------------------------------------------------

double TMesh::Dy_C(const TTMask &AM, int AValidNodes)
{

  return (VN(BC_B)?((AM.c()-AM.b())*c_dvdy):0.0);

}

//----------------------------------------------------------------------------------------

double TMesh::Dy_R(const TTRMask &AM, int AValidNodes)
{
  return (VN(BC_B)?((AM.c()-AM.b())*c_dvdy):0.0);
}

//----------------------------------------------------------------------------------------

double TMesh::TMesh::Dy_C(const TBMask &AM, int AValidNodes)
{

  return (VN(CT_B)?((AM.t()-AM.c())*c_dvdy):0.0);

}

//----------------------------------------------------------------------------------------

TIntCMask TMesh::CM_Dy(const TTMask &AM, int AValidNodes)
{
  return (VN(TRR_B)?((AM.tr()-AM.r())*c_tdvdy):0.0);
}

//----------------------------------------------------------------------------------------

TIntTMask TMesh::TM_Dx(const TTRMask &AM, int AValidNodes=ALL_B)
{
  return (VN(TTR_B)?((AM.tr()-AM.t())*c_rdudx):0.0);
}

//----------------------------------------------------------------------------------------

double TMesh::Dx_T(const TTRMask &AM, int AValidNodes=ALL_B)
{
  return (VN(LC_B)?((AM.c()-AM.l())*c_dudx):0.0);
}

//----------------------------------------------------------------------------------------

TIntRMask TMesh::RM_Dy(const TTRMask &AM, int AValidNodes=ALL_B)
{
  return (VN(TRR_B)?((AM.tr()-AM.r())*c_tdvdy):0.0);
}

//----------------------------------------------------------------------------------------

TVector TMesh::GradFc(const TMask &AMask, double AAlfa, const TVector &AV, int AValidNodes)
// Desc: This function returns the gradient of a F-kind variable at the center of the cell. AAlfa parameter
//       must be between [0..1] and controls how much the upwind scheme is used in the interpolation
//       (0 -> not used, 1 -> only upwind). AV is the velocity vector.
{

  TVector Res;
  double AL, AR, AT, AB;
  
  // This coefficients take account of how much importance have the derivatives at one side of the cell and at
  // the opposite one. If AAlfa=0.0, both derivatives have the same importance, and the derivative at the center 
  // is just the mean of both. If AAlfa is 1.0, the one at the side where the flux comes is set at the center.
  if (AV.x>=0) AL=AAlfa+1.0; else AL=1.0-AAlfa;
  if (AV.y>=0) AB=AAlfa+1.0; else AB=1.0-AAlfa;
  AR=2.0-AL;
  AT=2.0-AB;

  Res.x=0.5*(AR*(VN(R_B)?((AMask.r()-AMask.c())*c_dvdxc_r):0.0)+
             AL*(VN(L_B)?((AMask.c()-AMask.l())*c_dvdxc_l):0.0));

  Res.y=0.5*(AT*(VN(T_B)?((AMask.t()-AMask.c())*c_dudyc_t):0.0)+
             AB*(VN(B_B)?((AMask.c()-AMask.b())*c_dudyc_b):0.0));


  return Res;

}

//----------------------------------------------------------------------------------------