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

Application
    setBubble

Description
    Calculates the volume fraction field for a two-dimensional circular bubble.
    For the cells, which are cut by the interface, the volume fraction is
    obtained analytically and has a value between zero and one.

Usage
    \b setBubble

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

#include "argList.H"
#include "timeSelector.H"
#include "fvMesh.H"
#include "volFields.H"
#include "systemDict.H"

using namespace Foam;

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


//- Check, whether a point is inside the circle or not, in projection into given
//- plane pl
bool checkInside(const vector& p1, const scalar& r, const vector& pl)
{
    const scalar epsilon = (p1 & pl)/magSqr(pl);

    if (mag(p1 - epsilon*pl) <= r)
    {
        return true;
    }
    else
    {
        return false;
    }
}


//- Check whether two points, p1 and p2, are identical, when projected into
//- given plane pl
bool equivalent
(
    const vector& p1,
    const vector& p2,
    const vector& pl,
    const scalar& tolerance = small
)
{
    if (mag((p2 - p1) ^ pl) < tolerance)
    {
        return true;
    }
    else
    {
        return false;
    }
}


//- Check whether two points, p1 and p2, are located in the same given plane, pl
bool samePlane
(
    const vector& p1,
    const vector& p2,
    const vector& pl,
    const scalar& tolerance = small
)
{
    if (mag((p2 - p1) & pl) < tolerance)
    {
        return true;
    }
    else
    {
        return false;
    }
}


//- Return distance between points, p1 and p2, in direction normal to plane, pl
scalar planeDist
(
    const vector& p1,
    const vector& p2,
    const vector& pl
)
{
    return mag((p2 - p1) & pl)/mag(pl);
}


//- Return 3D radius for given point p1, in plane, pl, and 2D radius, r
scalar radius3D
(
    const vector& p1,
    const vector& pl,
    const scalar& r
)
{
    const scalar delta = planeDist(p1, vector(0, 0, 0), pl);

    return Foam::sqrt(sqr(r) + sqr(delta));
}


//- Calculate angle between two vectors, p1 and p2 by inner product
Foam::scalar vecAngle
(
    const vector& p1,
    const vector& p2,
    const scalar& tolerance = small
)
{
    const scalar cosPhi = (p1 & p2)/(max(mag(p1)*mag(p2), tolerance));

    // Enforce bounding between -1 and 1
    return Foam::acos(min(max(cosPhi, scalar(-1)), scalar(1)));
}


//- Get point of cut between edge and interface from two edge points
//- Use condition: p1 + epsilon*(p2 - p1) = mag(r)
//- and quadratic formula solution for epsilon
vector cutPoint
(
    const vector& p1,
    const vector& p2,
    const scalar& r,
    const scalar& tolerance = small
)
{
    // Solve a*epsilon^2 + b*epsilon + c = 0
    const double a = magSqr(p1) - 2.*(p1 & p2) + magSqr(p2);
    const double b = 2.*((p1 & p2) - magSqr(p1));
    const double c = magSqr(p1) - sqr(r);

    const double epsilon1 = ( -b + Foam::sqrt(sqr(b) - 4.*a*c))/(2.*a);
    const double epsilon2 = ( -b - Foam::sqrt(sqr(b) - 4.*a*c))/(2.*a);

    // Decide, which solution to choose
    if
    (
        (epsilon1 >= - tolerance && epsilon1 <= 1. + tolerance)
     &&!(epsilon2 >= - tolerance && epsilon2 <= 1. + tolerance)
    )
    {
        return p1 + min(max(epsilon1, Zero), scalar(1.0))*(p2 - p1);
    }
    else if
    (
       !(epsilon1 >= - tolerance && epsilon1 <= 1. + tolerance)
     && (epsilon2 >= - tolerance && epsilon2 <= 1. + tolerance)
    )
    {
        return p1 + min(max(epsilon2, Zero), scalar(1.0))*(p2 - p1);
    }
    else if
    (
        (epsilon1 >= - tolerance && epsilon1 <= 1. + tolerance)
     && (epsilon2 >= - tolerance && epsilon2 <= 1. + tolerance)
    )
    {
        if (mag(epsilon1 - epsilon2) < tolerance)
        {
            return p1 + min(max(epsilon1, Zero), scalar(1.0))*(p2 - p1);
        }
        else
        {
            WarningInFunction
                << "Two interface cuts found between points "
                << p1 << " and " << p2 << "."
                << nl << "    epsilon1 = " << epsilon1
                << nl << "    epsilon2 = " << epsilon2
                << nl << "    This case is not accounted for yet."
                << exit(FatalError);

            return vector(0, 0, 0);
        }
    }
    else
    {
        FatalErrorInFunction
            << "No interface cut found between points "
            << p1 << " and " << p2 << "." << nl
            << "Epsilon1 = " << epsilon1 << " and" << nl
            << "Epsilon2 = " << epsilon2 << "."
            << exit(FatalError);

        return vector(0, 0, 0);
    }
}


//- Return the surface area of a general triangle
Foam::scalar areaTriangle
(
    const vector& p1,
    const vector& p2,
    const vector& p3
)
{
    const scalar a = mag(p2 - p1);
    const scalar b = mag(p3 - p2);
    const scalar c = mag(p1 - p3);

    const scalar s = 0.5*(a + b + c);

    return Foam::sqrt(s*(s - a)*(s - b)*(s - c));
}


//- Return the surface area of a general quadrangle
Foam::scalar areaQuadrangle
(
    const vector& p1,
    const vector& p2,
    const vector& p3,
    const vector& p4
)
{
    const vector& diag1 = p3 - p1;
    const vector& diag2 = p4 - p2;

    const scalar angle = vecAngle(diag1, diag2);

    return 0.5*mag(diag1)*mag(diag2)*Foam::sin(angle);
}


//- Return the surface area of a segment of a circle
Foam::scalar areaSegment(const vector& p1, const vector& p2)
{
    const scalar r = mag(p1);
    const scalar alpha = vecAngle(p1, p2);

    return 0.5*sqr(r)*(alpha - Foam::sin(alpha));
}


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

int main(int argc, char *argv[])
{
    argList::noParallel();      // utility runs only serial

    #include "addDictOption.H"
    #include "addRegionOption.H"
    timeSelector::addOptions();
    #include "setRootCase.H"
    #include "createTime.H"

    const instantList Times = timeSelector::select0(runTime, args);

    #include "createRegionMeshNoChangers.H"
    #include "createFields.H"

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

    const cellList& cells = mesh.cells();
    const scalarField& cellVols = mesh.V();
    const pointField points(mesh.points() - centre);
    const List<edge> edges = mesh.edges();
    const volVectorField centres
    (
        mesh.C() - dimensionedVector("centre", dimLength, centre)
    );

    //----------------------------------CELLS----------------------------------
    forAll(cells, i)
    {
        const labelList& cellPoints = mesh.cellPoints()[i];
        const labelList& cellEdges = mesh.cellEdges()[i];

        //-------------------------------POINTS--------------------------------
        List<point> inside;
        List<point> inside2D;

        // Select relevant points inside circle in selected plane
        forAll(cellPoints, j)
        {
            point vertex = points[cellPoints[j]];

            // Check, if inside circle
            if (checkInside(vertex, radius, plane))
            {
                inside.append(vertex);

                if (inside2D.size() == 0)
                {
                    // Selected plane defined by first point
                    inside2D.append(vertex);
                }
                else if (samePlane(vertex, inside2D[0], plane))
                {
                    inside2D.append(vertex);
                }
            }
        }

        //---------------------------------FULL--------------------------------
        if (inside.size() == cellPoints.size())
        {
            // Cell completely inside circle
            alpha[i] = 1.0;
        }
        //--------------------------------EMPTY--------------------------------
        else if (inside.size() == 0)
        {
            // Cell completely outside circle
            alpha[i] = 0.0;
        }
        //---------------------------------HALF--------------------------------
        else if ((inside.size() > 0) && (inside.size() < cellPoints.size()))
        {
            // Cell cut by circle
            scalar gasVol = 0.;
            scalar depth = 0.;

            scalar numCutPoints = 0;

            point cutPoint1;
            point cutPoint2;

            point nextToCutPoint1;
            point nextToCutPoint2;

            //------------------------------EDGE-------------------------------
            // Get cut points and neighbouring corner points inside circle
            forAll(cellEdges, j)
            {
                const edge edge_ = edges[cellEdges[j]];

                const point& vertex1 = points[edge_.start()];
                const point& vertex2 = points[edge_.end()];

                point insidePoint;

                // Reject if points are not in same plane
                if (!samePlane(vertex1, vertex2, plane))
                {
                    if (equivalent(vertex1, vertex2, plane))
                    {
                        // Save third dimension for volume calculation
                        depth = mag(vertex2 - vertex1);
                    }

                    continue;
                }

                // Select single valid point inside circle
                bool isInside2D1 = false;
                bool isInside2D2 = false;
                for (int k = 0; k < inside2D.size(); k++)
                {
                    if (vertex1 == inside2D[k])
                    {
                        isInside2D1 = true;
                    }

                    if (vertex2 == inside2D[k])
                    {
                        isInside2D2 = true;
                    }
                }

                // Only point 1 is inside circle
                if (isInside2D1 && !isInside2D2)
                {
                    insidePoint = vertex1;
                }
                // Only point 2 is inside circle
                else if (!isInside2D1 && isInside2D2)
                {
                    insidePoint = vertex2;
                }
                // Both points or none of them are inside circle
                else
                {
                    continue;
                }

                // Get interface location at edge
                const scalar radiusStar = radius3D(vertex1, plane, radius);
                point edgeCutPoint =
                    cutPoint(vertex1, vertex2, radiusStar, 1.e-8);

                // Save 2 interface cut points and 2 neighbouring points in cell
                if (numCutPoints == 0)
                {
                    cutPoint1 = edgeCutPoint;
                    nextToCutPoint1 = insidePoint;
                }
                else if (numCutPoints == 1)
                {
                    cutPoint2 = edgeCutPoint;
                    nextToCutPoint2 = insidePoint;
                }
                else
                {
                    FatalErrorInFunction
                        << "Found more than " << numCutPoints
                        << " interface cut points in cell " << i << "."
                        << exit(FatalError);
                }

                numCutPoints += 1;
            }

            //-----------------------------VOLUME------------------------------
            if (inside2D.size() == 1)
            {
                // Interface exactly at corner point
                if (cutPoint1 == cutPoint2)
                {
                    gasVol = 0.;
                }
                // One corner inside circle
                else
                {
                    gasVol =
                        (
                            areaSegment(cutPoint1, cutPoint2)
                          + areaTriangle(nextToCutPoint1, cutPoint1, cutPoint2)
                        )
                       *depth;
                }
            }
            // Two corners inside circle
            else if (inside2D.size() == 2)
            {
                gasVol =
                    (
                        areaSegment(cutPoint1, cutPoint2)
                      + areaQuadrangle
                        (
                            nextToCutPoint1,
                            cutPoint1,
                            cutPoint2,
                            nextToCutPoint2
                        )
                    )
                   *depth;
            }
            // Three corners inside circle
            else if (inside2D.size() == 3)
            {
                point centralPoint;

                // Select point, which is not adjacent to cutPoints
                for (int k = 0; k < inside2D.size(); k++)
                {
                    if
                    (
                        (inside2D[k] != nextToCutPoint1)
                     && (inside2D[k] != nextToCutPoint2)
                    )
                    {
                        centralPoint = inside2D[k];

                        break;
                    }
                }

                gasVol =
                    (
                        areaSegment(cutPoint1, cutPoint2)
                      + areaQuadrangle
                        (
                            nextToCutPoint1,
                            cutPoint1,
                            cutPoint2,
                            nextToCutPoint2
                        )
                      + areaTriangle
                        (
                            nextToCutPoint1,
                            nextToCutPoint2,
                            centralPoint
                        )
                    )
                   *depth;
            }
            else
            {
                FatalErrorInFunction
                    << "Found " << inside2D.size()
                    << " points inside circle in cell " << i << "."
                    << exit(FatalError);
            }

            // Volume fraction
            alpha[i] = gasVol/cellVols[i];
        }
        else
        {
            FatalErrorInFunction
                << "Found " << inside.size()
                << " points inside the circle for cell " << i
                << ", but there are only " << cellPoints.size()
                << " cellPoints."
                << exit(FatalError);
        }
    }

    // Invert volume fraction field
    if (invert)
    {
        alpha = 1 - alpha;
    }

    Info<< "Writing field " << fieldName << nl << endl;
    alpha.write();

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

    return 0;
}


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