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

#include "CArrayMasks.h"
#include <math.h>
#include "CMathMisc.h"
#include "CMisc.h"
#include <assert.h>
#include "CVars.h"
#include <stdlib.h>
#include <stdio.h>

#ifdef __INTEL_COMPILER
#define PYM cbuildppematrixsoil
#else 
#define PYM cbuildppematrixsoil_
#endif

extern "C"  void _FORTRAN PYM (_INT AWaveMakerType,
                                 	   _DOUBLE AT, 
					   _DOUBLE ADelT, 
					   _DOUBLE ARhoF, 
					   _DOUBLE APSat, 
                                 	   _DOUBLE AFreeSurf, 
					   _DOUBLE ASource,
                                 	   _INT AKL, 
					   _INT AKT, 
                                 	   _INT AKR, 
					   _INT AKB,
                                 	   _INT AISourceS, 
					   _INT AISourceE, 
					   _INT AJSourceS, 
					   _INT AJSourceE,
                                 	   _DOUBLE_ARRAY AXCoefs, 
					   _DOUBLE_ARRAY AYCoefs, 
                                 	   _DOUBLE_ARRAY AXi,     
					   _DOUBLE_ARRAY AYj,
                                 	   _DOUBLE_ARRAY ADelX,   
					   _DOUBLE_ARRAY ADelY,
                                 	   _DOUBLE_ARRAY AAC,     
					   _DOUBLE_ARRAY AAR,     
					   _DOUBLE_ARRAY AAT, 
                                 	   _DOUBLE_ARRAY ADX,     
					   _DOUBLE_ARRAY ADY,
        	                           _DOUBLE_ARRAY AU,      
					   _DOUBLE_ARRAY AV, 
                                 	   _INT Aim1,             
					   _INT Ajm1,
        		                   _DOUBLE_ARRAY APorosity, 
                                 	   _DOUBLE_ARRAY AF, 
                                 	   _DOUBLE_ARRAY AP,      
					   _DOUBLE_ARRAY APn, 
                                 	   _DOUBLE_ARRAY ACVol,   
					   _DOUBLE_ARRAY AGc, 
                                 	   _INT_ARRAY Axnpc,
                                 	   _DOUBLE_ARRAY AA, 
					   _DOUBLE_ARRAY AB, 
					   _DOUBLE_ARRAY ABF, 
					   _DOUBLE_ARRAY AD, 
					   _DOUBLE_ARRAY ADF, 
                                 	   _DOUBLE_ARRAY AForcingTerm,
					   _DOUBLE AKw, _DOUBLE Asaturation, _DOUBLE Apermeability, _DOUBLE Aporo, _DOUBLE Az_soil)
{
/*!=======================================================================
! Purpose(s): 
!
! Author(s):  Javier Lopez Lara (jav.lopez@unican.es)
!             Gabriel Barajas   (barajasg@unican.es)
!=======================================================================*/

  int IMax=Aim1+1;
  int JMax=Ajm1+1;
  int MaxSize=(Aim1-1)*(Ajm1-1);
  int EqCount=0;
  double a, b, bf, d, df, ForcingTerm;


  // Creates the mesh object.
  TMesh Mesh(IMax, JMax, AXCoefs, AYCoefs, AXi, AYj, ADelX, ADelY);
  

  // Creates the masks.
  TMask     wAC(Mesh, AAC);
  TMask     wAR(Mesh, AAR);
  TMask     wAT(Mesh, AAT);
  TMask     wU(Mesh, AU);
  TMask     wV(Mesh, AV);
  TMask     wF(Mesh, AF);
  TMask     wP(Mesh, AP);
  TMask     wPn(Mesh, APn);
  TMask     wCellVol(Mesh, ACVol);
  TMask     wPor(Mesh, APorosity);
    
  TIntMask  wXPorIdx(Mesh, Axnpc);
  // javi
  TIntMask  wPorIdx(Mesh, Axnpc);

  TXMask  wDX(Mesh, ADelX);
  TYMask  wDY(Mesh, ADelY);
   
  
  
  Mesh.DefineArea(2, Aim1, 2, Ajm1);


  Mesh.First();

  do
  {
  	int i=Mesh.I;
  	int j=Mesh.J;
  	int ij=Mesh.IJ;
  	


    // PREVIOUS CALCULATIONS.
    // ----------------------
    // These magnitudes will be used later.
    // Compute fluid density at the centers of the sides of the cell.
  	double rhoL=ARhoF*(wDX.l()*wF.c()+wDX.c()*wF.l())/wDX.lc().Sum();
  	double rhoR=ARhoF*(wDX.r()*wF.c()+wDX.c()*wF.r())/wDX.cr().Sum();
  	double rhoT=ARhoF*(wDY.t()*wF.c()+wDY.c()*wF.t())/wDY.ct().Sum();
  	double rhoB=ARhoF*(wDY.b()*wF.c()+wDY.c()*wF.b())/wDY.bc().Sum();

  	// Compute fluid density at the center of the cell.
  	double Rho=wF.c()*ARhoF;


    // PPE Matrix can be diagonallized if pressure is zero in void cells. For doing this, we first
    // substract PSat from it.
    *wP.C=wPn.c()-APSat;


    // Computes velocity divergence.
    double DivV;
    if (wAC.c()>=1.0e-06)
    // If the cell is open to flow...
    {
      double dutdxc = (wAR.lc()*wU.lc()).Delta()/wDX.c();
      double dvtdyc = (wAT.bc()*wV.bc()).Delta()/wDY.c();
    	DivV=(dutdxc+dvtdyc);
    }
    else
    {
      DivV=0.0;	
    }

	//javi put this
  	//if (wPor.c()<1.0) DivV=0.0;	
	//javi put this
  	
  	// LAPLACIAN METRICS.
  	// ----------------------------
  	// Compute the coefficients of PPE matrix. These coefficients are multiplied by cell area.
    double Alfl=wAR.l()*2.0*wDY.c()/(wDX.lc().Sum()*rhoL);
  	double Alfr=wAR.c()*2.0*wDY.c()/(wDX.cr().Sum()*rhoR);
    double Gamt=wAT.c()*2.0*wDX.c()/(wDY.ct().Sum()*rhoT);
    double Gamb=wAT.b()*2.0*wDX.c()/(wDY.bc().Sum()*rhoB);

    //if (j==2 || i==2 || j==JMax-1 || i==IMax-1) DivV=0;

  	// INITIAL FORCING TERM.
  	// ----------------------------
    // Computes forcing term.
    ForcingTerm=-wCellVol.c()*DivV/ADelT;
    

    // POROUS MEDIA.
    // -------------
    // Modifies matrix coefficients to bring porous media into account.
    // We always assume non-lineal friction.
    

     int PorousCells=wXPorIdx.c();
     if (PorousCells>0)
    // If the cell has porous material at any of the sides...
     {
      // If a cell side is open to flux and there is porosity associated, modifies
      // matrix coefficients to take account of it.

      // Gets the porosity index at every side of the cell.
      int PIdxL, PIdxR, PIdxT, PIdxB;
      IdxUnpack(PorousCells, PIdxL, PIdxR, PIdxT, PIdxB);

      if (wAR.l()>OBSTACLE && PIdxL>0) Alfl=2.0*wDY.c()/(wDX.lc().Sum()*rhoL*AGc[PIdxL]);
      if (wAR.c()>OBSTACLE && PIdxR>0) Alfr=2.0*wDY.c()/(wDX.cr().Sum()*rhoR*AGc[PIdxR]);
      if (wAT.c()>OBSTACLE && PIdxT>0) Gamt=2.0*wDX.c()/(wDY.ct().Sum()*rhoT*AGc[PIdxT]);
      if (wAT.b()>OBSTACLE && PIdxB>0) Gamb=2.0*wDX.c()/(wDY.bc().Sum()*rhoB*AGc[PIdxB]);
 
     } //if
     
    a = Alfr+Alfl+Gamb+Gamt;
    b =-Alfr;
    bf=-Alfl;
    d =-Gamt;
    df=-Gamb;

	

    if (ARhoF*wF.l()<AFreeSurf) bf=0.0;
    if (ARhoF*wF.r()<AFreeSurf) b=0.0;
    if (ARhoF*wF.b()<AFreeSurf) df=0.0;
    if (ARhoF*wF.t()<AFreeSurf) d=0.0;  


/*if (rhoL<AFreeSurf) bf=0.0;
if (rhoR<AFreeSurf) b=0.0;
if (rhoB<AFreeSurf) df=0.0;
if (rhoT<AFreeSurf) d=0.0;*/

// JAVI PUT THIS
  	//if (wPor1.c()<1.0 && (ADelT!=AT))
	if (j<=Az_soil && (ADelT!=AT) && (wAC.c()>=OBSTACLE))
    {	
	 // ADDTIONAL TERMS FOR NON-SATURATED MECHANICS.
  	 // ----------------------------		
	 // Computes forcing term.	   
	 double beta = 1/AKw+(1-Asaturation)/wPn.c(); 	  	      
	 ForcingTerm = 9.8*Aporo*wPn.c()*beta*wCellVol.c()/(ADelT*Apermeability);
	 // Computes additional term for the diagonal matrix.	   
       a = a + 9.8*Aporo*beta*wCellVol.c()/(ADelT*Apermeability);
	  } 
// JAVI PUT THIS 


    

  	// Void cells are those empty of fluid.
    bool VoidCell=(ARhoF*wF.c()<AFreeSurf) || (wAC.c()<OBSTACLE);


    if (VoidCell)
    // If it is a void cell...
    {
      // For void cells, pressure is always Psat. We are solving for P'=P-Psat everywhere, so we can set P'
      // to zero in void cells.
      a=1.0;
      b=0.0;
      bf=0.0;
      d=0.0;
      df=0.0;
      ForcingTerm=0.0;
    }


    // WAVEMAKER-
    // ----------
    // Adds wave maker to the forcing term.
    if (AWaveMakerType==wtSourceRegion && i>=AISourceS && i<=AISourceE && j>=AJSourceS && j<=AJSourceE)
    // If we have a source function and cell is inside source rectangle...
    {
    	ForcingTerm=ForcingTerm+wAC.c()*wCellVol.c()*ASource;
    } //if


    
    
    // BOUNDARY CONDITIONS.
    // --------------------
    // Internal obstacle boundaries (enforce Neumann conditions).
    if (i>2 && i<IMax-1 && j>2 && j<JMax-1)
    // If cell is in the interior...
    {       if (wAC.c()>=OBSTACLE)
       // If cell is open to flow...
       {

       	  if (wAC.t()<OBSTACLE) 
		  {
			  a=a+d;  d=0.0;
			  if (j<=Az_soil && ADelT!=AT) ForcingTerm = ForcingTerm - 9.8*wCellVol.c()/wDY.c();
		  }
       	  if (wAC.b()<OBSTACLE) 
		  {
			  a=a+df; df=0.0;
			  if (j<=Az_soil && ADelT!=AT) ForcingTerm = ForcingTerm + 9.8*wCellVol.c()/wDY.c();  
		  }
       	  if (wAC.l()<OBSTACLE) { a=a+bf; bf=0.0;		  }
       	  if (wAC.r()<OBSTACLE) { a=a+b;  b=0.0;		  }

       } //if
       
    } 
    else
    {

      if (i==2 && AKL!=bcPeriodic)
      // Left boundary...
      {
    	  if (AKL==bcAppliedPressure) ForcingTerm=ForcingTerm-bf*wP.l();
                               else a=a+bf;
        bf=0;
		  

      }

      if (j==JMax-1)
      // Top boundary...
      {
    	  if (AKT==bcAppliedPressure) ForcingTerm=ForcingTerm-d*wP.c();
                               else a=a+d;
        d=0;		  

      }

      if (i==IMax-1 && AKR!=bcPeriodic)
      // Right boundary...
      {
    	  if (AKR==bcAppliedPressure) ForcingTerm=ForcingTerm-b*wP.c();
                               else a=a+b;
        b=0;
	  

      } 

      if (j==2)
      // Bottom boundary...
      {
    	  if (AKB==bcAppliedPressure) ForcingTerm=ForcingTerm-df*wP.b();
		  else {a=a+df;
		  df=0;
		  if (j<=Az_soil && ADelT!=AT) ForcingTerm = ForcingTerm + 9.8*wCellVol.c()/wDY.c();
      }
      

      } //if

    } //if
    


    // MATRIX ASSEMBLING.
    // ------------------

    *AA++=a;
    *AB++=b;
    *AD++=d;
    *ABF++=bf;
    *ADF++=df;
    *AForcingTerm++=ForcingTerm;


          





  } while(Mesh.Next());

}