/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  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/>.

Application
    test

Description
    Finite volume method test code.

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

#include "fvCFD.H"
#include "syncTools.H"
#include "Random.H"

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

template<class Type>
void checkSwap
(
    const GeometricField<Type, fvPatchField, volMesh>& fld
)
{
    const fvMesh& mesh = fld.mesh();

    List<Type> bVals(mesh.nFaces()-mesh.nInternalFaces(), pTraits<Type>::zero);
    forAll(fld.boundaryField(), patchI)
    {
        const polyPatch& pp = mesh.boundaryMesh()[patchI];
        const fvPatchField<Type>& pfld = fld.boundaryField()[patchI];

        const tmp<Field<Type> > pintFld = pfld.patchInternalField();

        forAll(pfld, i)
        {
            if (pp.coupled())
            {
                label meshFaceI = pp.start()+i;
                bVals[meshFaceI-mesh.nInternalFaces()] = pintFld()[i];
            }
        }
    }

    syncTools::swapBoundaryFaceList(mesh, bVals);

    forAll(fld.boundaryField(), patchI)
    {
        const polyPatch& pp = mesh.boundaryMesh()[patchI];
        const fvPatchField<Type>& pfld = fld.boundaryField()[patchI];

        Pout<< "Patch:" << pp.name() << endl;

        forAll(pfld, i)
        {
            if (pp.coupled())
            {
                label meshFaceI = pp.start()+i;
                const Type& bVal = bVals[meshFaceI-mesh.nInternalFaces()];

                if (bVal != pfld[i])
                {
                    Pout<< "face:" << meshFaceI << nl
                        << "    val:" << pfld[i] << nl
                        << "    swp:" << bVal << nl
                        << endl;
                }
            }
        }
    }
}


int main(int argc, char *argv[])
{

#   include "setRootCase.H"

#   include "createTime.H"
#   include "createMesh.H"

    volScalarField fx("fx", pow(mesh.C().component(vector::X), 2));
    fx.write();


    // Get randomised velocity field

    volVectorField U
    (
        IOobject
        (
            "U2",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        mesh,
        dimensionedVector("one", dimVelocity, vector::zero)
    );

    Random rndGen(0);
    forAll(U.internalField(), cellI)
    {
        rndGen.randomise(U[cellI]);
        U[cellI] -= vector(0.5, 0.5, 0.5);
    }
    U.correctBoundaryConditions();
    U.write();


    // Calculate flux

    surfaceScalarField phi
    (
        IOobject
        (
            "phi",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        linearInterpolate(U) & mesh.Sf()
    );


    // Calculate divergence

    volScalarField fld(fvc::div(phi, fx));
    checkSwap(fld);



//    volScalarField gradx4(fvc::grad(fx)().component(vector::X));
//    gradx4.write();

    //volVectorField curlC(fvc::curl(1.0*mesh.C()));
    //curlC.write();

    /*
    surfaceScalarField xf(mesh.Cf().component(vector::X));
    surfaceScalarField xf4(pow(xf, 4));

    for (int i=1; i<xf4.size()-1; i++)
    {
        scalar gradx4a = (xf4[i] - xf4[i-1])/(xf[i] - xf[i-1]);
        Info<< (gradx4a - gradx4[i])/gradx4a << endl;
    }
    */

    Info<< "end" << endl;
}


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

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  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/>.

Application
    scalarTransportFoam

Description
    Solves a transport equation for a passive scalar

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

#include "fvCFD.H"
#include "simpleControl.H"

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

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"
    #include "createMesh.H"
    #include "createFields.H"

    simpleControl simple(mesh);

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

    Info<< "\nCalculating scalar transport\n" << endl;

    #include "CourantNo.H"

    while (simple.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        while (simple.correctNonOrthogonal())
        {
            {
                tmp<volVectorField> tgradT(fvc::grad(T));
                volVectorField& gradT = tgradT();
                const volVectorField::GeometricBoundaryField& bfld =
                    gradT.boundaryField();

                Pout<< "--- Coupled patchFields before ---" << endl;
                forAll(bfld, patchI)
                {
                    if (bfld[patchI].coupled())
                    {
                        Pout<< "Patch:" << bfld[patchI].patch().name() << nl
                            << "internal:" << bfld[patchI].patchInternalField()
                            << nl
                            << "neighb:" << bfld[patchI].patchNeighbourField()
                            << endl;
                    }
                }




                gradT.correctBoundaryConditions();

                Pout<< "--- Coupled patchFields after ---" << endl;
                forAll(bfld, patchI)
                {
                    if (bfld[patchI].coupled())
                    {
                        Pout<< "Patch:" << bfld[patchI].patch().name() << nl
                            << "internal:" << bfld[patchI].patchInternalField()
                            << nl
                            << "neighb:" << bfld[patchI].patchNeighbourField()
                            << endl;
                    }
                }
                Pout<< nl << endl;
            }


            solve
            (
                fvm::ddt(T)
              + fvm::div(phi, T)
              - fvm::laplacian(DT, T)
            );
        }

        runTime.write();
    }

    Info<< "End\n" << endl;

    return 0;
}


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