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

overset modifications: allow overset pacthes overlap, allow fringe faces walk, mass conservation update

See merge request Development/openfoam!555

applications/solvers/basic/laplacianFoam/overLaplacianDyMFoam/Make/options

@@ -3,6 +3,7 @@ EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
     -I$(LIB_SRC)/overset/lnInclude \
+    -I$(LIB_SRC)/overset/include/lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \

applications/solvers/basic/laplacianFoam/overLaplacianDyMFoam/createFields.H

@@ -31,3 +31,25 @@
     Info<< "Reading diffusivity DT\n" << endl;
 
     dimensionedScalar DT("DT", dimViscosity, transportProperties);
+
+    bool oversetPatchErrOutput =
+        simple.dict().getOrDefault("oversetPatchErrOutput", false);
+
+    // Dummy phi for oversetPatchErrOutput
+    tmp<surfaceScalarField> tdummyPhi;
+    if (oversetPatchErrOutput)
+    {
+        tdummyPhi = tmp<surfaceScalarField>::New
+        (
+            IOobject
+            (
+                "dummyPhi",
+                runTime.timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimensionedScalar(dimless, Zero)
+        );
+    }

applications/solvers/basic/laplacianFoam/overLaplacianDyMFoam/overLaplacianDyMFoam.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2016-2017 OpenCFD Ltd.
+    Copyright (C) 2016-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -58,6 +58,7 @@ Description
 #include "fvOptions.H"
 #include "simpleControl.H"
 #include "dynamicFvMesh.H"
+#include "oversetPatchPhiErr.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -99,6 +100,11 @@ int main(int argc, char *argv[])
             fvOptions.constrain(TEqn);
             TEqn.solve();
             fvOptions.correct(T);
+
+            if (oversetPatchErrOutput)
+            {
+                oversetPatchPhiErr(TEqn, tdummyPhi.ref());
+            }
         }
 
         #include "write.H"

applications/solvers/basic/potentialFoam/overPotentialFoam/overPotentialFoam.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017 OpenCFD Ltd
+    Copyright (C) 2017-2022 OpenCFD Ltd
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -149,7 +149,6 @@ int main(int argc, char *argv[])
 
     mesh.update();
 
-    surfaceScalarField faceMask(localMin<scalar>(mesh).interpolate(cellMask));
 
     // Since solver contains no time loop it would never execute
     // function objects so do it ourselves

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

@@ -1,4 +0,0 @@
-bool ddtCorr
-(
-    pimple.dict().getOrDefault("ddtCorr", true)
-);

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

@@ -69,6 +69,8 @@ mesh.setFluxRequired(p.name());
 
 // Mask field for zeroing out contributions on hole cells
 #include "createCellMask.H"
+// Create bool field with interpolated cells
+#include "createInterpolatedCells.H"
 
 Info<< "Creating turbulence model\n" << endl;
 autoPtr<compressible::turbulenceModel> turbulence

applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/overRhoPimpleDyMFoam.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2016-2017 OpenCFD Ltd.
+    Copyright (C) 2016-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -43,7 +43,6 @@ Description
 #include "dynamicFvMesh.H"
 #include "fluidThermo.H"
 #include "turbulentFluidThermoModel.H"
-#include "bound.H"
 #include "pimpleControl.H"
 #include "pressureControl.H"
 #include "CorrectPhi.H"
@@ -89,10 +88,8 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readControls.H"
         #include "readDyMControls.H"
 
-
         // Store divrhoU from the previous mesh so that it can be mapped
         // and used in correctPhi to ensure the corrected phi has the
         // same divergence
@@ -128,7 +125,6 @@ int main(int argc, char *argv[])
         {
             if (pimple.firstIter() || moveMeshOuterCorrectors)
             {
-
                 // Do any mesh changes
                 mesh.update();
 
@@ -137,52 +133,22 @@ int main(int argc, char *argv[])
                     MRF.update();
 
                     #include "setCellMask.H"
-
+                    #include "setInterpolatedCells.H"
-                    const surfaceScalarField faceMaskOld
+                    #include "correctRhoPhiFaceMask.H"
-                    (
-                        localMin<scalar>(mesh).interpolate(cellMask.oldTime())
-                    );
-
-                    // Zero Uf on old faceMask (H-I)
-                    rhoUf() *= faceMaskOld;
-
-                    surfaceVectorField rhoUfint(fvc::interpolate(rho*U));
-
-                    // Update Uf and phi on new C-I faces
-                    rhoUf() += (1-faceMaskOld)*rhoUfint;
-
-                    // Update Uf boundary
-                    forAll(rhoUf().boundaryField(), patchI)
-                    {
-                        rhoUf().boundaryFieldRef()[patchI] =
-                            rhoUfint.boundaryField()[patchI];
-                    }
-
-                    // Calculate absolute flux from the mapped surface velocity
-                    phi = mesh.Sf() & rhoUf();
 
                     if (correctPhi)
                     {
+                        // Corrects flux on separated regions
                         #include "correctPhi.H"
                     }
 
-                    // Zero phi on current H-I
-                    const surfaceScalarField faceMask
-                    (
-                        localMin<scalar>(mesh).interpolate(cellMask)
-                    );
-
-                    phi *= faceMask;
-                    U   *= cellMask;
-
                      // Make the fluxes relative to the mesh-motion
                     fvc::makeRelative(phi, rho, U);
 
-                }
+                    if (checkMeshCourantNo)
-
+                    {
-                if (checkMeshCourantNo)
+                        #include "meshCourantNo.H"
-                {
+                    }
-                    #include "meshCourantNo.H"
                 }
             }
 

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

@@ -25,17 +25,6 @@ surfaceScalarField phiHbyA
     fvc::interpolate(rho)*fvc::flux(HbyA)
 );
 
-if (ddtCorr)
-{
-    surfaceScalarField faceMaskOld
-    (
-        localMin<scalar>(mesh).interpolate(cellMask.oldTime())
-    );
-
-    phiHbyA +=
-        faceMaskOld*MRF.zeroFilter(rhorAUf*fvc::ddtCorr(rho, U, phi, rhoUf));
-}
-
 fvc::makeRelative(phiHbyA, rho, U);
 MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
@@ -134,8 +123,4 @@ if (thermo.dpdt())
     }
 }
 
-surfaceScalarField faceMask
-(
-    localMin<scalar>(mesh).interpolate(cellMask)
-);
 phi *= faceMask;

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

@@ -1,9 +0,0 @@
-#include "readTimeControls.H"
-
-correctPhi = pimple.dict().getOrDefault("correctPhi", false);
-
-checkMeshCourantNo =
-    pimple.dict().getOrDefault("checkMeshCourantNo", false);
-
-
-ddtCorr = pimple.dict().getOrDefault("ddtCorr", true);

applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/createFields.H

@@ -124,3 +124,6 @@ dimensionedScalar initialMass("initialMass", fvc::domainIntegrate(rho));
 
 // Mask field for zeroing out contributions on hole cells
 #include "createCellMask.H"
+
+// Create bool field with interpolated cells
+#include "createInterpolatedCells.H"

applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/overBuoyantPimpleDyMFoam.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -50,6 +50,7 @@ Description
 #include "CorrectPhi.H"
 #include "cellCellStencilObject.H"
 #include "localMin.H"
+#include "oversetAdjustPhi.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -86,9 +87,6 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readTimeControls.H"
-
-        #include "readControls.H"
         #include "readDyMControls.H"
 
         #include "compressibleCourantNo.H"
@@ -128,45 +126,14 @@ int main(int argc, char *argv[])
                     MRF.update();
 
                     #include "setCellMask.H"
-
+                    #include "setInterpolatedCells.H"
-                    const surfaceScalarField faceMaskOld
+                    #include "correctRhoPhiFaceMask.H"
-                    (
-                        localMin<scalar>(mesh).interpolate(cellMask.oldTime())
-                    );
-
-                    // Zero Uf on old faceMask (H-I)
-                    rhoUf() *= faceMaskOld;
-
-                    //fvc::correctRhoUf(rhoUfint, rho, U, phi);
-                    surfaceVectorField rhoUfint(fvc::interpolate(rho*U));
-
-                    // Update Uf and phi on new C-I faces
-                    rhoUf() += (1-faceMaskOld)*rhoUfint;
-
-                    // Update Uf boundary
-                    forAll(rhoUf().boundaryField(), patchI)
-                    {
-                        rhoUf().boundaryFieldRef()[patchI] =
-                            rhoUfint.boundaryField()[patchI];
-                    }
-
-                    // Calculate absolute flux from the mapped surface velocity
-                    phi = mesh.Sf() & rhoUf();
 
                     if (correctPhi)
                     {
                         #include "correctPhi.H"
                     }
 
-                    // Zero phi on current H-I
-                    const surfaceScalarField faceMask
-                    (
-                        localMin<scalar>(mesh).interpolate(cellMask)
-                    );
-
-                    phi *= faceMask;
-                    U   *= cellMask;
-
                      // Make the fluxes relative to the mesh-motion
                     fvc::makeRelative(phi, rho, U);
                 }

applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/pEqn.H

@@ -21,17 +21,14 @@ surfaceScalarField phiHbyA
     fvc::flux(rho*HbyA) + phig
 );
 
-if (ddtCorr)
+if (adjustFringe)
 {
-    surfaceScalarField faceMaskOld
+    fvc::makeRelative(phiHbyA,rho, U);
-    (
+    oversetAdjustPhi(phiHbyA, U);
-        localMin<scalar>(mesh).interpolate(cellMask.oldTime())
+    fvc::makeAbsolute(phiHbyA,rho, U);
-    );
-
-    phiHbyA +=
-        faceMaskOld*MRF.zeroFilter(rhorAUf*fvc::ddtCorr(rho, U, phi));
 }
 
+
 MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
 // Update the pressure BCs to ensure flux consistency
@@ -122,8 +119,4 @@ if (thermo.dpdt())
     }
 }
 
-surfaceScalarField faceMask
-(
-    localMin<scalar>(mesh).interpolate(cellMask)
-);
 phi *= faceMask;

applications/solvers/heatTransfer/buoyantPimpleFoam/overBuoyantPimpleDyMFoam/readControls.H

@@ -1,9 +0,0 @@
-#include "readTimeControls.H"
-
-correctPhi = pimple.dict().getOrDefault("correctPhi", false);
-
-checkMeshCourantNo =
-    pimple.dict().getOrDefault("checkMeshCourantNo", false);
-
-
-ddtCorr = pimple.dict().getOrDefault("ddtCorr", true);

applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/CourantNo.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011 OpenFOAM Foundation
-    Copyright (C) 2016 OpenCFD Ltd.
+    Copyright (C) 2016-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -37,7 +37,10 @@ scalar meanCoNum = 0.0;
 
 if (mesh.nInternalFaces())
 {
-    surfaceScalarField phiMask(localMin<scalar>(mesh).interpolate(cellMask));
+    surfaceScalarField phiMask
+    (
+        localMin<scalar>(mesh).interpolate(cellMask + interpolatedCells)
+    );
 
     scalarField sumPhi(fvc::surfaceSum(mag(phiMask*phi))().internalField());
 

applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/UEqn.H

@@ -1,5 +1,4 @@
 // Solve the Momentum equation
-
 MRF.correctBoundaryVelocity(U);
 
 tmp<fvVectorMatrix> tUEqn

applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/createControls.H

@@ -1,26 +0,0 @@
-#include "createTimeControls.H"
-
-bool correctPhi
-(
-    pimple.dict().getOrDefault("correctPhi", false)
-);
-
-bool checkMeshCourantNo
-(
-    pimple.dict().getOrDefault("checkMeshCourantNo", false)
-);
-
-bool massFluxInterpolation
-(
-    pimple.dict().getOrDefault("massFluxInterpolation", false)
-);
-
-bool adjustFringe
-(
-    pimple.dict().getOrDefault("oversetAdjustPhi", false)
-);
-
-bool ddtCorr
-(
-    pimple.dict().getOrDefault("ddtCorr", true)
-);

applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/interpolatedFaces.H

@@ -1,273 +0,0 @@
-// Interpolation used
-interpolationCellPoint<vector> UInterpolator(HbyA);
-
-// Determine faces on outside of interpolated cells
-bitSet isOwnerInterpolatedFace(mesh.nInternalFaces());
-bitSet isNeiInterpolatedFace(mesh.nInternalFaces());
-
-// Determine donor cells
-labelListList donorCell(mesh.nInternalFaces());
-
-scalarListList weightCellCells(mesh.nInternalFaces());
-
-// Interpolated HbyA faces
-vectorField UIntFaces(mesh.nInternalFaces(), Zero);
-
-// Determine receptor neighbour cells
-labelList receptorNeigCell(mesh.nInternalFaces(), -1);
-
-{
-    const cellCellStencilObject& overlap = Stencil::New(mesh);
-    const labelList& cellTypes = overlap.cellTypes();
-    const labelIOList& zoneID = overlap.zoneID();
-
-    label nZones = gMax(zoneID)+1;
-    PtrList<fvMeshSubset> meshParts(nZones);
-    labelList nCellsPerZone(nZones, Zero);
-
-    // A mesh subset for each zone
-    forAll(meshParts, zonei)
-    {
-        meshParts.set
-        (
-            zonei,
-            // Select cells where the zoneID == zonei
-            new fvMeshSubset(mesh, zonei, zoneID)
-        );
-    }
-
-    for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
-    {
-        label ownType = cellTypes[mesh.faceOwner()[faceI]];
-        label neiType = cellTypes[mesh.faceNeighbour()[faceI]];
-        if
-        (
-            ownType == cellCellStencil::INTERPOLATED
-            && neiType == cellCellStencil::CALCULATED
-        )
-        {
-            isOwnerInterpolatedFace.set(faceI);
-
-            const vector& fc = mesh.faceCentres()[faceI];
-
-            for (label zoneI = 0; zoneI < nZones; zoneI++)
-            {
-                if (zoneI != zoneID[mesh.faceOwner()[faceI]])
-                {
-                    const fvMesh& partMesh = meshParts[zoneI].subMesh();
-                    const labelList& cellMap = meshParts[zoneI].cellMap();
-                    label cellI = partMesh.findCell(fc);
-
-                    if (cellI != -1)
-                    {
-                         // Determine weights
-                        labelList stencil(partMesh.cellCells()[cellI]);
-
-                        stencil.append(cellI);
-
-                        label st = stencil.size();
-
-                        donorCell[faceI].setSize(st);
-
-                        weightCellCells[faceI].setSize(st);
-
-                        scalarField weights(st);
-
-                        forAll(stencil, i)
-                        {
-                            scalar d = mag
-                            (
-                                partMesh.cellCentres()[stencil[i]]
-                              - fc
-                            );
-                            weights[i] = 1.0/d;
-                            donorCell[faceI][i] = cellMap[stencil[i]];
-                        }
-                        weights /= sum(weights);
-
-                        weightCellCells[faceI] = weights;
-
-                        forAll(stencil, i)
-                        {
-                            UIntFaces[faceI] +=
-                                weightCellCells[faceI][i]
-                               *UInterpolator.interpolate
-                                (
-                                    fc,
-                                    donorCell[faceI][i]
-                                );
-                        }
-
-                        break;
-                    }
-                }
-            }
-
-            receptorNeigCell[faceI] = mesh.faceNeighbour()[faceI];
-        }
-        else if
-        (
-            ownType == cellCellStencil::CALCULATED
-            && neiType == cellCellStencil::INTERPOLATED
-        )
-        {
-            isNeiInterpolatedFace.set(faceI);
-
-            const vector& fc = mesh.faceCentres()[faceI];
-            for (label zoneI = 0; zoneI < nZones; zoneI++)
-            {
-                if (zoneI != zoneID[mesh.faceNeighbour()[faceI]])
-                {
-                    const fvMesh& partMesh = meshParts[zoneI].subMesh();
-                    const labelList& cellMap = meshParts[zoneI].cellMap();
-                    label cellI = partMesh.findCell(fc);
-
-                    if (cellI != -1)
-                    {
-                        // Determine weights
-                        labelList stencil(partMesh.cellCells()[cellI]);
-
-                        stencil.append(cellI);
-
-                        label st = stencil.size();
-
-                        donorCell[faceI].setSize(st);
-
-                        weightCellCells[faceI].setSize(st);
-
-                        scalarField weights(st);
-
-                        forAll(stencil, i)
-                        {
-                            scalar d = mag
-                            (
-                                partMesh.cellCentres()[stencil[i]]
-                              - fc
-                            );
-                            weights[i] = 1.0/d;
-                            donorCell[faceI][i] = cellMap[stencil[i]];
-                        }
-                        weights /= sum(weights);
-
-                        weightCellCells[faceI] = weights;
-
-                        forAll(stencil, i)
-                        {
-                            UIntFaces[faceI] +=
-                                weightCellCells[faceI][i]
-                               *UInterpolator.interpolate
-                                (
-                                    fc,
-                                    donorCell[faceI][i]
-                                );
-                        }
-
-                        break;
-                    }
-                }
-            }
-
-            receptorNeigCell[faceI] = mesh.faceOwner()[faceI];
-        }
-    }
-}
-
-// contravariant U
-vectorField U1Contrav(mesh.nInternalFaces(), Zero);
-
-surfaceVectorField faceNormals(mesh.Sf()/mesh.magSf());
-
-forAll(isNeiInterpolatedFace, faceI)
-{
-    label cellId = -1;
-    if (isNeiInterpolatedFace.test(faceI))
-    {
-        cellId = mesh.faceNeighbour()[faceI];
-    }
-    else if (isOwnerInterpolatedFace.test(faceI))
-    {
-        cellId = mesh.faceOwner()[faceI];
-    }
-
-    if (cellId != -1)
-    {
-        const vector& n = faceNormals[faceI];
-        vector n1(Zero);
-
-        // 2-D cases
-        if (mesh.nSolutionD() == 2)
-        {
-            for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
-            {
-                if (mesh.geometricD()[cmpt] == -1)
-                {
-                    switch (cmpt)
-                    {
-                        case vector::X:
-                        {
-                            n1 = vector(0, n.z(), -n.y());
-                            break;
-                        }
-
-                        case vector::Y:
-                        {
-                            n1 = vector(n.z(), 0, -n.x());
-                            break;
-                        }
-
-                        case vector::Z:
-                        {
-                            n1 = vector(n.y(), -n.x(), 0);
-                            break;
-                        }
-                    }
-                }
-            }
-        }
-        else if (mesh.nSolutionD() == 3)
-        {
-            //Determine which is the primary direction
-            if (mag(n.x()) > mag(n.y()) && mag(n.x()) > mag(n.z()))
-            {
-                n1 = vector(n.y(), -n.x(), 0);
-            }
-            else if (mag(n.y()) > mag(n.z()))
-            {
-                n1 = vector(0, n.z(), -n.y());
-            }
-            else
-            {
-                n1 = vector(-n.z(), 0, n.x());
-            }
-        }
-        n1.normalise();
-
-        const vector n2 = normalised(n ^ n1);
-
-        tensor rot =
-            tensor
-            (
-               n.x() ,n.y(), n.z(),
-               n1.x() ,n1.y(), n1.z(),
-               n2.x() ,n2.y(), n2.z()
-            );
-
-//         tensor rot =
-//             tensor
-//             (
-//                n  & x ,n  & y, n  & z,
-//                n1 & x ,n1 & y, n1 & z,
-//                n2 & x ,n2 & y, n2 & z
-//             );
-
-        U1Contrav[faceI].x() =
-           2*transform(rot, UIntFaces[faceI]).x()
-           - transform(rot, HbyA[receptorNeigCell[faceI]]).x();
-
-        U1Contrav[faceI].y() = transform(rot, HbyA[cellId]).y();
-
-        U1Contrav[faceI].z() = transform(rot, HbyA[cellId]).z();
-
-        HbyA[cellId] = transform(inv(rot), U1Contrav[faceI]);
-    }
-}

applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/overPimpleDyMFoam.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2018 OpenCFD Ltd.
+    Copyright (C) 2016-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -46,12 +46,9 @@ Description
 #include "fvOptions.H"
 
 #include "cellCellStencilObject.H"
-#include "zeroGradientFvPatchFields.H"
 #include "localMin.H"
-#include "interpolationCellPoint.H"
-#include "transform.H"
-#include "fvMeshSubset.H"
 #include "oversetAdjustPhi.H"
+#include "oversetPatchPhiErr.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -68,10 +65,9 @@ int main(int argc, char *argv[])
     #include "setRootCaseLists.H"
     #include "createTime.H"
     #include "createDynamicFvMesh.H"
+    #include "createDyMControls.H"
     #include "initContinuityErrs.H"
 
-    pimpleControl pimple(mesh);
-
     #include "createFields.H"
     #include "createUf.H"
     #include "createMRF.H"
@@ -88,7 +84,9 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readControls.H"
+        #include "readDyMControls.H"
+        #include "readOversetDyMControls.H"
+
         #include "CourantNo.H"
 
         #include "setDeltaT.H"
@@ -97,45 +95,20 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        bool changed = mesh.update();
+        mesh.update();
 
-        if (changed)
+        if (mesh.changing())
         {
             #include "setCellMask.H"
             #include "setInterpolatedCells.H"
+            #include "correctPhiFaceMask.H"
 
-            surfaceScalarField faceMaskOld
+            fvc::makeRelative(phi, U);
-            (
-                localMin<scalar>(mesh).interpolate(cellMask.oldTime())
-            );
 
-            // Zero Uf on old faceMask (H-I)
+            if (checkMeshCourantNo)
-            Uf *= faceMaskOld;
+            {
-            // Update Uf and phi on new C-I faces
+                #include "meshCourantNo.H"
-            Uf += (1-faceMaskOld)*fvc::interpolate(U);
+            }
-            phi = mesh.Sf() & Uf;
-
-            // Zero phi on current H-I
-            surfaceScalarField faceMask
-            (
-                localMin<scalar>(mesh).interpolate(cellMask)
-            );
-            phi *= faceMask;
-        }
-
-
-        if (mesh.changing() && correctPhi)
-        {
-            // Calculate absolute flux from the mapped surface velocity
-            #include "correctPhi.H"
-        }
-
-        // Make the flux relative to the mesh motion
-        fvc::makeRelative(phi, U);
-
-        if (mesh.changing() && checkMeshCourantNo)
-        {
-            #include "meshCourantNo.H"
         }
 
         // --- Pressure-velocity PIMPLE corrector loop

applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/pEqn.H

@@ -1,36 +1,11 @@
-// Option 1: interpolate rAU, do not block out rAU on blocked cells
 volScalarField rAU("rAU", 1.0/UEqn.A());
-mesh.interpolate(rAU);
 
-// Option 2: do not interpolate rAU but block out rAU
-//surfaceScalarField rAUf("rAUf", fvc::interpolate(blockedCells*rAU));
-
-
-// Option 3: do not interpolate rAU but zero out rAUf on faces on holes
-// But what about:
-//
-//   H
-// H I C C C C
-//   H
-//
 surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
 volVectorField H("H", UEqn.H());
 
 volVectorField HbyA("HbyA", U);
 HbyA = constrainHbyA(rAU*H, U, p);
 
-if (massFluxInterpolation)
-{
-    #include "interpolatedFaces.H"
-}
-
-if (runTime.outputTime())
-{
-    H.write();
-    rAU.write();
-    HbyA.write();
-}
-
 if (pimple.nCorrPISO() <= 1)
 {
     tUEqn.clear();
@@ -38,33 +13,16 @@ if (pimple.nCorrPISO() <= 1)
 
 phiHbyA = fvc::flux(HbyA);
 
-if (ddtCorr)
-{
-    surfaceScalarField faceMaskOld
-    (
-        localMin<scalar>(mesh).interpolate(cellMask.oldTime())
-    );
-    phiHbyA += rAUf*faceMaskOld*fvc::ddtCorr(U, Uf);
-}
-
 MRF.makeRelative(phiHbyA);
 
-// WIP
-if (p.needReference())
-{
-    fvc::makeRelative(phiHbyA, U);
-    adjustPhi(phiHbyA, U, p);
-    fvc::makeAbsolute(phiHbyA, U);
-}
-
-
 if (adjustFringe)
 {
     fvc::makeRelative(phiHbyA, U);
-    oversetAdjustPhi(phiHbyA, U);
+    oversetAdjustPhi(phiHbyA, U, zoneIdMass);
     fvc::makeAbsolute(phiHbyA, U);
 }
 
+
 while (pimple.correctNonOrthogonal())
 {
     fvScalarMatrix pEqn
@@ -79,27 +37,26 @@ while (pimple.correctNonOrthogonal())
     if (pimple.finalNonOrthogonalIter())
     {
         phi = phiHbyA - pEqn.flux();
-        // option 2:
+        pEqn.relax();
-        // rAUf*fvc::snGrad(p)*mesh.magSf();
+        U =
+            cellMask*
+            (
+                HbyA
+              - rAU*fvc::reconstruct((pEqn.flux())/rAUf)
+            );
+        U.correctBoundaryConditions();
+        fvOptions.correct(U);
+    }
+
+    if (oversetPatchErrOutput)
+    {
+        oversetPatchPhiErr(pEqn, phiHbyA);
     }
 }
 
-
+// Excludes error in interpolated/hole cells
 #include "continuityErrs.H"
 
-// Explicitly relax pressure for momentum corrector
-p.relax();
-volVectorField gradP(fvc::grad(p));
-
-// Option 2: zero out velocity on blocked out cells
-//U = HbyA - rAU*cellMask*gradP;
-// Option 3: zero out velocity on blocked out cells
-// This is needed for the scalar Eq (k,epsilon, etc)
-// which can use U as source term
-U = cellMask*(HbyA - rAU*gradP);
-U.correctBoundaryConditions();
-
-fvOptions.correct(U);
 
 {
     Uf = fvc::interpolate(U);
@@ -109,9 +66,4 @@ fvOptions.correct(U);
 
 // Make the fluxes relative to the mesh motion
 fvc::makeRelative(phi, U);
-
-surfaceScalarField faceMask
-(
-    localMin<scalar>(mesh).interpolate(cellMask)
-);
 phi *= faceMask;

applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/readControls.H

@@ -1,10 +0,0 @@
-#include "readTimeControls.H"
-
-correctPhi = pimple.dict().getOrDefault("correctPhi", false);
-
-checkMeshCourantNo = pimple.dict().getOrDefault("checkMeshCourantNo", false);
-
-massFluxInterpolation =
-    pimple.dict().getOrDefault("massFluxInterpolation", false);
-
-ddtCorr = pimple.dict().getOrDefault("ddtCorr", true);

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

@@ -24,7 +24,3 @@ bool adjustFringe
 (
     simple.dict().getOrDefault("oversetAdjustPhi", false)
 );
-bool massFluxInterpolation
-(
-    simple.dict().getOrDefault("massFluxInterpolation", false)
-);

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

@@ -1,18 +1,10 @@
 {
-    surfaceScalarField faceMask(localMin<scalar>(mesh).interpolate(cellMask));
-
     volScalarField rAU(1.0/UEqn.A());
     surfaceScalarField rAUf("rAUf", faceMask*fvc::interpolate(rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = constrainHbyA(cellMask*rAU*UEqn.H(), U, p);
 
-    //mesh.interpolate(HbyA);
-    if (massFluxInterpolation)
-    {
-        #include "interpolatedFaces.H"
-    }
-
     tUEqn.clear();
 
     surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA));

applications/solvers/multiphase/compressibleInterFoam/overCompressibleInterDyMFoam/createFields.H

@@ -129,10 +129,5 @@ compressibleInterPhaseTransportModel turbulence
 // Mask field for zeroing out contributions on hole cells
 #include "createCellMask.H"
 
-surfaceScalarField faceMask
-(
-    localMin<scalar>(mesh).interpolate(cellMask)
-);
-
 // Create bool field with interpolated cells
 #include "createInterpolatedCells.H"

applications/solvers/multiphase/compressibleInterFoam/overCompressibleInterDyMFoam/overCompressibleInterDyMFoam.C

@@ -80,9 +80,6 @@ int main(int argc, char *argv[])
     #include "setRootCaseLists.H"
     #include "createTime.H"
     #include "createDynamicFvMesh.H"
-    pimpleControl pimple(mesh);
-
-    #include "createTimeControls.H"
     #include "createDyMControls.H"
     #include "createFields.H"
 
@@ -109,7 +106,7 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readControls.H"
+        #include "readDyMControls.H"
 
         if (LTS)
         {
@@ -154,40 +151,10 @@ int main(int argc, char *argv[])
                     // Update cellMask field for blocking out hole cells
                     #include "setCellMask.H"
                     #include "setInterpolatedCells.H"
-
+                    #include "correctPhiFaceMask.H"
-                    faceMask =
-                        localMin<scalar>(mesh).interpolate(cellMask.oldTime());
-
-                    // Zero Uf on old faceMask (H-I)
-                    Uf *= faceMask;
-
-                    const surfaceVectorField Uint(fvc::interpolate(U));
-                    // Update Uf and phi on new C-I faces
-                    Uf += (1-faceMask)*Uint;
-
-                    // Update Uf boundary
-                    forAll(Uf.boundaryField(), patchI)
-                    {
-                        Uf.boundaryFieldRef()[patchI] =
-                            Uint.boundaryField()[patchI];
-                    }
-
-                    phi = mesh.Sf() & Uf;
-
-                    // Correct phi on individual regions
-                    if (correctPhi)
-                    {
-                         #include "correctPhi.H"
-                    }
 
                     mixture.correct();
 
-                    // Zero phi on current H-I
-                    faceMask = localMin<scalar>(mesh).interpolate(cellMask);
-
-                    phi *= faceMask;
-                    U *= cellMask;
-
                     // Make the flux relative to the mesh motion
                     fvc::makeRelative(phi, U);
 
@@ -202,10 +169,6 @@ int main(int argc, char *argv[])
             #include "alphaControls.H"
             #include "compressibleAlphaEqnSubCycle.H"
 
-            const surfaceScalarField faceMask
-            (
-                localMin<scalar>(mesh).interpolate(cellMask)
-            );
             rhoPhi *= faceMask;
 
             turbulence.correctPhasePhi();

applications/solvers/multiphase/compressibleInterFoam/overCompressibleInterDyMFoam/pEqn.H

@@ -10,19 +10,6 @@
         fvc::flux(HbyA)
     );
 
-    if (ddtCorr)
-    {
-        surfaceScalarField faceMaskOld
-        (
-            localMin<scalar>(mesh).interpolate(cellMask.oldTime())
-        );
-        phiHbyA +=
-            MRF.zeroFilter
-            (
-                fvc::interpolate(rho*rAU)*faceMaskOld*fvc::ddtCorr(U, Uf)
-            );
-    }
-
     MRF.makeRelative(phiHbyA);
 
     surfaceScalarField phig

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

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011 OpenFOAM Foundation
+    Copyright (C) 2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -36,7 +36,10 @@ scalar meanCoNum = 0.0;
 
 if (mesh.nInternalFaces())
 {
-    surfaceScalarField phiMask(localMin<scalar>(mesh).interpolate(cellMask));
+    surfaceScalarField phiMask
+    (
+        localMin<scalar>(mesh).interpolate(cellMask + interpolatedCells)
+    );
 
     scalarField sumPhi
     (

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

@@ -7,3 +7,5 @@ volScalarField::Internal divU
   ? fvc::div(phiCN() + mesh.phi())
   : fvc::div(phiCN())
 );
+
+divU *= interpolatedCells*cellMask;

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

@@ -1,30 +0,0 @@
-bool correctPhi
-(
-    pimple.dict().getOrDefault("correctPhi", true)
-);
-
-bool checkMeshCourantNo
-(
-    pimple.dict().getOrDefault("checkMeshCourantNo", false)
-);
-
-bool moveMeshOuterCorrectors
-(
-    pimple.dict().getOrDefault("moveMeshOuterCorrectors", false)
-);
-
-
-bool massFluxInterpolation
-(
-    pimple.dict().getOrDefault("massFluxInterpolation", false)
-);
-
-bool adjustFringe
-(
-    pimple.dict().getOrDefault("oversetAdjustPhi", false)
-);
-
-bool ddtCorr
-(
-    pimple.dict().getOrDefault("ddtCorr", true)
-);

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

@@ -40,8 +40,12 @@ volVectorField U
     nonInt.insert("HbyA");
     nonInt.insert("grad(p_rgh)");
     nonInt.insert("nHat");
+    nonInt.insert("surfaceIntegrate(nHatf)");
+    nonInt.insert("surfaceIntegrate(phi+meshPhi)");
     nonInt.insert("surfaceIntegrate(phi)");
     nonInt.insert("surfaceIntegrate(phiHbyA)");
+    nonInt.insert("surfaceSum(((S|magSf)*S)");
+    nonInt.insert("surfaceIntegrate(((rAUf*magSf)*snGradCorr(p_rgh)))");
     nonInt.insert("cellMask");
     nonInt.insert("cellDisplacement");
     nonInt.insert("interpolatedCells");

applications/solvers/multiphase/interFoam/overInterDyMFoam/overInterDyMFoam.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2014 OpenFOAM Foundation
-    Copyright (C) 2016-2017 OpenCFD Ltd.
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -49,14 +49,11 @@ Description
 #include "turbulentTransportModel.H"
 #include "pimpleControl.H"
 #include "fvOptions.H"
-#include "CorrectPhi.H"
 #include "fvcSmooth.H"
 #include "cellCellStencilObject.H"
 #include "localMin.H"
-#include "interpolationCellPoint.H"
-#include "transform.H"
-#include "fvMeshSubset.H"
 #include "oversetAdjustPhi.H"
+#include "oversetPatchPhiErr.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -76,8 +73,7 @@ int main(int argc, char *argv[])
     #include "createTime.H"
     #include "createDynamicFvMesh.H"
     #include "initContinuityErrs.H"
-    pimpleControl pimple(mesh);
+
-    #include "createTimeControls.H"
     #include "createDyMControls.H"
     #include "createFields.H"
     #include "createAlphaFluxes.H"
@@ -97,11 +93,8 @@ int main(int argc, char *argv[])
         dimensionedScalar("rAUf", dimTime/rho.dimensions(), 1.0)
     );
 
-    if (correctPhi)
-    {
-        #include "correctPhi.H"
-    }
     #include "createUf.H"
+    #include "createControls.H"
 
     #include "setCellMask.H"
     #include "setInterpolatedCells.H"
@@ -119,7 +112,8 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readControls.H"
+        #include "readDyMControls.H"
+        #include "readOversetDyMControls.H"
 
         if (LTS)
         {
@@ -158,50 +152,17 @@ int main(int argc, char *argv[])
                         talphaPhi1Corr0.clear();
                     }
 
-                    gh = (g & mesh.C()) - ghRef;
-                    ghf = (g & mesh.Cf()) - ghRef;
-
                     // Update cellMask field for blocking out hole cells
                     #include "setCellMask.H"
                     #include "setInterpolatedCells.H"
+                    #include "correctPhiFaceMask.H"
 
-                    const surfaceScalarField faceMaskOld
+                    gh = (g & mesh.C()) - ghRef;
-                    (
+                    ghf = (g & mesh.Cf()) - ghRef;
-                        localMin<scalar>(mesh).interpolate(cellMask.oldTime())
-                    );
 
-                    // Zero Uf on old faceMask (H-I)
-                    Uf *= faceMaskOld;
-
-                    const surfaceVectorField Uint(fvc::interpolate(U));
-                    // Update Uf and phi on new C-I faces
-                    Uf += (1-faceMaskOld)*Uint;
-
-                    // Update Uf boundary
-                    forAll(Uf.boundaryField(), patchI)
-                    {
-                        Uf.boundaryFieldRef()[patchI] =
-                            Uint.boundaryField()[patchI];
-                    }
-
-                    phi = mesh.Sf() & Uf;
-
-                    // Correct phi on individual regions
-                    if (correctPhi)
-                    {
-                         #include "correctPhi.H"
-                    }
 
                     mixture.correct();
 
-                    // Zero phi on current H-I
-                    const surfaceScalarField faceMask
-                    (
-                        localMin<scalar>(mesh).interpolate(cellMask)
-                    );
-                    phi *= faceMask;
-                    U   *= cellMask;
-
                     // Make the flux relative to the mesh motion
                     fvc::makeRelative(phi, U);
 
@@ -213,14 +174,14 @@ int main(int argc, char *argv[])
                 }
             }
 
+            if (adjustFringe)
+            {
+                oversetAdjustPhi(phi, U, zoneIdMass);
+            }
 
             #include "alphaControls.H"
             #include "alphaEqnSubCycle.H"
 
-            const surfaceScalarField faceMask
-            (
-                localMin<scalar>(mesh).interpolate(cellMask)
-            );
             rhoPhi *= faceMask;
 
             mixture.correct();

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

@@ -1,64 +1,31 @@
 {
     rAU = 1.0/UEqn.A();
-    //mesh.interpolate(rAU);
-
-    surfaceScalarField faceMask(localMin<scalar>(mesh).interpolate(cellMask));
 
     surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
 
     volVectorField H("H", UEqn.H());
 
     volVectorField HbyA("HbyA", U);
-    //HbyA = rAU*UEqn.H();
-    HbyA = constrainHbyA(rAU*H, U, p_rgh);
 
-    if (massFluxInterpolation)
+    HbyA = constrainHbyA(rAU*H, U, p_rgh);
-    {
-        #include "interpolatedFaces.H"
-    }
 
     surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA));
 
-    if (ddtCorr)
-    {
-        surfaceScalarField faceMaskOld
-        (
-            localMin<scalar>(mesh).interpolate(cellMask.oldTime())
-        );
-        phiHbyA +=
-            fvc::interpolate(rho*rAU)*faceMaskOld*fvc::ddtCorr(U, Uf);
-    }
     MRF.makeRelative(phiHbyA);
 
-    if (p_rgh.needReference())
-    {
-        fvc::makeRelative(phiHbyA, U);
-        adjustPhi(phiHbyA, U, p_rgh);
-        fvc::makeAbsolute(phiHbyA, U);
-    }
-
-    if (adjustFringe)
-    {
-        fvc::makeRelative(phiHbyA, U);
-        oversetAdjustPhi(phiHbyA, U);
-        fvc::makeAbsolute(phiHbyA, U);
-    }
-
     surfaceScalarField phig
     (
         (
             mixture.surfaceTensionForce()
-          - ghf*fvc::snGrad(rho)
+          - ghf*fvc::snGrad(cellMask*rho)
-        )*faceMask*rAUf*mesh.magSf()
+        )*rAUf*faceMask*mesh.magSf()
     );
 
     phiHbyA += phig;
 
     if (adjustFringe)
     {
-        fvc::makeRelative(phiHbyA, U);
+        oversetAdjustPhi(phiHbyA, U, zoneIdMass);
-        oversetAdjustPhi(phiHbyA, U);
-        fvc::makeAbsolute(phiHbyA, U);
     }
 
     // Update the pressure BCs to ensure flux consistency
@@ -90,6 +57,10 @@
             U.correctBoundaryConditions();
             fvOptions.correct(U);
         }
+        if (oversetPatchErrOutput)
+        {
+            oversetPatchPhiErr(p_rghEqn, phiHbyA);
+        }
     }
 
     #include "continuityErrs.H"

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

@@ -1,16 +0,0 @@
-#include "readTimeControls.H"
-
-correctPhi = pimple.dict().getOrDefault("correctPhi", false);
-
-checkMeshCourantNo =
-    pimple.dict().getOrDefault("checkMeshCourantNo", false);
-
-moveMeshOuterCorrectors =
-    pimple.dict().getOrDefault("moveMeshOuterCorrectors", false);
-
-massFluxInterpolation =
-    pimple.dict().getOrDefault("massFluxInterpolation", false);
-
-ddtCorr = pimple.dict().getOrDefault("ddtCorr", true);
-
-adjustFringe =  pimple.dict().getOrDefault("oversetAdjustPhi", false);

applications/solvers/multiphase/interPhaseChangeFoam/overInterPhaseChangeDyMFoam/correctPhi.H

@@ -1,11 +0,0 @@
-CorrectPhi
-(
-    U,
-    phi,
-    p_rgh,
-    surfaceScalarField("rAUf", fvc::interpolate(rAU)),
-    divU,
-    pimple
-);
-
-#include "continuityErrs.H"

applications/solvers/multiphase/interPhaseChangeFoam/overInterPhaseChangeDyMFoam/createFields.H

@@ -149,10 +149,5 @@ surfaceScalarField alphaPhi10
 // Mask field for zeroing out contributions on hole cells
 #include "createCellMask.H"
 
-surfaceScalarField faceMask
-(
-    localMin<scalar>(mesh).interpolate(cellMask)
-);
-
 // Create bool field with interpolated cells
 #include "createInterpolatedCells.H"

applications/solvers/multiphase/interPhaseChangeFoam/overInterPhaseChangeDyMFoam/overInterPhaseChangeDyMFoam.C

@@ -82,9 +82,7 @@ int main(int argc, char *argv[])
     #include "setRootCaseLists.H"
     #include "createTime.H"
     #include "createDynamicFvMesh.H"
-    pimpleControl pimple(mesh);
 
-    #include "createTimeControls.H"
     #include "createDyMControls.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
@@ -118,7 +116,7 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readControls.H"
+        #include "readDyMControls.H"
 
         // Store divU from the previous mesh so that it can be mapped
         // and used in correctPhi to ensure the corrected phi has the
@@ -153,40 +151,10 @@ int main(int argc, char *argv[])
                      // Update cellMask field for blocking out hole cells
                     #include "setCellMask.H"
                     #include "setInterpolatedCells.H"
-
+                    #include "correctPhiFaceMask.H"
-                    faceMask =
-                        localMin<scalar>(mesh).interpolate(cellMask.oldTime());
-
-
-                    // Zero Uf on old faceMask (H-I)
-                    Uf *= faceMask;
-
-                    const surfaceVectorField Uint(fvc::interpolate(U));
-                    // Update Uf and phi on new C-I faces
-                    Uf += (1-faceMask)*Uint;
-
-                    // Update Uf boundary
-                    forAll(Uf.boundaryField(), patchI)
-                    {
-                        Uf.boundaryFieldRef()[patchI] =
-                            Uint.boundaryField()[patchI];
-                    }
-
-                    phi = mesh.Sf() & Uf;
-
-                    if (correctPhi)
-                    {
-                        #include "correctPhi.H"
-                    }
 
                     mixture->correct();
 
-                    // Zero phi on current H-I
-                    faceMask = localMin<scalar>(mesh).interpolate(cellMask);
-
-                    phi *= faceMask;
-                    U   *= cellMask;
-
                     // Make the flux relative to the mesh motion
                     fvc::makeRelative(phi, U);
                 }
@@ -214,10 +182,6 @@ int main(int argc, char *argv[])
             mixture->correct();
 
             #include "alphaEqnSubCycle.H"
-            const surfaceScalarField faceMask
-            (
-                localMin<scalar>(mesh).interpolate(cellMask)
-            );
             rhoPhi *= faceMask;
 
             interface.correct();

applications/solvers/multiphase/interPhaseChangeFoam/overInterPhaseChangeDyMFoam/pEqn.H

@@ -8,16 +8,6 @@
         fvc::flux(HbyA)
     );
 
-
-    if (ddtCorr)
-    {
-        surfaceScalarField faceMaskOld
-        (
-            localMin<scalar>(mesh).interpolate(cellMask.oldTime())
-        );
-        phiHbyA += faceMaskOld*fvc::ddtCorr(U, Uf);
-    }
-
     if (p_rgh.needReference())
     {
         fvc::makeRelative(phiHbyA, U);

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/pairGAMGAgglomeration/pairGAMGAgglomerate.C

@@ -36,6 +36,18 @@ void Foam::pairGAMGAgglomeration::agglomerate
     const scalarField& faceWeights
 )
 {
+    const label nFaces = mesh.lduAddr().lowerAddr().size();
+    if (faceWeights.size() != nFaces)
+    {
+        FatalErrorInFunction
+            << "Supplied number of face weights " << faceWeights.size()
+            << " does not correspond to the number of faces " << nFaces << endl
+            << "This may be because of using a geometry-based"
+            << " agglomeration method instead of a matrix-based one"
+            << exit(FatalError);
+    }
+
+
     // Start geometric agglomeration from the given faceWeights
     scalarField* faceWeightsPtr = const_cast<scalarField*>(&faceWeights);
 

src/fvOptions/Make/options

@@ -13,7 +13,8 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/overset/lnInclude
 
 LIB_LIBS = \
     -lfiniteVolume \
@@ -26,4 +27,5 @@ LIB_LIBS = \
     -lturbulenceModels \
     -lincompressibleTurbulenceModels \
     -lcompressibleTurbulenceModels \
-    -lreactionThermophysicalModels
+    -lreactionThermophysicalModels \
+    -loverset

src/fvOptions/cellSetOption/cellSetOption.C

@@ -30,6 +30,7 @@ License
 #include "cellSet.H"
 #include "cellBitSet.H"
 #include "volFields.H"
+#include "cellCellStencilObject.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -53,6 +54,7 @@ Foam::fv::cellSetOption::selectionModeTypeNames_
     { selectionModeType::smPoints, "points" },
     { selectionModeType::smCellSet, "cellSet" },
     { selectionModeType::smCellZone, "cellZone" },
+    { selectionModeType::smCellType, "cellType" }
 });
 
 
@@ -97,6 +99,10 @@ void Foam::fv::cellSetOption::setSelection(const dictionary& dict)
             }
             break;
         }
+        case smCellType:
+        {
+            break;
+        }
         default:
         {
             FatalErrorInFunction
@@ -234,6 +240,21 @@ void Foam::fv::cellSetOption::setCellSelection()
             }
             break;
         }
+        case smCellType:
+        {
+            labelHashSet selectedCells;
+            const cellCellStencilObject& overlap = Stencil::New(mesh_);
+            const labelList& cellTypes = overlap.cellTypes();
+            forAll(cellTypes, celli)
+            {
+                if (cellTypes[celli] == cellCellStencil::POROUS)
+                {
+                    selectedCells.insert(celli);
+                }
+                cells_ = selectedCells.sortedToc();
+            }
+            break;
+        }
         default:
         {
             FatalErrorInFunction
@@ -300,6 +321,7 @@ bool Foam::fv::cellSetOption::isActive()
             (
                 selectionMode_ == smGeometric
              || selectionMode_ == smPoints
+             || selectionMode_ == smCellType
             )
             {
                 // Geometric selection mode(s)

src/fvOptions/cellSetOption/cellSetOption.H

@@ -162,7 +162,8 @@ public:
             smCellSet,      //!< "cellSet"
             smCellZone,     //!< "cellZone"
             smPoints,       //!< "points"
-            smGeometric     //!< "geometric"
+            smGeometric,    //!< "geometric"
+            smCellType      //!< "overset type cells"
         };
 
         //- List of selection mode type names

src/fvOptions/constraints/derived/velocityDampingConstraint/velocityDampingConstraint.C

@@ -30,6 +30,7 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "fvMatrix.H"
 #include "volFields.H"
+#include "cellCellStencilObject.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -67,14 +68,14 @@ void Foam::fv::velocityDampingConstraint::addDamping(fvMatrix<vector>& eqn)
     label nDamped(0);
     const label nTotCells(returnReduce(cells_.size(), sumOp<label>()));
 
-    for (label celli : cells_)
+    for (const label celli : cells_)
     {
         const scalar magU = mag(U[celli]);
         if (magU > UMax_)
         {
             const scalar scale = sqr(Foam::cbrt(vol[celli]));
 
-            diag[celli] += scale*(magU-UMax_);
+            diag[celli] += C_*scale*(magU-UMax_);
 
             ++nDamped;
         }
@@ -105,7 +106,8 @@ Foam::fv::velocityDampingConstraint::velocityDampingConstraint
     const fvMesh& mesh
 )
 :
-    fv::cellSetOption(name, modelType, dict, mesh)
+    fv::cellSetOption(name, modelType, dict, mesh),
+    C_(1)
 {
     read(dict);
 }
@@ -135,6 +137,8 @@ bool Foam::fv::velocityDampingConstraint::read(const dictionary& dict)
     {
         coeffs_.readEntry("UMax", UMax_);
 
+        coeffs_.readIfPresent("C", C_);
+
         if (!coeffs_.readIfPresent("UNames", fieldNames_))
         {
             fieldNames_.resize(1);

src/fvOptions/constraints/derived/velocityDampingConstraint/velocityDampingConstraint.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2015 OpenFOAM Foundation
-    Copyright (C) 2015-2020 OpenCFD Ltd.
+    Copyright (C) 2015-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -62,6 +62,7 @@ Usage
 
         // Optional entries (runtime modifiable)
         UNames          (<Uname1> <Uname2> ... <UnameN>);
+        C               <scalar>;
 
         // Conditional optional entries (runtime modifiable)
 
@@ -80,6 +81,7 @@ Usage
       UMax      | Maximum velocity magnitude            | scalar | yes   | -
       UNames    | Names of operand velocity fields      | wordList | no  | -
       U         | Name of operand velocity field        | word   | cndtnl | U
+      C         | Model factor                          | scalar | no    | 1
     \endtable
 
     The inherited entries are elaborated in:
@@ -127,6 +129,9 @@ protected:
         //- Maximum velocity magnitude
         scalar UMax_;
 
+        //- Model factor
+        scalar C_;
+
 
     // Protected Member Functions
 

src/overset/Make/files

@@ -8,7 +8,10 @@ cellCellStencil/trackingInverseDistance/voxelMeshSearch.C
 cellCellStencil/trackingInverseDistance/trackingInverseDistanceCellCellStencil.C
 cellCellStencil/leastSquares/leastSquaresCellCellStencil.C
 
-dynamicOversetFvMesh/dynamicOversetFvMesh.C
+oversetFvMesh = oversetFvMesh
+$(oversetFvMesh)/oversetFvMeshBase.C
+$(oversetFvMesh)/dynamicOversetFvMesh/dynamicOversetFvMesh.C
+$(oversetFvMesh)/staticOversetFvMesh/staticOversetFvMesh.C
 
 fvMeshPrimitiveLduAddressing/fvMeshPrimitiveLduAddressing.C
 
@@ -23,8 +26,15 @@ oversetAdjustPhi/oversetAdjustPhi.C
 
 regionsToCell/regionsToCell.C
 
-
 lduPrimitiveProcessorInterface/GAMG/calculatedProcessorGAMGInterface.C
 lduPrimitiveProcessorInterface/GAMG/calculatedProcessorGAMGInterfaceField.C
 
+oversetCoupledPolyPatch/oversetLduInterfaceField/oversetLduInterfaceField.C
+oversetCoupledPolyPatch/oversetLduInterface/oversetLduInterface.C
+
+oversetCoupledPolyPatch/oversetGAMGInterface/oversetGAMGInterface.C
+oversetCoupledPolyPatch/oversetGAMGInterfaceField/oversetGAMGInterfaceField.C
+
+oversetPatchPhiErr/oversetPatchPhiErr.C
+
 LIB = $(FOAM_LIBBIN)/liboverset

src/overset/Make/options

@@ -1,7 +1,3 @@
-/* Add to EXE_INC as required:
-   -DFULLDEBUG -O0 -g
- */
-
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/fileFormats/lnInclude \

src/overset/cellCellStencil/cellCellStencil/cellCellStencil.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2021 OpenCFD Ltd.
+    Copyright (C) 2017-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -30,6 +30,7 @@ License
 #include "volFields.H"
 #include "syncTools.H"
 #include "globalIndex.H"
+#include "oversetFvPatchFields.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -48,6 +49,8 @@ Foam::cellCellStencil::cellTypeNames_
     { cellType::CALCULATED, "calculated" },
     { cellType::INTERPOLATED, "interpolated" },
     { cellType::HOLE, "hole" },
+    { cellType::SPECIAL, "special" },
+    { cellType::POROUS, "porous" },
 });
 
 
@@ -96,6 +99,17 @@ Foam::cellCellStencil::~cellCellStencil()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+Foam::word Foam::cellCellStencil::baseName(const word& name)
+{
+    if (name.ends_with("_0"))
+    {
+        return baseName(name.substr(0, name.size() - 2));
+    }
+
+    return name;
+}
+
+
 const Foam::labelIOList& Foam::cellCellStencil::zoneID(const fvMesh& mesh)
 {
     if (!mesh.foundObject<labelIOList>("zoneID"))
@@ -128,9 +142,9 @@ const Foam::labelIOList& Foam::cellCellStencil::zoneID(const fvMesh& mesh)
             ),
             mesh
         );
-        forAll(volZoneID, cellI)
+        forAll(volZoneID, celli)
         {
-            zoneID[cellI] = label(volZoneID[cellI]);
+            zoneID[celli] = label(volZoneID[celli]);
         }
 
         zoneIDPtr->store();
@@ -296,5 +310,512 @@ void Foam::cellCellStencil::globalCellCells
     }
 }
 
+void Foam::cellCellStencil::seedCell
+(
+    const label celli,
+    const scalar wantedFraction,
+    bitSet& isFront,
+    scalarField& fraction
+) const
+{
+    const cell& cFaces = mesh_.cells()[celli];
+    forAll(cFaces, i)
+    {
+        const label nbrFacei = cFaces[i];
+        if (fraction[nbrFacei] < wantedFraction)
+        {
+            fraction[nbrFacei] = wantedFraction;
+            isFront.set(nbrFacei);
+        }
+    }
+}
+
+
+void Foam::cellCellStencil::setUpFront
+(
+    const labelList& allCellTypes,
+    bitSet& isFront
+) const
+{
+    const labelList& own = mesh_.faceOwner();
+    const labelList& nei = mesh_.faceNeighbour();
+
+    for (label facei = 0; facei < mesh_.nInternalFaces(); ++facei)
+    {
+        const label ownType = allCellTypes[own[facei]];
+        const label neiType = allCellTypes[nei[facei]];
+        if
+        (
+            (ownType == HOLE && neiType != HOLE)
+         || (ownType != HOLE && neiType == HOLE)
+        )
+        {
+            isFront.set(facei);
+        }
+    }
+
+    labelList nbrCellTypes;
+    syncTools::swapBoundaryCellList(mesh_, allCellTypes, nbrCellTypes);
+
+    for
+    (
+        label facei = mesh_.nInternalFaces();
+        facei < mesh_.nFaces();
+        facei++
+    )
+    {
+        const label ownType = allCellTypes[own[facei]];
+        const label neiType = nbrCellTypes[facei-mesh_.nInternalFaces()];
+
+        if
+        (
+            (ownType == HOLE && neiType != HOLE)
+         || (ownType != HOLE && neiType == HOLE)
+        )
+        {
+            isFront.set(facei);
+        }
+    }
+}
+
+
+void Foam::cellCellStencil::setUpFrontOnOversetPatch
+(
+    const labelList& allCellTypes,
+    bitSet& isFront
+) const
+{
+    const fvBoundaryMesh& fvm = mesh_.boundary();
+    // 'overset' patches
+    forAll(fvm, patchi)
+    {
+        if (isA<oversetFvPatch>(fvm[patchi]))
+        {
+            const labelList& fc = fvm[patchi].faceCells();
+            forAll(fc, i)
+            {
+                const 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.set(fvm[patchi].start()+i);
+                }
+            }
+        }
+    }
+}
+
+
+void Foam::cellCellStencil::walkFront
+(
+    const globalIndex& globalCells,
+    const label layerRelax,
+    const labelListList& allStencil,
+    labelList& allCellTypes,
+    scalarField& allWeight,
+    const scalarList& compactCellVol,
+    const labelListList& compactStencil,
+    const labelList& zoneID,
+    const label holeLayers,
+    const label useLayer
+) const
+{
+    if (useLayer > holeLayers)
+    {
+        FatalErrorInFunction<< "useLayer: " << useLayer
+            << " is larger than : " <<  holeLayers
+            << abort(FatalError);
+    }
+
+    if (useLayer == 0)
+    {
+        FatalErrorInFunction<< "useLayer: " << useLayer
+            << " can not be zero."
+            << abort(FatalError);
+    }
+
+    // Current front
+    bitSet isFront(mesh_.nFaces());
+
+    // List of cellTYpes
+    DynamicList<labelList> dataSet;
+    // List of cellTypes for increasing layers
+    DynamicList<scalarField> dAlllWeight;
+    // List of average volumen ration interpolated/donor
+    DynamicList<scalar> daverageVolRatio;
+
+    // Counting holes of different layers
+    DynamicList<label> dHoles;
+
+    const scalarField& V = mesh_.V();
+
+    // Current layer
+    labelField nLayer(allCellTypes.size(), 0);
+
+    for (label currLayer = 1; currLayer < holeLayers+1; ++currLayer)
+    {
+        // Set up original front
+        isFront = false;
+
+        setUpFront(allCellTypes, isFront);
+
+        labelList allCellTypesWork(allCellTypes);
+
+        bitSet isFrontWork(isFront);
+        label nCurrLayer = currLayer;
+
+        while (nCurrLayer > 1 && returnReduce(isFrontWork.any(), orOp<bool>()))
+        {
+            bitSet newIsFront(mesh_.nFaces());
+            forAll(isFrontWork, facei)
+            {
+                if (isFrontWork.test(facei))
+                {
+                    const label own = mesh_.faceOwner()[facei];
+
+                    if (allCellTypesWork[own] != HOLE)
+                    {
+                        allCellTypesWork[own] = HOLE;
+                        newIsFront.set(mesh_.cells()[own]);
+                    }
+
+                    if (mesh_.isInternalFace(facei))
+                    {
+                        const label nei = mesh_.faceNeighbour()[facei];
+                        if (allCellTypesWork[nei] != HOLE)
+                        {
+                            allCellTypesWork[nei] = HOLE;
+                            newIsFront.set(mesh_.cells()[nei]);
+                        }
+                    }
+                }
+            }
+            syncTools::syncFaceList(mesh_, newIsFront, orEqOp<unsigned int>());
+
+            isFrontWork.transfer(newIsFront);
+
+            nCurrLayer--;
+        }
+
+
+        if ((debug&2) && (mesh_.time().outputTime()))
+        {
+            tmp<volScalarField> tfld
+            (
+                createField
+                (
+                    mesh_,
+                    "allCellTypesWork_Holes" + name(currLayer),
+                    allCellTypesWork
+                )
+            );
+            tfld().write();
+        }
+
+        if (currLayer == 1)
+        {
+            setUpFrontOnOversetPatch(allCellTypes, isFront);
+        }
+
+        if (currLayer > 1)
+        {
+            isFront = false;
+            setUpFrontOnOversetPatch(allCellTypesWork, isFront);
+            setUpFront(allCellTypesWork, isFront);
+        }
+
+        // Current interpolation fraction
+        scalarField fraction(mesh_.nFaces(), Zero);
+
+        // Ratio between Inter/donor
+        scalarField volRatio(allCellTypes.size(), Zero);
+
+        forAll(isFront, facei)
+        {
+            if (isFront.test(facei))
+            {
+                fraction[facei] = 1.0;
+            }
+        }
+
+        scalarField allWeightWork(allCellTypes.size(), Zero);
+        bitSet nHoles(allCellTypes.size());
+
+        while (returnReduce(isFront.any(), orOp<bool>()))
+        {
+            // Interpolate cells on front
+            bitSet newIsFront(mesh_.nFaces());
+            scalarField newFraction(fraction);
+
+            forAll(isFront, facei)
+            {
+                if (isFront.test(facei))
+                {
+                    const label own = mesh_.faceOwner()[facei];
+
+                    if (allCellTypesWork[own] != HOLE)
+                    {
+                        if (allWeightWork[own] < fraction[facei])
+                        {
+                            // Cell wants to become interpolated (if sufficient
+                            // stencil, otherwise becomes hole)
+                            if (allStencil[own].size())
+                            {
+                                nLayer[own] = currLayer;
+
+                                allWeightWork[own] = fraction[facei];
+                                allCellTypesWork[own] = INTERPOLATED;
+
+                                const label donorId = compactStencil[own][0];
+
+                                volRatio[own] = V[own]/compactCellVol[donorId];
+
+                                seedCell
+                                (
+                                    own,
+                                    fraction[facei]-layerRelax,
+                                    newIsFront,
+                                    newFraction
+                                );
+
+                                nHoles[own] = true;
+                            }
+                            else
+                            {
+                                allWeightWork[own] = 0.0;
+                                allCellTypesWork[own] = HOLE;
+                                // Add faces of cell as new front
+                                seedCell
+                                (
+                                    own,
+                                    1.0,
+                                    newIsFront,
+                                    newFraction
+                                );
+                            }
+                        }
+                    }
+                    if (mesh_.isInternalFace(facei))
+                    {
+                        label nei = mesh_.faceNeighbour()[facei];
+                        if (allCellTypesWork[nei] != HOLE)
+                        {
+                            if (allWeightWork[nei] < fraction[facei])
+                            {
+                                if (allStencil[nei].size())
+                                {
+                                    nLayer[nei] = currLayer;
+
+                                    allWeightWork[nei] = fraction[facei];
+                                    allCellTypesWork[nei] = INTERPOLATED;
+
+                                    const label donorId =  compactStencil[nei][0];
+
+                                    volRatio[nei] = V[nei]/compactCellVol[donorId];
+
+                                    seedCell
+                                    (
+                                        nei,
+                                        fraction[facei]-layerRelax,
+                                        newIsFront,
+                                        newFraction
+                                    );
+                                }
+                                else
+                                {
+                                    allWeightWork[nei] = 0.0;
+                                    allCellTypesWork[nei] = HOLE;
+                                    nHoles[nei] = true;
+                                    seedCell
+                                    (
+                                        nei,
+                                        1.0,
+                                        newIsFront,
+                                        newFraction
+                                    );
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+
+            syncTools::syncFaceList(mesh_, newIsFront, orEqOp<unsigned int>());
+            syncTools::syncFaceList(mesh_, newFraction, maxEqOp<scalar>());
+
+            isFront.transfer(newIsFront);
+            fraction.transfer(newFraction);
+        }
+
+        if ((debug&2) && (mesh_.time().outputTime()))
+        {
+            tmp<volScalarField> tfld
+            (
+                createField
+                (
+                    mesh_,
+                    "allCellTypesWork_Layers" + name(currLayer),
+                    allCellTypesWork
+                )
+            );
+            tfld().write();
+        }
+
+        dAlllWeight.append(allWeightWork);
+        dataSet.append(allCellTypesWork);
+
+        // Counting interpolated cells
+        label count = 1;
+        forAll (volRatio, i)
+        {
+            if (volRatio[i] > 0)
+            {
+                count++;
+            }
+        }
+        label nCount(count);
+        reduce(nCount, sumOp<label>());
+
+        const scalar aveVol = mag(gSum(volRatio));
+
+        daverageVolRatio.append(aveVol/nCount);
+
+        // Check holes number. A sudden increase occurs when the walk leaks
+        // out from the obstacle
+        label nTotalHoles(nHoles.count());
+        reduce(nTotalHoles, sumOp<label>());
+        dHoles.append(nTotalHoles);
+
+        // Check the increase between this layer and the last one
+        // if over 50% breaks the layer loop.
+        if
+        (
+            (currLayer > 1)
+          & (nTotalHoles > 2.0*dHoles[currLayer - 1])
+        )
+        {
+            break;
+        }
+    }
+
+
+    if ((debug&2) && (mesh_.time().outputTime()))
+    {
+        tmp<volScalarField> tfld
+        (
+            createField(mesh_, "walkFront_layers", nLayer)
+        );
+        tfld().write();
+    }
+
+    // Try to find the best averageVolRatio the further from the initial set
+    // As this one is next to HOLES
+    scalarList averageVolRatio;
+    averageVolRatio.transfer(daverageVolRatio);
+    label minVolId = findMin(averageVolRatio);
+
+    if (useLayer > -1)
+    {
+        Info<< nl << " Number of layers : " << averageVolRatio.size() << nl
+            << " Average volumetric ratio : " << averageVolRatio << nl
+            << " Number of holes cells : " << dHoles << nl
+            << " Using layer : " << useLayer <<  nl << endl;
+    }
+
+    if (useLayer != -1)
+    {
+        minVolId = useLayer - 1;
+    }
+
+    allCellTypes.transfer(dataSet[minVolId]);
+    allWeight.transfer(dAlllWeight[minVolId]);
+}
+
+
+// * * * * * * * * * * * * * Ostream operator  * * * * * * * * * * * * * * * //
+
+template<>
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const InfoProxy<cellCellStencil>& ic
+)
+{
+    const cellCellStencil& e = ic.t_;
+
+    const labelUList& cellTypes = e.cellTypes();
+    const labelUList& interpolationCells = e.interpolationCells();
+    const labelListList& cellStencil = e.cellStencil();
+
+    labelList nCells(cellCellStencil::count(4, cellTypes));
+
+    label nInvalidInterpolated = 0;
+    label nLocal = 0;
+    label nMixed = 0;
+    label nRemote = 0;
+    forAll(interpolationCells, i)
+    {
+        label celli = interpolationCells[i];
+        const labelList& slots = cellStencil[celli];
+
+        if (slots.empty())
+        {
+            nInvalidInterpolated++;
+        }
+
+        bool hasLocal = false;
+        bool hasRemote = false;
+
+        forAll(slots, sloti)
+        {
+            if (slots[sloti] >= cellTypes.size())
+            {
+                hasRemote = true;
+            }
+            else
+            {
+                hasLocal = true;
+            }
+        }
+
+        if (hasRemote)
+        {
+            if (!hasLocal)
+            {
+                nRemote++;
+            }
+            else
+            {
+                nMixed++;
+            }
+        }
+        else if (hasLocal)
+        {
+            nLocal++;
+        }
+    }
+    reduce(nLocal, sumOp<label>());
+    reduce(nMixed, sumOp<label>());
+    reduce(nRemote, sumOp<label>());
+    reduce(nInvalidInterpolated, sumOp<label>());
+
+    Info<< "Overset analysis : nCells : "
+        << returnReduce(cellTypes.size(), sumOp<label>()) << nl
+        << incrIndent
+        << indent << "calculated   : " << nCells[cellCellStencil::CALCULATED]
+        << nl
+        << indent << "interpolated : " << nCells[cellCellStencil::INTERPOLATED]
+        << " (interpolated from local:" << nLocal
+        << "  mixed local/remote:" << nMixed
+        << "  remote:" << nRemote
+        << "  special:" << nInvalidInterpolated << ")" << nl
+        << indent << "hole         : " << nCells[cellCellStencil::HOLE] << nl
+        << indent << "special      : " << nCells[cellCellStencil::SPECIAL] << nl
+        << decrIndent << endl;
+
+    return os;
+}
 
 // ************************************************************************* //

src/overset/cellCellStencil/cellCellStencil/cellCellStencil.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2019 OpenCFD Ltd.
+    Copyright (C) 2017-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -74,7 +74,9 @@ public:
     {
         CALCULATED = 0,     // normal operation
         INTERPOLATED = 1,   // interpolated
-        HOLE = 2            // hole
+        HOLE = 2,           // hole
+        SPECIAL = 3,        // hole that is donor
+        POROUS = 4          // tag for special treatment due to lack of donors
     };
 
 
@@ -91,12 +93,12 @@ protected:
         //- Set of fields that should not be interpolated
         wordHashSet nonInterpolatedFields_;
 
+        //- Dictionary of motion control parameters
+        const dictionary dict_;
+
 
     // Protected Member Functions
 
-        //- Count occurrences (in parallel)
-        static labelList count(const label size, const labelUList& lst);
-
         //- Helper: create volScalarField for postprocessing.
         template<class Type>
         static tmp<volScalarField> createField
@@ -106,6 +108,9 @@ protected:
             const UList<Type>&
         );
 
+        //- Helper: strip off trailing _0
+        static word baseName(const word& name);
+
 
 private:
 
@@ -213,6 +218,9 @@ public:
             return zoneID(mesh_);
         }
 
+        //- Count occurrences (in parallel)
+        static labelList count(const label size, const labelUList& lst);
+
         //- Helper: create cell-cell addressing in global numbering
         static void globalCellCells
         (
@@ -223,8 +231,89 @@ public:
             labelListList& cellCells,
             pointListList& cellCellCentres
         );
+
+        //- Interpolation of acceptor cells from donor cells
+        template<class T>
+        static void interpolate
+        (
+            Field<T>& psi,
+            const fvMesh& mesh,
+            const cellCellStencil& overlap,
+            const List<scalarList>& wghts
+        );
+
+        //- Explicit interpolation of acceptor cells from donor cells. Looks
+        //- up cellCellStencil.
+        template<class T>
+        void interpolate(const fvMesh& mesh, Field<T>& psi) const;
+
+        //- Explicit interpolation of acceptor cells from donor cells with
+        //- boundary condition handling
+        template<class GeoField>
+        void interpolate(GeoField& psi) const;
+
+        //- Explicit interpolation of all registered fields. No boundary
+        //- conditions. Excludes selected fields (and their old-time fields)
+        template<class GeoField>
+        void interpolate
+        (
+            const fvMesh& mesh,
+            const wordHashSet& suppressed
+        ) const;
+
+        //- Version of correctBoundaryConditions that excludes 'overset' bcs
+        template<class GeoField, class SuppressBC>
+        static void correctBoundaryConditions(GeoField& psi);
+
+        //- Return info proxy.
+        //- Used to print stencil information to a stream
+        InfoProxy<cellCellStencil> info() const
+        {
+            return *this;
+        }
+
+        //- Surround holes with layer(s) of interpolated cells
+        void walkFront
+        (
+            const globalIndex& globalCells,
+            const label layerRelax,
+            const labelListList& allStencil,
+            labelList& allCellTypes,
+            scalarField& allWeight,
+            const scalarList& compactCellVol,
+            const labelListList& compactStencil,
+            const labelList& zoneID,
+            const label holeLayers,
+            const label useLayer
+        ) const;
+
+        //- Set up front using allCellTypes
+        void setUpFront
+        (
+            const labelList& allCellTypes,
+            bitSet& isFront
+
+        ) const;
+
+        //- Set up front on overset patches
+        void setUpFrontOnOversetPatch
+        (
+            const labelList& allCellTypes,
+            bitSet& isFront
+        ) const;
+
+        //- Seed faces of cell with wantedFraction (if higher than current)
+        void seedCell
+        (
+            const label cellI,
+            const scalar wantedFraction,
+            bitSet& isFront,
+            scalarField& fraction
+        ) const;
 };
 
+template<>
+Ostream& operator<<(Ostream& os, const InfoProxy<cellCellStencil>& ic);
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/overset/cellCellStencil/cellCellStencil/cellCellStencilTemplates.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2019 OpenCFD Ltd.
+    Copyright (C) 2018-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -29,6 +29,126 @@ License
 
 // * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * * //
 
+template<class T>
+void Foam::cellCellStencil::interpolate
+(
+    Field<T>& psi,
+    const fvMesh& mesh,
+    const cellCellStencil& overlap,
+    const List<scalarList>& wghts
+)
+{
+    const labelListList& stencil = overlap.cellStencil();
+
+    if (stencil.size() != mesh.nCells())
+    {
+        return;
+    }
+
+    const mapDistribute& map = overlap.cellInterpolationMap();
+    const labelList& cellIDs = overlap.interpolationCells();
+    const scalarList& factor = overlap.cellInterpolationWeight();
+
+    Field<T> work(psi);
+    map.mapDistributeBase::distribute(work, UPstream::msgType()+1);
+
+    forAll(cellIDs, i)
+    {
+        const label celli = cellIDs[i];
+
+        const scalarList& w = wghts[celli];
+        const labelList& nbrs = stencil[celli];
+        const scalar f = factor[celli];
+
+        if (nbrs.size() == 0 && f != 0.0)
+        {
+            FatalErrorInFunction << "problem: cell:" << celli
+                << " at:" << mesh.cellCentres()[celli]
+                << " type:" << overlap.cellTypes()[celli]
+                << " stencil:" << nbrs
+                << " factor:" << f << exit(FatalError);
+        }
+
+        T s(pTraits<T>::zero);
+        forAll(nbrs, nbrI)
+        {
+            s += w[nbrI]*work[nbrs[nbrI]];
+        }
+
+        psi[celli] = (1.0-f)*psi[celli] + f*s;
+    }
+}
+
+
+template<class T>
+void Foam::cellCellStencil::interpolate(const fvMesh& mesh, Field<T>& psi) const
+{
+    const cellCellStencil& overlap = *this;
+
+    interpolate
+    (
+        psi,
+        mesh,
+        overlap,
+        overlap.cellInterpolationWeights()
+    );
+}
+
+
+template<class GeoField>
+void Foam::cellCellStencil::interpolate(GeoField& psi) const
+{
+    interpolate
+    (
+        psi.mesh(),
+        psi.primitiveFieldRef()
+    );
+    psi.correctBoundaryConditions();
+}
+
+
+template<class GeoField>
+void Foam::cellCellStencil::interpolate
+(
+    const fvMesh& mesh,
+    const wordHashSet& suppressed
+) const
+{
+    const cellCellStencil& overlap = *this;
+
+    auto flds(mesh.lookupClass<GeoField>());
+    for (auto fldPtr : flds)
+    {
+        const word& name = fldPtr->name();
+        if (!suppressed.found(baseName(name)))
+        {
+            if (debug)
+            {
+                Pout<< "cellCellStencil::interpolate: interpolating : "
+                    << name << endl;
+            }
+
+            auto& fld = const_cast<GeoField&>(*fldPtr);
+
+            interpolate
+            (
+                fld.primitiveFieldRef(),
+                mesh,
+                overlap,
+                overlap.cellInterpolationWeights()
+            );
+        }
+        else
+        {
+            if (debug)
+            {
+                Pout<< "cellCellStencil::interpolate: skipping : " << name
+                    << endl;
+            }
+        }
+    }
+}
+
 template<class Type>
 Foam::tmp<Foam::volScalarField>
 Foam::cellCellStencil::createField
@@ -38,23 +158,20 @@ Foam::cellCellStencil::createField
     const UList<Type>& psi
 )
 {
-    tmp<volScalarField> tfld
+    auto tfld = tmp<volScalarField>::New
     (
-        new volScalarField
+        IOobject
         (
-            IOobject
+            name,
-            (
+            mesh.time().timeName(),
-                name,
-                mesh.time().timeName(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
             mesh,
-            dimensionedScalar(dimless, Zero),
+            IOobject::NO_READ,
-            zeroGradientFvPatchScalarField::typeName
+            IOobject::NO_WRITE,
-        )
+            false
+        ),
+        mesh,
+        dimensionedScalar(dimless, Zero),
+        zeroGradientFvPatchScalarField::typeName
     );
     volScalarField& fld = tfld.ref();
 
@@ -66,4 +183,41 @@ Foam::cellCellStencil::createField
 }
 
 
+template<class GeoField, class SuppressBC>
+void Foam::cellCellStencil::correctBoundaryConditions
+(
+    GeoField& psi
+)
+{
+    // Version of GeoField::correctBoundaryConditions that exclude evaluation
+    // of oversetFvPatchFields
+
+    typename GeoField::Boundary& bfld = psi.boundaryFieldRef();
+
+    const label nReq = Pstream::nRequests();
+
+    forAll(bfld, patchi)
+    {
+        if (!isA<SuppressBC>(bfld[patchi]))
+        {
+            bfld[patchi].initEvaluate(Pstream::commsTypes::nonBlocking);
+        }
+    }
+
+    // Block for any outstanding requests
+    if (Pstream::parRun())
+    {
+        Pstream::waitRequests(nReq);
+    }
+
+    forAll(bfld, patchi)
+    {
+        if (!isA<SuppressBC>(bfld[patchi]))
+        {
+            bfld[patchi].evaluate(Pstream::commsTypes::nonBlocking);
+        }
+    }
+}
+
+
 // ************************************************************************* //

src/overset/cellCellStencil/cellVolumeWeight/cellVolumeWeightCellCellStencil.C

@@ -37,7 +37,8 @@ License
 #include "oversetFvPatch.H"
 #include "zeroGradientFvPatchFields.H"
 #include "syncTools.H"
-#include "dynamicOversetFvMesh.H"
+#include "oversetFvMeshBase.H"
+#include "oversetFvPatchField.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -56,150 +57,6 @@ Foam::cellCellStencils::cellVolumeWeight::defaultOverlapTolerance_ = 1e-6;
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::cellCellStencils::cellVolumeWeight::walkFront
-(
-    const scalar layerRelax,
-    labelList& allCellTypes,
-    scalarField& allWeight
-) const
-{
-    // Current front
-    bitSet isFront(mesh_.nFaces());
-    // unused: bitSet doneCell(mesh_.nCells());
-
-    const fvBoundaryMesh& fvm = mesh_.boundary();
-
-
-    // 'overset' patches
-
-    forAll(fvm, patchI)
-    {
-        if (isA<oversetFvPatch>(fvm[patchI]))
-        {
-            //Pout<< "Storing faces on patch " << fvm[patchI].name() << endl;
-            forAll(fvm[patchI], i)
-            {
-                isFront.set(fvm[patchI].start()+i);
-            }
-        }
-    }
-
-
-    // Outside of 'hole' region
-    {
-        const labelList& own = mesh_.faceOwner();
-        const labelList& nei = mesh_.faceNeighbour();
-
-        for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
-        {
-            label ownType = allCellTypes[own[faceI]];
-            label neiType = allCellTypes[nei[faceI]];
-            if
-            (
-                 (ownType == HOLE && neiType != HOLE)
-              || (ownType != HOLE && neiType == HOLE)
-            )
-            {
-                //Pout<< "Front at face:" << faceI
-                //    << " at:" << mesh_.faceCentres()[faceI] << endl;
-                isFront.set(faceI);
-            }
-        }
-
-        labelList nbrCellTypes;
-        syncTools::swapBoundaryCellList(mesh_, allCellTypes, nbrCellTypes);
-
-        for
-        (
-            label faceI = mesh_.nInternalFaces();
-            faceI < mesh_.nFaces();
-            faceI++
-        )
-        {
-            label ownType = allCellTypes[own[faceI]];
-            label neiType = nbrCellTypes[faceI-mesh_.nInternalFaces()];
-
-            if
-            (
-                 (ownType == HOLE && neiType != HOLE)
-              || (ownType != HOLE && neiType == HOLE)
-            )
-            {
-                //Pout<< "Front at coupled face:" << faceI
-                //    << " at:" << mesh_.faceCentres()[faceI] << endl;
-                isFront.set(faceI);
-            }
-        }
-    }
-
-
-
-    // Current interpolation fraction
-    scalar fraction = 1.0;
-
-    while (fraction > SMALL && returnReduceOr(isFront.any()))
-    {
-        // Interpolate cells on front
-
-        Info<< "Front : fraction:" << fraction
-            << " size:" << returnReduce(isFront.count(), sumOp<label>())
-            << endl;
-
-        bitSet newIsFront(mesh_.nFaces());
-        forAll(isFront, faceI)
-        {
-            if (isFront.test(faceI))
-            {
-                label own = mesh_.faceOwner()[faceI];
-                if (allCellTypes[own] != HOLE)
-                {
-                    if (allWeight[own] < fraction)
-                    {
-                        allWeight[own] = fraction;
-
-                        //if (debug)
-                        //{
-                        //    Pout<< "    setting cell "
-                        //        << mesh_.cellCentres()[own]
-                        //        << " to " << fraction << endl;
-                        //}
-                        allCellTypes[own] = INTERPOLATED;
-                        newIsFront.set(mesh_.cells()[own]);
-                    }
-                }
-                if (mesh_.isInternalFace(faceI))
-                {
-                    label nei = mesh_.faceNeighbour()[faceI];
-                    if (allCellTypes[nei] != HOLE)
-                    {
-                        if (allWeight[nei] < fraction)
-                        {
-                            allWeight[nei] = fraction;
-
-                            //if (debug)
-                            //{
-                            //    Pout<< "    setting cell "
-                            //        << mesh_.cellCentres()[nei]
-                            //        << " to " << fraction << endl;
-                            //}
-
-                            allCellTypes[nei] = INTERPOLATED;
-                            newIsFront.set(mesh_.cells()[nei]);
-                        }
-                    }
-                }
-            }
-        }
-
-        syncTools::syncFaceList(mesh_, newIsFront, orEqOp<unsigned int>());
-
-        isFront.transfer(newIsFront);
-
-        fraction -= layerRelax;
-    }
-}
-
-
 void Foam::cellCellStencils::cellVolumeWeight::findHoles
 (
     const globalIndex& globalCells,
@@ -213,13 +70,6 @@ void Foam::cellCellStencils::cellVolumeWeight::findHoles
     const labelList& own = mesh.faceOwner();
     const labelList& nei = mesh.faceNeighbour();
 
-
-    // The input cellTypes will be
-    // - HOLE           : cell part covered by other-mesh patch
-    // - INTERPOLATED   : cell fully covered by other-mesh patch
-    //                    or next to 'overset' patch
-    // - CALCULATED     : otherwise
-    //
     // so we start a walk from our patches and any cell we cannot reach
     // (because we walk is stopped by other-mesh patch) is a hole.
 
@@ -503,6 +353,7 @@ void Foam::cellCellStencils::cellVolumeWeight::interpolatePatchTypes
             {
                 // 'patch' overrides -1 and 'other'
                 result[cellI] = PATCH;
+                break;
             }
             else if (result[cellI] == -1)
             {
@@ -589,27 +440,26 @@ void Foam::cellCellStencils::cellVolumeWeight::combineCellTypes
                 // Patch-patch interaction... For now disable always
                 allCellTypes[cellI] = HOLE;
                 validDonors = false;
-
                 // Alternative is to look at the amount of overlap but this
                 // is not very robust
-                //if (allCellTypes[cellI] != HOLE)
+//                 if (allCellTypes[cellI] != HOLE)
-                //{
+//                 {
-                //    scalar overlapVol = sum(weights[subCellI]);
+//                    scalar overlapVol = sum(weights[subCellI]);
-                //    scalar v = mesh_.V()[cellI];
+//                    scalar v = mesh_.V()[cellI];
-                //    if (overlapVol < (1.0-overlapTolerance_)*v)
+//                    if (overlapVol < (1.0-overlapTolerance_)*v)
-                //    {
+//                    {
-                //        //Pout<< "** Patch overlap:" << cellI
+//                        //Pout<< "** Patch overlap:" << cellI
-                //        //    << " at:" << mesh_.cellCentres()[cellI] << endl;
+//                        //    << " at:" << mesh_.cellCentres()[cellI] << endl;
-                //        allCellTypes[cellI] = HOLE;
+//                        allCellTypes[cellI] = HOLE;
-                //        validDonors = false;
+//                        validDonors = false;
-                //    }
+//                    }
-                //}
+//                 }
             }
             break;
 
             case OVERSET:
             {
-                validDonors = false;
+                validDonors = true;
             }
             break;
         }
@@ -685,6 +535,10 @@ Foam::cellCellStencils::cellVolumeWeight::cellVolumeWeight
         mesh_,
         dimensionedScalar(dimless, Zero),
         zeroGradientFvPatchScalarField::typeName
+    ),
+    allowInterpolatedDonors_
+    (
+        dict.getOrDefault("allowInterpolatedDonors", true)
     )
 {
     // Protect local fields from interpolation
@@ -701,6 +555,9 @@ Foam::cellCellStencils::cellVolumeWeight::cellVolumeWeight
     // Read zoneID
     this->zoneID();
 
+    overlapTolerance_ =
+        dict_.getOrDefault("overlapTolerance", defaultOverlapTolerance_);
+
     // Read old-time cellTypes
     IOobject io
     (
@@ -755,7 +612,6 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
     }
     Pstream::listCombineReduce(nCellsPerZone, plusEqOp<label>());
 
-
     Info<< typeName << " : detected " << nZones
         << " mesh regions" << nl << endl;
 
@@ -777,7 +633,6 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
     Info<< decrIndent;
 
 
-
     // Current best guess for cells. Includes best stencil. Weights should
     // add up to volume.
     labelList allCellTypes(mesh_.nCells(), CALCULATED);
@@ -802,6 +657,15 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
         markPatchCells(partMesh, partCellMap, allPatchTypes);
     }
 
+    if ((debug&2) && (mesh_.time().outputTime()))
+    {
+        tmp<volScalarField> tfld
+        (
+            createField(mesh_, "allPatchTypes", allPatchTypes)
+        );
+        tfld().write();
+    }
+
 
     labelList nCells(count(3, allPatchTypes));
     Info<< nl
@@ -815,7 +679,6 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
 
     globalIndex globalCells(mesh_.nCells());
 
-
     for (label srcI = 0; srcI < meshParts.size()-1; srcI++)
     {
         const fvMesh& srcMesh = meshParts[srcI].subMesh();
@@ -837,7 +700,6 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
                 false                   // do not normalise
             );
 
-
             {
                 // Get tgt patch types on src mesh
                 labelList interpolatedTgtPatchTypes(srcMesh.nCells(), -1);
@@ -863,6 +725,8 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
                         mapper.tgtMap()->distribute(tgtGlobalCells);
                     }
                 }
+
+
                 combineCellTypes
                 (
                     srcI,
@@ -930,13 +794,20 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
     }
 
 
-    if (debug)
+    if ((debug&2) && (mesh_.time().outputTime()))
     {
         tmp<volScalarField> tfld
         (
             createField(mesh_, "allCellTypes", allCellTypes)
         );
         tfld().write();
+
+        tmp<volScalarField> tdonors
+        (
+            createField(mesh_, "allDonorID", allDonorID)
+        );
+        tdonors().write();
+
     }
 
 
@@ -958,8 +829,6 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
                     }
                     else
                     {
-                        //Pout<< "Holeing interpolated cell:" << cellI
-                        //    << " at:" << mesh_.cellCentres()[cellI] << endl;
                         allCellTypes[cellI] = HOLE;
                         allWeights[cellI].clear();
                         allStencil[cellI].clear();
@@ -970,36 +839,97 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
         }
     }
 
-/*
+
-    // Knock out cell with insufficient interpolation weights
+    if ((debug&2) && (mesh_.time().outputTime()))
-    forAll(allCellTypes, cellI)
-    {
-        if (allCellTypes[cellI] == INTERPOLATED)
-        {
-            scalar v = mesh_.V()[cellI];
-            scalar overlapVol = sum(allWeights[cellI]);
-            if (overlapVol < (1.0-overlapTolerance_)*v)
-            {
-                //Pout<< "Holeing cell:" << cellI
-                //    << " at:" << mesh_.cellCentres()[cellI] << endl;
-                allCellTypes[cellI] = HOLE;
-                allWeights[cellI].clear();
-                allStencil[cellI].clear();
-            }
-        }
-    }
-*/
-    if (debug)
     {
         tmp<volScalarField> tfld
         (
             createField(mesh_, "allCellTypes_patch", allCellTypes)
         );
-        //tfld.ref().correctBoundaryConditions();
         tfld().write();
+
+        labelList stencilSize(mesh_.nCells());
+        forAll(allStencil, celli)
+        {
+            stencilSize[celli] = allStencil[celli].size();
+        }
+        tmp<volScalarField> tfldStencil
+        (
+            createField(mesh_, "allStencil_patch", stencilSize)
+        );
+        tfldStencil().write();
     }
 
 
+    // Mark unreachable bits
+    findHoles(globalCells, mesh_, zoneID, allStencil, allCellTypes);
+
+    if ((debug&2) && (mesh_.time().outputTime()))
+    {
+        tmp<volScalarField> tfld
+        (
+            createField(mesh_, "allCellTypes_hole", allCellTypes)
+        );
+        tfld().write();
+
+        labelList stencilSize(mesh_.nCells());
+        forAll(allStencil, celli)
+        {
+            stencilSize[celli] = allStencil[celli].size();
+        }
+        tmp<volScalarField> tfldStencil
+        (
+            createField(mesh_, "allStencil_hole", stencilSize)
+        );
+        tfldStencil().write();
+    }
+
+
+    // Add buffer interpolation layer around holes
+    scalarField allWeight(mesh_.nCells(), Zero);
+
+    labelListList compactStencil(allStencil);
+    List<Map<label>> compactStencilMap;
+    mapDistribute map(globalCells, compactStencil, compactStencilMap);
+
+    scalarList compactCellVol(mesh_.V());
+    map.distribute(compactCellVol);
+
+    walkFront
+    (
+        globalCells,
+        layerRelax,
+        allStencil,
+        allCellTypes,
+        allWeight,
+        compactCellVol,
+        compactStencil,
+        zoneID,
+        dict_.getOrDefault("holeLayers", 1),
+        dict_.getOrDefault("useLayer", -1)
+    );
+
+
+    if ((debug&2) && (mesh_.time().outputTime()))
+    {
+        tmp<volScalarField> tfld
+        (
+            createField(mesh_, "allCellTypes_front", allCellTypes)
+        );
+        tfld().write();
+
+        labelList stencilSize(mesh_.nCells());
+        forAll(allStencil, celli)
+        {
+            stencilSize[celli] = allStencil[celli].size();
+        }
+        tmp<volScalarField> tfldStencil
+        (
+            createField(mesh_, "allStencil_front", stencilSize)
+        );
+        tfldStencil().write();
+    }
+
     // Check previous iteration cellTypes_ for any hole->calculated changes
     // If so set the cell either to interpolated (if there are donors) or
     // holes (if there are no donors). Note that any interpolated cell might
@@ -1035,53 +965,29 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
     }
 
 
-    // Mark unreachable bits
+    labelList compactCellTypes(allCellTypes);
-    findHoles(globalCells, mesh_, zoneID, allStencil, allCellTypes);
+    map.distribute(compactCellTypes);
 
-    if (debug)
+    label nHoleDonors = 0;
-    {
-        tmp<volScalarField> tfld
-        (
-            createField(mesh_, "allCellTypes_hole", allCellTypes)
-        );
-        //tfld.ref().correctBoundaryConditions();
-        tfld().write();
-    }
-
-
-    // Add buffer interpolation layer around holes
-    scalarField allWeight(mesh_.nCells(), Zero);
-    walkFront(layerRelax, allCellTypes, allWeight);
-
-    if (debug)
-    {
-        tmp<volScalarField> tfld
-        (
-            createField(mesh_, "allCellTypes_front", allCellTypes)
-        );
-        //tfld.ref().correctBoundaryConditions();
-        tfld().write();
-    }
-
-    // Normalise weights, Clear storage
     forAll(allCellTypes, cellI)
     {
         if (allCellTypes[cellI] == INTERPOLATED)
         {
-            if (allWeight[cellI] < SMALL || allStencil[cellI].size() == 0)
+            const labelList& slots = compactStencil[cellI];
+            forAll (slots, subCellI)
             {
-                //Pout<< "Clearing cell:" << cellI
+                if
-                //    << " at:" << mesh_.cellCentres()[cellI] << endl;
+                (
-                allWeights[cellI].clear();
+                    compactCellTypes[slots[0]] == HOLE
-                allStencil[cellI].clear();
+                 ||
-                allWeight[cellI] = 0.0;
+                    (
-            }
+                       !allowInterpolatedDonors_
-            else
+                     && compactCellTypes[slots[0]] == INTERPOLATED
-            {
+                    )
-                scalar s = sum(allWeights[cellI]);
+                )
-                forAll(allWeights[cellI], i)
                 {
-                    allWeights[cellI][i] /= s;
+                    allWeights[cellI][subCellI] = 0;
+                    nHoleDonors++;
                 }
             }
         }
@@ -1091,112 +997,52 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
             allStencil[cellI].clear();
         }
     }
+    reduce(nHoleDonors, sumOp<label>());
 
 
+    // Normalize weights
+    forAll(allCellTypes, cellI)
+    {
+        if (allCellTypes[cellI] == INTERPOLATED)
+        {
+            const scalar s = sum(allWeights[cellI]);
+
+            if (s < SMALL)
+            {
+                allCellTypes[cellI] = POROUS;
+                allWeights[cellI].clear();
+                allStencil[cellI].clear();
+            }
+            else
+            {
+                forAll(allWeights[cellI], i)
+                {
+                    allWeights[cellI][i] /= s;
+                }
+            }
+        }
+    }
+
     // Write to volField for debugging
-    if (debug)
+    if ((debug&2) && (mesh_.time().outputTime()))
     {
-        volScalarField patchTypes
+        if ((debug&2) && (mesh_.time().outputTime()))
-        (
-            IOobject
-            (
-                "patchTypes",
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            mesh_,
-            dimensionedScalar(dimless, Zero),
-            zeroGradientFvPatchScalarField::typeName
-        );
-
-        forAll(patchTypes.internalField(), cellI)
         {
-            patchTypes[cellI] = allPatchTypes[cellI];
+            tmp<volScalarField> tfld
-        }
-        //patchTypes.correctBoundaryConditions();
-        dynamicOversetFvMesh::correctBoundaryConditions
-        <
-            volScalarField,
-            oversetFvPatchField<scalar>
-        >(patchTypes.boundaryFieldRef(), false);
-        patchTypes.write();
-    }
-    if (debug)
-    {
-        volScalarField volTypes
-        (
-            IOobject
             (
-                "cellTypes",
+                createField(mesh_, "allCellTypes_final", allCellTypes)
-                mesh_.time().timeName(),
+            );
-                mesh_,
+            tfld().write();
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            mesh_,
-            dimensionedScalar(dimless, Zero),
-            zeroGradientFvPatchScalarField::typeName
-        );
-
-        forAll(volTypes.internalField(), cellI)
-        {
-            volTypes[cellI] = allCellTypes[cellI];
         }
-        //volTypes.correctBoundaryConditions();
-        dynamicOversetFvMesh::correctBoundaryConditions
-        <
-            volScalarField,
-            oversetFvPatchField<scalar>
-        >(volTypes.boundaryFieldRef(), false);
-        volTypes.write();
     }
 
 
-//  // Check previous iteration cellTypes_ for any hole->calculated changes
-//  {
-//      label nCalculated = 0;
-//
-//      forAll(cellTypes_, celli)
-//      {
-//          if (allCellTypes[celli] == CALCULATED && cellTypes_[celli] == HOLE)
-//          {
-//              if (allStencil[celli].size() == 0)
-//              {
-//                  FatalErrorInFunction
-//                      << "Cell:" << celli
-//                      << " at:" << mesh_.cellCentres()[celli]
-//                      << " zone:" << zoneID[celli]
-//                      << " changed from hole to calculated"
-//                      << " but there is no donor"
-//                      << exit(FatalError);
-//              }
-//              else
-//              {
-//                  allCellTypes[celli] = INTERPOLATED;
-//                  nCalculated++;
-//              }
-//          }
-//      }
-//
-//      if (debug)
-//      {
-//          Pout<< "Detected " << nCalculated << " cells changing from hole"
-//              << " to calculated. Changed these to interpolated"
-//              << endl;
-//      }
-//  }
-
-
     cellTypes_.transfer(allCellTypes);
     cellStencil_.transfer(allStencil);
     cellInterpolationWeights_.transfer(allWeights);
     cellInterpolationWeight_.transfer(allWeight);
     //cellInterpolationWeight_.correctBoundaryConditions();
-    dynamicOversetFvMesh::correctBoundaryConditions
+    oversetFvMeshBase::correctBoundaryConditions
     <
         volScalarField,
         oversetFvPatchField<scalar>
@@ -1220,7 +1066,7 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
     );
 
     // Dump interpolation stencil
-    if (debug)
+    if ((debug&2) && (mesh_.time().outputTime()))
     {
         // Dump weight
         cellInterpolationWeight_.instance() = mesh_.time().timeName();
@@ -1247,15 +1093,16 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
 
             forAll(slots, i)
             {
-                const point& donorCc = cc[slots[i]];
+                if (cellInterpolationWeights_[cellI][slots[i]] > 0)
-                const point& accCc = mesh_.cellCentres()[cellI];
+                {
-
+                    const point& donorCc = cc[slots[i]];
-                str.writeLine(accCc, 0.1*accCc+0.9*donorCc);
+                    const point& accCc = mesh_.cellCentres()[cellI];
+                    str.writeLine(accCc, 0.1*accCc+0.9*donorCc);
+                }
             }
         }
     }
 
-
     {
         labelList nCells(count(3, cellTypes_));
         Info<< "Overset analysis : nCells : "
@@ -1267,7 +1114,6 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
             << decrIndent << endl;
     }
 
-    // Tbd: detect if anything changed. Most likely it did!
     return true;
 }
 

src/overset/cellCellStencil/cellVolumeWeight/cellVolumeWeightCellCellStencil.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2018 OpenCFD Ltd.
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -90,15 +90,11 @@ protected:
         //- Amount of interpolation
         volScalarField cellInterpolationWeight_;
 
+        //- Allow interpolared as donors
+        const bool allowInterpolatedDonors_;
 
    // Protected Member Functions
 
-        void walkFront
-        (
-            const scalar layerRelax,
-            labelList& allCellTypes,
-            scalarField& allWeight
-        ) const;
 
         //- Find cells next to cells of type PATCH
         void findHoles

src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2020 OpenCFD Ltd.
+    Copyright (C) 2017-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -69,6 +69,15 @@ class inverseDistance
 :
     public cellCellStencil
 {
+    // Private Member Functions
+
+        //- No copy construct
+        inverseDistance(const inverseDistance&) = delete;
+
+        //- No copy assignment
+        void operator=(const inverseDistance&) = delete;
+
+
 protected:
 
     // Protected data
@@ -76,6 +85,12 @@ protected:
         //- Dictionary of motion control parameters
         const dictionary dict_;
 
+        //- Allow holes as donors
+        const bool allowHoleDonors_;
+
+        //- Allow interpolared as donors
+        const bool allowInterpolatedDonors_;
+
         //- Small perturbation vector for geometric tests
         vector smallVec_;
 
@@ -97,6 +112,13 @@ protected:
         //- Amount of interpolation
         volScalarField cellInterpolationWeight_;
 
+        //- Data for a set of interpolated/donor set
+        struct interpolatedDonorSet
+        {
+            scalar averVolRatio;
+            labelList cellTypes;
+        };
+
 
    // Protected Member Functions
 
@@ -217,6 +239,7 @@ protected:
             const boolList& isBlockedFace,
             labelList& cellRegion
         ) const;
+
         autoPtr<globalIndex> compactedRegionSplit
         (
             const fvMesh& mesh,
@@ -244,28 +267,18 @@ protected:
             scalarField& fraction
         ) const;
 
-        //- Surround holes with layer(s) of interpolated cells
-        void walkFront
-        (
-            const scalar layerRelax,
-            const labelListList& allStencil,
-            labelList& allCellTypes,
-            scalarField& allWeight
-        ) const;
-
         //- Create stencil starting from the donor containing the acceptor
-        virtual void createStencil(const globalIndex&);
+        //- allowHoleDonors : allow donors of type HOLE (otherwise are filtered
+        //- out)
+        virtual void createStencil
+        (
+            const globalIndex&,
+            const bool allowHoleDonors
+        );
 
-
+        //- Make holes next to live ones type SPECIAL and include in
-private:
+        //- interpolation stencil
-
+        void holeExtrapolationStencil(const globalIndex& globalCells);
-    // Private Member Functions
-
-        //- No copy construct
-        inverseDistance(const inverseDistance&) = delete;
-
-        //- No copy assignment
-        void operator=(const inverseDistance&) = delete;
 
 
 public:

src/overset/cellCellStencil/trackingInverseDistance/trackingInverseDistanceCellCellStencil.C

@@ -37,6 +37,7 @@ License
 #include "treeBoundBoxList.H"
 #include "voxelMeshSearch.H"
 #include "dynamicOversetFvMesh.H"
+#include "waveMethod.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -337,17 +338,24 @@ void Foam::cellCellStencils::trackingInverseDistance::markDonors
 
     const fvMesh& srcMesh = meshParts_[srcI].subMesh();
     const labelList& srcCellMap = meshParts_[srcI].cellMap();
-    const voxelMeshSearch& meshSearch = meshSearches[srcI];
     const fvMesh& tgtMesh = meshParts_[tgtI].subMesh();
     const pointField& tgtCc = tgtMesh.cellCentres();
     const labelList& tgtCellMap = meshParts_[tgtI].cellMap();
 
     // 1. do processor-local src/tgt overlap
     {
+        labelList tgtToSrcAddr;
+        waveMethod::calculate(tgtMesh, srcMesh, tgtToSrcAddr);
+
         forAll(tgtCellMap, tgtCelli)
         {
-            label srcCelli = meshSearch.findCell(tgtCc[tgtCelli]);
+            const label srcCelli = tgtToSrcAddr[tgtCelli];
-            if (srcCelli != -1 && allCellTypes[srcCellMap[srcCelli]] != HOLE)
+
+            if
+            (
+                srcCelli != -1
+             && allCellTypes[tgtCellMap[tgtCelli]] != HOLE
+            )
             {
                 label celli = tgtCellMap[tgtCelli];
 
@@ -444,8 +452,9 @@ void Foam::cellCellStencils::trackingInverseDistance::markDonors
         labelList donors(samples.size(), -1);
         forAll(samples, sampleI)
         {
-            label srcCelli = meshSearch.findCell(samples[sampleI]);
+            const point& sample = samples[sampleI];
-            if (srcCelli != -1 && allCellTypes[srcCellMap[srcCelli]] != HOLE)
+            label srcCelli = srcMesh.findCell(sample, polyMesh::CELL_TETS);
+            if (srcCelli != -1)
             {
                 donors[sampleI] = globalCells_.toGlobal(srcCellMap[srcCelli]);
             }
@@ -736,7 +745,7 @@ bool Foam::cellCellStencils::trackingInverseDistance::update()
     }
 
 
-    if (debug)
+    if (debug & 2)
     {
         forAll(patchParts, zonei)
         {
@@ -832,13 +841,19 @@ bool Foam::cellCellStencils::trackingInverseDistance::update()
     DebugInfo<< FUNCTION_NAME << " : Determined holes and donor-acceptors"
         << endl;
 
-    if (debug)
+    if ((debug & 2) && (mesh_.time().outputTime()))
     {
         tmp<volScalarField> tfld
         (
             createField(mesh_, "allCellTypes", allCellTypes)
         );
         tfld().write();
+
+        tmp<volScalarField> tallDonorID
+        (
+            createField(mesh_, "allDonorID", allDonorID)
+        );
+        tallDonorID().write();
     }
 
 
@@ -867,15 +882,86 @@ bool Foam::cellCellStencils::trackingInverseDistance::update()
     }
     DebugInfo<< FUNCTION_NAME << " : Removed bad donors" << endl;
 
-    if (debug)
+    if ((debug & 2) && (mesh_.time().outputTime()))
     {
         tmp<volScalarField> tfld
         (
             createField(mesh_, "allCellTypes_patch", allCellTypes)
         );
         tfld().write();
+
+        tmp<volScalarField> tfldOld
+        (
+            createField(mesh_, "allCellTypes_old", cellTypes_)
+        );
+        tfldOld().write();
     }
 
+    // Mark unreachable bits
+    findHoles(globalCells_, mesh_, zoneID, allStencil, allCellTypes);
+    DebugInfo<< FUNCTION_NAME << " : Flood-filled holes" << endl;
+
+    if ((debug & 2) && (mesh_.time().outputTime()))
+    {
+        tmp<volScalarField> tfld
+        (
+            createField(mesh_, "allCellTypes_hole", allCellTypes)
+        );
+        tfld().write();
+
+        labelList stencilSize(mesh_.nCells());
+        forAll(allStencil, celli)
+        {
+            stencilSize[celli] = allStencil[celli].size();
+        }
+        tmp<volScalarField> tfldsten
+        (
+            createField(mesh_, "allStencil_hole", stencilSize)
+        );
+        tfldsten().write();
+    }
+
+    // Update allStencil with new fill HOLES
+    forAll(allCellTypes, celli)
+    {
+        if (allCellTypes[celli] == HOLE && allStencil[celli].size())
+        {
+            allStencil[celli].clear();
+        }
+    }
+
+    // Add buffer interpolation layer(s) around holes
+    scalarField allWeight(mesh_.nCells(), Zero);
+
+    labelListList compactStencil(allStencil);
+    List<Map<label>> compactStencilMap;
+    mapDistribute map(globalCells, compactStencil, compactStencilMap);
+
+    scalarList compactCellVol(mesh_.V());
+    map.distribute(compactCellVol);
+
+    walkFront
+    (
+        globalCells,
+        layerRelax,
+        allStencil,
+        allCellTypes,
+        allWeight,
+        compactCellVol,
+        compactStencil,
+        zoneID,
+        dict_.getOrDefault("holeLayers", 1),
+        dict_.getOrDefault("useLayer", -1)
+    );
+
+    if ((debug & 2) && (mesh_.time().outputTime()))
+    {
+        tmp<volScalarField> tfld
+        (
+            createField(mesh_, "allCellTypes_front", allCellTypes)
+        );
+        tfld().write();
+    }
 
     // Check previous iteration cellTypes_ for any hole->calculated changes
     // If so set the cell either to interpolated (if there are donors) or
@@ -911,48 +997,47 @@ bool Foam::cellCellStencils::trackingInverseDistance::update()
     }
 
 
-    // Mark unreachable bits
-    findHoles(globalCells_, mesh_, zoneID, allStencil, allCellTypes);
-    DebugInfo<< FUNCTION_NAME << " : Flood-filled holes" << endl;
-
-    if (debug)
-    {
-        tmp<volScalarField> tfld
-        (
-            createField(mesh_, "allCellTypes_hole", allCellTypes)
-        );
-        tfld().write();
-    }
-
-    // Add buffer interpolation layer(s) around holes
-    scalarField allWeight(mesh_.nCells(), Zero);
-    walkFront(layerRelax, allStencil, allCellTypes, allWeight);
-    DebugInfo<< FUNCTION_NAME << " : Implemented layer relaxation" << endl;
-
-    if (debug)
-    {
-        tmp<volScalarField> tfld
-        (
-            createField(mesh_, "allCellTypes_front", allCellTypes)
-        );
-        tfld().write();
-    }
-
-
     // Convert cell-cell addressing to stencil in compact notation
 
     cellTypes_.transfer(allCellTypes);
     cellStencil_.setSize(mesh_.nCells());
     cellInterpolationWeights_.setSize(mesh_.nCells());
     DynamicList<label> interpolationCells;
+
+    labelList compactCellTypes(cellTypes_);
+    map.distribute(compactCellTypes);
+
+    label nHoleDonors = 0;
     forAll(cellTypes_, celli)
     {
         if (cellTypes_[celli] == INTERPOLATED)
         {
-            cellStencil_[celli].transfer(allStencil[celli]);
+            const labelList& slots = compactStencil[celli];
-            cellInterpolationWeights_[celli].setSize(1);
+            if (slots.size())
-            cellInterpolationWeights_[celli][0] = 1.0;
+            {
-            interpolationCells.append(celli);
+                if
+                (
+                    compactCellTypes[slots[0]] == HOLE
+                 ||
+                    (
+                       !allowInterpolatedDonors_
+                     && compactCellTypes[slots[0]] == INTERPOLATED
+                    )
+                )
+                {
+                    cellTypes_[celli] = POROUS;
+                    cellStencil_[celli].clear();
+                    cellInterpolationWeights_[celli].clear();
+                    nHoleDonors++;
+                }
+                else
+                {
+                    cellStencil_[celli].transfer(allStencil[celli]);
+                    cellInterpolationWeights_[celli].setSize(1);
+                    cellInterpolationWeights_[celli][0] = 1.0;
+                    interpolationCells.append(celli);
+                }
+            }
         }
         else
         {
@@ -960,6 +1045,11 @@ bool Foam::cellCellStencils::trackingInverseDistance::update()
             cellInterpolationWeights_[celli].clear();
         }
     }
+
+    reduce(nHoleDonors, sumOp<label>());
+
+    DebugInfo<< FUNCTION_NAME << "nHole Donors : " << nHoleDonors << endl;
+
     interpolationCells_.transfer(interpolationCells);
 
     List<Map<label>> compactMap;
@@ -973,16 +1063,18 @@ bool Foam::cellCellStencils::trackingInverseDistance::update()
         )
     );
     cellInterpolationWeight_.transfer(allWeight);
-    //cellInterpolationWeight_.correctBoundaryConditions();
+    oversetFvMeshBase::correctBoundaryConditions
-    dynamicOversetFvMesh::correctBoundaryConditions
     <
         volScalarField,
         oversetFvPatchField<scalar>
     >(cellInterpolationWeight_.boundaryFieldRef(), false);
 
 
-    if (debug & 2)
+    if ((debug & 2) && (mesh_.time().outputTime()))
     {
+        // Dump mesh
+        mesh_.time().write();
+
         // Dump stencil
         mkDir(mesh_.time().timePath());
         OBJstream str(mesh_.time().timePath()/"injectionStencil.obj");
@@ -1010,45 +1102,23 @@ bool Foam::cellCellStencils::trackingInverseDistance::update()
 
 
     // Extend stencil to get inverse distance weighted neighbours
-    createStencil(globalCells);
+    createStencil(globalCells, allowHoleDonors_);
     DebugInfo<< FUNCTION_NAME << " : Extended stencil" << endl;
 
+    // Optional: convert hole cells next to non-hole cells into
+    // interpolate-from-neighbours (of cell type SPECIAL)
+    if (allowHoleDonors_)
+    {
+        holeExtrapolationStencil(globalCells_);
+    }
 
-    if (debug & 2)
+
+    if ((debug & 2) && (mesh_.time().outputTime()))
     {
         // Dump weight
         cellInterpolationWeight_.instance() = mesh_.time().timeName();
         cellInterpolationWeight_.write();
 
-        // Dump cell types
-        volScalarField volTypes
-        (
-            IOobject
-            (
-                "cellTypes",
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            mesh_,
-            dimensionedScalar(dimless, Zero),
-            zeroGradientFvPatchScalarField::typeName
-        );
-
-        forAll(volTypes.internalField(), celli)
-        {
-            volTypes[celli] = cellTypes_[celli];
-        }
-        //volTypes.correctBoundaryConditions();
-        dynamicOversetFvMesh::correctBoundaryConditions
-        <
-            volScalarField,
-            oversetFvPatchField<scalar>
-        >(volTypes.boundaryFieldRef(), false);
-        volTypes.write();
-
         // Dump stencil
         mkDir(mesh_.time().timePath());
         OBJstream str(mesh_.time().timePath()/"stencil.obj");

src/overset/include/correctPhiFaceMask.H

@@ -0,0 +1,44 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Description
+    Correct phi on new faces C-I faces
+
+\*---------------------------------------------------------------------------*/
+
+surfaceScalarField faceMaskOld
+(
+    localMin<scalar>(mesh).interpolate(cellMask.oldTime())
+);
+
+// Zero Uf on old faceMask (H-I)
+Uf *= faceMaskOld;
+// Update Uf and phi on new C-I faces
+Uf += (1-faceMaskOld)*fvc::interpolate(U);
+phi = mesh.Sf() & Uf;
+
+phi *= faceMask;
+
+// ************************************************************************* //

src/overset/include/correctRhoPhiFaceMask.H

@@ -0,0 +1,55 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Description
+    Correct phi on new faces C-I faces
+
+\*---------------------------------------------------------------------------*/
+
+const surfaceScalarField faceMaskOld
+(
+    localMin<scalar>(mesh).interpolate(cellMask.oldTime())
+);
+
+// Zero Uf on old faceMask (H-I)
+rhoUf() *= faceMaskOld;
+
+surfaceVectorField rhoUfint(fvc::interpolate(rho*U));
+
+// Update Uf and phi on new C-I faces
+rhoUf() += (1-faceMaskOld)*rhoUfint;
+
+// Update Uf boundary
+forAll(rhoUf().boundaryField(), patchI)
+{
+    rhoUf().boundaryFieldRef()[patchI] = rhoUfint.boundaryField()[patchI];
+}
+
+// Calculate absolute flux from the mapped surface velocity
+phi = mesh.Sf() & rhoUf();
+
+phi *= faceMask;
+
+// ************************************************************************* //

src/overset/include/createCellMask.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2017 OpenCFD Ltd.
+    Copyright (C) 2016-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -50,6 +50,11 @@ volScalarField cellMask
     zeroGradientFvPatchScalarField::typeName
 );
 
+surfaceScalarField faceMask
+(
+    localMin<scalar>(mesh).interpolate(cellMask)
+);
+
 #include "setCellMask.H"
 
 // ************************************************************************* //

src/overset/include/createControls.H

@@ -0,0 +1,9 @@
+bool adjustFringe
+(
+    pimple.dict().getOrDefault("oversetAdjustPhi", false)
+);
+
+bool oversetPatchErrOutput
+(
+    pimple.dict().getOrDefault("oversetPatchErrOutput", false)
+);

src/overset/include/readOversetDyMControls.H

@@ -0,0 +1,4 @@
+adjustFringe = pimple.dict().getOrDefault("oversetAdjustPhi", false);
+label zoneIdMass =  pimple.dict().getOrDefault("adjustZone", -1);
+oversetPatchErrOutput =
+    pimple.dict().getOrDefault("oversetPatchErrOutput", false);

src/overset/include/setCellMask.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2017 OpenCFD Ltd.
+    Copyright (C) 2016-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -35,12 +35,19 @@ Description
     cellMask.primitiveFieldRef() = 1.0;
     forAll(cellMask, cellI)
     {
-        if (cellTypes[cellI] == cellCellStencil::HOLE)
+        if
+        (
+            cellTypes[cellI] == cellCellStencil::HOLE
+         || cellTypes[cellI] == cellCellStencil::SPECIAL
+        )
         {
             cellMask[cellI] = 0.0;
         }
     }
     cellMask.correctBoundaryConditions();
+
+    faceMask = localMin<scalar>(mesh).interpolate(cellMask);
 }
 
+
 // ************************************************************************* //

src/overset/oversetAdjustPhi/oversetAdjustPhi.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2018 OpenCFD Ltd.
+    Copyright (C) 2017-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -37,64 +37,62 @@ License
 bool Foam::oversetAdjustPhi
 (
     surfaceScalarField& phi,
-    const volVectorField& U
+    const volVectorField& U,
+    const label zoneId
 )
 {
     const fvMesh& mesh = U.mesh();
-
     const cellCellStencilObject& overlap = Stencil::New(mesh);
     const labelList& cellTypes = overlap.cellTypes();
     const labelList& zoneID = overlap.zoneID();
-    label nZones = gMax(zoneID)+1;
-
-
 
     // Pass1: accumulate all fluxes, calculate correction factor
 
-    scalarField massIn(nZones, Zero);
+    scalar massIn = 0;
-    scalarField adjustableMassOut(nZones, Zero);
+    scalar massOut = 0;
-
-    surfaceScalarField::Boundary& bphi =
-        phi.boundaryFieldRef();
 
+    surfaceScalarField::Boundary& bphi = phi.boundaryFieldRef();
 
     // Check all faces on the outside of interpolated cells
     const labelUList& own = mesh.owner();
     const labelUList& nei = mesh.neighbour();
+    forAll(own, facei)
     {
-        forAll(own, facei)
+        const label zonei = zoneID[own[facei]];
+        const label ownType = cellTypes[own[facei]];
+        const label neiType = cellTypes[nei[facei]];
+
+        const bool ownCalc =
+            (ownType == cellCellStencil::CALCULATED)
+         && (neiType == cellCellStencil::INTERPOLATED);
+
+        const bool neiCalc =
+            (ownType == cellCellStencil::INTERPOLATED)
+         && (neiType == cellCellStencil::CALCULATED);
+
+        const bool ownOrCalc = (ownCalc || neiCalc);
+
+        if
+        (
+            (ownOrCalc && (zonei == zoneId))
+         || (ownOrCalc && (zoneId == -1))
+        )
         {
-            label zonei = zoneID[own[facei]];   // note:own and nei in same zone
+            // Calculate flux w.r.t. calculated cell
+            scalar flux = phi[facei];
 
-            label ownType = cellTypes[own[facei]];
+            if (ownCalc)
-            label neiType = cellTypes[nei[facei]];
-
-            bool ownCalc =
-                (ownType == cellCellStencil::CALCULATED)
-             && (neiType == cellCellStencil::INTERPOLATED);
-
-            bool neiCalc =
-                (ownType == cellCellStencil::INTERPOLATED)
-             && (neiType == cellCellStencil::CALCULATED);
-
-
-            if (ownCalc || neiCalc)
             {
-                // Calculate flux w.r.t. calculated cell
+                flux = -flux;
-                scalar flux = phi[facei];
+            }
-                if (neiCalc)
-                {
-                    flux = -flux;
-                }
 
-                if (flux < 0.0)
+            if (flux < 0.0)
-                {
+            {
-                    massIn[zonei] -= flux;
+                massIn -= flux;
-                }
+            }
-                else
+            else
-                {
+            {
-                    adjustableMassOut[zonei] += flux;
+                massOut += flux;
-                }
             }
         }
     }
@@ -103,131 +101,95 @@ bool Foam::oversetAdjustPhi
     // Check all coupled faces on the outside of interpolated cells
     labelList neiCellTypes;
     syncTools::swapBoundaryCellList(mesh, cellTypes, neiCellTypes);
+
+    forAll(bphi, patchi)
     {
-        forAll(bphi, patchi)
+        const fvPatchVectorField& Up = U.boundaryField()[patchi];
+        const fvsPatchScalarField& phip = bphi[patchi];
+        const labelUList& fc = Up.patch().faceCells();
+
+        const label start = Up.patch().start();
+
+        forAll(fc, i)
         {
-            const fvPatchVectorField& Up = U.boundaryField()[patchi];
+            const label facei = start + i;
-            const fvsPatchScalarField& phip = bphi[patchi];
+            const label celli = fc[i];
-            const labelUList& fc = Up.patch().faceCells();
+            const label ownType = cellTypes[celli];
+            const label neiType = neiCellTypes[facei - mesh.nInternalFaces()];
 
-            label start = Up.patch().start();
+            const label zonei = zoneID[celli];
 
-            forAll(fc, i)
+            const bool ownCalc =
+                (ownType == cellCellStencil::CALCULATED)
+             && (neiType == cellCellStencil::INTERPOLATED);
+
+            if (ownCalc && (zonei == zoneId))
             {
-                label facei = start+i;
+                // Calculate flux w.r.t. calculated cell
-                label celli = fc[i];
+                scalar flux = phip[i];
-                label ownType = cellTypes[celli];
-                label neiType = neiCellTypes[facei-mesh.nInternalFaces()];
 
-                bool ownCalc =
+                if (flux < 0.0)
-                    (ownType == cellCellStencil::CALCULATED)
-                 && (neiType == cellCellStencil::INTERPOLATED);
-
-
-                if (ownCalc)
                 {
-                    // Calculate flux w.r.t. calculated cell
+                    massIn -= flux;
-                    scalar flux = phip[i];
+                }
-
+                else
-                    if (flux < 0.0)
+                {
-                    {
+                    massOut += flux;
-                        massIn[zoneID[celli]] -= flux;
-                    }
-                    else
-                    {
-                        adjustableMassOut[zoneID[celli]] += flux;
-                    }
                 }
             }
         }
     }
+    reduce(massIn, sumOp<scalar>());
+    reduce(massOut, sumOp<scalar>());
 
-    // Calculate the total flux in the domain, used for normalisation
+    const scalar massCorr = massIn/(massOut + SMALL);
-    scalar totalFlux = VSMALL + sum(mag(phi)).value();
 
-    forAll(massIn, zonei)
+    if (fv::debug)
     {
-        reduce(massIn[zonei], sumOp<scalar>());
+        Info<< "Zone                    : " << zoneId << nl
-        reduce(adjustableMassOut[zonei], sumOp<scalar>());
+            << "mass outflow            : " << massOut << nl
+            << "mass inflow             : " << massIn << nl
+            << "correction factor       : " << massCorr << endl;
     }
 
 
-    scalarField massCorr(nZones, 1.0);
-
-    forAll(massIn, zonei)
-    {
-        scalar magAdjustableMassOut = mag(adjustableMassOut[zonei]);
-
-        if
-        (
-            magAdjustableMassOut > VSMALL
-         && magAdjustableMassOut/totalFlux > SMALL
-        )
-        {
-            massCorr[zonei] = massIn[zonei]/adjustableMassOut[zonei];
-        }
-        else if (mag(massIn[zonei])/totalFlux > 1e-8)
-        {
-            WarningInFunction
-                << "Continuity error cannot be removed by adjusting the"
-                   " flow at fringe faces.\n    Please check the cell types"
-                << " from the overset analysis."
-                << nl
-                << "Zone                    : " << zonei << nl
-                << "Total flux              : " << totalFlux << nl
-                << "Specified mass inflow   : " << massIn[zonei] << nl
-                << "Adjustable mass outflow : " << adjustableMassOut[zonei]
-                << nl << endl;
-        }
-
-
-        if (fv::debug)
-        {
-            Info<< "Zone                    : " << zonei << nl
-                << "Total flux              : " << totalFlux << nl
-                << "Specified mass inflow   : " << massIn[zonei] << nl
-                << "Adjustable mass outflow : " << adjustableMassOut[zonei]
-                << nl
-                << "Correction factor       : " << massCorr[zonei] << nl
-                << endl;
-        }
-    }
-
-
-
     // Pass2: adjust fluxes
-
     forAll(own, facei)
     {
-        label zonei = zoneID[own[facei]];   // note:own and nei in same zone
+        const label zonei = zoneID[own[facei]];  // own and nei in same zone
 
-        label ownType = cellTypes[own[facei]];
+        const label ownType = cellTypes[own[facei]];
-        label neiType = cellTypes[nei[facei]];
+        const label neiType = cellTypes[nei[facei]];
 
-        bool ownCalc =
+        const bool ownCalc =
             (ownType == cellCellStencil::CALCULATED)
          && (neiType == cellCellStencil::INTERPOLATED);
 
-        bool neiCalc =
+        const bool neiCalc =
             (ownType == cellCellStencil::INTERPOLATED)
          && (neiType == cellCellStencil::CALCULATED);
 
-        if (ownCalc || neiCalc)
+        const bool ownOrCalc = (ownCalc || neiCalc);
+
+        if
+        (
+            (ownOrCalc && (zonei == zoneId)) || (ownOrCalc && (zoneId == -1))
+        )
         {
-            // Calculate flux w.r.t. calculated cell
             scalar flux = phi[facei];
-            if (neiCalc)
+
+            if (ownCalc)
             {
                 flux = -flux;
             }
 
             if (flux < 0.0)
             {
-                phi[facei] /= Foam::sqrt(massCorr[zonei]);
+                phi[facei] /= Foam::sqrt(massCorr);
             }
             else
             {
-                phi[facei] *= Foam::sqrt(massCorr[zonei]);
+                phi[facei] *= Foam::sqrt(massCorr);
             }
         }
     }
@@ -238,28 +200,31 @@ bool Foam::oversetAdjustPhi
         fvsPatchScalarField& phip = bphi[patchi];
         const labelUList& fc = Up.patch().faceCells();
 
-        label start = Up.patch().start();
+        const label start = Up.patch().start();
 
         forAll(fc, i)
         {
-            label facei = start+i;
+            const label facei = start + i;
-            label celli = fc[i];
+            const label celli = fc[i];
-            label zonei = zoneID[celli];   // note:own and nei in same zone
+            const label zonei = zoneID[celli];   // note:own and nei in same zone
-            label ownType = cellTypes[celli];
+            const label ownType = cellTypes[celli];
-            label neiType = neiCellTypes[facei-mesh.nInternalFaces()];
+            const label neiType = neiCellTypes[facei - mesh.nInternalFaces()];
 
-            bool ownCalc =
+            const bool ownCalc =
                 (ownType == cellCellStencil::CALCULATED)
              && (neiType == cellCellStencil::INTERPOLATED);
 
-            bool neiCalc =
+            const bool neiCalc =
                 (ownType == cellCellStencil::INTERPOLATED)
              && (neiType == cellCellStencil::CALCULATED);
 
-            if (ownCalc || neiCalc)
+            const bool ownOrCalc = (ownCalc || neiCalc);
+
+            if (ownOrCalc && (zonei == zoneId))
             {
                 // Calculate flux w.r.t. calculated cell
                 scalar flux = phip[i];
+
                 if (neiCalc)
                 {
                     flux = -flux;
@@ -267,18 +232,64 @@ bool Foam::oversetAdjustPhi
 
                 if (flux < 0.0)
                 {
-                    phip[i] /= Foam::sqrt(massCorr[zonei]);
+                    // phip[i] /= Foam::sqrt(massCorr[zonei]);
                 }
                 else
                 {
-                    phip[i] *= Foam::sqrt(massCorr[zonei]);
+                    phip[i] *= massCorr;
                 }
             }
         }
     }
 
-    return sum(mag(massIn))/totalFlux < SMALL
+
-        && sum(mag(adjustableMassOut))/totalFlux < SMALL;
+    // Check correction
+    if (fv::debug)
+    {
+        scalar massOutCheck = 0;
+        scalar massInCheck = 0;
+
+        forAll(own, facei)
+        {
+            const label zonei = zoneID[own[facei]];
+            const label ownType = cellTypes[own[facei]];
+            const label neiType = cellTypes[nei[facei]];
+
+            const bool ownCalc =
+                (ownType == cellCellStencil::CALCULATED)
+             && (neiType == cellCellStencil::INTERPOLATED);
+
+            const bool neiCalc =
+                (ownType == cellCellStencil::INTERPOLATED)
+             && (neiType == cellCellStencil::CALCULATED);
+
+            const bool ownOrCalc = (ownCalc || neiCalc);
+
+            if
+            (
+                (ownOrCalc && (zonei == zoneId))||(ownOrCalc && (zoneId == -1))
+            )
+            {
+                scalar flux = phi[facei];
+
+                if (ownCalc)
+                {
+                    flux = -flux;
+                }
+
+                if (flux < 0.0)
+                {
+                    massInCheck -= flux;
+                }
+                else
+                {
+                    massOutCheck += flux;
+                }
+            }
+        }
+    }
+
+    return true;
 }
 
 

src/overset/oversetAdjustPhi/oversetAdjustPhi.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017 OpenCFD Ltd.
+    Copyright (C) 2017-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -51,7 +51,8 @@ namespace Foam
 bool oversetAdjustPhi
 (
     surfaceScalarField& phi,
-    const volVectorField& U
+    const volVectorField& U,
+    const label zoneId = -1
 );
 
 } // End namespace Foam

src/overset/oversetCoupledPolyPatch/oversetGAMGInterface/oversetGAMGInterface.C

@@ -0,0 +1,131 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021-2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "oversetGAMGInterface.H"
+#include "addToRunTimeSelectionTable.H"
+
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(oversetGAMGInterface, 0);
+    addToRunTimeSelectionTable
+    (
+        GAMGInterface,
+        oversetGAMGInterface,
+        lduInterface
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::oversetGAMGInterface::oversetGAMGInterface
+(
+    const label index,
+    const lduInterfacePtrsList& coarseInterfaces,
+    const lduInterface& fineInterface,
+    const labelField& localRestrictAddressing,
+    const labelField& neighbourRestrictAddressing,
+    const label fineLevelIndex,
+    const label coarseComm
+)
+:
+    GAMGInterface
+    (
+        index,
+        coarseInterfaces
+    )
+{
+    // Problem:
+    // - we only want the oversetGAMGInterfaceField
+    // - but this needs a oversetGAMGInterface of at least size 1
+    //   (since faceRestrictAddressing cannot have -1 or so in it - it
+    //    can only have >= 0 i.e. we cannot agglomerate into nothing)
+    // - note that the end result is not used - the (current) corresponding
+    //   oversetGAMGInterfaceField has dummy updateInterfaceMatrix.
+
+    // Construct face agglomeration from cell agglomeration. Code same as
+    // in e.g. cyclicAMIGAMGInterface.
+    {
+        // From coarse face to cell
+        DynamicList<label> dynFaceCells(localRestrictAddressing.size());
+
+        // From face to coarse face
+        DynamicList<label> dynFaceRestrictAddressing
+        (
+            localRestrictAddressing.size()
+        );
+
+        Map<label> masterToCoarseFace(localRestrictAddressing.size());
+
+        for (const label curMaster : localRestrictAddressing)
+        {
+            const auto iter = masterToCoarseFace.cfind(curMaster);
+
+            if (iter.found())
+            {
+                // Already have coarse face
+                dynFaceRestrictAddressing.append(iter.val());
+            }
+            else
+            {
+                // New coarse face
+                const label coarseI = dynFaceCells.size();
+
+                dynFaceRestrictAddressing.append(coarseI);
+                dynFaceCells.append(curMaster);
+                masterToCoarseFace.insert(curMaster, coarseI);
+            }
+        }
+
+        faceCells_.transfer(dynFaceCells);
+        faceRestrictAddressing_.transfer(dynFaceRestrictAddressing);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::oversetGAMGInterface::~oversetGAMGInterface()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::labelField> Foam::oversetGAMGInterface::internalFieldTransfer
+(
+    const Pstream::commsTypes commsType,
+    const labelUList& iF
+) const
+{
+    return tmp<labelField>::New(iF);
+}
+
+
+// ************************************************************************* //

src/overset/oversetCoupledPolyPatch/oversetGAMGInterface/oversetGAMGInterface.H

@@ -0,0 +1,122 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::oversetGAMGInterface
+
+Description
+    GAMG agglomerated cyclic AMI interface.
+
+SourceFiles
+    oversetGAMGInterface.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef oversetGAMGInterface_H
+#define oversetGAMGInterface_H
+
+#include "GAMGInterface.H"
+#include "oversetLduInterface.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                  Class oversetGAMGInterface Declaration
+\*---------------------------------------------------------------------------*/
+
+class oversetGAMGInterface
+:
+    public GAMGInterface,
+    virtual public oversetLduInterface
+{
+    // Private Member Functions
+
+        //- No copy construct
+        oversetGAMGInterface(const oversetGAMGInterface&) = delete;
+
+        //- No copy assignment
+        void operator=(const oversetGAMGInterface&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("overset");
+
+
+    // Constructors
+
+        //- Construct from fine level interface,
+        //- local and neighbour restrict addressing
+        oversetGAMGInterface
+        (
+            const label index,
+            const lduInterfacePtrsList& coarseInterfaces,
+            const lduInterface& fineInterface,
+            const labelField& restrictAddressing,
+            const labelField& neighbourRestrictAddressing,
+            const label fineLevelIndex,
+            const label coarseComm
+        );
+
+
+    //- Destructor
+    virtual ~oversetGAMGInterface();
+
+
+    // Member Functions
+
+        // Interface transfer functions
+
+            //- Transfer and return internal field adjacent to the interface
+            virtual tmp<labelField> internalFieldTransfer
+            (
+                const Pstream::commsTypes commsType,
+                const labelUList& iF
+            ) const;
+
+
+        // I/O
+
+            //- Write to stream
+            virtual void write(Ostream&) const
+            {
+                NotImplemented;
+            }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/overset/oversetCoupledPolyPatch/oversetGAMGInterfaceField/oversetGAMGInterfaceField.C

@@ -0,0 +1,97 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "oversetGAMGInterfaceField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "lduMatrix.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(oversetGAMGInterfaceField, 0);
+    addToRunTimeSelectionTable
+    (
+        GAMGInterfaceField,
+        oversetGAMGInterfaceField,
+        lduInterface
+    );
+    addToRunTimeSelectionTable
+    (
+        GAMGInterfaceField,
+        oversetGAMGInterfaceField,
+        lduInterfaceField
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::oversetGAMGInterfaceField::oversetGAMGInterfaceField
+(
+    const GAMGInterface& GAMGCp,
+    const lduInterfaceField& fineInterface
+)
+:
+    GAMGInterfaceField(GAMGCp, fineInterface)
+{}
+
+
+Foam::oversetGAMGInterfaceField::oversetGAMGInterfaceField
+(
+    const GAMGInterface& GAMGCp,
+    const bool doTransform,
+    const int rank
+)
+:
+    GAMGInterfaceField(GAMGCp, doTransform, rank)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::oversetGAMGInterfaceField::~oversetGAMGInterfaceField()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::oversetGAMGInterfaceField::updateInterfaceMatrix
+(
+    solveScalarField& result,
+    const bool add,
+    const lduAddressing& lduAddr,
+    const label patchId,
+    const solveScalarField& psiInternal,
+    const scalarField& coeffs,
+    const direction cmpt,
+    const Pstream::commsTypes
+) const
+{}
+
+
+// ************************************************************************* //

src/overset/oversetCoupledPolyPatch/oversetGAMGInterfaceField/oversetGAMGInterfaceField.H

@@ -0,0 +1,125 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::oversetGAMGInterfaceField
+
+Description
+    GAMG agglomerated
+
+SourceFiles
+    oversetGAMGInterfaceField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef oversetGAMGInterfaceField_H
+#define oversetGAMGInterfaceField_H
+
+#include "GAMGInterfaceField.H"
+#include "oversetGAMGInterface.H"
+#include "oversetLduInterfaceField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                 Class oversetGAMGInterfaceField Declaration
+\*---------------------------------------------------------------------------*/
+
+class oversetGAMGInterfaceField
+:
+    public GAMGInterfaceField,
+    virtual public oversetLduInterfaceField
+{
+    // Private Member Functions
+
+        //- No copy construct
+        oversetGAMGInterfaceField
+        (
+            const oversetGAMGInterfaceField&
+        ) = delete;
+
+        //- No copy assignment
+        void operator=(const oversetGAMGInterfaceField&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("overset");
+
+
+    // Constructors
+
+        //- Construct from GAMG interface and fine level interface field
+        oversetGAMGInterfaceField
+        (
+            const GAMGInterface& GAMGCp,
+            const lduInterfaceField& fineInterfaceField
+        );
+
+        //- Construct from GAMG interface and fine level interface field
+        oversetGAMGInterfaceField
+        (
+            const GAMGInterface& GAMGCp,
+            const bool doTransform,
+            const int rank
+        );
+
+
+    //- Destructor
+    virtual ~oversetGAMGInterfaceField();
+
+
+    // Member Functions
+
+        // Interface matrix update
+
+            //- Update result field based on interface functionality
+            virtual void updateInterfaceMatrix
+            (
+                solveScalarField& result,
+                const bool add,
+                const lduAddressing& lduAddr,
+                const label patchId,
+                const solveScalarField& psiInternal,
+                const scalarField& coeffs,
+                const direction cmpt,
+                const Pstream::commsTypes commsType
+            ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/overset/oversetCoupledPolyPatch/oversetLduInterface/oversetLduInterface.C

@@ -0,0 +1,44 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+
+\*---------------------------------------------------------------------------*/
+
+#include "oversetLduInterface.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(oversetLduInterface, 0);
+}
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::oversetLduInterface::~oversetLduInterface()
+{}
+
+
+// ************************************************************************* //

src/overset/oversetCoupledPolyPatch/oversetLduInterface/oversetLduInterface.H

@@ -0,0 +1,79 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::oversetLduInterface
+
+Description
+    An abstract base class for overset coupled interfaces.
+
+SourceFiles
+    oversetLduInterface.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef oversetLduInterface_H
+#define oversetLduInterface_H
+
+#include "typeInfo.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class oversetLduInterface Declaration
+\*---------------------------------------------------------------------------*/
+
+class oversetLduInterface
+{
+public:
+
+    //- Runtime type information
+    TypeName("oversetLduInterface");
+
+
+    // Constructors
+
+        //- Construct null
+        oversetLduInterface()
+        {}
+
+
+    //- Destructor
+    virtual ~oversetLduInterface();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/overset/oversetCoupledPolyPatch/oversetLduInterfaceField/oversetLduInterfaceField.C

@@ -0,0 +1,45 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+
+\*---------------------------------------------------------------------------*/
+
+#include "oversetLduInterfaceField.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(oversetLduInterfaceField, 0);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::oversetLduInterfaceField::~oversetLduInterfaceField()
+{}
+
+
+// ************************************************************************* //

src/overset/oversetCoupledPolyPatch/oversetLduInterfaceField/oversetLduInterfaceField.H

@@ -0,0 +1,80 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::oversetLduInterfaceField
+
+Description
+    Abstract base class for overset coupled interface fields.
+
+SourceFiles
+    oversetLduInterfaceField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef oversetLduInterfaceField_H
+#define oversetLduInterfaceField_H
+
+#include "primitiveFieldsFwd.H"
+#include "typeInfo.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                  Class oversetLduInterfaceField Declaration
+\*---------------------------------------------------------------------------*/
+
+class oversetLduInterfaceField
+{
+public:
+
+    //- Runtime type information
+    TypeName("oversetLduInterfaceField");
+
+
+    // Constructors
+
+        //- Construct given coupled patch
+        oversetLduInterfaceField()
+        {}
+
+
+    //- Destructor
+    virtual ~oversetLduInterfaceField();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/overset/oversetFvMesh/dynamicOversetFvMesh/dynamicOversetFvMesh.C

@@ -0,0 +1,86 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "dynamicOversetFvMesh.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(dynamicOversetFvMesh, 0);
+    addToRunTimeSelectionTable(dynamicFvMesh, dynamicOversetFvMesh, IOobject);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::dynamicOversetFvMesh::dynamicOversetFvMesh
+(
+    const IOobject& io,
+    const bool doInit
+)
+:
+    dynamicMotionSolverListFvMesh(io, doInit),
+    oversetFvMeshBase(static_cast<const fvMesh&>(*this))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::dynamicOversetFvMesh::~dynamicOversetFvMesh()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::dynamicOversetFvMesh::update()
+{
+    if (!dynamicMotionSolverListFvMesh::update())
+    {
+        return false;
+    }
+
+    oversetFvMeshBase::update();
+
+    return true;
+}
+
+
+bool Foam::dynamicOversetFvMesh::writeObject
+(
+    IOstreamOption streamOpt,
+    const bool valid
+) const
+{
+    bool ok = dynamicMotionSolverListFvMesh::writeObject(streamOpt, valid);
+    ok = oversetFvMeshBase::writeObject(streamOpt, valid) && ok;
+    return ok;
+}
+
+
+// ************************************************************************* //

src/overset/oversetFvMesh/dynamicOversetFvMesh/dynamicOversetFvMesh.H

@@ -0,0 +1,226 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::dynamicOversetFvMesh
+
+Description
+    dynamicFvMesh with support for overset meshes.
+
+SourceFiles
+    dynamicOversetFvMesh.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef dynamicOversetFvMesh_H
+#define dynamicOversetFvMesh_H
+
+#include "dynamicMotionSolverListFvMesh.H"
+#include "oversetFvMeshBase.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                   Class dynamicOversetFvMesh Declaration
+\*---------------------------------------------------------------------------*/
+
+class dynamicOversetFvMesh
+:
+    public dynamicMotionSolverListFvMesh,
+    public oversetFvMeshBase
+{
+    // Private Member Functions
+
+        //- No copy construct
+        dynamicOversetFvMesh(const dynamicOversetFvMesh&) = delete;
+
+        //- No copy assignment
+        void operator=(const dynamicOversetFvMesh&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("dynamicOversetFvMesh");
+
+
+    // Constructors
+
+        //- Construct from IOobject
+        dynamicOversetFvMesh(const IOobject& io,  const bool doInit=true);
+
+
+    //- Destructor
+    virtual ~dynamicOversetFvMesh();
+
+
+    // Member Functions
+
+        //- Override ldu addressing
+        virtual const lduAddressing& lduAddr() const
+        {
+            return oversetFvMeshBase::lduAddr();
+        }
+
+        //- Override ldu addressing
+        virtual lduInterfacePtrsList interfaces() const
+        {
+            return oversetFvMeshBase::interfaces();
+        }
+
+
+        // Overset
+
+            // Explicit interpolation
+
+                virtual void interpolate(scalarField& psi) const
+                {
+                    Stencil::New(*this).interpolate<scalar>(*this, psi);
+                }
+
+                virtual void interpolate(vectorField& psi) const
+                {
+                    Stencil::New(*this).interpolate<vector>(*this, psi);
+                }
+
+                virtual void interpolate(sphericalTensorField& psi) const
+                {
+                    Stencil::New
+                    (
+                        *this
+                    ).interpolate<sphericalTensor>(*this, psi);
+                }
+
+                virtual void interpolate(symmTensorField& psi) const
+                {
+                    Stencil::New(*this).interpolate<symmTensor>(*this, psi);
+                }
+
+                virtual void interpolate(tensorField& psi) const
+                {
+                    Stencil::New(*this).interpolate<tensor>(*this, psi);
+                }
+
+                virtual void interpolate(volScalarField& psi) const
+                {
+                    Stencil::New(*this).interpolate<volScalarField>(psi);
+                }
+
+                virtual void interpolate(volVectorField& psi) const
+                {
+                    Stencil::New(*this).interpolate<volVectorField>(psi);
+                }
+
+                virtual void interpolate(volSphericalTensorField& psi) const
+                {
+                    Stencil::New(*this).interpolate
+                    <
+                        volSphericalTensorField
+                    >(psi);
+                }
+
+                virtual void interpolate(volSymmTensorField& psi) const
+                {
+                    Stencil::New(*this).interpolate<volSymmTensorField>(psi);
+                }
+
+                virtual void interpolate(volTensorField& psi) const
+                {
+                    Stencil::New(*this).interpolate<volTensorField>(psi);
+                }
+
+
+            // Implicit interpolation (matrix manipulation)
+
+                //- Solve returning the solution statistics given convergence
+                //- tolerance. Use the given solver controls
+                virtual SolverPerformance<scalar> solve
+                (
+                    fvMatrix<scalar>& m,
+                    const dictionary& dict
+                ) const
+                {
+                    return oversetFvMeshBase::solveOverset<scalar>(m, dict);
+                }
+
+                //- Solve returning the solution statistics given convergence
+                //- tolerance. Use the given solver controls
+                virtual SolverPerformance<vector> solve
+                (
+                    fvMatrix<vector>& m,
+                    const dictionary& dict
+                ) const
+                {
+                    return oversetFvMeshBase::solveOverset<vector>(m, dict);
+                }
+
+                //- Solve returning the solution statistics given convergence
+                //- tolerance. Use the given solver controls
+                virtual SolverPerformance<symmTensor> solve
+                (
+                    fvMatrix<symmTensor>& m,
+                    const dictionary& dict
+                ) const
+                {
+                    return oversetFvMeshBase::solveOverset<symmTensor>
+                    (m, dict);
+                }
+
+                //- Solve returning the solution statistics given convergence
+                //- tolerance. Use the given solver controls
+                virtual SolverPerformance<tensor> solve
+                (
+                    fvMatrix<tensor>& m,
+                    const dictionary& dict
+                ) const
+                {
+                    return oversetFvMeshBase::solveOverset<tensor>(m, dict);
+                }
+
+
+        //- Update the mesh for both mesh motion and topology change
+        virtual bool update();
+
+        //- Write using given format, version and compression
+        virtual bool writeObject
+        (
+            IOstreamOption streamOpt,
+            const bool valid
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/overset/oversetFvMesh/oversetFvMeshBase.C

@@ -25,28 +25,27 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "dynamicOversetFvMesh.H"
+#include "oversetFvMeshBase.H"
 #include "addToRunTimeSelectionTable.H"
 #include "cellCellStencilObject.H"
 #include "zeroGradientFvPatchFields.H"
 #include "lduPrimitiveProcessorInterface.H"
 #include "globalIndex.H"
+#include "GAMGAgglomeration.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    defineTypeNameAndDebug(dynamicOversetFvMesh, 0);
+    defineTypeNameAndDebug(oversetFvMeshBase, 0);
-    addToRunTimeSelectionTable(dynamicFvMesh, dynamicOversetFvMesh, IOobject);
-    addToRunTimeSelectionTable(dynamicFvMesh, dynamicOversetFvMesh, doInit);
 }
 
 
 // * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * * //
 
-bool Foam::dynamicOversetFvMesh::updateAddressing() const
+bool Foam::oversetFvMeshBase::updateAddressing() const
 {
-    const cellCellStencilObject& overlap = Stencil::New(*this);
+    const cellCellStencilObject& overlap = Stencil::New(mesh_);
 
     // The (processor-local part of the) stencil determines the local
     // faces to add to the matrix. tbd: parallel
@@ -54,7 +53,7 @@ bool Foam::dynamicOversetFvMesh::updateAddressing() const
 
     // Get the base addressing
     //const lduAddressing& baseAddr = dynamicMotionSolverFvMesh::lduAddr();
-    const lduAddressing& baseAddr = dynamicFvMesh::lduAddr();
+    const lduAddressing& baseAddr = mesh_.fvMesh::lduAddr();
 
     // Add to the base addressing
     labelList lowerAddr;
@@ -89,13 +88,12 @@ bool Foam::dynamicOversetFvMesh::updateAddressing() const
 
     if (debug)
     {
-        Pout<< "dynamicOversetFvMesh::update() : extended addressing from"
+        Pout<< "oversetFvMeshBase::update() : extended addressing from"
             << " nFaces:" << baseAddr.lowerAddr().size()
             << " to nFaces:" << lowerAddr.size()
             << " nExtraFaces:" << nExtraFaces << endl;
     }
 
-
     // Now for the tricky bits. We want to hand out processor faces according
     // to the localFaceCells/remoteFaceCells. Ultimately we need
     // per entry in stencil:
@@ -162,7 +160,6 @@ bool Foam::dynamicOversetFvMesh::updateAddressing() const
     //      proc1. We first want the pair : proc1 receive and proc2 send
     //      and then the pair : proc1 send, proc2 receive
 
-
     labelList procToInterface(Pstream::nProcs(), -1);
 
     forAll(procOwner, proci)
@@ -266,32 +263,31 @@ bool Foam::dynamicOversetFvMesh::updateAddressing() const
         {
             if (interface != -1)
             {
-                interface = procToInterface[interface]+boundary().size();
+                interface = procToInterface[interface]+mesh_.boundary().size();
             }
         }
     }
 
-
     // Get addressing and interfaces of all interfaces
 
-
     UPtrList<const labelUList> patchAddr;
     {
-        const fvBoundaryMesh& fvp = boundary();
+        const fvBoundaryMesh& fvp = mesh_.boundary();
 
         patchAddr.setSize(fvp.size() + remoteStencilInterfaces_.size());
 
         //allInterfaces_ = dynamicMotionSolverFvMesh::interfaces();
-        allInterfaces_ = dynamicFvMesh::interfaces();
+        allInterfaces_ = mesh_.fvMesh::interfaces();
         allInterfaces_.setSize(patchAddr.size());
 
         forAll(fvp, patchi)
         {
             patchAddr.set(patchi, &fvp[patchi].faceCells());
         }
+
         forAll(remoteStencilInterfaces_, i)
         {
-            const label patchi = fvp.size()+i;
+            label patchi = fvp.size()+i;
             const lduPrimitiveProcessorInterface& pp =
                 remoteStencilInterfaces_[i];
 
@@ -305,6 +301,7 @@ bool Foam::dynamicOversetFvMesh::updateAddressing() const
             allInterfaces_.set(patchi, &pp);
         }
     }
+
     const lduSchedule ps
     (
         lduPrimitiveMesh::nonBlockingSchedule<processorLduInterface>
@@ -317,7 +314,7 @@ bool Foam::dynamicOversetFvMesh::updateAddressing() const
     (
         new fvMeshPrimitiveLduAddressing
         (
-            nCells(),
+            mesh_.nCells(),
             std::move(lowerAddr),
             std::move(upperAddr),
             patchAddr,
@@ -336,22 +333,22 @@ bool Foam::dynamicOversetFvMesh::updateAddressing() const
             << " upper:" << addr.upperAddr().size() << endl;
 
         // Using lduAddressing::patch
-        forAll(patchAddr, patchI)
+        forAll(patchAddr, patchi)
         {
-            Pout<< "    " << patchI << "\tpatchAddr:"
+            Pout<< "    " << patchi << "\tpatchAddr:"
-                << addr.patchAddr(patchI).size()
+                << addr.patchAddr(patchi).size()
                 << endl;
         }
 
         // Using interfaces
-        const lduInterfacePtrsList& iFaces = allInterfaces_;
+        const lduInterfacePtrsList& iFaces = mesh_.interfaces();
         Pout<< "Adapted interFaces:" << iFaces.size() << endl;
-        forAll(iFaces, patchI)
+        forAll(iFaces, patchi)
         {
-            if (iFaces.set(patchI))
+            if (iFaces.set(patchi))
             {
-                Pout<< "    " << patchI << "\tinterface:"
+                Pout<< "    " << patchi << "\tinterface:"
-                    << iFaces[patchI].type() << endl;
+                    << iFaces[patchi].type() << endl;
             }
         }
     }
@@ -360,7 +357,7 @@ bool Foam::dynamicOversetFvMesh::updateAddressing() const
 }
 
 
-Foam::scalar Foam::dynamicOversetFvMesh::cellAverage
+Foam::scalar Foam::oversetFvMeshBase::cellAverage
 (
     const labelList& types,
     const labelList& nbrTypes,
@@ -370,16 +367,16 @@ Foam::scalar Foam::dynamicOversetFvMesh::cellAverage
     bitSet& isFront
 ) const
 {
-    const labelList& own = faceOwner();
+    const labelList& own = mesh_.faceOwner();
-    const labelList& nei = faceNeighbour();
+    const labelList& nei = mesh_.faceNeighbour();
-    const cell& cFaces = cells()[celli];
+    const cell& cFaces = mesh_.cells()[celli];
 
     scalar avg = 0.0;
     label n = 0;
     label nFront = 0;
     for (const label facei : cFaces)
     {
-        if (isInternalFace(facei))
+        if (mesh_.isInternalFace(facei))
         {
             label nbrCelli = (own[facei] == celli ? nei[facei] : own[facei]);
             if (norm[nbrCelli] == -GREAT)
@@ -399,7 +396,7 @@ Foam::scalar Foam::dynamicOversetFvMesh::cellAverage
         }
         else
         {
-            if (nbrNorm[facei-nInternalFaces()] == -GREAT)
+            if (nbrNorm[facei-mesh_.nInternalFaces()] == -GREAT)
             {
                 if (isFront.set(facei))
                 {
@@ -408,7 +405,7 @@ Foam::scalar Foam::dynamicOversetFvMesh::cellAverage
             }
             else
             {
-                avg += nbrNorm[facei-nInternalFaces()];
+                avg += nbrNorm[facei-mesh_.nInternalFaces()];
                 n++;
             }
         }
@@ -425,13 +422,13 @@ Foam::scalar Foam::dynamicOversetFvMesh::cellAverage
 }
 
 
-void Foam::dynamicOversetFvMesh::writeAgglomeration
+void Foam::oversetFvMeshBase::writeAgglomeration
 (
     const GAMGAgglomeration& agglom
 ) const
 {
-    labelList cellToCoarse(identity(nCells()));
+    labelList cellToCoarse(identity(mesh_.nCells()));
-    labelListList coarseToCell(invertOneToMany(nCells(), cellToCoarse));
+    labelListList coarseToCell(invertOneToMany(mesh_.nCells(), cellToCoarse));
 
     // Write initial agglomeration
     {
@@ -440,13 +437,13 @@ void Foam::dynamicOversetFvMesh::writeAgglomeration
             IOobject
             (
                 "agglomeration",
-                this->time().timeName(),
+                mesh_.time().timeName(),
-                *this,
+                mesh_,
                 IOobject::NO_READ,
                 IOobject::NO_WRITE,
                 false
             ),
-            *this,
+            mesh_,
             dimensionedScalar(dimless, Zero)
         );
         scalarField& fld = scalarAgglomeration.primitiveFieldRef();
@@ -463,7 +460,7 @@ void Foam::dynamicOversetFvMesh::writeAgglomeration
         scalarAgglomeration.write();
 
         Info<< "Writing initial cell distribution to "
-            << this->time().timeName() << endl;
+            << mesh_.time().timeName() << endl;
     }
 
 
@@ -495,13 +492,13 @@ void Foam::dynamicOversetFvMesh::writeAgglomeration
                 IOobject
                 (
                     "agglomeration_" + Foam::name(level),
-                    this->time().timeName(),
+                    mesh_.time().timeName(),
-                    *this,
+                    mesh_,
                     IOobject::NO_READ,
                     IOobject::NO_WRITE,
                     false
                 ),
-                *this,
+                mesh_,
                 dimensionedScalar(dimless, Zero)
             );
             scalarField& fld = scalarAgglomeration.primitiveFieldRef();
@@ -526,52 +523,30 @@ void Foam::dynamicOversetFvMesh::writeAgglomeration
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::dynamicOversetFvMesh::dynamicOversetFvMesh
+Foam::oversetFvMeshBase::oversetFvMeshBase(const fvMesh& mesh, bool doInit)
-(
-    const IOobject& io,
-    const bool doInit
-)
 :
-    dynamicMotionSolverListFvMesh(io, doInit)
+    mesh_(mesh),
+    active_(false)
 {
-    if (doInit)
-    {
-        init(false);    // do not initialise lower levels
-    }
-}
-
-
-bool Foam::dynamicOversetFvMesh::init(const bool doInit)
-{
-    if (doInit)
-    {
-        dynamicMotionSolverListFvMesh::init(doInit);
-    }
-
-    active_ = false;
-
     // Load stencil (but do not update)
-    (void)Stencil::New(*this, false);
+    (void)Stencil::New(mesh_, false);
-
-    // Assume something changed
-    return true;
 }
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-Foam::dynamicOversetFvMesh::~dynamicOversetFvMesh()
+Foam::oversetFvMeshBase::~oversetFvMeshBase()
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-const Foam::lduAddressing& Foam::dynamicOversetFvMesh::lduAddr() const
+const Foam::lduAddressing& Foam::oversetFvMeshBase::lduAddr() const
 {
     if (!active_)
     {
         //return dynamicMotionSolverFvMesh::lduAddr();
-        return dynamicFvMesh::lduAddr();
+        return mesh_.fvMesh::lduAddr();
     }
     if (!lduPtr_)
     {
@@ -582,12 +557,11 @@ const Foam::lduAddressing& Foam::dynamicOversetFvMesh::lduAddr() const
 }
 
 
-Foam::lduInterfacePtrsList Foam::dynamicOversetFvMesh::interfaces() const
+Foam::lduInterfacePtrsList Foam::oversetFvMeshBase::interfaces() const
 {
     if (!active_)
     {
-        //return dynamicMotionSolverFvMesh::interfaces();
+        return mesh_.fvMesh::interfaces();
-        return dynamicFvMesh::interfaces();
     }
     if (!lduPtr_)
     {
@@ -599,7 +573,7 @@ Foam::lduInterfacePtrsList Foam::dynamicOversetFvMesh::interfaces() const
 
 
 const Foam::fvMeshPrimitiveLduAddressing&
-Foam::dynamicOversetFvMesh::primitiveLduAddr() const
+Foam::oversetFvMeshBase::primitiveLduAddr() const
 {
     if (!lduPtr_)
     {
@@ -611,10 +585,8 @@ Foam::dynamicOversetFvMesh::primitiveLduAddr() const
 }
 
 
-bool Foam::dynamicOversetFvMesh::update()
+bool Foam::oversetFvMeshBase::update()
 {
-    //if (dynamicMotionSolverFvMesh::update())
-    if (dynamicMotionSolverListFvMesh::update())
     {
         // Calculate the local extra faces for the interpolation. Note: could
         // let demand-driven lduAddr() trigger it but just to make sure.
@@ -631,55 +603,43 @@ bool Foam::dynamicOversetFvMesh::update()
 }
 
 
-Foam::word Foam::dynamicOversetFvMesh::baseName(const word& name)
+bool Foam::oversetFvMeshBase::interpolateFields()
 {
-    if (name.ends_with("_0"))
+    const cellCellStencilObject& overlap = Stencil::New(mesh_);
-    {
-        return baseName(name.substr(0, name.size()-2));
-    }
 
-    return name;
-}
-
-
-bool Foam::dynamicOversetFvMesh::interpolateFields()
-{
     // Add the stencil suppression list
-    wordHashSet suppressed(Stencil::New(*this).nonInterpolatedFields());
+    wordHashSet suppressed(overlap.nonInterpolatedFields());
 
     // Use whatever the solver has set up as suppression list
     const dictionary* dictPtr
     (
-        this->schemesDict().findDict("oversetInterpolationSuppressed")
+        mesh_.schemesDict().findDict("oversetInterpolationSuppressed")
     );
     if (dictPtr)
     {
         suppressed.insert(dictPtr->toc());
     }
 
-    interpolate<volScalarField>(suppressed);
+    overlap.interpolate<volScalarField>(mesh_, suppressed);
-    interpolate<volVectorField>(suppressed);
+    overlap.interpolate<volVectorField>(mesh_, suppressed);
-    interpolate<volSphericalTensorField>(suppressed);
+    overlap.interpolate<volSphericalTensorField>(mesh_, suppressed);
-    interpolate<volSymmTensorField>(suppressed);
+    overlap.interpolate<volSymmTensorField>(mesh_, suppressed);
-    interpolate<volTensorField>(suppressed);
+    overlap.interpolate<volTensorField>(mesh_, suppressed);
 
     return true;
 }
 
 
-
+bool Foam::oversetFvMeshBase::writeObject
-bool Foam::dynamicOversetFvMesh::writeObject
 (
     IOstreamOption streamOpt,
     const bool valid
 ) const
 {
-    //bool ok = dynamicMotionSolverFvMesh::writeObject(streamOpt, valid);
-    bool ok = dynamicMotionSolverListFvMesh::writeObject(streamOpt, valid);
-
     // For postprocessing : write cellTypes and zoneID
+    bool ok = false;
     {
-        const cellCellStencilObject& overlap = Stencil::New(*this);
+        const cellCellStencilObject& overlap = Stencil::New(mesh_);
 
         const labelUList& cellTypes = overlap.cellTypes();
 
@@ -688,13 +648,13 @@ bool Foam::dynamicOversetFvMesh::writeObject
             IOobject
             (
                 "cellTypes",
-                this->time().timeName(),
+                mesh_.time().timeName(),
-                *this,
+                mesh_,
                 IOobject::NO_READ,
                 IOobject::NO_WRITE,
                 false
             ),
-            *this,
+            mesh_,
             dimensionedScalar(dimless, Zero),
             zeroGradientFvPatchScalarField::typeName
         );
@@ -712,18 +672,18 @@ bool Foam::dynamicOversetFvMesh::writeObject
             IOobject
             (
                 "zoneID",
-                this->time().timeName(),
+                mesh_.time().timeName(),
-                *this,
+                mesh_,
                 IOobject::NO_READ,
                 IOobject::NO_WRITE,
                 false
             ),
-            *this,
+            mesh_,
             dimensionedScalar(dimless, Zero),
             zeroGradientFvPatchScalarField::typeName
         );
 
-        const cellCellStencilObject& overlap = Stencil::New(*this);
+        const cellCellStencilObject& overlap = Stencil::New(mesh_);
         const labelIOList& zoneID = overlap.zoneID();
 
         forAll(zoneID, cellI)
@@ -735,7 +695,7 @@ bool Foam::dynamicOversetFvMesh::writeObject
     }
     if (debug)
     {
-        const cellCellStencilObject& overlap = Stencil::New(*this);
+        const cellCellStencilObject& overlap = Stencil::New(mesh_);
         const labelIOList& zoneID = overlap.zoneID();
         const labelListList& cellStencil = overlap.cellStencil();
 
@@ -744,7 +704,7 @@ bool Foam::dynamicOversetFvMesh::writeObject
         overlap.cellInterpolationMap().distribute(donorZoneID);
 
         // Get remote cellCentres
-        pointField cc(C());
+        pointField cc(mesh_.C());
         overlap.cellInterpolationMap().distribute(cc);
 
         volScalarField volDonorZoneID
@@ -752,13 +712,13 @@ bool Foam::dynamicOversetFvMesh::writeObject
             IOobject
             (
                 "donorZoneID",
-                this->time().timeName(),
+                mesh_.time().timeName(),
-                *this,
+                mesh_,
                 IOobject::NO_READ,
                 IOobject::NO_WRITE,
                 false
             ),
-            *this,
+            mesh_,
             dimensionedScalar("minOne", dimless, scalar(-1)),
             zeroGradientFvPatchScalarField::typeName
         );
@@ -774,7 +734,7 @@ bool Foam::dynamicOversetFvMesh::writeObject
                     if (donorZoneID[stencil[i]] != volDonorZoneID[cellI])
                     {
                         WarningInFunction << "Mixed donor meshes for cell "
-                            << cellI << " at " << C()[cellI]
+                            << cellI << " at " << mesh_.C()[cellI]
                             << " donors:" << UIndirectList<point>(cc, stencil)
                             << endl;
                         volDonorZoneID[cellI] = -2;
@@ -784,7 +744,15 @@ bool Foam::dynamicOversetFvMesh::writeObject
         }
         //- Do not correctBoundaryConditions since re-interpolates!
         //volDonorZoneID.correctBoundaryConditions();
-        volDonorZoneID.writeObject(streamOpt, valid);
+        cellCellStencil::correctBoundaryConditions
+        <
+            volScalarField,
+            oversetFvPatchField<scalar>
+        >
+        (
+            volDonorZoneID
+        );
+        ok = volDonorZoneID.writeObject(streamOpt, valid);
     }
 
     return ok;

src/overset/oversetFvMesh/oversetFvMeshBase.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2015-2020 OpenCFD Ltd.
+    Copyright (C) 2015-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -24,45 +24,53 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::dynamicOversetFvMesh
+    Foam::oversetFvMeshBase
 
 Description
-    dynamicFvMesh with support for overset meshes.
+    Support for overset functionality.
 
 SourceFiles
-    dynamicOversetFvMesh.C
+    oversetFvMeshBase.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef dynamicOversetFvMesh_H
+#ifndef oversetFvMeshBase_H
-#define dynamicOversetFvMesh_H
+#define oversetFvMeshBase_H
 
-#include "dynamicMotionSolverListFvMesh.H"
+#include "fvMeshPrimitiveLduAddressing.H"
-#include "labelIOList.H"
 #include "fvMeshPrimitiveLduAddressing.H"
 #include "lduInterfaceFieldPtrsList.H"
-#include "volFields.H"
+#include "volFieldsFwd.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
-class mapDistribute;
 class lduPrimitiveProcessorInterface;
 class GAMGAgglomeration;
 
 /*---------------------------------------------------------------------------*\
-                   Class dynamicOversetFvMesh Declaration
+                   Class oversetFvMeshBase Declaration
 \*---------------------------------------------------------------------------*/
 
-class dynamicOversetFvMesh
+class oversetFvMeshBase
-:
-    public dynamicMotionSolverListFvMesh
 {
+    // Private Member Functions
+
+        //- No copy construct
+        oversetFvMeshBase(const oversetFvMeshBase&) = delete;
+
+        //- No copy assignment
+        void operator=(const oversetFvMeshBase&) = delete;
+
+
 protected:
 
-    // Protected data
+    // Protected Data
+
+        //- Reference to mesh
+        const fvMesh& mesh_;
 
         //- Select base addressing (false) or locally stored extended
         //  lduAddressing (true)
@@ -72,12 +80,11 @@ protected:
         mutable autoPtr<fvMeshPrimitiveLduAddressing> lduPtr_;
 
         //- Added (processor)lduInterfaces for remote bits of stencil.
-        //PtrList<const lduInterface> remoteStencilInterfaces_;
         mutable PtrList<const lduPrimitiveProcessorInterface>
             remoteStencilInterfaces_;
 
         //- Interfaces for above mesh. Contains both original and
-        //  above added processorLduInterfaces
+        //- above added processorLduInterfaces
         mutable lduInterfacePtrsList allInterfaces_;
 
         //- Corresponding faces (in above lduPtr) to the stencil
@@ -105,44 +112,29 @@ protected:
             const lduInterfacePtrsList&
         ) const;
 
-        //- Explicit interpolation of acceptor cells from donor cells
-        template<class T>
-        void interpolate(Field<T>& psi) const;
-
-        //- Explicit interpolation of acceptor cells from donor cells with
-        //  boundary condition handling
-        template<class GeoField>
-        void interpolate(GeoField& psi) const;
-
-        //- Helper: strip off trailing _0
-        static word baseName(const word& name);
-
-        //- Explicit interpolation of all registered fields. Excludes
-        //  selected fields (and their old-time fields)
-        template<class GeoField>
-        void interpolate(const wordHashSet& suppressed);
-
         //- Freeze values at holes
         //template<class Type>
         //void freezeHoles(fvMatrix<Type>&) const;
 
-        //- Get scalar interfaces of certain type
+         //- Scale coefficient depending on cell type
-        //template<class GeoField, class PatchType>
-        //lduInterfaceFieldPtrsList scalarInterfaces(const GeoField& psi) const;
-
-        //- Determine normalisation for interpolation. This equals the
-        //  original diagonal except stabilised for zero diagonals (possible
-        //  in hole cells)
         template<class Type>
-        tmp<scalarField> normalisation(const fvMatrix<Type>& m) const;
+        void scaleConnection
-
+        (
-        //- Add interpolation to matrix (coefficients)
+            Field<Type>& coeffs,
-        template<class Type>
+            const labelUList& types,
-        void addInterpolation(fvMatrix<Type>&, const scalarField& norm) const;
+            const scalarList& factor,
+            const bool setHoleCellValue,
+            const label celli,
+            const label facei
+        ) const;
 
         //- Solve given dictionary with settings
         template<class Type>
-        SolverPerformance<Type> solve(fvMatrix<Type>&, const dictionary&) const;
+        SolverPerformance<Type> solveOverset
+        (
+            fvMatrix<Type>& m,
+            const dictionary&
+        ) const;
 
         //- Debug: correct coupled bc
         template<class GeoField>
@@ -166,36 +158,24 @@ protected:
         ) const;
 
 
-private:
-
-    // Private Member Functions
-
-        //- No copy construct
-        dynamicOversetFvMesh(const dynamicOversetFvMesh&) = delete;
-
-        //- No copy assignment
-        void operator=(const dynamicOversetFvMesh&) = delete;
-
-
 public:
 
     //- Runtime type information
-    TypeName("dynamicOversetFvMesh");
+    TypeName("oversetFvMeshBase");
 
 
     // Constructors
 
         //- Construct from IOobject
-        dynamicOversetFvMesh(const IOobject& io, const bool doInit=true);
+        oversetFvMeshBase(const fvMesh& mesh, const bool doInit=true);
 
 
     //- Destructor
-    virtual ~dynamicOversetFvMesh();
+    virtual ~oversetFvMeshBase();
 
 
     // Member Functions
 
-
         // Extended addressing
 
             //- Return extended ldu addressing
@@ -231,12 +211,14 @@ public:
                 {
                     DebugInfo<< "Switching to extended addressing with nFaces:"
                         << primitiveLduAddr().lowerAddr().size()
+                        << " nInterfaces:" << allInterfaces_.size()
                         << endl;
                 }
                 else
                 {
                     DebugInfo<< "Switching to base addressing with nFaces:"
-                        << fvMesh::lduAddr().lowerAddr().size()
+                        << mesh_.fvMesh::lduAddr().lowerAddr().size()
+                        << " nInterfaces:" << mesh_.fvMesh::interfaces().size()
                         << endl;
                 }
             }
@@ -244,108 +226,17 @@ public:
 
         // Overset
 
-            // Explicit interpolation
+            //- Manipulate the matrix to add the interpolation/set hole
+            //  values
+            template<class Type>
+            void addInterpolation
+            (
+                fvMatrix<Type>& m,
+                const scalarField& normalisation,
+                const bool setHoleCellValue,
+                const Type& holeCellValue
+            ) const;
 
-                virtual void interpolate(scalarField& psi) const
-                {
-                    interpolate<scalar>(psi);
-                }
-
-                virtual void interpolate(vectorField& psi) const
-                {
-                    interpolate<vector>(psi);
-                }
-
-                virtual void interpolate(sphericalTensorField& psi) const
-                {
-                    interpolate<sphericalTensor>(psi);
-                }
-
-                virtual void interpolate(symmTensorField& psi) const
-                {
-                    interpolate<symmTensor>(psi);
-                }
-
-                virtual void interpolate(tensorField& psi) const
-                {
-                    interpolate<tensor>(psi);
-                }
-
-                virtual void interpolate(volScalarField& psi) const
-                {
-                    interpolate<volScalarField>(psi);
-                }
-
-                virtual void interpolate(volVectorField& psi) const
-                {
-                    interpolate<volVectorField>(psi);
-                }
-
-                virtual void interpolate(volSphericalTensorField& psi) const
-                {
-                    interpolate<volSphericalTensorField>(psi);
-                }
-
-                virtual void interpolate(volSymmTensorField& psi) const
-                {
-                    interpolate<volSymmTensorField>(psi);
-                }
-
-                virtual void interpolate(volTensorField& psi) const
-                {
-                    interpolate<volTensorField>(psi);
-                }
-
-
-            // Implicit interpolation (matrix manipulation)
-
-                //- Solve returning the solution statistics given convergence
-                //  tolerance. Use the given solver controls
-                virtual SolverPerformance<scalar> solve
-                (
-                    fvMatrix<scalar>& m,
-                    const dictionary& dict
-                ) const
-                {
-                    return solve<scalar>(m, dict);
-                }
-
-                //- Solve returning the solution statistics given convergence
-                //  tolerance. Use the given solver controls
-                virtual SolverPerformance<vector> solve
-                (
-                    fvMatrix<vector>& m,
-                    const dictionary& dict
-                ) const
-                {
-                    return solve<vector>(m, dict);
-                }
-
-                //- Solve returning the solution statistics given convergence
-                //  tolerance. Use the given solver controls
-                virtual SolverPerformance<symmTensor> solve
-                (
-                    fvMatrix<symmTensor>& m,
-                    const dictionary& dict
-                ) const
-                {
-                    return solve<symmTensor>(m, dict);
-                }
-
-                //- Solve returning the solution statistics given convergence
-                //  tolerance. Use the given solver controls
-                virtual SolverPerformance<tensor> solve
-                (
-                    fvMatrix<tensor>& m,
-                    const dictionary& dict
-                ) const
-                {
-                    return solve<tensor>(m, dict);
-                }
-
-
-        //- Initialise all non-demand-driven data
-        virtual bool init(const bool doInit);
 
         //- Update the mesh for both mesh motion and topology change
         virtual bool update();
@@ -372,6 +263,12 @@ public:
             typename GeoField::Boundary& bfld,
             const bool typeOnly
         );
+
+        //- Determine normalisation for interpolation. This equals the
+        //- original diagonal except stabilised for zero diagonals (possible
+        //- in hole cells)
+        template<class Type>
+        tmp<scalarField> normalisation(const fvMatrix<Type>& m) const;
 };
 
 
@@ -382,7 +279,7 @@ public:
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #ifdef NoRepository
-#   include "dynamicOversetFvMeshTemplates.C"
+#   include "oversetFvMeshBaseTemplates.C"
 #endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/overset/oversetFvMesh/oversetFvMeshBaseTemplates.C

@@ -27,7 +27,6 @@ License
 
 #include "volFields.H"
 #include "fvMatrix.H"
-#include "cellCellStencilObject.H"
 #include "oversetFvPatchField.H"
 #include "calculatedProcessorFvPatchField.H"
 #include "lduInterfaceFieldPtrsList.H"
@@ -36,84 +35,46 @@ License
 
 // * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * * //
 
-template<class T>
+template<class Type>
-void Foam::dynamicOversetFvMesh::interpolate(Field<T>& psi) const
+void Foam::oversetFvMeshBase::scaleConnection
+(
+    Field<Type>& coeffs,
+    const labelUList& types,
+    const scalarList& factor,
+    const bool setHoleCellValue,
+    const label celli,
+    const label facei
+) const
 {
-    const cellCellStencil& overlap = Stencil::New(*this);
+    const label cType = types[celli];
-    const labelListList& stencil = overlap.cellStencil();
+    const scalar f = factor[celli];
 
-    if (stencil.size() != nCells())
+    if (cType == cellCellStencil::INTERPOLATED)
     {
-        return;
+        coeffs[facei] *= 1.0-f;
     }
-
+    else if (cType == cellCellStencil::HOLE)
-    const mapDistribute& map = overlap.cellInterpolationMap();
-    const List<scalarList>& wghts = overlap.cellInterpolationWeights();
-    const labelList& cellIDs = overlap.interpolationCells();
-    const scalarList& factor = overlap.cellInterpolationWeight();
-
-    Field<T> work(psi);
-    map.mapDistributeBase::distribute(work, UPstream::msgType()+1);
-
-    forAll(cellIDs, i)
     {
-        label celli = cellIDs[i];
+        // Disconnect hole cell from influence of neighbour
-
+        coeffs[facei] = pTraits<Type>::zero;
-        const scalarList& w = wghts[celli];
-        const labelList& nbrs = stencil[celli];
-        const scalar f = factor[celli];
-
-        T s(pTraits<T>::zero);
-        forAll(nbrs, nbrI)
-        {
-            s += w[nbrI]*work[nbrs[nbrI]];
-        }
-        //Pout<< "Interpolated value:" << s << endl;
-        //T oldPsi = psi[celli];
-        psi[celli] = (1.0-f)*psi[celli] + f*s;
-        //Pout<< "psi was:" << oldPsi << " now:" << psi[celli] << endl;
     }
-}
+    else if (cType == cellCellStencil::SPECIAL)
-
-
-template<class GeoField>
-void Foam::dynamicOversetFvMesh::interpolate(GeoField& psi) const
-{
-    interpolate(psi.primitiveFieldRef());
-    psi.correctBoundaryConditions();
-}
-
-
-template<class GeoField>
-void Foam::dynamicOversetFvMesh::interpolate(const wordHashSet& suppressed)
-{
-    auto flds(this->lookupClass<GeoField>());
-    for (auto fldPtr : flds)
     {
-        const word& name = fldPtr->name();
+        if (setHoleCellValue)
-        if (!suppressed.found(baseName(name)))
         {
-            if (debug)
+            // Behave like hole
-            {
+            coeffs[facei] = pTraits<Type>::zero;
-                Pout<< "dynamicOversetFvMesh::interpolate: interpolating : "
-                    << name << endl;
-            }
-            interpolate(fldPtr->primitiveFieldRef());
         }
         else
         {
-            if (debug)
+            // Behave like interpolated
-            {
+            coeffs[facei] *= 1.0-f;
-                Pout<< "dynamicOversetFvMesh::interpolate: skipping : " << name
-                    << endl;
-            }
         }
     }
 }
 
-
 template<class GeoField, class PatchType>
-void Foam::dynamicOversetFvMesh::correctBoundaryConditions
+void Foam::oversetFvMeshBase::correctBoundaryConditions
 (
     typename GeoField::Boundary& bfld,
     const bool typeOnly
@@ -133,7 +94,7 @@ void Foam::dynamicOversetFvMesh::correctBoundaryConditions
     // Wait for outstanding requests
     if (commsType == UPstream::commsTypes::nonBlocking)
     {
-        UPstream::waitRequests(startOfRequests);
+        Pstream::waitRequests(startOfRequests);
     }
 
     forAll(bfld, patchi)
@@ -147,7 +108,7 @@ void Foam::dynamicOversetFvMesh::correctBoundaryConditions
 
 
 template<class Type>
-Foam::tmp<Foam::scalarField> Foam::dynamicOversetFvMesh::normalisation
+Foam::tmp<Foam::scalarField> Foam::oversetFvMeshBase::normalisation
 (
     const fvMatrix<Type>& m
 ) const
@@ -165,7 +126,7 @@ Foam::tmp<Foam::scalarField> Foam::dynamicOversetFvMesh::normalisation
     {
         forAll(internalCoeffs, patchi)
         {
-            const labelUList& fc = lduAddr().patchAddr(patchi);
+            const labelUList& fc = mesh_.lduAddr().patchAddr(patchi);
             const Field<Type>& intCoeffs = internalCoeffs[patchi];
             const scalarField cmptCoeffs(intCoeffs.component(cmpt));
             forAll(fc, i)
@@ -199,9 +160,9 @@ Foam::tmp<Foam::scalarField> Foam::dynamicOversetFvMesh::normalisation
     {
         // Walk out the norm across hole cells
 
-        const labelList& own = faceOwner();
+        const labelList& own = mesh_.faceOwner();
-        const labelList& nei = faceNeighbour();
+        const labelList& nei = mesh_.faceNeighbour();
-        const cellCellStencilObject& overlap = Stencil::New(*this);
+        const cellCellStencilObject& overlap = Stencil::New(mesh_);
         const labelUList& types = overlap.cellTypes();
 
         label nHoles = 0;
@@ -215,8 +176,8 @@ Foam::tmp<Foam::scalarField> Foam::dynamicOversetFvMesh::normalisation
             }
         }
 
-        bitSet isFront(nFaces());
+        bitSet isFront(mesh_.nFaces());
-        for (label facei = 0; facei < nInternalFaces(); facei++)
+        for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
         {
             label ownType = types[own[facei]];
             label neiType = types[nei[facei]];
@@ -230,11 +191,16 @@ Foam::tmp<Foam::scalarField> Foam::dynamicOversetFvMesh::normalisation
             }
         }
         labelList nbrTypes;
-        syncTools::swapBoundaryCellList(*this, types, nbrTypes);
+        syncTools::swapBoundaryCellList(mesh_, types, nbrTypes);
-        for (label facei = nInternalFaces(); facei < nFaces(); facei++)
+        for
+        (
+            label facei = mesh_.nInternalFaces();
+            facei < mesh_.nFaces();
+            ++facei
+        )
         {
-            label ownType = types[own[facei]];
+            const label ownType = types[own[facei]];
-            label neiType = nbrTypes[facei-nInternalFaces()];
+            const label neiType = nbrTypes[facei-mesh_.nInternalFaces()];
             if
             (
                 (ownType == cellCellStencil::HOLE)
@@ -249,9 +215,9 @@ Foam::tmp<Foam::scalarField> Foam::dynamicOversetFvMesh::normalisation
         while (true)
         {
             scalarField nbrNorm;
-            syncTools::swapBoundaryCellList(*this, extrapolatedNorm, nbrNorm);
+            syncTools::swapBoundaryCellList(mesh_, extrapolatedNorm, nbrNorm);
 
-            bitSet newIsFront(nFaces());
+            bitSet newIsFront(mesh_.nFaces());
             scalarField newNorm(extrapolatedNorm);
 
             label nChanged = 0;
@@ -273,7 +239,7 @@ Foam::tmp<Foam::scalarField> Foam::dynamicOversetFvMesh::normalisation
                 }
                 if
                 (
-                    isInternalFace(facei)
+                    mesh_.isInternalFace(facei)
                  && extrapolatedNorm[nei[facei]] == -GREAT
                 )
                 {
@@ -299,7 +265,7 @@ Foam::tmp<Foam::scalarField> Foam::dynamicOversetFvMesh::normalisation
             // Transfer new front
             extrapolatedNorm.transfer(newNorm);
             isFront.transfer(newIsFront);
-            syncTools::syncFaceList(*this, isFront, maxEqOp<unsigned int>());
+            syncTools::syncFaceList(mesh_, isFront, maxEqOp<unsigned int>());
         }
 
 
@@ -324,16 +290,17 @@ Foam::tmp<Foam::scalarField> Foam::dynamicOversetFvMesh::normalisation
 
 
 template<class Type>
-void Foam::dynamicOversetFvMesh::addInterpolation
+void Foam::oversetFvMeshBase::addInterpolation
 (
     fvMatrix<Type>& m,
-    const scalarField& normalisation
+    const scalarField& normalisation,
+    const bool setHoleCellValue,
+    const Type& holeCellValue
 ) const
 {
-    const cellCellStencilObject& overlap = Stencil::New(*this);
+    const cellCellStencilObject& overlap = Stencil::New(mesh_);
     const List<scalarList>& wghts = overlap.cellInterpolationWeights();
     const labelListList& stencil = overlap.cellStencil();
-    const labelList& cellIDs = overlap.interpolationCells();
     const scalarList& factor = overlap.cellInterpolationWeight();
     const labelUList& types = overlap.cellTypes();
 
@@ -350,8 +317,6 @@ void Foam::dynamicOversetFvMesh::addInterpolation
     const lduAddressing& addr = lduAddr();
     const labelUList& upperAddr = addr.upperAddr();
     const labelUList& lowerAddr = addr.lowerAddr();
-    const labelUList& ownerStartAddr = addr.ownerStartAddr();
-    const labelUList& losortAddr = addr.losortAddr();
     const lduInterfacePtrsList& interfaces = allInterfaces_;
 
     if (!isA<fvMeshPrimitiveLduAddressing>(addr))
@@ -362,7 +327,6 @@ void Foam::dynamicOversetFvMesh::addInterpolation
     }
 
 
-
     // 1. Adapt lduMatrix for additional faces and new ordering
     upper.setSize(upperAddr.size(), 0.0);
     inplaceReorder(reverseFaceMap_, upper);
@@ -371,7 +335,7 @@ void Foam::dynamicOversetFvMesh::addInterpolation
 
 
     //const label nOldInterfaces = dynamicMotionSolverFvMesh::interfaces().size();
-    const label nOldInterfaces = dynamicFvMesh::interfaces().size();
+    const label nOldInterfaces = mesh_.fvMesh::interfaces().size();
 
 
     if (interfaces.size() > nOldInterfaces)
@@ -443,7 +407,6 @@ void Foam::dynamicOversetFvMesh::addInterpolation
     // fvMatrix:correction, both of which are outside the linear solver.
 
 
-
     // Clear out existing connections on cells to be interpolated
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Note: could avoid doing the zeroing of the new faces since these
@@ -451,6 +414,11 @@ void Foam::dynamicOversetFvMesh::addInterpolation
 
     forAll(upperAddr, facei)
     {
+        const label l = lowerAddr[facei];
+        scaleConnection(upper, types, factor, setHoleCellValue, l, facei);
+        const label u = upperAddr[facei];
+        scaleConnection(lower, types, factor, setHoleCellValue, u, facei);
+        /*
         if (types[upperAddr[facei]] == cellCellStencil::INTERPOLATED)
         {
             // Disconnect upper from lower
@@ -463,10 +431,21 @@ void Foam::dynamicOversetFvMesh::addInterpolation
             label celli = lowerAddr[facei];
             upper[facei] *= 1.0-factor[celli];
         }
+        */
     }
 
     for (label patchi = 0; patchi < nOldInterfaces; ++patchi)
     {
+        const labelUList& fc = addr.patchAddr(patchi);
+        Field<Type>& bCoeffs = m.boundaryCoeffs()[patchi];
+        Field<Type>& iCoeffs = m.internalCoeffs()[patchi];
+        forAll(fc, i)
+        {
+            scaleConnection(bCoeffs, types, factor, setHoleCellValue, fc[i], i);
+
+            scaleConnection(iCoeffs, types, factor, setHoleCellValue, fc[i], i);
+        }
+        /*
         const labelUList& fc = addr.patchAddr(patchi);
         Field<Type>& intCoeffs = m.internalCoeffs()[patchi];
         Field<Type>& bouCoeffs = m.boundaryCoeffs()[patchi];
@@ -487,6 +466,7 @@ void Foam::dynamicOversetFvMesh::addInterpolation
                 }
             }
         }
+        */
     }
 
 
@@ -497,8 +477,88 @@ void Foam::dynamicOversetFvMesh::addInterpolation
     // Do hole cells. Note: maybe put into interpolationCells() loop above?
     forAll(types, celli)
     {
-        if (types[celli] == cellCellStencil::HOLE)
+        if
+        (
+            types[celli] == cellCellStencil::HOLE
+         || (setHoleCellValue && types[celli] == cellCellStencil::SPECIAL)
+        )
         {
+            const Type wantedValue
+            (
+                setHoleCellValue
+              ? holeCellValue
+              : m.psi()[celli]
+            );
+            diag[celli] = normalisation[celli];
+            source[celli] = normalisation[celli]*wantedValue;
+        }
+        else if
+        (
+            types[celli] == cellCellStencil::INTERPOLATED
+         || (!setHoleCellValue && types[celli] == cellCellStencil::SPECIAL)
+        )
+        {
+            const scalar f = factor[celli];
+            const scalarList& w = wghts[celli];
+            const labelList& nbrs = stencil[celli];
+            const labelList& nbrFaces = stencilFaces_[celli];
+            const labelList& nbrPatches = stencilPatches_[celli];
+
+            diag[celli] *= (1.0-f);
+            source[celli] *= (1.0-f);
+
+            forAll(nbrs, nbri)
+            {
+                const label patchi = nbrPatches[nbri];
+                const label facei = nbrFaces[nbri];
+
+                if (patchi == -1)
+                {
+                    const label nbrCelli = nbrs[nbri];
+                    // Add the coefficients
+                    const scalar s = normalisation[celli]*f*w[nbri];
+
+                    scalar& u = upper[facei];
+                    scalar& l = lower[facei];
+                    if (celli < nbrCelli)
+                    {
+                        diag[celli] += s;
+                        u += -s;
+                    }
+                    else
+                    {
+                        diag[celli] += s;
+                        l += -s;
+                    }
+                }
+                else
+                {
+                    // Patch face. Store in boundaryCoeffs. Note sign change.
+                    //const label globalCelli = globalCellIDs[nbrs[nbri]];
+                    //const label proci =
+                    //    globalNumbering.whichProcID(globalCelli);
+                    //const label remoteCelli =
+                    //    globalNumbering.toLocal(proci, globalCelli);
+                    //
+                    //Pout<< "for cell:" << celli
+                    //    << " need weight from remote slot:" << nbrs[nbri]
+                    //    << " proc:" << proci << " remote cell:" << remoteCelli
+                    //    << " patch:" << patchi
+                    //    << " patchFace:" << facei
+                    //    << " weight:" << w[nbri]
+                    //    << endl;
+
+                    const scalar s = normalisation[celli]*f*w[nbri];
+                    m.boundaryCoeffs()[patchi][facei] += pTraits<Type>::one*s;
+                    m.internalCoeffs()[patchi][facei] += pTraits<Type>::one*s;
+
+                    // Note: do NOT add to diagonal - this is in the
+                    //       internalCoeffs and gets added to the diagonal
+                    //       inside fvMatrix::solve
+                }
+            }
+        }
+            /*
             label startLabel = ownerStartAddr[celli];
             label endLabel = ownerStartAddr[celli + 1];
 
@@ -518,7 +578,7 @@ void Foam::dynamicOversetFvMesh::addInterpolation
 
             diag[celli] = normalisation[celli];
             source[celli] = normalisation[celli]*m.psi()[celli];
-        }
+            */
     }
 
 
@@ -530,7 +590,7 @@ void Foam::dynamicOversetFvMesh::addInterpolation
     //}
     //overlap.cellInterpolationMap().distribute(globalCellIDs);
 
-
+/*
     forAll(cellIDs, i)
     {
         label celli = cellIDs[i];
@@ -636,11 +696,12 @@ void Foam::dynamicOversetFvMesh::addInterpolation
             //}
         }
     }
+*/
 }
 
 
 template<class Type>
-Foam::SolverPerformance<Type> Foam::dynamicOversetFvMesh::solve
+Foam::SolverPerformance<Type> Foam::oversetFvMeshBase::solveOverset
 (
     fvMatrix<Type>& m,
     const dictionary& dict
@@ -665,17 +726,17 @@ Foam::SolverPerformance<Type> Foam::dynamicOversetFvMesh::solve
     {
         if (debug)
         {
-            Pout<< "dynamicOversetFvMesh::solve() :"
+            Pout<< "oversetFvMeshBase::solveOverset() :"
                 << " bypassing matrix adjustment for field " << m.psi().name()
                 << endl;
         }
         //return dynamicMotionSolverFvMesh::solve(m, dict);
-        return dynamicFvMesh::solve(m, dict);
+        return mesh_.fvMesh::solve(m, dict);
     }
 
     if (debug)
     {
-        Pout<< "dynamicOversetFvMesh::solve() :"
+        Pout<< "oversetFvMeshBase::solveOverset() :"
             << " adjusting matrix for interpolation for field "
             << m.psi().name() << endl;
     }
@@ -683,20 +744,20 @@ Foam::SolverPerformance<Type> Foam::dynamicOversetFvMesh::solve
     // Calculate stabilised diagonal as normalisation for interpolation
     const scalarField norm(normalisation(m));
 
-    if (debug)
+    if (debug && mesh_.time().outputTime())
     {
         volScalarField scale
         (
             IOobject
             (
                 m.psi().name() + "_normalisation",
-                this->time().timeName(),
+                mesh_.time().timeName(),
-                *this,
+                mesh_,
                 IOobject::NO_READ,
                 IOobject::NO_WRITE,
                 false
             ),
-            *this,
+            mesh_,
             dimensionedScalar(dimless, Zero)
         );
         scale.ref().field() = norm;
@@ -709,7 +770,7 @@ Foam::SolverPerformance<Type> Foam::dynamicOversetFvMesh::solve
 
         if (debug)
         {
-            Pout<< "dynamicOversetFvMesh::solve() :"
+            Pout<< "oversetFvMeshBase::solveOverset() :"
                 << " writing matrix normalisation for field " << m.psi().name()
                 << " to " << scale.name() << endl;
         }
@@ -725,9 +786,15 @@ Foam::SolverPerformance<Type> Foam::dynamicOversetFvMesh::solve
     scalarField oldLower(m.lower());
     FieldField<Field, Type> oldInt(m.internalCoeffs());
     FieldField<Field, Type> oldBou(m.boundaryCoeffs());
+
+    Field<Type> oldSource(m.source());
+    scalarField oldDiag(m.diag());
+
     const label oldSize = bpsi.size();
 
-    addInterpolation(m, norm);
+    // Insert the interpolation into the matrix (done inside
+    // oversetFvPatchField<Type>::manipulateMatrix)
+    m.boundaryManipulate(bpsi);
 
     // Swap psi values so added patches have patchNeighbourField
     correctBoundaryConditions<GeoField, calculatedProcessorFvPatchField<Type>>
@@ -736,18 +803,9 @@ Foam::SolverPerformance<Type> Foam::dynamicOversetFvMesh::solve
         true
     );
 
-
-    // Print a bit
-    //write(Pout, m, lduAddr(), interfaces());
-    //{
-    //   const fvSolution& sol = static_cast<const fvSolution&>(*this);
-    //    const dictionary& pDict = sol.subDict("solvers").subDict("p");
-    //    writeAgglomeration(GAMGAgglomeration::New(m, pDict));
-    //}
-
     // Use lower level solver
     //SolverPerformance<Type> s(dynamicMotionSolverFvMesh::solve(m, dict));
-    SolverPerformance<Type> s(dynamicFvMesh::solve(m, dict));
+    SolverPerformance<Type> s(mesh_.fvMesh::solve(m, dict));
 
     // Restore boundary field
     bpsi.setSize(oldSize);
@@ -755,9 +813,45 @@ Foam::SolverPerformance<Type> Foam::dynamicOversetFvMesh::solve
     // Restore matrix
     m.upper().transfer(oldUpper);
     m.lower().transfer(oldLower);
+
+    m.source().transfer(oldSource);
+    m.diag().transfer(oldDiag);
+
     m.internalCoeffs().transfer(oldInt);
     m.boundaryCoeffs().transfer(oldBou);
 
+
+    const cellCellStencilObject& overlap = Stencil::New(mesh_);
+    const labelUList& types = overlap.cellTypes();
+    const labelList& own = mesh_.faceOwner();
+    const labelList& nei = mesh_.faceNeighbour();
+
+    auto& psi =
+        const_cast<GeometricField<Type, fvPatchField, volMesh>&>(m.psi());
+
+    // Mirror hole cell values next to calculated cells
+    for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
+    {
+        const label ownType = types[own[facei]];
+        const label neiType = types[nei[facei]];
+
+        if
+        (
+            ownType == cellCellStencil::HOLE
+         && neiType == cellCellStencil::CALCULATED)
+        {
+            psi[own[facei]] = psi[nei[facei]];
+        }
+        else if
+        (
+            ownType == cellCellStencil::CALCULATED
+         && neiType == cellCellStencil::HOLE
+        )
+        {
+             psi[nei[facei]] = psi[own[facei]];
+        }
+    }
+
     // Switch to original addressing
     active(false);
 
@@ -766,7 +860,7 @@ Foam::SolverPerformance<Type> Foam::dynamicOversetFvMesh::solve
 
 
 template<class Type>
-void Foam::dynamicOversetFvMesh::write
+void Foam::oversetFvMeshBase::write
 (
     Ostream& os,
     const fvMatrix<Type>& m,
@@ -893,7 +987,7 @@ void Foam::dynamicOversetFvMesh::write
 
 
 template<class GeoField>
-void Foam::dynamicOversetFvMesh::correctCoupledBoundaryConditions(GeoField& fld)
+void Foam::oversetFvMeshBase::correctCoupledBoundaryConditions(GeoField& fld)
 {
     typename GeoField::Boundary& bfld = fld.boundaryFieldRef();
 
@@ -912,7 +1006,7 @@ void Foam::dynamicOversetFvMesh::correctCoupledBoundaryConditions(GeoField& fld)
     // Wait for outstanding requests
     if (commsType == UPstream::commsTypes::nonBlocking)
     {
-        UPstream::waitRequests(startOfRequests);
+        Pstream::waitRequests(startOfRequests);
     }
 
     forAll(bfld, patchi)
@@ -927,7 +1021,7 @@ void Foam::dynamicOversetFvMesh::correctCoupledBoundaryConditions(GeoField& fld)
 
 
 template<class GeoField>
-void Foam::dynamicOversetFvMesh::checkCoupledBC(const GeoField& fld)
+void Foam::oversetFvMeshBase::checkCoupledBC(const GeoField& fld)
 {
     Pout<< "** starting checking of " << fld.name() << endl;
 

src/overset/oversetFvMesh/staticOversetFvMesh/staticOversetFvMesh.C

@@ -0,0 +1,84 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "staticOversetFvMesh.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(staticOversetFvMesh, 0);
+    addToRunTimeSelectionTable(dynamicFvMesh, staticOversetFvMesh, IOobject);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::staticOversetFvMesh::staticOversetFvMesh(const IOobject& io)
+:
+    staticFvMesh(io),
+    oversetFvMeshBase(static_cast<const fvMesh&>(*this))
+{
+    oversetFvMeshBase::update();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::staticOversetFvMesh::~staticOversetFvMesh()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::staticOversetFvMesh::update()
+{
+    if (!staticFvMesh::update())
+    {
+        return false;
+    }
+
+    oversetFvMeshBase::update();
+
+    return true;
+}
+
+
+bool Foam::staticOversetFvMesh::writeObject
+(
+    IOstreamOption streamOpt,
+    const bool valid
+) const
+{
+    bool ok = staticFvMesh::writeObject(streamOpt, valid);
+    ok = oversetFvMeshBase::writeObject(streamOpt, valid) && ok;
+    return ok;
+}
+
+
+// ************************************************************************* //

src/overset/oversetFvMesh/staticOversetFvMesh/staticOversetFvMesh.H

@@ -0,0 +1,225 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::staticOversetFvMesh
+
+Description
+    fvMesh with support for overset meshes.
+
+SourceFiles
+    staticOversetFvMesh.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef staticOversetFvMesh_H
+#define staticOversetFvMesh_H
+
+#include "staticFvMesh.H"
+#include "oversetFvMeshBase.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                   Class staticOversetFvMesh Declaration
+\*---------------------------------------------------------------------------*/
+
+class staticOversetFvMesh
+:
+    public staticFvMesh,
+    public oversetFvMeshBase
+{
+    // Private Member Functions
+
+        //- No copy construct
+        staticOversetFvMesh(const staticOversetFvMesh&) = delete;
+
+        //- No copy assignment
+        void operator=(const staticOversetFvMesh&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("staticOversetFvMesh");
+
+
+    // Constructors
+
+        //- Construct from IOobject
+        staticOversetFvMesh(const IOobject& io);
+
+
+    //- Destructor
+    virtual ~staticOversetFvMesh();
+
+
+    // Member Functions
+
+        //- Override ldu addressing
+        virtual const lduAddressing& lduAddr() const
+        {
+            return oversetFvMeshBase::lduAddr();
+        }
+
+        //- Override ldu addressing
+        virtual lduInterfacePtrsList interfaces() const
+        {
+            return oversetFvMeshBase::interfaces();
+        }
+
+
+        // Overset
+
+            // Explicit interpolation
+
+                virtual void interpolate(scalarField& psi) const
+                {
+                    Stencil::New(*this).interpolate<scalar>(*this, psi);
+                }
+
+                virtual void interpolate(vectorField& psi) const
+                {
+                    Stencil::New(*this).interpolate<vector>(*this, psi);
+                }
+
+                virtual void interpolate(sphericalTensorField& psi) const
+                {
+                    Stencil::New
+                    (
+                        *this
+                    ).interpolate<sphericalTensor>(*this, psi);
+                }
+
+                virtual void interpolate(symmTensorField& psi) const
+                {
+                    Stencil::New(*this).interpolate<symmTensor>(*this, psi);
+                }
+
+                virtual void interpolate(tensorField& psi) const
+                {
+                    Stencil::New(*this).interpolate<tensor>(*this, psi);
+                }
+
+                virtual void interpolate(volScalarField& psi) const
+                {
+                    Stencil::New(*this).interpolate<volScalarField>(psi);
+                }
+
+                virtual void interpolate(volVectorField& psi) const
+                {
+                    Stencil::New(*this).interpolate<volVectorField>(psi);
+                }
+
+                virtual void interpolate(volSphericalTensorField& psi) const
+                {
+                    Stencil::New(*this).interpolate
+                    <
+                        volSphericalTensorField
+                    >(psi);
+                }
+
+                virtual void interpolate(volSymmTensorField& psi) const
+                {
+                    Stencil::New(*this).interpolate<volSymmTensorField>(psi);
+                }
+
+                virtual void interpolate(volTensorField& psi) const
+                {
+                    Stencil::New(*this).interpolate<volTensorField>(psi);
+                }
+
+
+            // Implicit interpolation (matrix manipulation)
+
+                //- Solve returning the solution statistics given convergence
+                //- tolerance. Use the given solver controls
+                virtual SolverPerformance<scalar> solve
+                (
+                    fvMatrix<scalar>& m,
+                    const dictionary& dict
+                ) const
+                {
+                    return oversetFvMeshBase::solveOverset<scalar>(m, dict);
+                }
+
+                //- Solve returning the solution statistics given convergence
+                //- tolerance. Use the given solver controls
+                virtual SolverPerformance<vector> solve
+                (
+                    fvMatrix<vector>& m,
+                    const dictionary& dict
+                ) const
+                {
+                    return oversetFvMeshBase::solveOverset<vector>(m, dict);
+                }
+
+                //- Solve returning the solution statistics given convergence
+                //- tolerance. Use the given solver controls
+                virtual SolverPerformance<symmTensor> solve
+                (
+                    fvMatrix<symmTensor>& m,
+                    const dictionary& dict
+                ) const
+                {
+                    return oversetFvMeshBase::solveOverset<symmTensor>(m, dict);
+                }
+
+                //- Solve returning the solution statistics given convergence
+                //- tolerance. Use the given solver controls
+                virtual SolverPerformance<tensor> solve
+                (
+                    fvMatrix<tensor>& m,
+                    const dictionary& dict
+                ) const
+                {
+                    return oversetFvMeshBase::solveOverset<tensor>(m, dict);
+                }
+
+
+        //- Update the mesh for both mesh motion and topology change
+        virtual bool update();
+
+        //- Write using given format, version and compression
+        virtual bool writeObject
+        (
+            IOstreamOption streamOpt,
+            const bool valid
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/overset/oversetPatchPhiErr/oversetPatchPhiErr.C

@@ -0,0 +1,58 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "oversetPatchPhiErr.H"
+#include "oversetFvPatchField.H"
+#include "volFields.H"
+#include "volFieldsFwd.H"
+
+// * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
+
+void Foam::oversetPatchPhiErr
+(
+    const fvScalarMatrix& m,
+    const surfaceScalarField& phi
+)
+{
+    const volScalarField::Boundary& bm = m.psi().boundaryField();
+
+    forAll(bm, patchi)
+    {
+        const auto& fvp = bm[patchi];
+
+        if (isA<oversetFvPatchField<scalar>>(fvp))
+        {
+            const oversetFvPatchField<scalar>& op =
+                refCast<const oversetFvPatchField<scalar>>(fvp);
+
+            op.fringeFlux(m, phi);
+        }
+    }
+}
+
+
+// ************************************************************************* //

src/overset/oversetPatchPhiErr/oversetPatchPhiErr.H

@@ -0,0 +1,62 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+InNamespace
+    Foam
+
+Description
+    oversetPatchPhiErr
+
+SourceFiles
+    oversetPatchPhiErr.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef oversetPatchPhiErr_H
+#define oversetPatchPhiErr_H
+
+#include "volFieldsFwd.H"
+#include "surfaceFieldsFwd.H"
+#include "fvScalarMatrix.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+//- Prints out sum of m.flux and phi over the fringe faces
+void oversetPatchPhiErr
+(
+    const fvScalarMatrix& m,
+    const surfaceScalarField& phi
+);
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/overset/oversetPolyPatch/oversetFvPatch.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2019 OpenCFD Ltd.
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -41,4 +41,65 @@ namespace Foam
 }
 
 
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::labelField> Foam::oversetFvPatch::interfaceInternalField
+(
+    const labelUList& internalData
+) const
+{
+    // Return all internal values
+    return patchInternalField(internalData);
+}
+
+
+Foam::tmp<Foam::labelField> Foam::oversetFvPatch::interfaceInternalField
+(
+    const labelUList& internalData,
+    const labelUList& faceCells
+) const
+{
+    return patchInternalField(internalData, faceCells);
+}
+
+
+void Foam::oversetFvPatch::initInternalFieldTransfer
+(
+    const Pstream::commsTypes commsType,
+    const labelUList& iF
+) const
+{}
+
+
+void Foam::oversetFvPatch::initInternalFieldTransfer
+(
+    const Pstream::commsTypes commsType,
+    const labelUList& iF,
+    const labelUList& faceCells
+) const
+{}
+
+
+Foam::tmp<Foam::labelField> Foam::oversetFvPatch::internalFieldTransfer
+(
+    const Pstream::commsTypes commsType,
+    const labelUList& iF
+) const
+{
+    return tmp<labelField>::New(iF);
+}
+
+
+Foam::tmp<Foam::labelField> Foam::oversetFvPatch::internalFieldTransfer
+(
+    const Pstream::commsTypes commsType,
+    const labelUList& iF,
+    const labelUList& nbrFaceCells
+) const
+{
+    NotImplemented
+    return tmp<labelField>::New(iF);
+}
+
+
 // ************************************************************************* //

src/overset/oversetPolyPatch/oversetFvPatch.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2019 OpenCFD Ltd.
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -37,6 +37,9 @@ SourceFiles
 #ifndef oversetFvPatch_H
 #define oversetFvPatch_H
 
+#include "coupledFvPatch.H"
+#include "oversetLduInterface.H"
+#include "lduInterface.H"
 #include "fvPatch.H"
 #include "oversetPolyPatch.H"
 #include "fvBoundaryMesh.H"
@@ -52,9 +55,11 @@ namespace Foam
 
 class oversetFvPatch
 :
-    public fvPatch
+    public fvPatch,
+    public lduInterface,
+    public oversetLduInterface
 {
-    // Private data
+    // Private Data
 
         const oversetPolyPatch& oversetPolyPatch_;
 
@@ -75,7 +80,7 @@ public:
         {}
 
 
-    // Member functions
+    // Member Functions
 
         // Access
 
@@ -85,9 +90,6 @@ public:
                 return oversetPolyPatch_;
             }
 
-
-        // Access
-
             //- Return faceCell addressing
             virtual const labelUList& faceCells() const
             {
@@ -99,6 +101,56 @@ public:
             {
                 return oversetPolyPatch_.master();
             }
+
+
+        // Interface Transfer Functions
+
+            //- Return the values of the given internal data adjacent to
+            //- the interface as a field
+            virtual tmp<labelField> interfaceInternalField
+            (
+                const labelUList& internalData
+            ) const;
+
+            //- Return the values of the given internal data adjacent to
+            //- the interface as a field using faceCell mapping
+            virtual tmp<labelField> interfaceInternalField
+            (
+                const labelUList& internalData,
+                const labelUList& faceCells
+            ) const;
+
+            //- Initialise transfer of internal field adjacent to the interface
+            virtual void initInternalFieldTransfer
+            (
+                const Pstream::commsTypes commsType,
+                const labelUList& iF
+            ) const;
+
+            //- Initialise transfer of internal field adjacent to the interface
+            //- using faceCells mapping
+            virtual void initInternalFieldTransfer
+            (
+                const Pstream::commsTypes commsType,
+                const labelUList& iF,
+                const labelUList& faceCells
+            ) const;
+
+            //- Transfer and return internal field adjacent to the interface
+            virtual tmp<labelField> internalFieldTransfer
+            (
+                const Pstream::commsTypes commsType,
+                const labelUList& iF
+            ) const;
+
+            //- Transfer and return internal field adjacent to the interface
+            //- using nbrFaceCells mapping
+            virtual tmp<labelField> internalFieldTransfer
+            (
+                const Pstream::commsTypes commsType,
+                const labelUList& iF,
+                const labelUList& nbrFaceCells
+            ) const;
 };
 
 

src/overset/oversetPolyPatch/oversetFvPatchField.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2018 OpenCFD Ltd.
+    Copyright (C) 2016-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,7 +28,8 @@ License
 #include "volFields.H"
 #include "cellCellStencil.H"
 #include "cellCellStencilObject.H"
-#include "dynamicOversetFvMesh.H"
+#include "oversetFvMeshBase.H"
+#include "syncTools.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -39,8 +40,16 @@ Foam::oversetFvPatchField<Type>::oversetFvPatchField
     const DimensionedField<Type, volMesh>& iF
 )
 :
-    zeroGradientFvPatchField<Type>(p, iF),
+    coupledFvPatchField<Type>(p, iF),
-    oversetPatch_(refCast<const oversetFvPatch>(p))
+    oversetPatch_(refCast<const oversetFvPatch>(p)),
+    setHoleCellValue_(false),
+    fluxCorrection_(false),
+    interpolateHoleCellValue_(false),
+    holeCellValue_(pTraits<Type>::min),
+    fringeUpperCoeffs_(),
+    fringeLowerCoeffs_(),
+    fringeFaces_(),
+    zoneId_(-1)
 {}
 
 
@@ -53,8 +62,16 @@ Foam::oversetFvPatchField<Type>::oversetFvPatchField
     const fvPatchFieldMapper& mapper
 )
 :
-    zeroGradientFvPatchField<Type>(ptf, p, iF, mapper),
+    coupledFvPatchField<Type>(ptf, p, iF, mapper),
-    oversetPatch_(refCast<const oversetFvPatch>(p))
+    oversetPatch_(refCast<const oversetFvPatch>(p)),
+    setHoleCellValue_(ptf.setHoleCellValue_),
+    fluxCorrection_(ptf.fluxCorrection_),
+    interpolateHoleCellValue_(ptf.interpolateHoleCellValue_),
+    holeCellValue_(ptf.holeCellValue_),
+    fringeUpperCoeffs_(ptf.fringeUpperCoeffs_),
+    fringeLowerCoeffs_(ptf.fringeLowerCoeffs_),
+    fringeFaces_(ptf.fringeFaces_),
+    zoneId_(ptf.zoneId_)
 {}
 
 
@@ -66,9 +83,45 @@ Foam::oversetFvPatchField<Type>::oversetFvPatchField
     const dictionary& dict
 )
 :
-    zeroGradientFvPatchField<Type>(p, iF, dict),
+    coupledFvPatchField<Type>(p, iF, dict, false),
-    oversetPatch_(refCast<const oversetFvPatch>(p, dict))
+    oversetPatch_(refCast<const oversetFvPatch>(p, dict)),
-{}
+    setHoleCellValue_(dict.getOrDefault("setHoleCellValue", false)),
+    fluxCorrection_
+    (
+        dict.getOrDefaultCompat
+        (
+            "fluxCorrection",
+            {{"massCorrection", 2206}},
+            false
+        )
+    ),
+    interpolateHoleCellValue_
+    (
+        dict.getOrDefault("interpolateHoleCellValue", false)
+    ),
+    holeCellValue_
+    (
+        setHoleCellValue_
+      ? dict.get<Type>("holeCellValue")
+      : pTraits<Type>::min
+    ),
+    fringeUpperCoeffs_(),
+    fringeLowerCoeffs_(),
+    fringeFaces_(),
+    zoneId_(dict.getOrDefault<label>("zone", -1))
+{
+    if (dict.found("value"))
+    {
+        Field<Type>::operator=
+        (
+            Field<Type>("value", dict, p.size())
+        );
+    }
+    else
+    {
+        Field<Type>::operator=(this->patchInternalField());
+    }
+}
 
 
 template<class Type>
@@ -77,8 +130,16 @@ Foam::oversetFvPatchField<Type>::oversetFvPatchField
     const oversetFvPatchField<Type>& ptf
 )
 :
-    zeroGradientFvPatchField<Type>(ptf),
+    coupledFvPatchField<Type>(ptf),
-    oversetPatch_(ptf.oversetPatch_)
+    oversetPatch_(ptf.oversetPatch_),
+    setHoleCellValue_(ptf.setHoleCellValue_),
+    fluxCorrection_(ptf.fluxCorrection_),
+    interpolateHoleCellValue_(ptf.interpolateHoleCellValue_),
+    holeCellValue_(ptf.holeCellValue_),
+    fringeUpperCoeffs_(ptf.fringeUpperCoeffs_),
+    fringeLowerCoeffs_(ptf.fringeLowerCoeffs_),
+    fringeFaces_(ptf.fringeFaces_),
+    zoneId_(ptf.zoneId_)
 {}
 
 
@@ -89,13 +150,454 @@ Foam::oversetFvPatchField<Type>::oversetFvPatchField
     const DimensionedField<Type, volMesh>& iF
 )
 :
-    zeroGradientFvPatchField<Type>(ptf, iF),
+    coupledFvPatchField<Type>(ptf, iF),
-    oversetPatch_(ptf.oversetPatch_)
+    oversetPatch_(ptf.oversetPatch_),
+    setHoleCellValue_(ptf.setHoleCellValue_),
+    fluxCorrection_(ptf.fluxCorrection_),
+    interpolateHoleCellValue_(ptf.interpolateHoleCellValue_),
+    holeCellValue_(ptf.holeCellValue_),
+    fringeUpperCoeffs_(ptf.fringeUpperCoeffs_),
+    fringeLowerCoeffs_(ptf.fringeLowerCoeffs_),
+    fringeFaces_(ptf.fringeFaces_),
+    zoneId_(ptf.zoneId_)
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class Type>
+void Foam::oversetFvPatchField<Type>::storeFringeCoefficients
+(
+    const fvMatrix<Type>& matrix
+)
+{
+    const fvMesh& mesh = this->internalField().mesh();
+
+    const cellCellStencilObject& overlap = Stencil::New(mesh);
+    const labelList& cellTypes = overlap.cellTypes();
+    const labelList& zoneID = overlap.zoneID();
+    const labelUList& own = mesh.owner();
+    const labelUList& nei = mesh.neighbour();
+
+    const polyBoundaryMesh& boundaryMesh = mesh.boundaryMesh();
+
+    label fringesFaces = 0;
+
+    forAll(own, facei)
+    {
+        const label zonei = zoneID[own[facei]];
+
+        const label ownType = cellTypes[own[facei]];
+        const label neiType = cellTypes[nei[facei]];
+
+        const bool ownCalc =
+            (ownType == cellCellStencil::CALCULATED)
+         && (neiType == cellCellStencil::INTERPOLATED);
+
+        const bool neiCalc =
+            (ownType == cellCellStencil::INTERPOLATED)
+         && (neiType == cellCellStencil::CALCULATED);
+
+        const bool ownNei = (ownCalc || neiCalc);
+
+        if
+        (
+            (ownNei && (zonei == zoneId_)) || (ownNei && (zoneId_ == -1))
+        )
+        {
+            fringesFaces++;
+        }
+    }
+
+    const fvPatchList& patches = mesh.boundary();
+
+    labelList neiCellTypes;
+    syncTools::swapBoundaryCellList(mesh, cellTypes, neiCellTypes);
+    {
+        forAll(patches, patchi)
+        {
+            const fvPatch& curPatch = patches[patchi];
+
+            const labelUList& fc = curPatch.faceCells();
+
+            const label start = curPatch.start();
+
+            forAll(fc, i)
+            {
+                const label facei = start + i;
+                const label celli = fc[i];
+                const label ownType = cellTypes[celli];
+                const label neiType = neiCellTypes[facei-mesh.nInternalFaces()];
+
+                const label zonei = zoneID[celli];
+
+                const bool ownCalc =
+                    (ownType == cellCellStencil::CALCULATED)
+                 && (neiType == cellCellStencil::INTERPOLATED);
+
+
+                if (ownCalc && (zonei == zoneId_))
+                {
+                    fringesFaces++;
+                }
+            }
+        }
+    }
+    fringeUpperCoeffs_.setSize(fringesFaces, Zero);
+    fringeLowerCoeffs_.setSize(fringesFaces, Zero);
+    fringeFaces_.setSize(fringesFaces, -1);
+
+    const scalarField& upper = matrix.upper();
+    const scalarField& lower = matrix.lower();
+
+    fringesFaces = 0;
+    forAll(own, facei)
+    {
+        const label zonei = zoneID[own[facei]];
+
+        const label ownType = cellTypes[own[facei]];
+        const label neiType = cellTypes[nei[facei]];
+
+        const bool ownCalc =
+            (ownType == cellCellStencil::CALCULATED)
+         && (neiType == cellCellStencil::INTERPOLATED);
+
+        const bool neiCalc =
+            (ownType == cellCellStencil::INTERPOLATED)
+         && (neiType == cellCellStencil::CALCULATED);
+
+        const bool ownNei = (ownCalc || neiCalc);
+
+        if
+        (
+            (ownNei && (zonei == zoneId_)) || (ownNei && (zoneId_ == -1))
+        )
+        {
+            fringeUpperCoeffs_[fringesFaces] = upper[facei];
+            fringeLowerCoeffs_[fringesFaces] = lower[facei];
+            fringeFaces_[fringesFaces] = facei;
+            fringesFaces++;
+        }
+    }
+
+    forAll(boundaryMesh, patchi)
+    {
+        const polyPatch& p = boundaryMesh[patchi];
+
+        if (isA<coupledPolyPatch>(p))
+        {
+            const labelUList& fc = p.faceCells();
+            const label start = p.start();
+
+            forAll(fc, i)
+            {
+                const label facei = start + i;
+                const label celli = fc[i];
+                const label ownType = cellTypes[celli];
+                const label neiType = neiCellTypes[facei-mesh.nInternalFaces()];
+
+                const label zonei = zoneID[celli];
+
+                const bool ownCalc =
+                    (ownType == cellCellStencil::CALCULATED)
+                 && (neiType == cellCellStencil::INTERPOLATED);
+
+                const bool neiCalc =
+                    (neiType == cellCellStencil::CALCULATED)
+                 && (ownType == cellCellStencil::INTERPOLATED);
+
+                if ((ownCalc||neiCalc)  && (zonei == zoneId_))
+                {
+                    fringeLowerCoeffs_[fringesFaces] =
+                        component
+                        (
+                            matrix.internalCoeffs()[patchi][facei],
+                            0
+                        );
+
+                    fringeUpperCoeffs_[fringesFaces] =
+                        component
+                        (
+                            matrix.boundaryCoeffs()[patchi][facei],
+                            0
+                        );
+
+                    fringeFaces_[fringesFaces] = facei;
+
+                    fringesFaces++;
+                }
+            }
+        }
+    }
+}
+
+
+template<class Type>
+void Foam::oversetFvPatchField<Type>::fringeFlux
+(
+    const fvMatrix<Type>& matrix,
+    const surfaceScalarField& phi
+) const
+{
+    scalar massIn = 0;
+    scalar phiIn = 0;
+
+    const Field<Type>& psi = matrix.psi();
+
+    const scalarField& upper = matrix.upper();
+    const scalarField& lower = matrix.lower();
+
+    if (this->oversetPatch_.master())
+    {
+        const fvMesh& mesh = this->internalField().mesh();
+        const cellCellStencilObject& overlap = Stencil::New(mesh);
+        const labelList& cellTypes = overlap.cellTypes();
+        const labelList& zoneID = overlap.zoneID();
+
+        // Check all faces on the outside of interpolated cells
+        const labelUList& own = mesh.owner();
+        const labelUList& nei = mesh.neighbour();
+
+        label fringesFaces = 0;
+        forAll(own, facei)
+        {
+            const label zonei = zoneID[own[facei]];
+
+            const label ownType = cellTypes[own[facei]];
+            const label neiType = cellTypes[nei[facei]];
+
+            const bool ownCalc =
+                (ownType == cellCellStencil::CALCULATED)
+             && (neiType == cellCellStencil::INTERPOLATED);
+
+            const bool neiCalc =
+                (ownType == cellCellStencil::INTERPOLATED)
+             && (neiType == cellCellStencil::CALCULATED);
+
+            const bool ownNei = (ownCalc || neiCalc);
+
+            if
+            (
+                (ownNei && (zonei == zoneId_)) || (ownNei && (zoneId_ == -1))
+            )
+            {
+                const label fringei = fringeFaces_[fringesFaces];
+
+                // Get fringe upper/lower coeffs
+                //const scalar& ufc = fringeUpperCoeffs_[fringesFaces];
+                //const scalar& lfc = fringeLowerCoeffs_[fringesFaces];
+
+                const scalar& ufc = upper[fringei];
+                const scalar& lfc = lower[fringei];
+
+                const Type curFlux =
+                    ufc*psi[nei[fringei]] - lfc*psi[own[fringei]];
+
+                if (neiCalc)
+                {
+                    phiIn -= phi[fringei];
+                    massIn -= curFlux;
+                }
+                else
+                {
+                    phiIn += phi[fringei];
+                    massIn += curFlux;
+                }
+                fringesFaces++;
+            }
+        }
+    }
+
+    reduce(massIn, sumOp<scalar>());
+    reduce(phiIn, sumOp<scalar>());
+
+    Info << " gSum(phi) on fringes " << phiIn << endl;
+    Info << " gSum(p.flux) on fringes " << massIn << endl;
+}
+
+
+template<class Type>
+void Foam::oversetFvPatchField<Type>::adjustPsi
+(
+    solveScalarField& psi,
+    const lduAddressing& lduAddr,
+    solveScalarField& result
+) const
+{
+    const fvMesh& mesh = this->internalField().mesh();
+
+    const cellCellStencilObject& overlap = Stencil::New(mesh);
+    const labelList& cellTypes = overlap.cellTypes();
+    const labelList& zoneID = overlap.zoneID();
+
+    // Pass-1: accumulate all fluxes, calculate correction factor
+
+    scalarField interpolatedCoeffs(fringeUpperCoeffs_.size(), Zero);
+
+    // Options for scaling corrections
+    scalar massIn = 0;
+    scalar offDiagCoeffs = 0;
+    labelField facePerCell(cellTypes.size(), 0);
+
+    // Check all faces on the outside of interpolated cells
+    const labelUList& own = mesh.owner();
+    const labelUList& nei = mesh.neighbour();
+
+    label fringesFaces = 0;
+    {
+        forAll(own, facei)
+        {
+            const label zonei = zoneID[own[facei]];
+
+            const label ownType = cellTypes[own[facei]];
+            const label neiType = cellTypes[nei[facei]];
+
+            const bool ownCalc =
+                (ownType == cellCellStencil::CALCULATED)
+             && (neiType == cellCellStencil::INTERPOLATED);
+
+            const bool neiCalc =
+                (ownType == cellCellStencil::INTERPOLATED)
+             && (neiType == cellCellStencil::CALCULATED);
+
+            const bool ownNei = (ownCalc || neiCalc);
+
+            if
+            (
+                (ownNei && (zonei == zoneId_)) || (ownNei && (zoneId_ == -1))
+            )
+            {
+                // Get fringe upper/lower coeffs
+                const scalar& ufc = fringeUpperCoeffs_[fringesFaces];
+                const scalar& lfc = fringeLowerCoeffs_[fringesFaces];
+
+                const scalar curFlux =
+                    ufc*psi[nei[facei]] - lfc*psi[own[facei]];
+
+                if (neiCalc) // interpolated is owner
+                {
+                    massIn -= curFlux;
+                    offDiagCoeffs += lfc;
+                    facePerCell[own[facei]]++;
+                }
+                else
+                {
+                    massIn += curFlux;
+                    offDiagCoeffs += ufc;
+                    facePerCell[nei[facei]]++;
+                }
+                fringesFaces++;
+            }
+        }
+    }
+
+    scalarField weights(facePerCell.size(), scalar(1));
+    forAll (weights, celli)
+    {
+        if (facePerCell[celli] > 1)
+        {
+            weights[celli] = scalar(1)/facePerCell[celli];
+        }
+    }
+
+    // Check all coupled faces on the outside of interpolated cells
+    labelList neiCellTypes;
+    syncTools::swapBoundaryCellList(mesh, cellTypes, neiCellTypes);
+
+    const fvPatchList& boundaryMesh = mesh.boundary();
+
+    forAll(boundaryMesh, patchi)
+    {
+        const polyPatch& p = mesh.boundaryMesh()[patchi];
+
+        if (isA<coupledPolyPatch>(p))
+        {
+            const auto& coupledPatch = refCast<const coupledPolyPatch>(p);
+
+            const labelUList& fc = p.faceCells();
+            const label start = p.start();
+
+            forAll(fc, i)
+            {
+                const label facei = start + i;
+                const label celli = fc[i];
+                const label ownType = cellTypes[celli];
+                const label neiType = neiCellTypes[facei-mesh.nInternalFaces()];
+
+                const label zonei = zoneID[celli];
+
+                const bool ownCalc =
+                    (ownType == cellCellStencil::CALCULATED)
+                 && (neiType == cellCellStencil::INTERPOLATED);
+
+                const bool neiCalc =
+                    (neiType == cellCellStencil::CALCULATED)
+                 && (ownType == cellCellStencil::INTERPOLATED);
+
+                if ((ownCalc||neiCalc)  && (zonei == zoneId_))
+                {
+                    const scalar psiOwn = psi[celli];
+                    const scalar& lfc = fringeLowerCoeffs_[fringesFaces];
+                    const scalar curFlux = lfc*psiOwn;
+
+                    if (ownCalc)
+                    {
+                        massIn -= curFlux;
+
+                        if (coupledPatch.owner())
+                        {
+                            offDiagCoeffs -= lfc;
+                        }
+
+                        fringesFaces++;
+                    }
+                    else
+                    {
+                        massIn += curFlux;
+
+                        if (coupledPatch.owner())
+                        {
+                            offDiagCoeffs -= lfc;
+                        }
+
+                        fringesFaces++;
+                    }
+                }
+            }
+        }
+    }
+
+    reduce(massIn, sumOp<scalar>());
+    reduce(offDiagCoeffs, sumOp<scalar>());
+
+    scalar psiCorr = -massIn/offDiagCoeffs;
+
+    forAll (cellTypes, celli)
+    {
+        const bool bInter = (cellTypes[celli] == cellCellStencil::INTERPOLATED);
+        const label zonei = zoneID[celli];
+        if
+        (
+            (bInter && (zonei == zoneId_)) ||(bInter && (zoneId_ == -1))
+        )
+        {
+            psi[celli] += psiCorr;
+        }
+    }
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>> Foam::oversetFvPatchField<Type>::
+patchNeighbourField() const
+{
+    return tmp<Field<Type>>
+    (
+        new Field<Type>(this->size(), Zero)
+    );
+}
+
+
 template<class Type>
 void Foam::oversetFvPatchField<Type>::initEvaluate
 (
@@ -205,25 +707,374 @@ void Foam::oversetFvPatchField<Type>::initEvaluate
                     Info<< "Interpolating non-suppressed field " << fldName
                         << endl;
                 }
-                mesh.interpolate
+
+                // 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.
+
+                const cellCellStencil& overlap = Stencil::New(mesh);
+
+                Field<Type>& fld =
+                    const_cast<Field<Type>&>(this->primitiveField());
+
+                // tbd: different weights for different variables ...
+                cellCellStencil::interpolate
                 (
-                    const_cast<Field<Type>&>
+                    fld,
-                    (
+                    mesh,
-                        this->primitiveField()
+                    overlap,
-                    )
+                    overlap.cellInterpolationWeights()
                 );
+
+                if (this->setHoleCellValue_)
+                {
+                    const labelUList& types = overlap.cellTypes();
+                    label nConstrained = 0;
+                    forAll(types, celli)
+                    {
+                        const label cType = types[celli];
+                        if
+                        (
+                            cType == cellCellStencil::HOLE
+                         || cType == cellCellStencil::SPECIAL
+                        )
+                        {
+                            fld[celli] = this->holeCellValue_;
+                            nConstrained++;
+                        }
+                    }
+
+                    if (debug)
+                    {
+                        Pout<< FUNCTION_NAME << " field:" << fldName
+                            << " patch:" << this->oversetPatch_.name()
+                            << " set:" << nConstrained << " cells to:"
+                            << this->holeCellValue_ << endl;
+                    }
+                }
             }
         }
     }
+
+    coupledFvPatchField<Type>::initEvaluate(commsType);
+}
+
+
+template<class Type>
+void Foam::oversetFvPatchField<Type>::manipulateMatrix
+(
+    fvMatrix<Type>& matrix
+)
+{
+    if (this->manipulatedMatrix())
+    {
+        return;
+    }
+
+    const oversetFvPatch& ovp = this->oversetPatch_;
+
+    if (ovp.master())
+    {
+        if (fluxCorrection_ || (debug & 2))
+        {
+            storeFringeCoefficients(matrix);
+        }
+
+        // Trigger interpolation
+        const fvMesh& mesh = this->internalField().mesh();
+        const word& fldName = this->internalField().name();
+
+        // Try to find out if the solve routine comes from the unadapted mesh
+        // TBD. This should be cleaner.
+        if (&mesh.lduAddr() == &mesh.fvMesh::lduAddr())
+        {
+            return;
+        }
+
+        if (debug)
+        {
+            Pout<< FUNCTION_NAME << " field:" << fldName
+                << " patch:" << ovp.name() << endl;
+        }
+
+
+        // Calculate stabilised diagonal as normalisation for interpolation
+        const scalarField norm
+        (
+            dynamic_cast<const oversetFvMeshBase&>(mesh).normalisation(matrix)
+        );
+
+        const cellCellStencil& overlap = Stencil::New(mesh);
+        const labelUList& types = overlap.cellTypes();
+        const labelListList& stencil = overlap.cellStencil();
+
+        dynamic_cast<const oversetFvMeshBase&>(mesh).addInterpolation
+        (
+            matrix,
+            norm,
+            this->setHoleCellValue_,
+            this->holeCellValue_
+        );
+
+        if (debug)
+        {
+            pointField allCs(mesh.cellCentres());
+            const mapDistribute& map = overlap.cellInterpolationMap();
+            map.distribute(allCs, false, UPstream::msgType()+1);
+
+            // Make sure we don't interpolate from a hole
+
+            scalarField marker(this->primitiveField().size(), 0);
+            forAll(types, celli)
+            {
+                if (types[celli] == cellCellStencil::HOLE)
+                {
+                    marker[celli] = 1.0;
+                }
+            }
+            cellCellStencil::interpolate
+            (
+                marker,
+                mesh,
+                overlap,
+                overlap.cellInterpolationWeights()
+            );
+
+            forAll(marker, celli)
+            {
+                if
+                (
+                    types[celli] == cellCellStencil::INTERPOLATED
+                 && marker[celli] > SMALL
+                )
+                {
+                    WarningInFunction
+                        << " field:" << fldName
+                        << " patch:" << ovp.name()
+                        << " found:" << celli
+                        << " at:" << mesh.cellCentres()[celli]
+                        << " donorSlots:" << stencil[celli]
+                        << " at:"
+                        << UIndirectList<point>(allCs, stencil[celli])
+                        << " amount-of-hole:" << marker[celli]
+                        << endl;
+                }
+            }
+
+            // Make sure we don't have matrix coefficients for interpolated
+            // or hole cells
+
+            const lduAddressing& addr = mesh.lduAddr();
+            const labelUList& upperAddr = addr.upperAddr();
+            const labelUList& lowerAddr = addr.lowerAddr();
+            const scalarField& lower = matrix.lower();
+            const scalarField& upper = matrix.upper();
+
+            forAll(lowerAddr, facei)
+            {
+                const label l = lowerAddr[facei];
+                const bool lHole = (types[l] == cellCellStencil::HOLE);
+                const label u = upperAddr[facei];
+                const bool uHole = (types[u] == cellCellStencil::HOLE);
+
+                if
+                (
+                    (lHole && upper[facei] != 0.0)
+                 || (uHole && lower[facei] != 0.0)
+                )
+                {
+                    FatalErrorInFunction
+                        << "Hole-neighbouring face:" << facei
+                        << " lower:" << l
+                        << " type:" << types[l]
+                        << " coeff:" << lower[facei]
+                        << " upper:" << upperAddr[facei]
+                        << " type:" << types[u]
+                        << " coeff:" << upper[facei]
+                        << exit(FatalError);
+                }
+
+
+                // Empty donor list: treat like hole but still allow to
+                // influence neighbouring domains
+                const bool lEmpty =
+                (
+                    types[l] == cellCellStencil::INTERPOLATED
+                 && stencil[l].empty()
+                );
+                const bool uEmpty =
+                (
+                    types[u] == cellCellStencil::INTERPOLATED
+                 && stencil[u].empty()
+                );
+
+                if
+                (
+                    (lEmpty && upper[facei] != 0.0)
+                 || (uEmpty && lower[facei] != 0.0)
+                )
+                {
+                    FatalErrorInFunction
+                        << "Still connected face:" << facei << " lower:" << l
+                        << " type:" << types[l]
+                        << " coeff:" << lower[facei]
+                        << " upper:" << u
+                        << " type:" << types[u]
+                        << " coeff:" << upper[facei]
+                        << exit(FatalError);
+                }
+            }
+
+            forAll(matrix.internalCoeffs(), patchi)
+            {
+                const labelUList& fc = addr.patchAddr(patchi);
+                const Field<Type>& bouCoeffs = matrix.boundaryCoeffs()[patchi];
+
+                forAll(fc, i)
+                {
+                    const label celli = fc[i];
+                    const bool lHole = (types[celli] == cellCellStencil::HOLE);
+
+                    if (lHole && bouCoeffs[i] != pTraits<Type>::zero)
+                    {
+                        FatalErrorInFunction
+                            << "Patch:" << patchi
+                            << " patchFace:" << i
+                            << " lower:" << celli
+                            << " type:" << types[celli]
+                            << " bouCoeff:" << bouCoeffs[i]
+                            << exit(FatalError);
+                    }
+
+                    // Check whether this is influenced by neighbouring domains
+                    const bool lEmpty =
+                    (
+                        types[celli] == cellCellStencil::INTERPOLATED
+                     && stencil[celli].empty()
+                    );
+
+                    if (lEmpty && bouCoeffs[i] != pTraits<Type>::zero)
+                    {
+                        FatalErrorInFunction
+                            << "Patch:" << patchi
+                            << " patchFace:" << i
+                            << " lower:" << celli
+                            << " type:" << types[celli]
+                            << " bouCoeff:" << bouCoeffs[i]
+                            << exit(FatalError);
+                    }
+                }
+            }
+
+
+            // Make sure that diagonal is non-zero. Note: should add
+            // boundaryCoeff ...
+            const FieldField<Field, Type>& internalCoeffs =
+                matrix.internalCoeffs();
+
+            for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; ++cmpt)
+            {
+                // Replacement for m.addBoundaryDiag(norm, cmpt);
+                scalarField diag(matrix.diag());
+                forAll(internalCoeffs, patchi)
+                {
+                    const labelUList& fc = addr.patchAddr(patchi);
+                    const Field<Type>& intCoeffs = internalCoeffs[patchi];
+                    const scalarField cmptCoeffs(intCoeffs.component(cmpt));
+                    forAll(fc, i)
+                    {
+                        diag[fc[i]] += cmptCoeffs[i];
+                    }
+                }
+
+                forAll(diag, celli)
+                {
+                    if (mag(diag[celli]) < SMALL)
+                    {
+                        FatalErrorInFunction
+                            << "Patch:" << ovp.name()
+                            << " cell:" << celli
+                            << " at:" << mesh.cellCentres()[celli]
+                            << " diag:" << diag[celli]
+                            << exit(FatalError);
+                    }
+                }
+            }
+        }
+    }
+
+     fvPatchField<Type>::manipulateMatrix(matrix);
+}
+
+
+template<class Type>
+void Foam::oversetFvPatchField<Type>::initInterfaceMatrixUpdate
+(
+    solveScalarField& result,
+    const bool add,
+    const lduAddressing& lduAddr,
+    const label interfacei,
+    const solveScalarField& psiInternal,
+    const scalarField& coeffs,
+    const direction cmpt,
+    const Pstream::commsTypes commsType
+) const
+{
+}
+
+
+template<class Type>
+void Foam::oversetFvPatchField<Type>::updateInterfaceMatrix
+(
+    solveScalarField& result,
+    const bool add,
+    const lduAddressing& lduAddr,
+    const label patchId,
+    const solveScalarField& psiInternal,
+    const scalarField& coeffs,
+    const direction,
+    const Pstream::commsTypes commsType
+) const
+{
+    auto& psi = const_cast<solveScalarField&>(psiInternal);
+
+    if (fluxCorrection_ && this->oversetPatch_.master())
+    {
+        adjustPsi(psi, lduAddr, result);
+    }
 }
 
 
 template<class Type>
 void Foam::oversetFvPatchField<Type>::write(Ostream& os) const
 {
-    zeroGradientFvPatchField<Type>::write(os);
+    coupledFvPatchField<Type>::write(os);
-    // Make sure to write the value for ease of postprocessing.
+
-    this->writeEntry("value", os);
+    if (this->setHoleCellValue_)
+    {
+        os.writeEntry("setHoleCellValue", setHoleCellValue_);
+        os.writeEntry("holeCellValue", holeCellValue_);
+        os.writeEntryIfDifferent
+        (
+            "interpolateHoleCellValue",
+            false,
+            interpolateHoleCellValue_
+        );
+    }
+    os.writeEntryIfDifferent
+    (
+        "fluxCorrection",
+        false,
+        fluxCorrection_
+    );
+    os.writeEntryIfDifferent
+    (
+        "zone",
+        label(-1),
+        zoneId_
+    );
 }
 
 

src/overset/oversetPolyPatch/oversetFvPatchField.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2017 OpenCFD Ltd.
+    Copyright (C) 2016-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -42,7 +42,8 @@ SourceFiles
 #define oversetFvPatchField_H
 
 #include "oversetFvPatch.H"
-#include "zeroGradientFvPatchField.H"
+#include "coupledFvPatchFields.H"
+#include "oversetLduInterfaceField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -56,11 +57,12 @@ namespace Foam
 template<class Type>
 class oversetFvPatchField
 :
-    public zeroGradientFvPatchField<Type>
+    public oversetLduInterfaceField,
+    public coupledFvPatchField<Type>
 {
 protected:
 
-    // Protected data
+    // Protected Data
 
         //- Local reference cast into the overset patch
         const oversetFvPatch& oversetPatch_;
@@ -68,6 +70,32 @@ protected:
         //- Master patch ID
         mutable label masterPatchID_;
 
+        // Hole cell controls
+
+            //- Flag to set hole cell values
+            const bool setHoleCellValue_;
+
+            //- Flag to correct fluxes
+            const bool fluxCorrection_;
+
+            //- Flag to interpolate hole cell values (from nearby non-hole cell)
+            const bool interpolateHoleCellValue_;
+
+            //- Hole cell value
+            const Type holeCellValue_;
+
+            //- Fringe upper coefficients
+            mutable scalarField fringeUpperCoeffs_;
+
+            //- Fringe lower coefficients
+            mutable scalarField fringeLowerCoeffs_;
+
+            //- Fringe faces
+            mutable labelField fringeFaces_;
+
+            //- Zone to sum flux for mass conservation
+            label zoneId_;
+
 
 public:
 
@@ -133,13 +161,147 @@ public:
         }
 
 
-    // Member functions
+    // Member Functions
 
-        // Evaluation functions
+        // Coupled and adjust flux
+
+            //- Return neighbour field.  Dummy
+            virtual tmp<Field<Type> > patchNeighbourField() const;
+
+            //- Adjust psi for mass correction. Requires storeFringeCoefficients
+            //  to have been called before
+            void adjustPsi
+            (
+                solveScalarField& psi,
+                const lduAddressing& lduAddr,
+                solveScalarField& result
+            ) const;
+
+            //- Store fringe coefficients and faces
+            void storeFringeCoefficients(const fvMatrix<Type>& matrix);
+
+            //- Calculate patch flux (helper function). Requires
+            //  storeFringeCoefficients to have been called before
+            void fringeFlux
+            (
+                 const fvMatrix<Type>& m,
+                 const surfaceScalarField& phi
+            ) const;
+
+
+        // Evaluation
 
             //- Initialise the evaluation of the patch field
             virtual void initEvaluate(const Pstream::commsTypes commsType);
 
+            //- Return the matrix diagonal coefficients corresponding to the
+            //- evaluation of the value of this patchField with given weights
+            virtual tmp<Field<Type> > valueInternalCoeffs
+            (
+                const tmp<scalarField>&
+            ) const
+            {
+                return tmp<Field<Type>>
+                (
+                    new Field<Type>(this->size(), Zero)
+                );
+            }
+
+            //- Return the matrix source coefficients corresponding to the
+            //- evaluation of the value of this patchField with given weights
+            virtual tmp<Field<Type> > valueBoundaryCoeffs
+            (
+                const tmp<scalarField>&
+            ) const
+            {
+                return tmp<Field<Type>>
+                (
+                    new Field<Type>(this->size(), Zero)
+                );
+            }
+
+            //- Return the matrix diagonal coefficients corresponding to the
+            //- evaluation of the gradient of this patchField
+            tmp<Field<Type> > gradientInternalCoeffs() const
+            {
+                return tmp<Field<Type>>
+                (
+                    new Field<Type>(this->size(), Zero)
+                );
+            }
+
+            //- Return the matrix source coefficients corresponding to the
+            //- evaluation of the gradient of this patchField
+            tmp<Field<Type> > gradientBoundaryCoeffs() const
+            {
+                return tmp<Field<Type>>
+                (
+                    new Field<Type>(this->size(), Zero)
+                );
+            }
+
+            //- Manipulate matrix
+            virtual void manipulateMatrix(fvMatrix<Type>& matrix);
+
+
+        // Coupled interface functionality
+
+            //- Update result field based on interface functionality
+            virtual void updateInterfaceMatrix
+            (
+                solveScalarField& result,
+                const bool add,
+                const lduAddressing& lduAddr,
+                const label patchId,
+                const solveScalarField& psiInternal,
+                const scalarField& coeffs,
+                const direction,
+                const Pstream::commsTypes commsType
+            ) const;
+
+            //- Update result field based on interface functionality
+            virtual void updateInterfaceMatrix
+            (
+                Field<Type>&,
+                const bool add,
+                const lduAddressing& lduAddr,
+                const label patchId,
+                const Field<Type>&,
+                const scalarField&,
+                const Pstream::commsTypes commsType
+            ) const
+            {
+                NotImplemented;
+            }
+
+            //- Initialise neighbour matrix update
+            virtual void initInterfaceMatrixUpdate
+            (
+                Field<Type>&,
+                const bool add,
+                const lduAddressing&,
+                const label interfacei,
+                const Field<Type>&,
+                const scalarField&,
+                const Pstream::commsTypes commsType
+            ) const
+            {
+                NotImplemented;
+            }
+
+            virtual void initInterfaceMatrixUpdate
+            (
+                solveScalarField& result,
+                const bool add,
+                const lduAddressing&,
+                const label interfacei,
+                const solveScalarField& psiInternal,
+                const scalarField& coeffs,
+                const direction cmpt,
+                const Pstream::commsTypes commsType
+            ) const;
+
+
         // I-O
 
             //- Write

tutorials/basic/overLaplacianDyMFoam/heatTransferFluxConservation1D/0.orig/T

@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2206                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      T;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 1 0 0 0];
+
+internalField   uniform 273;
+
+boundaryField
+{
+    #includeEtc "caseDicts/setConstraintTypes"
+
+    "(walls|hole)"
+    {
+        type            zeroGradient;
+    }
+
+    free
+    {
+        type            overset;
+        fluxCorrection  true;
+        value           uniform 300;
+    }
+
+    left
+    {
+        type            fixedValue;
+        value           uniform 300;
+    }
+
+    right
+    {
+        type            fixedValue;
+        value           uniform 273;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/basic/overLaplacianDyMFoam/heatTransferFluxConservation1D/0.orig/cellDisplacement

@@ -0,0 +1,43 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2206                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    object      cellDisplacement;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 0 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    free
+    {
+        patchType       overset;
+        type            cellMotion;
+        value           uniform (0 0 0);
+    }
+
+    "(walls|hole|left1|right1)"
+    {
+        type            cellMotion;
+        value           uniform (0 0 0);
+    }
+
+    frontAndBack
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/basic/overLaplacianDyMFoam/heatTransferFluxConservation1D/0.orig/pointDisplacement

@@ -0,0 +1,43 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2206                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       pointVectorField;
+    object      pointDisplacement;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 0 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    ".*"
+    {
+        type            uniformFixedValue;
+        uniformValue    (0 0 0);
+    }
+
+    "(free|hole)"
+    {
+        patchType       overset;
+        type            uniformFixedValue;
+        uniformValue    table
+        (
+            (0.0         (0 0 0))
+            (1.0         (0.2 0 0))
+            (2.0         (0 0 0))
+        );
+    }
+}
+
+
+// ************************************************************************* //

tutorials/basic/overLaplacianDyMFoam/heatTransferFluxConservation1D/0.orig/zoneID

@@ -0,0 +1,32 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2206                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      zoneID;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    #includeEtc "caseDicts/setConstraintTypes"
+
+    ".*"
+    {
+        type        zeroGradient;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/basic/overLaplacianDyMFoam/heatTransferFluxConservation1D/Allclean

@@ -0,0 +1,8 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/CleanFunctions      # Tutorial clean functions
+#------------------------------------------------------------------------------
+
+cleanCase0
+
+#------------------------------------------------------------------------------

tutorials/basic/overLaplacianDyMFoam/heatTransferFluxConservation1D/Allrun

@@ -0,0 +1,10 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
+#------------------------------------------------------------------------------
+
+./Allrun.pre
+
+runApplication $(getApplication)
+
+#------------------------------------------------------------------------------

tutorials/basic/overLaplacianDyMFoam/heatTransferFluxConservation1D/Allrun.pre

@@ -0,0 +1,14 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
+#------------------------------------------------------------------------------
+
+runApplication blockMesh
+
+runApplication topoSet
+
+restore0Dir
+
+runApplication setFields
+
+#------------------------------------------------------------------------------

tutorials/basic/overLaplacianDyMFoam/heatTransferFluxConservation1D/README.md

@@ -0,0 +1,17 @@
+## Implicit flux conservation test
+
+The implicit flux conservation is switched on in the field (here T) overset BC.
+
+The keyword is used as:
+
+    free
+    {
+        type            overset;
+        fluxCorrection  true;
+        value           uniform 300;
+    }
+
+A special net fringe flux is output in this tutorial.
+
+Switch `fluxCorrection` to `false` to evaluate the flux correction.
+

tutorials/basic/overLaplacianDyMFoam/heatTransferFluxConservation1D/constant/dynamicMeshDict

@@ -0,0 +1,29 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2206                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      dynamicMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+motionSolverLibs (fvMotionSolvers overset);
+
+solver          displacementLaplacian;
+
+displacementLaplacianCoeffs
+{
+    diffusivity     uniform 1;
+}
+
+dynamicFvMesh       staticOversetFvMesh;
+
+
+// ************************************************************************* //

tutorials/basic/overLaplacianDyMFoam/heatTransferFluxConservation1D/constant/transportProperties

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

tutorials/basic/overLaplacianDyMFoam/heatTransferFluxConservation1D/system/blockMeshDict

@@ -0,0 +1,103 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2206                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+scale 1;
+
+vertices
+(
+    ( 0.00  0.0  0)
+    ( 1.00  0.0  0)
+    ( 1.00  1.0  0)
+    ( 0.00  1.0  0)
+    ( 0.00  0.0  1)
+    ( 1.00  0.0  1)
+    ( 1.00  1.0  1)
+    ( 0.00  1.0  1)
+
+//- 0 degrees
+    ( 0.5  0  0)
+    ( 1.5  0  0)
+    ( 1.5  1  0)
+    ( 0.5  1  0)
+    ( 0.5  0  1)
+    ( 1.5  0  1)
+    ( 1.5  1  1)
+    ( 0.5  1  1)
+);
+
+blocks
+(
+    hex (0 1 2 3 4 5 6 7) (10 10 1) simpleGrading (1 1 1)
+    hex (8 9 10 11 12 13 14 15) movingZone (7 7 1) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+    free
+    {
+        type    overset;
+        faces
+        (
+            (2 6 5 1)
+            (8 12 15 11)
+        );
+    }
+    walls
+    {
+        type    symmetry;
+        faces
+        (
+            (3 7 6 2)
+            (1 5 4 0)
+            (11 15 14 10)
+            ( 9 13 12  8)
+        );
+    }
+    frontAndBack
+    {
+        type    empty;
+        faces
+        (
+            (0 3 2 1)
+            (4 5 6 7)
+            ( 8 11 10  9)
+            (12 13 14 15)
+        );
+    }
+    left
+    {
+        type    patch;
+        faces
+        (
+            (0 4 7 3)
+        );
+    }
+    right
+    {
+        type    patch;
+        faces
+        (
+            (10 14 13 9)
+        );
+    }
+);
+
+
+// ************************************************************************* //

tutorials/basic/overLaplacianDyMFoam/heatTransferFluxConservation1D/system/controlDict

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

tutorials/basic/overLaplacianDyMFoam/heatTransferFluxConservation1D/system/fvSchemes

@@ -0,0 +1,81 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2206                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         steadyState;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+}
+
+divSchemes
+{
+    default         none;
+    div(phi,U)      bounded Gauss limitedLinearV 1;
+    div(phi,k)      bounded Gauss limitedLinear 1;
+    div(phi,epsilon) bounded Gauss limitedLinear 1;
+    div(phi,R)      bounded Gauss limitedLinear 1;
+    div(R)          Gauss linear;
+    div(phi,nuTilda) bounded Gauss limitedLinear 1;
+    div((nuEff*dev(T(grad(U))))) Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+fluxRequired
+{
+    default         no;
+    T;
+}
+
+oversetInterpolation
+{
+    // Interpolation scheme to use for overset calculation
+    method      inverseDistance;
+
+    // The inverseDistance method uses a 'voxel' like search structure.
+    // Optionally specify the extent and number of divisions n.
+    // Note that it will allocate an array of nx*ny*nz. If not specified:
+    // - searchBox          : local mesh bounding box
+    // - searchBoxDivisions : root (2D) or cube-root(3D) of number of cells
+    searchBox           (-0.05 -0.02 0)(1.52 1.05 1);
+    searchBoxDivisions  (120 120 1);
+    tolerance           1e-13;
+}
+
+oversetInterpolationSuppressed
+{
+    // Any additional fields that require overset interpolation
+}
+
+
+// ************************************************************************* //

tutorials/basic/overLaplacianDyMFoam/heatTransferFluxConservation1D/system/fvSolution

@@ -0,0 +1,66 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2206                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    T
+    {
+        tolerance       1e-10;
+        relTol          0;
+
+        solver          PBiCGStab;
+        preconditioner  DILU;
+
+        //// Note: GAMG not supported with fluxCorrection switched on
+        //solver          GAMG;
+        //smoother        GaussSeidel;
+        //nCellsInCoarsestLevel   1;
+        //// Important: overset adds faces so use matrix-based agglomerator
+        //agglomerator        algebraicPair;
+    }
+
+    cellDisplacement
+    {
+        solver          PCG;
+        preconditioner  DIC;
+        tolerance       1e-06;
+        relTol          0;
+        maxIter         100;
+    }
+}
+
+SIMPLE
+{
+    nNonOrthogonalCorrectors 0;
+    oversetPatchErrOutput    true;
+}
+
+relaxationFactors
+{
+    fields
+    {
+        p               0.3;
+    }
+    equations
+    {
+        U               0.7;
+        k               0.7;
+        omega           0.7;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/basic/overLaplacianDyMFoam/heatTransferFluxConservation1D/system/setFieldsDict

@@ -0,0 +1,46 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2206                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      setFieldsDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+defaultFieldValues
+(
+    volScalarFieldValue zoneID 123
+);
+
+regions
+(
+    cellToCell
+    {
+        set c0;
+
+        fieldValues
+        (
+            volScalarFieldValue zoneID 0
+        );
+    }
+
+    cellToCell
+    {
+        set c1;
+
+        fieldValues
+        (
+            volScalarFieldValue zoneID 1
+        );
+    }
+);
+
+
+// ************************************************************************* //

tutorials/basic/overLaplacianDyMFoam/heatTransferFluxConservation1D/system/topoSetDict

@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2206                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      topoSetDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+actions
+(
+    {
+        name    c0;
+        type    cellSet;
+        action  new;
+        source  regionToCell;
+        sourceInfo
+        {
+            insidePoints ((0.001 0.001 0.001));
+        }
+    }
+
+    {
+        name    c1;
+        type    cellSet;
+        action  new;
+        source  cellToCell;
+        sourceInfo
+        {
+            set c0;
+        }
+    }
+
+    {
+        name    c1;
+        type    cellSet;
+        action  invert;
+    }
+);
+
+
+// ************************************************************************* //

tutorials/compressible/overRhoPimpleDyMFoam/twoSimpleRotors/system/fvSolution

@@ -85,12 +85,9 @@ solvers
 PIMPLE
 {
     momentumPredictor   false;
-    correctPhi          true;
     nOuterCorrectors    1;
     nCorrectors         4;
     nNonOrthogonalCorrectors 0;
-    ddtCorr             false;
-
 }
 
 relaxationFactors

tutorials/heatTransfer/overBuoyantPimpleDyMFoam/movingBox/system/fvSolution

@@ -85,13 +85,10 @@ solvers
 PIMPLE
 {
     momentumPredictor   true;
-    correctPhi          false;
-    checkMeshCourantNo  yes;
     nOuterCorrectors    1;
     nCorrectors         4;
     nNonOrthogonalCorrectors 0;
-    ddtCorr             true;
+    pRefPoint           (0.0001 0.0001 0.001);
-    pRefCell            0;
     pRefValue           1e5;
 }
 

tutorials/incompressible/overPimpleDyMFoam/cylinder/cylinderAndBackground/system/fvSolution

@@ -73,13 +73,9 @@ solvers
 PIMPLE
 {
     momentumPredictor   false;
-    correctPhi          false;  //true;
-    oversetAdjustPhi    true;
     nOuterCorrectors    1;
     nCorrectors         2;
     nNonOrthogonalCorrectors 0;
-
-    ddtCorr             true;
 }
 
 relaxationFactors