/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 1991-2005 OpenCFD Ltd.
     \\/     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 2 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, write to the Free Software Foundation,
    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA

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

#include "extrudedMesh.H"
#include "wallPolyPatch.H"
#include "meshTools.H"
#include "ListOps.H"
#include "OFstream.H"

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

bool Foam::extrudedMesh::sameOrder(const face& f, const edge& e)
{
    label i = findIndex(f, e[0]);

    label nextI = (i == f.size()-1 ? 0 : i+1);

    return f[nextI] == e[1];
}


template
<
    class Face,
    template<class> class FaceList,
    class PointField,
    class PointExtruder
>
Foam::pointField Foam::extrudedMesh::extrudedPoints
(
    const PrimitivePatch<Face, FaceList, PointField>& extrudePatch,
    const label nLayers,
    const PointExtruder& extruder
)
{
    const pointField& surfacePoints = extrudePatch.localPoints();
    const vectorField& surfaceNormals = extrudePatch.pointNormals();

    pointField ePoints((nLayers + 1)*surfacePoints.size());

    forAll(surfacePoints, i)
    {
        ePoints[i] = surfacePoints[i];
    }

    for (label layer=0; layer<nLayers; layer++)
    {
        label offset = (layer + 1)*surfacePoints.size();

        forAll(surfacePoints, i)
        {
            ePoints[offset + i] =
                extruder
                (
                    surfacePoints[i],
                    surfaceNormals[i],
                    nLayers,
                    layer
                );
        }
    }

    return ePoints;
}


template<class Face, template<class> class FaceList, class PointField>
Foam::faceList Foam::extrudedMesh::extrudedFaces
(
    const PrimitivePatch<Face, FaceList, PointField>& extrudePatch,
    const label nLayers
)
{
    const pointField& surfacePoints = extrudePatch.localPoints();
    const List<face>& surfaceFaces = extrudePatch.localFaces();
    const edgeList& surfaceEdges = extrudePatch.edges();
    const label nInternalEdges = extrudePatch.nInternalEdges();

    label nFaces = 
        (nLayers + 1)*surfaceFaces.size() + nLayers*surfaceEdges.size();

    faceList eFaces(nFaces);

    labelList quad(4);
    label facei = 0;

    // Internal faces
    for (label layer=0; layer<nLayers; layer++)
    {
        label currentLayerOffset = layer*surfacePoints.size();
        label nextLayerOffset = currentLayerOffset + surfacePoints.size();

        // Side faces from layer to layer+1
        for (label i = 0; i < nInternalEdges; i++)
        {
            quad[0] = surfaceEdges[i][1] + currentLayerOffset;
            quad[1] = surfaceEdges[i][0] + currentLayerOffset;
            quad[2] = surfaceEdges[i][0] + nextLayerOffset;
            quad[3] = surfaceEdges[i][1] + nextLayerOffset;

            eFaces[facei++] = face(quad);
        }

        // Faces between layer and layer+1
        if (layer < nLayers-1)
        {
            forAll(surfaceFaces, i)
            {
                eFaces[facei++] =
                    face
                    (
                        surfaceFaces[i].reverseFace()
                      + nextLayerOffset
                    );
            }
        }
    }

    // External side faces
    for (label layer=0; layer<nLayers; layer++)
    {
        label currentLayerOffset = layer*surfacePoints.size();
        label nextLayerOffset = currentLayerOffset + surfacePoints.size();

        // Side faces across layer
        for (label i = nInternalEdges; i < surfaceEdges.size(); i++)
        {
            const edge& e = surfaceEdges[i];
            quad[0] = e[1] + currentLayerOffset;
            quad[1] = e[0] + currentLayerOffset;
            quad[2] = e[0] + nextLayerOffset;
            quad[3] = e[1] + nextLayerOffset;

            label ownerFace = extrudePatch.edgeFaces()[i][0];

            if (!sameOrder(surfaceFaces[ownerFace], e))
            {
                reverse(quad);
            }

            eFaces[facei++] = face(quad);
        }
    }

    // Top faces
    forAll(surfaceFaces, i)
    {
        eFaces[facei++] =  face(surfaceFaces[i]);
    }

    // Bottom faces
    forAll(surfaceFaces, i)
    {
        eFaces[facei++] =
            face
            (
                surfaceFaces[i].reverseFace()
             +  nLayers*surfacePoints.size()
            );
    }

    return eFaces;
}


template<class Face, template<class> class FaceList, class PointField>
Foam::cellList Foam::extrudedMesh::extrudedCells
(
    const PrimitivePatch<Face, FaceList, PointField>& extrudePatch,
    const label nLayers
)
{
    const List<face>& surfaceFaces = extrudePatch.localFaces();
    const edgeList& surfaceEdges = extrudePatch.edges();
    const label nInternalEdges = extrudePatch.nInternalEdges();

    cellList eCells(nLayers*surfaceFaces.size());

    // Size the cells
    forAll(surfaceFaces, i)
    {
        const face& f = surfaceFaces[i];

        for (label layer=0; layer<nLayers; layer++)
        {
            eCells[i+layer*surfaceFaces.size()].setSize(f.size() + 2);
        }
    }

    // Current face count per cell.
    labelList nCellFaces(eCells.size(), 0);


    label facei = 0;

    for (label layer=0; layer<nLayers; layer++)
    {
        // Side faces from layer to layer+1
        for (label i = 0; i < nInternalEdges; i++)
        {
            // Get patch faces using edge
            const labelList& edgeFaces = extrudePatch.edgeFaces()[i];

            // Get cells on this layer
            label cell0 = layer*surfaceFaces.size() + edgeFaces[0];
            label cell1 = layer*surfaceFaces.size() + edgeFaces[1];

            eCells[cell0][nCellFaces[cell0]++] = facei;
            eCells[cell1][nCellFaces[cell1]++] = facei;

            facei++;
        }

        // Faces between layer and layer+1
        if (layer < nLayers-1)
        {
            forAll(surfaceFaces, i)
            {
                label cell0 = layer*surfaceFaces.size() + i;
                label cell1 = (layer+1)*surfaceFaces.size() + i;

                eCells[cell0][nCellFaces[cell0]++] = facei;
                eCells[cell1][nCellFaces[cell1]++] = facei;

                facei++;
            }
        }
    }

    // External side faces
    for (label layer=0; layer<nLayers; layer++)
    {
        // Side faces across layer
        for (label i = nInternalEdges; i < surfaceEdges.size(); i++)
        {
            // Get patch faces using edge
            const labelList& edgeFaces = extrudePatch.edgeFaces()[i];

            // Get cells on this layer
            label cell0 = layer*surfaceFaces.size() + edgeFaces[0];

            eCells[cell0][nCellFaces[cell0]++] = facei;

            facei++;
        }
    }

    // Top faces
    forAll(surfaceFaces, i)
    {
        eCells[i][nCellFaces[i]++] = facei;

        facei++;
    }

    // Bottom faces
    forAll(surfaceFaces, i)
    {
        label cell0 = (nLayers-1)*surfaceFaces.size() + i;

        eCells[cell0][nCellFaces[cell0]++] = facei;

        facei++;
    }

    return eCells;
}


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

//- Construct from the primitivePatch to extrude
template
<
    class Face,
    template<class> class FaceList,
    class PointField,
    class PointExtruder
>
Foam::extrudedMesh::extrudedMesh
(
    const IOobject& io,
    const PrimitivePatch<Face, FaceList, PointField>& extrudePatch,
    const label nLayers,
    const PointExtruder& extruder
)
:
    polyMesh
    (
        io,
        extrudedPoints(extrudePatch, nLayers, extruder),
        extrudedFaces(extrudePatch, nLayers),
        extrudedCells(extrudePatch, nLayers)
    )
{
    List<polyPatch*> patches(3);

    label facei = nInternalFaces();

    label sz = nLayers*(extrudePatch.nEdges() - extrudePatch.nInternalEdges());

    patches[0] = new wallPolyPatch
    (
        "sides",
        sz,
        facei,
        0,
        boundaryMesh()
    );

    facei += sz;

    patches[1] = new polyPatch
    (
        "originalPatch",
        extrudePatch.size(),
        facei,
        1,
        boundaryMesh()
    );

    facei += extrudePatch.size();

    patches[2] = new polyPatch
    (
        "otherSide",
        extrudePatch.size(),
        facei,
        2,
        boundaryMesh()
    );

    addPatches(patches);
}


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


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