//-----------------------------------------------------------------------------
// NastDiffQuot2dDonorCell.cpp
//-----------------------------------------------------------------------------
//
//  Copyright (C) 1998 Technische Universitaet Muenchen, Germany
//  written by Bernhard Brueck
//
//  This file is part of Nast++
//
//-----------------------------------------------------------------------------
// CNastDiffQuot2dDonorCell ist eine Implementierung der Differenzenquotienten
// mit 
//-----------------------------------------------------------------------------
//  Changes:
//

#include "NastDiffQuot2dDonorCell.h"
#include "NastDebug.h"
#include "NastStaggeredGrid2d.h"
#include "NastGeometry2d.h"


//-----------------------------------------------------------------------------
//                    Konstruktor + Destruktor
//-----------------------------------------------------------------------------

//  Konstruktor
CNastDiffQuot2dDonorCell::CNastDiffQuot2dDonorCell( double gamma )
    :m_rdx (0),
     m_rdy (0),
     m_rdx2(0),
     m_rdy2(0),
     m_r4dx(0),
     m_r4dy(0),
     m_gamma(gamma),
     m_pGrid(0)
{
    NAST_ASSERT( gamma >= 0 );
}

//  Destruktor
CNastDiffQuot2dDonorCell::~CNastDiffQuot2dDonorCell()
{
    NAST_ASSERT_VALID( this );
    // auf Gitter und Geometrie sind nur Referenzen vorhanden,
    // und ansonsten gibt es nichts zum Freigeben
}

//-------------------------------------------------------------------------
//                          Memberfunktionen
//-------------------------------------------------------------------------
void 
CNastDiffQuot2dDonorCell::useGrid( CNastStaggeredGrid2d &grid )
{
    NAST_ASSERT_VALID( &grid );
    
    // die Konstanten berechnen
    m_rdx =  1 / grid.dx();
    m_rdy =  1 / grid.dy();
    m_rdx2 = 1 / (grid.dx() * grid.dx());
    m_rdy2 = 1 / (grid.dy() * grid.dy());
    m_r4dx = 1 / (4 * grid.dx());
    m_r4dy = 1 / (4 * grid.dy());

    // und das Gitter merken
    m_pGrid = &grid;
}


//-------------------------------------------------------------------------
//                          Diskretisierung der Werte
//-------------------------------------------------------------------------
//
//    |                     |                     |                     |
//  --+---------------------+---------------------+---------------------+--
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//  U(i-1,j+1)          U(i-1,j+1)                |                     |
//    |                     |                     |                     | 
//    |                     |                     |                     | 
//    |                     |                     |                     |
//    |                     |                     |                     |
//  --+---------------------+-------V(i,j+1)------+---------------------+--
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//  U(i-1,j)             U(i,j)      i,j      U(i+1,j)                  |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//  --+-------V(i-1,j)------+--------V(i,j)-------+-------V(i+1,j)------+--     
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                 U(i-1,j)               U(i+1,j-1)               |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//  --+---------------------+-------V(i,j-1)------+------V(i+1,j-1)-----+-- 
//    |                     |                     |                     |         

//      ^ 
//      |  
//      Y                            CNastCell2d::NORTH
//      |   
//  ----+--X-->    CNastCell2d::WEST  CNastCell2d::THIS   CNastCell2d::EAST
//      |               
//      |                             CNastCell2d::SOUTH


//-------------------------------------------------------------------------
// Differenzenquotienten
//-------------------------------------------------------------------------

//  --+---------------------+---------------------+---------------------+--
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     | 
//    |                     |                     |                     | 
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//  --+---------------------+---------------------+---------------------+--
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     U         i,j         U                     |
//    |                     |                     |                     | 
//    |                     |                     |                     | 
//    |                     |                     |                     |
//    |                     |         ***         |                     |
//  --+----------V----------U'--------*V*---------U'---------V----------+--     
//    |         v_w         |         ***         |         v_e         |
//    |                     |         v_m         |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     U                     U                     |     
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//  --+---------------------+---------------------+---------------------+-- 
//
//-----------------------------------------------------------------------------
// !! Achtung die Indizierung in den Zellen ist geandert
//-----------------------------------------------------------------------------
// DUVDX = ((U[i][j]+U[i][j+1])*(V[i][j]+V[i+1][j])+
// 	 m_gamma*fabs(U[i][j]+U[i][j+1])*(V[i][j]-V[i+1][j])-
// 	 (U[i-1][j]+U[i-1][j+1])*(V[i-1][j]+V[i][j])-
// 	 m_gamma*fabs(U[i-1][j]+U[i-1][j+1])*(V[i-1][j]-V[i][j]))
//      /(4.0*delx);
//-----------------------------------------------------------------------------
//
//  Der Differenzenquotient wird nur an den Kanten zwischen zwei Fluidzellen
//  berechnet. Wenn einer der beiden seitlichen V-Werte ausserhalb des 
//  Fluids liegt wird der Wert von CNastStaggeredGrid2d durch eine Anfrage
//  an die Geometrie besimmt. Deshalb darf auf die Werte nicht direkt
//  zugegriffen werden, sondern sie muessen ueber westV, eastV erfragt werden.
//

double CNastDiffQuot2dDonorCell::diff_duv_dx( int i, int j ) const
{ 
    //  nur an Kante zwischen zwei Fluidzellen
    NAST_ASSERT( grid()->cell(i,j).isFluid( CNastCell2d::THIS | 
        				    CNastCell2d::SOUTH ));

    const double v_w = grid()->westV( i, j );	   // westlicher V-Wert
    const double v_m = grid()->V    ( i, j );	   // mittlerer  V-Wert
    const double v_e = grid()->eastV( i, j );	   // oestlicher V-Wert

    const double uSum1  = grid()->U( i+1, j ) + grid()->U( i+1, j-1 );
    const double uSum2  = grid()->U( i  , j ) + grid()->U( i  , j-1 );

    const double vSum1  = v_m + v_e;
    const double vDiff1 = v_m - v_e;

    const double vSum2  = v_w + v_m;
    const double vDiff2 = v_w - v_m;
    
    return ( uSum1 * vSum1 + fabs(uSum1) * vDiff1 * m_gamma -
             uSum2 * vSum2 - fabs(uSum2) * vDiff2 * m_gamma) * m_r4dx;
}


//                          |                     |                     |
//  --+---------------------+---------------------+---------------------+--
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                    U u_n                  |                     | 
//    |                     |                     |                     | 
//    |                     |                     |                     | 
//    |                     |                     |                     |
//    |                     |                     |                     |
//    +----------V----------V'---------V----------+---------------------+-     
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                    ***                    |                     |
//    |                    *U* u_m    i,j         |                     |
//    |                    ***                    |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//  --+----------V----------V'---------V----------+---------------------+-
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     U u_s                 |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//  --+---------------------+---------------------+---------------------+--
//    |                     |                     |                     |         
//
//-----------------------------------------------------------------------------
// !! Achtung die Indizierung in den Zellen ist geandert
//-----------------------------------------------------------------------------
// DUVDY = ((V[i][j]+V[i+1][j])*(U[i][j]+U[i][j+1])+
// 	 m_gamma*fabs(V[i][j]+V[i+1][j])*(U[i][j]-U[i][j+1])-
// 	 (V[i][j-1]+V[i+1][j-1])*(U[i][j-1]+U[i][j])-
// 	 m_gamma*fabs(V[i][j-1]+V[i+1][j-1])*(U[i][j-1]-U[i][j]))
//                   /(4.0*dely);
//-----------------------------------------------------------------------------
//
//  Hier ist es moeglich, das die beiden aeusseren U Werte ausserhalb des 
//  Fluids liegen
//
double CNastDiffQuot2dDonorCell::diff_duv_dy( int i, int j ) const
{
    NAST_ASSERT( grid()->cell(i,j).isFluid( CNastCell2d::THIS | 
        				    CNastCell2d::WEST ));

    const double u_n = grid()->northU( i, j ); // noerdlicher U-Wort
    const double u_m = grid()->U     ( i, j ); // mittlerer U-Wert
    const double u_s = grid()->southU( i, j ); // suedlicher U-Wert

    const double vSum1  = grid()->V( i-1, j+1 ) + grid()->V( i, j+1 );
    const double vSum2  = grid()->V( i-1, j   ) + grid()->V( i, j   );

    const double uSum1  = u_m + u_n;
    const double uDiff1 = u_m - u_n;

    const double uSum2  = u_s + u_m;
    const double uDiff2 = u_s - u_m;

    return ( vSum1 * uSum1 + fabs(vSum1) * uDiff1 * m_gamma -
             vSum2 * uSum2 - fabs(vSum2) * uDiff2 * m_gamma) * m_r4dy;
}


//    |                     |                     |                     |
//  --+---------------------+---------------------+---------------------+--
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     | 
//    |                     |                     |                     | 
//    |                     |                     |                     |
//    |                     |                     |                     |
//  --+---------------------+---------------------+---------------------+--     
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                    ***                    |                     |
//    U                    *U*        i,j         U                     |
//    |                    ***                    |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//  --+---------------------+---------------------+---------------------+-- 
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//  --+---------------------+---------------------+---------------------+--
//    |                     |                     |                     |         
//
//-----------------------------------------------------------------------------
// !! Achtung die Indizierung in den Zellen ist geandert
//-----------------------------------------------------------------------------
// DU2DX = ((U[i][j]+U[i+1][j])*(U[i][j]+U[i+1][j])+
// 	 m_gamma*fabs(U[i][j]+U[i+1][j])*(U[i][j]-U[i+1][j])-
// 	 (U[i-1][j]+U[i][j])*(U[i-1][j]+U[i][j])-
// 	 m_gamma*fabs(U[i-1][j]+U[i][j])*(U[i-1][j]-U[i][j]))
//     /(4.0*delx);
//-----------------------------------------------------------------------------

double CNastDiffQuot2dDonorCell::diff_du2_dx( int i, int j ) const
{
    NAST_ASSERT( grid()->cell(i,j).isFluid( CNastCell2d::THIS | 
        				   CNastCell2d::WEST ));

    double uSum1  = grid()->U( i,   j ) + grid()->U( i+1, j );
    double uDiff1 = grid()->U( i,   j ) - grid()->U( i+1, j );
    double uSum2  = grid()->U( i-1, j ) + grid()->U( i,   j );
    double uDiff2 = grid()->U( i-1, j ) - grid()->U( i,   j );
    
    return ( uSum1 * uSum1 + fabs(uSum1) * uDiff1 * m_gamma -
             uSum2 * uSum2 + fabs(uSum2) * uDiff2 * m_gamma ) * m_r4dx;
}


//    |                     |                     |                     |
//  --+---------------------+---------------------+---------------------+--
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//  --+---------------------+----------V----------+---------------------+--
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |         i,j         |                     |     
//    |                     |                     |                     |
//    |                     |                     |                     | 
//    |                     |                     |                     | 
//    |                     |                     |                     |
//    |                     |         ***         |                     |
//  --+---------------------+---------*V*---------+---------------------+--     
//    |                     |         ***         |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//  --+---------------------+----------V----------+---------------------+-- 
//    |                     |                     |                     |         
//
//-----------------------------------------------------------------------------
// !! Achtung die Indizierung in den Zellen ist geandert
//-----------------------------------------------------------------------------
// DV2DY = ((V[i][j]+V[i][j+1])*(V[i][j]+V[i][j+1])+
// 	 m_gamma*fabs(V[i][j]+V[i][j+1])*(V[i][j]-V[i][j+1])-
// 	 (V[i][j-1]+V[i][j])*(V[i][j-1]+V[i][j])-
// 	 m_gamma*fabs(V[i][j-1]+V[i][j])*(V[i][j-1]-V[i][j]))
//     /(4.0*dely);
//-----------------------------------------------------------------------------

double CNastDiffQuot2dDonorCell::diff_dv2_dy( int i, int j ) const
{
    NAST_ASSERT( grid()->cell(i,j).isFluid( CNastCell2d::THIS | 
        				   CNastCell2d::SOUTH ));
        
    double vSum1  = grid()->V( i, j ) + grid()->V( i, j+1 );
    double vDiff1 = grid()->V( i, j ) - grid()->V( i, j+1 );
    
    double vSum2  = grid()->V( i, j-1 ) + grid()->V( i, j );
    double vDiff2 = grid()->V( i, j-1 ) - grid()->V( i, j );
    
    return ( vSum1 * vSum1 + fabs(vSum1) * vDiff1 * m_gamma -
             vSum2 * vSum2 + fabs(vSum2) * vDiff2 * m_gamma ) * m_r4dy;
}


//    |                     |                     |                     |
//  --+---------------------+---------------------+---------------------+--
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     U                     |
//    |                     |                     |                     | 
//    |                     |                     |                     | 
//    |                     |                     |                     |
//    |                     |                     |                     |
//  --+---------------------+---------------------+---------------------+--     
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                    ***                    |                     |
//    U                    *U*       i,j          U                     |
//    |                    ***                    |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//  --+---------------------+---------------------+---------------------+-- 
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     U                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//  --+---------------------+---------------------+---------------------+--
//    |                     |                     |                     |         
//
//-----------------------------------------------------------------------------
// LAPLU = (U[i+1][j]-2.0*U[i][j]+U[i-1][j])/delx/delx+
//         (U[i][j+1]-2.0*U[i][j]+U[i][j-1])/dely/dely;
//-----------------------------------------------------------------------------

double CNastDiffQuot2dDonorCell::diff_lapl_u( int i, int j ) const
{
    NAST_ASSERT( grid()->cell(i,j).isFluid( CNastCell2d::THIS | 
        				    CNastCell2d::WEST ));

    const double u_n = grid()->northU( i, j );	// noerdlicher U-Wert
    const double u_m = grid()->U     ( i, j );	// mittlerer U-Wert
    const double u_s = grid()->southU( i, j );	// suedlicher U-Wert

    return
        (grid()->U(i-1, j ) - 2*u_m  + grid()->U(i+1, j )) * m_rdx2 +
        (       u_n         - 2*u_m  +         u_s       ) * m_rdy2;
}

//    |                     |                     |                     |
//  --+---------------------+---------------------+---------------------+--
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//  --+---------------------+----------V----------+---------------------+--
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |        i,j          |                     |
//    |                     |                     |                     | 
//    |                     |                     |                     | 
//    |                     |                     |                     |
//    |                     |         ***         |                     |
//  --+----------V----------+---------*V*---------+----------V----------+--     
//    |                     |         ***         |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |     
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//    |                     |                     |                     |
//  --+---------------------+----------V----------+---------------------+-- 
//    |                     |                     |                     |
//
//-----------------------------------------------------------------------------
// LAPLV = (V[i+1][j]-2.0*V[i][j]+V[i-1][j])/delx/delx+
// (V[i][j+1]-2.0*V[i][j]+V[i][j-1])/dely/dely;
//-----------------------------------------------------------------------------

double CNastDiffQuot2dDonorCell::diff_lapl_v( int i, int j ) const
{
    NAST_ASSERT( grid()->cell(i,j).isFluid( CNastCell2d::THIS | 
        				    CNastCell2d::SOUTH ));

    const double v_w = grid()->westV( i, j ); // westlicher V-Wert
    const double v_m = grid()->V    ( i, j ); // mittlerer V-Wert
    const double v_e = grid()->eastV( i, j ); // suedlicher V-Wert

    return
        (        v_w        -  2 * v_m  +       v_e       ) * m_rdx2 +
        (grid()->V(i, j-1 ) -  2 * v_m  + grid()->V(i,j+1)) * m_rdy2;
}



// ----------------------------------------------------------------------------
//                                    Debug
// ----------------------------------------------------------------------------

void 
CNastDiffQuot2dDonorCell::debugDump( CNastDumpContext &dumpContext ) const
{
    CNastDiffQuot2d::debugDump( dumpContext );
    dumpContext << "CNastDiffQuot2dDonorCell" << "\n";
    dumpContext << "\tgamma = " << m_gamma << "\n";
}

void 
CNastDiffQuot2dDonorCell::assertValid() const
{
    // zuerst einmal AssertValid der Basisklasse aufrufen
    CNastDiffQuot2d::assertValid(); 

    //  falls noch kein Gitter vorhanden ist braucht nichts
    //  weiter getestet werden. Ansonsten muessen die Konstanten
    //  stimmen
    if( m_pGrid )
    {
        NAST_ASSERT( m_pGrid );	// ist ueberhaupt schon ein Gitter angegeben ?
        
        // sind die Konstanten alle gesetzt
        
        NAST_ASSERT( m_rdx  > 0);
        NAST_ASSERT( m_rdy  > 0);
        NAST_ASSERT( m_rdx2 > 0);
        NAST_ASSERT( m_rdy2 > 0);
        NAST_ASSERT( m_r4dx > 0);
        NAST_ASSERT( m_r4dy > 0);
    }
}


