/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
     \\/     M anipulation  | Copyright (C) 2017 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 "ListOps.H"
#include "point.H"
#include "Field.H"

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

template<class PointList>
Foam::label Foam::mergePoints
(
    const PointList& points,
    const scalar mergeTol,
    const bool verbose,
    labelList& pointMap,
    typename PointList::const_reference origin
)
{
    typedef typename PointList::value_type point_type;

    const label nPoints = points.size();

    // Create an old to new point mapping array
    pointMap.setSize(nPoints);
    pointMap = -1;

    if (!nPoints)
    {
        return 0;
    }

    point_type compareOrigin = origin;
    if (origin == point_type::max)
    {
        // Use average of input points to define a comparison origin.
        // Same as sum(points)/nPoints, but handles different list types
        compareOrigin = points[0];
        for (label pointi=1; pointi < nPoints; ++pointi)
        {
            compareOrigin += points[pointi];
        }
        compareOrigin /= nPoints;
    }

    // We're comparing distance squared to origin first.
    // Say if starting from two close points:
    //     x, y, z
    //     x+mergeTol, y+mergeTol, z+mergeTol
    // Then the magSqr of both will be
    //     x^2+y^2+z^2
    //     x^2+y^2+z^2 + 2*mergeTol*(x+z+y) + mergeTol^2*...
    // so the difference will be 2*mergeTol*(x+y+z)

    const scalar mergeTolSqr = Foam::sqr(mergeTol);

    // Sort points by magSqr
    List<scalar> magSqrDist(nPoints);
    forAll(points, pointi)
    {
        magSqrDist[pointi] = magSqr(points[pointi] - compareOrigin);
    }
    labelList order;
    Foam::sortedOrder(magSqrDist, order);


    Field<scalar> sortedTol(nPoints);
    forAll(order, sortI)
    {
        const point_type& pt = points[order[sortI]];

        // Use scalar precision
        sortedTol[sortI] =
            2*mergeTol*
            (
                mag(scalar(pt.x() - compareOrigin.x()))
              + mag(scalar(pt.y() - compareOrigin.y()))
              + mag(scalar(pt.z() - compareOrigin.z()))
            );
    }

    label newPointi = 0;

    // Handle 0th point separately (is always unique)
    label pointi = order[0];
    pointMap[pointi] = newPointi++;

    for (label sortI = 1; sortI < order.size(); ++sortI)
    {
        // Get original point index
        const label pointi = order[sortI];
        const scalar mag2 = magSqrDist[order[sortI]];

        // Convert to scalar precision
        // NOTE: not yet using point_type template parameter
        const point pt
        (
            scalar(points[pointi].x()),
            scalar(points[pointi].y()),
            scalar(points[pointi].z())
        );


        // Compare to previous points to find equal one.
        label equalPointi = -1;

        for
        (
            label prevSortI = sortI - 1;
            prevSortI >= 0
         && (mag(magSqrDist[order[prevSortI]] - mag2) <= sortedTol[sortI]);
            --prevSortI
        )
        {
            const label prevPointi = order[prevSortI];

            // Convert to scalar precision
            // NOTE: not yet using point_type template parameter
            const point prevPt
            (
                scalar(points[prevPointi].x()),
                scalar(points[prevPointi].y()),
                scalar(points[prevPointi].z())
            );

            if (magSqr(pt - prevPt) <= mergeTolSqr)
            {
                // Found match.
                equalPointi = prevPointi;

                break;
            }
        }


        if (equalPointi != -1)
        {
            // Same coordinate as equalPointi. Map to same new point.
            pointMap[pointi] = pointMap[equalPointi];

            if (verbose)
            {
                Pout<< "Foam::mergePoints : Merging points "
                    << pointi << " and " << equalPointi
                    << " with coordinates:" << points[pointi]
                    << " and " << points[equalPointi]
                    << endl;
            }
        }
        else
        {
            // Differs. Store new point.
            pointMap[pointi] = newPointi++;
        }
    }

    return newPointi;
}


template<class PointList>
bool Foam::mergePoints
(
    const PointList& points,
    const scalar mergeTol,
    const bool verbose,
    labelList& pointMap,
    List<typename PointList::value_type>& newPoints,
    typename PointList::const_reference origin
)
{
    const label nUnique = mergePoints
    (
        points,
        mergeTol,
        verbose,
        pointMap,
        origin
    );

    newPoints.setSize(nUnique);
    forAll(pointMap, pointi)
    {
        newPoints[pointMap[pointi]] = points[pointi];
    }

    return (nUnique != points.size());
}


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