/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
-------------------------------------------------------------------------------
    Copyright (C) 2011-2016 OpenFOAM Foundation
-------------------------------------------------------------------------------
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 "volgpuFields.H"
#include "surfacegpuFields.H"
#include "gpufvcGrad.H"
#include "coupledFvPatchgpuFields.H"
#include "gpusurfaceInterpolate.H"

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

namespace Foam
{
	template<class Type,class PhiLimiter>
	struct PhiSchemeLimiterFunctor{
		const PhiLimiter& limiter;
		PhiSchemeLimiterFunctor(const PhiLimiter& _limiter): limiter(_limiter) {}
		__host__ __device__
		scalar operator ()(const thrust::tuple<scalar,scalar,Type,Type,vector,scalar>& t){
			return limiter.limiter(thrust::get<0>(t),
			                       thrust::get<1>(t),
			                       thrust::get<2>(t),
			                       thrust::get<3>(t),
			                       thrust::get<4>(t),
			                       thrust::get<5>(t));
		}
	};
}

template<class Type, class PhiLimiter>
Foam::tmp<Foam::surfaceScalargpuField>
Foam::gpuPhiScheme<Type, PhiLimiter>::limiter
(
    const GeometricgpuField<Type, fvPatchgpuField, gpuvolMesh>& phi
) const
{
    const gpufvMesh& mesh = this->mesh();

    tmp<surfaceScalargpuField> tLimiter
    (
        new surfaceScalargpuField
        (
            IOobject
            (
                "PhiLimiter",
                mesh.time().timeName(),
                mesh.hostmesh()
            ),
            mesh,
            dimless
        )
    );
    surfaceScalargpuField& Limiter = tLimiter.ref();

    const surfaceScalargpuField& CDweights = mesh.gpusurfaceInterpolation::weights();

    const surfaceVectorgpuField& Sf = mesh.Sf();
    const surfaceScalargpuField& magSf = mesh.magSf();

    const labelgpuList& owner = mesh.owner();
    const labelgpuList& neighbour = mesh.neighbour();

    tmp<surfaceScalargpuField> tUflux = this->faceFlux_;

    if (this->faceFlux_.dimensions() == dimDensity*dimVelocity*dimArea)
    {
        const volScalargpuField& rho =
            phi.db().objectRegistry::template lookupObject<volScalargpuField>
            ("rho");

        tUflux = this->faceFlux_/fvc::interpolate(rho);
    }
    else if (this->faceFlux_.dimensions() != dimVelocity*dimArea)
    {
        FatalErrorInFunction
            << "dimensions of faceFlux are not correct"
            << exit(FatalError);
    }

    const surfaceScalargpuField& Uflux = tUflux();

    scalargpuField& pLimiter = Limiter.primitiveFieldRef();

    thrust::transform(thrust::make_zip_iterator(thrust::make_tuple(CDweights.field().begin(),
                                                                   Uflux.field().begin(),
                                                                   thrust::make_permutation_iterator(phi.field().begin(),owner.begin()),
                                                                   thrust::make_permutation_iterator(phi.field().begin(),neighbour.begin()),
                                                                   Sf.field().begin(),
                                                                   magSf.field().begin())),
                     thrust::make_zip_iterator(thrust::make_tuple(CDweights.field().end(),
                                                                   Uflux.field().end(),
                                                                   thrust::make_permutation_iterator(phi.field().begin(),owner.end()),
                                                                   thrust::make_permutation_iterator(phi.field().begin(),neighbour.end()),
                                                                   Sf.field().end(),
                                                                   magSf.field().end())),
                      pLimiter.begin(),
                      PhiSchemeLimiterFunctor<Type,PhiLimiter>(*this));

    surfaceScalargpuField::Boundary& bLimiter =
        Limiter.boundaryFieldRef();

    forAll(bLimiter, patchi)
    {
        scalargpuField& pLimiter = bLimiter[patchi];

        if (bLimiter[patchi].coupled())
        {
            const scalargpuField& pCDweights = CDweights.boundaryField()[patchi];
            const vectorgpuField& pSf = Sf.boundaryField()[patchi];
            const scalargpuField& pmagSf = magSf.boundaryField()[patchi];
            const scalargpuField& pFaceFlux = Uflux.boundaryField()[patchi];

            const gpuField<Type> pphiP
            (
                phi.boundaryField()[patchi].patchInternalField()
            );
            const gpuField<Type> pphiN
            (
                phi.boundaryField()[patchi].patchNeighbourField()
            );

            thrust::transform(thrust::make_zip_iterator(thrust::make_tuple(pCDweights.begin(),
                                                                   pFaceFlux.begin(),
                                                                   pphiP.begin(),
                                                                   pphiN.begin(),
                                                                   pSf.begin(),
                                                                   pmagSf.begin())),
                               thrust::make_zip_iterator(thrust::make_tuple(pCDweights.end(),
                                                                   pFaceFlux.end(),
                                                                   pphiP.end(),
                                                                   pphiN.end(),
                                                                   pSf.end(),
                                                                   pmagSf.end())),
                                pLimiter.begin(),
                                   PhiSchemeLimiterFunctor<Type,PhiLimiter>(*this));
        }
        else
        {
            pLimiter = 1.0;
        }
    }

    return tLimiter;
}


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