ENH: fused: 'fused' variants of explicit discretisation

src/Allwmake

@@ -55,6 +55,7 @@ wmake $targetType meshTools
 
 wmake $targetType finiteArea
 wmake $targetType finiteVolume
+wmake $targetType fused/finiteVolume
 
 wmake $targetType mesh/blockMesh
 wmake $targetType mesh/extrudeModel  # Requires: blockMesh

src/fused/finiteVolume/Make/files

@@ -0,0 +1,7 @@
+fusedGaussLaplacianSchemes.C
+fusedLeastSquaresGrads.C
+fusedGaussDivSchemes.C
+fusedGaussConvectionSchemes.C
+fusedGaussGrads.C
+
+LIB = $(FOAM_LIBBIN)/libfusedFiniteVolume

src/fused/finiteVolume/Make/options

@@ -0,0 +1,7 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
+
+LIB_LIBS = \
+    -lfiniteVolume

src/fused/finiteVolume/fusedGaussConvectionScheme.C

@@ -0,0 +1,265 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2016 OpenFOAM Foundation
+    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 "fusedGaussConvectionScheme.H"
+#include "fvcSurfaceIntegrate.H"
+#include "fvMatrices.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace fv
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+const surfaceInterpolationScheme<Type>&
+fusedGaussConvectionScheme<Type>::interpScheme() const
+{
+    return tinterpScheme_();
+}
+
+
+template<class Type>
+tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
+fusedGaussConvectionScheme<Type>::interpolate
+(
+    const surfaceScalarField&,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+) const
+{
+    return tinterpScheme_().interpolate(vf);
+}
+
+
+template<class Type>
+tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
+fusedGaussConvectionScheme<Type>::flux
+(
+    const surfaceScalarField& faceFlux,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+) const
+{
+    return faceFlux*interpolate(faceFlux, vf);
+}
+
+
+template<class Type>
+tmp<fvMatrix<Type>>
+fusedGaussConvectionScheme<Type>::fvmDiv
+(
+    const surfaceScalarField& faceFlux,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+) const
+{
+    DebugPout<< "fusedGaussConvectionScheme<Type>::fvmDiv on " << vf.name()
+        << " with flux " << faceFlux.name() << endl;
+
+    tmp<surfaceScalarField> tweights = tinterpScheme_().weights(vf);
+    const surfaceScalarField& weights = tweights();
+
+    tmp<fvMatrix<Type>> tfvm
+    (
+        new fvMatrix<Type>
+        (
+            vf,
+            faceFlux.dimensions()*vf.dimensions()
+        )
+    );
+    fvMatrix<Type>& fvm = tfvm.ref();
+
+    //fvm.lower() = -weights.primitiveField()*faceFlux.primitiveField();
+    multiplySubtract
+    (
+        fvm.lower(),
+        weights.primitiveField(),
+        faceFlux.primitiveField()
+    );
+
+    //fvm.upper() = fvm.lower() + faceFlux.primitiveField();
+    add(fvm.upper(), fvm.lower(), faceFlux.primitiveField());
+
+    fvm.negSumDiag();
+
+    forAll(vf.boundaryField(), patchi)
+    {
+        const fvPatchField<Type>& psf = vf.boundaryField()[patchi];
+        const fvsPatchScalarField& patchFlux = faceFlux.boundaryField()[patchi];
+        const fvsPatchScalarField& pw = weights.boundaryField()[patchi];
+
+        auto& intCoeffs = fvm.internalCoeffs()[patchi];
+        auto& bouCoeffs = fvm.boundaryCoeffs()[patchi];
+
+        //fvm.internalCoeffs()[patchi] = patchFlux*psf.valueInternalCoeffs(pw);
+        multiply(intCoeffs, patchFlux, psf.valueInternalCoeffs(pw)());
+
+        //fvm.boundaryCoeffs()[patchi] = -patchFlux*psf.valueBoundaryCoeffs(pw);
+        multiply(bouCoeffs, patchFlux, psf.valueBoundaryCoeffs(pw)());
+        bouCoeffs.negate();
+    }
+
+    if (tinterpScheme_().corrected())
+    {
+        fvm += fvc::surfaceIntegrate(faceFlux*tinterpScheme_().correction(vf));
+    }
+
+    return tfvm;
+}
+
+
+template<class Type>
+tmp<GeometricField<Type, fvPatchField, volMesh>>
+fusedGaussConvectionScheme<Type>::fvcDiv
+(
+    const surfaceScalarField& faceFlux,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+) const
+{
+    DebugPout<< "fusedGaussConvectionScheme<Type>::fvcDiv on " << vf.name()
+        << " with flux " << faceFlux.name() << endl;
+
+    typedef GeometricField<Type, fvPatchField, volMesh> FieldType;
+
+    const fvMesh& mesh = vf.mesh();
+
+    tmp<FieldType> tConvection
+    (
+        new FieldType
+        (
+            IOobject
+            (
+                "convection(" + faceFlux.name() + ',' + vf.name() + ')',
+                vf.instance(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimensioned<Type>
+            (
+                faceFlux.dimensions()
+               *vf.dimensions()
+               /dimVol,
+                Zero
+            ),
+            fvPatchFieldBase::extrapolatedCalculatedType()
+        )
+    );
+
+    if (this->tinterpScheme_().corrected())
+    {
+        const auto tfaceCorr(this->tinterpScheme_().correction(vf));
+        auto& faceCorr = tfaceCorr();
+
+        const auto interpolator = [&]
+        (
+            const vector& area,
+            const scalar lambda,
+
+            const Type& ownVal,
+            const Type& neiVal,
+
+            const scalar& faceVal,
+            const Type& correction
+
+        ) -> Type
+        {
+            return faceVal*((lambda*(ownVal - neiVal) + neiVal) + correction);
+        };
+
+        fvc::surfaceSum
+        (
+            // interpolation factors for volume field
+            this->tinterpScheme_().weights(vf),
+
+            // volume field(s)
+            vf,
+
+            // surface field(s)
+            faceFlux,
+            faceCorr,
+
+            // operation
+            interpolator,
+
+            tConvection.ref(),
+            false
+        );
+    }
+    else
+    {
+        const auto interpolator = [&]
+        (
+            const vector& area,
+            const scalar lambda,
+            const Type& ownVal,
+            const Type& neiVal,
+
+            const scalar& faceVal
+        ) -> Type
+        {
+            return faceVal*(lambda*(ownVal - neiVal) + neiVal);
+        };
+
+        fvc::surfaceSum
+        (
+            tinterpScheme_().weights(vf),
+
+            vf,
+
+            faceFlux,
+
+            interpolator,
+            tConvection.ref(),
+            false
+        );
+    }
+
+    tConvection.ref().primitiveFieldRef() /= mesh.Vsc();
+
+    tConvection.ref().correctBoundaryConditions();
+
+    return tConvection;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/fused/finiteVolume/fusedGaussConvectionScheme.H

@@ -0,0 +1,162 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2016 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+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/>.
+
+Class
+    Foam::fv::fusedGaussConvectionScheme
+
+Group
+    grpFvConvectionSchemes
+
+Description
+    Basic second-order convection using face-gradients and Gauss' theorem.
+
+SourceFiles
+    fusedGaussConvectionScheme.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef fusedGaussConvectionScheme_H
+#define fusedGaussConvectionScheme_H
+
+#include "convectionScheme.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace fv
+{
+
+/*---------------------------------------------------------------------------*\
+                 Class fusedGaussConvectionScheme Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class fusedGaussConvectionScheme
+:
+    public fv::convectionScheme<Type>
+{
+    // Private data
+
+        tmp<surfaceInterpolationScheme<Type>> tinterpScheme_;
+
+
+    // Private Member Functions
+
+        //- No copy construct
+        fusedGaussConvectionScheme(const fusedGaussConvectionScheme&) = delete;
+
+        //- No copy assignment
+        void operator=(const fusedGaussConvectionScheme&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("fusedGauss");
+
+
+    // Constructors
+
+        //- Construct from flux and interpolation scheme
+        fusedGaussConvectionScheme
+        (
+            const fvMesh& mesh,
+            const surfaceScalarField& faceFlux,
+            const tmp<surfaceInterpolationScheme<Type>>& scheme
+        )
+        :
+            convectionScheme<Type>(mesh, faceFlux),
+            tinterpScheme_(scheme)
+        {}
+
+        //- Construct from flux and Istream
+        fusedGaussConvectionScheme
+        (
+            const fvMesh& mesh,
+            const surfaceScalarField& faceFlux,
+            Istream& is
+        )
+        :
+            convectionScheme<Type>(mesh, faceFlux),
+            tinterpScheme_
+            (
+                surfaceInterpolationScheme<Type>::New(mesh, faceFlux, is)
+            )
+        {}
+
+
+    // Member Functions
+
+        const surfaceInterpolationScheme<Type>& interpScheme() const;
+
+        tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> interpolate
+        (
+            const surfaceScalarField&,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        ) const;
+
+        tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> flux
+        (
+            const surfaceScalarField&,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        ) const;
+
+        tmp<fvMatrix<Type>> fvmDiv
+        (
+            const surfaceScalarField&,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        ) const;
+
+        tmp<GeometricField<Type, fvPatchField, volMesh>> fvcDiv
+        (
+            const surfaceScalarField&,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "fusedGaussConvectionScheme.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/fused/finiteVolume/fusedGaussConvectionSchemes.C

@@ -0,0 +1,35 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2016 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+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 "fusedGaussConvectionScheme.H"
+#include "fvMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makeFvConvectionScheme(fusedGaussConvectionScheme)
+
+// ************************************************************************* //

src/fused/finiteVolume/fusedGaussDivScheme.C

@@ -0,0 +1,155 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2016 OpenFOAM Foundation
+    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 "fusedGaussDivScheme.H"
+#include "fvcSurfaceIntegrate.H"
+#include "fvMatrices.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace fv
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+tmp
+<
+    GeometricField
+    <typename innerProduct<vector, Type>::type, fvPatchField, volMesh>
+>
+fusedGaussDivScheme<Type>::fvcDiv
+(
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    typedef typename innerProduct<vector, Type>::type DivType;
+    typedef GeometricField<DivType, fvPatchField, volMesh> DivFieldType;
+
+    const fvMesh& mesh = vf.mesh();
+
+    DebugPout<< "fusedGaussDivScheme<Type>::fvcDiv on " << vf.name()
+         << endl;
+
+    tmp<DivFieldType> tDiv
+    (
+        new DivFieldType
+        (
+            IOobject
+            (
+                "div(" + vf.name() + ')',
+                vf.instance(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimensioned<DivType>(vf.dimensions()/dimLength, Zero),
+            fvPatchFieldBase::extrapolatedCalculatedType()
+        )
+    );
+    DivFieldType& div = tDiv.ref();
+
+    //fvcDiv(vf, this->tinterpScheme_().weights(vf), div);
+
+    if (this->tinterpScheme_().corrected())
+    {
+        const auto tfaceCorr(this->tinterpScheme_().correction(vf));
+        auto& faceCorr = tfaceCorr();
+
+        const auto dotInterpolate = [&]
+        (
+            const vector& area,
+            const scalar lambda,
+
+            const Type& ownVal,
+            const Type& neiVal,
+
+            const Type& correction
+
+        ) -> DivType
+        {
+            return area & ((lambda*(ownVal - neiVal) + neiVal) + correction);
+        };
+
+        fvc::surfaceSum
+        (
+            this->tinterpScheme_().weights(vf),
+            vf,
+            faceCorr,
+            dotInterpolate,
+            div,
+            false
+        );
+    }
+    else
+    {
+        const auto dotInterpolate = [&]
+        (
+            const vector& area,
+            const scalar lambda,
+            const Type& ownVal,
+            const Type& neiVal
+        ) -> DivType
+        {
+            return area & (lambda*(ownVal - neiVal) + neiVal);
+        };
+
+        fvc::surfaceSum
+        (
+            this->tinterpScheme_().weights(vf),
+            vf,
+            dotInterpolate,
+            div,
+            false
+        );
+    }
+
+    tDiv.ref().primitiveFieldRef() /= mesh.V();
+
+    tDiv.ref().correctBoundaryConditions();
+
+    return tDiv;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/fused/finiteVolume/fusedGaussDivScheme.H

@@ -0,0 +1,125 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2016 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+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/>.
+
+Class
+    Foam::fv::fusedGaussDivScheme
+
+Group
+    grpFvDivSchemes
+
+Description
+    Variant of gaussDiv that avoids intermediates
+
+SourceFiles
+    fusedGaussDivScheme.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef fusedGaussDivScheme_H
+#define fusedGaussDivScheme_H
+
+#include "divScheme.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace fv
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class fusedGaussDivScheme Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class fusedGaussDivScheme
+:
+    public fv::divScheme<Type>
+{
+    // Private Member Functions
+
+        //- No copy construct
+        fusedGaussDivScheme(const fusedGaussDivScheme&) = delete;
+
+        //- No copy assignment
+        void operator=(const fusedGaussDivScheme&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("fusedGauss");
+
+
+    // Constructors
+
+        //- Construct null
+        fusedGaussDivScheme(const fvMesh& mesh)
+        :
+            divScheme<Type>(mesh)
+        {}
+
+        //- Construct from Istream
+        fusedGaussDivScheme(const fvMesh& mesh, Istream& is)
+        :
+            divScheme<Type>(mesh, is)
+        {}
+
+
+    // Member Functions
+
+        tmp
+        <
+            GeometricField
+            <typename innerProduct<vector, Type>::type, fvPatchField, volMesh>
+        > fvcDiv
+        (
+            const GeometricField<Type, fvPatchField, volMesh>&
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "fusedGaussDivScheme.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/fused/finiteVolume/fusedGaussDivSchemes.C

@@ -0,0 +1,35 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2015 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+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 "fusedGaussDivScheme.H"
+#include "fvMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makeFvDivScheme(fusedGaussDivScheme)
+
+// ************************************************************************* //

src/fused/finiteVolume/fusedGaussGrad.C

@@ -0,0 +1,263 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2018-2021 OpenCFD Ltd.
+    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 "fusedGaussGrad.H"
+#include "extrapolatedCalculatedFvPatchField.H"
+#include "fvcSurfaceOps.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::tmp
+<
+    Foam::GeometricField
+    <
+        typename Foam::outerProduct<Foam::vector, Type>::type,
+        Foam::fvPatchField,
+        Foam::volMesh
+    >
+>
+Foam::fv::fusedGaussGrad<Type>::gradf
+(
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& ssf,
+    const word& name
+)
+{
+    typedef typename outerProduct<vector, Type>::type GradType;
+    typedef GeometricField<GradType, fvPatchField, volMesh> GradFieldType;
+
+    const fvMesh& mesh = ssf.mesh();
+
+    DebugPout<< "fusedGaussGrad<Type>::gradf on " << ssf.name()
+        << " with name " << name << endl;
+
+    tmp<GradFieldType> tgGrad
+    (
+        new GradFieldType
+        (
+            IOobject
+            (
+                name,
+                ssf.instance(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimensioned<GradType>(ssf.dimensions()/dimLength, Zero),
+            fvPatchFieldBase::extrapolatedCalculatedType()
+        )
+    );
+    GradFieldType& gGrad = tgGrad.ref();
+
+    const labelUList& owner = mesh.owner();
+    const labelUList& neighbour = mesh.neighbour();
+    const vectorField& Sf = mesh.Sf();
+
+    Field<GradType>& igGrad = gGrad;
+    const Field<Type>& issf = ssf;
+
+    forAll(owner, facei)
+    {
+        const GradType Sfssf = Sf[facei]*issf[facei];
+
+        igGrad[owner[facei]] += Sfssf;
+        igGrad[neighbour[facei]] -= Sfssf;
+    }
+
+    forAll(mesh.boundary(), patchi)
+    {
+        const labelUList& pFaceCells =
+            mesh.boundary()[patchi].faceCells();
+
+        const vectorField& pSf = mesh.Sf().boundaryField()[patchi];
+
+        const fvsPatchField<Type>& pssf = ssf.boundaryField()[patchi];
+
+        forAll(mesh.boundary()[patchi], facei)
+        {
+            igGrad[pFaceCells[facei]] += pSf[facei]*pssf[facei];
+        }
+    }
+
+    igGrad /= mesh.V();
+
+    gGrad.correctBoundaryConditions();
+
+    return tgGrad;
+}
+
+
+template<class Type>
+Foam::tmp
+<
+    Foam::GeometricField
+    <
+        typename Foam::outerProduct<Foam::vector, Type>::type,
+        Foam::fvPatchField,
+        Foam::volMesh
+    >
+>
+Foam::fv::fusedGaussGrad<Type>::calcGrad
+(
+    const GeometricField<Type, fvPatchField, volMesh>& vf,
+    const word& name
+) const
+{
+    typedef typename outerProduct<vector, Type>::type GradType;
+    typedef GeometricField<GradType, fvPatchField, volMesh> GradFieldType;
+
+    const fvMesh& mesh = vf.mesh();
+
+    DebugPout<< "fusedGaussGrad<Type>::calcGrad on " << vf.name()
+        << " with name " << name << endl;
+
+    tmp<GradFieldType> tgGrad
+    (
+        new GradFieldType
+        (
+            IOobject
+            (
+                name,
+                vf.instance(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimensioned<GradType>(vf.dimensions()/dimLength, Zero),
+            fvPatchFieldBase::extrapolatedCalculatedType()
+        )
+    );
+    GradFieldType& gGrad = tgGrad.ref();
+
+    if (this->tinterpScheme_().corrected())
+    {
+        const auto tfaceCorr(this->tinterpScheme_().correction(vf));
+        auto& faceCorr = tfaceCorr();
+
+        DebugPout<< "fusedGaussGrad<Type>::calcGrad corrected interpScheme "
+            << this->tinterpScheme_().type() << endl;
+
+        const auto interpolate = [&]
+        (
+            const vector& area,
+            const scalar lambda,
+
+            const Type& ownVal,
+            const Type& neiVal,
+
+            const Type& correction
+
+        ) -> GradType
+        {
+            return area*((lambda*(ownVal - neiVal) + neiVal) + correction);
+        };
+
+        fvc::surfaceSum
+        (
+            this->tinterpScheme_().weights(vf),
+            vf,
+            faceCorr,
+            interpolate,
+            gGrad,
+            false
+        );
+    }
+    else
+    {
+        DebugPout<< "fusedGaussGrad<Type>::calcGrad uncorrected interpScheme "
+            << this->tinterpScheme_().type() << endl;
+
+        const auto interpolate = [&]
+        (
+            const vector& area,
+            const scalar lambda,
+            const Type& ownVal,
+            const Type& neiVal
+        ) -> GradType
+        {
+            return area*(lambda*(ownVal - neiVal) + neiVal);
+        };
+
+        fvc::surfaceSum
+        (
+            tinterpScheme_().weights(vf),
+            vf,
+            interpolate,
+            gGrad,
+            false
+        );
+    }
+
+    gGrad.primitiveFieldRef() /= mesh.V();
+
+    gGrad.correctBoundaryConditions();
+
+    correctBoundaryConditions(vf, gGrad);
+
+    return tgGrad;
+}
+
+
+template<class Type>
+template<class GradType>
+void Foam::fv::fusedGaussGrad<Type>::correctBoundaryConditions
+(
+    const GeometricField<Type, fvPatchField, volMesh>& vf,
+    GeometricField<GradType, fvPatchField, volMesh>& gGrad
+)
+{
+    DebugPout<< "fusedGaussGrad<Type>::correctBoundaryConditions on "
+        << vf.name() << " with gGrad " << gGrad.name() << endl;
+
+    const fvMesh& mesh = vf.mesh();
+    auto& gGradbf = gGrad.boundaryFieldRef();
+
+    forAll(vf.boundaryField(), patchi)
+    {
+        if (!vf.boundaryField()[patchi].coupled())
+        {
+            const auto& pSf = mesh.Sf().boundaryField()[patchi];
+            const auto tsnGrad(vf.boundaryField()[patchi].snGrad());
+            const auto& snGrad = tsnGrad();
+            auto& pgrad = gGradbf[patchi];
+
+            forAll(pgrad, facei)
+            {
+                const vector n(pSf[facei]/mag(pSf[facei]));
+                const Type uncorrectSnGrad(n & pgrad[facei]);
+                pgrad[facei] += n*(snGrad[facei] - uncorrectSnGrad);
+            }
+        }
+    }
+}
+
+
+// ************************************************************************* //

src/fused/finiteVolume/fusedGaussGrad.H

@@ -0,0 +1,177 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2018-2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::fv::fusedGaussGrad
+
+Group
+    grpFvGradSchemes
+
+Description
+    Basic second-order gradient scheme using face-interpolation
+    and Gauss' theorem.
+
+SourceFiles
+    fusedGaussGrad.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_fusedGaussGrad_H
+#define Foam_fusedGaussGrad_H
+
+#include "gradScheme.H"
+#include "surfaceInterpolationScheme.H"
+#include "linear.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace fv
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class fusedGaussGrad Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class fusedGaussGrad
+:
+    public fv::gradScheme<Type>
+{
+    // Private Data
+
+        //- Interpolation scheme
+        tmp<surfaceInterpolationScheme<Type>> tinterpScheme_;
+
+
+    // Private Member Functions
+
+        //- No copy construct
+        fusedGaussGrad(const fusedGaussGrad&) = delete;
+
+        //- No copy assignment
+        void operator=(const fusedGaussGrad&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("fusedGauss");
+
+
+    // Constructors
+
+        //- Construct from mesh
+        fusedGaussGrad(const fvMesh& mesh)
+        :
+            gradScheme<Type>(mesh),
+            tinterpScheme_(new linear<Type>(mesh))
+        {}
+
+        //- Construct from mesh and Istream
+        fusedGaussGrad(const fvMesh& mesh, Istream& is)
+        :
+            gradScheme<Type>(mesh),
+            tinterpScheme_(nullptr)
+        {
+            if (is.eof())
+            {
+                tinterpScheme_.reset
+                (
+                    new linear<Type>(mesh)
+                );
+            }
+            else
+            {
+                tinterpScheme_.reset
+                (
+                    surfaceInterpolationScheme<Type>::New(mesh, is)
+                );
+            }
+        }
+
+
+    // Member Functions
+
+        //- Return the gradient of the given field
+        //- calculated using Gauss' theorem on the given surface field
+        static
+        tmp
+        <
+            GeometricField
+            <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>
+        > gradf
+        (
+            const GeometricField<Type, fvsPatchField, surfaceMesh>&,
+            const word& name
+        );
+
+        //- Return the gradient of the given field to the gradScheme::grad
+        //- for optional caching
+        virtual tmp
+        <
+            GeometricField
+            <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>
+        > calcGrad
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& vsf,
+            const word& name
+        ) const;
+
+        //- Correct the boundary values of the gradient using the patchField
+        //- snGrad functions
+        template<class GradType>
+        static void correctBoundaryConditions
+        (
+            const GeometricField<Type, fvPatchField, volMesh>&,
+            GeometricField<GradType, fvPatchField, volMesh>&
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "fusedGaussGrad.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/fused/finiteVolume/fusedGaussGrads.C

@@ -0,0 +1,35 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2015 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+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 "fvMesh.H"
+#include "fusedGaussGrad.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makeFvGradScheme(fusedGaussGrad)
+
+// ************************************************************************* //

src/fused/finiteVolume/fusedGaussLaplacianScheme.C

@@ -0,0 +1,581 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2016 OpenFOAM Foundation
+    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 "fusedGaussLaplacianScheme.H"
+#include "fvcSurfaceOps.H"
+#include "surfaceInterpolate.H"
+#include "fvcDiv.H"
+#include "fvcGrad.H"
+#include "fvMatrices.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace fv
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type, class GType>
+tmp<fvMatrix<Type>>
+fusedGaussLaplacianScheme<Type, GType>::fvmLaplacianUncorrected
+(
+    const surfaceScalarField& gammaMagSf,
+    const surfaceScalarField& deltaCoeffs,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    DebugPout
+        << "fusedGaussLaplacianScheme<Type, GType>::fvmLaplacianUncorrected on "
+        << vf.name()
+        << " with gammaMagSf " << gammaMagSf.name()
+        << " with deltaCoeffs " << deltaCoeffs.name()
+        << endl;
+
+    tmp<fvMatrix<Type>> tfvm
+    (
+        new fvMatrix<Type>
+        (
+            vf,
+            deltaCoeffs.dimensions()*gammaMagSf.dimensions()*vf.dimensions()
+        )
+    );
+    fvMatrix<Type>& fvm = tfvm.ref();
+
+    //fvm.upper() = deltaCoeffs.primitiveField()*gammaMagSf.primitiveField();
+    multiply
+    (
+        fvm.upper(),
+        deltaCoeffs.primitiveField(),
+        gammaMagSf.primitiveField()
+    );
+    fvm.negSumDiag();
+
+    forAll(vf.boundaryField(), patchi)
+    {
+        const fvPatchField<Type>& pvf = vf.boundaryField()[patchi];
+        const fvsPatchScalarField& pGamma = gammaMagSf.boundaryField()[patchi];
+        const fvsPatchScalarField& pDeltaCoeffs =
+            deltaCoeffs.boundaryField()[patchi];
+
+        auto& intCoeffs = fvm.internalCoeffs()[patchi];
+        auto& bouCoeffs = fvm.boundaryCoeffs()[patchi];
+
+        if (pvf.coupled())
+        {
+            //intCoeffs = pGamma*pvf.gradientInternalCoeffs(pDeltaCoeffs);
+            multiply
+            (
+                intCoeffs,
+                pGamma,
+                pvf.gradientInternalCoeffs(pDeltaCoeffs)()
+            );
+            //bouCoeffs = -pGamma*pvf.gradientBoundaryCoeffs(pDeltaCoeffs);
+            multiply
+            (
+                bouCoeffs,
+                pGamma,
+                pvf.gradientBoundaryCoeffs(pDeltaCoeffs)()
+            );
+            bouCoeffs.negate();
+        }
+        else
+        {
+            //intCoeffs = pGamma*pvf.gradientInternalCoeffs();
+            multiply
+            (
+                intCoeffs,
+                pGamma,
+                pvf.gradientInternalCoeffs()()
+            );
+            //bouCoeffs = -pGamma*pvf.gradientBoundaryCoeffs();
+            multiply
+            (
+                bouCoeffs,
+                pGamma,
+                pvf.gradientBoundaryCoeffs()()
+            );
+            bouCoeffs.negate();
+        }
+    }
+
+    return tfvm;
+}
+
+
+template<class Type, class GType>
+tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
+fusedGaussLaplacianScheme<Type, GType>::gammaSnGradCorr
+(
+    const surfaceVectorField& SfGammaCorr,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    const fvMesh& mesh = this->mesh();
+
+    DebugPout<< "fusedGaussLaplacianScheme<Type, GType>::gammaSnGradCorr on "
+        << vf.name() << " with SfGammCorr " << SfGammaCorr.name() << endl;
+
+    tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> tgammaSnGradCorr
+    (
+        new GeometricField<Type, fvsPatchField, surfaceMesh>
+        (
+            IOobject
+            (
+                "gammaSnGradCorr("+vf.name()+')',
+                vf.instance(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            SfGammaCorr.dimensions()
+           *vf.dimensions()*mesh.deltaCoeffs().dimensions()
+        )
+    );
+    auto& gammaSnGradCorr = tgammaSnGradCorr.ref();
+    gammaSnGradCorr.oriented() = SfGammaCorr.oriented();
+
+    for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
+    {
+        gammaSnGradCorr.replace
+        (
+            cmpt,
+            fvc::dotInterpolate(SfGammaCorr, fvc::grad(vf.component(cmpt)))
+        );
+    }
+
+    return tgammaSnGradCorr;
+}
+
+
+/*
+template<class Type, class GType>
+void fusedGaussLaplacianScheme<Type, GType>::gradComponent
+(
+    const surfaceScalarField& weights,
+    const GeometricField<Type, fvPatchField, volMesh>& vf,
+    const direction cmpt,
+    GeometricField<Type, fvPatchField, volMesh>& gGrad
+)
+{
+    gGrad = dimensioned<Type>(vf.dimensions()/dimLength, Zero);
+
+    // Calculate grad of vf.component(cmpt)
+    fvc::GaussOp
+    (
+        vf,
+        weights,
+        componentInterpolate<Type>(cmpt),
+        gGrad
+    );
+}
+template<class Type, class GType>
+tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
+fusedGaussLaplacianScheme<Type, GType>::gammaSnGradCorr
+(
+    const surfaceScalarField& weights,
+    const GeometricField<GType, fvPatchField, volMesh>& gamma,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    const fvMesh& mesh = this->mesh();
+
+    tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> tgammaSnGradCorr
+    (
+        new GeometricField<Type, fvsPatchField, surfaceMesh>
+        (
+            IOobject
+            (
+                "gammaSnGradCorr("+vf.name()+')',
+                vf.instance(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            gamma.dimensions()
+           *vf.dimensions()*mesh.deltaCoeffs().dimensions()
+        )
+    );
+    auto& gammaSnGradCorr = tgammaSnGradCorr.ref();
+
+    gammaSnGradCorr.oriented() = gamma.oriented();
+
+
+    GeometricField<Type, fvPatchField, volMesh> gradCmptFld
+    (
+        IOobject
+        (
+            vf.name() + ".component()",
+            vf.instance(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        mesh,
+        vf.dimensions()
+    );
+
+    for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
+    {
+        // Calculate fvc::grad(vf.component(cmpt)) into gradCmptFld
+        gradComponent
+        (
+            weights,
+            vf,
+            cmpt,
+
+            gradCmptFld
+        );
+
+        //gammaSnGradCorr.replace
+        //(
+        //    cmpt,
+        //    fvc::dotInterpolate(SfGammaCorr, gradCmptFld)
+        //);
+
+
+        fvc::interpolate
+        (
+            weights,
+
+            gradCmptFld,        // fvc::grad(vf.component(cmpt))
+
+            gamma,              // weight field
+
+            tanInterpolate<Type, GType>(cmpt),
+
+            gammaSnGradCorr
+        );
+    }
+
+    return tgammaSnGradCorr;
+}
+*/
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type, class GType>
+tmp<GeometricField<Type, fvPatchField, volMesh>>
+fusedGaussLaplacianScheme<Type, GType>::fvcLaplacian
+(
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> FieldType;
+    //typedef GeometricField<Type, fvsPatchField, surfaceMesh> SurfaceFieldType;
+
+    tmp<FieldType> tresult
+    (
+        new FieldType
+        (
+            IOobject
+            (
+                "laplacian(" + vf.name() + ')',
+                vf.instance(),
+                vf.mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            vf.mesh(),
+            dimensioned<Type>(vf.dimensions()/dimArea, Zero),
+            fvPatchFieldBase::extrapolatedCalculatedType()
+        )
+    );
+    FieldType& result = tresult.ref();
+
+    DebugPout
+        << "fusedGaussLaplacianScheme<Type, GType>::fvcLaplacian on "
+        << vf.name()
+        << " to generate " << result.name() << endl;
+
+    // Note: cannot use fvc::GaussOp since specialised handling on boundary.
+    // Maybe bypass for processor boundaries?
+
+    const auto tdeltaCoeffs(this->tsnGradScheme_().deltaCoeffs(vf));
+    const auto& deltaCoeffs = tdeltaCoeffs();
+
+
+    const fvMesh& mesh = vf.mesh();
+
+    if (this->tsnGradScheme_().corrected())
+    {
+        // Problem when instantiating for tensors - outerProduct not defined.
+        // scalar/vector specialisations in *Schemes.C file.
+
+        FatalErrorInFunction<< "Corrected snGrad not supported for field "
+            << vf.name() << exit(FatalError);
+
+        //typedef typename outerProduct<vector, Type>::type GradType;
+        //typedef GeometricField<GradType, fvPatchField, volMesh> GradFieldType;
+        //// Calculate sn gradient
+        //tmp<SurfaceFieldType> tfaceGrad
+        //(
+        //    new SurfaceFieldType
+        //    (
+        //        IOobject
+        //        (
+        //            "snGradCorr("+vf.name()+')',
+        //            vf.instance(),
+        //            mesh,
+        //            IOobject::NO_READ,
+        //            IOobject::NO_WRITE
+        //        ),
+        //        mesh,
+        //        vf.dimensions()
+        //    )
+        //);
+        //
+        //{
+        //    // Calculate gradient
+        //    tmp<GradFieldType> tgGrad
+        //    (
+        //        gradScheme<Type>::New
+        //        (
+        //            mesh,
+        //            mesh.gradScheme("grad(" + vf.name() + ')')
+        //        )().grad(vf, "grad(" + vf.name() + ')')
+        //    );
+        //    const auto& gGrad = tgGrad();
+        //
+        //    // Doing a dotinterpolate with nonOrthCorrectionVectors
+        //    const auto dotInterpolate = [&]
+        //    (
+        //        const vector& area,
+        //        const scalar lambda,
+        //
+        //        const GradType& ownVal,
+        //        const GradType& neiVal,
+        //
+        //        const vector& dotVector,
+        //
+        //        Type& result
+        //    )
+        //    {
+        //        result = dotVector&(lambda*(ownVal - neiVal) + neiVal);
+        //    };
+        //
+        //    fvc::interpolate
+        //    (
+        //        mesh.surfaceInterpolation::weights(),   // linear interp
+        //        gGrad,                          // volume field
+        //        mesh.nonOrthCorrectionVectors(),// surface multiplier
+        //        dotInterpolate,
+        //        tfaceGrad.ref()
+        //    );
+        //}
+        //const auto& faceGrad = tfaceGrad();
+        //
+        //const auto snGrad = [&]
+        //(
+        //    const vector& Sf,
+        //    const scalar dc,
+        //    const Type& ownVal,
+        //    const Type& neiVal,
+        //    const Type& correction
+        //) -> Type
+        //{
+        //    const auto snGrad(dc*(neiVal-ownVal) + correction);
+        //    return mag(Sf)*snGrad;
+        //};
+        //
+        //fvc::surfaceSnSum
+        //(
+        //    deltaCoeffs,
+        //    vf,
+        //    faceGrad,   // face-based addition
+        //    snGrad,
+        //    result,
+        //    false       // avoid boundary evaluation until volume division
+        //);
+    }
+    else
+    {
+        const auto snGrad = [&]
+        (
+            const vector& Sf,
+            const scalar dc,
+            const Type& ownVal,
+            const Type& neiVal
+        ) -> Type
+        {
+            const auto snGrad(dc*(neiVal-ownVal));
+            return mag(Sf)*snGrad;
+        };
+
+        fvc::surfaceSnSum
+        (
+            deltaCoeffs,
+            vf,
+            snGrad,
+            result,
+            false       // avoid boundary evaluation until volume division
+        );
+    }
+
+    result.primitiveFieldRef() /= mesh.V();
+    result.correctBoundaryConditions();
+
+    return tresult;
+}
+
+
+template<class Type, class GType>
+tmp<GeometricField<Type, fvPatchField, volMesh>>
+fusedGaussLaplacianScheme<Type, GType>::fvcLaplacian
+(
+    const GeometricField<GType, fvPatchField, volMesh>& gamma,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    DebugPout
+        << "fusedGaussLaplacianScheme<Type, GType>::fvcLaplacian on "
+        << vf.name() << " with gamma " << gamma.name() << endl;
+
+    return fvcLaplacian(this->tinterpGammaScheme_().interpolate(gamma)(), vf);
+}
+
+
+template<class Type, class GType>
+tmp<fvMatrix<Type>>
+fusedGaussLaplacianScheme<Type, GType>::fvmLaplacian
+(
+    const GeometricField<GType, fvsPatchField, surfaceMesh>& gamma,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    DebugPout<< "fusedGaussLaplacianScheme<Type, GType>::fvmLaplacian on "
+        << vf.name() << " with gamma " << gamma.name() << endl;
+
+    const fvMesh& mesh = this->mesh();
+
+    const surfaceVectorField Sn(mesh.Sf()/mesh.magSf());
+    const surfaceVectorField SfGamma(mesh.Sf() & gamma);
+    const GeometricField<scalar, fvsPatchField, surfaceMesh> SfGammaSn
+    (
+        SfGamma & Sn
+    );
+    const surfaceVectorField SfGammaCorr(SfGamma - SfGammaSn*Sn);
+
+    tmp<fvMatrix<Type>> tfvm = fvmLaplacianUncorrected
+    (
+        SfGammaSn,
+        this->tsnGradScheme_().deltaCoeffs(vf),
+        vf
+    );
+    fvMatrix<Type>& fvm = tfvm.ref();
+
+    tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> tfaceFluxCorrection
+        = gammaSnGradCorr(SfGammaCorr, vf);
+
+    if (this->tsnGradScheme_().corrected())
+    {
+        tfaceFluxCorrection.ref() +=
+            SfGammaSn*this->tsnGradScheme_().correction(vf);
+    }
+
+    fvm.source() -= mesh.V()*fvc::div(tfaceFluxCorrection())().primitiveField();
+
+    if (mesh.fluxRequired(vf.name()))
+    {
+        fvm.faceFluxCorrectionPtr(tfaceFluxCorrection.ptr());
+    }
+
+    return tfvm;
+}
+
+
+template<class Type, class GType>
+tmp<fvMatrix<Type>>
+fusedGaussLaplacianScheme<Type, GType>::fvmLaplacian
+(
+    const GeometricField<GType, fvPatchField, volMesh>& gamma,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    DebugPout<< "fusedGaussLaplacianScheme<Type, GType>::fvmLaplacian on "
+        << vf.name() << " with gamma " << gamma.name() << endl;
+    return fvmLaplacian(this->tinterpGammaScheme_().interpolate(gamma)(), vf);
+}
+
+
+template<class Type, class GType>
+tmp<GeometricField<Type, fvPatchField, volMesh>>
+fusedGaussLaplacianScheme<Type, GType>::fvcLaplacian
+(
+    const GeometricField<GType, fvsPatchField, surfaceMesh>& gamma,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    DebugPout<< "fusedGaussLaplacianScheme<Type, GType>::fvcLaplacian on "
+        << vf.name() << " with gamma " << gamma.name() << endl;
+
+    const fvMesh& mesh = this->mesh();
+
+    const surfaceVectorField Sn(mesh.Sf()/mesh.magSf());
+    const surfaceVectorField SfGamma(mesh.Sf() & gamma);
+    const GeometricField<scalar, fvsPatchField, surfaceMesh> SfGammaSn
+    (
+        SfGamma & Sn
+    );
+    const surfaceVectorField SfGammaCorr(SfGamma - SfGammaSn*Sn);
+
+    tmp<GeometricField<Type, fvPatchField, volMesh>> tLaplacian
+    (
+        fvc::div
+        (
+            SfGammaSn*this->tsnGradScheme_().snGrad(vf)
+          + gammaSnGradCorr(SfGammaCorr, vf)
+        )
+    );
+
+    tLaplacian.ref().rename
+    (
+        "laplacian(" + gamma.name() + ',' + vf.name() + ')'
+    );
+
+    return tLaplacian;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/fused/finiteVolume/fusedGaussLaplacianScheme.H

@@ -0,0 +1,264 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2019 OpenCFD Ltd.
+    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/>.
+
+Class
+    Foam::fv::GaussLaplacianScheme
+
+Group
+    grpFvLaplacianSchemes
+
+Description
+    Variant of GaussLaplacian that avoids intermediate fields
+
+SourceFiles
+    fusedGaussLaplacianScheme.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef fusedGaussLaplacianScheme_H
+#define fusedGaussLaplacianScheme_H
+
+#include "laplacianScheme.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace fv
+{
+
+/*---------------------------------------------------------------------------*\
+                  Class fusedGaussLaplacianScheme Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type, class GType>
+class fusedGaussLaplacianScheme
+:
+    public fv::laplacianScheme<Type, GType>
+{
+    // Private Member Functions
+
+        //- See gaussLaplacianScheme
+        tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> gammaSnGradCorr
+        (
+            const surfaceVectorField& SfGammaCorr,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        );
+
+        //- No copy construct
+        fusedGaussLaplacianScheme(const fusedGaussLaplacianScheme&)
+            = delete;
+
+        //- No copy assignment
+        void operator=(const fusedGaussLaplacianScheme&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("fusedGauss");
+
+
+    // Constructors
+
+        //- Construct null
+        fusedGaussLaplacianScheme(const fvMesh& mesh)
+        :
+            laplacianScheme<Type, GType>(mesh)
+        {}
+
+        //- Construct from Istream
+        fusedGaussLaplacianScheme(const fvMesh& mesh, Istream& is)
+        :
+            laplacianScheme<Type, GType>(mesh, is)
+        {}
+
+        //- Construct from mesh, interpolation and snGradScheme schemes
+        fusedGaussLaplacianScheme
+        (
+            const fvMesh& mesh,
+            const tmp<surfaceInterpolationScheme<GType>>& igs,
+            const tmp<snGradScheme<Type>>& sngs
+        )
+        :
+            laplacianScheme<Type, GType>(mesh, igs, sngs)
+        {}
+
+
+    //- Destructor
+    virtual ~fusedGaussLaplacianScheme() = default;
+
+
+    // Member Functions
+
+        static tmp<fvMatrix<Type>> fvmLaplacianUncorrected
+        (
+            const surfaceScalarField& gammaMagSf,
+            const surfaceScalarField& deltaCoeffs,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        );
+
+        ////- Helper: grad on component. Move to GaussGrad? Template on
+        ////  - cell-access operator?
+        ////  - cell-to-face operator?
+        ////  - face-to-cell operator?
+        //static void gradComponent
+        //(
+        //    const surfaceScalarField& weights,
+        //    const GeometricField<Type, fvPatchField, volMesh>& vf,
+        //    const direction cmpt,
+        //    GeometricField<Type, fvPatchField, volMesh>& gGrad
+        //);
+
+        virtual tmp<GeometricField<Type, fvPatchField, volMesh>> fvcLaplacian
+        (
+            const GeometricField<Type, fvPatchField, volMesh>&
+        );
+
+        virtual tmp<fvMatrix<Type>> fvmLaplacian
+        (
+            const GeometricField<GType, fvsPatchField, surfaceMesh>&,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        );
+
+        virtual tmp<GeometricField<Type, fvPatchField, volMesh>> fvcLaplacian
+        (
+            const GeometricField<GType, fvsPatchField, surfaceMesh>&,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        );
+
+
+    // New: override laplacianScheme::fvc/mLaplacian for scalar gamma and
+    //      scalar,vector volField
+
+        virtual tmp<fvMatrix<Type>> fvmLaplacian
+        (
+            const GeometricField<GType, fvPatchField, volMesh>&,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        );
+
+        virtual tmp<GeometricField<Type, fvPatchField, volMesh>> fvcLaplacian
+        (
+            const GeometricField<GType, fvPatchField, volMesh>&,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        );
+};
+
+
+// Use macros to emulate partial-specialisation of the Laplacian functions
+// for scalar diffusivity gamma
+
+#define defineFvmLaplacianScalarGamma(Type)                                    \
+                                                                               \
+template<>                                                                     \
+tmp<fvMatrix<Type>> fusedGaussLaplacianScheme<Type, scalar>::fvmLaplacian      \
+(                                                                              \
+    const GeometricField<scalar, fvsPatchField, surfaceMesh>&,                 \
+    const GeometricField<Type, fvPatchField, volMesh>&                         \
+);                                                                             \
+                                                                               \
+template<>                                                                     \
+tmp<GeometricField<Type, fvPatchField, volMesh>>                               \
+fusedGaussLaplacianScheme<Type, scalar>::fvcLaplacian                          \
+(                                                                              \
+    const GeometricField<scalar, fvsPatchField, surfaceMesh>&,                 \
+    const GeometricField<Type, fvPatchField, volMesh>&                         \
+);
+
+defineFvmLaplacianScalarGamma(scalar);
+defineFvmLaplacianScalarGamma(vector);
+defineFvmLaplacianScalarGamma(sphericalTensor);
+defineFvmLaplacianScalarGamma(symmTensor);
+defineFvmLaplacianScalarGamma(tensor);
+
+// Unweighted laplacian
+template<>
+tmp<GeometricField<scalar, fvPatchField, volMesh>>
+fusedGaussLaplacianScheme<scalar, scalar>::fvcLaplacian
+(
+    const GeometricField<scalar, fvPatchField, volMesh>&
+);
+template<>
+tmp<GeometricField<vector, fvPatchField, volMesh>>
+fusedGaussLaplacianScheme<vector, scalar>::fvcLaplacian
+(
+    const GeometricField<vector, fvPatchField, volMesh>&
+);
+
+// Weighted laplacian
+template<>
+tmp<GeometricField<scalar, fvPatchField, volMesh>>
+fusedGaussLaplacianScheme<scalar, scalar>::fvcLaplacian
+(
+    const GeometricField<scalar, fvPatchField, volMesh>&,
+    const GeometricField<scalar, fvPatchField, volMesh>&
+);
+template<>
+tmp<GeometricField<vector, fvPatchField, volMesh>>
+fusedGaussLaplacianScheme<vector, scalar>::fvcLaplacian
+(
+    const GeometricField<scalar, fvPatchField, volMesh>&,
+    const GeometricField<vector, fvPatchField, volMesh>&
+);
+template<>
+tmp<fvMatrix<scalar>>
+fusedGaussLaplacianScheme<scalar, scalar>::fvmLaplacian
+(
+    const GeometricField<scalar, fvPatchField, volMesh>&,
+    const GeometricField<scalar, fvPatchField, volMesh>&
+);
+template<>
+tmp<fvMatrix<vector>>
+fusedGaussLaplacianScheme<vector, scalar>::fvmLaplacian
+(
+    const GeometricField<scalar, fvPatchField, volMesh>&,
+    const GeometricField<vector, fvPatchField, volMesh>&
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "fusedGaussLaplacianScheme.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/fused/finiteVolume/fusedGaussLaplacianSchemes.C

@@ -0,0 +1,866 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2016 OpenFOAM Foundation
+    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 "fusedGaussLaplacianScheme.H"
+#include "fvMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makeFvLaplacianScheme(fusedGaussLaplacianScheme)
+
+#define declareFvmLaplacianScalarGamma(Type)                                   \
+                                                                               \
+template<>                                                                     \
+Foam::tmp<Foam::fvMatrix<Foam::Type>>                                          \
+Foam::fv::fusedGaussLaplacianScheme<Foam::Type, Foam::scalar>::                \
+fvmLaplacian                                                                   \
+(                                                                              \
+    const GeometricField<scalar, fvsPatchField, surfaceMesh>& gamma,           \
+    const GeometricField<Type, fvPatchField, volMesh>& vf                      \
+)                                                                              \
+{                                                                              \
+    DebugPout<< "fusedGaussLaplacianScheme::fvmLaplacian on " << vf.name()     \
+        << " with scalar gamma " << gamma.name() << endl;                      \
+                                                                               \
+    const fvMesh& mesh = this->mesh();                                         \
+                                                                               \
+    GeometricField<scalar, fvsPatchField, surfaceMesh> gammaMagSf              \
+    (                                                                          \
+        gamma*mesh.magSf()                                                     \
+    );                                                                         \
+                                                                               \
+    tmp<fvMatrix<Type>> tfvm = fvmLaplacianUncorrected                         \
+    (                                                                          \
+        gammaMagSf,                                                            \
+        this->tsnGradScheme_().deltaCoeffs(vf),                                \
+        vf                                                                     \
+    );                                                                         \
+    fvMatrix<Type>& fvm = tfvm.ref();                                          \
+                                                                               \
+    if (this->tsnGradScheme_().corrected())                                    \
+    {                                                                          \
+        if (mesh.fluxRequired(vf.name()))                                      \
+        {                                                                      \
+            fvm.faceFluxCorrectionPtr() = std::make_unique                     \
+            <                                                                  \
+                GeometricField<Type, fvsPatchField, surfaceMesh>               \
+            >                                                                  \
+            (                                                                  \
+                gammaMagSf*this->tsnGradScheme_().correction(vf)               \
+            );                                                                 \
+                                                                               \
+            fvm.source() -=                                                    \
+                mesh.V()*                                                      \
+                fvc::div                                                       \
+                (                                                              \
+                    *fvm.faceFluxCorrectionPtr()                               \
+                )().primitiveField();                                          \
+        }                                                                      \
+        else                                                                   \
+        {                                                                      \
+            fvm.source() -=                                                    \
+                mesh.V()*                                                      \
+                fvc::div                                                       \
+                (                                                              \
+                    gammaMagSf*this->tsnGradScheme_().correction(vf)           \
+                )().primitiveField();                                          \
+        }                                                                      \
+    }                                                                          \
+                                                                               \
+    return tfvm;                                                               \
+}                                                                              \
+                                                                               \
+                                                                               \
+template<>                                                                     \
+Foam::tmp<Foam::GeometricField<Foam::Type, Foam::fvPatchField, Foam::volMesh>> \
+Foam::fv::fusedGaussLaplacianScheme<Foam::Type, Foam::scalar>::                \
+fvcLaplacian                                                                   \
+(                                                                              \
+    const GeometricField<scalar, fvsPatchField, surfaceMesh>& gamma,           \
+    const GeometricField<Type, fvPatchField, volMesh>& vf                      \
+)                                                                              \
+{                                                                              \
+    DebugPout<< "fvcLaplacian on " << vf.name()                                \
+        << " with scalar gamma " << gamma.name() << endl;                      \
+                                                                               \
+    const fvMesh& mesh = this->mesh();                                         \
+                                                                               \
+    tmp<GeometricField<Type, fvPatchField, volMesh>> tLaplacian                \
+    (                                                                          \
+        fvc::div(gamma*this->tsnGradScheme_().snGrad(vf)*mesh.magSf())         \
+    );                                                                         \
+                                                                               \
+    tLaplacian.ref().rename                                                    \
+    (                                                                          \
+        "laplacian(" + gamma.name() + ',' + vf.name() + ')'                    \
+    );                                                                         \
+                                                                               \
+    return tLaplacian;                                                         \
+}
+
+
+template<>
+Foam::tmp
+<
+    Foam::GeometricField<Foam::scalar, Foam::fvPatchField, Foam::volMesh>
+>
+Foam::fv::fusedGaussLaplacianScheme<Foam::scalar, Foam::scalar>::fvcLaplacian
+(
+    const GeometricField<scalar, fvPatchField, volMesh>& gamma,
+    const GeometricField<scalar, fvPatchField, volMesh>& vf
+)
+{
+    typedef scalar Type;
+    typedef GeometricField<Type, fvPatchField, volMesh> FieldType;
+    typedef GeometricField<Type, fvsPatchField, surfaceMesh> SurfaceFieldType;
+    typedef typename outerProduct<vector, Type>::type GradType;
+    typedef GeometricField<GradType, fvPatchField, volMesh> GradFieldType;
+
+    DebugPout
+        << "fusedGaussLaplacianScheme<scalar, scalar>::fvcLaplacian"
+        << " on " << vf.name() << " with gamma " << gamma.name() << endl;
+
+    const fvMesh& mesh = vf.mesh();
+
+    tmp<FieldType> tresult
+    (
+        new FieldType
+        (
+            IOobject
+            (
+                "laplacian(" + vf.name() + ')',
+                vf.instance(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimensioned<Type>
+            (
+                gamma.dimensions()*vf.dimensions()/dimArea, Zero
+            ),
+            fvPatchFieldBase::extrapolatedCalculatedType()
+        )
+    );
+    FieldType& result = tresult.ref();
+
+    const auto tweights(this->tinterpGammaScheme_().weights(gamma));
+    const auto& weights = tweights();
+    const auto tdeltaCoeffs(this->tsnGradScheme_().deltaCoeffs(vf));
+    const auto& deltaCoeffs = tdeltaCoeffs();
+
+    if (this->tsnGradScheme_().corrected())
+    {
+        // Calculate sn gradient
+        tmp<SurfaceFieldType> tfaceGrad
+        (
+            new SurfaceFieldType
+            (
+                IOobject
+                (
+                    "snGradCorr("+vf.name()+')',
+                    vf.instance(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                mesh,
+                vf.dimensions()
+            )
+        );
+
+        {
+            // Calculate gradient
+            tmp<GradFieldType> tgGrad
+            (
+                gradScheme<Type>::New
+                (
+                    mesh,
+                    mesh.gradScheme("grad(" + vf.name() + ')')
+                )().grad(vf, "grad(" + vf.name() + ')')
+            );
+            const auto& gGrad = tgGrad();
+
+            // Doing a dotinterpolate with nonOrthCorrectionVectors
+            const auto dotInterpolate = [&]
+            (
+                const vector& area,
+                const scalar lambda,
+
+                const GradType& ownVal,
+                const GradType& neiVal,
+
+                const vector& dotVector,
+
+                Type& result
+            )
+            {
+                result = dotVector&(lambda*(ownVal - neiVal) + neiVal);
+            };
+
+            fvc::interpolate
+            (
+                mesh.surfaceInterpolation::weights(),   // linear interpolation
+                gGrad,                          // volume field
+                mesh.nonOrthCorrectionVectors(),// surface multiplier
+                dotInterpolate,
+                tfaceGrad.ref()
+            );
+        }
+        const auto& faceGrad = tfaceGrad();
+
+        const auto snGrad = [&]
+        (
+            const vector& Sf,
+
+            const scalar weight,
+            const scalar ownGamma,
+            const scalar neiGamma,
+
+            const scalar dc,
+            const Type& ownVal,
+            const Type& neiVal,
+
+            const Type& correction
+        ) -> Type
+        {
+            const auto snGrad(dc*(neiVal-ownVal) + correction);
+            const scalar faceGamma(weight*(ownGamma-neiGamma)+neiGamma);
+            return mag(Sf)*faceGamma*snGrad;
+        };
+
+        fvc::surfaceSnSum
+        (
+            weights,    // gamma weights
+            gamma,
+
+            deltaCoeffs,
+            vf,
+
+            faceGrad,   // face-based addition
+
+            snGrad,
+
+            result,
+            false       // avoid boundary evaluation until volume division
+        );
+    }
+    else
+    {
+        const auto snGrad = [&]
+        (
+            const vector& Sf,
+
+            const scalar weight,
+            const scalar ownGamma,
+            const scalar neiGamma,
+
+            const scalar dc,
+            const Type& ownVal,
+            const Type& neiVal
+        ) -> Type
+        {
+            const auto snGrad(dc*(neiVal-ownVal));
+            const scalar faceGamma(weight*(ownGamma-neiGamma)+neiGamma);
+            return mag(Sf)*faceGamma*snGrad;
+        };
+
+        fvc::surfaceSnSum
+        (
+            weights,
+            gamma,
+
+            deltaCoeffs,
+            vf,
+
+            snGrad,
+
+            result,
+            false       // avoid boundary evaluation until volume division
+        );
+    }
+
+    result.primitiveFieldRef() /= mesh.V();
+    result.correctBoundaryConditions();
+
+    return tresult;
+}
+
+
+template<>
+Foam::tmp
+<
+    Foam::GeometricField<Foam::vector, Foam::fvPatchField, Foam::volMesh>
+>
+Foam::fv::fusedGaussLaplacianScheme<Foam::vector, Foam::scalar>::fvcLaplacian
+(
+    const GeometricField<scalar, fvPatchField, volMesh>& gamma,
+    const GeometricField<vector, fvPatchField, volMesh>& vf
+)
+{
+    DebugPout
+        << "fusedGaussLaplacianScheme<vector, scalar>::fvcLaplacian"
+        << " on " << vf.name() << " with gamma " << gamma.name() << endl;
+
+    typedef vector Type;
+    typedef GeometricField<Type, fvPatchField, volMesh> FieldType;
+    typedef GeometricField<Type, fvsPatchField, surfaceMesh> SurfaceFieldType;
+    typedef typename outerProduct<vector, Type>::type GradType;
+    typedef GeometricField<GradType, fvPatchField, volMesh> GradFieldType;
+
+    const fvMesh& mesh = vf.mesh();
+
+    tmp<FieldType> tresult
+    (
+        new FieldType
+        (
+            IOobject
+            (
+                "laplacian(" + vf.name() + ')',
+                vf.instance(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimensioned<Type>
+            (
+                gamma.dimensions()*vf.dimensions()/dimArea, Zero
+            ),
+            fvPatchFieldBase::extrapolatedCalculatedType()
+        )
+    );
+    FieldType& result = tresult.ref();
+
+    const auto tweights(this->tinterpGammaScheme_().weights(gamma));
+    const auto& weights = tweights();
+    const auto tdeltaCoeffs(this->tsnGradScheme_().deltaCoeffs(vf));
+    const auto& deltaCoeffs = tdeltaCoeffs();
+
+    if (this->tsnGradScheme_().corrected())
+    {
+        // Calculate sn gradient
+        tmp<SurfaceFieldType> tfaceGrad
+        (
+            new SurfaceFieldType
+            (
+                IOobject
+                (
+                    "snGradCorr("+vf.name()+')',
+                    vf.instance(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                mesh,
+                vf.dimensions()
+            )
+        );
+
+        {
+            // Calculate gradient
+            tmp<GradFieldType> tgGrad
+            (
+                gradScheme<Type>::New
+                (
+                    mesh,
+                    mesh.gradScheme("grad(" + vf.name() + ')')
+                )().grad(vf, "grad(" + vf.name() + ')')
+            );
+            const auto& gGrad = tgGrad();
+
+            // Doing a dotinterpolate with nonOrthCorrectionVectors
+            const auto dotInterpolate = [&]
+            (
+                const vector& area,
+                const scalar lambda,
+
+                const GradType& ownVal,
+                const GradType& neiVal,
+
+                const vector& dotVector,
+
+                Type& result
+            )
+            {
+                result = dotVector&(lambda*(ownVal - neiVal) + neiVal);
+            };
+
+            fvc::interpolate
+            (
+                mesh.surfaceInterpolation::weights(),   // linear interpolation
+                gGrad,                          // volume field
+                mesh.nonOrthCorrectionVectors(),// surface multiplier
+                dotInterpolate,
+                tfaceGrad.ref()
+            );
+        }
+        const auto& faceGrad = tfaceGrad();
+
+        const auto snGrad = [&]
+        (
+            const vector& Sf,
+
+            const scalar weight,
+            const scalar ownGamma,
+            const scalar neiGamma,
+
+            const scalar dc,
+            const Type& ownVal,
+            const Type& neiVal,
+
+            const Type& correction
+        ) -> Type
+        {
+            const auto snGrad(dc*(neiVal-ownVal) + correction);
+            const scalar faceGamma(weight*(ownGamma-neiGamma)+neiGamma);
+            return mag(Sf)*faceGamma*snGrad;
+        };
+
+        fvc::surfaceSnSum
+        (
+            weights,    // gamma weights
+            gamma,
+
+            deltaCoeffs,
+            vf,
+
+            faceGrad,   // face-based addition
+
+            snGrad,
+
+            result,
+            false       // avoid boundary evaluation until volume division
+        );
+    }
+    else
+    {
+        const auto snGrad = [&]
+        (
+            const vector& Sf,
+
+            const scalar weight,
+            const scalar ownGamma,
+            const scalar neiGamma,
+
+            const scalar dc,
+            const Type& ownVal,
+            const Type& neiVal
+        ) -> Type
+        {
+            const auto snGrad(dc*(neiVal-ownVal));
+            const scalar faceGamma(weight*(ownGamma-neiGamma)+neiGamma);
+            return mag(Sf)*faceGamma*snGrad;
+        };
+
+        fvc::surfaceSnSum
+        (
+            weights,    // gamma weights
+            gamma,
+
+            deltaCoeffs,
+            vf,
+
+            snGrad,
+
+            result,
+            false       // avoid boundary evaluation until volume division
+        );
+    }
+
+    result.primitiveFieldRef() /= mesh.V();
+    result.correctBoundaryConditions();
+
+    return tresult;
+}
+
+
+template<>
+Foam::tmp<Foam::fvMatrix<Foam::scalar>>
+Foam::fv::fusedGaussLaplacianScheme<Foam::scalar, Foam::scalar>::fvmLaplacian
+(
+    const GeometricField<scalar, fvPatchField, volMesh>& gamma,
+    const GeometricField<scalar, fvPatchField, volMesh>& vf
+)
+{
+    // TBD
+    DebugPout
+        << "fusedGaussLaplacianScheme<scalar, scalar>::fvmLaplacian"
+        << " on " << vf.name() << " with gamma " << gamma.name() << endl;
+    return fvmLaplacian(this->tinterpGammaScheme_().interpolate(gamma)(), vf);
+}
+
+
+template<>
+Foam::tmp<Foam::fvMatrix<Foam::vector>>
+Foam::fv::fusedGaussLaplacianScheme<Foam::vector, Foam::scalar>::fvmLaplacian
+(
+    const GeometricField<scalar, fvPatchField, volMesh>& gamma,
+    const GeometricField<vector, fvPatchField, volMesh>& vf
+)
+{
+    // TBD
+    DebugPout
+        << "fusedGaussLaplacianScheme<vector, scalar>::fvmLaplacian"
+        << " on " << vf.name() << " with gamma " << gamma.name() << endl;
+    return fvmLaplacian(this->tinterpGammaScheme_().interpolate(gamma)(), vf);
+}
+
+
+template<>
+Foam::tmp
+<
+    Foam::GeometricField<Foam::scalar, Foam::fvPatchField, Foam::volMesh>
+>
+Foam::fv::fusedGaussLaplacianScheme<Foam::scalar, Foam::scalar>::fvcLaplacian
+(
+    const GeometricField<scalar, fvPatchField, volMesh>& vf
+)
+{
+    typedef scalar Type;
+
+    typedef GeometricField<Type, fvPatchField, volMesh> FieldType;
+    typedef GeometricField<Type, fvsPatchField, surfaceMesh> SurfaceFieldType;
+    typedef typename outerProduct<vector, Type>::type GradType;
+    typedef GeometricField<GradType, fvPatchField, volMesh> GradFieldType;
+
+    const fvMesh& mesh = vf.mesh();
+
+    tmp<FieldType> tresult
+    (
+        new FieldType
+        (
+            IOobject
+            (
+                "laplacian(" + vf.name() + ')',
+                vf.instance(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimensioned<Type>(vf.dimensions()/dimArea, Zero),
+            fvPatchFieldBase::extrapolatedCalculatedType()
+        )
+    );
+    FieldType& result = tresult.ref();
+
+    DebugPout
+        << "fusedGaussLaplacianScheme<scalar, GType>::fvcLaplacian on "
+        << vf.name()
+        << " to generate " << result.name() << endl;
+
+
+    const auto tdeltaCoeffs(this->tsnGradScheme_().deltaCoeffs(vf));
+    const auto& deltaCoeffs = tdeltaCoeffs();
+
+
+    if (this->tsnGradScheme_().corrected())
+    {
+        // Calculate sn gradient
+        tmp<SurfaceFieldType> tfaceGrad
+        (
+            new SurfaceFieldType
+            (
+                IOobject
+                (
+                    "snGradCorr("+vf.name()+')',
+                    vf.instance(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                mesh,
+                vf.dimensions()
+            )
+        );
+
+        {
+            // Calculate gradient
+            tmp<GradFieldType> tgGrad
+            (
+                gradScheme<Type>::New
+                (
+                    mesh,
+                    mesh.gradScheme("grad(" + vf.name() + ')')
+                )().grad(vf, "grad(" + vf.name() + ')')
+            );
+            const auto& gGrad = tgGrad();
+
+            // Doing a dotinterpolate with nonOrthCorrectionVectors
+            const auto dotInterpolate = [&]
+            (
+                const vector& area,
+                const scalar lambda,
+
+                const GradType& ownVal,
+                const GradType& neiVal,
+
+                const vector& dotVector,
+
+                Type& result
+            )
+            {
+                result = dotVector&(lambda*(ownVal - neiVal) + neiVal);
+            };
+
+            fvc::interpolate
+            (
+                mesh.surfaceInterpolation::weights(),   // linear interpolation
+                gGrad,                          // volume field
+                mesh.nonOrthCorrectionVectors(),// surface multiplier
+                dotInterpolate,
+                tfaceGrad.ref()
+            );
+        }
+        const auto& faceGrad = tfaceGrad();
+
+        const auto snGrad = [&]
+        (
+            const vector& Sf,
+            const scalar dc,
+            const Type& ownVal,
+            const Type& neiVal,
+            const Type& correction
+        ) -> Type
+        {
+            const auto snGrad(dc*(neiVal-ownVal) + correction);
+            return mag(Sf)*snGrad;
+        };
+
+        fvc::surfaceSnSum
+        (
+            deltaCoeffs,
+            vf,
+            faceGrad,   // face-based addition
+            snGrad,
+            result,
+            false       // avoid boundary evaluation until volume division
+        );
+    }
+    else
+    {
+        const auto snGrad = [&]
+        (
+            const vector& Sf,
+            const scalar dc,
+            const Type& ownVal,
+            const Type& neiVal
+        ) -> Type
+        {
+            const auto snGrad(dc*(neiVal-ownVal));
+            return mag(Sf)*snGrad;
+        };
+
+        fvc::surfaceSnSum
+        (
+            deltaCoeffs,
+            vf,
+            snGrad,
+            result,
+            false       // avoid boundary evaluation until volume division
+        );
+    }
+
+    result.primitiveFieldRef() /= mesh.V();
+    result.correctBoundaryConditions();
+
+    return tresult;
+}
+
+
+template<>
+Foam::tmp
+<
+    Foam::GeometricField<Foam::vector, Foam::fvPatchField, Foam::volMesh>
+>
+Foam::fv::fusedGaussLaplacianScheme<Foam::vector, Foam::scalar>::fvcLaplacian
+(
+    const GeometricField<vector, fvPatchField, volMesh>& vf
+)
+{
+    typedef vector Type;
+
+    typedef GeometricField<Type, fvPatchField, volMesh> FieldType;
+    typedef GeometricField<Type, fvsPatchField, surfaceMesh> SurfaceFieldType;
+    typedef typename outerProduct<vector, Type>::type GradType;
+    typedef GeometricField<GradType, fvPatchField, volMesh> GradFieldType;
+
+    const fvMesh& mesh = vf.mesh();
+
+    tmp<FieldType> tresult
+    (
+        new FieldType
+        (
+            IOobject
+            (
+                "laplacian(" + vf.name() + ')',
+                vf.instance(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimensioned<Type>(vf.dimensions()/dimArea, Zero),
+            fvPatchFieldBase::extrapolatedCalculatedType()
+        )
+    );
+    FieldType& result = tresult.ref();
+
+    DebugPout
+        << "fusedGaussLaplacianScheme<vector, GType>::fvcLaplacian on "
+        << vf.name()
+        << " to generate " << result.name() << endl;
+
+
+    const auto tdeltaCoeffs(this->tsnGradScheme_().deltaCoeffs(vf));
+    const auto& deltaCoeffs = tdeltaCoeffs();
+
+
+    if (this->tsnGradScheme_().corrected())
+    {
+        // Calculate sn gradient
+        tmp<SurfaceFieldType> tfaceGrad
+        (
+            new SurfaceFieldType
+            (
+                IOobject
+                (
+                    "snGradCorr("+vf.name()+')',
+                    vf.instance(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                mesh,
+                vf.dimensions()
+            )
+        );
+
+        {
+            // Calculate gradient
+            tmp<GradFieldType> tgGrad
+            (
+                gradScheme<Type>::New
+                (
+                    mesh,
+                    mesh.gradScheme("grad(" + vf.name() + ')')
+                )().grad(vf, "grad(" + vf.name() + ')')
+            );
+            const auto& gGrad = tgGrad();
+
+            // Doing a dotinterpolate with nonOrthCorrectionVectors
+            const auto dotInterpolate = [&]
+            (
+                const vector& area,
+                const scalar lambda,
+
+                const GradType& ownVal,
+                const GradType& neiVal,
+
+                const vector& dotVector,
+
+                Type& result
+            )
+            {
+                result = dotVector&(lambda*(ownVal - neiVal) + neiVal);
+            };
+
+            fvc::interpolate
+            (
+                mesh.surfaceInterpolation::weights(),   // linear interpolation
+                gGrad,                          // volume field
+                mesh.nonOrthCorrectionVectors(),// surface multiplier
+                dotInterpolate,
+                tfaceGrad.ref()
+            );
+        }
+        const auto& faceGrad = tfaceGrad();
+
+        const auto snGrad = [&]
+        (
+            const vector& Sf,
+            const scalar dc,
+            const Type& ownVal,
+            const Type& neiVal,
+            const Type& correction
+        ) -> Type
+        {
+            const auto snGrad(dc*(neiVal-ownVal) + correction);
+            return mag(Sf)*snGrad;
+        };
+
+        fvc::surfaceSnSum
+        (
+            deltaCoeffs,
+            vf,
+            faceGrad,   // face-based addition
+            snGrad,
+            result,
+            false       // avoid boundary evaluation until volume division
+        );
+    }
+    else
+    {
+        const auto snGrad = [&]
+        (
+            const vector& Sf,
+            const scalar dc,
+            const Type& ownVal,
+            const Type& neiVal
+        ) -> Type
+        {
+            const auto snGrad(dc*(neiVal-ownVal));
+            return mag(Sf)*snGrad;
+        };
+
+        fvc::surfaceSnSum
+        (
+            deltaCoeffs,
+            vf,
+            snGrad,
+            result,
+            false       // avoid boundary evaluation until volume division
+        );
+    }
+
+    result.primitiveFieldRef() /= mesh.V();
+    result.correctBoundaryConditions();
+
+    return tresult;
+}
+
+
+declareFvmLaplacianScalarGamma(scalar);
+declareFvmLaplacianScalarGamma(vector);
+declareFvmLaplacianScalarGamma(sphericalTensor);
+declareFvmLaplacianScalarGamma(symmTensor);
+declareFvmLaplacianScalarGamma(tensor);
+
+
+// ************************************************************************* //

src/fused/finiteVolume/fusedLeastSquaresGrad.C

@@ -0,0 +1,115 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2018-2021 OpenCFD Ltd.
+    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 "fusedLeastSquaresGrad.H"
+#include "leastSquaresVectors.H"
+#include "gaussGrad.H"
+#include "fvMesh.H"
+#include "volMesh.H"
+#include "surfaceMesh.H"
+#include "GeometricField.H"
+#include "extrapolatedCalculatedFvPatchField.H"
+#include "fvcSurfaceOps.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::tmp
+<
+    Foam::GeometricField
+    <
+        typename Foam::outerProduct<Foam::vector, Type>::type,
+        Foam::fvPatchField,
+        Foam::volMesh
+    >
+>
+Foam::fv::fusedLeastSquaresGrad<Type>::calcGrad
+(
+    const GeometricField<Type, fvPatchField, volMesh>& vf,
+    const word& name
+) const
+{
+    typedef typename outerProduct<vector, Type>::type GradType;
+    typedef GeometricField<GradType, fvPatchField, volMesh> GradFieldType;
+
+    const fvMesh& mesh = vf.mesh();
+
+    DebugPout<< "fusedLeastSquaresGrad<Type>::calcGrad on " << vf.name()
+        << " with name " << name << endl;
+
+    tmp<GradFieldType> tlsGrad
+    (
+        new GradFieldType
+        (
+            IOobject
+            (
+                name,
+                vf.instance(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimensioned<GradType>(vf.dimensions()/dimLength, Zero),
+            fvPatchFieldBase::extrapolatedCalculatedType()
+        )
+    );
+    GradFieldType& lsGrad = tlsGrad.ref();
+
+    // Get reference to least square vectors
+    const leastSquaresVectors& lsv = leastSquaresVectors::New(mesh);
+    const surfaceVectorField& ownLs = lsv.pVectors();
+    const surfaceVectorField& neiLs = lsv.nVectors();
+
+    const auto differenceOp = [&]
+    (
+        const vector& area,
+        const Type& ownVal,
+        const Type& neiVal
+    ) -> Type
+    {
+        return neiVal - ownVal;
+    };
+
+    fvc::surfaceOp
+    (
+        vf,
+        ownLs,
+        neiLs,
+        differenceOp,
+        lsGrad
+    );
+
+    gaussGrad<Type>::correctBoundaryConditions(vf, lsGrad);
+
+    return tlsGrad;
+}
+
+
+// ************************************************************************* //

src/fused/finiteVolume/fusedLeastSquaresGrad.H

@@ -0,0 +1,129 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2018-2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::fv::fusedLeastSquaresGrad
+
+Group
+    grpFvGradSchemes
+
+Description
+    Second-order gradient scheme using least-squares.
+
+SourceFiles
+    leastSquaresGrad.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef fusedLeastSquaresGrad_H
+#define fusedLeastSquaresGrad_H
+
+#include "gradScheme.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace fv
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class fusedLeastSquaresGrad Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class fusedLeastSquaresGrad
+:
+    public fv::gradScheme<Type>
+{
+    // Private Member Functions
+
+        //- No copy construct
+        fusedLeastSquaresGrad(const fusedLeastSquaresGrad&) = delete;
+
+        //- No copy assignment
+        void operator=(const fusedLeastSquaresGrad&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("fusedLeastSquares");
+
+
+    // Constructors
+
+        //- Construct from mesh
+        fusedLeastSquaresGrad(const fvMesh& mesh)
+        :
+            gradScheme<Type>(mesh)
+        {}
+
+        //- Construct from Istream
+        fusedLeastSquaresGrad(const fvMesh& mesh, Istream&)
+        :
+            gradScheme<Type>(mesh)
+        {}
+
+
+    // Member Functions
+
+        //- Return the gradient of the given field to the gradScheme::grad
+        //- for optional caching
+        virtual tmp
+        <
+            GeometricField
+            <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>
+        > calcGrad
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& vsf,
+            const word& name
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "fusedLeastSquaresGrad.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/fused/finiteVolume/fusedLeastSquaresGrads.C

@@ -0,0 +1,35 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2015 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+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 "fvMesh.H"
+#include "fusedLeastSquaresGrad.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makeFvGradScheme(fusedLeastSquaresGrad)
+
+// ************************************************************************* //

src/fused/finiteVolume/fvcSurfaceOps.H

@@ -0,0 +1,291 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+InNamespace
+    Foam::fvc
+
+Description
+    Surface integrate surfaceField creating a volField.
+    Surface sum a surfaceField creating a volField.
+
+SourceFiles
+    fvcSurfaceOps.C
+
+\*---------------------------------------------------------------------------*/
+
+
+#ifndef fvcSurfaceOps_H
+#define fvcSurfaceOps_H
+
+#include "primitiveFieldsFwd.H"
+#include "volFieldsFwd.H"
+#include "surfaceFieldsFwd.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+//- In-place operations on Fields. Add to FieldFunctions.C ?
+#define INPLACE_PRODUCT_OPERATOR(product, CombineOp, Op, OpFunc)               \
+                                                                               \
+template<class Type1, class Type2>                                             \
+void OpFunc                                                                    \
+(                                                                              \
+    Field<typename product<Type1, Type2>::type>& result,                       \
+    const UList<Type1>& f1,                                                    \
+    const UList<Type2>& f2                                                     \
+)                                                                              \
+{                                                                              \
+    typedef typename product<Type1, Type2>::type resultType;                   \
+    TFOR_ALL_F_OP_F_OP_F                                                       \
+    (resultType, result, CombineOp, Type1, f1, Op, Type2, f2)                  \
+}
+
+
+INPLACE_PRODUCT_OPERATOR(outerProduct, +=, *, multiplyAdd)
+INPLACE_PRODUCT_OPERATOR(outerProduct, -=, *, multiplySubtract)
+
+#undef INPLACE_PRODUCT_OPERATOR
+
+
+
+/*---------------------------------------------------------------------------*\
+                     Namespace fvc functions Declaration
+\*---------------------------------------------------------------------------*/
+
+namespace fvc
+{
+    // Interpolation
+
+        //- Interpolate to face (using cop) and additional face field
+        template<class Type, class FType, class ResultType, class CellToFaceOp>
+        void interpolate
+        (
+            const surfaceScalarField& lambdas,
+            const GeometricField<Type, fvPatchField, volMesh>& vf,
+            const GeometricField<FType, fvsPatchField, surfaceMesh>& sf,
+            const CellToFaceOp& cop,
+            GeometricField<ResultType, fvsPatchField, surfaceMesh>& result
+        );
+
+        //- Interpolate to face (using cop)
+        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
+        );
+
+
+    // Interpolation and accumulation
+
+        //- Interpolate to face (using cop) and accumulate.
+        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 = true
+        );
+
+        //- Interpolate to face (using cop) and accumulate. Additional
+        //- face field
+        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>& sf,
+            const CellToFaceOp& cop,
+            GeometricField<ResultType, fvPatchField, volMesh>& result,
+            const bool doCorrectBoundaryConditions = true
+        );
+
+        //- Interpolate to face (using cop) and accumulate. Additional
+        //- face fields
+        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 = true
+        );
+
+        //- Interpolate to face (using cop) and apply Gauss. Note: uses V(),
+        //  not Vsc()
+        template<class Type, class ResultType, class CellToFaceOp>
+        void GaussOp
+        (
+            const surfaceScalarField& lambdas,
+            const GeometricField<Type, fvPatchField, volMesh>& vf,
+            const CellToFaceOp& cop,
+            GeometricField<ResultType, fvPatchField, volMesh>& result
+        );
+
+
+    // Difference and accumulation
+
+        //- sum of snGrad
+        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
+        );
+        //- sum of snGrad with additional surface field
+        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
+        );
+
+        //- sum of snGrad with additional (interpolated) volField
+        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
+        );
+
+        //- sum of snGrad with additional (interpolated) volfields
+        template
+        <
+            class Type,
+            class GType0,
+            class GType1,
+            class ResultType,
+            class CellToFaceOp
+        >
+        void surfaceSnSum
+        (
+            const surfaceScalarField& weights,
+
+            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
+        );
+
+        //- sum of snGrad with additional surface field
+        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
+        );
+
+
+    // Other
+
+        //- Interpolate to face (using cop) and apply distribution vectors
+        template<class Type, class ResultType, class CellToFaceOp>
+        void surfaceOp
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& vf,
+            const surfaceVectorField& ownLs,
+            const surfaceVectorField& neiLs,
+            const CellToFaceOp& cop,
+            GeometricField<ResultType, fvPatchField, volMesh>& result
+        );
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "fvcSurfaceOps.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

tutorials/incompressible/simpleFoam/pitzDaily_fused/0/U

@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2406                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           uniform (10 0 0);
+    }
+
+    outlet
+    {
+        type            zeroGradient;
+    }
+
+    upperWall
+    {
+        type            noSlip;
+    }
+
+    lowerWall
+    {
+        type            noSlip;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/incompressible/simpleFoam/pitzDaily_fused/0/epsilon

@@ -0,0 +1,53 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2406                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      epsilon;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -3 0 0 0 0];
+
+internalField   uniform 14.855;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           uniform 14.855;
+    }
+
+    outlet
+    {
+        type            zeroGradient;
+    }
+
+    upperWall
+    {
+        type            epsilonWallFunction;
+        value           uniform 14.855;
+    }
+
+    lowerWall
+    {
+        type            epsilonWallFunction;
+        value           uniform 14.855;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/incompressible/simpleFoam/pitzDaily_fused/0/k

@@ -0,0 +1,53 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2406                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      k;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 0.375;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           uniform 0.375;
+    }
+
+    outlet
+    {
+        type            zeroGradient;
+    }
+
+    upperWall
+    {
+        type            kqRWallFunction;
+        value           uniform 0.375;
+    }
+
+    lowerWall
+    {
+        type            kqRWallFunction;
+        value           uniform 0.375;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/incompressible/simpleFoam/pitzDaily_fused/0/nuTilda

@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2406                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      nuTilda;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           uniform 0;
+    }
+
+    outlet
+    {
+        type            zeroGradient;
+    }
+
+    upperWall
+    {
+        type            zeroGradient;
+    }
+
+    lowerWall
+    {
+        type            zeroGradient;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/incompressible/simpleFoam/pitzDaily_fused/0/nut

@@ -0,0 +1,54 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2406                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      nut;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+
+    outlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+
+    upperWall
+    {
+        type            nutkWallFunction;
+        value           uniform 0;
+    }
+
+    lowerWall
+    {
+        type            nutkWallFunction;
+        value           uniform 0;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/incompressible/simpleFoam/pitzDaily_fused/0/omega

@@ -0,0 +1,53 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2406                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      omega;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 -1 0 0 0 0];
+
+internalField   uniform 440.15;
+
+boundaryField
+{
+    inlet
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+
+    outlet
+    {
+        type            zeroGradient;
+    }
+
+    upperWall
+    {
+        type            omegaWallFunction;
+        value           $internalField;
+    }
+
+    lowerWall
+    {
+        type            omegaWallFunction;
+        value           $internalField;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/incompressible/simpleFoam/pitzDaily_fused/0/p

@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2406                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            zeroGradient;
+    }
+
+    outlet
+    {
+        type            fixedValue;
+        value           uniform 0;
+    }
+
+    upperWall
+    {
+        type            zeroGradient;
+    }
+
+    lowerWall
+    {
+        type            zeroGradient;
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/incompressible/simpleFoam/pitzDaily_fused/Allrun

@@ -0,0 +1,22 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
+#------------------------------------------------------------------------------
+
+runApplication blockMesh
+
+runApplication decomposePar
+
+echo "Updating fvSchemes to use Gauss"
+sed "s/GAUSS/Gauss/g" system/fvSchemes.template > system/fvSchemes
+
+runParallel -s Gauss $(getApplication)
+
+
+echo "Updating fvSchemes to use fusedGauss"
+sed "s/GAUSS/fusedGauss/g" system/fvSchemes.template > system/fvSchemes
+
+runParallel -s fusedGauss $(getApplication)
+
+
+#------------------------------------------------------------------------------

tutorials/incompressible/simpleFoam/pitzDaily_fused/README

@@ -0,0 +1,4 @@
+Testing fused discretisation
+============================
+
+Runs both standard Gauss and 'fused' Gauss

tutorials/incompressible/simpleFoam/pitzDaily_fused/constant/transportProperties

@@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2406                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+transportModel  Newtonian;
+
+nu              1e-05;
+
+
+// ************************************************************************* //

tutorials/incompressible/simpleFoam/pitzDaily_fused/constant/turbulenceProperties

@@ -0,0 +1,31 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2406                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType      RAS;
+
+RAS
+{
+    // Tested with kEpsilon, realizableKE, kOmega, kOmegaSST,
+    // ShihQuadraticKE, LienCubicKE.
+    RASModel        kEpsilon;
+
+    turbulence      on;
+
+    printCoeffs     on;
+}
+
+
+// ************************************************************************* //

tutorials/incompressible/simpleFoam/pitzDaily_fused/system/blockMeshDict

@@ -0,0 +1,153 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2406                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+scale   0.001;
+
+vertices
+(
+    (-20.6 0 -0.5)
+    (-20.6 25.4 -0.5)
+    (0 -25.4 -0.5)
+    (0 0 -0.5)
+    (0 25.4 -0.5)
+    (206 -25.4 -0.5)
+    (206 0 -0.5)
+    (206 25.4 -0.5)
+    (290 -16.6 -0.5)
+    (290 0 -0.5)
+    (290 16.6 -0.5)
+
+    (-20.6 0 0.5)
+    (-20.6 25.4 0.5)
+    (0 -25.4 0.5)
+    (0 0 0.5)
+    (0 25.4 0.5)
+    (206 -25.4 0.5)
+    (206 0 0.5)
+    (206 25.4 0.5)
+    (290 -16.6 0.5)
+    (290 0 0.5)
+    (290 16.6 0.5)
+);
+
+negY
+(
+    (2 4 1)
+    (1 3 0.3)
+);
+
+posY
+(
+    (1 4 2)
+    (2 3 4)
+    (2 4 0.25)
+);
+
+posYR
+(
+    (2 1 1)
+    (1 1 0.25)
+);
+
+
+blocks
+(
+    hex (0 3 4 1 11 14 15 12)
+    (18 30 1)
+    simpleGrading (0.5 $posY 1)
+
+    hex (2 5 6 3 13 16 17 14)
+    (180 27 1)
+    edgeGrading (4 4 4 4 $negY 1 1 $negY 1 1 1 1)
+
+    hex (3 6 7 4 14 17 18 15)
+    (180 30 1)
+    edgeGrading (4 4 4 4 $posY $posYR $posYR $posY 1 1 1 1)
+
+    hex (5 8 9 6 16 19 20 17)
+    (25 27 1)
+    simpleGrading (2.5 1 1)
+
+    hex (6 9 10 7 17 20 21 18)
+    (25 30 1)
+    simpleGrading (2.5 $posYR 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+    inlet
+    {
+        type patch;
+        faces
+        (
+            (0 1 12 11)
+        );
+    }
+    outlet
+    {
+        type patch;
+        faces
+        (
+            (8 9 20 19)
+            (9 10 21 20)
+        );
+    }
+    upperWall
+    {
+        type wall;
+        faces
+        (
+            (1 4 15 12)
+            (4 7 18 15)
+            (7 10 21 18)
+        );
+    }
+    lowerWall
+    {
+        type wall;
+        faces
+        (
+            (0 3 14 11)
+            (3 2 13 14)
+            (2 5 16 13)
+            (5 8 19 16)
+        );
+    }
+    frontAndBack
+    {
+        type empty;
+        faces
+        (
+            (0 3 4 1)
+            (2 5 6 3)
+            (3 6 7 4)
+            (5 8 9 6)
+            (6 9 10 7)
+            (11 14 15 12)
+            (13 16 17 14)
+            (14 17 18 15)
+            (16 19 20 17)
+            (17 20 21 18)
+        );
+    }
+);
+
+
+// ************************************************************************* //

tutorials/incompressible/simpleFoam/pitzDaily_fused/system/controlDict

@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2406                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+libs            (fusedFiniteVolume);
+
+application     simpleFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         200;
+
+deltaT          1;
+
+writeControl    timeStep;
+
+writeInterval   100000;
+
+purgeWrite      0;
+
+writeFormat     binary;
+
+writePrecision  16;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable true;
+
+// ************************************************************************* //

tutorials/incompressible/simpleFoam/pitzDaily_fused/system/decomposeParDict

@@ -0,0 +1,25 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2406                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    note        "mesh decomposition control dictionary";
+    object      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//- The total number of domains (mandatory)
+numberOfSubdomains  5;
+
+//- The decomposition method (mandatory)
+method          hierarchical;
+n               (5 1 1);
+
+// ************************************************************************* //

tutorials/incompressible/simpleFoam/pitzDaily_fused/system/fvSchemes

@@ -0,0 +1,63 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2406                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         steadyState;
+}
+
+gradSchemes
+{
+    default         fusedGauss linear;
+}
+
+divSchemes
+{
+    default         none;
+
+    div(phi,U)      bounded fusedGauss linearUpwind grad(U);
+
+    turbulence      bounded fusedGauss limitedLinear 1;
+    div(phi,k)      $turbulence;
+    div(phi,epsilon) $turbulence;
+    div(phi,omega)  $turbulence;
+
+    div(nonlinearStress) fusedGauss linear;
+    div((nuEff*dev2(T(grad(U))))) fusedGauss linear;
+}
+
+laplacianSchemes
+{
+    default         fusedGauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+wallDist
+{
+    method          meshWave;
+}
+
+
+// ************************************************************************* //

tutorials/incompressible/simpleFoam/pitzDaily_fused/system/fvSchemes.template

@@ -0,0 +1,63 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2406                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         steadyState;
+}
+
+gradSchemes
+{
+    default         GAUSS linear;
+}
+
+divSchemes
+{
+    default         none;
+
+    div(phi,U)      bounded GAUSS linearUpwind grad(U);
+
+    turbulence      bounded GAUSS limitedLinear 1;
+    div(phi,k)      $turbulence;
+    div(phi,epsilon) $turbulence;
+    div(phi,omega)  $turbulence;
+
+    div(nonlinearStress) GAUSS linear;
+    div((nuEff*dev2(T(grad(U))))) GAUSS linear;
+}
+
+laplacianSchemes
+{
+    default         GAUSS linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+wallDist
+{
+    method          meshWave;
+}
+
+
+// ************************************************************************* //

tutorials/incompressible/simpleFoam/pitzDaily_fused/system/fvSolution

@@ -0,0 +1,59 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2406                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    p
+    {
+        solver          GAMG;
+        tolerance       1e-06;
+        relTol          0.1;
+        smoother        GaussSeidel;
+    }
+
+    "(U|k|epsilon|omega|f|v2)"
+    {
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-05;
+        relTol          0.1;
+    }
+}
+
+SIMPLE
+{
+    nNonOrthogonalCorrectors 0;
+    consistent      yes;
+
+    residualControl
+    {
+        p               1e-2;
+        U               1e-3;
+        "(k|epsilon|omega|f|v2)" 1e-3;
+    }
+}
+
+relaxationFactors
+{
+    equations
+    {
+        U               0.9; // 0.9 is more stable but 0.95 more convergent
+        ".*"            0.9; // 0.9 is more stable but 0.95 more convergent
+    }
+}
+
+
+// ************************************************************************* //