fvcSurfaceOps.C

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
-------------------------------------------------------------------------------
    Copyright (C) 2024 M.Janssens
-------------------------------------------------------------------------------
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 "fvcSurfaceOps.H"
#include "fvMesh.H"
#include "fusedGaussLaplacianScheme.H"
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
namespace Foam
{
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
namespace fvc
{
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
template<class Type, class ResultType, class CellToFaceOp>
void surfaceSum
(
    const surfaceScalarField& lambdas,
    const GeometricField<Type, fvPatchField, volMesh>& vf,
    const CellToFaceOp& cop,
    GeometricField<ResultType, fvPatchField, volMesh>& result,
    const bool doCorrectBoundaryConditions
)
{
    const fvMesh& mesh = vf.mesh();
    const auto& Sf = mesh.Sf();
    const auto& P = mesh.owner();
    const auto& N = mesh.neighbour();
    const auto& pbm = mesh.boundaryMesh();
 
    const auto& vfi = vf.primitiveField();
    auto& sfi = result.primitiveFieldRef();
 
    // See e.g. surfaceInterpolationScheme<Type>::dotInterpolate
 
    // Fetch patchNeighbourField or value onto end of vf
    bitSet isCoupled(vf.boundaryField().size());
 
    const label oldSize = vf.constCast().boundaryEvaluate
    (
        [&](const auto& pfld, UList<Type>& slice)
        {
            if (pfld.coupled())
            {
                pfld.patchNeighbourField(slice);
                isCoupled.set(pfld.patch().index());
            }
            else
            {
                SubList<Type>(slice, pfld.size()) = pfld;
            }
        }
    );
 
    // Internal field
    const auto& Sfi = Sf.primitiveField();
    const auto& lambda = lambdas.primitiveField();
 
    #if 0
    {
        //- Example of cell-loop instead of face-loop
        const auto& cells = mesh.cells();
        const label nCells = mesh.nCells();
        const label nIntFaces = mesh.nInternalFaces();
        const auto& patchID = pbm.patchID();
        for (label celli = 0; celli < nCells; celli++)
        {
            const auto& cFaces = cells[celli];
            auto& res = sfi[celli];
            res = Zero;
            for (const label facei : cFaces)
            {
                if (facei < nIntFaces)
                {
                    ResultType faceVal
                    (
                        cop
                        (
                            Sfi[facei],
                            lambda[facei],
                            vfi[P[facei]],
                            vfi[N[facei]]
                        )
                    );
                    if (P[facei] != celli)
                    {
                        faceVal = -faceVal;
                    }
                    sfi[celli] += faceVal;
                }
                else
                {
                    const label patchi = patchID[facei-nIntFaces];
                    const label patchFacei = facei-pbm[patchi].start();
                    const label offset = nCells+pbm[patchi].offset();
                    const auto& pSf = Sf.boundaryField()[patchi];
                    const auto& pLambda = lambdas.boundaryField()[patchi];
 
                    if (isCoupled(patchi))
                    {
                        // Interpolate between owner-side and neighbour-side values
                        const ResultType faceVal
                        (
                            cop
                            (
                                pSf[patchFacei],
                                pLambda[patchFacei],
                                vfi[celli],
                                vfi[offset+patchFacei]
                            )
                        );
                        sfi[celli] += faceVal;
                    }
                    else if (pSf.size())
                    {
                        // Use patch value only
                        const ResultType faceVal
                        (
                            cop
                            (
                                pSf[patchFacei],
                                scalar(1.0),
                                vfi[offset+patchFacei],
                                pTraits<Type>::zero  // not used
                            )
                        );
                        sfi[celli] += faceVal;
                    }
                }
            }
        }
    }
    #else
    {
        for (label facei=0; facei<P.size(); facei++)
        {
            const label ownCelli = P[facei];
            const label neiCelli = N[facei];
 
            const ResultType faceVal
            (
                cop
                (
                    Sfi[facei],
                    lambda[facei],
                    vfi[ownCelli],
                    vfi[neiCelli]
                )
            );
            sfi[ownCelli] += faceVal;
            sfi[neiCelli] -= faceVal;
        }
 
        // Boundary field
        forAll(mesh.boundary(), patchi)
        {
            const auto& pFaceCells = mesh.boundary()[patchi].faceCells();
            const label offset = mesh.nCells()+pbm[patchi].offset();
            const auto& pSf = Sf.boundaryField()[patchi];
            const auto& pLambda = lambdas.boundaryField()[patchi];
 
            if (isCoupled(patchi))
            {
                //auto tpnf(pvf.patchNeighbourField());
                //auto& pnf = tpnf();
 
                for (label facei=0; facei<pFaceCells.size(); facei++)
                {
                    // Interpolate between owner-side and neighbour-side values
                    const ResultType faceVal
                    (
                        cop
                        (
                            pSf[facei],
                            pLambda[facei],
                            vfi[pFaceCells[facei]],
                            //pnf[facei]
                            vfi[offset+facei]
                        )
                    );
                    sfi[pFaceCells[facei]] += faceVal;
                }
            }
            else
            {
                for (label facei=0; facei<pFaceCells.size(); facei++)
                {
                    // Use patch value only
                    const ResultType faceVal
                    (
                        cop
                        (
                            pSf[facei],
                            scalar(1.0),
                            vfi[offset+facei],
                            pTraits<Type>::zero  // not used
                        )
                    );
                    sfi[pFaceCells[facei]] += faceVal;
                }
            }
        }
    }
    #endif
 
    // Restore original size
    vf.constCast().resize(oldSize);
 
    if (doCorrectBoundaryConditions)
    {
        result.correctBoundaryConditions();
    }
}
 
 
template<class Type, class FType, class ResultType, class CellToFaceOp>
void surfaceSum
(
    const surfaceScalarField& lambdas,
    const GeometricField<Type, fvPatchField, volMesh>& vf,
    const GeometricField<FType, fvsPatchField, surfaceMesh>& sadd,
    const CellToFaceOp& cop,
    GeometricField<ResultType, fvPatchField, volMesh>& result,
    const bool doCorrectBoundaryConditions
)
{
    const fvMesh& mesh = vf.mesh();
    const auto& Sf = mesh.Sf();
    const auto& P = mesh.owner();
    const auto& N = mesh.neighbour();
 
    const auto& vfi = vf.primitiveField();
    auto& sfi = result.primitiveFieldRef();
 
    // See e.g. surfaceInterpolationScheme<Type>::dotInterpolate
 
    // Internal field
    {
        const auto& Sfi = Sf.primitiveField();
        const auto& lambda = lambdas.primitiveField();
        const auto& saddi = sadd.primitiveField();
 
        for (label facei=0; facei<P.size(); facei++)
        {
            const label ownCelli = P[facei];
            const label neiCelli = N[facei];
 
            const ResultType faceVal
            (
                cop
                (
                    Sfi[facei],
                    lambda[facei],
                    vfi[ownCelli],
                    vfi[neiCelli],
 
                    saddi[facei]        // additional face value
                )
            );
            sfi[ownCelli] += faceVal;
            sfi[neiCelli] -= faceVal;
        }
    }
 
 
    // Boundary field
    {
        // Fetch patchNeighbourField onto end of GeoField
        const label oldSize = vf.constCast().boundaryEvaluate
        (
            [](const auto& pfld, UList<Type>& slice)
            {
                if (pfld.coupled())
                {
                    pfld.patchNeighbourField(slice);
                }
            }
        );
 
        forAll(mesh.boundary(), patchi)
        {
            const auto& pFaceCells = mesh.boundary()[patchi].faceCells();
            const auto& pSf = Sf.boundaryField()[patchi];
            const auto& pvf = vf.boundaryField()[patchi];
            const auto& pLambda = lambdas.boundaryField()[patchi];
            const auto& psadd = sadd.boundaryField()[patchi];
 
            if (pvf.coupled())
            {
                //auto tpnf(pvf.patchNeighbourField());
                //auto& pnf = tpnf();
 
                const label offset = mesh.nCells()+pvf.patch().offset();
 
                for (label facei=0; facei<pFaceCells.size(); facei++)
                {
                    // Interpolate between owner-side and neighbour-side values
                    const ResultType faceVal
                    (
                        cop
                        (
                            pSf[facei],
                            pLambda[facei],
                            vfi[pFaceCells[facei]],
                            //pnf[facei],
                            vfi[offset+facei],
                            psadd[facei]
                        )
                    );
                    sfi[pFaceCells[facei]] += faceVal;
                }
            }
            else
            {
                for (label facei=0; facei<pFaceCells.size(); facei++)
                {
                    // Use patch value only
                    const ResultType faceVal
                    (
                        cop
                        (
                            pSf[facei],
                            scalar(1.0),
                            pvf[facei],
                            pTraits<Type>::zero,  // not used
 
                            psadd[facei]
                        )
                    );
                    sfi[pFaceCells[facei]] += faceVal;
                }
            }
        }
 
        // Restore original size
        vf.constCast().resize(oldSize);
    }
 
    if (doCorrectBoundaryConditions)
    {
        result.correctBoundaryConditions();
    }
}
 
 
template
<
    class Type,
    class FType0,
    class FType1,
    class ResultType,
    class CellToFaceOp
>
void surfaceSum
(
    const surfaceScalarField& lambdas,
    const GeometricField<Type, fvPatchField, volMesh>& vf,
 
    const GeometricField<FType0, fvsPatchField, surfaceMesh>& sf0,
    const GeometricField<FType1, fvsPatchField, surfaceMesh>& sf1,
 
    const CellToFaceOp& cop,
    GeometricField<ResultType, fvPatchField, volMesh>& result,
    const bool doCorrectBoundaryConditions
)
{
    const fvMesh& mesh = vf.mesh();
    const auto& Sf = mesh.Sf();
    const auto& P = mesh.owner();
    const auto& N = mesh.neighbour();
 
    const auto& vfi = vf.primitiveField();
    auto& resulti = result.primitiveFieldRef();
 
    // See e.g. surfaceInterpolationScheme<Type>::dotInterpolate
 
    // Internal field
    {
        const auto& Sfi = Sf.primitiveField();
        const auto& lambda = lambdas.primitiveField();
        const auto& sf0i = sf0.primitiveField();
        const auto& sf1i = sf1.primitiveField();
 
        for (label facei=0; facei<P.size(); facei++)
        {
            const label ownCelli = P[facei];
            const label neiCelli = N[facei];
 
            const ResultType faceVal
            (
                cop
                (
                    Sfi[facei],
 
                    lambda[facei],
                    vfi[ownCelli],
                    vfi[neiCelli],
 
                    sf0i[facei],        // additional face value
                    sf1i[facei]         // additional face value
                )
            );
            resulti[ownCelli] += faceVal;
            resulti[neiCelli] -= faceVal;
        }
    }
 
 
    // Boundary field
    {
        // Fetch patchNeighbourField onto end of GeoField
        const label oldSize = vf.constCast().boundaryEvaluate
        (
            [](const auto& pfld, UList<Type>& slice)
            {
                if (pfld.coupled())
                {
                    pfld.patchNeighbourField(slice);
                }
            }
        );
 
        forAll(mesh.boundary(), patchi)
        {
            const auto& pFaceCells = mesh.boundary()[patchi].faceCells();
            const auto& pSf = Sf.boundaryField()[patchi];
            const auto& pvf = vf.boundaryField()[patchi];
            const auto& pLambda = lambdas.boundaryField()[patchi];
            const auto& psf0 = sf0.boundaryField()[patchi];
            const auto& psf1 = sf1.boundaryField()[patchi];
 
            if (pvf.coupled())
            {
                //auto tpnf(pvf.patchNeighbourField());
                //auto& pnf = tpnf();
 
                const label offset = mesh.nCells()+pvf.patch().offset();
 
                for (label facei=0; facei<pFaceCells.size(); facei++)
                {
                    // Interpolate between owner-side and neighbour-side values
                    const ResultType faceVal
                    (
                        cop
                        (
                            pSf[facei],
 
                            pLambda[facei],
                            vfi[pFaceCells[facei]],
                            //pnf[facei],
                            vfi[offset+facei],
 
                            psf0[facei],
                            psf1[facei]
                        )
                    );
 
                    resulti[pFaceCells[facei]] += faceVal;
                }
            }
            else
            {
                for (label facei=0; facei<pFaceCells.size(); facei++)
                {
                    // Use patch value only
                    const ResultType faceVal
                    (
                        cop
                        (
                            pSf[facei],
                            scalar(1.0),
                            pvf[facei],
                            pTraits<Type>::zero,  // not used
 
                            psf0[facei],
                            psf1[facei]
                        )
                    );
 
                    resulti[pFaceCells[facei]] += faceVal;
                }
            }
        }
 
        // Restore original size
        vf.constCast().resize(oldSize);
    }
 
    if (doCorrectBoundaryConditions)
    {
        result.correctBoundaryConditions();
    }
}
 
 
template<class Type, class ResultType, class CombineOp>
void GaussOp
(
    const surfaceScalarField& lambdas,
    const GeometricField<Type, fvPatchField, volMesh>& vf,
    const CombineOp& cop,
    GeometricField<ResultType, fvPatchField, volMesh>& result
)
{
    const fvMesh& mesh = vf.mesh();
 
    // Sum contributions
    surfaceSum(lambdas, vf, cop, result, false);
 
    auto& sfi = result.primitiveFieldRef();
    sfi /= mesh.V();
 
    result.correctBoundaryConditions();
}
 
 
template<class Type, class ResultType, class CombineOp>
void surfaceOp
(
    const GeometricField<Type, fvPatchField, volMesh>& vf,
    const surfaceVectorField& ownLs,
    const surfaceVectorField& neiLs,
    const CombineOp& cop,
    GeometricField<ResultType, fvPatchField, volMesh>& result
)
{
    const fvMesh& mesh = vf.mesh();
    const auto& Sf = mesh.Sf();
    const auto& P = mesh.owner();
    const auto& N = mesh.neighbour();
 
    const auto& vfi = vf.primitiveField();
    auto& sfi = result.primitiveFieldRef();
 
    // Internal field
    {
        const auto& Sfi = Sf.primitiveField();
 
        for (label facei=0; facei<P.size(); facei++)
        {
            const label ownCelli = P[facei];
            const label neiCelli = N[facei];
 
            const Type faceVal
            (
                cop
                (
                    Sfi[facei],     // needed?
                    vfi[ownCelli],
                    vfi[neiCelli]
                )
            );
            sfi[ownCelli] += ownLs[facei]*faceVal;
            sfi[neiCelli] -= neiLs[facei]*faceVal;
        }
    }
 
 
    // Boundary field
    {
        // Fetch patchNeighbourField onto end of GeoField
        const label oldSize = vf.constCast().boundaryEvaluate
        (
            [](const auto& pfld, UList<Type>& slice)
            {
                if (pfld.coupled())
                {
                    pfld.patchNeighbourField(slice);
                }
            }
        );
 
        forAll(mesh.boundary(), patchi)
        {
            const auto& pFaceCells = mesh.boundary()[patchi].faceCells();
            const auto& pSf = Sf.boundaryField()[patchi];
            const auto& pvf = vf.boundaryField()[patchi];
            const auto& pOwnLs = ownLs.boundaryField()[patchi];
 
            if (pvf.coupled())
            {
                //auto tpnf(pvf.patchNeighbourField());
                //auto& pnf = tpnf();
 
                const label offset = mesh.nCells()+pvf.patch().offset();
 
                for (label facei=0; facei<pFaceCells.size(); facei++)
                {
                    const Type faceVal
                    (
                        cop
                        (
                            pSf[facei], // needed?
                            vfi[pFaceCells[facei]],
                            //pnf[facei]
                            vfi[offset+facei]
                        )
                    );
                    sfi[pFaceCells[facei]] += pOwnLs[facei]*faceVal;
                }
            }
            else
            {
                for (label facei=0; facei<pFaceCells.size(); facei++)
                {
                    const Type faceVal
                    (
                        cop
                        (
                            pSf[facei],
                            vfi[pFaceCells[facei]],
                            pvf[facei]
                        )
                    );
                    sfi[pFaceCells[facei]] += pOwnLs[facei]*faceVal;
                }
            }
        }
 
        // Restore original size
        vf.constCast().resize(oldSize);
    }
 
    result.correctBoundaryConditions();
}
 
 
template<class Type, class ResultType, class CellToFaceOp>
void surfaceSnSum
(
    const surfaceScalarField& deltaCoeffs,
    const GeometricField<Type, fvPatchField, volMesh>& vf,
 
    const CellToFaceOp& cop,
 
    GeometricField<ResultType, fvPatchField, volMesh>& result,
    const bool doCorrectBoundaryConditions
)
{
    const fvMesh& mesh = vf.mesh();
    const auto& Sf = mesh.Sf();
    const auto& P = mesh.owner();
    const auto& N = mesh.neighbour();
 
    const auto& vfi = vf.primitiveField();
    auto& sfi = result.primitiveFieldRef();
 
    // Internal field
    {
        const auto& Sfi = Sf.primitiveField();
        const auto& dc = deltaCoeffs.primitiveField();
 
        for (label facei=0; facei<P.size(); facei++)
        {
            const label ownCelli = P[facei];
            const label neiCelli = N[facei];
 
            const ResultType faceVal
            (
                cop
                (
                    Sfi[facei],         // area vector
                    dc[facei],          // delta coefficients
                    vfi[ownCelli],
                    vfi[neiCelli]
                )
            );
            sfi[ownCelli] += faceVal;
            sfi[neiCelli] -= faceVal;
        }
    }
 
 
    // Boundary field
    {
        const label oldSize = vf.constCast().boundaryEvaluate
        (
            [](const auto& pfld, UList<Type>& slice)
            {
                if (pfld.coupled())
                {
                    pfld.patchNeighbourField(slice);
                }
                else
                {
                    pfld.snGrad(slice);
                }
            }
        );
 
        forAll(mesh.boundary(), patchi)
        {
            const auto& pFaceCells = mesh.boundary()[patchi].faceCells();
            const auto& pSf = Sf.boundaryField()[patchi];
            const auto& pvf = vf.boundaryField()[patchi];
            const auto& pdc = deltaCoeffs.boundaryField()[patchi];
 
            const label offset = mesh.nCells()+pvf.patch().offset();
 
            if (pvf.coupled())
            {
                //auto tpnf(pvf.patchNeighbourField());
                //auto& pnf = tpnf();
 
                for (label facei=0; facei<pFaceCells.size(); facei++)
                {
                    const label ownCelli = pFaceCells[facei];
 
                    // Interpolate between owner-side and neighbour-side values
                    const ResultType faceVal
                    (
                        cop
                        (
                            pSf[facei],         // area vector
                            pdc[facei],
                            vfi[ownCelli],
                            //pnf[facei]
                            vfi[offset+facei]
                        )
                    );
                    sfi[ownCelli] += faceVal;
                }
            }
            else
            {
                //auto tpnf(pvf.snGrad());
                //auto& pnf = tpnf();
 
                for (label facei=0; facei<pFaceCells.size(); facei++)
                {
                    const label ownCelli = pFaceCells[facei];
 
                    // Use patch value only
                    const ResultType faceVal
                    (
                        cop
                        (
                            pSf[facei],             // area vector
                            scalar(1.0),            // use 100% of pnf
                            pTraits<Type>::zero,
                            //pnf[facei]
                            vfi[offset+facei]
                        )
                    );
                    sfi[ownCelli] += faceVal;
                }
            }
        }
 
        // Restore original size
        vf.constCast().resize(oldSize);
    }
 
    if (doCorrectBoundaryConditions)
    {
        result.correctBoundaryConditions();
    }
}
 
 
template<class Type, class ResultType, class CellToFaceOp>
void surfaceSnSum
(
    const surfaceScalarField& deltaCoeffs,
    const GeometricField<Type, fvPatchField, volMesh>& vf,
    const GeometricField<Type, fvsPatchField, surfaceMesh>& sadd,
 
    const CellToFaceOp& cop,
 
    GeometricField<ResultType, fvPatchField, volMesh>& result,
    const bool doCorrectBoundaryConditions
)
{
    const fvMesh& mesh = vf.mesh();
    const auto& Sf = mesh.Sf();
    const auto& P = mesh.owner();
    const auto& N = mesh.neighbour();
 
    const auto& vfi = vf.primitiveField();
    auto& sfi = result.primitiveFieldRef();
 
    // Internal field
    {
        const auto& Sfi = Sf.primitiveField();
        const auto& dc = deltaCoeffs.primitiveField();
        const auto& saddi = sadd.primitiveField();
 
        for (label facei=0; facei<P.size(); facei++)
        {
            const label ownCelli = P[facei];
            const label neiCelli = N[facei];
 
            const ResultType faceVal
            (
                cop
                (
                    Sfi[facei],         // area vector
                    dc[facei],          // delta coefficients
                    vfi[ownCelli],
                    vfi[neiCelli],
                    saddi[facei]        // face addition value
                )
            );
            sfi[ownCelli] += faceVal;
            sfi[neiCelli] -= faceVal;
        }
    }
 
 
    // Boundary field
    {
        // Fetch patchNeighbourField onto end of GeoField
        const label oldSize = vf.constCast().boundaryEvaluate
        (
            [](const auto& pfld, UList<Type>& slice)
            {
                if (pfld.coupled())
                {
                    pfld.patchNeighbourField(slice);
                }
                else
                {
                    pfld.snGrad(slice);
                }
            }
        );
 
        forAll(mesh.boundary(), patchi)
        {
            const auto& pFaceCells = mesh.boundary()[patchi].faceCells();
            const auto& pSf = Sf.boundaryField()[patchi];
            const auto& pvf = vf.boundaryField()[patchi];
            const auto& pdc = deltaCoeffs.boundaryField()[patchi];
            const auto& psadd = sadd.boundaryField()[patchi];
 
            const label offset = mesh.nCells()+pvf.patch().offset();
 
            if (pvf.coupled())
            {
                //auto tpnf(pvf.patchNeighbourField());
                //auto& pnf = tpnf();
 
                for (label facei=0; facei<pFaceCells.size(); facei++)
                {
                    const label ownCelli = pFaceCells[facei];
 
                    // Interpolate between owner-side and neighbour-side values
                    const ResultType faceVal
                    (
                        cop
                        (
                            pSf[facei],         // area vector
                            pdc[facei],
                            vfi[ownCelli],
                            //pnf[facei],
                            vfi[offset+facei],
                            psadd[facei]
                        )
                    );
                    sfi[ownCelli] += faceVal;
                }
            }
            else
            {
                //auto tpnf(pvf.snGrad());
                //auto& pnf = tpnf();
 
                for (label facei=0; facei<pFaceCells.size(); facei++)
                {
                    const label ownCelli = pFaceCells[facei];
 
                    // Use patch value only
                    const ResultType faceVal
                    (
                        cop
                        (
                            pSf[facei],             // area vector
                            scalar(1.0),            // use 100% of pnf
                            pTraits<Type>::zero,
                            //pnf[facei],
                            vfi[offset+facei],
                            psadd[facei]
                        )
                    );
                    sfi[ownCelli] += faceVal;
                }
            }
        }
 
        // Restore original size
        vf.constCast().resize(oldSize);
    }
 
    if (doCorrectBoundaryConditions)
    {
        result.correctBoundaryConditions();
    }
}
 
 
template<class Type, class GType, class ResultType, class CellToFaceOp>
void surfaceSnSum
(
    const surfaceScalarField& gammaWeights,
    const GeometricField<GType, fvPatchField, volMesh>& gamma,
 
    const surfaceScalarField& deltaCoeffs,
    const GeometricField<Type, fvPatchField, volMesh>& vf,
 
    const CellToFaceOp& cop,
 
    GeometricField<ResultType, fvPatchField, volMesh>& result,
    const bool doCorrectBoundaryConditions
)
{
    const fvMesh& mesh = vf.mesh();
    const auto& Sf = mesh.Sf();
    const auto& P = mesh.owner();
    const auto& N = mesh.neighbour();
 
    const auto& gammai = gamma.primitiveField();
    const auto& vfi = vf.primitiveField();
    auto& sfi = result.primitiveFieldRef();
 
    // Internal field
    {
        const auto& Sfi = Sf.primitiveField();
        const auto& weights = gammaWeights.primitiveField();
        const auto& dc = deltaCoeffs.primitiveField();
 
        for (label facei=0; facei<P.size(); facei++)
        {
            const label ownCelli = P[facei];
            const label neiCelli = N[facei];
 
            const ResultType faceVal
            (
                cop
                (
                    Sfi[facei],         // area vector
 
                    weights[facei],     // interpolation weights
                    gammai[ownCelli],
                    gammai[neiCelli],
 
                    dc[facei],          // delta coefficients
                    vfi[ownCelli],
                    vfi[neiCelli]
                )
            );
            sfi[ownCelli] += faceVal;
            sfi[neiCelli] -= faceVal;
        }
    }
 
 
    // Boundary field
    {
        // Fetch patchNeighbourField onto end of GeoField
        const label oldSize = vf.constCast().boundaryEvaluate
        (
            [](const auto& pfld, UList<Type>& slice)
            {
                if (pfld.coupled())
                {
                    pfld.patchNeighbourField(slice);
                }
                else
                {
                    pfld.snGrad(slice);
                }
            }
        );
 
        (void) gamma.constCast().boundaryEvaluate
        (
            [](const auto& pfld, UList<GType>& slice)
            {
                if (pfld.coupled())
                {
                    pfld.patchNeighbourField(slice);
                }
            }
        );
 
        forAll(mesh.boundary(), patchi)
        {
            const auto& pFaceCells = mesh.boundary()[patchi].faceCells();
            const auto& pSf = Sf.boundaryField()[patchi];
            const auto& pvf = vf.boundaryField()[patchi];
            const auto& pdc = deltaCoeffs.boundaryField()[patchi];
            const auto& pweights = gammaWeights.boundaryField()[patchi];
            const auto& pgamma = gamma.boundaryField()[patchi];
 
            const label offset = mesh.nCells()+pvf.patch().offset();
 
            if (pvf.coupled())
            {
                //auto tpnf(pvf.patchNeighbourField());
                //auto& pnf = tpnf();
                //auto tgammanf(pgamma.patchNeighbourField());
                //auto& gammanf = tgammanf();
 
                for (label facei=0; facei<pFaceCells.size(); facei++)
                {
                    const label ownCelli = pFaceCells[facei];
 
                    // Interpolate between owner-side and neighbour-side values
                    const ResultType faceVal
                    (
                        cop
                        (
                            pSf[facei],         // area vector
 
                            pweights[facei],
                            gammai[ownCelli],
                            //gammanf[facei],
                            gammai[offset+facei],
 
                            pdc[facei],
                            vfi[ownCelli],
                            //pnf[facei]
                            vfi[offset+facei]
                        )
                    );
                    sfi[ownCelli] += faceVal;
                }
            }
            else
            {
                //auto tpnf(pvf.snGrad());
                //auto& pnf = tpnf();
 
                for (label facei=0; facei<pFaceCells.size(); facei++)
                {
                    const label ownCelli = pFaceCells[facei];
 
                    // Use patch value only
                    const ResultType faceVal
                    (
                        cop
                        (
                            pSf[facei],             // area vector
 
                            scalar(1.0),
                            pgamma[facei],
                            pTraits<GType>::zero,   // not used
 
 
                            scalar(1.0),            // use 100% of pnf
                            pTraits<Type>::zero,
                            //pnf[facei]
                            vfi[offset+facei]
                        )
                    );
                    sfi[ownCelli] += faceVal;
                }
            }
        }
 
        // Restore original size
        vf.constCast().resize(oldSize);
        gamma.constCast().resize(oldSize);
    }
 
    if (doCorrectBoundaryConditions)
    {
        result.correctBoundaryConditions();
    }
}
 
 
template<class Type, class GType, class ResultType, class CellToFaceOp>
void surfaceSnSum
(
    const surfaceScalarField& gammaWeights,
    const GeometricField<GType, fvPatchField, volMesh>& gamma,
 
    const surfaceScalarField& deltaCoeffs,
    const GeometricField<Type, fvPatchField, volMesh>& vf,
 
    const GeometricField<Type, fvsPatchField, surfaceMesh>& sadd,
 
    const CellToFaceOp& cop,
 
    GeometricField<ResultType, fvPatchField, volMesh>& result,
    const bool doCorrectBoundaryConditions
)
{
    const fvMesh& mesh = vf.mesh();
    const auto& Sf = mesh.Sf();
    const auto& P = mesh.owner();
    const auto& N = mesh.neighbour();
 
    const auto& gammai = gamma.primitiveField();
    const auto& vfi = vf.primitiveField();
    auto& sfi = result.primitiveFieldRef();
 
    // Internal field
    {
        const auto& Sfi = Sf.primitiveField();
        const auto& weights = gammaWeights.primitiveField();
        const auto& dc = deltaCoeffs.primitiveField();
        const auto& saddi = sadd.primitiveField();
 
        for (label facei=0; facei<P.size(); facei++)
        {
            const label ownCelli = P[facei];
            const label neiCelli = N[facei];
 
            const ResultType faceVal
            (
                cop
                (
                    Sfi[facei],         // area vector
 
                    weights[facei],     // interpolation weights
                    gammai[ownCelli],
                    gammai[neiCelli],
 
                    dc[facei],          // delta coefficients
                    vfi[ownCelli],
                    vfi[neiCelli],
 
                    saddi[facei]        // face addition value
                )
            );
            sfi[ownCelli] += faceVal;
            sfi[neiCelli] -= faceVal;
        }
    }
 
 
    // Boundary field
    {
        // Fetch patchNeighbourField onto end of GeoField
        const label oldSize = vf.constCast().boundaryEvaluate
        (
            [](const auto& pfld, UList<Type>& slice)
            {
                if (pfld.coupled())
                {
                    pfld.patchNeighbourField(slice);
                }
                else
                {
                    pfld.snGrad(slice);
                }
            }
        );
 
        (void) gamma.constCast().boundaryEvaluate
        (
            [](const auto& pfld, UList<GType>& slice)
            {
                if (pfld.coupled())
                {
                    pfld.patchNeighbourField(slice);
                }
            }
        );
 
        forAll(mesh.boundary(), patchi)
        {
            const auto& pFaceCells = mesh.boundary()[patchi].faceCells();
            const auto& pSf = Sf.boundaryField()[patchi];
            const auto& pvf = vf.boundaryField()[patchi];
            const auto& pdc = deltaCoeffs.boundaryField()[patchi];
            const auto& pweights = gammaWeights.boundaryField()[patchi];
            const auto& pgamma = gamma.boundaryField()[patchi];
            const auto& psadd = sadd.boundaryField()[patchi];
 
            const label offset = mesh.nCells()+pvf.patch().offset();
 
            if (pvf.coupled())
            {
                for (label facei=0; facei<pFaceCells.size(); facei++)
                {
                    const label ownCelli = pFaceCells[facei];
 
                    // Interpolate between owner-side and neighbour-side values
                    const ResultType faceVal
                    (
                        cop
                        (
                            pSf[facei],         // area vector
 
                            pweights[facei],
                            gammai[ownCelli],
                            //gammanf[facei],
                            gammai[offset+facei],
 
                            pdc[facei],
                            vfi[ownCelli],
                            //pnf[facei],
                            vfi[offset+facei],
 
                            psadd[facei]
                        )
                    );
                    sfi[ownCelli] += faceVal;
                }
            }
            else
            {
                //auto tpnf(pvf.snGrad());
                //auto& pnf = tpnf();
 
                for (label facei=0; facei<pFaceCells.size(); facei++)
                {
                    const label ownCelli = pFaceCells[facei];
 
                    // Use patch value only
                    const ResultType faceVal
                    (
                        cop
                        (
                            pSf[facei],             // area vector
 
                            scalar(1.0),
                            pgamma[facei],
                            pTraits<GType>::zero,   // not used
 
 
                            scalar(1.0),            // use 100% of pnf
                            pTraits<Type>::zero,
                            //pnf[facei],
                            vfi[offset+facei],
 
                            psadd[facei]
                        )
                    );
                    sfi[ownCelli] += faceVal;
                }
            }
        }
 
        // Restore original size
        vf.constCast().resize(oldSize);
        gamma.constCast().resize(oldSize);
    }
 
    if (doCorrectBoundaryConditions)
    {
        result.correctBoundaryConditions();
    }
}
 
 
template
<
    class Type,
    class GType0,
    class GType1,
    class ResultType,
    class CellToFaceOp
>
void surfaceSnSum
(
    const surfaceScalarField& gammaWeights,
    const GeometricField<GType0, fvPatchField, volMesh>& gamma0,
    const GeometricField<GType1, fvPatchField, volMesh>& gamma1,
 
    const surfaceScalarField& deltaCoeffs,
    const GeometricField<Type, fvPatchField, volMesh>& vf,
 
    const CellToFaceOp& cop,
 
    GeometricField<ResultType, fvPatchField, volMesh>& result,
    const bool doCorrectBoundaryConditions
)
{
    const fvMesh& mesh = vf.mesh();
    const auto& Sf = mesh.Sf();
    const auto& P = mesh.owner();
    const auto& N = mesh.neighbour();
 
    const auto& gamma0i = gamma0.primitiveField();
    const auto& gamma1i = gamma1.primitiveField();
    const auto& vfi = vf.primitiveField();
    auto& sfi = result.primitiveFieldRef();
 
    // Internal field
    {
        const auto& Sfi = Sf.primitiveField();
        const auto& weights = gammaWeights.primitiveField();
        const auto& dc = deltaCoeffs.primitiveField();
 
        for (label facei=0; facei<P.size(); facei++)
        {
            const label ownCelli = P[facei];
            const label neiCelli = N[facei];
 
            const ResultType faceVal
            (
                cop
                (
                    Sfi[facei],         // area vector
 
                    weights[facei],     // interpolation weights
 
                    gamma0i[ownCelli],
                    gamma0i[neiCelli],
 
                    gamma1i[ownCelli],
                    gamma1i[neiCelli],
 
                    dc[facei],          // delta coefficients
                    vfi[ownCelli],
                    vfi[neiCelli]
                )
            );
            sfi[ownCelli] += faceVal;
            sfi[neiCelli] -= faceVal;
        }
    }
 
 
    // Boundary field
    {
        // Fetch patchNeighbourField onto end of GeoField
        const label oldSize = vf.constCast().boundaryEvaluate
        (
            [](const auto& pfld, UList<Type>& slice)
            {
                if (pfld.coupled())
                {
                    pfld.patchNeighbourField(slice);
                }
                else
                {
                    pfld.snGrad(slice);
                }
            }
        );
 
        (void) gamma0.constCast().boundaryEvaluate
        (
            [](const auto& pfld, UList<GType0>& slice)
            {
                if (pfld.coupled())
                {
                    pfld.patchNeighbourField(slice);
                }
            }
        );
 
        (void) gamma1.constCast().boundaryEvaluate
        (
            [](const auto& pfld, UList<GType1>& slice)
            {
                if (pfld.coupled())
                {
                    pfld.patchNeighbourField(slice);
                }
            }
        );
 
        forAll(mesh.boundary(), patchi)
        {
            const auto& pFaceCells = mesh.boundary()[patchi].faceCells();
            const auto& pSf = Sf.boundaryField()[patchi];
            const auto& pvf = vf.boundaryField()[patchi];
            const auto& pdc = deltaCoeffs.boundaryField()[patchi];
            const auto& pweights = gammaWeights.boundaryField()[patchi];
            const auto& pgamma0 = gamma0.boundaryField()[patchi];
            const auto& pgamma1 = gamma1.boundaryField()[patchi];
 
            const label offset = mesh.nCells()+pvf.patch().offset();
 
            if (pvf.coupled())
            {
                //auto tpnf(pvf.patchNeighbourField());
                //auto& pnf = tpnf();
                //auto tgamma0nf(pgamma0.patchNeighbourField());
                //auto& gamma0nf = tgamma0nf();
                //auto tgamma1nf(pgamma1.patchNeighbourField());
                //auto& gamma1nf = tgamma1nf();
 
                for (label facei=0; facei<pFaceCells.size(); facei++)
                {
                    const label ownCelli = pFaceCells[facei];
 
                    // Interpolate between owner-side and neighbour-side values
                    const ResultType faceVal
                    (
                        cop
                        (
                            pSf[facei],         // area vector
 
                            pweights[facei],
 
                            gamma0i[ownCelli],
                            gamma0i[offset+facei],  //gamma0nf[facei],
 
                            gamma1i[ownCelli],
                            gamma1i[offset+facei],  //gamma1nf[facei],
 
                            pdc[facei],
                            vfi[ownCelli],
                            vfi[offset+facei]       //pnf[facei]
                        )
                    );
                    sfi[ownCelli] += faceVal;
                }
            }
            else
            {
                //auto tpnf(pvf.snGrad());
                //auto& pnf = tpnf();
                for (label facei=0; facei<pFaceCells.size(); facei++)
                {
                    const label ownCelli = pFaceCells[facei];
 
                    // Use patch value only
                    const ResultType faceVal
                    (
                        cop
                        (
                            pSf[facei],             // area vector
 
                            scalar(1.0),
 
                            pgamma0[facei],
                            pTraits<GType0>::zero,  // not used
 
                            pgamma1[facei],
                            pTraits<GType1>::zero,  // not used
 
                            scalar(1.0),            // use 100% of pnf
                            pTraits<Type>::zero,
                            vfi[offset+facei]       //pnf[facei]
                        )
                    );
                    sfi[ownCelli] += faceVal;
                }
            }
        }
 
        // Restore original size
        vf.constCast().resize(oldSize);
        gamma0.constCast().resize(oldSize);
        gamma1.constCast().resize(oldSize);
    }
 
    if (doCorrectBoundaryConditions)
    {
        result.correctBoundaryConditions();
    }
}
 
 
template<class Type, class FType, class ResultType, class CellToFaceOp>
void interpolate
(
    const surfaceScalarField& weights,
    const GeometricField<Type, fvPatchField, volMesh>& vf,
    const GeometricField<FType, fvsPatchField, surfaceMesh>& sf,
    const CellToFaceOp& cop,
    GeometricField<ResultType, fvsPatchField, surfaceMesh>& result
)
{
    const fvMesh& mesh = vf.mesh();
    const auto& Sf = mesh.Sf();
    const auto& P = mesh.owner();
    const auto& N = mesh.neighbour();
 
    const auto& vfi = vf.primitiveField();
 
    // Internal field
    {
        const auto& Sfi = Sf.primitiveField();
        const auto& weight = weights.primitiveField();
        const auto& sfi = sf.primitiveField();
 
        auto& resulti = result.primitiveFieldRef();
 
        for (label facei=0; facei<P.size(); facei++)
        {
            const label ownCelli = P[facei];
            const label neiCelli = N[facei];
 
            cop
            (
                Sfi[facei],
 
                weight[facei],
                vfi[ownCelli],
                vfi[neiCelli],
 
                sfi[facei],
 
                resulti[facei]
            );
        }
    }
 
 
    // Boundary field
    {
        // Fetch patchNeighbourField onto end of GeoField
        const label oldSize = vf.constCast().boundaryEvaluate
        (
            [](const auto& pfld, UList<Type>& slice)
            {
                if (pfld.coupled())
                {
                    pfld.patchNeighbourField(slice);
                }
            }
        );
 
        forAll(mesh.boundary(), patchi)
        {
            const auto& pFaceCells = mesh.boundary()[patchi].faceCells();
            const auto& pSf = Sf.boundaryField()[patchi];
            const auto& pvf = vf.boundaryField()[patchi];
            const auto& pweight = weights.boundaryField()[patchi];
            const auto& psf = sf.boundaryField()[patchi];
            auto& presult = result.boundaryFieldRef()[patchi];
 
            const label offset = mesh.nCells()+pvf.patch().offset();
 
            if (pvf.coupled())
            {
                //auto tpnf(pvf.patchNeighbourField());
                //auto& pnf = tpnf();
 
                for (label facei=0; facei<pFaceCells.size(); facei++)
                {
                    // Interpolate between owner-side and neighbour-side values
                    cop
                    (
                        pSf[facei],
 
                        pweight[facei],
                        vfi[pFaceCells[facei]],
                        vfi[offset+facei],  //pnf[facei],
 
                        psf[facei],
 
                        presult[facei]
                    );
                }
            }
            else
            {
                for (label facei=0; facei<pFaceCells.size(); facei++)
                {
                    // Use patch value only
                    cop
                    (
                        pSf[facei],
 
                        scalar(1.0),
                        pvf[facei],
                        pTraits<Type>::zero,    // not used
 
                        psf[facei],
 
                        presult[facei]
                    );
                }
            }
        }
 
        // Restore original size
        vf.constCast().resize(oldSize);
    }
}
 
 
template
<
    class Type0,
    class Type1,
    class ResultType,
    class CellToFaceOp
>
void interpolate
(
    const surfaceScalarField& weights,
    const GeometricField<Type0, fvPatchField, volMesh>& vf0,
    const GeometricField<Type1, fvPatchField, volMesh>& vf1,
    const CellToFaceOp& cop,
    GeometricField<ResultType, fvsPatchField, surfaceMesh>& result
)
{
    const fvMesh& mesh = vf0.mesh();
    const auto& Sf = mesh.Sf();
    const auto& P = mesh.owner();
    const auto& N = mesh.neighbour();
 
    const auto& vf0i = vf0.primitiveField();
    const auto& vf1i = vf1.primitiveField();
 
    // Internal field
    {
        const auto& Sfi = Sf.primitiveField();
        const auto& weight = weights.primitiveField();
 
        auto& resulti = result.primitiveFieldRef();
 
        for (label facei=0; facei<P.size(); facei++)
        {
            const label ownCelli = P[facei];
            const label neiCelli = N[facei];
 
            cop
            (
                Sfi[facei],
 
                weight[facei],
 
                vf0i[ownCelli],
                vf0i[neiCelli],
 
                vf1i[ownCelli],
                vf1i[neiCelli],
 
                resulti[facei]
            );
        }
    }
 
 
    // Boundary field
    {
        // Fetch patchNeighbourField onto end of GeoField
        const label oldSize = vf0.constCast().boundaryEvaluate
        (
            [](const auto& pfld, UList<Type0>& slice)
            {
                if (pfld.coupled())
                {
                    pfld.patchNeighbourField(slice);
                }
            }
        );
 
        // Fetch patchNeighbourField onto end of GeoField
        (void) vf1.constCast().boundaryEvaluate
        (
            [](const auto& pfld, UList<Type1>& slice)
            {
                if (pfld.coupled())
                {
                    pfld.patchNeighbourField(slice);
                }
            }
        );
 
        forAll(mesh.boundary(), patchi)
        {
            const auto& pFaceCells = mesh.boundary()[patchi].faceCells();
            const auto& pSf = Sf.boundaryField()[patchi];
            const auto& pvf0 = vf0.boundaryField()[patchi];
            const auto& pvf1 = vf1.boundaryField()[patchi];
            const auto& pweight = weights.boundaryField()[patchi];
            auto& presult = result.boundaryFieldRef()[patchi];
 
            const label offset = mesh.nCells()+pvf0.patch().offset();
 
            if (pvf0.coupled() || pvf1.coupled())
            {
                //auto tpnf0(pvf0.patchNeighbourField());
                //auto& pnf0 = tpnf0();
                //auto tpnf1(pvf1.patchNeighbourField());
                //auto& pnf1 = tpnf1();
 
                for (label facei=0; facei<pFaceCells.size(); facei++)
                {
                    // Interpolate between owner-side and neighbour-side values
                    cop
                    (
                        pSf[facei],
 
                        pweight[facei],
 
                        vf0i[pFaceCells[facei]],
                        vf0i[offset+facei], //pnf0[facei],
 
                        vf1i[pFaceCells[facei]],
                        vf1i[offset+facei], //pnf1[facei],
 
                        presult[facei]
                    );
                }
            }
            else
            {
                for (label facei=0; facei<pFaceCells.size(); facei++)
                {
                    // Use patch value only
                    cop
                    (
                        pSf[facei],
 
                        scalar(1.0),
 
                        pvf0[facei],
                        pTraits<Type0>::zero, // not used
 
                        pvf1[facei],
                        pTraits<Type1>::zero, // not used
 
                        presult[facei]
                    );
                }
            }
        }
 
        // Restore original size
        vf0.constCast().resize(oldSize);
        vf1.constCast().resize(oldSize);
    }
}
 
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
} // End namespace fvc
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
} // End namespace Foam
 
// ************************************************************************* //