From 86da88cdfd1a22ef1ba86df52e5f37eac9b43503 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Wed, 25 Sep 2019 12:40:23 +0100
Subject: [PATCH] ENH: Added relaxed nonOrthogonal laplacian corrector scheme

---
 src/finiteVolume/Make/files                   |   1 +
 .../relaxedNonOrthoGaussLaplacianScheme.C     | 276 ++++++++++++++++++
 .../relaxedNonOrthoGaussLaplacianScheme.H     | 210 +++++++++++++
 .../relaxedNonOrthoGaussLaplacianSchemes.C    | 140 +++++++++
 4 files changed, 627 insertions(+)
 create mode 100644 src/finiteVolume/finiteVolume/laplacianSchemes/relaxedNonOrthoGaussLaplacianScheme/relaxedNonOrthoGaussLaplacianScheme.C
 create mode 100644 src/finiteVolume/finiteVolume/laplacianSchemes/relaxedNonOrthoGaussLaplacianScheme/relaxedNonOrthoGaussLaplacianScheme.H
 create mode 100644 src/finiteVolume/finiteVolume/laplacianSchemes/relaxedNonOrthoGaussLaplacianScheme/relaxedNonOrthoGaussLaplacianSchemes.C

diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files
index 37cfc5f1664..8b184cfaf18 100644
--- a/src/finiteVolume/Make/files
+++ b/src/finiteVolume/Make/files
@@ -467,6 +467,7 @@ $(convectionSchemes)/boundedConvectionScheme/boundedConvectionSchemes.C
 laplacianSchemes = finiteVolume/laplacianSchemes
 $(laplacianSchemes)/laplacianScheme/laplacianSchemes.C
 $(laplacianSchemes)/gaussLaplacianScheme/gaussLaplacianSchemes.C
+$(laplacianSchemes)/relaxedNonOrthoGaussLaplacianScheme/relaxedNonOrthoGaussLaplacianSchemes.C
 
 finiteVolume/fvc/fvcFlux.C
 finiteVolume/fvc/fvcMeshPhi.C
diff --git a/src/finiteVolume/finiteVolume/laplacianSchemes/relaxedNonOrthoGaussLaplacianScheme/relaxedNonOrthoGaussLaplacianScheme.C b/src/finiteVolume/finiteVolume/laplacianSchemes/relaxedNonOrthoGaussLaplacianScheme/relaxedNonOrthoGaussLaplacianScheme.C
new file mode 100644
index 00000000000..d0732bf1241
--- /dev/null
+++ b/src/finiteVolume/finiteVolume/laplacianSchemes/relaxedNonOrthoGaussLaplacianScheme/relaxedNonOrthoGaussLaplacianScheme.C
@@ -0,0 +1,276 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "relaxedNonOrthoGaussLaplacianScheme.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>>
+relaxedNonOrthoGaussLaplacianScheme<Type, GType>::fvmLaplacianUncorrected
+(
+    const surfaceScalarField& gammaMagSf,
+    const surfaceScalarField& deltaCoeffs,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    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();
+    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];
+
+        if (pvf.coupled())
+        {
+            fvm.internalCoeffs()[patchi] =
+                pGamma*pvf.gradientInternalCoeffs(pDeltaCoeffs);
+            fvm.boundaryCoeffs()[patchi] =
+               -pGamma*pvf.gradientBoundaryCoeffs(pDeltaCoeffs);
+        }
+        else
+        {
+            fvm.internalCoeffs()[patchi] = pGamma*pvf.gradientInternalCoeffs();
+            fvm.boundaryCoeffs()[patchi] = -pGamma*pvf.gradientBoundaryCoeffs();
+        }
+    }
+
+    return tfvm;
+}
+
+
+template<class Type, class GType>
+tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
+relaxedNonOrthoGaussLaplacianScheme<Type, GType>::gammaSnGradCorr
+(
+    const surfaceVectorField& SfGammaCorr,
+    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,
+            SfGammaCorr.dimensions()
+           *vf.dimensions()*mesh.deltaCoeffs().dimensions()
+        )
+    );
+
+    for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
+    {
+        tgammaSnGradCorr.ref().replace
+        (
+            cmpt,
+            fvc::dotInterpolate(SfGammaCorr, fvc::grad(vf.component(cmpt)))
+        );
+    }
+
+    return tgammaSnGradCorr;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type, class GType>
+tmp<GeometricField<Type, fvPatchField, volMesh>>
+relaxedNonOrthoGaussLaplacianScheme<Type, GType>::fvcLaplacian
+(
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    const fvMesh& mesh = this->mesh();
+
+    tmp<GeometricField<Type, fvPatchField, volMesh>> tLaplacian
+    (
+        fvc::div(this->tsnGradScheme_().snGrad(vf)*mesh.magSf())
+    );
+
+    tLaplacian.ref().rename("laplacian(" + vf.name() + ')');
+
+    return tLaplacian;
+}
+
+
+template<class Type, class GType>
+tmp<fvMatrix<Type>>
+relaxedNonOrthoGaussLaplacianScheme<Type, GType>::fvmLaplacian
+(
+    const GeometricField<GType, fvsPatchField, surfaceMesh>& gamma,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    const fvMesh& mesh = this->mesh();
+
+    typedef GeometricField<Type, fvsPatchField, surfaceMesh> SType;
+
+    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<SType> tfaceFluxCorrection = gammaSnGradCorr(SfGammaCorr, vf);
+
+    if (this->tsnGradScheme_().corrected())
+    {
+        tfaceFluxCorrection.ref() +=
+            SfGammaSn*this->tsnGradScheme_().correction(vf);
+    }
+
+    const word corrName(tfaceFluxCorrection().name());
+
+    tmp<SType> trelaxedCorrection(new SType(tfaceFluxCorrection()));
+
+    const word oldName(corrName + "_0");
+    const scalar relax(vf.mesh().equationRelaxationFactor(oldName));
+
+    const objectRegistry& obr = vf.db();
+    if (obr.foundObject<SType>(oldName))
+    {
+        SType& oldCorrection = obr.lookupObjectRef<SType>(oldName);
+
+        trelaxedCorrection.ref() *= relax;
+        trelaxedCorrection.ref() += (1.0-relax)*oldCorrection;
+
+        oldCorrection = tfaceFluxCorrection;
+    }
+    else
+    {
+        SType* s = new SType(oldName, tfaceFluxCorrection);
+        s->store();
+    }
+
+    fvm.source() -=
+        mesh.V()
+       *fvc::div
+        (
+            trelaxedCorrection()
+        )().primitiveField();
+
+    if (mesh.fluxRequired(vf.name()))
+    {
+        fvm.faceFluxCorrectionPtr() = trelaxedCorrection.ptr();
+    }
+
+    return tfvm;
+}
+
+
+template<class Type, class GType>
+tmp<GeometricField<Type, fvPatchField, volMesh>>
+relaxedNonOrthoGaussLaplacianScheme<Type, GType>::fvcLaplacian
+(
+    const GeometricField<GType, fvsPatchField, surfaceMesh>& gamma,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    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
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/finiteVolume/laplacianSchemes/relaxedNonOrthoGaussLaplacianScheme/relaxedNonOrthoGaussLaplacianScheme.H b/src/finiteVolume/finiteVolume/laplacianSchemes/relaxedNonOrthoGaussLaplacianScheme/relaxedNonOrthoGaussLaplacianScheme.H
new file mode 100644
index 00000000000..ff4a7f6c67a
--- /dev/null
+++ b/src/finiteVolume/finiteVolume/laplacianSchemes/relaxedNonOrthoGaussLaplacianScheme/relaxedNonOrthoGaussLaplacianScheme.H
@@ -0,0 +1,210 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019 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::relaxedNonOrthoGaussLaplacianScheme
+
+Description
+    Basic second-order laplacian using face-gradients and Gauss' theorem.
+
+Usage
+    Minimal example by using \c system/fvSchemes:
+    \verbatim
+    laplacianSchemes
+    {
+        laplacian(<term>)    relaxedNonOrthoGauss <other options>;
+    }
+    \endverbatim
+
+    and by using \c system/fvSolution:
+    \verbatim
+    relaxationFactors
+    {
+        equations
+        {
+            <term>           <relaxation factor>;
+        }
+    }
+    \endverbatim
+
+SourceFiles
+    relaxedNonOrthoGaussLaplacianScheme.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef relaxedNonOrthoGaussLaplacianScheme_H
+#define relaxedNonOrthoGaussLaplacianScheme_H
+
+#include "laplacianScheme.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace fv
+{
+
+/*---------------------------------------------------------------------------*\
+             Class relaxedNonOrthoGaussLaplacianScheme Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type, class GType>
+class relaxedNonOrthoGaussLaplacianScheme
+:
+    public fv::laplacianScheme<Type, GType>
+{
+    // Private Member Functions
+
+        tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> gammaSnGradCorr
+        (
+            const surfaceVectorField& SfGammaCorr,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        );
+
+        //- No copy construct
+        relaxedNonOrthoGaussLaplacianScheme
+        (
+            const relaxedNonOrthoGaussLaplacianScheme&
+        ) = delete;
+
+        //- No copy assignment
+        void operator=(const relaxedNonOrthoGaussLaplacianScheme&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("relaxedNonOrthoGauss");
+
+
+    // Constructors
+
+        //- Construct null
+        relaxedNonOrthoGaussLaplacianScheme(const fvMesh& mesh)
+        :
+            laplacianScheme<Type, GType>(mesh)
+        {}
+
+        //- Construct from Istream
+        relaxedNonOrthoGaussLaplacianScheme(const fvMesh& mesh, Istream& is)
+        :
+            laplacianScheme<Type, GType>(mesh, is)
+        {}
+
+        //- Construct from mesh, interpolation and snGradScheme schemes
+        relaxedNonOrthoGaussLaplacianScheme
+        (
+            const fvMesh& mesh,
+            const tmp<surfaceInterpolationScheme<GType>>& igs,
+            const tmp<snGradScheme<Type>>& sngs
+        )
+        :
+            laplacianScheme<Type, GType>(mesh, igs, sngs)
+        {}
+
+
+    //- Destructor
+    virtual ~relaxedNonOrthoGaussLaplacianScheme() = default;
+
+
+    // Member Functions
+
+        static tmp<fvMatrix<Type>> fvmLaplacianUncorrected
+        (
+            const surfaceScalarField& gammaMagSf,
+            const surfaceScalarField& deltaCoeffs,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        );
+
+        tmp<GeometricField<Type, fvPatchField, volMesh>> fvcLaplacian
+        (
+            const GeometricField<Type, fvPatchField, volMesh>&
+        );
+
+        tmp<fvMatrix<Type>> fvmLaplacian
+        (
+            const GeometricField<GType, fvsPatchField, surfaceMesh>&,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        );
+
+        tmp<GeometricField<Type, fvPatchField, volMesh>> fvcLaplacian
+        (
+            const GeometricField<GType, fvsPatchField, surfaceMesh>&,
+            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>>                                                            \
+relaxedNonOrthoGaussLaplacianScheme<Type, scalar>::fvmLaplacian                \
+(                                                                              \
+    const GeometricField<scalar, fvsPatchField, surfaceMesh>&,                 \
+    const GeometricField<Type, fvPatchField, volMesh>&                         \
+);                                                                             \
+                                                                               \
+template<>                                                                     \
+tmp<GeometricField<Type, fvPatchField, volMesh>>                               \
+relaxedNonOrthoGaussLaplacianScheme<Type, scalar>::fvcLaplacian                \
+(                                                                              \
+    const GeometricField<scalar, fvsPatchField, surfaceMesh>&,                 \
+    const GeometricField<Type, fvPatchField, volMesh>&                         \
+);
+
+
+defineFvmLaplacianScalarGamma(scalar);
+defineFvmLaplacianScalarGamma(vector);
+defineFvmLaplacianScalarGamma(sphericalTensor);
+defineFvmLaplacianScalarGamma(symmTensor);
+defineFvmLaplacianScalarGamma(tensor);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "relaxedNonOrthoGaussLaplacianScheme.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/finiteVolume/laplacianSchemes/relaxedNonOrthoGaussLaplacianScheme/relaxedNonOrthoGaussLaplacianSchemes.C b/src/finiteVolume/finiteVolume/laplacianSchemes/relaxedNonOrthoGaussLaplacianScheme/relaxedNonOrthoGaussLaplacianSchemes.C
new file mode 100644
index 00000000000..7aa0027ec98
--- /dev/null
+++ b/src/finiteVolume/finiteVolume/laplacianSchemes/relaxedNonOrthoGaussLaplacianScheme/relaxedNonOrthoGaussLaplacianSchemes.C
@@ -0,0 +1,140 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "relaxedNonOrthoGaussLaplacianScheme.H"
+#include "fvMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makeFvLaplacianScheme(relaxedNonOrthoGaussLaplacianScheme)
+
+#define declareFvmLaplacianScalarGamma(Type)                                   \
+                                                                               \
+template<>                                                                     \
+Foam::tmp<Foam::fvMatrix<Foam::Type>>                                          \
+Foam::fv::relaxedNonOrthoGaussLaplacianScheme<Foam::Type, Foam::scalar>::      \
+fvmLaplacian                                                                   \
+(                                                                              \
+    const GeometricField<scalar, fvsPatchField, surfaceMesh>& gamma,           \
+    const GeometricField<Type, fvPatchField, volMesh>& vf                      \
+)                                                                              \
+{                                                                              \
+    const fvMesh& mesh = this->mesh();                                         \
+                                                                               \
+    typedef GeometricField<Type, fvsPatchField, surfaceMesh> SType;            \
+                                                                               \
+    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())                                    \
+    {                                                                          \
+        tmp<SType> tCorr(this->tsnGradScheme_().correction(vf));               \
+        const word corrName(tCorr().name());                                   \
+        tmp<SType> tfaceFluxCorrection(gammaMagSf*tCorr);                      \
+                                                                               \
+        tmp<SType> trelaxedCorrection(new SType(tfaceFluxCorrection()));       \
+                                                                               \
+        const word oldName(corrName + "_0");                                   \
+        const scalar relax(vf.mesh().equationRelaxationFactor(corrName));      \
+        const objectRegistry& obr = vf.db();                                   \
+        if (obr.foundObject<SType>(oldName))                                   \
+        {                                                                      \
+            SType& oldCorrection = obr.lookupObjectRef<SType>(oldName);        \
+            trelaxedCorrection.ref() *= relax;                                 \
+            trelaxedCorrection.ref() += (1.0-relax)*oldCorrection;             \
+                                                                               \
+            oldCorrection = trelaxedCorrection();                              \
+        }                                                                      \
+        else                                                                   \
+        {                                                                      \
+            SType* s = new SType(oldName, tfaceFluxCorrection);                \
+            s->store();                                                        \
+        }                                                                      \
+                                                                               \
+        tmp<Field<Type>> tcorr                                                 \
+        (                                                                      \
+            mesh.V()                                                           \
+           *fvc::div                                                           \
+            (                                                                  \
+                trelaxedCorrection()                                           \
+            )().primitiveField()                                               \
+        );                                                                     \
+                                                                               \
+        fvm.source() -= tcorr();                                               \
+                                                                               \
+        if (mesh.fluxRequired(vf.name()))                                      \
+        {                                                                      \
+            fvm.faceFluxCorrectionPtr() = trelaxedCorrection.ptr();            \
+        }                                                                      \
+    }                                                                          \
+                                                                               \
+    return tfvm;                                                               \
+}                                                                              \
+                                                                               \
+                                                                               \
+template<>                                                                     \
+Foam::tmp<Foam::GeometricField<Foam::Type, Foam::fvPatchField, Foam::volMesh>> \
+Foam::fv::relaxedNonOrthoGaussLaplacianScheme<Foam::Type, Foam::scalar>::fvcLaplacian         \
+(                                                                              \
+    const GeometricField<scalar, fvsPatchField, surfaceMesh>& gamma,           \
+    const GeometricField<Type, fvPatchField, volMesh>& vf                      \
+)                                                                              \
+{                                                                              \
+    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;                                                         \
+}
+
+
+declareFvmLaplacianScalarGamma(scalar);
+declareFvmLaplacianScalarGamma(vector);
+declareFvmLaplacianScalarGamma(sphericalTensor);
+declareFvmLaplacianScalarGamma(symmTensor);
+declareFvmLaplacianScalarGamma(tensor);
+
+
+// ************************************************************************* //
-- 
GitLab