/*---------------------------------------------------------------------------*\
  ==  == ====== ====   ====    |
                    \\     ||  | Multiphase Code Repository by HZDR
  ======   //   ||  || ===//   | Website: https://doi.org/10.14278/rodare.767
  ||  ||  //    ||  // || \\   | License: GPL-3.0-or-later
  ==  == ====== ====   ==  ==  |
-------------------------------------------------------------------------------
License
    This file is part of the Multiphase Code Repository by HZDR.

    Copyright (C) 2025 by Helmholtz-Zentrum Dresden-Rossendorf e.V. (HZDR),
    Website: https://hzdr.de

    Multiphase Code Repository by HZDR is based on the free software for
    computational fluid dynamics (CFD) from the OpenFOAM Foundation.
    Copyright (C) 2025 by OpenFOAM Foundation, Website: https://openfoam.org

    If you are interested in which files are original OpenFOAM Foundation
    files, which OpenFOAM Foundation files were modified, and which files were
    newly created, see FILES.md.

    Multiphase Code Repository by HZDR 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.

    Multiphase Code Repository by HZDR 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 the Multiphase Code Repository by HZDR. If not, see
    <http://www.gnu.org/licenses/>.

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

#include "clusterTrackingModel.H"
#include "phaseSystem.H"
#include "fvcGrad.H"
#include "surfaceInterpolate.H"

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

namespace Foam
{
    defineTypeNameAndDebug(clusterTrackingModel, 0);
    defineRunTimeSelectionTable(clusterTrackingModel, dictionary);
}


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

Foam::clusterTrackingModel::clusterTrackingModel
(
    const dictionary& dict,
    const phaseInterface& interface,
    const bool registerObject
)
:
    regIOobject
    (
        IOobject
        (
            IOobject::groupName(typeName, interface.name()),
            interface.mesh().time().name(),
            interface.mesh(),
            IOobject::NO_READ,
            IOobject::NO_WRITE,
            registerObject
        )
    ),
    mesh_(interface.mesh()),
    twoDimensional_(false),
    thickness_(1),
    compressedPhase_
    (
        mesh_.lookupObject<phaseModel>
        (
            IOobject::groupName("alpha", dict.lookup("compressedPhase"))
        )
    ),
    emergingPhase_
    (
        mesh_.lookupObject<phaseModel>
        (
            IOobject::groupName("alpha", dict.lookup("emergingPhase"))
        )
    ),
    cellsPerDiameterThreshold_
    (
        dict.lookupOrDefault<scalar>
        (
            "cellsPerDiameterThreshold",
            dimless,
            8
        )
    ),
    maxIter_
    (
        dict.lookupOrDefault<label>
        (
            "maxIter",
            100
        )
    ),
    particleVolumeTol_
    (
        dict.lookupOrDefault<label>
        (
            "particleVolumeTol",
            10
        )
    ),
    diameter_(0),
    nClusters_(0),
    nActiveClusters_(0),
    clusterCentre_(0),
    clusterRadius_(0),
    initialRadius_(0),
    centreProximityFieldName_(dict.lookup("proximityFieldName")),
    centreProximity_
    (
        IOobject
        (
            centreProximityFieldName_,
            mesh_.time().name(),
            mesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh_,
        dimensionedScalar("zero", dimless, Zero)
    ),
    colour_
    (
        IOobject
        (
            IOobject::groupName("clusterColour", interface.name()),
            mesh_.time().name(),
            mesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh_,
        dimensionedScalar("zero", dimless, Zero)
    ),
    overResolutionIndicator_
    (
        IOobject
        (
            IOobject::groupName
            (
                "overResolutionIndicator", compressedPhase_.name()
            ),
            mesh_.time().name(),
            mesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh_,
        dimensionedScalar("zero", dimless, Zero)
    )
{

    // Define initial cluster radius
    initialRadius_ = mag(mesh_.bounds().span());

    // Check dimension of domain
    if (mesh_.nGeometricD() == 2)
    {
        Info<< "Case is 2D" << endl;

        twoDimensional_ = true;

        vector plane(0,0,0);
        if (mesh_.geometricD().x() == -1)
        {
            plane = vector(1,0,0);
            Info<< "Empty direction is x" << endl;
        }
        else if (mesh_.geometricD().y() == -1)
        {
            plane = vector(0,1,0);
            Info<< "Third direction is y" << endl;
        }
        else if (mesh_.geometricD().z() == -1)
        {
            plane = vector(0,0,1);
            Info<< "Third direction is z" << endl;
        }

        const cell& faces = mesh_.cells()[0];
        const labelList points = faces.labels(mesh_.faces());
        pointField pLocal(points.size(), vector::zero);
        forAll(points, i)
        {
            const label& pointI = points[i];
            pLocal[i] = mesh_.points()[pointI];
        }

        thickness_ = Foam::max(pLocal & plane) - Foam::min(pLocal & plane);
    }
    else if (mesh_.nGeometricD() == 3)
    {
        twoDimensional_ = false;
    }
    else
    {
        FatalErrorInFunction
            << "Not implemented for 1D."
            << exit(FatalError);
    }
}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

Foam::clusterTrackingModel::~clusterTrackingModel()
{}


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

Foam::phaseInterfaceKey Foam::clusterTrackingModel::key() const
{
    return phaseInterfaceKey(compressedPhase_, emergingPhase_);
}


Foam::scalar Foam::clusterTrackingModel::diameter() const
{
    return diameter_;
}


const Foam::phaseModel& Foam::clusterTrackingModel::compressedPhase() const
{
    return compressedPhase_;
}


Foam::label Foam::clusterTrackingModel::updateClusters()
{
    label nInvalidClusters(0);

    // Find particle diameter
    diameter_ = gMax(volScalarField(compressedPhase_.d()));

    // Target particle volume
    const scalar particleVolume
    (
        twoDimensional_
      ? 0.25*constant::mathematical::pi*sqr(diameter_)*thickness_
      : constant::mathematical::pi/6.0*pow3(diameter_)
    );

    const scalarField cellWidth
    (
        twoDimensional_
      ? sqrt(mesh_.V()/thickness_)
      : cbrt(mesh_.V())
    );

    overResolutionIndicator_.primitiveFieldRef() =
        pos0(diameter_*1.001 - cellsPerDiameterThreshold_*cellWidth);

    overResolutionIndicator_.correctBoundaryConditions();

    // No over-resolved mesh regions
    if (!gSum(overResolutionIndicator_))
    {
        this->discard();
        return 0;
    }

    const scalarField alphaHybrid
    (
        overResolutionIndicator_*compressedPhase_ + emergingPhase_
    );

    // Calculate number of resolvable structures
    const label nClustersOld(nClusters_);

    nClusters_ = floor(gSum(alphaHybrid*mesh_.V())/particleVolume);

    const bool nClustersChanged(nClusters_ != nClustersOld);

    if (debug)
    {
        Info<< "Cluster tracking: potential structures = "
            << gSum(alphaHybrid*mesh_.V())/particleVolume
            << endl;
    }

    // Detect clusters based on recent cluster properties
    if (nClusters_)
    {
        // Re-init clusters if number of clusters changes
        if (nClustersChanged)
        {
            // Reset cluster properties
            clusterCentre_.clear();
            clusterRadius_.clear();

            const vector deltaLineEquidistant
            (
                1/float(nClusters_ + 1)*mesh_.bounds().span()
            );

            for (int i = 0; i < nClusters_; i++)
            {
                // Initialise cluster centre positions along diagonal,
                // add small offsets to break symmetry
                vector pointPos
                (
                    mesh_.bounds().min() + (i + 1)*deltaLineEquidistant
                );
                const scalar offset(0.1*mag(deltaLineEquidistant));
                const label dir(i%3); // cycle offset direction
                pointPos[dir] += offset; // add offset
                // Make sure to be within domain
                pointPos[dir] = min(pointPos[dir], mesh_.bounds().max()[dir]);
                pointPos[dir] = max(pointPos[dir], mesh_.bounds().min()[dir]);

                clusterCentre_.insert(i, pointPos);
                clusterRadius_.insert(i, initialRadius_);
            }

            if (debug)
            {
                forAll(clusterCentre_, i)
                {
                    Info<< "Cluster initialised " << i + 1
                        << " " << clusterCentre_[i]
                        << " " << clusterRadius_[i] << endl;
                }
            }
        }

        // Quantities reused during and after iteration
        List<label> cellsInCluster(nClusters_);
        List<scalar> materialInCluster(nClusters_);

        const scalarField hybridVolume(alphaHybrid*mesh_.V());
        const vectorField hybridVolumeC(hybridVolume*mesh_.C());

        // Adjust clusters iteratively
        for (int k = 1; k <= maxIter_; k++)
        {
            if (debug)
            {
                Info << "Cluster iteration " << k << endl;
            }

            // Find closest cluster
            centreProximity_ = 0;
            colour_ = 0;

            forAll(clusterCentre_, i)
            {
                // Measure proximity of cells to this cluster
                const scalarField proximity
                (
                    posPart
                    (
                        1
                      - magSqr
                        (
                            mesh_.C()*dimensionedScalar(inv(dimLength), 1)
                          - clusterCentre_[i]
                        )
                       /sqr(clusterRadius_[i])
                    )
                );

                // Check if this cluster is closer than previous clusters
                const scalarField delta(proximity - centreProximity_);

                // Overwrite proximity if this cluster is closer
                centreProximity_.primitiveFieldRef() += posPart(delta);

                // Overwrite colour if this cluster is closer
                colour_.primitiveFieldRef() =
                    pos(delta)*(i+1)
                  + neg0(delta)*colour_.primitiveFieldRef();
            }

            // Update cluster positions and radii
            bool clustersAdjusted(false);

            forAll(clusterCentre_, i)
            {
                // Mark cells in cluster (cluster i  has colour i+1)
                const scalarField mask(neg0(mag(i+1 - colour_)));

                cellsInCluster[i] = gSum(mask);

                materialInCluster[i] = gSum(mask*hybridVolume);

                if (cellsInCluster[i])
                {
                    // Some material in cluster
                    if (materialInCluster[i] > 1e-9*particleVolume)
                    {
                        const scalar volumeError
                        (
                            materialInCluster[i] - particleVolume
                        );

                        if
                        (
                            fabs(volumeError) > 0.01*particleVolume
                         || clusterRadius_[i] > 5*diameter_
                        )
                        {
                            clustersAdjusted = true;

                            // Update cluster position and cluster size
                            clusterCentre_[i] =
                                gSum(mask*hybridVolumeC)/materialInCluster[i];

                            clusterRadius_[i] *= (1 - 0.01*sign(volumeError));
                        }
                    }
                    else // No material in cluster
                    {
                        // Simple cluster position update
                        clusterCentre_[i] =
                            gSum(mask*mesh_.V()*mesh_.C())
                           /gSum(mask*mesh_.V());
                    }
                }

                if (debug)
                {
                    if (cellsInCluster[i])
                    {
                        Info<< " Cluster " << i+1
                            << " Vol actual/target = "
                            << materialInCluster[i]/particleVolume
                            << ", continuous fraction = "
                            << gSum(mask*emergingPhase_*mesh_.V())
                              /particleVolume
                            << endl;
                    }
                    else
                    {
                        Info<< " Cluster " << i+1 << " is empty" << endl;
                    }
                }
            }

            if (!clustersAdjusted)
            {
                break;
            }
        }

        // Handle invalid clusters
        forAll(clusterCentre_, i)
        {
            if (debug)
            {
                Info<< "Check cluster " << i+1 << endl;
            }

            bool valid(true);

            // Mark cells in cluster (cluster i  has colour i+1)
            const scalarField mask(neg0(mag(i+1 - colour_)));

            // Not enough material contained (no cells or no material)
            if (materialInCluster[i] < 0.99*particleVolume)
            {
                valid = false;
                if (debug)
                {
                    Info<< "> Passive (not enough material)" << endl;
                }
            }
            else // Enough material contained
            {
                const scalar continuousVolumeInCluster
                (
                    gSum(mask*emergingPhase_*mesh_.V())
                );

                const scalar clusterVolume(gSum(mask*mesh_.V()));

                const scalar averageCellVolumeInCluster
                (
                    clusterVolume/cellsInCluster[i]
                );

                // Cluster saturated with continuous material
                if
                (
                    fabs(continuousVolumeInCluster - particleVolume)
                    <= particleVolumeTol_*averageCellVolumeInCluster
                )
                {
                    valid = false;
                    if (debug)
                    {
                        Info<< "> Passive (saturated)" << endl;
                    }
                }
                // Cluster too huge to describe a particle (compactness)
                else if (materialInCluster[i] < 0.01*clusterVolume)
                {
                    valid = false;
                    if (debug)
                    {
                        Info<< "> Passive (residual cluster)" << endl;
                    }
                }
                else // Cluster valid
                {
                    continue;
                }
            }

            // Count invalid clusters and set corresponding output to zero
            if (!valid)
            {
                nInvalidClusters++;
                centreProximity_.primitiveFieldRef() *= (1 - mask);
            }
        }

        centreProximity_.correctBoundaryConditions();
    }

    // Remaining clusters
    nActiveClusters_ = nClusters_ - nInvalidClusters;

    Info<< "Cluster tracking: " << compressedPhase_.name()
        << " structures active/all = " << nActiveClusters_ << "/"
        << nClusters_  << endl;

    if (debug && mesh_.time().writeTime())
    {
        overResolutionIndicator_.write();
        colour_.write();
        centreProximity_.write();
    }

    return nActiveClusters_;
}


bool Foam::clusterTrackingModel::isActive()
{
    return (nActiveClusters_ > 0);
}


void Foam::clusterTrackingModel::discard()
{
    if (nActiveClusters_)
    {
        centreProximity_ = 0;
        centreProximity_.correctBoundaryConditions();

        colour_ = 0;
        colour_.correctBoundaryConditions();
    }

    nActiveClusters_ = 0;
}


bool Foam::clusterTrackingModel::writeData(Ostream& os) const
{
    return os.good();
}


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