Concentrate harmonic interpolation

src/finiteVolume/fields/fvPatchFields/constraint/regionCoupling/regionCouplingFvPatchField.C

@@ -29,7 +29,7 @@ Author
 
 #include "regionCouplingFvPatchField.H"
 #include "symmTransformField.H"
-#include "magLongDelta.H"
+#include "harmonic.H"
 #include "volFields.H"
 #include "surfaceFields.H"
 
@@ -40,111 +40,6 @@ namespace Foam
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-template<class Type>
-tmp<scalarField>
-regionCouplingFvPatchField<Type>::weights
-(
-    const Field<Type>& fOwn,
-    const Field<Type>& fNei
-) const
-{
-    // Implement weights-based stabilised harmonic interpolation using
-    // magnitude of type
-    // Algorithm:
-    // 1) calculate magnitude of internal field and neighbour field
-    // 2) calculate harmonic mean magnitude
-    // 3) express harmonic mean magnitude as: mean = w*mOwn + (1 - w)*mNei
-    // 4) Based on above, calculate w = (mean - mNei)/(mOwn - mNei)
-    // 5) Use weights to interpolate values
-
-    tmp<scalarField> tweights(new scalarField(fOwn.size(), 0.5));
-    scalarField& weights = tweights();
-
-    // Larger small for complex arithmetic accuracy
-    const scalar kSmall = 1000*SMALL;
-
-# if 0
-    // Hrv's treatment
-    scalarField mOwn = mag(fOwn);
-    scalarField mNei = mag(fNei);
-    scalarField mean = 2*(mOwn*mNei)/(mOwn + mNei);
-
-    scalar den;
-
-    forAll (weights, faceI)
-    {
-        den = (mNei[faceI] - mOwn[faceI]);
-
-        // Note: complex arithmetic requires extra accuracy
-        // This is a division of two close subtractions
-        // HJ, 28/Sep/2011
-        if (mag(den) > kSmall)
-        {
-            // Limit weights for round-off safety
-            weights[faceI] =
-                Foam::max(0, Foam::min((mNei[faceI] - mean[faceI])/den, 1));
-        }
-        else
-        {
-            // Use 0.5 weights
-        }
-    }
-
-# else
-
-    // Henrik's treatment
-    const fvPatch& p = this->patch();
-
-    // Note: for interpolation, work with face fields, to allow wall-corrected
-    // diffusivity (eg wall functions) to operate correctly.
-    // HJ, 28/Sep/2011
-
-    // Mag long deltas are identical on both sides.  HJ, 28/Sep/2011
-    const magLongDelta& mld = magLongDelta::New(p.boundaryMesh().mesh());
-
-    scalarField magPhiOwn = mag(fOwn);
-    scalarField magPhiNei = mag(fNei);
-
-    const scalarField& pWeights = p.weights();
-    const scalarField& pDeltaCoeffs = p.deltaCoeffs();
-    const scalarField& pLongDelta = mld.magDelta(p.index());
-
-    forAll (weights, faceI)
-    {
-        scalar mOwn = magPhiOwn[faceI]/(1 - pWeights[faceI]);
-        scalar mNei = magPhiNei[faceI]/pWeights[faceI];
-
-        scalar den = magPhiNei[faceI] - magPhiOwn[faceI];
-
-        // Note: complex arithmetic requires extra accuracy
-        // This is a division of two close subtractions
-        // HJ, 28/Sep/2011
-        if (mag(den) > kSmall)
-        {
-            scalar mean = mOwn*mNei/
-                (
-                    (mOwn + mNei)*
-                    pLongDelta[faceI]*
-                    pDeltaCoeffs[faceI]
-                );
-
-            // Limit weights for round-off safety
-            weights[faceI] =
-                Foam::max(0, Foam::min((magPhiNei[faceI] - mean)/den, 1));
-        }
-        else
-        {
-            // Use 0.5 weights
-            //weights[faceI] = 0.5;
-        }
-    }
-
-#endif
-
-    return tweights;
-}
-
-
 template<class Type>
 const Foam::Field<Type>& regionCouplingFvPatchField<Type>::originalPatchField() const
 {
@@ -152,7 +47,6 @@ const Foam::Field<Type>& regionCouplingFvPatchField<Type>::originalPatchField()
     {
         // Store original field for symmetric evaluation
         // Henrik Rusche, Aug/2011
-        Info << "store original field" << endl;
 
         originalPatchField_ = *this;
         curTimeIndex_ = this->db().time().timeIndex();
@@ -394,7 +288,9 @@ void regionCouplingFvPatchField<Type>::initEvaluate
     );
 
     // Do interpolation
-    scalarField weights = this->weights(fOwn, fNei);
+    harmonic<Type> interp(this->patch().boundaryMesh().mesh());
+
+    scalarField weights = interp.weights(fOwn, fNei, this->patch());
 
     Field<Type>::operator=(weights*fOwn + (1.0 - weights)*fNei);
 
@@ -445,7 +341,9 @@ void regionCouplingFvPatchField<Type>::updateCoeffs()
     Field<Type> fNei = this->patchNeighbourField();
 
     // Do interpolation
-    scalarField weights = this->weights(fOwn, fNei);
+    harmonic<Type> interp(this->patch().boundaryMesh().mesh());
+
+    scalarField weights = interp.weights(fOwn, fNei, this->patch());
 
     Field<Type>::operator=(weights*fOwn + (1.0 - weights)*fNei);
 

src/finiteVolume/fields/fvPatchFields/constraint/regionCoupling/regionCouplingFvPatchField.H

@@ -88,13 +88,6 @@ protected:
             return matrixUpdateBuffer_;
         }
 
-        //- Calculate interpolation weights
-        tmp<scalarField> weights
-        (
-            const Field<Type>& fOwn,
-            const Field<Type>& fNei
-        ) const;
-
 
 public:
 

src/finiteVolume/interpolation/surfaceInterpolation/schemes/harmonic/harmonic.C

@@ -37,4 +37,5 @@ namespace Foam
     makeSurfaceInterpolationScheme(harmonic)
 }
 
+
 // ************************************************************************* //

src/finiteVolume/interpolation/surfaceInterpolation/schemes/harmonic/harmonic.H

@@ -1,4 +1,4 @@
-/*---------------------------------------------------------------------------*\
+/*---------------------------------------------------------------------------* \
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
@@ -108,135 +108,16 @@ public:
         virtual tmp<surfaceScalarField> weights
         (
             const GeometricField<Type, fvPatchField, volMesh>& phi
-        ) const
+        ) const;
-        {
-            tmp<surfaceScalarField> tw
-            (
-                new surfaceScalarField
-                (
-                    IOobject
-                    (
-                        "harmonicWeightingFactors" + phi.name(),
-                        this->mesh().time().timeName(),
-                        this->mesh()
-                    ),
-                    this->mesh() ,
-                    dimless
-                )
-            );
 
-            const surfaceScalarField& deltaCoeffs = this->mesh().deltaCoeffs();
+        //- Return the interpolation weighting factors for a patch
-            const surfaceScalarField& weights = this->mesh().weights();
+        tmp<scalarField> weights
+        (
+            const Field<Type>& fOwn,
+            const Field<Type>& fNei,
+            const fvPatch& patch
+        ) const;
 
-            const magLongDelta& mld = magLongDelta::New(this->mesh());
-            const surfaceScalarField& longDelta = mld.magDelta();
-
-            surfaceScalarField& w = tw();
-
-            const unallocLabelList& owner = this->mesh().owner();
-            const unallocLabelList& neighbour = this->mesh().neighbour();
-
-            scalarField magPhi = mag(phi);
-
-            scalarField& wIn = w.internalField();
-
-            // Larger small for complex arithmetic accuracy
-            const scalar kSmall = 1000*SMALL;
-
-            // Calculate internal weights using field magnitude
-            forAll (owner, faceI)
-            {
-                label own = owner[faceI];
-                label nei = neighbour[faceI];
-
-                scalar mOwn = magPhi[own]/(1 - weights[faceI]);
-                scalar mNei = magPhi[nei]/weights[faceI];
-
-                scalar den = magPhi[nei] - magPhi[own];
-
-                scalar mean = mOwn*mNei/
-                    ((mOwn + mNei)*longDelta[faceI]*deltaCoeffs[faceI]);
-
-                // Note: complex arithmetic requires extra accuracy
-                // This is a division of two close subtractions
-                // HJ, 28/Sep/2011
-                if (mag(den) > kSmall)
-                {
-                    // Limit weights for round-off safety
-                    wIn[faceI] =
-                        Foam::max(0, Foam::min((magPhi[nei] - mean)/den, 1));
-                }
-                else
-                {
-                    wIn[faceI] = 0.5;
-                }
-            }
-
-            forAll (phi.boundaryField(), pi)
-            {
-                fvsPatchScalarField& wp = w.boundaryField()[pi];
-
-                const fvPatchField<Type>& pPhi = phi.boundaryField()[pi];
-
-                if (pPhi.coupled())
-                {
-                    const scalarField& pWeights = weights.boundaryField()[pi];
-
-                    scalarField magPhiOwn = mag(pPhi.patchInternalField());
-                    scalarField magPhiNei = mag(pPhi.patchNeighbourField());
-
-                    const scalarField& pDeltaCoeffs =
-                        deltaCoeffs.boundaryField()[pi];
-
-                    const scalarField& pLongDelta = mld.magDelta(pi);
-
-                    // Calculate internal weights using field magnitude
-                    forAll (pPhi, faceI)
-                    {
-                        scalar mOwn = magPhiOwn[faceI]/(1 - pWeights[faceI]);
-                        scalar mNei = magPhiNei[faceI]/pWeights[faceI];
-
-                        scalar den = magPhiNei[faceI] - magPhiOwn[faceI];
-
-                        // Note: complex arithmetic requires extra accuracy
-                        // This is a division of two close subtractions
-                        // HJ, 28/Sep/2011
-                        if (mag(den) > kSmall)
-                        {
-                            scalar mean = mOwn*mNei/
-                                (
-                                    (mOwn + mNei)*
-                                    pLongDelta[faceI]*
-                                    pDeltaCoeffs[faceI]
-                                );
-
-                            // Limit weights for round-off safety
-                            wp[faceI] =
-                                Foam::max
-                                (
-                                    0,
-                                    Foam::min
-                                    (
-                                        (magPhiNei[faceI] - mean)/den,
-                                        1
-                                    )
-                                );
-                        }
-                        else
-                        {
-                            wp[faceI] = 0.5;
-                        }
-                    }
-                }
-                else
-                {
-                    // Boundary weights for uncoupled patches are 1
-                    wp = 1;
-                }
-            }
-
-            return tw;
-        }
 };
 
 
@@ -246,6 +127,12 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#ifdef NoRepository
+#   include "harmonicTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //

src/finiteVolume/interpolation/surfaceInterpolation/schemes/harmonic/harmonicTemplates.C

@@ -0,0 +1,201 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+    Harmonic-mean differencing scheme class.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::tmp<Foam::surfaceScalarField> Foam::harmonic<Type>::weights
+(
+    const GeometricField<Type, fvPatchField, volMesh>& phi
+) const
+{
+    // Implement weights-based stabilised harmonic interpolation using
+    // magnitude of type
+    // Algorithm:
+    // 1) calculate magnitude of internal field and neighbour field
+    // 2) calculate harmonic mean magnitude
+    // 3) express harmonic mean magnitude as: mean = w*mOwn + (1 - w)*mNei
+    // 4) Based on above, calculate w = (mean - mNei)/(mOwn - mNei)
+    // 5) Use weights to interpolate values
+
+    tmp<surfaceScalarField> tw
+    (
+        new surfaceScalarField
+        (
+            IOobject
+            (
+                "harmonicWeightingFactors" + phi.name(),
+                this->mesh().time().timeName(),
+                this->mesh()
+            ),
+            this->mesh() ,
+            dimless
+        )
+    );
+
+    const surfaceScalarField& deltaCoeffs = this->mesh().deltaCoeffs();
+    const surfaceScalarField& weights = this->mesh().weights();
+
+    const magLongDelta& mld = magLongDelta::New(this->mesh());
+    const surfaceScalarField& longDelta = mld.magDelta();
+
+    surfaceScalarField& w = tw();
+
+    const unallocLabelList& owner = this->mesh().owner();
+    const unallocLabelList& neighbour = this->mesh().neighbour();
+
+    scalarField magPhi = mag(phi);
+
+    scalarField& wIn = w.internalField();
+
+    // Larger small for complex arithmetic accuracy
+    const scalar kSmall = 1000*SMALL;
+
+    // Calculate internal weights using field magnitude
+    forAll (owner, faceI)
+    {
+        label own = owner[faceI];
+        label nei = neighbour[faceI];
+
+        scalar mOwn = magPhi[own]/(1 - weights[faceI]);
+        scalar mNei = magPhi[nei]/weights[faceI];
+
+        scalar den = magPhi[nei] - magPhi[own];
+
+        scalar mean = mOwn*mNei/
+            ((mOwn + mNei)*longDelta[faceI]*deltaCoeffs[faceI]);
+
+        // Note: complex arithmetic requires extra accuracy
+        // This is a division of two close subtractions
+        // HJ, 28/Sep/2011
+        if (mag(den) > kSmall)
+        {
+            // Limit weights for round-off safety
+            wIn[faceI] =
+                Foam::max(0, Foam::min((magPhi[nei] - mean)/den, 1));
+        }
+        else
+        {
+            wIn[faceI] = 0.5;
+        }
+    }
+
+    forAll (phi.boundaryField(), pi)
+    {
+        fvsPatchScalarField& wp = w.boundaryField()[pi];
+
+        const fvPatchField<Type>& pPhi = phi.boundaryField()[pi];
+
+        if (pPhi.coupled())
+        {
+            wp = this->weights
+            (
+                pPhi.patchInternalField(),
+                pPhi.patchNeighbourField(),
+                pPhi.patch()
+            );
+        }
+        else
+        {
+            // Boundary weights for uncoupled patches are 1
+            wp = 1;
+        }
+    }
+
+    return tw;
+}
+
+
+template<class Type>
+Foam::tmp<Foam::scalarField> Foam::harmonic<Type>::weights
+(
+    const Field<Type>& fOwn,
+    const Field<Type>& fNei,
+    const fvPatch& patch
+) const
+{
+    // Implement weights-based stabilised harmonic interpolation using
+    // magnitude of type
+    // Algorithm:
+    // 1) calculate magnitude of internal field and neighbour field
+    // 2) calculate harmonic mean magnitude
+    // 3) express harmonic mean magnitude as: mean = w*mOwn + (1 - w)*mNei
+    // 4) Based on above, calculate w = (mean - mNei)/(mOwn - mNei)
+    // 5) Use weights to interpolate values
+
+    tmp<scalarField> tweights(new scalarField(fOwn.size(), 0.5));
+    scalarField& weights = tweights();
+
+    // Larger small for complex arithmetic accuracy
+    const scalar kSmall = 1000*SMALL;
+
+    // Mag long deltas are identical on both sides.  HJ, 28/Sep/2011
+    const magLongDelta& mld = magLongDelta::New(this->mesh());
+
+    scalarField magPhiOwn = mag(fOwn);
+    scalarField magPhiNei = mag(fNei);
+
+    const scalarField& pWeights = patch.weights();
+    const scalarField& pDeltaCoeffs = patch.deltaCoeffs();
+    const scalarField& pLongDelta = mld.magDelta(patch.index());
+
+    forAll (weights, faceI)
+    {
+        scalar mOwn = magPhiOwn[faceI]/(1 - pWeights[faceI]);
+        scalar mNei = magPhiNei[faceI]/pWeights[faceI];
+
+        scalar den = magPhiNei[faceI] - magPhiOwn[faceI];
+
+        // Note: complex arithmetic requires extra accuracy
+        // This is a division of two close subtractions
+        // HJ, 28/Sep/2011
+        if (mag(den) > kSmall)
+        {
+            scalar mean = mOwn*mNei/
+                (
+                    (mOwn + mNei)*
+                    pLongDelta[faceI]*
+                    pDeltaCoeffs[faceI]
+                );
+
+            // Limit weights for round-off safety
+            weights[faceI] =
+                Foam::max(0, Foam::min((magPhiNei[faceI] - mean)/den, 1));
+        }
+        else
+        {
+            // Use 0.5 weights
+        }
+    }
+
+    return tweights;
+}
+
+
+// ************************************************************************* //