//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
//  Source File: QRelations
//  Author:      jcm
//  Date:        Sat - Nov 18, 2000
//  Namespace:   conejo
//  Framework:   Tecolote
//  Copyright:   Los Alamos National Laboratory
//               Full Copyright=$(TECOLOTE_ROOT)/Doc/Copyright
//  RCS_VERSION_ID: $Id:
//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::


#ifndef  __conejo_QRelations_t_hh
#define  __conejo_QRelations_t_hh

#include "PhysicsBaseClasses_src/HelperClasses/MMField.hh"
#include "Demo_src/Model/CompatibleHydro/QRelations.hh"
#include "Evaluator/ScalarCode.h"

//#define ENTER(a)
#define ENTER(a) tecout << "Entering " << a << endl;

namespace conejo
{
   using namespace TecFramework;
   using namespace poomalote;
   using namespace PhysicsBaseClasses;
   using namespace Hydrodynamics;
   using namespace std;

   const Real spokeCutoff = 1.0e-12; // We get this value from Ed's code

   // $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
   template <class Traits>
   BEGIN_PERSISTENT(QRelations<Traits>)
      PERSISTENT( linearQ,  "LinearQ" )
      PERSISTENT( quadQ,    "QuadQ"   )
   END_PERSISTENT

   //$ linearQ : Real - linear coefficient for Q
   //
   //$ quadQ   : Real - quadratic coefficient for Q
   //
   // $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

    //-----------------------------------------------------------------------------
    // Test function
    //-----------------------------------------------------------------------------
    template<class Traits>
    struct EdgeQPressureInfo
    {
      void scalarCodeInfo(ScalarCodeInfo &info) const
      {
        info.arguments(5);

        info.write(0, true );
        info.write(1, false);
        info.write(2, false);
        info.write(3, false);
        info.write(4, false);
//        info.write(5, false);
//        info.write(6, false);
 //       info.write(7, false);
        info.useGuards(0, false);
        info.useGuards(1, true);
        info.useGuards(2, true);
        info.useGuards(3, true);
        info.useGuards(4, true);
//        info.useGuards(5, true);
//        info.useGuards(6, true);
//        info.useGuards(7, true);

        info.dimensions(Traits::Dim);

        for (int i = 0; i < Traits::Dim; ++i)
        {
          info.lowerExtent(i) = 2;
          info.upperExtent(i) = 2;
        }
      }
    };

    template<class Traits>
    struct ScalarEdgeQPressure
      : public EdgeQPressureInfo<Traits>
    {
      // Typedefs
      FIELD_TYPEDEFS(Traits)

      ScalarEdgeQPressure(const Real& inLinearQ )
        : EdgeQPressureInfo<Traits>(),
          linearQ(inLinearQ)
      {
      }

      template<class F1, class F2, class F3, class F4, class F5>
      void operator()(const F1& EdgeQPressure,     const F2& EdgeGammaConstant, const F3& EdgeSoundSpeed,
                      const F4& EdgeVelocity,      const F5& EdgePsiLimiter,    const Loc<Traits::Dim> &loc)
      {
		if( EdgePsiLimiter(loc) < eps ) {
		    EdgeQPressure(loc) = 0.0;
		    return;
		}

		Real edgeVelocityMagnitude = sqrt(dot(EdgeVelocity(loc),EdgeVelocity(loc)));
		
		EdgeQPressure(loc) = edgeVelocityMagnitude * EdgePsiLimiter(loc) *
		                     (EdgeGammaConstant(loc) * edgeVelocityMagnitude +
		                      sqrt( linearQ * linearQ * EdgeSoundSpeed(loc) * EdgeSoundSpeed(loc) +
		                      EdgeGammaConstant(loc) * EdgeGammaConstant(loc) * edgeVelocityMagnitude * edgeVelocityMagnitude));
      }

    private:
      Real linearQ;
    };

	//======================================================================
	// Constructor -- QRelations<Traits>::QRelations
	//======================================================================

	template <class Traits>
	QRelations<Traits>::QRelations( DataDirectory* pDataDir, const string& inName )
	 :  CompatibleRelations(pDataDir,inName),
	    Old(pDataDir->strictGet<DataDirectory>("CompatibleHydroOld")),
	    linearQ(1.0),
	    quadQ(1.0)
	{
        VectorField& EdgeLength   = DataDir.get<VectorField>( "EdgeLength",   Mesh.getField<VectorField>    ( AllEdge() ) );
	    for(int d=0;d<Dim;++d) {
	        Interval<1> CVert(CoarseVert[d]);
	        RDomain.push_back(CoarseVert);
	        RDomain[d][d] = Interval<1>(CVert.first(),    CVert.last() - 2);
	        LDomain.push_back(CoarseVert);
	        LDomain[d][d] = Interval<1>(CVert.first() + 1,CVert.last() - 1);
	        Loc<Dim> offset(0); offset[d] = 1;
	        RightEdgeNgbr.push_back(RDomain[d]);
	        RightEdgeNgbr[d] += offset;
	        LeftEdgeNgbr.push_back(LDomain[d]);
	        LeftEdgeNgbr[d] -= offset;
	        CoarseEdges.push_back(Range<Dim>(EdgeLength[d].domain()));
	        LowerVert.push_back(CoarseEdges[d]);
	        UpperVert.push_back(CoarseEdges[d] + offset);
	        EdgeNgbr.push_back( Range<Dim>(EdgeLength[d].domain() + offset) );
	    }
	}          // end constructor

	//======================================================================
	// Destructor -- QRelations<Traits>::~QRelations
	//======================================================================

	template <class Traits>
	QRelations<Traits>::~QRelations( ) { }

	//======================================================================
	// Function -- QRelations<Traits>::createRelations
	//======================================================================

	template <class Traits>
	void QRelations<Traits>::createRelations( ) {
		ENTER("QRelations::createRelations");
		
        gammaConst = 2.6666666666666667 * 0.25 * quadQ;

        // Dt is an Independent Constant Field used in Relationships
        ScalarField& Dt                = DataDir.get<ScalarField>( "Dt",                 Mesh.getField<ScalarField>    ( Vert   () ) );

        // -------------------------------------------------------------------------------------------------------------------------
        // Input Fields to QRelations
        // -------------------------------------------------------------------------------------------------------------------------
        ScalarField& OldDt             = Old.get<ScalarField>    ( "Dt",                 Mesh.getField<ScalarField>    ( Vert   () ) );
        VectorField& OldVelocity       = Old.get<VectorField>    ( "Velocity" ,          Mesh.getField<VectorField>    ( Vert   () ) );
        VectorField& OldEdgeLength     = Old.get<VectorField>    ( "EdgeLength" ,        FineMesh.getField<VectorField>( AllEdge() ) );
	    VectorField& OldEdgeVelocity   = Old.get<VectorField>    ( "EdgeVelocity",       Mesh.getField<VectorField>    ( AllEdge() ) );
        ScalarField& Volume            = DataDir.get<ScalarField>( "Volume",             Mesh.getField<ScalarField>    ( Cell   () ) );
        ScalarField& SubVolume         = DataDir.get<ScalarField>( "SubVolume",          FineMesh.getField<ScalarField>( Cell   () ) );
        VectorField& SubFaceAreas      = DataDir.get<VectorField>( "SubFaceAreas" ,      FineMesh.getField<VectorField>( AllFace() ) );
        VectorField& EdgeLength        = DataDir.get<VectorField>( "EdgeLength" ,        FineMesh.getField<VectorField>( AllEdge() ) );
        VectorField& Velocity          = DataDir.get<VectorField>( "Velocity",           Mesh.getField<VectorField>    ( Vert   () ) );
        ScalarField& CellDensity       = DataDir.get<ScalarField>( "CellDensity",        Mesh.getField<ScalarField>    ( Cell   () ) );
        ScalarField& CellSoundSpeedSq  = DataDir.get<ScalarField>( "CellSoundSpeedSq",   Mesh.getField<ScalarField>    ( Cell   () ) );
        ScalarField& CellGammaConstant = DataDir.get<ScalarField>( "CellGammaConstant",  Mesh.getField<ScalarField>    ( Cell   () ) );
        CellGammaConstant = gammaConst; // Should go away when we use more sophisticated EOSs.
        // -------------------------------------------------------------------------------------------------------------------------
		// Store the sqrt because it is used to calculate the VertSoundSpeed
        // -------------------------------------------------------------------------------------------------------------------------
        ScalarField& CellSoundSpeed    = DataDir.get<ScalarField>( "CellSoundSpeed",     Mesh.getField<ScalarField>    ( Cell   () ) );
		
        // -------------------------------------------------------------------------------------------------------------------------
		// Vertex Fields -- (weighted) sums of SubCell Fields
        // -------------------------------------------------------------------------------------------------------------------------
	    ScalarField& VertVolume        = DataDir.get<ScalarField>( "VertVolume",         Mesh.getField<ScalarField>    ( Vert   () ) );
//	    ScalarField& WeightedVertVol   = DataDir.get<ScalarField>( "WeightedVertVol",    Mesh.getField<ScalarField>    ( Vert   () ) );
	    ScalarField& VertDensity       = DataDir.get<ScalarField>( "VertDensity",        Mesh.getField<ScalarField>    ( Vert   () ) );
        ScalarField& VertSoundSpeed    = DataDir.get<ScalarField>( "VertSoundSpeed",     Mesh.getField<ScalarField>    ( Vert   () ) );
	    ScalarField& VertGammaConstant = DataDir.get<ScalarField>( "VertGammaConstant",  Mesh.getField<ScalarField>    ( Vert   () ) );

        // -------------------------------------------------------------------------------------------------------------------------
        // Edge Fields that hold Q-related quantities
        // -------------------------------------------------------------------------------------------------------------------------
	    VectorField& EdgeVelocity      = DataDir.get<VectorField>( "EdgeVelocity",       Mesh.getField<VectorField>    ( AllEdge() ) );
	    ScalarField& EdgeDensity       = DataDir.get<ScalarField>( "EdgeDensity",        Mesh.getField<ScalarField>    ( AllEdge() ) );
	    ScalarField& EdgeQPressure     = DataDir.get<ScalarField>( "EdgeQPressure",      Mesh.getField<ScalarField>    ( AllEdge() ) );
	    ScalarField& EdgeSoundSpeed    = DataDir.get<ScalarField>( "EdgeSoundSpeed",     Mesh.getField<ScalarField>    ( AllEdge() ) );
	    ScalarField& EdgeGammaConstant = DataDir.get<ScalarField>( "EdgeGammaConstant",  Mesh.getField<ScalarField>    ( AllEdge() ) );
	    ScalarField& EdgeQTmpMax       = DataDir.get<ScalarField>( "EdgeQTmpMax",        Mesh.getField<ScalarField>    ( AllEdge() ) );
	    ScalarField& SpokeDVolDt       = DataDir.get<ScalarField>( "SpokeDVolDt",        FineMesh.getField<ScalarField>( AllFace() ) );
	    ScalarField& SpokeQSwitch      = DataDir.get<ScalarField>( "SpokeQSwitch",       FineMesh.getField<ScalarField>( AllFace() ) );

	
        ScalarField& RightLimiterRatio = DataDir.get<ScalarField>( "RightLimiterRatio",  Mesh.getField<ScalarField>    ( AllEdge() ) );
        ScalarField& LeftLimiterRatio  = DataDir.get<ScalarField>( "LeftLimiterRatio",   Mesh.getField<ScalarField>    ( AllEdge() ) );
	    ScalarField& EdgePsiLimiter    = DataDir.get<ScalarField>( "EdgePsiLimiter",     Mesh.getField<ScalarField>    ( AllEdge() ) );
	
        // -------------------------------------------------------------------------------------------------------------------------
        // Output Fields from QRelations
        // -------------------------------------------------------------------------------------------------------------------------
        ScalarField& CellQ             = DataDir.get<ScalarField>( "CellQ",              Mesh.getField<ScalarField>    ( Cell   () ) );
        VectorField& SubForceQ         = DataDir.get<VectorField>( "SubForceQ",          FineMesh.getField<VectorField>( Cell   () ) );
        ScalarField& SubPressureQMod   = DataDir.get<ScalarField>( "SubPressureQMod",    FineMesh.getField<ScalarField>( Cell   () ) );

		Let::NewRelation( *(pParent)this, &Parent::calcSqrt,                 CellSoundSpeed,      CellSoundSpeedSq );
		Let::NewRelation( *(pParent)this, &Parent::calcWeightedVertAvg,      VertSoundSpeed,      CellSoundSpeed,    SubVolume,        VertVolume );
		Let::NewRelation( *(pParent)this, &Parent::calcWeightedVertAvg,      VertGammaConstant,   CellGammaConstant, SubVolume,        VertVolume );

		Let::NewRelation( *this, &QRelations<Traits>::calcEdgeVelocity,      EdgeVelocity,        OldVelocity );
		Let::NewRelation( *this, &QRelations<Traits>::calcEdgeDensity,       EdgeDensity,         VertDensity );
		Let::NewRelation( *this, &QRelations<Traits>::calcEdgeGammaConstant, EdgeGammaConstant,   VertGammaConstant );
		Let::NewRelation( *this, &QRelations<Traits>::calcEdgeQPressure,     EdgeQPressure,       EdgeGammaConstant, EdgeSoundSpeed,   EdgeVelocity, EdgePsiLimiter );
		Let::NewRelation( *this, &QRelations<Traits>::calcEdgeSoundSpeed,    EdgeSoundSpeed,      VertSoundSpeed );
		Let::NewRelation( *this, &QRelations<Traits>::calcEdgeQTmpMax,       EdgeQTmpMax,         EdgeDensity,       EdgeQPressure  );
		Let::NewRelation( *this, &QRelations<Traits>::calcSpokeDVolDt,       SpokeDVolDt,         EdgeVelocity,      SubFaceAreas   );
        Let::NewRelation( *this, &QRelations<Traits>::calcSpokeQSwitch,      SpokeQSwitch,        SpokeDVolDt,       Volume,            OldDt );
		
		Let::NewRelation( *this, &QRelations<Traits>::calcRightLimiterRatio, RightLimiterRatio,   OldEdgeLength,     OldEdgeVelocity );
		Let::NewRelation( *this, &QRelations<Traits>::calcLeftLimiterRatio,  LeftLimiterRatio,    OldEdgeLength,     OldEdgeVelocity );
		Let::NewRelation( *this, &QRelations<Traits>::calcEdgePsiLimiter,    EdgePsiLimiter,      RightLimiterRatio, LeftLimiterRatio, OldEdgeLength, OldEdgeVelocity, OldDt );

		Let::NewRelation( *this, &QRelations<Traits>::calcCellQ,             CellQ,               EdgeQTmpMax,       SpokeQSwitch      );
		Let::NewRelation( *this, &QRelations<Traits>::calcSubForceQ,         SubForceQ,           EdgeQTmpMax,       EdgeVelocity,     SpokeDVolDt,   SpokeQSwitch );
		Let::NewRelation( *this, &QRelations<Traits>::calcSubPressureQMod,   SubPressureQMod,     EdgeQPressure,     SpokeQSwitch      );
		
	}
	

	//======================================================================
	// Function -- QRelations<Traits>::calcEdgeVelocity
	//======================================================================

	template <class Traits>
	void QRelations<Traits>::calcEdgeVelocity ( const VectorField&  EdgeVelocity,
	                    	   			        const VectorField&  OldVelocity)
	{
		ENTER("QRelations::calcEdgeVelocity");
		int d = getEdgeDirection(EdgeVelocity);
		EdgeVelocity(CoarseEdges[d]) = OldVelocity(EdgeNgbr[d]) - OldVelocity(CoarseEdges[d]);
	}           // end function calcEdgeVelocity
	

	//======================================================================
	// Function -- QRelations<Traits>::calcEdgeDensity
	//======================================================================

	template <class Traits>
	void QRelations<Traits>::calcEdgeDensity ( const ScalarField&  EdgeDensity,
	                    	   			       const ScalarField&  VertDensity)
	{
		ENTER("QRelations::calcEdgeDensity");
		
		int d = getEdgeDirection(EdgeDensity);
		EdgeDensity(CoarseEdges[d]) = 2.0 * ( VertDensity(CoarseEdges[d]) * VertDensity(EdgeNgbr[d]) ) /
		                                    ( VertDensity(CoarseEdges[d]) + VertDensity(EdgeNgbr[d]) );
	    ScalarField& VertVolume        = DataDir.get<ScalarField>( "VertVolume",         Mesh.getField<ScalarField>    ( Vert   () ) );
	    ScalarField& CellDensity       = DataDir.strictGet<ScalarField>( "CellDensity" );
	}           // end function calcEdgeDensity
	
	//======================================================================
	// Function -- QRelations<Traits>::calcEdgeGammaConstant
	//======================================================================

	template <class Traits>
	void QRelations<Traits>::calcEdgeGammaConstant ( const ScalarField&  EdgeGammaConstant,
	                    	   			             const ScalarField&  VertGammaConstant)
	{
		ENTER("QRelations::calcEdgeGammaConstant");
		
		int d = getEdgeDirection(EdgeGammaConstant);
		EdgeGammaConstant(CoarseEdges[d]) = 0.5 * ( VertGammaConstant(CoarseEdges[d]) + VertGammaConstant(EdgeNgbr[d]) );
	}           // end function calcEdgeGammaConstant
	
	
	//======================================================================
	// Function -- QRelations<Traits>::calcEdgeQPressure
	//======================================================================

	template <class Traits>
	void QRelations<Traits>::calcEdgeQPressure ( const ScalarField&  EdgeQPressure,
	                    	   			         const ScalarField&  EdgeGammaConstant,
	                    	   			         const ScalarField&  EdgeSoundSpeed,
	                    	   			         const VectorField&  EdgeVelocity,
 	                    	   			         const ScalarField&  EdgePsiLimiter )
	{
		ENTER("QRelations::calcEdgeQPressure");
		
		int d = getEdgeDirection(EdgeQPressure);
		EdgeGammaConstant[d].update();
		EdgeSoundSpeed[d].update();
		EdgeVelocity[d].update();
		EdgePsiLimiter[d].update();
		
		ScalarEdgeQPressure<Traits> sEdgeQPressure(linearQ);
        ScalarCode<ScalarEdgeQPressure<Traits> > scEdgeQPressure(sEdgeQPressure);

        scEdgeQPressure(EdgeQPressure,EdgeGammaConstant[d], EdgeSoundSpeed[d],EdgeVelocity[d],EdgePsiLimiter[d]);
	}           // end function EdgeQPressure

	//======================================================================
	// Function -- QRelations<Traits>::calcEdgeSoundSpeed
	//======================================================================

	template <class Traits>
	void QRelations<Traits>::calcEdgeSoundSpeed ( const ScalarField&  EdgeSoundSpeed,
	                    	   			          const ScalarField&  VertSoundSpeed )
	{
		ENTER("QRelations::calcEdgeSoundSpeed");
		
		int d = getEdgeDirection(EdgeSoundSpeed);
		VertSoundSpeed.update();
		EdgeSoundSpeed(CoarseEdges[d]) = min( VertSoundSpeed(CoarseEdges[d]),VertSoundSpeed(EdgeNgbr[d]) );
	}           // end function calcEdgeSoundSpeed


	//======================================================================
	// Function -- QRelations<Traits>::calcEdgeQTmpMax
	//======================================================================

	template <class Traits>
	void QRelations<Traits>::calcEdgeQTmpMax ( const ScalarField&  EdgeQTmpMax,
	                    	   			       const ScalarField&  EdgeDensity,
	                    	   			       const ScalarField&  EdgeQPressure )
	{
		ENTER("QRelations::calcEdgeQTmpMax");
		int d = getEdgeDirection(EdgeQTmpMax);
		EdgeQTmpMax = EdgeDensity[d] * EdgeQPressure[d];
	}           // end function calcEdgeQTmpMax
	
	//======================================================================
	// Function -- QRelations<Traits>::calcSpokeDVolDt
	//======================================================================

	template <class Traits>
	void QRelations<Traits>::calcSpokeDVolDt ( const ScalarField&  SpokeDVolDt,
	                    	   			       const VectorField&  EdgeVelocity,
	                    	   			       const VectorField&  SubFaceAreas)
	{
		ENTER("QRelations::calcSpokeDVolDt");
		int d = getFaceDirection(SpokeDVolDt);
		EdgeVelocity[d].update();
		SubFaceAreas[d].update();
	   	for(int edg=0;edg<nEdgesPerDimension;++edg) {
		    SpokeDVolDt(Spoke[d][edg]) = dot(EdgeVelocity[d](CellEdge[d][edg]),SubFaceAreas[d](Spoke[d][edg]));
		}
	}           // end function calcSpokeDVolDt
	
    //======================================================================
	// Function -- QRelations<Traits>::calcSpokeQSwitch
	//======================================================================

	template <class Traits>
	void QRelations<Traits>::calcSpokeQSwitch ( const ScalarField&  SpokeQSwitch,
	                    	   			   const ScalarField&  SpokeDVolDt,
	                    	   			   const ScalarField&  Volume,
	                    	   			   const ScalarField&  Dt )
	{
		ENTER("QRelations::calcSpokeQSwitch");
		int d = getFaceDirection(SpokeQSwitch);
	   	for(int edg=0;edg<nEdgesPerDimension;++edg) {
		    SpokeQSwitch(Spoke[d][edg]) = (-Dt(CoarseCell) * SpokeDVolDt[d](Spoke[d][edg]))/Volume(CoarseCell);
		}
		SpokeQSwitch = where(SpokeQSwitch > spokeCutoff, 1.0, 0.0);
	}           // end function calcSpokeQSwitch
	
	//======================================================================
	// Function -- QRelations<Traits>::calcRightLimiterRatio
	//======================================================================

	template <class Traits>
    void QRelations<Traits>::calcRightLimiterRatio ( const ScalarField&  RightLimiterRatio,
                                                     const VectorField&  EdgeLength,
                                                     const VectorField&  EdgeVelocity)
	{
		ENTER("QRelations::calcRightLimiterRatio");
		int d = getEdgeDirection(RightLimiterRatio);
        EdgeLength[d].update();
        EdgeVelocity[d].update();
		RightLimiterRatio = 1.0;

		RightLimiterRatio(RDomain[d]) = dot(EdgeLength[d](RightEdgeNgbr[d]),
		                                    EdgeLength[d](RDomain[d]))       *
                                        dot(EdgeVelocity[d](RDomain[d]),
		                                    EdgeVelocity[d](RDomain[d]));

		RightLimiterRatio(RDomain[d]) = where(RightLimiterRatio(RDomain[d]) > numeric_limits<Real>::epsilon(),
                                                (dot(EdgeVelocity[d](RightEdgeNgbr[d]),
                                                     EdgeVelocity[d](RDomain[d]))       *
                                                 dot(EdgeLength[d](RDomain[d]),
                                                     EdgeLength[d](RDomain[d])) ) /
                                               RightLimiterRatio(RDomain[d]),
                                               1.0 );

        //-------------------------------------------------------------
        //  if ( d == 0 ) {
        //  } else if ( d == 1) {
        //     RightLimiterRatio -= 1.0;
        //     tecout << "RightLimiterRatio[1] = " << endl;
        //     tecout << RightLimiterRatio << endl;
        //     RightLimiterRatio += 1.0;
        //  }
        //-------------------------------------------------------------


	}           // end function calcRightLimiterRatio

	//======================================================================
	// Function -- QRelations<Traits>::calcLeftLimiterRatio
	//======================================================================

	template <class Traits>
    void QRelations<Traits>::calcLeftLimiterRatio ( const ScalarField&  LeftLimiterRatio,
                                                    const VectorField&  EdgeLength,
                                                    const VectorField&  EdgeVelocity)
	{
		ENTER("QRelations::calcLeftLimiterRatio");
		int d = getEdgeDirection(LeftLimiterRatio);
        EdgeLength[d].update();
        EdgeVelocity[d].update();
		LeftLimiterRatio = 1.0;

		LeftLimiterRatio(LDomain[d]) = dot(EdgeLength[d](LeftEdgeNgbr[d]),
		                                    EdgeLength[d](LDomain[d]))       *
                                        dot(EdgeVelocity[d](LDomain[d]),
		                                    EdgeVelocity[d](LDomain[d]));

		LeftLimiterRatio(LDomain[d]) = where(LeftLimiterRatio(LDomain[d]) > numeric_limits<Real>::epsilon(),
                                              (dot(EdgeVelocity[d](LeftEdgeNgbr[d]),
                                                    EdgeVelocity[d](LDomain[d]))       *
                                               dot(EdgeLength[d](LDomain[d]),
                                                    EdgeLength[d](LDomain[d])) ) /
                                               LeftLimiterRatio(LDomain[d]),
                                               1.0 );

        //-------------------------------------------------------------
        // if ( d == 0 ) {
        // } else if ( d == 1) {
        //    LeftLimiterRatio -= 1.0;
        //    tecout << "LeftLimiterRatio[1] = " << endl;
        //    tecout << LeftLimiterRatio << endl;
        //    LeftLimiterRatio += 1.0;
        // }
        //-------------------------------------------------------------

	}           // end function calcLeftLimiterRatio

	//======================================================================
	// Function -- QRelations<Traits>::calcEdgePsiLimiter
	//======================================================================

	template <class Traits>
    void QRelations<Traits>::calcEdgePsiLimiter ( const ScalarField&  EdgePsiLimiter,
                                                  const ScalarField&  RightLimiterRatio,
                                                  const ScalarField&  LeftLimiterRatio,
                                                  const VectorField&  EdgeLength,
                                                  const VectorField&  EdgeVelocity,
                                                  const ScalarField&  Dt)
	{
		ENTER("QRelations::calcEdgePsiLimiter");

		int d = getEdgeDirection(EdgePsiLimiter);

        //------------------------------------------------------------------------
        //  Typically the EdgePsiLimiter is limited between 0 and 1,
        //  For this algorithm, 0 means do not have any artificial viscosity
        //  along the edge.
        //
        //  (In Randy Christensen's paper, EdgePsiLimiter is
        //  PsiLimiter = 1.0 - EdgePsiLimiter, so that when he calculates the
        //  Q, there is a factor of ( 1.0 - PsiLimiter) in the expression.  This
        //  means his PsiLimiter has a different meaning than the one used here.
        //  Don M )
        //------------------------------------------------------------------------
        RightLimiterRatio[d].update();
        LeftLimiterRatio[d].update();
        EdgeLength[d].update();
        EdgeVelocity[d].update();
		EdgePsiLimiter = 1.0 - max( 0.0, min( min( 2.0 * RightLimiterRatio[d],
		                                       2.0 * LeftLimiterRatio[d]),
                                          min( 0.5 * (RightLimiterRatio[d] +
		                                              LeftLimiterRatio[d]),
		                                       1.0 ) ) );
		


        //------------------------------------------------------------------------
        // The where statement turns off the limiter when:
        //    1.  The limiter is non zero because of round off
        //    2.  The edge is degenerate (the vertex points coincide)
        //    3.  The change in the velocity along the edge is too small for
        //           anything to happen during the time step.
        //------------------------------------------------------------------------
        EdgePsiLimiter = where( (EdgePsiLimiter > numeric_limits<Real>::epsilon() &&
		                      dot(EdgeLength[d],EdgeLength[d])    >
		                     numeric_limits<Real>::epsilon()             &&
		                      dot(EdgeVelocity[d],EdgeVelocity[d]) * Dt * Dt  >
		                     numeric_limits<Real>::epsilon() *
		                      dot(EdgeLength[d],EdgeLength[d]) ),
		                   EdgePsiLimiter, 0.0 );

//---------------------------------------------------------------------------------------------------------
       if ( d == 0 ) {
//        //   tecout << "eps = " << eps <<endl;
//        //   tecout << "numeric_limits<Real>::epsilon() = " << numeric_limits<Real>::epsilon() << endl;
//           tecout << " Dt = " << Dt << endl;
//           tecout << "EdgePsiLimiter[0] = " << endl;
//           tecout << EdgePsiLimiter << endl;
//           tecout << "EdgeVelocity[0] = " << endl;
//           tecout << EdgeVelocity[0] << endl;
        } else if ( d == 1) {
//        //   tecout << "RightLimiterRatio[1] = " << endl;
//        //   tecout << RightLimiterRatio[1] << endl;
//        //   tecout << "LeftLimiterRatio[1] = " << endl;
//        //   tecout << LeftLimiterRatio[1] << endl;
//           tecout << "EdgePsiLimiter[1] = " << endl;
//           tecout << EdgePsiLimiter << endl;
//           tecout << "EdgeVelocity[1] = " << endl;   
//           tecout << EdgeVelocity[1] << endl;
        }
//---------------------------------------------------------------------------------------------------------

	}           // end function calcEdgePsiLimiter

	//======================================================================
	// Function -- QRelations<Traits>::calcCellQ
	//======================================================================

	template <class Traits>
	void QRelations<Traits>::calcCellQ ( const ScalarField&  CellQ,
	                    	   			 const ScalarField&  EdgeQTmpMax,
                                         const ScalarField&  SpokeQSwitch )
	{
		ENTER("QRelations::calcCellQ");
		EdgeQTmpMax.update();
		SpokeQSwitch.update();
		CellQ = 0.0;
		for ( int d = 0; d < Dim; ++d ) {
		    for ( int edg = 0; edg < nEdgesPerDimension; ++edg ) {
		        CellQ(CoarseCell) = max( CellQ(CoarseCell), EdgeQTmpMax[d](CellEdge[d][edg]) * SpokeQSwitch[d](Spoke[d][edg]) );
		    }
		}

//----------------------------------------
//         tecout << " in calcCellQ: CellQ = " << endl;
//         tecout << CellQ << endl;
//----------------------------------------

	}           // end function calcCellQ

	//======================================================================
	// Function -- QRelations<Traits>::calcSubForceQ
	//======================================================================

	template <class Traits>
	void QRelations<Traits>::calcSubForceQ ( const VectorField&  SubForceQ,
	                    	   			     const ScalarField&  EdgeQTmpMax,
	                    	   			     const VectorField&  EdgeVelocity,
	                    	   			     const ScalarField&  SpokeDVolDt,
                                             const ScalarField&  SpokeQSwitch )
	{
		ENTER("QRelations::calcSubForceQ");
	    SubForceQ = 0.0;
		for ( int d = 0; d < Dim; ++d ) {
		    Loc<Dim> offset(0); offset[d] = 1;
		    for ( int edg = 0; edg < nEdgesPerDimension; ++edg ) {
		        SubForceQ(Spoke[d][edg] - offset) -=(EdgeVelocity[d](CellEdge[d][edg]) *
		                                             SpokeDVolDt[d] (Spoke[d][edg])    *
		                                             EdgeQTmpMax[d] (CellEdge[d][edg]) *
                                                     SpokeQSwitch[d](Spoke[d][edg]) )  /
                                                    (dot(EdgeVelocity[d](CellEdge[d][edg]),EdgeVelocity[d](CellEdge[d][edg]))+eps);

		        SubForceQ(Spoke[d][edg])          +=(EdgeVelocity[d](CellEdge[d][edg]) *
		                                             SpokeDVolDt[d] (Spoke[d][edg])    *
		                                             EdgeQTmpMax[d] (CellEdge[d][edg]) *
                                                     SpokeQSwitch[d](Spoke[d][edg]) )  /
                                                    (dot(EdgeVelocity[d](CellEdge[d][edg]),EdgeVelocity[d](CellEdge[d][edg]))+eps);
		    }
		}

//--------------------------------------------
//         tecout << " in calcSubForceQ: SubForceQ = " << endl;
//         tecout << SubForceQ << endl;
//--------------------------------------------

	}           // end function calcSubForceQ

	//======================================================================
	// Function -- QRelations<Traits>::calcSubPressureQMod
	//======================================================================

	template <class Traits>
	void QRelations<Traits>::calcSubPressureQMod ( const ScalarField&  SubPressureQMod,
	                    	   			           const ScalarField&  EdgeQPressure,
                                                   const ScalarField&  SpokeQSwitch )
	{
		ENTER("QRelations::calcSubPressureQMod");
		
	    SubPressureQMod = 0.0;
		for ( int d = 0; d < Dim; ++d ) {
		    Loc<Dim> offset(0); offset[d] = 1;
		    for ( int edg = 0; edg < nEdgesPerDimension; ++edg ) {
		        SubPressureQMod(Spoke[d][edg] - offset) += EdgeQPressure[d](CellEdge[d][edg]) *
                                                           SpokeQSwitch[d](Spoke[d][edg]);

		        SubPressureQMod(Spoke[d][edg])          += EdgeQPressure[d](CellEdge[d][edg]) *
                                                           SpokeQSwitch[d](Spoke[d][edg]);

		    }
		}
	}           // end function calcSubPressureQMod


}          // end namespace conejo

#endif     // end shroud __conejo_QRelations_t_hh