/*====================================================================================
! Purpose(s): 
!
! Contain(s): cstress
!
! Author(s):  Javier Lopez Lara (jav.lopez@unican.es)
!             Gabriel Barajas   (barajasg@unican.es)
!====================================================================================*/
#include "CArrayMasks.h"
#include "CVars.h"
#include "CMisc.h"
#include "CErrors.h"
#include "math.h"
#include <stdlib.h>
#include <stdio.h>
#include "CMathMisc.h"


#include "CCalcFricVel.h"
#include "CSetFricVel.h"

#ifdef __INTEL_COMPILER
#define PYM cstress
#else 
#define PYM cstress_
#endif

extern "C"  void _FORTRAN PYM ( _DOUBLE_ARRAY AAC, _DOUBLE_ARRAY AAR, _DOUBLE_ARRAY AAT,  
	                       	     _DOUBLE_ARRAY AU, _DOUBLE_ARRAY AV, _DOUBLE_ARRAY AXNut, 
		                     _DOUBLE_ARRAY ATauxx, _DOUBLE_ARRAY ATauyy, _DOUBLE_ARRAY ATauxy,  
                        	     _DOUBLE_ARRAY ATauxxTurb, _DOUBLE_ARRAY ATauyyTurb, _DOUBLE_ARRAY ATauxyTurb, 
                        	     _DOUBLE_ARRAY APorosity, _DOUBLE_ARRAY AK, _DOUBLE_ARRAY AEp, _DOUBLE_ARRAY AF, 
                        	     _DOUBLE_ARRAY AUtmp, _DOUBLE_ARRAY AVtmp,                          
                        	     _DOUBLE_ARRAY AXCoefs, _DOUBLE_ARRAY AYCoefs, _DOUBLE_ARRAY AXi, _DOUBLE_ARRAY AYj, _DOUBLE_ARRAY ADelX, _DOUBLE_ARRAY ADelY,	 
	                             _DOUBLE AMu, _DOUBLE ANu, _DOUBLE Asd, _INT Anks, 
                        	     _DOUBLE ARhoF, _INT AKb, _INT AKt, _INT AKl, _INT AKr, _INT AIMax, _INT AJMax, 
		                     _INT_ARRAY ANf, _INT Aim1, _INT Ajm1, _INT AKEpModel, _INT Anrough) 
{
/*!=======================================================================
! Purpose(s): 
!
! Author(s):  Javier Lopez Lara (jav.lopez@unican.es)
!             Gabriel Barajas   (barajasg@unican.es)
!=======================================================================*/

  const double MIN_NUT = 1e-20;
  const double MIN_EP  = 1.0e-32;

  double C1, C2, C3, Cd;
  double UxC,  UyC,  VxC,  VyC;
  double UxTR, UyTR, VxTR, VyTR;
  double Smax, DmaxC, DmaxTR;
  double xk3xep2c, xnutc, xnutr, xnutru, xxl1;
  double xxl1turb, xxnl1, xxnl2, xxnl3, xyl1, xyl1turb, xynl2, xynl3;
  double yyl1, yyl1turb, yynl1, yynl2, yynl3;

  double provK0, provK1, provK2, provK3, provK4, provK5, provK6, provK7, provK8, provK9;


  bool KbRigid = (AKb==bcRigidFreeSlip || AKb==bcRigidNoSlip);
  bool KtRigid = (AKt==bcRigidFreeSlip || AKt==bcRigidNoSlip);
  bool KlRigid = (AKl==bcRigidFreeSlip || AKl==bcRigidNoSlip);
  bool KrRigid = (AKr==bcRigidFreeSlip || AKr==bcRigidNoSlip);

  // Is body surface rough?
  bool IsRoughSurface=(Anrough==1);

  // Roughness.
  double nks=Anks*1.0;
  double Ks=nks*Asd;

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

  // Creates the masks.
  TCMask      wAC        (Mesh, AAC);
  TRMask      wAR        (Mesh, AAR);
  TTMask      wAT        (Mesh, AAT);
  TRMask      wU         (Mesh, AU);
  TTMask      wV         (Mesh, AV);
  TCMask      wK         (Mesh, AK);
  TCMask      wEp        (Mesh, AEp);

  TCMask      wF         (Mesh, AF);
  TCMask      wPorosity  (Mesh, APorosity);
  TCMask      wNuT       (Mesh, AXNut);
  TIntMask    wNf        (Mesh, ANf);
  TCMask      wTauxx     (Mesh, ATauxx);
  TTRMask     wTauxy     (Mesh, ATauxy);
  TCMask      wTauyy     (Mesh, ATauyy);
  TCMask      wTauxxTurb (Mesh, ATauxxTurb);
  TTRMask     wTauxyTurb (Mesh, ATauxyTurb);
  TCMask      wTauyyTurb (Mesh, ATauyyTurb);

  TRMask      wrPorosity (Mesh);
  TTMask      wtPorosity (Mesh);

  TRMask      wUPor      (Mesh);
  TTMask      wVPor      (Mesh);

  TCMask      wUPorx     (Mesh);
  TBLMask     wUPory     (Mesh);
  TBLMask     wVPorx     (Mesh);
  TCMask      wVPory     (Mesh);

  TCMask      wK3_Ep2Por2(Mesh);

  Mesh.DefineArea(1, AIMax, 1, AJMax);

  Mesh.First2();
  do
  {

    // Finds the nodes that will be considered in various derivatives.
    int OpenToFlowCells=wAC.NodesGreaterThan(1e-6);
    int ValidNodes=wNf.NodesNotEqualTo(ctEmptyCell);

    ValidNodes=ValidNodes & (wF.NodesGreaterThan(1e-6));

    wrPorosity = wPorosity.MCopyToR();
    wtPorosity = wPorosity.MCopyToT();

    wUPor = wU/wrPorosity;
    wVPor = wV/wtPorosity;

    wK3_Ep2Por2 = Cube(wK)/(Sqr(wEp*wPorosity));

    // Evaluate du/dx, du/dy, dv/dx and dv/dy at its natural points.
    wUPorx = CM_Dx (wUPor, ValidNodes);
    wUPory = BLM_Dy(wUPor, ValidNodes);
    wVPorx = BLM_Dx(wVPor, ValidNodes);
    wVPory = CM_Dy (wVPor, ValidNodes);

    if (Mesh.IsInDefinedArea())
    // If the cell is inside the area of calculus...
    {
      UxC = wUPorx.c();
      UyC = wUPory.IC();
      VxC = wVPorx.IC();
      VyC = wVPory.c();

      UxTR = wUPorx.ITR();
      UyTR = wUPory.tr();   // <- tr() in a BL mask is equivalent to c() in a TR mask.
      VxTR = wVPorx.tr();
      VyTR = wVPory.ITR();

      if (AKEpModel==keNoTurbulence)
      // If turbulence is not coinsidered...
      {
        if ((OpenToFlowCells & C_B)==C_B)
        // If cell is open to flow...
        {
          *wTauxx.C = 2.0*AMu*UxC;
          *wTauyy.C = 2.0*AMu*VyC;
          *wTauxy.C = AMu*(UyTR+VxTR);
        }
        else
        {
          *wTauxx.C = 0.0;
          *wTauyy.C = 0.0; 
          *wTauxy.C = 0.0;
        }

      }
      else
      // If turbulence is coinsidered...
      {
        // when kemodel=1 & 4, turbulence production is calculated by the
        // product of stress and strain. It is noted the stress here
        // contains Reynolds stress only based on the derivation. The
        // inclusion of the molecular stress may contaminate results when
        // turbulence is low, though it is insignificant for large
        // turbulence flow. Thus, it is more accurate to introduce Reynolds
        // stress (turbulence) and total stress (mean flow) separatedly.
        // in certain case, very small epsilon (e.g., 1.0e-160) may occur during
        // computation; the third power will make it become zero.
        if (fabs(wEp.c())<MIN_EP)
        // if Ep is really small...
        {
          *wTauxxTurb.C=0.0;
          *wTauxyTurb.C=0.0;
          *wTauyyTurb.C=0.0;

          *wTauxx.C=0.0;
          *wTauxy.C=0.0;
          *wTauyy.C=0.0;
        } 
        else
        // if Ep is big enough...
        {

          switch (AKEpModel)
          {
      	    case keLinearEddyVisc: //standard k-e model

                printf("Error, Linear Eddy Visc turbulence model is not yet implemented", __FILE__, __LINE__);
              break;


      	    case keNonLinearEddyVisc:

    	
              if ((OpenToFlowCells & C_B)==C_B)
              // If cell is open to flow...
              {
                // Calculate the type of obstacle at the cell.
                TVector ObsSrfce;
                double CDist, TRDist;
                bool IsObsIntrfce=IsObstacleInterface(Mesh, wAC, wAR, wAT, 
                                                      KlRigid, KrRigid, KtRigid, KbRigid, ObsSrfce, CDist, TRDist);

                double K_Ep=wK.c()/wEp.c();

                // Accounts for cells near the obstacles.
                if (IsObsIntrfce)
                // If cell is near an obstacle...
                {
          	      // The derivative of the component of the velocity parallel to the surface in the
          	      // direction perpendicular to it is given by the logaritmic law.
                	
          	      // If we call n the unitary outward normal vector to the surface ("ObsSrfce" in the code) and 
          	      // t a unitary tangential one, then, the derivative of the parallel component of vector V along the
          	      // direction perpendicular to the surface can be written as:
          	      //
          	      //  n\B7grad(t\B7V)
                	
                  // Calculate a unitary vector tangent to the outward normal from the solid surface.
                  // We choose the one that is 90\BA clockwise from the normal one. The direction of this
                  // vector is important because we store the frictional velocity with its sign.
                  TVector ObsSrfceTan=NormalVector(ObsSrfce);                  
                  TVector VelC(wUPor.IC(), wVPor.IC());
                  
                  // Find the projection of V at the center of the cell along the direction parallel to the surface.
                  double VC=DotProduct(ObsSrfceTan, VelC);
                  
                  // Computes frictional velocity at the center of the cell.
                  double FricVC=0.0;
                  if (wAC.c()>OBSTACLE) FricVC=CalcFrictionVelocity(E, VC, CDist, ANu, Kappa, IsRoughSurface, Ks);
                  
                  // Modifies the gradient tensor in order to include the logarithmic law profile.
                  // This routine is completely general and works for any angle of the surface.
                  CSetFricVel(UxC, UyC, VxC, VyC, FricVC, CDist, ObsSrfce);

                  // Does exactly the same thing for the TR point.
                  TVector VelTR(wUPor.ITR(), wVPor.ITR());
                  double VTR=DotProduct(ObsSrfceTan, VelTR);

                  double FricVTR=0.0;
                  if (wAC.c()>OBSTACLE) FricVTR=CalcFrictionVelocity(E, VTR, TRDist, ANu, Kappa, IsRoughSurface, Ks);

                  CSetFricVel(UxTR, UyTR, VxTR, VyTR, FricVTR, TRDist, ObsSrfce);
                }

                // Compute some coefficients.
                DmaxC  = K_Ep*max2(fabs(UxC+UyC),   fabs(VxC+VyC));

                C1 =  1.0/(185.2+3.0*DmaxC*DmaxC);
                C2 = -1.0/( 58.5+2.0*DmaxC*DmaxC);
                C3 =  1.0/(370.4+3.0*DmaxC*DmaxC);

                // Computes some previous values.

                double Ux2C=Sqr(UxC);
                double Uy2C=Sqr(UyC);
                double Vx2C=Sqr(VxC);
                double Vy2C=Sqr(VyC);
                double UyVxC=UyC*VxC;
                double Ux2Uy2Vx2Vy2C=Ux2C+Uy2C+Vx2C+Vy2C;

                // Calculate taus.
                *wTauxxTurb.C = ARhoF*(2.0*wNuT.c()*UxC - (2.0/3.0)*wK.c()/wPorosity.c() +
                                      wK3_Ep2Por2.c()*(C1*(2.0*(Ux2C+UyVxC) - (2.0/3.0)*(Ux2C+2.0*UyVxC+Vy2C)) +
                                                       C2*(Ux2C+Uy2C-(1.0/3.0)*Ux2Uy2Vx2Vy2C) +
                                                       C3*(Ux2C+Vx2C-(1.0/3.0)*Ux2Uy2Vx2Vy2C)));

                *wTauyyTurb.C = ARhoF*(2.0*wNuT.c()*VyC - (2.0/3.0)*wK.c()/wPorosity.c() +
                                      wK3_Ep2Por2.c()*(C1*(2.0*(Vy2C+UyVxC) - (2.0/3.0)*(Ux2C+2.0*UyVxC+Vy2C)) +
                                                       C2*(Vx2C+Vy2C-(1.0/3.0)*Ux2Uy2Vx2Vy2C) +
                                                       C3*(Uy2C+Vy2C-(1.0/3.0)*Ux2Uy2Vx2Vy2C)));

                *wTauxx.C = 2.0*AMu*UxC + wTauxxTurb.c();
                *wTauyy.C = 2.0*AMu*VyC + wTauyyTurb.c();


              
                // find some magnitudes at the TR corner.
                double NuTtr, K3_Ep2Por2tr;
                int ValidEpNodes=wEp.NodesGreaterThan(MIN_EP);
                NuTtr=wNuT.ITR(ValidEpNodes);
                K3_Ep2Por2tr=wK3_Ep2Por2.ITR(ValidEpNodes);
                
                DmaxTR = K_Ep*max2(fabs(UxTR+UyTR), fabs(VxTR+VyTR));
                C2 = -1.0/( 58.5+2.0*DmaxTR*DmaxTR);
                C3 =  1.0/(370.4+3.0*DmaxTR*DmaxTR);

                *wTauxyTurb.C = ARhoF*(NuTtr*(UyTR+VxTR) +
                                      K3_Ep2Por2tr*(C2*(UxTR*VxTR+UyTR*VyTR) +
                                                    C3*(UxTR*UyTR+VxTR*VyTR)));

                *wTauxy.C = AMu*(UyTR+VxTR) + wTauxyTurb.c();
                
              }
              else
              // If cell is not open to flow...
              {
                // Values of stresses on closed cells are assumed to be zero.
                *wTauxx.C = 0.0;
                *wTauyy.C = 0.0; 
                *wTauxy.C = 0.0;
                *wTauxxTurb.C = 0.0;
                *wTauyyTurb.C = 0.0; 
                *wTauxyTurb.C = 0.0; 
              }

              break; //keNonLinearEddyVisc

          } //switch

        } //if

      } //if

    } //if

  } while(Mesh.Next2());

}