/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "solidParticleCloud.H"

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

namespace Foam
{
    defineTemplateTypeNameAndDebug(Cloud<solidParticle>, 0);
}

// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //


bool Foam::solidParticle::move( trackingData& td,
                                const scalar trackTime
                              )
{
    td.switchProcessor = false;
    td.keepParticle = true;

    const polyBoundaryMesh& pbMesh = mesh_.boundaryMesh();

    scalar tEnd = (1.0 - stepFraction())*trackTime;
    scalar dtMax = tEnd;

    while (td.keepParticle && !td.switchProcessor && tEnd > SMALL)
    {
        if (debug)
        {
            Info<< "Time = " << mesh_.time().timeName()
                << " trackTime = " << trackTime
                << " tEnd = " << tEnd
                << " steptFraction() = " << stepFraction() << endl;
        }

        // set the lagrangian time-step
        scalar dt = min(dtMax, tEnd);

        // remember which cell the parcel is in
        // since this will change if a face is hit
        label cellI = cell();
  
        //- Track particle to a given position and returns 1.0 if the
        //  trajectory is completed without hitting a face otherwise
        //  stops at the face and returns the fraction of the trajectory
        //  completed.
        //  on entry 'stepFraction()' should be set to the fraction of the
        //  time-step at which the tracking starts.
//        Pout << ".";
        dt *= trackToFace(position() + dt*U_, td);

/*        if(dt<=0) 
        {
            Pout << "S";
        }*/

        tEnd -= dt;
        stepFraction() = 1.0 - tEnd/trackTime;
  
        cellPointWeight cpw(mesh_, position(), cellI, face());
        vector Uc = td.UInterp().interpolate(cpw);
      
        // Add in here the interpolation for p, p_rgh, alpha, which we use in ParaFoam to colour the dots
        p_ = td.pInterp().interpolate(cpw);
        p_rgh_ = td.p_rghInterp().interpolate(cpw);
        alpha_ = td.alphaInterp().interpolate(cpw);
        // -----------
        scalar rhoc = td.rhoInterp().interpolate(cpw);    // density of fluid surrounding 
        scalar rhop  = td.cloud().rhop();
   
        // Lets assume that if the particle are smaller than zero, then we merely want to track the flow
        // So we'll just skip the following
        if(d_>0)
        {
        	if (debug) Info << "Real Particles" <<endl;
            scalar nuc  = td.nuInterp().interpolate(cpw);

            // Seems to be based on: http://www.freestudy.co.uk/fluid%20mechanics/saet3a.pdf
            scalar magUr = mag(Uc - U_);
            
            scalar ReFunc = 1.0;
            scalar Re = magUr*d_/nuc;   // where nu is kinematic viscosity

            if (Re > 0.01)
            {
                // This looks like a coefficient from "Clift et al, 1978", where they give
                // Cd = 24 * (1+0.15 Re^0.687) / Re
                ReFunc = 1.0 + 0.15*pow(Re, 0.687);
            }

            // Compute the drag coefficient
            scalar Dc = (24.0*nuc/d_)*ReFunc*(3.0/4.0)*(rhoc/(d_*rhop));

            // Compute new velocity of particle. Initial velocity + acceleration from drag
            // This is a rearrangement to get U1, for
            //   U1 = U0 + dt*Dc*(Uc-U1) + dt*(1-rhoc/rhop)*g

            U_ = (U_ + dt*(Dc*Uc + (1.0 - rhoc/rhop)*td.g()))/(1.0 + dt*Dc);
        }
        else
        {
        	if (debug) Info << "Trace particles" <<endl;

            // As we are merely tracking the stream flow, but we want to retain acceleration due to 
            // gravity on particle as they go into the air. But I'm assuming these particle reach terminal 
            // velocity immediately, i.e. there is no retained acceleration from a previous time step. 
            U_ = Uc + (1.0 - rhoc/rhop)*td.g()*dt;
        }

        if (onBoundary() && td.keepParticle)
        {
            if (isA<processorPolyPatch>(pbMesh[patch(face())]))
            {
                td.switchProcessor = true;
            }
            else 
            {
                // Particle has reached the edge of the domain, so it must vanish
                td.keepParticle=false;
            }
        }
    }

    return td.keepParticle;
}


bool Foam::solidParticle::hitPatch
(
    const polyPatch&,
    trackingData&,
    const label,
    const scalar,
    const tetIndices&
)
{
    return false;
}


void Foam::solidParticle::hitProcessorPatch
(
    const processorPolyPatch&,
    trackingData& td
)
{
    td.switchProcessor = true;
}


void Foam::solidParticle::hitWallPatch(
                                        const wallPolyPatch& wpp,
                                        trackingData& td,
                                        const tetIndices& tetIs
                                      )
{
    vector nw = tetIs.faceTri(mesh_).normal();
    nw /= mag(nw);

    scalar Un = U_ & nw;
    vector Ut = U_ - Un*nw;

    if (Un > 0)
    {
        U_ -= (1.0 + td.cloud().e())*Un*nw;
    }

    U_ -= td.cloud().mu()*Ut;
}


void Foam::solidParticle::hitPatch( const polyPatch&,
                                    trackingData& td )
{
    td.keepParticle = false;
}


void Foam::solidParticle::transformProperties (const tensor& T)
{
    particle::transformProperties(T);
    U_ = transform(T, U_);
}


void Foam::solidParticle::transformProperties(const vector& separation)
{
    particle::transformProperties(separation);
}


Foam::scalar Foam::solidParticle::wallImpactDistance(const vector&) const
{
    if(d_<=0) return 0.0005;  // Some small number
    return 0.5*d_;
}


// ************************************************************************* //