Cumulative development by Vuko Vukcevic.

applications/solvers/overset/icoDyMOversetFoam/Make/files

@@ -0,0 +1,3 @@
+icoDyMOversetFoam.C  
+
+EXE = $(FOAM_APPBIN)/icoDyMOversetFoam

applications/solvers/overset/icoDyMOversetFoam/Make/options

@@ -0,0 +1,14 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/overset/oversetMesh/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools \
+    -ldynamicFvMesh \
+    -ldynamicMesh \
+    -loversetMesh \
+    -llduSolvers

applications/solvers/overset/icoDyMOversetFoam/UEqn.H

@@ -0,0 +1,12 @@
+// Momentum equation
+fvVectorMatrix UEqn
+(
+    fvm::ddt(U)
+  + fvm::div(phi, U)
+  - fvm::laplacian(nu, U)
+);
+
+if (piso.momentumPredictor())
+{
+    solve(UEqn == -fvc::grad(p));
+}

applications/solvers/overset/icoDyMOversetFoam/correctPhi.H

@@ -0,0 +1,57 @@
+{
+    wordList pcorrTypes(p.boundaryField().types());
+
+    for (label i = 0; i < p.boundaryField().size(); i++)
+    {
+        if (p.boundaryField()[i].fixesValue())
+        {
+            pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
+        }
+    }
+
+    volScalarField pcorr
+    (
+        IOobject
+        (
+            "pcorr",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("pcorr", p.dimensions(), 0),
+        pcorrTypes
+    );
+
+    // Initialise flux with interpolated velocity
+    phi = faceOversetMask*(fvc::interpolate(U) & mesh.Sf());
+
+    // Adjust overset fluxes
+    oversetAdjustPhi(phi, U); // Fringe flux adjustment
+    globalOversetAdjustPhi(phi, U, pcorr); // Global flux adjustment
+
+    while (piso.correctNonOrthogonal())
+    {
+        fvScalarMatrix pcorrEqn
+        (
+            fvm::laplacian(rAU, pcorr) == fvc::div(phi)
+        );
+
+        // Adjust non-orthogonal fluxes if necessary
+        om.correctNonOrthoFluxes(pcorrEqn, U);
+
+        pcorrEqn.setReference(pRefCell, pRefValue);
+        pcorrEqn.solve();
+
+        if (piso.finalNonOrthogonalIter())
+        {
+            phi -= pcorrEqn.flux();
+        }
+
+        // Fluxes are corrected to absolute velocity and further corrected
+        // later.  HJ, 6/Feb/2009
+    }
+
+#   include "oversetContinuityErrs.H"
+}

applications/solvers/overset/icoDyMOversetFoam/createControls.H

@@ -0,0 +1,11 @@
+#include "createTimeControls.H"
+
+bool correctPhi
+(
+    piso.dict().lookupOrDefault("correctPhi", false)
+);
+
+bool checkMeshCourantNo
+(
+    piso.dict().lookupOrDefault("checkMeshCourantNo", false)
+);

applications/solvers/overset/icoDyMOversetFoam/createFields.H

@@ -0,0 +1,79 @@
+    Info<< "Reading transportProperties\n" << endl;
+
+    IOdictionary transportProperties
+    (
+        IOobject
+        (
+            "transportProperties",
+            runTime.constant(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    );
+
+    dimensionedScalar nu
+    (
+        transportProperties.lookup("nu")
+    );
+
+    Info<< "Reading field p\n" << endl;
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "Reading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+#   include "createPhi.H"
+
+    // Set reference cell taking into account whether implicit or explicit
+    // overset aglorithm is used for field p.
+    const oversetMesh& om = oversetMesh::New(mesh);
+    label pRefCell = 0;
+    scalar pRefValue = 0.0;
+    om.setRefCell(p, piso.dict(), pRefCell, pRefValue);
+
+    mesh.schemesDict().setFluxRequired(p.name());
+    if (correctPhi)
+    {
+        mesh.schemesDict().setFluxRequired("pcorr");
+    }
+
+    Info<< "Reading field rAU if present\n" << endl;
+    volScalarField rAU
+    (
+        IOobject
+        (
+            "rAU",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("dt", dimTime, 1.0),
+        zeroGradientFvPatchScalarField::typeName
+    );
+
+#   include "writeOversetMasks.H"

applications/solvers/overset/icoDyMOversetFoam/icoDyMOversetFoam.C

@@ -0,0 +1,167 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    icoDyMOversetFoam
+
+Description
+    Transient solver for incompressible, laminar flow of Newtonian fluids
+    with dynamic mesh and overset mesh support.
+
+Author
+    Hrvoje Jasak, Wikki Ltd.  All rights reserved
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "dynamicFvMesh.H"
+#include "oversetMesh.H"
+#include "oversetFvPatchFields.H"
+#include "oversetAdjustPhi.H"
+#include "globalOversetAdjustPhi.H"
+#include "pisoControl.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+#   include "setRootCase.H"
+#   include "createTime.H"
+#   include "createDynamicFvMesh.H"
+
+    pisoControl piso(mesh);
+
+#   include "initContinuityErrs.H"
+#   include "initTotalVolume.H"
+#   include "createControls.H"
+#   include "createFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+#       include "readControls.H"
+#       include "checkTotalVolume.H"
+
+        // Make the fluxes absolute
+        fvc::makeAbsolute(phi, U);
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        bool meshChanged = mesh.update();
+        reduce(meshChanged, orOp<bool>());
+
+#       include "createOversetMasks.H"
+
+        if (correctPhi && (mesh.moving() || meshChanged))
+        {
+#           include "correctPhi.H"
+        }
+
+        // Make the fluxes relative to the mesh motion
+        fvc::makeRelative(phi, U);
+
+#       include "oversetCourantNo.H"
+#       include "setDeltaT.H"
+
+        if (mesh.moving() && checkMeshCourantNo)
+        {
+#           include "meshCourantNo.H"
+        }
+
+#       include "UEqn.H"
+
+        // --- PISO loop
+
+        while (piso.correct())
+        {
+            p.boundaryField().updateCoeffs();
+
+            rAU = 1.0/UEqn.A();
+            oversetFvPatchScalarField::oversetInterpolate(rAU); // Overset update
+
+            U = rAU*UEqn.H();
+            oversetFvPatchVectorField::oversetInterpolate(U); // Overset update
+
+            phi = fvc::interpolate(U) & mesh.Sf();
+
+            // Adjust overset fluxes
+            oversetAdjustPhi(phi, U); // Fringe flux adjustment
+            globalOversetAdjustPhi(phi, U, p); // Global flux adjustment
+
+            while (piso.correctNonOrthogonal())
+            {
+                fvScalarMatrix pEqn
+                (
+                    fvm::laplacian(rAU, p) == fvc::div(phi)
+                );
+
+                // Adjust non-orthogonal fringe fluxes if necessary
+                om.correctNonOrthoFluxes(pEqn, U);
+
+                pEqn.setReference(pRefCell, pRefValue);
+                pEqn.solve
+                (
+                    mesh.solutionDict().solver(p.select(piso.finalInnerIter()))
+                );
+
+                if (piso.finalNonOrthogonalIter())
+                {
+                    phi -= pEqn.flux();
+                }
+
+                // Perform overset interpolation (after flux reconstruction)
+                oversetFvPatchScalarField::oversetInterpolate(p);
+            }
+
+#           include "oversetContinuityErrs.H"
+
+            U -= rAU*fvc::grad(p);
+            U.correctBoundaryConditions();
+            oversetFvPatchVectorField::oversetInterpolate(U); // Overset update
+            // Note: if implicit fringe is switched on, we are doing the
+            // interpolation twice (once in correctBoundaryConditions and once
+            // in oversetInterpolate). Reorganize. VV, 4/Oct/2016.
+
+            // Make the fluxes relative to the mesh motion
+            fvc::makeRelative(phi, U);
+        }
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return(0);
+}
+
+
+// ************************************************************************* //

applications/solvers/overset/icoDyMOversetFoam/readControls.H

@@ -0,0 +1,5 @@
+#include "readTimeControls.H"
+
+correctPhi = piso.dict().lookupOrDefault("correctPhi", false);
+
+checkMeshCourantNo = piso.dict().lookupOrDefault("checkMeshCourantNo", false);

applications/solvers/overset/icoOversetFoam/Make/files

@@ -0,0 +1,3 @@
+icoOversetFoam.C  
+
+EXE = $(FOAM_APPBIN)/icoOversetFoam

applications/solvers/overset/icoOversetFoam/Make/options

@@ -0,0 +1,10 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/triSurface/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/overset/oversetMesh/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -loversetMesh \
+    -llduSolvers

applications/solvers/overset/icoOversetFoam/createFields.H

@@ -0,0 +1,74 @@
+    Info<< "Reading transportProperties\n" << endl;
+
+    IOdictionary transportProperties
+    (
+        IOobject
+        (
+            "transportProperties",
+            runTime.constant(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    );
+
+    dimensionedScalar nu
+    (
+        transportProperties.lookup("nu")
+    );
+
+    Info<< "Reading field p\n" << endl;
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "Reading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+#   include "createPhi.H"
+
+    // Set reference cell taking into account whether implicit or explicit
+    // overset aglorithm is used for field p.
+    label pRefCell = 0;
+    scalar pRefValue = 0.0;
+    om.setRefCell(p, piso.dict(), pRefCell, pRefValue);
+
+    mesh.schemesDict().setFluxRequired(p.name());
+
+    Info<< "Reading field rAU if present\n" << endl;
+    volScalarField rAU
+    (
+        IOobject
+        (
+            "rAU",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("dt", dimTime, 1.0),
+        zeroGradientFvPatchScalarField::typeName
+    );
+
+#   include "writeOversetMasks.H"

applications/solvers/overset/icoOversetFoam/icoOversetFoam.C

@@ -0,0 +1,150 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    icoOversetFoam
+
+Description
+    Transient solver for incompressible, laminar flow of Newtonian fluids
+    with overset mesh support.
+
+Author
+    Hrvoje Jasak, Wikki Ltd.  All rights reserved
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "pisoControl.H"
+
+#include "oversetMesh.H"
+#include "oversetFvPatchFields.H"
+#include "oversetAdjustPhi.H"
+#include "globalOversetAdjustPhi.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+#   include "setRootCase.H"
+#   include "createTime.H"
+#   include "createMesh.H"
+
+    pisoControl piso(mesh);
+
+#   include "createOversetMasks.H"
+#   include "createFields.H"
+#   include "initContinuityErrs.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.loop())
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+#       include "oversetCourantNo.H"
+
+        // Momentum equation
+        fvVectorMatrix UEqn
+        (
+            fvm::ddt(U)
+          + fvm::div(phi, U)
+          - fvm::laplacian(nu, U)
+        );
+
+        if (piso.momentumPredictor())
+        {
+            solve(UEqn == -fvc::grad(p));
+        }
+
+        // --- PISO loop
+
+        while (piso.correct())
+        {
+            p.boundaryField().updateCoeffs();
+
+            rAU = 1.0/UEqn.A();
+            oversetFvPatchScalarField::oversetInterpolate(rAU); // Overset update
+
+            U = rAU*UEqn.H();
+            oversetFvPatchVectorField::oversetInterpolate(U); // Overset update
+
+            phi = fvc::interpolate(U) & mesh.Sf();
+
+            // Adjust overset fluxes
+            oversetAdjustPhi(phi, U); // Fringe flux adjustment
+            globalOversetAdjustPhi(phi, U, p); // Global flux adjustment
+
+            while (piso.correctNonOrthogonal())
+            {
+                fvScalarMatrix pEqn
+                (
+                    fvm::laplacian(rAU, p) == fvc::div(phi)
+                );
+
+                // Adjust non-orthogonal fringe fluxes if necessary
+                om.correctNonOrthoFluxes(pEqn, U);
+
+                pEqn.setReference(pRefCell, pRefValue);
+                pEqn.solve
+                (
+                    mesh.solutionDict().solver(p.select(piso.finalInnerIter()))
+                );
+
+                if (piso.finalNonOrthogonalIter())
+                {
+                    phi -= pEqn.flux();
+                }
+
+                // Perform overset interpolation (after flux reconstruction)
+                oversetFvPatchScalarField::oversetInterpolate(p);
+            }
+
+#           include "oversetContinuityErrs.H"
+
+            U -= rAU*fvc::grad(p);
+            U.correctBoundaryConditions();
+            oversetFvPatchVectorField::oversetInterpolate(U); // Overset update
+            // Note: if implicit fringe is switched on, we are doing the
+            // interpolation twice (once in correctBoundaryConditions and once
+            // in oversetInterpolate). Reorganize. VV, 4/Oct/2016.
+
+            // Make the fluxes relative to the mesh motion
+            fvc::makeRelative(phi, U);
+        }
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return(0);
+}
+
+
+// ************************************************************************* //

applications/solvers/overset/interOversetFoam/Make/files

@@ -0,0 +1,3 @@
+interOversetFoam.C
+
+EXE = $(FOAM_APPBIN)/interOversetFoam

applications/solvers/overset/interOversetFoam/Make/options

@@ -0,0 +1,21 @@
+EXE_INC = \
+    -I$(LIB_SRC)/turbulenceModels \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/triSurface/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/overset/oversetMesh/lnInclude
+
+EXE_LIBS = \
+    -linterfaceProperties \
+    -lincompressibleTransportModels \
+    -lincompressibleTurbulenceModel \
+    -lincompressibleRASModels \
+    -lincompressibleLESModels \
+    -lfiniteVolume \
+    -loversetMesh \
+    -llduSolvers

applications/solvers/overset/interOversetFoam/UEqn.H

@@ -0,0 +1,33 @@
+    surfaceScalarField muEff
+    (
+        "muEff",
+        twoPhaseProperties.muf()
+      + fvc::interpolate(rho*turbulence->nut())
+    );
+
+    fvVectorMatrix UEqn
+    (
+        fvm::ddt(rho, U)
+      + fvm::div(rhoPhi, U)
+      - fvm::laplacian(muEff, U)
+      - (fvc::grad(U) & fvc::grad(muEff))
+    );
+
+    UEqn.relax();
+
+    if (pimple.momentumPredictor())
+    {
+        solve
+        (
+            UEqn
+         ==
+            fvc::reconstruct
+            (
+                (
+                    fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
+                  - ghf*fvc::snGrad(rho)
+                  - fvc::snGrad(pd)
+                ) * mesh.magSf()
+            )
+        );
+    }

applications/solvers/overset/interOversetFoam/alphaEqn.H

@@ -0,0 +1,60 @@
+{
+    word alphaScheme("div(phi,alpha)");
+    word alpharScheme("div(phirb,alpha)");
+
+    // New formulation of phir: Co condition
+    // Note faceOversetMask.  HJ, 16/Jun/2015
+    surfaceScalarField phir
+    (
+        "phir",
+        faceOversetMask*interface.cAlpha()*interface.nHatf()*
+        min
+        (
+            scalar(0.5)/    // This is ref compression Co number for cAlpha = 1
+            (
+                mesh.surfaceInterpolation::deltaCoeffs()*
+                mesh.time().deltaT()
+            ),
+            mag(phi/mesh.magSf())
+          + dimensionedScalar("small", dimVelocity, 1e-3)
+        )
+    );
+
+    fvScalarMatrix alpha1Eqn
+    (
+        fvm::ddt(alpha1)
+      + fvm::div(phi, alpha1, alphaScheme)
+      + fvm::div
+        (
+            -fvc::flux(-phir, scalar(1) - alpha1, alpharScheme),
+            alpha1,
+            alpharScheme
+        )
+    );
+
+    alpha1Eqn.relax();
+    alpha1Eqn.solve();
+
+    // Update mass flux
+    rhoPhi = alpha1Eqn.flux()*(rho1 - rho2) + phi*rho2;
+
+    Info<< "Liquid phase volume fraction = "
+        << alpha1.weightedAverage(mesh.V()).value()
+        << "  Min(alpha1) = " << min(cellOversetMask*alpha1).value()
+        << "  Max(alpha1) = " << max(cellOversetMask*alpha1).value()
+        << endl;
+
+    // Limit alpha to handle overset cells
+    alpha1.max(0);
+    alpha1.min(1);
+
+    // Update of interface and density
+    interface.correct();
+    twoPhaseProperties.correct();
+
+    // Calculate density using limited alpha1
+    rho = twoPhaseProperties.rho();
+
+    // Parallel update in snGrad.  HJ, 8/Dec/2012
+    rho.correctBoundaryConditions();
+}

applications/solvers/overset/interOversetFoam/correctPhi.H

@@ -0,0 +1,60 @@
+{
+    // Set flux required
+    mesh.schemesDict().setFluxRequired("pcorr");
+
+    wordList pcorrTypes(p.boundaryField().types());
+
+    for (label i = 0; i < p.boundaryField().size(); i++)
+    {
+        if (p.boundaryField()[i].fixesValue())
+        {
+            pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
+        }
+    }
+
+    volScalarField pcorr
+    (
+        IOobject
+        (
+            "pcorr",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("pcorr", p.dimensions(), 0),
+        pcorrTypes
+    );
+
+    // Initialise flux with interpolated velocity
+    phi = faceOversetMask*(fvc::interpolate(U) & mesh.Sf());
+
+    // Adjust overset fluxes
+    oversetAdjustPhi(phi, U); // Fringe flux adjustment
+    globalOversetAdjustPhi(phi, U, pcorr); // Global flux adjustment
+
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pcorrEqn
+        (
+            fvm::laplacian(rAU, pcorr) == fvc::div(phi)
+        );
+
+        // Adjust non-orthogonal fluxes if necessary
+        om.correctNonOrthoFluxes(pcorrEqn, U);
+
+        pcorrEqn.setReference(pRefCell, pRefValue);
+        pcorrEqn.solve();
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi -= pcorrEqn.flux();
+        }
+
+        // Fluxes are corrected to absolute velocity and further corrected
+        // later.  HJ, 6/Feb/2009
+    }
+
+#   include "oversetContinuityErrs.H"
+}

applications/solvers/overset/interOversetFoam/createFields.H

@@ -0,0 +1,137 @@
+    Info<< "Reading field pd\n" << endl;
+    volScalarField pd
+    (
+        IOobject
+        (
+            "pd",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "Reading field alpha1\n" << endl;
+    volScalarField alpha1
+    (
+        IOobject
+        (
+            "alpha1",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "Reading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+#   include "createPhi.H"
+
+
+    Info<< "Reading transportProperties\n" << endl;
+    twoPhaseMixture twoPhaseProperties(U, phi, "alpha1");
+
+    const dimensionedScalar& rho1 = twoPhaseProperties.rho1();
+    const dimensionedScalar& rho2 = twoPhaseProperties.rho2();
+
+
+    // Need to store rho for ddt(rho, U)
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT
+        ),
+        alpha1*rho1 + (scalar(1) - alpha1)*rho2,
+        alpha1.boundaryField().types()
+    );
+    rho.oldTime();
+
+
+    // Mass flux
+    // Initialisation does not matter because rhoPhi is reset after the
+    // alpha1 solution before it is used in the U equation.
+    surfaceScalarField rhoPhi
+    (
+        IOobject
+        (
+            "rho*phi",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        rho1*phi
+    );
+
+
+    Info<< "Calculating field g.h\n" << endl;
+    volScalarField gh("gh", g & mesh.C());
+    surfaceScalarField ghf("gh", g & mesh.Cf());
+
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        pd + rho*gh
+    );
+
+    // Set reference cell taking into account whether implicit or explicit
+    // overset aglorithm is used for field p.
+    label pRefCell = 0;
+    scalar pRefValue = 0.0;
+    om.setRefCell(p, pimple.dict(), pRefCell, pRefValue);
+
+    mesh.schemesDict().setFluxRequired(pd.name());
+    mesh.schemesDict().setFluxRequired(alpha1.name());
+
+    // Construct interface from alpha1 distribution
+    interfaceProperties interface(alpha1, U, twoPhaseProperties);
+
+     // Construct incompressible turbulence model
+    autoPtr<incompressible::turbulenceModel> turbulence
+    (
+        incompressible::turbulenceModel::New(U, phi, twoPhaseProperties)
+    );
+
+    Info<< "Reading field rAU if present\n" << endl;
+    volScalarField rAU
+    (
+        IOobject
+        (
+            "rAU",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        runTime.deltaT()/rho1,
+        zeroGradientFvPatchScalarField::typeName
+    );
+
+#   include "writeOversetMasks.H"

applications/solvers/overset/interOversetFoam/interOversetFoam.C

@@ -0,0 +1,125 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    interOversetFoam
+
+Description
+    Solver for 2 incompressible, isothermal immiscible fluids using a VOF
+    (volume of fluid) phase-fraction based interface capturing approach,
+    with overset mesh support
+
+    The momentum and other fluid properties are of the "mixture" and a single
+    momentum equation is solved.
+
+    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
+
+    For a two-fluid approach see twoPhaseEulerFoam.
+
+Author
+    Hrvoje Jasak, Wikki Ltd.  All rights reserved
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "interfaceProperties.H"
+#include "twoPhaseMixture.H"
+#include "turbulenceModel.H"
+#include "pimpleControl.H"
+
+#include "oversetMesh.H"
+#include "oversetFvPatchFields.H"
+#include "oversetAdjustPhi.H"
+#include "globalOversetAdjustPhi.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+#   include "setRootCase.H"
+#   include "createTime.H"
+#   include "createMesh.H"
+
+    pimpleControl pimple(mesh);
+
+#   include "readGravitationalAcceleration.H"
+#   include "initContinuityErrs.H"
+#   include "createOversetMasks.H"
+#   include "createFields.H"
+#   include "createTimeControls.H"
+#   include "correctPhi.H"
+#   include "CourantNo.H"
+#   include "setInitialDeltaT.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+#       include "readTimeControls.H"
+#       include "oversetCourantNo.H"
+#       include "setDeltaT.H"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // Pressure-velocity corrector
+        while (pimple.loop())
+        {
+            twoPhaseProperties.correct();
+
+#           include "alphaEqn.H"
+
+#           include "UEqn.H"
+
+            // --- PISO loop
+            while (pimple.correct())
+            {
+#               include "pdEqn.H"
+            }
+
+#           include "oversetContinuityErrs.H"
+
+#           include "limitU.H"
+
+            p = pd + cellOversetMask*rho*gh;
+
+            turbulence->correct();
+        }
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/overset/interOversetFoam/limitU.H

@@ -0,0 +1,29 @@
+{
+    scalar limitMagU = readScalar(pimple.dict().lookup("limitMagU"));
+
+    volScalarField magU(mag(U));
+
+    scalar maxMagU = max(magU).value();
+
+    Info<< "mag(U): max: " << maxMagU
+        << " avg: " << magU.weightedAverage(mesh.V()).value();
+
+    if (maxMagU > limitMagU)
+    {
+        U.internalField() *=
+            neg(magU.internalField() - limitMagU)
+          + pos(magU.internalField() - limitMagU)*
+            limitMagU/(magU.internalField() + SMALL);
+        U.correctBoundaryConditions();
+        oversetFvPatchVectorField::oversetInterpolate(U); // Overset update
+        // Note: if implicit fringe is switched on, we are doing the
+        // interpolation twice (once in correctBoundaryConditions and once in
+        // oversetInterpolate). Reorganize. VV, 4/Oct/2016.
+
+        Info << " ...limiting" << endl;
+    }
+    else
+    {
+        Info << endl;
+    }
+}

applications/solvers/overset/interOversetFoam/pdEqn.H

@@ -0,0 +1,65 @@
+{
+    pd.boundaryField().updateCoeffs();
+
+    rAU = 1.0/UEqn.A();
+    oversetFvPatchScalarField::oversetInterpolate(rAU); // Overset update
+    surfaceScalarField rAUf = fvc::interpolate(rAU);
+
+    U = rAU*UEqn.H();
+    oversetFvPatchVectorField::oversetInterpolate(U); // Overset update
+
+    surfaceScalarField phiU
+    (
+        "phiU",
+        fvc::interpolate(U) & mesh.Sf()
+    );
+
+    phi = phiU
+      + (
+            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
+          - ghf*fvc::snGrad(rho)
+        )*rAUf*mesh.magSf();
+
+    // Adjust overset fluxes
+    oversetAdjustPhi(phi, U); // Fringe flux adjustment
+    globalOversetAdjustPhi(phi, U, p); // Global flux adjustment
+
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pdEqn
+        (
+            fvm::laplacian(rAUf, pd) == fvc::div(phi)
+        );
+
+        // Adjust non-orthogonal fringe fluxes if necessary
+        om.correctNonOrthoFluxes(pdEqn, U);
+
+        pdEqn.setReference(pRefCell, pRefValue);
+
+        pdEqn.solve
+        (
+            mesh.solutionDict().solver(pd.select(pimple.finalInnerIter()))
+        );
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi -= pdEqn.flux();
+        }
+
+        // Perform overset interpolation (after flux reconstruction)
+        oversetFvPatchScalarField::oversetInterpolate(p);
+    }
+
+    // Explicitly relax pressure for momentum corrector
+    if (!pimple.finalIter())
+    {
+        p.relax();
+    }
+
+    U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
+    U.correctBoundaryConditions();
+    oversetFvPatchVectorField::oversetInterpolate(U); // Overset update
+    // Note: if implicit fringe is switched on, we are doing the interpolation
+    // twice (once in correctBoundaryConditions and once in oversetInterpolate)
+    // Reorganize. VV, 4/Oct/2016.
+}

applications/solvers/overset/laplacianOversetFoam/Make/files

@@ -0,0 +1,3 @@
+laplacianOversetFoam.C
+
+EXE = $(FOAM_APPBIN)/laplacianOversetFoam

applications/solvers/overset/laplacianOversetFoam/Make/options

@@ -0,0 +1,8 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/overset/oversetMesh/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -loversetMesh \
+    -llduSolvers

applications/solvers/overset/laplacianOversetFoam/createFields.H

@@ -0,0 +1,37 @@
+    Info<< "Reading field T\n" << endl;
+
+    volScalarField T
+    (
+        IOobject
+        (
+            "T",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+
+    Info<< "Reading transportProperties\n" << endl;
+
+    IOdictionary transportProperties
+    (
+        IOobject
+        (
+            "transportProperties",
+            runTime.constant(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    );
+
+
+    Info<< "Reading diffusivity DT\n" << endl;
+
+    dimensionedScalar DT
+    (
+        transportProperties.lookup("DT")
+    );

applications/solvers/overset/laplacianOversetFoam/laplacianOversetFoam.C

@@ -0,0 +1,99 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    laplaciantOversetFoam
+
+Description
+    Solves a simple Laplace equation, e.g. for thermal diffusion in a solid
+    with support for overset meshes.
+
+    Experimental: no overset mesh changes required
+
+Author
+    Hrvoje Jasak, Wikki Ltd.  All rights reserved
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "simpleControl.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+
+#   include "setRootCase.H"
+
+#   include "createTime.H"
+#   include "createMesh.H"
+
+    simpleControl simple(mesh);
+
+#   include "createFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nCalculating temperature distribution\n" << endl;
+
+    while (simple.loop())
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        T.correctBoundaryConditions();
+        Pout<< "T: " << T.internalField() << endl;
+        return 0;
+
+        fvScalarMatrix TEqn
+        (
+            fvm::laplacian(DT, T)
+        );
+        Pout<< "TEqn: " << TEqn << endl;
+//         TEqn.solve();
+
+        volScalarField residual
+        (
+            "residual",
+            T
+        );
+
+        residual.internalField() = TEqn.residual();
+//         residual.boundaryField() == 0;
+        residual.write();
+
+        Info<< "residual " << gSumMag(residual.internalField()) << endl;
+
+#       include "write.H"
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/overset/laplacianOversetFoam/write.H

@@ -0,0 +1,46 @@
+    if (runTime.outputTime())
+    {
+        volVectorField gradT = fvc::grad(T);
+
+        volScalarField gradTx
+        (
+            IOobject
+            (
+                "gradTx",
+                runTime.timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            gradT.component(vector::X)
+        );
+
+        volScalarField gradTy
+        (
+            IOobject
+            (
+                "gradTy",
+                runTime.timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            gradT.component(vector::Y)
+        );
+
+        volScalarField gradTz
+        (
+            IOobject
+            (
+                "gradTz",
+                runTime.timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            gradT.component(vector::Z)
+        );
+
+
+        runTime.write();
+    }

applications/solvers/overset/pimpleDyMOversetFoam/Make/files

@@ -0,0 +1,3 @@
+pimpleDyMOversetFoam.C
+
+EXE = $(FOAM_APPBIN)/pimpleDyMOversetFoam

applications/solvers/overset/pimpleDyMOversetFoam/Make/options

@@ -0,0 +1,25 @@
+EXE_INC = \
+    -I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/overset/oversetMesh/lnInclude \
+    -I$(LIB_SRC)/overset/oversetDynamicFvMesh/lnInclude
+
+EXE_LIBS = \
+    -ldynamicFvMesh \
+    -ltopoChangerFvMesh \
+    -ldynamicMesh \
+    -lmeshTools \
+    -lincompressibleTransportModels \
+    -lincompressibleTurbulenceModel \
+    -lincompressibleRASModels \
+    -lincompressibleLESModels \
+    -lfiniteVolume \
+    -llduSolvers \
+    -L$(MESQUITE_LIB_DIR) -lmesquite \
+    -loversetMesh \
+    -loversetDynamicFvMesh

applications/solvers/overset/pimpleDyMOversetFoam/UEqn.H

@@ -0,0 +1,19 @@
+    fvVectorMatrix UEqn
+    (
+        fvm::ddt(U)
+      + fvm::div(phi, U)
+      + turbulence->divDevReff()
+    );
+
+    UEqn.relax
+    (
+        mesh.solutionDict().equationRelaxationFactor
+        (
+            U.select(pimple.finalIter())
+        )
+    );
+
+    if (pimple.momentumPredictor())
+    {
+        solve(UEqn == -fvc::grad(p));
+    }

applications/solvers/overset/pimpleDyMOversetFoam/correctPhi.H

@@ -0,0 +1,57 @@
+{
+    wordList pcorrTypes(p.boundaryField().types());
+
+    for (label i = 0; i < p.boundaryField().size(); i++)
+    {
+        if (p.boundaryField()[i].fixesValue())
+        {
+            pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
+        }
+    }
+
+    volScalarField pcorr
+    (
+        IOobject
+        (
+            "pcorr",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("pcorr", p.dimensions(), 0),
+        pcorrTypes
+    );
+
+    // Initialise flux with interpolated velocity
+    phi = faceOversetMask*(fvc::interpolate(U) & mesh.Sf());
+
+    // Adjust overset fluxes
+    oversetAdjustPhi(phi, U); // Fringe flux adjustment
+    globalOversetAdjustPhi(phi, U, pcorr); // Global flux adjustment
+
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pcorrEqn
+        (
+            fvm::laplacian(rAU, pcorr) == fvc::div(phi)
+        );
+
+        // Adjust non-orthogonal fluxes if necessary
+        om.correctNonOrthoFluxes(pcorrEqn, U);
+
+        pcorrEqn.setReference(pRefCell, pRefValue);
+        pcorrEqn.solve();
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi -= pcorrEqn.flux();
+        }
+
+        // Fluxes are corrected to absolute velocity and further corrected
+        // later.  HJ, 6/Feb/2009
+    }
+
+#   include "oversetContinuityErrs.H"
+}

applications/solvers/overset/pimpleDyMOversetFoam/createControls.H

@@ -0,0 +1,11 @@
+#include "createTimeControls.H"
+
+bool correctPhi
+(
+    pimple.dict().lookupOrDefault("correctPhi", false)
+);
+
+bool checkMeshCourantNo
+(
+    pimple.dict().lookupOrDefault("checkMeshCourantNo", false)
+);

applications/solvers/overset/pimpleDyMOversetFoam/createFields.H

@@ -0,0 +1,68 @@
+    Info<< "Reading field p\n" << endl;
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+
+    Info<< "Reading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+#   include "createPhi.H"
+
+    // Set reference cell taking into account whether implicit or explicit
+    // overset aglorithm is used for field p.
+    const oversetMesh& om = oversetMesh::New(mesh);
+    label pRefCell = 0;
+    scalar pRefValue = 0.0;
+    om.setRefCell(p, pimple.dict(), pRefCell, pRefValue);
+
+    mesh.schemesDict().setFluxRequired(p.name());
+    if (correctPhi)
+    {
+        mesh.schemesDict().setFluxRequired("pcorr");
+    }
+
+    singlePhaseTransportModel laminarTransport(U, phi);
+
+    autoPtr<incompressible::turbulenceModel> turbulence
+    (
+        incompressible::turbulenceModel::New(U, phi, laminarTransport)
+    );
+
+    Info<< "Reading field rAU if present\n" << endl;
+    volScalarField rAU
+    (
+        IOobject
+        (
+            "rAU",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        runTime.deltaT(),
+        zeroGradientFvPatchScalarField::typeName
+    );
+
+#   include "writeOversetMasks.H"

applications/solvers/overset/pimpleDyMOversetFoam/pEqn.H

@@ -0,0 +1,76 @@
+{
+    p.boundaryField().updateCoeffs();
+
+    rAU = 1.0/UEqn.A();
+    oversetFvPatchScalarField::oversetInterpolate(rAU); // Overset update
+    surfaceScalarField rAUf = fvc::interpolate(rAU);
+
+    U = rAU*UEqn.H();
+    oversetFvPatchVectorField::oversetInterpolate(U); // Overset update
+
+    surfaceScalarField phiU
+    (
+        "phiU",
+        fvc::interpolate(U) & mesh.Sf()
+    );
+
+    phi = phiU;
+
+    // Adjust overset fluxes
+    oversetAdjustPhi(phi, U); // Fringe flux adjustment
+    globalOversetAdjustPhi(phi, U, p); // Global flux adjustment
+
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pEqn
+        (
+            fvm::laplacian(rAUf, p) == fvc::div(phi)
+        );
+
+        // Adjust non-orthogonal fringe fluxes if necessary
+        om.correctNonOrthoFluxes(pEqn, U);
+
+        pEqn.setReference(pRefCell, pRefValue);
+
+        if
+        (
+            pimple.corrPISO() == pimple.nCorrPISO()
+         && pimple.finalNonOrthogonalIter()
+        )
+        {
+            pEqn.solve(mesh.solutionDict().solver("pFinal"));
+        }
+        else
+        {
+            pEqn.solve();
+        }
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi -= pEqn.flux();
+        }
+
+        // Perform overset interpolation (after flux reconstruction)
+        oversetFvPatchScalarField::oversetInterpolate(p);
+    }
+
+#   include "oversetContinuityErrs.H"
+
+    // Explicitly relax pressure for momentum corrector
+    if (!pimple.finalIter())
+    {
+        p.relax();
+    }
+
+#   include "movingMeshContinuityErrs.H"
+
+    U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
+    U.correctBoundaryConditions();
+    oversetFvPatchVectorField::oversetInterpolate(U); // Overset update
+    // Note: if implicit fringe is switched on, we are doing the interpolation
+    // twice (once in correctBoundaryConditions and once in oversetInterpolate)
+    // Reorganize. VV, 4/Oct/2016.
+
+    // Make the fluxes relative to the mesh motion
+    fvc::makeRelative(phi, U);
+}

applications/solvers/overset/pimpleDyMOversetFoam/pimpleDyMOversetFoam.C

@@ -0,0 +1,129 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    pimpleDyMOversetFoam.C
+
+Description
+    Transient solver for incompressible, flow of Newtonian fluids
+    with with dynamic mesh and overset mesh using the PIMPLE
+    (merged PISO-SIMPLE) algorithm.
+
+    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "singlePhaseTransportModel.H"
+#include "turbulenceModel.H"
+#include "dynamicFvMesh.H"
+#include "pimpleControl.H"
+
+#include "oversetMesh.H"
+#include "oversetFvPatchFields.H"
+#include "oversetAdjustPhi.H"
+#include "globalOversetAdjustPhi.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+#   include "setRootCase.H"
+
+#   include "createTime.H"
+#   include "createDynamicFvMesh.H"
+
+    pimpleControl pimple(mesh);
+
+#   include "initContinuityErrs.H"
+#   include "createControls.H"
+#   include "createFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+#       include "readControls.H"
+
+        // Make the fluxes absolute
+        fvc::makeAbsolute(phi, U);
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        bool meshChanged = mesh.update();
+        reduce(meshChanged, orOp<bool>());
+
+#       include "createOversetMasks.H"
+#       include "volContinuity.H"
+
+        if (checkMeshCourantNo)
+        {
+#           include "meshCourantNo.H"
+        }
+
+        // Mesh motion update
+        if (correctPhi && meshChanged)
+        {
+            // Fluxes will be corrected to absolute velocity
+            // HJ, 6/Feb/2009
+#           include "correctPhi.H"
+        }
+
+        // Make the fluxes relative to the mesh motion
+        fvc::makeRelative(phi, U);
+
+#       include "oversetCourantNo.H"
+#       include "setDeltaT.H"
+
+        // --- PIMPLE loop
+        while (pimple.loop())
+        {
+#           include "UEqn.H"
+
+            // --- PISO loop
+            while (pimple.correct())
+            {
+#               include "pEqn.H"
+            }
+
+            turbulence->correct();
+        }
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/overset/pimpleDyMOversetFoam/readControls.H

@@ -0,0 +1,5 @@
+#include "readTimeControls.H"
+
+correctPhi = pimple.dict().lookupOrDefault("correctPhi", false);
+
+checkMeshCourantNo = pimple.dict().lookupOrDefault("checkMeshCourantNo", false);

applications/solvers/overset/potentialDyMOversetFoam/Make/files

@@ -0,0 +1,3 @@
+potentialDyMOversetFoam.C
+
+EXE = $(FOAM_APPBIN)/potentialDyMOversetFoam

applications/solvers/overset/potentialDyMOversetFoam/Make/options

@@ -0,0 +1,18 @@
+EXE_INC = \
+    -I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/overset/oversetMesh/lnInclude \
+    -I$(LIB_SRC)/overset/oversetDynamicFvMesh/lnInclude
+
+EXE_LIBS = \
+    -ldynamicFvMesh \
+    -ldynamicMesh \
+    -lengine \
+    -lmeshTools \
+    -lfiniteVolume \
+    -loversetMesh \
+    -loversetDynamicFvMesh \
+    -llduSolvers \
+    -L$(MESQUITE_LIB_DIR) -lmesquite

applications/solvers/overset/potentialDyMOversetFoam/createFields.H

@@ -0,0 +1,66 @@
+    Info<< "Reading field p\n" << endl;
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "Reading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    surfaceScalarField phi
+    (
+        IOobject
+        (
+            "phi",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        fvc::interpolate(U) & mesh.Sf()
+    );
+
+    // Write overset div phi for visualisation purposes
+    volScalarField oversetDivPhi
+    (
+        IOobject
+        (
+            "oversetDivPhi",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        fvc::div(phi)
+    );
+
+    // Set reference cell taking into account whether implicit or explicit
+    // overset algorithm is used for field p.
+    const oversetMesh& om = oversetMesh::New(mesh);
+    label pRefCell = 0;
+    scalar pRefValue = 0.0;
+    om.setRefCell(p, piso.dict(), pRefCell, pRefValue);
+
+    mesh.schemesDict().setFluxRequired(p.name());
+
+    // Write overset masks
+#   include "writeOversetMasks.H"

applications/solvers/overset/potentialDyMOversetFoam/potentialDyMOversetFoam.C

@@ -0,0 +1,190 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    potentialDyMOversetFoam
+
+Description
+    Transient solver for potential flow with dynamic overset mesh.
+
+Author
+    Hrvoje Jasak, Wikki Ltd.  All rights reserved.
+    Vuko Vukcevic, FMENA Zagreb. All rights reserved.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "dynamicFvMesh.H"
+#include "pisoControl.H"
+
+#include "oversetMesh.H"
+#include "oversetFvPatchFields.H"
+#include "oversetAdjustPhi.H"
+#include "globalOversetAdjustPhi.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::validOptions.insert("reconstructU", "");
+
+#   include "setRootCase.H"
+#   include "createTime.H"
+#   include "createDynamicFvMesh.H"
+
+    pisoControl piso(mesh);
+
+#   include "createFields.H"
+#   include "initTotalVolume.H"
+#   include "initContinuityErrs.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.loop())
+    {
+#       include "checkTotalVolume.H"
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        bool meshChanged = mesh.update();
+        reduce(meshChanged, orOp<bool>());
+
+#       include "createOversetMasks.H"
+
+        // Update moving wall velocity boundary condition and calculate the flux
+        U.correctBoundaryConditions();
+
+        // Forced overset update: make sure the overset interpolation is
+        // performed regardless of whether the coupledFringe is specified.
+        oversetFvPatchVectorField::oversetInterpolate(U);
+
+        phi == (linearInterpolate(U) & mesh.Sf());
+
+        // Resetting pressure field
+        p.internalField() = 0;
+
+#       include "volContinuity.H"
+#       include "meshCourantNo.H"
+
+        // Solve potential flow equations
+
+        // Adjust fluxes
+        oversetAdjustPhi(phi, U); // Fringe flux adjustment
+        globalOversetAdjustPhi(phi, U, p); // Global flux adjustment
+
+        while (piso.correctNonOrthogonal())
+        {
+            p.storePrevIter();
+
+            Info<< "Initial flux contour continuity error = "
+                << mag(sum(phi.boundaryField()))
+                << endl;
+
+            fvScalarMatrix pEqn
+            (
+                fvm::laplacian
+                (
+                    dimensionedScalar
+                    (
+                        "1",
+                        dimTime/p.dimensions()*dimensionSet(0, 2, -2, 0, 0),
+                        1
+                    ),
+                    p
+                )
+             ==
+                fvc::div(phi)
+            );
+
+            // Adjust non-orthogonal fringe fluxes if necessary
+            om.correctNonOrthoFluxes(pEqn, U);
+
+            pEqn.setReference(pRefCell, pRefValue);
+            pEqn.solve();
+
+            if (piso.finalNonOrthogonalIter())
+            {
+                phi -= pEqn.flux();
+#               include "oversetContinuityErrs.H"
+            }
+            else
+            {
+                p.relax();
+            }
+
+            // Perform overset interpolation (after flux reconstruction)
+            oversetFvPatchScalarField::oversetInterpolate(p);
+        }
+
+        // Update div phi field for visualisation purposes
+        oversetDivPhi = cellOversetMask*fvc::div(phi);
+
+        if (args.optionFound("reconstructU"))
+        {
+            U = fvc::reconstruct(phi);
+            U.correctBoundaryConditions();
+        }
+
+        Info<< "Interpolated U error = "
+            << (
+                   sqrt
+                   (
+                       sum
+                       (
+                           sqr
+                           (
+                               faceOversetMask*
+                               (
+                                   (fvc::interpolate(U) & mesh.Sf())
+                                 - phi
+                               )
+                           )
+                       )
+                   )/sum(mesh.magSf())
+               ).value()
+            << endl;
+
+        // Calculate velocity magnitude
+        {
+            volScalarField magU = mag(U);
+
+            Info<< "mag(U): max: " << gMax(magU.internalField())
+                << " min: " << gMin(magU.internalField()) << endl;
+        }
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return(0);
+}
+
+
+// ************************************************************************* //

applications/solvers/overset/potentialOversetFoam/Make/files

@@ -0,0 +1,3 @@
+potentialOversetFoam.C
+
+EXE = $(FOAM_APPBIN)/potentialOversetFoam

applications/solvers/overset/potentialOversetFoam/Make/options

@@ -0,0 +1,10 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/triSurface/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/overset/oversetMesh/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -loversetMesh \
+    -llduSolvers

applications/solvers/overset/potentialOversetFoam/createFields.H

@@ -0,0 +1,56 @@
+    Info<< "Reading field p\n" << endl;
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh
+    );
+
+    // Do not reset p
+
+    Info<< "Reading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    if (args.optionFound("resetU"))
+    {
+        U.internalField() = vector::zero;
+    }
+
+    surfaceScalarField phi
+    (
+        IOobject
+        (
+            "phi",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        fvc::interpolate(U) & mesh.Sf()
+    );
+
+
+    // Set reference cell taking into account whether implicit or explicit
+    // overset algorithm is used for field p.
+    label pRefCell = 0;
+    scalar pRefValue = 0.0;
+    om.setRefCell(p, simple.dict(), pRefCell, pRefValue);
+
+    mesh.schemesDict().setFluxRequired(p.name());

applications/solvers/overset/potentialOversetFoam/potentialOversetFoam.C

@@ -0,0 +1,262 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    potentialOversetFoam
+
+Description
+    Potential flow solver with overset mesh support.
+
+Author
+    Hrvoje Jasak, Wikki Ltd.  All rights reserved
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "simpleControl.H"
+
+#include "oversetMesh.H"
+#include "oversetFvPatchFields.H"
+#include "oversetAdjustPhi.H"
+#include "globalOversetAdjustPhi.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::validOptions.insert("resetU", "");
+    argList::validOptions.insert("writep", "");
+
+#   include "setRootCase.H"
+
+#   include "createTime.H"
+#   include "createMesh.H"
+
+    simpleControl simple(mesh);
+
+#   include "createOversetMasks.H"
+#   include "createFields.H"
+
+#   include "writeOversetMasks.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< nl << "Calculating potential flow" << endl;
+
+    // Do correctors over the complete set
+    while (simple.correctNonOrthogonal())
+    {
+        phi = (linearInterpolate(U) & mesh.Sf());
+
+        // Adjust overset fluxes
+        oversetAdjustPhi(phi, U); // Fringe flux adjustment
+        globalOversetAdjustPhi(phi, U, p); // Global flux adjustment
+
+        Info<< "Initial flux contour continuity error = "
+            << mag(sum(phi.boundaryField()))
+            << endl;
+
+        p.storePrevIter();
+
+        fvScalarMatrix pEqn
+        (
+            fvm::laplacian
+            (
+                dimensionedScalar
+                (
+                    "1",
+                    dimTime/p.dimensions()*dimensionSet(0, 2, -2, 0, 0),
+                    1
+                ),
+                p
+            )
+         ==
+            fvc::div(phi)
+        );
+
+        // Adjust non-orthogonal fringe fluxes if necessary
+        om.correctNonOrthoFluxes(pEqn, U);
+
+        pEqn.setReference(pRefCell, pRefValue);
+        pEqn.solve();
+
+        // Correct the flux
+        phi -= pEqn.flux();
+
+        volScalarField divFlux
+        (
+            "divFlux",
+            cellOversetMask*fvc::div(phi)
+        );
+        divFlux.write();
+
+        // Perform overset interpolation (after flux reconstruction)
+        oversetFvPatchScalarField::oversetInterpolate(p);
+
+        if (!simple.finalNonOrthogonalIter())
+        {
+            p.relax();
+        }
+
+        Info<< "p min " << gMin(p.internalField())
+            << " max " << gMax(p.internalField())
+            << " masked min "
+            << gMin(cellOversetMask.internalField()*p.internalField())
+            << " max "
+            << gMax(cellOversetMask.internalField()*p.internalField())
+            << endl;
+
+        Info<< "continuity error = "
+            << mag
+               (
+                    cellOversetMask*fvc::div(phi)
+               )().weightedAverage(mesh.V()).value()
+            << endl;
+
+        Info<< "Contour continuity error = "
+            << mag(sum(phi.boundaryField()))
+            << endl;
+
+        U = fvc::reconstruct(phi);
+        U.correctBoundaryConditions();
+        oversetFvPatchVectorField::oversetInterpolate(U); // Overset update
+        // Note: if implicit fringe is switched on, we are doing the
+        // interpolation twice (once in correctBoundaryConditions and once
+        // in oversetInterpolate). Reorganize. VV, 4/Oct/2016.
+
+        Info<< "Interpolated U error = "
+            << (
+                   sqrt
+                   (
+                       sum
+                       (
+                           sqr
+                           (
+                               faceOversetMask*
+                               (
+                                   (fvc::interpolate(U) & mesh.Sf())
+                                 - phi
+                               )
+                           )
+                       )
+                   )/sum(mesh.magSf())
+               ).value()
+            << endl;
+    }
+
+    // Calculate velocity magnitude
+    {
+        volScalarField magU = cellOversetMask*mag(U);
+
+        Info << "Masked mag(U): max: " << gMax(magU.internalField())
+            << " min: " << gMin(magU.internalField()) << endl;
+    }
+
+    // Force the write
+    U.write();
+    p.write();
+    phi.write();
+    cellOversetMask.write();
+
+    if (args.optionFound("writep"))
+    {
+        // Find reference patch
+        label refPatch = -1;
+        scalar maxMagU = 0;
+
+        // Go through all velocity patches and find the one that fixes
+        // velocity to the largest value
+
+        forAll (U.boundaryField(), patchI)
+        {
+            const fvPatchVectorField& Upatch = U.boundaryField()[patchI];
+
+            if (Upatch.fixesValue())
+            {
+                // Calculate mean velocity
+                scalar u = sum(mag(Upatch));
+                label patchSize = Upatch.size();
+
+                reduce(u, sumOp<scalar>());
+                reduce(patchSize, sumOp<label>());
+
+                if (patchSize > 0)
+                {
+                    scalar curMag = u/patchSize;
+
+                    if (curMag > maxMagU)
+                    {
+                        refPatch = patchI;
+
+                        maxMagU = curMag;
+                    }
+                }
+            }
+        }
+
+        if (refPatch > -1)
+        {
+            // Calculate reference pressure
+            const fvPatchVectorField& Upatch = U.boundaryField()[refPatch];
+            const fvPatchScalarField& pPatch = p.boundaryField()[refPatch];
+
+            scalar patchE = sum(mag(pPatch + 0.5*magSqr(Upatch)));
+            label patchSize = Upatch.size();
+
+            reduce(patchE, sumOp<scalar>());
+            reduce(patchSize, sumOp<label>());
+
+            scalar e = patchE/patchSize;
+
+
+            Info<< "Using reference patch " << refPatch
+                << " with mag(U) = " << maxMagU
+                << " p + 0.5*U^2 = " << e << endl;
+
+            p.internalField() = e - 0.5*magSqr(U.internalField());
+            p.correctBoundaryConditions();
+            oversetFvPatchScalarField::oversetInterpolate(p); // Overset update
+            // Note: if implicit fringe is switched on, we are doing the
+            // interpolation twice (once in correctBoundaryConditions and once
+            // in oversetInterpolate). Reorganize. VV, 4/Oct/2016.
+        }
+        else
+        {
+            Info<< "No reference patch found.  Writing potential function"
+                << endl;
+        }
+
+        p.write();
+    }
+
+    Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+        << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+        << nl << endl;
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/overset/scalarTransportOversetFoam/Make/files

@@ -0,0 +1,3 @@
+scalarTransportOversetFoam.C
+
+EXE = $(FOAM_APPBIN)/scalarTransportOversetFoam

applications/solvers/overset/scalarTransportOversetFoam/Make/options

@@ -0,0 +1,8 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/overset/oversetMesh/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -loversetMesh \
+    -llduSolvers

applications/solvers/overset/scalarTransportOversetFoam/createFields.H

@@ -0,0 +1,75 @@
+    Info<< "Reading field T\n" << endl;
+
+    volScalarField T
+    (
+        IOobject
+        (
+            "T",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "Reading field V\n" << endl;
+
+    volVectorField V
+    (
+        IOobject
+        (
+            "V",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+
+    Info<< "Reading field U\n" << endl;
+
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+
+    Info<< "Reading transportProperties\n" << endl;
+
+    IOdictionary transportProperties
+    (
+        IOobject
+        (
+            "transportProperties",
+            runTime.constant(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    );
+
+
+    Info<< "Reading diffusivity DT\n" << endl;
+
+    dimensionedScalar DT
+    (
+        transportProperties.lookup("DT")
+    );
+
+    dimensionedScalar DV
+    (
+        transportProperties.lookup("DV")
+    );
+
+#   include "createPhi.H"

applications/solvers/overset/scalarTransportOversetFoam/scalarTransportOversetFoam.C

@@ -0,0 +1,120 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    scalarTransportOversetFoam
+
+Description
+    Solves a transport equation for a passive scalar with support for
+    overset meshes.
+
+    Experimental: no overset mesh changes required
+
+Author
+    Hrvoje Jasak, Wikki Ltd.  All rights reserved
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "simpleControl.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+#   include "setRootCase.H"
+
+#   include "createTime.H"
+#   include "createMesh.H"
+
+    simpleControl simple(mesh);
+
+#   include "createFields.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nCalculating scalar transport\n" << endl;
+
+#   include "CourantNo.H"
+
+    while (simple.loop())
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        fvScalarMatrix TEqn
+        (
+            fvm::ddt(T)
+          + fvm::div(phi, T)
+          - fvm::laplacian(DT, T)
+        );
+
+        TEqn.solve();
+
+//         volScalarField Tresidual
+//         (
+//             "Tresidual",
+//             T
+//         );
+//         Tresidual.writeOpt() = IOobject::AUTO_WRITE;
+
+//         Tresidual.internalField() = TEqn.residual();
+//         Tresidual.boundaryField() == 0;
+
+//         Info<< "residual " << gSumMag(Tresidual.internalField()) << endl;
+
+        fvVectorMatrix VEqn
+        (
+            fvm::ddt(V)
+          + fvm::div(phi, V)
+          - fvm::laplacian(DV, V)
+        );
+
+        VEqn.solve();
+
+//         return 0;
+
+//         VEqn.solve();
+
+//         volVectorField Vresidual
+//         (
+//             "Vresidual",
+//             V
+//         );
+//         Vresidual.writeOpt() = IOobject::AUTO_WRITE;
+
+//         Vresidual.internalField() = VEqn.residual();
+//         Vresidual.boundaryField() == vector::zero;
+
+//         Info<< "residual " << gSumMag(Vresidual.internalField()) << endl;
+
+        runTime.write();
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/overset/simpleMRFOversetFoam/Make/files

@@ -0,0 +1,3 @@
+simpleMRFOversetFoam.C 
+
+EXE = $(FOAM_APPBIN)/simpleMRFOversetFoam

applications/solvers/overset/simpleMRFOversetFoam/Make/options

@@ -0,0 +1,16 @@
+EXE_INC = \
+    -I$(LIB_SRC)/turbulenceModels \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/triSurface/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/overset/oversetMesh/lnInclude
+
+EXE_LIBS = \
+    -lincompressibleRASModels \
+    -lincompressibleTransportModels \
+    -lfiniteVolume \
+    -loversetMesh \
+    -llduSolvers

applications/solvers/overset/simpleMRFOversetFoam/UEqn.H

@@ -0,0 +1,14 @@
+    // Solve the Momentum equation
+    tmp<fvVectorMatrix> UEqn
+    (
+        fvm::div(phi, U)
+      + turbulence->divDevReff()
+    );
+    mrfZones.addCoriolis(UEqn());
+
+    UEqn().relax();
+
+    solve
+    (
+        UEqn() == -fvc::grad(p)
+    );

applications/solvers/overset/simpleMRFOversetFoam/createFields.H

@@ -0,0 +1,83 @@
+    Info << "Reading field p\n" << endl;
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "Reading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+#   include "createPhi.H"
+
+    // Set reference cell taking into account whether implicit or explicit
+    // overset aglorithm is used for field p.
+    label pRefCell = 0;
+    scalar pRefValue = 0.0;
+    om.setRefCell(p, simple.dict(), pRefCell, pRefValue);
+
+    mesh.schemesDict().setFluxRequired(p.name());
+
+
+    singlePhaseTransportModel laminarTransport(U, phi);
+
+    autoPtr<incompressible::RASModel> turbulence
+    (
+        incompressible::RASModel::New(U, phi, laminarTransport)
+    );
+
+    // Create MRF zones
+    MRFZones mrfZones(mesh);
+    mrfZones.correctBoundaryVelocity(U);
+
+    // Create Urel as a permanent field to make it available for on-the-fly
+    // post-processing operations
+    volVectorField Urel
+    (
+        IOobject
+        (
+            "Urel",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        U
+    );
+    mrfZones.relativeVelocity(Urel);
+
+    Info<< "Reading field rAU if present\n" << endl;
+    volScalarField rAU
+    (
+        IOobject
+        (
+            "rAU",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("dt", dimTime, 1.0),
+        zeroGradientFvPatchScalarField::typeName
+    );
+
+#   include "writeOversetMasks.H"

applications/solvers/overset/simpleMRFOversetFoam/pEqn.H

@@ -0,0 +1,59 @@
+{
+    p.boundaryField().updateCoeffs();
+
+    rAU = 1.0/UEqn().A();
+    oversetFvPatchScalarField::oversetInterpolate(rAU); // Overset update
+
+    // Update boundary velocity for consistency with the flux
+    // This is needed for a ramped MRF
+    mrfZones.correctBoundaryVelocity(U);
+
+    U = rAU*UEqn().H();
+    oversetFvPatchVectorField::oversetInterpolate(U); // Overset update
+
+    UEqn.clear();
+    phi = (fvc::interpolate(U) & mesh.Sf());
+    mrfZones.relativeFlux(phi);
+
+    // Adjust overset fluxes
+    oversetAdjustPhi(phi, U); // Fringe flux adjustment
+    globalOversetAdjustPhi(phi, U, p); // Global flux adjustment
+
+    // Non-orthogonal pressure corrector loop
+    while (simple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pEqn
+        (
+            fvm::laplacian(rAU, p) == fvc::div(phi)
+        );
+
+        // Adjust non-orthogonal fringe fluxes if necessary
+        om.correctNonOrthoFluxes(pEqn, U);
+
+        pEqn.setReference(pRefCell, pRefValue);
+
+        // Retain the residual from the first iteration
+        pEqn.solve();
+
+        if (simple.finalNonOrthogonalIter())
+        {
+            phi -= pEqn.flux();
+        }
+
+        // Perform overset interpolation (after flux reconstruction)
+        oversetFvPatchScalarField::oversetInterpolate(p);
+    }
+
+#   include "oversetContinuityErrs.H"
+
+    // Explicitly relax pressure for momentum corrector
+    p.relax();
+
+    // Momentum corrector
+    U -= rAU*fvc::grad(p);
+    U.correctBoundaryConditions();
+    oversetFvPatchVectorField::oversetInterpolate(U); // Overset update
+    // Note: if implicit fringe is switched on, we are doing the interpolation
+    // twice (once in correctBoundaryConditions and once in oversetInterpolate)
+    // Reorganize. VV, 4/Oct/2016.
+}

applications/solvers/overset/simpleMRFOversetFoam/simpleMRFOversetFoam.C

@@ -0,0 +1,94 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    simpleMRFOversetFoam
+
+Description
+    Steady-state solver for incompressible, turbulent flow
+    with overset mesh support and MRF regions.
+
+Author
+    Vuko Vukcevic, Wikki Ltd.  All rights reserved.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "singlePhaseTransportModel.H"
+#include "RASModel.H"
+#include "simpleControl.H"
+#include "MRFZones.H"
+
+#include "oversetMesh.H"
+#include "oversetFvPatchFields.H"
+#include "oversetAdjustPhi.H"
+#include "globalOversetAdjustPhi.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+#   include "setRootCase.H"
+#   include "createTime.H"
+#   include "createMesh.H"
+
+    simpleControl simple(mesh);
+
+#   include "createOversetMasks.H"
+#   include "createFields.H"
+#   include "initContinuityErrs.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (simple.loop())
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // Pressure-velocity SIMPLE corrector
+        {
+#           include "UEqn.H"
+#           include "pEqn.H"
+        }
+
+        // Calculate relative velocity
+        Urel == U;
+        mrfZones.relativeVelocity(Urel);
+
+        turbulence->correct();
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/overset/simpleOversetFoam/Make/files

@@ -0,0 +1,3 @@
+simpleOversetFoam.C 
+
+EXE = $(FOAM_APPBIN)/simpleOversetFoam

applications/solvers/overset/simpleOversetFoam/Make/options

@@ -0,0 +1,16 @@
+EXE_INC = \
+    -I$(LIB_SRC)/turbulenceModels \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/triSurface/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/overset/oversetMesh/lnInclude
+
+EXE_LIBS = \
+    -lincompressibleRASModels \
+    -lincompressibleTransportModels \
+    -lfiniteVolume \
+    -loversetMesh \
+    -llduSolvers

applications/solvers/overset/simpleOversetFoam/UEqn.H

@@ -0,0 +1,13 @@
+    // Solve the Momentum equation
+    tmp<fvVectorMatrix> UEqn
+    (
+        fvm::div(phi, U)
+      + turbulence->divDevReff()
+    );
+
+    UEqn().relax();
+
+    solve
+    (
+        UEqn() == -fvc::grad(p)
+    );

applications/solvers/overset/simpleOversetFoam/createFields.H

@@ -0,0 +1,63 @@
+    Info << "Reading field p\n" << endl;
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "Reading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+#   include "createPhi.H"
+
+    // Set reference cell taking into account whether implicit or explicit
+    // overset aglorithm is used for field p.
+    label pRefCell = 0;
+    scalar pRefValue = 0.0;
+    om.setRefCell(p, simple.dict(), pRefCell, pRefValue);
+
+    mesh.schemesDict().setFluxRequired(p.name());
+
+
+    singlePhaseTransportModel laminarTransport(U, phi);
+
+    autoPtr<incompressible::RASModel> turbulence
+    (
+        incompressible::RASModel::New(U, phi, laminarTransport)
+    );
+
+    Info<< "Reading field rAU if present\n" << endl;
+    volScalarField rAU
+    (
+        IOobject
+        (
+            "rAU",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("dt", dimTime, 1.0),
+        zeroGradientFvPatchScalarField::typeName
+    );
+
+#   include "writeOversetMasks.H"

applications/solvers/overset/simpleOversetFoam/pEqn.H

@@ -0,0 +1,54 @@
+{
+    p.boundaryField().updateCoeffs();
+
+    rAU = 1.0/UEqn().A();
+    oversetFvPatchScalarField::oversetInterpolate(rAU); // Overset update
+
+    U = rAU*UEqn().H();
+    oversetFvPatchVectorField::oversetInterpolate(U); // Overset update
+
+    UEqn.clear();
+    phi = (fvc::interpolate(U) & mesh.Sf());
+
+    // Adjust overset fluxes
+    oversetAdjustPhi(phi, U); // Fringe flux adjustment
+    globalOversetAdjustPhi(phi, U, p); // Global flux adjustment
+
+    // Non-orthogonal pressure corrector loop
+    while (simple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pEqn
+        (
+            fvm::laplacian(rAU, p) == fvc::div(phi)
+        );
+
+        // Adjust non-orthogonal fringe fluxes if necessary
+        om.correctNonOrthoFluxes(pEqn, U);
+
+        pEqn.setReference(pRefCell, pRefValue);
+
+        // Retain the residual from the first iteration
+        pEqn.solve();
+
+        if (simple.finalNonOrthogonalIter())
+        {
+            phi -= pEqn.flux();
+        }
+
+        // Perform overset interpolation (after flux reconstruction)
+        oversetFvPatchScalarField::oversetInterpolate(p);
+    }
+
+#   include "oversetContinuityErrs.H"
+
+    // Explicitly relax pressure for momentum corrector
+    p.relax();
+
+    // Momentum corrector
+    U -= rAU*fvc::grad(p);
+    U.correctBoundaryConditions();
+    oversetFvPatchVectorField::oversetInterpolate(U); // Overset update
+    // Note: if implicit fringe is switched on, we are doing the interpolation
+    // twice (once in correctBoundaryConditions and once in oversetInterpolate)
+    // Reorganize. VV, 4/Oct/2016.
+}

applications/solvers/overset/simpleOversetFoam/simpleOversetFoam.C

@@ -0,0 +1,89 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    simpleOversetFoam
+
+Description
+    Steady-state solver for incompressible, turbulent flow
+    with overset mesh support.
+
+Author
+    Hrvoje Jasak, Wikki Ltd.  All rights reserved
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "singlePhaseTransportModel.H"
+#include "RASModel.H"
+#include "simpleControl.H"
+
+#include "oversetMesh.H"
+#include "oversetFvPatchFields.H"
+#include "oversetAdjustPhi.H"
+#include "globalOversetAdjustPhi.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+#   include "setRootCase.H"
+#   include "createTime.H"
+#   include "createMesh.H"
+
+    simpleControl simple(mesh);
+
+#   include "createOversetMasks.H"
+#   include "createFields.H"
+#   include "initContinuityErrs.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (simple.loop())
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // Pressure-velocity SIMPLE corrector
+        {
+#           include "UEqn.H"
+#           include "pEqn.H"
+        }
+
+        turbulence->correct();
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/utilities/mesh/manipulation/deformedGeom/deformedGeom.C

@@ -53,8 +53,6 @@ int main(int argc, char *argv[])
 
     instantList timeDirs = timeSelector::select0(runTime, args);
 
-    volPointInterpolation pInterp(mesh);
-
     pointField zeroPoints(mesh.points());
 
     forAll(timeDirs, timeI)
@@ -73,6 +71,9 @@ int main(int argc, char *argv[])
             IOobject::MUST_READ
         );
 
+        // Create volPointInterpolation object after the mesh has been read
+        volPointInterpolation pInterp(mesh);
+
         // Check U exists
         if (Uheader.headerOk())
         {

applications/utilities/overset/calcOverset/Make/files

@@ -0,0 +1,3 @@
+calcOverset.C 
+
+EXE = $(FOAM_APPBIN)/calcOverset

applications/utilities/overset/calcOverset/Make/options

@@ -0,0 +1,16 @@
+EXE_INC = \
+    -I$(LIB_SRC)/turbulenceModels \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/triSurface/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/overset/oversetMesh/lnInclude
+
+EXE_LIBS = \
+    -lincompressibleRASModels \
+    -lincompressibleTransportModels \
+    -lfiniteVolume \
+    -loversetMesh \
+    -llduSolvers

applications/utilities/overset/calcOverset/calcOverset.C

@@ -0,0 +1,87 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+    calcOverset
+
+Description
+    Calculates overset addressing and everything necessary for overset
+    interpolation. Used for parallel scaling tests of overset assembly.
+
+Author
+    Vuko Vukcevic, FMENA Zagreb.  All rights reserved.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "oversetMesh.H"
+#include "oversetFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+#   include "setRootCase.H"
+#   include "createTime.H"
+#   include "createMesh.H"
+
+    volVectorField cellCentres
+    (
+        IOobject
+        (
+            "cellCentres",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh.C(),
+        "zeroGradient"
+    );
+
+    volVectorField origCellCentres("origCellCentres", cellCentres);
+    origCellCentres.write();
+
+    // Create CPU time object
+    cpuTime oversetTime;
+
+    // Make sure everything necessary for overset interpolation has been
+    // calculated by calling a single overset interpolation
+    oversetFvPatchVectorField::oversetInterpolate(cellCentres);
+
+    Info<< "Elapsed time for overset assembly and single interpolation: "
+        << oversetTime.cpuTimeIncrement() << endl;
+
+    // Write interpolated cell centres
+    cellCentres.write();
+
+#   include "writeOversetMasks.H"
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/utilities/overset/oversetRegionIndex/Make/files

@@ -0,0 +1,3 @@
+oversetRegionIndex.C
+
+EXE = $(FOAM_APPBIN)/oversetRegionIndex

applications/utilities/overset/oversetRegionIndex/Make/options

@@ -0,0 +1,7 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools

applications/utilities/overset/oversetRegionIndex/oversetRegionIndex.C

@@ -0,0 +1,82 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+    oversetRegionID
+
+Description
+    Visualise region split in overset mesh
+
+Author
+    Hrvoje Jasak, Wikki Ltd.  All rights reserved.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "regionSplit.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+#   include "setRootCase.H"
+#   include "createTime.H"
+#   include "createMesh.H"
+
+    volScalarField regionIndex
+    (
+        IOobject
+        (
+            "regionIndex",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimless, 0)
+    );
+    scalarField& rIn = regionIndex.internalField();
+
+    // Create region split.  Region index depends only on mesh ordering
+    regionSplit rs(mesh);
+
+    forAll (rIn, cellI)
+    {
+        rIn[cellI] = rs[cellI];
+    }
+
+    Info<< "Number of regions: " << rs.nRegions() << endl;
+
+    Info << "Write overset region split ... ";
+
+    regionIndex.write();
+
+    Info << "done" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

src/Allwmake

@@ -81,6 +81,9 @@ wmake libso immersedBoundary/immersedBoundaryTurbulence
 wmake libso immersedBoundary/immersedBoundaryForce
 wmake libso immersedBoundary/immersedBoundaryDynamicMesh
 
+wmake libso overset/oversetMesh
+wmake libso overset/oversetDynamicFvMesh
+
 ( cd cudaSolvers ; ./Allwmake )
 
 # ----------------------------------------------------------------- end-of-file

src/dynamicMesh/dynamicFvMesh/Make/files

@@ -10,7 +10,6 @@ dynamicMotionSolverFvMesh/dynamicMotionSolverFvMesh.C
 subsetMotionSolverFvMesh/subsetMotionSolverFvMesh.C
 dynamicInkJetFvMesh/dynamicInkJetFvMesh.C
 dynamicRefineFvMesh/dynamicRefineFvMesh.C
-dynamicRefinePolyFvMesh/dynamicRefinePolyFvMesh.C
 
 mixerGgiFvMesh/mixerGgiFvMesh.C
 turboFvMesh/turboFvMesh.C

src/dynamicMesh/dynamicFvMesh/dynamicRefinePolyFvMesh/dynamicMeshDict

@@ -1,77 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | foam-extend: Open Source CFD                    |
-|  \\    /   O peration     | Version:     4.0                                |
-|   \\  /    A nd           | Web:         http://www.foam-extend.org         |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    object      dynamicMeshDict;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-//dynamicFvMeshLib   "libtopoChangerFvMesh.so";
-dynamicFvMesh      dynamicRefinePolyFvMesh;
-//staticFvMesh;
-
-mixerFvMeshCoeffs
-{
-    coordinateSystem
-    {
-        type            cylindrical;
-        origin          (0 0 0);
-        axis            (0 0 1);
-        direction       (1 0 0);
-    }
-
-    rpm             10;
-
-    slider
-    {
-        inside      insideSlider;
-        outside     outsideSlider;
-    }
-}
-
-// Refinement
-dynamicRefineFvMeshCoeffs
-{
-    // Refine every refineInterval timesteps
-    refineInterval 3;
-
-    // Maximum refinement level (starts from 0)
-    maxRefinement 2;
-
-    // Maximum cell limit (approximate)
-    maxCells 1000000;
-
-    // volScalarField to base refinement on
-    field gamma;
-
-    // Which cells to un/refine: based on point values (simple averaging).
-    // - refine pointCells of point value inbetween minLevel..maxLevel
-    // - unrefine pointCells that are within nBufferLayers of points marked
-    //   for refinement.
-    minLevel 0.01;
-    maxLevel 0.99;
-    nBufferLayers 1;
-
-    // Newly introduced patch points optionally get projected onto a surface
-    //projectSurfaces ("fixedWalls4.stl");
-    //projectPatches (fixedWalls);
-    // Maximum project distance
-    //projectDistance 1;
-
-    // Fluxes to adapt. For newly created faces or split faces the flux
-    // gets estimated from an interpolated volVectorField ('velocity')
-    // First is name of the flux to adapt, second is velocity that will
-    // be interpolated and inner-producted with the face area vector.
-    correctFluxes ((phi U));
-}
-
-
-// ************************************************************************* //

src/dynamicMesh/dynamicFvMesh/dynamicRefinePolyFvMesh/dynamicRefinePolyFvMesh.H

@@ -1,215 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | foam-extend: Open Source CFD
-   \\    /   O peration     | Version:     4.0
-    \\  /    A nd           | Web:         http://www.foam-extend.org
-     \\/     M anipulation  | For copyright notice see file Copyright
--------------------------------------------------------------------------------
-License
-    This file is part of foam-extend.
-
-    foam-extend 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.
-
-    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::dynamicRefinePolyFvMesh
-
-Description
-    A fvMesh with built-in refinement of arbitrary polyhedral cells.
-
-    Determines which cells to refine/unrefine and does all in update().
-
-SourceFiles
-    dynamicRefinePolyFvMesh.C
-
-Author
-    Vuko Vukcevic, Wikki Ltd.  All rights reserved
-
-Notes
-    Generalisation of dynamicRefineMesh for polyhedral cells
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef dynamicRefinePolyFvMesh_H
-#define dynamicRefinePolyFvMesh_H
-
-#include "dynamicFvMesh.H"
-#include "polyRef.H"
-#include "PackedBoolList.H"
-#include "Switch.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                      Class dynamicRefinePolyFvMesh Declaration
-\*---------------------------------------------------------------------------*/
-
-class dynamicRefinePolyFvMesh
-:
-    public dynamicFvMesh
-{
-
-    // Private data
-
-        // Helper variables to enable switching between a single and multiple
-        // mesh motion updates within a time step (if update() is called more
-        // than once in a single time step)
-
-            //- Switch for single motion update (true by default)
-            Switch singleMotionUpdate_;
-
-            //- Helper varaible: current time index
-            label curTimeIndex_;
-
-
-protected:
-
-        //- Mesh cutting engine
-        polyRef meshCutter_;
-
-        //- Dump cellLevel for postprocessing
-        Switch dumpLevel_;
-
-        //- Fluxes to map
-        List<Pair<word> > correctFluxes_;
-
-        //- Number of refinement/unrefinement steps done so far.
-        label nRefinementIterations_;
-
-
-    // Private Member Functions
-
-        //- Count set/unset elements in packedlist.
-        static label count(const PackedBoolList&, const unsigned int);
-
-        //- Read the projection parameters from dictionary
-        void readDict();
-
-
-        //- Refine cells. Update mesh and fields.
-        autoPtr<mapPolyMesh> refine(const labelList&);
-
-        //- Unrefine cells. Gets passed in centre points of cells to combine.
-        autoPtr<mapPolyMesh> unrefine(const labelList&);
-
-
-        // Selection of cells to un/refine
-
-            //- Get per cell max of connected point
-            scalarField maxPointField(const scalarField&) const;
-
-            //- Get point min of connected cell
-            scalarField minCellField(const volScalarField&) const;
-
-            //- Simple (non-parallel) interpolation by averaging
-            scalarField cellToPoint(const scalarField& vFld) const;
-
-            //- Calculate error (= -1 by default or distance from inbetween
-            //  levels
-            scalarField error
-            (
-                const scalarField& fld,
-                const scalar minLevel,
-                const scalar maxLevel
-            ) const;
-
-            //- Select candidate cells for refinement
-            virtual void selectRefineCandidates
-            (
-                const scalar lowerRefineLevel,
-                const scalar upperRefineLevel,
-                const scalarField& vFld,
-                PackedBoolList& candidateCell
-            ) const;
-
-            //- Subset candidate cells for refinement
-            virtual labelList selectRefineCells
-            (
-                const label maxCells,
-                const label maxRefinement,
-                const PackedBoolList& candidateCell
-            ) const;
-
-            //- Select points that can be unrefined
-            virtual labelList selectUnrefinePoints
-            (
-                const scalar unrefineLevel,
-                const PackedBoolList& markedCell,
-                const scalarField& pFld
-            ) const;
-
-            //- Extend markedCell with cell-face-cell
-            void extendMarkedCells(PackedBoolList& markedCell) const;
-
-
-private:
-
-        //- Disallow default bitwise copy construct
-        dynamicRefinePolyFvMesh(const dynamicRefinePolyFvMesh&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const dynamicRefinePolyFvMesh&);
-
-public:
-
-    //- Runtime type information
-    TypeName("dynamicRefinePolyFvMesh");
-
-
-    // Constructors
-
-        //- Construct from IOobject
-        explicit dynamicRefinePolyFvMesh(const IOobject& io);
-
-
-    // Destructor
-
-        virtual ~dynamicRefinePolyFvMesh();
-
-
-    // Member Functions
-
-        //- Direct access to the refinement engine
-        const polyRef& meshCutter() const
-        {
-            return meshCutter_;
-        }
-
-        //- Update the mesh for both mesh motion and topology change
-        virtual bool update();
-
-
-    // Writing
-
-        //- Write using given format, version and compression
-        virtual bool writeObject
-        (
-            IOstream::streamFormat fmt,
-            IOstream::versionNumber ver,
-            IOstream::compressionType cmp
-        ) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/dynamicMesh/dynamicMesh/Make/files

@@ -37,6 +37,9 @@ $(slidingInterface)/decoupleSlidingInterface.C
 repatchCoverage = $(polyMeshModifiers)/repatchCoverage
 $(repatchCoverage)/repatchCoverage.C
 
+polyhedralRefinement = $(polyMeshModifiers)/polyhedralRefinement
+$(polyhedralRefinement)/polyhedralRefinement.C
+
 polyTopoChange/polyTopoChange/polyTopoChange.C
 polyTopoChange/polyTopoChange/actions/topoAction/topoActions.C
 
@@ -52,7 +55,6 @@ motionSolver/motionSolver.C
 refinementData/refinementData.C
 refinementData/refinementDistanceData.C
 refinementData/refinementHistory.C
-refinementData/polyRefinementHistory.C
 
 directTopoChange/directTopoChange/directTopoChange.C
 
@@ -61,7 +63,6 @@ $(directActions)/addPatchCellLayer.C
 $(directActions)/edgeCollapser.C
 $(directActions)/faceCollapser.C
 $(directActions)/hexRef8.C
-$(directActions)/polyRef.C
 $(directActions)/removeCells.C
 $(directActions)/removeFaces.C
 $(directActions)/removePoints.C

src/dynamicMesh/dynamicMesh/directTopoChange/directTopoChange/directActions/polyRef.H

@@ -1,558 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | foam-extend: Open Source CFD
-   \\    /   O peration     | Version:     4.0
-    \\  /    A nd           | Web:         http://www.foam-extend.org
-     \\/     M anipulation  | For copyright notice see file Copyright
--------------------------------------------------------------------------------
-License
-    This file is part of foam-extend.
-
-    foam-extend 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.
-
-    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::polyRef
-
-Description
-    Refinement of (split) polyhedra using directTopoChange.
-
-    Each polyhedron is split by the following procedure:
-    1. Adding points at the edge centre, face centre and cell centre,
-    2. Adding cells n cells where n is the number of points of the cell,
-    3. Splitting each face into multiple faces going from:
-       existing corner point -> new edge centre point -> new face centre
-       point -> other new edge centre point (sharing the same corner point)
-    4. Adding internal faces going from:
-       new edge centre point -> new face centre point -> new cell centre
-       point -> other new face centre point (sharing the same edge)
-
-SourceFiles
-    polyRef.C
-
-Author
-    Vuko Vukcevic, Wikki Ltd.  All rights reserved.
-
-Notes
-    Generalisation of hexRef8 for polyhedral cells
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef polyRef_H
-#define polyRef_H
-
-#include "labelIOList.H"
-#include "face.H"
-#include "HashSet.H"
-#include "DynamicList.H"
-#include "primitivePatch.H"
-#include "removeFaces.H"
-#include "polyRefinementHistory.H"
-#include "PackedList.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// Forward declaration of classes
-class polyMesh;
-class polyPatch;
-class directTopoChange;
-class mapPolyMesh;
-class mapDistributePolyMesh;
-
-/*---------------------------------------------------------------------------*\
-                           Class polyRef Declaration
-\*---------------------------------------------------------------------------*/
-
-class polyRef
-{
-    // Private data
-
-        //- Reference to underlying mesh
-        const polyMesh& mesh_;
-
-        //- Per cell the refinement level
-        labelIOList cellLevel_;
-
-        //- Per point the refinement level
-        labelIOList pointLevel_;
-
-        //- Typical edge length between unrefined points
-        const scalar level0Edge_;
-
-        //- Refinement history
-        polyRefinementHistory history_;
-
-        //- Face remover engine
-        removeFaces faceRemover_;
-
-        //- Level of saved points
-        Map<label> savedPointLevel_;
-
-        //- Level of saved cells
-        Map<label> savedCellLevel_;
-
-
-    // Private Member Functions
-
-        //- Reorder according to map
-        static void reorder
-        (
-            const labelList& map,
-            const label len,
-            const label null,
-            labelList& elems
-        );
-
-        //- Get patch and zone info
-        void getFaceInfo
-        (
-            const label faceI,
-            label& patchID,
-            label& zoneID,
-            label& zoneFlip
-        ) const;
-
-        //- Adds a face on top of existing faceI. Reverses if nessecary
-        label addFace
-        (
-            directTopoChange& meshMod,
-            const label faceI,
-            const face& newFace,
-            const label own,
-            const label nei
-        ) const;
-
-        //- Adds internal face from point. No checks on reversal
-        label addInternalFace
-        (
-            directTopoChange& meshMod,
-            const label meshFaceI,
-            const label meshPointI,
-            const face& newFace,
-            const label own,
-            const label nei
-        ) const;
-
-        //- Modifies existing faceI for either new owner/neighbour or new face
-        //  points. Reverses if nessecary
-        void modFace
-        (
-            directTopoChange& meshMod,
-            const label faceI,
-            const face& newFace,
-            const label own,
-            const label nei
-        ) const;
-
-        scalar getLevel0EdgeLength() const;
-
-        //- Get cell added to point of cellI (if any)
-        label getAnchorCell
-        (
-            const labelListList& cellAnchorPoints,
-            const labelListList& cellAddedCells,
-            const label cellI,
-            const label faceI,
-            const label pointI
-        ) const;
-
-        //- Get new owner and neighbour (in unspecified order) of pointI
-        //  on faceI
-        void getFaceNeighbours
-        (
-            const labelListList& cellAnchorPoints,
-            const labelListList& cellAddedCells,
-            const label faceI,
-            const label pointI,
-
-            label& own,
-            label& nei
-        ) const;
-
-
-        //- Get index of point with minimum pointlevel
-        label findMinLevel(const labelList& f) const;
-
-        //- Get index of point with maximum pointlevel
-        label findMaxLevel(const labelList& f) const;
-
-        //- Count number of vertices <= anchorLevel
-        label countAnchors(const labelList&, const label) const;
-
-        //- Find index of point with wantedLevel, starting from fp
-        label findLevel
-        (
-            const face& f,
-            const label startFp,
-            const bool searchForward,
-            const label wantedLevel
-        ) const;
-
-        ////- Print levels of list of points.
-        //void printLevels(Ostream&, const labelList&) const;
-
-        //- debug: check orientation of added internal face
-        static void checkInternalOrientation
-        (
-            directTopoChange& meshMod,
-            const label cellI,
-            const label faceI,
-            const point& ownPt,
-            const point& neiPt,
-            const face& newFace
-        );
-
-        //- debug: check orientation of new boundary face
-        static void checkBoundaryOrientation
-        (
-            directTopoChange& meshMod,
-            const label cellI,
-            const label faceI,
-            const point& ownPt,
-            const point& boundaryPt,
-            const face& newFace
-        );
-
-        //- If p0 and p1 are existing vertices check if edge is split and insert
-        //  splitPoint
-        void insertEdgeSplit
-        (
-            const labelList& edgeMidPoint,
-            const label p0,
-            const label p1,
-            dynamicLabelList& verts
-        ) const;
-
-        //- Store in maps correspondence from midpoint to anchors and faces
-        label storeMidPointInfo
-        (
-            const labelListList& cellAnchorPoints,
-            const labelListList& cellAddedCells,
-            const labelList& cellMidPoint,
-            const labelList& edgeMidPoint,
-            const label cellI,
-            const label faceI,
-            const bool faceOrder,
-            const label midPointI,
-            const label anchorPointI,
-            const label faceMidPointI,
-
-            Map<edge>& midPointToAnchors,
-            Map<edge>& midPointToFaceMids,
-            directTopoChange& meshMod
-        ) const;
-
-        //- Create all internal faces from an unsplit face
-        void createInternalFromSplitFace
-        (
-            const labelListList& cellAnchorPoints,
-            const labelListList& cellAddedCells,
-            const labelList& cellMidPoint,
-            const labelList& faceMidPoint,
-            const labelList& edgeMidPoint,
-            const label cellI,
-            const label faceI,
-
-            Map<edge>& midPointToAnchors,
-            Map<edge>& midPointToFaceMids,
-            directTopoChange& meshMod,
-            label& nFacesAdded
-        ) const;
-
-        //- Create all internal faces to split cellI into n cells where n is the
-        //  number of cell points
-        void createInternalFaces
-        (
-            const labelListList& cellAnchorPoints,
-            const labelListList& cellAddedCells,
-            const labelList& cellMidPoint,
-            const labelList& faceMidPoint,
-            const labelList& faceAnchorLevel,
-            const labelList& edgeMidPoint,
-            const label cellI,
-            directTopoChange& meshMod
-        ) const;
-
-        //- Store vertices from startFp upto face split point.
-        //  Used when splitting face into n faces where n is the number of
-        //  points in a face (or number of edges)
-        void walkFaceToMid
-        (
-            const labelList& edgeMidPoint,
-            const label cLevel,
-            const label faceI,
-            const label startFp,
-            dynamicLabelList& faceVerts
-        ) const;
-
-        //- Same as walkFaceToMid but now walk back
-        void walkFaceFromMid
-        (
-            const labelList& edgeMidPoint,
-            const label cLevel,
-            const label faceI,
-            const label startFp,
-            dynamicLabelList& faceVerts
-        ) const;
-
-        //- Updates refineCell so consistent 2:1 refinement. Returns local
-        //  number of cells changed
-        label faceConsistentRefinement
-        (
-            const bool maxSet,
-            PackedList<1>& refineCell
-        ) const;
-
-        //- Check wanted refinement for 2:1 consistency
-        void checkWantedRefinementLevels(const labelList&) const;
-
-
-        // Copy control
-
-            //- Disallow default bitwise copy construct
-            polyRef(const polyRef&);
-
-            //- Disallow default bitwise assignment
-            void operator=(const polyRef&);
-
-
-public:
-
-    //- Runtime type information
-    ClassName("polyRef");
-
-
-    // Constructors
-
-        //- Construct from mesh, read_if_present refinement data
-        //  (from write below)
-        polyRef(const polyMesh& mesh);
-
-        //- Construct from mesh and un/refinement data
-        polyRef
-        (
-            const polyMesh& mesh,
-            const labelList& cellLevel,
-            const labelList& pointLevel,
-            const polyRefinementHistory& history
-        );
-
-        //- Construct from mesh and refinement data.
-        polyRef
-        (
-            const polyMesh& mesh,
-            const labelList& cellLevel,
-            const labelList& pointLevel
-        );
-
-
-    // Member Functions
-
-        // Access
-
-            const labelIOList& cellLevel() const
-            {
-                return cellLevel_;
-            }
-
-            const labelIOList& pointLevel() const
-            {
-                return pointLevel_;
-            }
-
-            const polyRefinementHistory& history() const
-            {
-                return history_;
-            }
-
-            //- Typical edge length between unrefined points
-            scalar level0EdgeLength() const
-            {
-                return level0Edge_;
-            }
-
-        // Refinement
-
-            // Get cell level such that the face has the most points that are
-            // <= level (at least three points)
-            label getAnchorLevel(const label faceI) const;
-
-            //- Helper: get points of a cell without using cellPoints addressing
-            labelList cellPoints(const label cellI) const;
-
-            //- Given valid mesh and current cell level and proposed
-            //  cells to refine calculate any clashes (due to 2:1) and return
-            //  ok list of cells to refine.
-            //  Either adds cells to refine to set (maxSet = true) or
-            //  removes cells to refine (maxSet = false)
-            labelList consistentRefinement
-            (
-                const labelList& cellsToRefine,
-                const bool maxSet
-            ) const;
-
-            //- Like consistentRefinement but slower:
-            //  - specify number of cells between consecutive refinement levels
-            //    (consistentRefinement equivalent to 1)
-            //  - specify max level difference between point-connected cells.
-            //    (-1 to disable). Note that with normal 2:1 limitation
-            //    (maxFaceDiff=1) there can be 8:1 size difference across point
-            //    connected cells so maxPointDiff allows you to make that less.
-            //  cellsToRefine : cells we're thinking about refining. It will
-            //                  extend this set. All refinement levels will be
-            //                  at least maxFaceDiff layers thick.
-            //  facesToCheck  : additional faces where to implement the
-            //                  maxFaceDiff thickness (usually only boundary
-            //                  faces)
-            labelList consistentSlowRefinement
-            (
-                const label maxFaceDiff,
-                const labelList& cellsToRefine,
-                const labelList& facesToCheck,
-                const label maxPointDiff,
-                const labelList& pointsToCheck
-            ) const;
-
-            //- Like consistentSlowRefinement but uses different meshWave
-            //  (proper distance instead of toplogical count). No point checks
-            //  yet
-            labelList consistentSlowRefinement2
-            (
-                const label maxFaceDiff,
-                const labelList& cellsToRefine,
-                const labelList& facesToCheck
-            ) const;
-
-            //- Insert refinement. All selected cells will be split into n cells
-            //  where n is the number of points per cell.
-            //  Returns per element in cells the n cells they were split into.
-            //  Guarantees that the 0th element is the original cell label.
-            //  Mapping:
-            //  -split cells: n new ones get added from original
-            //  -split faces: original gets modified; n - 1 new ones get added
-            //                from original
-            //  -added internal faces: added from original cell face (if
-            //   that was internal) or created out-of-nothing (so will not
-            //   get mapped!). Note: could make this inflate from point but
-            //   that will allocate interpolation.
-            //  -points added to split edge: added from edge start()
-            //  -midpoints added: added from cellPoints[0].
-            labelListList setRefinement
-            (
-                const labelList& cells,
-                directTopoChange&
-            );
-
-            //- Update local numbering for changed mesh
-            void updateMesh(const mapPolyMesh&);
-
-
-            // Restoring : is where other processes delete and reinsert data.
-            //             These callbacks allow this to restore the cellLevel
-            //             and pointLevel for reintroduced points.
-            //             Is not related to undoing my refinement
-
-                //- Signal points/face/cells for which to store data
-                void storeData
-                (
-                    const labelList& pointsToStore,
-                    const labelList& facesToStore,
-                    const labelList& cellsToStore
-                );
-
-                //- Update local numbering + undo
-                //  Data to restore given as new pointlabel + stored pointlabel
-                //  (i.e. what was in pointsToStore)
-                void updateMesh
-                (
-                    const mapPolyMesh&,
-                    const Map<label>& pointsToRestore,
-                    const Map<label>& facesToRestore,
-                    const Map<label>& cellsToRestore
-                );
-
-
-            //- Update local numbering for subsetted mesh.
-            //  Gets new-to-old maps. Not compatible with unrefinement.
-            void subset
-            (
-                const labelList& pointMap,
-                const labelList& faceMap,
-                const labelList& cellMap
-            );
-
-            //- Update local numbering for mesh redistribution
-            void distribute(const mapDistributePolyMesh&);
-
-            //- Debug: Check coupled mesh for correctness
-            void checkMesh() const;
-
-            //- Debug: Check 2:1 consistency across faces.
-            //  maxPointDiff==-1 : only check 2:1 across faces
-            //  maxPointDiff!=-1 : check point-connected cells.
-            void checkRefinementLevels
-            (
-                const label maxPointDiff,
-                const labelList& pointsToCheck
-            ) const;
-
-
-        // Unrefinement (undoing refinement, not arbitrary coarsening)
-
-            //- Return the points at the centre of top-level split cells
-            //  that can be unsplit.
-            labelList getSplitPoints() const;
-
-            //- Given proposed splitPoints to unrefine according to calculate
-            //  any clashes (due to 2:1) and return ok list of points to
-            //  unrefine. Either adds points to refine to set (maxSet = true)
-            //  or removes points to refine (maxSet = false)
-            labelList consistentUnrefinement
-            (
-                const labelList& pointsToUnrefine,
-                const bool maxSet
-            ) const;
-
-            //- Remove some refinement. Needs to be supplied output of
-            //  consistentUnrefinement. Only call if undoable set.
-            //  All n pointCells of a split point will be combined into
-            //  the lowest numbered cell of those n.
-            void setUnrefinement
-            (
-                const labelList& splitPointLabels,
-                directTopoChange&
-            );
-
-
-        // Write
-
-            // Set instance for mesh files
-            void setInstance(const fileName& inst);
-
-            //- Force writing refinement+history to polyMesh directory.
-            bool write() const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/dynamicMesh/dynamicMesh/directTopoChange/directTopoChange/directActions/removeFaces.H

@@ -33,6 +33,15 @@ Description
 SourceFiles
     removeFaces.C
 
+Notes
+    All classes (hexRef8, undoableMeshCutter) that have a private member object
+    of this class (removeFaces) should be rewritten such that they are
+    polyMeshModifier classes. In order to allow for smooth transition between
+    using the directTopoChange engine towards polyTopoChange engine, the
+    member functions that actually do the refinement in this class are templated
+    on a Type which can either be polyTopoChange engine or directTopoChange
+    since they provide the same interface. VV, 16/Jan/2018.
+
 \*---------------------------------------------------------------------------*/
 
 #ifndef removeFaces_H
@@ -64,11 +73,11 @@ class removeFaces
 {
     // Private data
 
-        //- Reference to mesh
+        //- Const reference to mesh
         const polyMesh& mesh_;
 
         //- Cosine of angles between boundary faces. Boundary faces can be
-        //  merged only if angle between faces > minCos.
+        //  merged only if angle between faces is greater than minCos
         const scalar minCos_;
 
 
@@ -81,10 +90,12 @@ class removeFaces
             const label cellI,
             const label oldRegion,
             const label newRegion,
+
             labelList& cellRegion
         ) const;
 
-        //- Changes region of connected set of faces
+        //- Changes region of connected set of faces. Returns number of changed
+        //  faces
         label changeFaceRegion
         (
             const labelList& cellRegion,
@@ -92,11 +103,12 @@ class removeFaces
             const labelList& nFacesPerEdge,
             const label faceI,
             const label newRegion,
+
             labelList& faceRegion
         ) const;
 
-        //- Get all affected faces (including faces marked for removal)
+        //- Find and return all affected faces (including faces marked for removal)
-        boolList getFacesAffected
+        Xfer<boolList> affectedFaces
         (
             const labelList& cellRegion,
             const labelList& cellRegionMaster,
@@ -115,31 +127,28 @@ class removeFaces
                 const fileName&
             );
 
-            //- Merge faceLabels into single face.
+            //- Merge faceLabels into single face
+            template<class TopoChangeEngine>
             void mergeFaces
             (
                 const labelList& cellRegion,
                 const labelList& cellRegionMaster,
                 const labelHashSet& pointsToRemove,
                 const labelList& faceLabels,
-                directTopoChange& meshMod
+
+                TopoChangeEngine& ref
             ) const;
 
-            //- Get patch, zone info for faceI
+            //- Return face with all pointsToRemove removed
-            void getFaceInfo
+            face filterFace
             (
-                const label faceI,
+                const labelHashSet& pointsToRemove,
-                label& patchID,
+                const label
-                label& zoneID,
-                label& zoneFlip
             ) const;
 
-            //- Return face with all pointsToRemove removed.
-            face filterFace(const labelHashSet& pointsToRemove, const label)
-             const;
-
             //- Wrapper for meshMod.modifyFace. Reverses face if own>nei.
-            void modFace
+            template<class TopoChangeEngine>
+            void modifyFace
             (
                 const face& f,
                 const label masterFaceID,
@@ -151,7 +160,7 @@ class removeFaces
                 const label zoneID,
                 const bool zoneFlip,
 
-                directTopoChange& meshMod
+                TopoChangeEngine& ref
             ) const;
 
 
@@ -178,7 +187,7 @@ public:
 
     // Member Functions
 
-        //- Given set of faces to pierce calculates:
+        //- Given a set of faces to pierce calculates:
         //  - region for connected cells
         //  - mastercell for each region. This is the lowest numbered cell
         //    of all cells that get merged.
@@ -188,26 +197,30 @@ public:
         label compatibleRemoves
         (
             const labelList& inPiercedFaces,
+
             labelList& cellRegion,
             labelList& cellRegionMaster,
             labelList& outPiercedFaces
         ) const;
 
 
-        //- Play commands into directTopoChange to remove faces.
+        //- Play commands into TopoChangeEngine to remove faces.
+        template<class TopoChangeEngine>
         void setRefinement
         (
             const labelList& piercedFaces,
             const labelList& cellRegion,
             const labelList& cellRegionMaster,
-            directTopoChange&
+            TopoChangeEngine& ref
         ) const;
 
-        //- Force recalculation of locally stored data on topological change
+        //- Force recalculation of locally stored data on topological
+        //  change. Does nothing
         void updateMesh(const mapPolyMesh&)
         {}
 
-        //- Force recalculation of locally stored data for mesh distribution
+        //- Force recalculation of locally stored data for mesh
+        //  distribution. Does nothing.
         void distribute(const mapDistributePolyMesh&)
         {}
 };
@@ -217,6 +230,10 @@ public:
 
 } // End namespace Foam
 
+#ifdef NoRepository
+#   include "removeFacesTemplates.C"
+#endif
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #endif

src/dynamicMesh/dynamicMesh/polyMeshModifiers/polyhedralRefinement/polyhedralRefinement.H

@@ -0,0 +1,478 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::polyhedralRefinement
+
+Description
+    Isotropic refinement of polyhedral cells using the mesh modifier engine
+
+    Each polyhedral cell is split by the following procedure:
+    1. Adding points at the edge centre, face centre and cell centre,
+    2. Adding cells n cells where n is the number of points of the cell,
+    3. Splitting each face into multiple faces going from:
+       existing corner point -> new edge centre point -> new face centre
+       point -> other new edge centre point (sharing the same corner point)
+    4. Adding internal faces going from:
+       new edge centre point -> new face centre point -> new cell centre
+       point -> other new face centre point (sharing the same edge)
+
+SourceFiles
+    polyhedralRefinement.C
+
+Author
+    Vuko Vukcevic, Wikki Ltd.  All rights reserved.
+
+Notes
+    Generalisation of hexRef8 for polyhedral cells and refactorisation using
+    polyMesh modifier engine.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef polyhedralRefinement_H
+#define polyhedralRefinement_H
+
+#include "polyMeshModifier.H"
+#include "labelIOList.H"
+#include "removeFaces.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class polyhedralRefinement Declaration
+\*---------------------------------------------------------------------------*/
+
+class polyhedralRefinement
+:
+    public polyMeshModifier
+{
+public:
+
+    // Public enumeration for refinement status
+    enum refinementStatus
+    {
+        UNREFINED = -1,
+        UNCHANGED = 0,
+        REFINED = 1
+    };
+
+
+private:
+
+    // Private data
+
+        //- Reference to polyMesh for easy access in helper functions
+        const polyMesh& mesh_;
+
+
+        // Refinement control and handling
+
+            //- List of cells to refine in this time step
+            mutable labelList cellsToRefine_;
+
+            //- List of split point labels to unrefine in this time step
+            mutable labelList splitPointsToUnrefine_;
+
+            //- Cell refinement level
+            mutable labelIOList cellLevel_;
+
+            //- Point refinement level
+            mutable labelIOList pointLevel_;
+
+            //- Helper list for original (old) cells that will be refined or
+            //  unrefined. The list is updated in setPolyhedralRefinement and
+            //  setPolyhedralUnfereinement and is used in updateMesh to update
+            //  the cellLevel on the new mesh.
+            //  Values stored in the list are:
+            //  a) UNREFINED = -1 = cell is unrefined
+            //  b) UNCHANGED =  0 = cell is untouched
+            //  c) REFINED   = +1 = cell is refined
+            mutable labelList refinementLevelIndicator_;
+
+            //- Typical edge length between unrefined points
+            scalar level0EdgeLength_;
+
+            //- Face remover engine
+            mutable removeFaces faceRemover_;
+
+            //- Maximum number of cells in the mesh. Note: not strictly enforced
+            label maxCells_;
+
+            //- Maximum number of refinement levels for a given cell
+            label maxRefinementLevel_;
+
+            //- Switch whether to use point based consistency on refinement
+            Switch pointBasedConsistency_;
+
+            //- Number of buffer layers for refinement
+            label nBufferLayers_;
+
+
+    // Private Member Functions
+
+        // Helper functions
+
+            //- Get least cell level such that the face has at least three
+            //  points smaller than the level
+            label getAnchorLevel(const label faceI) const;
+
+            //- Calculate level0EdgeLength_ (constructor helper)
+            void calcLevel0EdgeLength();
+
+            //- Set file instance for cellLevel_ and pointLevel_
+            void setInstance(const fileName& inst) const;
+
+            //- Extend marked cells across faces
+            void extendMarkedCellsAcrossFaces(boolList& markedCell) const;
+
+            //- Extend marked cells across points
+            void extendMarkedCellsAcrossPoints(boolList& markedCell) const;
+
+
+        // Global topology modification functions (operate on whole polyMesh)
+        // and directly use local topology modification functions below
+
+            //- Refine polyhedral cells. All cellsToRefine_ cells will be split
+            //  into n cells where n is the number of points per cell.
+            //  Guarantees that the 0th element is the original cell label.
+            //  Mapping:
+            //  -split cells: n new ones get added from original
+            //  -split faces: original gets modified; n - 1 new ones get added
+            //                from original
+            //  -added internal faces: added from original cell face (if
+            //   that was internal) or created out-of-nothing (so will not
+            //   get mapped!). Note: could make this inflate from point but
+            //   that will allocate interpolation.
+            //  -points added to split edge: added from edge start()
+            //  -midpoints added: added from cellPoints[0].
+            void setPolyhedralRefinement(polyTopoChange& ref) const;
+
+            //- Remove some of the previously performed refinement. Uses
+            //  splitPointLabels_ to determine the unrefinement.
+            //  All n pointCells of a split point will be combined into
+            //  the lowest numbered cell of those n.
+            void setPolyhedralUnrefinement(polyTopoChange& ref) const;
+
+
+        // Local topology modification functions (operate on cells/faces)
+
+            //- Adds a face on top of existing faceI. Reverses if nessecary
+            label addFace
+            (
+                polyTopoChange& ref,
+                const label faceI,
+                const face& newFace,
+                const label own,
+                const label nei
+            ) const;
+
+            //- Adds internal face from point. No checks on reversal
+            label addInternalFace
+            (
+                polyTopoChange& ref,
+                const label meshFaceI,
+                const label meshPointI,
+                const face& newFace,
+                const label own,
+                const label nei
+            ) const;
+
+            //- Modifies existing faceI for either new owner/neighbour or new face
+            //  points. Reverses if nessecary
+            void modifyFace
+            (
+                polyTopoChange& ref,
+                const label faceI,
+                const face& newFace,
+                const label own,
+                const label nei
+            ) const;
+
+            //- Create all internal faces of split cellI into n cells where n is the
+            //  number of cell points
+            void createInternalFaces
+            (
+                const labelListList& cellAnchorPoints,
+                const labelListList& cellAddedCells,
+                const labelList& cellMidPoint,
+                const labelList& faceMidPoint,
+                const labelList& faceAnchorLevel,
+                const labelList& edgeMidPoint,
+                const label cellI,
+                polyTopoChange& ref
+            ) const;
+
+
+        // Topological change helper functions
+
+            //- Get cell added to point of cellI (if any)
+            label getAnchorCell
+            (
+                const labelListList& cellAnchorPoints,
+                const labelListList& cellAddedCells,
+                const label cellI,
+                const label faceI,
+                const label pointI
+            ) const;
+
+            //- Set new owner and neighbour (in unspecified order) of pointI
+            //  on faceI
+            void setNewFaceNeighbours
+            (
+                const labelListList& cellAnchorPoints,
+                const labelListList& cellAddedCells,
+                const label faceI,
+                const label pointI,
+
+                label& own,
+                label& nei
+            ) const;
+
+            //- Store vertices from startFp up to face split point.
+            //  Used when splitting face into n faces where n is the number of
+            //  points in a face (or number of edges)
+            void walkFaceToMid
+            (
+                const labelList& edgeMidPoint,
+                const label cLevel,
+                const label faceI,
+                const label startFp,
+                dynamicLabelList& faceVerts
+            ) const;
+
+            //- Same as walkFaceToMid but now walk back
+            void walkFaceFromMid
+            (
+                const labelList& edgeMidPoint,
+                const label cLevel,
+                const label faceI,
+                const label startFp,
+                dynamicLabelList& faceVerts
+            ) const;
+
+            //- Get index of point with minimum point level
+            label findMinLevel(const labelList& f) const;
+
+            //- Get index of point with maximum point level
+            label findMaxLevel(const labelList& f) const;
+
+            //- Count number of vertices <= anchorLevel for a given face
+            label countAnchors
+            (
+                const labelList& f,
+                const label anchorLevel
+            ) const;
+
+            //- Find index of point with wantedLevel, starting from fp
+            label findLevel
+            (
+                const face& f,
+                const label startFp,
+                const bool searchForward,
+                const label wantedLevel
+            ) const;
+
+            //- Store in maps correspondence from midpoint to anchors and
+            //  faces. Used when creating internal faces
+            label storeMidPointInfo
+            (
+                const labelListList& cellAnchorPoints,
+                const labelListList& cellAddedCells,
+                const labelList& cellMidPoint,
+                const labelList& edgeMidPoint,
+                const label cellI,
+                const label faceI,
+                const bool faceOrder,
+                const label midPointI,
+                const label anchorPointI,
+                const label faceMidPointI,
+
+                Map<edge>& midPointToAnchors,
+                Map<edge>& midPointToFaceMids,
+                polyTopoChange& ref
+            ) const;
+
+            //- If p0 and p1 are existing vertices check if edge is split and insert
+            //  splitPoint. Used with storing mid point
+            void insertEdgeSplit
+            (
+                const labelList& edgeMidPoint,
+                const label p0,
+                const label p1,
+                dynamicLabelList& verts
+            ) const;
+
+
+        // Debug functions
+
+            //- Check orientation of added internal face
+            void checkInternalOrientation
+            (
+                polyTopoChange& ref,
+                const label cellI,
+                const label faceI,
+                const point& ownPt,
+                const point& neiPt,
+                const face& newFace
+            ) const;
+
+            //- Check orientation of a new boundary face
+            void checkBoundaryOrientation
+            (
+                polyTopoChange& ref,
+                const label cellI,
+                const label faceI,
+                const point& ownPt,
+                const point& boundaryPt,
+                const face& newFace
+            ) const;
+
+
+        // Refinement/unrefinement consistency checks
+
+            //- Updates cellsToRefine such that a face consistent 2:1 refinement
+            //  is obtained. Returns local number of cells changed
+            label faceConsistentRefinement(boolList& cellsToRefine) const;
+
+            //- Updates cellsToRefine such that a point consistent 4:1 refinement
+            //  is obtained. Returns local number of cells changed
+            label pointConsistentRefinement(boolList& cellsToRefine) const;
+
+            //- Updates cellsToUnrefine such that a face consistent 2:1
+            //  unrefinement is obtained. Returns local number of cells changed
+            label faceConsistentUnrefinement(boolList& cellsToUnrefine) const;
+
+            //- Updates cellsToUnrefine such that a point consistent 4:1
+            //  unrefinement is obtained. Returns local number of cells changed
+            label pointConsistentUnrefinement(boolList& cellsToUnrefine) const;
+
+
+        // Copy control
+
+            //- Disallow default bitwise copy construct
+            polyhedralRefinement(const polyhedralRefinement&);
+
+            //- Disallow default bitwise assignment
+            void operator=(const polyhedralRefinement&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("polyhedralRefinement");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        polyhedralRefinement
+        (
+            const word& name,
+            const dictionary& dict,
+            const label index,
+            const polyTopoChanger& mme
+        );
+
+
+    // Destructor
+
+        virtual ~polyhedralRefinement();
+
+
+    // Member Functions
+
+        // Access
+
+            const labelIOList& cellLevel() const
+            {
+                return cellLevel_;
+            }
+
+            const labelIOList& pointLevel() const
+            {
+                return pointLevel_;
+            }
+
+            //- Typical edge length between unrefined points
+            scalar level0EdgeLength() const
+            {
+                return level0EdgeLength_;
+            }
+
+
+        // Edit
+
+            //- Set cells to refine given a list of refinement
+            //  candidates. Refinement candidates are extended within the
+            //  function due to possible 4:1 conflicts and specified number of
+            //  buffer layers.
+            //  Note: must be called BEFORE setSplitPointsToUnrefine
+            void setCellsToRefine(const labelList& refinementCellCandidates);
+
+            //- Set split points to unrefine given a list of all mesh points
+            //  that are candidates for unrefinement. Split points are
+            //  determined as a subset of unrefinement candidates, avoiding
+            //  splitting points of cells that are going to be refined at the
+            //  same time and ensuring consistent unrefinement.
+            //  Note: must be called AFTER setCellsToRefine
+            void setSplitPointsToUnrefine
+            (
+                const labelList& unrefinementPointCandidates
+            );
+
+
+        // Inherited interface from polyMeshModifier
+
+            //- Check for topology change
+            virtual bool changeTopology() const;
+
+            //- Insert the polyhedral refinement/unrefinement into the
+            //  topological change
+            virtual void setRefinement(polyTopoChange&) const;
+
+            //- Modify motion points to comply with the topological change
+            virtual void modifyMotionPoints(pointField& motionPoints) const;
+
+            //- Force recalculation of locally stored data on topological change
+            virtual void updateMesh(const mapPolyMesh&);
+
+            //- Write
+            virtual void write(Ostream&) const;
+
+            //- Write dictionary
+            virtual void writeDict(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/dynamicMesh/dynamicMesh/polyTopoChange/polyTopoChange/actions/removeObject/polyRemoveCell.H

@@ -56,6 +56,7 @@ class polyRemoveCell
         //- Merge cell ID or -1
         label mergeCellID_;
 
+
 public:
 
     // Static data members

src/dynamicMesh/dynamicMesh/refinementData/polyRefinementHistory.H

@@ -1,378 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | foam-extend: Open Source CFD
-   \\    /   O peration     | Version:     4.0
-    \\  /    A nd           | Web:         http://www.foam-extend.org
-     \\/     M anipulation  | For copyright notice see file Copyright
--------------------------------------------------------------------------------
-License
-    This file is part of foam-extend.
-
-    foam-extend 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.
-
-    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::polyRefinementHistory
-
-Description
-    All refinement history. Used in unrefinement.
-
-    - visibleCells: valid for the current mesh and contains per cell -1
-      (cell unrefined) or an index into splitCells_.
-    - splitCells: for every split contains the parent (also index into
-      splitCells) and optionally a subsplit as n indices into splitCells, where
-      n is the number of points of a cell. Note that the numbers in splitCells
-      are not cell labels, they are purely indices into splitCells.
-
-    E.g. 2 cells, cell 1 gets refined so end up with 9 cells:
-    @verbatim
-        // splitCells
-        9
-        (
-        -1 (1 2 3 4 5 6 7 8)
-        0 0()
-        0 0()
-        0 0()
-        0 0()
-        0 0()
-        0 0()
-        0 0()
-        0 0()
-        )
-
-        // visibleCells
-        9(-1 1 2 3 4 5 6 7 8)
-    @endverbatim
-
-
-    So cell0 (visibleCells=-1) is unrefined.
-    Cells 1-8 have all valid splitCells entries which are:
-      - parent:0
-      - subsplits:0()
-
-    The parent 0 refers back to the splitcell entries.
-
-
-SourceFiles
-    polyRefinementHistory.C
-
-Author
-    Vuko Vukcevic, Wikki Ltd.  All rights reserved
-
-Notes
-    Generalisation of refinementHistory for polyhedral cells
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef polyRefinementHistory_H
-#define polyRefinementHistory_H
-
-#include "DynamicList.H"
-#include "dynamicLabelList.H"
-#include "labelList.H"
-#include "FixedList.H"
-#include "SLList.H"
-#include "autoPtr.H"
-#include "regIOobject.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// Forward declaration of classes
-class mapPolyMesh;
-class mapDistributePolyMesh;
-
-/*---------------------------------------------------------------------------*\
-                     Class polyRefinementHistory Declaration
-\*---------------------------------------------------------------------------*/
-
-class polyRefinementHistory
-:
-    public regIOobject
-{
-public:
-
-    class splitPolyCell
-    {
-    public:
-
-        // Index to original splitCell this cell was refined off from
-        // -1: top level cell
-        // -2: free splitCell (so should also be in freeSplitCells_)
-        label parent_;
-
-        //- cells this cell was refined into
-        autoPtr<labelList> addedCellsPtr_;
-
-        //- Construct null (parent = -1)
-        splitPolyCell();
-
-        //- Construct from parent
-        splitPolyCell(const label parent);
-
-        //- Construct from Istream
-        splitPolyCell(Istream& is);
-
-        //- Construct as deep copy
-        splitPolyCell(const splitPolyCell&);
-
-        //- Copy operator since autoPtr otherwise 'steals' storage.
-        void operator=(const splitPolyCell& s)
-        {
-            // Check for assignment to self
-            if (this == &s)
-            {
-                FatalErrorIn
-                (
-                    "splitPolyCell::operator=(const Foam::splitPolyCell&)"
-                )   << "Attempted assignment to self"
-                    << abort(FatalError);
-            }
-
-            parent_ = s.parent_;
-
-            addedCellsPtr_.reset
-            (
-                s.addedCellsPtr_.valid()
-              ? new labelList(s.addedCellsPtr_())
-              : NULL
-            );
-        }
-
-        bool operator==(const splitPolyCell& s) const
-        {
-            if (addedCellsPtr_.valid() != s.addedCellsPtr_.valid())
-            {
-                return false;
-            }
-            else if (parent_ != s.parent_)
-            {
-                return false;
-            }
-            else if (addedCellsPtr_.valid())
-            {
-                return addedCellsPtr_() == s.addedCellsPtr_();
-            }
-            else
-            {
-                return true;
-            }
-        }
-
-        bool operator!=(const splitPolyCell& s) const
-        {
-            return !operator==(s);
-        }
-
-        friend Istream& operator>>(Istream&, splitPolyCell&);
-        friend Ostream& operator<<(Ostream&, const splitPolyCell&);
-    };
-
-
-private:
-
-    TypeName("polyRefinementHistory");
-
-    // Private data
-
-        //- Storage for splitCells
-        DynamicList<splitPolyCell> splitCells_;
-
-        //- Unused indices in splitCells
-        dynamicLabelList freeSplitCells_;
-
-        //- Currently visible cells. Indices into splitCells.
-        labelList visibleCells_;
-
-
-    // Private member functions
-
-        //- Debug write
-        static void writeEntry
-        (
-            const List<splitPolyCell>&,
-            const splitPolyCell&
-        );
-
-        //- Debug write
-        static void writeDebug
-        (
-            const labelList&,
-            const List<splitPolyCell>&
-        );
-
-        //- Check consistency of structure, i.e. indices into splitCells_.
-        void checkIndices() const;
-
-        //- Allocate a splitCell. Return index in splitCells_.
-        label allocateSplitCell
-        (
-            const label parent,
-            const label i,
-            const label nCells
-        );
-
-        //- Free a splitCell.
-        void freeSplitCell(const label index);
-
-        //- Mark entry in splitCells. Recursively mark its parent and subs.
-        void markSplit
-        (
-            const label,
-            labelList& oldToNew,
-            DynamicList<splitPolyCell>&
-        ) const;
-
-        void countProc
-        (
-            const label index,
-            const label newProcNo,
-            labelList& splitCellProc,
-            labelList& splitCellNum
-        ) const;
-
-public:
-
-    // Constructors
-
-        //- Construct (read) given an IOobject
-        polyRefinementHistory(const IOobject&);
-
-        //- Construct (read) or construct null
-        polyRefinementHistory
-        (
-            const IOobject&,
-            const List<splitPolyCell>& splitCells,
-            const labelList& visibleCells
-        );
-
-        //- Construct (read) or construct from initial number of cells
-        //  (all visible)
-        polyRefinementHistory(const IOobject&, const label nCells);
-
-        //- Construct as copy
-        polyRefinementHistory(const IOobject&, const polyRefinementHistory&);
-
-        //- Construct from Istream
-        polyRefinementHistory(const IOobject&, Istream&);
-
-
-    // Member Functions
-
-        //- Per cell in the current mesh (i.e. visible) either -1 (unrefined)
-        //  or an index into splitCells.
-        const labelList& visibleCells() const
-        {
-            return visibleCells_;
-        }
-
-        //- Storage for splitPolyCells.
-        const DynamicList<splitPolyCell>& splitCells() const
-        {
-            return splitCells_;
-        }
-
-        //- Cache of unused indices in splitCells
-        const dynamicLabelList& freeSplitCells() const
-        {
-            return freeSplitCells_;
-        }
-
-        //- Is there unrefinement history. Note that this will fall over if
-        //  there are 0 cells in the mesh. But this gives problems with
-        //  lots of other programs anyway.
-        bool active() const
-        {
-            return visibleCells_.size() > 0;
-        }
-
-        //- Get parent of cell
-        label parentIndex(const label cellI) const
-        {
-            label index = visibleCells_[cellI];
-
-            if (index < 0)
-            {
-                FatalErrorIn("polyRefinementHistory::parentIndex(const label)")
-                    << "Cell " << cellI << " is not visible"
-                    << abort(FatalError);
-            }
-            return splitCells_[index].parent_;
-        }
-
-        //- Store splitting of cell into n (number of points for a cell)
-        void storeSplit
-        (
-            const label cellI,
-            const labelList& addedCells
-        );
-
-        //- Store combining n cells into master
-        void combineCells
-        (
-            const label masterCellI,
-            const labelList& combinedCells
-        );
-
-        //- Update numbering for mesh changes
-        void updateMesh(const mapPolyMesh&);
-
-        //- Update numbering for subsetting
-        void subset
-        (
-            const labelList& pointMap,
-            const labelList& faceMap,
-            const labelList& cellMap
-        );
-
-        //- Update local numbering for mesh redistribution.
-        //  Can only distribute clusters sent across in one go; cannot
-        //  handle parts recombined in multiple passes.
-        void distribute(const mapDistributePolyMesh&);
-
-        //- Compact splitCells_. Removes all freeSplitCells_ elements.
-        void compact();
-
-        //- Extend/shrink storage. additional visibleCells_ elements get
-        //  set to -1.
-        void resize(const label nCells);
-
-        //- Debug write
-        void writeDebug() const;
-
-
-        //- ReadData function required for regIOobject read operation
-        virtual bool readData(Istream&);
-
-        //- WriteData function required for regIOobject write operation
-        virtual bool writeData(Ostream&) const;
-
-
-    // IOstream Operators
-
-        friend Istream& operator>>(Istream&, polyRefinementHistory&);
-        friend Ostream& operator<<(Ostream&, const polyRefinementHistory&);
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/dynamicMesh/topoChangerFvMesh/Make/files

@@ -11,4 +11,8 @@ mixerFvMesh/mixerFvMesh.C
 multiMixerFvMesh/mixerRotor.C
 multiMixerFvMesh/multiMixerFvMesh.C
 
+dynamicPolyRefinementFvMesh/dynamicPolyRefinementFvMesh.C
+dynamicPolyRefinementFvMesh/refinementSelection/refinementSelection/refinementSelection.C
+dynamicPolyRefinementFvMesh/refinementSelection/fieldBoundsRefinement/fieldBoundsRefinement.C
+
 LIB = $(FOAM_LIBBIN)/libtopoChangerFvMesh

src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicMeshDict

@@ -0,0 +1,58 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | foam-extend: Open Source CFD                    |
+|  \\    /   O peration     | Version:     4.0                                |
+|   \\  /    A nd           | Web:         http://www.foam-extend.org         |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      dynamicMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dynamicFvMesh   dynamicPolyRefinementFvMesh;
+
+dynamicPolyRefinementFvMeshCoeffs
+{
+    // Dynamic mesh procedure controls
+    refineInterval   1;
+
+    // Refinement selection criteria
+    refinementSelection
+    {
+        // Refines all cells with 0.001 < alpha < 0.999, otherwise unrefines
+        // previously refined cells
+        type        fieldBoundsRefinement;
+        fieldName   alpha;
+        lowerBound  0.001;
+        upperBound  0.999;
+
+        // Whether to use cell-point-cell smoothing for selecting refinement
+        // candidates. Off by default
+        cellPointCellSmoothing off;
+    }
+
+    // Polyhedral refinement engine controls
+    active yes;
+
+    // Maximum number of cells to allow (not strictly controlled)
+    maxCells             2000000;
+
+    // Maximum refinement level
+    maxRefinementLevel   3;
+
+    // Number of buffer layers between refinement levels
+    nBufferLayers        1;
+
+    // Whether to use point based consistency check. Needed when one allows more
+    // than 2 refinement levels (automatically switched on in that case)
+    pointBasedRefinement yes;
+}
+
+
+// ************************************************************************* //

src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicPolyRefinementFvMesh.C

@@ -0,0 +1,210 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "dynamicPolyRefinementFvMesh.H"
+#include "addToRunTimeSelectionTable.H"
+#include "refinementSelection.H"
+#include "polyhedralRefinement.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(dynamicPolyRefinementFvMesh, 0);
+
+addToRunTimeSelectionTable
+(
+    topoChangerFvMesh,
+    dynamicPolyRefinementFvMesh,
+    IOobject
+);
+
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::dynamicPolyRefinementFvMesh::readDict()
+{
+    // Read and check refinement interval
+    refineInterval_ = readLabel(refinementDict_.lookup("refineInterval"));
+    if (refineInterval_ < 1)
+    {
+        FatalErrorIn("dynamicPolyRefinementFvMesh::readDict()")
+            << "Illegal refineInterval found: " << refineInterval_ << nl
+            << "The refineInterval controls the refinement/unrefinement"
+            << " trigerring within a certain time step and should be > 0"
+            << exit(FatalError);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::dynamicPolyRefinementFvMesh::dynamicPolyRefinementFvMesh
+(
+    const IOobject& io
+)
+:
+    topoChangerFvMesh(io),
+    refinementDict_
+    (
+        IOdictionary
+        (
+            IOobject
+            (
+                "dynamicMeshDict",
+                time().constant(),
+                *this,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE,
+                false
+            )
+        ).subDict(typeName + "Coeffs")
+    ),
+    refineInterval_(readLabel(refinementDict_.lookup("refineInterval"))),
+    curTimeIndex_(-1),
+
+    refinementSelectionPtr_(refinementSelection::New(*this, refinementDict_))
+{
+    // Add the topology modifier engine
+    Info<< "Adding polyhedralRefinement topology modifier" << endl;
+
+    topoChanger_.setSize(1);
+    topoChanger_.set
+    (
+        0,
+        new polyhedralRefinement
+        (
+            "polyhedralRefinement",
+            refinementDict_,
+            0,
+            topoChanger_
+        )
+    );
+
+    // Write mesh and modifiers
+    topoChanger_.writeOpt() = IOobject::AUTO_WRITE;
+    topoChanger_.write();
+    write();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::dynamicPolyRefinementFvMesh::~dynamicPolyRefinementFvMesh()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::dynamicPolyRefinementFvMesh::update()
+{
+    // Re-read the data from dictionary for on-the-fly changes
+    readDict();
+
+    // Performing refinement/unrefinement when:
+    // 1. We are at the first time step
+    // 2. Time step is a multiplier of specified refineInterval
+    // 3. Only once per time step
+    if
+    (
+        time().timeIndex() > 0
+     && time().timeIndex() % refineInterval_ == 0
+     && curTimeIndex_ < time().timeIndex()
+    )
+    {
+        // Update current time index to skip multiple topo changes per single
+        // time step
+        curTimeIndex_ = time().timeIndex();
+
+        // Get reference to polyhedralRefinement polyMeshModifier
+        polyhedralRefinement& polyRefModifier =
+            refCast<polyhedralRefinement>(topoChanger_[0]);
+
+        // Get refinement candidates from refinement selection algorithm. Note:
+        // return type is Xfer<labelList> so there's no copying
+        const labelList refCandidates
+        (
+            refinementSelectionPtr_->refinementCellCandidates()
+        );
+
+        // Set cells to refine. Note: polyhedralRefinement ensures that face and
+        // point consistent refinement is performed
+        polyRefModifier.setCellsToRefine(refCandidates);
+
+        // Get unrefinement point candidates from refinement selection
+        // algorithm. Note: return type is Xfer<labelList> so there's no copying
+        const labelList unrefCandidates
+        (
+            refinementSelectionPtr_->unrefinementPointCandidates()
+        );
+
+        // Set split points to unrefine around.
+        // Notes:
+        // 1. polyhedralRefinement ensures that only a consistent set of split
+        //    points is used for unrefinement
+        // 2. Must be called after polyhedralRefinement::setCellsToRefine
+        polyRefModifier.setSplitPointsToUnrefine(unrefCandidates);
+
+        // Activate the polyhedral refinement engine if there are some cells to
+        // refine or there are some split points to unrefine around
+        bool enableTopoChange =
+            !refCandidates.empty() || !unrefCandidates.empty();
+
+        // Note: must enable topo change for all processors since face and point
+        // consistent refinement must be ensured across coupled patches
+        reduce(enableTopoChange, orOp<bool>());
+
+        if (enableTopoChange)
+        {
+            polyRefModifier.enable();
+        }
+        else
+        {
+            polyRefModifier.disable();
+        }
+
+        // Perform refinement and unrefinement in one go
+        autoPtr<mapPolyMesh> topoChangeMap = topoChanger_.changeMesh();
+
+        // Output cell balance
+        Info<< "Successfully performed polyhedral refinement. "
+            << "Changed from "
+            << returnReduce(topoChangeMap->nOldCells(), sumOp<label>())
+            << " to "
+            << returnReduce(topoChangeMap->cellMap().size(), sumOp<label>())
+            << " cells." << endl;
+
+        return topoChangeMap->morphing();
+    }
+
+    // No refinement/unrefinement at this time step. Return false
+    return false;
+}
+
+
+// ************************************************************************* //

src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicPolyRefinementFvMesh.H

@@ -0,0 +1,124 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::dynamicPolyRefinementFvMesh
+
+Description
+    Adaptive mesh refinement for arbitrary polyhedral cells.
+
+SourceFiles
+    dynamicPolyRefinementFvMesh.C
+
+Author
+    Vuko Vukcevic, Wikki Ltd.  All rights reserved
+
+Notes
+    Generalisation and refactorisation of dynamicRefineMesh for polyhedral cells
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef dynamicPolyRefinementFvMesh_H
+#define dynamicPolyRefinementFvMesh_H
+
+#include "topoChangerFvMesh.H"
+#include "refinementSelection.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                Class dynamicPolyRefinementFvMesh Declaration
+\*---------------------------------------------------------------------------*/
+
+class dynamicPolyRefinementFvMesh
+:
+    public topoChangerFvMesh
+{
+    // Private data
+
+        //- Dictionary containing dynamic mesh coefficient controls
+        dictionary refinementDict_;
+
+        //- Refinement interval
+        label refineInterval_;
+
+        //- Current time index (helper variable to skip multiple topo changes in
+        //  a single time step)
+        label curTimeIndex_;
+
+        //- Refinement selection algorithm that has two jobs:
+        //  1. Selects cells to refine from all cells
+        //  2. Selects split points that are ok to refine based on given set of
+        //     all split points
+        autoPtr<refinementSelection> refinementSelectionPtr_;
+
+
+    // Private Member Functions
+
+        //- Helper function for reading the dictionary and updating the data
+        void readDict();
+
+
+        // Copy control
+
+            //- Disallow default bitwise copy construct
+            dynamicPolyRefinementFvMesh(const dynamicPolyRefinementFvMesh&);
+
+            //- Disallow default bitwise assignment
+            void operator=(const dynamicPolyRefinementFvMesh&);
+
+public:
+
+    //- Runtime type information
+    TypeName("dynamicPolyRefinementFvMesh");
+
+
+    // Constructors
+
+        //- Construct from IOobject
+        explicit dynamicPolyRefinementFvMesh(const IOobject& io);
+
+
+    //- Destructor
+    virtual ~dynamicPolyRefinementFvMesh();
+
+
+    // Member Functions
+
+        //- Update the mesh for topology change
+        virtual bool update();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/refinementSelection/fieldBoundsRefinement/fieldBoundsRefinement.C

@@ -0,0 +1,187 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Author
+    Vuko Vukcevic, Wikki Ltd.  All rights reserved.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fieldBoundsRefinement.H"
+#include "addToRunTimeSelectionTable.H"
+#include "volFields.H"
+#include "pointFields.H"
+#include "volPointInterpolation.H"
+#include "pointVolInterpolation.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(fieldBoundsRefinement, 0);
+addToRunTimeSelectionTable
+(
+    refinementSelection,
+    fieldBoundsRefinement,
+    dictionary
+);
+
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::fieldBoundsRefinement::fieldBoundsRefinement
+(
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+:
+    refinementSelection(mesh, dict),
+    fieldName_(coeffDict().lookup("fieldName")),
+    lowerBound_(readScalar(coeffDict().lookup("lowerBound"))),
+    upperBound_(readScalar(coeffDict().lookup("upperBound"))),
+    cellPointCellSmoothing_
+    (
+        coeffDict().lookupOrDefault<Switch>("cellPointCellSmoothing", false)
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor* * * * * * * * * * * * * * * * //
+
+Foam::fieldBoundsRefinement::~fieldBoundsRefinement()
+{}
+
+
+// * * * * * * * * * * * * Public Member Functions * * * * * * * * * * * * * //
+
+Foam::Xfer<Foam::labelList>
+Foam::fieldBoundsRefinement::refinementCellCandidates() const
+{
+    // Get the field
+    const volScalarField& vField =
+        mesh().lookupObject<volScalarField>(fieldName_);
+
+    // Create temporary for the field (sharing the reference to volume field)
+    tmp<volScalarField> tvf(vField);
+
+    // Use cell-point-cell interpolation to smooth out the field
+    if (cellPointCellSmoothing_)
+    {
+        // Create volume to point interpolation object
+        const volPointInterpolation& vpi = volPointInterpolation::New(mesh());
+
+        // Interpolate from cell centres to points
+        pointScalarField pField(vpi.interpolate(vField));
+
+        // Create point to volume interpolation object
+        const pointMesh& pMesh = pointMesh::New(mesh());
+        const pointVolInterpolation pvi(pMesh, mesh());
+
+        // Interpolate from points back to cell centres
+        tmp<volScalarField> tInterpolatedVf = pvi.interpolate(pField);
+
+        // Assign to the temporary object
+        tvf = tInterpolatedVf;
+    }
+
+    // Get const reference to (possibly smoothed volume field)
+    const volScalarField& vf = tvf();
+
+    // Create storage for collection of cells. Assume that one in five cells
+    // will be refined to prevent excessive resizing.
+    dynamicLabelList refinementCandidates(mesh().nCells()/5);
+
+    // Loop through internal field and collect cells to refine
+    const scalarField& vfIn = vf.internalField();
+
+    forAll (vfIn, cellI)
+    {
+        // Get current cell value
+        const scalar& cellValue = vfIn[cellI];
+
+        if (cellValue > lowerBound_ && cellValue < upperBound_)
+        {
+            // Cell value is within the bounds, append cell for potential
+            // refinement
+            refinementCandidates.append(cellI);
+        }
+    }
+
+    // Print out some information
+    Info<< "Selection algorithm " << type() << " selected "
+        << returnReduce(refinementCandidates.size(), sumOp<label>())
+        << " cells as refinement candidates."
+        << endl;
+
+    // Return the list in the Xfer container to prevent copying
+    return refinementCandidates.xfer();
+}
+
+
+Foam::Xfer<Foam::labelList>
+Foam::fieldBoundsRefinement::unrefinementPointCandidates() const
+{
+    // Get the field
+    const volScalarField& vField =
+        mesh().lookupObject<volScalarField>(fieldName_);
+
+    // Create volume to point interpolation object
+    const volPointInterpolation& vpi = volPointInterpolation::New(mesh());
+
+    // Interpolate the volume field from cell centres to points and get the
+    // internal point field
+    pointScalarField pField(vpi.interpolate(vField));
+    const scalarField& pFieldIn = pField.internalField();
+
+    // Create storage for collection of candidates. Assume that one in ten
+    // mesh points will be unrefined to prevent excessive resizing
+    dynamicLabelList unrefinementCandidates(mesh().nPoints()/10);
+
+    // Loop through all split points and select candidates to unrefine
+    forAll (pField, pointI)
+    {
+        // Get point value
+        const scalar pointValue = pFieldIn[pointI];
+
+        if (pointValue > upperBound_ || pointValue < lowerBound_)
+        {
+            // Point value is outside of bounds, append point for potential
+            // unrefinement
+            unrefinementCandidates.append(pointI);
+        }
+    }
+
+    // Print out some information
+    Info<< "Selection algorithm " << type() << " selected "
+        << returnReduce(unrefinementCandidates.size(), sumOp<label>())
+        << " points as unrefinement candidates."
+        << endl;
+
+    // Return the list in the Xfer container to prevent copying
+    return unrefinementCandidates.xfer();
+}
+
+
+// ************************************************************************* //

src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/refinementSelection/fieldBoundsRefinement/fieldBoundsRefinement.H

@@ -0,0 +1,121 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::fieldBoundsRefinement
+
+Description
+    Selection of refinement cells based on a given scalar field bounds:
+    if a given (scalar) field is between lowerBound and upperBound, cell will
+    be marked for refinement.
+
+SourceFiles
+    fieldBoundsRefinement.C
+
+Author
+    Vuko Vukcevic, Wikki Ltd.  All rights reserved.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef fieldBoundsRefinement_H
+#define fieldBoundsRefinement_H
+
+#include "refinementSelection.H"
+#include "Switch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                   Class fieldBoundsRefinement Declaration
+\*---------------------------------------------------------------------------*/
+
+class fieldBoundsRefinement
+:
+    public refinementSelection
+{
+    // Private data
+
+        //- Name of the field to refine
+        word fieldName_;
+
+        //- Lower bound for the field
+        scalar lowerBound_;
+
+        //- Upper bound for the field
+        scalar upperBound_;
+
+        //- Whether to use cell-point-cell interpolation to smooth out the field
+        //  before selection. Switched off by default
+        Switch cellPointCellSmoothing_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        fieldBoundsRefinement(const fieldBoundsRefinement&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const fieldBoundsRefinement&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("fieldBoundsRefinement");
+
+
+    // Constructors
+
+        //- Construct from components
+        fieldBoundsRefinement(const fvMesh& mesh, const dictionary& dict);
+
+
+    //- Destructor
+    virtual ~fieldBoundsRefinement();
+
+
+    // Member Functions
+
+        // Selection of refinement/unrefinement candidates
+
+            //- Return transferable list of cells to refine
+            virtual Xfer<labelList> refinementCellCandidates() const;
+
+            //- Return transferable list of split points to unrefine
+            virtual Xfer<labelList> unrefinementPointCandidates() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/refinementSelection/refinementSelection/refinementSelection.C

@@ -0,0 +1,99 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Author
+    Vuko Vukcevic, Wikki Ltd.  All rights reserved.
+
+\*---------------------------------------------------------------------------*/
+
+#include "refinementSelection.H"
+#include "fvMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(refinementSelection, 0);
+defineRunTimeSelectionTable(refinementSelection, dictionary);
+
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::refinementSelection::refinementSelection
+(
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+:
+    mesh_(mesh),
+    coeffDict_(dict.subDict("refinementSelection"))
+{}
+
+
+// * * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::refinementSelection> Foam::refinementSelection::New
+(
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+{
+    // Get subdictionary from the dictionary
+    const dictionary coeffDict(dict.subDict("refinementSelection"));
+
+    // Get the name of the desired refinement selection algorithm
+    const word refinementSelectionTypeName(coeffDict.lookup("type"));
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(refinementSelectionTypeName);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+            (
+            "refinementSelection::refinementSelection::New\n"
+            "(\n"
+            "    const fvMesh& mesh,\n"
+            "    const dictionary& dict\n"
+            ")"
+            )   << "Unknown refinementSelection type "
+                << refinementSelectionTypeName << endl << endl
+                << "Valid refinementSelection types are :" << endl
+                << dictionaryConstructorTablePtr_->toc()
+                << exit(FatalError);
+    }
+
+    return autoPtr<refinementSelection>(cstrIter()(mesh, dict));
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor* * * * * * * * * * * * * * * * //
+
+Foam::refinementSelection::~refinementSelection()
+{}
+
+
+// ************************************************************************* //

src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/refinementSelection/refinementSelection/refinementSelection.H

@@ -0,0 +1,173 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::refinementSelection
+
+Description
+    An abstract base class for providing an interface for refinement selection
+    algorithms. The interface provides two functionalities:
+
+    1. Selects candidate cells to refine based on chosen criteria through
+       refinementCellCandidates() member function. Returns a list of all cells
+       that satisfy given criteria,
+    2. Selects candidate split points to unrefine based on chosen criteria
+       through unrefinementPointCandidates() member function. Returns a list of
+       all points (not just split points!) that satisfy the criteria.
+
+    Note: Here, we do not check that the unrefinement point candidates do not
+    clash the points of cells marked as candidates for refinement. This is taken
+    care of in polyhedralRefinement class.
+
+SourceFiles
+    refinementSelection.C
+
+Author
+    Vuko Vukcevic, Wikki Ltd.  All rights reserved.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef refinementSelection_H
+#define refinementSelection_H
+
+#include "dictionary.H"
+#include "Xfer.H"
+#include "labelList.H"
+#include "autoPtr.H"
+#include "runTimeSelectionTables.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class fvMesh;
+
+/*---------------------------------------------------------------------------*\
+                   Class refinementSelection Declaration
+\*---------------------------------------------------------------------------*/
+
+class refinementSelection
+{
+    // Private data
+
+        //- Const reference to fvMesh
+        const fvMesh& mesh_;
+
+        //- Refinement selection dictionary (subdictionary of
+        //  dynamicPolyRefinementFvMeshCoeffs dictionary)
+        const dictionary coeffDict_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        refinementSelection(const refinementSelection&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const refinementSelection&);
+
+
+protected:
+
+    // Protected member functions
+
+        // Access functions for derived classes
+
+            //- Const access to the fvMesh
+            const fvMesh& mesh() const
+            {
+                return mesh_;
+            }
+
+            //- Const access to the coefficient dictionary
+            const dictionary& coeffDict() const
+            {
+                return coeffDict_;
+            }
+
+
+public:
+
+    //- Runtime type information
+    TypeName("refinementSelection");
+
+
+    // Declare run-time constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            refinementSelection,
+            dictionary,
+            (
+                const fvMesh& mesh,
+                const dictionary& dict
+            ),
+            (mesh, dict)
+        );
+
+
+    // Constructors
+
+        //- Construct from components
+        refinementSelection(const fvMesh& mesh, const dictionary& dict);
+
+
+    // Selectors
+
+        //- Return an autoPtr to the selected refinementSelection
+        static autoPtr<refinementSelection> New
+        (
+            const fvMesh& mesh,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~refinementSelection();
+
+
+    // Member Functions
+
+        // Selection of refinement/unrefinement candidates
+
+            //- Return transferable list of refinement cell candidates
+            virtual Xfer<labelList> refinementCellCandidates() const = 0;
+
+            //- Return transferable list of unrefinement split point candidates
+            virtual Xfer<labelList> unrefinementPointCandidates() const = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteArea/faMesh/faMesh.C

@@ -262,7 +262,7 @@ Foam::faMesh::faMesh(const polyMesh& pMesh)
     // Calculate the geometry for the patches (transformation tensors etc.)
     boundary_.calcGeometry();
 
-    if (isFile(pMesh.time().timePath()/pMesh.dbDir()/"S0"))
+    if (isFile(pMesh.time().timePath()/"S0"))
     {
         S0Ptr_ = new DimensionedField<scalar, areaMesh>
         (
@@ -811,7 +811,7 @@ Foam::faMesh::faMesh
     // Calculate the geometry for the patches (transformation tensors etc.)
     boundary_.calcGeometry();
 
-    if (isFile(mesh().time().timePath()/mesh().dbDir()/"S0"))
+    if (isFile(mesh().time().timePath()/"S0"))
     {
         S0Ptr_ = new DimensionedField<scalar, areaMesh>
         (
@@ -1292,7 +1292,7 @@ bool Foam::faMesh::movePoints() const
             if (debug)
             {
                 InfoIn("bool faMesh::movePoints() const")
-                    << "Creating old cell volumes." << endl;
+                    << "Creating old face areas." << endl;
             }
 
             S0Ptr_ = new DimensionedField<scalar, areaMesh>

src/finiteArea/faMesh/faMeshUpdate.C

@@ -237,7 +237,6 @@ void Foam::faMesh::mapOldAreas(const faMeshMapper& mapper) const
             }
         }
     }
-
 }
 
 

src/finiteArea/include/createFaMesh.H

@@ -1,2 +1,2 @@
 
-    faMesh aMesh(mesh);
+    const faMesh& aMesh = faMesh::New(mesh);

src/finiteVolume/fields/fvPatchFields/constraint/ggi/ggiFvPatchField.C

@@ -178,17 +178,28 @@ tmp<Field<Type> > ggiFvPatchField<Type>::patchNeighbourField() const
 
     if (ggiPatch_.bridgeOverlap())
     {
-        // Symmetry treatment used for overlap
+        // Use mirrored neighbour field for interpolation. Note: mirroring needs
-        const vectorField nHat = this->patch().nf();
+        // to take into account the weights, i.e. how "far" we are actually
-
+        // mirroring. VV, 19/Jan/2018.
-        // Use mirrored neighbour field for interpolation
-        // HJ, 21/Jan/2009
         const Field<Type> bridgeField =
-            transform(I - 2.0*sqr(nHat), this->patchInternalField());
+        transform
+        (
+            (I - sqr(this->patch().nf())/(1.0 - ggiPatch_.fvPatch::weights())),
+            this->patchInternalField()
+        );
 
-        // Note: bridging now takes into account fully uncovered and partially
+//        if (pTraits<Type>::rank == 0)
-        // covered faces. VV, 18/Oct/2017.
+//        {
-        ggiPatch_.bridge(bridgeField, pnf);
+//            // Scale the field for scalars to ensure conservative and consistent
+//            // flux on both sides
+//            ggiPatch_.scaleForPartialCoverage(bridgeField, pnf);
+//        }
+//        else
+        {
+            // Bridge the field for higher order tensors to correctly take into
+            // account mirroring
+            ggiPatch_.bridge(bridgeField, pnf);
+        }
     }
 
     return tpnf;
@@ -259,15 +270,19 @@ void ggiFvPatchField<Type>::initInterfaceMatrixUpdate
 
     if (ggiPatch_.bridgeOverlap())
     {
-        // Note: will not work properly for types with rank > 0 (everything
+        // Note: this implicit treatment does not really work implicitly for
-        // above scalar) if the symmetry plane is not aligned with one of the
+        // types with rank > 0 (everything above scalar) if the symmetry plane
-        // coordinate axes. VV, 18/Oct/2017.
+        // is not aligned with one of the coordinate axes. VV, 18/Oct/2017.
+
+        // Use mirrored neighbour field for interpolation. Note: mirroring needs
+        // to take into account the weights, i.e. how "far" we are actually
+        // mirroring. VV, 19/Jan/2018.
         const scalarField bridgeField =
-            transform
+        transform
-            (
+        (
-                I - 2.0*sqr(this->patch().nf()),
+            (I - sqr(this->patch().nf())/(1.0 - ggiPatch_.fvPatch::weights())),
-                ggiPatch_.patchInternalField(psiInternal)
+            ggiPatch_.patchInternalField(psiInternal)
-            );
+        );
 
         ggiPatch_.bridge(bridgeField, pnf);
     }
@@ -335,15 +350,15 @@ void ggiFvPatchField<Type>::initInterfaceMatrixUpdate
 
     if (ggiPatch_.bridgeOverlap())
     {
-        // Note: will not work properly for types with rank > 0 (everything
+        // Use mirrored neighbour field for interpolation. Note: mirroring needs
-        // above scalar) if the symmetry plane is not aligned with one of the
+        // to take into account the weights, i.e. how "far" we are actually
-        // coordinate axes. VV, 18/Oct/2017.
+        // mirroring. VV, 19/Jan/2018.
         const Field<Type> bridgeField =
-            transform
+        transform
-            (
+        (
-                I - 2.0*sqr(this->patch().nf()),
+            (I - sqr(this->patch().nf())/(1.0 - ggiPatch_.fvPatch::weights())),
-                ggiPatch_.patchInternalField(psiInternal)
+            ggiPatch_.patchInternalField(psiInternal)
-            );
+        );
 
         ggiPatch_.bridge(bridgeField, pnf);
     }

src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C

@@ -660,7 +660,35 @@ CoEulerDdtScheme<Type>::fvcDdtConsistentPhiCorr
     const surfaceScalarField& rAUf
 )
 {
-    return (mesh().Sf() & faceU.oldTime())*rAUf*CofrDeltaT();
+    tmp<fluxFieldType> toldTimeFlux =
+        (mesh().Sf() & faceU.oldTime())*rAUf*CofrDeltaT();
+
+    if (mesh().moving())
+    {
+        // Mesh is moving, need to take into account the ratio between old and
+        // current cell volumes
+        volScalarField V0ByV
+        (
+            IOobject
+            (
+                "V0ByV",
+                mesh().time().timeName(),
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            dimensionedScalar("one", dimless, 1.0),
+            zeroGradientFvPatchScalarField::typeName
+        );
+        V0ByV.internalField() = mesh().V0()/mesh().V();
+        V0ByV.correctBoundaryConditions();
+
+        // Correct the flux with interpolated volume ratio
+        toldTimeFlux() *= fvc::interpolate(V0ByV);
+    }
+
+    return toldTimeFlux;
 }
 
 

src/finiteVolume/finiteVolume/ddtSchemes/CrankNicolsonDdtScheme/CrankNicolsonDdtScheme.C

@@ -1202,15 +1202,38 @@ CrankNicolsonDdtScheme<Type>::fvcDdtConsistentPhiCorr
         );
     }
     
-    return
+    // Calculate old time flux
-        rAUf*
+    fluxFieldType oldTimeFlux =
+        rAUf*rDtCoef_(faceUDdt0)*(mesh().Sf() & faceU.oldTime());
+
+    if (mesh().moving())
+    {
+        // Mesh is moving, need to take into account the ratio between old and
+        // current cell volumes
+        volScalarField V0ByV
         (
-            mesh().Sf()
+            IOobject
-          & (
+            (
-                rDtCoef_(faceUDdt0)*faceU.oldTime()
+                "V0ByV",
-              + offCentre_(faceUDdt0())
+                mesh().time().timeName(),
-            )
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            dimensionedScalar("one", dimless, 1.0),
+            zeroGradientFvPatchScalarField::typeName
         );
+        V0ByV.internalField() = mesh().V0()/mesh().V();
+        V0ByV.correctBoundaryConditions();
+
+        // Correct the flux with interpolated volume ratio
+        oldTimeFlux *= fvc::interpolate(V0ByV);
+    }
+
+    return
+        oldTimeFlux
+      + rAUf*rDtCoef_(faceUDdt0)*(mesh().Sf() & offCentre_(faceUDdt0()));
 }
 
 

src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C

@@ -522,7 +522,35 @@ EulerDdtScheme<Type>::fvcDdtConsistentPhiCorr
     const surfaceScalarField& rAUf
 )
 {
-    return (mesh().Sf() & faceU.oldTime())*rAUf/mesh().time().deltaT();
+    tmp<fluxFieldType> toldTimeFlux =
+        (mesh().Sf() & faceU.oldTime())*rAUf/mesh().time().deltaT();
+
+    if (mesh().moving())
+    {
+        // Mesh is moving, need to take into account the ratio between old and
+        // current cell volumes
+        volScalarField V0ByV
+        (
+            IOobject
+            (
+                "V0ByV",
+                mesh().time().timeName(),
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            dimensionedScalar("one", dimless, 1.0),
+            zeroGradientFvPatchScalarField::typeName
+        );
+        V0ByV.internalField() = mesh().V0()/mesh().V();
+        V0ByV.correctBoundaryConditions();
+
+        // Correct the flux with interpolated volume ratio
+        toldTimeFlux() *= fvc::interpolate(V0ByV);
+    }
+
+    return toldTimeFlux;
 }
 
 

src/finiteVolume/finiteVolume/ddtSchemes/SLTSDdtScheme/SLTSDdtScheme.C

@@ -665,7 +665,35 @@ SLTSDdtScheme<Type>::fvcDdtConsistentPhiCorr
     const surfaceScalarField& rAUf
 )
 {
-    return (mesh().Sf() & faceU.oldTime())*rAUf*fvc::interpolate(SLrDeltaT());
+    tmp<fluxFieldType> toldTimeFlux =
+        (mesh().Sf() & faceU.oldTime())*rAUf*fvc::interpolate(SLrDeltaT());
+
+    if (mesh().moving())
+    {
+        // Mesh is moving, need to take into account the ratio between old and
+        // current cell volumes
+        volScalarField V0ByV
+        (
+            IOobject
+            (
+                "V0ByV",
+                mesh().time().timeName(),
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            dimensionedScalar("one", dimless, 1.0),
+            zeroGradientFvPatchScalarField::typeName
+        );
+        V0ByV.internalField() = mesh().V0()/mesh().V();
+        V0ByV.correctBoundaryConditions();
+
+        // Correct the flux with interpolated volume ratio
+        toldTimeFlux() *= fvc::interpolate(V0ByV);
+    }
+
+    return toldTimeFlux;
 }
 
 

src/finiteVolume/finiteVolume/ddtSchemes/backwardDdtScheme/backwardDdtScheme.C

@@ -732,18 +732,66 @@ backwardDdtScheme<Type>::fvcDdtConsistentPhiCorr
 
     const scalar rDeltaT = 1.0/deltaT;
 
+    // Note: minus sign in gamma coefficient so we can simply add the fluxes
+    // together at the end
     const dimensionedScalar beta("beta", dimless/dimTime, coefft0*rDeltaT);
     const dimensionedScalar gamma("gamma", dimless/dimTime, -coefft00*rDeltaT);
 
-    return
+    // Calculate old and old-old flux contributions
-        rAUf*
+    fluxFieldType oldTimeFlux =
+        beta*rAUf*(mesh().Sf() & faceU.oldTime());
+    fluxFieldType oldOldTimeFlux =
+        gamma*rAUf*(mesh().Sf() & faceU.oldTime().oldTime());
+
+    if (mesh().moving())
+    {
+        // Mesh is moving, need to take into account the ratio between old and
+        // current cell volumes for old flux contribution
+        volScalarField V0ByV
         (
-            mesh().Sf()
+            IOobject
-          & (
+            (
-                beta*faceU.oldTime()
+                "V0ByV",
-              + gamma*faceU.oldTime().oldTime()
+                mesh().time().timeName(),
-            )
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            dimensionedScalar("one", dimless, 1.0),
+            zeroGradientFvPatchScalarField::typeName
         );
+        V0ByV.internalField() = mesh().V0()/mesh().V();
+        V0ByV.correctBoundaryConditions();
+
+        // Correct old time flux contribution
+        oldTimeFlux *= fvc::interpolate(V0ByV);
+
+
+        // Also need to take into account the ratio between old-old and current
+        // cell volumes for old-old time flux contribution
+        volScalarField V00ByV
+        (
+            IOobject
+            (
+                "V00ByV",
+                mesh().time().timeName(),
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            dimensionedScalar("one", dimless, 1.0),
+            zeroGradientFvPatchScalarField::typeName
+        );
+        V00ByV.internalField() = mesh().V00()/mesh().V();
+        V00ByV.correctBoundaryConditions();
+
+        // Correct old-old time flux contribution
+        oldOldTimeFlux *= fvc::interpolate(V00ByV);
+    }
+
+    return oldTimeFlux + oldOldTimeFlux;
 }
 
 

src/finiteVolume/finiteVolume/ddtSchemes/steadyInertialDdtScheme/steadyInertialDdtScheme.C

@@ -667,7 +667,35 @@ steadyInertialDdtScheme<Type>::fvcDdtConsistentPhiCorr
     const surfaceScalarField& rAUf
 )
 {
-    return (mesh().Sf() & faceU.oldTime())*rAUf*CofrDeltaT();
+    tmp<fluxFieldType> toldTimeFlux =
+        (mesh().Sf() & faceU.oldTime())*rAUf*CofrDeltaT();
+
+    if (mesh().moving())
+    {
+        // Mesh is moving, need to take into account the ratio between old and
+        // current cell volumes
+        volScalarField V0ByV
+        (
+            IOobject
+            (
+                "V0ByV",
+                mesh().time().timeName(),
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            dimensionedScalar("one", dimless, 1.0),
+            zeroGradientFvPatchScalarField::typeName
+        );
+        V0ByV.internalField() = mesh().V0()/mesh().V();
+        V0ByV.correctBoundaryConditions();
+
+        // Correct the flux with interpolated volume ratio
+        toldTimeFlux() *= fvc::interpolate(V0ByV);
+    }
+
+    return toldTimeFlux;
 }