Merge branch 'feature-overset-lsq' into 'develop'

ENH: overset: new solvers, new stencil

See merge request Development/OpenFOAM-plus!180

applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/pEqn.H

@@ -7,8 +7,8 @@ surfaceScalarField faceMask(localMin<scalar>(mesh).interpolate(cellMask));
 
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", faceMask*fvc::interpolate(rho*rAU));
-volVectorField HbyA("HbyA", constrainHbyA(rAU*UEqn.H(), U, p));
+volVectorField HbyA("HbyA", U);
-//mesh.interpolate(HbyA);
+HbyA = constrainHbyA(rAU*UEqn.H(), U, p);
 
 if (pimple.nCorrPISO() <= 1)
 {

applications/solvers/incompressible/simpleFoam/overSimpleFoam/UEqn.H

@@ -16,6 +16,9 @@
 
     fvOptions.constrain(UEqn);
 
+    if (simple.momentumPredictor())
+    {
         solve(UEqn == -cellMask*fvc::grad(p));
+    }
 
     fvOptions.correct(U);

applications/solvers/incompressible/simpleFoam/overSimpleFoam/pEqn.H

@@ -5,7 +5,7 @@
     surfaceScalarField rAUf("rAUf", faceMask*fvc::interpolate(rAU));
 
     volVectorField HbyA("HbyA", U);
-    HbyA = constrainHbyA(rAU*UEqn.H(), U, p);
+    HbyA = constrainHbyA(cellMask*rAU*UEqn.H(), U, p);
 
     //mesh.interpolate(HbyA);
     if (massFluxInterpolation)

applications/solvers/multiphase/interCondensatingEvaporatingFoam/createFields.H

@@ -138,3 +138,4 @@ volScalarField pDivU
     mesh,
     dimensionedScalar("pDivU", p.dimensions()/dimTime, 0)
 );
+

applications/solvers/multiphase/interFoam/overInterDyMFoam/readControls.H

@@ -1,6 +1,6 @@
 #include "readTimeControls.H"
 
-correctPhi = pimple.dict().lookupOrDefault<Switch>("correctPhi", true);
+correctPhi = pimple.dict().lookupOrDefault<Switch>("correctPhi", false);
 
 checkMeshCourantNo =
     pimple.dict().lookupOrDefault<Switch>("checkMeshCourantNo", false);

applications/test/rigidBodyDynamics/pendulum/pendulum

@@ -32,6 +32,7 @@ bodies
         type            sphere;
         mass            1;
         radius          0.05;
+        centreOfMass    (0 10 0);
         transform       (1 0 0 0 1 0 0 0 1) (0 -1 0);
         mergeWith       hinge;
     }

applications/test/rigidBodyDynamics/pendulumAndSpring/pendulumAndSpring

@@ -14,6 +14,7 @@ bodies
         parent          root;
         mass            6e4;
         radius          0.01;
+        centreOfMass    (0 0 0);
         transform       (1 0 0 0 1 0 0 0 1) (0 0 0);
         joint
         {
@@ -37,6 +38,7 @@ bodies
         parent          M;
         mass            6e3;
         radius          0.01;
+        centreOfMass    (0 0 0);
         transform       (1 0 0 0 1 0 0 0 1) (0 -5 0);
         mergeWith       M;
     }

src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H

@@ -233,10 +233,10 @@ public:
         void reorder(const labelUList&, const bool validBoundary);
 
         //- writeData member function required by regIOobject
-        bool writeData(Ostream&) const;
+        virtual bool writeData(Ostream&) const;
 
         //- Write using given format, version and form uncompression
-        bool writeObject
+        virtual bool writeObject
         (
             IOstream::streamFormat fmt,
             IOstream::versionNumber ver,

src/finiteVolume/fvMesh/wallDist/patchDistMethods/Poisson/PoissonPatchDistMethod.C

@@ -101,7 +101,18 @@ bool Foam::patchDistMethods::Poisson::correct
     volVectorField gradyPsi(fvc::grad(yPsi));
     volScalarField magGradyPsi(mag(gradyPsi));
 
-    y = sqrt(magSqr(gradyPsi) + 2*yPsi) - magGradyPsi;
+    // Need to stabilise the y for overset meshes since the holed cells
+    // keep the initial value (0.0) so the gradient of that will be
+    // zero as well. Turbulence models do not like zero wall distance.
+    y = max
+    (
+        sqrt(magSqr(gradyPsi) + 2*yPsi) - magGradyPsi,
+        dimensionedScalar("smallY", dimLength, SMALL)
+    );
+
+    // For overset: enforce smooth y field (yPsi is smooth, magGradyPsi is
+    // not)
+    mesh_.interpolate(y);
 
     // Need to stabilise the y for overset meshes since the holed cells
     // keep the initial value (0.0) so the gradient of that will be

src/overset/Make/files

@@ -6,6 +6,7 @@ cellCellStencil/inverseDistance/waveMethod.C
 cellCellStencil/inverseDistance/meshToMeshData.C
 cellCellStencil/trackingInverseDistance/voxelMeshSearch.C
 cellCellStencil/trackingInverseDistance/trackingInverseDistanceCellCellStencil.C
+cellCellStencil/leastSquares/leastSquaresCellCellStencil.C
 
 dynamicOversetFvMesh/dynamicOversetFvMesh.C
 
@@ -21,6 +22,9 @@ oversetPolyPatch/oversetGAMGInterfaceField.C
 oversetPolyPatch/oversetPointPatch.C
 oversetPolyPatch/oversetPointPatchFields.C
 
+oversetPolyPatch/semiImplicitOversetFvPatchFields.C
+oversetPolyPatch/semiImplicitOversetGAMGInterfaceField.C
+
 oversetAdjustPhi/oversetAdjustPhi.C
 
 regionsToCell/regionsToCell.C

src/overset/cellCellStencil/cellCellStencil/cellCellStencil.C

@@ -94,22 +94,22 @@ Foam::cellCellStencil::~cellCellStencil()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-const Foam::labelIOList& Foam::cellCellStencil::zoneID() const
+const Foam::labelIOList& Foam::cellCellStencil::zoneID(const fvMesh& mesh)
 {
-    if (!mesh_.foundObject<labelIOList>("zoneID"))
+    if (!mesh.foundObject<labelIOList>("zoneID"))
     {
         labelIOList* zoneIDPtr = new labelIOList
         (
             IOobject
             (
                 "zoneID",
-                mesh_.facesInstance(),
+                mesh.facesInstance(),
                 polyMesh::meshSubDir,
-                mesh_,
+                mesh,
                 IOobject::NO_READ,
                 IOobject::NO_WRITE
             ),
-            mesh_.nCells()
+            mesh.nCells()
         );
         labelIOList& zoneID = *zoneIDPtr;
 
@@ -118,13 +118,13 @@ const Foam::labelIOList& Foam::cellCellStencil::zoneID() const
             IOobject
             (
                 "zoneID",
-                mesh_.time().findInstance(mesh_.dbDir(), "zoneID"),
+                mesh.time().findInstance(mesh.dbDir(), "zoneID"),
-                mesh_,
+                mesh,
                 IOobject::MUST_READ,
                 IOobject::NO_WRITE,
                 false
             ),
-            mesh_
+            mesh
         );
         forAll(volZoneID, cellI)
         {
@@ -133,7 +133,7 @@ const Foam::labelIOList& Foam::cellCellStencil::zoneID() const
 
         zoneIDPtr->store();
     }
-    return mesh_.lookupObject<labelIOList>("zoneID");
+    return mesh.lookupObject<labelIOList>("zoneID");
 }
 
 

src/overset/cellCellStencil/cellCellStencil/cellCellStencil.H

@@ -168,12 +168,27 @@ public:
         //  calculated, 1 = use interpolated)
         virtual const scalarList& cellInterpolationWeight() const = 0;
 
+        //- Calculate weights for a single acceptor
+        virtual void stencilWeights
+        (
+            const point& sample,
+            const pointList& donorCcs,
+            scalarList& weights
+        ) const = 0;
+
         //- Helper: is stencil fully local
         bool localStencil(const labelUList&) const;
 
         //- Helper: get reference to registered zoneID. Loads volScalarField
         //  if not registered.
-        const labelIOList& zoneID() const;
+        static const labelIOList& zoneID(const fvMesh&);
+
+        //- Helper: get reference to registered zoneID. Loads volScalarField
+        //  if not registered.
+        const labelIOList& zoneID() const
+        {
+            return zoneID(mesh_);
+        }
 
         //- Helper: create cell-cell addressing in global numbering
         static void globalCellCells

src/overset/cellCellStencil/cellCellStencil/cellCellStencilObject.H

@@ -154,6 +154,17 @@ public:
         {
             return stencilPtr_().cellInterpolationWeight();
         }
+
+        //- Calculate weights for a single acceptor
+        virtual void stencilWeights
+        (
+            const point& sample,
+            const pointList& donorCcs,
+            scalarList& weights
+        ) const
+        {
+            stencilPtr_().stencilWeights(sample, donorCcs, weights);
+        }
 };
 
 

src/overset/cellCellStencil/cellVolumeWeight/cellVolumeWeightCellCellStencil.C

@@ -1100,4 +1100,39 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
 }
 
 
+void Foam::cellCellStencils::cellVolumeWeight::stencilWeights
+(
+    const point& sample,
+    const pointList& donorCcs,
+    scalarList& weights
+) const
+{
+    // Inverse-distance weighting
+
+    weights.setSize(donorCcs.size());
+    scalar sum = 0.0;
+    forAll(donorCcs, i)
+    {
+        scalar d = mag(sample-donorCcs[i]);
+
+        if (d > ROOTVSMALL)
+        {
+            weights[i] = 1.0/d;
+            sum += weights[i];
+        }
+        else
+        {
+            // Short circuit
+            weights = 0.0;
+            weights[i] = 1.0;
+            return;
+        }
+    }
+    forAll(weights, i)
+    {
+        weights[i] /= sum;
+    }
+}
+
+
 // ************************************************************************* //

src/overset/cellCellStencil/cellVolumeWeight/cellVolumeWeightCellCellStencil.H

@@ -227,6 +227,15 @@ public:
         {
             return cellInterpolationWeight_;
         }
+
+        //- Calculate inverse distance weights for a single acceptor. Revert
+        //  to inverse distance (so not consistent with volume overlap!)
+        virtual void stencilWeights
+        (
+            const point& sample,
+            const pointList& donorCcs,
+            scalarList& weights
+        ) const;
 };
 
 

src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.C

@@ -155,6 +155,7 @@ void Foam::cellCellStencils::inverseDistance::fill
 void Foam::cellCellStencils::inverseDistance::markBoundaries
 (
     const fvMesh& mesh,
+    const vector& smallVec,
 
     const boundBox& bb,
     const labelVector& nDivs,
@@ -188,6 +189,9 @@ void Foam::cellCellStencils::inverseDistance::markBoundaries
 
                 // Mark in voxel mesh
                 boundBox faceBb(pp.points(), pp[i], false);
+                faceBb.min() -= smallVec;
+                faceBb.max() += smallVec;
+
                 if (bb.overlaps(faceBb))
                 {
                     fill(patchTypes, bb, nDivs, faceBb, patchCellType::PATCH);
@@ -213,6 +217,9 @@ void Foam::cellCellStencils::inverseDistance::markBoundaries
 
                 // Mark in voxel mesh
                 boundBox faceBb(pp.points(), pp[i], false);
+                faceBb.min() -= smallVec;
+                faceBb.max() += smallVec;
+
                 if (bb.overlaps(faceBb))
                 {
                     fill(patchTypes, bb, nDivs, faceBb, patchCellType::OVERSET);
@@ -333,6 +340,10 @@ void Foam::cellCellStencils::inverseDistance::markPatchesAsHoles
             forAll(tgtCellMap, tgtCelli)
             {
                 label celli = tgtCellMap[tgtCelli];
+                treeBoundBox cBb(cellBb(mesh_, celli));
+                cBb.min() -= smallVec_;
+                cBb.max() += smallVec_;
+
                 if
                 (
                     overlaps
@@ -340,7 +351,7 @@ void Foam::cellCellStencils::inverseDistance::markPatchesAsHoles
                         srcPatchBb,
                         zoneDivs,
                         srcPatchTypes,
-                        cellBb(mesh_, celli),
+                        cBb,
                         patchCellType::PATCH
                     )
                 )
@@ -399,6 +410,9 @@ void Foam::cellCellStencils::inverseDistance::markPatchesAsHoles
                 forAll(tgtCellMap, tgtCelli)
                 {
                     label celli = tgtCellMap[tgtCelli];
+                    treeBoundBox cBb(cellBb(mesh_, celli));
+                    cBb.min() -= smallVec_;
+                    cBb.max() += smallVec_;
                     if
                     (
                         overlaps
@@ -406,7 +420,7 @@ void Foam::cellCellStencils::inverseDistance::markPatchesAsHoles
                             srcPatchBb,
                             zoneDivs,
                             srcPatchTypes,
-                            cellBb(mesh_, celli),
+                            cBb,
                             patchCellType::PATCH
                         )
                     )
@@ -501,7 +515,9 @@ void Foam::cellCellStencils::inverseDistance::markDonors
         label celli = tgtCellMap[tgtCelli];
         if (allStencil[celli].empty())
         {
-            const treeBoundBox subBb(cellBb(mesh_, celli));
+            treeBoundBox subBb(cellBb(mesh_, celli));
+            subBb.min() -= smallVec_;
+            subBb.max() += smallVec_;
 
             forAll(srcOverlapProcs, i)
             {
@@ -1194,12 +1210,20 @@ void Foam::cellCellStencils::inverseDistance::walkFront
     {
         if (isA<oversetFvPatch>(fvm[patchI]))
         {
-            forAll(fvm[patchI], i)
+            const labelList& fc = fvm[patchI].faceCells();
+            forAll(fc, i)
             {
+                label cellI = fc[i];
+                if (allCellTypes[cellI] == INTERPOLATED)
+                {
+                    // Note that acceptors might have been marked hole if
+                    // there are no donors in which case we do not want to
+                    // walk this out. This is an extreme situation.
                     isFront[fvm[patchI].start()+i] = true;
                 }
             }
         }
+    }
 
 
     // Outside of 'hole' region
@@ -1353,12 +1377,12 @@ void Foam::cellCellStencils::inverseDistance::walkFront
 }
 
 
-void Foam::cellCellStencils::inverseDistance::calcStencilWeights
+void Foam::cellCellStencils::inverseDistance::stencilWeights
 (
     const point& sample,
     const pointList& donorCcs,
     scalarList& weights
-)
+) const
 {
     // Inverse-distance weighting
 
@@ -1500,7 +1524,7 @@ void Foam::cellCellStencils::inverseDistance::createStencil
         {
             label cellI = donorCells[i];
             const pointList& donorCentres = donorCellCentres[cellI];
-            calcStencilWeights
+            stencilWeights
             (
                 samples[cellI],
                 donorCentres,
@@ -1565,6 +1589,7 @@ Foam::cellCellStencils::inverseDistance::inverseDistance
 :
     cellCellStencil(mesh),
     dict_(dict),
+    smallVec_(vector::zero),
     cellTypes_(labelList(mesh.nCells(), CALCULATED)),
     interpolationCells_(0),
     cellInterpolationMap_(),
@@ -1605,6 +1630,9 @@ bool Foam::cellCellStencils::inverseDistance::update()
 {
     scalar layerRelax(dict_.lookupOrDefault("layerRelax", 1.0));
 
+    scalar tol = dict_.lookupOrDefault("tolerance", 1e-10);
+    smallVec_ = mesh_.bounds().span()*tol;
+
     const labelIOList& zoneID = this->zoneID();
 
     label nZones = gMax(zoneID)+1;
@@ -1745,6 +1773,7 @@ bool Foam::cellCellStencils::inverseDistance::update()
         markBoundaries
         (
             meshParts[zoneI].subMesh(),
+            smallVec_,
 
             patchBb[zoneI][Pstream::myProcNo()],
             patchDivisions[zoneI],
@@ -1758,7 +1787,7 @@ bool Foam::cellCellStencils::inverseDistance::update()
 
     // Print a bit
     {
-        Info<< typeName << " : detected " << nZones
+        Info<< type() << " : detected " << nZones
             << " mesh regions" << endl;
         Info<< incrIndent;
         forAll(nCellsPerZone, zoneI)
@@ -1862,6 +1891,16 @@ bool Foam::cellCellStencils::inverseDistance::update()
                     }
                     else
                     {
+                        // If there are no donors we can either set the
+                        // acceptor cell to 'hole' i.e. nothing gets calculated
+                        // if the acceptor cells go outside the donor meshes or
+                        // we can set it to 'calculated' to have something
+                        // like an acmi type behaviour where only covered
+                        // acceptor cell use the interpolation and non-covered
+                        // become calculated. Uncomment below line. Note: this
+                        // should be switchable?
+                        //allCellTypes[cellI] = CALCULATED;
+
                         allCellTypes[cellI] = HOLE;
                     }
                 }
@@ -1915,7 +1954,7 @@ bool Foam::cellCellStencils::inverseDistance::update()
         // Dump stencil
         mkDir(mesh_.time().timePath());
         OBJstream str(mesh_.time().timePath()/"injectionStencil.obj");
-        Pout<< typeName << " : dumping injectionStencil to "
+        Pout<< type() << " : dumping injectionStencil to "
             << str.name() << endl;
         pointField cc(mesh_.cellCentres());
         cellInterpolationMap_.distribute(cc);
@@ -1973,7 +2012,7 @@ bool Foam::cellCellStencils::inverseDistance::update()
         // Dump stencil
         mkDir(mesh_.time().timePath());
         OBJstream str(mesh_.time().timePath()/"stencil.obj");
-        Pout<< typeName << " : dumping to " << str.name() << endl;
+        Pout<< type() << " : dumping to " << str.name() << endl;
         pointField cc(mesh_.cellCentres());
         cellInterpolationMap_.distribute(cc);
 

src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.H

@@ -74,6 +74,9 @@ protected:
         //- Dictionary of motion control parameters
         const dictionary dict_;
 
+        //- Small perturbation vector for geometric tests
+        vector smallVec_;
+
         //- Per cell the cell type
         labelList cellTypes_;
 
@@ -143,6 +146,8 @@ protected:
         static void markBoundaries
         (
             const fvMesh& mesh,
+            const vector& smallVec,
+
             const boundBox& bb,
             const labelVector& nDivs,
             PackedList<2>& patchTypes,
@@ -237,15 +242,6 @@ protected:
             scalarField& allWeight
         ) const;
 
-
-        //- Calculate inverse distance weights
-        static void calcStencilWeights
-        (
-            const point& sample,
-            const pointList& donorCcs,
-            scalarList& weights
-        );
-
         //- Create stencil starting from the donor containing the acceptor
         virtual void createStencil(const globalIndex&);
 
@@ -319,6 +315,14 @@ public:
         {
             return cellInterpolationWeight_;
         }
+
+        //- Calculate inverse distance weights for a single acceptor
+        virtual void stencilWeights
+        (
+            const point& sample,
+            const pointList& donorCcs,
+            scalarList& weights
+        ) const;
 };
 
 

src/overset/cellCellStencil/leastSquares/leastSquaresCellCellStencil.C

@@ -0,0 +1,177 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify i
+    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 "leastSquaresCellCellStencil.H"
+#include "addToRunTimeSelectionTable.H"
+#include "SVD.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace cellCellStencils
+{
+    defineTypeNameAndDebug(leastSquares, 0);
+    addToRunTimeSelectionTable(cellCellStencil, leastSquares, mesh);
+}
+}
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::cellCellStencils::leastSquares::stencilWeights
+(
+    const point& sample,
+    const pointList& donorCcs,
+    scalarList& weights
+) const
+{
+    // Implicit least squares weighting
+    // Number of donors
+    label nD = donorCcs.size();
+
+    weights.setSize(nD);
+
+    // List for distance vectors and LSQ weights
+    List<vector> d(nD);
+    scalarList LSQw(nD);
+
+    // Sum of weights
+    scalar W = 0;
+
+    // Sum of weighted distance vectors
+    vector dw = vector::zero;
+
+    RectangularMatrix<scalar> A(nD, 3);
+
+    bool shortC = false;
+
+    // Compute distance vectors and fill rectangular matrix
+    forAll(donorCcs, j)
+    {
+        // Neighbour weights
+        d[j] = donorCcs[j] - sample;
+
+        // Check for short-circuiting if zero distance
+        // is detected with respect to any donor
+        if (mag(d[j]) < ROOTVSMALL)
+        {
+            shortC = true;
+            break;
+        }
+
+        LSQw[j] = 1.0/magSqr(d[j]);
+
+        // T matrix
+        vector wd = LSQw[j]*d[j];
+        A[j][0] = wd.x();
+        A[j][1] = wd.y();
+        A[j][2] = wd.z();
+
+        // Sum of weighted distance vectors
+        dw += wd;
+
+        // Sum of weights
+        W += LSQw[j];
+    }
+
+    if (!shortC)
+    {
+        // Use Singular Value Decomposition to avoid problems
+        // with 1D, 2D stencils
+        SVD svd(A.T()*A, SMALL);
+
+        // Least squares vectors
+        RectangularMatrix<scalar> ATAinvAT(svd.VSinvUt()*A.T());
+
+        scalar saveDiag(W);
+
+        // Diagonal coefficient
+        for (label i = 0; i < 3; i++)
+        {
+            // Get row
+            scalarList Row(UList<scalar>(ATAinvAT[i], nD));
+
+            forAll(donorCcs, j)
+            {
+                W -= Row[j]*LSQw[j]*dw[i];
+            }
+        }
+
+        if (mag(W) < SMALL*mag(saveDiag))
+        {
+            shortC = true;
+        }
+        else
+        {
+            // Compute final neighbour weights with  additional scaling
+            forAll(donorCcs, j)
+            {
+                weights[j] =
+                (
+                    LSQw[j]
+                - ATAinvAT[0][j]*LSQw[j]*dw[0]
+                - ATAinvAT[1][j]*LSQw[j]*dw[1]
+                - ATAinvAT[2][j]*LSQw[j]*dw[2]
+                )
+            /W;
+            }
+        }
+    }
+
+    if (shortC)
+    {
+        // Matrix ill conditioned. Use straight injection from central
+        // donor.
+        weights = 0.0;
+        weights[0] = 1.0;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::cellCellStencils::leastSquares::leastSquares
+(
+    const fvMesh& mesh,
+    const dictionary& dict,
+    const bool doUpdate
+)
+:
+    inverseDistance(mesh, dict, false)
+{
+    if (doUpdate)
+    {
+        update();
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::cellCellStencils::leastSquares::~leastSquares()
+{}
+
+
+// ************************************************************************* //

src/overset/cellCellStencil/leastSquares/leastSquaresCellCellStencil.H

@@ -0,0 +1,107 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     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 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::cellCellStencils::leastSquares
+
+Description
+    Least-squares-weighted interpolation stencil.
+
+    Base machinery is similar to inverse distance interpolation stencil
+    but weights minimize error in LSQ sense recovering exact solution
+    for linear solution problems. Gradient and values are found
+    simultaneously.
+
+SourceFiles
+    leastSquaresCellCellStencil.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef cellCellStencils_leastSquares_H
+#define cellCellStencils_leastSquares_H
+
+#include "inverseDistanceCellCellStencil.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+namespace cellCellStencils
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class leastSquares Declaration
+\*---------------------------------------------------------------------------*/
+
+class leastSquares
+:
+    public inverseDistance
+{
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        leastSquares(const leastSquares&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const leastSquares&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("leastSquares");
+
+
+    // Constructors
+
+        //- Construct from fvMesh
+        leastSquares(const fvMesh&, const dictionary&, const bool);
+
+
+    //- Destructor
+    virtual ~leastSquares();
+
+
+    // Member Functions
+
+        //- Calculate lsq weights for single acceptor
+        virtual void stencilWeights
+        (
+            const point& sample,
+            const pointList& donorCcs,
+            scalarList& weights
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace cellCellStencils
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/overset/cellCellStencil/trackingInverseDistance/trackingInverseDistanceCellCellStencil.C

@@ -51,6 +51,7 @@ namespace cellCellStencils
 void Foam::cellCellStencils::trackingInverseDistance::markBoundaries
 (
     const fvMesh& mesh,
+    const vector& smallVec,
 
     const boundBox& bb,
     const labelVector& nDivs,
@@ -84,6 +85,9 @@ void Foam::cellCellStencils::trackingInverseDistance::markBoundaries
 
                 // Mark in voxel mesh
                 boundBox faceBb(pp.points(), pp[i], false);
+                faceBb.min() -= smallVec;
+                faceBb.max() += smallVec;
+
                 if (bb.overlaps(faceBb))
                 {
                     voxelMeshSearch::fill
@@ -116,6 +120,8 @@ void Foam::cellCellStencils::trackingInverseDistance::markBoundaries
 
                 // Mark in voxel mesh
                 boundBox faceBb(pp.points(), pp[i], false);
+                faceBb.min() -= smallVec;
+                faceBb.max() += smallVec;
                 if (bb.overlaps(faceBb))
                 {
                     voxelMeshSearch::fill
@@ -165,6 +171,9 @@ void Foam::cellCellStencils::trackingInverseDistance::markPatchesAsHoles
             forAll(tgtCellMap, tgtCelli)
             {
                 label celli = tgtCellMap[tgtCelli];
+                boundBox cBb(allPoints, allCellPoints[celli], false);
+                cBb.min() -= smallVec_;
+                cBb.max() += smallVec_;
 
                 if
                 (
@@ -172,7 +181,7 @@ void Foam::cellCellStencils::trackingInverseDistance::markPatchesAsHoles
                     (
                         srcPatchBb,
                         srcDivs,
-                        boundBox(allPoints, allCellPoints[celli], false),
+                        cBb,
                         srcPatchTypes,
                         static_cast<unsigned int>(patchCellType::PATCH)
                     )
@@ -231,13 +240,17 @@ void Foam::cellCellStencils::trackingInverseDistance::markPatchesAsHoles
                 forAll(tgtCellMap, tgtCelli)
                 {
                     label celli = tgtCellMap[tgtCelli];
+                    boundBox cBb(allPoints, allCellPoints[celli], false);
+                    cBb.min() -= smallVec_;
+                    cBb.max() += smallVec_;
+
                     if
                     (
                         voxelMeshSearch::overlaps
                         (
                             srcPatchBb,
                             srcDivs,
-                            boundBox(allPoints, allCellPoints[celli], false),
+                            cBb,
                             srcPatchTypes,
                             static_cast<unsigned int>(patchCellType::PATCH)
                         )
@@ -628,6 +641,7 @@ bool Foam::cellCellStencils::trackingInverseDistance::update()
         markBoundaries
         (
             meshParts_[zonei].subMesh(),
+            smallVec_,
 
             patchBb[zonei][Pstream::myProcNo()],
             patchDivisions[zonei],

src/overset/cellCellStencil/trackingInverseDistance/trackingInverseDistanceCellCellStencil.H

@@ -77,6 +77,7 @@ protected:
         static void markBoundaries
         (
             const fvMesh& mesh,
+            const vector& smallVec,
 
             const boundBox& bb,
             const labelVector& nDivs,

src/overset/dynamicOversetFvMesh/dynamicOversetFvMesh.C

@@ -228,7 +228,10 @@ Foam::dynamicOversetFvMesh::dynamicOversetFvMesh(const IOobject& io)
 :
     dynamicMotionSolverFvMesh(io),
     active_(false)
-{}
+{
+    // Force loading zoneID field before time gets incremented
+    (void)cellCellStencil::zoneID(*this);
+}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

src/overset/dynamicOversetFvMesh/dynamicOversetFvMesh.H

@@ -38,6 +38,8 @@ SourceFiles
 #include "dynamicMotionSolverFvMesh.H"
 #include "labelIOList.H"
 #include "fvMeshPrimitiveLduAddressing.H"
+#include "lduInterfaceFieldPtrsList.H"
+#include "volFields.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -100,6 +102,14 @@ protected:
         template<class GeoField>
         void interpolate(GeoField& psi) const;
 
+        //- Freeze values at holes
+        template<class Type>
+        void freezeHoles(fvMatrix<Type>&) const;
+
+        //- Get scalar interfaces of certain type
+        template<class GeoField, class PatchType>
+        lduInterfaceFieldPtrsList scalarInterfaces(const GeoField& psi) const;
+
         //- Add interpolation to matrix (coefficients)
         template<class Type>
         void addInterpolation(fvMatrix<Type>&) const;
@@ -112,25 +122,6 @@ protected:
         template<class GeoField>
         static void correctCoupledBoundaryConditions(GeoField& fld);
 
-        //- Use extended addressing
-        void active(const bool f) const
-        {
-            active_ = f;
-
-            if (active_)
-            {
-                DebugInfo<< "Switching to extended addressing with nFaces:"
-                    << primitiveLduAddr().lowerAddr().size()
-                    << endl;
-            }
-            else
-            {
-                DebugInfo<< "Switching to base addressing with nFaces:"
-                    << fvMesh::lduAddr().lowerAddr().size()
-                    << endl;
-            }
-        }
-
 
 private:
 
@@ -164,12 +155,6 @@ public:
 
         // Extended addressing
 
-            //- Use extended addressing
-            bool active() const
-            {
-                return active_;
-            }
-
             //- Return extended ldu addressing
             const fvMeshPrimitiveLduAddressing& primitiveLduAddr() const;
 
@@ -177,6 +162,37 @@ public:
             //  primitiveLduAddr
             virtual const lduAddressing& lduAddr() const;
 
+            //- Return old to new face addressing
+            const labelList& reverseFaceMap() const
+            {
+                return reverseFaceMap_;
+            }
+
+            //- Return true if using extended addressing
+            bool active() const
+            {
+                return active_;
+            }
+
+            //- Enable/disable extended addressing
+            void active(const bool f) const
+            {
+                active_ = f;
+
+                if (active_)
+                {
+                    DebugInfo<< "Switching to extended addressing with nFaces:"
+                        << primitiveLduAddr().lowerAddr().size()
+                        << endl;
+                }
+                else
+                {
+                    DebugInfo<< "Switching to base addressing with nFaces:"
+                        << fvMesh::lduAddr().lowerAddr().size()
+                        << endl;
+                }
+            }
+
 
         // Overset
 

src/overset/dynamicOversetFvMesh/dynamicOversetFvMeshTemplates.C

@@ -26,6 +26,7 @@ License
 #include "volFields.H"
 #include "fvMatrix.H"
 #include "cellCellStencilObject.H"
+#include "oversetFvPatchField.H"
 
 // * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * * //
 
@@ -77,41 +78,34 @@ void Foam::dynamicOversetFvMesh::interpolate(GeoField& psi) const
 }
 
 
-//template<class T>
+template<class GeoField, class PatchType>
-//void Foam::dynamicOversetFvMesh::setReference(fvMatrix<T>& m) const
+Foam::lduInterfaceFieldPtrsList
-//{
+Foam::dynamicOversetFvMesh::scalarInterfaces(const GeoField& psi) const
-//    const mapDistribute& map = cellInterpolationMap();
+{
-//    const List<scalarList>& wghts = cellInterpolationWeights();
+    const typename GeoField::Boundary& bpsi = psi.boundaryField();
-//    const labelListList& stencil = cellStencil();
+
-//    const labelList& cellIDs = interpolationCells();
+    lduInterfaceFieldPtrsList interfaces(bpsi.size());
-//    const scalarList& factor = cellInterpolationWeight();
+
-//
+    forAll(interfaces, patchi)
-//    Field<T> work(m.psi());
+    {
-//    map.mapDistributeBase::distribute(work, UPstream::msgType()+1);
+        if
-//
+        (
-//    forAll(cellIDs, i)
+            isA<lduInterfaceField>(bpsi[patchi])
-//    {
+         && isA<PatchType>(bpsi[patchi])
-//        label celli = cellIDs[i];
+        )
-//
+        {
-//        const scalarList& w = wghts[celli];
+            interfaces.set
-//        const labelList& nbrs = stencil[celli];
+            (
-//        const scalar f = factor[celli];
+                patchi,
-//
+                &refCast<const lduInterfaceField>
-//        //Pout<< "Interpolating " << celli << " from values:" << endl;
+                (
-//        T s(pTraits<T>::zero);
+                    bpsi[patchi]
-//        forAll(nbrs, nbrI)
+                )
-//        {
+            );
-//            //Pout<< "    " << work[nbrs[nbrI]]
+        }
-//            //    << " from slot " << nbrs[nbrI] << endl;
+    }
-//            s += w[nbrI]*work[nbrs[nbrI]];
+    return interfaces;
-//        }
+}
-//        //Pout<< "Interpolated value:" << s << endl;
-//        const T oldPsi = m.psi()[celli];
-//        m.setReference(celli, (1.0-f)*oldPsi + f*s, true);
-//        Pout<< "psi was:" << oldPsi << " now forcing to:"
-//            << (1.0-f)*oldPsi + f*s << endl;
-//    }
-//}
 
 
 template<class Type>
@@ -185,6 +179,7 @@ void Foam::dynamicOversetFvMesh::addInterpolation(fvMatrix<Type>& m) const
         }
     }
 
+
     forAll(m.internalCoeffs(), patchI)
     {
         const labelUList& fc = addr.patchAddr(patchI);
@@ -209,6 +204,76 @@ void Foam::dynamicOversetFvMesh::addInterpolation(fvMatrix<Type>& m) const
         }
     }
 
+    // NOTE: For a non-scalar matrix, in the function fvMatrix::solveSegregated
+    // it uses updateInterfaceMatrix to correct for remote source contributions.
+    // This unfortunately also triggers the overset interpolation which then
+    // produces a non-zero source for interpolated cells.
+    // The procedure below calculates this contribution 'correctionSource'
+    // and substracts it from the source later in order to compensate.
+    Field<Type> correctionSource(diag.size(), pTraits<Type>::zero);
+
+    if (pTraits<Type>::nComponents > 1)
+    {
+        typedef GeometricField<Type, fvPatchField, volMesh> GeoField;
+
+        typename pTraits<Type>::labelType validComponents
+        (
+            this->template validComponents<Type>()
+        );
+
+        // Get all the overset lduInterfaces
+        lduInterfaceFieldPtrsList interfaces
+        (
+            scalarInterfaces<GeoField, oversetFvPatchField<Type>>
+            (
+                m.psi()
+            )
+        );
+
+        const Field<Type>& psiInt = m.psi().primitiveField();
+
+        for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
+        {
+            if (component(validComponents, cmpt) != -1)
+            {
+                const scalarField psiCmpt(psiInt.component(cmpt));
+                scalarField corrSourceCmpt(correctionSource.component(cmpt));
+
+                forAll(interfaces, patchi)
+                {
+                    if (interfaces.set(patchi))
+                    {
+                        interfaces[patchi].initInterfaceMatrixUpdate
+                        (
+                            corrSourceCmpt,
+                            true,
+                            psiCmpt,
+                            m.boundaryCoeffs()[patchi].component(cmpt),
+                            cmpt,
+                            Pstream::defaultCommsType
+                        );
+                    }
+                }
+                forAll(interfaces, patchi)
+                {
+                    if (interfaces.set(patchi))
+                    {
+                        interfaces[patchi].updateInterfaceMatrix
+                        (
+                            corrSourceCmpt,
+                            true,
+                            psiCmpt,
+                            m.boundaryCoeffs()[patchi].component(cmpt),
+                            cmpt,
+                            Pstream::defaultCommsType
+                        );
+                    }
+                }
+                correctionSource.replace(cmpt, corrSourceCmpt);
+            }
+        }
+    }
+
 
 
     // Modify matrix
@@ -271,12 +336,14 @@ void Foam::dynamicOversetFvMesh::addInterpolation(fvMatrix<Type>& m) const
         {
             // Create interpolation stencil. Leave any non-local contributions
             //  to be done by oversetFvPatchField updateMatrixInterface
+            // 'correctionSource' is only applied for vector and higher
+            // fvMatrix; it is zero for scalar matrices (see above).
 
             diag[celli] *= (1.0-f);
             diag[celli] += normalisation*f;
 
             source[celli] *= (1.0-f);
-            source[celli] += normalisation*f*pTraits<Type>::zero; // dummy
+            source[celli] -= correctionSource[celli];
 
             //Pout<< "Interpolating " << celli << " with factor:" << f
             //    << " new diag:" << diag[celli]
@@ -328,6 +395,39 @@ Foam::SolverPerformance<Type> Foam::dynamicOversetFvMesh::solve
     const dictionary& dict
 ) const
 {
+    // Check we're running with bcs that can handle implicit matrix manipulation
+    const typename GeometricField<Type, fvPatchField, volMesh>::Boundary bpsi =
+        m.psi().boundaryField();
+
+    bool hasOverset = false;
+    forAll(bpsi, patchi)
+    {
+        if (isA<oversetFvPatchField<Type>>(bpsi[patchi]))
+        {
+            hasOverset = true;
+            break;
+        }
+    }
+
+    if (!hasOverset)
+    {
+        if (debug)
+        {
+            Pout<< "dynamicOversetFvMesh::solve() :"
+                << " bypassing matrix adjustment for field " << m.psi().name()
+                << endl;
+        }
+        return dynamicMotionSolverFvMesh::solve(m, dict);
+    }
+
+    if (debug)
+    {
+        Pout<< "dynamicOversetFvMesh::solve() :"
+            << " adjusting matrix for interpolation for field "
+            << m.psi().name() << endl;
+    }
+
+
     // Switch to extended addressing (requires mesh::update() having been
     // called)
     active(true);
@@ -346,8 +446,7 @@ Foam::SolverPerformance<Type> Foam::dynamicOversetFvMesh::solve
     // Use lower level solver
     SolverPerformance<Type> s(dynamicMotionSolverFvMesh::solve(m, dict));
 
-    // Reset matrix (shuffle faces into original order)
+    // Restore matrix
-    //resetMatrix(m);
     m.upper().transfer(oldUpper);
     m.lower().transfer(oldLower);
     m.internalCoeffs().transfer(oldInt);

src/overset/oversetPolyPatch/oversetFvPatchField.C

@@ -24,6 +24,8 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "volFields.H"
+#include "cellCellStencil.H"
+#include "dynamicOversetFvMesh.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -34,9 +36,7 @@ Foam::oversetFvPatchField<Type>::oversetFvPatchField
     const DimensionedField<Type, volMesh>& iF
 )
 :
-    LduInterfaceField<Type>(refCast<const oversetFvPatch>(p)),
+    semiImplicitOversetFvPatchField<Type>(p, iF)
-    zeroGradientFvPatchField<Type>(p, iF),
-    oversetPatch_(refCast<const oversetFvPatch>(p))
 {}
 
 
@@ -49,21 +49,8 @@ Foam::oversetFvPatchField<Type>::oversetFvPatchField
     const fvPatchFieldMapper& mapper
 )
 :
-    LduInterfaceField<Type>(refCast<const oversetFvPatch>(p)),
+    semiImplicitOversetFvPatchField<Type>(ptf, p, iF, mapper)
-    zeroGradientFvPatchField<Type>(ptf, p, iF, mapper),
+{}
-    oversetPatch_(refCast<const oversetFvPatch>(p))
-{
-    if (!isA<oversetFvPatch>(this->patch()))
-    {
-        FatalErrorInFunction
-            << "    patch type '" << p.type()
-            << "' not constraint type '" << typeName << "'"
-            << "\n    for patch " << p.name()
-            << " of field " << this->internalField().name()
-            << " in file " << this->internalField().objectPath()
-            << exit(FatalIOError);
-    }
-}
 
 
 template<class Type>
@@ -74,26 +61,8 @@ Foam::oversetFvPatchField<Type>::oversetFvPatchField
     const dictionary& dict
 )
 :
-    LduInterfaceField<Type>(refCast<const oversetFvPatch>(p)),
+    semiImplicitOversetFvPatchField<Type>(p, iF, dict)
-    zeroGradientFvPatchField<Type>(p, iF, dict),
+{}
-    oversetPatch_(refCast<const oversetFvPatch>(p))
-{
-    if (!isA<oversetFvPatch>(p))
-    {
-        FatalIOErrorInFunction(dict)
-            << "    patch type '" << p.type()
-            << "' not constraint type '" << typeName << "'"
-            << "\n    for patch " << p.name()
-            << " of field " << this->internalField().name()
-            << " in file " << this->internalField().objectPath()
-            << exit(FatalIOError);
-    }
-
-    if (!dict.found("value") && this->coupled())
-    {
-        this->evaluate(Pstream::commsTypes::blocking);
-    }
-}
 
 
 template<class Type>
@@ -102,9 +71,7 @@ Foam::oversetFvPatchField<Type>::oversetFvPatchField
     const oversetFvPatchField<Type>& ptf
 )
 :
-    LduInterfaceField<Type>(ptf.oversetPatch_),
+    semiImplicitOversetFvPatchField<Type>(ptf)
-    zeroGradientFvPatchField<Type>(ptf),
-    oversetPatch_(ptf.oversetPatch_)
 {}
 
 
@@ -115,22 +82,19 @@ Foam::oversetFvPatchField<Type>::oversetFvPatchField
     const DimensionedField<Type, volMesh>& iF
 )
 :
-    LduInterfaceField<Type>(ptf.oversetPatch_),
+    semiImplicitOversetFvPatchField<Type>(ptf, iF)
-    zeroGradientFvPatchField<Type>(ptf, iF),
-    oversetPatch_(ptf.oversetPatch_)
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-
 template<class Type>
 void Foam::oversetFvPatchField<Type>::initEvaluate
 (
     const Pstream::commsTypes commsType
 )
 {
-    if (oversetPatch_.master())
+    if (this->oversetPatch_.master())
     {
         // Trigger interpolation
         const fvMesh& mesh = this->internalField().mesh();
@@ -139,24 +103,14 @@ void Foam::oversetFvPatchField<Type>::initEvaluate
 
         if (&mesh.lduAddr() != &mesh.fvMesh::lduAddr())
         {
-            // Running extended addressing. Interpolate always.
+            // Running extended addressing. Flux correction already done
+            // in the linear solver (through the initUpdateInterfaceMatrix
+            // method below)
             if (debug)
             {
-                Info<< "Interpolating solved-for field " << fldName << endl;
+                Info<< "Skipping overset interpolation for solved-for field "
+                    << fldName << endl;
             }
-
-            // Interpolate without bc update (since would trigger infinite
-            // recursion back to oversetFvPatchField<Type>::evaluate)
-            // The only problem is bcs that use the new cell values
-            // (e.g. zeroGradient, processor). These need to appear -after-
-            // the 'overset' bc.
-            mesh.interpolate
-            (
-                const_cast<Field<Type>&>
-                (
-                    this->primitiveField()
-                )
-            );
         }
         else if
         (
@@ -206,31 +160,6 @@ void Foam::oversetFvPatchField<Type>::initEvaluate
 }
 
 
-template<class Type>
-void Foam::oversetFvPatchField<Type>::evaluate(const Pstream::commsTypes)
-{
-    if (!this->updated())
-    {
-        this->updateCoeffs();
-    }
-
-    zeroGradientFvPatchField<Type>::evaluate();
-}
-
-
-template<class Type>
-void Foam::oversetFvPatchField<Type>::initInterfaceMatrixUpdate
-(
-    scalarField& result,
-    const bool add,
-    const scalarField& psiInternal,
-    const scalarField& coeffs,
-    const direction cmpt,
-    const Pstream::commsTypes commsType
-) const
-{}
-
-
 template<class Type>
 void Foam::oversetFvPatchField<Type>::updateInterfaceMatrix
 (
@@ -243,7 +172,10 @@ void Foam::oversetFvPatchField<Type>::updateInterfaceMatrix
 ) const
 {
     // Add remote values
-    if (oversetPatch_.master())
+
+    const oversetFvPatch& ovp = this->oversetPatch_;
+
+    if (ovp.master())
     {
         const fvMesh& mesh = this->patch().boundaryMesh().mesh();
 
@@ -251,11 +183,10 @@ void Foam::oversetFvPatchField<Type>::updateInterfaceMatrix
         // TBD. This should be cleaner.
         if (&mesh.lduAddr() == &mesh.fvMesh::lduAddr())
         {
-            //Pout<< "** not running extended addressing..." << endl;
             return;
         }
 
-        const labelListList& stencil = oversetPatch_.stencil();
+        const labelListList& stencil = ovp.stencil();
 
         if (stencil.size() != psiInternal.size())
         {
@@ -263,12 +194,11 @@ void Foam::oversetFvPatchField<Type>::updateInterfaceMatrix
                 << " stencil:" << stencil.size() << exit(FatalError);
         }
 
-        const mapDistribute& map = oversetPatch_.cellInterpolationMap();
+        const mapDistribute& map = ovp.cellInterpolationMap();
-        const List<scalarList>& wghts =
+        const List<scalarList>& wghts = ovp.cellInterpolationWeights();
-            oversetPatch_.cellInterpolationWeights();
+        const labelList& cellIDs = ovp.interpolationCells();
-        const labelList& cellIDs = oversetPatch_.interpolationCells();
+        const scalarList& factor = ovp.cellInterpolationWeight();
-        const scalarList& factor = oversetPatch_.cellInterpolationWeight();
+        const scalarField& normalisation = ovp.normalisation();
-        const scalarField& normalisation = oversetPatch_.normalisation();
 
 
         // Since we're inside initEvaluate/evaluate there might be processor
@@ -306,41 +236,4 @@ void Foam::oversetFvPatchField<Type>::updateInterfaceMatrix
 }
 
 
-template<class Type>
-void Foam::oversetFvPatchField<Type>::initInterfaceMatrixUpdate
-(
-    Field<Type>&,
-    const bool add,
-    const Field<Type>&,
-    const scalarField&,
-    const Pstream::commsTypes commsType
-) const
-{
-    NotImplemented;
-}
-
-
-template<class Type>
-void Foam::oversetFvPatchField<Type>::updateInterfaceMatrix
-(
-    Field<Type>& result,
-    const bool add,
-    const Field<Type>& psiInternal,
-    const scalarField&,
-    const Pstream::commsTypes
-) const
-{
-    NotImplemented;
-}
-
-
-template<class Type>
-void Foam::oversetFvPatchField<Type>::write(Ostream& os) const
-{
-    zeroGradientFvPatchField<Type>::write(os);
-    // Make sure to write the value for ease of postprocessing.
-    this->writeEntry("value", os);
-}
-
-
 // ************************************************************************* //

src/overset/oversetPolyPatch/oversetFvPatchField.H

@@ -28,8 +28,9 @@ Group
     grpCoupledBoundaryConditions
 
 Description
-
+    Boundary condition for use on overset patches. To be run in combination
-SeeAlso
+    with special dynamicFvMesh type that inserts local contributions into
+    the matrix. This boundary condition adds the remote contributions.
 
 SourceFiles
     oversetFvPatchField.C
@@ -39,9 +40,7 @@ SourceFiles
 #ifndef oversetFvPatchField_H
 #define oversetFvPatchField_H
 
-#include "oversetFvPatch.H"
+#include "semiImplicitOversetFvPatchField.H"
-#include "zeroGradientFvPatchField.H"
-#include "LduInterfaceField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -55,18 +54,8 @@ namespace Foam
 template<class Type>
 class oversetFvPatchField
 :
-    public LduInterfaceField<Type>,
+    public semiImplicitOversetFvPatchField<Type>
-    public zeroGradientFvPatchField<Type>
 {
-    // Private data
-
-        //- Local reference cast into the overset patch
-        const oversetFvPatch& oversetPatch_;
-
-        //- Master patch ID
-        mutable label masterPatchID_;
-
-
 public:
 
     //- Runtime type information
@@ -133,34 +122,11 @@ public:
 
     // Member functions
 
-        // Access
-
-            //- Return local reference cast into the cyclic AMI patch
-            const oversetFvPatch& oversetPatch() const
-            {
-                return oversetPatch_;
-            }
-
-
         // Evaluation functions
 
             //- Initialise the evaluation of the patch field
             virtual void initEvaluate(const Pstream::commsTypes commsType);
 
-            //- Evaluate the patch field
-            virtual void evaluate(const Pstream::commsTypes commsType);
-
-            //- Initialise neighbour matrix update
-            virtual void initInterfaceMatrixUpdate
-            (
-                scalarField&,
-                const bool add,
-                const scalarField&,
-                const scalarField&,
-                const direction,
-                const Pstream::commsTypes commsType
-            ) const;
-
             //- Update result field based on interface functionality
             virtual void updateInterfaceMatrix
             (
@@ -171,32 +137,6 @@ public:
                 const direction cmpt,
                 const Pstream::commsTypes commsType
             ) const;
-
-            //- Initialise neighbour matrix update
-            virtual void initInterfaceMatrixUpdate
-            (
-                Field<Type>&,
-                const bool add,
-                const Field<Type>&,
-                const scalarField&,
-                const Pstream::commsTypes commsType
-            ) const;
-
-            //- Update result field based on interface functionality
-            virtual void updateInterfaceMatrix
-            (
-                Field<Type>&,
-                const bool add,
-                const Field<Type>&,
-                const scalarField&,
-                const Pstream::commsTypes commsType
-            ) const;
-
-
-        // I-O
-
-            //- Write
-            virtual void write(Ostream& os) const;
 };
 
 

src/overset/oversetPolyPatch/oversetGAMGInterfaceField.H

@@ -52,12 +52,16 @@ class oversetGAMGInterfaceField
 :
     public GAMGInterfaceField
 {
-    // Private data
+protected:
+
+    // Protected data
 
         //- Local reference cast into the interface
         const oversetGAMGInterface& oversetInterface_;
 
 
+private:
+
     // Private Member Functions
 
         //- Disallow default bitwise copy construct

src/overset/oversetPolyPatch/semiImplicitOversetFvPatchField.C

@@ -0,0 +1,272 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     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 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 "volFields.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::semiImplicitOversetFvPatchField<Type>::semiImplicitOversetFvPatchField
+(
+    const fvPatch& p,
+    const DimensionedField<Type, volMesh>& iF
+)
+:
+    LduInterfaceField<Type>(refCast<const oversetFvPatch>(p)),
+    zeroGradientFvPatchField<Type>(p, iF),
+    oversetPatch_(refCast<const oversetFvPatch>(p))
+{}
+
+
+template<class Type>
+Foam::semiImplicitOversetFvPatchField<Type>::semiImplicitOversetFvPatchField
+(
+    const semiImplicitOversetFvPatchField<Type>& ptf,
+    const fvPatch& p,
+    const DimensionedField<Type, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    LduInterfaceField<Type>(refCast<const oversetFvPatch>(p)),
+    zeroGradientFvPatchField<Type>(ptf, p, iF, mapper),
+    oversetPatch_(refCast<const oversetFvPatch>(p))
+{
+    if (!isA<oversetFvPatch>(this->patch()))
+    {
+        FatalErrorInFunction
+            << "    patch type '" << p.type()
+            << "' not constraint type '" << typeName << "'"
+            << "\n    for patch " << p.name()
+            << " of field " << this->internalField().name()
+            << " in file " << this->internalField().objectPath()
+            << exit(FatalIOError);
+    }
+}
+
+
+template<class Type>
+Foam::semiImplicitOversetFvPatchField<Type>::semiImplicitOversetFvPatchField
+(
+    const fvPatch& p,
+    const DimensionedField<Type, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    LduInterfaceField<Type>(refCast<const oversetFvPatch>(p)),
+    zeroGradientFvPatchField<Type>(p, iF, dict),
+    oversetPatch_(refCast<const oversetFvPatch>(p))
+{
+    if (!isA<oversetFvPatch>(p))
+    {
+        FatalIOErrorInFunction(dict)
+            << "    patch type '" << p.type()
+            << "' not constraint type '" << typeName << "'"
+            << "\n    for patch " << p.name()
+            << " of field " << this->internalField().name()
+            << " in file " << this->internalField().objectPath()
+            << exit(FatalIOError);
+    }
+
+    if (!dict.found("value") && this->coupled())
+    {
+        this->evaluate(Pstream::commsTypes::blocking);
+    }
+}
+
+
+template<class Type>
+Foam::semiImplicitOversetFvPatchField<Type>::semiImplicitOversetFvPatchField
+(
+    const semiImplicitOversetFvPatchField<Type>& ptf
+)
+:
+    LduInterfaceField<Type>(ptf.oversetPatch_),
+    zeroGradientFvPatchField<Type>(ptf),
+    oversetPatch_(ptf.oversetPatch_)
+{}
+
+
+template<class Type>
+Foam::semiImplicitOversetFvPatchField<Type>::semiImplicitOversetFvPatchField
+(
+    const semiImplicitOversetFvPatchField<Type>& ptf,
+    const DimensionedField<Type, volMesh>& iF
+)
+:
+    LduInterfaceField<Type>(ptf.oversetPatch_),
+    zeroGradientFvPatchField<Type>(ptf, iF),
+    oversetPatch_(ptf.oversetPatch_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::semiImplicitOversetFvPatchField<Type>::evaluate
+(
+    const Pstream::commsTypes
+)
+{
+    if (debug)
+    {
+        Pout<< FUNCTION_NAME << " field " <<  this->internalField().name()
+            << " patch " << this->patch().name() << endl;
+    }
+
+
+    if (!this->updated())
+    {
+        this->updateCoeffs();
+    }
+
+    zeroGradientFvPatchField<Type>::evaluate();
+}
+
+
+template<class Type>
+void Foam::semiImplicitOversetFvPatchField<Type>::initEvaluate
+(
+    const Pstream::commsTypes commsType
+)
+{
+    if (debug)
+    {
+        Pout<< FUNCTION_NAME << " field " <<  this->internalField().name()
+            << " patch " << this->patch().name() << endl;
+    }
+
+    if (this->oversetPatch().master())
+    {
+        // Trigger interpolation
+        const fvMesh& mesh = this->internalField().mesh();
+        //const dictionary& fvSchemes = mesh.schemesDict();
+        const word& fldName = this->internalField().name();
+
+        if (&mesh.lduAddr() != &mesh.fvMesh::lduAddr())
+        {
+            // Running extended addressing so called from within fvMatrix::solve
+            FatalErrorInFunction
+                << "Running extended addressing is not allowed when solving "
+                << fldName
+                << " Please choose a dynamicFvMesh without matrix adaptation"
+                << exit(FatalError);
+        }
+        else
+        {
+            if (debug)
+            {
+                Info<< "Interpolating field " << fldName << endl;
+            }
+
+            // Interpolate without bc update (since would trigger infinite
+            // recursion back to
+            // semiImplicitOversetFvPatchField<Type>::evaluate)
+            // The only problem is bcs that use the new cell values
+            // (e.g. zeroGradient, processor). These need to appear -after-
+            // the 'overset' bc.
+            mesh.interpolate
+            (
+                const_cast<Field<Type>&>
+                (
+                    this->primitiveField()
+                )
+            );
+        }
+    }
+}
+
+
+template<class Type>
+void Foam::semiImplicitOversetFvPatchField<Type>::updateInterfaceMatrix
+(
+    scalarField& result,
+    const bool add,
+    const scalarField& psiInternal,
+    const scalarField& coeffs,
+    const direction cmpt,
+    const Pstream::commsTypes
+) const
+{
+    if (debug)
+    {
+        Pout<< FUNCTION_NAME << " field " <<  this->internalField().name()
+            << " patch " << this->patch().name() << endl;
+    }
+
+    const oversetFvPatch& ovp = this->oversetPatch();
+
+    // Set all interpolated cells
+    if (ovp.master())
+    {
+        const labelListList& stencil = ovp.stencil();
+
+        if (stencil.size() != psiInternal.size())
+        {
+            FatalErrorInFunction << "psiInternal:" << psiInternal.size()
+                << " stencil:" << stencil.size() << exit(FatalError);
+        }
+
+        const mapDistribute& map = ovp.cellInterpolationMap();
+        const List<scalarList>& wghts = ovp.cellInterpolationWeights();
+        const labelList& cellIDs = ovp.interpolationCells();
+        //const scalarList& factor = ovp.cellInterpolationWeight();
+
+        // Since we're inside initEvaluate/evaluate there might be processor
+        // comms underway. Change the tag we use.
+        scalarField work(psiInternal);
+        map.mapDistributeBase::distribute(work, UPstream::msgType()+1);
+
+        forAll(cellIDs, i)
+        {
+            label celli = cellIDs[i];
+
+            const scalarList& w = wghts[celli];
+            const labelList& nbrs = stencil[celli];
+
+            //scalar f = factor[celli];
+
+            scalar s(0.0);
+            forAll(nbrs, nbrI)
+            {
+                s += w[nbrI]*work[nbrs[nbrI]];
+            }
+
+            //Pout<< "cell:" << celli << " interpolated value:" << s << endl;
+            result[celli] = s;  //(1.0-f)*result[celli] + f*s;
+        }
+    }
+}
+
+
+template<class Type>
+void Foam::semiImplicitOversetFvPatchField<Type>::write(Ostream& os) const
+{
+    zeroGradientFvPatchField<Type>::write(os);
+    // Make sure to write the value for ease of postprocessing.
+    this->writeEntry("value", os);
+}
+
+
+// ************************************************************************* //

src/overset/oversetPolyPatch/semiImplicitOversetFvPatchField.H

@@ -0,0 +1,207 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     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 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::semiImplicitOversetFvPatchField
+
+Group
+    grpCoupledBoundaryConditions
+
+Description
+    Boundary condition for use on overset patches. Implements update of
+    interpolated cells inside the linear solvers.
+
+SeeAlso
+
+SourceFiles
+    semiImplicitOversetFvPatchField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef semiImplicitOversetFvPatchField_H
+#define semiImplicitOversetFvPatchField_H
+
+#include "oversetFvPatch.H"
+#include "zeroGradientFvPatchField.H"
+#include "LduInterfaceField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                 Class semiImplicitOversetFvPatchField Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class semiImplicitOversetFvPatchField
+:
+    public LduInterfaceField<Type>,
+    public zeroGradientFvPatchField<Type>
+{
+protected:
+
+    // Protected data
+
+        //- Local reference cast into the overset patch
+        const oversetFvPatch& oversetPatch_;
+
+        //- Master patch ID
+        mutable label masterPatchID_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("semiImplicitOverset");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        semiImplicitOversetFvPatchField
+        (
+            const fvPatch&,
+            const DimensionedField<Type, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        semiImplicitOversetFvPatchField
+        (
+            const fvPatch&,
+            const DimensionedField<Type, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given semiImplicitOversetFvPatchField onto a
+        //  new patch
+        semiImplicitOversetFvPatchField
+        (
+            const semiImplicitOversetFvPatchField<Type>&,
+            const fvPatch&,
+            const DimensionedField<Type, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        semiImplicitOversetFvPatchField
+        (
+            const semiImplicitOversetFvPatchField<Type>&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchField<Type> > clone() const
+        {
+            return tmp<fvPatchField<Type> >
+            (
+                new semiImplicitOversetFvPatchField<Type>(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        semiImplicitOversetFvPatchField
+        (
+            const semiImplicitOversetFvPatchField<Type>&,
+            const DimensionedField<Type, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchField<Type> > clone
+        (
+            const DimensionedField<Type, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchField<Type> >
+            (
+                new semiImplicitOversetFvPatchField<Type>(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        // Access
+
+            //- Return local reference cast into the cyclic AMI patch
+            const oversetFvPatch& oversetPatch() const
+            {
+                return oversetPatch_;
+            }
+
+
+        // Evaluation functions
+
+            //- Evaluate the patch field
+            virtual void evaluate(const Pstream::commsTypes commsType);
+
+             //- Initialise the evaluation of the patch field
+            virtual void initEvaluate(const Pstream::commsTypes commsType);
+
+            //- Update result field based on interface functionality
+            virtual void updateInterfaceMatrix
+            (
+                scalarField& result,
+                const bool add,
+                const scalarField& psiInternal,
+                const scalarField& coeffs,
+                const direction cmpt,
+                const Pstream::commsTypes commsType
+            ) const;
+
+            //- Update result field based on interface functionality
+            virtual void updateInterfaceMatrix
+            (
+                Field<Type>&,
+                const bool add,
+                const Field<Type>&,
+                const scalarField&,
+                const Pstream::commsTypes commsType
+            ) const
+            {
+                NotImplemented;
+            }
+
+
+        // I-O
+
+            //- Write
+            virtual void write(Ostream& os) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "semiImplicitOversetFvPatchField.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/overset/oversetPolyPatch/semiImplicitOversetFvPatchFields.C

@@ -0,0 +1,43 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016-2017 OpenCFD Ltd.
+     \\/     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 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 "semiImplicitOversetFvPatchFields.H"
+#include "addToRunTimeSelectionTable.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+makePatchFields(semiImplicitOverset);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/overset/oversetPolyPatch/semiImplicitOversetFvPatchFields.H

@@ -0,0 +1,49 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     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 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef semiImplicitOversetFvPatchFields_H
+#define semiImplicitOversetFvPatchFields_H
+
+#include "semiImplicitOversetFvPatchField.H"
+#include "fieldTypes.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeFieldTypedefs(semiImplicitOverset);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/overset/oversetPolyPatch/semiImplicitOversetGAMGInterfaceField.C

@@ -0,0 +1,122 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     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 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 "semiImplicitOversetGAMGInterfaceField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "lduMatrix.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(semiImplicitOversetGAMGInterfaceField, 0);
+    addToRunTimeSelectionTable
+    (
+        GAMGInterfaceField,
+        semiImplicitOversetGAMGInterfaceField,
+        lduInterfaceField
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::semiImplicitOversetGAMGInterfaceField::
+semiImplicitOversetGAMGInterfaceField
+(
+    const GAMGInterface& GAMGCp,
+    const lduInterfaceField& fineInterface
+)
+:
+    oversetGAMGInterfaceField(GAMGCp, fineInterface)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::semiImplicitOversetGAMGInterfaceField::
+~semiImplicitOversetGAMGInterfaceField()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::semiImplicitOversetGAMGInterfaceField::updateInterfaceMatrix
+(
+    scalarField& result,
+    const bool add,
+    const scalarField& psiInternal,
+    const scalarField& coeffs,
+    const direction cmpt,
+    const Pstream::commsTypes commsType
+) const
+{
+    // Add remote values
+    if (oversetInterface_.master())
+    {
+        const labelListList& stencil = oversetInterface_.stencil();
+
+        if (stencil.size() != psiInternal.size())
+        {
+            FatalErrorInFunction << "psiInternal:" << psiInternal.size()
+                << " stencil:" << stencil.size() << exit(FatalError);
+        }
+
+        const mapDistribute& map = oversetInterface_.cellInterpolationMap();
+        const List<scalarList>& wghts =
+            oversetInterface_.cellInterpolationWeights();
+        const labelList& cellIDs = oversetInterface_.interpolationCells();
+        const scalarList& factor = oversetInterface_.cellInterpolationWeight();
+        const scalarField& normalisation = oversetInterface_.normalisation();
+
+        scalarField work(psiInternal);
+        map.mapDistributeBase::distribute(work, UPstream::msgType()+1);
+
+        forAll(cellIDs, i)
+        {
+            label celli = cellIDs[i];
+
+            const scalarList& w = wghts[celli];
+            const labelList& nbrs = stencil[celli];
+            const scalar f = factor[celli];
+
+            scalar s(0.0);
+            forAll(nbrs, nbrI)
+            {
+                label slotI = nbrs[nbrI];
+                s += w[nbrI]*work[slotI];
+            }
+            if (add)
+            {
+                s = -1.0*s;
+            }
+
+            result[celli] += normalisation[celli]*f*s;
+        }
+    }
+}
+
+
+// ************************************************************************* //

src/overset/oversetPolyPatch/semiImplicitOversetGAMGInterfaceField.H

@@ -0,0 +1,108 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     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 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::semiImplicitOversetGAMGInterfaceField
+
+Description
+    GAMG agglomerated processor interface field.
+
+SourceFiles
+    processorGAMGInterfaceField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef semiImplicitOversetGAMGInterfaceField_H
+#define semiImplicitOversetGAMGInterfaceField_H
+
+#include "oversetGAMGInterfaceField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+              Class semiImplicitOversetGAMGInterfaceField Declaration
+\*---------------------------------------------------------------------------*/
+
+class semiImplicitOversetGAMGInterfaceField
+:
+    public oversetGAMGInterfaceField
+{
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        semiImplicitOversetGAMGInterfaceField
+        (
+            const semiImplicitOversetGAMGInterfaceField&
+        );
+
+        //- Disallow default bitwise assignment
+        void operator=(const semiImplicitOversetGAMGInterfaceField&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("semiImplicitOverset");
+
+
+    // Constructors
+
+        //- Construct from GAMG interface and fine level interface field
+        semiImplicitOversetGAMGInterfaceField
+        (
+            const GAMGInterface& GAMGCp,
+            const lduInterfaceField& fineInterface
+        );
+
+
+    //- Destructor
+    virtual ~semiImplicitOversetGAMGInterfaceField();
+
+
+    // Member Functions
+
+        //- Update result field based on interface functionality
+        virtual void updateInterfaceMatrix
+        (
+            scalarField& result,
+            const bool add,
+            const scalarField& psiInternal,
+            const scalarField& coeffs,
+            const direction cmpt,
+            const Pstream::commsTypes commsType
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

tutorials/incompressible/overPimpleDyMFoam/simpleRotor/system/fvSchemes

@@ -57,6 +57,9 @@ snGradSchemes
 oversetInterpolation
 {
     method              inverseDistance;
+
+    searchBox           (0 0 0)(0.01 0.01 0.01);
+    searchBoxDivisions  (100 100 1);
 }
 
 fluxRequired

tutorials/incompressible/overSimpleFoam/aeroFoil/background_overset/system/decomposeParDict

@@ -16,11 +16,11 @@ FoamFile
 
 numberOfSubdomains 8;
 
-method          hierarchical;
+method          scotch; //hierarchical;
 
 coeffs
 {
-    n           (8 1 1);
+    n           (4 2 1);
     //delta       0.001; // default=0.001
     //order       xyz;   // default=xzy
 }

tutorials/incompressible/overSimpleFoam/aeroFoil/background_overset/system/fvSolution

@@ -60,6 +60,7 @@ solvers
 
 SIMPLE
 {
+    momentumPredictor   true;
     nOuterCorrectors    1;
     nCorrectors         2;
     nNonOrthogonalCorrectors 0;

tutorials/multiphase/overInterDyMFoam/floatingBody/background/system/fvSolution

@@ -29,8 +29,8 @@ solvers
 
     "alpha.water.*"
     {
-        nAlphaCorr      2;
+        nAlphaCorr      3;
-        nAlphaSubCycles 1;
+        nAlphaSubCycles 2;
         cAlpha          1;
         icAlpha         0;
 
@@ -84,6 +84,10 @@ PIMPLE
 
     ddtCorr                     yes;
     correctPhi                 no;
+    massFluxInterpolation      no;
+
+    moveMeshOuterCorrectors no;
+    turbOnFinalIterOnly no;
 
     oversetAdjustPhi    no;