TUT: Add and update overset tutorials

1) Deleting ddtCorr option
2) Deleting fluxInterpolation flux option
3) Reducing control readers files

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

@@ -149,7 +149,7 @@ int main(int argc, char *argv[])
 
     mesh.update();
 
-    surfaceScalarField faceMask(localMin<scalar>(mesh).interpolate(cellMask));
+    //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

@@ -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

@@ -86,9 +86,6 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readTimeControls.H"
-
-        #include "readControls.H"
         #include "readDyMControls.H"
 
         #include "compressibleCourantNo.H"
@@ -128,45 +125,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,16 +21,6 @@ surfaceScalarField phiHbyA
     fvc::flux(rho*HbyA) + phig
 );
 
-if (ddtCorr)
-{
-    surfaceScalarField faceMaskOld
-    (
-        localMin<scalar>(mesh).interpolate(cellMask.oldTime())
-    );
-
-    phiHbyA +=
-        faceMaskOld*MRF.zeroFilter(rhorAUf*fvc::ddtCorr(rho, U, phi));
-}
 
 MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
@@ -122,8 +112,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

@@ -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

@@ -46,11 +46,7 @@ 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"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -68,10 +64,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 +83,7 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readControls.H"
+        #include "readDyMControls.H"
         #include "CourantNo.H"
 
         #include "setDeltaT.H"
@@ -97,45 +92,26 @@ 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
+            if (correctPhi)
-            (
+            {
-                localMin<scalar>(mesh).interpolate(cellMask.oldTime())
+                // Corrects flux on separated regions
-            );
+                #include "correctPhi.H"
+            }
 
-            // Zero Uf on old faceMask (H-I)
+            fvc::makeRelative(phi, U);
-            Uf *= faceMaskOld;
-            // Update Uf and phi on new C-I faces
-            Uf += (1-faceMaskOld)*fvc::interpolate(U);
-            phi = mesh.Sf() & Uf;
 
-            // Zero phi on current H-I
+            if (checkMeshCourantNo)
-            surfaceScalarField faceMask
+            {
-            (
+                #include "meshCourantNo.H"
-                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,29 +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();
@@ -38,18 +20,14 @@ if (pimple.nCorrPISO() <= 1)
 
 phiHbyA = fvc::flux(HbyA);
 
-if (ddtCorr)
+if (runTime.outputTime())
 {
-    surfaceScalarField faceMaskOld
+    volScalarField divPhiHbyA("divPhiHbyA", fvc::div(phiHbyA));
-    (
+    divPhiHbyA.write();
-        localMin<scalar>(mesh).interpolate(cellMask.oldTime())
-    );
-    phiHbyA += rAUf*faceMaskOld*fvc::ddtCorr(U, Uf);
 }
 
 MRF.makeRelative(phiHbyA);
 
-// WIP
 if (p.needReference())
 {
     fvc::makeRelative(phiHbyA, U);
@@ -57,7 +35,8 @@ if (p.needReference())
     fvc::makeAbsolute(phiHbyA, U);
 }
 
-
+// WIP: To adjust phi on fringe faces to help mass
+// conservation
 if (adjustFringe)
 {
     fvc::makeRelative(phiHbyA, U);
@@ -79,27 +58,29 @@ 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);
     }
 }
 
-
+// Excludes error in interpolated/hole cells
 #include "continuityErrs.H"
 
-// Explicitly relax pressure for momentum corrector
+if (runTime.outputTime())
-p.relax();
+{
-volVectorField gradP(fvc::grad(p));
+    volVectorField Ucorrected("Ucorrected", U);
+    Ucorrected.write();
 
-// Option 2: zero out velocity on blocked out cells
+    volScalarField divPhiCon("divPhiCon", fvc::div(phi));
-//U = HbyA - rAU*cellMask*gradP;
+    divPhiCon.write();
-// 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 +90,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,5 +1,5 @@
 {
-    surfaceScalarField faceMask(localMin<scalar>(mesh).interpolate(cellMask));
+    //surfaceScalarField faceMask(localMin<scalar>(mesh).interpolate(cellMask));
 
     volScalarField rAU(1.0/UEqn.A());
     surfaceScalarField rAUf("rAUf", faceMask*fvc::interpolate(rAU));
@@ -7,12 +7,6 @@
     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,25 +151,7 @@ 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)
@@ -182,12 +161,6 @@ int main(int argc, char *argv[])
 
                     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 +175,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

@@ -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/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/overInterDyMFoam.C

@@ -49,13 +49,9 @@ 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"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -76,8 +72,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"
@@ -102,6 +97,7 @@ int main(int argc, char *argv[])
         #include "correctPhi.H"
     }
     #include "createUf.H"
+    #include "createControls.H"
 
     #include "setCellMask.H"
     #include "setInterpolatedCells.H"
@@ -119,7 +115,7 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readControls.H"
+        #include "readDyMControls.H"
 
         if (LTS)
         {
@@ -164,27 +160,7 @@ int main(int argc, char *argv[])
                     // Update cellMask field for blocking out hole cells
                     #include "setCellMask.H"
                     #include "setInterpolatedCells.H"
-
+                    #include "correctPhiFaceMask.H"
-                    const surfaceScalarField faceMaskOld
-                    (
-                        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)
@@ -194,17 +170,8 @@ int main(int argc, char *argv[])
 
                     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);
-
                 }
 
                 if (mesh.changing() && checkMeshCourantNo)
@@ -213,14 +180,9 @@ int main(int argc, char *argv[])
                 }
             }
 
-
             #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,8 +1,5 @@
 {
     rAU = 1.0/UEqn.A();
-    //mesh.interpolate(rAU);
-
-    surfaceScalarField faceMask(localMin<scalar>(mesh).interpolate(cellMask));
 
     surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
 
@@ -12,22 +9,8 @@
     //HbyA = rAU*UEqn.H();
     HbyA = constrainHbyA(rAU*H, U, p_rgh);
 
-    if (massFluxInterpolation)
-    {
-        #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())
@@ -37,13 +20,6 @@
         fvc::makeAbsolute(phiHbyA, U);
     }
 
-    if (adjustFringe)
-    {
-        fvc::makeRelative(phiHbyA, U);
-        oversetAdjustPhi(phiHbyA, U);
-        fvc::makeAbsolute(phiHbyA, U);
-    }
-
     surfaceScalarField phig
     (
         (

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/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,26 +151,7 @@ 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)
                     {
@@ -181,12 +160,6 @@ int main(int argc, char *argv[])
 
                     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);
                 }

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/Allwmake

@@ -88,6 +88,8 @@ wmakeLnInclude -u regionFaModels
 wmakeLnInclude -u faOptions
 regionModels/Allwmake $targetType $*
 
+wmakeLnInclude -u overset
+
 wmake $targetType fvOptions
 wmake $targetType faOptions
 wmake $targetType fvMotionSolver

src/fvOptions/Make/options

@@ -10,7 +10,8 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/solidThermo/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 \
@@ -20,4 +21,5 @@ LIB_LIBS = \
     -lsolidThermo \
     -lturbulenceModels \
     -lincompressibleTurbulenceModels \
-    -lcompressibleTurbulenceModels
+    -lcompressibleTurbulenceModels \
+    -loverset

src/fvOptions/cellSetOption/cellSetOption.C

@@ -28,6 +28,7 @@ License
 
 #include "cellSetOption.H"
 #include "volFields.H"
+#include "cellCellStencilObject.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -50,6 +51,7 @@ Foam::fv::cellSetOption::selectionModeTypeNames_
     { selectionModeType::smCellSet, "cellSet" },
     { selectionModeType::smCellZone, "cellZone" },
     { selectionModeType::smAll, "all" },
+    { selectionModeType::smCellType, "cellType" },
 });
 
 
@@ -78,6 +80,11 @@ void Foam::fv::cellSetOption::setSelection(const dictionary& dict)
         {
             break;
         }
+        case smCellType:
+        {
+            break;
+
+        }
         default:
         {
             FatalErrorInFunction
@@ -181,6 +188,21 @@ void Foam::fv::cellSetOption::setCellSet()
             cells_ = identity(mesh_.nCells());
             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
@@ -235,7 +257,7 @@ bool Foam::fv::cellSetOption::isActive()
                 // Force printing of new set volume
                 V_ = -GREAT;
             }
-            else if (selectionMode_ == smPoints)
+            else if ((selectionMode_ == smPoints) || (selectionMode_ == smCellType))
             {
                 // This is the only geometric selection mode
                 setCellSet();

src/fvOptions/cellSetOption/cellSetOption.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2017-2020 OpenCFD Ltd.
+    Copyright (C) 2017-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -129,7 +129,8 @@ public:
             smPoints,
             smCellSet,
             smCellZone,
-            smAll
+            smAll,
+            smCellType
         };
 
         //- 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 * * * * * * * * * * * * * //
 
@@ -64,14 +65,14 @@ void Foam::fv::velocityDampingConstraint::addDamping(fvMatrix<vector>& eqn)
 
     label nDamped = 0;
 
-    forAll(U, cellI)
+    for (const label celli : cells_)
     {
-        const scalar magU = mag(U[cellI]);
+        const scalar magU = mag(U[celli]);
         if (magU > UMax_)
         {
-            const scalar scale = sqr(Foam::cbrt(vol[cellI]));
+            const scalar scale = sqr(Foam::cbrt(vol[celli]));
 
-            diag[cellI] += scale*(magU-UMax_);
+            diag[celli] += C_*scale*(magU-UMax_);
 
             ++nDamped;
         }
@@ -96,7 +97,8 @@ Foam::fv::velocityDampingConstraint::velocityDampingConstraint
     const fvMesh& mesh
 )
 :
-    fv::cellSetOption(name, modelType, dict, mesh)
+    fv::cellSetOption(name, modelType, dict, mesh),
+    C_(1)
 {
     read(dict);
 }
@@ -126,6 +128,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-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -127,6 +127,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
 

src/overset/cellCellStencil/cellCellStencil/cellCellStencil.C

@@ -48,6 +48,7 @@ Foam::cellCellStencil::cellTypeNames_
     { cellType::CALCULATED, "calculated" },
     { cellType::INTERPOLATED, "interpolated" },
     { cellType::HOLE, "hole" },
+    { cellType::SPECIAL, "special" },
 });
 
 
@@ -96,6 +97,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"))
@@ -297,4 +309,88 @@ void Foam::cellCellStencil::globalCellCells
 }
 
 
+// * * * * * * * * * * * * * 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

@@ -74,7 +74,9 @@ public:
     {
         CALCULATED = 0,     // normal operation
         INTERPOLATED = 1,   // interpolated
-        HOLE = 2            // hole
+        HOLE = 2,           // hole
+        SPECIAL = 3,         // hole that changes
+        POROUS = 4
     };
 
 
@@ -94,8 +96,6 @@ protected:
 
     // Protected Member Functions
 
-        //- Count occurrences (in parallel)
-        static labelList count(const label size, const labelUList& lst);
 
         //- Helper: create volScalarField for postprocessing.
         template<class Type>
@@ -106,6 +106,9 @@ protected:
             const UList<Type>&
         );
 
+        //- Helper: strip off trailing _0
+        static word baseName(const word& name);
+
 
 private:
 
@@ -213,6 +216,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 +229,50 @@ 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;
+        }
 };
 
+template<>
+Ostream& operator<<(Ostream& os, const InfoProxy<cellCellStencil>& ic);
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/overset/cellCellStencil/cellCellStencil/cellCellStencilTemplates.C

@@ -29,6 +29,128 @@ 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)
+    {
+        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]];
+        }
+        //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;
+    }
+}
+
+
+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
@@ -66,4 +188,46 @@ Foam::cellCellStencil::createField
 }
 
 
+template<class GeoField, class SuppressBC>
+void Foam::cellCellStencil::correctBoundaryConditions
+(
+    GeoField& psi
+    //const labelHashSet& suppress
+)
+{
+    // Version of GeoField::correctBoundaryConditions that exclude evaluation
+    // of oversetFvPatchFields
+
+    typename GeoField::Boundary& bfld = psi.boundaryFieldRef();
+
+    const label nReq = Pstream::nRequests();
+
+    forAll(bfld, patchi)
+    {
+        //if (!suppress.found(patchi))
+        //if (!isA<oversetFvPatchField<GeoField::value_type>(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 (!suppress.found(patchi))
+        //if (!isA<oversetFvPatchField<GeoField::value_type>(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,11 +57,33 @@ Foam::cellCellStencils::cellVolumeWeight::defaultOverlapTolerance_ = 1e-6;
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+void Foam::cellCellStencils::cellVolumeWeight::seedCell
+(
+    const label cellI,
+    const scalar wantedFraction,
+    bitSet& isFront,
+    scalarField& fraction
+) const
+{
+    const cell& cFaces = mesh_.cells()[cellI];
+    forAll(cFaces, i)
+    {
+        label nbrFacei = cFaces[i];
+        if (fraction[nbrFacei] < wantedFraction)
+        {
+            fraction[nbrFacei] = wantedFraction;
+            isFront.set(nbrFacei);
+        }
+    }
+}
+
+
 void Foam::cellCellStencils::cellVolumeWeight::walkFront
 (
     const scalar layerRelax,
     labelList& allCellTypes,
-    scalarField& allWeight
+    scalarField& allWeight,
+    const labelList& zoneID
 ) const
 {
     // Current front
@@ -69,17 +92,20 @@ void Foam::cellCellStencils::cellVolumeWeight::walkFront
 
     const fvBoundaryMesh& fvm = mesh_.boundary();
 
-
     // 'overset' patches
-
     forAll(fvm, patchI)
     {
         if (isA<oversetFvPatch>(fvm[patchI]))
         {
+            const labelList& fc = fvm[patchI].faceCells();
             //Pout<< "Storing faces on patch " << fvm[patchI].name() << endl;
-            forAll(fvm[patchI], i)
+            forAll(fc, i)
             {
-                isFront.set(fvm[patchI].start()+i);
+                label cellI = fc[i];
+                if (allCellTypes[cellI] == INTERPOLATED)
+                {
+                    isFront.set(fvm[patchI].start()+i);
+                }
             }
         }
     }
@@ -96,8 +122,9 @@ void Foam::cellCellStencils::cellVolumeWeight::walkFront
             label neiType = allCellTypes[nei[faceI]];
             if
             (
-                 (ownType == HOLE && neiType != HOLE)
+                 ((ownType == HOLE && neiType != HOLE)
-              || (ownType != HOLE && neiType == HOLE)
+              || (ownType != HOLE && neiType == HOLE))
+              //&& (zoneID[own[faceI]] == 0 && zoneID[nei[faceI]] == 0)
             )
             {
                 //Pout<< "Front at face:" << faceI
@@ -121,8 +148,9 @@ void Foam::cellCellStencils::cellVolumeWeight::walkFront
 
             if
             (
-                 (ownType == HOLE && neiType != HOLE)
+                 ((ownType == HOLE && neiType != HOLE)
-              || (ownType != HOLE && neiType == HOLE)
+              || (ownType != HOLE && neiType == HOLE))
+             // && (zoneID[own[faceI]] == 0 && zoneID[nei[faceI]] == 0)
             )
             {
                 //Pout<< "Front at coupled face:" << faceI
@@ -505,6 +533,7 @@ void Foam::cellCellStencils::cellVolumeWeight::interpolatePatchTypes
             {
                 // 'patch' overrides -1 and 'other'
                 result[cellI] = PATCH;
+                break;
             }
             else if (result[cellI] == -1)
             {
@@ -591,27 +620,27 @@ 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 = false;
+                validDonors = true;
             }
             break;
         }
@@ -687,6 +716,10 @@ Foam::cellCellStencils::cellVolumeWeight::cellVolumeWeight
         mesh_,
         dimensionedScalar(dimless, Zero),
         zeroGradientFvPatchScalarField::typeName
+    ),
+    allowInterpolatedDonors_
+    (
+        dict.getOrDefault<bool>("allowInterpolatedDonors", true)
     )
 {
     // Protect local fields from interpolation
@@ -703,6 +736,9 @@ Foam::cellCellStencils::cellVolumeWeight::cellVolumeWeight
     // Read zoneID
     this->zoneID();
 
+    overlapTolerance_ =
+        dict_.getOrDefault("overlapTolerance", defaultOverlapTolerance_);
+
     // Read old-time cellTypes
     IOobject io
     (
@@ -758,7 +794,6 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
     Pstream::listCombineGather(nCellsPerZone, plusEqOp<label>());
     Pstream::listCombineScatter(nCellsPerZone);
 
-
     Info<< typeName << " : detected " << nZones
         << " mesh regions" << nl << endl;
 
@@ -780,7 +815,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);
@@ -805,6 +839,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
@@ -818,7 +861,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();
@@ -840,7 +882,6 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
                 false                   // do not normalise
             );
 
-
             {
                 // Get tgt patch types on src mesh
                 labelList interpolatedTgtPatchTypes(srcMesh.nCells(), -1);
@@ -866,6 +907,8 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
                         mapper.tgtMap()->distribute(tgtGlobalCells);
                     }
                 }
+
+
                 combineCellTypes
                 (
                     srcI,
@@ -933,13 +976,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();
+
     }
 
 
@@ -961,8 +1011,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();
@@ -973,26 +1021,32 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
         }
     }
 
-/*
+
     // Knock out cell with insufficient interpolation weights
-    forAll(allCellTypes, cellI)
+//     forAll(allCellTypes, cellI)
-    {
+//     {
-        if (allCellTypes[cellI] == INTERPOLATED)
+//         if (allCellTypes[cellI] == INTERPOLATED)
-        {
+//         {
-            scalar v = mesh_.V()[cellI];
+//             scalar v = mesh_.V()[cellI];
-            scalar overlapVol = sum(allWeights[cellI]);
+//             scalar overlapVol = sum(allWeights[cellI]);
-            if (overlapVol < (1.0-overlapTolerance_)*v)
+//             if (overlapVol < (1.0-overlapTolerance_)*v)
-            {
+//             {
-                //Pout<< "Holeing cell:" << cellI
+//                 //Pout<< "Holeing cell:" << cellI
-                //    << " at:" << mesh_.cellCentres()[cellI] << endl;
+//                 //    << " at:" << mesh_.cellCentres()[cellI] << endl;
-                allCellTypes[cellI] = HOLE;
+//                 allCellTypes[cellI] = HOLE;
-                allWeights[cellI].clear();
+//                 allWeights[cellI].clear();
-                allStencil[cellI].clear();
+//                 allStencil[cellI].clear();
-            }
+//             }
-        }
+//         }
-    }
+//         // Knocks out HOLE stencil
-*/
+//         if (allCellTypes[cellI] == HOLE || allCellTypes[cellI] == CALCULATED)
-    if (debug)
+//         {
+//             allWeights[cellI].clear();
+//             allStencil[cellI].clear();
+//         }
+//     }
+
+    if ((debug&2) && (mesh_.time().outputTime()))
     {
         tmp<volScalarField> tfld
         (
@@ -1000,6 +1054,42 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
         );
         //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.ref().correctBoundaryConditions();
+        tfld().write();
+
+        labelList stencilSize(mesh_.nCells());
+        forAll(allStencil, celli)
+        {
+            stencilSize[celli] = allStencil[celli].size();
+        }
+        tmp<volScalarField> tfldStencil
+        (
+            createField(mesh_, "allStencil_hole", stencilSize)
+        );
+        tfldStencil().write();
     }
 
 
@@ -1038,25 +1128,53 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
     }
 
 
-    // Mark unreachable bits
-    findHoles(globalCells, mesh_, zoneID, allStencil, allCellTypes);
-
-    if (debug)
-    {
-        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);
+    walkFront(layerRelax, allCellTypes, allWeight, zoneID);
 
-    if (debug)
+    // Knock out cell with insufficient interpolation weights
+//      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();
+//             }
+//         }
+//         // Knocks out HOLE stencil
+//         if (allCellTypes[cellI] == HOLE || allCellTypes[cellI] == CALCULATED)
+//         {
+//             allWeights[cellI].clear();
+//             allStencil[cellI].clear();
+//         }
+//     }
+
+    // Try to make calculated cells that are nrb with patches
+    //
+//     forAll(allPatchTypes, cellI)
+//     {
+//         if (allPatchTypes[cellI] == patchCellType::PATCH)
+//         {
+//             if (allCellTypes[cellI] == INTERPOLATED)
+//             {
+//                 allCellTypes[cellI] = CALCULATED;
+//                 allWeights[cellI].clear();
+//                 allStencil[cellI].clear();
+//                 allWeight[cellI] = 0.0;
+//             }
+//         }
+//     }
+
+
+
+    if ((debug&2) && (mesh_.time().outputTime()))
     {
         tmp<volScalarField> tfld
         (
@@ -1064,27 +1182,99 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
         );
         //tfld.ref().correctBoundaryConditions();
         tfld().write();
+
+        labelList stencilSize(mesh_.nCells());
+        forAll(allStencil, celli)
+        {
+            stencilSize[celli] = allStencil[celli].size();
+        }
+        tmp<volScalarField> tfldStencil
+        (
+            createField(mesh_, "allStencil_front", stencilSize)
+        );
+        tfldStencil().write();
     }
 
     // Normalise weights, Clear storage
+//     forAll(allCellTypes, cellI)
+//     {
+//         if (allCellTypes[cellI] == INTERPOLATED)
+//         {
+//             if (allWeight[cellI] < SMALL || allStencil[cellI].size() == 0)
+//             {
+//                 //Pout<< "Clearing cell:" << cellI
+//                 //    << " at:" << mesh_.cellCentres()[cellI] << endl;
+//                 allWeights[cellI].clear();
+//                 allStencil[cellI].clear();
+//                 allWeight[cellI] = 0.0;
+//             }
+//             else
+//             {
+//                 scalar s = sum(allWeights[cellI]);
+//                 forAll(allWeights[cellI], i)
+//                 {
+//                     allWeights[cellI][i] /= s;
+//                 }
+//             }
+//         }
+//         else
+//         {
+//             allWeights[cellI].clear();
+//             allStencil[cellI].clear();
+//         }
+//     }
+
+    // Eliminate weights with HOLE donors
+    // Clear storage
+//     forAll(allCellTypes, cellI)
+//     {
+//         if (allCellTypes[cellI] == INTERPOLATED)
+//         {
+//             forAll (allStencil[cellI], subCellI)
+//             {
+//                 label subCell = allStencil[cellI][subCellI];
+//
+//                 if (allCellTypes[subCell] == HOLE)
+//                 {
+//                     allWeights[cellI][subCellI] = 0;
+//                 }
+//             }
+//         }
+//         else
+//         {
+//             allWeights[cellI].clear();
+//             allStencil[cellI].clear();
+//         }
+//     }
+
+
+    labelListList compactStencil(allStencil);
+    List<Map<label>> compactStencilMap;
+    mapDistribute map(globalCells, compactStencil, compactStencilMap);
+
+    labelList compactCellTypes(allCellTypes);
+    map.distribute(compactCellTypes);
+
+    label nHoleDonors = 0;
     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++;
                 }
             }
         }
@@ -1095,10 +1285,60 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
         }
     }
 
+    reduce(nHoleDonors, sumOp<label>());
+
+    Pout << "nHole Donors : " << nHoleDonors << endl;
+
+    // Normalize weights
+    forAll(allCellTypes, cellI)
+    {
+        if (allCellTypes[cellI] == INTERPOLATED)
+        {
+            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()))
     {
+         if ((debug&2) && (mesh_.time().outputTime()))
+        {
+            tmp<volScalarField> tfld
+            (
+                createField(mesh_, "allCellTypes_final", allCellTypes)
+            );
+            //tfld.ref().correctBoundaryConditions();
+            tfld().write();
+        }
+
+
+        labelList stencilSize(mesh_.nCells());
+        forAll(allStencil, celli)
+        {
+            stencilSize[celli] = allStencil[celli].size();
+        }
+        tmp<volScalarField> tfldStencil
+        (
+            createField(mesh_, "allStencilFinal", stencilSize)
+        );
+        tfldStencil().write();
+
         volScalarField patchTypes
         (
             IOobject
@@ -1120,14 +1360,14 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
             patchTypes[cellI] = allPatchTypes[cellI];
         }
         //patchTypes.correctBoundaryConditions();
-        dynamicOversetFvMesh::correctBoundaryConditions
+        oversetFvMeshBase::correctBoundaryConditions
         <
             volScalarField,
             oversetFvPatchField<scalar>
         >(patchTypes.boundaryFieldRef(), false);
         patchTypes.write();
     }
-    if (debug)
+    if ((debug&2) && (mesh_.time().outputTime()))
     {
         volScalarField volTypes
         (
@@ -1150,7 +1390,7 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
             volTypes[cellI] = allCellTypes[cellI];
         }
         //volTypes.correctBoundaryConditions();
-        dynamicOversetFvMesh::correctBoundaryConditions
+        oversetFvMeshBase::correctBoundaryConditions
         <
             volScalarField,
             oversetFvPatchField<scalar>
@@ -1159,47 +1399,12 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
     }
 
 
-//  // 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>
@@ -1223,7 +1428,7 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
     );
 
     // Dump interpolation stencil
-    if (debug)
+    if ((debug&2) && (mesh_.time().outputTime()))
     {
         // Dump weight
         cellInterpolationWeight_.instance() = mesh_.time().timeName();
@@ -1250,10 +1455,12 @@ 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.write(linePointRef(accCc, 0.1*accCc+0.9*donorCc));
+                    const point& accCc = mesh_.cellCentres()[cellI];
+                    str.write(linePointRef(accCc, 0.1*accCc+0.9*donorCc));
+                }
             }
         }
     }

src/overset/cellCellStencil/cellVolumeWeight/cellVolumeWeightCellCellStencil.H

@@ -90,14 +90,35 @@ 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 labelList& zoneID
+//         ) const;
+
+        //- Seed faces of cell with wantedFraction (if higher than current)
+        void seedCell
+        (
+            const label cellI,
+            const scalar wantedFraction,
+            bitSet& isFront,
+            scalarField& fraction
+        ) const;
+
+
         void walkFront
         (
             const scalar layerRelax,
             labelList& allCellTypes,
-            scalarField& allWeight
+            scalarField& allWeight,
+            const labelList& zoneID
         ) const;
 
         //- Find cells next to cells of type PATCH

src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.C

@@ -39,9 +39,12 @@ License
 #include "waveMethod.H"
 
 #include "regionSplit.H"
-#include "dynamicOversetFvMesh.H"
+#include "oversetFvPatchFields.H"
-//#include "minData.H"
+#include "topoDistanceData.H"
-//#include "FaceCellWave.H"
+#include "FaceCellWave.H"
+
+#include "OBJstream.H"
+#include "uindirectPrimitivePatch.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -510,7 +513,12 @@ void Foam::cellCellStencils::inverseDistance::markDonors
         forAll(tgtCellMap, tgtCelli)
         {
             label srcCelli = tgtToSrcAddr[tgtCelli];
-            if (srcCelli != -1 && allCellTypes[srcCellMap[srcCelli]] != HOLE)
+            if
+            (
+                srcCelli != -1
+             //&& allCellTypes[srcCellMap[srcCelli]] != HOLE
+             && allCellTypes[tgtCellMap[tgtCelli]] != HOLE
+            )
             {
                 label celli = tgtCellMap[tgtCelli];
 
@@ -553,7 +561,6 @@ void Foam::cellCellStencils::inverseDistance::markDonors
     }
 
 
-
     // Indices of tgtcells to send over to each processor
     List<DynamicList<label>> tgtSendCells(Pstream::nProcs());
     forAll(srcOverlapProcs, i)
@@ -577,7 +584,10 @@ void Foam::cellCellStencils::inverseDistance::markDonors
                 label procI = srcOverlapProcs[i];
                 if (subBb.overlaps(srcBbs[procI]))
                 {
-                    tgtSendCells[procI].append(tgtCelli);
+                    //if (allCellTypes[tgtCellMap[tgtCelli]] != HOLE)
+                    {
+                        tgtSendCells[procI].append(tgtCelli);
+                    }
                 }
             }
         }
@@ -610,7 +620,7 @@ void Foam::cellCellStencils::inverseDistance::markDonors
         {
             const point& sample = samples[sampleI];
             label srcCelli = srcMesh.findCell(sample, polyMesh::CELL_TETS);
-            if (srcCelli != -1 && allCellTypes[srcCellMap[srcCelli]] != HOLE)
+            if ((srcCelli != -1))// && (allCellTypes[srcCellMap[srcCelli]] != HOLE))
             {
                 donors[sampleI] = globalCells.toGlobal(srcCellMap[srcCelli]);
             }
@@ -618,7 +628,7 @@ void Foam::cellCellStencils::inverseDistance::markDonors
 
         // Use same pStreamBuffers to send back.
         UOPstream os(procI, pBufs);
-        os << donors;
+        os << donors;// << srcCellis;
     }
     pBufs.finishedSends();
 
@@ -643,18 +653,20 @@ void Foam::cellCellStencils::inverseDistance::markDonors
             if (globalDonor != -1)
             {
                 label celli = tgtCellMap[cellIDs[donorI]];
-
-                // TBD: check for multiple donors. Maybe better one? For
-                //      now check 'nearer' mesh
-                if (betterDonor(tgtI, allDonor[celli], srcI))
                 {
-                    allStencil[celli].setSize(1);
+                    // TBD: check for multiple donors. Maybe better one? For
-                    allStencil[celli][0] = globalDonor;
+                    //      now check 'nearer' mesh
-                    allDonor[celli] = srcI;
+                    if (betterDonor(tgtI, allDonor[celli], srcI))
+                    {
+                        allStencil[celli].setSize(1);
+                        allStencil[celli][0] = globalDonor;
+                        allDonor[celli] = srcI;
+                    }
                 }
             }
         }
     }
+
 }
 
 
@@ -1211,15 +1223,18 @@ void Foam::cellCellStencils::inverseDistance::findHoles
 
 
     // See which regions have not been visited (regionType == 1)
+    label count = 0;
     forAll(cellRegion, cellI)
     {
         label type = regionType[cellRegion[cellI]];
         if (type == 1 && cellTypes[cellI] != HOLE)
         {
             cellTypes[cellI] = HOLE;
+            count++;
         }
     }
-    DebugInfo<< FUNCTION_NAME << " : Finished hole flood filling" << endl;
+    DebugInfo<< FUNCTION_NAME << " : Finished hole flood filling : " << count
+        << endl;
 }
 
 
@@ -1249,7 +1264,8 @@ void Foam::cellCellStencils::inverseDistance::walkFront
     const scalar layerRelax,
     const labelListList& allStencil,
     labelList& allCellTypes,
-    scalarField& allWeight
+    scalarField& allWeight,
+    const labelList& zoneID
 ) const
 {
     // Current front
@@ -1259,7 +1275,6 @@ void Foam::cellCellStencils::inverseDistance::walkFront
 
 
     // 'overset' patches
-
     forAll(fvm, patchI)
     {
         if (isA<oversetFvPatch>(fvm[patchI]))
@@ -1291,8 +1306,9 @@ void Foam::cellCellStencils::inverseDistance::walkFront
             label neiType = allCellTypes[nei[faceI]];
             if
             (
-                 (ownType == HOLE && neiType != HOLE)
+                 ((ownType == HOLE && neiType != HOLE)
-              || (ownType != HOLE && neiType == HOLE)
+              || (ownType != HOLE && neiType == HOLE))
+             // && (zoneID[own[faceI]] == 0 && zoneID[nei[faceI]] == 0)
             )
             {
                 //Pout<< "Front at face:" << faceI
@@ -1316,8 +1332,9 @@ void Foam::cellCellStencils::inverseDistance::walkFront
 
             if
             (
-                 (ownType == HOLE && neiType != HOLE)
+                 ((ownType == HOLE && neiType != HOLE)
-              || (ownType != HOLE && neiType == HOLE)
+              || (ownType != HOLE && neiType == HOLE))
+             // && (zoneID[own[faceI]] == 0 && zoneID[nei[faceI]] == 0)
             )
             {
                 //Pout<< "Front at coupled face:" << faceI
@@ -1466,10 +1483,301 @@ void Foam::cellCellStencils::inverseDistance::stencilWeights
 }
 
 
-void Foam::cellCellStencils::inverseDistance::createStencil
+void Foam::cellCellStencils::inverseDistance::holeExtrapolationStencil
 (
     const globalIndex& globalCells
 )
+{
+    // Detects holes that are used for interpolation. If so
+    // - make type interpolated
+    // - add interpolation to nearest non-hole cell(s)
+
+    const labelList& owner = mesh_.faceOwner();
+    const labelList& neighbour = mesh_.faceNeighbour();
+
+
+    // Get hole-cells that are used in an interpolation stencil
+    boolList isDonor(cellInterpolationMap().constructSize(), false);
+    label nHoleDonors = 0;
+    {
+        for (const label celli : interpolationCells_)
+        {
+            const labelList& slots = cellStencil_[celli];
+            UIndirectList<bool>(isDonor, slots) = true;
+        }
+
+        mapDistributeBase::distribute<bool, orEqOp<bool>, flipOp>
+        (
+            Pstream::commsTypes::nonBlocking,
+            List<labelPair>(),
+            mesh_.nCells(),
+            cellInterpolationMap().constructMap(),
+            false,
+            cellInterpolationMap().subMap(),
+            false,
+            isDonor,
+            false,
+            orEqOp<bool>(),
+            flipOp()                               // negateOp
+            //false                                   // nullValue
+        );
+
+
+        nHoleDonors = 0;
+        forAll(cellTypes_, celli)
+        {
+            if (cellTypes_[celli] == cellCellStencil::HOLE && isDonor[celli])
+            {
+                nHoleDonors++;
+            }
+        }
+        reduce(nHoleDonors, sumOp<label>());
+        if (debug)
+        {
+            Pout<< "Detected " << nHoleDonors
+                << " hole cells that are used as donors" << endl;
+        }
+    }
+
+
+    if (nHoleDonors)
+    {
+        // Convert map and stencil back to global indices by 'interpolating'
+        // global cells. Note: since we're only adding interpolations
+        // (HOLE->SPECIAL) this could probably be done more efficiently.
+
+        labelList globalCellIDs(identity(mesh_.nCells()));
+        globalCells.inplaceToGlobal(Pstream::myProcNo(), globalCellIDs);
+        cellInterpolationMap().distribute(globalCellIDs);
+
+        forAll(interpolationCells_, i)
+        {
+            label cellI = interpolationCells_[i];
+            const labelList& slots = cellStencil_[cellI];
+            cellStencil_[cellI] = UIndirectList<label>(globalCellIDs, slots)();
+        }
+
+
+        labelList nbrTypes;
+        syncTools::swapBoundaryCellList(mesh_, cellTypes_, nbrTypes);
+
+        label nSpecialNear = 0;
+        label nSpecialFar = 0;
+
+
+        // 1. Find hole cells next to live cells
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        {
+            List<pointList> donorCcs(mesh_.nCells());
+            labelListList donorCells(mesh_.nCells());
+
+            for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
+            {
+                label own = owner[facei];
+                bool ownHole = (cellTypes_[own] == cellCellStencil::HOLE);
+                label nbr = neighbour[facei];
+                bool nbrHole = (cellTypes_[nbr] == cellCellStencil::HOLE);
+                if (isDonor[own] && ownHole && !nbrHole)
+                {
+                    donorCells[own].append(globalCells.toGlobal(nbr));
+                    donorCcs[own].append(mesh_.cellCentres()[nbr]);
+                }
+                else if (isDonor[nbr] && nbrHole && !ownHole)
+                {
+                    donorCells[nbr].append(globalCells.toGlobal(own));
+                    donorCcs[nbr].append(mesh_.cellCentres()[own]);
+                }
+            }
+            labelList globalCellIDs(identity(mesh_.nCells()));
+            globalCells.inplaceToGlobal(Pstream::myProcNo(), globalCellIDs);
+            labelList nbrCells;
+            syncTools::swapBoundaryCellList(mesh_, globalCellIDs, nbrCells);
+            pointField nbrCc;
+            syncTools::swapBoundaryCellList(mesh_, mesh_.cellCentres(), nbrCc);
+            forAll(nbrTypes, bFacei)
+            {
+                label facei = bFacei+mesh_.nInternalFaces();
+                label own = owner[facei];
+                bool ownHole = (cellTypes_[own] == cellCellStencil::HOLE);
+                bool nbrHole = (nbrTypes[bFacei] == cellCellStencil::HOLE);
+                if (isDonor[own] && ownHole && !nbrHole)
+                {
+                    donorCells[own].append(nbrCells[bFacei]);
+                    donorCcs[own].append(nbrCc[bFacei]);
+                }
+            }
+
+            forAll(donorCcs, celli)
+            {
+                if (donorCcs[celli].size())
+                {
+                    cellStencil_[celli] = std::move(donorCells[celli]);
+                    stencilWeights
+                    (
+                        mesh_.cellCentres()[celli],
+                        donorCcs[celli],
+                        cellInterpolationWeights_[celli]
+                    );
+                    cellTypes_[celli] = SPECIAL;
+                    interpolationCells_.append(celli);
+                    cellInterpolationWeight_[celli] = 1.0;
+                    nSpecialNear++;
+                }
+            }
+        }
+
+
+
+
+        // 2. Walk to find (topologically) nearest 'live' cell
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        // This will do the remainder but will only find a single donor
+
+        // Field on cells and faces.
+        List<topoDistanceData<int>> cellData(mesh_.nCells());
+        List<topoDistanceData<int>> faceData(mesh_.nFaces());
+
+        int dummyTrackData(0);
+
+        {
+            // Mark all non-hole cells to disable any wave across them
+            forAll(cellTypes_, celli)
+            {
+                label cType = cellTypes_[celli];
+                if
+                (
+                    cType == cellCellStencil::CALCULATED
+                 || cType == cellCellStencil::INTERPOLATED
+                )
+                {
+                    cellData[celli] = topoDistanceData<int>(labelMin, 0);
+                }
+            }
+
+            DynamicList<label> seedFaces(nHoleDonors);
+            DynamicList<topoDistanceData<int>> seedData(nHoleDonors);
+
+            for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
+            {
+                label own = owner[facei];
+                bool ownLive =
+                (
+                    cellTypes_[own] == cellCellStencil::CALCULATED
+                 || cellTypes_[own] == cellCellStencil::INTERPOLATED
+                );
+
+                label nbr = neighbour[facei];
+                bool nbrLive =
+                (
+                    cellTypes_[nbr] == cellCellStencil::CALCULATED
+                 || cellTypes_[nbr] == cellCellStencil::INTERPOLATED
+                );
+
+                if (ownLive && !nbrLive)
+                {
+                    seedFaces.append(facei);
+                    const label globalOwn = globalCells.toGlobal(own);
+                    seedData.append(topoDistanceData<int>(globalOwn, 0));
+                }
+                else if (!ownLive && nbrLive)
+                {
+                    seedFaces.append(facei);
+                    const label globalNbr = globalCells.toGlobal(nbr);
+                    seedData.append(topoDistanceData<int>(globalNbr, 0));
+                }
+            }
+
+            forAll(nbrTypes, bFacei)
+            {
+                label facei = bFacei+mesh_.nInternalFaces();
+                label own = owner[facei];
+                bool ownLive =
+                (
+                    cellTypes_[own] == cellCellStencil::CALCULATED
+                 || cellTypes_[own] == cellCellStencil::INTERPOLATED
+                );
+                bool nbrLive =
+                (
+                    nbrTypes[bFacei] == cellCellStencil::CALCULATED
+                 || nbrTypes[bFacei] == cellCellStencil::INTERPOLATED
+                );
+                if (ownLive && !nbrLive)
+                {
+                    seedFaces.append(facei);
+                    const label globalOwn = globalCells.toGlobal(own);
+                    seedData.append(topoDistanceData<int>(globalOwn, 0));
+                }
+            }
+
+            // Propagate information inwards
+            FaceCellWave<topoDistanceData<int>> distanceCalc
+            (
+                mesh_,
+                seedFaces,
+                seedData,
+                faceData,
+                cellData,
+                mesh_.globalData().nTotalCells()+1,
+                dummyTrackData
+            );
+        }
+
+
+        // Add the new donor ones
+        forAll(cellData, celli)
+        {
+            if
+            (
+                cellTypes_[celli] == cellCellStencil::HOLE
+             && isDonor[celli]
+             && cellData[celli].valid(dummyTrackData)
+            )
+            {
+                cellStencil_[celli].setSize(1);
+                cellStencil_[celli] = cellData[celli].data();   // global cellID
+                cellInterpolationWeights_[celli].setSize(1);
+                cellInterpolationWeights_[celli] = 1.0;
+                cellTypes_[celli] = SPECIAL;        // handled later on
+                interpolationCells_.append(celli);
+                cellInterpolationWeight_[celli] = 1.0;
+                nSpecialFar++;
+            }
+        }
+
+        if (debug&2)
+        {
+            reduce(nSpecialNear, sumOp<label>());
+            reduce(nSpecialFar, sumOp<label>());
+
+            Pout<< "Detected " << nHoleDonors
+                << " hole cells that are used as donors of which" << nl
+                << "    next to live cells : " << nSpecialNear << nl
+                << "    other              : " << nSpecialFar << nl
+                << endl;
+        }
+
+
+        // Re-do the mapDistribute
+        List<Map<label>> compactMap;
+        cellInterpolationMap_.reset
+        (
+            new mapDistribute
+            (
+                globalCells,
+                cellStencil_,
+                compactMap
+            )
+        );
+    }
+}
+
+
+void Foam::cellCellStencils::inverseDistance::createStencil
+(
+    const globalIndex& globalCells,
+    const bool allowHoleDonors
+)
 {
     // Send cell centre back to donor
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -1489,15 +1797,17 @@ void Foam::cellCellStencils::inverseDistance::createStencil
     const vector greatPoint(GREAT, GREAT, GREAT);
 
     boolList isValidDonor(mesh_.nCells(), true);
-    forAll(cellTypes_, celli)
+    if (!allowHoleDonors)
     {
-        if (cellTypes_[celli] == HOLE)
+        forAll(cellTypes_, celli)
         {
-            isValidDonor[celli] = false;
+            if (cellTypes_[celli] == HOLE)
+            {
+                isValidDonor[celli] = false;
+            }
         }
     }
 
-
     // Has acceptor been handled already?
     bitSet doneAcceptor(interpolationCells_.size());
 
@@ -1516,22 +1826,26 @@ void Foam::cellCellStencils::inverseDistance::createStencil
                 const point& cc = mesh_.cellCentres()[cellI];
                 const labelList& slots = cellStencil_[cellI];
 
-                if (slots.size() != 1)
+                //- Note: empty slots can happen for interpolated cells with
+                //        bad/insufficient donors. These are handled later on.
+                if (slots.size() > 1)
                 {
                     FatalErrorInFunction<< "Problem:" << slots
                         << abort(FatalError);
                 }
-
+                else if (slots.size())
-                forAll(slots, slotI)
                 {
-                    label elemI = slots[slotI];
+                    forAll(slots, slotI)
-                    //Pout<< "    acceptor:" << cellI
+                    {
-                    //    << " at:" << mesh_.cellCentres()[cellI]
+                        label elemI = slots[slotI];
-                    //    << " global:" << globalCells.toGlobal(cellI)
+                        //Pout<< "    acceptor:" << cellI
-                    //    << " found in donor:" << elemI << endl;
+                        //    << " at:" << mesh_.cellCentres()[cellI]
-                    minMagSqrEqOp<point>()(samples[elemI], cc);
+                        //    << " global:" << globalCells.toGlobal(cellI)
+                        //    << " found in donor:" << elemI << endl;
+                        minMagSqrEqOp<point>()(samples[elemI], cc);
+                    }
+                    nSamples++;
                 }
-                nSamples++;
             }
         }
 
@@ -1614,25 +1928,27 @@ void Foam::cellCellStencils::inverseDistance::createStencil
                 label cellI = interpolationCells_[i];
                 const labelList& slots = cellStencil_[cellI];
 
-                if (slots.size() != 1)
+                if (slots.size() > 1)
                 {
                     FatalErrorInFunction << "Problem:" << slots
                         << abort(FatalError);
                 }
-
+                else if (slots.size() == 1)
-                label slotI = slots[0];
-
-                // Important: check if the stencil is actually for this cell
-                if (samples[slotI] == mesh_.cellCentres()[cellI])
                 {
-                    cellStencil_[cellI].transfer(donorCellCells[slotI]);
+                    label slotI = slots[0];
-                    cellInterpolationWeights_[cellI].transfer
+
-                    (
+                    // Important: check if the stencil is actually for this cell
-                        donorWeights[slotI]
+                    if (samples[slotI] == mesh_.cellCentres()[cellI])
-                    );
+                    {
-                    // Mark cell as being done so it does not get sent over
+                        cellStencil_[cellI].transfer(donorCellCells[slotI]);
-                    // again.
+                        cellInterpolationWeights_[cellI].transfer
-                    doneAcceptor.set(i);
+                        (
+                            donorWeights[slotI]
+                        );
+                        // Mark cell as being done so it does not get sent over
+                        // again.
+                        doneAcceptor.set(i);
+                    }
                 }
             }
         }
@@ -1663,6 +1979,11 @@ Foam::cellCellStencils::inverseDistance::inverseDistance
 :
     cellCellStencil(mesh),
     dict_(dict),
+    allowHoleDonors_(dict.getOrDefault<bool>("allowHoleDonors", false)),
+    allowInterpolatedDonors_
+    (
+        dict.getOrDefault<bool>("allowInterpolatedDonors", true)
+    ),
     smallVec_(Zero),
     cellTypes_(labelList(mesh.nCells(), CALCULATED)),
     interpolationCells_(0),
@@ -2001,13 +2322,19 @@ bool Foam::cellCellStencils::inverseDistance::update()
         }
     }
 
-    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();
     }
 
     // Use the patch types and weights to decide what to do
@@ -2020,23 +2347,12 @@ bool Foam::cellCellStencils::inverseDistance::update()
                 case OVERSET:
                 {
                     // Require interpolation. See if possible.
-
                     if (allStencil[cellI].size())
                     {
                         allCellTypes[cellI] = INTERPOLATED;
                     }
                     else
                     {
-                        // If there are no donors we can either set the
-                        // acceptor cell to 'hole' i.e. nothing gets calculated
-                        // if the acceptor cells go outside the donor meshes or
-                        // we can set it to 'calculated' to have something
-                        // like an acmi type behaviour where only covered
-                        // acceptor cell use the interpolation and non-covered
-                        // become calculated. Uncomment below line. Note: this
-                        // should be switchable?
-                        //allCellTypes[cellI] = CALCULATED;
-
                         allCellTypes[cellI] = HOLE;
                     }
                 }
@@ -2044,15 +2360,55 @@ bool Foam::cellCellStencils::inverseDistance::update()
         }
     }
 
-    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);
+
+
+    if ((debug&2) && (mesh_.time().outputTime()))
+    {
+        tmp<volScalarField> tfld
+        (
+            createField(mesh_, "allCellTypes_hole", allCellTypes)
+        );
+        tfld().write();
+    }
+    if ((debug&2) && (mesh_.time().outputTime()))
+    {
+        labelList stencilSize(mesh_.nCells());
+        forAll(allStencil, celli)
+        {
+            stencilSize[celli] = allStencil[celli].size();
+        }
+        tmp<volScalarField> tfld
+        (
+            createField(mesh_, "allStencil_hole", stencilSize)
+        );
+        tfld().write();
+    }
+
+     // Update allStencil with new fill HOLES
+    forAll(allCellTypes, celli)
+    {
+        if (allCellTypes[celli] == HOLE && allStencil[celli].size())
+        {
+            allStencil[celli].clear();
+        }
+    }
 
     // Check previous iteration cellTypes_ for any hole->calculated changes
     // If so set the cell either to interpolated (if there are donors) or
@@ -2087,38 +2443,12 @@ bool Foam::cellCellStencils::inverseDistance::update()
         }
     }
 
-
-    // Mark unreachable bits
-    findHoles(globalCells, mesh_, zoneID, allStencil, allCellTypes);
-
-    if (debug)
-    {
-        tmp<volScalarField> tfld
-        (
-            createField(mesh_, "allCellTypes_hole", allCellTypes)
-        );
-        tfld().write();
-    }
-    if (debug)
-    {
-        labelList stencilSize(mesh_.nCells());
-        forAll(allStencil, celli)
-        {
-            stencilSize[celli] = allStencil[celli].size();
-        }
-        tmp<volScalarField> tfld
-        (
-            createField(mesh_, "allStencil_hole", stencilSize)
-        );
-        tfld().write();
-    }
-
-
     // Add buffer interpolation layer(s) around holes
     scalarField allWeight(mesh_.nCells(), Zero);
-    walkFront(layerRelax, allStencil, allCellTypes, allWeight);
+    walkFront(layerRelax, allStencil, allCellTypes, allWeight, zoneID);
 
-    if (debug)
+
+    if ((debug&2) && (mesh_.time().outputTime()))
     {
         tmp<volScalarField> tfld
         (
@@ -2127,21 +2457,32 @@ bool Foam::cellCellStencils::inverseDistance::update()
         tfld().write();
     }
 
-
     // Convert cell-cell addressing to stencil in compact notation
-
+    // Clean any potential INTERPOLATED with HOLE donors.
+    // This is not eliminated in the front walk as the HOLE cells might
+    // leak when inset region is overlapping backgroung mesh
     cellTypes_.transfer(allCellTypes);
     cellStencil_.setSize(mesh_.nCells());
     cellInterpolationWeights_.setSize(mesh_.nCells());
     DynamicList<label> interpolationCells;
+
+/*
     forAll(cellTypes_, cellI)
     {
         if (cellTypes_[cellI] == INTERPOLATED)
         {
-            cellStencil_[cellI].transfer(allStencil[cellI]);
+            if (cellTypes_[allStencil[cellI][0]] != HOLE)
-            cellInterpolationWeights_[cellI].setSize(1);
+            {
-            cellInterpolationWeights_[cellI][0] = 1.0;
+                cellStencil_[cellI].transfer(allStencil[cellI]);
-            interpolationCells.append(cellI);
+                cellInterpolationWeights_[cellI].setSize(1);
+                cellInterpolationWeights_[cellI][0] = 1.0;                 interpolationCells.append(cellI);
+            }
+            else
+            {
+                cellTypes_[cellI] = POROUS;//CALCULATED;
+                cellStencil_[cellI].clear();
+                cellInterpolationWeights_[cellI].clear();
+            }
         }
         else
         {
@@ -2149,22 +2490,76 @@ bool Foam::cellCellStencils::inverseDistance::update()
             cellInterpolationWeights_[cellI].clear();
         }
     }
-    interpolationCells_.transfer(interpolationCells);
+*/
 
+    labelListList compactStencil(allStencil);
+    List<Map<label>> compactStencilMap;
+    mapDistribute map(globalCells, compactStencil, compactStencilMap);
+
+    labelList compactCellTypes(cellTypes_);
+    map.distribute(compactCellTypes);
+
+    label nHoleDonors = 0;
+    forAll(cellTypes_, celli)
+    {
+        if (cellTypes_[celli] == INTERPOLATED)
+        {
+            const labelList& slots = compactStencil[celli];
+            if (slots.size())
+            {
+                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
+        {
+            cellStencil_[celli].clear();
+            cellInterpolationWeights_[celli].clear();
+        }
+    }
+
+    reduce(nHoleDonors, sumOp<label>());
+
+    DebugInfo<< FUNCTION_NAME << "nHole Donors : " << nHoleDonors << endl;
+
+    // Reset map with new cellStencil
     List<Map<label>> compactMap;
     cellInterpolationMap_.reset
     (
         new mapDistribute(globalCells, cellStencil_, compactMap)
     );
+
+    interpolationCells_.transfer(interpolationCells);
+
     cellInterpolationWeight_.transfer(allWeight);
-    dynamicOversetFvMesh::correctBoundaryConditions
+    oversetFvMeshBase::correctBoundaryConditions
     <
         volScalarField,
         oversetFvPatchField<scalar>
     >(cellInterpolationWeight_.boundaryFieldRef(), false);
 
 
-    if (debug&2)
+    if ((debug&2) && (mesh_.time().outputTime()))
     {
         // Dump mesh
         mesh_.time().write();
@@ -2194,11 +2589,18 @@ bool Foam::cellCellStencils::inverseDistance::update()
 
 
     // Extend stencil to get inverse distance weighted neighbours
-    createStencil(globalCells);
+    createStencil(globalCells, allowHoleDonors_);
 
-
+    // Optional: convert hole cells next to non-hole cells into
-    if (debug&2)
+    // interpolate-from-neighbours (of cell type SPECIAL)
+    if (allowHoleDonors_)
     {
+        holeExtrapolationStencil(globalCells);
+    }
+
+    if ((debug&2) && (mesh_.time().outputTime()))
+    {
+
         // Dump weight
         cellInterpolationWeight_.instance() = mesh_.time().timeName();
         cellInterpolationWeight_.write();
@@ -2235,7 +2637,7 @@ bool Foam::cellCellStencils::inverseDistance::update()
             (
                 createField(mesh_, "maxMagWeight", maxMagWeight)
             );
-            dynamicOversetFvMesh::correctBoundaryConditions
+            oversetFvMeshBase::correctBoundaryConditions
             <
                 volScalarField,
                 oversetFvPatchField<scalar>
@@ -2250,7 +2652,7 @@ bool Foam::cellCellStencils::inverseDistance::update()
                 createField(mesh_, "cellTypes", cellTypes_)
             );
             //tfld.ref().correctBoundaryConditions();
-            dynamicOversetFvMesh::correctBoundaryConditions
+            oversetFvMeshBase::correctBoundaryConditions
             <
                 volScalarField,
                 oversetFvPatchField<scalar>

src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.H

@@ -76,6 +76,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_;
 
@@ -217,6 +223,8 @@ protected:
             const boolList& isBlockedFace,
             labelList& cellRegion
         ) const;
+
+
         autoPtr<globalIndex> compactedRegionSplit
         (
             const fvMesh& mesh,
@@ -250,11 +258,22 @@ protected:
             const scalar layerRelax,
             const labelListList& allStencil,
             labelList& allCellTypes,
-            scalarField& allWeight
+            scalarField& allWeight,
+            const labelList& zoneID
         ) 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
+        //  interpolation stencil
+        void holeExtrapolationStencil(const globalIndex& globalCells);
 
 
 private:
@@ -286,6 +305,7 @@ public:
 
     // Member Functions
 
+
         //- Update stencils. Return false if nothing changed.
         virtual bool update();
 

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,27 @@ void Foam::cellCellStencils::trackingInverseDistance::markDonors
 
     const fvMesh& srcMesh = meshParts_[srcI].subMesh();
     const labelList& srcCellMap = meshParts_[srcI].cellMap();
-    const voxelMeshSearch& meshSearch = meshSearches[srcI];
+    //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]);
+            //label srcCelli = meshSearch.findCell(tgtCc[tgtCelli]);
-            if (srcCelli != -1 && allCellTypes[srcCellMap[srcCelli]] != HOLE)
+            label srcCelli = tgtToSrcAddr[tgtCelli];
+
+            if
+            (
+                srcCelli != -1
+             //&& allCellTypes[srcCellMap[srcCelli]] != HOLE
+             && allCellTypes[tgtCellMap[tgtCelli]] != HOLE
+            )
             {
                 label celli = tgtCellMap[tgtCelli];
 
@@ -444,8 +455,10 @@ void Foam::cellCellStencils::trackingInverseDistance::markDonors
         labelList donors(samples.size(), -1);
         forAll(samples, sampleI)
         {
-            label srcCelli = meshSearch.findCell(samples[sampleI]);
+//             label srcCelli = meshSearch.findCell(samples[sampleI]);
-            if (srcCelli != -1 && allCellTypes[srcCellMap[srcCelli]] != HOLE)
+            const point& sample = samples[sampleI];
+            label srcCelli = srcMesh.findCell(sample, polyMesh::CELL_TETS);
+            if (srcCelli != -1)// && allCellTypes[srcCellMap[srcCelli]] != HOLE)
             {
                 donors[sampleI] = globalCells_.toGlobal(srcCellMap[srcCelli]);
             }
@@ -738,7 +751,7 @@ bool Foam::cellCellStencils::trackingInverseDistance::update()
     }
 
 
-    if (debug)
+    if (debug&2)
     {
         forAll(patchParts, zonei)
         {
@@ -834,13 +847,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();
     }
 
 
@@ -869,15 +888,65 @@ 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();
+        }
+    }
+
+    // Try to make calculated cells that are nrb with patches
+//     forAll(allPatchTypes, cellI)
+//     {
+//         if (allPatchTypes[cellI] == patchCellType::PATCH)
+//         {
+//             if (allCellTypes[cellI] == INTERPOLATED)
+//             {
+//                 allCellTypes[cellI] = CALCULATED;
+//             }
+//         }
+//     }
 
     // Check previous iteration cellTypes_ for any hole->calculated changes
     // If so set the cell either to interpolated (if there are donors) or
@@ -913,25 +982,14 @@ 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);
+    walkFront(layerRelax, allStencil, allCellTypes, allWeight, zoneID);
     DebugInfo<< FUNCTION_NAME << " : Implemented layer relaxation" << endl;
 
-    if (debug)
+
+    if ((debug&2) && (mesh_.time().outputTime()))
     {
         tmp<volScalarField> tfld
         (
@@ -947,14 +1005,71 @@ bool Foam::cellCellStencils::trackingInverseDistance::update()
     cellStencil_.setSize(mesh_.nCells());
     cellInterpolationWeights_.setSize(mesh_.nCells());
     DynamicList<label> interpolationCells;
+    /*
+    forAll(cellTypes_, cellI)
+    {
+        if (cellTypes_[cellI] == INTERPOLATED)
+        {
+            if (cellTypes_[allStencil[cellI][0]] != HOLE)
+            {
+                cellStencil_[cellI].transfer(allStencil[cellI]);
+                cellInterpolationWeights_[cellI].setSize(1);
+                cellInterpolationWeights_[cellI][0] = 1.0;
+                interpolationCells.append(cellI);
+            }
+            else
+            {
+                cellTypes_[cellI] = POROUS;
+                cellStencil_[cellI].clear();
+                cellInterpolationWeights_[cellI].clear();
+            }
+        }
+        else
+        {
+            cellStencil_[cellI].clear();
+            cellInterpolationWeights_[cellI].clear();
+        }
+    }
+    */
+
+    labelListList compactStencil(allStencil);
+    List<Map<label>> compactStencilMap;
+    mapDistribute map(globalCells, compactStencil, compactStencilMap);
+
+    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
         {
@@ -962,6 +1077,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;
@@ -975,16 +1095,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");
@@ -1012,45 +1134,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) 2016-2017 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,56 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2016-2017 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

@@ -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,5 @@
+bool adjustFringe
+(
+    pimple.dict().getOrDefault("oversetAdjustPhi", false)
+);
+

src/overset/include/setCellMask.H

@@ -35,12 +35,18 @@ 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/oversetFvMesh/dynamicOversetFvMesh/dynamicOversetFvMesh.C

@@ -0,0 +1,80 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014-2018 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify i
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "dynamicOversetFvMesh.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(dynamicOversetFvMesh, 0);
+    addToRunTimeSelectionTable(dynamicFvMesh, dynamicOversetFvMesh, IOobject);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::dynamicOversetFvMesh::dynamicOversetFvMesh(const IOobject& io)
+:
+    dynamicMotionSolverFvMesh(io),
+    oversetFvMeshBase(static_cast<const fvMesh&>(*this))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::dynamicOversetFvMesh::~dynamicOversetFvMesh()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::dynamicOversetFvMesh::update()
+{
+    if (dynamicMotionSolverFvMesh::update())
+    {
+        oversetFvMeshBase::update();
+
+        return true;
+    }
+
+    return false;
+}
+
+
+bool Foam::dynamicOversetFvMesh::writeObject
+(
+    IOstreamOption streamOpt,
+    const bool valid
+) const
+{
+    bool ok = dynamicMotionSolverFvMesh::writeObject(streamOpt, valid);
+    ok = oversetFvMeshBase::writeObject(streamOpt, valid) && ok;
+    return ok;
+}
+
+
+// ************************************************************************* //

src/overset/oversetFvMesh/dynamicOversetFvMesh/dynamicOversetFvMesh.H

@@ -0,0 +1,224 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015-2018 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::dynamicOversetFvMesh
+
+Description
+    dynamicFvMesh with support for overset meshes.
+
+SourceFiles
+    dynamicOversetFvMesh.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef dynamicOversetFvMesh_H
+#define dynamicOversetFvMesh_H
+
+#include "dynamicMotionSolverFvMesh.H"
+#include "oversetFvMeshBase.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                   Class dynamicOversetFvMesh Declaration
+\*---------------------------------------------------------------------------*/
+
+class dynamicOversetFvMesh
+:
+    public dynamicMotionSolverFvMesh,
+    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);
+
+
+    //- 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/dynamicOversetFvMesh/dynamicOversetFvMeshTemplates.C

@@ -2,10 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
+    \\  /    A nd           | Copyright (C) 2014-2017 OpenCFD Ltd.
      \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2014-2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -32,10 +30,9 @@ License
 #include "calculatedProcessorFvPatchField.H"
 #include "lduInterfaceFieldPtrsList.H"
 #include "processorFvPatch.H"
-#include "syncTools.H"
 
 // * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * * //
-
+/*
 template<class T>
 void Foam::dynamicOversetFvMesh::interpolate(Field<T>& psi) const
 {
@@ -151,17 +148,18 @@ Foam::tmp<Foam::scalarField> Foam::dynamicOversetFvMesh::normalisation
     const fvMatrix<Type>& m
 ) const
 {
-    // Determine normalisation. This is normally the original diagonal plus
+    // Determine normalisation. This is normally the original diagonal.
-    // remote contributions. This needs to be stabilised for hole cells
+    // This needs to be stabilised for hole cells
     // which can have a zero diagonal. Assume that if any component has
     // a non-zero diagonal the cell does not need stabilisation.
     tmp<scalarField> tnorm(tmp<scalarField>::New(m.diag()));
     scalarField& norm = tnorm.ref();
 
-    // Add remote coeffs to duplicate behaviour of fvMatrix::addBoundaryDiag
+    // Add boundary coeffs to duplicate behaviour of fvMatrix
     const FieldField<Field, Type>& internalCoeffs = m.internalCoeffs();
     for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
     {
+        //m.addBoundaryDiag(norm, cmpt);
         forAll(internalCoeffs, patchi)
         {
             const labelUList& fc = lduAddr().patchAddr(patchi);
@@ -174,151 +172,17 @@ Foam::tmp<Foam::scalarField> Foam::dynamicOversetFvMesh::normalisation
         }
     }
 
-    // Count number of problematic cells
-    label nZeroDiag = 0;
     forAll(norm, celli)
     {
-        const scalar& n = norm[celli];
+        scalar& n = norm[celli];
-        if (magSqr(n) < sqr(SMALL))
+        if (mag(n) < SMALL)
         {
-            //Pout<< "For field " << m.psi().name()
+            n = 1.0; //?
-            //    << " have diagonal " << n << " for cell " << celli
-            //    << " at:" << cellCentres()[celli] << endl;
-            nZeroDiag++;
         }
-    }
+        else
-
-    reduce(nZeroDiag, sumOp<label>());
-
-    if (debug)
-    {
-        Pout<< "For field " << m.psi().name() << " have zero diagonals for "
-            << nZeroDiag << " cells" << endl;
-    }
-
-    if (nZeroDiag > 0)
-    {
-        // Walk out the norm across hole cells
-
-        const labelList& own = faceOwner();
-        const labelList& nei = faceNeighbour();
-        const cellCellStencilObject& overlap = Stencil::New(*this);
-        const labelUList& types = overlap.cellTypes();
-
-        label nHoles = 0;
-        scalarField extrapolatedNorm(norm);
-        forAll(types, celli)
         {
-            if (types[celli] == cellCellStencil::HOLE)
+            // Restore original diagonal
-            {
+            n = m.diag()[celli];
-                extrapolatedNorm[celli] = -GREAT;
-                nHoles++;
-            }
-        }
-
-        bitSet isFront(nFaces());
-        for (label facei = 0; facei < nInternalFaces(); facei++)
-        {
-            label ownType = types[own[facei]];
-            label neiType = types[nei[facei]];
-            if
-            (
-                (ownType == cellCellStencil::HOLE)
-             != (neiType == cellCellStencil::HOLE)
-            )
-            {
-                isFront.set(facei);
-            }
-        }
-        labelList nbrTypes;
-        syncTools::swapBoundaryCellList(*this, types, nbrTypes);
-        for (label facei = nInternalFaces(); facei < nFaces(); facei++)
-        {
-            label ownType = types[own[facei]];
-            label neiType = nbrTypes[facei-nInternalFaces()];
-            if
-            (
-                (ownType == cellCellStencil::HOLE)
-             != (neiType == cellCellStencil::HOLE)
-            )
-            {
-                isFront.set(facei);
-            }
-        }
-
-
-        while (true)
-        {
-            scalarField nbrNorm;
-            syncTools::swapBoundaryCellList(*this, extrapolatedNorm, nbrNorm);
-
-            bitSet newIsFront(nFaces());
-            scalarField newNorm(extrapolatedNorm);
-
-            label nChanged = 0;
-            for (const label facei : isFront)
-            {
-                if (extrapolatedNorm[own[facei]] == -GREAT)
-                {
-                    // Average owner cell, add faces to newFront
-                    newNorm[own[facei]] = cellAverage
-                    (
-                        types,
-                        nbrTypes,
-                        extrapolatedNorm,
-                        nbrNorm,
-                        own[facei],
-                        newIsFront
-                    );
-                    nChanged++;
-                }
-                if
-                (
-                    isInternalFace(facei)
-                 && extrapolatedNorm[nei[facei]] == -GREAT
-                )
-                {
-                    // Average nei cell, add faces to newFront
-                    newNorm[nei[facei]] = cellAverage
-                    (
-                        types,
-                        nbrTypes,
-                        extrapolatedNorm,
-                        nbrNorm,
-                        nei[facei],
-                        newIsFront
-                    );
-                    nChanged++;
-                }
-            }
-
-            reduce(nChanged, sumOp<label>());
-            if (nChanged == 0)
-            {
-                break;
-            }
-
-            // Transfer new front
-            extrapolatedNorm.transfer(newNorm);
-            isFront.transfer(newIsFront);
-            syncTools::syncFaceList(*this, isFront, maxEqOp<unsigned int>());
-        }
-
-
-        forAll(norm, celli)
-        {
-            scalar& n = norm[celli];
-            if (magSqr(n) < sqr(SMALL))
-            {
-                //Pout<< "For field " << m.psi().name()
-                //    << " for cell " << celli
-                //    << " at:" << cellCentres()[celli]
-                //    << " have norm " << n
-                //    << " have extrapolated norm " << extrapolatedNorm[celli]
-                //    << endl;
-                // Override the norm
-                n = extrapolatedNorm[celli];
-            }
         }
     }
     return tnorm;
@@ -372,9 +236,7 @@ void Foam::dynamicOversetFvMesh::addInterpolation
     inplaceReorder(reverseFaceMap_, lower);
 
 
-    //const label nOldInterfaces = dynamicMotionSolverFvMesh::interfaces().size();
+    const label nOldInterfaces = dynamicMotionSolverFvMesh::interfaces().size();
-    const label nOldInterfaces = dynamicFvMesh::interfaces().size();
-
 
     if (interfaces.size() > nOldInterfaces)
     {
@@ -671,8 +533,7 @@ Foam::SolverPerformance<Type> Foam::dynamicOversetFvMesh::solve
                 << " bypassing matrix adjustment for field " << m.psi().name()
                 << endl;
         }
-        //return dynamicMotionSolverFvMesh::solve(m, dict);
+        return dynamicMotionSolverFvMesh::solve(m, dict);
-        return dynamicFvMesh::solve(m, dict);
     }
 
     if (debug)
@@ -685,38 +546,6 @@ Foam::SolverPerformance<Type> Foam::dynamicOversetFvMesh::solve
     // Calculate stabilised diagonal as normalisation for interpolation
     const scalarField norm(normalisation(m));
 
-    if (debug)
-    {
-        volScalarField scale
-        (
-            IOobject
-            (
-                m.psi().name() + "_normalisation",
-                this->time().timeName(),
-                *this,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            *this,
-            dimensionedScalar(dimless, Zero)
-        );
-        scale.ref().field() = norm;
-        correctBoundaryConditions
-        <
-            volScalarField,
-            oversetFvPatchField<scalar>
-        >(scale.boundaryFieldRef(), false);
-        scale.write();
-
-        if (debug)
-        {
-            Pout<< "dynamicOversetFvMesh::solve() :"
-                << " writing matrix normalisation for field " << m.psi().name()
-                << " to " << scale.name() << endl;
-        }
-    }
-
 
     // Switch to extended addressing (requires mesh::update() having been
     // called)
@@ -748,8 +577,7 @@ Foam::SolverPerformance<Type> Foam::dynamicOversetFvMesh::solve
     //}
 
     // Use lower level solver
-    //SolverPerformance<Type> s(dynamicMotionSolverFvMesh::solve(m, dict));
+    SolverPerformance<Type> s(dynamicMotionSolverFvMesh::solve(m, dict));
-    SolverPerformance<Type> s(dynamicFvMesh::solve(m, dict));
 
     // Restore boundary field
     bpsi.setSize(oldSize);
@@ -892,14 +720,14 @@ void Foam::dynamicOversetFvMesh::write
 
     os  << "** End of Matrix **" << endl;
 }
-
+*/
-
+/*
 template<class GeoField>
 void Foam::dynamicOversetFvMesh::correctCoupledBoundaryConditions(GeoField& fld)
 {
     typename GeoField::Boundary& bfld = fld.boundaryFieldRef();
 
-    const label nReq = Pstream::nRequests();
+    label nReq = Pstream::nRequests();
 
     forAll(bfld, patchi)
     {
@@ -956,6 +784,6 @@ void Foam::dynamicOversetFvMesh::checkCoupledBC(const GeoField& fld)
     }
     Pout<< "** end of " << fld.name() << endl;
 }
-
+*/
 
 // ************************************************************************* //

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,7 +88,7 @@ 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;
@@ -266,7 +265,7 @@ bool Foam::dynamicOversetFvMesh::updateAddressing() const
         {
             if (interface != -1)
             {
-                interface = procToInterface[interface]+boundary().size();
+                interface = procToInterface[interface]+mesh_.boundary().size();
             }
         }
     }
@@ -277,12 +276,12 @@ bool Foam::dynamicOversetFvMesh::updateAddressing() const
 
     List<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)
@@ -317,7 +316,7 @@ bool Foam::dynamicOversetFvMesh::updateAddressing() const
     (
         new fvMeshPrimitiveLduAddressing
         (
-            nCells(),
+            mesh_.nCells(),
             std::move(lowerAddr),
             std::move(upperAddr),
             patchAddr,
@@ -344,7 +343,7 @@ bool Foam::dynamicOversetFvMesh::updateAddressing() const
         }
 
         // Using interfaces
-        const lduInterfacePtrsList& iFaces = allInterfaces_;
+        const lduInterfacePtrsList& iFaces = mesh_.interfaces();
         Pout<< "Adapted interFaces:" << iFaces.size() << endl;
         forAll(iFaces, patchI)
         {
@@ -360,7 +359,7 @@ bool Foam::dynamicOversetFvMesh::updateAddressing() const
 }
 
 
-Foam::scalar Foam::dynamicOversetFvMesh::cellAverage
+Foam::scalar Foam::oversetFvMeshBase::cellAverage
 (
     const labelList& types,
     const labelList& nbrTypes,
@@ -370,16 +369,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 +398,7 @@ Foam::scalar Foam::dynamicOversetFvMesh::cellAverage
         }
         else
         {
-            if (nbrNorm[facei-nInternalFaces()] == -GREAT)
+            if (nbrNorm[facei-mesh_.nInternalFaces()] == -GREAT)
             {
                 if (isFront.set(facei))
                 {
@@ -408,7 +407,7 @@ Foam::scalar Foam::dynamicOversetFvMesh::cellAverage
             }
             else
             {
-                avg += nbrNorm[facei-nInternalFaces()];
+                avg += nbrNorm[facei-mesh_.nInternalFaces()];
                 n++;
             }
         }
@@ -425,13 +424,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 +439,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 +462,7 @@ void Foam::dynamicOversetFvMesh::writeAgglomeration
         scalarAgglomeration.write();
 
         Info<< "Writing initial cell distribution to "
-            << this->time().timeName() << endl;
+            << mesh_.time().timeName() << endl;
     }
 
 
@@ -496,13 +495,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();
@@ -527,52 +526,53 @@ 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
+//     if (doInit)
-    return true;
+//     {
+//         init(false);    // do not initialise lower levels
+//     }
 }
 
 
+
+// bool Foam::oversetFvMeshBase::init(const bool doInit)
+// {
+//     if (doInit)
+//     {
+//         dynamicMotionSolverListFvMesh::init(doInit);
+//     }
+//
+//     active_ = false;
+//
+//     // Load stencil (but do not update)
+//     (void)Stencil::New(*this, 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_)
     {
@@ -583,12 +583,13 @@ 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 dynamicFvMesh::interfaces();
+        //return dynamicFvMesh::interfaces();
+        return mesh_.fvMesh::interfaces();
     }
     if (!lduPtr_)
     {
@@ -600,7 +601,7 @@ Foam::lduInterfacePtrsList Foam::dynamicOversetFvMesh::interfaces() const
 
 
 const Foam::fvMeshPrimitiveLduAddressing&
-Foam::dynamicOversetFvMesh::primitiveLduAddr() const
+Foam::oversetFvMeshBase::primitiveLduAddr() const
 {
     if (!lduPtr_)
     {
@@ -612,10 +613,10 @@ Foam::dynamicOversetFvMesh::primitiveLduAddr() const
 }
 
 
-bool Foam::dynamicOversetFvMesh::update()
+bool Foam::oversetFvMeshBase::update()
 {
     //if (dynamicMotionSolverFvMesh::update())
-    if (dynamicMotionSolverListFvMesh::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.
@@ -632,55 +633,48 @@ 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::dynamicOversetFvMesh::writeObject
+bool Foam::oversetFvMeshBase::writeObject
 (
     IOstreamOption streamOpt,
     const bool valid
 ) const
 {
     //bool ok = dynamicMotionSolverFvMesh::writeObject(streamOpt, valid);
-    bool ok = dynamicMotionSolverListFvMesh::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();
 
@@ -689,13 +683,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
         );
@@ -713,18 +707,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)
@@ -736,7 +730,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();
 
@@ -745,7 +739,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
@@ -753,13 +747,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
         );
@@ -775,7 +769,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;
@@ -785,7 +779,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

@@ -24,46 +24,45 @@ 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 "labelIOList.H"
+#include "fvMeshPrimitiveLduAddressing.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
 {
 protected:
 
     // Protected data
 
+        const fvMesh& mesh_;
+
         //- Select base addressing (false) or locally stored extended
         //  lduAddressing (true)
         mutable bool active_;
@@ -105,22 +104,8 @@ protected:
             const lduInterfacePtrsList&
         ) const;
 
-        //- Explicit interpolation of acceptor cells from donor cells
+//         //- Helper: strip off trailing _0
-        template<class T>
+//         static word baseName(const word& name);
-        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>
@@ -130,19 +115,30 @@ protected:
         //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;
-
         //- Add interpolation to matrix (coefficients)
+        //template<class Type>
+        //void addInterpolation(fvMatrix<Type>&, const scalarField& norm) const;
+
+
+         //- Scale coefficient depending on cell type
         template<class Type>
-        void addInterpolation(fvMatrix<Type>&, const scalarField& norm) const;
+        void scaleConnection
+        (
+            Field<Type>& coeffs,
+            const labelUList& types,
+            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>
@@ -171,26 +167,26 @@ private:
     // Private Member Functions
 
         //- No copy construct
-        dynamicOversetFvMesh(const dynamicOversetFvMesh&) = delete;
+        oversetFvMeshBase(const oversetFvMeshBase&) = delete;
 
         //- No copy assignment
-        void operator=(const dynamicOversetFvMesh&) = delete;
+        void operator=(const oversetFvMeshBase&) = 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
@@ -231,12 +227,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 +242,22 @@ public:
 
         // Overset
 
-            // Explicit interpolation
 
-                virtual void interpolate(scalarField& psi) const
+                //- Manipulate the matrix to add the interpolation/set hole
-                {
+                //  values
-                    interpolate<scalar>(psi);
+                template<class Type>
-                }
+                void addInterpolation
-
-                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,
+                    fvMatrix<Type>& m,
-                    const dictionary& dict
+                    const scalarField& normalisation,
-                ) const
+                    const bool setHoleCellValue,
-                {
+                    const Type& holeCellValue
-                    return solve<scalar>(m, dict);
+                ) const;
-                }
 
-                //- 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);
+        //virtual bool init(const bool doInit);
 
         //- Update the mesh for both mesh motion and topology change
         virtual bool update();
@@ -372,6 +284,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 +300,7 @@ public:
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #ifdef NoRepository
-#   include "dynamicOversetFvMeshTemplates.C"
+#   include "oversetFvMeshBaseTemplates.C"
 #endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/overset/oversetFvMesh/staticOversetFvMesh/staticOversetFvMesh.C

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

@@ -28,7 +28,7 @@ License
 #include "volFields.H"
 #include "cellCellStencil.H"
 #include "cellCellStencilObject.H"
-#include "dynamicOversetFvMesh.H"
+#include "oversetFvMeshBase.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -40,7 +40,10 @@ Foam::oversetFvPatchField<Type>::oversetFvPatchField
 )
 :
     zeroGradientFvPatchField<Type>(p, iF),
-    oversetPatch_(refCast<const oversetFvPatch>(p))
+    oversetPatch_(refCast<const oversetFvPatch>(p)),
+    setHoleCellValue_(false),
+    interpolateHoleCellValue_(false),
+    holeCellValue_(pTraits<Type>::min)
 {}
 
 
@@ -54,7 +57,10 @@ Foam::oversetFvPatchField<Type>::oversetFvPatchField
 )
 :
     zeroGradientFvPatchField<Type>(ptf, p, iF, mapper),
-    oversetPatch_(refCast<const oversetFvPatch>(p))
+    oversetPatch_(refCast<const oversetFvPatch>(p)),
+    setHoleCellValue_(ptf.setHoleCellValue_),
+    interpolateHoleCellValue_(ptf.interpolateHoleCellValue_),
+    holeCellValue_(ptf.holeCellValue_)
 {}
 
 
@@ -67,7 +73,18 @@ Foam::oversetFvPatchField<Type>::oversetFvPatchField
 )
 :
     zeroGradientFvPatchField<Type>(p, iF, dict),
-    oversetPatch_(refCast<const oversetFvPatch>(p, dict))
+    oversetPatch_(refCast<const oversetFvPatch>(p, dict)),
+    setHoleCellValue_(dict.lookupOrDefault<bool>("setHoleCellValue", false)),
+    interpolateHoleCellValue_
+    (
+        dict.lookupOrDefault<bool>("interpolateHoleCellValue", false)
+    ),
+    holeCellValue_
+    (
+        setHoleCellValue_
+      ? dict.get<Type>("holeCellValue")
+      : pTraits<Type>::min
+    )
 {}
 
 
@@ -78,7 +95,10 @@ Foam::oversetFvPatchField<Type>::oversetFvPatchField
 )
 :
     zeroGradientFvPatchField<Type>(ptf),
-    oversetPatch_(ptf.oversetPatch_)
+    oversetPatch_(ptf.oversetPatch_),
+    setHoleCellValue_(ptf.setHoleCellValue_),
+    interpolateHoleCellValue_(ptf.interpolateHoleCellValue_),
+    holeCellValue_(ptf.holeCellValue_)
 {}
 
 
@@ -90,7 +110,10 @@ Foam::oversetFvPatchField<Type>::oversetFvPatchField
 )
 :
     zeroGradientFvPatchField<Type>(ptf, iF),
-    oversetPatch_(ptf.oversetPatch_)
+    oversetPatch_(ptf.oversetPatch_),
+    setHoleCellValue_(ptf.setHoleCellValue_),
+    interpolateHoleCellValue_(ptf.interpolateHoleCellValue_),
+    holeCellValue_(ptf.holeCellValue_)
 {}
 
 
@@ -205,6 +228,53 @@ void Foam::oversetFvPatchField<Type>::initEvaluate
                     Info<< "Interpolating non-suppressed field " << fldName
                         << endl;
                 }
+
+                // Interpolate without bc update (since would trigger infinite
+                // recursion back to oversetFvPatchField<Type>::evaluate)
+                // The only problem is bcs that use the new cell values
+                // (e.g. zeroGradient, processor). These need to appear -after-
+                // the 'overset' bc.
+
+                const cellCellStencil& overlap = Stencil::New(mesh);
+
+                Field<Type>& fld = const_cast<Field<Type>&>(this->primitiveField());
+
+                // tbd: different weights for different variables ...
+                cellCellStencil::interpolate
+                (
+                    fld,
+                    mesh,
+                    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;
+                    }
+                }
+                /*
                 mesh.interpolate
                 (
                     const_cast<Field<Type>&>
@@ -212,9 +282,260 @@ void Foam::oversetFvPatchField<Type>::initEvaluate
                         this->primitiveField()
                     )
                 );
+                */
             }
         }
     }
+
+    zeroGradientFvPatchField<Type>::initEvaluate(commsType);
+}
+
+
+template<class Type>
+void Foam::oversetFvPatchField<Type>::manipulateMatrix
+(
+    fvMatrix<Type>& matrix
+)
+{
+    //const word& fldName = this->internalField().name();
+
+    if (this->manipulatedMatrix())
+    {
+        return;
+    }
+
+    const oversetFvPatch& ovp = this->oversetPatch_;
+
+    if (ovp.master())
+    {
+        // 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())
+        {
+            //Pout<< FUNCTION_NAME << "SKIPPING MANIP field:" << fldName
+            //    << " patch:" << ovp.name() << endl;
+            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
+                )
+                {
+                    //FatalErrorInFunction
+                    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]
+                        //<< exit(FatalError);
+                        << 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>& intCoeffs =
+                //    matrix.internalCoeffs()[patchi];
+                const Field<Type>& bouCoeffs = matrix.boundaryCoeffs()[patchi];
+                forAll(fc, i)
+                {
+                    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 I am 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);
+                    }
+                }
+            }
+        }
+    }
+
+    zeroGradientFvPatchField<Type>::manipulateMatrix(matrix);
 }
 
 
@@ -222,6 +543,17 @@ template<class Type>
 void Foam::oversetFvPatchField<Type>::write(Ostream& os) const
 {
     zeroGradientFvPatchField<Type>::write(os);
+    if (this->setHoleCellValue_)
+    {
+        os.writeEntry("setHoleCellValue", setHoleCellValue_);
+        os.writeEntry("holeCellValue", holeCellValue_);
+        os.writeEntryIfDifferent
+        (
+            "interpolateHoleCellValue",
+            false,
+            interpolateHoleCellValue_
+        );
+    }
     // Make sure to write the value for ease of postprocessing.
     this->writeEntry("value", os);
 }

src/overset/oversetPolyPatch/oversetFvPatchField.H

@@ -68,6 +68,17 @@ protected:
         //- Master patch ID
         mutable label masterPatchID_;
 
+        // Hole cell controls
+
+            //- Set hole cell value
+            const bool setHoleCellValue_;
+
+            //- Interpolate hole cell value (from nearby non-hole cell)
+            const bool interpolateHoleCellValue_;
+
+            //- Hole cell value
+            const Type holeCellValue_;
+
 
 public:
 
@@ -140,6 +151,9 @@ public:
             //- Initialise the evaluation of the patch field
             virtual void initEvaluate(const Pstream::commsTypes commsType);
 
+            //- Manipulate matrix
+            virtual void manipulateMatrix(fvMatrix<Type>& matrix);
+
         // I-O
 
             //- Write

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

tutorials/incompressible/overPimpleDyMFoam/simpleRotor/system/fvSolution

@@ -70,13 +70,10 @@ solvers
 PIMPLE
 {
     momentumPredictor   no;
-    correctPhi          no;
     nOuterCorrectors    1;
     nCorrectors         2;
     nNonOrthogonalCorrectors 0;
 
-    ddtCorr             false;
-
     pRefPoint           (0.0001 0.0001 0.001);
     pRefValue           0.0;
 }

tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/system/fvSchemes

@@ -53,12 +53,13 @@ snGradSchemes
 
 oversetInterpolation
 {
-    method          cellVolumeWeight;
+    //method          cellVolumeWeight;
 
     // Faster but less accurate
-    //method              trackingInverseDistance;
+    method              trackingInverseDistance;
-    //searchBox           (0 0 0)(0.02 0.01 0.01);
+    searchBox           (0 0 0)(0.02 0.01 0.01);
-    //searchBoxDivisions  3{(64 64 1)};
+    searchBoxDivisions  3{(64 64 1)};
+    allowHoleDonors     false;
 }
 
 fluxRequired

tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/system/fvSolution

@@ -71,13 +71,10 @@ solvers
 PIMPLE
 {
     momentumPredictor   false;
-    correctPhi          no;
     nOuterCorrectors    1;
     nCorrectors         3;
     nNonOrthogonalCorrectors 0;
 
-    ddtCorr             true;
-
     pRefPoint           (0.0001 0.0001 0.001);
     pRefValue           0.0;
 }

tutorials/multiphase/overCompressibleInterDyMFoam/compressibleTwoSimpleRotors/constant/turbulenceProperties

@@ -15,7 +15,7 @@ FoamFile
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 simulationType      RAS;
-
+density             variable;
 RAS
 {
     RASModel        kEpsilon;

tutorials/multiphase/overCompressibleInterDyMFoam/compressibleTwoSimpleRotors/system/controlDict

@@ -87,24 +87,6 @@ functions
              (0.0009999 0.0015 0.003)
         );
     }
-
-    alphaVol
-    {
-        libs            (utilityFunctionObjects);
-        type            coded;
-        name            alphaVolume;
-        writeControl    timeStep;
-        writeInterval   10;
-
-        codeWrite
-        #{
-            const auto& alpha =
-                mesh().lookupObject<volScalarField>("alpha.water");
-
-            Info<< "Alpha volume = " << alpha.weightedAverage(mesh().Vsc())
-                << endl;
-        #};
-    }
 }
 
 

tutorials/multiphase/overCompressibleInterDyMFoam/compressibleTwoSimpleRotors/system/fvSolution

@@ -96,7 +96,6 @@ PIMPLE
     nOuterCorrectors    2;
     nCorrectors         3;
     nNonOrthogonalCorrectors 0;
-    ddtCorr             false;
 }
 
 relaxationFactors

tutorials/multiphase/overInterDyMFoam/twoSquaresOutDomain/0.orig/U

@@ -0,0 +1,55 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2012                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    #includeEtc "caseDicts/setConstraintTypes"
+
+    walls
+    {
+        type            uniformFixedValue;
+        uniformValue    (0 0 0);
+    }
+
+    hole
+    {
+        type            movingWallVelocity;
+        value           uniform (0 0 0);
+    }
+
+    outlet
+    {
+        type            uniformFixedValue;
+        uniformValue    (0 0 0);
+    }
+
+    inlet
+    {
+        type            uniformFixedValue;
+        uniformValue     (0 0 0);
+    }
+
+    overset
+    {
+        type            overset;
+    }
+}
+
+// ************************************************************************* //

tutorials/multiphase/overInterDyMFoam/twoSquaresOutDomain/0.orig/alpha.water

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

tutorials/multiphase/overInterDyMFoam/twoSquaresOutDomain/0.orig/epsilon

@@ -0,0 +1,37 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2012                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      epsilon;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [ 0 2 -3 0 0 0 0 ];
+
+internalField   uniform 0.1;
+
+boundaryField
+{
+    #includeEtc "caseDicts/setConstraintTypes"
+
+    overset
+    {
+        type            overset;
+    }
+
+    "(walls|hole|inlet|outlet)"
+    {
+        type            epsilonWallFunction;
+        value           $internalField;
+    }
+}
+
+// ************************************************************************* //

tutorials/multiphase/overInterDyMFoam/twoSquaresOutDomain/0.orig/k

@@ -0,0 +1,37 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2012                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      k;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [ 0 2 -2 0 0 0 0 ];
+
+internalField   uniform 0.01;
+
+boundaryField
+{
+    #includeEtc "caseDicts/setConstraintTypes"
+
+    overset
+    {
+        type            overset;
+    }
+
+    "(walls|hole|inlet|outlet)"
+    {
+        type            kqRWallFunction;
+        value           $internalField;
+    }
+}
+
+// ************************************************************************* //

tutorials/multiphase/overInterDyMFoam/twoSquaresOutDomain/0.orig/nut

@@ -0,0 +1,37 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2012                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      nut;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [ 0 2 -1 0 0 0 0 ];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    #includeEtc "caseDicts/setConstraintTypes"
+
+    overset
+    {
+        type            overset;
+    }
+
+    "(walls|hole|inlet|outlet)"
+    {
+        type            nutkWallFunction;
+        value           $internalField;
+    }
+}
+
+// ************************************************************************* //

tutorials/multiphase/overInterDyMFoam/twoSquaresOutDomain/0.orig/p

@@ -0,0 +1,37 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2012                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    #includeEtc "caseDicts/setConstraintTypes"
+
+    "(walls|hole|outlet|inlet)"
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    overset
+    {
+        type            overset;
+    }
+}
+
+// ************************************************************************* //

tutorials/multiphase/overInterDyMFoam/twoSquaresOutDomain/0.orig/p_rgh

@@ -0,0 +1,36 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2012                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p_rgh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    #includeEtc "caseDicts/setConstraintTypes"
+
+    oversetPatch
+    {
+        type            overset;
+    }
+
+    "(walls|hole|outlet|inlet)"
+    {
+        type            fixedFluxPressure;
+    }
+}
+
+// ************************************************************************* //

tutorials/multiphase/overInterDyMFoam/twoSquaresOutDomain/0.orig/pointDisplacement

@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2012                                 |
+|   \\  /    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);
+    }
+
+    hole
+    {
+        type            zeroGradient;
+    }
+
+    overset
+    {
+        patchType       overset;
+        type            zeroGradient;
+    }
+}
+
+// ************************************************************************* //

tutorials/multiphase/overInterDyMFoam/twoSquaresOutDomain/0.orig/zoneID

@@ -0,0 +1,36 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2012                                 |
+|   \\  /    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"
+
+    overset
+    {
+        type            overset;
+    }
+
+    ".*"
+    {
+        type            zeroGradient;
+    }
+}
+
+// ************************************************************************* //

tutorials/multiphase/overInterDyMFoam/twoSquaresOutDomain/Allclean

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

tutorials/multiphase/overInterDyMFoam/twoSquaresOutDomain/Allrun

@@ -0,0 +1,18 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
+#------------------------------------------------------------------------------
+
+canCompile || exit 0    # Dynamic code
+
+./Allrun.pre
+
+# Serial
+#runApplication $(getApplication)
+
+# Parallel
+runApplication decomposePar -cellDist
+runParallel  $(getApplication)
+runApplication reconstructPar
+
+#------------------------------------------------------------------------------

tutorials/multiphase/overInterDyMFoam/twoSquaresOutDomain/Allrun.pre

@@ -0,0 +1,21 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
+#------------------------------------------------------------------------------
+
+runApplication blockMesh
+
+# Select cellSets
+runApplication -s 1 topoSet
+
+runApplication subsetMesh box -patch hole -overwrite
+
+# Select cellSets
+runApplication -s 2 topoSet
+
+restore0Dir
+
+# Use cellSets to write zoneID
+runApplication setFields
+
+#------------------------------------------------------------------------------

tutorials/multiphase/overInterDyMFoam/twoSquaresOutDomain/constant/dynamicMeshDict

@@ -0,0 +1,86 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2012                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      dynamicMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dynamicFvMesh       dynamicOversetFvMesh;
+
+solver          multiSolidBodyMotionSolver;
+
+multiSolidBodyMotionSolverCoeffs
+{
+    movingZone1
+    {
+        solidBodyMotionFunction multiMotion;//linearMotio;
+        multiMotionCoeffs
+        {
+            one
+            {
+                solidBodyMotionFunction linearMotion;
+                origin      (0.005 0.005 0.005);
+                axis        (0 0 1);
+                omega       0;//100.0;
+                velocity    (-0.1 -0.1 0);
+            }
+            two
+            {
+                solidBodyMotionFunction rotatingMotion;
+                origin      (0.005 0.005 0.005);
+                axis        (0 0 1);
+                omega       20.0;
+                velocity    (-0.1 0 0);
+            }
+        }
+    }
+
+    movingZone2
+    {
+        solidBodyMotionFunction linearMotion;
+        linearMotionCoeffs
+        {
+            velocity    (0.1 0 0);
+            origin      (0.009 0.005 0.005);
+            axis        (0 0 1);
+            omega       -100.0;
+        }
+    }
+}
+
+
+// multiSolidBodyMotionSolverCoeffs
+// {
+//     movingZone1
+//     {
+//         solidBodyMotionFunction rotatingMotion;
+//         rotatingMotionCoeffs
+//         {
+//             origin      (0.005 0.005 0.005);
+//             axis        (0 0 1);
+//             omega       20.0;
+//         }
+//     }
+//
+//     movingZone2
+//     {
+//         solidBodyMotionFunction rotatingMotion;
+//         rotatingMotionCoeffs
+//         {
+//             origin      (0.013 0.005 0.005);
+//             axis        (0 0 1);
+//             omega       -20.0;
+//         }
+//     }
+// }
+
+// ************************************************************************* //

tutorials/multiphase/overInterDyMFoam/twoSquaresOutDomain/constant/fvOptions

@@ -0,0 +1,25 @@
+/*--------------------------------*- C++ -*---------------------------------*\
+|    //\\          A            |                                            |
+|   //   \\        X            | Axion:    GM's End-2-End CFD Solution      |
+|  //      \\      I            | Version:  0.2                              |
+| //         \\    O            | Web:      axion.gm.com                     |
+| ==============   N            |                                            |
+\*--------------------------------------------------------------------------*/
+
+FoamFile
+{
+	object         fvOptions;
+	class          dictionary;
+	format         ascii;
+	version        2.0;
+	location       "system";
+}
+
+limitU
+{
+    type            velocityDampingConstraint;
+    active          true;
+    selectionMode   cellType;
+    UMax        0;
+    C           4;
+}

tutorials/multiphase/overInterDyMFoam/twoSquaresOutDomain/constant/g

@@ -0,0 +1,20 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2012                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       uniformDimensionedVectorField;
+    object      g;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -2 0 0 0 0];
+value           ( 0 -9.81 0);
+
+// ************************************************************************* //

tutorials/multiphase/overInterDyMFoam/twoSquaresOutDomain/constant/transportProperties

@@ -0,0 +1,35 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2012                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+phases (water air);
+
+water
+{
+    transportModel  Newtonian;
+    nu              1e-06;
+    rho             998.2;
+}
+
+air
+{
+    transportModel  Newtonian;
+    nu              1.48e-05;
+    rho             1;
+}
+
+sigma               0.0;
+
+// ************************************************************************* //

tutorials/multiphase/overInterDyMFoam/twoSquaresOutDomain/constant/turbulenceProperties

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

tutorials/multiphase/overInterDyMFoam/twoSquaresOutDomain/system/blockMeshDict

@@ -0,0 +1,148 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2012                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+scale   0.01;
+
+vertices
+(
+    ( 0.00  0.0  0)
+    ( 2.00  0.0  0)
+    ( 2.00  1.0  0)
+    ( 0.00  1.0  0)
+    ( 0.00  0.0  1)
+    ( 2.00  0.0  1)
+    ( 2.00  1.0  1)
+    ( 0.00  1.0  1)
+
+    // movingZone1
+    ( 0.15  0.35  0)
+    ( 0.85  0.35  0)
+    ( 0.85  0.65  0)
+    ( 0.15  0.65  0)
+    ( 0.15  0.35  1)
+    ( 0.85  0.35  1)
+    ( 0.85  0.65  1)
+    ( 0.15  0.65  1)
+
+    // movingZone2
+    ( 1.15  0.15  0)
+    ( 1.45  0.15  0)
+    ( 1.45  0.85  0)
+    ( 1.15  0.85  0)
+    ( 1.15  0.15  1)
+    ( 1.45  0.15  1)
+    ( 1.45  0.85  1)
+    ( 1.15  0.85  1)
+);
+
+blocks
+(
+    hex (0 1 2 3 4 5 6 7) (140 70 1) simpleGrading (1 1 1)
+
+    hex (8 9 10 11 12 13 14 15) movingZone1 (60 24 1) simpleGrading (1 1 1)
+
+    hex (16 17 18 19 20 21 22 23) movingZone2 (24 60 1) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+    overset1
+    {
+        type        overset;
+        faces
+        (
+            ( 8 12 15 11)
+            (10 14 13  9)
+            (11 15 14 10)
+            ( 9 13 12  8)
+        );
+    }
+    overset2
+    {
+        type        overset;
+        faces
+        (
+            (16 20 23 19)
+            (18 22 21 17)
+            (19 23 22 18)
+            (17 21 20 16)
+        );
+    }
+    walls
+    {
+        type wall;
+        faces
+        (
+            (3 7 6 2)
+            (1 5 4 0)
+        );
+    }
+    inlet
+    {
+        type wall;
+        faces
+        (
+            (0 4 7 3)
+        );
+    }
+    outlet
+    {
+        type wall;
+        faces
+        (
+            (2 6 5 1)
+        );
+    }
+    // Populated by subsetMesh
+    hole
+    {
+        type wall;
+        faces ();
+    }
+    frontAndBack
+    {
+        type empty;
+        faces
+        (
+            (0 3 2 1)
+            (4 5 6 7)
+        );
+    }
+    frontAndBack1
+    {
+        type empty;
+        faces
+        (
+            ( 8 11 10  9)
+            (12 13 14 15)
+        );
+    }
+    frontAndBack2
+    {
+        type empty;
+        faces
+        (
+            (16 19 18 17)
+            (20 21 22 23)
+        );
+    }
+);
+
+// ************************************************************************* //

tutorials/multiphase/overInterDyMFoam/twoSquaresOutDomain/system/controlDict

@@ -0,0 +1,116 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2012                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+libs            (overset fvMotionSolvers);
+
+DebugSwitches
+{
+    overset                 0;
+    dynamicOversetFvMesh    0;
+    cellVolumeWeight        0;
+}
+
+application     overInterDyMFoam;
+
+startFrom       startTime;
+
+startTime       0.0;
+
+stopAt          endTime;
+
+endTime         0.08;
+
+deltaT          0.001;
+
+writeControl    adjustable;
+
+writeInterval   0.01;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  12;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable yes;
+
+adjustTimeStep  yes;
+
+maxCo           1.5;
+
+maxAlphaCo      2.0;
+
+maxDeltaT       1;
+
+
+functions
+{
+    probes
+    {
+        type            probes;
+        libs            (sampling);
+
+        // Name of the directory for probe data
+        name            probes;
+
+        // Write at same frequency as fields
+        writeControl    timeStep;
+        writeInterval   1;
+
+        // Fields to be probed
+        fields          (p U);
+
+        // Optional: interpolation scheme to use (default is cell)
+        interpolationScheme cell;
+
+        probeLocations
+        (
+             (0.0009999 0.0015 0.003)
+        );
+    }
+
+    alphaVol
+    {
+        // Mandatory entries (unmodifiable)
+        type            volFieldValue;
+        libs            (fieldFunctionObjects);
+
+        // Mandatory entries (runtime modifiable)
+        fields          (alpha.water);
+        operation       volIntegrate;
+        regionType      all;
+
+        // Optional entries (runtime modifiable)
+        postOperation   none;
+        //weightField     alpha1;
+
+        // Optional (inherited) entries
+        // Write at same frequency as fields
+        writeControl    timeStep;
+        writeInterval   1;
+        writeFields     false;
+        log             true;
+
+    }
+}
+
+// ************************************************************************* //

tutorials/multiphase/overInterDyMFoam/twoSquaresOutDomain/system/decomposeParDict

@@ -0,0 +1,26 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2012                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains  3;
+
+method          hierarchical;
+
+coeffs
+{
+    n           (3 1 1);
+}
+
+// ************************************************************************* //

tutorials/multiphase/overInterDyMFoam/twoSquaresOutDomain/system/fvSchemes

@@ -0,0 +1,82 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2012                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+}
+
+divSchemes
+{
+    default         none;
+
+    div(rhoPhi,U)   Gauss upwind;
+    div(U)          Gauss linear;
+
+    div(phi,alpha)  Gauss vanLeer;
+    div(phirb,alpha) Gauss linear;
+
+    div(phi,alpha.water) Gauss upwind;
+
+    div(phi,epsilon) Gauss upwind;
+    div(phi,k)       Gauss upwind;
+
+    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
+    div((nuEff*dev2(T(grad(U)))))  Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+oversetInterpolation
+{
+    method             trackingInverseDistance;
+    searchBox           (-0.01 -0.01 0)(0.03 0.02 0.01);
+    searchBoxDivisions  3{(264 264 1)};//(264 264 1);//
+    allowInterpolatedDonors false;
+}
+
+fluxRequired
+{
+    default         no;
+    pcorr           ;
+    p               ;
+}
+
+oversetInterpolationSuppressed
+{
+    grad(p_rgh);
+    surfaceIntegrate(phiHbyA);
+}
+
+// ************************************************************************* //

tutorials/multiphase/overInterDyMFoam/twoSquaresOutDomain/system/fvSolution

@@ -0,0 +1,115 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2012                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    cellDisplacement
+    {
+        solver          PCG;
+        preconditioner  DIC;
+
+        tolerance       1e-06;
+        relTol          0;
+        maxIter         100;
+    }
+
+    "alpha.water.*"
+    {
+        nAlphaCorr      3;
+        nAlphaSubCycles 2;
+        cAlpha          1;
+        icAlpha         0;
+        scAlpha         0;
+
+        MULESCorr       yes;
+        nLimiterIter    5;
+        alphaApplyPrevCorr  no;
+
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-8;
+        relTol          0;
+    }
+
+    p_rgh
+    {
+        solver          PBiCGStab;
+        preconditioner  DILU;
+        tolerance       1e-12;
+        relTol          0.01;
+        maxIter         500;
+    }
+
+    p_rghFinal
+    {
+        $p_rgh;
+        relTol          0;
+    }
+
+    pcorr
+    {
+        $p;
+        solver          PCG;
+        preconditioner  DIC;
+    }
+
+    pcorrFinal
+    {
+        $pcorr;
+        relTol          0;
+    }
+
+    "(U|k|epsilon)"
+    {
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-6;
+        relTol          0.01;
+        maxIter         200;
+        minIter         1;
+    }
+
+    "(U|k|epsilon)Final"
+    {
+        $U;
+        relTol          0;
+    }
+}
+
+PIMPLE
+{
+    momentumPredictor   no;
+    correctPhi          no;
+    nOuterCorrectors    2;
+    nCorrectors         3;
+    nNonOrthogonalCorrectors 0;
+
+    pRefPoint           (0.0002 0.0099 0.001);
+    pRefValue           0.0;
+}
+
+relaxationFactors
+{
+   fields
+   {
+   }
+   equations
+   {
+       ".*"             1;
+   }
+}
+
+// ************************************************************************* //

tutorials/multiphase/overInterDyMFoam/twoSquaresOutDomain/system/setFieldsDict

@@ -0,0 +1,67 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2012                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      setFieldsDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+defaultFieldValues
+(
+    volScalarFieldValue zoneID 123
+    volScalarFieldValue alpha.water 0
+);
+
+regions
+(
+    // Set cell values
+    // (does zerogradient on boundaries)
+    cellToCell
+    {
+        set c0;
+
+        fieldValues
+        (
+            volScalarFieldValue zoneID 0
+        );
+    }
+
+    cellToCell
+    {
+        set c1;
+
+        fieldValues
+        (
+            volScalarFieldValue zoneID 1
+        );
+    }
+
+    cellToCell
+    {
+        set c2;
+
+        fieldValues
+        (
+            volScalarFieldValue zoneID 2
+        );
+    }
+
+    boxToCell
+    {
+        box ( -100 -100 -100 ) ( 100 0.005 100 );
+        fieldValues
+        (
+            volScalarFieldValue alpha.water 1
+        );
+    }
+);
+
+// ************************************************************************* //

tutorials/multiphase/overInterDyMFoam/twoSquaresOutDomain/system/topoSetDict

@@ -0,0 +1,94 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2012                                 |
+|   \\  /    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  regionsToCell;
+        insidePoints ((0.001 0.001 0.001));
+    }
+
+    {
+        name    c1;
+        type    cellSet;
+        action  new;
+        source  cellToCell;
+        set     c0;
+    }
+
+    {
+        name    c1;
+        type    cellSet;
+        action  invert;
+    }
+
+    {
+        name    c2;
+        type    cellSet;
+        action  new;
+        source  regionsToCell;
+        insidePoints ((0.0116 0.00151 0.001));
+        set     c1;
+    }
+
+    {
+        name    c1;
+        type    cellSet;
+        action  subtract;
+        source  cellToCell;
+        set     c2;
+    }
+
+    // Select box to remove from region 1 and 2
+    {
+        name    box;
+        type    cellSet;
+        action  new;
+        source  cellToCell;
+        set     c1;
+    }
+
+    {
+        name    box;
+        type    cellSet;
+        action  add;
+        source  cellToCell;
+        set     c2;
+    }
+
+    {
+        name    box;
+        type    cellSet;
+        action  subset;
+        source  boxToCell;
+        boxes
+        (
+            (0.0025 0.0045 -100)(0.0075 0.0055 100)
+            (0.0125 0.0025 -100)(0.0135 0.0075 100)
+        );
+    }
+
+    {
+        name    box;
+        type    cellSet;
+        action  invert;
+    }
+);
+
+// ************************************************************************* //

tutorials/multiphase/overInterPhaseChangeDyMFoam/twoSimpleRotors/system/controlDict

@@ -31,7 +31,7 @@ startTime       0.0;
 
 stopAt          endTime;
 
-endTime         0.2;
+endTime         0.1;
 
 deltaT          0.00001;
 
@@ -87,24 +87,6 @@ functions
              (0.0009999 0.0015 0.003)
         );
     }
-
-    alphaVol
-    {
-        libs            (utilityFunctionObjects);
-        type            coded;
-        name            alphaVolume;
-        writeControl    timeStep;
-        writeInterval   10;
-
-        codeWrite
-        #{
-            const auto& alpha =
-                mesh().lookupObject<volScalarField>("alpha.water");
-
-            Info<< "Alpha volume = " << alpha.weightedAverage(mesh().Vsc())
-                << endl;
-        #};
-    }
 }
 
 // ************************************************************************* //

tutorials/multiphase/overInterPhaseChangeDyMFoam/twoSimpleRotors/system/fvSolution

@@ -63,8 +63,8 @@ solvers
     pcorr
     {
         $p;
-        solver          PCG;
+        solver          PBiCGStab;
-        preconditioner  DIC;
+        preconditioner  DILU;
     }
 
     pcorrFinal
@@ -98,10 +98,8 @@ PIMPLE
     nCorrectors         3;
     nNonOrthogonalCorrectors 0;
 
-    ddtCorr             false;
-
     pRefPoint           (0.0002 0.0099 0.001);
-    pRefValue           1e5; //Not used
+    pRefValue           1e5;
 }
 
 relaxationFactors