Merge branch 'develop' of develop.openfoam.com:Development/OpenFOAM-plus into develop

Conflicts:
	tutorials/basic/overLaplacianDyMFoam/heatTransfer/0.orig/T
	tutorials/basic/overLaplacianDyMFoam/heatTransfer/0.orig/zoneID

applications/solvers/multiphase/interIsoFoam/Make/files

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

applications/solvers/multiphase/interIsoFoam/Make/options

@@ -0,0 +1,20 @@
+EXE_INC = \
+    -I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/transportModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
+
+EXE_LIBS = \
+    -limmiscibleIncompressibleTwoPhaseMixture \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
+    -lfiniteVolume \
+    -lfvOptions \
+    -lmeshTools \
+    -lsampling

applications/solvers/multiphase/interIsoFoam/UEqn.H

@@ -0,0 +1,33 @@
+    MRF.correctBoundaryVelocity(U);
+
+    fvVectorMatrix UEqn
+    (
+        fvm::ddt(rho, U) + fvm::div(rhoPhi, U)
+      + MRF.DDt(rho, U)
+      + turbulence->divDevRhoReff(rho, U)
+     ==
+        fvOptions(rho, U)
+    );
+
+    UEqn.relax();
+
+    fvOptions.constrain(UEqn);
+
+    if (pimple.momentumPredictor())
+    {
+        solve
+        (
+            UEqn
+         ==
+            fvc::reconstruct
+            (
+                (
+                    mixture.surfaceTensionForce()
+                  - ghf*fvc::snGrad(rho)
+                  - fvc::snGrad(p_rgh)
+                ) * mesh.magSf()
+            )
+        );
+
+        fvOptions.correct(U);
+    }

applications/solvers/multiphase/interIsoFoam/alphaControls.H

@@ -0,0 +1,3 @@
+const dictionary& alphaControls = mesh.solverDict(alpha1.name());
+
+label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")));

applications/solvers/multiphase/interIsoFoam/alphaCourantNo.H

@@ -0,0 +1,57 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Global
+    alphaCourantNo
+
+Description
+    Calculates and outputs the mean and maximum Courant Numbers.
+
+\*---------------------------------------------------------------------------*/
+
+scalar maxAlphaCo
+(
+    readScalar(runTime.controlDict().lookup("maxAlphaCo"))
+);
+
+scalar alphaCoNum = 0.0;
+scalar meanAlphaCoNum = 0.0;
+
+if (mesh.nInternalFaces())
+{
+    scalarField sumPhi
+    (
+        mixture.nearInterface()().primitiveField()
+       *fvc::surfaceSum(mag(phi))().primitiveField()
+    );
+
+    alphaCoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
+
+    meanAlphaCoNum =
+        0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();
+}
+
+Info<< "Interface Courant Number mean: " << meanAlphaCoNum
+    << " max: " << alphaCoNum << endl;
+
+// ************************************************************************* //

applications/solvers/multiphase/interIsoFoam/alphaEqn.H

@@ -0,0 +1,35 @@
+// If there are more than one outer corrector, we use a mixture of old and
+// new U and phi for propagating alpha1 in all but the first outer iteration
+if (!pimple.firstIter())
+{
+    // We are recalculating alpha1 from the its old time value
+    alpha1 = alpha1.oldTime();
+    // Temporarily storing new U and phi values in prevIter storage
+    U.storePrevIter();
+    phi.storePrevIter();
+
+    // Overwriting new U and phi values with mixture of old and new values
+    phi = 0.5*(phi + phi.oldTime());
+    U = 0.5*(U + U.oldTime());
+}
+
+// Update alpha1
+advector.advect();
+
+// Update rhoPhi
+rhoPhi = advector.getRhoPhi(rho1, rho2);
+
+alpha2 = 1.0 - alpha1;
+
+if (!pimple.firstIter())
+{
+    // Restoring new U and phi values temporarily saved in prevIter() above
+    U = U.prevIter();
+    phi = phi.prevIter();
+}
+
+Info<< "Phase-1 volume fraction = "
+    << alpha1.weightedAverage(mesh.Vsc()).value()
+    << "  Min(" << alpha1.name() << ") = " << min(alpha1).value()
+    << "  Max(" << alpha1.name() << ") - 1 = " << max(alpha1).value() - 1
+    << endl;

applications/solvers/multiphase/interIsoFoam/alphaEqnSubCycle.H

@@ -0,0 +1,35 @@
+if (nAlphaSubCycles > 1)
+{
+    dimensionedScalar totalDeltaT = runTime.deltaT();
+    surfaceScalarField rhoPhiSum
+    (
+        IOobject
+        (
+            "rhoPhiSum",
+            runTime.timeName(),
+            mesh
+        ),
+        mesh,
+        dimensionedScalar("0", rhoPhi.dimensions(), 0)
+    );
+
+    tmp<volScalarField> trSubDeltaT;
+
+    for
+    (
+        subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
+        !(++alphaSubCycle).end();
+    )
+    {
+        #include "alphaEqn.H"
+        rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
+    }
+
+    rhoPhi = rhoPhiSum;
+}
+else
+{
+    #include "alphaEqn.H"
+}
+
+rho == alpha1*rho1 + alpha2*rho2;

applications/solvers/multiphase/interIsoFoam/correctPhi.H

@@ -0,0 +1,11 @@
+CorrectPhi
+(
+    U,
+    phi,
+    p_rgh,
+    dimensionedScalar("rAUf", dimTime/rho.dimensions(), 1),
+    geometricZeroField(),
+    pimple
+);
+
+#include "continuityErrs.H"

applications/solvers/multiphase/interIsoFoam/createFields.H

@@ -0,0 +1,123 @@
+Info<< "Reading field p_rgh\n" << endl;
+volScalarField p_rgh
+(
+    IOobject
+    (
+        "p_rgh",
+        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;
+immiscibleIncompressibleTwoPhaseMixture mixture(U, phi);
+
+volScalarField& alpha1(mixture.alpha1());
+volScalarField& alpha2(mixture.alpha2());
+
+const dimensionedScalar& rho1 = mixture.rho1();
+const dimensionedScalar& rho2 = mixture.rho2();
+
+
+// Need to store rho for ddt(rho, U)
+volScalarField rho
+(
+    IOobject
+    (
+        "rho",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT
+    ),
+    alpha1*rho1 + alpha2*rho2
+);
+rho.oldTime();
+
+
+// Mass flux
+surfaceScalarField rhoPhi
+(
+    IOobject
+    (
+        "rhoPhi",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    fvc::interpolate(rho)*phi
+);
+
+
+// Construct incompressible turbulence model
+autoPtr<incompressible::turbulenceModel> turbulence
+(
+    incompressible::turbulenceModel::New(U, phi, mixture)
+);
+
+
+#include "readGravitationalAcceleration.H"
+#include "readhRef.H"
+#include "gh.H"
+
+
+volScalarField p
+(
+    IOobject
+    (
+        "p",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    p_rgh + rho*gh
+);
+
+label pRefCell = 0;
+scalar pRefValue = 0.0;
+setRefCell
+(
+    p,
+    p_rgh,
+    pimple.dict(),
+    pRefCell,
+    pRefValue
+);
+
+if (p_rgh.needReference())
+{
+    p += dimensionedScalar
+    (
+        "p",
+        p.dimensions(),
+        pRefValue - getRefCellValue(p, pRefCell)
+    );
+    p_rgh = p - rho*gh;
+}
+
+mesh.setFluxRequired(p_rgh.name());
+mesh.setFluxRequired(alpha1.name());
+
+#include "createMRF.H"
+#include "createIsoAdvection.H"

applications/solvers/multiphase/interIsoFoam/createIsoAdvection.H

@@ -0,0 +1,10 @@
+// Defining isoAdvection
+isoAdvection advector(alpha1, phi, U);
+
+bool frozenFlow = pimple.dict().lookupOrDefault<bool>("frozenFlow", false);
+if (frozenFlow)
+{
+    Info<< "Employing frozen-flow assumption: "
+        << "pressure-velocity system will not be updated"
+        << endl;
+}

applications/solvers/multiphase/interIsoFoam/interIsoFoam.C

@@ -0,0 +1,143 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+                isoAdvector | Copyright (C) 2016 DHI
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    interIsoFoam
+
+Group
+    grpMultiphaseSolvers
+
+Description
+    Solver derived from interFoam for 2 incompressible, isothermal immiscible
+    fluids using the iso-advector phase-fraction based interface capturing
+    approach.
+
+    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.
+
+    Reference:
+        \verbatim
+            Roenby, J., Bredmose, H. and Jasak, H. (2016).
+            A computational method for sharp interface advection
+            Royal Society Open Science, 3
+            doi 10.1098/rsos.160405
+        \endverbatim
+
+    isoAdvector code supplied by Johan Roenby, DHI (2016)
+
+\*---------------------------------------------------------------------------*/
+
+#include "isoAdvection.H"
+#include "fvCFD.H"
+#include "subCycle.H"
+#include "immiscibleIncompressibleTwoPhaseMixture.H"
+#include "turbulentTransportModel.H"
+#include "pimpleControl.H"
+#include "fvOptions.H"
+#include "CorrectPhi.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "postProcess.H"
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createControl.H"
+    #include "createTimeControls.H"
+    #include "initContinuityErrs.H"
+    #include "createFields.H"
+    #include "createFvOptions.H"
+    #include "correctPhi.H"
+
+    turbulence->validate();
+
+    #include "readTimeControls.H"
+    #include "CourantNo.H"
+    #include "setInitialDeltaT.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readTimeControls.H"
+
+        #include "CourantNo.H"
+        #include "alphaCourantNo.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
+        {
+            #include "alphaControls.H"
+            #include "alphaEqnSubCycle.H"
+
+            mixture.correct();
+
+            if (frozenFlow)
+            {
+                continue;
+            }
+
+            #include "UEqn.H"
+
+            // --- Pressure corrector loop
+            while (pimple.correct())
+            {
+                #include "pEqn.H"
+            }
+
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/interIsoFoam/pEqn.H

@@ -0,0 +1,67 @@
+{
+    volScalarField rAU("rAU", 1.0/UEqn.A());
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
+
+    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
+
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        fvc::flux(HbyA)
+      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
+    );
+    MRF.makeRelative(phiHbyA);
+    adjustPhi(phiHbyA, U, p_rgh);
+
+    surfaceScalarField phig
+    (
+        (
+            mixture.surfaceTensionForce()
+          - ghf*fvc::snGrad(rho)
+        )*rAUf*mesh.magSf()
+//          - ghf*(fvc::grad(rho) & mesh.Sf())*rAUf
+    );
+
+    phiHbyA += phig;
+
+    // Update the pressure BCs to ensure flux consistency
+    constrainPressure(p_rgh, U, phiHbyA, rAUf, MRF);
+
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix p_rghEqn
+        (
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
+        );
+
+        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
+
+        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi = phiHbyA - p_rghEqn.flux();
+
+            p_rgh.relax();
+
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
+            U.correctBoundaryConditions();
+            fvOptions.correct(U);
+        }
+    }
+
+    #include "continuityErrs.H"
+
+    p == p_rgh + rho*gh;
+
+    if (p_rgh.needReference())
+    {
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pRefCell)
+        );
+        p_rgh = p - rho*gh;
+    }
+}

applications/solvers/multiphase/interIsoFoam/setDeltaT.H

@@ -0,0 +1,53 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Global
+    setDeltaT
+
+Description
+    Reset the timestep to maintain a constant maximum courant Number.
+    Reduction of time-step is immediate, but increase is damped to avoid
+    unstable oscillations.
+
+\*---------------------------------------------------------------------------*/
+
+if (adjustTimeStep)
+{
+    scalar maxDeltaTFact =
+        min(maxCo/(CoNum + SMALL), maxAlphaCo/(alphaCoNum + SMALL));
+
+    scalar deltaTFact = min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
+
+    runTime.setDeltaT
+    (
+        min
+        (
+            deltaTFact*runTime.deltaTValue(),
+            maxDeltaT
+        )
+    );
+
+    Info<< "deltaT = " <<  runTime.deltaTValue() << endl;
+}
+
+// ************************************************************************* //

applications/test/parallelOverset/heatTransfer/0.orig/T

@@ -20,8 +20,7 @@ internalField   uniform 273;
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
+    #includeEtc "caseDicts/setConstraintTypes"
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
 
     "(walls|hole)"
     {

applications/test/parallelOverset/heatTransfer/0.orig/zoneID

@@ -21,8 +21,7 @@ internalField   uniform 0;
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
+    #includeEtc "caseDicts/setConstraintTypes"
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
 
     overset
     {

applications/utilities/mesh/generation/snappyHexMesh/meshQualityDict

@@ -15,7 +15,7 @@ FoamFile
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 // Include defaults parameters from master dictionary
-#include "$WM_PROJECT_DIR/etc/caseDicts/meshQualityDict"
+#includeEtc "caseDicts/meshQualityDict"
 
 
 // ************************************************************************* //

applications/utilities/preProcessing/setAlphaField/Make/files

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

applications/utilities/preProcessing/setAlphaField/Make/options

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

applications/utilities/preProcessing/setAlphaField/setAlphaField.C

@@ -0,0 +1,191 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+                isoAdvector | Copyright (C) 2016-2017 DHI
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    setAlphaField
+
+Description
+    Uses isoCutCell to create a volume fraction field from either a cylinder,
+    a sphere or a plane.
+
+    Original code supplied by Johan Roenby, DHI (2016)
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "isoCutFace.H"
+#include "isoCutCell.H"
+#include "Enum.H"
+#include "mathematicalConstants.H"
+
+using namespace Foam::constant;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+class shapeSelector
+{
+    public:
+
+        enum class shapeType
+        {
+            PLANE,
+            SPHERE,
+            CYLINDER,
+            SIN
+        };
+
+    static const Foam::Enum<shapeType> shapeTypeNames;
+};
+
+const Foam::Enum<shapeSelector::shapeType> shapeSelector::shapeTypeNames
+{
+    { shapeSelector::shapeType::PLANE, "plane" },
+    { shapeSelector::shapeType::SPHERE, "sphere" },
+    { shapeSelector::shapeType::CYLINDER, "cylinder" },
+    { shapeSelector::shapeType::SIN, "sin" },
+};
+
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+
+    Info<< "Reading setAlphaFieldDict\n" << endl;
+
+    IOdictionary dict
+    (
+        IOobject
+        (
+            "setAlphaFieldDict",
+            runTime.system(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    );
+
+    const shapeSelector::shapeType surfType
+    (
+        shapeSelector::shapeTypeNames.read(dict.lookup("type"))
+    );
+    const vector centre(dict.lookup("centre"));
+    const word fieldName(dict.lookup("field"));
+
+    Info<< "Reading field " << fieldName << "\n" << endl;
+    volScalarField alpha1
+    (
+        IOobject
+        (
+            fieldName,
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    scalar f0 = 0.0;
+    scalarField f(mesh.points().size());
+
+    Info<< "Processing type '" << shapeSelector::shapeTypeNames[surfType]
+        << "'" << endl;
+
+    switch (surfType)
+    {
+        case shapeSelector::shapeType::PLANE:
+        {
+            const vector direction(dict.lookup("direction"));
+
+            f = -(mesh.points() - centre) & (direction/mag(direction));
+            f0 = 0.0;
+            break;
+        }
+        case shapeSelector::shapeType::SPHERE:
+        {
+            const scalar radius(readScalar(dict.lookup("radius")));
+
+            f = -mag(mesh.points() - centre);
+            f0 = -radius;
+            break;
+        }
+        case shapeSelector::shapeType::CYLINDER:
+        {
+            const scalar radius(readScalar(dict.lookup("radius")));
+            const vector direction(dict.lookup("direction"));
+
+            f = -sqrt
+            (
+                sqr(mag(mesh.points() - centre))
+              - sqr(mag((mesh.points() - centre) & direction))
+            );
+            f0 = -radius;
+            break;
+        }
+        case shapeSelector::shapeType::SIN:
+        {
+            const scalar period(readScalar(dict.lookup("period")));
+            const scalar amplitude(readScalar(dict.lookup("amplitude")));
+            const vector up(dict.lookup("up"));
+            const vector direction(dict.lookup("direction"));
+
+            const scalarField xx
+            (
+                (mesh.points() - centre) & direction/mag(direction)
+            );
+            const scalarField zz((mesh.points() - centre) & up/mag(up));
+
+            f = amplitude*Foam::sin(2*mathematical::pi*xx/period) - zz;
+            f0 = 0;
+            break;
+        }
+    }
+
+
+    // Define function on mesh points and isovalue
+
+    // Calculating alpha1 volScalarField from f = f0 isosurface
+    isoCutCell icc(mesh, f);
+    icc.volumeOfFluid(alpha1, f0);
+
+    // Writing volScalarField alpha1
+    ISstream::defaultPrecision(18);
+    alpha1.write();
+
+    Info<< nl << "Phase-1 volume fraction = "
+        << alpha1.weightedAverage(mesh.Vsc()).value()
+        << "  Min(" << alpha1.name() << ") = " << min(alpha1).value()
+        << "  Max(" << alpha1.name() << ") - 1 = " << max(alpha1).value() - 1
+        << nl << endl;
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/utilities/preProcessing/setAlphaField/setAlphaFieldDict

@@ -0,0 +1,24 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+field           alpha.water;
+type            cylinder;
+radius          0.25;
+direction       (0 1 0);
+centre          (0.5 0 0.5);
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -234,25 +234,16 @@ void Foam::epsilonWallFunctionFvPatchScalarField::calculate
     {
         const label celli = patch.faceCells()[facei];
 
-        const scalar yPlus = Cmu25*y[facei]*sqrt(k[celli])/nuw[facei];
-
         const scalar w = cornerWeights[facei];
 
-        if (yPlus > yPlusLam_)
+        epsilon0[celli] += w*Cmu75*pow(k[celli], 1.5)/(kappa_*y[facei]);
-        {
-            epsilon0[celli] += w*Cmu75*pow(k[celli], 1.5)/(kappa_*y[facei]);
 
-            G0[celli] +=
+        G0[celli] +=
-                w
+            w
-               *(nutw[facei] + nuw[facei])
+            *(nutw[facei] + nuw[facei])
-               *magGradUw[facei]
+            *magGradUw[facei]
-               *Cmu25*sqrt(k[celli])
+            *Cmu25*sqrt(k[celli])
-               /(kappa_*y[facei]);
+            /(kappa_*y[facei]);
-        }
-        else
-        {
-            epsilon0[celli] += w*2.0*k[celli]*nuw[facei]/sqr(y[facei]);
-        }
     }
 }
 

src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -238,8 +238,6 @@ void omegaWallFunctionFvPatchScalarField::calculate
     {
         const label celli = patch.faceCells()[facei];
 
-        const scalar yPlus = Cmu25*y[facei]*sqrt(k[celli])/nuw[facei];
-
         const scalar w = cornerWeights[facei];
 
         const scalar omegaVis = 6*nuw[facei]/(beta1_*sqr(y[facei]));
@@ -257,27 +255,17 @@ void omegaWallFunctionFvPatchScalarField::calculate
             omega0[celli] += w*sqrt(sqr(omegaVis) + sqr(omegaLog));
         }
 
-        if (yPlus > yPlusLam_)
+        if (!blended_)
         {
-            if (!blended_)
+            omega0[celli] += w*omegaLog;
-            {
+        }
-                omega0[celli] += w*omegaLog;
-            }
 
-            G0[celli] +=
+        G0[celli] +=
-                w
+            w
-               *(nutw[facei] + nuw[facei])
+            *(nutw[facei] + nuw[facei])
-               *magGradUw[facei]
+            *magGradUw[facei]
-               *Cmu25*sqrt(k[celli])
+            *Cmu25*sqrt(k[celli])
-               /(kappa_*y[facei]);
+            /(kappa_*y[facei]);
-        }
-        else
-        {
-            if (!blended_)
-            {
-                omega0[celli] += w*omegaVis;
-            }
-        }
     }
 }
 

src/finiteVolume/Make/files

@@ -238,6 +238,9 @@ fvMatrices/fvMatrices.C
 fvMatrices/fvScalarMatrix/fvScalarMatrix.C
 fvMatrices/solvers/MULES/MULES.C
 fvMatrices/solvers/MULES/CMULES.C
+fvMatrices/solvers/isoAdvection/isoCutCell/isoCutCell.C
+fvMatrices/solvers/isoAdvection/isoCutFace/isoCutFace.C
+fvMatrices/solvers/isoAdvection/isoAdvection/isoAdvection.C
 fvMatrices/solvers/GAMGSymSolver/GAMGAgglomerations/faceAreaPairGAMGAgglomeration/faceAreaPairGAMGAgglomeration.C
 
 interpolation = interpolation/interpolation

src/finiteVolume/fvMatrices/solvers/isoAdvection/isoAdvection/isoAdvection.H

@@ -0,0 +1,388 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016-2017 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+                isoAdvector | Copyright (C) 2016-2017 DHI
+-------------------------------------------------------------------------------
+
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::isoAdvection
+
+Description
+    Calculates the new VOF (alpha) field after time step dt given the initial
+    VOF field and a velocity field U and face fluxes phi. The fluid transport
+    calculation is based on an idea of using isosurfaces to estimate the
+    internal distribution of fluid in cells and advecting such isosurfaces
+    across the mesh faces with the velocity field interpolated to the
+    isosurfaces.
+
+    Reference:
+        \verbatim
+            Roenby, J., Bredmose, H. and Jasak, H. (2016).
+            A computational method for sharp interface advection
+            Royal Society Open Science, 3
+            doi 10.1098/rsos.160405
+        \endverbatim
+
+    Original code supplied by Johan Roenby, DHI (2016)
+
+SourceFiles
+    isoAdvection.C
+    isoAdvectionTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef isoAdvection_H
+#define isoAdvection_H
+
+#include "fvMesh.H"
+#include "volFieldsFwd.H"
+#include "surfaceFields.H"
+#include "className.H"
+#include "isoCutCell.H"
+#include "isoCutFace.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class isoAdvection Declaration
+\*---------------------------------------------------------------------------*/
+
+class isoAdvection
+{
+    // Private data types
+
+        typedef DynamicList<label> DynamicLabelList;
+        typedef DynamicList<scalar> DynamicScalarList;
+        typedef DynamicList<vector> DynamicVectorList;
+        typedef DynamicList<point> DynamicPointList;
+
+
+    // Private data
+
+        //- Reference to mesh
+        const fvMesh& mesh_;
+
+        //- Dictionary for isoAdvection controls
+        const dictionary dict_;
+
+        //- VOF field
+        volScalarField& alpha1_;
+
+        //- Often used reference to alpha1 internal field
+        scalarField& alpha1In_;
+
+        //- Reference to flux field
+        const surfaceScalarField& phi_;
+
+        //- Reference to velocity field
+        const volVectorField& U_;
+
+        //- Face volumetric water transport
+        surfaceScalarField dVf_;
+
+        //- Time spent performing interface advection
+        scalar advectionTime_;
+
+
+        // Point interpolation data
+
+            //- VOF field interpolated to mesh points
+            scalarField ap_;
+
+
+        // Switches and tolerances. Tolerances need to go into toleranceSwitches
+
+            //- Number of alpha bounding steps
+            label nAlphaBounds_;
+
+            //- Tolerance for search of isoFace giving specified VOF value
+            scalar vof2IsoTol_;
+
+            //- Tolerance for marking of surface cells:
+            //  Those with surfCellTol_ < alpha1 < 1 - surfCellTol_
+            scalar surfCellTol_;
+
+            //- Switch controlling whether to use isoface normals for interface
+            //  orientation (default corresponding to false) to base it on
+            //  a smoothed gradient of alpha calculation (giving better results
+            //  on tri on tet meshes).
+            bool gradAlphaBasedNormal_;
+
+            //- Print isofaces in a <case>/isoFaces/isoFaces_#N.vtk files.
+            //  Intended for debugging
+            bool writeIsoFacesToFile_;
+
+        // Cell and face cutting
+
+            //- List of surface cells
+            DynamicLabelList surfCells_;
+
+            //- Cell cutting object
+            isoCutCell isoCutCell_;
+
+            //- Face cutting object
+            isoCutFace isoCutFace_;
+
+            //- Bool list for cells that have been touched by the bounding step
+            boolList cellIsBounded_;
+
+            //- True for all surface cells and their neighbours
+            boolList checkBounding_;
+
+            //- Storage for boundary faces downwind to a surface cell
+            DynamicLabelList bsFaces_;
+
+            //- Storage for boundary surface iso face centre
+            DynamicVectorList bsx0_;
+
+            //- Storage for boundary surface iso face normal
+            DynamicVectorList bsn0_;
+
+            //- Storage for boundary surface iso face speed
+            DynamicScalarList bsUn0_;
+
+            //- Storage for boundary surface iso value
+            DynamicScalarList bsf0_;
+
+
+        // Additional data for parallel runs
+
+            //- List of processor patch labels
+            DynamicLabelList procPatchLabels_;
+
+            //- For each patch if it is a processor patch this is a list of the
+            //  face labels on this patch that are downwind to a surface cell.
+            //  For non-processor patches the list will be empty.
+            List<DynamicLabelList> surfaceCellFacesOnProcPatches_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        isoAdvection(const isoAdvection&);
+
+        //- Disallow default bitwise copy assignment
+        void operator=(const isoAdvection&);
+
+
+        // Advection functions
+
+            //- For each face calculate volumetric face transport during dt
+            void timeIntegratedFlux();
+
+            //- Calculate volumetric face transport during dt given the isoFace
+            //  data provided as input for face facei
+            scalar timeIntegratedFaceFlux
+            (
+                const label facei,
+                const vector& x0,
+                const vector& n0,
+                const scalar Un0,
+                const scalar f0,
+                const scalar dt,
+                const scalar phi,
+                const scalar magSf
+            );
+
+            //- For a given cell return labels of faces fluxing out of this cell
+            //  (based on sign of phi)
+            void setDownwindFaces
+            (
+                const label celli,
+                DynamicLabelList& downwindFaces
+            ) const;
+
+            // Limit fluxes
+            void limitFluxes();
+
+            // Bound fluxes
+            void boundFromAbove
+            (
+                const scalarField& alpha1,
+                surfaceScalarField& dVfcorrected,
+                DynamicLabelList& correctedFaces
+            );
+
+            //- Given the face volume transport dVf calculates the total volume
+            //  leaving a given cell. Note: cannot use dVf member because
+            //  netFlux is called also for corrected dVf
+            scalar netFlux
+            (
+                const surfaceScalarField& dVf,
+                const label celli
+            ) const;
+
+            //- Determine if a cell is a surface cell
+            bool isASurfaceCell(const label celli) const
+            {
+                return
+                (
+                    surfCellTol_ < alpha1In_[celli]
+                 && alpha1In_[celli] < 1 - surfCellTol_
+                );
+            }
+
+            //- Clear out isoFace data
+            void clearIsoFaceData()
+            {
+                surfCells_.clear();
+                bsFaces_.clear();
+                bsx0_.clear();
+                bsn0_.clear();
+                bsUn0_.clear();
+                bsf0_.clear();
+            }
+
+        // Face value functions needed for random face access where the face
+        // can be either internal or boundary face
+
+            //- Return face value for a given Geometric surface field
+            template<typename Type>
+            Type faceValue
+            (
+                const GeometricField<Type, fvsPatchField, surfaceMesh>& f,
+                const label facei
+            ) const;
+
+            //- Set face value for a given Geometric surface field
+            template<typename Type>
+            void setFaceValue
+            (
+                GeometricField<Type, fvsPatchField, surfaceMesh>& f,
+                const label facei,
+                const Type& value
+            ) const;
+
+
+        // Parallel run handling functions
+
+            //- Synchronize dVf across processor boundaries using upwind value
+            void syncProcPatches
+            (
+                surfaceScalarField& dVf,
+                const surfaceScalarField& phi
+            );
+
+            //- Check if the face is on processor patch and append it to the
+            //  list of surface cell faces on processor patches
+            void checkIfOnProcPatch(const label facei);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("isoAdvection");
+
+    //- Constructors
+
+        //- Construct given alpha, phi and velocity field. Note: phi should be
+        // divergence free up to a sufficient tolerance
+        isoAdvection
+        (
+            volScalarField& alpha1,
+            const surfaceScalarField& phi,
+            const volVectorField& U
+        );
+
+
+    //- Destructor
+    virtual ~isoAdvection()
+    {}
+
+
+    // Member functions
+
+        //- Advect the free surface. Updates alpha field, taking into account
+        //  multiple calls within a single time step.
+        void advect();
+
+        //- Apply the bounding based on user inputs
+        void applyBruteForceBounding();
+
+        // Access functions
+
+            //- Return alpha field
+            const volScalarField& alpha() const
+            {
+                return alpha1_;
+            }
+
+            //- Return the controls dictionary
+            const dictionary& dict() const
+            {
+                return dict_;
+            }
+
+            //- Return cellSet of surface cells
+            void writeSurfaceCells() const;
+
+            //- Return cellSet of bounded cells
+            void writeBoundedCells() const;
+
+            //- Return mass flux
+            tmp<surfaceScalarField> getRhoPhi
+            (
+                const dimensionedScalar rho1,
+                const dimensionedScalar rho2
+            ) const
+            {
+                return tmp<surfaceScalarField>
+                (
+                    new surfaceScalarField
+                    (
+                        "rhoPhi",
+                        (rho1 - rho2)*dVf_/mesh_.time().deltaT() + rho2*phi_
+                    )
+                );
+            }
+
+            scalar advectionTime() const
+            {
+                return advectionTime_;
+            }
+
+            //- Write isoface points to .obj file
+            void writeIsoFaces
+            (
+                const DynamicList<List<point>>& isoFacePts
+            ) const;
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "isoAdvectionTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/fvMatrices/solvers/isoAdvection/isoAdvection/isoAdvectionTemplates.C

@@ -0,0 +1,111 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016-2017 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+                isoAdvector | Copyright (C) 2016-2017 DHI
+-------------------------------------------------------------------------------
+
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "isoAdvection.H"
+
+// ************************************************************************* //
+
+template<typename Type>
+Type Foam::isoAdvection::faceValue
+(
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& f,
+    const label facei
+) const
+{
+    if (mesh_.isInternalFace(facei))
+    {
+        return f.primitiveField()[facei];
+    }
+    else
+    {
+        const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
+
+        // Boundary face. Find out which face of which patch
+        const label patchi = pbm.patchID()[facei - mesh_.nInternalFaces()];
+
+        if (patchi < 0 || patchi >= pbm.size())
+        {
+            FatalErrorInFunction
+                << "Cannot find patch for face " << facei
+                << abort(FatalError);
+        }
+
+        // Handle empty patches
+        const polyPatch& pp = pbm[patchi];
+        if (isA<emptyPolyPatch>(pp) || pp.empty())
+        {
+            return pTraits<Type>::zero;
+        }
+
+        const label patchFacei = pp.whichFace(facei);
+        return f.boundaryField()[patchi][patchFacei];
+    }
+}
+
+
+template<typename Type>
+void Foam::isoAdvection::setFaceValue
+(
+    GeometricField<Type, fvsPatchField, surfaceMesh>& f,
+    const label facei,
+    const Type& value
+) const
+{
+    if (mesh_.isInternalFace(facei))
+    {
+        f.primitiveFieldRef()[facei] = value;
+    }
+    else
+    {
+        const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
+
+        // Boundary face. Find out which face of which patch
+        const label patchi = pbm.patchID()[facei - mesh_.nInternalFaces()];
+
+        if (patchi < 0 || patchi >= pbm.size())
+        {
+            FatalErrorInFunction
+                << "Cannot find patch for face " << facei
+                << abort(FatalError);
+        }
+
+        // Handle empty patches
+        const polyPatch& pp = pbm[patchi];
+        if (isA<emptyPolyPatch>(pp) || pp.empty())
+        {
+            return;
+        }
+
+        const label patchFacei = pp.whichFace(facei);
+
+        f.boundaryFieldRef()[patchi][patchFacei] = value;
+    }
+}
+
+
+// ************************************************************************* //

src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutCell/isoCutCell.C

@@ -0,0 +1,709 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+                isoAdvector | Copyright (C) 2016 DHI
+-------------------------------------------------------------------------------
+
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "isoCutCell.H"
+#include "scalarMatrices.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+int Foam::isoCutCell::debug = 0;
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::isoCutCell::isoCutCell(const fvMesh& mesh, scalarField& f)
+:
+    mesh_(mesh),
+    cellI_(-1),
+    f_(f),
+    isoValue_(0),
+    isoCutFace_(isoCutFace(mesh_, f_)),
+    isoCutFaces_(10),
+    isoCutFacePoints_(10),
+    isoCutFaceCentres_(10),
+    isoCutFaceAreas_(10),
+    isoFaceEdges_(10),
+    isoFacePoints_(10),
+    isoFaceCentre_(vector::zero),
+    isoFaceArea_(vector::zero),
+    subCellCentre_(vector::zero),
+    subCellVolume_(-10),
+    VOF_(-10),
+    fullySubFaces_(10),
+    cellStatus_(-1),
+    subCellCentreAndVolumeCalculated_(false),
+    isoFaceCentreAndAreaCalculated_(false)
+{
+    clearStorage();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::isoCutCell::calcSubCellCentreAndVolume()
+{
+    if (cellStatus_ == 0)
+    {
+        subCellCentre_ = vector::zero;
+        subCellVolume_ = 0.0;
+
+        // Estimate the approximate cell centre as the average of face centres
+        label nCellFaces(1 + isoCutFaceCentres_.size() + fullySubFaces_.size());
+        vector cEst = isoFaceCentre_ + sum(isoCutFaceCentres_);
+        forAll(fullySubFaces_, facei)
+        {
+            cEst += mesh_.faceCentres()[fullySubFaces_[facei]];
+        }
+        cEst /= scalar(nCellFaces);
+
+
+        // Contribution to subcell centre and volume from isoface
+        const scalar pyr3Vol0 =
+            max(mag(isoFaceArea_ & (isoFaceCentre_ - cEst)), VSMALL);
+
+        // Calculate face-pyramid centre
+        const vector pc0 = 0.75*isoFaceCentre_ + 0.25*cEst;
+
+        // Accumulate volume-weighted face-pyramid centre
+        subCellCentre_ += pyr3Vol0*pc0;
+
+        // Accumulate face-pyramid volume
+        subCellVolume_ += pyr3Vol0;
+
+        // Contribution to subcell centre and volume from cut faces
+        forAll(isoCutFaceCentres_, facei)
+        {
+            // Calculate 3*face-pyramid volume
+            scalar pyr3Vol =
+                max
+                (
+                    mag
+                    (
+                        isoCutFaceAreas_[facei]
+                      & (isoCutFaceCentres_[facei] - cEst)
+                    ),
+                    VSMALL
+                );
+
+            // Calculate face-pyramid centre
+            vector pc = 0.75*isoCutFaceCentres_[facei] + 0.25*cEst;
+
+            // Accumulate volume-weighted face-pyramid centre
+            subCellCentre_ += pyr3Vol*pc;
+
+            // Accumulate face-pyramid volume
+            subCellVolume_ += pyr3Vol;
+        }
+
+        // Contribution to subcell centre and volume from fully submerged faces
+        forAll(fullySubFaces_, i)
+        {
+            const label facei = fullySubFaces_[i];
+            const point& fCentre = mesh_.faceCentres()[facei];
+            const vector& fArea = mesh_.faceAreas()[facei];
+
+            // Calculate 3*face-pyramid volume
+            scalar pyr3Vol = max(mag(fArea & (fCentre - cEst)), VSMALL);
+
+            // Calculate face-pyramid centre
+            vector pc = 0.75*fCentre + 0.25*cEst;
+
+            // Accumulate volume-weighted face-pyramid centre
+            subCellCentre_ += pyr3Vol*pc;
+
+            // Accumulate face-pyramid volume
+            subCellVolume_ += pyr3Vol;
+        }
+
+        subCellCentre_ /= subCellVolume_;
+        subCellVolume_ /= scalar(3);
+        VOF_ = subCellVolume_/mesh_.cellVolumes()[cellI_];
+
+        subCellCentreAndVolumeCalculated_ = true;
+
+        if (debug)
+        {
+            vector sumSf = isoFaceArea_;
+            scalar sumMagSf = mag(isoFaceArea_);
+            forAll(isoCutFaceCentres_, facei)
+            {
+                sumSf += isoCutFaceAreas_[facei];
+                sumMagSf += mag(isoCutFaceAreas_[facei]);
+            }
+            forAll(fullySubFaces_, facei)
+            {
+                sumSf += mesh_.faceAreas()[fullySubFaces_[facei]];
+                sumMagSf += mag(isoCutFaceAreas_[facei]);
+            }
+            if (mag(sumSf) > 1e-10)
+            {
+                Pout<< "Warninig: mag(sumSf)/magSumSf = "
+                    << mag(sumSf)/sumMagSf << " for surface cell"
+                    << cellI_ << endl;
+            }
+        }
+    }
+    else if (cellStatus_ == 1)
+    {
+        // Cell fully above isosurface
+        subCellCentre_ = vector::zero;
+        subCellVolume_ = 0;
+        VOF_ = 0;
+    }
+    else if (cellStatus_ == -1)
+    {
+        // Cell fully below isosurface
+        subCellCentre_ = mesh_.cellCentres()[cellI_];
+        subCellVolume_ = mesh_.cellVolumes()[cellI_];
+        VOF_ = 1;
+    }
+}
+
+
+void Foam::isoCutCell::calcIsoFaceCentreAndArea()
+{
+    // Initial guess of face centre from edge points
+    point fCentre = vector::zero;
+    label nEdgePoints = 0;
+    forAll(isoFaceEdges_, ei)
+    {
+        DynamicList<point>& edgePoints = isoFaceEdges_[ei];
+        forAll(edgePoints, pi)
+        {
+            fCentre += edgePoints[pi];
+            nEdgePoints++;
+        }
+    }
+
+    if (nEdgePoints > 0)
+    {
+        fCentre /= nEdgePoints;
+    }
+    else
+    {
+        DebugPout << "Warning: nEdgePoints = 0 for cell " << cellI_ << endl;
+    }
+
+    vector sumN = vector::zero;
+    scalar sumA = 0;
+    vector sumAc = vector::zero;
+
+    forAll(isoFaceEdges_, ei)
+    {
+        const DynamicList<point>& edgePoints = isoFaceEdges_[ei];
+        const label nPoints = edgePoints.size();
+        for (label pi = 0; pi < nPoints-1; pi++)
+        {
+            const point& nextPoint = edgePoints[pi + 1];
+
+            vector c = edgePoints[pi] + nextPoint + fCentre;
+            vector n = (nextPoint - edgePoints[pi])^(fCentre - edgePoints[pi]);
+            scalar a = mag(n);
+
+            // Edges may have different orientation
+            sumN += Foam::sign(n & sumN)*n;
+            sumA += a;
+            sumAc += a*c;
+        }
+    }
+
+    // This is to deal with zero-area faces. Mark very small faces
+    // to be detected in e.g., processorPolyPatch.
+    if (sumA < ROOTVSMALL)
+    {
+        isoFaceCentre_ = fCentre;
+        isoFaceArea_ = vector::zero;
+    }
+    else
+    {
+        isoFaceCentre_ = sumAc/sumA/scalar(3);
+        isoFaceArea_ = 0.5*sumN;
+    }
+
+
+    // Check isoFaceArea_ direction and change if not pointing out of subcell
+    if ((isoFaceArea_ & (isoFaceCentre_ - subCellCentre())) < 0)
+    {
+        isoFaceArea_ *= (-1);
+    }
+
+    isoFaceCentreAndAreaCalculated_ = true;
+}
+
+
+void Foam::isoCutCell::calcIsoFacePointsFromEdges()
+{
+    DebugPout
+        << "Enter calcIsoFacePointsFromEdges() with isoFaceArea_ = "
+        << isoFaceArea_ << " and isoFaceCentre_ = " << isoFaceCentre_
+        << " and isoFaceEdges_ = " << isoFaceEdges_ << endl;
+
+    // Defining local coordinates with zhat along isoface normal and xhat from
+    // isoface centre to first point in isoFaceEdges_
+    const vector zhat = isoFaceArea_/mag(isoFaceArea_);
+    vector xhat = isoFaceEdges_[0][0] - isoFaceCentre_;
+    xhat = (xhat - (xhat & zhat)*zhat);
+    xhat /= mag(xhat);
+    vector yhat = zhat^xhat;
+    yhat /= mag(yhat);
+
+    DebugPout << "Calculated local coordinates" << endl;
+
+    // Calculating isoface point angles in local coordinates
+    DynamicList<point> unsortedIsoFacePoints(3*isoFaceEdges_.size());
+    DynamicList<scalar> unsortedIsoFacePointAngles(3*isoFaceEdges_.size());
+    forAll(isoFaceEdges_, ei)
+    {
+        const DynamicList<point>& edgePoints = isoFaceEdges_[ei];
+        forAll(edgePoints, pi)
+        {
+            const point& p = edgePoints[pi];
+            unsortedIsoFacePoints.append(p);
+            unsortedIsoFacePointAngles.append
+            (
+                Foam::atan2
+                (
+                    ((p - isoFaceCentre_) & yhat),
+                    ((p - isoFaceCentre_) & xhat)
+                )
+            );
+        }
+    }
+
+    DebugPout<< "Calculated isoFace point angles" << endl;
+
+    // Sorting isoface points by angle and inserting into isoFacePoints_
+    labelList order(unsortedIsoFacePointAngles.size());
+    Foam::sortedOrder(unsortedIsoFacePointAngles, order);
+    isoFacePoints_.append(unsortedIsoFacePoints[order[0]]);
+    for (label pi = 1; pi < order.size(); pi++)
+    {
+        if
+        (
+            mag
+            (
+                unsortedIsoFacePointAngles[order[pi]]
+               -unsortedIsoFacePointAngles[order[pi-1]]
+            ) > 1e-8
+        )
+        {
+            isoFacePoints_.append(unsortedIsoFacePoints[order[pi]]);
+        }
+    }
+
+    DebugPout<< "Sorted isoface points by angle" << endl;
+}
+
+
+Foam::label Foam::isoCutCell::calcSubCell
+(
+    const label celli,
+    const scalar isoValue
+)
+{
+    // Populate isoCutFaces_, isoCutFacePoints_, fullySubFaces_, isoFaceCentres_
+    // and isoFaceArea_.
+
+    clearStorage();
+    cellI_ = celli;
+    isoValue_ = isoValue;
+    const cell& c = mesh_.cells()[celli];
+
+    forAll(c, fi)
+    {
+        const label facei = c[fi];
+
+        const label faceStatus = isoCutFace_.calcSubFace(facei, isoValue_);
+
+        if (faceStatus == 0)
+        {
+            // Face is cut
+            isoCutFacePoints_.append(isoCutFace_.subFacePoints());
+            isoCutFaceCentres_.append(isoCutFace_.subFaceCentre());
+            isoCutFaceAreas_.append(isoCutFace_.subFaceArea());
+            isoFaceEdges_.append(isoCutFace_.surfacePoints());
+        }
+        else if (faceStatus == -1)
+        {
+            // Face fully below
+            fullySubFaces_.append(facei);
+        }
+    }
+
+    if (isoCutFacePoints_.size())
+    {
+        // Cell cut at least at one face
+        cellStatus_ = 0;
+        calcIsoFaceCentreAndArea();
+    }
+    else if (fullySubFaces_.empty())
+    {
+        // Cell fully above isosurface
+        cellStatus_ = 1;
+    }
+    else
+    {
+        // Cell fully below isosurface
+        cellStatus_ = -1;
+    }
+
+    return cellStatus_;
+}
+
+
+const Foam::point& Foam::isoCutCell::subCellCentre()
+{
+    if (!subCellCentreAndVolumeCalculated_)
+    {
+        calcSubCellCentreAndVolume();
+    }
+
+    return subCellCentre_;
+}
+
+
+Foam::scalar Foam::isoCutCell::subCellVolume()
+{
+    if (!subCellCentreAndVolumeCalculated_)
+    {
+        calcSubCellCentreAndVolume();
+    }
+
+    return subCellVolume_;
+}
+
+
+const Foam::DynamicList<Foam::point>& Foam::isoCutCell::isoFacePoints()
+{
+    if (cellStatus_ == 0 && isoFacePoints_.size() == 0)
+    {
+        calcIsoFacePointsFromEdges();
+    }
+
+    return isoFacePoints_;
+}
+
+
+const Foam::point& Foam::isoCutCell::isoFaceCentre()
+{
+    if (!isoFaceCentreAndAreaCalculated_)
+    {
+        calcIsoFaceCentreAndArea();
+    }
+
+    return isoFaceCentre_;
+}
+
+
+const Foam::vector& Foam::isoCutCell::isoFaceArea()
+{
+    if (!isoFaceCentreAndAreaCalculated_)
+    {
+        calcIsoFaceCentreAndArea();
+    }
+
+    return isoFaceArea_;
+}
+
+
+Foam::scalar Foam::isoCutCell::volumeOfFluid()
+{
+    if (!subCellCentreAndVolumeCalculated_)
+    {
+        calcSubCellCentreAndVolume();
+    }
+
+    return VOF_;
+}
+
+
+Foam::scalar Foam::isoCutCell::isoValue() const
+{
+    return isoValue_;
+}
+
+
+void Foam::isoCutCell::clearStorage()
+{
+    cellI_ = -1;
+    isoValue_ = 0;
+    isoCutFace_.clearStorage();
+    isoCutFaces_.clear();
+    isoCutFacePoints_.clear();
+    isoCutFaceCentres_.clear();
+    isoCutFaceAreas_.clear();
+    isoFaceEdges_.clear();
+    isoFacePoints_.clear();
+    isoFaceCentre_ = vector::zero;
+    isoFaceArea_ = vector::zero;
+    subCellCentre_ = vector::zero;
+    subCellVolume_ = -10;
+    VOF_ = -10;
+    fullySubFaces_.clear();
+    cellStatus_ = -1;
+    subCellCentreAndVolumeCalculated_ = false;
+    isoFaceCentreAndAreaCalculated_ = false;
+}
+
+
+Foam::label Foam::isoCutCell::vofCutCell
+(
+    const label celli,
+    const scalar alpha1,
+    const scalar tol,
+    const label maxIter
+)
+{
+    DebugInFunction
+        << "vofCutCell for cell " << celli << " with alpha1 = "
+        << alpha1 << " ------" << endl;
+
+    // Finding cell vertex extrema values
+    const labelList& pLabels = mesh_.cellPoints(celli);
+    scalarField fvert(pLabels.size());
+    forAll(pLabels, pi)
+    {
+        fvert[pi] = f_[pLabels[pi]];
+    }
+    labelList order(fvert.size());
+    sortedOrder(fvert, order);
+    scalar f1 = fvert[order.first()];
+    scalar f2 = fvert[order.last()];
+
+    DebugPout << "fvert = " << fvert << ", and order = " << order << endl;
+
+    // Handling special case where method is handed an almost full/empty cell
+    if (alpha1 < tol)
+    {
+        return calcSubCell(celli, f2);
+    }
+    else if (1 - alpha1 < tol)
+    {
+        return calcSubCell(celli, f1);
+    }
+
+    // Finding the two vertices inbetween which the isovalue giving alpha1 lies
+    label L1 = 0;
+    label L2 = fvert.size() - 1;
+    scalar a1 = 1;
+    scalar a2 = 0;
+    scalar L3, f3, a3;
+
+    while (L2 - L1 > 1)
+    {
+        L3 = round(0.5*(L1 + L2));
+        f3 = fvert[order[L3]];
+        calcSubCell(celli, f3);
+        a3 = volumeOfFluid();
+        if (a3 > alpha1)
+        {
+            L1 = L3; f1 = f3; a1 = a3;
+        }
+        else if (a3 < alpha1)
+        {
+            L2 = L3; f2 = f3; a2 = a3;
+        }
+    }
+
+    if (mag(f1 - f2) < 10*SMALL)
+    {
+        DebugPout<< "Warning: mag(f1 - f2) < 10*SMALL" << endl;
+        return calcSubCell(celli, f1);
+    }
+
+    if (mag(a1 - a2) < tol)
+    {
+        DebugPout<< "Warning: mag(a1 - a2) < tol for cell " << celli << endl;
+        return calcSubCell(celli, 0.5*(f1 + f2));
+    }
+
+    // Now we know that a(f) = alpha1 is to be found on the f interval
+    // [f1, f2], i.e. alpha1 will be in the interval [a2,a1]
+    DebugPout
+        << "L1 = " << L1 << ", f1 = " << f1 << ", a1 = " << a1 << nl
+        << "L2 = " << L2 << ", f2 = " << f2  << ", a2 = " << a2 << endl;
+
+
+    // Finding coefficients in 3 deg polynomial alpha(f) from 4 solutions
+
+    // Finding 2 additional points on 3 deg polynomial
+    f3 = f1 + (f2 - f1)/scalar(3);
+    calcSubCell(celli, f3);
+    a3 = volumeOfFluid();
+
+    scalar f4 = f1 + 2*(f2 - f1)/3;
+    calcSubCell(celli, f4);
+    scalar a4 = volumeOfFluid();
+
+    // Building and solving Vandermonde matrix equation
+    scalarField a(4), f(4), C(4);
+    {
+        a[0] = a1, a[1] = a3, a[2] = a4, a[3] = a2;
+        f[0] = 0, f[1] = (f3-f1)/(f2-f1), f[2] = (f4-f1)/(f2-f1), f[3] = 1;
+        scalarSquareMatrix M(4);
+        forAll(f, i)
+        {
+            forAll(f, j)
+            {
+                M[i][j] = pow(f[i], 3 - j);
+            }
+        }
+
+        // C holds the 4 polynomial coefficients
+        C = a;
+        LUsolve(M, C);
+    }
+
+    // Finding root with Newton method
+
+    f3 = f[1]; a3 = a[1];
+    label nIter = 0;
+    scalar res = mag(a3 - alpha1);
+    while (res > tol && nIter < 10*maxIter)
+    {
+        f3 -=
+            (C[0]*pow3(f3) + C[1]*sqr(f3) + C[2]*f3 + C[3] - alpha1)
+           /(3*C[0]*sqr(f3) + 2*C[1]*f3 + C[2]);
+        a3 = C[0]*pow3(f3) + C[1]*sqr(f3) + C[2]*f3 + C[3];
+        res = mag(a3 - alpha1);
+        nIter++;
+    }
+    // Scaling back to original range
+    f3 = f3*(f2 - f1) + f1;
+
+    // Check result
+    calcSubCell(celli, f3);
+    const scalar VOF = volumeOfFluid();
+    res = mag(VOF - alpha1);
+
+    if (res > tol)
+    {
+        DebugPout
+            << "Newton obtained f3 = " << f3 << " and a3 = " << a3
+            << " with mag(a3-alpha1) = " << mag(a3-alpha1)
+            << " but calcSubCell(celli,f3) gives VOF  = " << VOF << nl
+            << "M(f)*C = a with " << nl
+            << "f_scaled = " << f << nl
+            << "f = " << f*(f2 - f1) + f1 << nl
+            << "a = " << a << nl
+            << "C = " << C << endl;
+    }
+    else
+    {
+        DebugPout<< "Newton did the job" << endl;
+        return cellStatus_;
+    }
+
+    // If tolerance not met use the secant method  with f3 as a hopefully very
+    // good initial guess to crank res the last piece down below tol
+
+    scalar x2 = f3;
+    scalar g2 = VOF - alpha1;
+    scalar x1 = max(1e-3*(f2 - f1), 100*SMALL);
+    x1 = min(max(x1, f1), f2);
+    calcSubCell(celli, x1);
+    scalar g1 = volumeOfFluid() - alpha1;
+
+    nIter = 0;
+    scalar g0(0), x0(0);
+    while (res > tol && nIter < maxIter && g1 != g2)
+    {
+        x0 = (x2*g1 - x1*g2)/(g1 - g2);
+        calcSubCell(celli, x0);
+        g0 = volumeOfFluid() - alpha1;
+        res = mag(g0);
+        x2 = x1; g2 = g1;
+        x1 = x0; g1 = g0;
+        nIter++;
+    }
+
+    if (debug)
+    {
+        if (res < tol)
+        {
+            Pout<< "Bisection finished the job in " << nIter << " iterations."
+                << endl;
+        }
+        else
+        {
+            Pout<< "Warning: Bisection not converged " << endl;
+            Pout<< "Leaving vofCutCell with f3 = " << f3 << " giving a3 = "
+                << a3 << " so alpha1 - a3 = " << alpha1 - a3 << endl;
+        }
+    }
+
+    return cellStatus_;
+}
+
+
+void Foam::isoCutCell::volumeOfFluid
+(
+    volScalarField& alpha1,
+    const scalar f0
+)
+{
+    // Setting internal field
+    scalarField& alphaIn = alpha1;
+    forAll(alphaIn, celli)
+    {
+        const label cellStatus = calcSubCell(celli, f0);
+        if (cellStatus != 1)
+        {
+            // If cell not entirely above isosurface
+            alphaIn[celli] = volumeOfFluid();
+        }
+    }
+
+    // Setting boundary alpha1 values
+    forAll(mesh_.boundary(), patchi)
+    {
+        if (mesh_.boundary()[patchi].size() > 0)
+        {
+            const label start = mesh_.boundary()[patchi].patch().start();
+            scalarField& alphap = alpha1.boundaryFieldRef()[patchi];
+            const scalarField& magSfp = mesh_.magSf().boundaryField()[patchi];
+
+            forAll(alphap, patchFacei)
+            {
+                const label facei = patchFacei + start;
+                const label faceStatus = isoCutFace_.calcSubFace(facei, f0);
+
+                if (faceStatus != 1)
+                {
+                    // Face not entirely above isosurface
+                    alphap[patchFacei] =
+                        mag(isoCutFace_.subFaceArea())/magSfp[patchFacei];
+                }
+            }
+        }
+    }
+}
+
+
+// ************************************************************************* //

src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutCell/isoCutCell.H

@@ -0,0 +1,200 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+                isoAdvector | Copyright (C) 2016 DHI
+-------------------------------------------------------------------------------
+
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::isoCutCell
+
+Description
+    Class for cutting a cell, celli, of an fvMesh, mesh_, at its intersection
+    with an isosurface defined by the mesh point values f_ and the isovalue,
+    isoValue_.
+
+    Reference:
+        \verbatim
+            Roenby, J., Bredmose, H. and Jasak, H. (2016).
+            A computational method for sharp interface advection
+            Royal Society Open Science, 3
+            doi 10.1098/rsos.160405
+        \endverbatim
+
+    Original code supplied by Johan Roenby, DHI (2016)
+
+SourceFiles
+    isoCutCell.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef isoCutCell_H
+#define isoCutCell_H
+
+#include "fvMesh.H"
+#include "volFieldsFwd.H"
+#include "isoCutFace.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class isoCutCell Declaration
+\*---------------------------------------------------------------------------*/
+
+class isoCutCell
+{
+    // Private data
+
+        //- Mesh whose cells and faces to cut at their intersection with an
+        //  isosurface.
+        const fvMesh& mesh_;
+
+        //- Cell to cut
+        label cellI_;
+
+        //- Isofunction values at mesh points. f_size() = mesh_.nPoints().
+        scalarField& f_;
+
+        //- Isovalue used to cut cell
+        scalar isoValue_;
+
+        //- An isoCutFace object to get access to its face cutting functionality
+        isoCutFace isoCutFace_;
+
+        //- List of face labels of isoCutFaces
+        DynamicList<label> isoCutFaces_;
+
+        //- List of point lists each defining an isoCutFace
+        DynamicList<DynamicList<point>> isoCutFacePoints_;
+
+        //- List of face centres for isoCutFaces
+        DynamicList<point> isoCutFaceCentres_;
+
+        //- List of face area vectors for isoCutFaces
+        DynamicList<vector> isoCutFaceAreas_;
+
+        //- Storage for subFace edges belonging to isoFace
+        DynamicList<DynamicList<point>> isoFaceEdges_;
+
+        //- Points constituting the cell-isosurface intersection (isoface)
+        DynamicList<point> isoFacePoints_;
+
+        //- Face centre of the isoface
+        point isoFaceCentre_;
+
+        //- Face normal of the isoface by convention pointing from high to low
+        //  values (i.e. opposite of the gradient vector).
+        vector isoFaceArea_;
+
+        //- Cell centre of the subcell of celli which is "fully submerged", i.e.
+        //  where the function value is higher than the isoValue_
+        point subCellCentre_;
+
+        //- Volume of fully submerged subcell
+        scalar subCellVolume_;
+
+        //- Volume of Fluid for celli (subCellVolume_/mesh_.V()[celli])
+        scalar VOF_;
+
+        //- List of fully submerged faces
+        DynamicList<label> fullySubFaces_;
+
+        //- A cell status label taking one of the values:
+        //
+        //  - -1: cell is fully below the isosurface
+        //  -  0: cell is cut
+        //  - +1: cell is fully above the isosurface
+        label cellStatus_;
+
+        //- Boolean telling if subcell centre and volume have been calculated
+        bool subCellCentreAndVolumeCalculated_;
+
+        //- Boolean telling if isoface centre and area have been calculated
+        bool isoFaceCentreAndAreaCalculated_;
+
+
+    // Private Member Functions
+
+            void calcSubCellCentreAndVolume();
+
+            void calcIsoFaceCentreAndArea();
+
+            void calcIsoFacePointsFromEdges();
+
+
+public:
+
+    // Constructors
+
+        //- Construct from fvMesh and a scalarField
+        //  Length of scalarField should equal number of mesh points
+        isoCutCell(const fvMesh&, scalarField&);
+
+    // Static data
+
+        static int debug;
+
+
+    // Member functions
+
+        label calcSubCell(const label celli, const scalar isoValue);
+
+        const point& subCellCentre();
+
+        scalar subCellVolume();
+
+        const DynamicList<point>& isoFacePoints();
+
+        const point& isoFaceCentre();
+
+        const vector& isoFaceArea();
+
+        scalar volumeOfFluid();
+
+        scalar isoValue() const;
+
+        void clearStorage();
+
+        label vofCutCell
+        (
+            const label celli,
+            const scalar alpha1,
+            const scalar tol,
+            const label maxIter
+        );
+
+        void volumeOfFluid(volScalarField& alpha1, const scalar f0);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutFace/isoCutFace.C

@@ -0,0 +1,629 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+                isoAdvector | Copyright (C) 2016 DHI
+-------------------------------------------------------------------------------
+
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "isoCutFace.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+int Foam::isoCutFace::debug = 0;
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::isoCutFace::isoCutFace
+(
+    const fvMesh& mesh,
+    scalarField& f
+)
+:
+    mesh_(mesh),
+    f_(f),
+    firstEdgeCut_(-1),
+    lastEdgeCut_(-1),
+    firstFullySubmergedPoint_(-1),
+    nFullySubmergedPoints_(0),
+    subFaceCentre_(vector::zero),
+    subFaceArea_(vector::zero),
+    subFacePoints_(10),
+    surfacePoints_(4),
+    subFaceCentreAndAreaIsCalculated_(false)
+{
+    clearStorage();
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::isoCutFace::calcSubFaceCentreAndArea()
+{
+    const label nPoints = subFacePoints_.size();
+
+    // If the face is a triangle, do a direct calculation for efficiency
+    // and to avoid round-off error-related problems
+    if (nPoints == 3)
+    {
+        subFaceCentre_ = sum(subFacePoints_)/scalar(3);
+        subFaceArea_ =
+            0.5
+           *(
+                (subFacePoints_[1] - subFacePoints_[0])
+               ^(subFacePoints_[2] - subFacePoints_[0])
+            );
+    }
+    else if (nPoints > 0)
+    {
+        vector sumN = vector::zero;
+        scalar sumA = 0.0;
+        vector sumAc = vector::zero;
+        const point fCentre = sum(subFacePoints_)/scalar(nPoints);
+
+        for (label pi = 0; pi < nPoints; pi++)
+        {
+            const point& nextPoint = subFacePoints_[subFacePoints_.fcIndex(pi)];
+
+            vector c = subFacePoints_[pi] + nextPoint + fCentre;
+            vector n =
+                (nextPoint - subFacePoints_[pi])^(fCentre - subFacePoints_[pi]);
+            scalar a = magSqr(n);
+
+            sumN += n;
+            sumA += a;
+            sumAc += a*c;
+        }
+
+        // This is to deal with zero-area faces. Mark very small faces
+        // to be detected in e.g., processorPolyPatch.
+        if (sumA < ROOTVSMALL)
+        {
+            subFaceCentre_ = fCentre;
+            subFaceArea_ = vector::zero;
+        }
+        else
+        {
+            subFaceCentre_ = (1.0/3.0)*sumAc/sumA;
+            subFaceArea_ = 0.5*sumN;
+        }
+    }
+
+    subFaceCentreAndAreaIsCalculated_ = true;
+}
+
+
+Foam::label Foam::isoCutFace::calcSubFace
+(
+    const scalar isoValue,
+    const pointField& points,
+    const scalarField& f,
+    const labelList& pLabels
+)
+{
+    // Face status set to one of the values:
+    //  -1: face is fully below the isosurface
+    //   0: face is cut, i.e. has values larger and smaller than isoValue
+    //  +1: face is fully above the isosurface
+    label faceStatus;
+
+    label pl1 = pLabels[0];
+    scalar f1 = f[pl1];
+
+    // If vertex values are very close to isoValue lift them slightly to avoid
+    // dealing with the many special cases of a face being touched either at a
+    // single point, along an edge, or the entire face being on the surface.
+    if (mag(f1 - isoValue) < 10*SMALL)
+    {
+        f1 += sign(f1 - isoValue)*10*SMALL;
+    }
+
+    // Finding cut edges, the point along them where they are cut, and all fully
+    // submerged face points.
+    forAll(pLabels, pi)
+    {
+        label pl2 = pLabels[pLabels.fcIndex(pi)];
+        scalar f2 = f[pl2];
+        if (mag(f2 - isoValue) < 10*SMALL)
+        {
+            f2 += sign(f2 - isoValue)*10*SMALL;
+        }
+
+        if (f1 > isoValue)
+        {
+            nFullySubmergedPoints_ += 1;
+
+            if (f2 < isoValue)
+            {
+                lastEdgeCut_ = (isoValue - f1)/(f2 - f1);
+            }
+        }
+        else if (f1 < isoValue && f2 > isoValue)
+        {
+            if (firstFullySubmergedPoint_ == -1)
+            {
+                firstFullySubmergedPoint_ = pLabels.fcIndex(pi);
+
+                firstEdgeCut_ = (isoValue - f1)/(f2 - f1);
+            }
+            else
+            {
+                if (debug)
+                {
+                    const face fl(pLabels);
+
+                    WarningInFunction
+                        << "More than two face cuts for face " << fl
+                        << endl;
+
+                    Pout<< "Face values: f-isoValue = " << endl;
+                    forAll(fl, fpi)
+                    {
+                        Pout<< f[fl[fpi]] - isoValue << " ";
+                    }
+                    Pout<< " " << endl;
+                }
+            }
+        }
+        pl1 = pl2;
+        f1 = f2;
+    }
+
+    if (firstFullySubmergedPoint_ != -1)
+    {
+        // Face is cut
+        faceStatus = 0;
+        subFacePoints(points, pLabels);
+    }
+    else if (f1 < isoValue)
+    {
+        // Face entirely above isosurface
+        faceStatus = 1;
+    }
+    else
+    {
+        // Face entirely below isosurface
+        faceStatus = -1;
+    }
+
+    return faceStatus;
+}
+
+
+void Foam::isoCutFace::subFacePoints
+(
+    const pointField& points,
+    const labelList& pLabels
+)
+{
+    const label nPoints = pLabels.size();
+
+    surfacePoints(points, pLabels);
+
+    forAll(surfacePoints_, pi)
+    {
+        subFacePoints_.append(surfacePoints_[pi]);
+    }
+
+    for (label pi = 0; pi < nFullySubmergedPoints_; pi++)
+    {
+        subFacePoints_.append
+        (
+            points[pLabels[(firstFullySubmergedPoint_ + pi) % nPoints]]
+        );
+    }
+}
+
+
+void Foam::isoCutFace::surfacePoints
+(
+    const pointField& points,
+    const labelList& pLabels
+)
+{
+    const label nPoints = pLabels.size();
+
+    const label n = firstFullySubmergedPoint_ + nFullySubmergedPoints_;
+
+    label pl1 = pLabels[(n - 1) % nPoints];
+
+    label pl2 = pLabels[n % nPoints];
+
+    surfacePoints_.append
+    (
+        points[pl1] + lastEdgeCut_*(points[pl2] - points[pl1])
+    );
+
+    pl1 = pLabels[(firstFullySubmergedPoint_ - 1 + nPoints) % nPoints];
+    pl2 = pLabels[firstFullySubmergedPoint_];
+
+    surfacePoints_.append
+    (
+        points[pl1] + firstEdgeCut_*(points[pl2] - points[pl1])
+    );
+}
+
+
+// * * * * * * * * * * * Public Member Functions  * * * * * * * * * * * * * //
+
+Foam::label Foam::isoCutFace::calcSubFace
+(
+    const label faceI,
+    const scalar isoValue
+)
+{
+    clearStorage();
+    const labelList& pLabels = mesh_.faces()[faceI];
+    const pointField& points = mesh_.points();
+    return calcSubFace(isoValue, points, f_, pLabels);
+}
+
+
+Foam::label Foam::isoCutFace::calcSubFace
+(
+    const pointField& points,
+    const scalarField& f,
+    const scalar isoValue
+)
+{
+    clearStorage();
+    const labelList pLabels(identity(f.size()));
+    return calcSubFace(isoValue, points, f, pLabels);
+}
+
+
+const Foam::point& Foam::isoCutFace::subFaceCentre()
+{
+    if (!subFaceCentreAndAreaIsCalculated_)
+    {
+        calcSubFaceCentreAndArea();
+    }
+    return subFaceCentre_;
+}
+
+
+const Foam::vector& Foam::isoCutFace::subFaceArea()
+{
+    if (!subFaceCentreAndAreaIsCalculated_)
+    {
+        calcSubFaceCentreAndArea();
+    }
+    return subFaceArea_;
+}
+
+
+const Foam::DynamicList<Foam::point>& Foam::isoCutFace::subFacePoints() const
+{
+    return subFacePoints_;
+}
+
+
+const Foam::DynamicList<Foam::point>& Foam::isoCutFace::surfacePoints() const
+{
+    return surfacePoints_;
+}
+
+
+void Foam::isoCutFace::clearStorage()
+{
+    firstEdgeCut_ = -1;
+    lastEdgeCut_ = -1;
+    firstFullySubmergedPoint_ = -1;
+    nFullySubmergedPoints_ = 0;
+    subFaceCentre_ = vector::zero;
+    subFaceArea_ = vector::zero;
+    subFacePoints_.clear();
+    surfacePoints_.clear();
+    subFaceCentreAndAreaIsCalculated_ = false;
+}
+
+
+Foam::scalar Foam::isoCutFace::timeIntegratedArea
+(
+    const pointField& fPts,
+    const scalarField& pTimes,
+    const scalar dt,
+    const scalar magSf,
+    const scalar Un0
+)
+{
+    // Initialise time integrated area returned by this function
+    scalar tIntArea = 0.0;
+
+    // Finding ordering of vertex points
+    labelList order(pTimes.size());
+    sortedOrder(pTimes, order);
+    const scalar firstTime = pTimes[order.first()];
+    const scalar lastTime = pTimes[order.last()];
+
+    // Dealing with case where face is not cut by surface during time interval
+    // [0,dt] because face was already passed by surface
+    if (lastTime <= 0)
+    {
+        // If all face cuttings were in the past and cell is filling up (Un0>0)
+        // then face must be full during whole time interval
+        tIntArea = magSf*dt*pos(Un0);
+        return tIntArea;
+    }
+
+    // Dealing with case where face is not cut by surface during time interval
+    // [0, dt] because dt is too small for surface to reach closest face point
+    if (firstTime >= dt)
+    {
+        // If all cuttings are in the future but non of them within [0,dt] then
+        // if cell is filling up (Un0 > 0) face must be empty during whole time
+        // interval
+        tIntArea = magSf*dt*(1 - pos(Un0));
+        return tIntArea;
+    }
+
+    // If we reach this point in the code some part of the face will be swept
+    // during [tSmall, dt-tSmall]. However, it may be the case that there are no
+    // vertex times within the interval. This will happen sometimes for small
+    // time steps where both the initial and the final face-interface
+    // intersection line (FIIL) will be along the same two edges.
+
+    // Face-interface intersection line (FIIL) to be swept across face
+    DynamicList<point> FIIL(3);
+    // Submerged area at beginning of each sub time interval time
+    scalar initialArea = 0.0;
+    //Running time keeper variable for the integration process
+    scalar time = 0.0;
+
+    // Special treatment of first sub time interval
+    if (firstTime > 0)
+    {
+        // If firstTime > 0 the face is uncut in the time interval
+        // [0, firstTime] and hence fully submerged in fluid A or B.
+        // If Un0 > 0 cell is filling up and it must initially be empty.
+        // If Un0 < 0 cell must initially be full(y immersed in fluid A).
+        time = firstTime;
+        initialArea = magSf*(1.0 - pos(Un0));
+        tIntArea = initialArea*time;
+        cutPoints(fPts, pTimes, time, FIIL);
+    }
+    else
+    {
+        // If firstTime <= 0 then face is initially cut and we must
+        // calculate the initial submerged area and FIIL:
+        time = 0.0;
+        // Note: calcSubFace assumes well-defined 2-point FIIL!!!!
+        calcSubFace(fPts, -sign(Un0)*pTimes, time);
+        initialArea = mag(subFaceArea());
+        cutPoints(fPts, pTimes, time, FIIL);
+    }
+
+    // Making sorted array of all vertex times that are between max(0,firstTime)
+    // and dt and further than tSmall from the previous time.
+    DynamicList<scalar> sortedTimes(pTimes.size());
+    {
+        scalar prevTime = time;
+        const scalar tSmall = max(1e-6*dt, 10*SMALL);
+        forAll(order, ti)
+        {
+            const scalar timeI = pTimes[order[ti]];
+            if ( timeI > prevTime + tSmall && timeI <= dt)
+            {
+                sortedTimes.append(timeI);
+                prevTime = timeI;
+            }
+        }
+    }
+
+    // Sweeping all quadrilaterals corresponding to the intervals defined above
+    forAll(sortedTimes, ti)
+    {
+        const scalar newTime = sortedTimes[ti];
+        // New face-interface intersection line
+        DynamicList<point> newFIIL(3);
+        cutPoints(fPts, pTimes, newTime, newFIIL);
+
+        // quadrilateral area coefficients
+        scalar alpha = 0, beta = 0;
+        quadAreaCoeffs(FIIL, newFIIL, alpha, beta);
+        // Integration of area(t) = A*t^2+B*t from t = 0 to 1
+        tIntArea += (newTime - time)*
+            (initialArea + sign(Un0)*(alpha/3.0 + 0.5*beta));
+        // Adding quad area to submerged area
+        initialArea += sign(Un0)*(alpha + beta);
+
+        FIIL = newFIIL;
+        time = newTime;
+    }
+
+    // Special treatment of last time interval
+    if (lastTime > dt)
+    {
+        // FIIL will end up cutting the face at dt
+        // New face-interface intersection line
+        DynamicList<point> newFIIL(3);
+        cutPoints(fPts, pTimes, dt, newFIIL);
+
+        // quadrilateral area coefficients
+        scalar alpha = 0, beta = 0;
+        quadAreaCoeffs(FIIL, newFIIL, alpha, beta);
+        // Integration of area(t) = A*t^2+B*t from t = 0 to 1
+        tIntArea += (dt - time)*
+            (initialArea + sign(Un0)*(alpha/3.0 + 0.5*beta));
+    }
+    else
+    {
+        // FIIL will leave the face at lastTime and face will be fully in fluid
+        // A or fluid B in the time interval from lastTime to dt.
+        tIntArea += magSf*(dt - lastTime)*pos(Un0);
+    }
+
+    return tIntArea;
+}
+
+
+void Foam::isoCutFace::cutPoints
+(
+    const pointField& pts,
+    const scalarField& f,
+    const scalar f0,
+    DynamicList<point>& cutPoints
+)
+{
+    const label nPoints = pts.size();
+    scalar f1(f[0]);
+
+    // Snapping vertex value to f0 if very close (needed for 2D cases)
+    if (mag(f1 - f0) < 10*SMALL)
+    {
+        f1 = f0;
+    }
+
+    forAll(pts, pi)
+    {
+        label pi2 = (pi + 1) % nPoints;
+        scalar f2 = f[pi2];
+
+        // Snapping vertex value
+        if (mag(f2 - f0) < 10*SMALL)
+        {
+            f2 = f0;
+        }
+
+        if ((f1 < f0 && f2 > f0) || (f1 > f0 && f2 < f0))
+        {
+            const scalar s = (f0 - f1)/(f2 - f1);
+            cutPoints.append(pts[pi] + s*(pts[pi2] - pts[pi]));
+        }
+        else if (f1 == f0)
+        {
+            cutPoints.append(pts[pi]);
+        }
+        f1 = f2;
+    }
+
+    if (cutPoints.size() > 2)
+    {
+        WarningInFunction << "cutPoints = " << cutPoints << " for pts = " << pts
+            << ", f - f0 = " << f - f0 << " and f0 = " << f0 << endl;
+    }
+}
+
+
+void Foam::isoCutFace::quadAreaCoeffs
+(
+    const DynamicList<point>& pf0,
+    const DynamicList<point>& pf1,
+    scalar& alpha,
+    scalar& beta
+) const
+{
+    // Number of points in provided face-interface intersection lines
+    const label np0 = pf0.size();
+    const label np1 = pf1.size();
+
+    // quad area coeffs such that area(t) = alpha*t^2 + beta*t.
+    // With time interval normalised, we have full quadArea = alpha + beta
+    // and time integrated quad area = alpha/3 + beta/2;
+    alpha = 0.0;
+    beta = 0.0;
+
+    if (np0 && np1)
+    {
+        // Initialising quadrilateral vertices A, B, C and D
+        vector A(pf0[0]);
+        vector C(pf1[0]);
+        vector B(pf0[0]);
+        vector D(pf1[0]);
+
+        if (np0 == 2)
+        {
+            B = pf0[1];
+        }
+        else if (np0 > 2)
+        {
+            WarningInFunction << "Vertex face was cut at pf0 = " << pf0 << endl;
+        }
+
+        if (np1 == 2)
+        {
+            D = pf1[1];
+        }
+        else if (np1 > 2)
+        {
+            WarningInFunction << "Vertex face was cut at pf1 = " << pf1 << endl;
+        }
+
+        // Swapping pf1 points if pf0 and pf1 point in same general direction
+        // (because we want a quadrilateral ABCD where pf0 = AB and pf1 = CD)
+        if (((B - A) & (D - C)) > 0)
+        {
+            vector tmp = D;
+            D = C;
+            C = tmp;
+        }
+
+        // Defining local coordinates (xhat, yhat) for area integration of swept
+        // quadrilateral ABCD such that A = (0,0), B = (Bx,0), C = (Cx,Cy) and
+        // D = (Dx,Dy) with Cy = 0 and Dy > 0.
+
+        const scalar Bx = mag(B - A);
+
+        vector xhat(vector::zero);
+        if (Bx > 10*SMALL)
+        {
+            // If |AB| > 0 ABCD we use AB to define xhat
+            xhat = (B - A)/mag(B - A);
+        }
+        else if (mag(C - D) > 10*SMALL)
+        {
+            // If |AB| ~ 0 ABCD is a triangle ACD and we use CD for xhat
+            xhat = (C - D)/mag(C - D);
+        }
+        else
+        {
+            return;
+        }
+
+        // Defining vertical axis in local coordinates
+        vector yhat = D - A;
+        yhat -= ((yhat & xhat)*xhat);
+
+        if (mag(yhat) > 10*SMALL)
+        {
+            yhat /= mag(yhat);
+
+            const scalar Cx = (C - A) & xhat;
+            const scalar Cy = mag((C - A) & yhat);
+            const scalar Dx = (D - A) & xhat;
+            const scalar Dy = mag((D - A) & yhat);
+
+            // area = ((Cx - Bx)*Dy - Dx*Cy)/6.0 + 0.25*Bx*(Dy + Cy);
+            alpha = 0.5*((Cx - Bx)*Dy - Dx*Cy);
+            beta = 0.5*Bx*(Dy + Cy);
+        }
+    }
+    else
+    {
+        WarningInFunction
+            << "Vertex face was cut at " << pf0 << " and at " << pf1 << endl;
+    }
+}
+
+
+// ************************************************************************* //

src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutFace/isoCutFace.H

@@ -0,0 +1,200 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+                isoAdvector | Copyright (C) 2016 DHI
+-------------------------------------------------------------------------------
+
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::isoCutFace
+
+Description
+    Class for cutting a face, faceI, of an fvMesh, mesh_, at its intersection
+    with an isosurface defined by the mesh point values f_ and the isovalue,
+    isoValue_.
+
+    Reference:
+        \verbatim
+            Roenby, J., Bredmose, H. and Jasak, H. (2016).
+            A computational method for sharp interface advection
+            Royal Society Open Science, 3
+            doi 10.1098/rsos.160405
+        \endverbatim
+
+    Original code supplied by Johan Roenby, DHI (2016)
+
+SourceFiles
+    isoCutFace.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef isoCutFace_H
+#define isoCutFace_H
+
+#include "fvMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class isoCutFaces Declaration
+\*---------------------------------------------------------------------------*/
+
+class isoCutFace
+{
+    // Private data
+
+        //- Mesh whose cells and faces to cut at their intersection with an
+        //  isoface
+        const fvMesh& mesh_;
+
+        //- Isofunction values at mesh points. f_size() = mesh_.nPoints()
+        scalarField& f_;
+
+        //- Point along first cut edge where isosurface cuts edge
+        scalar firstEdgeCut_;
+
+        //- Point along last cut edge where isosurface cuts edge
+        scalar lastEdgeCut_;
+
+        //- Index in mesh_.faces()[faceI_] of first fully submerged (f > f0)
+        //  face point
+        label firstFullySubmergedPoint_;
+
+        //- Index in mesh_.faces()[faceI_] of last fully submerged (f > f0)
+        //  face point
+        label nFullySubmergedPoints_;
+
+        //- Storage for centre of subface
+        point subFaceCentre_;
+
+        //- Storage for area vector of subface
+        vector subFaceArea_;
+
+        //- Storage for subFacePoints
+        DynamicList<point> subFacePoints_;
+
+        //- Storage for subFacePoints
+        DynamicList<point> surfacePoints_;
+
+        //- Boolean telling if subface centre and area have been calculated
+        bool subFaceCentreAndAreaIsCalculated_;
+
+
+    // Private Member Functions
+
+        void calcSubFaceCentreAndArea();
+
+        //- Calculate cut points along edges of face with values f[pLabels]
+        //  Returns the face status, where:
+        //  -1: face is fully below the isosurface
+        //   0: face is cut, i.e. has values larger and smaller than isoValue
+        //  +1: face is fully above the isosurface
+        label calcSubFace
+        (
+            const scalar isoValue,
+            const pointField& points,
+            const scalarField& f,
+            const labelList& pLabels
+        );
+
+        void subFacePoints(const pointField& points, const labelList& pLabels);
+
+        void surfacePoints(const pointField& points, const labelList& pLabels);
+
+
+public:
+
+    // Constructors
+
+        //- Construct from fvMesh and a scalarField
+        //  Length of scalarField should equal number of mesh points
+        isoCutFace(const fvMesh& mesh, scalarField& f);
+
+    // Static data
+
+        static int debug;
+
+
+    // Member functions
+
+        //- Calculate cut points along edges of faceI
+        label calcSubFace(const label faceI, const scalar isoValue);
+
+        //- Calculate cut points along edges of face with values f
+        label calcSubFace
+        (
+            const pointField& points,
+            const scalarField& f,
+            const scalar isoValue
+        );
+
+        const point& subFaceCentre();
+
+        const vector& subFaceArea();
+
+        const DynamicList<point>& subFacePoints() const;
+
+        const DynamicList<point>& surfacePoints() const;
+
+        void clearStorage();
+
+        //- Calculate time integrated area for a face
+        scalar timeIntegratedArea
+        (
+            const pointField& fPts,
+            const scalarField& pTimes,
+            const scalar dt,
+            const scalar magSf,
+            const scalar Un0
+        );
+
+        void cutPoints
+        (
+            const pointField& pts,
+            const scalarField& f,
+            const scalar f0,
+            DynamicList<point>& cutPoints
+        );
+
+
+        void quadAreaCoeffs
+        (
+            const DynamicList<point>& pf0,
+            const DynamicList<point>& pf1,
+            scalar& alpha,
+            scalar& beta
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/functionObjects/field/Make/files

@@ -18,6 +18,8 @@ nearWallFields/findCellParticleCloud.C
 processorField/processorField.C
 readFields/readFields.C
 
+setFlow/setFlow.C
+
 streamLine/streamLine.C
 streamLine/streamLineBase.C
 streamLine/streamLineParticle.C

src/functionObjects/field/externalCoupled/externalCoupled.C

@@ -906,7 +906,7 @@ bool Foam::functionObjects::externalCoupled::read(const dictionary& dict)
 
     timeOut_ = dict.lookupOrDefault("timeOut", 100*waitInterval_);
     stateEnd_ =
-        stateEndNames_.lookupOrDefault("stateEnd", dict, stateEnd::REMOVE);
+        stateEndNames_.lookupOrDefault("stateEnd", dict, stateEnd::DONE);
 
 
     // Get names of all fvMeshes (and derived types)

src/functionObjects/field/externalCoupled/externalCoupled.H

@@ -32,12 +32,13 @@ Description
     an external application. The coupling is through plain text files
     where OpenFOAM boundary data is read/written as one line per face
     (data from all processors collated):
-
+    \verbatim
         # Patch: <patch name>
         <fld1> <fld2> .. <fldn>             //face0
         <fld1> <fld2> .. <fldn>             //face1
         ..
         <fld1> <fld2> .. <fldn>             //faceN
+    \endverbatim
 
     where the actual entries depend on the bc type:
     - mixed: value, snGrad, refValue, refGrad, valueFraction
@@ -47,29 +48,33 @@ Description
     These text files are located in a user specified communications directory
     which gets read/written on the master processor only. In the
     communications directory the structure will be
-
+    \verbatim
         <regionsName>/<patchGroup>/<fieldName>.[in|out]
+    \endverbatim
 
     (where regionsName is either the name of a single region or a composite
     of multiple region names)
 
     At start-up, the boundary creates a lock file, i.e..
-
+    \verbatim
         OpenFOAM.lock
+    \endverbatim
 
     ... to signal the external source to wait. During the functionObject
     execution the boundary values are written to files (one per region,
     per patch(group), per field), e.g.
-
+    \verbatim
         <regionsName>/<patchGroup>/<fieldName>.out
+    \endverbatim
 
     The lock file is then removed, instructing the external source to take
     control of the program execution. When ready, the external program
     should create the return values, e.g. to files
-
+    \verbatim
         <regionsName>/<patchGroup>/<fieldName>.in
+    \endverbatim
 
-    ... and then re-instate the lock file. The functionObject will then
+    ... and then reinstate the lock file. The functionObject will then
     read these values, apply them to the boundary conditions and pass
     program execution back to OpenFOAM.
 
@@ -98,21 +103,39 @@ Usage
     }
     \endverbatim
 
-    This reads/writes (on the master processor) the directory:
 
+    This reads/writes (on the master processor) the directory:
+    \verbatim
         comms/region0_region1/TPatchGroup/
+    \endverbatim
 
     with contents:
+    \verbatim
         patchPoints     (collected points)
         patchFaces      (collected faces)
         p.in            (input file of p, written by external application)
         T.out           (output file of T, written by OpenFOAM)
+    \endverbatim
 
     The patchPoints/patchFaces files denote the (collated) geometry
     which will be written if it does not exist yet or can be written as
     a preprocessing step using the createExternalCoupledPatchGeometry
     application.
 
+    The entries comprise:
+    \table
+        Property     | Description             | Required | Default value
+        type         | type name: externalCoupled | yes |
+        commsDir     | communication directory | yes |
+        waitInterval | wait interval in (s)    | no | 1
+        timeOut      | timeout in (s)          | no | 100*waitInterval
+        stateEnd     | Lockfile treatment on termination | no | done
+        initByExternal | initialization values supplied by external app | yes
+        calcFrequency | calculation frequency  | no | 1
+        regions  | the regions to couple | yes
+    \endtable
+
+
 SourceFiles
     externalCoupled.C
     externalCoupledTemplates.C

src/functionObjects/field/extractEulerianParticles/extractEulerianParticles/extractEulerianParticles.H

@@ -38,10 +38,10 @@ Usage
         ...
         faceZone        f0;
         nLocations      10;
-        alphaName       alpha.water;
+        alpha           alpha.water;
-        UName           U;
+        U               U;
-        rhoName         rho;
+        rho             rho;
-        phiName         phi;
+        phi             phi;
     }
     \endverbatim
 
@@ -50,10 +50,14 @@ Usage
         Property     | Description             | Required | Default value
         type         | type name: extractEulerianParticles | yes |
         faceZone     | Name of faceZone used as collection surface | yes |
-        nLocations   | Number of injection bins to generate  | yes |
+        alpha        | Name of phase indicator field | yes |
-        aplhaName    | Name of phase indicator field | yes   |
+        alphaThreshold | Threshold for alpha field | no | 0.1 |
-        rhoName      | Name of density field   |      yes |
+        nLocations   | Number of injection bins to generate  | no | 0 |
-        phiNane      | Name of flux field      |      yes |
+        U            | Name of velocity field  | no | U |
+        rho          | Name of density field   | no | rho |
+        phi          | Name of flux field      | no | phi |
+        minDiameter  | min diameter            | no | small |
+        maxDiameter  | max diameter            | no | great |
     \endtable
 
 SourceFiles

src/functionObjects/field/setFlow/setFlow.C

@@ -0,0 +1,448 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "setFlow.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+#include "fvcFlux.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvcSurfaceIntegrate.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+    defineTypeNameAndDebug(setFlow, 0);
+
+    addToRunTimeSelectionTable
+    (
+        functionObject,
+        setFlow,
+        dictionary
+    );
+}
+}
+
+
+const Foam::Enum<Foam::functionObjects::setFlow::modeType>
+Foam::functionObjects::setFlow::modeTypeNames
+{
+    { functionObjects::setFlow::modeType::FUNCTION, "function" },
+    { functionObjects::setFlow::modeType::ROTATION, "rotation" },
+    { functionObjects::setFlow::modeType::VORTEX2D, "vortex2D" },
+    { functionObjects::setFlow::modeType::VORTEX3D, "vortex3D" }
+};
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::functionObjects::setFlow::setPhi(const volVectorField& U)
+{
+    surfaceScalarField* phiptr =
+        mesh_.lookupObjectRefPtr<surfaceScalarField>(phiName_);
+
+    if (!phiptr)
+    {
+        return;
+    }
+
+    if (rhoName_ != "none")
+    {
+        const volScalarField* rhoptr =
+            mesh_.lookupObjectPtr<volScalarField>(rhoName_);
+
+        if (rhoptr)
+        {
+            const volScalarField& rho = *rhoptr;
+            *phiptr = fvc::flux(rho*U);
+        }
+        else
+        {
+            FatalErrorInFunction
+                << "Unable to find rho field'" << rhoName_
+                << "' in the mesh database.  Available fields are:"
+                << mesh_.names<volScalarField>()
+                << exit(FatalError);
+        }
+    }
+    else
+    {
+        *phiptr = fvc::flux(U);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::functionObjects::setFlow::setFlow
+(
+    const word& name,
+    const Time& runTime,
+    const dictionary& dict
+)
+:
+    fvMeshFunctionObject(name, runTime, dict),
+    UName_("U"),
+    rhoName_("none"),
+    phiName_("phi"),
+    mode_(modeType::FUNCTION),
+    reverseTime_(VGREAT),
+    scalePtr_(nullptr),
+    origin_(vector::zero),
+    R_(tensor::I),
+    omegaPtr_(nullptr),
+    velocityPtr_(nullptr)
+{
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::functionObjects::setFlow::~setFlow()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::functionObjects::setFlow::read(const dictionary& dict)
+{
+    if (fvMeshFunctionObject::read(dict))
+    {
+        Info<< name() << ":" << endl;
+        mode_ = modeTypeNames.read(dict.lookup("mode"));
+
+        Info<< "    operating mode: " << modeTypeNames[mode_] << endl;
+
+        if (dict.readIfPresent("U", UName_))
+        {
+            Info<< "    U field name: " << UName_ << endl;
+        }
+
+        if (dict.readIfPresent("rho", rhoName_))
+        {
+            Info<< "    rho field name: " << rhoName_ << endl;
+        }
+
+        if (dict.readIfPresent("phi", phiName_))
+        {
+            Info<< "    phi field name: " << phiName_ << endl;
+        }
+
+        if (dict.readIfPresent("reverseTime", reverseTime_))
+        {
+            Info<< "    reverse flow direction at time: " << reverseTime_
+                << endl;
+            reverseTime_ = mesh_.time().userTimeToTime(reverseTime_);
+        }
+
+        // Scaling is applied across all modes
+        scalePtr_ = Function1<scalar>::New("scale", dict);
+
+        switch (mode_)
+        {
+            case modeType::FUNCTION:
+            {
+                velocityPtr_ = Function1<vector>::New("velocity", dict);
+                break;
+            }
+            case modeType::ROTATION:
+            {
+                omegaPtr_ = Function1<scalar>::New("omega", dict);
+                dict.lookup("origin") >> origin_;
+                vector refDir(dict.lookup("refDir"));
+                refDir /= mag(refDir) + ROOTVSMALL;
+                vector axis(dict.lookup("axis"));
+                axis /= mag(axis) + ROOTVSMALL;
+                R_ = tensor(refDir, axis, refDir^axis);
+                break;
+            }
+            case modeType::VORTEX2D:
+            case modeType::VORTEX3D:
+            {
+                dict.lookup("origin") >> origin_;
+                vector refDir(dict.lookup("refDir"));
+                refDir /= mag(refDir) + ROOTVSMALL;
+                vector axis(dict.lookup("axis"));
+                axis /= mag(axis) + ROOTVSMALL;
+                R_ = tensor(refDir, axis, refDir^axis);
+                break;
+            }
+        }
+
+        Info<< endl;
+
+        return true;
+    }
+
+    return false;
+}
+
+
+bool Foam::functionObjects::setFlow::execute()
+{
+    volVectorField* Uptr = mesh_.lookupObjectRefPtr<volVectorField>(UName_);
+
+    surfaceScalarField* phiptr =
+        mesh_.lookupObjectRefPtr<surfaceScalarField>(phiName_);
+
+    Log << nl << name() << ":" << nl;
+
+    if (!Uptr || !phiptr)
+    {
+        Info<< "    Either field " << UName_ << " or " << phiName_
+            << " not found in the mesh database" << nl;
+
+        return true;
+    }
+
+    const scalar t = mesh_.time().timeOutputValue();
+
+    Log << "    setting " << UName_ << " and " << phiName_ << nl;
+
+    volVectorField& U = *Uptr;
+    surfaceScalarField& phi = *phiptr;
+
+    switch (mode_)
+    {
+        case modeType::FUNCTION:
+        {
+            const vector Uc = velocityPtr_->value(t);
+            U == dimensionedVector("Uc", dimVelocity, Uc);
+            U.correctBoundaryConditions();
+            setPhi(U);
+
+            break;
+        }
+        case modeType::ROTATION:
+        {
+            const volVectorField& C = mesh_.C();
+            const volVectorField d
+            (
+                typeName + ":d",
+                C - dimensionedVector("origin", dimLength, origin_)
+            );
+            const scalarField x(d.component(vector::X));
+            const scalarField z(d.component(vector::Z));
+
+            const scalar omega = omegaPtr_->value(t);
+            vectorField& Uc = U.primitiveFieldRef();
+            Uc.replace(vector::X, -omega*z);
+            Uc.replace(vector::Y, scalar(0));
+            Uc.replace(vector::Z, omega*x);
+
+            volVectorField::Boundary& Ubf = U.boundaryFieldRef();
+            forAll(Ubf, patchi)
+            {
+                vectorField& Uf = Ubf[patchi];
+                if (Uf.size())
+                {
+                    const vectorField& Cp = C.boundaryField()[patchi];
+                    const vectorField dp(Cp - origin_);
+                    const scalarField xp(dp.component(vector::X));
+                    const scalarField zp(dp.component(vector::Z));
+                    Uf.replace(vector::X, -omega*zp);
+                    Uf.replace(vector::Y, scalar(0));
+                    Uf.replace(vector::Z, omega*xp);
+                }
+            }
+
+            U = U & R_;
+            U.correctBoundaryConditions();
+            setPhi(U);
+
+            break;
+        }
+        case modeType::VORTEX2D:
+        {
+            const scalar pi = Foam::constant::mathematical::pi;
+
+            const volVectorField& C = mesh_.C();
+
+            const volVectorField d
+            (
+                typeName + ":d",
+                C - dimensionedVector("origin", dimLength, origin_)
+            );
+            const scalarField x(d.component(vector::X));
+            const scalarField z(d.component(vector::Z));
+
+            vectorField& Uc = U.primitiveFieldRef();
+            Uc.replace(vector::X, -sin(2*pi*z)*sqr(sin(pi*x)));
+            Uc.replace(vector::Y, scalar(0));
+            Uc.replace(vector::Z, sin(2*pi*x)*sqr(sin(pi*z)));
+
+            U = U & R_;
+
+            // Calculating incompressible flux based on stream function
+            // Note: R_ rotation not implemented in phi calculation
+            const scalarField xp(mesh_.points().component(0) - origin_[0]);
+            const scalarField yp(mesh_.points().component(1) - origin_[1]);
+            const scalarField zp(mesh_.points().component(2) - origin_[2]);
+            const scalarField psi((1.0/pi)*sqr(sin(pi*xp))*sqr(sin(pi*zp)));
+
+            scalarField& phic = phi.primitiveFieldRef();
+            forAll(phic, fi)
+            {
+                phic[fi] = 0;
+                const face& f = mesh_.faces()[fi];
+                const label nPoints = f.size();
+
+                forAll(f, fpi)
+                {
+                    const label p1 = f[fpi];
+                    const label p2 = f[(fpi + 1) % nPoints];
+                    phic[fi] += 0.5*(psi[p1] + psi[p2])*(yp[p2] - yp[p1]);
+                }
+            }
+
+            surfaceScalarField::Boundary& phibf = phi.boundaryFieldRef();
+            forAll(phibf, patchi)
+            {
+                scalarField& phif = phibf[patchi];
+                const label start = mesh_.boundaryMesh()[patchi].start();
+
+                forAll(phif, fi)
+                {
+                    phif[fi] = 0;
+                    const face& f = mesh_.faces()[start + fi];
+                    const label nPoints = f.size();
+
+                    forAll(f, fpi)
+                    {
+                        const label p1 = f[fpi];
+                        const label p2 = f[(fpi + 1) % nPoints];
+                        phif[fi] += 0.5*(psi[p1] + psi[p2])*(yp[p2] - yp[p1]);
+                    }
+                }
+            }
+
+            break;
+        }
+        case modeType::VORTEX3D:
+        {
+            const scalar pi = Foam::constant::mathematical::pi;
+            const volVectorField& C = mesh_.C();
+
+            const volVectorField d
+            (
+                typeName + ":d",
+                C - dimensionedVector("origin", dimLength, origin_)
+            );
+            const scalarField x(d.component(vector::X));
+            const scalarField y(d.component(vector::Y));
+            const scalarField z(d.component(vector::Z));
+
+            vectorField& Uc = U.primitiveFieldRef();
+            Uc.replace(vector::X, 2*sqr(sin(pi*x))*sin(2*pi*y)*sin(2*pi*z));
+            Uc.replace(vector::Y, -sin(2*pi*x)*sqr(sin(pi*y))*sin(2*pi*z));
+            Uc.replace(vector::Z, -sin(2*pi*x)*sin(2*pi*y)*sqr(sin(pi*z)));
+
+            U = U & R_;
+            U.correctBoundaryConditions();
+
+            // Calculating phi
+            // Note: R_ rotation not implemented in phi calculation
+            const vectorField Cf(mesh_.Cf().primitiveField() - origin_);
+            const scalarField Xf(Cf.component(vector::X));
+            const scalarField Yf(Cf.component(vector::Y));
+            const scalarField Zf(Cf.component(vector::Z));
+            vectorField Uf(Xf.size());
+            Uf.replace(0, 2*sqr(sin(pi*Xf))*sin(2*pi*Yf)*sin(2*pi*Zf));
+            Uf.replace(1, -sin(2*pi*Xf)*sqr(sin(pi*Yf))*sin(2*pi*Zf));
+            Uf.replace(2, -sin(2*pi*Xf)*sin(2*pi*Yf)*sqr(sin(pi*Zf)));
+
+            scalarField& phic = phi.primitiveFieldRef();
+            const vectorField& Sfc = mesh_.Sf().primitiveField();
+            phic = Uf & Sfc;
+
+            surfaceScalarField::Boundary& phibf = phi.boundaryFieldRef();
+            const surfaceVectorField::Boundary& Sfbf =
+                mesh_.Sf().boundaryField();
+            const surfaceVectorField::Boundary& Cfbf =
+                mesh_.Cf().boundaryField();
+
+            forAll(phibf, patchi)
+            {
+                scalarField& phif = phibf[patchi];
+                const vectorField& Sff = Sfbf[patchi];
+                const vectorField& Cff = Cfbf[patchi];
+                const scalarField xf(Cff.component(vector::X));
+                const scalarField yf(Cff.component(vector::Y));
+                const scalarField zf(Cff.component(vector::Z));
+                vectorField Uf(xf.size());
+                Uf.replace(0, 2*sqr(sin(pi*xf))*sin(2*pi*yf)*sin(2*pi*zf));
+                Uf.replace(1, -sin(2*pi*xf)*sqr(sin(pi*yf))*sin(2*pi*zf));
+                Uf.replace(2, -sin(2*pi*xf)*sin(2*pi*yf)*sqr(sin(pi*zf)));
+                phif = Uf & Sff;
+            }
+
+            break;
+        }
+    }
+
+    if (t > reverseTime_)
+    {
+        Log << "    flow direction: reverse" << nl;
+        U.negate();
+        phi.negate();
+    }
+
+    // Apply scaling
+    const scalar s = scalePtr_->value(t);
+    U *= s;
+    phi *= s;
+
+    U.correctBoundaryConditions();
+
+    const scalarField sumPhi(fvc::surfaceIntegrate(phi));
+    Log << "    Continuity error: max(mag(sum(phi))) = "
+        << gMax(mag(sumPhi)) << nl << endl;
+
+    return true;
+}
+
+
+bool Foam::functionObjects::setFlow::write()
+{
+    const auto& Uptr = mesh_.lookupObjectRefPtr<volVectorField>(UName_);
+    if (Uptr)
+    {
+        Uptr->write();
+    }
+
+    const auto& phiptr = mesh_.lookupObjectRefPtr<surfaceScalarField>(phiName_);
+    if (phiptr)
+    {
+        phiptr->write();
+    }
+
+    return true;
+}
+
+
+// ************************************************************************* //

src/functionObjects/field/setFlow/setFlow.H

@@ -0,0 +1,199 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::functionObjects::setFlow
+
+Group
+    grpFieldFunctionObjects
+
+Description
+    Provides options to set the velocity and flux fields as a function of time
+
+    Useful for testing advection behaviour of numerical schemes by e.g.
+    imposing solid body rotation, vortex flows.  All types include a scaling
+    Foam::Function1 type enabling the strength of the transformation to vary as
+    a function of time.
+
+Usage
+    Example of function object specification:
+    \verbatim
+    setFlow1
+    {
+        type            setFlow;
+        libs            ("libfieldFunctionObjects.so");
+        ...
+        mode            rotation;
+        scale           1;
+        reverseTime     1;
+        omega           6.28318530718;
+        origin          (0.5 0 0.5);
+        refDir          (1 0 0);
+        axis            (0 1 0);
+    }
+    \endverbatim
+
+    Where the entries comprise:
+    \table
+        Property     | Description                | Required | Default value
+        type         | Type name: setFlow         | yes      |
+        U            | Name of velocity field     | no       | U
+        rho          | Name of density field      | no       | none
+        phi          | Name of flux field         | no       | phi
+        mode         | operating mode - see below | yes      |
+    \endtable
+
+    Available \c mode types include:
+    - function
+    - rortation
+    - vortex2D
+    - vortex3D
+
+
+See also
+    Foam::functionObjects::fvMeshFunctionObject
+
+SourceFiles
+    setFlow.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef functionObjects_setFlow_H
+#define functionObjects_setFlow_H
+
+#include "fvMeshFunctionObject.H"
+#include "Function1.H"
+#include "Enum.H"
+#include "point.H"
+#include "volFieldsFwd.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class setFlow Declaration
+\*---------------------------------------------------------------------------*/
+
+class setFlow
+:
+    public fvMeshFunctionObject
+{
+    enum class modeType
+    {
+        FUNCTION,
+        ROTATION,
+        VORTEX2D,
+        VORTEX3D
+    };
+
+   static const Foam::Enum<modeType> modeTypeNames;
+
+
+   // Private Data
+
+        //- Name of vecloity field, default = "U"
+        word UName_;
+
+        //- Name of density field, default = "none"
+        word rhoName_;
+
+        //- Name of flux field, default = "phi"
+        word phiName_;
+
+        //- Operating mode
+        modeType mode_;
+
+        //- Reverse time
+        scalar reverseTime_;
+
+        //- Scaling function
+        autoPtr<Function1<scalar>> scalePtr_;
+
+
+        // Rotation
+
+            //- Origin
+            point origin_;
+
+            //- Rotation tensor for rotational mode
+            tensor R_;
+
+            //- Rotational speed function
+            autoPtr<Function1<scalar>> omegaPtr_;
+
+
+        // Function
+
+            //- Velocity function
+            autoPtr<Function1<vector>> velocityPtr_;
+
+
+    // Private Member Functions
+
+        //- Set the flux field
+        void setPhi(const volVectorField& U);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("setFlow");
+
+
+    // Constructors
+
+        //- Construct from Time and dictionary
+        setFlow
+        (
+            const word& name,
+            const Time& runTime,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~setFlow();
+
+
+    virtual bool read(const dictionary& dict);
+
+    virtual bool execute();
+
+    virtual bool write();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace functionObjects
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -20,7 +20,6 @@ internalField   uniform 273;
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
     #includeEtc "caseDicts/setConstraintTypes"
 
     "(walls|hole)"

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

@@ -20,14 +20,8 @@ internalField   uniform 0;
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
     #includeEtc "caseDicts/setConstraintTypes"
 
-    overset
-    {
-        type        overset;
-    }
-
     ".*"
     {
         type        zeroGradient;

tutorials/combustion/XiDyMFoam/annularCombustorTurbine/0.orig/Su

@@ -21,7 +21,7 @@ internalField   uniform 0.135;
 
 boundaryField
 {
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
+    #includeEtc "caseDicts/setConstraintTypes"
 
     inlet
     {

tutorials/combustion/XiDyMFoam/annularCombustorTurbine/0.orig/T

@@ -22,7 +22,7 @@ internalField   uniform 2400;
 
 boundaryField
 {
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
+    #includeEtc "caseDicts/setConstraintTypes"
 
     inlet
     {

tutorials/combustion/XiDyMFoam/annularCombustorTurbine/0.orig/Tu

@@ -21,7 +21,7 @@ internalField   uniform 573;
 
 boundaryField
 {
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
+    #includeEtc "caseDicts/setConstraintTypes"
 
     inlet
     {

tutorials/combustion/XiDyMFoam/annularCombustorTurbine/0.orig/U

@@ -21,7 +21,7 @@ internalField   uniform (4 0 0);
 
 boundaryField
 {
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
+    #includeEtc "caseDicts/setConstraintTypes"
 
     inlet
     {

tutorials/combustion/XiDyMFoam/annularCombustorTurbine/0.orig/Xi

@@ -21,7 +21,7 @@ internalField   uniform 1;
 
 boundaryField
 {
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
+    #includeEtc "caseDicts/setConstraintTypes"
 
     inlet
     {

tutorials/combustion/XiDyMFoam/annularCombustorTurbine/0.orig/alphat

@@ -21,7 +21,7 @@ internalField   uniform 0;
 
 boundaryField
 {
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
+    #includeEtc "caseDicts/setConstraintTypes"
 
     inlet
     {

tutorials/combustion/XiDyMFoam/annularCombustorTurbine/0.orig/b

@@ -21,7 +21,7 @@ internalField   uniform 0;
 
 boundaryField
 {
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
+    #includeEtc "caseDicts/setConstraintTypes"
 
     inlet
     {

tutorials/combustion/XiDyMFoam/annularCombustorTurbine/0.orig/ft

@@ -21,7 +21,7 @@ internalField   uniform 0.0834168596800118;
 
 boundaryField
 {
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
+    #includeEtc "caseDicts/setConstraintTypes"
 
     inlet
     {

tutorials/combustion/XiDyMFoam/annularCombustorTurbine/0.orig/k

@@ -21,7 +21,7 @@ internalField   uniform 6;
 
 boundaryField
 {
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
+    #includeEtc "caseDicts/setConstraintTypes"
 
     inlet
     {

tutorials/combustion/XiDyMFoam/annularCombustorTurbine/0.orig/nut

@@ -21,7 +21,7 @@ internalField   uniform 0;
 
 boundaryField
 {
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
+    #includeEtc "caseDicts/setConstraintTypes"
 
     inlet
     {

tutorials/combustion/XiDyMFoam/annularCombustorTurbine/0.orig/p

@@ -21,7 +21,7 @@ internalField   uniform 10e5;
 
 boundaryField
 {
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
+    #includeEtc "caseDicts/setConstraintTypes"
 
     inlet
     {

tutorials/combustion/XiDyMFoam/annularCombustorTurbine/0.orig/pPotential

@@ -21,7 +21,7 @@ internalField   uniform 0;
 
 boundaryField
 {
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
+    #includeEtc "caseDicts/setConstraintTypes"
 
     inlet
     {

tutorials/combustion/XiDyMFoam/annularCombustorTurbine/system/meshQualityDict

@@ -16,7 +16,7 @@ FoamFile
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 // Include defaults parameters from master dictionary
-#include "$WM_PROJECT_DIR/etc/caseDicts/meshQualityDict"
+#includeEtc "caseDicts/meshQualityDict"
 
 maxNonOrtho 55;
 

tutorials/compressible/rhoSimpleFoam/gasMixing/injectorPipe/system/meshQualityDict

@@ -16,7 +16,7 @@ FoamFile
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 // Include defaults parameters from master dictionary
-#include "$WM_PROJECT_DIR/etc/caseDicts/meshQualityDict"
+#includeEtc "caseDicts/meshQualityDict"
 
 maxNonOrtho 55;
 

tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/0/T

@@ -20,7 +20,6 @@ internalField   uniform 265;
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
     #includeEtc "caseDicts/setConstraintTypes"
 
     wall

tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/0/U

@@ -20,7 +20,6 @@ internalField   uniform (0 0 0);
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
     #includeEtc "caseDicts/setConstraintTypes"
 
     wall

tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/0/alphat

@@ -21,7 +21,6 @@ internalField   uniform 0;
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
     #includeEtc "caseDicts/setConstraintTypes"
 
     wall

tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/0/epsilon

@@ -21,7 +21,6 @@ internalField   uniform 0.01;
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
     #includeEtc "caseDicts/setConstraintTypes"
 
     wall

tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/0/k

@@ -21,7 +21,6 @@ internalField   uniform 0.1;
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
     #includeEtc "caseDicts/setConstraintTypes"
 
     wall

tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/0/nut

@@ -21,7 +21,6 @@ internalField   uniform 0;
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
     #includeEtc "caseDicts/setConstraintTypes"
 
     wall

tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/0/p

@@ -20,7 +20,6 @@ internalField   uniform 0;
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
     #includeEtc "caseDicts/setConstraintTypes"
 
     wall

tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/0/p_rgh

@@ -20,7 +20,6 @@ internalField   uniform 0;
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
     #includeEtc "caseDicts/setConstraintTypes"
 
     wall

tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.orig/U

@@ -21,7 +21,6 @@ internalField   uniform (0.1 0 0);
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
     #includeEtc "caseDicts/setConstraintTypes"
 
     floor

tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.orig/alphat

@@ -21,7 +21,6 @@ internalField   uniform 0;
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
     #includeEtc "caseDicts/setConstraintTypes"
 
     floor

tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.orig/epsilon

@@ -21,7 +21,6 @@ internalField   uniform 0.01;
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
     #includeEtc "caseDicts/setConstraintTypes"
 
     floor

tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.orig/k

@@ -21,7 +21,6 @@ internalField   uniform 0.1;
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
     #includeEtc "caseDicts/setConstraintTypes"
 
     floor

tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.orig/nut

@@ -21,7 +21,6 @@ internalField   uniform 0;
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
     #includeEtc "caseDicts/setConstraintTypes"
 
     floor

tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.orig/p

@@ -21,7 +21,6 @@ internalField   uniform 101325;
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
     #includeEtc "caseDicts/setConstraintTypes"
 
     floor

tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.orig/p_rgh

@@ -21,7 +21,6 @@ internalField   uniform 101325;
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
     #includeEtc "caseDicts/setConstraintTypes"
 
     floor

tutorials/incompressible/overPimpleDyMFoam/cylinder/cylinderAndBackground/0.orig/U

@@ -22,8 +22,7 @@ internalField   uniform $flowVelocity;
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
+    #includeEtc "caseDicts/setConstraintTypes"
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
 
     overset
     {

tutorials/incompressible/overPimpleDyMFoam/cylinder/cylinderAndBackground/0.orig/epsilon

@@ -22,8 +22,7 @@ internalField   uniform $turbulentEpsilon;
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
+    #includeEtc "caseDicts/setConstraintTypes"
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
 
     overset
     {

tutorials/incompressible/overPimpleDyMFoam/cylinder/cylinderAndBackground/0.orig/k

@@ -22,8 +22,7 @@ internalField   uniform $turbulentKE;
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
+    #includeEtc "caseDicts/setConstraintTypes"
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
 
     overset
     {

tutorials/incompressible/overPimpleDyMFoam/cylinder/cylinderAndBackground/0.orig/nut

@@ -20,8 +20,7 @@ internalField   uniform 0;
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
+    #includeEtc "caseDicts/setConstraintTypes"
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
 
     overset
     {

tutorials/incompressible/overPimpleDyMFoam/cylinder/cylinderAndBackground/0.orig/p

@@ -22,8 +22,7 @@ internalField   uniform $pressure;
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
+    #includeEtc "caseDicts/setConstraintTypes"
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
 
     overset
     {

tutorials/incompressible/overPimpleDyMFoam/cylinder/cylinderAndBackground/0.orig/zoneID

@@ -20,8 +20,7 @@ internalField   uniform 0;
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
+    #includeEtc "caseDicts/setConstraintTypes"
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
 
     overset
     {

tutorials/incompressible/overPimpleDyMFoam/simpleRotor/0.orig/U

@@ -20,8 +20,7 @@ internalField   uniform (0 0 0);
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
+    #includeEtc "caseDicts/setConstraintTypes"
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
 
     walls
     {

tutorials/incompressible/overPimpleDyMFoam/simpleRotor/0.orig/epsilon

@@ -21,8 +21,7 @@ internalField   uniform 200;
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
+    #includeEtc "caseDicts/setConstraintTypes"
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
 
     "(walls|hole)"
     {

tutorials/incompressible/overPimpleDyMFoam/simpleRotor/0.orig/k

@@ -20,8 +20,7 @@ internalField   uniform 0.2;
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
+    #includeEtc "caseDicts/setConstraintTypes"
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
 
     "(walls|hole)"
     {

tutorials/incompressible/overPimpleDyMFoam/simpleRotor/0.orig/nuTilda

@@ -20,8 +20,7 @@ internalField   uniform 0;
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
+    #includeEtc "caseDicts/setConstraintTypes"
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
 
     overset
     {

tutorials/incompressible/overPimpleDyMFoam/simpleRotor/0.orig/nut

@@ -20,8 +20,7 @@ internalField   uniform 0;
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
+    #includeEtc "caseDicts/setConstraintTypes"
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
 
     "(walls|hole)"
     {

tutorials/incompressible/overPimpleDyMFoam/simpleRotor/0.orig/p

@@ -20,8 +20,7 @@ internalField   uniform 0;
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
+    #includeEtc "caseDicts/setConstraintTypes"
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
 
     "(walls|hole)"
     {

tutorials/incompressible/overPimpleDyMFoam/simpleRotor/0.orig/zoneID

@@ -20,8 +20,7 @@ internalField   uniform 0;
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
+    #includeEtc "caseDicts/setConstraintTypes"
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
 
     ".*"
     {

tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/0.orig/U

@@ -20,8 +20,7 @@ internalField   uniform (0.1 0 0);
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
+    #includeEtc "caseDicts/setConstraintTypes"
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
 
     walls
     {

tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/0.orig/p

@@ -20,8 +20,7 @@ internalField   uniform 0;
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
+    #includeEtc "caseDicts/setConstraintTypes"
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
 
     "(walls|hole)"
     {

tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/0.orig/zoneID

@@ -20,8 +20,7 @@ internalField   uniform 0;
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
+    #includeEtc "caseDicts/setConstraintTypes"
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
 
     overset
     {

tutorials/incompressible/overSimpleFoam/aeroFoil/background_overset/0.orig/U

@@ -22,8 +22,7 @@ internalField   uniform $flowVelocity;
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
+    #includeEtc "caseDicts/setConstraintTypes"
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
 
     walls
     {

tutorials/incompressible/overSimpleFoam/aeroFoil/background_overset/0.orig/epsilon

@@ -22,8 +22,7 @@ internalField   uniform $turbulentEpsilon;
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
+    #includeEtc "caseDicts/setConstraintTypes"
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
 
     wall
     {

tutorials/incompressible/overSimpleFoam/aeroFoil/background_overset/0.orig/k

@@ -22,8 +22,7 @@ internalField   uniform $turbulentKE;
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
+    #includeEtc "caseDicts/setConstraintTypes"
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
 
     wall
     {

tutorials/incompressible/overSimpleFoam/aeroFoil/background_overset/0.orig/nuTilda

@@ -20,8 +20,7 @@ internalField   uniform 0;
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
+    #includeEtc "caseDicts/setConstraintTypes"
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
 
     ".*"
     {

tutorials/incompressible/overSimpleFoam/aeroFoil/background_overset/0.orig/nut

@@ -20,8 +20,7 @@ internalField   uniform 0;
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
+    #includeEtc "caseDicts/setConstraintTypes"
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
 
     wall
     {

tutorials/incompressible/overSimpleFoam/aeroFoil/background_overset/0.orig/omega

@@ -22,8 +22,7 @@ internalField   uniform $turbulentOmega;
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
+    #includeEtc "caseDicts/setConstraintTypes"
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
 
     wall
     {

tutorials/incompressible/overSimpleFoam/aeroFoil/background_overset/0.orig/p

@@ -22,8 +22,7 @@ internalField   uniform $pressure;
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
+    #includeEtc "caseDicts/setConstraintTypes"
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
 
     wall
     {

tutorials/incompressible/overSimpleFoam/aeroFoil/background_overset/0.orig/zoneID

@@ -20,8 +20,7 @@ internalField   uniform 0;
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
+    #includeEtc "caseDicts/setConstraintTypes"
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
 
     ".*"
     {

tutorials/incompressible/pimpleDyMFoam/oscillatingInletPeriodicAMI2D/0/U

@@ -21,7 +21,7 @@ internalField   uniform (0 0 0);
 
 boundaryField
 {
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
+    #includeEtc "caseDicts/setConstraintTypes"
 
     inlet
     {

tutorials/incompressible/pimpleDyMFoam/oscillatingInletPeriodicAMI2D/0/epsilon

@@ -21,7 +21,7 @@ internalField   uniform 1.8e-3;
 
 boundaryField
 {
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
+    #includeEtc "caseDicts/setConstraintTypes"
 
     inlet
     {

tutorials/incompressible/pimpleDyMFoam/oscillatingInletPeriodicAMI2D/0/k

@@ -21,7 +21,7 @@ internalField   uniform 3.75e-3;
 
 boundaryField
 {
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
+    #includeEtc "caseDicts/setConstraintTypes"
 
     inlet
     {

tutorials/incompressible/pimpleDyMFoam/oscillatingInletPeriodicAMI2D/0/nut

@@ -21,7 +21,7 @@ internalField   uniform 0;
 
 boundaryField
 {
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
+    #includeEtc "caseDicts/setConstraintTypes"
 
     inlet
     {

tutorials/incompressible/pimpleDyMFoam/oscillatingInletPeriodicAMI2D/0/p

@@ -21,7 +21,7 @@ internalField   uniform 0;
 
 boundaryField
 {
-    #include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
+    #includeEtc "caseDicts/setConstraintTypes"
 
     inlet
     {

tutorials/incompressible/pimpleDyMFoam/propeller/0.orig/U

@@ -21,7 +21,6 @@ internalField   uniform (0 0 0);
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
     #includeEtc "caseDicts/setConstraintTypes"
 
     inlet

tutorials/incompressible/pimpleDyMFoam/propeller/0.orig/epsilon

@@ -21,7 +21,6 @@ internalField   uniform 0.0495;
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
     #includeEtc "caseDicts/setConstraintTypes"
 
     inlet

tutorials/incompressible/pimpleDyMFoam/propeller/0.orig/k

@@ -21,7 +21,6 @@ internalField   uniform 0.06;
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
     #includeEtc "caseDicts/setConstraintTypes"
 
     inlet

tutorials/incompressible/pimpleDyMFoam/propeller/0.orig/nut

@@ -21,7 +21,6 @@ internalField   uniform 0;
 
 boundaryField
 {
-    //- Set patchGroups for constraint patches
     #includeEtc "caseDicts/setConstraintTypes"
 
     inlet