diff --git a/applications/solvers/multiphase/interIsoFoam/Make/files b/applications/solvers/multiphase/interIsoFoam/Make/files
new file mode 100644
index 0000000000000000000000000000000000000000..10f61a688541da0b2ba1d37a7aa98ab384ccbd2a
--- /dev/null
+++ b/applications/solvers/multiphase/interIsoFoam/Make/files
@@ -0,0 +1,3 @@
+interIsoFoam.C
+
+EXE = $(FOAM_APPBIN)/interIsoFoam
diff --git a/applications/solvers/multiphase/interIsoFoam/Make/options b/applications/solvers/multiphase/interIsoFoam/Make/options
new file mode 100644
index 0000000000000000000000000000000000000000..cfe4a920729ed8bd1497f89de334caf51b85a312
--- /dev/null
+++ b/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
diff --git a/applications/solvers/multiphase/interIsoFoam/UEqn.H b/applications/solvers/multiphase/interIsoFoam/UEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..77d1dcd83e8c665fa8127bfadb8d59a94222bc61
--- /dev/null
+++ b/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);
+ }
diff --git a/applications/solvers/multiphase/interIsoFoam/alphaControls.H b/applications/solvers/multiphase/interIsoFoam/alphaControls.H
new file mode 100644
index 0000000000000000000000000000000000000000..92a8fdf9510302a9c6d032a2c5c6962acf66059a
--- /dev/null
+++ b/applications/solvers/multiphase/interIsoFoam/alphaControls.H
@@ -0,0 +1,4 @@
+const dictionary& alphaControls = mesh.solverDict(alpha1.name());
+
+label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")));
+
diff --git a/applications/solvers/multiphase/interIsoFoam/alphaCourantNo.H b/applications/solvers/multiphase/interIsoFoam/alphaCourantNo.H
new file mode 100644
index 0000000000000000000000000000000000000000..a084155c8ae70e1e0b4a6cc9ae9b6f5fba828235
--- /dev/null
+++ b/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;
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/interIsoFoam/alphaEqn.H b/applications/solvers/multiphase/interIsoFoam/alphaEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..cefabfda7354a26c6745131a630333c66596b83f
--- /dev/null
+++ b/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;
diff --git a/applications/solvers/multiphase/interIsoFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/interIsoFoam/alphaEqnSubCycle.H
new file mode 100644
index 0000000000000000000000000000000000000000..8f0af80e0da0de71599babef39d407eca5138621
--- /dev/null
+++ b/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;
diff --git a/applications/solvers/multiphase/interIsoFoam/correctPhi.H b/applications/solvers/multiphase/interIsoFoam/correctPhi.H
new file mode 100644
index 0000000000000000000000000000000000000000..9afcd58a66efaf0aa09bae9ecb5c39d095bf7099
--- /dev/null
+++ b/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"
diff --git a/applications/solvers/multiphase/interIsoFoam/createFields.H b/applications/solvers/multiphase/interIsoFoam/createFields.H
new file mode 100644
index 0000000000000000000000000000000000000000..8c3a05e3fab2c528971f94d24f3acbab990fb26a
--- /dev/null
+++ b/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"
diff --git a/applications/solvers/multiphase/interIsoFoam/createIsoAdvection.H b/applications/solvers/multiphase/interIsoFoam/createIsoAdvection.H
new file mode 100644
index 0000000000000000000000000000000000000000..889af7228f5b80fe69037560c5ee945763698575
--- /dev/null
+++ b/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;
+}
diff --git a/applications/solvers/multiphase/interIsoFoam/interIsoFoam.C b/applications/solvers/multiphase/interIsoFoam/interIsoFoam.C
new file mode 100644
index 0000000000000000000000000000000000000000..1514c5c0cb73516f47f81cdbd66871693a8c6587
--- /dev/null
+++ b/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;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/interIsoFoam/pEqn.H b/applications/solvers/multiphase/interIsoFoam/pEqn.H
new file mode 100644
index 0000000000000000000000000000000000000000..17c3ae9748fe754acc6dbb7971f9b19e15937b46
--- /dev/null
+++ b/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;
+ }
+}
diff --git a/applications/solvers/multiphase/interIsoFoam/setDeltaT.H b/applications/solvers/multiphase/interIsoFoam/setDeltaT.H
new file mode 100644
index 0000000000000000000000000000000000000000..9cc860c032eed2e9a46138518589b83f1e8aed19
--- /dev/null
+++ b/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;
+}
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/setAlphaField/Make/files b/applications/utilities/preProcessing/setAlphaField/Make/files
new file mode 100755
index 0000000000000000000000000000000000000000..23145059e2763a92c69ee3e5b45429d91adb1404
--- /dev/null
+++ b/applications/utilities/preProcessing/setAlphaField/Make/files
@@ -0,0 +1,3 @@
+setAlphaField.C
+
+EXE = $(FOAM_APPBIN)/setAlphaField
diff --git a/applications/utilities/preProcessing/setAlphaField/Make/options b/applications/utilities/preProcessing/setAlphaField/Make/options
new file mode 100755
index 0000000000000000000000000000000000000000..b4ad3d8bdff2f2475467775f28a6c62effa80575
--- /dev/null
+++ b/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
diff --git a/applications/utilities/preProcessing/setAlphaField/setAlphaField.C b/applications/utilities/preProcessing/setAlphaField/setAlphaField.C
new file mode 100644
index 0000000000000000000000000000000000000000..6fe8001e31d440e8a8d4e188c80308cd84d5be49
--- /dev/null
+++ b/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;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/setAlphaField/setAlphaFieldDict b/applications/utilities/preProcessing/setAlphaField/setAlphaFieldDict
new file mode 100644
index 0000000000000000000000000000000000000000..6bbb616b549e174c1360681c03405f0a628bcf8d
--- /dev/null
+++ b/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);
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files
index 1087c3e378540fe4988f956fed8b379fe1b177fe..ca9bf31ff6f9ab606cca79fa955659661e65db54 100644
--- a/src/finiteVolume/Make/files
+++ b/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
diff --git a/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoAdvection/isoAdvection.C b/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoAdvection/isoAdvection.C
new file mode 100644
index 0000000000000000000000000000000000000000..33f63c2b6cfae32d3ac4dd2d5bcc3219d28f0fed
--- /dev/null
+++ b/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoAdvection/isoAdvection.C
@@ -0,0 +1,1165 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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"
+#include "volFields.H"
+#include "interpolationCellPoint.H"
+#include "volPointInterpolation.H"
+#include "fvcSurfaceIntegrate.H"
+#include "fvcGrad.H"
+#include "upwind.H"
+#include "cellSet.H"
+#include "meshTools.H"
+#include "OBJstream.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+ defineTypeNameAndDebug(isoAdvection, 0);
+}
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::isoAdvection::isoAdvection
+(
+ volScalarField& alpha1,
+ const surfaceScalarField& phi,
+ const volVectorField& U
+)
+:
+ // General data
+ mesh_(alpha1.mesh()),
+ dict_(mesh_.solverDict(alpha1.name())),
+ alpha1_(alpha1),
+ alpha1In_(alpha1.ref()),
+ phi_(phi),
+ U_(U),
+ dVf_
+ (
+ IOobject
+ (
+ "dVf_",
+ mesh_.time().timeName(),
+ mesh_,
+ IOobject::NO_READ,
+ IOobject::NO_WRITE
+ ),
+ mesh_,
+ dimensionedScalar("zero", dimVol, 0)
+ ),
+ advectionTime_(0),
+
+ // Interpolation data
+ ap_(mesh_.nPoints()),
+
+ // Tolerances and solution controls
+ nAlphaBounds_(dict_.lookupOrDefault<label>("nAlphaBounds", 3)),
+ vof2IsoTol_(dict_.lookupOrDefault<scalar>("vof2IsoTol", 1e-8)),
+ surfCellTol_(dict_.lookupOrDefault<scalar>("surfCellTol", 1e-8)),
+ gradAlphaBasedNormal_
+ (
+ dict_.lookupOrDefault<bool>("gradAlphaNormal", false)
+ ),
+ writeIsoFacesToFile_
+ (
+ dict_.lookupOrDefault<bool>("writeIsoFaces", false)
+ ),
+
+ // Cell cutting data
+ surfCells_(label(0.2*mesh_.nCells())),
+ isoCutCell_(mesh_, ap_),
+ isoCutFace_(mesh_, ap_),
+ cellIsBounded_(mesh_.nCells(), false),
+ checkBounding_(mesh_.nCells(), false),
+ bsFaces_(label(0.2*(mesh_.nFaces() - mesh_.nInternalFaces()))),
+ bsx0_(bsFaces_.size()),
+ bsn0_(bsFaces_.size()),
+ bsUn0_(bsFaces_.size()),
+ bsf0_(bsFaces_.size()),
+
+ // Parallel run data
+ procPatchLabels_(mesh_.boundary().size()),
+ surfaceCellFacesOnProcPatches_(0)
+{
+ isoCutCell::debug = debug;
+ isoCutFace::debug = debug;
+
+ // Prepare lists used in parallel runs
+ if (Pstream::parRun())
+ {
+ // Force calculation of required demand driven data (else parallel
+ // communication may crash)
+ mesh_.cellCentres();
+ mesh_.cellVolumes();
+ mesh_.faceCentres();
+ mesh_.faceAreas();
+ mesh_.magSf();
+ mesh_.boundaryMesh().patchID();
+ mesh_.cellPoints();
+ mesh_.cellCells();
+ mesh_.cells();
+
+ // Get boundary mesh and resize the list for parallel comms
+ const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+
+ surfaceCellFacesOnProcPatches_.resize(patches.size());
+
+ // Append all processor patch labels to the list
+ forAll(patches, patchi)
+ {
+ if
+ (
+ isA<processorPolyPatch>(patches[patchi])
+ && patches[patchi].size() > 0
+ )
+ {
+ procPatchLabels_.append(patchi);
+ }
+ }
+ }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+void Foam::isoAdvection::timeIntegratedFlux()
+{
+ // Get time step
+ const scalar dt = mesh_.time().deltaTValue();
+
+ // Create object for interpolating velocity to isoface centres
+ interpolationCellPoint<vector> UInterp(U_);
+
+ // For each downwind face of each surface cell we "isoadvect" to find dVf
+ label nSurfaceCells = 0;
+
+ // Clear out the data for re-use and reset list containing information
+ // whether cells could possibly need bounding
+ clearIsoFaceData();
+ checkBounding_ = false;
+
+ // Get necessary references
+ const scalarField& phiIn = phi_.primitiveField();
+ const scalarField& magSfIn = mesh_.magSf().primitiveField();
+ scalarField& dVfIn = dVf_.primitiveFieldRef();
+
+ // Get necessary mesh data
+ const labelListList& cellPoints = mesh_.cellPoints();
+ const labelListList& cellCells = mesh_.cellCells();
+ const cellList& cellFaces = mesh_.cells();
+ const labelList& own = mesh_.faceOwner();
+ const labelList& nei = mesh_.faceNeighbour();
+ const vectorField& cellCentres = mesh_.cellCentres();
+ const pointField& points = mesh_.points();
+
+ // Storage for isoFace points. Only used if writeIsoFacesToFile_
+ DynamicList<List<point>> isoFacePts;
+
+ // Interpolating alpha1 cell centre values to mesh points (vertices)
+ ap_ = volPointInterpolation::New(mesh_).interpolate(alpha1_);
+
+ vectorField gradAlpha(mesh_.nPoints(), vector::zero);
+ if (gradAlphaBasedNormal_)
+ {
+ // Calculate gradient of alpha1 and interpolate to vertices
+ volVectorField gradA("gradA", fvc::grad(alpha1_));
+ gradAlpha = volPointInterpolation::New(mesh_).interpolate(gradA);
+ }
+
+ // Loop through cells
+ forAll(alpha1In_, celli)
+ {
+ if (isASurfaceCell(celli))
+ {
+ // This is a surface cell, increment counter, append and mark cell
+ nSurfaceCells++;
+ surfCells_.append(celli);
+ checkBounding_[celli] = true;
+
+ DebugInfo
+ << "\n------------ Cell " << celli << " with alpha1 = "
+ << alpha1In_[celli] << " and 1-alpha1 = "
+ << 1.0 - alpha1In_[celli] << " ------------"
+ << endl;
+
+ // Calculate isoFace centre x0, normal n0 at time t
+ label maxIter = 100; // NOTE: make it a debug switch
+
+ const labelList& cp = cellPoints[celli];
+ scalarField ap_org(cp.size(), 0);
+ if (gradAlphaBasedNormal_)
+ {
+ // Calculating smoothed alpha gradient in surface cell in order
+ // to use it as the isoface orientation.
+ vector smoothedGradA = vector::zero;
+ const point& cellCentre = cellCentres[celli];
+ scalar wSum = 0;
+ forAll(cp, pointI)
+ {
+ point vertex = points[cp[pointI]];
+ scalar w = 1.0/mag(vertex - cellCentre);
+ wSum += w;
+ smoothedGradA += w*gradAlpha[cp[pointI]];
+ }
+ smoothedGradA /= wSum;
+
+ // Temporarily overwrite the interpolated vertex alpha values in
+ // ap_ with the vertex-cell centre distance along smoothedGradA.
+ forAll(ap_org, vi)
+ {
+ ap_org[vi] = ap_[cp[vi]];
+ const point& vertex = points[cp[vi]];
+ ap_[cp[vi]] =
+ (
+ (vertex - cellCentre)
+ & (smoothedGradA/mag(smoothedGradA))
+ );
+ }
+ }
+
+ // Calculate cell status (-1: cell is fully below the isosurface, 0:
+ // cell is cut, 1: cell is fully above the isosurface)
+ label cellStatus = isoCutCell_.vofCutCell
+ (
+ celli,
+ alpha1In_[celli],
+ vof2IsoTol_,
+ maxIter
+ );
+
+ if (gradAlphaBasedNormal_)
+ {
+ // Restoring ap_ by putting the original values back into it.
+ forAll(ap_org, vi)
+ {
+ ap_[cp[vi]] = ap_org[vi];
+ }
+ }
+
+ // Cell is cut
+ if (cellStatus == 0)
+ {
+ const scalar f0 = isoCutCell_.isoValue();
+ const point& x0 = isoCutCell_.isoFaceCentre();
+ vector n0 = isoCutCell_.isoFaceArea();
+ n0 /= (mag(n0));
+
+ if (writeIsoFacesToFile_ && mesh_.time().writeTime())
+ {
+ isoFacePts.append(isoCutCell_.isoFacePoints());
+ }
+
+ // Get the speed of the isoface by interpolating velocity and
+ // dotting it with isoface normal
+ const scalar Un0 = UInterp.interpolate(x0, celli) & n0;
+
+ DebugInfo
+ << "calcIsoFace gives initial surface: \nx0 = " << x0
+ << ", \nn0 = " << n0 << ", \nf0 = " << f0 << ", \nUn0 = "
+ << Un0 << endl;
+
+ // Estimate time integrated flux through each downwind face
+ // Note: looping over all cell faces - in reduced-D, some of
+ // these faces will be on empty patches
+ const cell& celliFaces = cellFaces[celli];
+ forAll(celliFaces, fi)
+ {
+ const label facei = celliFaces[fi];
+
+ if (mesh_.isInternalFace(facei))
+ {
+ bool isDownwindFace = false;
+ label otherCell = -1;
+
+ if (celli == own[facei])
+ {
+ if (phiIn[facei] > 10*SMALL)
+ {
+ isDownwindFace = true;
+ }
+
+ otherCell = nei[facei];
+ }
+ else
+ {
+ if (phiIn[facei] < -10*SMALL)
+ {
+ isDownwindFace = true;
+ }
+
+ otherCell = own[facei];
+ }
+
+ if (isDownwindFace)
+ {
+ dVfIn[facei] = timeIntegratedFaceFlux
+ (
+ facei,
+ x0,
+ n0,
+ Un0,
+ f0,
+ dt,
+ phiIn[facei],
+ magSfIn[facei]
+ );
+ }
+
+ // We want to check bounding of neighbour cells to
+ // surface cells as well:
+ checkBounding_[otherCell] = true;
+
+ // Also check neighbours of neighbours.
+ // Note: consider making it a run time selectable
+ // extension level (easily done with recursion):
+ // 0 - only neighbours
+ // 1 - neighbours of neighbours
+ // 2 - ...
+ const labelList& nNeighbourCells = cellCells[otherCell];
+ forAll(nNeighbourCells, ni)
+ {
+ checkBounding_[nNeighbourCells[ni]] = true;
+ }
+ }
+ else
+ {
+ bsFaces_.append(facei);
+ bsx0_.append(x0);
+ bsn0_.append(n0);
+ bsUn0_.append(Un0);
+ bsf0_.append(f0);
+
+ // Note: we must not check if the face is on the
+ // processor patch here.
+ }
+ }
+ }
+ }
+ }
+
+ if (writeIsoFacesToFile_ && mesh_.time().writeTime())
+ {
+ writeIsoFaces(isoFacePts);
+ }
+
+ // Get references to boundary fields
+ const polyBoundaryMesh& boundaryMesh = mesh_.boundaryMesh();
+ const surfaceScalarField::Boundary& phib = phi_.boundaryField();
+ const surfaceScalarField::Boundary& magSfb = mesh_.magSf().boundaryField();
+ surfaceScalarField::Boundary& dVfb = dVf_.boundaryFieldRef();
+ const label nInternalFaces = mesh_.nInternalFaces();
+
+ // Loop through boundary surface faces
+ forAll(bsFaces_, i)
+ {
+ // Get boundary face index (in the global list)
+ const label facei = bsFaces_[i];
+ const label patchi = boundaryMesh.patchID()[facei - nInternalFaces];
+ const label start = boundaryMesh[patchi].start();
+
+ if (phib[patchi].size())
+ {
+ const label patchFacei = facei - start;
+ const scalar phiP = phib[patchi][patchFacei];
+
+ if (phiP > 10*SMALL)
+ {
+ const scalar magSf = magSfb[patchi][patchFacei];
+
+ dVfb[patchi][patchFacei] = timeIntegratedFaceFlux
+ (
+ facei,
+ bsx0_[i],
+ bsn0_[i],
+ bsUn0_[i],
+ bsf0_[i],
+ dt,
+ phiP,
+ magSf
+ );
+
+ // Check if the face is on processor patch and append it to
+ // the list if necessary
+ checkIfOnProcPatch(facei);
+ }
+ }
+ }
+
+ Info<< "Number of isoAdvector surface cells = "
+ << returnReduce(nSurfaceCells, sumOp<label>()) << endl;
+}
+
+
+Foam::scalar Foam::isoAdvection::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
+)
+{
+ // Treating rare cases where isoface normal is not calculated properly
+ if (mag(n0) < 0.5)
+ {
+ scalar alphaf = 0;
+ scalar waterInUpwindCell = 0;
+
+ if (phi > 10*SMALL || !mesh_.isInternalFace(facei))
+ {
+ const label upwindCell = mesh_.faceOwner()[facei];
+ alphaf = alpha1In_[upwindCell];
+ waterInUpwindCell = alphaf*mesh_.cellVolumes()[upwindCell];
+ }
+ else
+ {
+ const label upwindCell = mesh_.faceNeighbour()[facei];
+ alphaf = alpha1In_[upwindCell];
+ waterInUpwindCell = alphaf*mesh_.cellVolumes()[upwindCell];
+ }
+
+ if (debug)
+ {
+ WarningInFunction
+ << "mag(n0) = " << mag(n0)
+ << " so timeIntegratedFlux calculates dVf from upwind"
+ << " cell alpha value: " << alphaf << endl;
+ }
+
+ return min(alphaf*phi*dt, waterInUpwindCell);
+ }
+
+
+ // Find sorted list of times where the isoFace will arrive at face points
+ // given initial position x0 and velocity Un0*n0
+
+ // Get points for this face
+ const face& f = mesh_.faces()[facei];
+ const pointField fPts(f.points(mesh_.points()));
+ const label nPoints = fPts.size();
+
+ scalarField pTimes(fPts.size());
+ if (mag(Un0) > 10*SMALL) // Note: tolerances
+ {
+ // Here we estimate time of arrival to the face points from their normal
+ // distance to the initial surface and the surface normal velocity
+
+ pTimes = ((fPts - x0) & n0)/Un0;
+
+ scalar dVf = 0;
+
+ // Check if pTimes changes direction more than twice when looping face
+ label nShifts = 0;
+ forAll(pTimes, pi)
+ {
+ const label oldEdgeSign =
+ sign(pTimes[(pi + 1) % nPoints] - pTimes[pi]);
+ const label newEdgeSign =
+ sign(pTimes[(pi + 2) % nPoints] - pTimes[(pi + 1) % nPoints]);
+
+ if (newEdgeSign != oldEdgeSign)
+ {
+ nShifts++;
+ }
+ }
+
+ if (nShifts == 2)
+ {
+ dVf =
+ phi/magSf
+ *isoCutFace_.timeIntegratedArea(fPts, pTimes, dt, magSf, Un0);
+ }
+ else if (nShifts > 2)
+ {
+ // Triangle decompose the face
+ pointField fPts_tri(3);
+ scalarField pTimes_tri(3);
+ fPts_tri[0] = mesh_.faceCentres()[facei];
+ pTimes_tri[0] = ((fPts_tri[0] - x0) & n0)/Un0;
+ for (label pi = 0; pi < nPoints; pi++)
+ {
+ fPts_tri[1] = fPts[pi];
+ pTimes_tri[1] = pTimes[pi];
+ fPts_tri[2] = fPts[(pi + 1) % nPoints];
+ pTimes_tri[2] = pTimes[(pi + 1) % nPoints];
+ const scalar magSf_tri =
+ mag
+ (
+ 0.5
+ *(fPts_tri[2] - fPts_tri[0])
+ ^(fPts_tri[1] - fPts_tri[0])
+ );
+ const scalar phi_tri = phi*magSf_tri/magSf;
+ dVf +=
+ phi_tri
+ /magSf_tri
+ *isoCutFace_.timeIntegratedArea
+ (
+ fPts_tri,
+ pTimes_tri,
+ dt,
+ magSf_tri,
+ Un0
+ );
+ }
+ }
+ else
+ {
+ if (debug)
+ {
+ WarningInFunction
+ << "Warning: nShifts = " << nShifts << " on face " << facei
+ << " with pTimes = " << pTimes << " owned by cell "
+ << mesh_.faceOwner()[facei] << endl;
+ }
+ }
+
+ return dVf;
+ }
+ else
+ {
+ // Un0 is almost zero and isoFace is treated as stationary
+ isoCutFace_.calcSubFace(facei, f0);
+ const scalar alphaf = mag(isoCutFace_.subFaceArea()/magSf);
+
+ if (debug)
+ {
+ WarningInFunction
+ << "Un0 is almost zero (" << Un0
+ << ") - calculating dVf on face " << facei
+ << " using subFaceFraction giving alphaf = " << alphaf
+ << endl;
+ }
+
+ return phi*dt*alphaf;
+ }
+}
+
+
+void Foam::isoAdvection::setDownwindFaces
+(
+ const label celli,
+ DynamicLabelList& downwindFaces
+) const
+{
+ DebugInFunction << endl;
+
+ // Get necessary mesh data and cell information
+ const labelList& own = mesh_.faceOwner();
+ const cellList& cells = mesh_.cells();
+ const cell& c = cells[celli];
+
+ downwindFaces.clear();
+
+ // Check all faces of the cell
+ forAll(c, fi)
+ {
+ // Get face and corresponding flux
+ const label facei = c[fi];
+ const scalar phi = faceValue(phi_, facei);
+
+ if (own[facei] == celli)
+ {
+ if (phi > 10*SMALL)
+ {
+ downwindFaces.append(facei);
+ }
+ }
+ else if (phi < -10*SMALL)
+ {
+ downwindFaces.append(facei);
+ }
+ }
+
+ downwindFaces.shrink();
+}
+
+
+void Foam::isoAdvection::limitFluxes()
+{
+ DebugInFunction << endl;
+
+ // Get time step size
+ const scalar dt = mesh_.time().deltaT().value();
+
+// scalarField alphaNew = alpha1In_ - fvc::surfaceIntegrate(dVf_);
+ const scalar aTol = 1.0e-12; // Note: tolerances
+ const scalar maxAlphaMinus1 = 1; // max(alphaNew - 1);
+ const scalar minAlpha = -1; // min(alphaNew);
+ const label nUndershoots = 20; // sum(neg(alphaNew + aTol));
+ const label nOvershoots = 20; // sum(pos(alphaNew - 1 - aTol));
+ cellIsBounded_ = false;
+
+ // Loop number of bounding steps
+ for (label n = 0; n < nAlphaBounds_; n++)
+ {
+ Info<< "isoAdvection: bounding iteration " << n + 1 << endl;
+
+ if (maxAlphaMinus1 > aTol) // Note: tolerances
+ {
+ DebugInfo << "Bound from above... " << endl;
+
+// scalarField& dVfcorrected = dVf_.primitiveFieldRef();
+
+ surfaceScalarField dVfcorrected("dVfcorrected", dVf_);
+ DynamicList<label> correctedFaces(3*nOvershoots);
+ boundFromAbove(alpha1In_, dVfcorrected, correctedFaces);
+
+ forAll(correctedFaces, fi)
+ {
+ label facei = correctedFaces[fi];
+
+ // Change to treat boundaries consistently
+ setFaceValue(dVf_, facei, faceValue(dVfcorrected, facei));
+ }
+
+ syncProcPatches(dVf_, phi_);
+ }
+
+ if (minAlpha < -aTol) // Note: tolerances
+ {
+ DebugInfo << "Bound from below... " << endl;
+
+ scalarField alpha2(1.0 - alpha1In_);
+ surfaceScalarField dVfcorrected
+ (
+ "dVfcorrected",
+ phi_*dimensionedScalar("dt", dimTime, dt) - dVf_
+ );
+// dVfcorrected -= dVf_; // phi_ and dVf_ have same sign and dVf_ is
+ // the portion of phi_*dt that is water.
+ // If phi_ > 0 then dVf_ > 0 and mag(phi_*dt-dVf_) < mag(phi_*dt) as
+ // it should.
+ // If phi_ < 0 then dVf_ < 0 and mag(phi_*dt-dVf_) < mag(phi_*dt) as
+ // it should.
+ DynamicList<label> correctedFaces(3*nUndershoots);
+ boundFromAbove(alpha2, dVfcorrected, correctedFaces);
+ forAll(correctedFaces, fi)
+ {
+ label facei = correctedFaces[fi];
+
+ // Change to treat boundaries consistently
+ scalar phi = faceValue(phi_, facei);
+ scalar dVcorr = faceValue(dVfcorrected, facei);
+ setFaceValue(dVf_, facei, phi*dt - dVcorr);
+ }
+
+ syncProcPatches(dVf_, phi_);
+ }
+
+ if (debug)
+ {
+ // Check if still unbounded
+ scalarField alphaNew(alpha1In_ - fvc::surfaceIntegrate(dVf_)());
+ label maxAlphaMinus1 = max(alphaNew - 1);
+ scalar minAlpha = min(alphaNew);
+ label nUndershoots = sum(neg(alphaNew + aTol));
+ label nOvershoots = sum(pos(alphaNew - 1 - aTol));
+ Info<< "After bounding number " << n + 1 << " of time "
+ << mesh_.time().value() << ":" << endl;
+ Info<< "nOvershoots = " << nOvershoots << " with max(alphaNew-1) = "
+ << maxAlphaMinus1 << " and nUndershoots = " << nUndershoots
+ << " with min(alphaNew) = " << minAlpha << endl;
+ }
+ }
+}
+
+
+void Foam::isoAdvection::boundFromAbove
+(
+ const scalarField& alpha1,
+ surfaceScalarField& dVf,
+ DynamicList<label>& correctedFaces
+)
+{
+ DebugInFunction << endl;
+
+ correctedFaces.clear();
+ scalar aTol = 10*SMALL; // Note: tolerances
+
+ const scalarField& meshV = mesh_.cellVolumes();
+ const scalar dt = mesh_.time().deltaTValue();
+
+ DynamicList<label> downwindFaces(10);
+ DynamicList<label> facesToPassFluidThrough(downwindFaces.size());
+ DynamicList<scalar> dVfmax(downwindFaces.size());
+ DynamicList<scalar> phi(downwindFaces.size());
+
+ // Loop through alpha cell centred field
+ forAll(alpha1, celli)
+ {
+ if (checkBounding_[celli])
+ {
+ const scalar Vi = meshV[celli];
+ scalar alpha1New = alpha1[celli] - netFlux(dVf, celli)/Vi;
+ scalar alphaOvershoot = alpha1New - 1.0;
+ scalar fluidToPassOn = alphaOvershoot*Vi;
+ label nFacesToPassFluidThrough = 1;
+
+ bool firstLoop = true;
+
+ // First try to pass surplus fluid on to neighbour cells that are
+ // not filled and to which dVf < phi*dt
+ while (alphaOvershoot > aTol && nFacesToPassFluidThrough > 0)
+ {
+ DebugInfo
+ << "\n\nBounding cell " << celli
+ << " with alpha overshooting " << alphaOvershoot
+ << endl;
+
+ facesToPassFluidThrough.clear();
+ dVfmax.clear();
+ phi.clear();
+
+ cellIsBounded_[celli] = true;
+
+ // Find potential neighbour cells to pass surplus phase to
+ setDownwindFaces(celli, downwindFaces);
+
+ scalar dVftot = 0;
+ nFacesToPassFluidThrough = 0;
+
+ forAll(downwindFaces, fi)
+ {
+ const label facei = downwindFaces[fi];
+ const scalar phif = faceValue(phi_, facei);
+ const scalar dVff = faceValue(dVf, facei);
+ const scalar maxExtraFaceFluidTrans = mag(phif*dt - dVff);
+
+ // dVf has same sign as phi and so if phi>0 we have
+ // mag(phi_[facei]*dt) - mag(dVf[facei]) = phi_[facei]*dt
+ // - dVf[facei]
+ // If phi < 0 we have mag(phi_[facei]*dt) -
+ // mag(dVf[facei]) = -phi_[facei]*dt - (-dVf[facei]) > 0
+ // since mag(dVf) < phi*dt
+ DebugInfo
+ << "downwindFace " << facei
+ << " has maxExtraFaceFluidTrans = "
+ << maxExtraFaceFluidTrans << endl;
+
+ if (maxExtraFaceFluidTrans/Vi > aTol)
+ {
+// if (maxExtraFaceFluidTrans/Vi > aTol &&
+// mag(dVfIn[facei])/Vi > aTol) //Last condition may be
+// important because without this we will flux through uncut
+// downwind faces
+ facesToPassFluidThrough.append(facei);
+ phi.append(phif);
+ dVfmax.append(maxExtraFaceFluidTrans);
+ dVftot += mag(phif*dt);
+ }
+ }
+
+ DebugInfo
+ << "\nfacesToPassFluidThrough: "
+ << facesToPassFluidThrough << ", dVftot = "
+ << dVftot << " m3 corresponding to dalpha = "
+ << dVftot/Vi << endl;
+
+ forAll(facesToPassFluidThrough, fi)
+ {
+ const label facei = facesToPassFluidThrough[fi];
+ scalar fluidToPassThroughFace =
+ fluidToPassOn*mag(phi[fi]*dt)/dVftot;
+
+ nFacesToPassFluidThrough +=
+ pos(dVfmax[fi] - fluidToPassThroughFace);
+
+ fluidToPassThroughFace =
+ min(fluidToPassThroughFace, dVfmax[fi]);
+
+ scalar dVff = faceValue(dVf, facei);
+ dVff += sign(phi[fi])*fluidToPassThroughFace;
+ setFaceValue(dVf, facei, dVff);
+
+ if (firstLoop)
+ {
+ checkIfOnProcPatch(facei);
+ correctedFaces.append(facei);
+ }
+ }
+
+ firstLoop = false;
+ alpha1New = alpha1[celli] - netFlux(dVf, celli)/Vi;
+ alphaOvershoot = alpha1New - 1.0;
+ fluidToPassOn = alphaOvershoot*Vi;
+
+ DebugInfo
+ << "\nNew alpha for cell " << celli << ": "
+ << alpha1New << endl;
+ }
+ }
+ }
+
+ DebugInfo << "correctedFaces = " << correctedFaces << endl;
+}
+
+
+Foam::scalar Foam::isoAdvection::netFlux
+(
+ const surfaceScalarField& dVf,
+ const label celli
+) const
+{
+ scalar dV = 0;
+
+ // Get face indices
+ const cell& c = mesh_.cells()[celli];
+
+ // Get mesh data
+ const labelList& own = mesh_.faceOwner();
+
+ forAll(c, fi)
+ {
+ const label facei = c[fi];
+ const scalar dVff = faceValue(dVf, facei);
+
+ if (own[facei] == celli)
+ {
+ dV += dVff;
+ }
+ else
+ {
+ dV -= dVff;
+ }
+ }
+
+ return dV;
+}
+
+
+void Foam::isoAdvection::syncProcPatches
+(
+ surfaceScalarField& dVf,
+ const surfaceScalarField& phi
+)
+{
+ const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+
+ if (Pstream::parRun())
+ {
+ PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking);
+
+ // Send
+ forAll(procPatchLabels_, i)
+ {
+ const label patchi = procPatchLabels_[i];
+
+ const processorPolyPatch& procPatch =
+ refCast<const processorPolyPatch>(patches[patchi]);
+
+ UOPstream toNbr(procPatch.neighbProcNo(), pBufs);
+ const scalarField& pFlux = dVf.boundaryField()[patchi];
+
+ const List<label>& surfCellFacesOnProcPatch =
+ surfaceCellFacesOnProcPatches_[patchi];
+
+ const UIndirectList<scalar> dVfPatch
+ (
+ pFlux,
+ surfCellFacesOnProcPatch
+ );
+
+ toNbr << surfCellFacesOnProcPatch << dVfPatch;
+ }
+
+ pBufs.finishedSends();
+
+
+ // Receive and combine
+ forAll(procPatchLabels_, patchLabeli)
+ {
+ const label patchi = procPatchLabels_[patchLabeli];
+
+ const processorPolyPatch& procPatch =
+ refCast<const processorPolyPatch>(patches[patchi]);
+
+ UIPstream fromNeighb(procPatch.neighbProcNo(), pBufs);
+ List<label> faceIDs;
+ List<scalar> nbrdVfs;
+
+ fromNeighb >> faceIDs >> nbrdVfs;
+
+ if (debug)
+ {
+ Pout<< "Received at time = " << mesh_.time().value()
+ << ": surfCellFacesOnProcPatch = " << faceIDs << nl
+ << "Received at time = " << mesh_.time().value()
+ << ": dVfPatch = " << nbrdVfs << endl;
+ }
+
+ // Combine fluxes
+ scalarField& localFlux = dVf.boundaryFieldRef()[patchi];
+
+ forAll(faceIDs, i)
+ {
+ const label facei = faceIDs[i];
+ localFlux[facei] = - nbrdVfs[i];
+ if (debug && mag(localFlux[facei] + nbrdVfs[i]) > 10*SMALL)
+ {
+ Pout<< "localFlux[facei] = " << localFlux[facei]
+ << " and nbrdVfs[i] = " << nbrdVfs[i]
+ << " for facei = " << facei << endl;
+ }
+ }
+ }
+
+ if (debug)
+ {
+ // Write out results for checking
+ forAll(procPatchLabels_, patchLabeli)
+ {
+ const label patchi = procPatchLabels_[patchLabeli];
+ const scalarField& localFlux = dVf.boundaryField()[patchi];
+ Pout<< "time = " << mesh_.time().value() << ": localFlux = "
+ << localFlux << endl;
+ }
+ }
+
+ // Reinitialising list used for minimal parallel communication
+ forAll(surfaceCellFacesOnProcPatches_, patchi)
+ {
+ surfaceCellFacesOnProcPatches_[patchi].clear();
+ }
+ }
+}
+
+
+void Foam::isoAdvection::checkIfOnProcPatch(const label facei)
+{
+ if (!mesh_.isInternalFace(facei))
+ {
+ const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
+ const label patchi = pbm.patchID()[facei - mesh_.nInternalFaces()];
+
+ if (isA<processorPolyPatch>(pbm[patchi]) && pbm[patchi].size())
+ {
+ const label patchFacei = pbm[patchi].whichFace(facei);
+ surfaceCellFacesOnProcPatches_[patchi].append(patchFacei);
+ }
+ }
+}
+
+
+void Foam::isoAdvection::advect()
+{
+ DebugInFunction << endl;
+
+ scalar advectionStartTime = mesh_.time().elapsedCpuTime();
+
+ // Initialising dVf with upwind values
+ // i.e. phi[facei]*alpha1[upwindCell[facei]]*dt
+ dVf_ = upwind<scalar>(mesh_, phi_).flux(alpha1_)*mesh_.time().deltaT();
+
+ // Do the isoAdvection on surface cells
+ timeIntegratedFlux();
+
+ // Synchronize processor patches
+ syncProcPatches(dVf_, phi_);
+
+ // Adjust dVf for unbounded cells
+ limitFluxes();
+
+ // Advect the free surface
+ alpha1_ -= fvc::surfaceIntegrate(dVf_);
+ alpha1_.correctBoundaryConditions();
+
+ // Apply non-conservative bounding mechanisms (clipping and snapping)
+ // Note: We should be able to write out alpha before this is done!
+ applyBruteForceBounding();
+
+ // Write surface cell set and bound cell set if required by user
+ writeSurfaceCells();
+ writeBoundedCells();
+
+ advectionTime_ += (mesh_.time().elapsedCpuTime() - advectionStartTime);
+}
+
+
+void Foam::isoAdvection::applyBruteForceBounding()
+{
+ bool alpha1Changed = false;
+
+ scalar snapAlphaTol = dict_.lookupOrDefault<scalar>("snapTol", 0);
+ if (snapAlphaTol > 0)
+ {
+ alpha1_ =
+ alpha1_
+ *pos(alpha1_ - snapAlphaTol)
+ *neg(alpha1_ - (1.0 - snapAlphaTol))
+ + pos(alpha1_ - (1.0 - snapAlphaTol));
+
+ alpha1Changed = true;
+ }
+
+ bool clip = dict_.lookupOrDefault<bool>("clip", true);
+ if (clip)
+ {
+ alpha1_ = min(1.0, max(0.0, alpha1_));
+ alpha1Changed = true;
+ }
+
+ if (alpha1Changed)
+ {
+ alpha1_.correctBoundaryConditions();
+ }
+}
+
+
+void Foam::isoAdvection::writeSurfaceCells() const
+{
+ if (dict_.lookupOrDefault<bool>("writeSurfCells", false))
+ {
+ cellSet cSet
+ (
+ IOobject
+ (
+ "surfCells",
+ mesh_.time().timeName(),
+ mesh_,
+ IOobject::NO_READ
+ )
+ );
+
+ forAll(surfCells_, i)
+ {
+ cSet.insert(surfCells_[i]);
+ }
+
+ cSet.write();
+ }
+}
+
+
+void Foam::isoAdvection::writeBoundedCells() const
+{
+ if (dict_.lookupOrDefault<bool>("writeBoundedCells", false))
+ {
+ cellSet cSet
+ (
+ IOobject
+ (
+ "boundedCells",
+ mesh_.time().timeName(),
+ mesh_,
+ IOobject::NO_READ
+ )
+ );
+
+ forAll(cellIsBounded_, i)
+ {
+ if (cellIsBounded_[i])
+ {
+ cSet.insert(i);
+ }
+ }
+
+ cSet.write();
+ }
+}
+
+
+void Foam::isoAdvection::writeIsoFaces
+(
+ const DynamicList<List<point>>& faces
+) const
+{
+ // Writing isofaces to obj file for inspection, e.g. in paraview
+ const fileName dirName
+ (
+ Pstream::parRun() ?
+ mesh_.time().path()/".."/"isoFaces"
+ : mesh_.time().path()/"isoFaces"
+ );
+ const string fName
+ (
+ "isoFaces_" + Foam::name("%012d", mesh_.time().timeIndex())
+ );
+
+ if (Pstream::parRun())
+ {
+ // Collect points from all the processors
+ List<DynamicList<List<point>>> allProcFaces(Pstream::nProcs());
+ allProcFaces[Pstream::myProcNo()] = faces;
+ Pstream::gatherList(allProcFaces);
+
+ if (Pstream::master())
+ {
+ mkDir(dirName);
+ OBJstream os(dirName/fName + ".obj");
+ Info<< nl << "isoAdvection: writing iso faces to file: "
+ << os.name() << nl << endl;
+
+ face f;
+ forAll(allProcFaces, proci)
+ {
+ const DynamicList<List<point>>& procFacePts =
+ allProcFaces[proci];
+
+ forAll(procFacePts, i)
+ {
+ const List<point>& facePts = procFacePts[i];
+
+ if (facePts.size() != f.size())
+ {
+ f = face(identity(facePts.size()));
+ }
+
+ os.write(f, facePts, false);
+ }
+ }
+ }
+ }
+ else
+ {
+ mkDir(dirName);
+ OBJstream os(dirName/fName + ".obj");
+ Info<< nl << "isoAdvection: writing iso faces to file: "
+ << os.name() << nl << endl;
+
+ face f;
+ forAll(faces, i)
+ {
+ const List<point>& facePts = faces[i];
+
+ if (facePts.size() != f.size())
+ {
+ f = face(identity(facePts.size()));
+ }
+
+ os.write(f, facePts, false);
+ }
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoAdvection/isoAdvection.H b/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoAdvection/isoAdvection.H
new file mode 100644
index 0000000000000000000000000000000000000000..13899054fe872dd0e730283cf9dc462fe1f33a40
--- /dev/null
+++ b/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
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoAdvection/isoAdvectionTemplates.C b/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoAdvection/isoAdvectionTemplates.C
new file mode 100644
index 0000000000000000000000000000000000000000..baec97c65409d0b24dc3c6844ebeb4e73875814a
--- /dev/null
+++ b/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;
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutCell/isoCutCell.C b/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutCell/isoCutCell.C
new file mode 100644
index 0000000000000000000000000000000000000000..e4309ef48a06d1990bb8e31b5b8c67a5629d745f
--- /dev/null
+++ b/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];
+ }
+ }
+ }
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutCell/isoCutCell.H b/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutCell/isoCutCell.H
new file mode 100644
index 0000000000000000000000000000000000000000..171a4715a4b9d6a22598fa787c82193b268b84d7
--- /dev/null
+++ b/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
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutFace/isoCutFace.C b/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutFace/isoCutFace.C
new file mode 100644
index 0000000000000000000000000000000000000000..8eedcc0c81f11e46cec344ad814c665c2aa3dced
--- /dev/null
+++ b/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;
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutFace/isoCutFace.H b/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutFace/isoCutFace.H
new file mode 100644
index 0000000000000000000000000000000000000000..6a5e19013c36c3cd5e77f8e77dd28be5a8f77457
--- /dev/null
+++ b/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
+
+// ************************************************************************* //
diff --git a/src/functionObjects/field/Make/files b/src/functionObjects/field/Make/files
index 14959cbfa56448faca26c228c702cdc6e275db05..27a5790b3be84132689653badae9dfb88b790aac 100644
--- a/src/functionObjects/field/Make/files
+++ b/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
diff --git a/src/functionObjects/field/setFlow/setFlow.C b/src/functionObjects/field/setFlow/setFlow.C
new file mode 100644
index 0000000000000000000000000000000000000000..0d61bbff3f450b4f96c584bf8ebd9346397bd68e
--- /dev/null
+++ b/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;
+}
+
+
+// ************************************************************************* //
diff --git a/src/functionObjects/field/setFlow/setFlow.H b/src/functionObjects/field/setFlow/setFlow.H
new file mode 100644
index 0000000000000000000000000000000000000000..26425d8092781a2b526ab94b49de0bb05d1142bc
--- /dev/null
+++ b/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
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/damBreak/0.orig/U b/tutorials/multiphase/interIsoFoam/damBreak/0.orig/U
new file mode 100644
index 0000000000000000000000000000000000000000..5b146420656b9de88de6fa2f0ce7b1e30ff14228
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/damBreak/0.orig/U
@@ -0,0 +1,51 @@
+/*--------------------------------*- 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 volVectorField;
+ location "0";
+ object U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions [0 1 -1 0 0 0 0];
+
+internalField uniform (0 0 0);
+
+boundaryField
+{
+ leftWall
+ {
+ type fixedValue;
+ value uniform (0 0 0);
+ }
+ rightWall
+ {
+ type fixedValue;
+ value uniform (0 0 0);
+ }
+ lowerWall
+ {
+ type fixedValue;
+ value uniform (0 0 0);
+ }
+ atmosphere
+ {
+ type pressureInletOutletVelocity;
+ value uniform (0 0 0);
+ }
+ defaultFaces
+ {
+ type empty;
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/damBreak/0.orig/alpha.water b/tutorials/multiphase/interIsoFoam/damBreak/0.orig/alpha.water
new file mode 100644
index 0000000000000000000000000000000000000000..c0a881b0282d6910fa722a9e0002913bac5e76aa
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/damBreak/0.orig/alpha.water
@@ -0,0 +1,51 @@
+/*--------------------------------*- 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 volScalarField;
+ object alpha.water;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions [0 0 0 0 0 0 0];
+
+internalField uniform 0;
+
+boundaryField
+{
+ leftWall
+ {
+ type zeroGradient;
+ }
+
+ rightWall
+ {
+ type zeroGradient;
+ }
+
+ lowerWall
+ {
+ type zeroGradient;
+ }
+
+ atmosphere
+ {
+ type inletOutlet;
+ inletValue uniform 0;
+ value uniform 0;
+ }
+
+ defaultFaces
+ {
+ type empty;
+ }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/damBreak/0.orig/p_rgh b/tutorials/multiphase/interIsoFoam/damBreak/0.orig/p_rgh
new file mode 100644
index 0000000000000000000000000000000000000000..55339d401e4f3b2abbfa166e581e8d06d07c3b86
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/damBreak/0.orig/p_rgh
@@ -0,0 +1,59 @@
+/*--------------------------------*- 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 volScalarField;
+ object p_rgh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions [1 -1 -2 0 0 0 0];
+
+internalField uniform 0;
+
+boundaryField
+{
+ leftWall
+ {
+ type fixedFluxPressure;
+ value uniform 0;
+ }
+
+ rightWall
+ {
+ type fixedFluxPressure;
+ value uniform 0;
+ }
+
+ lowerWall
+ {
+ type fixedFluxPressure;
+ value uniform 0;
+ }
+
+ atmosphere
+ {
+ type totalPressure;
+ p0 uniform 0;
+ U U;
+ phi phi;
+ rho rho;
+ psi none;
+ gamma 1;
+ value uniform 0;
+ }
+
+ defaultFaces
+ {
+ type empty;
+ }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/damBreak/Allclean b/tutorials/multiphase/interIsoFoam/damBreak/Allclean
new file mode 100755
index 0000000000000000000000000000000000000000..513c366d4e453d40bf10eae880a1071f3e8bcc48
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/damBreak/Allclean
@@ -0,0 +1,9 @@
+#!/bin/sh
+cd ${0%/*} || exit 1 # Run from this directory
+
+# Source tutorial clean functions
+. $WM_PROJECT_DIR/bin/tools/CleanFunctions
+
+cleanCase
+
+\rm -rf 0 isoFaces
diff --git a/tutorials/multiphase/interIsoFoam/damBreak/Allrun b/tutorials/multiphase/interIsoFoam/damBreak/Allrun
new file mode 100755
index 0000000000000000000000000000000000000000..73a4aa44fc904dd62d6bb71affb53b394bce85c7
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/damBreak/Allrun
@@ -0,0 +1,12 @@
+#!/bin/sh
+cd ${0%/*} || exit 1 # Run from this directory
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+restore0Dir
+
+runApplication blockMesh
+runApplication setFields
+
+runApplication $(getApplication)
diff --git a/tutorials/multiphase/interIsoFoam/damBreak/Allrun-parallel b/tutorials/multiphase/interIsoFoam/damBreak/Allrun-parallel
new file mode 100755
index 0000000000000000000000000000000000000000..729fedf62cb4974d370296dc161b6bc756f728a8
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/damBreak/Allrun-parallel
@@ -0,0 +1,13 @@
+#!/bin/sh
+cd ${0%/*} || exit 1 # Run from this directory
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+restore0Dir
+
+runApplication blockMesh
+runApplication decomposePar
+runParallel setFields
+
+runParallel $(getApplication)
diff --git a/tutorials/multiphase/interIsoFoam/damBreak/constant/dynamicMeshDict b/tutorials/multiphase/interIsoFoam/damBreak/constant/dynamicMeshDict
new file mode 100644
index 0000000000000000000000000000000000000000..0d234e3050053c7590aa9b6089f7c448c2932b97
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/damBreak/constant/dynamicMeshDict
@@ -0,0 +1,21 @@
+/*--------------------------------*- 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 "constant";
+ object dynamicMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dynamicFvMesh staticFvMesh;
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/damBreak/constant/g b/tutorials/multiphase/interIsoFoam/damBreak/constant/g
new file mode 100644
index 0000000000000000000000000000000000000000..efd2b0b9ea08489348efd846c2fc22a831e71e7f
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/damBreak/constant/g
@@ -0,0 +1,22 @@
+/*--------------------------------*- 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 uniformDimensionedVectorField;
+ location "constant";
+ object g;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions [0 1 -2 0 0 0 0];
+value ( 0 -9.81 0 );
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/damBreak/constant/transportProperties b/tutorials/multiphase/interIsoFoam/damBreak/constant/transportProperties
new file mode 100644
index 0000000000000000000000000000000000000000..5bed7c36b6175c13b98c1b485ea678bf74c52c2c
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/damBreak/constant/transportProperties
@@ -0,0 +1,38 @@
+/*--------------------------------*- 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 "constant";
+ object transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+phases (water air);
+
+water
+{
+ transportModel Newtonian;
+ nu [0 2 -1 0 0 0 0] 1e-06;
+ rho [1 -3 0 0 0 0 0] 1000;
+}
+
+air
+{
+ transportModel Newtonian;
+ nu [0 2 -1 0 0 0 0] 1.48e-05;
+ rho [1 -3 0 0 0 0 0] 1;
+}
+
+
+sigma sigma [ 1 0 -2 0 0 0 0 ] 0.07;
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/damBreak/constant/turbulenceProperties b/tutorials/multiphase/interIsoFoam/damBreak/constant/turbulenceProperties
new file mode 100644
index 0000000000000000000000000000000000000000..5eec04267266e7fd15e7701a875d683d5e658cd9
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/damBreak/constant/turbulenceProperties
@@ -0,0 +1,21 @@
+/*--------------------------------*- 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 "constant";
+ object turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType laminar;
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/damBreak/system/blockMeshDict b/tutorials/multiphase/interIsoFoam/damBreak/system/blockMeshDict
new file mode 100644
index 0000000000000000000000000000000000000000..9570d8dd9c064c3b00c28335ce57fa40360bb23a
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/damBreak/system/blockMeshDict
@@ -0,0 +1,108 @@
+/*--------------------------------*- 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;
+ object blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 0.146;
+
+vertices
+(
+ (0 0 0)
+ (2 0 0)
+ (2.16438 0 0)
+ (4 0 0)
+ (0 0.32876 0)
+ (2 0.32876 0)
+ (2.16438 0.32876 0)
+ (4 0.32876 0)
+ (0 4 0)
+ (2 4 0)
+ (2.16438 4 0)
+ (4 4 0)
+ (0 0 0.1)
+ (2 0 0.1)
+ (2.16438 0 0.1)
+ (4 0 0.1)
+ (0 0.32876 0.1)
+ (2 0.32876 0.1)
+ (2.16438 0.32876 0.1)
+ (4 0.32876 0.1)
+ (0 4 0.1)
+ (2 4 0.1)
+ (2.16438 4 0.1)
+ (4 4 0.1)
+);
+
+blocks
+(
+ hex (0 1 5 4 12 13 17 16) (46 16 1) simpleGrading (1 1 1)
+ hex (2 3 7 6 14 15 19 18) (38 16 1) simpleGrading (1 1 1)
+ hex (4 5 9 8 16 17 21 20) (46 84 1) simpleGrading (1 1 1)
+ hex (5 6 10 9 17 18 22 21) (8 84 1) simpleGrading (1 1 1)
+ hex (6 7 11 10 18 19 23 22) (38 84 1) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+ leftWall
+ {
+ type wall;
+ faces
+ (
+ (0 12 16 4)
+ (4 16 20 8)
+ );
+ }
+ rightWall
+ {
+ type wall;
+ faces
+ (
+ (7 19 15 3)
+ (11 23 19 7)
+ );
+ }
+ lowerWall
+ {
+ type wall;
+ faces
+ (
+ (0 1 13 12)
+ (1 5 17 13)
+ (5 6 18 17)
+ (2 14 18 6)
+ (2 3 15 14)
+ );
+ }
+ atmosphere
+ {
+ type patch;
+ faces
+ (
+ (8 20 21 9)
+ (9 21 22 10)
+ (10 22 23 11)
+ );
+ }
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/damBreak/system/controlDict b/tutorials/multiphase/interIsoFoam/damBreak/system/controlDict
new file mode 100644
index 0000000000000000000000000000000000000000..4ac22e31c5494405190f4b9107258a750b926ed6
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/damBreak/system/controlDict
@@ -0,0 +1,56 @@
+/*--------------------------------*- 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 controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application interIsoFoam;
+
+startFrom latestTime;
+
+startTime 0;
+
+stopAt endTime;
+
+endTime 2;
+
+deltaT 0.001;
+
+writeControl adjustableRunTime;
+
+writeInterval 0.02;
+
+purgeWrite 0;
+
+writeFormat ascii;
+
+writePrecision 6;
+
+writeCompression uncompressed;
+
+timeFormat general;
+
+timePrecision 6;
+
+runTimeModifiable yes;
+
+adjustTimeStep yes;
+
+maxCo 10;
+maxAlphaCo 0.5;
+
+maxDeltaT 1;
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/damBreak/system/decomposeParDict b/tutorials/multiphase/interIsoFoam/damBreak/system/decomposeParDict
new file mode 100644
index 0000000000000000000000000000000000000000..b2904624d351c464ff314a3e2e57f6dfde1176b9
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/damBreak/system/decomposeParDict
@@ -0,0 +1,45 @@
+/*--------------------------------*- 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 decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains 4;
+
+method simple;
+
+simpleCoeffs
+{
+ n ( 2 2 1 );
+ delta 0.001;
+}
+
+hierarchicalCoeffs
+{
+ n ( 1 1 1 );
+ delta 0.001;
+ order xyz;
+}
+
+manualCoeffs
+{
+ dataFile "";
+}
+
+distributed no;
+
+roots ( );
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/damBreak/system/fvSchemes b/tutorials/multiphase/interIsoFoam/damBreak/system/fvSchemes
new file mode 100644
index 0000000000000000000000000000000000000000..dad2988ef4ba45b5b1b56ffd7e3ffcfcd0f9d557
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/damBreak/system/fvSchemes
@@ -0,0 +1,60 @@
+/*--------------------------------*- 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 fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+ default Euler;
+}
+
+gradSchemes
+{
+ default Gauss linear;
+}
+
+divSchemes
+{
+ div(rhoPhi,U) Gauss limitedLinearV 1;
+ div(phi,alpha) Gauss vanLeer;
+ div(phirb,alpha) Gauss interfaceCompression;
+ div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
+}
+
+laplacianSchemes
+{
+ default Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+ default linear;
+}
+
+snGradSchemes
+{
+ default corrected;
+}
+
+fluxRequired
+{
+ default no;
+ p_rgh;
+ pcorr;
+ alpha.water;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/damBreak/system/fvSolution b/tutorials/multiphase/interIsoFoam/damBreak/system/fvSolution
new file mode 100644
index 0000000000000000000000000000000000000000..0a189fa124167b12162d6724c99922ee62ecfb0a
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/damBreak/system/fvSolution
@@ -0,0 +1,76 @@
+/*--------------------------------*- 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;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+ "alpha.water.*"
+ {
+ isoFaceTol 1e-6;
+ surfCellTol 1e-6;
+ nAlphaBounds 3;
+ snapTol 1e-12;
+ clip true;
+
+ nAlphaSubCycles 1;
+ cAlpha 1;
+ }
+
+ "pcorr.*"
+ {
+ solver PCG;
+ preconditioner DIC;
+ tolerance 1e-10;
+ relTol 0;
+ }
+
+ p_rgh
+ {
+ solver GAMG;
+ smoother DICGaussSeidel;
+ tolerance 1e-07;
+ relTol 0.05;
+ }
+
+ p_rghFinal
+ {
+ $p_rgh;
+ tolerance 1e-07;
+ relTol 0;
+ }
+
+ U
+ {
+ solver PBiCGStab;
+ preconditioner DILU;
+ tolerance 1e-06;
+ relTol 0;
+ }
+}
+
+
+PIMPLE
+{
+ momentumPredictor no;
+ nCorrectors 3;
+ nOuterCorrectors 1;
+ nNonOrthogonalCorrectors 0;
+ pRefCell 0;
+ pRefValue 0;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/damBreak/system/setFieldsDict b/tutorials/multiphase/interIsoFoam/damBreak/system/setFieldsDict
new file mode 100644
index 0000000000000000000000000000000000000000..0605f3124706bc60c1ac1fc06ec9be3ab11f590c
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/damBreak/system/setFieldsDict
@@ -0,0 +1,36 @@
+/*--------------------------------*- 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 setFieldsDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+defaultFieldValues
+(
+ volScalarFieldValue alpha.water 0
+);
+
+regions
+(
+ boxToCell
+ {
+ box (0 0 -1) (0.1461 0.292 1);
+ fieldValues
+ (
+ volScalarFieldValue alpha.water 1
+ );
+ }
+);
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/discInConstantFlow/0.orig/U b/tutorials/multiphase/interIsoFoam/discInConstantFlow/0.orig/U
new file mode 100644
index 0000000000000000000000000000000000000000..c3c591f385b59e0ee6d5bc69900b304dd20bcc9c
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/discInConstantFlow/0.orig/U
@@ -0,0 +1,39 @@
+/*--------------------------------*- 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 volVectorField;
+ location "0";
+ object U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions [0 1 -1 0 0 0 0];
+
+internalField uniform (1 0 0.5);
+
+boundaryField
+{
+ rim
+ {
+ type zeroGradient;
+ }
+ front
+ {
+ type empty;
+ }
+ back
+ {
+ type empty;
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/discInConstantFlow/0.orig/alpha.water b/tutorials/multiphase/interIsoFoam/discInConstantFlow/0.orig/alpha.water
new file mode 100644
index 0000000000000000000000000000000000000000..5f274cc29edf48ca0fcb95482dae6269a5bff834
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/discInConstantFlow/0.orig/alpha.water
@@ -0,0 +1,37 @@
+/*--------------------------------*- 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 volScalarField;
+ object alpha.water;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions [0 0 0 0 0 0 0];
+
+internalField uniform 0;
+
+boundaryField
+{
+ rim
+ {
+ type zeroGradient;
+ }
+ front
+ {
+ type empty;
+ }
+ back
+ {
+ type empty;
+ }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/discInConstantFlow/0.orig/p_rgh b/tutorials/multiphase/interIsoFoam/discInConstantFlow/0.orig/p_rgh
new file mode 100644
index 0000000000000000000000000000000000000000..e4231ad4d6a4b7ba2e6dcb8c12411f193536ad46
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/discInConstantFlow/0.orig/p_rgh
@@ -0,0 +1,37 @@
+/*--------------------------------*- 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 volScalarField;
+ object p_rgh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions [1 -1 -2 0 0 0 0];
+
+internalField uniform 0;
+
+boundaryField
+{
+ rim
+ {
+ type zeroGradient;
+ }
+ front
+ {
+ type empty;
+ }
+ back
+ {
+ type empty;
+ }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/discInConstantFlow/Allclean b/tutorials/multiphase/interIsoFoam/discInConstantFlow/Allclean
new file mode 100755
index 0000000000000000000000000000000000000000..78449efd83660c6911b1e93e950765b750428078
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/discInConstantFlow/Allclean
@@ -0,0 +1,8 @@
+#!/bin/sh
+cd ${0%/*} || exit 1 # run from this directory
+
+# Source tutorial clean functions
+. $WM_PROJECT_DIR/bin/tools/CleanFunctions
+
+cleanCase
+\rm -rf 0
diff --git a/tutorials/multiphase/interIsoFoam/discInConstantFlow/Allrun b/tutorials/multiphase/interIsoFoam/discInConstantFlow/Allrun
new file mode 100755
index 0000000000000000000000000000000000000000..47d57eb4f7754a3c3dff129a5ce08642683e70b6
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/discInConstantFlow/Allrun
@@ -0,0 +1,11 @@
+cd ${0%/*} || exit 1 # Run from this directory
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+restore0Dir
+
+runApplication blockMesh
+runApplication setAlphaField
+
+runApplication $(getApplication)
diff --git a/tutorials/multiphase/interIsoFoam/discInConstantFlow/constant/g b/tutorials/multiphase/interIsoFoam/discInConstantFlow/constant/g
new file mode 100644
index 0000000000000000000000000000000000000000..efd2b0b9ea08489348efd846c2fc22a831e71e7f
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/discInConstantFlow/constant/g
@@ -0,0 +1,22 @@
+/*--------------------------------*- 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 uniformDimensionedVectorField;
+ location "constant";
+ object g;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions [0 1 -2 0 0 0 0];
+value ( 0 -9.81 0 );
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/discInConstantFlow/constant/transportProperties b/tutorials/multiphase/interIsoFoam/discInConstantFlow/constant/transportProperties
new file mode 100644
index 0000000000000000000000000000000000000000..3f872a9f54fdddb3ce7aacba5134e0f70d26cf99
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/discInConstantFlow/constant/transportProperties
@@ -0,0 +1,37 @@
+/*--------------------------------*- 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 "constant";
+ object transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+phases (water air);
+
+water
+{
+ transportModel Newtonian;
+ nu 0;
+ rho 1000;
+}
+
+air
+{
+ transportModel Newtonian;
+ nu 0;
+ rho 1;
+}
+
+sigma 0;
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/discInConstantFlow/constant/turbulenceProperties b/tutorials/multiphase/interIsoFoam/discInConstantFlow/constant/turbulenceProperties
new file mode 100644
index 0000000000000000000000000000000000000000..5eec04267266e7fd15e7701a875d683d5e658cd9
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/discInConstantFlow/constant/turbulenceProperties
@@ -0,0 +1,21 @@
+/*--------------------------------*- 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 "constant";
+ object turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType laminar;
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/discInConstantFlow/system/blockMeshDict b/tutorials/multiphase/interIsoFoam/discInConstantFlow/system/blockMeshDict
new file mode 100644
index 0000000000000000000000000000000000000000..27ef06566d1da00a450252d9bb6574f162e9df24
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/discInConstantFlow/system/blockMeshDict
@@ -0,0 +1,85 @@
+/*--------------------------------*- 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;
+ object blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 1;
+
+L 5;
+nx 100;
+
+y1 -.5;
+y2 .5;
+ny 1;
+
+H 3;
+nz 60;
+
+vertices
+(
+ (0 $y1 0)
+ ($L $y1 0)
+ ($L $y2 0)
+ (0 $y2 0)
+ (0 $y1 $H)
+ ($L $y1 $H)
+ ($L $y2 $H)
+ (0 $y2 $H)
+);
+
+blocks
+(
+ hex (0 1 2 3 4 5 6 7) ($nx $ny $nz) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+ rim
+ {
+ type patch;
+ faces
+ (
+ (4 5 6 7)
+ (0 4 7 3)
+ (1 2 6 5)
+ (0 3 2 1)
+ );
+ }
+ front
+ {
+ type empty;
+ faces
+ (
+ (0 1 5 4)
+ );
+ }
+ back
+ {
+ type empty;
+ faces
+ (
+ (2 3 7 6)
+ );
+ }
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/discInConstantFlow/system/controlDict b/tutorials/multiphase/interIsoFoam/discInConstantFlow/system/controlDict
new file mode 100644
index 0000000000000000000000000000000000000000..27c2889db194da104aa429c0ad0f5f08a6adedcc
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/discInConstantFlow/system/controlDict
@@ -0,0 +1,60 @@
+/*--------------------------------*- 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 controlDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application interIsoFoam;
+
+startFrom latestTime;
+
+startTime 0.0;
+
+stopAt endTime;
+
+endTime 4;
+
+writeControl adjustableRunTime;
+
+writeInterval 0.1;
+
+deltaT 1e-3;
+
+purgeWrite 0;
+
+writeFormat ascii;
+
+writePrecision 14;
+
+writeCompression uncompressed;
+
+timeFormat general;
+
+timePrecision 6;
+
+graphFormat raw;
+
+runTimeModifiable no;
+
+adjustTimeStep yes;
+
+maxCo 1e6;
+maxAlphaCo 0.5;
+
+maxDeltaT 1e6;
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/discInConstantFlow/system/fvSchemes b/tutorials/multiphase/interIsoFoam/discInConstantFlow/system/fvSchemes
new file mode 100644
index 0000000000000000000000000000000000000000..cc9aa4323c354e3606de9d4ca6bb79b838a1ba4d
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/discInConstantFlow/system/fvSchemes
@@ -0,0 +1,60 @@
+/*--------------------------------*- 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 fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+ default Euler;
+}
+
+gradSchemes
+{
+ default Gauss linear;
+}
+
+divSchemes
+{
+ div(rho*phi,U) Gauss limitedLinearV 1;
+ div(phi,alpha) Gauss vanLeer;
+ div(phirb,alpha) Gauss interfaceCompression;
+ div((muEff*dev(T(grad(U))))) Gauss linear;
+}
+
+laplacianSchemes
+{
+ default Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+ default linear;
+}
+
+snGradSchemes
+{
+ default corrected;
+}
+
+fluxRequired
+{
+ default no;
+ p_rgh;
+ pcorr;
+ alpha1;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/discInConstantFlow/system/fvSolution b/tutorials/multiphase/interIsoFoam/discInConstantFlow/system/fvSolution
new file mode 100644
index 0000000000000000000000000000000000000000..444dba14b893c1097b855cfc92929e6c70035c6b
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/discInConstantFlow/system/fvSolution
@@ -0,0 +1,44 @@
+/*--------------------------------*- 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;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+ "alpha.*"
+ {
+ isoFaceTol 1e-8;
+ surfCellTol 1e-8;
+ nAlphaBounds 3;
+ snapTol 0;
+ clip false;
+ cAlpha 1;
+ nAlphaSubCycles 1;
+ }
+}
+
+
+PIMPLE
+{
+ frozenFlow yes;
+ momentumPredictor no;
+ nCorrectors -1;
+ nNonOrthogonalCorrectors -1;
+ pRefCell 0;
+ pRefValue 0;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/discInConstantFlow/system/setAlphaFieldDict b/tutorials/multiphase/interIsoFoam/discInConstantFlow/system/setAlphaFieldDict
new file mode 100644
index 0000000000000000000000000000000000000000..6bbb616b549e174c1360681c03405f0a628bcf8d
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/discInConstantFlow/system/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);
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/0.orig/U b/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/0.orig/U
new file mode 100644
index 0000000000000000000000000000000000000000..990863a30b3007a76fa337f6298d4ee849dad03b
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/0.orig/U
@@ -0,0 +1,39 @@
+/*--------------------------------*- 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 volVectorField;
+ location "0";
+ object U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions [0 1 -1 0 0 0 0];
+
+internalField uniform (0 0 0);
+
+boundaryField
+{
+ rim
+ {
+ type zeroGradient;
+ }
+ front
+ {
+ type empty;
+ }
+ back
+ {
+ type empty;
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/0.orig/alpha.water b/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/0.orig/alpha.water
new file mode 100644
index 0000000000000000000000000000000000000000..5b041ff5f7f40d24b56c7b7acbbe179e92ef7b55
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/0.orig/alpha.water
@@ -0,0 +1,38 @@
+/*--------------------------------*- 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 volScalarField;
+ object alpha;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions [0 0 0 0 0 0 0];
+
+internalField uniform 0;
+
+boundaryField
+{
+ rim
+ {
+ type fixedValue;
+ value uniform 0;
+ }
+ front
+ {
+ type empty;
+ }
+ back
+ {
+ type empty;
+ }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/0.orig/p_rgh b/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/0.orig/p_rgh
new file mode 100644
index 0000000000000000000000000000000000000000..e4231ad4d6a4b7ba2e6dcb8c12411f193536ad46
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/0.orig/p_rgh
@@ -0,0 +1,37 @@
+/*--------------------------------*- 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 volScalarField;
+ object p_rgh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions [1 -1 -2 0 0 0 0];
+
+internalField uniform 0;
+
+boundaryField
+{
+ rim
+ {
+ type zeroGradient;
+ }
+ front
+ {
+ type empty;
+ }
+ back
+ {
+ type empty;
+ }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/Allclean b/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/Allclean
new file mode 100755
index 0000000000000000000000000000000000000000..c11c2183476730adbcb8206b57df5f367f29be34
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/Allclean
@@ -0,0 +1,7 @@
+#!/bin/bash
+
+# Source tutorial clean functions
+. $WM_PROJECT_DIR/bin/tools/CleanFunctions
+
+cleanCase
+\rm -rf 0
diff --git a/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/Allrun b/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/Allrun
new file mode 100755
index 0000000000000000000000000000000000000000..474424491f37b168313472d53fd3c819656c1d4c
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/Allrun
@@ -0,0 +1,17 @@
+#!/bin/sh
+cd ${0%/*} || exit 1 # Run from this directory
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+restore0Dir
+
+runApplication blockMesh
+runApplication -s 1 topoSet
+runApplication -s 1 refineMesh -overwrite
+runApplication -s 2 topoSet -dict system/topoSetDict2
+runApplication -s 2 refineMesh -overwrite
+
+runApplication setAlphaField
+
+runApplication $(getApplication)
diff --git a/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/constant/g b/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/constant/g
new file mode 100644
index 0000000000000000000000000000000000000000..efd2b0b9ea08489348efd846c2fc22a831e71e7f
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/constant/g
@@ -0,0 +1,22 @@
+/*--------------------------------*- 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 uniformDimensionedVectorField;
+ location "constant";
+ object g;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions [0 1 -2 0 0 0 0];
+value ( 0 -9.81 0 );
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/constant/transportProperties b/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/constant/transportProperties
new file mode 100644
index 0000000000000000000000000000000000000000..6e38d0134f3347395c5bee4d2bee64a625bdf8a8
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/constant/transportProperties
@@ -0,0 +1,86 @@
+/*--------------------------------*- 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 "constant";
+ object transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+DT DT [ 0 2 -1 0 0 0 0 ] 0;
+
+//Input for OF-2.3.0 and later
+phases (water air);
+
+water
+{
+ transportModel Newtonian;
+ nu [0 2 -1 0 0 0 0] 1e-06;
+ rho [1 -3 0 0 0 0 0] 1000;
+}
+
+air
+{
+ transportModel Newtonian;
+ nu [0 2 -1 0 0 0 0] 1.48e-05;
+ rho [1 -3 0 0 0 0 0] 1;
+}
+
+//Input for OF-2.2.0 and earlier
+
+phase1
+{
+ transportModel Newtonian;
+ nu nu [ 0 2 -1 0 0 0 0 ] 0;
+ rho rho [ 1 -3 0 0 0 0 0 ] 1000;
+ CrossPowerLawCoeffs
+ {
+ nu0 nu0 [ 0 2 -1 0 0 0 0 ] 1e-06;
+ nuInf nuInf [ 0 2 -1 0 0 0 0 ] 1e-06;
+ m m [ 0 0 1 0 0 0 0 ] 1;
+ n n [ 0 0 0 0 0 0 0 ] 0;
+ }
+
+ BirdCarreauCoeffs
+ {
+ nu0 nu0 [ 0 2 -1 0 0 0 0 ] 0.0142515;
+ nuInf nuInf [ 0 2 -1 0 0 0 0 ] 1e-06;
+ k k [ 0 0 1 0 0 0 0 ] 99.6;
+ n n [ 0 0 0 0 0 0 0 ] 0.1003;
+ }
+}
+
+phase2
+{
+ transportModel Newtonian;
+ nu nu [ 0 2 -1 0 0 0 0 ] 0;
+ rho rho [ 1 -3 0 0 0 0 0 ] 1;
+ CrossPowerLawCoeffs
+ {
+ nu0 nu0 [ 0 2 -1 0 0 0 0 ] 1e-06;
+ nuInf nuInf [ 0 2 -1 0 0 0 0 ] 1e-06;
+ m m [ 0 0 1 0 0 0 0 ] 1;
+ n n [ 0 0 0 0 0 0 0 ] 0;
+ }
+
+ BirdCarreauCoeffs
+ {
+ nu0 nu0 [ 0 2 -1 0 0 0 0 ] 0.0142515;
+ nuInf nuInf [ 0 2 -1 0 0 0 0 ] 1e-06;
+ k k [ 0 0 1 0 0 0 0 ] 99.6;
+ n n [ 0 0 0 0 0 0 0 ] 0.1003;
+ }
+}
+
+sigma sigma [ 1 0 -2 0 0 0 0 ] 0.0;
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/constant/turbulenceProperties b/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/constant/turbulenceProperties
new file mode 100644
index 0000000000000000000000000000000000000000..5eec04267266e7fd15e7701a875d683d5e658cd9
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/constant/turbulenceProperties
@@ -0,0 +1,21 @@
+/*--------------------------------*- 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 "constant";
+ object turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType laminar;
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/system/blockMeshDict b/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/system/blockMeshDict
new file mode 100644
index 0000000000000000000000000000000000000000..b8c7f146193d157efa7490069e27ae8ecc2927dc
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/system/blockMeshDict
@@ -0,0 +1,85 @@
+/*--------------------------------*- 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;
+ object blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 1;
+
+L 1;
+nx 100;
+
+y1 -.01;
+y2 .01;
+ny 1;
+
+H 1;
+nz 100;
+
+vertices
+(
+ (0 $y1 0)
+ ($L $y1 0)
+ ($L $y2 0)
+ (0 $y2 0)
+ (0 $y1 $H)
+ ($L $y1 $H)
+ ($L $y2 $H)
+ (0 $y2 $H)
+);
+
+blocks
+(
+ hex (0 1 2 3 4 5 6 7) ($nx $ny $nz) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+ rim
+ {
+ type patch;
+ faces
+ (
+ (4 5 6 7)
+ (0 4 7 3)
+ (1 2 6 5)
+ (0 3 2 1)
+ );
+ }
+ front
+ {
+ type empty;
+ faces
+ (
+ (0 1 5 4)
+ );
+ }
+ back
+ {
+ type empty;
+ faces
+ (
+ (2 3 7 6)
+ );
+ }
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/system/controlDict b/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/system/controlDict
new file mode 100644
index 0000000000000000000000000000000000000000..e72e9490476d1b9e263057853aa56f41d58d2161
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/system/controlDict
@@ -0,0 +1,81 @@
+/*--------------------------------*- 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 controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application interIsoFoam;
+
+startFrom latestTime;
+
+startTime 0.0;
+
+stopAt endTime;
+
+endTime 8;
+
+writeControl adjustableRunTime;
+
+writeInterval 0.1;
+
+deltaT 0.002;
+
+purgeWrite 0;
+
+writeFormat binary;
+
+writePrecision 14;
+
+writeCompression uncompressed;
+
+timeFormat general;
+
+timePrecision 6;
+
+graphFormat raw;
+
+runTimeModifiable no;
+
+adjustTimeStep no;
+
+maxCo 1e10;
+maxAlphaCo 0.5;
+
+maxDeltaT 0.05;
+
+functions
+{
+ setVelocity
+ {
+ type setFlow;
+ libs ("libfieldFunctionObjects.so");
+ writeControl writeTime;
+ mode vortex2D;
+ scale sine;
+ scaleCoeffs
+ {
+ amplitude 1;
+ frequency 0.0625; // = 1/16
+ t0 -4; // want cos -> time shift = -(pi/2)/(2 pi f)
+ scale 1;
+ level 0;
+ }
+ origin (0 0 0);
+ refDir (1 0 0);
+ axis (0 1 0);
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/system/decomposeParDict b/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/system/decomposeParDict
new file mode 100644
index 0000000000000000000000000000000000000000..02204f4d2d9aa9dc5c673422e5d32a7779951d60
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/system/decomposeParDict
@@ -0,0 +1,61 @@
+/*---------------------------------------------------------------------------*\
+| ========= | |
+| \\ / 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 decomposeParDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+numberOfSubdomains 4;
+
+method simple;
+
+simpleCoeffs
+{
+ n (2 1 2);
+ delta 0.001;
+}
+
+hierarchicalCoeffs
+{
+ n (1 1 1);
+ delta 0.001;
+ order xyz;
+}
+
+metisCoeffs
+{
+ processorWeights
+ (
+ 1
+ 1
+ 1
+ 1
+ );
+}
+
+manualCoeffs
+{
+ dataFile "";
+}
+
+distributed no;
+
+roots
+(
+);
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/system/fvSchemes b/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/system/fvSchemes
new file mode 100644
index 0000000000000000000000000000000000000000..4927912a7f76cadc3bfd5d332c4d3dbde7bce0b4
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/system/fvSchemes
@@ -0,0 +1,51 @@
+/*--------------------------------*- 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 fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+ default Euler;
+}
+
+gradSchemes
+{
+ default Gauss linear;
+}
+
+divSchemes
+{
+ div(rhoPhi,U) Gauss limitedLinearV 1;
+ div(phi,alpha) Gauss vanLeer;
+ div(phirb,alpha) Gauss linear;
+}
+
+laplacianSchemes
+{
+ default Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+ default linear;
+}
+
+snGradSchemes
+{
+ default corrected;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/system/fvSolution b/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/system/fvSolution
new file mode 100644
index 0000000000000000000000000000000000000000..25e58664db2b2804627b7b8c0c88cfa24319ea1a
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/system/fvSolution
@@ -0,0 +1,43 @@
+/*--------------------------------*- 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;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+ "alpha.*"
+ {
+ isoFaceTol 1e-8;
+ surfCellTol 1e-8;
+ nAlphaBounds 3;
+ snapTol 0;
+ clip false;
+
+ nAlphaSubCycles 1;
+ cAlpha 1;
+ }
+}
+
+PIMPLE
+{
+ frozenFlow yes;
+ momentumPredictor no;
+ nCorrectors -1;
+ nNonOrthogonalCorrectors -1;
+ pRefCell 0;
+ pRefValue 0;
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/system/refineMeshDict b/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/system/refineMeshDict
new file mode 100644
index 0000000000000000000000000000000000000000..e3efa074d9f446aeb98ad981687b55b89e1b9744
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/system/refineMeshDict
@@ -0,0 +1,62 @@
+/*--------------------------------*- 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;
+ object refineMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Cells to refine; name of cell set
+set c1;
+
+// Type of coordinate system:
+// - global : coordinate system same for every cell. Usually aligned with
+// x,y,z axis. Specify in globalCoeffs section below.
+// - patchLocal : coordinate system different for every cell. Specify in
+// patchLocalCoeffs section below.
+coordinateSystem global;
+//coordinateSystem patchLocal;
+
+// .. and its coefficients. x,y in this case. (normal direction is calculated
+// as tan1^tan2)
+globalCoeffs
+{
+ tan1 (1 0 0);
+ tan2 (0 1 0);
+}
+
+patchLocalCoeffs
+{
+ patch outside; // Normal direction is facenormal of zero'th face of patch
+ tan1 (1 0 0);
+}
+
+// List of directions to refine
+directions
+(
+ tan1
+// tan2
+ normal
+);
+
+// Whether to use hex topology. This will
+// - if patchLocal: all cells on selected patch should be hex
+// - split all hexes in 2x2x2 through the middle of edges.
+useHexTopology true;
+
+// Cut purely geometric (will cut hexes through vertices) or take topology
+// into account. Incompatible with useHexTopology
+geometricCut false;
+
+// Write meshes from intermediate steps
+writeMesh false;
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/system/setAlphaFieldDict b/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/system/setAlphaFieldDict
new file mode 100644
index 0000000000000000000000000000000000000000..f4424dabcbbb45872abfb089bc6a017d2a9412bb
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/system/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.15;
+direction (0 1 0);
+centre (0.5 0 0.75);
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/system/topoSetDict b/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/system/topoSetDict
new file mode 100644
index 0000000000000000000000000000000000000000..bcd7ded2824cc05a8ff8f01e140c627d931fa420
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/system/topoSetDict
@@ -0,0 +1,37 @@
+/*--------------------------------*- 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;
+ object topoSetDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+actions
+(
+ {
+ name c1;
+ type cellSet;
+ action new;
+ source boxToCell;
+ sourceInfo
+ {
+ boxes
+ (
+ (.1 -1e10 0.08) (.9 1e10 0.2)
+ (.75 -1e10 0.2) (.92 1e10 0.9)
+ (.1 -1e10 0.2) (.25 1e10 0.7)
+ );
+ }
+ }
+);
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/system/topoSetDict2 b/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/system/topoSetDict2
new file mode 100644
index 0000000000000000000000000000000000000000..04f70800238a9621a8ef29aec361340a4c43d528
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/discInReversedVortexFlow/system/topoSetDict2
@@ -0,0 +1,35 @@
+/*--------------------------------*- 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;
+ object topoSetDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+actions
+(
+ {
+ name c1;
+ type cellSet;
+ action new;
+ source boxToCell;
+ sourceInfo
+ {
+ boxes
+ (
+ (.4 -1e10 0.09) (.7 1e10 0.15)
+ );
+ }
+ }
+);
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/0.orig/U b/tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/0.orig/U
new file mode 100644
index 0000000000000000000000000000000000000000..990863a30b3007a76fa337f6298d4ee849dad03b
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/0.orig/U
@@ -0,0 +1,39 @@
+/*--------------------------------*- 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 volVectorField;
+ location "0";
+ object U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions [0 1 -1 0 0 0 0];
+
+internalField uniform (0 0 0);
+
+boundaryField
+{
+ rim
+ {
+ type zeroGradient;
+ }
+ front
+ {
+ type empty;
+ }
+ back
+ {
+ type empty;
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/0.orig/alpha.water b/tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/0.orig/alpha.water
new file mode 100644
index 0000000000000000000000000000000000000000..739e3e2dd2bb4d3f54100477e6a87064e184c0cf
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/0.orig/alpha.water
@@ -0,0 +1,38 @@
+/*--------------------------------*- 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 volScalarField;
+ object alpha.water;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions [0 0 0 0 0 0 0];
+
+internalField uniform 0;
+
+boundaryField
+{
+ rim
+ {
+ type fixedValue;
+ value uniform 0;
+ }
+ front
+ {
+ type empty;
+ }
+ back
+ {
+ type empty;
+ }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/0.orig/p_rgh b/tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/0.orig/p_rgh
new file mode 100644
index 0000000000000000000000000000000000000000..e4231ad4d6a4b7ba2e6dcb8c12411f193536ad46
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/0.orig/p_rgh
@@ -0,0 +1,37 @@
+/*--------------------------------*- 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 volScalarField;
+ object p_rgh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions [1 -1 -2 0 0 0 0];
+
+internalField uniform 0;
+
+boundaryField
+{
+ rim
+ {
+ type zeroGradient;
+ }
+ front
+ {
+ type empty;
+ }
+ back
+ {
+ type empty;
+ }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/Allclean b/tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/Allclean
new file mode 100755
index 0000000000000000000000000000000000000000..97f2707db81d7086eac366adedadf01c2a642142
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/Allclean
@@ -0,0 +1,9 @@
+#!/bin/sh
+cd ${0%/*} || exit 1 # run from this directory
+
+# Source tutorial clean functions
+. $WM_PROJECT_DIR/bin/tools/CleanFunctions
+
+cleanCase
+
+\rm -rf 0
diff --git a/tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/Allrun b/tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/Allrun
new file mode 100755
index 0000000000000000000000000000000000000000..317234a300a109eef71b899d681e149c68e34df8
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/Allrun
@@ -0,0 +1,19 @@
+#!/bin/sh
+cd ${0%/*} || exit 1 # run from this directory
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+application=$(getApplication)
+
+restore0Dir
+
+runApplication blockMesh
+
+runApplication -s postProcess ${application} -postProcess -time 0
+
+runApplication setAlphaField
+
+runApplication setFields
+
+runApplication ${application}
diff --git a/tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/constant/g b/tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/constant/g
new file mode 100644
index 0000000000000000000000000000000000000000..efd2b0b9ea08489348efd846c2fc22a831e71e7f
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/constant/g
@@ -0,0 +1,22 @@
+/*--------------------------------*- 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 uniformDimensionedVectorField;
+ location "constant";
+ object g;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions [0 1 -2 0 0 0 0];
+value ( 0 -9.81 0 );
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/constant/transportProperties b/tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/constant/transportProperties
new file mode 100644
index 0000000000000000000000000000000000000000..3f872a9f54fdddb3ce7aacba5134e0f70d26cf99
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/constant/transportProperties
@@ -0,0 +1,37 @@
+/*--------------------------------*- 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 "constant";
+ object transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+phases (water air);
+
+water
+{
+ transportModel Newtonian;
+ nu 0;
+ rho 1000;
+}
+
+air
+{
+ transportModel Newtonian;
+ nu 0;
+ rho 1;
+}
+
+sigma 0;
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/constant/turbulenceProperties b/tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/constant/turbulenceProperties
new file mode 100644
index 0000000000000000000000000000000000000000..5eec04267266e7fd15e7701a875d683d5e658cd9
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/constant/turbulenceProperties
@@ -0,0 +1,21 @@
+/*--------------------------------*- 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 "constant";
+ object turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType laminar;
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/system/blockMeshDict b/tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/system/blockMeshDict
new file mode 100644
index 0000000000000000000000000000000000000000..b8c7f146193d157efa7490069e27ae8ecc2927dc
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/system/blockMeshDict
@@ -0,0 +1,85 @@
+/*--------------------------------*- 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;
+ object blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 1;
+
+L 1;
+nx 100;
+
+y1 -.01;
+y2 .01;
+ny 1;
+
+H 1;
+nz 100;
+
+vertices
+(
+ (0 $y1 0)
+ ($L $y1 0)
+ ($L $y2 0)
+ (0 $y2 0)
+ (0 $y1 $H)
+ ($L $y1 $H)
+ ($L $y2 $H)
+ (0 $y2 $H)
+);
+
+blocks
+(
+ hex (0 1 2 3 4 5 6 7) ($nx $ny $nz) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+ rim
+ {
+ type patch;
+ faces
+ (
+ (4 5 6 7)
+ (0 4 7 3)
+ (1 2 6 5)
+ (0 3 2 1)
+ );
+ }
+ front
+ {
+ type empty;
+ faces
+ (
+ (0 1 5 4)
+ );
+ }
+ back
+ {
+ type empty;
+ faces
+ (
+ (2 3 7 6)
+ );
+ }
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/system/controlDict b/tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/system/controlDict
new file mode 100644
index 0000000000000000000000000000000000000000..7cdc121f32e424c34649836eee576b2e8ff9971b
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/system/controlDict
@@ -0,0 +1,76 @@
+/*--------------------------------*- 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 controlDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application interIsoFoam;
+
+startFrom latestTime;
+
+startTime 0.0;
+
+stopAt endTime;
+
+endTime 2;
+
+writeControl adjustableRunTime;
+
+writeInterval 0.05;
+
+deltaT 0.001;
+
+purgeWrite 0;
+
+writeFormat ascii;
+
+writePrecision 14;
+
+writeCompression uncompressed;
+
+timeFormat general;
+
+timePrecision 6;
+
+graphFormat raw;
+
+runTimeModifiable no;
+
+adjustTimeStep no;
+
+maxCo 1e6;
+maxAlphaCo 0.5;
+
+maxDeltaT 0.05;
+
+functions
+{
+ setVelocity
+ {
+ type setFlow;
+ libs ("libfieldFunctionObjects.so");
+ writeControl writeTime;
+ mode rotation;
+ scale 1;
+ reverseTime 1;
+ omega 6.28318530718;
+ origin (0.5 0 0.5);
+ refDir (1 0 0);
+ axis (0 1 0);
+ }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/system/fvSchemes b/tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/system/fvSchemes
new file mode 100644
index 0000000000000000000000000000000000000000..9d9a3c05434d8a6469e7edbaeaf8fbbdc06d7b6b
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/system/fvSchemes
@@ -0,0 +1,52 @@
+/*--------------------------------*- 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 fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+ default Euler;
+}
+
+gradSchemes
+{
+ default Gauss linear;
+}
+
+divSchemes
+{
+ div(rhoPhi,U) Gauss limitedLinearV 1;
+ div(phi,alpha) Gauss vanLeer;
+ div(phirb,alpha) Gauss linear;
+ div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
+}
+
+laplacianSchemes
+{
+ default Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+ default linear;
+}
+
+snGradSchemes
+{
+ default corrected;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/system/fvSolution b/tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/system/fvSolution
new file mode 100644
index 0000000000000000000000000000000000000000..84b3554cf8ee15e7f8289aede6a83be74646868e
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/system/fvSolution
@@ -0,0 +1,43 @@
+/*--------------------------------*- 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;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+ "alpha.*"
+ {
+ isoFaceTol 1e-8;
+ surfCellTol 1e-8;
+ snapTol 0;
+ nAlphaBounds 3;
+ clip false;
+ nAlphaSubCycles 1;
+ cAlpha 1;
+ }
+}
+
+
+PIMPLE
+{
+ frozenFlow yes;
+ momentumPredictor no;
+ nCorrectors -1;
+ nNonOrthogonalCorrectors -1;
+ pRefCell 0;
+ pRefValue 0;
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/system/setAlphaFieldDict b/tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/system/setAlphaFieldDict
new file mode 100644
index 0000000000000000000000000000000000000000..ac15050b1fd4229accbcb7f63b6fbf7c053eecf0
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/system/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.15;
+direction (0 1 0);
+centre (0.5 0 0.75);
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/system/setFieldsDict b/tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/system/setFieldsDict
new file mode 100644
index 0000000000000000000000000000000000000000..f933d6fa14c2e1eba7e781e8d98b247f9eb6c311
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/notchedDiscInSolidBodyRotation/system/setFieldsDict
@@ -0,0 +1,31 @@
+/*--------------------------------*- 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 setFieldsDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+regions
+(
+ boxToCell
+ {
+ box (.47 -1 0) (.53 1 .85);
+ fieldValues
+ (
+ volScalarFieldValue alpha.water 0
+ );
+ }
+);
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/0.orig/U b/tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/0.orig/U
new file mode 100644
index 0000000000000000000000000000000000000000..ca61b972e3253f557dfcfa55ffa7e37a47d00b25
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/0.orig/U
@@ -0,0 +1,31 @@
+/*--------------------------------*- 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 volVectorField;
+ location "0";
+ object U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions [0 1 -1 0 0 0 0];
+
+internalField uniform (0 0 0);
+
+boundaryField
+{
+ sides
+ {
+ type zeroGradient;
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/0.orig/alpha.water b/tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/0.orig/alpha.water
new file mode 100644
index 0000000000000000000000000000000000000000..0ad7555652f799391390d97845cee551e2257559
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/0.orig/alpha.water
@@ -0,0 +1,30 @@
+/*--------------------------------*- 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 volScalarField;
+ object alpha.water;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions [0 0 0 0 0 0 0];
+
+internalField uniform 0;
+
+boundaryField
+{
+ sides
+ {
+ type fixedValue;
+ value uniform 0;
+ }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/0.orig/p_rgh b/tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/0.orig/p_rgh
new file mode 100644
index 0000000000000000000000000000000000000000..b2e6dadb332131cf66256ef9674c80cbc8b71d4d
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/0.orig/p_rgh
@@ -0,0 +1,30 @@
+/*--------------------------------*- 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 volScalarField;
+ object p_rgh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions [1 -1 -2 0 0 0 0];
+
+internalField uniform 0;
+
+boundaryField
+{
+ sides
+ {
+ type fixedValue;
+ value uniform 0;
+ }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/Allclean b/tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/Allclean
new file mode 100755
index 0000000000000000000000000000000000000000..c11c2183476730adbcb8206b57df5f367f29be34
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/Allclean
@@ -0,0 +1,7 @@
+#!/bin/bash
+
+# Source tutorial clean functions
+. $WM_PROJECT_DIR/bin/tools/CleanFunctions
+
+cleanCase
+\rm -rf 0
diff --git a/tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/Allrun b/tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/Allrun
new file mode 100755
index 0000000000000000000000000000000000000000..f21f9b5925b1ad60cf46ff83a8db1bef052efdff
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/Allrun
@@ -0,0 +1,13 @@
+#!/bin/sh
+cd ${0%/*} || exit 1 # Run from this directory
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+restore0Dir
+
+runApplication blockMesh
+
+runApplication setAlphaField
+runApplication decomposePar
+runParallel $(getApplication)
diff --git a/tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/constant/g b/tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/constant/g
new file mode 100644
index 0000000000000000000000000000000000000000..efd2b0b9ea08489348efd846c2fc22a831e71e7f
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/constant/g
@@ -0,0 +1,22 @@
+/*--------------------------------*- 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 uniformDimensionedVectorField;
+ location "constant";
+ object g;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions [0 1 -2 0 0 0 0];
+value ( 0 -9.81 0 );
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/constant/transportProperties b/tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/constant/transportProperties
new file mode 100644
index 0000000000000000000000000000000000000000..6b4453c9f7c94fc633cfa41f2d9cbcca258e294c
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/constant/transportProperties
@@ -0,0 +1,37 @@
+/*--------------------------------*- 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 "constant";
+ object transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+phases (water air);
+
+water
+{
+ transportModel Newtonian;
+ nu [0 2 -1 0 0 0 0] 1e-06;
+ rho [1 -3 0 0 0 0 0] 1000;
+}
+
+air
+{
+ transportModel Newtonian;
+ nu [0 2 -1 0 0 0 0] 1.48e-05;
+ rho [1 -3 0 0 0 0 0] 1;
+}
+
+sigma sigma [ 1 0 -2 0 0 0 0 ] 0.0;
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/constant/turbulenceProperties b/tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/constant/turbulenceProperties
new file mode 100644
index 0000000000000000000000000000000000000000..5eec04267266e7fd15e7701a875d683d5e658cd9
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/constant/turbulenceProperties
@@ -0,0 +1,21 @@
+/*--------------------------------*- 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 "constant";
+ object turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType laminar;
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/system/blockMeshDict b/tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/system/blockMeshDict
new file mode 100644
index 0000000000000000000000000000000000000000..a04acb0bb2c9fb52521c57f38602eed565e85392
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/system/blockMeshDict
@@ -0,0 +1,71 @@
+/*--------------------------------*- 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;
+ object blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 1;
+
+x1 0;
+x2 1;
+y1 0;
+y2 1;
+z1 0;
+z2 1;
+nx 64;
+ny 64;
+nz 64;
+
+vertices
+(
+ ($x1 $y1 $z1)
+ ($x2 $y1 $z1)
+ ($x2 $y2 $z1)
+ ($x1 $y2 $z1)
+ ($x1 $y1 $z2)
+ ($x2 $y1 $z2)
+ ($x2 $y2 $z2)
+ ($x1 $y2 $z2)
+);
+
+blocks
+(
+ hex (0 1 2 3 4 5 6 7) ($nx $ny $nz) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+ sides
+ {
+ type patch;
+ faces
+ (
+ (4 5 6 7)
+ (0 4 7 3)
+ (1 2 6 5)
+ (0 3 2 1)
+ (0 1 5 4)
+ (2 3 7 6)
+ );
+ }
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/system/controlDict b/tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/system/controlDict
new file mode 100644
index 0000000000000000000000000000000000000000..7de4706a29c4793dfec5b079f6f5012ddf9e2e32
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/system/controlDict
@@ -0,0 +1,81 @@
+/*--------------------------------*- 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 controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application interIsoFoam;
+
+startFrom startTime;
+
+startTime 0.0;
+
+stopAt endTime;
+
+endTime 3.0;
+
+writeControl adjustableRunTime;
+
+writeInterval 0.1;
+
+deltaT 0.0025;
+
+purgeWrite 0;
+
+writeFormat ascii;
+
+writePrecision 14;
+
+writeCompression uncompressed;
+
+timeFormat general;
+
+timePrecision 6;
+
+graphFormat raw;
+
+runTimeModifiable yes;
+
+adjustTimeStep no;
+
+maxCo 1e10;
+maxAlphaCo 0.5;
+
+maxDeltaT 0.2;
+
+functions
+{
+ setVelocity
+ {
+ type setFlow;
+ libs ("libfieldFunctionObjects.so");
+ writeControl writeTime;
+ mode vortex3D;
+ scale sine;
+ scaleCoeffs
+ {
+ amplitude 1;
+ frequency 0.16666;// = 1/6
+ t0 -1.5; // want cos -> time shift = -(pi/2)/(2 pi f)
+ scale 1;
+ level 0;
+ }
+ origin (0 0 0);
+ refDir (1 0 0);
+ axis (0 1 0);
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/system/decomposeParDict b/tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/system/decomposeParDict
new file mode 100644
index 0000000000000000000000000000000000000000..b8a6f06f4f58af0561a4a1d43652fbd6ba348b8c
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/system/decomposeParDict
@@ -0,0 +1,31 @@
+/*---------------------------------------------------------------------------*\
+| ========= | |
+| \\ / 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 decomposeParDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+numberOfSubdomains 8;
+
+method simple;
+
+simpleCoeffs
+{
+ n (2 2 2);
+ delta 0.001;
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/system/fvSchemes b/tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/system/fvSchemes
new file mode 100644
index 0000000000000000000000000000000000000000..4927912a7f76cadc3bfd5d332c4d3dbde7bce0b4
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/system/fvSchemes
@@ -0,0 +1,51 @@
+/*--------------------------------*- 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 fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+ default Euler;
+}
+
+gradSchemes
+{
+ default Gauss linear;
+}
+
+divSchemes
+{
+ div(rhoPhi,U) Gauss limitedLinearV 1;
+ div(phi,alpha) Gauss vanLeer;
+ div(phirb,alpha) Gauss linear;
+}
+
+laplacianSchemes
+{
+ default Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+ default linear;
+}
+
+snGradSchemes
+{
+ default corrected;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/system/fvSolution b/tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/system/fvSolution
new file mode 100644
index 0000000000000000000000000000000000000000..d92af9b76b62894885fc8417345c218513c46ff7
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/system/fvSolution
@@ -0,0 +1,43 @@
+/*--------------------------------*- 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;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+ "alpha.*"
+ {
+ isoFaceTol 1e-8;
+ surfCellTol 1e-8;
+ snapTol 0;
+ nAlphaBounds 4;
+ clip false;
+ nAlphaSubCycles 1;
+ cAlpha 1;
+ }
+}
+
+PIMPLE
+{
+ frozenFlow yes;
+ momentumPredictor no;
+ nCorrectors -1;
+ nNonOrthogonalCorrectors -1;
+ pRefCell 0;
+ pRefValue 0;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/system/setAlphaFieldDict b/tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/system/setAlphaFieldDict
new file mode 100644
index 0000000000000000000000000000000000000000..8990547c1247a42fdacf7dcc24d0aa9b81062b86
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/sphereInReversedVortexFlow/system/setAlphaFieldDict
@@ -0,0 +1,23 @@
+/*--------------------------------*- 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 sphere;
+radius 0.15;
+centre (0.35 0.35 0.35);
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/standingWave/0.orig/U b/tutorials/multiphase/interIsoFoam/standingWave/0.orig/U
new file mode 100644
index 0000000000000000000000000000000000000000..c1140f952980372997a5200b0f04756df4b5494d
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/standingWave/0.orig/U
@@ -0,0 +1,51 @@
+/*--------------------------------*- 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 volVectorField;
+ location "0";
+ object U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions [0 1 -1 0 0 0 0];
+
+internalField uniform (0 0 0);
+
+boundaryField
+{
+ left
+ {
+ type slip;
+ }
+ right
+ {
+ type slip;
+ }
+ top
+ {
+ type slip;
+ }
+ bottom
+ {
+ type slip;
+ }
+ front
+ {
+ type empty;
+ }
+ back
+ {
+ type empty;
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/standingWave/0.orig/alpha.water b/tutorials/multiphase/interIsoFoam/standingWave/0.orig/alpha.water
new file mode 100644
index 0000000000000000000000000000000000000000..c425019d1254d3ff0c0ec66b9a962debaef339ba
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/standingWave/0.orig/alpha.water
@@ -0,0 +1,49 @@
+/*--------------------------------*- 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 volScalarField;
+ object alpha.water;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions [0 0 0 0 0 0 0];
+
+internalField uniform 0;
+
+boundaryField
+{
+ left
+ {
+ type zeroGradient;
+ }
+ right
+ {
+ type zeroGradient;
+ }
+ top
+ {
+ type zeroGradient;
+ }
+ bottom
+ {
+ type zeroGradient;
+ }
+ front
+ {
+ type empty;
+ }
+ back
+ {
+ type empty;
+ }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/standingWave/0.orig/p_rgh b/tutorials/multiphase/interIsoFoam/standingWave/0.orig/p_rgh
new file mode 100644
index 0000000000000000000000000000000000000000..4b7ce2d53863bd5584a5c8cdfbe3e7a92253c262
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/standingWave/0.orig/p_rgh
@@ -0,0 +1,53 @@
+/*--------------------------------*- 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 volScalarField;
+ object p_rgh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions [1 -1 -2 0 0 0 0];
+
+internalField uniform 0;
+
+boundaryField
+{
+ left
+ {
+ type fixedFluxPressure;
+ value uniform 0;
+ }
+ right
+ {
+ type fixedFluxPressure;
+ value uniform 0;
+ }
+ top
+ {
+ type fixedFluxPressure;
+ value uniform 0;
+ }
+ bottom
+ {
+ type fixedFluxPressure;
+ value uniform 0;
+ }
+ front
+ {
+ type empty;
+ }
+ back
+ {
+ type empty;
+ }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/standingWave/Allclean b/tutorials/multiphase/interIsoFoam/standingWave/Allclean
new file mode 100755
index 0000000000000000000000000000000000000000..9738fda663b36f826c5f2ddb759bda3d4a0510de
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/standingWave/Allclean
@@ -0,0 +1,9 @@
+#!/bin/sh
+cd ${0%/*} || exit 1 # Run from this directory
+
+# Source tutorial clean functions
+. $WM_PROJECT_DIR/bin/tools/CleanFunctions
+
+cleanCase
+
+\rm -rf 0
diff --git a/tutorials/multiphase/interIsoFoam/standingWave/Allrun b/tutorials/multiphase/interIsoFoam/standingWave/Allrun
new file mode 100755
index 0000000000000000000000000000000000000000..3f7a4cb538a521d93128f436ef60a8166df90528
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/standingWave/Allrun
@@ -0,0 +1,17 @@
+#!/bin/sh
+cd ${0%/*} || exit 1 # Run from this directory
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+restore0Dir
+
+runApplication blockMesh
+runApplication -s 1 topoSet -dict system/topoSetDict1
+runApplication -s 1 refineMesh -dict system/refineMeshDict1 -overwrite
+runApplication -s 2 topoSet -dict system/topoSetDict2
+runApplication -s 2 refineMesh -dict system/refineMeshDict2 -overwrite
+
+runApplication setAlphaField
+
+runApplication $(getApplication)
diff --git a/tutorials/multiphase/interIsoFoam/standingWave/constant/dynamicMeshDict b/tutorials/multiphase/interIsoFoam/standingWave/constant/dynamicMeshDict
new file mode 100644
index 0000000000000000000000000000000000000000..0d234e3050053c7590aa9b6089f7c448c2932b97
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/standingWave/constant/dynamicMeshDict
@@ -0,0 +1,21 @@
+/*--------------------------------*- 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 "constant";
+ object dynamicMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dynamicFvMesh staticFvMesh;
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/standingWave/constant/g b/tutorials/multiphase/interIsoFoam/standingWave/constant/g
new file mode 100644
index 0000000000000000000000000000000000000000..195dbea3f3cac083e26f182b8ae44f0074c8486c
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/standingWave/constant/g
@@ -0,0 +1,22 @@
+/*--------------------------------*- 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 uniformDimensionedVectorField;
+ location "constant";
+ object g;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions [0 1 -2 0 0 0 0];
+value ( 0 0 -9.81 );
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/standingWave/constant/transportProperties b/tutorials/multiphase/interIsoFoam/standingWave/constant/transportProperties
new file mode 100644
index 0000000000000000000000000000000000000000..f7055c5a380fdefb3602e3d1a1da4e72dd28a81f
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/standingWave/constant/transportProperties
@@ -0,0 +1,37 @@
+/*--------------------------------*- 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 "constant";
+ object transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+phases (water air);
+
+water
+{
+ transportModel Newtonian;
+ nu [0 2 -1 0 0 0 0] 1e-06;
+ rho [1 -3 0 0 0 0 0] 1000;
+}
+
+air
+{
+ transportModel Newtonian;
+ nu [0 2 -1 0 0 0 0] 1.48e-05;
+ rho [1 -3 0 0 0 0 0] 1;
+}
+
+sigma sigma [ 1 0 -2 0 0 0 0 ] 0.00;
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/standingWave/constant/turbulenceProperties b/tutorials/multiphase/interIsoFoam/standingWave/constant/turbulenceProperties
new file mode 100644
index 0000000000000000000000000000000000000000..5eec04267266e7fd15e7701a875d683d5e658cd9
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/standingWave/constant/turbulenceProperties
@@ -0,0 +1,21 @@
+/*--------------------------------*- 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 "constant";
+ object turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType laminar;
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/standingWave/system/blockMeshDict b/tutorials/multiphase/interIsoFoam/standingWave/system/blockMeshDict
new file mode 100644
index 0000000000000000000000000000000000000000..0330ca8d5fb083923f95f5dee83d24318eddf06f
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/standingWave/system/blockMeshDict
@@ -0,0 +1,106 @@
+/*--------------------------------*- 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;
+ object blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 1;
+
+L 1;
+nx 50;
+
+y1 -.01;
+y2 .01;
+ny 1;
+
+H 1;
+nz 50;
+
+vertices
+(
+ (0 $y1 0)
+ ($L $y1 0)
+ ($L $y2 0)
+ (0 $y2 0)
+ (0 $y1 $H)
+ ($L $y1 $H)
+ ($L $y2 $H)
+ (0 $y2 $H)
+);
+
+blocks
+(
+ hex (0 1 2 3 4 5 6 7) ($nx $ny $nz) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+ left
+ {
+ type patch;
+ faces
+ (
+ (0 4 7 3)
+ );
+ }
+ right
+ {
+ type patch;
+ faces
+ (
+ (1 2 6 5)
+ );
+ }
+ top
+ {
+ type patch;
+ faces
+ (
+ (4 5 6 7)
+ );
+ }
+ bottom
+ {
+ type patch;
+ faces
+ (
+ (0 3 2 1)
+ );
+ }
+ front
+ {
+ type empty;
+ faces
+ (
+ (0 1 5 4)
+ );
+ }
+ back
+ {
+ type empty;
+ faces
+ (
+ (2 3 7 6)
+ );
+ }
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/standingWave/system/controlDict b/tutorials/multiphase/interIsoFoam/standingWave/system/controlDict
new file mode 100644
index 0000000000000000000000000000000000000000..152a0075faff6a04ed6e9c2cc7daf984c8899905
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/standingWave/system/controlDict
@@ -0,0 +1,56 @@
+/*--------------------------------*- 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 controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application interIsoFoam;
+
+startFrom startTime;
+
+startTime 0;
+
+stopAt endTime;
+
+endTime 10;
+
+deltaT 0.001;
+
+writeControl adjustableRunTime;
+
+writeInterval .1;
+
+purgeWrite 0;
+
+writeFormat ascii;
+
+writePrecision 6;
+
+writeCompression uncompressed;
+
+timeFormat general;
+
+timePrecision 6;
+
+runTimeModifiable yes;
+
+adjustTimeStep yes;
+
+maxCo 0.5;
+maxAlphaCo 0.5;
+
+maxDeltaT 1;
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/standingWave/system/decomposeParDict b/tutorials/multiphase/interIsoFoam/standingWave/system/decomposeParDict
new file mode 100644
index 0000000000000000000000000000000000000000..28839d76e1c3259378a8a35330871cb3f7c06947
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/standingWave/system/decomposeParDict
@@ -0,0 +1,45 @@
+/*--------------------------------*- 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 decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains 1;
+
+method simple;
+
+simpleCoeffs
+{
+ n ( 2 1 2 );
+ delta 0.001;
+}
+
+hierarchicalCoeffs
+{
+ n ( 1 1 1 );
+ delta 0.001;
+ order xyz;
+}
+
+manualCoeffs
+{
+ dataFile "";
+}
+
+distributed no;
+
+roots ( );
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/standingWave/system/fvSchemes b/tutorials/multiphase/interIsoFoam/standingWave/system/fvSchemes
new file mode 100644
index 0000000000000000000000000000000000000000..afd820abeaaacf3b759ab784af5561b7ecf276f2
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/standingWave/system/fvSchemes
@@ -0,0 +1,59 @@
+/*--------------------------------*- 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 fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+ default CrankNicolson 0.5;
+}
+
+gradSchemes
+{
+ default Gauss linear;
+}
+
+divSchemes
+{
+ div(rhoPhi,U) Gauss limitedLinearV 1;
+ div(phi,alpha) Gauss vanLeer;
+ div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
+}
+
+laplacianSchemes
+{
+ default Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+ default linear;
+}
+
+snGradSchemes
+{
+ default corrected;
+}
+
+fluxRequired
+{
+ default no;
+ p_rgh;
+ pcorr;
+ alpha.water;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/standingWave/system/fvSolution b/tutorials/multiphase/interIsoFoam/standingWave/system/fvSolution
new file mode 100644
index 0000000000000000000000000000000000000000..993e09cc7cadf7c4c14f3cabe8e75bdef1eba3ba
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/standingWave/system/fvSolution
@@ -0,0 +1,75 @@
+/*--------------------------------*- 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;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+
+ "alpha.*"
+ {
+ isoFaceTol 1e-6;
+ surfCellTol 1e-6;
+ nAlphaBounds 2;
+ snapTol 1e-12;
+ clip true;
+
+ nAlphaSubCycles 1;
+ cAlpha 1;
+ }
+
+ "pcorr.*"
+ {
+ solver PCG;
+ preconditioner DIC;
+ tolerance 1e-10;
+ relTol 0;
+ }
+
+ p_rgh
+ {
+ solver PCG;
+ preconditioner DIC;
+ tolerance 1e-07;
+ relTol 0.05;
+ }
+
+ p_rghFinal
+ {
+ $p_rgh;
+ tolerance 1e-07;
+ relTol 0;
+ }
+
+ U
+ {
+ solver PBiCG;
+ preconditioner DILU;
+ tolerance 1e-06;
+ relTol 0;
+ }
+}
+
+PIMPLE
+{
+ momentumPredictor no;
+ nCorrectors 3;
+ nOuterCorrectors 1;
+ nNonOrthogonalCorrectors 0;
+ pRefCell 0;
+ pRefValue 0;
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/standingWave/system/refineMeshDict1 b/tutorials/multiphase/interIsoFoam/standingWave/system/refineMeshDict1
new file mode 100644
index 0000000000000000000000000000000000000000..e3efa074d9f446aeb98ad981687b55b89e1b9744
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/standingWave/system/refineMeshDict1
@@ -0,0 +1,62 @@
+/*--------------------------------*- 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;
+ object refineMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Cells to refine; name of cell set
+set c1;
+
+// Type of coordinate system:
+// - global : coordinate system same for every cell. Usually aligned with
+// x,y,z axis. Specify in globalCoeffs section below.
+// - patchLocal : coordinate system different for every cell. Specify in
+// patchLocalCoeffs section below.
+coordinateSystem global;
+//coordinateSystem patchLocal;
+
+// .. and its coefficients. x,y in this case. (normal direction is calculated
+// as tan1^tan2)
+globalCoeffs
+{
+ tan1 (1 0 0);
+ tan2 (0 1 0);
+}
+
+patchLocalCoeffs
+{
+ patch outside; // Normal direction is facenormal of zero'th face of patch
+ tan1 (1 0 0);
+}
+
+// List of directions to refine
+directions
+(
+ tan1
+// tan2
+ normal
+);
+
+// Whether to use hex topology. This will
+// - if patchLocal: all cells on selected patch should be hex
+// - split all hexes in 2x2x2 through the middle of edges.
+useHexTopology true;
+
+// Cut purely geometric (will cut hexes through vertices) or take topology
+// into account. Incompatible with useHexTopology
+geometricCut false;
+
+// Write meshes from intermediate steps
+writeMesh false;
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/standingWave/system/refineMeshDict2 b/tutorials/multiphase/interIsoFoam/standingWave/system/refineMeshDict2
new file mode 100644
index 0000000000000000000000000000000000000000..eaf14b04ed670d05661e22597158426414cbf8d0
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/standingWave/system/refineMeshDict2
@@ -0,0 +1,62 @@
+/*--------------------------------*- 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;
+ object refineMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Cells to refine; name of cell set
+set c2;
+
+// Type of coordinate system:
+// - global : coordinate system same for every cell. Usually aligned with
+// x,y,z axis. Specify in globalCoeffs section below.
+// - patchLocal : coordinate system different for every cell. Specify in
+// patchLocalCoeffs section below.
+coordinateSystem global;
+//coordinateSystem patchLocal;
+
+// .. and its coefficients. x,y in this case. (normal direction is calculated
+// as tan1^tan2)
+globalCoeffs
+{
+ tan1 (1 0 0);
+ tan2 (0 1 0);
+}
+
+patchLocalCoeffs
+{
+ patch outside; // Normal direction is facenormal of zero'th face of patch
+ tan1 (1 0 0);
+}
+
+// List of directions to refine
+directions
+(
+ tan1
+// tan2
+ normal
+);
+
+// Whether to use hex topology. This will
+// - if patchLocal: all cells on selected patch should be hex
+// - split all hexes in 2x2x2 through the middle of edges.
+useHexTopology true;
+
+// Cut purely geometric (will cut hexes through vertices) or take topology
+// into account. Incompatible with useHexTopology
+geometricCut false;
+
+// Write meshes from intermediate steps
+writeMesh false;
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/standingWave/system/setAlphaFieldDict b/tutorials/multiphase/interIsoFoam/standingWave/system/setAlphaFieldDict
new file mode 100644
index 0000000000000000000000000000000000000000..a500bbd2be3c8d94e1e77fa2f42bf96971d3dfe7
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/standingWave/system/setAlphaFieldDict
@@ -0,0 +1,26 @@
+/*--------------------------------*- 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 sin;
+direction (1 0 0);
+up (0 0 1);
+centre (-0.5 0 0.5);
+period 2;
+amplitude 0.05;
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/standingWave/system/topoSetDict1 b/tutorials/multiphase/interIsoFoam/standingWave/system/topoSetDict1
new file mode 100644
index 0000000000000000000000000000000000000000..d6ac6f57be3c775531ab3414b69733e87233ede5
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/standingWave/system/topoSetDict1
@@ -0,0 +1,32 @@
+/*--------------------------------*- 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;
+ object topoSetDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+actions
+(
+ {
+ name c1;
+ type cellSet;
+ action new;
+ source boxToCell;
+ sourceInfo
+ {
+ box (-1e10 -1e10 0.3) (1e10 1e10 0.7);
+ }
+ }
+);
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interIsoFoam/standingWave/system/topoSetDict2 b/tutorials/multiphase/interIsoFoam/standingWave/system/topoSetDict2
new file mode 100644
index 0000000000000000000000000000000000000000..1971415ba2270b8ff089a3e80e72f1946706e235
--- /dev/null
+++ b/tutorials/multiphase/interIsoFoam/standingWave/system/topoSetDict2
@@ -0,0 +1,32 @@
+/*--------------------------------*- 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;
+ object topoSetDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+actions
+(
+ {
+ name c2;
+ type cellSet;
+ action new;
+ source boxToCell;
+ sourceInfo
+ {
+ box (-1e10 -1e10 0.4) (1e10 1e10 0.6);
+ }
+ }
+);
+
+// ************************************************************************* //