/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
-------------------------------------------------------------------------------
    Copyright (C) 2011-2016 OpenFOAM Foundation
    Copyright (C) 2019 OpenCFD Ltd.
-------------------------------------------------------------------------------
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 "WallLocalSpringSliderDashpot.H"

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

template<class CloudType>
void Foam::WallLocalSpringSliderDashpot<CloudType>::findMinMaxProperties
(
    scalar& rMin,
    scalar& rhoMax,
    scalar& UMagMax
) const
{
    rMin = VGREAT;
    rhoMax = -VGREAT;
    UMagMax = -VGREAT;

    for (const typename CloudType::parcelType& p : this->owner())
    {
        // Finding minimum diameter to avoid excessive arithmetic

        scalar dEff = p.d();

        if (useEquivalentSize_)
        {
            dEff *= cbrt(p.nParticle()*volumeFactor_);
        }

        rMin = min(dEff, rMin);

        rhoMax = max(p.rho(), rhoMax);

        UMagMax = max
        (
            mag(p.U()) + mag(p.omega())*dEff/2,
            UMagMax
        );
    }

    // Transform the minimum diameter into minimum radius
    //     rMin = dMin/2

    rMin /= 2.0;
}


template<class CloudType>
void Foam::WallLocalSpringSliderDashpot<CloudType>::evaluateWall
(
    typename CloudType::parcelType& p,
    const point& site,
    const WallSiteData<vector>& data,
    scalar pREff,
    bool cohesion
) const
{
    // wall patch index
    label wPI = patchMap_[data.patchIndex()];

    // data for this patch
    scalar Estar = Estar_[wPI];
    scalar Gstar = Gstar_[wPI];
    scalar alpha = alpha_[wPI];
    scalar b = b_[wPI];
    scalar mu = mu_[wPI];
    scalar cohesionEnergyDensity = cohesionEnergyDensity_[wPI];
    cohesion = cohesion && cohesion_[wPI];

    vector r_PW = p.position() - site;

    vector U_PW = p.U() - data.wallData();

    scalar r_PW_mag = mag(r_PW);

    scalar normalOverlapMag = max(pREff - r_PW_mag, 0.0);

    vector rHat_PW = r_PW/(r_PW_mag + VSMALL);

    scalar kN = (4.0/3.0)*sqrt(pREff)*Estar;

    scalar etaN = alpha*sqrt(p.mass()*kN)*pow025(normalOverlapMag);

    vector fN_PW =
        rHat_PW
       *(kN*pow(normalOverlapMag, b) - etaN*(U_PW & rHat_PW));

    // Cohesion force, energy density multiplied by the area of wall/particle
    // overlap
    if (cohesion)
    {
        fN_PW +=
           -cohesionEnergyDensity
           *mathematical::pi*(sqr(pREff) - sqr(r_PW_mag))
           *rHat_PW;
    }

    p.f() += fN_PW;

    vector USlip_PW =
        U_PW - (U_PW & rHat_PW)*rHat_PW
      + (p.omega() ^ (pREff*-rHat_PW));

    scalar deltaT = this->owner().mesh().time().deltaTValue();

    vector& tangentialOverlap_PW =
        p.collisionRecords().matchWallRecord(-r_PW, pREff).collisionData();

    tangentialOverlap_PW += USlip_PW*deltaT;

    scalar tangentialOverlapMag = mag(tangentialOverlap_PW);

    if (tangentialOverlapMag > VSMALL)
    {
        scalar kT = 8.0*sqrt(pREff*normalOverlapMag)*Gstar;

        scalar etaT = etaN;

        // Tangential force
        vector fT_PW;

        if (kT*tangentialOverlapMag > mu*mag(fN_PW))
        {
            // Tangential force greater than sliding friction,
            // particle slips

            fT_PW = -mu*mag(fN_PW)*USlip_PW/mag(USlip_PW);

            tangentialOverlap_PW = Zero;
        }
        else
        {
            fT_PW =
                -kT*tangentialOverlapMag
               *tangentialOverlap_PW/tangentialOverlapMag
              - etaT*USlip_PW;
        }

        p.f() += fT_PW;

        p.torque() += (pREff*-rHat_PW) ^ fT_PW;
    }
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

template<class CloudType>
Foam::WallLocalSpringSliderDashpot<CloudType>::WallLocalSpringSliderDashpot
(
    const dictionary& dict,
    CloudType& cloud
)
:
    WallModel<CloudType>(dict, cloud, typeName),
    Estar_(),
    Gstar_(),
    alpha_(),
    b_(),
    mu_(),
    cohesionEnergyDensity_(),
    cohesion_(),
    patchMap_(),
    maxEstarIndex_(-1),
    collisionResolutionSteps_
    (
        this->coeffDict().getScalar("collisionResolutionSteps")
    ),
    volumeFactor_(1.0),
    useEquivalentSize_(Switch(this->coeffDict().lookup("useEquivalentSize")))
{
    if (useEquivalentSize_)
    {
        this->coeffDict().readEntry("volumeFactor", volumeFactor_);
    }

    scalar pNu = this->owner().constProps().poissonsRatio();

    scalar pE = this->owner().constProps().youngsModulus();

    const polyMesh& mesh = cloud.mesh();

    const polyBoundaryMesh& bMesh = mesh.boundaryMesh();

    patchMap_.setSize(bMesh.size(), -1);

    DynamicList<label> wallPatchIndices;

    for (const polyPatch& pp : bMesh)
    {
        if (isA<wallPolyPatch>(pp))
        {
            wallPatchIndices.append(pp.index());
        }
    }

    label nWallPatches = wallPatchIndices.size();

    Estar_.setSize(nWallPatches);
    Gstar_.setSize(nWallPatches);
    alpha_.setSize(nWallPatches);
    b_.setSize(nWallPatches);
    mu_.setSize(nWallPatches);
    cohesionEnergyDensity_.setSize(nWallPatches);
    cohesion_.setSize(nWallPatches);

    scalar maxEstar = -GREAT;

    forAll(wallPatchIndices, wPI)
    {
        const dictionary& patchCoeffDict
        (
            this->coeffDict().subDict(bMesh[wallPatchIndices[wPI]].name())
        );

        patchMap_[wallPatchIndices[wPI]] = wPI;

        scalar nu = patchCoeffDict.get<scalar>("poissonsRatio");

        scalar E = patchCoeffDict.get<scalar>("youngsModulus");

        Estar_[wPI] = 1/((1 - sqr(pNu))/pE + (1 - sqr(nu))/E);

        Gstar_[wPI] = 1/(2*((2 + pNu - sqr(pNu))/pE + (2 + nu - sqr(nu))/E));

        alpha_[wPI] = patchCoeffDict.get<scalar>("alpha");

        b_[wPI] = patchCoeffDict.get<scalar>("b");

        mu_[wPI] = patchCoeffDict.get<scalar>("mu");

        cohesionEnergyDensity_[wPI] =
            patchCoeffDict.get<scalar>("cohesionEnergyDensity");

        cohesion_[wPI] = (mag(cohesionEnergyDensity_[wPI]) > VSMALL);

        if (Estar_[wPI] > maxEstar)
        {
            maxEstarIndex_ = wPI;

            maxEstar = Estar_[wPI];
        }
    }
}


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

template<class CloudType>
Foam::scalar Foam::WallLocalSpringSliderDashpot<CloudType>::pREff
(
    const typename CloudType::parcelType& p
) const
{
    if (useEquivalentSize_)
    {
        return p.d()/2*cbrt(p.nParticle()*volumeFactor_);
    }

    return p.d()/2;
}


template<class CloudType>
bool Foam::WallLocalSpringSliderDashpot<CloudType>::controlsTimestep() const
{
    return true;
}


template<class CloudType>
Foam::label Foam::WallLocalSpringSliderDashpot<CloudType>::nSubCycles() const
{
    if (!(this->owner().size()))
    {
        return 1;
    }

    scalar rMin;
    scalar rhoMax;
    scalar UMagMax;

    findMinMaxProperties(rMin, rhoMax, UMagMax);

    // Note:  pi^(7/5)*(5/4)^(2/5) = 5.429675
    scalar minCollisionDeltaT =
        5.429675
       *rMin
       *pow(rhoMax/(Estar_[maxEstarIndex_]*sqrt(UMagMax) + VSMALL), 0.4)
       /collisionResolutionSteps_;

    return ceil(this->owner().time().deltaTValue()/minCollisionDeltaT);
}


template<class CloudType>
void Foam::WallLocalSpringSliderDashpot<CloudType>::evaluateWall
(
    typename CloudType::parcelType& p,
    const List<point>& flatSitePoints,
    const List<WallSiteData<vector>>& flatSiteData,
    const List<point>& sharpSitePoints,
    const List<WallSiteData<vector>>& sharpSiteData
) const
{
    scalar pREff = this->pREff(p);

    forAll(flatSitePoints, siteI)
    {
        evaluateWall
        (
            p,
            flatSitePoints[siteI],
            flatSiteData[siteI],
            pREff,
            true
        );
    }

    forAll(sharpSitePoints, siteI)
    {
        // Treating sharp sites like flat sites, except suppress cohesion

        evaluateWall
        (
            p,
            sharpSitePoints[siteI],
            sharpSiteData[siteI],
            pREff,
            false
        );
    }
}


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