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

Class
    Foam::blendingMethods::psiSurf

Description
    This blending method utilises the AIAD algorithm based on the phase fraction
    and the GENTOP algorithm based on the gradient of the phase fraction.

    For details refer to
    \verbatim
        Gauss, F., & Porombka, P. (2015).
        Memo on the unified blending functions in AIAD and GENTOP.
        Helmholtz-Zentrum Dresden-Rossendorf,
        23. Sept. 2015, pp. 1-5.
    \endverbatim

    Blending factors
    \f[
        f_{1} = \phi_{1,surf}\,f_{1 \rightarrow 2}
    \f]
    \f[
        f_{2} = \phi_{2,surf}\,f_{2 \rightarrow 1}
    \f]

    Unified blending function
    \f[
        \Psi_{surf}=f_{1} - f_{2}
    \f]

    Note, that the blending factors could be stored during runtime by enabling
    the debug switch 'blendingMethod'. Calls to functions f1() and f2() may be
    enforced with a debug switch for a specific model type,
    e.g. 'BlendedInterfacialModel<dragModel> '.

    Blending based on phase fraction

    \f[
        f_{1 \rightarrow 2} = 1-3\,\tilde{\alpha}_1^2+2\,\tilde{\alpha}_1^3
    \f]
    \f[
        f_{2 \rightarrow 1} = 1-3\,\tilde{\alpha}_2^2+2\,\tilde{\alpha}_2^3
    \f]
    \f[
        \tilde{\alpha}_1 =
            min\left[max\left[\frac{1}{2\delta_{\alpha}}\left(
            \frac{\alpha_1}{\alpha_1+\alpha_2}-\alpha_{1,crit}+
            \delta_{1}\right),0\right],1\right]
    \f]
    \f[
        \tilde{\alpha}_2 =
            min\left[max\left[\frac{1}{2\delta_{\alpha}}\left(
            \frac{\alpha_2}{\alpha_1+\alpha_2}-\alpha_{2,crit}+
            \delta_{2}\right),0\right],1\right]
    \f]

    Blending based on gradient of phase fraction

    \f[
        \phi_{1,surf} =
            1-3\,\widetilde{\nabla\alpha_1}^2+2\,\widetilde{\nabla\alpha_1}^3
    \f]
    \f[
        \phi_{2,surf} =
            1-3\,\widetilde{\nabla\alpha_2}^2+2\,\widetilde{\nabla\alpha_2}^3
    \f]
    \f[
        \widetilde{\nabla\alpha_1} =
            min\left[max\left[\frac{1}{2\delta_{\nabla}}\left(
            \left|\nabla\alpha_1\right|-\nabla\alpha_{crit}+
            \delta_{\nabla}\right),0\right],1\right]
    \f]
    \f[
        \widetilde{\nabla\alpha_2} =
            min\left[max\left[\frac{1}{2\delta_{\nabla}}\left(
            \left|\nabla\alpha_2\right|-\nabla\alpha_{crit}+
            \delta_{\nabla}\right),0\right],1\right]
    \f]
    \f[
        \nabla\alpha_{crit} = \frac{1}{n\sqrt[3]{V}}
    \f]
    \f[
        \delta_{\nabla} = f\cdot\nabla\alpha_{crit}
    \f]

Usage
    \table
        Name | Description | Default
        critAlpha.air | centre of transition region for air | 0.3
        critAlpha.water | centre of transition region for water | 0.3
        deltaAlpha.air | half-width of transition region for air | 0.05
        deltaAlpha.water | half-width of transition region for water | 0.05
        n | number of grid cells covering the transition region | 5
        f | factor to determine half-width of transition region | 0.1
    \endtable

SourceFiles
    psiSurf.C

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

#ifndef psiSurf_H
#define psiSurf_H

#include "blendingMethod.H"

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

namespace Foam
{
namespace blendingMethods
{

/*---------------------------------------------------------------------------*\
                           Class psiSurf Declaration
\*---------------------------------------------------------------------------*/

class psiSurf
:
    public blendingMethod
{
    // Private Data

        //- Centre of transition region (phase fraction)
        HashTable<dimensionedScalar, word, word::hash>
            critAlpha_;

        //- Half-with of transition region (phase fraction)
        HashTable<dimensionedScalar, word, word::hash>
            deltaAlpha_;

        //- Number of grid cells covering the transition region (phase gradient)
        const label nCells_;

        //- Factor to determine half-width of transition region (phase gradient)
        const scalar fDelta_;


    // Private Member Functions

        inline tmp<volScalarField> limit01
        (
            const volScalarField& field
        ) const
        {
            return min(max(field, Zero), scalar(1.0));
        }

        inline tmp<volScalarField> blendingFunction
        (
            const volScalarField& field
        ) const
        {
            //return 0.5*(1.0+cos(Foam::constant::mathematical::pi*field));
            return 1.0 - 3.0*sqr(field) + 2.0*pow3(field);
        }

        tmp<volScalarField> f1in2
        (
            const phaseModel& phase1,
            const phaseModel& phase2
        ) const;

        Foam::tmp<volScalarField> phiSurf
        (
            const phaseModel& phase
        ) const;


protected:

    // Protected Member Functions

        //- Evaluate the blending function
        virtual tmp<volScalarField> fContinuous
        (
            const UPtrList<const volScalarField>& alphas,
            const label phaseSet,
            const label systemSet
        ) const;

        //- Evaluate the blending function. Filters out phases that cannot
        //  become continuous from the sets, then calls fContinuous
        virtual tmp<volScalarField> f
        (
            const UPtrList<const volScalarField>& alphas,
            const label phaseSet,
            const label systemSet
        ) const;


public:

    //- Runtime type information
    TypeNameNoDebug("psiSurf");


    // Constructors

        //- Construct from a dictionary and two phases
        psiSurf
        (
            const dictionary& dict,
            const phaseInterface& interface
        );


    //- Destructor
    ~psiSurf();


    // Member Functions

        //- Return whether or not a phase can be considered continuous
        virtual bool canBeContinuous(const label index) const;

        //- Return whether or not this interface can segregate
        virtual bool canSegregate() const;
};


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

} // End namespace blendingMethods
} // End namespace Foam

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

#endif

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