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

#ifndef CArrayMasksH
#define CArrayMasksH


// Class
#include <stdlib.h>

#define TOT_COEF 6

#define TL_B  1
#define T_B   2
#define TR_B  4
#define L_B   8
#define C_B   16
#define R_B   32
#define BL_B  64
#define B_B   128
#define BR_B  256
#define CTTRR_B (C_B | T_B | TR_B | R_B)
#define CBBRR_B (C_B | B_B | BR_B | R_B)
#define CTTLL_B (C_B | T_B | TL_B | L_B)
#define CBBLL_B (C_B | B_B | BL_B | L_B)

#define CLTRB_B (C_B | L_B | T_B | R_B | B_B)
#define LTRB_B  (L_B | T_B | R_B | B_B)
#define CLB_B   (C_B | L_B | B_B)
#define LC_B    (C_B | L_B)
#define BC_B    (C_B | B_B)
#define CR_B    (C_B | R_B)
#define CT_B    (C_B | T_B)

#define LR_B   (L_B | R_B)
#define TB_B   (T_B | B_B)
#define TRR_B  (TR_B | R_B)
#define TLL_B  (TL_B | L_B)
#define BLL_B  (BL_B | L_B)
#define TTR_B  (TR_B | T_B)
#define BBR_B  (BR_B | B_B)
#define BBL_B  (BL_B | B_B)
#define BRR_B  (BR_B | R_B)

#define ALL_B (C_B | L_B | T_B | R_B | B_B | BL_B | BR_B | TR_B | TL_B )

#define VN(Nodes) (((AValidNodes) & (Nodes))==(Nodes))

#define OF Offsets

const int MAXNUMCALCMASKS=30;


//----------------------------------------------------------------------------------------
//  Miscelaneous classes and functions.
//----------------------------------------------------------------------------------------

inline unsigned char _Eq(int AOp1, int AOp2)
// returns 1 if AOp1==AOp2 and 0 otherwise.
{
#ifdef _M_IX86aaaaaa
  // If processor is 8x86...
  register unsigned char res;
  asm
  {
    mov ecx, AOp1;
    xor ecx, AOp2;
    sete cl;
    mov res, cl;

  };
  return res;
#else
  return ((AOp1==AOp2)?1:0);
#endif

};

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

inline unsigned char _NEq(int AOp1, int AOp2)
// returns 1 if AOp1==AOp2 and 0 otherwise.
{
#ifdef _M_IX86aaaaaaa
  // If processor is 8x86...
  register unsigned char res;
  asm
  {
    mov eax, AOp1;
    xor eax, AOp2;
    setne al;
    mov res, al;
  };
  return res;
#else
  return ((AOp1!=AOp2)?1:0);
#endif


};

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

// Simple class that let us define a rectangle and test whether a point is inside it or not.
class TIntRect
{
private:
  int L, R, B, T;

public:
  TIntRect(int AL, int AR, int AT, int AB) {if (AL<=AR) {L=AL; R=AR;} else {L=AR; R=AL;}; if (AB<=AT) {T=AT; B=AB;} else {T=AB; B=AT;};};
  TIntRect(int *ARect) {if (ARect[0]<=ARect[1]) {L=ARect[0]; R=ARect[1];} else {L=ARect[1]; R=ARect[0];}; 
                        if (ARect[2]<=ARect[3]) {B=ARect[2]; T=ARect[3];} else {B=ARect[3]; T=ARect[2];};};
  
  void Grow(int AWidth) {L=L-AWidth; R=R+AWidth; T=T+AWidth; B=B-AWidth;};

  bool IsInside(int Ai, int Aj) {return (Ai>=L && Ai<=R && Aj>=B && Aj<=T);};
  bool IsInsideEx(int Ai, int Aj) {return (Ai>L && Ai<R && Aj>B && Aj<T);};

}; //TIntRect

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


// Valid nodes constructors.

inline int ValidNodes(bool AL, bool AR, bool AT, bool AB)
{ return (((AT&&AL)?TL_B:0) | (AT?T_B:0) | ((AT&&AR)?TR_B:0) |
          (AL?L_B:0       ) | (     C_B) | (AR?R_B:0       ) |   
          ((AB&&AL)?BL_B:0) | (AB?B_B:0) | ((AB&&AR)?BR_B:0)); };

inline int ValidNodes(bool ATL, bool AT, bool ATR, bool AL, bool AC, bool AR, bool ABL, bool AB, bool ABR)
{ return ((ATL?TL_B:0) | (AT?T_B:0) | (ATR?TR_B:0) |
          (AL?L_B:0  ) | (AC?C_B:0) | (AR?R_B:0) |
          (ABL?BL_B:0) | (AB?B_B:0) | (ABR?BR_B:0)); };

inline int BLValidNodesFromC(int AValidNodes)
{
  return ((VN(CTTRR_B)?TR_B:0) | (VN(CBBRR_B)?R_B:0) |
          (VN(CTTLL_B)?T_B:0 ) | (VN(CBBLL_B)?C_B:0) |
          (VN(TLL_B)?TL_B:0  ) | (VN(BBR_B)?BR_B:0 ) |
          (VN(BLL_B)?L_B:0   ) | (VN(BBL_B)?B_B:0  ) |
          (VN(BL_B)?BL_B:0));


}

// Forward declarations.
//----------------------------------------------------------------------------------------

class TMesh;

class TIntermedMask;
class TIntLMask;
class TIntRMask;
class TIntTMask;
class TIntBMask;
class TIntCMask;
class TIntTRMask;
class TIntBLMask;


class TLMask;
class TRMask;
class TTMask;
class TBMask;
class TCMask;
class TTRMask;
class TBLMask;

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


// This structure stores the indexes we have to add to the central element in order to obtain the
// surounded ones. Being "Offsets" variable global to the program, this means we cannot neast
// Mesh.First()-Mesh.Next() loops.
struct TOffsets
{
  int tl, t, tr;
  int  l,    r;
  int bl, b, br;
};

extern TOffsets Offsets;


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


//                                        Points of the cell where the different
//     Name of the points in a mask.            type of masks are defined.
//
//                                                       TTMask       TTRMask
//        -------------------                      ---------O--------O
//        |     |     |     |                      |                 |
//        | tl  |  t  |  tr |                      |                 |
//        |     |     |     |                      |                 |
//        ---------T---------                      |                 |
//        |     |     |     |                      |      TCMask     |
//        |  l  L  c  R  r  |               TLMask O        O        O TRMask
//        |     |     |     |                      |                 |
//        ---------B---------                      |                 |
//        |     |     |     |                      |                 |
//        | bl  |  b  |  br |                      |                 |
//        |     |     |     |                      |                 |
//        -------------------                      0--------O---------
//                                          TBRMask       TBMask
//




// Base for TMask class. It stores all the data needed by that class.
class TMaskBase
{
protected:
  double *Array;


public:
  TMesh  *Mesh;
  double *C;

};

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


// This class is used to work with pairs of values. We can operate with them and compute their
// sum, difference or mean.
class TPair
{
private:
  double r, l;

public:
  inline TPair(double AL, double AR) {r=AR;l=AL;};
  inline TPair(double *AL, double *AR) {r=*AR;l=*AL;};


  inline TPair operator*(const TPair &AOp2) {return TPair(l*AOp2.l, r*AOp2.r);};
  inline TPair operator/(const TPair &AOp2) {return TPair(l/AOp2.l, r/AOp2.r);};


  inline double Delta(void) {return r-l;};
  inline double Sum(void) {return r+l;};
  inline double Mean(void) {return 0.5*(r+l);};
};




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

class Tcrt
{
public:
  double t, tr;
  double c, r;


  inline Tcrt(double At, double Atr, double Ac, double Ar) : t(At), tr(Atr), c(Ac), r(Ar) {};
  inline Tcrt(TMaskBase *AMask)
    {double *C=AMask->C;t=*(C+OF.t);tr=*(C+OF.tr);;c=*C;r=*(C+OF.r);};


  inline Tcrt operator*(const Tcrt &AOp2) {return Tcrt(t*AOp2.t, tr*AOp2.tr, c*AOp2.c, r*AOp2.r);};
  inline Tcrt operator/(const Tcrt &AOp2) {return Tcrt(t/AOp2.t, tr/AOp2.tr, c/AOp2.c, r/AOp2.r);};
  inline Tcrt operator*(const double AOp2) {return Tcrt(t*AOp2, tr*AOp2, c*AOp2, r*AOp2);};
  friend inline Tcrt operator*(const double AOp1, const Tcrt &AOp2) {return Tcrt(AOp2.t*AOp1, AOp2.tr*AOp1, AOp2.c*AOp1, AOp2.r*AOp1);};

};


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


class Tbc_R 
{
public:
  double l, c, r;
  double bl, b, br;


  inline Tbc_R(double Al, double Ac, double Ar, double Abl, double Ab, double Abr) : l(Al), c(Ac), r(Ar), bl(Abl), b(Ab), br(Abr) {};
  inline Tbc_R(double *Al, double *Ac, double *Ar, double *Abl, double *Ab, double *Abr) : l(*Al), c(*Ac), r(*Ar), bl(*Abl), b(*Ab), br(*Abr) {};
  inline Tbc_R(TMaskBase *AMask)
    {double *C=AMask->C;l=*(C+OF.l);c=*C;r=*(C+OF.r);bl=*(C+OF.bl);b=*(C+OF.b);br=*(C+OF.br);};


  inline Tbc_R operator*(const Tbc_R &AOp2)
    {return Tbc_R(l*AOp2.l, c*AOp2.c, r*AOp2.r, bl*AOp2.bl, b*AOp2.b, br*AOp2.br);};
  inline Tbc_R operator/(const Tbc_R &AOp2)
    {return Tbc_R(l/AOp2.l, c/AOp2.c, r/AOp2.r, bl/AOp2.bl, b/AOp2.b, br/AOp2.br);};


};


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


class Tct_R
{
public:
  double tl, t, tr;
  double l, c, r;


  inline Tct_R(double Atl, double At, double Atr, double Al, double Ac, double Ar) : tl(Atl), t(At), tr(Atr), l(Al), c(Ac), r(Ar) {};
  inline Tct_R(double *Atl, double *At, double *Atr, double *Al, double *Ac, double *Ar) : tl(*Atl), t(*At), tr(*Atr), l(*Al), c(*Ac), r(*Ar) {};
  inline Tct_R(TMaskBase *AMask)
    {double *C=AMask->C;tl=*(C+OF.tl);t=*(C+OF.t);tr=*(C+OF.tr);l=*(C+OF.l);c=*C;r=*(C+OF.r);};

  inline Tct_R operator*(const Tct_R &AOp2)
    {return Tct_R(tl*AOp2.tl, t*AOp2.t, tr*AOp2.tr, l*AOp2.l, c*AOp2.c, r*AOp2.r);};
  inline Tct_R operator/(const Tct_R &AOp2)
    {return Tct_R(tl/AOp2.tl, t/AOp2.t, tr/AOp2.tr, l/AOp2.l, c/AOp2.c, r/AOp2.r);};

};

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


class Tlc_C
{
public:
  double tl, t;
  double l,  c;
  double bl, b;

  inline Tlc_C(double Atl, double At, double Al, double Ac, double Abl, double Ab) : tl(Atl), t(At), l(Al), c(Ac), bl(Abl), b(Ab) {};
  inline Tlc_C(double *Atl, double *At, double *Al, double *Ac, double *Abl, double *Ab) : tl(*Atl), t(*At), l(*Al), c(*Ac), bl(*Abl), b(*Ab) {};
  inline Tlc_C(TMaskBase *AMask)
    {double *C=AMask->C;tl=*(C+OF.tl);t=*(C+OF.t);l=*(C+OF.l);c=*C;bl=*(C+OF.bl);b=*(C+OF.b);};


  inline Tlc_C operator*(const Tlc_C &AOp2)
    {return Tlc_C(tl*AOp2.tl, t*AOp2.t, l*AOp2.l, c*AOp2.c, bl*AOp2.bl, b*AOp2.b);};
  inline Tlc_C operator/(const Tlc_C &AOp2)
    {return Tlc_C(tl/AOp2.tl, t/AOp2.t, l/AOp2.l, c/AOp2.c, bl/AOp2.bl, b/AOp2.b);};

};

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

// TVector struct defines a (x,y) vector. Its used in functions that accepts or returns a vector,
// like, for example, the one that returns the gradient of a magnitude in TMesh.
class TVector
{
public:

  double x, y;
  
  TVector(double AX, double AY) {x=AX; y=AY;};
  TVector(void) {};
};



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


// TMask class represents a 3x3 square window set at the position of the considered node.
// Usually, there must be defined a TMask object for each magnitude used in our computations.
// All TMask objects are associated to aa TMesh object that represents the grid. Mesh objects
// is in charge of sweeping the mask along the grid for every node.
// Conventionally, masks neme derives from array name with the addition of "w" prefix.
//
// Example:
// 
//   wU = TMask(Mesh, U);
//
// Means that mask wu is associated to array U, and Mesh is in charge of moving the mask for
// every node in the grid.
//
// As magnitudes can be defined at nodes, cellsides or cells, the same is applicable to mask.
// Comments usually refer to "cells", but it can be also "nodes" or "cellsides".
class TMask : public TMaskBase
{
protected:
  inline TIntermedMask Mul(const TMask &AOp2);
  inline TIntermedMask Div(const TMask &AOp2);
  inline TIntermedMask Add(const TMask &AOp2);
  inline TIntermedMask Sub(const TMask &AOp2);

  inline TIntermedMask Mul(const TIntermedMask &AOp2);
  inline TIntermedMask Div(const TIntermedMask &AOp2);
  inline TIntermedMask Add(const TIntermedMask &AOp2);
  inline TIntermedMask Sub(const TIntermedMask &AOp2);

  inline TIntermedMask Mul(const double AOp2);
  inline TIntermedMask Div(const double AOp2);
  inline TIntermedMask Add(const double AOp2);
  inline TIntermedMask Sub(const double AOp2);


  inline void Assign(const TMask &AOp2);
  inline void Assign(const TIntermedMask &AOp2);
  inline void Assign(const double AOp2);

public:

  TMask(TMesh &AMesh, double *AArray);
  TMask(TMesh &AMesh);
  
  inline void GoRight(void) {C++;};                   // Moves the mask 1 cell to the right.
  inline void CenterMask(int AIJ) {C=Array+AIJ;};     // Sets mask at the position of cell AIJ.
  
  // Finds what nodes of the mask are greater than AVal
  inline int NodesGreaterThan(double AVal) const    {return ((tl()>AVal)?TL_B:0)+((t()>AVal)?T_B:0)+((tr()>AVal)?TR_B:0)+((l()>AVal)?L_B:0)+((c()>AVal)?C_B:0)+((r()>AVal)?R_B:0)+((bl()>AVal)?BL_B:0)+((b()>AVal)?B_B:0)+((br()>AVal)?BR_B:0);};

  // Return the value at each cell of the mask.
  inline double tl(void) const {return *(C+OF.tl);};
  inline double t(void)  const {return *(C+OF.t);};
  inline double tr(void) const {return *(C+OF.tr);};
  inline double l(void)  const {return *(C+OF.l);};
  inline double c(void)  const {return *C;};
  inline double r(void)  const {return *(C+OF.r);};
  inline double bl(void) const {return *(C+OF.bl);};
  inline double b(void)  const {return *(C+OF.b);};
  inline double br(void) const {return *(C+OF.br);};
  
  // Sets the value at each cell of the mask.
  inline void Set_tl(double AV) const {*(C+OF.tl)=AV;};
  inline void Set_t(double AV)  const {*(C+OF.t)=AV;};
  inline void Set_tr(double AV) const {*(C+OF.tr)=AV;};
  inline void Set_l(double AV)  const {*(C+OF.l)=AV;};
  inline void Set_c(double AV)  const {*C=AV;};
  inline void Set_r(double AV)  const {*(C+OF.r)=AV;};
  inline void Set_bl(double AV) const {*(C+OF.bl)=AV;};
  inline void Set_b(double AV)  const {*(C+OF.b)=AV;};
  inline void Set_br(double AV) const {*(C+OF.br)=AV;};

  // Return pairs of values in the mask.
  inline TPair lc(void) const {return TPair(*(C+OF.l), *C);};
  inline TPair cr(void) const {return TPair(*C, *(C+OF.r));};
  inline TPair bc(void) const {return TPair(*(C+OF.b), *C);};
  inline TPair ct(void) const {return TPair(*C, *(C+OF.t));};

  inline TPair tl_t(void) const {return TPair(*(C+OF.tl), *(C+OF.t));};
  inline TPair bl_b(void) const {return TPair(*(C+OF.bl), *(C+OF.b));};
  inline TPair br_r(void) const {return TPair(*(C+OF.br), *(C+OF.r));};
  inline TPair bl_l(void) const {return TPair(*(C+OF.bl), *(C+OF.l));};
  inline TPair l_tl(void) const {return TPair(*(C+OF.l),  *(C+OF.tl));};
  inline TPair r_tr(void) const {return TPair(*(C+OF.r),  *(C+OF.tr));};
  

  // Return pairs of rows or columns in the mask.
  inline Tbc_R bc_R(void) const {return Tbc_R((TMaskBase *)this);};
  inline Tct_R ct_R(void) const {return Tct_R((TMaskBase *)this);};
  inline Tlc_C lc_C(void) const {return Tlc_C((TMaskBase *)this);};

  // Return 2x2 subsquares of the mask.
  inline Tcrt crt(void) const {return Tcrt((TMaskBase *)this);};
  

};


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

class TLMask : public TMask
{
public:
  TLMask(TMesh &AMesh, double *AArray) : TMask(AMesh, AArray) {};
  TLMask(TMesh &AMesh) : TMask(AMesh) {};

  TIntLMask operator+(const TIntLMask &AM);
  TIntLMask operator-(const TIntLMask &AM);
  TIntLMask operator*(const TIntLMask &AM);
  TIntLMask operator/(const TIntLMask &AM);
  void      operator=(const TIntLMask &AM);

  TIntLMask operator+(const TLMask &AM);
  TIntLMask operator-(const TLMask &AM);
  TIntLMask operator*(const TLMask &AM);
  TIntLMask operator/(const TLMask &AM);
  void      operator=(const TLMask &AM);

  TIntLMask operator+(const double AM);
  TIntLMask operator-(const double AM);
  TIntLMask operator*(const double AM);
  TIntLMask operator/(const double AM);
  void      operator=(const double AM);
  
  double    IC(void) {return 0.5*(c()+r());}
};

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

class TRMask : public TMask
{
public:
  TRMask(TMesh &AMesh, double *AArray) : TMask(AMesh, AArray) {};
  TRMask(TMesh &AMesh) : TMask(AMesh) {};

  TIntRMask operator+(const TIntRMask &AM);
  TIntRMask operator-(const TIntRMask &AM); 
  TIntRMask operator*(const TIntRMask &AM); 
  TIntRMask operator/(const TIntRMask &AM);
  void      operator=(const TIntRMask &AM);
                                            
  TIntRMask operator+(const TRMask &AM);
  TIntRMask operator-(const TRMask &AM); 
  TIntRMask operator*(const TRMask &AM); 
  TIntRMask operator/(const TRMask &AM);
  void      operator=(const TRMask &AM);

  TIntRMask operator+(const double AM);
  TIntRMask operator-(const double AM);
  TIntRMask operator*(const double AM);
  TIntRMask operator/(const double AM);
  void      operator=(const double AM);
  
  double    IC(void) {return 0.5*(l()+c());}
  double    ITR(void);
};

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


class TTMask : public TMask
{
public:
  TTMask(TMesh &AMesh, double *AArray) : TMask(AMesh, AArray) {};
  TTMask(TMesh &AMesh) : TMask(AMesh) {};

  TIntTMask operator+(const TIntTMask &AM);
  TIntTMask operator-(const TIntTMask &AM);
  TIntTMask operator*(const TIntTMask &AM);
  TIntTMask operator/(const TIntTMask &AM);
  void      operator=(const TIntTMask &AM);

  TIntTMask operator+(const TTMask &AM);
  TIntTMask operator-(const TTMask &AM);
  TIntTMask operator*(const TTMask &AM);
  TIntTMask operator/(const TTMask &AM);
  void      operator=(const TTMask &AM);

  TIntTMask operator+(const double AM);
  TIntTMask operator-(const double AM);
  TIntTMask operator*(const double AM);
  TIntTMask operator/(const double AM);
  void      operator=(const double AM);
  
  double    IC(void) {return 0.5*(b()+c());}
  double    ITR(void);
};

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


class TBMask : public TMask
{
public:
  TBMask(TMesh &AMesh, double *AArray) : TMask(AMesh, AArray) {};
  TBMask(TMesh &AMesh) : TMask(AMesh) {};

  TIntBMask operator+(const TIntBMask &AM);
  TIntBMask operator-(const TIntBMask &AM);
  TIntBMask operator*(const TIntBMask &AM);
  TIntBMask operator/(const TIntBMask &AM);
  void      operator=(const TIntBMask &AM);

  TIntBMask operator+(const TBMask &AM);
  TIntBMask operator-(const TBMask &AM);
  TIntBMask operator*(const TBMask &AM);
  TIntBMask operator/(const TBMask &AM);
  void      operator=(const TBMask &AM);

  TIntBMask operator+(const double AM);
  TIntBMask operator-(const double AM);
  TIntBMask operator*(const double AM);
  TIntBMask operator/(const double AM);
  void      operator=(const double AM);

  double    IC(void) {return 0.5*(c()+t());}
};

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


class TCMask : public TMask
{
public:
  TCMask(TMesh &AMesh, double *AArray) : TMask(AMesh, AArray) {};
  TCMask(TMesh &AMesh) : TMask(AMesh) {};

  TIntCMask operator+(const TIntCMask &AM);
  TIntCMask operator-(const TIntCMask &AM); 
  TIntCMask operator*(const TIntCMask &AM); 
  TIntCMask operator/(const TIntCMask &AM);
  void      operator=(const TIntCMask &AM);

  TIntCMask operator+(const TCMask &AM);
  TIntCMask operator-(const TCMask &AM);
  TIntCMask operator*(const TCMask &AM);
  TIntCMask operator/(const TCMask &AM);
  void      operator=(const TCMask &AM);

  TIntCMask operator+(const double AM);
  TIntCMask operator-(const double AM);
  TIntCMask operator*(const double AM);
  TIntCMask operator/(const double AM);
  void      operator=(const double AM);
  
  double    IR(void);
  double    IT(void);
  double    ITR(void);
  double    ITR(int AValidNodes);

  TIntLMask MIL(int AValidNodes=ALL_B);
  TIntBMask MIB(int AValidNodes=ALL_B);

  // Warning: These methods are used for deal with some tricks. In the future they should not exist.
  TIntRMask MCopyToR(void);
  TIntTMask MCopyToT(void);
};

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

class TTRMask : public TMask
{
public:
  TTRMask(TMesh &AMesh, double *AArray) : TMask(AMesh, AArray) {};
  TTRMask(TMesh &AMesh) : TMask(AMesh) {};

  TIntTRMask operator+(const TIntTRMask &AM);
  TIntTRMask operator-(const TIntTRMask &AM); 
  TIntTRMask operator*(const TIntTRMask &AM); 
  TIntTRMask operator/(const TIntTRMask &AM);
  void       operator=(const TIntTRMask &AM);

  TIntTRMask operator+(const TTRMask &AM);
  TIntTRMask operator-(const TTRMask &AM);
  TIntTRMask operator*(const TTRMask &AM);
  TIntTRMask operator/(const TTRMask &AM);
  void       operator=(const TTRMask &AM);
  
  TIntBLMask operator+(const TBLMask &AM);
  TIntBLMask operator-(const TBLMask &AM);
  TIntBLMask operator*(const TBLMask &AM);
  TIntBLMask operator/(const TBLMask &AM);

  TIntTRMask operator+(const double AM);
  TIntTRMask operator-(const double AM);
  TIntTRMask operator*(const double AM);
  TIntTRMask operator/(const double AM);
  void       operator=(const double AM);
         
  
  double    IC(void) {return 0.25*(c()+r()+b()+bl());}
};

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

class TBLMask : public TMask
{
public:
  TBLMask(TMesh &AMesh, double *AArray) : TMask(AMesh, AArray) {};
  TBLMask(TMesh &AMesh) : TMask(AMesh) {};

  TIntBLMask operator+(const TIntBLMask &AM);
  TIntBLMask operator-(const TIntBLMask &AM);
  TIntBLMask operator*(const TIntBLMask &AM);
  TIntBLMask operator/(const TIntBLMask &AM);
  void       operator=(const TIntBLMask &AM);

  TIntBLMask operator+(const TBLMask &AM);
  TIntBLMask operator-(const TBLMask &AM);
  TIntBLMask operator*(const TBLMask &AM);
  TIntBLMask operator/(const TBLMask &AM);
  void       operator=(const TBLMask &AM);
  
  TIntBLMask operator+(const TTRMask &AM);
  TIntBLMask operator-(const TTRMask &AM);
  TIntBLMask operator*(const TTRMask &AM);
  TIntBLMask operator/(const TTRMask &AM);
  
  TIntBLMask operator+(const TIntTRMask &AM);
  TIntBLMask operator-(const TIntTRMask &AM);
  TIntBLMask operator*(const TIntTRMask &AM);
  TIntBLMask operator/(const TIntTRMask &AM);

  TIntBLMask operator+(const double AM);
  TIntBLMask operator-(const double AM);
  TIntBLMask operator*(const double AM);
  TIntBLMask operator/(const double AM);
  void       operator=(const double AM);
  
  double    IC(void) {return 0.25*(c()+r()+t()+tr());}
  double    IC(int AValidNodes) {return 0.25*((VN(C_B)?c():0.0)+(VN(R_B)?r():0.0)+(VN(T_B)?t():0.0)+(VN(TR_B)?tr():0.0));}

  
};

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



class TIntermedMask
{
private:

protected:
  inline TIntermedMask Mul(const TMask &AOp2) {TIntermedMask Res(tr*AOp2.tr()); return Res;};
  inline TIntermedMask Div(const TMask &AOp2) {TIntermedMask Res(tr/AOp2.tr()); return Res;};
  inline TIntermedMask Add(const TMask &AOp2) {TIntermedMask Res(tr+AOp2.tr()); return Res;};
  inline TIntermedMask Sub(const TMask &AOp2) {TIntermedMask Res(tr-AOp2.tr()); return Res;};
         
  inline TIntermedMask Mul(const TIntermedMask &AOp2) {TIntermedMask Res(tr*AOp2.tr); return Res;};
  inline TIntermedMask Div(const TIntermedMask &AOp2) {TIntermedMask Res(tr/AOp2.tr); return Res;};
  inline TIntermedMask Add(const TIntermedMask &AOp2) {TIntermedMask Res(tr+AOp2.tr); return Res;};
  inline TIntermedMask Sub(const TIntermedMask &AOp2) {TIntermedMask Res(tr-AOp2.tr); return Res;};

  inline TIntermedMask Mul(const double AOp2) {TIntermedMask Res(tr*AOp2); return Res;};
  inline TIntermedMask Div(const double AOp2) {TIntermedMask Res(tr/AOp2); return Res;};
  inline TIntermedMask Add(const double AOp2) {TIntermedMask Res(tr+AOp2); return Res;};
  inline TIntermedMask Sub(const double AOp2) {TIntermedMask Res(tr-AOp2); return Res;};

public:
  inline TIntermedMask(double AInitialValue) {tr=AInitialValue;};
  inline TIntermedMask(void) {};


  double tr;


}; //TIntermedMask

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



class TIntLMask  : public TIntermedMask
{
friend class TLMask;
protected:
  inline TIntLMask (TIntermedMask AM) {*(TIntermedMask *)(this)=AM;};

public:
  inline TIntLMask(double AInitialValue) {tr=AInitialValue;};      



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

  inline TIntLMask operator+(const double AM)     {return Add(AM);};
  inline TIntLMask operator-(const double AM)     {return Sub(AM);};
  inline TIntLMask operator*(const double AM)     {return Mul(AM);};
  inline TIntLMask operator/(const double AM)     {return Div(AM);};
};

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

class TIntRMask  : public TIntermedMask
{ 
friend class TRMask;
protected:
  inline TIntRMask(TIntermedMask AM) {*(TIntermedMask *)(this)=AM;};
public:
  inline TIntRMask(double AInitialValue) {tr=AInitialValue;};


  inline TIntRMask operator+(const TIntRMask &AM) {return Add(AM);};
  inline TIntRMask operator-(const TIntRMask &AM) {return Sub(AM);};
  inline TIntRMask operator*(const TIntRMask &AM) {return Mul(AM);};
  inline TIntRMask operator/(const TIntRMask &AM) {return Div(AM);};
                                   
  inline TIntRMask operator+(const TRMask &AM)    {return Add(TMask(AM));};
  inline TIntRMask operator-(const TRMask &AM)    {return Sub(TMask(AM));};
  inline TIntRMask operator*(const TRMask &AM)    {return Mul(TMask(AM));};
  inline TIntRMask operator/(const TRMask &AM)    {return Div(TMask(AM));};

  inline TIntRMask operator+(const double AM)     {return Add(AM);};
  inline TIntRMask operator-(const double AM)     {return Sub(AM);};
  inline TIntRMask operator*(const double AM)     {return Mul(AM);};
  inline TIntRMask operator/(const double AM)     {return Div(AM);};
};

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

class TIntTMask  : public TIntermedMask 
{ 
friend class TTMask;
protected:
  inline TIntTMask(TIntermedMask AM) {*(TIntermedMask *)(this)=AM;};

public:
  inline TIntTMask(double AInitialValue) {tr=AInitialValue;};

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

  inline TIntTMask operator+(const double AM)     {return Add(AM);};
  inline TIntTMask operator-(const double AM)     {return Sub(AM);};
  inline TIntTMask operator*(const double AM)     {return Mul(AM);};
  inline TIntTMask operator/(const double AM)     {return Div(AM);};
};

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

class TIntBMask  : public TIntermedMask
{
friend class TBMask;
protected:
  inline TIntBMask(TIntermedMask AM) {*(TIntermedMask *)(this)=AM;};

public:
  inline TIntBMask(double AInitialValue) {tr=AInitialValue;};

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

  inline TIntBMask operator+(const double AM)     {return Add(AM);};
  inline TIntBMask operator-(const double AM)     {return Sub(AM);};
  inline TIntBMask operator*(const double AM)     {return Mul(AM);};
  inline TIntBMask operator/(const double AM)     {return Div(AM);};
};

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

class TIntCMask  : public TIntermedMask 
{ 
friend class TCMask;
protected:
  inline TIntCMask(TIntermedMask AM) {*(TIntermedMask *)(this)=AM;};

public:
  inline TIntCMask(double AInitialValue) {tr=AInitialValue;};
  
  // Creates a TIntBLMask from a TBLMask, if the cells are invalid, it set them to AVoidVal.
  inline TIntCMask(TCMask &AM, double AVoidVal, int AValidCells) {tr=(((AValidCells & TR_B)==TR_B)?AM.tr():AVoidVal);};


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

  inline TIntCMask operator+(const double AM)     {return Add(AM);};
  inline TIntCMask operator-(const double AM)     {return Sub(AM);};
  inline TIntCMask operator*(const double AM)     {return Mul(AM);};
  inline TIntCMask operator/(const double AM)     {return Div(AM);};
};


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

class TIntBLMask  : public TIntermedMask
{ 
friend class TBLMask;
protected:
  inline TIntBLMask(TIntermedMask AM) {*(TIntermedMask *)(this)=AM;};

public:
  inline TIntBLMask(double AInitialValue) {tr=AInitialValue;};
  // Creates a TIntBLMask from a TBLMask, if the cells surrounding a node are invalid, it set them to AVoidVal.
  // Note that AValidCells are cell-centered, not BL centered.
  inline TIntBLMask(TBLMask &AM, double AVoidVal, int AValidCells) {tr=((AValidCells & CTTRR_B)==CTTRR_B)?AM.tr():AVoidVal;};

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

  inline TIntBLMask operator+(const double AM)     {return Add(AM);};
  inline TIntBLMask operator-(const double AM)     {return Sub(AM);};
  inline TIntBLMask operator*(const double AM)     {return Mul(AM);};
  inline TIntBLMask operator/(const double AM)     {return Div(AM);};                                   


};

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

class TIntTRMask  : public TIntermedMask
{ 
friend class TTRMask;
protected:
  inline TIntTRMask(TIntermedMask AM) {*(TIntermedMask *)(this)=AM;};
public:
  inline TIntTRMask(double AInitialValue) {tr=AInitialValue;};
  inline TIntTRMask(TTRMask &AM, double AVoidVal, int AValidCells) {tr=((AValidCells & TR_B)==TR_B)?AM.tr():AVoidVal;};

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

  inline TIntTRMask operator+(const double AM)     {return Add(AM);};
  inline TIntTRMask operator-(const double AM)     {return Sub(AM);};
  inline TIntTRMask operator*(const double AM)     {return Mul(AM);};
  inline TIntTRMask operator/(const double AM)     {return Div(AM);};
  
};
//----------------------------------------------------------------------------------------
//----------------------------------------------------------------------------------------
//----------------------------------------------------------------------------------------

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

double R_Dx(const TRMask &AM, int AValidNodes=ALL_B);
double R_Dx(const TCMask &AM, int AValidNodes=ALL_B);
double T_Dx(const TTMask &AM, int AValidNodes=ALL_B);
double T_Dx(const TTRMask &AM, int AValidNodes=ALL_B);
double C_Dx(const TCMask &AM, int AValidNodes=ALL_B);
double C_Dx(const TCMask &AM, double AAlfa, const TVector &AVector, int AValidNodes=ALL_B);
double C_Dx(const TLMask &AM, int AValidNodes=ALL_B);


TIntCMask  CM_Dx(const TRMask &AM, int AValidNodes=ALL_B);
TIntBLMask BLM_Dx(const TTMask &AM, int AValidNodes=ALL_B);
TIntBLMask BLM_Dx(const TRMask &AM, int AValidNodes=ALL_B);
TIntLMask  LM_Dx(const TCMask &AM, int AValidNodes=ALL_B);
TIntTMask  TM_Dx(const TTRMask &AM, int AValidNodes=ALL_B);


double T_Dy(const TTMask &AM, int AValidNodes=ALL_B);
double T_Dy(const TCMask &AM, int AValidNodes=ALL_B);
double R_Dy(const TTMask &AM, int AValidNodes=ALL_B);
double R_Dy(const TTRMask &AM, int AValidNodes=ALL_B);
double C_Dy(const TCMask &AM, int AValidNodes=ALL_B);
double C_Dy(const TCMask &AM, double AAlfa, const TVector &AVector, int AValidNodes=ALL_B);
double C_Dy(const TBMask &AM, int AValidNodes=ALL_B);

TIntCMask  CM_Dy(const TTMask &AM, int AValidNodes=ALL_B);
TIntBLMask BLM_Dy(const TRMask &AM, int AValidNodes=ALL_B);
TIntBLMask BLM_Dy(const TTMask &AM, int AValidNodes=ALL_B);
TIntBMask  BM_Dy(const TCMask &AM, int AValidNodes=ALL_B);
TIntRMask  RM_Dy(const TTRMask &AM, int AValidNodes=ALL_B);
TIntCMask  CM_Dy(const TCMask &AM, double AAlfa, const TVector &AVector, int AValidNodes=ALL_B);

double C_Div(const TLMask &AMx, const TBMask &AMy, int AValidNodes=ALL_B);
double C_Div(const TCMask &AMx, const TCMask &AMy, int AValidNodes=ALL_B);


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



// TIntMask is similar to TMask but for integer arrays.
class TIntMask
{
protected:
  int *Array;



public:
  int *C;
  
  TIntMask(TMesh &AMesh, int *AArray);

  inline void GoRight(void) {C++;};        // Moves the mask 1 cell to the right.
  void CenterMask(int AIJ) {C=Array+AIJ;}; // Sets mask at the position of cell AIJ.


  // Return the value at each cell of the mask.
  inline int tl(void) const {return *(C+OF.tl);};
  inline int t(void)  const {return *(C+OF.t);};
  inline int tr(void) const {return *(C+OF.tr);};
  inline int l(void)  const {return *(C+OF.l);};
  inline int c(void)  const {return *C;};
  inline int r(void)  const {return *(C+OF.r);};
  inline int bl(void) const {return *(C+OF.bl);};
  inline int b(void)  const {return *(C+OF.b);};
  inline int br(void) const {return *(C+OF.br);};
  
  // Return an integer which bits represent what nodes are equal (or not equal) to AVal.
  inline int NodesEqualTo(int AVal) const    {return ((tl()==AVal)?TL_B:0)+((t()==AVal)?T_B:0)+((tr()==AVal)?TR_B:0)+((l()==AVal)?L_B:0)+((c()==AVal)?C_B:0)+((r()==AVal)?R_B:0)+((bl()==AVal)?BL_B:0)+((b()==AVal)?B_B:0)+((br()==AVal)?BR_B:0);};
  inline int NodesNotEqualTo(int AVal) const {return _NEq(tl(),AVal)+(_NEq(t(),AVal) << 1)+(_NEq(tr(),AVal) << 2)+(_NEq(l(),AVal) << 3)+(_NEq(c(),AVal) << 4)+(_NEq(r(),AVal) << 5)+(_NEq(bl(),AVal) << 6)+(_NEq(b(),AVal) << 7)+(_NEq(br(),AVal) << 8);};

};

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


// Similar to TMask but for arrays that only vary with X (usually X and DelX).
class TXMask
{
  double *Array;
  int Width;

public:
  double *C;


  TXMask(TMesh &AMesh, double *AArray);

  inline void GoRight(void) {C++;};  // Moves the mask 1 cell to the right.
  inline void CenterMask(int AI);    // Sets mask at the position of cell AIJ.
                                                                            

  // Return the value at each cell of the mask.
  inline double l(void) const {return *(C+OF.l);};
  inline double c(void) const {return *C;};
  inline double r(void) const {return *(C+OF.r);};

  // Return pairs of values in the mask.
  inline TPair lc(void) const {return TPair(*(C+OF.l), *C);};
  inline TPair cr(void) const {return TPair(*C, *(C+OF.r));};



};

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

// Similar to TMask but for arrays that only vary with Y (usually Y and DelY).
class TYMask
{
  double *Array;
  int Ot, Ob;
  int Height;

public:
  double *C;


  TYMask(TMesh &AMesh, double *AArray);

  inline void CenterMask(int AJ);  // Sets mask at the position of cell AIJ.

  // Return the value at each cell of the mask.
  inline double t(void) const {return *(C+Ot);};
  inline double c(void) const {return *C;};
  inline double b(void) const {return *(C+Ob);};

  // Return pairs of values in the mask.
  inline TPair bc(void) const {return TPair(*(C+Ob), *C);};
  inline TPair ct(void)const  {return TPair(*C, *(C+Ot));};



};


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


// TMesh class represents the grid where calculations are performed. This class let's do
// different things:
//
//   - It have the list of all the masks we are using, and let us move all these mask at
//     the same time along the grid. The way for doing so is shown in the following code:
//
//       Mesh.DefineArea(IMin, IMax, JMin, JMax);  // Defines the subarea we want to "sweep". If this
//                                                 // method is not called, we'll sweep the whole grid.
//       Mesh.First();          // Set all the masks at the first element.
//       do
//       {
//         // Our code.
//       } while(Mesh.Next());  // Repeat for every grid cell until all cells have been visited.
//
//   - It stores precalculated coefficients that let us calculate the magnitudes and it's spatial
//     derivatives at different points of the cells.
//
//   Some comments about notation:
//
//   In cobras model, magnitudes can be defined at different points of the cell. For example, VOF function
// is defined at the cell center, U is defined at the middle of the right side and V velocity at the middle
// of the top side.
//
// For notation purposes, we use these magnitudes F, U & V as eponym for all magnitude defined at these points.
// So all magnitudes defined at the center will be called F magnitudes, and the same for U and V.
//
// For instance, let's take a look to the method "DUDYc". This method takes an U-type magnitude and calculates
// its derivative (D) with respect to Y (DY) at the center of the cell (C).
//
// As another example, method "Ftr" returns the value of a F-type magnitude at the top right (tr) corner of the cell.
class TMesh
{
private:
  int Count;
  int CountYTCB;
  int CountXLCR;
  int CountInt;
  
  int BufferCount;
  

  TMask    *Masks[100];
  TYMask   *MasksYTCB[100];
  TXMask   *MasksXLCR[100];
  TIntMask *MasksInt[100];


  // Metric parameters.
  double YCoef2, YCoef3, YCoef4, YCoef5;
  double *FXC, *FYC, *FDelX, *FDelY;
  double *XCoefs, *YCoefs;

  bool IncludeMetricInfo;

  // Coefficients for calculating derivatives and interpolations over the grid.
  double c_dudyc_b, c_dudyc_t;
  double c_dvdxc_l, c_dvdxc_r;

  double c_dudxtr_t, c_dudxtr_c, c_dudxtr_tr, c_dudxtr_r;
  double c_dvdytr_tr, c_dvdytr_t, c_dvdytr_r, c_dvdytr_c;

  double c_dudx, c_dvdy, c_rdudx, c_tdvdy;

  double c_ftr_t, c_ftr_tr, c_ftr_c, c_ftr_r;


public:
  int I, J, IJ, Width, Height;
  int IMin, IMax, JMin, JMax;
  double XC, YC, DelX, DelY;

  TMesh(int AWidth, int AHeight, double *AXCoefs=NULL, double *AYCoefs=NULL, double *AXC=NULL, double *AYC=NULL, 
        double *ADX=NULL, double *ADY=NULL);
  ~TMesh(void);

  void AddMask(TMask *AMask)     {Masks[Count]=AMask;Count++;};
  void AddMask(TYMask *AMask)    {MasksYTCB[CountYTCB]=AMask;CountYTCB++;};
  void AddMask(TXMask *AMask)    {MasksXLCR[CountXLCR]=AMask;CountXLCR++;};
  void AddMask(TIntMask *AMask)  {MasksInt[CountInt]=AMask;CountInt++;};


  void CenterMasks(int AI, int AJ);
  void GoRight(void);

  void DefineArea(int AIMin, int AIMax, int AJMin, int AJMax) {IMin=AIMin;IMax=AIMax;JMin=AJMin;JMax=AJMax;};
  void First(void);
  bool Next(void);

  void CenterMasks2(int AI, int AJ);
  void GoRight2(void);

  void First2(void);
  bool Next2(void);
  double *NewBuffer(void);

  // Calculate derivatives of magnitudes of U, V or F types at points c, tr, etc.
  inline double DUDYc(double ATL, double AT, double AL, double AC, double ABL, double AB)
    {return (ATL+AT-AL-AC)*c_dudyc_t-(ABL+AB-AL-AC)*c_dudyc_b;};
  inline double DVDXc(double ABL, double AL, double AB, double AC, double ABR, double AR)
    {return (ABR+AR-AB-AC)*c_dvdxc_r-(ABL+AL-AB-AC)*c_dvdxc_l;};

  inline double DUDYc(double ATL, double AT, double AL, double AC, double ABL, double AB, int AValidNodes)
    {return (VN(T_B)?((ATL+AT-AL-AC)*c_dudyc_t):0.0)-(VN(B_B)?(ABL+AB-AL-AC)*c_dudyc_b:0.0);};
  inline double DVDXc(double ABL, double AL, double AB, double AC, double ABR, double AR, int AValidNodes)
    {return (VN(R_B)?((ABR+AR-AB-AC)*c_dvdxc_r):0.0)-(VN(L_B)?(ABL+AL-AB-AC)*c_dvdxc_l:0.0);};

  inline double DUDYc(const Tlc_C &Alc_C, int AValidNodes)
    {return (VN(T_B)?((Alc_C.tl+Alc_C.t-Alc_C.l-Alc_C.c)*c_dudyc_t):0.0)-(VN(B_B)?(Alc_C.bl+Alc_C.b-Alc_C.l-Alc_C.c)*c_dudyc_b:0.0);};
  inline double DVDXc(const Tbc_R &Abc_R, int AValidNodes)
    {return (VN(R_B)?((Abc_R.br+Abc_R.r-Abc_R.b-Abc_R.c)*c_dvdxc_r):0.0)-(VN(L_B)?(Abc_R.bl+Abc_R.l-Abc_R.b-Abc_R.c)*c_dvdxc_l:0.0);};

  inline double DUDYc(TMask &AMask, int AValidNodes) 
    {return (((AValidNodes & T_B)==T_B)?((AMask.tl()+AMask.t()-AMask.l()-AMask.c())*c_dudyc_t):0.0)-(((AValidNodes & B_B)==B_B)?(AMask.bl()+AMask.b()-AMask.l()-AMask.c())*c_dudyc_b:0.0);};


  inline double DUDYtr(double AC, double AT, int AValidNodes)
    {return VN(CT_B)?((AT-AC)*c_dudyc_t):0.0;};
  inline double DVDXtr(double AC, double AR, int AValidNodes)
    {return VN(CR_B)?((AR-AC)*c_dvdxc_r):0.0;};


  inline double DUDXtr(double ATL, double AL, double AT, double AC, double ATR, double AR)
    {return c_dudxtr_t*(AT-ATL)+c_dudxtr_c*(AC-AL)+c_dudxtr_tr*(ATR-AT)+c_dudxtr_r*(AR-AC);};
  inline double DVDYtr(double ATR, double AT, double AR, double AC, double ABR, double AB)
    {return c_dvdytr_t*(AT-AC) +c_dvdytr_c*(AC-AB)+c_dvdytr_tr*(ATR-AR)+c_dvdytr_r*(AR-ABR);};


  inline double DUDXtr(double ATL, double AL, double AT, double AC, double ATR, double AR, int AValidNodes)
    { return (VN(T_B)?(c_dudxtr_t*(AT-ATL)):(c_dudxtr_t*(AC-AL)))+
             c_dudxtr_c*(AC-AL)+
             (VN(TRR_B)?(c_dudxtr_tr*(ATR-AT)+c_dudxtr_r*(AR-AC)):0.0);};
  inline double DVDYtr(double ATR, double AT, double AR, double AC, double ABR, double AB, int AValidNodes)
    { return (VN(CT_B) ?(c_dvdytr_t* (AT-AC)) :0.0)+
             (VN(BC_B) ?(c_dvdytr_c* (AC-AB)) :0.0)+
             (VN(TRR_B)?(c_dvdytr_tr*(ATR-AR)):0.0)+
             (VN(BRR_B)?(c_dvdytr_r* (AR-ABR)):0.0);};


  inline double DUDXtr(const Tct_R &Act_R, int AValidNodes)
    { return (VN(T_B)?c_dudxtr_t*(Act_R.t-Act_R.tl):c_dudxtr_t*(Act_R.c-Act_R.l))+
             c_dudxtr_c*(Act_R.c-Act_R.l)+
             (VN(TRR_B)?(c_dudxtr_tr*(Act_R.tr-Act_R.t)+c_dudxtr_r*(Act_R.r-Act_R.c)):0.0);};
             
  inline double DUDXr(double AL, double AC, double AR, int AValidNodes)
    {return (VN(L_B)?(c_dudxtr_t+c_dudxtr_c)*(AC-AL):0.0)+(VN(R_B)?((c_dudxtr_tr+c_dudxtr_r)*(AR-AC)):0.0);};
    
  inline double DVDYt(double AB, double AC, double AT, int AValidNodes)
    {return (VN(B_B)?(c_dvdytr_t+c_dvdytr_c)*(AC-AB):0.0)+(VN(T_B)?((c_dvdytr_tr+c_dvdytr_r)*(AT-AC)):0.0);};



  // Derivatives of masks at points. DRx_C, means Derivative of a TRMask at the center (C) point.
  double C_Dx(const TCMask &AM, double AAlfa, const TVector &AVector, int AValidNodes);
  double Dx_C(const TRMask &AM, int AValidNodes);
  double Dx_C(const TLMask &AM, int AValidNodes);
  double Dx_T(const TTRMask &AM, int AValidNodes);
  
  TIntCMask CM_Dx(const TRMask &AM, int AValidNodes);

  double C_Dy(const TCMask &AM, double AAlfa, const TVector &AVector, int AValidNodes);
  double Dy_C(const TTMask &AM, int AValidNodes);
  double Dy_C(const TBMask &AM, int AValidNodes);
  double Dy_R(const TTRMask &AM, int AValidNodes);
  
  TIntCMask CM_Dy(const TTMask &AM, int AValidNodes);

  TIntTMask  TM_Dx(const TTRMask &AM, int AValidNodes);
  TIntRMask  RM_Dy(const TTRMask &AM, int AValidNodes);


  // This function returns the gradient of a F-kind variable at the center of the cell.
  TVector GradFc(const TMask &AMask, double AAlfa, const TVector &AV, int AValidNodes);

  // Calculate interpolation of magnitudes of U, V or F types at points c, tr, etc.
  inline double Ftr(double AT, double ATR, double AC, double AR) {return c_ftr_t*AT+c_ftr_tr*ATR+c_ftr_c*AC+c_ftr_r*AR;};
  inline double Ftr(const Tcrt &Acrt) {return c_ftr_t*Acrt.t+c_ftr_tr*Acrt.tr+c_ftr_c*Acrt.c+c_ftr_r*Acrt.r;};
  
  

  
  // Calculate interpolation of masks at points where they are not defined.
  inline double IR (const TCMask *AM) {return (c_ftr_t+c_ftr_c)*AM->c()+(c_ftr_tr+c_ftr_r)*AM->r();};
  inline double IT (const TCMask *AM) {return (c_ftr_r+c_ftr_c)*AM->c()+(c_ftr_tr+c_ftr_t)*AM->t();};
  inline double ITR(const TCMask *AM) {return c_ftr_t*AM->t()+c_ftr_c*AM->c()+c_ftr_tr*AM->tr()+c_ftr_r*AM->r();};
  inline double ITR(const TTMask *AM) {return (c_ftr_t+c_ftr_c)*AM->c()+(c_ftr_tr+c_ftr_r)*AM->r();};
  inline double ITR(const TRMask *AM) {return (c_ftr_r+c_ftr_c)*AM->c()+(c_ftr_tr+c_ftr_t)*AM->t();};

  double ITR(const TCMask *AM, int AValidNodes);

  inline TIntLMask MIL(const TCMask *AM, int AValidNodes=ALL_B)
    {return TIntLMask(VN(T_B|TR_B)?(c_ftr_t+c_ftr_c)*AM->t()+(c_ftr_tr+c_ftr_r)*AM->tr(): (VN(T_B)?AM->t():AM->tr()));};
  inline TIntBMask MIB(const TCMask *AM, int AValidNodes=ALL_B)
    {return TIntBMask(VN(R_B|TR_B)?(c_ftr_r+c_ftr_c)*AM->r()+(c_ftr_tr+c_ftr_t)*AM->tr(): (VN(R_B)?AM->r():AM->tr()));};



  // Indexes to the corresponding elements.
  inline int TLij(void) {return IJ+OF.tl;};
  inline int Tij(void)  {return IJ+OF.t;};
  inline int TRij(void) {return IJ+OF.tr;};
  inline int Lij(void)  {return IJ+OF.l;};
  inline int Cij(void)  {return IJ;};
  inline int Rij(void)  {return IJ+OF.r;};
  inline int BLij(void) {return IJ+OF.bl;};
  inline int Bij(void)  {return IJ+OF.b;};
  inline int BRij(void) {return IJ+OF.br;};
  
  // Returns true if the cell is inside the defined area.
  inline bool IsInDefinedArea(void) {return (I>=IMin && I<=IMax && J>=JMin && J<=JMax);};


};

//----------------------------------------------------------------------------------------
#endif