Merge branch 'master' of /home/dm4/OpenFOAM/OpenFOAM-dev

applications/solvers/multiphase/compressibleMultiphaseInterFoam/Allwclean

@@ -0,0 +1,8 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+set -x
+
+wclean libso multiphaseMixtureThermo
+wclean
+
+# ----------------------------------------------------------------- end-of-file

applications/solvers/multiphase/compressibleMultiphaseInterFoam/Allwmake

@@ -0,0 +1,8 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+set -x
+
+wmake libso multiphaseMixtureThermo
+wmake
+
+# ----------------------------------------------------------------- end-of-file

applications/solvers/multiphase/compressibleMultiphaseInterFoam/Make/files

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

applications/solvers/multiphase/compressibleMultiphaseInterFoam/Make/options

@@ -0,0 +1,17 @@
+EXE_INC = \
+    -ImultiphaseMixtureThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
+    -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS = \
+    -lmultiphaseMixtureThermo \
+    -lfluidThermophysicalModels \
+    -lspecie \
+    -linterfaceProperties \
+    -lcompressibleTurbulenceModel \
+    -lcompressibleRASModels \
+    -lcompressibleLESModels \
+    -lfiniteVolume

applications/solvers/multiphase/compressibleMultiphaseInterFoam/TEqn.H

@@ -0,0 +1,17 @@
+{
+    fvScalarMatrix TEqn
+    (
+        fvm::ddt(rho, T)
+      + fvm::div(multiphaseProperties.rhoPhi(), T)
+      - fvm::laplacian(multiphaseProperties.alphaEff(turbulence->mut()), T)
+      + (
+            fvc::div(fvc::absolute(phi, U), p)
+          + fvc::ddt(rho, K) + fvc::div(multiphaseProperties.rhoPhi(), K)
+        )*multiphaseProperties.rCv()
+    );
+
+    TEqn.relax();
+    TEqn.solve();
+
+    multiphaseProperties.correct();
+}

applications/solvers/multiphase/compressibleMultiphaseInterFoam/UEqn.H

@@ -0,0 +1,27 @@
+    fvVectorMatrix UEqn
+    (
+        fvm::ddt(rho, U)
+      + fvm::div(multiphaseProperties.rhoPhi(), U)
+      + turbulence->divDevRhoReff(U)
+    );
+
+    UEqn.relax();
+
+    if (pimple.momentumPredictor())
+    {
+        solve
+        (
+            UEqn
+         ==
+            fvc::reconstruct
+            (
+                (
+                    multiphaseProperties.surfaceTensionForce()
+                  - ghf*fvc::snGrad(rho)
+                  - fvc::snGrad(p_rgh)
+                ) * mesh.magSf()
+            )
+        );
+
+        K = 0.5*magSqr(U);
+    }

applications/solvers/multiphase/compressibleMultiphaseInterFoam/alphaCourantNo.H

@@ -0,0 +1,57 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 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
+    CourantNo
+
+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
+    (
+        multiphaseProperties.nearInterface()().internalField()
+      * fvc::surfaceSum(mag(phi))().internalField()
+    );
+
+    alphaCoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
+
+    meanAlphaCoNum =
+        0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();
+}
+
+Info<< "Interface Courant Number mean: " << meanAlphaCoNum
+    << " max: " << alphaCoNum << endl;
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseInterFoam.C

@@ -0,0 +1,111 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 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/>.
+
+Application
+    compressibleMultiphaseInterFoam
+
+Description
+    Solver for n compressible, non-isothermal immiscible fluids using a VOF
+    (volume of fluid) 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.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "multiphaseMixtureThermo.H"
+#include "turbulenceModel.H"
+#include "pimpleControl.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readGravitationalAcceleration.H"
+
+    pimpleControl pimple(mesh);
+
+    #include "readTimeControls.H"
+    #include "initContinuityErrs.H"
+    #include "createFields.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())
+        {
+            multiphaseProperties.solve();
+
+            solve(fvm::ddt(rho) + fvc::div(multiphaseProperties.rhoPhi()));
+
+            #include "UEqn.H"
+            #include "TEqn.H"
+
+            // --- Pressure corrector loop
+            while (pimple.correct())
+            {
+                #include "pEqn.H"
+            }
+
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleMultiphaseInterFoam/createFields.H

@@ -0,0 +1,71 @@
+    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<< "Constructing multiphaseMixtureThermo\n" << endl;
+    multiphaseMixtureThermo multiphaseProperties(U, phi);
+
+    volScalarField& p = multiphaseProperties.p();
+    volScalarField& T = multiphaseProperties.T();
+
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT
+        ),
+        multiphaseProperties.rho()
+    );
+
+    Info<< max(rho) << min(rho);
+
+    dimensionedScalar pMin(multiphaseProperties.lookup("pMin"));
+
+
+    Info<< "Calculating field g.h\n" << endl;
+    volScalarField gh("gh", g & mesh.C());
+    surfaceScalarField ghf("ghf", g & mesh.Cf());
+
+    // Construct compressible turbulence model
+    autoPtr<compressible::turbulenceModel> turbulence
+    (
+        compressible::turbulenceModel::New
+        (
+            rho,
+            U,
+            multiphaseProperties.rhoPhi(),
+            multiphaseProperties
+        )
+    );
+
+    Info<< "Creating field kinetic energy K\n" << endl;
+    volScalarField K("K", 0.5*magSqr(U));

applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/Make/files

@@ -0,0 +1,5 @@
+phaseModel/phaseModel.C
+alphaContactAngle/alphaContactAngleFvPatchScalarField.C
+multiphaseMixtureThermo.C
+
+LIB = $(FOAM_LIBBIN)/libmultiphaseMixtureThermo

applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/Make/options

@@ -0,0 +1,8 @@
+EXE_INC = \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+LIB_LIBS = \
+    -lfluidThermophysicalModels \
+    -lspecie \
+    -lfiniteVolume

applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/alphaContactAngle/alphaContactAngleFvPatchScalarField.C

@@ -0,0 +1,146 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "alphaContactAngleFvPatchScalarField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvPatchFieldMapper.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+alphaContactAngleFvPatchScalarField::interfaceThetaProps::interfaceThetaProps
+(
+    Istream& is
+)
+:
+    theta0_(readScalar(is)),
+    uTheta_(readScalar(is)),
+    thetaA_(readScalar(is)),
+    thetaR_(readScalar(is))
+{}
+
+
+Istream& operator>>
+(
+    Istream& is,
+    alphaContactAngleFvPatchScalarField::interfaceThetaProps& tp
+)
+{
+    is >> tp.theta0_ >> tp.uTheta_ >> tp.thetaA_ >> tp.thetaR_;
+    return is;
+}
+
+
+Ostream& operator<<
+(
+    Ostream& os,
+    const alphaContactAngleFvPatchScalarField::interfaceThetaProps& tp
+)
+{
+    os  << tp.theta0_ << token::SPACE
+        << tp.uTheta_ << token::SPACE
+        << tp.thetaA_ << token::SPACE
+        << tp.thetaR_;
+
+    return os;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    zeroGradientFvPatchScalarField(p, iF)
+{}
+
+
+alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField
+(
+    const alphaContactAngleFvPatchScalarField& gcpsf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    zeroGradientFvPatchScalarField(gcpsf, p, iF, mapper),
+    thetaProps_(gcpsf.thetaProps_)
+{}
+
+
+alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    zeroGradientFvPatchScalarField(p, iF),
+    thetaProps_(dict.lookup("thetaProperties"))
+{
+    evaluate();
+}
+
+
+alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField
+(
+    const alphaContactAngleFvPatchScalarField& gcpsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    zeroGradientFvPatchScalarField(gcpsf, iF),
+    thetaProps_(gcpsf.thetaProps_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void alphaContactAngleFvPatchScalarField::write(Ostream& os) const
+{
+    fvPatchScalarField::write(os);
+    os.writeKeyword("thetaProperties")
+        << thetaProps_ << token::END_STATEMENT << nl;
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeField
+(
+    fvPatchScalarField,
+    alphaContactAngleFvPatchScalarField
+);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/alphaContactAngle/alphaContactAngleFvPatchScalarField.H

@@ -0,0 +1,215 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 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/>.
+
+Class
+    Foam::alphaContactAngleFvPatchScalarField
+
+Description
+    Contact-angle boundary condition for multi-phase interface-capturing
+    simulations.  Used in conjuction with multiphaseMixture.
+
+SourceFiles
+    alphaContactAngleFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef alphaContactAngleFvPatchScalarField_H
+#define alphaContactAngleFvPatchScalarField_H
+
+#include "zeroGradientFvPatchFields.H"
+#include "multiphaseMixtureThermo.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class alphaContactAngleFvPatch Declaration
+\*---------------------------------------------------------------------------*/
+
+class alphaContactAngleFvPatchScalarField
+:
+    public zeroGradientFvPatchScalarField
+{
+public:
+
+    class interfaceThetaProps
+    {
+        //- Equilibrium contact angle
+        scalar theta0_;
+
+        //- Dynamic contact angle velocity scale
+        scalar uTheta_;
+
+        //- Limiting advancing contact angle
+        scalar thetaA_;
+
+        //- Limiting receeding contact angle
+        scalar thetaR_;
+
+
+    public:
+
+        // Constructors
+            interfaceThetaProps()
+            {}
+
+            interfaceThetaProps(Istream&);
+
+
+        // Member functions
+
+            //- Return the equilibrium contact angle theta0
+            scalar theta0(bool matched=true) const
+            {
+                if (matched) return theta0_;
+                else return 180.0 - theta0_;
+            }
+
+            //- Return the dynamic contact angle velocity scale
+            scalar uTheta() const
+            {
+                return uTheta_;
+            }
+
+            //- Return the limiting advancing contact angle
+            scalar thetaA(bool matched=true) const
+            {
+                if (matched) return thetaA_;
+                else return 180.0 - thetaA_;
+            }
+
+            //- Return the limiting receeding contact angle
+            scalar thetaR(bool matched=true) const
+            {
+                if (matched) return thetaR_;
+                else return 180.0 - thetaR_;
+            }
+
+
+        // IO functions
+
+            friend Istream& operator>>(Istream&, interfaceThetaProps&);
+            friend Ostream& operator<<(Ostream&, const interfaceThetaProps&);
+    };
+
+    typedef HashTable
+    <
+        interfaceThetaProps,
+        multiphaseMixtureThermo::interfacePair,
+        multiphaseMixtureThermo::interfacePair::hash
+    > thetaPropsTable;
+
+
+private:
+
+    // Private data
+
+        thetaPropsTable thetaProps_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("alphaContactAngle");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        alphaContactAngleFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        alphaContactAngleFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given alphaContactAngleFvPatchScalarField
+        //  onto a new patch
+        alphaContactAngleFvPatchScalarField
+        (
+            const alphaContactAngleFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new alphaContactAngleFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        alphaContactAngleFvPatchScalarField
+        (
+            const alphaContactAngleFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new alphaContactAngleFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        //- Return the contact angle properties
+        const thetaPropsTable& thetaProps() const
+        {
+            return thetaProps_;
+        }
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.H

@@ -0,0 +1,434 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 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/>.
+
+Class
+    Foam::multiphaseMixtureThermo
+
+Description
+
+SourceFiles
+    multiphaseMixtureThermo.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef multiphaseMixtureThermo_H
+#define multiphaseMixtureThermo_H
+
+#include "phaseModel.H"
+#include "PtrDictionary.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+#include "rhoThermo.H"
+#include "psiThermo.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class multiphaseMixtureThermo Declaration
+\*---------------------------------------------------------------------------*/
+
+class multiphaseMixtureThermo
+:
+    public psiThermo
+{
+public:
+
+    class interfacePair
+    :
+        public Pair<word>
+    {
+    public:
+
+        class hash
+        :
+            public Hash<interfacePair>
+        {
+        public:
+
+            hash()
+            {}
+
+            label operator()(const interfacePair& key) const
+            {
+                return word::hash()(key.first()) + word::hash()(key.second());
+            }
+        };
+
+
+        // Constructors
+
+            interfacePair()
+            {}
+
+            interfacePair(const word& alpha1Name, const word& alpha2Name)
+            :
+                Pair<word>(alpha1Name, alpha2Name)
+            {}
+
+            interfacePair(const phaseModel& alpha1, const phaseModel& alpha2)
+            :
+                Pair<word>(alpha1.name(), alpha2.name())
+            {}
+
+
+        // Friend Operators
+
+            friend bool operator==
+            (
+                const interfacePair& a,
+                const interfacePair& b
+            )
+            {
+                return
+                (
+                    ((a.first() == b.first()) && (a.second() == b.second()))
+                 || ((a.first() == b.second()) && (a.second() == b.first()))
+                );
+            }
+
+            friend bool operator!=
+            (
+                const interfacePair& a,
+                const interfacePair& b
+            )
+            {
+                return (!(a == b));
+            }
+    };
+
+
+private:
+
+    // Private data
+
+        //- Dictionary of phases
+        PtrDictionary<phaseModel> phases_;
+
+        const fvMesh& mesh_;
+        const volVectorField& U_;
+        const surfaceScalarField& phi_;
+
+        surfaceScalarField rhoPhi_;
+
+        volScalarField alphas_;
+
+        typedef HashTable<scalar, interfacePair, interfacePair::hash>
+            sigmaTable;
+
+        sigmaTable sigmas_;
+        dimensionSet dimSigma_;
+
+        //- Stabilisation for normalisation of the interface normal
+        const dimensionedScalar deltaN_;
+
+        //- Conversion factor for degrees into radians
+        static const scalar convertToRad;
+
+
+    // Private member functions
+
+        void calcAlphas();
+
+        void solveAlphas(const scalar cAlpha);
+
+        tmp<surfaceVectorField> nHatfv
+        (
+            const volScalarField& alpha1,
+            const volScalarField& alpha2
+        ) const;
+
+        tmp<surfaceScalarField> nHatf
+        (
+            const volScalarField& alpha1,
+            const volScalarField& alpha2
+        ) const;
+
+        void correctContactAngle
+        (
+            const phaseModel& alpha1,
+            const phaseModel& alpha2,
+            surfaceVectorField::GeometricBoundaryField& nHatb
+        ) const;
+
+        tmp<volScalarField> K
+        (
+            const phaseModel& alpha1,
+            const phaseModel& alpha2
+        ) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("multiphaseMixtureThermo");
+
+
+    // Constructors
+
+        //- Construct from components
+        multiphaseMixtureThermo
+        (
+            const volVectorField& U,
+            const surfaceScalarField& phi
+        );
+
+
+    //- Destructor
+    virtual ~multiphaseMixtureThermo()
+    {}
+
+
+    // Member Functions
+
+        //- Return the phases
+        const PtrDictionary<phaseModel>& phases() const
+        {
+            return phases_;
+        }
+
+        //- Return non-const access to the phases
+        PtrDictionary<phaseModel>& phases()
+        {
+            return phases_;
+        }
+
+        //- Return the velocity
+        const volVectorField& U() const
+        {
+            return U_;
+        }
+
+        //- Return the volumetric flux
+        const surfaceScalarField& phi() const
+        {
+            return phi_;
+        }
+
+        const surfaceScalarField& rhoPhi() const
+        {
+            return rhoPhi_;
+        }
+
+        //- Update properties
+        virtual void correct();
+
+        //- Update densities for given pressure change
+        void correctRho(const volScalarField& dp);
+
+        //- Return true if the equation of state is incompressible
+        //  i.e. rho != f(p)
+        virtual bool incompressible() const;
+
+        //- Return true if the equation of state is isochoric
+        //  i.e. rho = const
+        virtual bool isochoric() const;
+
+
+        // Access to thermodynamic state variables
+
+            //- Enthalpy/Internal energy [J/kg]
+            //  Non-const access allowed for transport equations
+            virtual volScalarField& he()
+            {
+                notImplemented("multiphaseMixtureThermo::he()");
+                return phases_[0]->thermo().he();
+            }
+
+            //- Enthalpy/Internal energy [J/kg]
+            virtual const volScalarField& he() const
+            {
+                notImplemented("multiphaseMixtureThermo::he() const");
+                return phases_[0]->thermo().he();
+            }
+
+            //- Enthalpy/Internal energy
+            //  for given pressure and temperature [J/kg]
+            virtual tmp<volScalarField> he
+            (
+                const volScalarField& p,
+                const volScalarField& T
+            ) const;
+
+            //- Enthalpy/Internal energy for cell-set [J/kg]
+            virtual tmp<scalarField> he
+            (
+                const scalarField& p,
+                const scalarField& T,
+                const labelList& cells
+            ) const;
+
+            //- Enthalpy/Internal energy for patch [J/kg]
+            virtual tmp<scalarField> he
+            (
+                const scalarField& p,
+                const scalarField& T,
+                const label patchi
+            ) const;
+
+            //- Chemical enthalpy [J/kg]
+            virtual tmp<volScalarField> hc() const;
+
+            //- Temperature from enthalpy/internal energy for cell-set
+            virtual tmp<scalarField> THE
+            (
+                const scalarField& h,
+                const scalarField& p,
+                const scalarField& T0,      // starting temperature
+                const labelList& cells
+            ) const;
+
+            //- Temperature from enthalpy/internal energy for patch
+            virtual tmp<scalarField> THE
+            (
+                const scalarField& h,
+                const scalarField& p,
+                const scalarField& T0,      // starting temperature
+                const label patchi
+            ) const;
+
+
+        // Fields derived from thermodynamic state variables
+
+            //- Density [kg/m^3]
+            virtual tmp<volScalarField> rho() const;
+
+            //- Heat capacity at constant pressure [J/kg/K]
+            virtual tmp<volScalarField> Cp() const;
+
+            //- Heat capacity at constant pressure for patch [J/kg/K]
+            virtual tmp<scalarField> Cp
+            (
+                const scalarField& p,
+                const scalarField& T,
+                const label patchi
+            ) const;
+
+            //- Heat capacity at constant volume [J/kg/K]
+            virtual tmp<volScalarField> Cv() const;
+
+            //- Heat capacity at constant volume for patch [J/kg/K]
+            virtual tmp<scalarField> Cv
+            (
+                const scalarField& p,
+                const scalarField& T,
+                const label patchi
+            ) const;
+
+            //- gamma = Cp/Cv []
+            virtual tmp<volScalarField> gamma() const;
+
+            //- gamma = Cp/Cv for patch []
+            virtual tmp<scalarField> gamma
+            (
+                const scalarField& p,
+                const scalarField& T,
+                const label patchi
+            ) const;
+
+            //- Heat capacity at constant pressure/volume [J/kg/K]
+            virtual tmp<volScalarField> Cpv() const;
+
+            //- Heat capacity at constant pressure/volume for patch [J/kg/K]
+            virtual tmp<scalarField> Cpv
+            (
+                const scalarField& p,
+                const scalarField& T,
+                const label patchi
+            ) const;
+
+            //- Heat capacity ratio []
+            virtual tmp<volScalarField> CpByCpv() const;
+
+            //- Heat capacity ratio for patch []
+            virtual tmp<scalarField> CpByCpv
+            (
+                const scalarField& p,
+                const scalarField& T,
+                const label patchi
+            ) const;
+
+
+        // Fields derived from transport state variables
+
+            //- Thermal diffusivity for temperature of mixture [J/m/s/K]
+            virtual tmp<volScalarField> kappa() const;
+
+            //- Thermal diffusivity of mixture for patch [J/m/s/K]
+            virtual tmp<scalarField> kappa
+            (
+                const label patchi
+            ) const;
+
+            //- Effective thermal diffusivity of mixture [J/m/s/K]
+            virtual tmp<volScalarField> kappaEff
+            (
+                const volScalarField& alphat
+            ) const;
+
+            //- Effective thermal diffusivity of mixture for patch [J/m/s/K]
+            virtual tmp<scalarField> kappaEff
+            (
+                const scalarField& alphat,
+                const label patchi
+            ) const;
+
+            //- Effective thermal diffusivity of mixture [J/m/s/K]
+            virtual tmp<volScalarField> alphaEff
+            (
+                const volScalarField& alphat
+            ) const;
+
+            //- Effective thermal diffusivity of mixture for patch [J/m/s/K]
+            virtual tmp<scalarField> alphaEff
+            (
+                const scalarField& alphat,
+                const label patchi
+            ) const;
+
+
+        //- Return the phase-averaged reciprocal Cv
+        tmp<volScalarField> rCv() const;
+
+        tmp<surfaceScalarField> surfaceTensionForce() const;
+
+        //- Indicator of the proximity of the interface
+        //  Field values are 1 near and 0 away for the interface.
+        tmp<volScalarField> nearInterface() const;
+
+        //- Solve for the mixture phase-fractions
+        void solve();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/phaseModel/phaseModel.C

@@ -0,0 +1,95 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "phaseModel.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::phaseModel::phaseModel
+(
+    const word& phaseName,
+    const volScalarField& p,
+    const volScalarField& T
+)
+:
+    volScalarField
+    (
+        IOobject
+        (
+            IOobject::groupName("alpha", phaseName),
+            p.mesh().time().timeName(),
+            p.mesh(),
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        p.mesh()
+    ),
+    name_(phaseName),
+    p_(p),
+    T_(T),
+    thermo_(NULL),
+    dgdt_
+    (
+        IOobject
+        (
+            IOobject::groupName("dgdt", phaseName),
+            p.mesh().time().timeName(),
+            p.mesh(),
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        p.mesh(),
+        dimensionedScalar("0", dimless/dimTime, 0)
+    )
+{
+    {
+        volScalarField Tp(IOobject::groupName("T", phaseName), T);
+        Tp.write();
+    }
+
+    thermo_ = rhoThermo::New(p.mesh(), phaseName);
+    thermo_->validate(phaseName, "e");
+
+    correct();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::phaseModel> Foam::phaseModel::clone() const
+{
+    notImplemented("phaseModel::clone() const");
+    return autoPtr<phaseModel>(NULL);
+}
+
+
+void Foam::phaseModel::correct()
+{
+    thermo_->he() = thermo_->he(p_, T_);
+    thermo_->correct();
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/phaseModel/phaseModel.H

@@ -0,0 +1,156 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 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/>.
+
+Class
+    Foam::phaseModel
+
+Description
+    Single incompressible phase derived from the phase-fraction.
+    Used as part of the multiPhaseMixture for interface-capturing multi-phase
+    simulations.
+
+SourceFiles
+    phaseModel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef phaseModel_H
+#define phaseModel_H
+
+#include "rhoThermo.H"
+#include "volFields.H"
+#include "dictionaryEntry.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class phaseModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class phaseModel
+:
+    public volScalarField
+{
+    // Private data
+
+        word name_;
+        const volScalarField& p_;
+        const volScalarField& T_;
+        autoPtr<rhoThermo> thermo_;
+        volScalarField dgdt_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        phaseModel
+        (
+            const word& phaseName,
+            const volScalarField& p,
+            const volScalarField& T
+        );
+
+        //- Return clone
+        autoPtr<phaseModel> clone() const;
+
+        //- Return a pointer to a new phaseModel created on freestore
+        //  from Istream
+        class iNew
+        {
+            const volScalarField& p_;
+            const volScalarField& T_;
+
+        public:
+
+            iNew
+            (
+                const volScalarField& p,
+                const volScalarField& T
+            )
+            :
+                p_(p),
+                T_(T)
+            {}
+
+            autoPtr<phaseModel> operator()(Istream& is) const
+            {
+                return autoPtr<phaseModel>(new phaseModel(is, p_, T_));
+            }
+        };
+
+
+    // Member Functions
+
+        const word& name() const
+        {
+            return name_;
+        }
+
+        const word& keyword() const
+        {
+            return name();
+        }
+
+        //- Return const-access to phase rhoThermo
+        const rhoThermo& thermo() const
+        {
+            return thermo_();
+        }
+
+        //- Return access to phase rhoThermo
+        rhoThermo& thermo()
+        {
+            return thermo_();
+        }
+
+        //- Return const-access to phase divergence
+        const volScalarField& dgdt() const
+        {
+            return dgdt_;
+        }
+
+        //- Return access to phase divergence
+        volScalarField& dgdt()
+        {
+            return dgdt_;
+        }
+
+        void correct();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H

@@ -0,0 +1,142 @@
+{
+    volScalarField rAU("rAU", 1.0/UEqn.A());
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
+
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();
+
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        (fvc::interpolate(HbyA) & mesh.Sf())
+      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
+    );
+
+    surfaceScalarField phig
+    (
+        (
+            multiphaseProperties.surfaceTensionForce()
+          - ghf*fvc::snGrad(rho)
+        )*rAUf*mesh.magSf()
+    );
+
+    phiHbyA += phig;
+
+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    (
+        p_rgh.boundaryField(),
+        (
+            phiHbyA.boundaryField()
+          - (mesh.Sf().boundaryField() & U.boundaryField())
+        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+    );
+
+    PtrList<fvScalarMatrix> p_rghEqnComps(multiphaseProperties.phases().size());
+
+    label phasei = 0;
+    forAllConstIter
+    (
+        PtrDictionary<phaseModel>,
+        multiphaseProperties.phases(),
+        phase
+    )
+    {
+        const rhoThermo& thermo = phase().thermo();
+        const volScalarField& rho = thermo.rho()();
+
+        p_rghEqnComps.set
+        (
+            phasei,
+            (
+                fvc::ddt(rho) + thermo.psi()*correction(fvm::ddt(p_rgh))
+              + fvc::div(phi, rho) - fvc::Sp(fvc::div(phi), rho)
+            ).ptr()
+        );
+
+        phasei++;
+    }
+
+    // Cache p_rgh prior to solve for density update
+    volScalarField p_rgh_0(p_rgh);
+
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix p_rghEqnIncomp
+        (
+            fvc::div(phiHbyA)
+          - fvm::laplacian(rAUf, p_rgh)
+        );
+
+        tmp<fvScalarMatrix> p_rghEqnComp;
+
+        phasei = 0;
+        forAllConstIter
+        (
+            PtrDictionary<phaseModel>,
+            multiphaseProperties.phases(),
+            phase
+        )
+        {
+            tmp<fvScalarMatrix> hmm
+            (
+                (max(phase(), scalar(0))/phase().thermo().rho())
+               *p_rghEqnComps[phasei]
+            );
+
+            if (phasei == 0)
+            {
+                p_rghEqnComp = hmm;
+            }
+            else
+            {
+                p_rghEqnComp() += hmm;
+            }
+
+            phasei++;
+        }
+
+        solve
+        (
+            p_rghEqnComp
+          + p_rghEqnIncomp,
+            mesh.solver(p_rgh.select(pimple.finalInnerIter()))
+        );
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            // p = max(p_rgh + multiphaseProperties.rho()*gh, pMin);
+            // p_rgh = p - multiphaseProperties.rho()*gh;
+
+            phasei = 0;
+            forAllIter
+            (
+                PtrDictionary<phaseModel>,
+                multiphaseProperties.phases(),
+                phase
+            )
+            {
+                phase().dgdt() =
+                    pos(phase())
+                  *(p_rghEqnComps[phasei] & p_rgh)/phase().thermo().rho();
+            }
+
+            phi = phiHbyA + p_rghEqnIncomp.flux();
+
+            U = HbyA
+              + rAU*fvc::reconstruct((phig + p_rghEqnIncomp.flux())/rAUf);
+            U.correctBoundaryConditions();
+        }
+    }
+
+    p = max(p_rgh + multiphaseProperties.rho()*gh, pMin);
+
+    // Update densities from change in p_rgh
+    multiphaseProperties.correctRho(p_rgh - p_rgh_0);
+    rho = multiphaseProperties.rho();
+
+    K = 0.5*magSqr(U);
+
+    Info<< "max(U) " << max(mag(U)).value() << endl;
+    Info<< "min(p_rgh) " << min(p_rgh).value() << endl;
+}

applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C

@@ -126,11 +126,25 @@ void Foam::multiphaseSystem::solveAlphas()
         }
 
         // Ensure that the flux at inflow BCs is preserved
-        phiAlphaCorr.boundaryField() = min
-        (
-            phase1.phi().boundaryField()*alpha1.boundaryField(),
-            phiAlphaCorr.boundaryField()
-        );
+        forAll(phiAlphaCorr.boundaryField(), patchi)
+        {
+            fvsPatchScalarField& phiAlphaCorrp =
+                phiAlphaCorr.boundaryField()[patchi];
+
+            if (!phiAlphaCorrp.coupled())
+            {
+                const scalarField& phi1p = phase1.phi().boundaryField()[patchi];
+                const scalarField& alpha1p = alpha1.boundaryField()[patchi];
+
+                forAll(phiAlphaCorrp, facei)
+                {
+                    if (phi1p[facei] < 0)
+                    {
+                        phiAlphaCorrp[facei] = alpha1p[facei]*phi1p[facei];
+                    }
+                }
+            }
+        }
 
         MULES::limit
         (

applications/solvers/multiphase/twoPhaseEulerFoam/EEqns.H

@@ -23,7 +23,7 @@
         )/rho1
 
       //***HGW- fvm::laplacian(alpha1*turbulence1->alphaEff(), he1)
-      - fvm::laplacian(alpha1*turbulence1->nuEff(), he1)
+      - fvm::laplacian(alpha1*phase1.turbulence().nuEff(), he1)
      ==
         heatTransferCoeff*(thermo2.T() - thermo1.T())/rho1
       + heatTransferCoeff*he1/Cpv1/rho1
@@ -46,7 +46,7 @@
         )/rho2
 
       //***HGW- fvm::laplacian(alpha2*turbulence2->alphaEff(), he2)
-      - fvm::laplacian(alpha2*turbulence2->nuEff(), he2)
+      - fvm::laplacian(alpha2*phase2.turbulence().nuEff(), he2)
      ==
         heatTransferCoeff*(thermo1.T() - thermo2.T())/rho2
       + heatTransferCoeff*he2/Cpv2/rho2

applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H

@@ -26,7 +26,7 @@ volScalarField dragCoeff(fluid.dragCoeff());
               - fvm::Sp(fvc::div(phi1), U1)
             )
 
-          + turbulence1->divDevReff(U1)
+          + phase1.turbulence().divDevReff(U1)
          ==
           - fvm::Sp(dragCoeff/rho1, U1)
           - alpha1*alpha2/rho1*(liftForce - fluid.Cvm()*rho2*DDtU2)
@@ -50,7 +50,7 @@ volScalarField dragCoeff(fluid.dragCoeff());
               + fvm::div(phi2, U2)
               - fvm::Sp(fvc::div(phi2), U2)
             )
-          + turbulence2->divDevReff(U2)
+          + phase2.turbulence().divDevReff(U2)
          ==
           - fvm::Sp(dragCoeff/rho2, U2)
           + alpha1*alpha2/rho2*(liftForce + fluid.Cvm()*rho2*DDtU1)

applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H

@@ -1,165 +0,0 @@
-{
-    word alphaScheme("div(phi," + alpha1.name() + ')');
-    word alpharScheme("div(phir," + alpha1.name() + ')');
-
-    alpha1.correctBoundaryConditions();
-
-    surfaceScalarField phic("phic", phi);
-    surfaceScalarField phir("phir", phi1 - phi2);
-
-    surfaceScalarField alpha1f(fvc::interpolate(max(alpha1, scalar(0))));
-
-    tmp<surfaceScalarField> pPrimeByA;
-
-    if (implicitPhasePressure)
-    {
-        pPrimeByA =
-            fvc::interpolate((1.0/rho1)*rAU1*turbulence1().pPrime())
-          + fvc::interpolate((1.0/rho2)*rAU2*turbulence2().pPrime());
-
-        surfaceScalarField phiP
-        (
-            pPrimeByA()*fvc::snGrad(alpha1, "bounded")*mesh.magSf()
-        );
-
-        phic += alpha1f*phiP;
-        phir += phiP;
-    }
-
-    for (int acorr=0; acorr<nAlphaCorr; acorr++)
-    {
-        volScalarField::DimensionedInternalField Sp
-        (
-            IOobject
-            (
-                "Sp",
-                runTime.timeName(),
-                mesh
-            ),
-            mesh,
-            dimensionedScalar("Sp", dgdt.dimensions(), 0.0)
-        );
-
-        volScalarField::DimensionedInternalField Su
-        (
-            IOobject
-            (
-                "Su",
-                runTime.timeName(),
-                mesh
-            ),
-            // Divergence term is handled explicitly to be
-            // consistent with the explicit transport solution
-            fvc::div(phi)*min(alpha1, scalar(1))
-        );
-
-        forAll(dgdt, celli)
-        {
-            if (dgdt[celli] > 0.0 && alpha1[celli] > 0.0)
-            {
-                Sp[celli] -= dgdt[celli]*alpha1[celli];
-                Su[celli] += dgdt[celli]*alpha1[celli];
-            }
-            else if (dgdt[celli] < 0.0 && alpha1[celli] < 1.0)
-            {
-                Sp[celli] += dgdt[celli]*(1.0 - alpha1[celli]);
-            }
-        }
-
-        dimensionedScalar totalDeltaT = runTime.deltaT();
-        if (nAlphaSubCycles > 1)
-        {
-            alphaPhi1 = dimensionedScalar("0", alphaPhi1.dimensions(), 0);
-        }
-
-        for
-        (
-            subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
-            !(++alphaSubCycle).end();
-        )
-        {
-            surfaceScalarField alphaPhic1
-            (
-                fvc::flux
-                (
-                    phic,
-                    alpha1,
-                    alphaScheme
-                )
-              + fvc::flux
-                (
-                    -fvc::flux(-phir, scalar(1) - alpha1, alpharScheme),
-                    alpha1,
-                    alpharScheme
-                )
-            );
-
-            // Ensure that the flux at inflow BCs is preserved
-            forAll(alphaPhic1.boundaryField(), patchi)
-            {
-                fvsPatchScalarField& alphaPhic1p =
-                    alphaPhic1.boundaryField()[patchi];
-
-                if (!alphaPhic1p.coupled())
-                {
-                    const scalarField& phi1p = phi1.boundaryField()[patchi];
-                    const scalarField& alpha1p = alpha1.boundaryField()[patchi];
-
-                    forAll(alphaPhic1p, facei)
-                    {
-                        if (phi1p[facei] < 0)
-                        {
-                            alphaPhic1p[facei] = alpha1p[facei]*phi1p[facei];
-                        }
-                    }
-                }
-            }
-
-            MULES::explicitSolve
-            (
-                geometricOneField(),
-                alpha1,
-                phi,
-                alphaPhic1,
-                Sp,
-                Su,
-                1,
-                0
-            );
-
-            if (nAlphaSubCycles > 1)
-            {
-                alphaPhi1 += (runTime.deltaT()/totalDeltaT)*alphaPhic1;
-            }
-            else
-            {
-                alphaPhi1 = alphaPhic1;
-            }
-        }
-
-        if (implicitPhasePressure)
-        {
-            fvScalarMatrix alpha1Eqn
-            (
-                fvm::ddt(alpha1) - fvc::ddt(alpha1)
-              - fvm::laplacian(alpha1f*pPrimeByA, alpha1, "bounded")
-            );
-
-            alpha1Eqn.relax();
-            alpha1Eqn.solve();
-
-            alphaPhi1 += alpha1Eqn.flux();
-        }
-
-        alphaPhi2 = phi - alphaPhi1;
-        alpha2 = scalar(1) - alpha1;
-
-        Info<< "Dispersed phase volume fraction = "
-            << alpha1.weightedAverage(mesh.V()).value()
-            << "  Min(alpha1) = " << min(alpha1).value()
-            << "  Max(alpha1) = " << max(alpha1).value()
-            << endl;
-    }
-}
-
-rho = fluid.rho();

applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H

@@ -10,19 +10,13 @@
 
     volVectorField& U1 = phase1.U();
     surfaceScalarField& phi1 = phase1.phi();
-    surfaceScalarField alphaPhi1
-    (
-        IOobject::groupName("alphaPhi", phase1.name()),
-        fvc::interpolate(alpha1)*phi1
-    );
+    surfaceScalarField& alphaPhi1 = phase1.phiAlpha();
 
     volVectorField& U2 = phase2.U();
     surfaceScalarField& phi2 = phase2.phi();
-    surfaceScalarField alphaPhi2
-    (
-        IOobject::groupName("alphaPhi", phase2.name()),
-        fvc::interpolate(alpha2)*phi2
-    );
+    surfaceScalarField& alphaPhi2 = phase2.phiAlpha();
+
+    surfaceScalarField& phi = fluid.phi();
 
     dimensionedScalar pMin
     (
@@ -55,19 +49,6 @@
         fluid.U()
     );
 
-    surfaceScalarField phi
-    (
-        IOobject
-        (
-            "phi",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        fluid.phi()
-    );
-
     volScalarField rho
     (
         IOobject
@@ -97,10 +78,6 @@
       - fvc::div(phi2)*U2
     );
 
-
-    Info<< "Calculating field g.h\n" << endl;
-    volScalarField gh("gh", g & mesh.C());
-
     volScalarField rAU1
     (
         IOobject
@@ -133,13 +110,6 @@
     scalar pRefValue = 0.0;
     setRefCell(p, mesh.solutionDict().subDict("PIMPLE"), pRefCell, pRefValue);
 
-
-    volScalarField dgdt
-    (
-        pos(alpha2)*fvc::div(phi)/max(alpha2, scalar(0.0001))
-    );
-
-
     Info<< "Creating field dpdt\n" << endl;
     volScalarField dpdt
     (
@@ -157,29 +127,3 @@
     Info<< "Creating field kinetic energy K\n" << endl;
     volScalarField K1(IOobject::groupName("K", phase1.name()), 0.5*magSqr(U1));
     volScalarField K2(IOobject::groupName("K", phase2.name()), 0.5*magSqr(U2));
-
-    autoPtr<PhaseIncompressibleTurbulenceModel<phaseModel> >
-    turbulence1
-    (
-        PhaseIncompressibleTurbulenceModel<phaseModel>::New
-        (
-            alpha1,
-            U1,
-            alphaPhi1,
-            phi1,
-            phase1
-        )
-    );
-
-    autoPtr<PhaseIncompressibleTurbulenceModel<phaseModel> >
-    turbulence2
-    (
-        PhaseIncompressibleTurbulenceModel<phaseModel>::New
-        (
-            alpha2,
-            U2,
-            alphaPhi2,
-            phi2,
-            phase2
-        )
-    );

applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H

@@ -31,17 +31,19 @@
     surfaceScalarField phiP1
     (
         "phiP1",
-        fvc::interpolate((1.0/rho1)*rAU1*turbulence1().pPrime())
+        fvc::interpolate((1.0/rho1)*rAU1*phase1.turbulence().pPrime())
        *fvc::snGrad(alpha1)*mesh.magSf()
     );
+    phiP1.boundaryField() == 0;
 
     // Phase-2 pressure flux (e.g. due to particle-particle pressure)
     surfaceScalarField phiP2
     (
         "phiP2",
-        fvc::interpolate((1.0/rho2)*rAU2*turbulence2().pPrime())
+        fvc::interpolate((1.0/rho2)*rAU2*phase2.turbulence().pPrime())
        *fvc::snGrad(alpha2)*mesh.magSf()
     );
+    phiP2.boundaryField() == 0;
 
     surfaceScalarField phiHbyA1
     (
@@ -126,9 +128,11 @@
         );
 
         pEqnComp1 =
-            fvc::ddt(rho1)
-          + fvc::div(phi1, rho1) - fvc::Sp(fvc::div(phi1), rho1)
-          + correction
+            (
+                fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1, rho1)
+              - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
+            )/rho1
+          + (alpha1/rho1)*correction
             (
                 psi1*fvm::ddt(p)
               + fvm::div(phid1, p) - fvm::Sp(fvc::div(phid1), p)
@@ -137,9 +141,11 @@
         pEqnComp1().relax();
 
         pEqnComp2 =
-            fvc::ddt(rho2)
-          + fvc::div(phi2, rho2) - fvc::Sp(fvc::div(phi2), rho2)
-          + correction
+            (
+                fvc::ddt(alpha2, rho2) + fvc::div(alphaPhi2, rho2)
+              - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
+            )/rho2
+          + (alpha2/rho2)*correction
             (
                 psi2*fvm::ddt(p)
               + fvm::div(phid2, p) - fvm::Sp(fvc::div(phid2), p)
@@ -150,12 +156,18 @@
     else
     {
         pEqnComp1 =
-            fvc::ddt(rho1) + psi1*correction(fvm::ddt(p))
-          + fvc::div(phi1, rho1) - fvc::Sp(fvc::div(phi1), rho1);
+            (
+                fvc::ddt(alpha1, rho1) + fvc::div(alphaPhi1, rho1)
+              - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
+            )/rho1
+          + (alpha1*psi1/rho1)*correction(fvm::ddt(p));
 
         pEqnComp2 =
-            fvc::ddt(rho2) + psi2*correction(fvm::ddt(p))
-          + fvc::div(phi2, rho2) - fvc::Sp(fvc::div(phi2), rho2);
+            (
+                fvc::ddt(alpha2, rho2) + fvc::div(alphaPhi2, rho2)
+              - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
+            )/rho2
+          + (alpha2*psi2/rho2)*correction(fvm::ddt(p));
     }
 
     // Cache p prior to solve for density update
@@ -171,11 +183,7 @@
 
         solve
         (
-            (
-                (alpha1/rho1)*pEqnComp1()
-              + (alpha2/rho2)*pEqnComp2()
-            )
-          + pEqnIncomp,
+            pEqnComp1() + pEqnComp2() + pEqnIncomp,
             mesh.solver(p.select(pimple.finalInnerIter()))
         );
 
@@ -199,10 +207,10 @@
 
             phi = alpha1f*phi1 + alpha2f*phi2;
 
-            dgdt =
+            fluid.dgdt() =
             (
-                pos(alpha2)*(pEqnComp2 & p)/rho2
-              - pos(alpha1)*(pEqnComp1 & p)/rho1
+                pos(alpha2)*(pEqnComp2 & p)/max(alpha2, scalar(1e-3))
+              - pos(alpha1)*(pEqnComp1 & p)/max(alpha1, scalar(1e-3))
             );
 
             p.relax();

applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H

@@ -1,7 +0,0 @@
-    #include "readTimeControls.H"
-    #include "alphaControls.H"
-
-    Switch implicitPhasePressure
-    (
-        alphaControls.lookupOrDefault<Switch>("implicitPhasePressure", false)
-    );

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C

@@ -31,13 +31,10 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "MULES.H"
-#include "subCycle.H"
-#include "rhoThermo.H"
 #include "twoPhaseSystem.H"
+#include "PhaseIncompressibleTurbulenceModel.H"
 #include "dragModel.H"
 #include "heatTransferModel.H"
-#include "PhaseIncompressibleTurbulenceModel.H"
 #include "pimpleControl.H"
 #include "IOMRFZoneList.H"
 #include "fixedFluxPressureFvPatchScalarField.H"
@@ -66,7 +63,7 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readTwoPhaseEulerFoamControls.H"
+        #include "readTimeControls.H"
         #include "CourantNos.H"
         #include "setDeltaT.H"
 
@@ -76,7 +73,10 @@ int main(int argc, char *argv[])
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {
-            #include "alphaEqn.H"
+            fluid.solve();
+            rho = fluid.rho();
+            fluid.correct();
+
             #include "EEqns.H"
             #include "UEqns.H"
 
@@ -90,8 +90,7 @@ int main(int argc, char *argv[])
 
             if (pimple.turbCorr())
             {
-                turbulence1->correct();
-                turbulence2->correct();
+                fluid.correctTurbulence();
             }
         }
 

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/Make/files

@@ -4,6 +4,12 @@ diameterModels/diameterModel/newDiameterModel.C
 diameterModels/constantDiameter/constantDiameter.C
 diameterModels/isothermalDiameter/isothermalDiameter.C
 
+diameterModels/IATE/IATE.C
+diameterModels/IATE/IATEsources/IATEsource/IATEsource.C
+diameterModels/IATE/IATEsources/wakeEntrainmentCoalescence/wakeEntrainmentCoalescence.C
+diameterModels/IATE/IATEsources/turbulentBreakUp/turbulentBreakUp.C
+diameterModels/IATE/IATEsources/randomCoalescence/randomCoalescence.C
+
 twoPhaseSystem.C
 
 LIB = $(FOAM_LIBBIN)/libcompressibleTwoPhaseSystem

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/Make/options

@@ -3,7 +3,10 @@ EXE_INC = \
     -I../interfacialModels/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/transportModels/incompressible/lnInclude
+    -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/phaseIncompressible/lnInclude
 
 LIB_LIBS = \
     -lincompressibleTransportModels \

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATE.C

@@ -0,0 +1,188 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "IATE.H"
+#include "IATEsource.H"
+#include "twoPhaseSystem.H"
+#include "fvmDdt.H"
+#include "fvmDiv.H"
+#include "fvmSup.H"
+#include "fvcDdt.H"
+#include "fvcDiv.H"
+#include "fvcAverage.H"
+#include "mathematicalConstants.H"
+#include "fundamentalConstants.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+    defineTypeNameAndDebug(IATE, 0);
+
+    addToRunTimeSelectionTable
+    (
+        diameterModel,
+        IATE,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::diameterModels::IATE::IATE
+(
+    const dictionary& diameterProperties,
+    const phaseModel& phase
+)
+:
+    diameterModel(diameterProperties, phase),
+    kappai_
+    (
+        IOobject
+        (
+            IOobject::groupName("kappai", phase.name()),
+            phase_.U().time().timeName(),
+            phase_.U().mesh(),
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        phase_.U().mesh()
+    ),
+    dMax_("dMax", dimLength, diameterProperties_.lookup("dMax")),
+    dMin_("dMin", dimLength, diameterProperties_.lookup("dMin")),
+    d_
+    (
+        IOobject
+        (
+            IOobject::groupName("d", phase.name()),
+            phase_.U().time().timeName(),
+            phase_.U().mesh(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        dsm()
+    ),
+    sources_
+    (
+        diameterProperties_.lookup("sources"),
+        IATEsource::iNew(*this)
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::diameterModels::IATE::~IATE()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField> Foam::diameterModels::IATE::dsm() const
+{
+    return max(6/max(kappai_, 6/dMax_), dMin_);
+}
+
+// Placeholder for the nucleation/condensation model
+// Foam::tmp<Foam::volScalarField> Foam::diameterModels::IATE::Rph() const
+// {
+//     const volScalarField& T = phase_.thermo().T();
+//     const volScalarField& p = phase_.thermo().p();
+//
+//     scalar A, B, C, sigma, vm, Rph;
+//
+//     volScalarField ps(1e5*pow(10, A - B/(T + C)));
+//     volScalarField Dbc
+//     (
+//         4*sigma*vm/(constant::physicoChemical::k*T*log(p/ps))
+//     );
+//
+//     return constant::mathematical::pi*sqr(Dbc)*Rph;
+// }
+
+void Foam::diameterModels::IATE::correct()
+{
+    // Initialise the accumulated source term to the dilatation effect
+    volScalarField R
+    (
+        (
+            (2.0/3.0)
+           /max
+            (
+                fvc::average(phase_ + phase_.oldTime()),
+                phase_.fluid().residualPhaseFraction()
+            )
+        )
+       *(fvc::ddt(phase_) + fvc::div(phase_.phiAlpha()))
+    );
+
+    // Accumulate the run-time selectable sources
+    forAll(sources_, j)
+    {
+        R -= sources_[j].R();
+    }
+
+    // Construct the interfacial curvature equation
+    fvScalarMatrix kappaiEqn
+    (
+        fvm::ddt(kappai_) + fvm::div(phase_.phi(), kappai_)
+      - fvm::Sp(fvc::div(phase_.phi()), kappai_)
+     ==
+      - fvm::SuSp(R, kappai_)
+    //+ Rph() // Omit the nucleation/condensation term
+    );
+
+    kappaiEqn.relax();
+    kappaiEqn.solve();
+
+    // Update the Sauter-mean diameter
+    d_ = dsm();
+}
+
+
+bool Foam::diameterModels::IATE::read(const dictionary& phaseProperties)
+{
+    diameterModel::read(phaseProperties);
+
+    diameterProperties_.lookup("dMax") >> dMax_;
+    diameterProperties_.lookup("dMin") >> dMin_;
+
+    // Re-create all the sources updating number, type and coefficients
+    PtrList<IATEsource>
+    (
+        diameterProperties_.lookup("sources"),
+        IATEsource::iNew(*this)
+    ).transfer(sources_);
+
+    return true;
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATE.H

@@ -0,0 +1,154 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 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/>.
+
+Class
+    Foam::diameterModels::IATE
+
+Description
+    IATE (Interfacial Area Transport Equation) bubble diameter model.
+
+    Solves for the interfacial curvature per unit volume of the phase rather
+    than interfacial area per unit volume to avoid stability issues relating to
+    the consistency requirements between the phase fraction and interfacial area
+    per unit volume.  In every other respect this model is as presented in the
+    paper:
+
+    \verbatim
+        "Development of Interfacial Area Transport Equation"
+        M. Ishii,
+        S. Kim,
+        J Kelly,
+        Nuclear Engineering and Technology, Vol.37 No.6 December 2005
+    \endverbatim
+
+SourceFiles
+    IATE.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef IATE_H
+#define IATE_H
+
+#include "diameterModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+
+// Forward declaration of classes
+class IATEsource;
+
+/*---------------------------------------------------------------------------*\
+                           Class IATE Declaration
+\*---------------------------------------------------------------------------*/
+
+class IATE
+:
+    public diameterModel
+{
+    // Private data
+
+        //- Interfacial curvature (alpha*interfacial area)
+        volScalarField kappai_;
+
+        //- Maximum diameter used for stabilisation in the limit kappai->0
+        dimensionedScalar dMax_;
+
+        //- Minimum diameter used for stabilisation in the limit kappai->inf
+        dimensionedScalar dMin_;
+
+        //- The Sauter-mean diameter of the phase
+        volScalarField d_;
+
+        //- IATE sources
+        PtrList<IATEsource> sources_;
+
+
+    // Private member functions
+
+        tmp<volScalarField> dsm() const;
+
+
+public:
+
+    friend class IATEsource;
+
+    //- Runtime type information
+    TypeName("IATE");
+
+
+    // Constructors
+
+        //- Construct from components
+        IATE
+        (
+            const dictionary& diameterProperties,
+            const phaseModel& phase
+        );
+
+
+    //- Destructor
+    virtual ~IATE();
+
+
+    // Member Functions
+
+        //- Return the interfacial curvature
+        const volScalarField& kappai() const
+        {
+            return kappai_;
+        }
+
+        //- Return the interfacial area
+        tmp<volScalarField> a() const
+        {
+            return phase_*kappai_;
+        }
+
+        //- Return the Sauter-mean diameter
+        virtual tmp<volScalarField> d() const
+        {
+            return d_;
+        }
+
+        //- Correct the diameter field
+        virtual void correct();
+
+        //- Read phaseProperties dictionary
+        virtual bool read(const dictionary& phaseProperties);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace diameterModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/IATEsource/IATEsource.C

@@ -0,0 +1,147 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "IATEsource.H"
+#include "twoPhaseSystem.H"
+#include "fvMatrix.H"
+#include "PhaseIncompressibleTurbulenceModel.H"
+#include "uniformDimensionedFields.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+    defineTypeNameAndDebug(IATEsource, 0);
+    defineRunTimeSelectionTable(IATEsource, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::diameterModels::IATEsource>
+Foam::diameterModels::IATEsource::New
+(
+    const word& type,
+    const IATE& iate,
+    const dictionary& dict
+)
+{
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(type);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "IATEsource::New"
+            "(const word& type, const IATE&, const dictionary&)"
+        )   << "Unknown IATE source type "
+            << type << nl << nl
+            << "Valid IATE source types : " << endl
+            << dictionaryConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<IATEsource>(cstrIter()(iate, dict));
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField> Foam::diameterModels::IATEsource::Ur() const
+{
+    const uniformDimensionedVectorField& g =
+        phase().U().db().lookupObject<uniformDimensionedVectorField>("g");
+
+    return
+        sqrt(2.0)
+       *pow025
+        (
+            fluid().sigma()*mag(g)
+           *(otherPhase().rho() - phase().rho())
+           /sqr(otherPhase().rho())
+        )
+       *pow(max(1 - phase(), scalar(0)), 1.75);
+}
+
+Foam::tmp<Foam::volScalarField> Foam::diameterModels::IATEsource::Ut() const
+{
+    return sqrt(2*otherPhase().turbulence().k());
+}
+
+Foam::tmp<Foam::volScalarField> Foam::diameterModels::IATEsource::Re() const
+{
+    return max(Ur()*phase().d()/otherPhase().nu(), scalar(1.0e-3));
+}
+
+Foam::tmp<Foam::volScalarField> Foam::diameterModels::IATEsource::CD() const
+{
+    const volScalarField Eo(this->Eo());
+    const volScalarField Re(this->Re());
+
+    return
+        max
+        (
+            min
+            (
+                (16/Re)*(1 + 0.15*pow(Re, 0.687)),
+                48/Re
+            ),
+            8*Eo/(3*(Eo + 4))
+        );
+}
+
+Foam::tmp<Foam::volScalarField> Foam::diameterModels::IATEsource::Mo() const
+{
+    const uniformDimensionedVectorField& g =
+        phase().U().db().lookupObject<uniformDimensionedVectorField>("g");
+
+    return
+        mag(g)*pow4(otherPhase().nu())*sqr(otherPhase().rho())
+       *(otherPhase().rho() - phase().rho())
+       /pow3(fluid().sigma());
+}
+
+Foam::tmp<Foam::volScalarField> Foam::diameterModels::IATEsource::Eo() const
+{
+    const uniformDimensionedVectorField& g =
+        phase().U().db().lookupObject<uniformDimensionedVectorField>("g");
+
+    return
+        mag(g)*sqr(phase().d())
+       *(otherPhase().rho() - phase().rho())
+       /fluid().sigma();
+}
+
+Foam::tmp<Foam::volScalarField> Foam::diameterModels::IATEsource::We() const
+{
+    return otherPhase().rho()*sqr(Ur())*phase().d()/fluid().sigma();
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/IATEsource/IATEsource.H

@@ -0,0 +1,192 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 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/>.
+
+Class
+    Foam::diameterModels::IATEsource
+
+Description
+    IATE (Interfacial Area Transport Equation) bubble diameter model
+    run-time selectable sources.
+
+SourceFiles
+    IATEsource.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef IATEsource_H
+#define IATEsource_H
+
+#include "IATE.H"
+#include "mathematicalConstants.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class IATEsource Declaration
+\*---------------------------------------------------------------------------*/
+
+class IATEsource
+{
+
+protected:
+
+    // Protected data
+
+        //- Reference to the IATE this source applies to
+        const IATE& iate_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("IATEsource");
+
+
+    // Declare run-time constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            IATEsource,
+            dictionary,
+            (
+                const IATE& iate,
+                const dictionary& dict
+            ),
+            (iate, dict)
+        );
+
+
+    //- Class used for the read-construction of
+    //  PtrLists of IATE sources
+    class iNew
+    {
+        const IATE& iate_;
+
+    public:
+
+        iNew(const IATE& iate)
+        :
+            iate_(iate)
+        {}
+
+        autoPtr<IATEsource> operator()(Istream& is) const
+        {
+            word type(is);
+            dictionary dict(is);
+            return IATEsource::New(type, iate_, dict);
+        }
+    };
+
+
+    // Constructors
+
+        IATEsource(const IATE& iate)
+        :
+            iate_(iate)
+        {}
+
+        autoPtr<IATEsource> clone() const
+        {
+            notImplemented("autoPtr<IATEsource> clone() const");
+            return autoPtr<IATEsource>(NULL);
+        }
+
+
+    // Selectors
+
+        static autoPtr<IATEsource> New
+        (
+            const word& type,
+            const IATE& iate,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~IATEsource()
+    {}
+
+
+    // Member Functions
+
+        const phaseModel& phase() const
+        {
+            return iate_.phase();
+        }
+
+        const twoPhaseSystem& fluid() const
+        {
+            return iate_.phase().fluid();
+        }
+
+        const phaseModel& otherPhase() const
+        {
+            return phase().otherPhase();
+        }
+
+        scalar phi() const
+        {
+            return 1.0/(36*constant::mathematical::pi);
+        }
+
+        //- Return the bubble relative velocity
+        tmp<volScalarField> Ur() const;
+
+        //- Return the bubble turbulent velocity
+        tmp<volScalarField> Ut() const;
+
+        //- Return the bubble Reynolds number
+        tmp<volScalarField> Re() const;
+
+        //- Return the bubble drag coefficient
+        tmp<volScalarField> CD() const;
+
+        //- Return the bubble Morton number
+        tmp<volScalarField> Mo() const;
+
+        //- Return the bubble Eotvos number
+        tmp<volScalarField> Eo() const;
+
+        //- Return the bubble Webber number
+        tmp<volScalarField> We() const;
+
+        virtual tmp<volScalarField> R() const = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace diameterModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/dummy/dummy.C

@@ -0,0 +1,66 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "dummy.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+namespace IATEsources
+{
+    defineTypeNameAndDebug(dummy, 0);
+    addToRunTimeSelectionTable(IATEsource, dummy, word);
+}
+}
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField>
+Foam::diameterModels::IATEsources::dummy::R() const
+{
+    return tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "R",
+                iate_.phase().U().time().timeName(),
+                iate_.phase().mesh()
+            ),
+            iate_.phase().U().mesh(),
+            dimensionedScalar("R", dimless/dimTime, 0)
+        )
+    );
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/dummy/dummy.H

@@ -0,0 +1,96 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 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/>.
+
+Class
+    Foam::diameterModels::IATEsources::dummy
+
+Description
+
+SourceFiles
+    dummy.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef dummy_H
+#define dummy_H
+
+#include "IATEsource.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+namespace IATEsources
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class dummy Declaration
+\*---------------------------------------------------------------------------*/
+
+class dummy
+:
+    public IATEsource
+{
+
+public:
+
+    //- Runtime type information
+    TypeName("dummy");
+
+    // Constructors
+
+        dummy
+        (
+            const word& name,
+            const IATE& iate,
+            const dictionary& dict
+        )
+        :
+            IATEsource(iate)
+        {}
+
+
+    //- Destructor
+    virtual ~dummy()
+    {}
+
+
+    // Member Functions
+
+        virtual tmp<volScalarField> R() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace IATEsources
+} // End namespace diameterModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/randomCoalescence/randomCoalescence.C

@@ -0,0 +1,109 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "randomCoalescence.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+namespace IATEsources
+{
+    defineTypeNameAndDebug(randomCoalescence, 0);
+    addToRunTimeSelectionTable(IATEsource, randomCoalescence, dictionary);
+}
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::diameterModels::IATEsources::randomCoalescence::
+randomCoalescence
+(
+    const IATE& iate,
+    const dictionary& dict
+)
+:
+    IATEsource(iate),
+    Crc_("Crc", dimless, dict.lookup("Crc")),
+    C_("C", dimless, dict.lookup("C")),
+    alphaMax_("alphaMax", dimless, dict.lookup("alphaMax"))
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField>
+Foam::diameterModels::IATEsources::randomCoalescence::R() const
+{
+    tmp<volScalarField> tR
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "R",
+                iate_.phase().U().time().timeName(),
+                iate_.phase().mesh()
+            ),
+            iate_.phase().U().mesh(),
+            dimensionedScalar("R", dimless/dimTime, 0)
+        )
+    );
+
+    volScalarField R = tR();
+
+    scalar Crc = Crc_.value();
+    scalar C = C_.value();
+    scalar alphaMax = alphaMax_.value();
+    volScalarField Ut(this->Ut());
+    const volScalarField& alpha = phase();
+    const volScalarField& kappai = iate_.kappai();
+    scalar cbrtAlphaMax = cbrt(alphaMax);
+
+    forAll(R, celli)
+    {
+        if (alpha[celli] < alphaMax - SMALL)
+        {
+            scalar cbrtAlphaMaxMAlpha = cbrtAlphaMax - cbrt(alpha[celli]);
+
+            R[celli] =
+                12*phi()*kappai[celli]*alpha[celli]
+               *Crc
+               *Ut[celli]
+               *(1 - exp(-C*cbrt(alpha[celli]*alphaMax)/cbrtAlphaMaxMAlpha))
+               /(cbrtAlphaMax*cbrtAlphaMaxMAlpha);
+        }
+    }
+
+    return tR;
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/randomCoalescence/randomCoalescence.H

@@ -0,0 +1,108 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 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/>.
+
+Class
+    Foam::diameterModels::IATEsources::randomCoalescence
+
+Description
+    Random coalescence IATE source as defined in paper:
+
+    \verbatim
+        "Development of Interfacial Area Transport Equation"
+        M. Ishii,
+        S. Kim,
+        J Kelly,
+        Nuclear Engineering and Technology, Vol.37 No.6 December 2005
+    \endverbatim
+
+
+SourceFiles
+    randomCoalescence.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef randomCoalescence_H
+#define randomCoalescence_H
+
+#include "IATEsource.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+namespace IATEsources
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class randomCoalescence Declaration
+\*---------------------------------------------------------------------------*/
+
+class randomCoalescence
+:
+    public IATEsource
+{
+    // Private data
+
+        dimensionedScalar Crc_;
+        dimensionedScalar C_;
+        dimensionedScalar alphaMax_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("randomCoalescence");
+
+    // Constructors
+
+        randomCoalescence
+        (
+            const IATE& iate,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~randomCoalescence()
+    {}
+
+
+    // Member Functions
+
+        virtual tmp<volScalarField> R() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace IATEsources
+} // End namespace diameterModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/turbulentBreakUp/turbulentBreakUp.C

@@ -0,0 +1,104 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "turbulentBreakUp.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+namespace IATEsources
+{
+    defineTypeNameAndDebug(turbulentBreakUp, 0);
+    addToRunTimeSelectionTable(IATEsource, turbulentBreakUp, dictionary);
+}
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::diameterModels::IATEsources::turbulentBreakUp::
+turbulentBreakUp
+(
+    const IATE& iate,
+    const dictionary& dict
+)
+:
+    IATEsource(iate),
+    Cti_("Cti", dimless, dict.lookup("Cti")),
+    WeCr_("WeCr", dimless, dict.lookup("WeCr"))
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField>
+Foam::diameterModels::IATEsources::turbulentBreakUp::R() const
+{
+    tmp<volScalarField> tR
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "R",
+                iate_.phase().U().time().timeName(),
+                iate_.phase().mesh()
+            ),
+            iate_.phase().U().mesh(),
+            dimensionedScalar("R", dimless/dimTime, 0)
+        )
+    );
+
+    volScalarField R = tR();
+
+    scalar Cti = Cti_.value();
+    scalar WeCr = WeCr_.value();
+    volScalarField Ut(this->Ut());
+    volScalarField We(this->We());
+    const volScalarField& d(iate_.d()());
+
+    forAll(R, celli)
+    {
+        if (We[celli] > WeCr)
+        {
+            R[celli] =
+                (1.0/3.0)
+               *Cti/d[celli]
+               *Ut[celli]
+               *sqrt(1 - WeCr/We[celli])
+               *exp(-WeCr/We[celli]);
+        }
+    }
+
+    return tR;
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/turbulentBreakUp/turbulentBreakUp.H

@@ -0,0 +1,106 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 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/>.
+
+Class
+    Foam::diameterModels::IATEsources::turbulentBreakUp
+
+Description
+    Turbulence-induced break-up IATE source as defined in paper:
+
+    \verbatim
+        "Development of Interfacial Area Transport Equation"
+        M. Ishii,
+        S. Kim,
+        J Kelly,
+        Nuclear Engineering and Technology, Vol.37 No.6 December 2005
+    \endverbatim
+
+SourceFiles
+    turbulentBreakUp.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef turbulentBreakUp_H
+#define turbulentBreakUp_H
+
+#include "IATEsource.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+namespace IATEsources
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class turbulentBreakUp Declaration
+\*---------------------------------------------------------------------------*/
+
+class turbulentBreakUp
+:
+    public IATEsource
+{
+    // Private data
+
+        dimensionedScalar Cti_;
+        dimensionedScalar WeCr_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("turbulentBreakUp");
+
+    // Constructors
+
+        turbulentBreakUp
+        (
+            const IATE& iate,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~turbulentBreakUp()
+    {}
+
+
+    // Member Functions
+
+        virtual tmp<volScalarField> R() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace IATEsources
+} // End namespace diameterModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/wakeEntrainmentCoalescence/wakeEntrainmentCoalescence.C

@@ -0,0 +1,72 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "wakeEntrainmentCoalescence.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+namespace IATEsources
+{
+    defineTypeNameAndDebug(wakeEntrainmentCoalescence, 0);
+    addToRunTimeSelectionTable
+    (
+        IATEsource,
+        wakeEntrainmentCoalescence,
+        dictionary
+    );
+}
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::diameterModels::IATEsources::wakeEntrainmentCoalescence::
+wakeEntrainmentCoalescence
+(
+    const IATE& iate,
+    const dictionary& dict
+)
+:
+    IATEsource(iate),
+    Cwe_("Cwe", dimless, dict.lookup("Cwe"))
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField>
+Foam::diameterModels::IATEsources::wakeEntrainmentCoalescence::R() const
+{
+    return (-12)*phi()*Cwe_*cbrt(CD())*iate_.a()*Ur();
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/wakeEntrainmentCoalescence/wakeEntrainmentCoalescence.H

@@ -0,0 +1,105 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 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/>.
+
+Class
+    Foam::diameterModels::IATEsources::wakeEntrainmentCoalescence
+
+Description
+    Bubble coalescence due to wake entrainment IATE source as defined in paper:
+
+    \verbatim
+        "Development of Interfacial Area Transport Equation"
+        M. Ishii,
+        S. Kim,
+        J Kelly,
+        Nuclear Engineering and Technology, Vol.37 No.6 December 2005
+    \endverbatim
+
+SourceFiles
+    wakeEntrainmentCoalescence.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef wakeEntrainmentCoalescence_H
+#define wakeEntrainmentCoalescence_H
+
+#include "IATEsource.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace diameterModels
+{
+namespace IATEsources
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class wakeEntrainmentCoalescence Declaration
+\*---------------------------------------------------------------------------*/
+
+class wakeEntrainmentCoalescence
+:
+    public IATEsource
+{
+    // Private data
+
+        dimensionedScalar Cwe_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("wakeEntrainmentCoalescence");
+
+    // Constructors
+
+        wakeEntrainmentCoalescence
+        (
+            const IATE& iate,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~wakeEntrainmentCoalescence()
+    {}
+
+
+    // Member Functions
+
+        virtual tmp<volScalarField> R() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace IATEsources
+} // End namespace diameterModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/constantDiameter/constantDiameter.C

@@ -48,12 +48,12 @@ namespace diameterModels
 
 Foam::diameterModels::constant::constant
 (
-    const dictionary& dict,
+    const dictionary& diameterProperties,
     const phaseModel& phase
 )
 :
-    diameterModel(dict, phase),
-    d_("d", dimLength, dict.lookup("d"))
+    diameterModel(diameterProperties, phase),
+    d_("d", dimLength, diameterProperties_.lookup("d"))
 {}
 
 
@@ -84,4 +84,14 @@ Foam::tmp<Foam::volScalarField> Foam::diameterModels::constant::d() const
 }
 
 
+bool Foam::diameterModels::constant::read(const dictionary& phaseProperties)
+{
+    diameterModel::read(phaseProperties);
+
+    diameterProperties_.lookup("d") >> d_;
+
+    return true;
+}
+
+
 // ************************************************************************* //

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/constantDiameter/constantDiameter.H

@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::constant
+    Foam::diameterModels::constant
 
 Description
     Constant dispersed-phase particle diameter model.
@@ -69,7 +69,7 @@ public:
         //- Construct from components
         constant
         (
-            const dictionary& dict,
+            const dictionary& diameterProperties,
             const phaseModel& phase
         );
 
@@ -80,7 +80,11 @@ public:
 
     // Member Functions
 
-        tmp<volScalarField> d() const;
+        //- Return the diameter as a field
+        virtual tmp<volScalarField> d() const;
+
+        //- Read diameterProperties dictionary
+        virtual bool read(const dictionary& diameterProperties);
 };
 
 

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/diameterModel/diameterModel.C

@@ -38,11 +38,11 @@ namespace Foam
 
 Foam::diameterModel::diameterModel
 (
-    const dictionary& dict,
+    const dictionary& diameterProperties,
     const phaseModel& phase
 )
 :
-    dict_(dict),
+    diameterProperties_(diameterProperties),
     phase_(phase)
 {}
 
@@ -53,4 +53,18 @@ Foam::diameterModel::~diameterModel()
 {}
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::diameterModel::correct()
+{}
+
+
+bool Foam::diameterModel::read(const dictionary& phaseProperties)
+{
+    diameterProperties_ = phaseProperties.subDict(type() + "Coeffs");
+
+    return true;
+}
+
+
 // ************************************************************************* //

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/diameterModel/diameterModel.H

@@ -36,12 +36,12 @@ SourceFiles
 #ifndef diameterModel_H
 #define diameterModel_H
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 #include "dictionary.H"
 #include "phaseModel.H"
 #include "runTimeSelectionTables.H"
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 namespace Foam
 {
 
@@ -51,11 +51,12 @@ namespace Foam
 
 class diameterModel
 {
+
 protected:
 
     // Protected data
 
-        const dictionary& dict_;
+        dictionary diameterProperties_;
         const phaseModel& phase_;
 
 
@@ -73,10 +74,10 @@ public:
             diameterModel,
             dictionary,
             (
-                const dictionary& dict,
+                const dictionary& diameterProperties,
                 const phaseModel& phase
             ),
-            (dict, phase)
+            (diameterProperties, phase)
         );
 
 
@@ -84,7 +85,7 @@ public:
 
         diameterModel
         (
-            const dictionary& dict,
+            const dictionary& diameterProperties,
             const phaseModel& phase
         );
 
@@ -97,15 +98,33 @@ public:
 
         static autoPtr<diameterModel> New
         (
-            const dictionary& dict,
+            const dictionary& diameterProperties,
             const phaseModel& phase
         );
 
 
     // Member Functions
 
+        //- Return the phase diameter properties dictionary
+        const dictionary& diameterProperties() const
+        {
+            return diameterProperties_;
+        }
+
+        //- Return the phase
+        const phaseModel& phase() const
+        {
+            return phase_;
+        }
+
         //- Return the phase mean diameter field
         virtual tmp<volScalarField> d() const = 0;
+
+        //- Correct the diameter field
+        virtual void correct();
+
+        //- Read phaseProperties dictionary
+        virtual bool read(const dictionary& phaseProperties) = 0;
 };
 
 

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/isothermalDiameter/isothermalDiameter.C

@@ -48,13 +48,13 @@ namespace diameterModels
 
 Foam::diameterModels::isothermal::isothermal
 (
-    const dictionary& dict,
+    const dictionary& diameterProperties,
     const phaseModel& phase
 )
 :
-    diameterModel(dict, phase),
-    d0_("d0", dimLength, dict.lookup("d0")),
-    p0_("p0", dimPressure, dict.lookup("p0"))
+    diameterModel(diameterProperties, phase),
+    d0_("d0", dimLength, diameterProperties_.lookup("d0")),
+    p0_("p0", dimPressure, diameterProperties_.lookup("p0"))
 {}
 
 
@@ -77,4 +77,15 @@ Foam::tmp<Foam::volScalarField> Foam::diameterModels::isothermal::d() const
 }
 
 
+bool Foam::diameterModels::isothermal::read(const dictionary& phaseProperties)
+{
+    diameterModel::read(phaseProperties);
+
+    diameterProperties_.lookup("d0") >> d0_;
+    diameterProperties_.lookup("p0") >> p0_;
+
+    return true;
+}
+
+
 // ************************************************************************* //

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/isothermalDiameter/isothermalDiameter.H

@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::isothermal
+    Foam::diameterModels::isothermal
 
 Description
     Isothermal dispersed-phase particle diameter model.
@@ -72,7 +72,7 @@ public:
         //- Construct from components
         isothermal
         (
-            const dictionary& dict,
+            const dictionary& diameterProperties,
             const phaseModel& phase
         );
 
@@ -83,7 +83,11 @@ public:
 
     // Member Functions
 
-        tmp<volScalarField> d() const;
+        //- Return the diameter field
+        virtual tmp<volScalarField> d() const;
+
+        //- Read phaseProperties dictionary
+        virtual bool read(const dictionary& phaseProperties);
 };
 
 

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.C

@@ -26,6 +26,8 @@ License
 #include "phaseModel.H"
 #include "twoPhaseSystem.H"
 #include "diameterModel.H"
+#include "fvMatrix.H"
+#include "PhaseIncompressibleTurbulenceModel.H"
 #include "dragModel.H"
 #include "heatTransferModel.H"
 #include "fixedValueFvPatchFields.H"
@@ -73,6 +75,17 @@ Foam::phaseModel::phaseModel
             IOobject::AUTO_WRITE
         ),
         fluid.mesh()
+    ),
+    phiAlpha_
+    (
+        IOobject
+        (
+            IOobject::groupName("alphaPhi", name_),
+            fluid.mesh().time().timeName(),
+            fluid.mesh()
+        ),
+        fluid.mesh(),
+        dimensionedScalar("0", dimensionSet(0, 3, -1, 0, 0), 0)
     )
 {
     thermo_->validate("phaseModel " + name_, "h", "e");
@@ -153,6 +166,16 @@ Foam::phaseModel::phaseModel
         phaseDict_,
         *this
     );
+
+    turbulence_ =
+        PhaseIncompressibleTurbulenceModel<phaseModel>::New
+        (
+            *this,
+            U_,
+            phiAlpha_,
+            phi(),
+            *this
+        );
 }
 
 
@@ -164,10 +187,39 @@ Foam::phaseModel::~phaseModel()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+const Foam::phaseModel& Foam::phaseModel::otherPhase() const
+{
+    return fluid_.otherPhase(*this);
+}
+
+
 Foam::tmp<Foam::volScalarField> Foam::phaseModel::d() const
 {
     return dPtr_().d();
 }
 
+Foam::PhaseIncompressibleTurbulenceModel<Foam::phaseModel>&
+Foam::phaseModel::turbulence()
+{
+    return turbulence_();
+}
+
+const Foam::PhaseIncompressibleTurbulenceModel<Foam::phaseModel>&
+Foam::phaseModel::turbulence() const
+{
+    return turbulence_();
+}
+
+void Foam::phaseModel::correct()
+{
+    return dPtr_->correct();
+}
+
+bool Foam::phaseModel::read(const dictionary& phaseProperties)
+{
+    phaseDict_ = phaseProperties.subDict(name_);
+    return dPtr_->read(phaseDict_);
+}
+
 
 // ************************************************************************* //

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/phaseModel/phaseModel.H

@@ -48,6 +48,9 @@ namespace Foam
 class twoPhaseSystem;
 class diameterModel;
 
+template<class Phase>
+class PhaseIncompressibleTurbulenceModel;
+
 
 /*---------------------------------------------------------------------------*\
                            Class phaseModel Declaration
@@ -74,12 +77,18 @@ class phaseModel
         //- Velocity
         volVectorField U_;
 
-        //- Fluxes
+        //- Volumetric flux of the phase
+        surfaceScalarField phiAlpha_;
+
+        //- Volumetric flux of the phase
         autoPtr<surfaceScalarField> phiPtr_;
 
         //- Diameter model
         autoPtr<diameterModel> dPtr_;
 
+        //- turbulence model
+        autoPtr<PhaseIncompressibleTurbulenceModel<phaseModel> > turbulence_;
+
 
 public:
 
@@ -99,19 +108,47 @@ public:
 
     // Member Functions
 
+        //- Return the name of this phase
+        const word& name() const
+        {
+            return name_;
+        }
+
         //- Return the twoPhaseSystem to which this phase belongs
         const twoPhaseSystem& fluid() const
         {
             return fluid_;
         }
 
-        const word& name() const
-        {
-            return name_;
-        }
+        //- Return the other phase in this two-phase system
+        const phaseModel& otherPhase() const;
 
+        //- Return the Sauter-mean diameter
         tmp<volScalarField> d() const;
 
+        //- Return the turbulence model
+        const PhaseIncompressibleTurbulenceModel<phaseModel>&
+            turbulence() const;
+
+        //- Return non-const access to the turbulence model
+        //  for correction
+        PhaseIncompressibleTurbulenceModel<phaseModel>&
+            turbulence();
+
+        //- Return the thermophysical model
+        const rhoThermo& thermo() const
+        {
+            return thermo_();
+        }
+
+        //- Return non-const access to the thermophysical model
+        //  for correction
+        rhoThermo& thermo()
+        {
+            return thermo_();
+        }
+
+        //- Return the laminar viscosity
         tmp<volScalarField> nu() const
         {
             return thermo_->nu();
@@ -123,57 +160,71 @@ public:
             return thermo_->nu(patchi);
         }
 
+        //- Return the thermal conductivity
         tmp<volScalarField> kappa() const
         {
             return thermo_->kappa();
         }
 
+        //- Return the specific heat capacity
         tmp<volScalarField> Cp() const
         {
             return thermo_->Cp();
         }
 
+        //- Return the density
         const volScalarField& rho() const
         {
             return thermo_->rho();
         }
 
-        const rhoThermo& thermo() const
-        {
-            return thermo_();
-        }
-
-        rhoThermo& thermo()
-        {
-            return thermo_();
-        }
-
+        //- Return the velocity
         const volVectorField& U() const
         {
             return U_;
         }
 
+        //- Return non-const access to the velocity
+        //  Used in the momentum equation
         volVectorField& U()
         {
             return U_;
         }
 
+        //- Return the volumetric flux
         const surfaceScalarField& phi() const
         {
             return phiPtr_();
         }
 
+        //- Return non-const access to the volumetric flux
         surfaceScalarField& phi()
         {
             return phiPtr_();
         }
 
-        //- Dummy correct
-        void correct()
-        {}
+        //- Return the volumetric flux of the phase
+        const surfaceScalarField& phiAlpha() const
+        {
+            return phiAlpha_;
+        }
 
-        //- Dummy read
-        bool read()
+        //- Return non-const access to the volumetric flux of the phase
+        surfaceScalarField& phiAlpha()
+        {
+            return phiAlpha_;
+        }
+
+        //- Correct the phase properties
+        //  other than the thermodynamics and turbulence
+        //  which have special treatment
+        void correct();
+
+        //- Read phaseProperties dictionary
+        virtual bool read(const dictionary& phaseProperties);
+
+        //- Dummy Read for transportModel
+        virtual bool read()
         {
             return true;
         }

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C

@@ -24,9 +24,19 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "twoPhaseSystem.H"
+#include "fvMatrix.H"
+#include "PhaseIncompressibleTurbulenceModel.H"
 #include "surfaceInterpolate.H"
-#include "fixedValueFvsPatchFields.H"
+#include "MULES.H"
+#include "subCycle.H"
+#include "fvcDdt.H"
+#include "fvcDiv.H"
+#include "fvcSnGrad.H"
+#include "fvcFlux.H"
 #include "fvcCurl.H"
+#include "fvmDdt.H"
+#include "fvmLaplacian.H"
+#include "fixedValueFvsPatchFields.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -63,6 +73,37 @@ Foam::twoPhaseSystem::twoPhaseSystem
         wordList(lookup("phases"))[1]
     ),
 
+    phi_
+    (
+        IOobject
+        (
+            "phi",
+            mesh.time().timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        this->calcPhi()
+    ),
+
+    dgdt_
+    (
+        IOobject
+        (
+            "dgdt",
+            mesh.time().timeName(),
+            mesh
+        ),
+        pos(phase2_)*fvc::div(phi_)/max(phase2_, scalar(0.0001))
+    ),
+
+    sigma_
+    (
+        "sigma",
+        dimensionSet(1, 0, -2, 0, 0),
+        lookup("sigma")
+    ),
+
     Cvm_
     (
         "Cvm",
@@ -170,7 +211,7 @@ Foam::tmp<Foam::volVectorField> Foam::twoPhaseSystem::U() const
 }
 
 
-Foam::tmp<Foam::surfaceScalarField> Foam::twoPhaseSystem::phi() const
+Foam::tmp<Foam::surfaceScalarField> Foam::twoPhaseSystem::calcPhi() const
 {
     return
         fvc::interpolate(phase1_)*phase1_.phi()
@@ -366,18 +407,242 @@ Foam::tmp<Foam::volScalarField> Foam::twoPhaseSystem::heatTransferCoeff() const
 }
 
 
+void Foam::twoPhaseSystem::solve()
+{
+    const Time& runTime = mesh_.time();
+
+    volScalarField& alpha1 = phase1_;
+    volScalarField& alpha2 = phase2_;
+
+    const surfaceScalarField& phi1 = phase1_.phi();
+    const surfaceScalarField& phi2 = phase2_.phi();
+
+    const dictionary& alphaControls = mesh_.solverDict
+    (
+        alpha1.name()
+    );
+
+    label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")));
+    label nAlphaCorr(readLabel(alphaControls.lookup("nAlphaCorr")));
+    Switch implicitPhasePressure
+    (
+        alphaControls.lookupOrDefault<Switch>("implicitPhasePressure", false)
+    );
+
+    word alphaScheme("div(phi," + alpha1.name() + ')');
+    word alpharScheme("div(phir," + alpha1.name() + ')');
+
+    alpha1.correctBoundaryConditions();
+
+
+    surfaceScalarField phic("phic", phi_);
+    surfaceScalarField phir("phir", phi1 - phi2);
+
+    surfaceScalarField alpha1f(fvc::interpolate(max(alpha1, scalar(0))));
+
+    tmp<surfaceScalarField> pPrimeByA;
+
+    if (implicitPhasePressure)
+    {
+        const volScalarField& rAU1 = mesh_.lookupObject<volScalarField>
+        (
+            IOobject::groupName("rAU", phase1_.name())
+        );
+        const volScalarField& rAU2 = mesh_.lookupObject<volScalarField>
+        (
+            IOobject::groupName("rAU", phase2_.name())
+        );
+
+        pPrimeByA =
+            fvc::interpolate((1.0/phase1_.rho())
+           *rAU1*phase1_.turbulence().pPrime())
+          + fvc::interpolate((1.0/phase2_.rho())
+           *rAU2*phase2_.turbulence().pPrime());
+
+        surfaceScalarField phiP
+        (
+            pPrimeByA()*fvc::snGrad(alpha1, "bounded")*mesh_.magSf()
+        );
+
+        phic += alpha1f*phiP;
+        phir += phiP;
+    }
+
+    for (int acorr=0; acorr<nAlphaCorr; acorr++)
+    {
+        volScalarField::DimensionedInternalField Sp
+        (
+            IOobject
+            (
+                "Sp",
+                runTime.timeName(),
+                mesh_
+            ),
+            mesh_,
+            dimensionedScalar("Sp", dgdt_.dimensions(), 0.0)
+        );
+
+        volScalarField::DimensionedInternalField Su
+        (
+            IOobject
+            (
+                "Su",
+                runTime.timeName(),
+                mesh_
+            ),
+            // Divergence term is handled explicitly to be
+            // consistent with the explicit transport solution
+            fvc::div(phi_)*min(alpha1, scalar(1))
+        );
+
+        forAll(dgdt_, celli)
+        {
+            if (dgdt_[celli] > 0.0 && alpha1[celli] > 0.0)
+            {
+                Sp[celli] -= dgdt_[celli]*alpha1[celli];
+                Su[celli] += dgdt_[celli]*alpha1[celli];
+            }
+            else if (dgdt_[celli] < 0.0 && alpha1[celli] < 1.0)
+            {
+                Sp[celli] += dgdt_[celli]*(1.0 - alpha1[celli]);
+            }
+        }
+
+        dimensionedScalar totalDeltaT = runTime.deltaT();
+        if (nAlphaSubCycles > 1)
+        {
+            phase1_.phiAlpha() =
+                dimensionedScalar("0", phase1_.phiAlpha().dimensions(), 0);
+        }
+
+        for
+        (
+            subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
+            !(++alphaSubCycle).end();
+        )
+        {
+            surfaceScalarField alphaPhic1
+            (
+                fvc::flux
+                (
+                    phic,
+                    alpha1,
+                    alphaScheme
+                )
+              + fvc::flux
+                (
+                    -fvc::flux(-phir, scalar(1) - alpha1, alpharScheme),
+                    alpha1,
+                    alpharScheme
+                )
+            );
+
+            // Ensure that the flux at inflow BCs is preserved
+            forAll(alphaPhic1.boundaryField(), patchi)
+            {
+                fvsPatchScalarField& alphaPhic1p =
+                    alphaPhic1.boundaryField()[patchi];
+
+                if (!alphaPhic1p.coupled())
+                {
+                    const scalarField& phi1p = phi1.boundaryField()[patchi];
+                    const scalarField& alpha1p = alpha1.boundaryField()[patchi];
+
+                    forAll(alphaPhic1p, facei)
+                    {
+                        if (phi1p[facei] < 0)
+                        {
+                            alphaPhic1p[facei] = alpha1p[facei]*phi1p[facei];
+                        }
+                    }
+                }
+            }
+
+            MULES::explicitSolve
+            (
+                geometricOneField(),
+                alpha1,
+                phi_,
+                alphaPhic1,
+                Sp,
+                Su,
+                1,
+                0
+            );
+
+            if (nAlphaSubCycles > 1)
+            {
+                phase1_.phiAlpha() += (runTime.deltaT()/totalDeltaT)*alphaPhic1;
+            }
+            else
+            {
+                phase1_.phiAlpha() = alphaPhic1;
+            }
+        }
+
+        if (implicitPhasePressure)
+        {
+            fvScalarMatrix alpha1Eqn
+            (
+                fvm::ddt(alpha1) - fvc::ddt(alpha1)
+              - fvm::laplacian(alpha1f*pPrimeByA, alpha1, "bounded")
+            );
+
+            alpha1Eqn.relax();
+            alpha1Eqn.solve();
+
+            phase1_.phiAlpha() += alpha1Eqn.flux();
+        }
+
+        phase2_.phiAlpha() = phi_ - phase1_.phiAlpha();
+        alpha2 = scalar(1) - alpha1;
+
+        Info<< alpha1.name() << " volume fraction = "
+            << alpha1.weightedAverage(mesh_.V()).value()
+            << "  Min(alpha1) = " << min(alpha1).value()
+            << "  Max(alpha1) = " << max(alpha1).value()
+            << endl;
+    }
+}
+
+
+void Foam::twoPhaseSystem::correct()
+{
+    phase1_.correct();
+    phase2_.correct();
+}
+
+
+void Foam::twoPhaseSystem::correctTurbulence()
+{
+    phase1_.turbulence().correct();
+    phase2_.turbulence().correct();
+}
+
+
 bool Foam::twoPhaseSystem::read()
 {
     if (regIOobject::read())
     {
         bool readOK = true;
 
-        readOK &= phase1_.read();
-        readOK &= phase2_.read();
+        readOK &= phase1_.read(*this);
+        readOK &= phase2_.read(*this);
 
+        lookup("sigma") >> sigma_;
         lookup("Cvm") >> Cvm_;
         lookup("Cl") >> Cl_;
 
+        // drag1_->read(*this);
+        // drag2_->read(*this);
+
+        // heatTransfer1_->read(*this);
+        // heatTransfer2_->read(*this);
+
+        lookup("dispersedPhase") >> dispersedPhase_;
+        lookup("residualPhaseFraction") >> residualPhaseFraction_;
+        lookup("residualSlip") >> residualSlip_;
+
         return readOK;
     }
     else

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H

@@ -74,6 +74,12 @@ class twoPhaseSystem
         phaseModel phase1_;
         phaseModel phase2_;
 
+        surfaceScalarField phi_;
+
+        volScalarField dgdt_;
+
+        dimensionedScalar sigma_;
+
         dimensionedScalar Cvm_;
         dimensionedScalar Cl_;
 
@@ -88,6 +94,10 @@ class twoPhaseSystem
         dimensionedScalar residualSlip_;
 
 
+        //- Return the mixture flux
+        tmp<surfaceScalarField> calcPhi() const;
+
+
 public:
 
     // Constructors
@@ -140,6 +150,28 @@ public:
             return phase2_;
         }
 
+        //- Return the mixture flux
+        const surfaceScalarField& phi() const
+        {
+            return phi_;
+        }
+
+        //- Return the mixture flux
+        surfaceScalarField& phi()
+        {
+            return phi_;
+        }
+
+        const volScalarField& dgdt() const
+        {
+            return dgdt_;
+        }
+
+        volScalarField& dgdt()
+        {
+            return dgdt_;
+        }
+
         const dragModel& drag1() const
         {
             return drag1_();
@@ -193,8 +225,11 @@ public:
         //- Return the mixture velocity
         tmp<volVectorField> U() const;
 
-        //- Return the mixture flux
-        tmp<surfaceScalarField> phi() const;
+        //- Return the surface tension coefficient
+        dimensionedScalar sigma() const
+        {
+            return sigma_;
+        }
 
         //- Return the virtual-mass coefficient
         dimensionedScalar Cvm() const
@@ -208,9 +243,14 @@ public:
             return Cl_;
         }
 
-        //- Dummy correct
-        void correct()
-        {}
+        //- Solve for the two-phase-fractions
+        void solve();
+
+        //- Correct two-phase properties other than turbulence
+        void correct();
+
+        //- Correct two-phase turbulence
+        void correctTurbulence();
 
         //- Read base phaseProperties dictionary
         bool read();

applications/test/HashingSpeed/Test-HashingSpeed.C

@@ -1016,7 +1016,6 @@ int i;
 // found somewhere in the 2nd addition
 uint32_t stroustrup (const char * s, int len) {
     uint32_t h;
-    int i;
 
     for (register int i=0; i < len; ++s, ++i)
     {

applications/test/List/Test-List.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -137,7 +137,7 @@ int main(int argc, char *argv[])
     Info<< "    std::list constructed from Foam::List: ";
 
     std::list<vector>::iterator it;
-    for (it=stlList.begin(); it != stlList.end(); it++)
+    for (it=stlList.begin(); it != stlList.end(); ++it)
     {
         Info<< *it << " ";
     }

applications/test/string/Test-string.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -79,7 +79,7 @@ int main(int argc, char *argv[])
     while (iter != iter2)
     {
         Info<< *iter;
-        iter++;
+        ++iter;
     }
     Info<< "<\n";
 

applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict

@@ -477,7 +477,9 @@ meshQualityControls
 //  0 : only write final meshes
 //  1 : write intermediate meshes
 //  2 : write volScalarField with cellLevel for postprocessing
-// 4 : write current intersections as .obj files
+//  4 : write current mesh intersections as .obj files
+//  8 : write information about explicit feature edge refinement
+// 16 : write information about layers
 debug 0;
 
 // Merge tolerance. Is fraction of overall bounding box of initial mesh.

applications/utilities/mesh/manipulation/checkMesh/checkTopology.C

@@ -1,3 +1,28 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 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/>.
+
+\*---------------------------------------------------------------------------*/
+
 #include "checkTopology.H"
 #include "polyMesh.H"
 #include "Time.H"
@@ -9,6 +34,8 @@
 #include "emptyPolyPatch.H"
 #include "processorPolyPatch.H"
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 Foam::label Foam::checkTopology
 (
     const polyMesh& mesh,
@@ -287,6 +314,37 @@ Foam::label Foam::checkTopology
                 rs
             );
             ctr.write();
+
+
+            // write cellSet for each region
+            PtrList<cellSet> cellRegions(rs.nRegions());
+            for (label i = 0; i < rs.nRegions(); i++)
+            {
+                cellRegions.set
+                (
+                    i,
+                    new cellSet
+                    (
+                        mesh,
+                        "region" + Foam::name(i),
+                        mesh.nCells()/100
+                    )
+                );
+            }
+
+            forAll(rs, i)
+            {
+                cellRegions[rs[i]].insert(i);
+            }
+
+            for (label i = 0; i < rs.nRegions(); i++)
+            {
+                Info<< "  <<Writing region " << i << " with "
+                    << returnReduce(cellRegions[i].size(), sumOp<scalar>())
+                    << " cells to cellSet " << cellRegions[i].name() << endl;
+
+                cellRegions[i].write();
+            }
         }
     }
 
@@ -437,3 +495,6 @@ Foam::label Foam::checkTopology
 
     return noFailedChecks;
 }
+
+
+// ************************************************************************* //

applications/utilities/mesh/manipulation/orientFaceZone/cavity/Allrun

@@ -1,22 +0,0 @@
-#!/bin/sh
-cd ${0%/*} || exit 1    # run from this directory
-
-# Source tutorial run functions
-. $WM_PROJECT_DIR/bin/tools/RunFunctions
-
-# Get application name
-application=`getApplication`
-
-runApplication blockMesh
-# Make random cell numbering to get faceZone orientation random
-runApplication renumberMesh -dict system/renumberMeshDict
-
-# Construct faceZone
-runApplication topoSet
-
-runApplication decomposePar -force -cellDist
-runParallel orientFaceZone 2 f0 '(10 0 0)'
-runApplication reconstructParMesh -latestTime -mergeTol 1e-6
-
-
-# ----------------------------------------------------------------- end-of-file

applications/utilities/mesh/manipulation/orientFaceZone/cavity/README.txt

@@ -1 +0,0 @@
-2013-09-16 To test orientFaceZone. Make random orientation. Orient afterwards.

applications/utilities/mesh/manipulation/orientFaceZone/cavity/system/decomposeParDict

@@ -1,135 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    note        "mesh decomposition control dictionary";
-    object      decomposeParDict;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-numberOfSubdomains  2;
-
-//- Keep owner and neighbour on same processor for faces in zones:
-// preserveFaceZones (heater solid1 solid3);
-
-//- Keep owner and neighbour on same processor for faces in patches:
-//  (makes sense only for cyclic patches)
-//preservePatches (cyclic_half0 cyclic_half1);
-
-//- Keep all of faceSet on a single processor. This puts all cells
-//  connected with a point, edge or face on the same processor.
-//  (just having face connected cells might not guarantee a balanced
-//  decomposition)
-// The processor can be -1 (the decompositionMethod chooses the processor
-// for a good load balance) or explicitly provided (upsets balance).
-//singleProcessorFaceSets ((f0 -1));
-
-
-//- Use the volScalarField named here as a weight for each cell in the
-//  decomposition.  For example, use a particle population field to decompose
-//  for a balanced number of particles in a lagrangian simulation.
-// weightField dsmcRhoNMean;
-
-//method          scotch;
-method          hierarchical;
-// method          simple;
-// method          metis;
-// method          manual;
-// method          multiLevel;
-// method          structured;  // does 2D decomposition of structured mesh
-
-multiLevelCoeffs
-{
-    // Decomposition methods to apply in turn. This is like hierarchical but
-    // fully general - every method can be used at every level.
-
-    level0
-    {
-        numberOfSubdomains  64;
-        //method simple;
-        //simpleCoeffs
-        //{
-        //    n           (2 1 1);
-        //    delta       0.001;
-        //}
-        method scotch;
-    }
-    level1
-    {
-        numberOfSubdomains  4;
-        method scotch;
-    }
-}
-
-// Desired output
-
-simpleCoeffs
-{
-    n           (2 1 1);
-    delta       0.001;
-}
-
-hierarchicalCoeffs
-{
-    n           (2 1 1);
-    delta       0.001;
-    order       xyz;
-}
-
-metisCoeffs
-{
- /*
-    processorWeights
-    (
-        1
-        1
-        1
-        1
-    );
-  */
-}
-
-scotchCoeffs
-{
-    //processorWeights
-    //(
-    //    1
-    //    1
-    //    1
-    //    1
-    //);
-    //writeGraph  true;
-    //strategy "b";
-}
-
-manualCoeffs
-{
-    dataFile    "decompositionData";
-}
-
-structuredCoeffs
-{
-    // Patches to do 2D decomposition on. Structured mesh only; cells have
-    // to be in 'columns' on top of patches.
-    patches     (bottomPatch);
-}
-
-//// Is the case distributed? Note: command-line argument -roots takes
-//// precedence
-//distributed     yes;
-//// Per slave (so nProcs-1 entries) the directory above the case.
-//roots
-//(
-//    "/tmp"
-//    "/tmp"
-//);
-
-// ************************************************************************* //

applications/utilities/mesh/manipulation/orientFaceZone/cavity/system/renumberMeshDict

@@ -1,110 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    note        "mesh renumbering dictionary";
-    object      renumberMeshDict;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-// Write maps from renumbered back to original mesh
-writeMaps true;
-
-// Optional entry: sort cells on coupled boundaries to last for use with
-// e.g. nonBlockingGaussSeidel.
-sortCoupledFaceCells false;
-
-// Optional entry: renumber on a block-by-block basis. It uses a
-// blockCoeffs dictionary to construct a decompositionMethod to do
-// a block subdivision) and then applies the renumberMethod to each
-// block in turn. This can be used in large cases to keep the blocks
-// fitting in cache with all the the cache misses bunched at the end.
-// This number is the approximate size of the blocks - this gets converted
-// to a number of blocks that is the input to the decomposition method.
-//blockSize 1000;
-
-// Optional entry: sort points into internal and boundary points
-//orderPoints false;
-
-
-
-//method          CuthillMcKee;
-//method          Sloan;
-//method          manual;
-method          random;
-//method          structured;
-//method          spring;
-//method          zoltan;             // only if compiled with zoltan support
-
-//CuthillMcKeeCoeffs
-//{
-//    // Reverse CuthillMcKee (RCM) or plain
-//    reverse true;
-//}
-
-manualCoeffs
-{
-    // In system directory: new-to-original (i.e. order) labelIOList
-    dataFile "cellMap";
-}
-
-
-// For extruded (i.e. structured in one direction) meshes
-structuredCoeffs
-{
-    // Patches that mesh was extruded from. These determine the starting
-    // layer of cells
-    patches (movingWall);
-    // Method to renumber the starting layer of cells
-    method  random;
-
-    // Renumber in columns (depthFirst) or in layers
-    depthFirst true;
-
-    // Optional: reverse ordering
-    //reverse false;
-}
-
-
-springCoeffs
-{
-    // Maximum jump of cell indices. Is fraction of number of cells
-    maxCo 0.01;
-
-    // Limit the amount of movement; the fraction maxCo gets decreased
-    // with every iteration
-    freezeFraction 0.999;
-
-    // Maximum number of iterations
-    maxIter 1000;
-}
-
-
-blockCoeffs
-{
-    method          scotch;
-    //method hierarchical;
-    //hierarchicalCoeffs
-    //{
-    //    n           (1 2 1);
-    //    delta       0.001;
-    //    order       xyz;
-    //}
-}
-
-
-zoltanCoeffs
-{
-    ORDER_METHOD    LOCAL_HSFC;
-}
-
-
-// ************************************************************************* //

applications/utilities/mesh/manipulation/orientFaceZone/cavity/system/topoSetDict

@@ -1,64 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  2.2.1                                 |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    object      topoSetDict;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-actions
-(
-    // Get 'T' of faces
-    {
-        name    f0;
-        type    faceSet;
-        action  new;
-        source  boxToFace;
-        sourceInfo
-        {
-            boxes
-            (
-                (0.0499 -100 -100)(0.0501 100 100)
-                (0.03 0.0499 -100)(0.05 0.0501 100)
-            );
-        }
-    }
-
-//    // Pick up some cells
-//    {
-//        name    c0;
-//        type    cellSet;
-//        action  new;
-//        source  boxToCell;
-//        sourceInfo
-//        {
-//            boxes
-//            (
-//                (-100 -100 -100)(0.05 0.05 100)
-//                (0.05 0.05 -100)(100 100 100)
-//            );
-//        }
-//    }
-
-    // Convert faceSet to faceZone
-    {
-        name    f0;
-        type    faceZoneSet;
-        action  new;
-        source  setToFaceZone;
-        sourceInfo
-        {
-            faceSet f0;
-        }
-    }
-);
-
-// ************************************************************************* //

applications/utilities/mesh/manipulation/topoSet/topoSetDict

@@ -368,6 +368,8 @@ FoamFile
 //    {
 //        faceSet f0;       // name of faceSet
 //        cellSet c0;       // name of cellSet of slave side
+//        flip    false;    // optional: flip the faceZone (so now the cellSet
+//                          //           is the master side)
 //    }
 //
 //    // Select based on surface. Orientation from normals on surface

applications/utilities/miscellaneous/foamInfoExec/foamInfoExec.C

@@ -46,8 +46,8 @@ int main(int argc, char *argv[])
     );
 
     argList::noBanner();
-    argList::noParallel();
     argList::addBoolOption("times", "list available times");
+    argList::addBoolOption("latestTime", "list last time");
     argList::addBoolOption
     (
         "keywords",
@@ -75,6 +75,15 @@ int main(int argc, char *argv[])
             Info<< times[i].name() << endl;
         }
     }
+    else if (args.optionFound("latestTime"))
+    {
+        instantList times
+        (
+            Foam::Time::findTimes(args.rootPath()/args.caseName())
+        );
+
+        Info<< times.last().name() << endl;
+    }
 
     if (args.optionFound("dict"))
     {

applications/utilities/parallelProcessing/decomposePar/decomposePar.C

@@ -99,6 +99,40 @@ Usage
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+const labelIOList& procAddressing
+(
+    const PtrList<fvMesh>& procMeshList,
+    const label procI,
+    const word& name,
+    PtrList<labelIOList>& procAddressingList
+)
+{
+    const fvMesh& procMesh = procMeshList[procI];
+
+    if (!procAddressingList.set(procI))
+    {
+        procAddressingList.set
+        (
+            procI,
+            new labelIOList
+            (
+                IOobject
+                (
+                    name,
+                    procMesh.facesInstance(),
+                    procMesh.meshSubDir,
+                    procMesh,
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE
+                )
+            )
+        );
+    }
+    return procAddressingList[procI];
+}
+
+
+
 int main(int argc, char *argv[])
 {
     argList::addNote
@@ -805,74 +839,29 @@ int main(int argc, char *argv[])
                 }
                 const fvMesh& procMesh = procMeshList[procI];
 
-
-                if (!faceProcAddressingList.set(procI))
-                {
-                    faceProcAddressingList.set
+                const labelIOList& faceProcAddressing = procAddressing
                 (
+                    procMeshList,
                     procI,
-                        new labelIOList
-                        (
-                            IOobject
-                            (
                     "faceProcAddressing",
-                                procMesh.facesInstance(),
-                                procMesh.meshSubDir,
-                                procMesh,
-                                IOobject::MUST_READ,
-                                IOobject::NO_WRITE
-                            )
-                        )
+                    faceProcAddressingList
                 );
-                }
-                const labelIOList& faceProcAddressing =
-                    faceProcAddressingList[procI];
 
-
-                if (!cellProcAddressingList.set(procI))
-                {
-                    cellProcAddressingList.set
+                const labelIOList& cellProcAddressing = procAddressing
                 (
+                    procMeshList,
                     procI,
-                        new labelIOList
-                        (
-                            IOobject
-                            (
                     "cellProcAddressing",
-                                procMesh.facesInstance(),
-                                procMesh.meshSubDir,
-                                procMesh,
-                                IOobject::MUST_READ,
-                                IOobject::NO_WRITE
-                            )
-                        )
+                    cellProcAddressingList
                 );
-                }
-                const labelIOList& cellProcAddressing =
-                    cellProcAddressingList[procI];
 
-
-                if (!boundaryProcAddressingList.set(procI))
-                {
-                    boundaryProcAddressingList.set
+                const labelIOList& boundaryProcAddressing = procAddressing
                 (
+                    procMeshList,
                     procI,
-                        new labelIOList
-                        (
-                            IOobject
-                            (
                     "boundaryProcAddressing",
-                                procMesh.facesInstance(),
-                                procMesh.meshSubDir,
-                                procMesh,
-                                IOobject::MUST_READ,
-                                IOobject::NO_WRITE
-                            )
-                        )
+                    boundaryProcAddressingList
                 );
-                }
-                const labelIOList& boundaryProcAddressing =
-                    boundaryProcAddressingList[procI];
 
 
                 // FV fields
@@ -959,27 +948,13 @@ int main(int argc, char *argv[])
                  || pointTensorFields.size()
                 )
                 {
-                    if (!pointProcAddressingList.set(procI))
-                    {
-                        pointProcAddressingList.set
+                    const labelIOList& pointProcAddressing = procAddressing
                     (
+                        procMeshList,
                         procI,
-                            new labelIOList
-                            (
-                                IOobject
-                                (
                         "pointProcAddressing",
-                                    procMesh.facesInstance(),
-                                    procMesh.meshSubDir,
-                                    procMesh,
-                                    IOobject::MUST_READ,
-                                    IOobject::NO_WRITE
-                                )
-                            )
+                        pointProcAddressingList
                     );
-                    }
-                    const labelIOList& pointProcAddressing =
-                        pointProcAddressingList[procI];
 
                     const pointMesh& procPMesh = pointMesh::New(procMesh);
 

applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C

@@ -38,6 +38,7 @@ License
 #include "cellSet.H"
 #include "faceSet.H"
 #include "pointSet.H"
+#include "uniformDimensionedFields.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -195,6 +196,60 @@ bool Foam::domainDecomposition::writeDecomposition(const bool decomposeSets)
     }
 
 
+    autoPtr<labelIOList> cellLevelPtr;
+    {
+        IOobject io
+        (
+            "cellLevel",
+            facesInstance(),
+            polyMesh::meshSubDir,
+            *this,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        );
+        if (io.headerOk())
+        {
+            Info<< "Reading hexRef8 data : " << io.name() << endl;
+            cellLevelPtr.reset(new labelIOList(io));
+        }
+    }
+    autoPtr<labelIOList> pointLevelPtr;
+    {
+        IOobject io
+        (
+            "pointLevel",
+            facesInstance(),
+            polyMesh::meshSubDir,
+            *this,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        );
+        if (io.headerOk())
+        {
+            Info<< "Reading hexRef8 data : " << io.name() << endl;
+            pointLevelPtr.reset(new labelIOList(io));
+        }
+    }
+    autoPtr<uniformDimensionedScalarField> level0EdgePtr;
+    {
+        IOobject io
+        (
+            "level0Edge",
+            facesInstance(),
+            polyMesh::meshSubDir,
+            *this,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        );
+        if (io.headerOk())
+        {
+            Info<< "Reading hexRef8 data : " << io.name() << endl;
+            level0EdgePtr.reset(new uniformDimensionedScalarField(io));
+        }
+    }
+
+
+
     label maxProcCells = 0;
     label totProcFaces = 0;
     label maxProcPatches = 0;
@@ -767,8 +822,28 @@ bool Foam::domainDecomposition::writeDecomposition(const bool decomposeSets)
 
         procMesh.write();
 
+        // Write points if pointsInstance differing from facesInstance
+        if (facesInstancePointsPtr_.valid())
+        {
+            pointIOField pointsInstancePoints
+            (
+                IOobject
+                (
+                    "points",
+                    pointsInstance(),
+                    polyMesh::meshSubDir,
+                    procMesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                xferMove(procPoints)
+            );
+            pointsInstancePoints.write();
+        }
 
 
+        // Decompose any sets
         if (decomposeSets)
         {
             forAll(cellSets, i)
@@ -813,25 +888,67 @@ bool Foam::domainDecomposition::writeDecomposition(const bool decomposeSets)
         }
 
 
-        // Write points if pointsInstance differing from facesInstance
-        if (facesInstancePointsPtr_.valid())
+        // hexRef8 data
+        if (cellLevelPtr.valid())
         {
-            pointIOField pointsInstancePoints
+            labelIOList
             (
                 IOobject
                 (
-                    "points",
-                    pointsInstance(),
+                    cellLevelPtr().name(),
+                    facesInstance(),
                     polyMesh::meshSubDir,
                     procMesh,
                     IOobject::NO_READ,
-                    IOobject::NO_WRITE,
-                    false
+                    IOobject::AUTO_WRITE
                 ),
-                xferMove(procPoints)
-            );
-            pointsInstancePoints.write();
+                UIndirectList<label>
+                (
+                    cellLevelPtr(),
+                    procCellAddressing_[procI]
+                )()
+            ).write();
         }
+        if (pointLevelPtr.valid())
+        {
+            labelIOList
+            (
+                IOobject
+                (
+                    pointLevelPtr().name(),
+                    facesInstance(),
+                    polyMesh::meshSubDir,
+                    procMesh,
+                    IOobject::NO_READ,
+                    IOobject::AUTO_WRITE
+                ),
+                UIndirectList<label>
+                (
+                    pointLevelPtr(),
+                    procPointAddressing_[procI]
+                )()
+            ).write();
+        }
+        if (level0EdgePtr.valid())
+        {
+            uniformDimensionedScalarField
+            (
+                IOobject
+                (
+                    level0EdgePtr().name(),
+                    facesInstance(),
+                    polyMesh::meshSubDir,
+                    procMesh,
+                    IOobject::NO_READ,
+                    IOobject::AUTO_WRITE
+                ),
+                level0EdgePtr()
+            ).write();
+        }
+
+
+
+        // Statistics
 
         Info<< endl
             << "Processor " << procI << nl

applications/utilities/postProcessing/dataConversion/foamToEnsight/cellSets.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -78,6 +78,26 @@ public:
             hexes(nCells),
             polys(nCells)
         {}
+
+
+    // Member Functions
+
+        void setSize(const label nCells)
+        {
+            nTets = 0;
+            nPyrs = 0;
+            nPrisms = 0;
+            nHexesWedges = 0;
+            nPolys = 0;
+
+            tets.setSize(nCells);
+            pyrs.setSize(nCells);
+            prisms.setSize(nCells);
+            wedges.setSize(nCells);
+            hexes.setSize(nCells);
+            polys.setSize(nCells);
+        }
+
 };
 
 

applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.C

@@ -47,7 +47,8 @@ License
 void Foam::ensightMesh::correct()
 {
     patchPartOffset_ = 2;
-    meshCellSets_ = mesh_.nCells();
+    meshCellSets_.setSize(mesh_.nCells());
+
     boundaryFaceSets_.setSize(mesh_.boundary().size());
     allPatchNames_.clear();
     patchNames_.clear();
@@ -193,8 +194,6 @@ void Foam::ensightMesh::correct()
     if (!noPatches_)
     {
         forAll(mesh_.boundary(), patchi)
-        {
-            if (mesh_.boundary()[patchi].size())
         {
             const polyPatch& p = mesh_.boundaryMesh()[patchi];
 
@@ -233,7 +232,6 @@ void Foam::ensightMesh::correct()
             polys.setSize(nPolys);
         }
     }
-    }
 
     forAll(allPatchNames_, patchi)
     {
@@ -241,14 +239,11 @@ void Foam::ensightMesh::correct()
         nFacePrimitives nfp;
 
         if (patchNames_.empty() || patchNames_.found(patchName))
-        {
-            if (mesh_.boundary()[patchi].size())
         {
             nfp.nTris   = boundaryFaceSets_[patchi].tris.size();
             nfp.nQuads  = boundaryFaceSets_[patchi].quads.size();
             nfp.nPolys  = boundaryFaceSets_[patchi].polys.size();
         }
-        }
 
         reduce(nfp.nTris, sumOp<label>());
         reduce(nfp.nQuads, sumOp<label>());
@@ -1148,6 +1143,7 @@ void Foam::ensightMesh::write
             if (nfp.nTris || nfp.nQuads || nfp.nPolys)
             {
                 const polyPatch& p = mesh_.boundaryMesh()[patchi];
+
                 const labelList& tris = boundaryFaceSets_[patchi].tris;
                 const labelList& quads = boundaryFaceSets_[patchi].quads;
                 const labelList& polys = boundaryFaceSets_[patchi].polys;

applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/PV3FoamReader/pqPV3FoamReaderPanel.cxx

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -441,8 +441,6 @@ void pqPV3FoamReaderPanel::IncludeZonesToggled()
 
 void pqPV3FoamReaderPanel::ExtrapolatePatchesToggled()
 {
-    vtkSMProperty* prop;
-
     vtkSMIntVectorProperty::SafeDownCast
     (
         this->proxy()->GetProperty("UiExtrapolatePatches")
@@ -454,8 +452,6 @@ void pqPV3FoamReaderPanel::ExtrapolatePatchesToggled()
 
 void pqPV3FoamReaderPanel::InterpolateVolFieldsToggled()
 {
-    vtkSMProperty* prop;
-
     vtkSMIntVectorProperty::SafeDownCast
     (
         this->proxy()->GetProperty("UiInterpolateVolFields")

src/OSspecific/POSIX/fileMonitor.C

@@ -133,7 +133,8 @@ namespace Foam
         //- initialise inotify
         inline fileMonitorWatcher(const bool useInotify, const label sz = 20)
         :
-            useInotify_(useInotify)
+            useInotify_(useInotify),
+            inotifyFd_(-1)
         {
             if (useInotify_)
             {

src/OpenFOAM/db/functionObjects/OutputFilterFunctionObject/OutputFilterFunctionObject.C

@@ -179,7 +179,7 @@ bool Foam::OutputFilterFunctionObject<OutputFilter>::execute
 template<class OutputFilter>
 bool Foam::OutputFilterFunctionObject<OutputFilter>::end()
 {
-    if (active())
+    if (enabled_)
     {
         if (!storeFilter_)
         {

src/OpenFOAM/fields/Fields/Field/FieldFunctions.C

@@ -358,6 +358,59 @@ Type sum(const UList<Type>& f)
 
 TMP_UNARY_FUNCTION(Type, sum)
 
+template<class Type>
+Type maxMagSqr(const UList<Type>& f)
+{
+    if (f.size())
+    {
+        Type Max(f[0]);
+        TFOR_ALL_S_OP_FUNC_F_S
+        (
+            Type,
+            Max,
+            =,
+            maxMagSqrOp<Type>(),
+            Type,
+            f,
+            Type,
+            Max
+        )
+        return Max;
+    }
+    else
+    {
+        return pTraits<Type>::min;
+    }
+}
+
+TMP_UNARY_FUNCTION(Type, maxMagSqr)
+
+template<class Type>
+Type minMagSqr(const UList<Type>& f)
+{
+    if (f.size())
+    {
+        Type Min(f[0]);
+        TFOR_ALL_S_OP_FUNC_F_S
+        (
+            Type,
+            Min,
+            =,
+            minMagSqrOp<Type>(),
+            Type,
+            f,
+            Type,
+            Min
+        )
+        return Min;
+    }
+    else
+    {
+        return pTraits<Type>::max;
+    }
+}
+
+TMP_UNARY_FUNCTION(Type, minMagSqr)
 
 template<class Type>
 scalar sumProd(const UList<Type>& f1, const UList<Type>& f2)
@@ -488,6 +541,8 @@ TMP_UNARY_FUNCTION(ReturnType, gFunc)
 G_UNARY_FUNCTION(Type, gMax, max, max)
 G_UNARY_FUNCTION(Type, gMin, min, min)
 G_UNARY_FUNCTION(Type, gSum, sum, sum)
+G_UNARY_FUNCTION(Type, gMaxMagSqr, maxMagSqr, maxMagSqr)
+G_UNARY_FUNCTION(Type, gMinMagSqr, minMagSqr, minMagSqr)
 G_UNARY_FUNCTION(scalar, gSumSqr, sumSqr, sum)
 G_UNARY_FUNCTION(scalar, gSumMag, sumMag, sum)
 G_UNARY_FUNCTION(Type, gSumCmptMag, sumCmptMag, sum)

src/OpenFOAM/fields/Fields/Field/FieldFunctions.H

@@ -171,6 +171,16 @@ Type sum(const UList<Type>& f);
 
 TMP_UNARY_FUNCTION(Type, sum)
 
+template<class Type>
+Type maxMagSqr(const UList<Type>& f);
+
+TMP_UNARY_FUNCTION(Type, maxMagSqr)
+
+template<class Type>
+Type minMagSqr(const UList<Type>& f);
+
+TMP_UNARY_FUNCTION(Type, minMagSqr)
+
 
 template<class Type>
 scalar sumProd(const UList<Type>& f1, const UList<Type>& f2);
@@ -208,6 +218,8 @@ TMP_UNARY_FUNCTION(ReturnType, gFunc)
 G_UNARY_FUNCTION(Type, gMax, max, max)
 G_UNARY_FUNCTION(Type, gMin, min, min)
 G_UNARY_FUNCTION(Type, gSum, sum, sum)
+G_UNARY_FUNCTION(Type, gMaxMagSqr, maxMagSqr, maxMagSqr)
+G_UNARY_FUNCTION(Type, gMinMagSqr, minMagSqr, minMagSqr)
 G_UNARY_FUNCTION(scalar, gSumSqr, sumSqr, sum)
 G_UNARY_FUNCTION(scalar, gSumMag, sumMag, sum)
 G_UNARY_FUNCTION(Type, gSumCmptMag, sumCmptMag, sum)

src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C

@@ -927,6 +927,19 @@ Foam::word Foam::GeometricField<Type, PatchField, GeoMesh>::select
 }
 
 
+template<class Type, template<class> class PatchField, class GeoMesh>
+void Foam::GeometricField<Type, PatchField, GeoMesh>::writeMinMax
+(
+    Ostream& os
+) const
+{
+    os  << "min/max(" << this->name() << ") = "
+        << Foam::min(*this).value() << ", "
+        << Foam::max(*this).value()
+        << endl;
+}
+
+
 // writeData member function required by regIOobject
 template<class Type, template<class> class PatchField, class GeoMesh>
 bool Foam::GeometricField<Type, PatchField, GeoMesh>::

src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -488,6 +488,9 @@ public:
         //  otherwise the standard parameters by returning the field name
         word select(bool final) const;
 
+        //- Helper function to write the min and max to an Ostream
+        void writeMinMax(Ostream& os) const;
+
 
     // Member function *this operators
 

src/OpenFOAM/fields/pointPatchFields/derived/timeVaryingUniformFixedValue/timeVaryingUniformFixedValuePointPatchField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -55,7 +55,9 @@ timeVaryingUniformFixedValuePointPatchField
 :
     fixedValuePointPatchField<Type>(ptf, p, iF, mapper),
     timeSeries_(ptf.timeSeries_)
-{}
+{
+    updateCoeffs();
+}
 
 
 template<class Type>

src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C

@@ -1963,7 +1963,7 @@ Foam::pointField Foam::globalMeshData::geometricSharedPoints() const
     pointField sharedPoints(mesh_.points(), sharedPointLabels());
 
     // Append from all processors
-    combineReduce(sharedPoints, plusEqOp<pointField>());
+    combineReduce(sharedPoints, ListPlusEqOp<pointField>());
 
     // Merge tolerance
     scalar tolDim = matchTol_ * mesh_.bounds().mag();

src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H

@@ -109,29 +109,6 @@ class globalMeshData
     public processorTopology
 {
 
-    // Private class
-
-        // To combineReduce a pointField. Just appends all lists.
-        template<class T>
-        class plusEqOp
-        {
-
-        public:
-
-            void operator()(T& x, const T& y) const
-            {
-                label n = x.size();
-
-                x.setSize(x.size() + y.size());
-
-                forAll(y, i)
-                {
-                    x[n++] = y[i];
-                }
-            }
-        };
-
-
     // Private data
 
         //- Reference to mesh
@@ -337,6 +314,29 @@ class globalMeshData
 
 public:
 
+    // Public class
+
+        // To combineReduce a List. Just appends all lists.
+        template<class T>
+        class ListPlusEqOp
+        {
+
+        public:
+
+            void operator()(T& x, const T& y) const
+            {
+                label n = x.size();
+
+                x.setSize(x.size() + y.size());
+
+                forAll(y, i)
+                {
+                    x[n++] = y[i];
+                }
+            }
+        };
+
+
     //- Runtime type information
     ClassName("globalMeshData");
 

src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H

@@ -251,27 +251,35 @@ public:
         {
         public:
 
-            template<class T>
+            template<class Type>
             void operator()
             (
                 const vectorTensorTransform& vt,
                 const bool forward,
-                List<T>& fld
+                List<Type>& fld
             ) const
             {
-                if (forward)
-                {
-                    transformList(vt.R(), fld);
+                const tensor T(forward ? vt.R() : vt.R().T());
+                transformList(T, fld);
             }
-                else
+
+            template<class Type>
+            void operator()
+            (
+                const vectorTensorTransform& vt,
+                const bool forward,
+                List<List<Type> >& flds
+            ) const
             {
-                    transformList(vt.R().T(), fld);
+                forAll(flds, i)
+                {
+                    operator()(vt, forward, flds[i]);
                 }
             }
 
             //- Transform patch-based field
-            template<class T>
-            void operator()(const coupledPolyPatch& cpp, UList<T>& fld) const
+            template<class Type>
+            void operator()(const coupledPolyPatch& cpp, UList<Type>& fld) const
             {
                 if (!cpp.parallel())
                 {
@@ -280,8 +288,8 @@ public:
             }
 
             //- Transform sparse field
-            template<class T, template<class> class Container>
-            void operator()(const coupledPolyPatch& cpp, Container<T>& map)
+            template<class Type, template<class> class Container>
+            void operator()(const coupledPolyPatch& cpp, Container<Type>& map)
             const
             {
                 if (!cpp.parallel())
@@ -313,6 +321,18 @@ public:
                     fld = vt.invTransformPosition(pfld);
                 }
             }
+            void operator()
+            (
+                const vectorTensorTransform& vt,
+                const bool forward,
+                List<List<point> >& flds
+            ) const
+            {
+                forAll(flds, i)
+                {
+                    operator()(vt, forward, flds[i]);
+                }
+            }
             //- Transform patch-based field
             void operator()(const coupledPolyPatch& cpp, pointField& fld) const
             {

src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsNormals.C

@@ -28,43 +28,6 @@ License
 #include "indirectPrimitivePatch.H"
 #include "globalMeshData.H"
 
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    //- Transformation
-    class listTransform
-    {
-    public:
-
-        void operator()
-        (
-            const vectorTensorTransform& vt,
-            const bool forward,
-            List<List<point> >& fld
-        ) const
-        {
-            const tensor T
-            (
-                forward
-              ? vt.R()
-              : vt.R().T()
-            );
-
-            forAll(fld, i)
-            {
-                List<point>& elems = fld[i];
-                forAll(elems, elemI)
-                {
-                    elems[elemI] = transform(T, elems[elemI]);
-                }
-            }
-        }
-    };
-}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template
@@ -141,7 +104,7 @@ Foam::PatchTools::pointNormals
     (
         transforms,
         pointFaceNormals,
-        listTransform()
+        mapDistribute::transform()
     );
 
 

src/dynamicMesh/layerAdditionRemoval/addCellLayer.C

@@ -161,6 +161,9 @@ void Foam::layerAdditionRemoval::addCellLayer
 
     forAll(mf, faceI)
     {
+        label cellI = mc[faceI];
+        label zoneI =  mesh.cellZones().whichZone(cellI);
+
         addedCells[faceI] =
             ref.setAction
             (
@@ -170,7 +173,7 @@ void Foam::layerAdditionRemoval::addCellLayer
                     -1,           // master edge
                     mf[faceI],    // master face
                     -1,           // master cell
-                    -1            // zone for cell
+                    zoneI         // zone for cell
                 )
             );
     }

src/dynamicMesh/polyTopoChange/polyTopoChange/removePoints.C

@@ -219,7 +219,7 @@ Foam::label Foam::removePoints::countPointUsage
 
     pointCanBeDeleted.setSize(mesh_.nPoints());
     pointCanBeDeleted = false;
-    label nDeleted = 0;
+    //label nDeleted = 0;
 
     forAll(edge0, pointI)
     {
@@ -243,14 +243,14 @@ Foam::label Foam::removePoints::countPointUsage
             if ((e0Vec & e1Vec) > minCos)
             {
                 pointCanBeDeleted[pointI] = true;
-                nDeleted++;
+                //nDeleted++;
             }
         }
         else if (edge0[pointI] == -1)
         {
             // point not used at all
             pointCanBeDeleted[pointI] = true;
-            nDeleted++;
+            //nDeleted++;
         }
     }
     edge0.clear();
@@ -300,6 +300,15 @@ Foam::label Foam::removePoints::countPointUsage
         true                // null value
     );
 
+    label nDeleted = 0;
+    forAll(pointCanBeDeleted, pointI)
+    {
+        if (pointCanBeDeleted[pointI])
+        {
+            nDeleted++;
+        }
+    }
+
     return returnReduce(nDeleted, sumOp<label>());
 }
 

src/finiteVolume/fields/fvPatchFields/derived/mappedField/mappedPatchFieldBase.C

@@ -165,21 +165,6 @@ mappedPatchFieldBase<Type>::sampleField() const
 }
 
 
-template<class Type>
-const interpolation<Type>& mappedPatchFieldBase<Type>::interpolator() const
-{
-    if (!interpolator_.valid())
-    {
-        interpolator_ = interpolation<Type>::New
-        (
-            interpolationScheme_,
-            sampleField()
-        );
-    }
-    return interpolator_();
-}
-
-
 template<class Type>
 tmp<Field<Type> > mappedPatchFieldBase<Type>::mappedField() const
 {
@@ -218,6 +203,14 @@ tmp<Field<Type> > mappedPatchFieldBase<Type>::mappedField() const
                     samples
                 );
 
+                autoPtr<interpolation<Type> > interpolator
+                (
+                    interpolation<Type>::New
+                    (
+                        interpolationScheme_,
+                        sampleField()
+                    )
+                );
                 const interpolation<Type>& interp = interpolator();
 
                 newValues.setSize(samples.size(), pTraits<Type>::max);

src/finiteVolume/fields/fvPatchFields/derived/mappedField/mappedPatchFieldBase.H

@@ -88,11 +88,6 @@ protected:
         //- Interpolation scheme to use for nearestcell mode
         word interpolationScheme_;
 
-        mutable autoPtr<interpolation<Type> > interpolator_;
-
-
-    // Protected Member Functions
-
 
 public:
 
@@ -149,9 +144,6 @@ public:
         //- Field to sample. Either on my or nbr mesh
         const GeometricField<Type, fvPatchField, volMesh>& sampleField() const;
 
-        //- Access the interpolation method
-        const interpolation<Type>& interpolator() const;
-
         //- Map sampleField onto *this patch
         virtual tmp<Field<Type> > mappedField() const;
 

src/finiteVolume/fvMatrices/fvScalarMatrix/fvScalarMatrix.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -36,6 +36,7 @@ SourceFiles
 #define fvScalarMatrix_H
 
 #include "fvMatrix.H"
+#include "fvMatricesFwd.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/finiteVolume/interpolation/interpolation/interpolationCellPoint/interpolationCellPoint.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,7 +35,15 @@ Foam::interpolationCellPoint<Type>::interpolationCellPoint
 )
 :
     interpolation<Type>(psi),
-    psip_(volPointInterpolation::New(psi.mesh()).interpolate(psi))
+    psip_
+    (
+        volPointInterpolation::New(psi.mesh()).interpolate
+        (
+            psi,
+            "volPointInterpolate(" + psi.name() + ')',
+            true        // use cache
+        )
+    )
 {
     // Uses cellPointWeight to do interpolation which needs tet decomposition
     (void)psi.mesh().tetBasePtIs();

src/finiteVolume/interpolation/interpolation/interpolationCellPointFace/interpolationCellPointFace.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -45,7 +45,15 @@ interpolationCellPointFace<Type>::interpolationCellPointFace
 )
 :
     interpolation<Type>(psi),
-    psip_(volPointInterpolation::New(psi.mesh()).interpolate(psi)),
+    psip_
+    (
+        volPointInterpolation::New(psi.mesh()).interpolate
+        (
+            psi,
+            "volPointInterpolate(" + psi.name() + ')',
+            true        // use cache
+        )
+    ),
     psis_(linearInterpolate(psi))
 {}
 

src/finiteVolume/interpolation/interpolation/interpolationPointMVC/interpolationPointMVC.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,7 +35,15 @@ Foam::interpolationPointMVC<Type>::interpolationPointMVC
 )
 :
     interpolation<Type>(psi),
-    psip_(volPointInterpolation::New(psi.mesh()).interpolate(psi))
+    psip_
+    (
+        volPointInterpolation::New(psi.mesh()).interpolate
+        (
+            psi,
+            "volPointInterpolate(" + psi.name() + ')',
+            true            // use cache
+        )
+    )
 {}
 
 

src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/LimitedScheme.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -28,33 +28,17 @@ License
 #include "fvcGrad.H"
 #include "coupledFvPatchFields.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
 template<class Type, class Limiter, template<class> class LimitFunc>
-Foam::tmp<Foam::surfaceScalarField>
-Foam::LimitedScheme<Type, Limiter, LimitFunc>::limiter
+void Foam::LimitedScheme<Type, Limiter, LimitFunc>::calcLimiter
 (
-    const GeometricField<Type, fvPatchField, volMesh>& phi
+    const GeometricField<Type, fvPatchField, volMesh>& phi,
+    surfaceScalarField& limiterField
 ) const
 {
     const fvMesh& mesh = this->mesh();
 
-    tmp<surfaceScalarField> tLimiter
-    (
-        new surfaceScalarField
-        (
-            IOobject
-            (
-                type() + "Limiter(" + phi.name() + ')',
-                mesh.time().timeName(),
-                mesh
-            ),
-            mesh,
-            dimless
-        )
-    );
-    surfaceScalarField& lim = tLimiter();
-
     tmp<GeometricField<typename Limiter::phiType, fvPatchField, volMesh> >
         tlPhi = LimitFunc<Type>()(phi);
 
@@ -73,7 +57,7 @@ Foam::LimitedScheme<Type, Limiter, LimitFunc>::limiter
 
     const vectorField& C = mesh.C();
 
-    scalarField& pLim = lim.internalField();
+    scalarField& pLim = limiterField.internalField();
 
     forAll(pLim, face)
     {
@@ -92,7 +76,8 @@ Foam::LimitedScheme<Type, Limiter, LimitFunc>::limiter
         );
     }
 
-    surfaceScalarField::GeometricBoundaryField& bLim = lim.boundaryField();
+    surfaceScalarField::GeometricBoundaryField& bLim =
+        limiterField.boundaryField();
 
     forAll(bLim, patchi)
     {
@@ -143,8 +128,80 @@ Foam::LimitedScheme<Type, Limiter, LimitFunc>::limiter
             pLim = 1.0;
         }
     }
+}
 
-    return tLimiter;
+
+// * * * * * * * * * * * * Public Member Functions  * * * * * * * * * * * * //
+
+template<class Type, class Limiter, template<class> class LimitFunc>
+Foam::tmp<Foam::surfaceScalarField>
+Foam::LimitedScheme<Type, Limiter, LimitFunc>::limiter
+(
+    const GeometricField<Type, fvPatchField, volMesh>& phi
+) const
+{
+    const fvMesh& mesh = this->mesh();
+
+    const word limiterFieldName(type() + "Limiter(" + phi.name() + ')');
+
+    if (this->mesh().cache("limiter"))
+    {
+        if (!mesh.foundObject<surfaceScalarField>(limiterFieldName))
+        {
+            surfaceScalarField* limiterField
+            (
+                new surfaceScalarField
+                (
+                    IOobject
+                    (
+                        limiterFieldName,
+                        mesh.time().timeName(),
+                        mesh,
+                        IOobject::NO_READ,
+                        IOobject::NO_WRITE
+                    ),
+                    mesh,
+                    dimless
+                )
+            );
+
+            mesh.objectRegistry::store(limiterField);
+        }
+
+        surfaceScalarField& limiterField =
+            const_cast<surfaceScalarField&>
+            (
+                mesh.lookupObject<surfaceScalarField>
+                (
+                    limiterFieldName
+                )
+            );
+
+        calcLimiter(phi, limiterField);
+
+        return limiterField;
+    }
+    else
+    {
+        tmp<surfaceScalarField> tlimiterField
+        (
+            new surfaceScalarField
+            (
+                IOobject
+                (
+                    limiterFieldName,
+                    mesh.time().timeName(),
+                    mesh
+                ),
+                mesh,
+                dimless
+            )
+        );
+
+        calcLimiter(phi, tlimiterField());
+
+        return tlimiterField;
+    }
 }
 
 

src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/LimitedScheme.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -68,6 +68,13 @@ class LimitedScheme
 {
     // Private Member Functions
 
+        //- Calculate the limiter
+        void calcLimiter
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& phi,
+            surfaceScalarField& limiterField
+        ) const;
+
         //- Disallow default bitwise copy construct
         LimitedScheme(const LimitedScheme&);
 

src/fvOptions/sources/derived/effectivenessHeatExchangerSource/effectivenessHeatExchangerSource.H

@@ -82,6 +82,41 @@ Description
     }
     \endverbatim
 
+    The effectiveness table is described in terms of the primary and secondary
+    mass flow rates.  For example, the table:
+
+                           secondary MFR
+                       |  0.1   0.2   0.3
+                  -----+-----------------
+                  0.02 |   A     B     C
+     primary MFR  0.04 |   D     E     F
+                  0.06 |   G     H     I
+
+
+    Is specified by the following:
+
+        (
+            0.02
+            (
+                (0.1    A)
+                (0.2    B)
+                (0.3    C)
+            ),
+            0.04
+            (
+                (0.1    D)
+                (0.2    E)
+                (0.3    F)
+            ),
+            0.06
+            (
+                (0.1    G)
+                (0.2    H)
+                (0.3    I)
+            )
+        );
+
+
 Note
 - the table with name "fileName" should have the same units as the
   secondary mass flow rate and kg/s for phi

src/fvOptions/sources/general/semiImplicitSource/SemiImplicitSource.C

@@ -189,7 +189,7 @@ void Foam::fv::SemiImplicitSource<Type>::addSup
 
     UIndirectList<scalar>(Sp, cells_) = injectionRate_[fieldI].second()/VDash_;
 
-    eqn += Su + fvm::Sp(Sp, psi);
+    eqn += Su + fvm::SuSp(Sp, psi);
 }
 
 

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.C

@@ -110,7 +110,7 @@ void Foam::cloudSolution::read()
     dict_.lookup("transient") >> transient_;
     dict_.lookup("coupled") >> coupled_;
     dict_.lookup("cellValueSourceCorrection") >> cellValueSourceCorrection_;
-    dict_.lookup("maxCo") >> maxCo_;
+    dict_.readIfPresent("maxCo", maxCo_);
 
     if (steadyState())
     {

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ConeInjection/ConeInjection.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -31,7 +31,7 @@ Description
       - list of injector positions and directions (along injection axes)
       - number of parcels to inject per injector
       - parcel velocities
-      - inner and outer cone angles
+      - inner and outer half-cone angles
     - Parcel diameters obtained by distribution model
 
 SourceFiles
@@ -87,10 +87,10 @@ class ConeInjection
         //- Parcel velocity magnitude relative to SOI [m/s]
         const TimeDataEntry<scalar> Umag_;
 
-        //- Inner cone angle relative to SOI [deg]
+        //- Inner half-cone angle relative to SOI [deg]
         const TimeDataEntry<scalar> thetaInner_;
 
-        //- Outer cone angle relative to SOI [deg]
+        //- Outer half-cone angle relative to SOI [deg]
         const TimeDataEntry<scalar> thetaOuter_;
 
         //- Parcel size distribution model

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ConeNozzleInjection/ConeNozzleInjection.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -32,7 +32,7 @@ Description
       - injector position
       - direction (along injection axis)
       - parcel flow rate
-      - inner and outer cone angles
+      - inner and outer half-cone angles
 
     - Parcel diameters obtained by size distribution model
 
@@ -134,10 +134,10 @@ private:
         //- Flow rate profile relative to SOI []
         const TimeDataEntry<scalar> flowRateProfile_;
 
-        //- Inner cone angle relative to SOI [deg]
+        //- Inner half-cone angle relative to SOI [deg]
         const TimeDataEntry<scalar> thetaInner_;
 
-        //- Outer cone angle relative to SOI [deg]
+        //- Outer half-cone angle relative to SOI [deg]
         const TimeDataEntry<scalar> thetaOuter_;
 
         //- Parcel size PDF model

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C

@@ -550,7 +550,7 @@ void Foam::InjectionModel<CloudType>::inject(TrackData& td)
                 if (cellI > -1)
                 {
                     // Lagrangian timestep
-                    scalar dt = time - timeInj;
+                    const scalar dt = time - timeInj;
 
                     // Apply corrections to position for 2-D cases
                     meshTools::constrainToMeshCentre(mesh, pos);
@@ -586,8 +586,11 @@ void Foam::InjectionModel<CloudType>::inject(TrackData& td)
                             pPtr->rho()
                         );
 
+                    const scalar mParcel0 = pPtr->nParticle()*pPtr->mass();
+
                     if (!pPtr->move(td, dt))
                     {
+                        massAdded += mParcel0;
                         delete pPtr;
                     }
                     else
@@ -595,7 +598,7 @@ void Foam::InjectionModel<CloudType>::inject(TrackData& td)
                         if (pPtr->nParticle() >= 1.0)
                         {
                             td.cloud().addParticle(pPtr);
-                            massAdded += pPtr->nParticle()*pPtr->mass();
+                            massAdded += mParcel0;
                             parcelsAdded++;
                         }
                         else

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C

@@ -1262,6 +1262,9 @@ void Foam::autoLayerDriver::calculateLayerThickness
             << setf(ios_base::left) << setw(maxPatchNameLen) << "-----"
             << setw(0) << " -----    ------ --------- -------" << endl;
 
+
+        const PackedBoolList isMasterPoint(syncTools::getMasterPoints(mesh));
+
         forAll(patchIDs, i)
         {
             label patchI = patchIDs[i];
@@ -1270,10 +1273,14 @@ void Foam::autoLayerDriver::calculateLayerThickness
 
             scalar sumThickness = 0;
             scalar sumNearWallThickness = 0;
+            label nMasterPoints = 0;
 
             forAll(meshPoints, patchPointI)
             {
-                label ppPointI = pp.meshPointMap()[meshPoints[patchPointI]];
+                label meshPointI = meshPoints[patchPointI];
+                if (isMasterPoint[meshPointI])
+                {
+                    label ppPointI = pp.meshPointMap()[meshPointI];
 
                     sumThickness += thickness[ppPointI];
                     sumNearWallThickness += layerParams.firstLayerThickness
@@ -1284,9 +1291,11 @@ void Foam::autoLayerDriver::calculateLayerThickness
                         thickness[ppPointI],
                         expansionRatio[ppPointI]
                     );
+                    nMasterPoints++;
+                }
             }
 
-            label totNPoints = returnReduce(meshPoints.size(), sumOp<label>());
+            label totNPoints = returnReduce(nMasterPoints, sumOp<label>());
 
             // For empty patches, totNPoints is 0.
             scalar avgThickness = 0;
@@ -1319,7 +1328,7 @@ void Foam::autoLayerDriver::calculateLayerThickness
 // Synchronize displacement among coupled patches.
 void Foam::autoLayerDriver::syncPatchDisplacement
 (
-    const motionSmoother& meshMover,
+    const indirectPrimitivePatch& pp,
     const scalarField& minThickness,
     pointField& patchDisp,
     labelList& patchNLayers,
@@ -1327,7 +1336,7 @@ void Foam::autoLayerDriver::syncPatchDisplacement
 ) const
 {
     const fvMesh& mesh = meshRefiner_.mesh();
-    const labelList& meshPoints = meshMover.patch().meshPoints();
+    const labelList& meshPoints = pp.meshPoints();
 
     label nChangedTotal = 0;
 
@@ -1437,9 +1446,9 @@ void Foam::autoLayerDriver::syncPatchDisplacement
         }
     }
 
-    Info<< "Prevented extrusion on "
-        << returnReduce(nChangedTotal, sumOp<label>())
-        << " coupled patch points during syncPatchDisplacement." << endl;
+    //Info<< "Prevented extrusion on "
+    //    << returnReduce(nChangedTotal, sumOp<label>())
+    //    << " coupled patch points during syncPatchDisplacement." << endl;
 }
 
 
@@ -1569,7 +1578,7 @@ void Foam::autoLayerDriver::getPatchDisplacement
     // Make sure displacement is equal on both sides of coupled patches.
     syncPatchDisplacement
     (
-        meshMover,
+        pp,
         minThickness,
         patchDisp,
         patchNLayers,
@@ -1767,7 +1776,7 @@ Foam::label Foam::autoLayerDriver::truncateDisplacement
     {
         syncPatchDisplacement
         (
-            meshMover,
+            pp,
             minThickness,
             patchDisp,
             patchNLayers,
@@ -2856,7 +2865,7 @@ void Foam::autoLayerDriver::addLayers
         // Make sure displacement is equal on both sides of coupled patches.
         syncPatchDisplacement
         (
-            meshMover,
+            pp,
             minThickness,
             patchDisp,
             patchNLayers,

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H

@@ -257,7 +257,7 @@ class autoLayerDriver
                 //- Synchronize displacement among coupled patches.
                 void syncPatchDisplacement
                 (
-                    const motionSmoother& meshMover,
+                    const indirectPrimitivePatch& pp,
                     const scalarField& minThickness,
                     pointField& patchDisp,
                     labelList& patchNLayers,
@@ -403,6 +403,7 @@ class autoLayerDriver
                 void smoothField
                 (
                     const motionSmoother& meshMover,
+                    const PackedBoolList& isMasterPoint,
                     const PackedBoolList& isMasterEdge,
                     const labelList& meshEdges,
                     const scalarField& fieldMin,
@@ -414,6 +415,7 @@ class autoLayerDriver
                 void smoothPatchNormals
                 (
                     const motionSmoother& meshMover,
+                    const PackedBoolList& isMasterPoint,
                     const PackedBoolList& isMasterEdge,
                     const labelList& meshEdges,
                     const label nSmoothDisp,
@@ -424,6 +426,7 @@ class autoLayerDriver
                 void smoothNormals
                 (
                     const label nSmoothDisp,
+                    const PackedBoolList& isMasterPoint,
                     const PackedBoolList& isMasterEdge,
                     const labelList& fixedPoints,
                     pointVectorField& normals
@@ -440,8 +443,11 @@ class autoLayerDriver
                 //  large feature angle
                 void handleFeatureAngleLayerTerminations
                 (
-                    const indirectPrimitivePatch& pp,
                     const scalar minCos,
+                    const PackedBoolList& isMasterPoint,
+                    const indirectPrimitivePatch& pp,
+                    const labelList& meshEdges,
+
                     List<extrudeMode>& extrudeStatus,
                     pointField& patchDisp,
                     labelList& patchNLayers,
@@ -453,10 +459,13 @@ class autoLayerDriver
                 // in the layer mesh and stop any layer growth at these points.
                 void findIsolatedRegions
                 (
-                    const indirectPrimitivePatch& pp,
-                    const PackedBoolList& isMasterEdge,
-                    const labelList& meshEdges,
                     const scalar minCosLayerTermination,
+                    const PackedBoolList& isMasterPoint,
+                    const PackedBoolList& isMasterEdge,
+                    const indirectPrimitivePatch& pp,
+                    const labelList& meshEdges,
+                    const scalarField& minThickness,
+
                     List<extrudeMode>& extrudeStatus,
                     pointField& patchDisp,
                     labelList& patchNLayers

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C

@@ -103,6 +103,7 @@ void Foam::autoLayerDriver::sumWeights
 void Foam::autoLayerDriver::smoothField
 (
     const motionSmoother& meshMover,
+    const PackedBoolList& isMasterPoint,
     const PackedBoolList& isMasterEdge,
     const labelList& meshEdges,
     const scalarField& fieldMin,
@@ -136,7 +137,7 @@ void Foam::autoLayerDriver::smoothField
             isMasterEdge,
             meshEdges,
             meshPoints,
-            pp.edges(),
+            edges,
             invSumWeight,
             field,
             average
@@ -162,10 +163,14 @@ void Foam::autoLayerDriver::smoothField
         // Do residual calculation every so often.
         if ((iter % 10) == 0)
         {
-            Info<< "    Iteration " << iter << "   residual "
-                <<  gSum(mag(field-average))
-                   /returnReduce(average.size(), sumOp<label>())
-                << endl;
+            scalar resid = meshRefinement::gAverage
+            (
+                meshMover.mesh(),
+                isMasterPoint,
+                meshPoints,
+                mag(field-average)()
+            );
+            Info<< "    Iteration " << iter << "   residual " << resid << endl;
         }
     }
 }
@@ -265,6 +270,7 @@ void Foam::autoLayerDriver::smoothField
 void Foam::autoLayerDriver::smoothPatchNormals
 (
     const motionSmoother& meshMover,
+    const PackedBoolList& isMasterPoint,
     const PackedBoolList& isMasterEdge,
     const labelList& meshEdges,
     const label nSmoothDisp,
@@ -299,7 +305,7 @@ void Foam::autoLayerDriver::smoothPatchNormals
             isMasterEdge,
             meshEdges,
             meshPoints,
-            pp.edges(),
+            edges,
             invSumWeight,
             normals,
             average
@@ -308,10 +314,14 @@ void Foam::autoLayerDriver::smoothPatchNormals
         // Do residual calculation every so often.
         if ((iter % 10) == 0)
         {
-            Info<< "    Iteration " << iter << "   residual "
-                <<  gSum(mag(normals-average))
-                   /returnReduce(average.size(), sumOp<label>())
-                << endl;
+            scalar resid = meshRefinement::gAverage
+            (
+                meshMover.mesh(),
+                isMasterPoint,
+                meshPoints,
+                mag(normals-average)()
+            );
+            Info<< "    Iteration " << iter << "   residual " << resid << endl;
         }
 
         // Transfer to normals vector field
@@ -330,6 +340,7 @@ void Foam::autoLayerDriver::smoothPatchNormals
 void Foam::autoLayerDriver::smoothNormals
 (
     const label nSmoothDisp,
+    const PackedBoolList& isMasterPoint,
     const PackedBoolList& isMasterEdge,
     const labelList& fixedPoints,
     pointVectorField& normals
@@ -372,8 +383,6 @@ void Foam::autoLayerDriver::smoothNormals
         invSumWeight
     );
 
-    Info<< "shrinkMeshDistance : Smoothing normals in interior ..." << endl;
-
     for (label iter = 0; iter < nSmoothDisp; iter++)
     {
         vectorField average(mesh.nPoints());
@@ -392,10 +401,13 @@ void Foam::autoLayerDriver::smoothNormals
         // Do residual calculation every so often.
         if ((iter % 10) == 0)
         {
-            Info<< "    Iteration " << iter << "   residual "
-                <<  gSum(mag(normals-average))
-                   /returnReduce(average.size(), sumOp<label>())
-                << endl;
+            scalar resid = meshRefinement::gAverage
+            (
+                mesh,
+                isMasterPoint,
+                mag(normals-average)()
+            );
+            Info<< "    Iteration " << iter << "   residual " << resid << endl;
         }
 
 
@@ -479,14 +491,18 @@ bool Foam::autoLayerDriver::isMaxEdge
 // large feature angle
 void Foam::autoLayerDriver::handleFeatureAngleLayerTerminations
 (
-    const indirectPrimitivePatch& pp,
     const scalar minCos,
+    const PackedBoolList& isMasterPoint,
+    const indirectPrimitivePatch& pp,
+    const labelList& meshEdges,
     List<extrudeMode>& extrudeStatus,
     pointField& patchDisp,
     labelList& patchNLayers,
     label& nPointCounter
 ) const
 {
+    const fvMesh& mesh = meshRefiner_.mesh();
+
     // Mark faces that have all points extruded
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -507,15 +523,60 @@ void Foam::autoLayerDriver::handleFeatureAngleLayerTerminations
     }
 
 
+
+    //label nOldPointCounter = nPointCounter;
+
     // Detect situation where two featureedge-neighbouring faces are partly or
     // not extruded and the edge itself is extruded. In this case unmark the
     // edge for extrusion.
 
-    forAll(pp.edgeFaces(), edgeI)
-    {
-        const labelList& eFaces = pp.edgeFaces()[edgeI];
 
-        if (eFaces.size() == 2)
+    List<List<point> > edgeFaceNormals(pp.nEdges());
+    List<List<bool> > edgeFaceExtrude(pp.nEdges());
+
+    const labelListList& edgeFaces = pp.edgeFaces();
+    const vectorField& faceNormals = pp.faceNormals();
+    const labelList& meshPoints = pp.meshPoints();
+
+    forAll(edgeFaces, edgeI)
+    {
+        const labelList& eFaces = edgeFaces[edgeI];
+
+        edgeFaceNormals[edgeI].setSize(eFaces.size());
+        edgeFaceExtrude[edgeI].setSize(eFaces.size());
+        forAll(eFaces, i)
+        {
+            label faceI = eFaces[i];
+            edgeFaceNormals[edgeI][i] = faceNormals[faceI];
+            edgeFaceExtrude[edgeI][i] = extrudedFaces[faceI];
+        }
+    }
+
+    syncTools::syncEdgeList
+    (
+        mesh,
+        meshEdges,
+        edgeFaceNormals,
+        globalMeshData::ListPlusEqOp<List<point> >(),   // combine operator
+        List<point>()               // null value
+    );
+
+    syncTools::syncEdgeList
+    (
+        mesh,
+        meshEdges,
+        edgeFaceExtrude,
+        globalMeshData::ListPlusEqOp<List<bool> >(),    // combine operator
+        List<bool>()                // null value
+    );
+
+
+    forAll(edgeFaceNormals, edgeI)
+    {
+        const List<point>& eFaceNormals = edgeFaceNormals[edgeI];
+        const List<bool>& eFaceExtrude = edgeFaceExtrude[edgeI];
+
+        if (eFaceNormals.size() == 2)
         {
             const edge& e = pp.edges()[edgeI];
             label v0 = e[0];
@@ -527,10 +588,10 @@ void Foam::autoLayerDriver::handleFeatureAngleLayerTerminations
              || extrudeStatus[v1] != NOEXTRUDE
             )
             {
-                if (!extrudedFaces[eFaces[0]] || !extrudedFaces[eFaces[1]])
+                if (!eFaceExtrude[0] || !eFaceExtrude[1])
                 {
-                    const vector& n0 = pp.faceNormals()[eFaces[0]];
-                    const vector& n1 = pp.faceNormals()[eFaces[1]];
+                    const vector& n0 = eFaceNormals[0];
+                    const vector& n1 = eFaceNormals[1];
 
                     if ((n0 & n1) < minCos)
                     {
@@ -544,9 +605,12 @@ void Foam::autoLayerDriver::handleFeatureAngleLayerTerminations
                                 extrudeStatus
                             )
                         )
+                        {
+                            if (isMasterPoint[meshPoints[v0]])
                             {
                                 nPointCounter++;
                             }
+                        }
                         if
                         (
                             unmarkExtrusion
@@ -557,6 +621,8 @@ void Foam::autoLayerDriver::handleFeatureAngleLayerTerminations
                                 extrudeStatus
                             )
                         )
+                        {
+                            if (isMasterPoint[meshPoints[v1]])
                             {
                                 nPointCounter++;
                             }
@@ -567,15 +633,22 @@ void Foam::autoLayerDriver::handleFeatureAngleLayerTerminations
         }
     }
 
+    //Info<< "Added "
+    //    << returnReduce(nPointCounter-nOldPointCounter, sumOp<label>())
+    //    << " point not to extrude." << endl;
+}
+
 
 // Find isolated islands (points, edges and faces and layer terminations)
 // in the layer mesh and stop any layer growth at these points.
 void Foam::autoLayerDriver::findIsolatedRegions
 (
-    const indirectPrimitivePatch& pp,
-    const PackedBoolList& isMasterEdge,
-    const labelList& meshEdges,
     const scalar minCosLayerTermination,
+    const PackedBoolList& isMasterPoint,
+    const PackedBoolList& isMasterEdge,
+    const indirectPrimitivePatch& pp,
+    const labelList& meshEdges,
+    const scalarField& minThickness,
     List<extrudeMode>& extrudeStatus,
     pointField& patchDisp,
     labelList& patchNLayers
@@ -594,8 +667,10 @@ void Foam::autoLayerDriver::findIsolatedRegions
         // large feature angle
         handleFeatureAngleLayerTerminations
         (
-            pp,
             minCosLayerTermination,
+            isMasterPoint,
+            pp,
+            meshEdges,
 
             extrudeStatus,
             patchDisp,
@@ -603,6 +678,15 @@ void Foam::autoLayerDriver::findIsolatedRegions
             nPointCounter
         );
 
+        syncPatchDisplacement
+        (
+            pp,
+            minThickness,
+            patchDisp,
+            patchNLayers,
+            extrudeStatus
+        );
+
 
         // Do not extrude from point where all neighbouring
         // faces are not grown
@@ -802,25 +886,18 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo
     // Predetermine mesh edges
     // ~~~~~~~~~~~~~~~~~~~~~~~
 
-    // Precalulate master edge (only relevant for shared edges)
-    PackedBoolList isMasterEdge(syncTools::getMasterEdges(mesh));
+    // Precalulate master point/edge (only relevant for shared points/edges)
+    const PackedBoolList isMasterPoint(syncTools::getMasterPoints(mesh));
+    const PackedBoolList isMasterEdge(syncTools::getMasterEdges(mesh));
     // Precalculate meshEdge per pp edge
-    labelList meshEdges(pp.nEdges());
-
-    forAll(meshEdges, patchEdgeI)
-    {
-        const edge& e = pp.edges()[patchEdgeI];
-
-        label v0 = pp.meshPoints()[e[0]];
-        label v1 = pp.meshPoints()[e[1]];
-        meshEdges[patchEdgeI] = meshTools::findEdge
+    const labelList meshEdges
+    (
+        pp.meshEdges
         (
             mesh.edges(),
-            mesh.pointEdges()[v0],
-            v0,
-            v1
+            mesh.pointEdges()
+        )
     );
-    }
 
 
     // Determine pointNormal
@@ -832,7 +909,7 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo
     if (debug&meshRefinement::MESH || debug&meshRefinement::LAYERINFO)
     {
         pointField meshPointNormals(mesh.nPoints(), point(1, 0, 0));
-        UIndirectList<point>(meshPointNormals, pp.meshPoints()) = pointNormals;
+        UIndirectList<point>(meshPointNormals, meshPoints) = pointNormals;
         meshRefinement::testSyncPointList
         (
             "pointNormals",
@@ -845,6 +922,7 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo
     smoothPatchNormals
     (
         meshMover,
+        isMasterPoint,
         isMasterEdge,
         meshEdges,
         nSmoothSurfaceNormals,
@@ -855,7 +933,7 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo
     if (debug&meshRefinement::MESH || debug&meshRefinement::LAYERINFO)
     {
         pointField meshPointNormals(mesh.nPoints(), point(1, 0, 0));
-        UIndirectList<point>(meshPointNormals, pp.meshPoints()) = pointNormals;
+        UIndirectList<point>(meshPointNormals, meshPoints) = pointNormals;
         meshRefinement::testSyncPointList
         (
             "smoothed pointNormals",
@@ -929,18 +1007,18 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo
     {
         pointField origin(pointWallDist.size());
         scalarField distSqr(pointWallDist.size());
-        scalarField passiveS(pointWallDist.size());
+        //NA scalarField passiveS(pointWallDist.size());
         pointField passiveV(pointWallDist.size());
         forAll(pointWallDist, pointI)
         {
             origin[pointI] = pointWallDist[pointI].origin();
             distSqr[pointI] = pointWallDist[pointI].distSqr();
-            passiveS[pointI] = pointWallDist[pointI].s();
+            //passiveS[pointI] = pointWallDist[pointI].s();
             passiveV[pointI] = pointWallDist[pointI].v();
         }
         meshRefinement::testSyncPointList("origin", mesh, origin);
         meshRefinement::testSyncPointList("distSqr", mesh, distSqr);
-        meshRefinement::testSyncPointList("passiveS", mesh, passiveS);
+        //meshRefinement::testSyncPointList("passiveS", mesh, passiveS);
         meshRefinement::testSyncPointList("passiveV", mesh, passiveV);
     }
 
@@ -1186,7 +1264,14 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo
     }
 
     // Smooth normal vectors. Do not change normals on pp.meshPoints
-    smoothNormals(nSmoothNormals, isMasterEdge, meshPoints, dispVec);
+    smoothNormals
+    (
+        nSmoothNormals,
+        isMasterPoint,
+        isMasterEdge,
+        meshPoints,
+        dispVec
+    );
 
     if (debug&meshRefinement::MESH || debug&meshRefinement::LAYERINFO)
     {
@@ -1288,25 +1373,18 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance
     const indirectPrimitivePatch& pp = meshMover.patch();
     const labelList& meshPoints = pp.meshPoints();
 
-    // Precalulate master edge (only relevant for shared edges)
-    PackedBoolList isMasterEdge(syncTools::getMasterEdges(mesh));
+    // Precalulate master points/edge (only relevant for shared points/edges)
+    const PackedBoolList isMasterPoint(syncTools::getMasterPoints(mesh));
+    const PackedBoolList isMasterEdge(syncTools::getMasterEdges(mesh));
     // Precalculate meshEdge per pp edge
-    labelList meshEdges(pp.nEdges());
-
-    forAll(meshEdges, patchEdgeI)
-    {
-        const edge& e = pp.edges()[patchEdgeI];
-
-        label v0 = pp.meshPoints()[e[0]];
-        label v1 = pp.meshPoints()[e[1]];
-        meshEdges[patchEdgeI] = meshTools::findEdge
+    const labelList meshEdges
+    (
+        pp.meshEdges
         (
             mesh.edges(),
-            mesh.pointEdges()[v0],
-            v0,
-            v1
+            mesh.pointEdges()
+        )
     );
-    }
 
 
     scalarField thickness(layerThickness.size());
@@ -1405,7 +1483,10 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance
 
                 patchDisp[patchPointI] = thickness[patchPointI]*n;
 
+                if (isMasterPoint[pointI])
+                {
                     numThicknessRatioExclude++;
+                }
 
                 if (str.valid())
                 {
@@ -1433,14 +1514,17 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance
         << numThicknessRatioExclude
         << " nodes where thickness to medial axis distance is large " << endl;
 
+
     // find points where layer growth isolated to a lone point, edge or face
 
     findIsolatedRegions
     (
-        pp,
-        isMasterEdge,
-        meshEdges,
         minCosLayerTermination,
+        isMasterPoint,
+        isMasterEdge,
+        pp,
+        meshEdges,
+        minThickness,
 
         extrudeStatus,
         patchDisp,
@@ -1461,6 +1545,7 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance
     smoothField
     (
         meshMover,
+        isMasterPoint,
         isMasterEdge,
         meshEdges,
         minThickness,
@@ -1608,9 +1693,13 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance
             // Do residual calculation every so often.
             if ((iter % 10) == 0)
             {
-                Info<< "    Iteration " << iter << "   residual "
-                    <<  gSum(mag(displacement-average))
-                       /returnReduce(average.size(), sumOp<label>())
+                scalar resid = meshRefinement::gAverage
+                (
+                    mesh,
+                    syncTools::getMasterPoints(mesh),
+                    mag(displacement-average)()
+                );
+                Info<< "    Iteration " << iter << "   residual " << resid
                     << endl;
             }
         }

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C

@@ -33,6 +33,7 @@ Description
 #include "fvMesh.H"
 #include "Time.H"
 #include "OFstream.H"
+#include "OBJstream.H"
 #include "mapPolyMesh.H"
 #include "pointEdgePoint.H"
 #include "PointEdgeWave.H"
@@ -909,6 +910,7 @@ void Foam::autoSnapDriver::detectNearSurfaces
     Info<< "Detecting near surfaces ..." << endl;
 
     const pointField& localPoints = pp.localPoints();
+    const labelList& meshPoints = pp.meshPoints();
     const refinementSurfaces& surfaces = meshRefiner_.surfaces();
     const fvMesh& mesh = meshRefiner_.mesh();
 
@@ -1193,15 +1195,33 @@ void Foam::autoSnapDriver::detectNearSurfaces
 
     const pointField avgCc(avgCellCentres(mesh, pp));
 
-    // Construct rays from localPoints to beyond cell centre
+    // Construct rays through localPoints to beyond cell centre
+    pointField start(pp.nPoints());
     pointField end(pp.nPoints());
     forAll(localPoints, pointI)
     {
         const point& pt = localPoints[pointI];
-        end[pointI] = pt + 2*(avgCc[pointI]-pt);
+        const vector d = 2*(avgCc[pointI]-pt);
+        start[pointI] = pt - d;
+        end[pointI] = pt + d;
     }
 
 
+    autoPtr<OBJstream> gapStr;
+    if (debug&meshRefinement::OBJINTERSECTIONS)
+    {
+        gapStr.reset
+        (
+            new OBJstream
+            (
+                mesh.time().path()
+              / "detectNearSurfaces_" + meshRefiner_.timeName() + ".obj"
+            )
+        );
+    }
+
+
+    const PackedBoolList isMasterPoint(syncTools::getMasterPoints(mesh));
     label nOverride = 0;
 
     // 1. All points to non-interface surfaces
@@ -1226,7 +1246,7 @@ void Foam::autoSnapDriver::detectNearSurfaces
         surfaces.findNearestIntersection
         (
             unzonedSurfaces,
-            localPoints,
+            start,
             end,
 
             surface1,
@@ -1248,44 +1268,73 @@ void Foam::autoSnapDriver::detectNearSurfaces
 
             bool override = false;
 
-            if (hit1[pointI].hit())
+            //if (hit1[pointI].hit())
+            //{
+            //    if
+            //    (
+            //        meshRefiner_.isGap
+            //        (
+            //            planarCos,
+            //            nearestPoint[pointI],
+            //            nearestNormal[pointI],
+            //            hit1[pointI].hitPoint(),
+            //            normal1[pointI]
+            //        )
+            //    )
+            //    {
+            //        disp[pointI] = hit1[pointI].hitPoint()-pt;
+            //        override = true;
+            //    }
+            //}
+            //if (hit2[pointI].hit())
+            //{
+            //    if
+            //    (
+            //        meshRefiner_.isGap
+            //        (
+            //            planarCos,
+            //            nearestPoint[pointI],
+            //            nearestNormal[pointI],
+            //            hit2[pointI].hitPoint(),
+            //            normal2[pointI]
+            //        )
+            //    )
+            //    {
+            //        disp[pointI] = hit2[pointI].hitPoint()-pt;
+            //        override = true;
+            //    }
+            //}
+
+            if (hit1[pointI].hit() && hit2[pointI].hit())
             {
                 if
                 (
                     meshRefiner_.isGap
                     (
                         planarCos,
-                        nearestPoint[pointI],
-                        nearestNormal[pointI],
                         hit1[pointI].hitPoint(),
-                        normal1[pointI]
-                    )
-                )
-                {
-                    disp[pointI] = hit1[pointI].hitPoint()-pt;
-                    override = true;
-                }
-            }
-            if (hit2[pointI].hit())
-            {
-                if
-                (
-                    meshRefiner_.isGap
-                    (
-                        planarCos,
-                        nearestPoint[pointI],
-                        nearestNormal[pointI],
+                        normal1[pointI],
                         hit2[pointI].hitPoint(),
                         normal2[pointI]
                     )
                 )
                 {
+                    // TBD: check if the attraction (to nearest) would attract
+                    // good enough and not override attraction
+
+                    if (gapStr.valid())
+                    {
+                        const point& intPt = hit2[pointI].hitPoint();
+                        gapStr().write(linePointRef(pt, intPt));
+                    }
+
+                    // Choose hit2 : nearest to end point (so inside the domain)
                     disp[pointI] = hit2[pointI].hitPoint()-pt;
                     override = true;
                 }
             }
 
-            if (override)
+            if (override && isMasterPoint[meshPoints[pointI]])
             {
                 nOverride++;
             }
@@ -1337,7 +1386,7 @@ void Foam::autoSnapDriver::detectNearSurfaces
             surfaces.findNearestIntersection
             (
                 surfacesToTest,
-                pointField(localPoints, zonePointIndices),
+                pointField(start, zonePointIndices),
                 pointField(end, zonePointIndices),
 
                 surface1,
@@ -1352,8 +1401,6 @@ void Foam::autoSnapDriver::detectNearSurfaces
             );
 
 
-            label nOverride = 0;
-
             forAll(hit1, i)
             {
                 label pointI = zonePointIndices[i];
@@ -1363,44 +1410,69 @@ void Foam::autoSnapDriver::detectNearSurfaces
 
                 bool override = false;
 
-                if (hit1[i].hit())
+                //if (hit1[i].hit())
+                //{
+                //    if
+                //    (
+                //        meshRefiner_.isGap
+                //        (
+                //            planarCos,
+                //            nearestPoint[pointI],
+                //            nearestNormal[pointI],
+                //            hit1[i].hitPoint(),
+                //            normal1[i]
+                //        )
+                //    )
+                //    {
+                //        disp[pointI] = hit1[i].hitPoint()-pt;
+                //        override = true;
+                //    }
+                //}
+                //if (hit2[i].hit())
+                //{
+                //    if
+                //    (
+                //        meshRefiner_.isGap
+                //        (
+                //            planarCos,
+                //            nearestPoint[pointI],
+                //            nearestNormal[pointI],
+                //            hit2[i].hitPoint(),
+                //            normal2[i]
+                //        )
+                //    )
+                //    {
+                //        disp[pointI] = hit2[i].hitPoint()-pt;
+                //        override = true;
+                //    }
+                //}
+
+                if (hit1[i].hit() && hit2[i].hit())
                 {
                     if
                     (
                         meshRefiner_.isGap
                         (
                             planarCos,
-                            nearestPoint[pointI],
-                            nearestNormal[pointI],
                             hit1[i].hitPoint(),
-                            normal1[i]
-                        )
-                    )
-                    {
-                        disp[pointI] = hit1[i].hitPoint()-pt;
-                        override = true;
-                    }
-                }
-                if (hit2[i].hit())
-                {
-                    if
-                    (
-                        meshRefiner_.isGap
-                        (
-                            planarCos,
-                            nearestPoint[pointI],
-                            nearestNormal[pointI],
+                            normal1[i],
                             hit2[i].hitPoint(),
                             normal2[i]
                         )
                     )
                     {
+                        if (gapStr.valid())
+                        {
+                            const point& intPt = hit2[i].hitPoint();
+                            gapStr().write(linePointRef(pt, intPt));
+                        }
+
                         disp[pointI] = hit2[i].hitPoint()-pt;
                         override = true;
                     }
                 }
 
-                if (override)
+                if (override && isMasterPoint[meshPoints[pointI]])
                 {
                     nOverride++;
                 }
@@ -1618,7 +1690,13 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
             scalarField magDisp(mag(patchDisp));
 
             Info<< "Wanted displacement : average:"
-                << gSum(magDisp)/returnReduce(patchDisp.size(), sumOp<label>())
+                <<  meshRefinement::gAverage
+                    (
+                        mesh,
+                        syncTools::getMasterPoints(mesh),
+                        pp.meshPoints(),
+                        magDisp
+                    )
                 << " min:" << gMin(magDisp)
                 << " max:" << gMax(magDisp) << endl;
         }
@@ -2476,25 +2554,30 @@ void Foam::autoSnapDriver::doSnap
                 adaptPatchIDs
             )
         );
-        indirectPrimitivePatch& pp = ppPtr();
+
 
         // Distance to attract to nearest feature on surface
-        const scalarField snapDist(calcSnapDistance(mesh, snapParams, pp));
+        const scalarField snapDist(calcSnapDistance(mesh, snapParams, ppPtr()));
 
 
         // Construct iterative mesh mover.
         Info<< "Constructing mesh displacer ..." << endl;
         Info<< "Using mesh parameters " << motionDict << nl << endl;
 
-        const pointMesh& pMesh = pointMesh::New(mesh);
-
-        motionSmoother meshMover
+        autoPtr<motionSmoother> meshMoverPtr
+        (
+            new motionSmoother
             (
                 mesh,
-            pp,
+                ppPtr(),
                 adaptPatchIDs,
-            meshRefinement::makeDisplacementField(pMesh, adaptPatchIDs),
+                meshRefinement::makeDisplacementField
+                (
+                    pointMesh::New(mesh),
+                    adaptPatchIDs
+                ),
                 motionDict
+            )
         );
 
 
@@ -2523,16 +2606,95 @@ void Foam::autoSnapDriver::doSnap
             snapParams,
             nInitErrors,
             baffles,
-            meshMover
+            meshMoverPtr()
         );
 
 
+
+        //- Only if in feature attraction mode:
+        // Nearest feature
+        vectorField patchAttraction;
+        // Constraints at feature
+        List<pointConstraint> patchConstraints;
+
+
         for (label iter = 0; iter < nFeatIter; iter++)
         {
             Info<< nl
                 << "Morph iteration " << iter << nl
                 << "-----------------" << endl;
 
+
+            //if (iter > 0 && iter == nFeatIter/2)
+            //{
+            //    Info<< "Splitting diagonal attractions" << endl;
+            //    const labelList& bFaces = ppPtr().addressing();
+            //
+            //    DynamicList<label> splitFaces(bFaces.size());
+            //    DynamicList<labelPair> splits(bFaces.size());
+            //
+            //    forAll(bFaces, faceI)
+            //    {
+            //        const labelPair split
+            //        (
+            //            findDiagonalAttraction
+            //            (
+            //                ppPtr(),
+            //                patchAttraction,
+            //                patchConstraints,
+            //                faceI
+            //            )
+            //        );
+            //
+            //        if (split != labelPair(-1, -1))
+            //        {
+            //            splitFaces.append(bFaces[faceI]);
+            //            splits.append(split);
+            //        }
+            //    }
+            //
+            //    autoPtr<mapPolyMesh> mapPtr = meshRefiner_.splitFaces
+            //    (
+            //        splitFaces,
+            //        splits
+            //    );
+            //
+            //    const labelList& faceMap = mapPtr().faceMap();
+            //    meshRefinement::updateList(faceMap, -1, duplicateFace);
+            //    const labelList& reverseFaceMap = mapPtr().reverseFaceMap();
+            //    forAll(baffles, i)
+            //    {
+            //        labelPair& baffle = baffles[i];
+            //        baffle.first() = reverseFaceMap[baffle.first()];
+            //        baffle.second() = reverseFaceMap[baffle.second()];
+            //    }
+            //
+            //    meshMoverPtr.clear();
+            //    ppPtr.clear();
+            //
+            //    ppPtr = meshRefinement::makePatch(mesh, adaptPatchIDs);
+            //    meshMoverPtr.reset
+            //    (
+            //        new motionSmoother
+            //        (
+            //            mesh,
+            //            ppPtr(),
+            //            adaptPatchIDs,
+            //            meshRefinement::makeDisplacementField
+            //            (
+            //                pointMesh::New(mesh),
+            //                adaptPatchIDs
+            //            ),
+            //            motionDict
+            //        )
+            //    );
+            //}
+
+
+            indirectPrimitivePatch& pp = ppPtr();
+            motionSmoother& meshMover = meshMoverPtr();
+
+
             // Calculate displacement at every patch point. Insert into
             // meshMover.
             // Calculate displacement at every patch point
@@ -2580,7 +2742,9 @@ void Foam::autoSnapDriver::doSnap
                     scalar(iter+1)/nFeatIter,
                     snapDist,
                     disp,
-                    meshMover
+                    meshMover,
+                    patchAttraction,
+                    patchConstraints
                 );
             }
 

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H

@@ -125,7 +125,9 @@ class autoSnapDriver
 
                 void smoothAndConstrain
                 (
+                    const PackedBoolList& isMasterEdge,
                     const indirectPrimitivePatch& pp,
+                    const labelList& meshEdges,
                     const List<pointConstraint>& constraints,
                     vectorField& disp
                 ) const;
@@ -196,6 +198,16 @@ class autoSnapDriver
                     List<pointConstraint>& patchConstraints
                 ) const;
 
+                //- Detect any diagonal attraction. Returns indices in face
+                //  or (-1, -1) if none
+                labelPair findDiagonalAttraction
+                (
+                    const indirectPrimitivePatch& pp,
+                    const vectorField& patchAttraction,
+                    const List<pointConstraint>& patchConstraints,
+                    const label faceI
+                ) const;
+
                 //- Return hit if on multiple points
                 pointIndexHit findMultiPatchPoint
                 (
@@ -371,7 +383,9 @@ class autoSnapDriver
                     const scalar featureAttract,
                     const scalarField& snapDist,
                     const vectorField& nearestDisp,
-                    motionSmoother& meshMover
+                    motionSmoother& meshMover,
+                    vectorField& patchAttraction,
+                    List<pointConstraint>& patchConstraints
                 ) const;
 
 

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C

@@ -38,50 +38,18 @@ License
 #include "treeDataPoint.H"
 #include "indexedOctree.H"
 #include "snapParameters.H"
+#include "PatchTools.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    class listTransform
-    {
-    public:
-
-        void operator()
-        (
-            const vectorTensorTransform& vt,
-            const bool forward,
-            List<List<point> >& fld
-        ) const
-        {
-            const tensor T
-            (
-                forward
-              ? vt.R()
-              : vt.R().T()
-            );
-
-            forAll(fld, i)
-            {
-                List<point>& elems = fld[i];
-                forAll(elems, elemI)
-                {
-                    elems[elemI] = transform(T, elems[elemI]);
-                }
-            }
-        }
-    };
-
     template<class T>
     class listPlusEqOp
     {
     public:
 
-        void operator()
-        (
-            List<T>& x,
-            const List<T>& y
-        ) const
+        void operator()(List<T>& x, const List<T>& y) const
         {
             label sz = x.size();
             x.setSize(sz+y.size());
@@ -159,7 +127,9 @@ bool Foam::autoSnapDriver::isFeaturePoint
 
 void Foam::autoSnapDriver::smoothAndConstrain
 (
+    const PackedBoolList& isMasterEdge,
     const indirectPrimitivePatch& pp,
+    const labelList& meshEdges,
     const List<pointConstraint>& constraints,
     vectorField& disp
 ) const
@@ -196,6 +166,10 @@ void Foam::autoSnapDriver::smoothAndConstrain
             if (nConstraints <= 1)
             {
                 forAll(pEdges, i)
+                {
+                    label edgeI = pEdges[i];
+
+                    if (isMasterEdge[meshEdges[edgeI]])
                     {
                         label nbrPointI = edges[pEdges[i]].otherVertex(pointI);
                         if (constraints[nbrPointI].first() >= nConstraints)
@@ -206,6 +180,7 @@ void Foam::autoSnapDriver::smoothAndConstrain
                     }
                 }
             }
+        }
 
         syncTools::syncPointList
         (
@@ -564,6 +539,9 @@ void Foam::autoSnapDriver::calcNearestFacePointProperties
 {
     const fvMesh& mesh = meshRefiner_.mesh();
 
+    const PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh));
+
+
     // For now just get all surrounding face data. Expensive - should just
     // store and sync data on coupled points only
     // (see e.g PatchToolsNormals.C)
@@ -583,7 +561,7 @@ void Foam::autoSnapDriver::calcNearestFacePointProperties
         forAll(pFaces, i)
         {
             label faceI = pFaces[i];
-            if (faceSurfaceGlobalRegion[faceI] != -1)
+            if (isMasterFace[faceI] && faceSurfaceGlobalRegion[faceI] != -1)
             {
                 nFaces++;
             }
@@ -605,7 +583,7 @@ void Foam::autoSnapDriver::calcNearestFacePointProperties
             label faceI = pFaces[i];
             label globalRegionI = faceSurfaceGlobalRegion[faceI];
 
-            if (globalRegionI != -1)
+            if (isMasterFace[faceI] && globalRegionI != -1)
             {
                 pNormals[nFaces] = faceSurfaceNormal[faceI];
                 pDisp[nFaces] = faceDisp[faceI];
@@ -694,7 +672,7 @@ void Foam::autoSnapDriver::calcNearestFacePointProperties
         pointFaceSurfNormals,
         listPlusEqOp<point>(),
         List<point>(),
-        listTransform()
+        mapDistribute::transform()
     );
     syncTools::syncPointList
     (
@@ -703,7 +681,7 @@ void Foam::autoSnapDriver::calcNearestFacePointProperties
         pointFaceDisp,
         listPlusEqOp<point>(),
         List<point>(),
-        listTransform()
+        mapDistribute::transform()
     );
     syncTools::syncPointList
     (
@@ -712,7 +690,7 @@ void Foam::autoSnapDriver::calcNearestFacePointProperties
         pointFaceCentres,
         listPlusEqOp<point>(),
         List<point>(),
-        listTransform()
+        mapDistribute::transformPosition()
     );
     syncTools::syncPointList
     (
@@ -722,6 +700,25 @@ void Foam::autoSnapDriver::calcNearestFacePointProperties
         listPlusEqOp<label>(),
         List<label>()
     );
+
+
+    // Sort the data according to the face centres. This is only so we get
+    // consistent behaviour serial and parallel.
+    labelList visitOrder;
+    forAll(pointFaceDisp, pointI)
+    {
+        List<point>& pNormals = pointFaceSurfNormals[pointI];
+        List<point>& pDisp = pointFaceDisp[pointI];
+        List<point>& pFc = pointFaceCentres[pointI];
+        labelList& pFid = pointFacePatchID[pointI];
+
+        sortedOrder(mag(pFc)(), visitOrder);
+
+        pNormals = List<point>(pNormals, visitOrder);
+        pDisp = List<point>(pDisp, visitOrder);
+        pFc = List<point>(pFc, visitOrder);
+        pFid = UIndirectList<label>(pFid, visitOrder);
+    }
 }
 
 
@@ -1360,6 +1357,70 @@ void Foam::autoSnapDriver::stringFeatureEdges
 }
 
 
+// If only two attractions and across face return the face indices
+Foam::labelPair Foam::autoSnapDriver::findDiagonalAttraction
+(
+    const indirectPrimitivePatch& pp,
+    const vectorField& patchAttraction,
+    const List<pointConstraint>& patchConstraints,
+    const label faceI
+) const
+{
+    const face& f = pp.localFaces()[faceI];
+    // For now just detect any attraction. Improve this to look at
+    // actual attraction position and orientation
+
+    labelPair attractIndices(-1, -1);
+
+    if (f.size() >= 4)
+    {
+        forAll(f, fp)
+        {
+            label pointI = f[fp];
+            if (patchConstraints[pointI].first() >= 2)
+            {
+                // Attract to feature edge or point
+                if (attractIndices[0] == -1)
+                {
+                    // First attraction. Store
+                    attractIndices[0] = fp;
+                }
+                else if (attractIndices[1] == -1)
+                {
+                    // Second attraction. Check if not consecutive to first
+                    // attraction
+                    label fp0 = attractIndices[0];
+                    if (f.fcIndex(fp0) == fp || f.fcIndex(fp) == fp0)
+                    {
+                        // Consecutive. Skip.
+                        attractIndices = labelPair(-1, -1);
+                        break;
+                    }
+                    else
+                    {
+                        attractIndices[1] = fp;
+                    }
+                }
+                else
+                {
+                    // More than two attractions. Skip.
+                    attractIndices = labelPair(-1, -1);
+                    break;
+                }
+            }
+        }
+
+
+        if (attractIndices[1] == -1)
+        {
+            // Found only one attraction. Skip.
+            attractIndices = labelPair(-1, -1);
+        }
+    }
+    return attractIndices;
+}
+
+
 Foam::pointIndexHit Foam::autoSnapDriver::findNearFeatureEdge
 (
     const indirectPrimitivePatch& pp,
@@ -1929,7 +1990,7 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
                 edgeFaceNormals,
                 listPlusEqOp<point>(),
                 List<point>(),
-                listTransform()
+                mapDistribute::transform()
             );
         }
 
@@ -2778,7 +2839,9 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
     const scalar featureAttract,
     const scalarField& snapDist,
     const vectorField& nearestDisp,
-    motionSmoother& meshMover
+    motionSmoother& meshMover,
+    vectorField& patchAttraction,
+    List<pointConstraint>& patchConstraints
 ) const
 {
     const Switch implicitFeatureAttraction = snapParams.implicitFeatureSnap();
@@ -2851,7 +2914,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // - faceSurfaceNormal
     // - faceDisp
-    // - faceCentres&faceNormal
+    // - faceCentres
     List<List<point> > pointFaceSurfNormals(pp.nPoints());
     List<List<point> > pointFaceDisp(pp.nPoints());
     List<List<point> > pointFaceCentres(pp.nPoints());
@@ -2884,10 +2947,11 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
     // here.
 
     // Nearest feature
-    vectorField patchAttraction(localPoints.size(), vector::zero);
+    patchAttraction.setSize(localPoints.size());
+    patchAttraction = vector::zero;
     // Constraints at feature
-    List<pointConstraint> patchConstraints(localPoints.size());
-
+    patchConstraints.setSize(localPoints.size());
+    patchConstraints = pointConstraint();
 
     if (implicitFeatureAttraction)
     {
@@ -2951,15 +3015,29 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
         patchConstraints
     );
 
+    const PackedBoolList isMasterPoint(syncTools::getMasterPoints(mesh));
+    {
+        vector avgPatchDisp = meshRefinement::gAverage
+        (
+            mesh,
+            isMasterPoint,
+            pp.meshPoints(),
+            patchDisp
+        );
+        vector avgPatchAttr = meshRefinement::gAverage
+        (
+            mesh,
+            isMasterPoint,
+            pp.meshPoints(),
+            patchAttraction
+        );
 
         Info<< "Attraction:" << endl
-        << "     linear   : max:" << gMax(patchDisp)
-        << " avg:" << gAverage(patchDisp)
-        << endl
-        << "     feature  : max:" << gMax(patchAttraction)
-        << " avg:" << gAverage(patchAttraction)
-        << endl;
-
+            << "     linear   : max:" << gMaxMagSqr(patchDisp)
+            << " avg:" << avgPatchDisp << endl
+            << "     feature  : max:" << gMaxMagSqr(patchAttraction)
+            << " avg:" << avgPatchAttr << endl;
+    }
 
     // So now we have:
     // - patchDisp          : point movement to go to nearest point on surface
@@ -2986,12 +3064,19 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
 
     // Count
     {
+        const labelList& meshPoints = pp.meshPoints();
+
+        label nMasterPoints = 0;
         label nPlanar = 0;
         label nEdge = 0;
         label nPoint = 0;
 
         forAll(patchConstraints, pointI)
         {
+            if (isMasterPoint[meshPoints[pointI]])
+            {
+                nMasterPoints++;
+
                 if (patchConstraints[pointI].first() == 1)
                 {
                     nPlanar++;
@@ -3005,18 +3090,19 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
                     nPoint++;
                 }
             }
+        }
 
-        label nTotPoints = returnReduce(pp.nPoints(), sumOp<label>());
+        reduce(nMasterPoints, sumOp<label>());
         reduce(nPlanar, sumOp<label>());
         reduce(nEdge, sumOp<label>());
         reduce(nPoint, sumOp<label>());
-        Info<< "Feature analysis : total points:"
-            << nTotPoints
+        Info<< "Feature analysis : total master points:"
+            << nMasterPoints
             << " attraction to :" << nl
             << "    feature point   : " << nPoint << nl
             << "    feature edge    : " << nEdge << nl
             << "    nearest surface : " << nPlanar << nl
-            << "    rest            : " << nTotPoints-nPoint-nEdge-nPlanar
+            << "    rest            : " << nMasterPoints-nPoint-nEdge-nPlanar
             << nl
             << endl;
     }
@@ -3035,11 +3121,29 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
 
     if (featureAttract < 1-0.001)
     {
+        const PackedBoolList isMasterEdge(syncTools::getMasterEdges(mesh));
+
+        const vectorField pointNormals
+        (
+            PatchTools::pointNormals
+            (
+                mesh,
+                pp
+            )
+        );
+        const labelList meshEdges
+        (
+            pp.meshEdges(mesh.edges(), mesh.pointEdges())
+        );
+
+
         // 1. Smoothed all displacement
         vectorField smoothedPatchDisp = patchDisp;
         smoothAndConstrain
         (
+            isMasterEdge,
             pp,
+            meshEdges,
             patchConstraints,
             smoothedPatchDisp
         );
@@ -3047,16 +3151,18 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
 
         // 2. Smoothed tangential component
         vectorField tangPatchDisp = patchDisp;
-        tangPatchDisp -= (pp.pointNormals() & patchDisp) * pp.pointNormals();
+        tangPatchDisp -= (pointNormals & patchDisp) * pointNormals;
         smoothAndConstrain
         (
+            isMasterEdge,
             pp,
+            meshEdges,
             patchConstraints,
             tangPatchDisp
         );
 
         // Re-add normal component
-        tangPatchDisp += (pp.pointNormals() & patchDisp) * pp.pointNormals();
+        tangPatchDisp += (pointNormals & patchDisp) * pointNormals;
 
         if (debug&meshRefinement::OBJINTERSECTIONS)
         {