From 6adba987a2cff3f54f8d5d9eac299ab3a1214876 Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Fri, 23 Aug 2013 13:38:33 +0100
Subject: [PATCH 01/20] ENH: Updated cloud for topology change in parallel

---
 .../clouds/Templates/KinematicCloud/KinematicCloud.C        | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
index acf9e8c235..7cab777680 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C
@@ -652,7 +652,11 @@ void Foam::KinematicCloud<CloudType>::scaleSources()
 template<class CloudType>
 void Foam::KinematicCloud<CloudType>::preEvolve()
 {
-    Info<< "\nSolving cloud " << this->name() << endl;
+    // force calculaion of mesh dimensions - needed for parallel runs
+    // with topology change due to lazy evaluation of valid mesh dimensions
+    label nGeometricD = mesh_.nGeometricD();
+
+    Info<< "\nSolving " << nGeometricD << "-D cloud " << this->name() << endl;
 
     this->dispersion().cacheFields(true);
     forces_.cacheFields(true);

From 7c3d2ec15eaf88998cbf61ac014617029f8f0018 Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Fri, 23 Aug 2013 14:19:32 +0100
Subject: [PATCH 02/20] ENH: Updated R utility to include compressible cases

---
 .../postProcessing/turbulence/R/Make/options  |  18 ++-
 .../utilities/postProcessing/turbulence/R/R.C | 141 +++++++++++++++---
 .../turbulence/R/createFields.H               |  22 ---
 3 files changed, 135 insertions(+), 46 deletions(-)
 delete mode 100644 applications/utilities/postProcessing/turbulence/R/createFields.H

diff --git a/applications/utilities/postProcessing/turbulence/R/Make/options b/applications/utilities/postProcessing/turbulence/R/Make/options
index a136b16661..27b70cae0a 100644
--- a/applications/utilities/postProcessing/turbulence/R/Make/options
+++ b/applications/utilities/postProcessing/turbulence/R/Make/options
@@ -1,13 +1,19 @@
 EXE_INC = \
-    -I$(LIB_SRC)/postProcessing/postCalc \
-    -I$(LIB_SRC)/turbulenceModels \
-    -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \
+    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/turbulenceModels \
+    -I$(LIB_SRC)/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
 EXE_LIBS = \
-    $(FOAM_LIBBIN)/postCalc.o \
-    -lincompressibleRASModels \
     -lincompressibleTransportModels \
+    -lincompressibleRASModels \
+    -lfluidThermophysicalModels \
+    -lspecie \
+    -lcompressibleRASModels \
     -lfiniteVolume \
-    -lgenericPatchFields
+    -lgenericPatchFields \
+    -lmeshTools \
+    -lsampling
diff --git a/applications/utilities/postProcessing/turbulence/R/R.C b/applications/utilities/postProcessing/turbulence/R/R.C
index b96bb03947..708cda7214 100644
--- a/applications/utilities/postProcessing/turbulence/R/R.C
+++ b/applications/utilities/postProcessing/turbulence/R/R.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
@@ -29,35 +29,140 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
-#include "calc.H"
 #include "fvCFD.H"
 #include "incompressible/singlePhaseTransportModel/singlePhaseTransportModel.H"
-#include "RASModel.H"
+#include "incompressible/turbulenceModel/turbulenceModel.H"
+
+#include "fluidThermo.H"
+#include "compressible/turbulenceModel/turbulenceModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
+void calcIncompressibleR
+(
+    const fvMesh& mesh,
+    const Time& runTime,
+    const volVectorField& U
+)
 {
-    #include "createFields.H"
+    #include "createPhi.H"
 
-    Info<< "\nCalculating the Reynolds Stress R\n" << endl;
+    singlePhaseTransportModel laminarTransport(U, phi);
 
-    volSymmTensorField R
+    autoPtr<incompressible::turbulenceModel> model
     (
-        IOobject
-        (
-            "R",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        RASModel->R()
+        incompressible::turbulenceModel::New(U, phi, laminarTransport)
     );
 
-    R.write();
+    Info<< "Writing R field" << nl << endl;
 
-    Info<< "End" << endl;
+    model->R()().write();
+}
+
+
+void calcCompressibleR
+(
+    const fvMesh& mesh,
+    const Time& runTime,
+    const volVectorField& U
+)
+{
+    IOobject rhoHeader
+    (
+        "rho",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    );
+
+    if (!rhoHeader.headerOk())
+    {
+        Info<< "    no " << rhoHeader.name() <<" field" << endl;
+        return;
+    }
+
+    Info<< "Reading field rho\n" << endl;
+    volScalarField rho(rhoHeader, mesh);
+
+    #include "compressibleCreatePhi.H"
+
+    autoPtr<fluidThermo> pThermo(fluidThermo::New(mesh));
+    fluidThermo& thermo = pThermo();
+
+    autoPtr<compressible::turbulenceModel> model
+    (
+        compressible::turbulenceModel::New
+        (
+            rho,
+            U,
+            phi,
+            thermo
+        )
+    );
+
+    Info<< "Writing R field" << nl << endl;
+
+    model->R()().write();
+}
+
+
+int main(int argc, char *argv[])
+{
+    timeSelector::addOptions();
+
+    #include "addRegionOption.H"
+
+    argList::addBoolOption
+    (
+        "compressible",
+        "calculate compressible R"
+    );
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+    instantList timeDirs = timeSelector::select0(runTime, args);
+    #include "createNamedMesh.H"
+
+    const bool compressible = args.optionFound("compressible");
+
+    forAll(timeDirs, timeI)
+    {
+        runTime.setTime(timeDirs[timeI], timeI);
+        Info<< "Time = " << runTime.timeName() << endl;
+
+        IOobject UHeader
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        );
+
+        if (UHeader.headerOk())
+        {
+            Info<< "Reading field " << UHeader.name() << nl << endl;
+            volVectorField U(UHeader, mesh);
+
+            if (compressible)
+            {
+                calcCompressibleR(mesh, runTime, U);
+            }
+            else
+            {
+                calcIncompressibleR(mesh, runTime, U);
+            }
+        }
+        else
+        {
+            Info<< "    no " << UHeader.name() << " field" << endl;
+        }
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
 }
 
 
diff --git a/applications/utilities/postProcessing/turbulence/R/createFields.H b/applications/utilities/postProcessing/turbulence/R/createFields.H
deleted file mode 100644
index 7c07f44f8b..0000000000
--- a/applications/utilities/postProcessing/turbulence/R/createFields.H
+++ /dev/null
@@ -1,22 +0,0 @@
-    Info<< "Reading field U\n" << endl;
-    volVectorField U
-    (
-        IOobject
-        (
-            "U",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    #include "createPhi.H"
-
-    singlePhaseTransportModel laminarTransport(U, phi);
-
-    autoPtr<incompressible::RASModel> RASModel
-    (
-        incompressible::RASModel::New(U, phi, laminarTransport)
-    );

From 2f04e9754b683472baac761ace02d111a4f69eef Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Tue, 27 Aug 2013 11:22:27 +0100
Subject: [PATCH 03/20] ENH: Updated error messages

---
 .../StandardWallInteraction/StandardWallInteraction.C | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/StandardWallInteraction/StandardWallInteraction.C b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/StandardWallInteraction/StandardWallInteraction.C
index d9b38daccf..fa86f67ca9 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/StandardWallInteraction/StandardWallInteraction.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/StandardWallInteraction/StandardWallInteraction.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
@@ -57,7 +57,7 @@ Foam::StandardWallInteraction<CloudType>::StandardWallInteraction
                 "StandardWallInteraction<CloudType>::StandardWallInteraction"
                 "("
                     "const dictionary&, "
-                    "CloudType& cloud"
+                    "CloudType&"
                 ")"
             )   << "Unknown interaction result type "
                 << interactionTypeName
@@ -174,10 +174,11 @@ bool Foam::StandardWallInteraction<CloudType>::correct
                 (
                     "bool StandardWallInteraction<CloudType>::correct"
                     "("
+                        "typename CloudType::parcelType&, "
                         "const polyPatch&, "
-                        "const label, "
-                        "bool&, "
-                        "vector&"
+                        "bool& keepParticle, "
+                        "const scalar, "
+                        "const tetIndices&"
                     ") const"
                 )   << "Unknown interaction type "
                     << this->interactionTypeToWord(interactionType_)

From 6a1b56732771973dfd1ee7ef7b8551959c59eeb7 Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Tue, 27 Aug 2013 18:11:13 +0100
Subject: [PATCH 04/20] ENH: Cloud tracking - handle case where tracking stalls
 due to stepFraction

---
 src/lagrangian/basic/particle/particle.C      |  4 ++-
 src/lagrangian/basic/particle/particle.H      |  3 ++
 .../KinematicParcel/KinematicParcel.C         | 30 ++++++++++++++-----
 3 files changed, 29 insertions(+), 8 deletions(-)

diff --git a/src/lagrangian/basic/particle/particle.C b/src/lagrangian/basic/particle/particle.C
index 8ffb988d45..c71559d042 100644
--- a/src/lagrangian/basic/particle/particle.C
+++ b/src/lagrangian/basic/particle/particle.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
@@ -34,6 +34,8 @@ const Foam::scalar Foam::particle::trackingCorrectionTol = 1e-5;
 
 const Foam::scalar Foam::particle::lambdaDistanceToleranceCoeff = 1e3*SMALL;
 
+const Foam::scalar Foam::particle::minStepFractionTol = 1e5*SMALL;
+
 namespace Foam
 {
     defineTypeNameAndDebug(particle, 0);
diff --git a/src/lagrangian/basic/particle/particle.H b/src/lagrangian/basic/particle/particle.H
index 15742cacc9..0275cc5fd1 100644
--- a/src/lagrangian/basic/particle/particle.H
+++ b/src/lagrangian/basic/particle/particle.H
@@ -307,6 +307,9 @@ public:
         //  for the denominator and numerator of lambda
         static const scalar lambdaDistanceToleranceCoeff;
 
+        //- Minimum stepFraction tolerance
+        static const scalar minStepFractionTol;
+
 
     // Constructors
 
diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
index 8cf9317c86..10bd38d62e 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
@@ -271,6 +271,8 @@ bool Foam::KinematicParcel<ParcelType>::move
     scalar tEnd = (1.0 - p.stepFraction())*trackTime;
     const scalar dtMax = tEnd;
 
+    bool moving = true;
+
     while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL)
     {
         // Apply correction to position for reduced-D cases
@@ -281,22 +283,36 @@ bool Foam::KinematicParcel<ParcelType>::move
         // Set the Lagrangian time-step
         scalar dt = min(dtMax, tEnd);
 
-        // Remember which cell the parcel is in since this will change if
-        // a face is hit
+        // Cache the parcel current cell as this will change if a face is hit
         const label cellI = p.cell();
 
         const scalar magU = mag(U_);
-        if (p.active() && magU > ROOTVSMALL)
+        if (p.active())
         {
             const scalar d = dt*magU;
             const scalar dCorr = min(d, maxCo*cbrt(V[cellI]));
-            dt *=
-                dCorr/d
-               *p.trackToFace(p.position() + dCorr*U_/magU, td);
+            dt *= dCorr/d;
+
+            if (moving && (magU > ROOTVSMALL))
+            {
+                dt *= p.trackToFace(p.position() + dCorr*U_/magU, td);
+            }
         }
 
         tEnd -= dt;
-        p.stepFraction() = 1.0 - tEnd/trackTime;
+
+        scalar newStepFraction = 1.0 - tEnd/trackTime;
+
+        if
+        (
+            mag(p.stepFraction() - newStepFraction)
+          < particle::minStepFractionTol
+        )
+        {
+            moving = false;
+        }
+
+        p.stepFraction() = newStepFraction;
 
         // Avoid problems with extremely small timesteps
         if (dt > ROOTVSMALL)

From d62efed988908822280370de0fb4b68a58ed6a31 Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Tue, 27 Aug 2013 18:11:35 +0100
Subject: [PATCH 05/20] ENH: CollidingParcel - refactored to remove code
 duplication

---
 .../CollidingParcel/CollidingParcel.C         | 72 +------------------
 1 file changed, 1 insertion(+), 71 deletions(-)

diff --git a/src/lagrangian/intermediate/parcels/Templates/CollidingParcel/CollidingParcel.C b/src/lagrangian/intermediate/parcels/Templates/CollidingParcel/CollidingParcel.C
index f04ac77a15..f422f7ea37 100644
--- a/src/lagrangian/intermediate/parcels/Templates/CollidingParcel/CollidingParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/CollidingParcel/CollidingParcel.C
@@ -69,13 +69,6 @@ bool Foam::CollidingParcel<ParcelType>::move
     typename TrackData::cloudType::parcelType& p =
         static_cast<typename TrackData::cloudType::parcelType&>(*this);
 
-    td.switchProcessor = false;
-    td.keepParticle = true;
-
-    const polyMesh& mesh = td.cloud().pMesh();
-    const polyBoundaryMesh& pbMesh = mesh.boundaryMesh();
-    const scalarField& V = mesh.cellVolumes();
-
     switch (td.part())
     {
         case TrackData::tpVelocityHalfStep:
@@ -92,70 +85,7 @@ bool Foam::CollidingParcel<ParcelType>::move
 
         case TrackData::tpLinearTrack:
         {
-            scalar tEnd = (1.0 - p.stepFraction())*trackTime;
-            const scalar dtMax = tEnd;
-
-            while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL)
-            {
-                // Apply correction to position for reduced-D cases
-                meshTools::constrainToMeshCentre(mesh, p.position());
-
-                const point start(p.position());
-
-                // Set the Lagrangian time-step
-                scalar dt = min(dtMax, tEnd);
-
-                // Remember which cell the parcel is in since this
-                // will change if a face is hit
-                const label cellI = p.cell();
-
-                const scalar magU = mag(p.U());
-                if (p.active() && magU > ROOTVSMALL)
-                {
-                    const scalar d = dt*magU;
-                    const scalar maxCo = td.cloud().solution().maxCo();
-                    const scalar dCorr = min(d, maxCo*cbrt(V[cellI]));
-                    dt *=
-                        dCorr/d
-                       *p.trackToFace(p.position() + dCorr*p.U()/magU, td);
-                }
-
-                tEnd -= dt;
-                p.stepFraction() = 1.0 - tEnd/trackTime;
-
-                // Avoid problems with extremely small timesteps
-                if (dt > ROOTVSMALL)
-                {
-                    // Update cell based properties
-                    p.setCellValues(td, dt, cellI);
-
-                    if (td.cloud().solution().cellValueSourceCorrection())
-                    {
-                        p.cellValueSourceCorrection(td, dt, cellI);
-                    }
-
-                    p.calc(td, dt, cellI);
-                }
-
-                if (p.onBoundary() && td.keepParticle)
-                {
-                    if (isA<processorPolyPatch>(pbMesh[p.patch(p.face())]))
-                    {
-                        td.switchProcessor = true;
-                    }
-                }
-
-                p.age() += dt;
-
-                td.cloud().functions().postMove
-                (
-                    p,
-                    cellI,
-                    dt,
-                    start,
-                    td.keepParticle
-                );
-            }
+            ParcelType::move(td, trackTime);
 
             break;
         }

From 14a309d4c29f01fcc95c778d7531a8e5e550bd37 Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Wed, 28 Aug 2013 09:16:31 +0100
Subject: [PATCH 06/20] ENH: cloudSolution - maxCo available for both steady
 state and transient

---
 .../Templates/KinematicCloud/cloudSolution/cloudSolution.C    | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.C b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.C
index fd7bf31706..c2bea38236 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.C
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.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
@@ -110,11 +110,11 @@ void Foam::cloudSolution::read()
     dict_.lookup("transient") >> transient_;
     dict_.lookup("coupled") >> coupled_;
     dict_.lookup("cellValueSourceCorrection") >> cellValueSourceCorrection_;
+    dict_.lookup("maxCo") >> maxCo_;
 
     if (steadyState())
     {
         dict_.lookup("calcFrequency") >> calcFrequency_;
-        dict_.lookup("maxCo") >> maxCo_;
         dict_.lookup("maxTrackTime") >> maxTrackTime_;
 
         if (coupled_)

From df2273f0c46b8f4fe0279fcd4888c079ea0d1842 Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Wed, 28 Aug 2013 10:28:10 +0100
Subject: [PATCH 07/20] ENH: Refactored code

---
 .../reactingParcelFilmFoam/reactingParcelFilmFoam.C       | 8 ++------
 1 file changed, 2 insertions(+), 6 deletions(-)

diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C b/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C
index 65cfcb4516..4067c125eb 100644
--- a/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C
+++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C
@@ -104,13 +104,9 @@ int main(int argc, char *argv[])
             }
 
             rho = thermo.rho();
+        }
 
-            runTime.write();
-        }
-        else
-        {
-            runTime.write();
-        }
+        runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"

From d7f67c78b1a35b15680e2a6b218ebf6cbf9c7ade Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Wed, 28 Aug 2013 10:44:09 +0100
Subject: [PATCH 08/20] ENH: Tutorial update

---
 .../les/oppositeBurningPanels/constant/reactingCloud1Properties  | 1 +
 .../les/smallPoolFire2D/constant/reactingCloud1Properties        | 1 +
 .../les/smallPoolFire3D/constant/reactingCloud1Properties        | 1 +
 .../simplifiedSiwek/constant/coalCloud1Properties                | 1 +
 .../simplifiedSiwek/constant/limestoneCloud1Properties           | 1 +
 .../hopper/hopperEmptying/constant/kinematicCloudProperties      | 1 +
 .../hopper/hopperInitialState/constant/kinematicCloudProperties  | 1 +
 .../cylinder/constant/reactingCloud1Properties                   | 1 +
 .../hotBoxes/constant/reactingCloud1Properties                   | 1 +
 .../splashPanel/constant/reactingCloud1Properties                | 1 +
 .../reactingParcelFoam/filter/constant/reactingCloud1Properties  | 1 +
 .../parcelInBox/constant/reactingCloud1Properties                | 1 +
 .../verticalChannel/constant/reactingCloud1Properties            | 1 +
 .../sprayFoam/aachenBomb/constant/sprayCloudProperties           | 1 +
 14 files changed, 14 insertions(+)

diff --git a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/reactingCloud1Properties b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/reactingCloud1Properties
index 5408eb7928..9228607e91 100644
--- a/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/reactingCloud1Properties
+++ b/tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/reactingCloud1Properties
@@ -21,6 +21,7 @@ solution
     coupled         yes;
     transient       yes;
     cellValueSourceCorrection yes;
+    maxCo           0.3;
 
     sourceTerms
     {
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/reactingCloud1Properties b/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/reactingCloud1Properties
index 5408eb7928..9228607e91 100644
--- a/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/reactingCloud1Properties
+++ b/tutorials/combustion/fireFoam/les/smallPoolFire2D/constant/reactingCloud1Properties
@@ -21,6 +21,7 @@ solution
     coupled         yes;
     transient       yes;
     cellValueSourceCorrection yes;
+    maxCo           0.3;
 
     sourceTerms
     {
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/reactingCloud1Properties b/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/reactingCloud1Properties
index 5408eb7928..9228607e91 100644
--- a/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/reactingCloud1Properties
+++ b/tutorials/combustion/fireFoam/les/smallPoolFire3D/constant/reactingCloud1Properties
@@ -21,6 +21,7 @@ solution
     coupled         yes;
     transient       yes;
     cellValueSourceCorrection yes;
+    maxCo           0.3;
 
     sourceTerms
     {
diff --git a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/coalCloud1Properties b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/coalCloud1Properties
index 7c71c8971e..e164395342 100644
--- a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/coalCloud1Properties
+++ b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/coalCloud1Properties
@@ -21,6 +21,7 @@ solution
     transient       yes;
     coupled         true;
     cellValueSourceCorrection on;
+    maxCo           0.3;
 
     sourceTerms
     {
diff --git a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/limestoneCloud1Properties b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/limestoneCloud1Properties
index 4b5b5393e4..36b287137d 100644
--- a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/limestoneCloud1Properties
+++ b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/limestoneCloud1Properties
@@ -21,6 +21,7 @@ solution
     coupled         true;
     transient       yes;
     cellValueSourceCorrection on;
+    maxCo           0.3;
 
     sourceTerms
     {
diff --git a/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperEmptying/constant/kinematicCloudProperties b/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperEmptying/constant/kinematicCloudProperties
index 06571ffdd1..150adb2953 100644
--- a/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperEmptying/constant/kinematicCloudProperties
+++ b/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperEmptying/constant/kinematicCloudProperties
@@ -21,6 +21,7 @@ solution
     coupled         false;
     transient       yes;
     cellValueSourceCorrection off;
+    maxCo           0.3;
 
     sourceTerms
     {
diff --git a/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperInitialState/constant/kinematicCloudProperties b/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperInitialState/constant/kinematicCloudProperties
index 874c386653..2fa48b6da0 100644
--- a/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperInitialState/constant/kinematicCloudProperties
+++ b/tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperInitialState/constant/kinematicCloudProperties
@@ -21,6 +21,7 @@ solution
     coupled         false;
     transient       yes;
     cellValueSourceCorrection off;
+    maxCo           0.3;
 
     interpolationSchemes
     {
diff --git a/tutorials/lagrangian/reactingParcelFilmFoam/cylinder/constant/reactingCloud1Properties b/tutorials/lagrangian/reactingParcelFilmFoam/cylinder/constant/reactingCloud1Properties
index c96b27a4b4..6e01c75332 100644
--- a/tutorials/lagrangian/reactingParcelFilmFoam/cylinder/constant/reactingCloud1Properties
+++ b/tutorials/lagrangian/reactingParcelFilmFoam/cylinder/constant/reactingCloud1Properties
@@ -21,6 +21,7 @@ solution
     coupled         no;
     transient       yes;
     cellValueSourceCorrection no;
+    maxCo           0.3;
 
     sourceTerms
     {
diff --git a/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/constant/reactingCloud1Properties b/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/constant/reactingCloud1Properties
index e6afd3b55d..4e8221412f 100644
--- a/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/constant/reactingCloud1Properties
+++ b/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/constant/reactingCloud1Properties
@@ -21,6 +21,7 @@ solution
     coupled         yes;
     transient       yes;
     cellValueSourceCorrection yes;
+    maxCo           0.3;
 
     sourceTerms
     {
diff --git a/tutorials/lagrangian/reactingParcelFilmFoam/splashPanel/constant/reactingCloud1Properties b/tutorials/lagrangian/reactingParcelFilmFoam/splashPanel/constant/reactingCloud1Properties
index 5652c17be9..e38561ee8f 100644
--- a/tutorials/lagrangian/reactingParcelFilmFoam/splashPanel/constant/reactingCloud1Properties
+++ b/tutorials/lagrangian/reactingParcelFilmFoam/splashPanel/constant/reactingCloud1Properties
@@ -21,6 +21,7 @@ solution
     coupled         no;
     transient       yes;
     cellValueSourceCorrection no;
+    maxCo           0.3;
 
     sourceTerms
     {
diff --git a/tutorials/lagrangian/reactingParcelFoam/filter/constant/reactingCloud1Properties b/tutorials/lagrangian/reactingParcelFoam/filter/constant/reactingCloud1Properties
index 2ba702faa8..2eba294927 100644
--- a/tutorials/lagrangian/reactingParcelFoam/filter/constant/reactingCloud1Properties
+++ b/tutorials/lagrangian/reactingParcelFoam/filter/constant/reactingCloud1Properties
@@ -21,6 +21,7 @@ solution
     coupled         true;
     transient       yes;
     cellValueSourceCorrection on;
+    maxCo           0.3;
 
     sourceTerms
     {
diff --git a/tutorials/lagrangian/reactingParcelFoam/parcelInBox/constant/reactingCloud1Properties b/tutorials/lagrangian/reactingParcelFoam/parcelInBox/constant/reactingCloud1Properties
index 41fe1084d4..1ce6679e09 100644
--- a/tutorials/lagrangian/reactingParcelFoam/parcelInBox/constant/reactingCloud1Properties
+++ b/tutorials/lagrangian/reactingParcelFoam/parcelInBox/constant/reactingCloud1Properties
@@ -21,6 +21,7 @@ solution
     coupled         false;
     transient       yes;
     cellValueSourceCorrection off;
+    maxCo           0.3;
 
     sourceTerms
     {
diff --git a/tutorials/lagrangian/reactingParcelFoam/verticalChannel/constant/reactingCloud1Properties b/tutorials/lagrangian/reactingParcelFoam/verticalChannel/constant/reactingCloud1Properties
index 70524de8a8..c6e386958c 100644
--- a/tutorials/lagrangian/reactingParcelFoam/verticalChannel/constant/reactingCloud1Properties
+++ b/tutorials/lagrangian/reactingParcelFoam/verticalChannel/constant/reactingCloud1Properties
@@ -21,6 +21,7 @@ solution
     coupled         true;
     transient       yes;
     cellValueSourceCorrection on;
+    maxCo           0.3;
 
     sourceTerms
     {
diff --git a/tutorials/lagrangian/sprayFoam/aachenBomb/constant/sprayCloudProperties b/tutorials/lagrangian/sprayFoam/aachenBomb/constant/sprayCloudProperties
index b0989cc8b9..4d9cfe2711 100644
--- a/tutorials/lagrangian/sprayFoam/aachenBomb/constant/sprayCloudProperties
+++ b/tutorials/lagrangian/sprayFoam/aachenBomb/constant/sprayCloudProperties
@@ -21,6 +21,7 @@ solution
     coupled         true;
     transient       yes;
     cellValueSourceCorrection on;
+    maxCo           0.3;
 
     sourceTerms
     {

From 668cd823b5c302ab70b6fe3afc9f86580cb812d2 Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Thu, 29 Aug 2013 17:24:30 +0100
Subject: [PATCH 09/20] BUG: corrected steadyParticleTracks utility

---
 .../steadyParticleTracks.C                    | 42 ++++++++----
 .../steadyParticleTracksTemplates.C           | 66 ++++++++++++-------
 .../steadyParticleTracksTemplates.H           |  4 +-
 3 files changed, 73 insertions(+), 39 deletions(-)

diff --git a/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracks.C b/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracks.C
index 5a6fad75a6..4d2ab943ac 100644
--- a/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracks.C
+++ b/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracks.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
@@ -184,6 +184,7 @@ int main(int argc, char *argv[])
             }
         }
 
+
         if (nTracks == 0)
         {
             Info<< "\n    No track data" << endl;
@@ -202,11 +203,13 @@ int main(int argc, char *argv[])
 
             // particle "age" property used to sort the tracks
             List<SortableList<scalar> > agePerTrack(nTracks);
+            List<List<label> > particleMap(nTracks);
 
             forAll(trackLengths, i)
             {
                 const label length = trackLengths[i];
                 agePerTrack[i].setSize(length);
+                particleMap[i].setSize(length);
             }
 
             // store the particle age per track
@@ -228,6 +231,7 @@ int main(int argc, char *argv[])
                     const label trackI = particleToTrack[i];
                     const label sampleI = trackSamples[trackI];
                     agePerTrack[trackI][sampleI] = age[i];
+                    particleMap[trackI][sampleI] = i;
                     trackSamples[trackI]++;
                 }
                 tage.clear();
@@ -251,21 +255,30 @@ int main(int argc, char *argv[])
                 Info<< "\n    Writing points" << endl;
 
                 {
-                    label offset = 0;
                     forAll(agePerTrack, i)
                     {
                         agePerTrack[i].sort();
+
                         const labelList& ids = agePerTrack[i].indices();
+                        labelList& particleIds = particleMap[i];
+
+                        {
+                            // update addressing
+                            List<label> sortedIds(ids);
+                            forAll(sortedIds, j)
+                            {
+                                sortedIds[j] = particleIds[ids[j]];
+                            }
+                            particleIds = sortedIds;
+                        }
 
                         forAll(ids, j)
                         {
-                            const label localId = offset + ids[j];
+                            const label localId = particleIds[j];
                             const vector& pos = particles[localId].position();
                             os  << pos.x() << ' ' << pos.y() << ' ' << pos.z()
                                 << nl;
                         }
-
-                        offset += trackLengths[i];
                     }
                 }
 
@@ -278,13 +291,14 @@ int main(int argc, char *argv[])
                 // Write ids of track points to file
                 {
                     label globalPtI = 0;
-                    forAll(agePerTrack, i)
+                    forAll(particleMap, i)
                     {
-                        os  << agePerTrack[i].size() << nl;
+                        os  << particleMap[i].size() << nl;
 
-                        forAll(agePerTrack[i], j)
+                        forAll(particleMap[i], j)
                         {
                             os  << ' ' << globalPtI++;
+
                             if (((j + 1) % 10 == 0) && (j != 0))
                             {
                                 os  << nl;
@@ -303,14 +317,14 @@ int main(int argc, char *argv[])
 
                 Info<< "\n    Processing fields" << nl << endl;
 
-                processFields<label>(os, agePerTrack, userFields, cloudObjs);
-                processFields<scalar>(os, agePerTrack, userFields, cloudObjs);
-                processFields<vector>(os, agePerTrack, userFields, cloudObjs);
+                processFields<label>(os, particleMap, userFields, cloudObjs);
+                processFields<scalar>(os, particleMap, userFields, cloudObjs);
+                processFields<vector>(os, particleMap, userFields, cloudObjs);
                 processFields<sphericalTensor>
-                    (os, agePerTrack, userFields, cloudObjs);
+                    (os, particleMap, userFields, cloudObjs);
                 processFields<symmTensor>
-                    (os, agePerTrack, userFields, cloudObjs);
-                processFields<tensor>(os, agePerTrack, userFields, cloudObjs);
+                    (os, particleMap, userFields, cloudObjs);
+                processFields<tensor>(os, particleMap, userFields, cloudObjs);
 
             }
         }
diff --git a/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracksTemplates.C b/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracksTemplates.C
index 89e499b17d..3dcc365c1d 100644
--- a/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracksTemplates.C
+++ b/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracksTemplates.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
@@ -57,39 +57,57 @@ tmp<Field<Type> > readParticleField
         return tmp<Field<Type> >(new Field<Type>(newField.xfer()));
     }
 
-    Info<< "error: cloud field name " << name << " not found" << endl;
+    FatalErrorIn
+    (
+        "template<class Type>"
+        "void readParticleField"
+        "("
+            "const word&, "
+            "const IOobjectList"
+        ")"
+    )
+        << "error: cloud field name " << name << " not found"
+        << abort(FatalError);
 
     return Field<Type>::null();
 }
 
 
 template<class Type>
-PtrList<List<Type> > readFields
+void readFields
 (
     PtrList<List<Type> >& values,
-    const List<word>& fields,
+    const List<word>& fieldNames,
     const IOobjectList& cloudObjs
 )
 {
     IOobjectList objects(cloudObjs.lookupClass(IOField<Type>::typeName));
 
-    label fieldI = 0;
-    forAllConstIter(IOobjectList, objects, iter)
+    forAll(fieldNames, j)
     {
-        const IOobject& obj = *iter();
-        forAll(fields, j)
+        const IOobject* obj = objects.lookup(fieldNames[j]);
+        if (obj != NULL)
         {
-            if (obj.name() == fields[j])
-            {
-                Info<< "        reading field " << obj.name() << endl;
-                IOField<Type> newField(obj);
-                values.set(fieldI++, new List<Type>(newField.xfer()));
-                break;
-            }
+            Info<< "        reading field " << fieldNames[j] << endl;
+            IOField<Type> newField(*obj);
+            values.set(j, new List<Type>(newField.xfer()));
+        }
+        else
+        {
+            FatalErrorIn
+            (
+                "template<class Type>"
+                "void readFields"
+                "("
+                    "PtrList<List<Type> >&, "
+                    "const List<word>&, "
+                    "const IOobjectList&"
+                ")"
+            )
+                << "Unable to read field " << fieldNames[j]
+                << abort(FatalError);
         }
     }
-
-    return values;
 }
 
 
@@ -109,7 +127,7 @@ void writeVTKFields
 (
     OFstream& os,
     const PtrList<List<Type> >& values,
-    const List<SortableList<scalar> >& agePerTrack,
+    const List<List<label> >& addr,
     const List<word>& fieldNames
 )
 {
@@ -121,9 +139,9 @@ void writeVTKFields
         os  << nl << fieldNames[fieldI] << ' ' << pTraits<Type>::nComponents
             << ' ' << values[fieldI].size() << " float" << nl;
         label offset = 0;
-        forAll(agePerTrack, trackI)
+        forAll(addr, trackI)
         {
-            const List<label> ids = agePerTrack[trackI].indices() + offset;
+            const List<label> ids(addr[trackI]);
 
             List<Type> data(UIndirectList<Type>(values[fieldI], ids));
             label nData = data.size() - 1;
@@ -149,7 +167,7 @@ template<class Type>
 void processFields
 (
     OFstream& os,
-    const List<SortableList<scalar> >& agePerTrack,
+    const List<List<label> >& addr,
     const List<word>& userFieldNames,
     const IOobjectList& cloudObjs
 )
@@ -176,12 +194,14 @@ void processFields
         (
             os,
             values,
-            agePerTrack,
-            fieldNames.xfer()
+            addr,
+            fieldNames
         );
     }
 }
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 } // End namespace Foam
 
 // ************************************************************************* //
diff --git a/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracksTemplates.H b/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracksTemplates.H
index 973e6da255..f3b5d48d2d 100644
--- a/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracksTemplates.H
+++ b/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracksTemplates.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
@@ -49,7 +49,7 @@ namespace Foam
     );
 
     template<class Type>
-    PtrList<List<Type> > readFields
+    void readFields
     (
         PtrList<List<Type> >& values,
         const List<word>& fields,

From b7187c6018c0039e232524708b08806c6b2a36dc Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Mon, 2 Sep 2013 22:50:22 +0100
Subject: [PATCH 10/20] Remove unnecessary library

---
 .../multiphase/multiphaseEulerFoam/multiphaseSystem/Make/options | 1 -
 1 file changed, 1 deletion(-)

diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/Make/options b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/Make/options
index ce7e7a2d85..68bcd9ef06 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/Make/options
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/Make/options
@@ -10,5 +10,4 @@ EXE_INC = \
 LIB_LIBS = \
     -linterfaceProperties \
     -lincompressibleTransportModels \
-    -lcompressibleMultiphaseEulerianInterfacialModels \
     -lfiniteVolume

From 3e12a89a041f728471846cbc7828bbd4d837ea8b Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Mon, 2 Sep 2013 22:51:15 +0100
Subject: [PATCH 11/20] Check mesh is 2D and scale resulting steam-function by
 the thickness

---
 .../streamFunction/streamFunction.C           | 23 +++++++++++++++----
 1 file changed, 18 insertions(+), 5 deletions(-)

diff --git a/applications/utilities/postProcessing/velocityField/streamFunction/streamFunction.C b/applications/utilities/postProcessing/velocityField/streamFunction/streamFunction.C
index aa56086c90..de9af38b87 100644
--- a/applications/utilities/postProcessing/velocityField/streamFunction/streamFunction.C
+++ b/applications/utilities/postProcessing/velocityField/streamFunction/streamFunction.C
@@ -42,14 +42,26 @@ Description
 int main(int argc, char *argv[])
 {
     timeSelector::addOptions();
-#   include "addRegionOption.H"
+    #include "addRegionOption.H"
 
-#   include "setRootCase.H"
-#   include "createTime.H"
+    #include "setRootCase.H"
+    #include "createTime.H"
 
     instantList timeDirs = timeSelector::select0(runTime, args);
 
-#   include "createNamedMesh.H"
+    #include "createNamedMesh.H"
+
+    label nD = mesh.nGeometricD();
+
+    if (nD != 2)
+    {
+        FatalErrorIn(args.executable())
+            << "Case is not 2D, stream-function cannot be computed"
+            << exit(FatalError);
+    }
+
+    const vector slabDir((Vector<label>::one - mesh.geometricD())/2);
+    scalar thickness = slabDir & mesh.bounds().span();
 
     const pointMesh& pMesh = pointMesh::New(mesh);
 
@@ -428,7 +440,6 @@ int main(int argc, char *argv[])
                      }
 
                      Info<< ".";
-//                     Info<< "One pass, n visited = " << nVisited << endl;
 
                      if (nVisited == nVisitedOld)
                      {
@@ -449,6 +460,8 @@ int main(int argc, char *argv[])
                 Info<< endl;
             } while (!finished);
 
+            // Normalise the stream-function by the 2D mesh thickness
+            streamFunction /= thickness;
             streamFunction.boundaryField() = 0.0;
             streamFunction.write();
         }

From e9cdbdef5e47c6e5e6ed5e1338b2b1a9d07bd81c Mon Sep 17 00:00:00 2001
From: laurence <laurence>
Date: Fri, 30 Aug 2013 12:18:15 +0100
Subject: [PATCH 12/20] ENH: foamyHexMesh: Allow surface point insertion on
 close surfaces

ENH: foamyHexMesh: do not point motion to cross surfaces
ENH: foamyHexMesh: do not allow perturbation to cross surface
ENH: foamyHexMesh: change inside/outside test to check for nearest intersection
BUG: foamyHexMesh: With surfaces of same priority do not overwrite cell size
ENH: foamyHexMesh: insert surface points on close parallel surfaces
BUG: foamyHexMesh: write to own processor file in parallel
ENH: foamyHexMesh: Insert surface points if internal points are exposed
ENH: foamyHexMesh: add hasBoundaryPoint() to indexedCell
REVERT: foamyHexMesh: do not check inside/outside for surface points
REVERT: foamyHexMesh: correct function argument list
REVERT: foamyHexMesh: correct function argument list
---
 .../controlMeshRefinement.C                   |  16 +-
 .../conformalVoronoiMesh.C                    |   9 +-
 .../conformalVoronoiMeshConformToSurface.C    | 183 ++++++++++++------
 .../conformalVoronoiMeshI.H                   |  35 ++--
 .../indexedCell/indexedCell.H                 |   1 +
 .../indexedCell/indexedCellI.H                |  13 ++
 .../conformationSurfaces.C                    |  18 +-
 .../rayShooting/rayShooting.C                 | 120 +++++++++---
 .../rayShooting/rayShooting.H                 |  10 +
 9 files changed, 289 insertions(+), 116 deletions(-)

diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cellShapeControl/controlMeshRefinement/controlMeshRefinement.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cellShapeControl/controlMeshRefinement/controlMeshRefinement.C
index b100d8eebb..8d954799a1 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cellShapeControl/controlMeshRefinement/controlMeshRefinement.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cellShapeControl/controlMeshRefinement/controlMeshRefinement.C
@@ -295,14 +295,14 @@ void Foam::controlMeshRefinement::initialMeshPopulation
                     shapeController_.minimumCellSize()
                 );
             }
-            else if (maxPriority == controlFunction.maxPriority())
-            {
-                vertices[vI].targetCellSize() = max
-                (
-                    min(sizes[vI], size),
-                    shapeController_.minimumCellSize()
-                );
-            }
+//            else if (maxPriority == controlFunction.maxPriority())
+//            {
+//                vertices[vI].targetCellSize() = max
+//                (
+//                    min(sizes[vI], size),
+//                    shapeController_.minimumCellSize()
+//                );
+//            }
             else
             {
                 vertices[vI].targetCellSize() = max
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
index e580febdd8..8224492921 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
@@ -1487,8 +1487,13 @@ void Foam::conformalVoronoiMesh::move()
     // Save displacements to file.
     if (foamyHexMeshControls().objOutput() && time().outputTime())
     {
-        Pout<< "Writing point displacement vectors to file." << endl;
-        OFstream str("displacements_" + runTime_.timeName() + ".obj");
+        Info<< "Writing point displacement vectors to file." << endl;
+        OFstream str
+        (
+            time().path()
+          + "displacements_" + runTime_.timeName()
+          + ".obj"
+        );
 
         for
         (
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
index db77f62f5f..f5a023592c 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
@@ -550,53 +550,130 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation()
             }
         }
 
+//        for
+//        (
+//            Delaunay::Finite_edges_iterator eit = finite_edges_begin();
+//            eit != finite_edges_end();
+//            ++eit
+//        )
+//        {
+//            Cell_handle c = eit->first;
+//            Vertex_handle vA = c->vertex(eit->second);
+//            Vertex_handle vB = c->vertex(eit->third);
+//
+//            if
+//            (
+//                vA->referred()
+//             || vB->referred()
+//            )
+//            {
+//                continue;
+//            }
+//
+//            if
+//            (
+//                (vA->internalPoint() && vA->referred())
+//             || (vB->internalPoint() && vB->referred())
+//            )
+//            {
+//                continue;
+//            }
+//
+//            if
+//            (
+//                (vA->internalPoint() && vB->externalBoundaryPoint())
+//             || (vB->internalPoint() && vA->externalBoundaryPoint())
+//             || (vA->internalBoundaryPoint() && vB->internalBoundaryPoint())
+//            )
+//            {
+//                pointIndexHitAndFeatureDynList surfaceIntersections(0.5*AtoV);
+//                pointIndexHit surfHit;
+//                label hitSurface;
+//
+//                geometryToConformTo_.findSurfaceNearest
+//                (
+//                    (
+//                        vA->internalPoint()
+//                      ? topoint(vA->point())
+//                      : topoint(vB->point())
+//                    ),
+//                    magSqr(topoint(vA->point()) - topoint(vB->point())),
+//                    surfHit,
+//                    hitSurface
+//                );
+//
+//                if (surfHit.hit())
+//                {
+//                    surfaceIntersections.append
+//                    (
+//                        pointIndexHitAndFeature(surfHit, hitSurface)
+//                    );
+//
+//                    addSurfaceAndEdgeHits
+//                    (
+//                        (
+//                            vA->internalPoint()
+//                          ? topoint(vA->point())
+//                          : topoint(vB->point())
+//                        ),
+//                        surfaceIntersections,
+//                        surfacePtReplaceDistCoeffSqr,
+//                        edgeSearchDistCoeffSqr,
+//                        surfaceHits,
+//                        featureEdgeHits,
+//                        surfaceToTreeShape,
+//                        edgeToTreeShape,
+//                        surfacePtToEdgePtDist,
+//                        false
+//                    );
+//                }
+//            }
+//        }
+
         for
         (
-            Delaunay::Finite_edges_iterator eit = finite_edges_begin();
-            eit != finite_edges_end();
-            ++eit
+            Delaunay::Finite_cells_iterator cit = finite_cells_begin();
+            cit != finite_cells_end();
+            ++cit
         )
         {
-            Cell_handle c = eit->first;
-            Vertex_handle vA = c->vertex(eit->second);
-            Vertex_handle vB = c->vertex(eit->third);
-
-            if
-            (
-                vA->referred()
-             || vB->referred()
-            )
+            label nInternal = 0;
+            for (label vI = 0; vI < 4; vI++)
             {
-                continue;
+                if (cit->vertex(vI)->internalPoint())
+                {
+                    nInternal++;
+                }
             }
 
-            if
-            (
-                (vA->internalPoint() && vA->referred())
-             || (vB->internalPoint() && vB->referred())
-            )
+            if (nInternal == 1 && cit->hasBoundaryPoint())
+            //if (cit->boundaryDualVertex() && !cit->hasReferredPoint())
             {
-                continue;
-            }
+                const Foam::point& pt = cit->dual();
+
+                const Foam::point cellCentre =
+                    topoint
+                    (
+                        CGAL::centroid
+                        (
+                            CGAL::Tetrahedron_3<baseK>
+                            (
+                                cit->vertex(0)->point(),
+                                cit->vertex(1)->point(),
+                                cit->vertex(2)->point(),
+                                cit->vertex(3)->point()
+                            )
+                        )
+                    );
 
-            if
-            (
-                (vA->internalPoint() && vB->externalBoundaryPoint())
-             || (vB->internalPoint() && vA->externalBoundaryPoint())
-            )
-            {
                 pointIndexHitAndFeatureDynList surfaceIntersections(0.5*AtoV);
                 pointIndexHit surfHit;
                 label hitSurface;
 
-                geometryToConformTo_.findSurfaceNearest
+                geometryToConformTo_.findSurfaceNearestIntersection
                 (
-                    (
-                        vA->internalPoint()
-                      ? topoint(vA->point())
-                      : topoint(vB->point())
-                    ),
-                    magSqr(topoint(vA->point()) - topoint(vB->point())),
+                    cellCentre,
+                    pt,
                     surfHit,
                     hitSurface
                 );
@@ -610,11 +687,7 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation()
 
                     addSurfaceAndEdgeHits
                     (
-                        (
-                            vA->internalPoint()
-                          ? topoint(vA->point())
-                          : topoint(vB->point())
-                        ),
+                        pt,
                         surfaceIntersections,
                         surfacePtReplaceDistCoeffSqr,
                         edgeSearchDistCoeffSqr,
@@ -1613,6 +1686,7 @@ void Foam::conformalVoronoiMesh::limitDisplacement
     // Do not allow infinite recursion
     if (callCount > 7)
     {
+        displacement = vector::zero;
         return;
     }
 
@@ -1661,6 +1735,7 @@ void Foam::conformalVoronoiMesh::limitDisplacement
             if (magSqr(pt - surfHit.hitPoint()) <= searchDistanceSqr)
             {
                 // Cannot limit displacement, point closer than tolerance
+                displacement = vector::zero;
                 return;
             }
         }
@@ -1748,17 +1823,19 @@ bool Foam::conformalVoronoiMesh::nearSurfacePoint
     const Foam::point& pt = pHit.first().hitPoint();
 
     pointIndexHit closePoint;
-
     const bool closeToSurfacePt = pointIsNearSurfaceLocation(pt, closePoint);
 
     if
     (
         closeToSurfacePt
-     && mag(pt - closePoint.hitPoint()) > pointPairDistance(pt)
+     && (
+            magSqr(pt - closePoint.hitPoint())
+          > sqr(pointPairDistance(pt))
+        )
     )
     {
-        const scalar cosAngle
-            = angleBetweenSurfacePoints(pt, closePoint.hitPoint());
+        const scalar cosAngle =
+            angleBetweenSurfacePoints(pt, closePoint.hitPoint());
 
         // @todo make this tolerance run-time selectable?
         if (cosAngle < searchAngleOppositeSurface)
@@ -1768,11 +1845,6 @@ bool Foam::conformalVoronoiMesh::nearSurfacePoint
 
             const scalar searchDist = targetCellSize(closePoint.hitPoint());
 
-            if (searchDist < SMALL)
-            {
-                Pout<< "WARNING: SMALL CELL SIZE" << endl;
-            }
-
             geometryToConformTo_.findSurfaceNearest
             (
                 closePoint.hitPoint(),
@@ -1789,15 +1861,15 @@ bool Foam::conformalVoronoiMesh::nearSurfacePoint
                 norm
             );
 
-            const vector nA = norm[0];
+            const vector& nA = norm[0];
 
             pointIndexHit oppositeHit;
             label oppositeSurfaceHit = -1;
 
             geometryToConformTo_.findSurfaceNearestIntersection
             (
-                closePoint.hitPoint() + SMALL*nA,
-                closePoint.hitPoint() + mag(pt - closePoint.hitPoint())*nA,
+                closePoint.hitPoint() + 0.5*pointPairDistance(pt)*nA,
+                closePoint.hitPoint() + 5*targetCellSize(pt)*nA,
                 oppositeHit,
                 oppositeSurfaceHit
             );
@@ -1808,21 +1880,10 @@ bool Foam::conformalVoronoiMesh::nearSurfacePoint
                 pHit.first() = oppositeHit;
                 pHit.second() = oppositeSurfaceHit;
 
-//                appendToSurfacePtTree(pHit.first().hitPoint());
-//                surfaceToTreeShape.append
-//                (
-//                    existingSurfacePtLocations_.size() - 1
-//                );
-
                 return !closeToSurfacePt;
             }
         }
     }
-    else
-    {
-//        appendToSurfacePtTree(pt);
-//        surfaceToTreeShape.append(existingSurfacePtLocations_.size() - 1);
-    }
 
     return closeToSurfacePt;
 }
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H
index 64534cb839..e9e03a6330 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H
@@ -230,13 +230,13 @@ inline void Foam::conformalVoronoiMesh::createPointPair
 {
     vector ppDistn = ppDist*n;
 
-    const Foam::point internalPt = surfPt - ppDistn;
-    const Foam::point externalPt = surfPt + ppDistn;
+//    const Foam::point internalPt = surfPt - ppDistn;
+//    const Foam::point externalPt = surfPt + ppDistn;
 
-    bool internalInside = geometryToConformTo_.inside(internalPt);
-    bool externalOutside = geometryToConformTo_.outside(externalPt);
+//    bool internalInside = geometryToConformTo_.inside(internalPt);
+//    bool externalOutside = geometryToConformTo_.outside(externalPt);
 
-    if (internalInside && externalOutside)
+//    if (internalInside && externalOutside)
     {
         pts.append
         (
@@ -259,15 +259,24 @@ inline void Foam::conformalVoronoiMesh::createPointPair
                 Pstream::myProcNo()
             )
         );
+
+//        if (ptPair)
+//        {
+//            ptPairs_.addPointPair
+//            (
+//                pts[pts.size() - 2].index(),
+//                pts[pts.size() - 1].index() // external 0 -> slave
+//            );
+//        }
     }
-    else
-    {
-        Info<< "Warning: point pair not inside/outside" << nl
-            << "    surfPt = " << surfPt << nl
-            << "    internal = " << internalPt << " " << internalInside << nl
-            << "    external = " << externalPt << " " << externalOutside
-            << endl;
-    }
+//    else
+//    {
+//        Info<< "Warning: point pair not inside/outside" << nl
+//            << "    surfPt = " << surfPt << nl
+//            << "    internal = " << internalPt << " " << internalInside << nl
+//            << "    external = " << externalPt << " " << externalOutside
+//            << endl;
+//    }
 }
 
 
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell/indexedCell.H b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell/indexedCell.H
index 7a95b036a9..82bbe43b95 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell/indexedCell.H
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell/indexedCell.H
@@ -183,6 +183,7 @@ public:
         inline bool hasSeedPoint() const;
 
         inline bool hasInternalPoint() const;
+        inline bool hasBoundaryPoint() const;
 
         inline bool hasConstrainedPoint() const;
 
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell/indexedCellI.H b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell/indexedCellI.H
index 8e035dbc7b..2f5d076138 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell/indexedCellI.H
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell/indexedCellI.H
@@ -228,6 +228,19 @@ inline bool CGAL::indexedCell<Gt, Cb>::hasInternalPoint() const
 }
 
 
+template<class Gt, class Cb>
+inline bool CGAL::indexedCell<Gt, Cb>::hasBoundaryPoint() const
+{
+    return
+    (
+        this->vertex(0)->boundaryPoint()
+     || this->vertex(1)->boundaryPoint()
+     || this->vertex(2)->boundaryPoint()
+     || this->vertex(3)->boundaryPoint()
+    );
+}
+
+
 template<class Gt, class Cb>
 inline bool CGAL::indexedCell<Gt, Cb>::hasConstrainedPoint() const
 {
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
index 0360cf0d83..73bd6038d9 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
@@ -683,14 +683,18 @@ Foam::Field<bool> Foam::conformationSurfaces::wellInOutSide
                 vector hitDir = info[0].rawPoint() - samplePts[i];
                 hitDir /= mag(hitDir) + SMALL;
 
-                if
+                pointIndexHit surfHit;
+                label hitSurface;
+
+                findSurfaceNearestIntersection
                 (
-                    findSurfaceAnyIntersection
-                    (
-                        samplePts[i],
-                        info[0].rawPoint() - 1e-3*mag(hitDir)*(hitDir)
-                    )
-                )
+                    samplePts[i],
+                    info[0].rawPoint(),
+                    surfHit,
+                    hitSurface
+                );
+
+                if (surfHit.hit() && hitSurface != surfaces_[s])
                 {
                     continue;
                 }
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.C
index ebecad88e1..66dc10737d 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.C
@@ -37,6 +37,90 @@ namespace Foam
 defineTypeNameAndDebug(rayShooting, 0);
 addToRunTimeSelectionTable(initialPointsMethod, rayShooting, dictionary);
 
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void rayShooting::splitLine
+(
+    const line<point, point>& l,
+    const scalar& pert,
+    DynamicList<Vb::Point>& initialPoints
+) const
+{
+    Foam::point midPoint(l.centre());
+    const scalar localCellSize(cellShapeControls().cellSize(midPoint));
+    const scalar lineLength(l.mag());
+
+    const scalar minDistFromSurfaceSqr
+    (
+        minimumSurfaceDistanceCoeffSqr_
+       *sqr(localCellSize)
+    );
+
+    if
+    (
+        magSqr(midPoint - l.start()) > minDistFromSurfaceSqr
+     && magSqr(midPoint - l.end()) > minDistFromSurfaceSqr
+    )
+    {
+        // Add extra points if line length is much bigger than local cell size
+//        if (lineLength > 4.0*localCellSize)
+//        {
+//            splitLine
+//            (
+//                line<point, point>(l.start(), midPoint),
+//                pert,
+//                initialPoints
+//            );
+//
+//            splitLine
+//            (
+//                line<point, point>(midPoint, l.end()),
+//                pert,
+//                initialPoints
+//            );
+//        }
+
+        if (randomiseInitialGrid_)
+        {
+            Foam::point newPt
+            (
+                midPoint.x() + pert*(rndGen().scalar01() - 0.5),
+                midPoint.y() + pert*(rndGen().scalar01() - 0.5),
+                midPoint.z() + pert*(rndGen().scalar01() - 0.5)
+            );
+
+            if
+            (
+                !geometryToConformTo().findSurfaceAnyIntersection
+                (
+                    midPoint,
+                    newPt
+                )
+            )
+            {
+                midPoint = newPt;
+            }
+            else
+            {
+                WarningIn
+                (
+                    "rayShooting::splitLine"
+                    "("
+                    "   const line<point,point>&,"
+                    "   const scalar&,"
+                    "   DynamicList<Vb::Point>&"
+                    ")"
+                )   << "Point perturbation crosses a surface. Not inserting."
+                    << endl;
+            }
+        }
+
+        initialPoints.append(toPoint(midPoint));
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 rayShooting::rayShooting
@@ -75,7 +159,7 @@ List<Vb::Point> rayShooting::initialPoints() const
     const searchableSurfaces& surfaces = geometryToConformTo().geometry();
     const labelList& surfacesToConformTo = geometryToConformTo().surfaces();
 
-    const scalar maxRayLength = surfaces.bounds().mag();
+    const scalar maxRayLength(surfaces.bounds().mag());
 
     // Initialise points list
     label initialPointsSize = 0;
@@ -94,8 +178,7 @@ List<Vb::Point> rayShooting::initialPoints() const
         const pointField& faceCentres = faceCentresTmp();
 
         Info<< "    Shoot rays from " << s.name() << nl
-            << "    nRays = " << faceCentres.size() << endl;
-
+            << "        nRays = " << faceCentres.size() << endl;
 
         forAll(faceCentres, fcI)
         {
@@ -110,9 +193,11 @@ List<Vb::Point> rayShooting::initialPoints() const
                 continue;
             }
 
-            const scalar pert =
+            const scalar pert
+            (
                 randomPerturbationCoeff_
-               *cellShapeControls().cellSize(fC);
+               *cellShapeControls().cellSize(fC)
+            );
 
             pointIndexHit surfHitStart;
             label hitSurfaceStart;
@@ -181,27 +266,12 @@ List<Vb::Point> rayShooting::initialPoints() const
                         }
                     }
 
-                    Foam::point midPoint(l.centre());
-
-                    const scalar minDistFromSurfaceSqr =
-                        minimumSurfaceDistanceCoeffSqr_
-                       *sqr(cellShapeControls().cellSize(midPoint));
-
-                    if (randomiseInitialGrid_)
-                    {
-                        midPoint.x() += pert*(rndGen().scalar01() - 0.5);
-                        midPoint.y() += pert*(rndGen().scalar01() - 0.5);
-                        midPoint.z() += pert*(rndGen().scalar01() - 0.5);
-                    }
-
-                    if
+                    splitLine
                     (
-                        magSqr(midPoint - l.start()) > minDistFromSurfaceSqr
-                     && magSqr(midPoint - l.end()) > minDistFromSurfaceSqr
-                    )
-                    {
-                        initialPoints.append(toPoint(midPoint));
-                    }
+                        l,
+                        pert,
+                        initialPoints
+                    );
                 }
             }
         }
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.H b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.H
index e66df62b7c..e5bd1936ff 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.H
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/rayShooting/rayShooting.H
@@ -61,6 +61,16 @@ private:
         scalar randomPerturbationCoeff_;
 
 
+    // Private Member Functions
+
+        void splitLine
+        (
+            const line<point, point>& l,
+            const scalar& pert,
+            DynamicList<Vb::Point>& initialPoints
+        ) const;
+
+
 public:
 
     //- Runtime type information

From e8985cbd2954382f123f7bcd5a4e1af1004f7e49 Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Fri, 30 Aug 2013 16:31:22 +0100
Subject: [PATCH 13/20] BUG: Corrected commit 6a1b567

---
 .../Templates/KinematicParcel/KinematicParcel.C       | 11 ++++-------
 1 file changed, 4 insertions(+), 7 deletions(-)

diff --git a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
index 10bd38d62e..b833962825 100644
--- a/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
+++ b/src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C
@@ -287,16 +287,13 @@ bool Foam::KinematicParcel<ParcelType>::move
         const label cellI = p.cell();
 
         const scalar magU = mag(U_);
-        if (p.active())
+        if (p.active() && moving && (magU > ROOTVSMALL))
         {
             const scalar d = dt*magU;
             const scalar dCorr = min(d, maxCo*cbrt(V[cellI]));
-            dt *= dCorr/d;
-
-            if (moving && (magU > ROOTVSMALL))
-            {
-                dt *= p.trackToFace(p.position() + dCorr*U_/magU, td);
-            }
+            dt *=
+                dCorr/d
+               *p.trackToFace(p.position() + dCorr*U_/magU, td);
         }
 
         tEnd -= dt;

From 94da7c46a5103cca1529fa8ac0f6d45f4a5a68ed Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Mon, 2 Sep 2013 09:58:28 +0100
Subject: [PATCH 14/20] STYLE: Minor code formatting

---
 .../layerAdditionRemoval.C                    | 32 ++++++++-----------
 1 file changed, 14 insertions(+), 18 deletions(-)

diff --git a/src/dynamicMesh/layerAdditionRemoval/layerAdditionRemoval.C b/src/dynamicMesh/layerAdditionRemoval/layerAdditionRemoval.C
index b712a51c7f..b052b4de86 100644
--- a/src/dynamicMesh/layerAdditionRemoval/layerAdditionRemoval.C
+++ b/src/dynamicMesh/layerAdditionRemoval/layerAdditionRemoval.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
@@ -96,6 +96,7 @@ void Foam::layerAdditionRemoval::checkDefinition()
     }
 }
 
+
 Foam::scalar Foam::layerAdditionRemoval::readOldThickness
 (
     const dictionary& dict
@@ -115,7 +116,6 @@ void Foam::layerAdditionRemoval::clearAddressing() const
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 Foam::layerAdditionRemoval::layerAdditionRemoval
 (
     const word& name,
@@ -140,7 +140,6 @@ Foam::layerAdditionRemoval::layerAdditionRemoval
 }
 
 
-// Construct from dictionary
 Foam::layerAdditionRemoval::layerAdditionRemoval
 (
     const word& name,
@@ -377,7 +376,7 @@ void Foam::layerAdditionRemoval::setRefinement(polyTopoChange& ref) const
         // Clear addressing.  This also resets the addition/removal data
         if (debug)
         {
-            Pout<< "layerAdditionRemoval::setRefinement(polyTopoChange& ref) "
+            Pout<< "layerAdditionRemoval::setRefinement(polyTopoChange&) "
                 << " for object " << name() << " : "
                 << "Clearing addressing after layer removal. " << endl;
         }
@@ -393,7 +392,7 @@ void Foam::layerAdditionRemoval::setRefinement(polyTopoChange& ref) const
         // Clear addressing.  This also resets the addition/removal data
         if (debug)
         {
-            Pout<< "layerAdditionRemoval::setRefinement(polyTopoChange& ref) "
+            Pout<< "layerAdditionRemoval::setRefinement(polyTopoChange&) "
                 << " for object " << name() << " : "
                 << "Clearing addressing after layer addition. " << endl;
         }
@@ -431,16 +430,14 @@ void Foam::layerAdditionRemoval::updateMesh(const mapPolyMesh&)
 
 void Foam::layerAdditionRemoval::setMinLayerThickness(const scalar t) const
 {
-    if
-    (
-        t < VSMALL
-     || maxLayerThickness_ < t
-    )
+    if (t < VSMALL || maxLayerThickness_ < t)
     {
         FatalErrorIn
         (
-            "void layerAdditionRemoval::setMinLayerThickness("
-            "const scalar t) const"
+            "void layerAdditionRemoval::setMinLayerThickness"
+            "("
+                "const scalar"
+            ") const"
         )   << "Incorrect layer thickness definition."
             << abort(FatalError);
     }
@@ -451,15 +448,14 @@ void Foam::layerAdditionRemoval::setMinLayerThickness(const scalar t) const
 
 void Foam::layerAdditionRemoval::setMaxLayerThickness(const scalar t) const
 {
-    if
-    (
-        t < minLayerThickness_
-    )
+    if (t < minLayerThickness_)
     {
         FatalErrorIn
         (
-            "void layerAdditionRemoval::setMaxLayerThickness("
-            "const scalar t) const"
+            "void layerAdditionRemoval::setMaxLayerThickness"
+            "("
+                "const scalar"
+            ") const"
         )   << "Incorrect layer thickness definition."
             << abort(FatalError);
     }

From 57019e885532d9c511bbc5731522bd4a9e1c7dad Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Mon, 2 Sep 2013 11:31:23 +0100
Subject: [PATCH 15/20] ENH: Added checking for mass fractions summing to zero

---
 .../multiComponentMixture/multiComponentMixture.C | 15 +++++++++++++--
 1 file changed, 13 insertions(+), 2 deletions(-)

diff --git a/src/thermophysicalModels/reactionThermo/mixtures/multiComponentMixture/multiComponentMixture.C b/src/thermophysicalModels/reactionThermo/mixtures/multiComponentMixture/multiComponentMixture.C
index dbdc804d8b..f10a1ce48b 100644
--- a/src/thermophysicalModels/reactionThermo/mixtures/multiComponentMixture/multiComponentMixture.C
+++ b/src/thermophysicalModels/reactionThermo/mixtures/multiComponentMixture/multiComponentMixture.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
@@ -49,7 +49,7 @@ const ThermoType& Foam::multiComponentMixture<ThermoType>::constructSpeciesData
 template<class ThermoType>
 void Foam::multiComponentMixture<ThermoType>::correctMassFractions()
 {
-    // It changes Yt patches to "calculated"
+    // Multiplication by 1.0 changes Yt patches to "calculated"
     volScalarField Yt("Yt", 1.0*Y_[0]);
 
     for (label n=1; n<Y_.size(); n++)
@@ -57,6 +57,17 @@ void Foam::multiComponentMixture<ThermoType>::correctMassFractions()
         Yt += Y_[n];
     }
 
+    if (mag(max(Yt).value()) < ROOTVSMALL)
+    {
+        FatalErrorIn
+        (
+            "void Foam::multiComponentMixture<ThermoType>::"
+            "correctMassFractions()"
+        )
+            << "Sum of mass fractions is zero for species " << this->species()
+            << exit(FatalError);
+    }
+
     forAll(Y_, n)
     {
         Y_[n] /= Yt;

From d9805740c70aad4bf5ba039e34c2717df265413d Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Mon, 2 Sep 2013 14:25:29 +0100
Subject: [PATCH 16/20] ENH: Tutorial update

---
 tutorials/incompressible/icoFoam/cavity/0/U | 13 ++-----------
 1 file changed, 2 insertions(+), 11 deletions(-)

diff --git a/tutorials/incompressible/icoFoam/cavity/0/U b/tutorials/incompressible/icoFoam/cavity/0/U
index 21b0d542ec..711702f987 100644
--- a/tutorials/incompressible/icoFoam/cavity/0/U
+++ b/tutorials/incompressible/icoFoam/cavity/0/U
@@ -22,17 +22,8 @@ boundaryField
 {
     movingWall
     {
-        type            uniformFixedValue;
-        uniformValue    csvFile;
-        csvFileCoeffs
-        {
-            nHeaderLine     1;
-            refColumn       0;
-            componentColumns (1 2 3);
-            mergeSeparators no;
-            fileName        "$FOAM_CASE/csvData";
-            separator       " ";
-        }
+        type            fixedValue;
+        value           uniform (1 0 0);
     }
 
     fixedWalls

From 203c79307151f85376a912c7566ad8b89b74d4f0 Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Mon, 2 Sep 2013 14:25:44 +0100
Subject: [PATCH 17/20] 
 tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D: Updated for new
 phase naming

---
 .../0.org/{alphamercury => alpha.air}         |  2 +-
 .../0.org/{alphas => alpha.mercury}           |  2 +-
 .../0.org/{alphaoil => alpha.oil}             |  2 +-
 .../0.org/{alphaair => alpha.water}           |  0
 .../mixerVessel2D/0.org/alphawater            | 44 -------------------
 .../mixerVessel2D/0/{alphaair => alpha.air}   |  2 +-
 .../0/{alphamercury => alpha.mercury}         |  2 +-
 .../mixerVessel2D/0/{alphaoil => alpha.oil}   |  2 +-
 .../0/{alphawater => alpha.water}             |  2 +-
 .../mixerVessel2D/constant/polyMesh/boundary  |  2 +
 .../constant/transportProperties              |  2 -
 .../mixerVessel2D/system/setFieldsDict        | 32 +++++++-------
 12 files changed, 25 insertions(+), 69 deletions(-)
 rename tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/{alphamercury => alpha.air} (97%)
 rename tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/{alphas => alpha.mercury} (97%)
 rename tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/{alphaoil => alpha.oil} (97%)
 rename tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/{alphaair => alpha.water} (100%)
 delete mode 100644 tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alphawater
 rename tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/{alphaair => alpha.air} (99%)
 rename tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/{alphamercury => alpha.mercury} (99%)
 rename tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/{alphaoil => alpha.oil} (99%)
 rename tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/{alphawater => alpha.water} (99%)

diff --git a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alphamercury b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alpha.air
similarity index 97%
rename from tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alphamercury
rename to tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alpha.air
index 0f3dcd4639..210222182a 100644
--- a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alphamercury
+++ b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alpha.air
@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    object      alpha1;
+    object      alpha.air;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alphas b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alpha.mercury
similarity index 97%
rename from tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alphas
rename to tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alpha.mercury
index 0f3dcd4639..cbb4647b47 100644
--- a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alphas
+++ b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alpha.mercury
@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    object      alpha1;
+    object      alpha.mercury;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alphaoil b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alpha.oil
similarity index 97%
rename from tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alphaoil
rename to tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alpha.oil
index 0f3dcd4639..e8b4d42497 100644
--- a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alphaoil
+++ b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alpha.oil
@@ -10,7 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    object      alpha1;
+    object      alpha.oil;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alphaair b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alpha.water
similarity index 100%
rename from tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alphaair
rename to tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alpha.water
diff --git a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alphawater b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alphawater
deleted file mode 100644
index 0f3dcd4639..0000000000
--- a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0.org/alphawater
+++ /dev/null
@@ -1,44 +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       volScalarField;
-    object      alpha1;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [0 0 0 0 0 0 0];
-
-internalField   uniform 0;
-
-boundaryField
-{
-    rotor
-    {
-        type            zeroGradient;
-    }
-
-    stator
-    {
-        type            zeroGradient;
-    }
-
-    front
-    {
-        type            empty;
-    }
-
-    back
-    {
-        type            empty;
-    }
-}
-
-// ************************************************************************* //
diff --git a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alphaair b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alpha.air
similarity index 99%
rename from tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alphaair
rename to tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alpha.air
index b25f228620..a53034d33d 100644
--- a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alphaair
+++ b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alpha.air
@@ -11,7 +11,7 @@ FoamFile
     format      ascii;
     class       volScalarField;
     location    "0";
-    object      alphaair;
+    object      alpha.air;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alphamercury b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alpha.mercury
similarity index 99%
rename from tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alphamercury
rename to tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alpha.mercury
index 0e6d124fd0..76002901ed 100644
--- a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alphamercury
+++ b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alpha.mercury
@@ -11,7 +11,7 @@ FoamFile
     format      ascii;
     class       volScalarField;
     location    "0";
-    object      alphamercury;
+    object      alpha.mercury;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alphaoil b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alpha.oil
similarity index 99%
rename from tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alphaoil
rename to tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alpha.oil
index 2311855110..56d642d0d9 100644
--- a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alphaoil
+++ b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alpha.oil
@@ -11,7 +11,7 @@ FoamFile
     format      ascii;
     class       volScalarField;
     location    "0";
-    object      alphaoil;
+    object      alpha.oil;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alphawater b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alpha.water
similarity index 99%
rename from tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alphawater
rename to tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alpha.water
index d61eb80ddf..766d8c1462 100644
--- a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alphawater
+++ b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/0/alpha.water
@@ -11,7 +11,7 @@ FoamFile
     format      ascii;
     class       volScalarField;
     location    "0";
-    object      alphawater;
+    object      alpha.water;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/constant/polyMesh/boundary b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/constant/polyMesh/boundary
index 292f25b806..188a0f0c58 100644
--- a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/constant/polyMesh/boundary
+++ b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/constant/polyMesh/boundary
@@ -32,12 +32,14 @@ FoamFile
     front
     {
         type            empty;
+        inGroups        1(empty);
         nFaces          3072;
         startFace       6336;
     }
     back
     {
         type            empty;
+        inGroups        1(empty);
         nFaces          3072;
         startFace       9408;
     }
diff --git a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/constant/transportProperties b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/constant/transportProperties
index e338328156..3cfdbabfaf 100644
--- a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/constant/transportProperties
+++ b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/constant/transportProperties
@@ -46,8 +46,6 @@ phases
      }
 );
 
-refPhase        air;
-
 sigmas
 (
     (air water) 0.07
diff --git a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/system/setFieldsDict b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/system/setFieldsDict
index d5dd8aca43..82d9ffc783 100644
--- a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/system/setFieldsDict
+++ b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/system/setFieldsDict
@@ -17,10 +17,10 @@ FoamFile
 
 defaultFieldValues
 (
-    volScalarFieldValue alphaair 1
-    volScalarFieldValue alphawater 0
-    volScalarFieldValue alphaoil 0
-    volScalarFieldValue alphamercury 0
+    volScalarFieldValue alpha.air 1
+    volScalarFieldValue alpha.water 0
+    volScalarFieldValue alpha.oil 0
+    volScalarFieldValue alpha.mercury 0
 );
 
 regions
@@ -30,10 +30,10 @@ regions
         box (0 0 -1) (1 1 1);
         fieldValues
         (
-            volScalarFieldValue alphawater 1
-            volScalarFieldValue alphaoil 0
-            volScalarFieldValue alphamercury 0
-            volScalarFieldValue alphaair 0
+            volScalarFieldValue alpha.water 1
+            volScalarFieldValue alpha.oil 0
+            volScalarFieldValue alpha.mercury 0
+            volScalarFieldValue alpha.air 0
         );
     }
     boxToCell
@@ -41,10 +41,10 @@ regions
         box (0 -1 -1) (1 0 1);
         fieldValues
         (
-            volScalarFieldValue alphawater 0
-            volScalarFieldValue alphaoil 1
-            volScalarFieldValue alphamercury 0
-            volScalarFieldValue alphaair 0
+            volScalarFieldValue alpha.water 0
+            volScalarFieldValue alpha.oil 1
+            volScalarFieldValue alpha.mercury 0
+            volScalarFieldValue alpha.air 0
         );
     }
     boxToCell
@@ -52,10 +52,10 @@ regions
         box (-1 -1 -1) (0 0 1);
         fieldValues
         (
-            volScalarFieldValue alphawater 0
-            volScalarFieldValue alphaoil 0
-            volScalarFieldValue alphamercury 1
-            volScalarFieldValue alphaair 0
+            volScalarFieldValue alpha.water 0
+            volScalarFieldValue alpha.oil 0
+            volScalarFieldValue alpha.mercury 1
+            volScalarFieldValue alpha.air 0
         );
     }
 );

From 604fb95172cd2eac932912cdc62913401f0cbadc Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Mon, 2 Sep 2013 18:32:05 +0100
Subject: [PATCH 18/20] BUG: Corrected Pilch-Erdman break-up model

---
 .../BreakupModel/PilchErdman/PilchErdman.C    | 75 ++++++++++++-------
 .../BreakupModel/PilchErdman/PilchErdman.H    | 22 +++++-
 2 files changed, 68 insertions(+), 29 deletions(-)

diff --git a/src/lagrangian/spray/submodels/BreakupModel/PilchErdman/PilchErdman.C b/src/lagrangian/spray/submodels/BreakupModel/PilchErdman/PilchErdman.C
index d2ca617bbd..e6863ad0db 100644
--- a/src/lagrangian/spray/submodels/BreakupModel/PilchErdman/PilchErdman.C
+++ b/src/lagrangian/spray/submodels/BreakupModel/PilchErdman/PilchErdman.C
@@ -36,7 +36,7 @@ Foam::PilchErdman<CloudType>::PilchErdman
 :
     BreakupModel<CloudType>(dict, owner, typeName),
     B1_(0.375),
-    B2_(0.236)
+    B2_(0.2274)
 {
     if (!this->defaultCoeffs(true))
     {
@@ -90,58 +90,79 @@ bool Foam::PilchErdman<CloudType>::update
     scalar& massChild
 )
 {
-    scalar semiMass = nParticle*pow3(d);
-    scalar We = 0.5*rhoc*sqr(Urmag)*d/sigma;
+    // Weber number - eq (1)
+    scalar We = rhoc*sqr(Urmag)*d/sigma;
+
+    // Ohnesorge number - eq (2)
     scalar Oh = mu/sqrt(rho*d*sigma);
 
-    scalar Wec = 6.0*(1.0 + 1.077*pow(Oh, 1.6));
+    // Critical Weber number - eq (5)
+    scalar Wec = 12.0*(1.0 + 1.077*pow(Oh, 1.6));
 
     if (We > Wec)
     {
-        // We > 1335, wave crest stripping
+        // We > 2670, wave crest stripping - eq (12)
         scalar taubBar = 5.5;
 
-        if (We < 1335)
+        if (We < 2670)
         {
-            if (We > 175.0)
+            if (We > 351)
             {
-                // sheet stripping
-                taubBar = 0.766*pow(2.0*We - 12.0, 0.25);
+                // sheet stripping - eq (11)
+                taubBar = 0.766*pow(We - 12.0, 0.25);
             }
-            else if (We > 22.0)
+            else if (We > 45)
             {
-                // Bag-and-stamen breakup
-                taubBar = 14.1*pow(2.0*We - 12.0, -0.25);
+                // bag-and-stamen breakup  - eq (10)
+                taubBar = 14.1*pow(We - 12.0, 0.25);
             }
-            else if (We > 9.0)
+            else if (We > 18)
             {
-                // Bag breakup
-                taubBar = 2.45*pow(2.0*We - 12.0, 0.25);
+                // bag breakup - eq (9)
+                taubBar = 2.45*pow(We - 12.0, 0.25);
             }
-            else if (We > 6.0)
+            else if (We > 12)
             {
-                // Vibrational breakup
-                taubBar = 6.0*pow(2.0*We - 12.0, -0.25);
+                // vibrational breakup - eq (8)
+                taubBar = 6.0*pow(We - 12.0, -0.25);
+            }
+            else
+            {
+                // no break-up
+                taubBar = GREAT;
             }
         }
 
         scalar rho12 = sqrt(rhoc/rho);
 
-        scalar Vd = Urmag*rho12*(B1_*taubBar + B2_*taubBar*taubBar);
+        // velocity of fragmenting drop - eq (20)
+        scalar Vd = Urmag*rho12*(B1_*taubBar + B2_*sqr(taubBar));
+
+        // maximum stable diameter - eq (33)
         scalar Vd1 = sqr(1.0 - Vd/Urmag);
         Vd1 = max(Vd1, SMALL);
-        scalar Ds = 2.0*Wec*sigma/(Vd1*rhoc*sqr(Urmag));
-        scalar A = Urmag*rho12/d;
+        scalar dStable = Wec*sigma/(Vd1*rhoc*sqr(Urmag));
 
-        scalar taub = taubBar/A;
+        if (d < dStable)
+        {
+            // droplet diameter already stable = no break-up
+            // - do not update d and nParticle
+            return false;
+        }
+        else
+        {
+            scalar semiMass = nParticle*pow3(d);
 
-        scalar frac = dt/taub;
+            // invert eq (3) to create a dimensional break-up time
+            scalar taub = taubBar*d/(Urmag*rho12);
 
-        // update the droplet diameter according to the rate eq. (implicitly)
-        d = (d + frac*Ds)/(1.0 + frac);
+            // update droplet diameter according to the rate eq (implicitly)
+            scalar frac = dt/taub;
+            d = (d + frac*dStable)/(1.0 + frac);
 
-        // correct the number of particles to conserve mass
-        nParticle = semiMass/pow3(d);
+            // correct the number of particles to conserve mass
+            nParticle = semiMass/pow3(d);
+        }
     }
 
     return false;
diff --git a/src/lagrangian/spray/submodels/BreakupModel/PilchErdman/PilchErdman.H b/src/lagrangian/spray/submodels/BreakupModel/PilchErdman/PilchErdman.H
index 830c1feeb5..09d1f3fe03 100644
--- a/src/lagrangian/spray/submodels/BreakupModel/PilchErdman/PilchErdman.H
+++ b/src/lagrangian/spray/submodels/BreakupModel/PilchErdman/PilchErdman.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
@@ -25,7 +25,7 @@ Class
     Foam::PilchErdman
 
 Description
-    secondary breakup model
+    Particle secondary breakup model, based on the reference:
 
     @verbatim
     Pilch, M. and Erdman, C.A.
@@ -35,6 +35,24 @@ Description
     Int. J. Multiphase Flows 13 (1987), 741-757
     @endverbatim
 
+    The droplet fragment velocity is described by the equation:
+
+    \f[
+        V_d = V sqrt(epsilon)(B1 T + B2 T^2)
+    \f]
+
+    Where:
+        V_d     : fragment velocity
+        V       : magnitude of the relative velocity
+        epsilon : density ratio (rho_carrier/rho_droplet)
+        T       : characteristic break-up time
+        B1, B2  : model input coefficients
+
+    The authors suggest that:
+        compressible flow   : B1 = 0.75*1.0; B2 = 3*0.116
+        incompressible flow : B1 = 0.75*0.5; B2 = 3*0.0758
+
+
 \*---------------------------------------------------------------------------*/
 
 #ifndef PilchErdman_H

From ef2dcd9ad145c495049fca26eb82c9ba6708c559 Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Tue, 3 Sep 2013 11:32:26 +0100
Subject: [PATCH 19/20] BUG: layer addition/removal - updated parallel
 behaviour

---
 .../layerAdditionRemoval.C                    | 179 ++++++++++--------
 .../layerAdditionRemoval.H                    |  16 +-
 2 files changed, 108 insertions(+), 87 deletions(-)

diff --git a/src/dynamicMesh/layerAdditionRemoval/layerAdditionRemoval.C b/src/dynamicMesh/layerAdditionRemoval/layerAdditionRemoval.C
index b052b4de86..f99612eed6 100644
--- a/src/dynamicMesh/layerAdditionRemoval/layerAdditionRemoval.C
+++ b/src/dynamicMesh/layerAdditionRemoval/layerAdditionRemoval.C
@@ -58,10 +58,8 @@ void Foam::layerAdditionRemoval::checkDefinition()
 {
     if (!faceZoneID_.active())
     {
-        FatalErrorIn
-        (
-            "void Foam::layerAdditionRemoval::checkDefinition()"
-        )   << "Master face zone named " << faceZoneID_.name()
+        FatalErrorIn("void Foam::layerAdditionRemoval::checkDefinition()")
+            << "Master face zone named " << faceZoneID_.name()
             << " cannot be found."
             << abort(FatalError);
     }
@@ -79,13 +77,15 @@ void Foam::layerAdditionRemoval::checkDefinition()
             << abort(FatalError);
     }
 
-    if (topoChanger().mesh().faceZones()[faceZoneID_.index()].empty())
+    label nFaces = topoChanger().mesh().faceZones()[faceZoneID_.index()].size();
+
+    reduce(nFaces, sumOp<label>());
+
+    if (nFaces == 0)
     {
-        FatalErrorIn
-        (
-            "void Foam::layerAdditionRemoval::checkDefinition()"
-        )   << "Face extrusion zone contains no faces. "
-            << " Please check your mesh definition."
+        FatalErrorIn("void Foam::layerAdditionRemoval::checkDefinition()")
+            << "Face extrusion zone contains no faces. "
+            << "Please check your mesh definition."
             << abort(FatalError);
     }
 
@@ -120,16 +120,18 @@ Foam::layerAdditionRemoval::layerAdditionRemoval
 (
     const word& name,
     const label index,
-    const polyTopoChanger& mme,
+    const polyTopoChanger& ptc,
     const word& zoneName,
     const scalar minThickness,
-    const scalar maxThickness
+    const scalar maxThickness,
+    const Switch thicknessFromVolume
 )
 :
-    polyMeshModifier(name, index, mme, true),
-    faceZoneID_(zoneName, mme.mesh().faceZones()),
+    polyMeshModifier(name, index, ptc, true),
+    faceZoneID_(zoneName, ptc.mesh().faceZones()),
     minLayerThickness_(minThickness),
     maxLayerThickness_(maxThickness),
+    thicknessFromVolume_(thicknessFromVolume),
     oldLayerThickness_(-1.0),
     pointsPairingPtr_(NULL),
     facesPairingPtr_(NULL),
@@ -145,13 +147,17 @@ Foam::layerAdditionRemoval::layerAdditionRemoval
     const word& name,
     const dictionary& dict,
     const label index,
-    const polyTopoChanger& mme
+    const polyTopoChanger& ptc
 )
 :
-    polyMeshModifier(name, index, mme, Switch(dict.lookup("active"))),
-    faceZoneID_(dict.lookup("faceZoneName"), mme.mesh().faceZones()),
+    polyMeshModifier(name, index, ptc, Switch(dict.lookup("active"))),
+    faceZoneID_(dict.lookup("faceZoneName"), ptc.mesh().faceZones()),
     minLayerThickness_(readScalar(dict.lookup("minLayerThickness"))),
     maxLayerThickness_(readScalar(dict.lookup("maxLayerThickness"))),
+    thicknessFromVolume_
+    (
+        dict.lookupOrDefault<Switch>("thicknessFromVolume", true)
+    ),
     oldLayerThickness_(readOldThickness(dict)),
     pointsPairingPtr_(NULL),
     facesPairingPtr_(NULL),
@@ -188,11 +194,13 @@ bool Foam::layerAdditionRemoval::changeTopology() const
     // Layer removal:
     //     When the min thickness falls below the threshold, trigger removal.
 
-    const faceZone& fz = topoChanger().mesh().faceZones()[faceZoneID_.index()];
+    const polyMesh& mesh = topoChanger().mesh();
+
+    const faceZone& fz = mesh.faceZones()[faceZoneID_.index()];
     const labelList& mc = fz.masterCells();
 
-    const scalarField& V = topoChanger().mesh().cellVolumes();
-    const vectorField& S = topoChanger().mesh().faceAreas();
+    const scalarField& V = mesh.cellVolumes();
+    const vectorField& S = mesh.faceAreas();
 
     if (min(V) < -VSMALL)
     {
@@ -205,63 +213,68 @@ bool Foam::layerAdditionRemoval::changeTopology() const
     scalar avgDelta = 0;
     scalar minDelta = GREAT;
     scalar maxDelta = 0;
+    label nDelta = 0;
 
-    forAll(fz, faceI)
+    if (thicknessFromVolume_)
     {
-        scalar curDelta = V[mc[faceI]]/mag(S[fz[faceI]]);
-        avgDelta += curDelta;
-        minDelta = min(minDelta, curDelta);
-        maxDelta = max(maxDelta, curDelta);
+        // Thickness calculated from cell volume/face area
+        forAll(fz, faceI)
+        {
+            scalar curDelta = V[mc[faceI]]/mag(S[fz[faceI]]);
+            avgDelta += curDelta;
+            minDelta = min(minDelta, curDelta);
+            maxDelta = max(maxDelta, curDelta);
+        }
+
+        nDelta = fz.size();
+    }
+    else
+    {
+        // Thickness calculated from edges on layer
+        const Map<label>& zoneMeshPointMap = fz().meshPointMap();
+
+        // Edges with only one point on zone
+        forAll(mc, faceI)
+        {
+            const cell& cFaces = mesh.cells()[mc[faceI]];
+            const edgeList cellEdges(cFaces.edges(mesh.faces()));
+
+            forAll(cellEdges, i)
+            {
+                const edge& e = cellEdges[i];
+
+                if (zoneMeshPointMap.found(e[0]))
+                {
+                    if (!zoneMeshPointMap.found(e[1]))
+                    {
+                        scalar curDelta = e.mag(mesh.points());
+                        avgDelta += curDelta;
+                        nDelta++;
+                        minDelta = min(minDelta, curDelta);
+                        maxDelta = max(maxDelta, curDelta);
+                    }
+                }
+                else
+                {
+                    if (zoneMeshPointMap.found(e[1]))
+                    {
+                        scalar curDelta = e.mag(mesh.points());
+                        avgDelta += curDelta;
+                        nDelta++;
+                        minDelta = min(minDelta, curDelta);
+                        maxDelta = max(maxDelta, curDelta);
+                    }
+                }
+            }
+        }
     }
 
-    avgDelta /= fz.size();
+    reduce(minDelta, minOp<scalar>());
+    reduce(maxDelta, maxOp<scalar>());
+    reduce(avgDelta, sumOp<scalar>());
+    reduce(nDelta, sumOp<label>());
 
-    ////MJ Alternative thickness determination
-    //{
-    //    // Edges on layer.
-    //    const Map<label>& zoneMeshPointMap = fz().meshPointMap();
-    //
-    //    label nDelta = 0;
-    //
-    //    // Edges with only one point on zone
-    //    const polyMesh& mesh = topoChanger().mesh();
-    //
-    //    forAll(mc, faceI)
-    //    {
-    //        const cell& cFaces = mesh.cells()[mc[faceI]];
-    //        const edgeList cellEdges(cFaces.edges(mesh.faces()));
-    //
-    //        forAll(cellEdges, i)
-    //        {
-    //            const edge& e = cellEdges[i];
-    //
-    //            if (zoneMeshPointMap.found(e[0]))
-    //            {
-    //                if (!zoneMeshPointMap.found(e[1]))
-    //                {
-    //                    scalar curDelta = e.mag(mesh.points());
-    //                    avgDelta += curDelta;
-    //                    nDelta++;
-    //                    minDelta = min(minDelta, curDelta);
-    //                    maxDelta = max(maxDelta, curDelta);
-    //                }
-    //            }
-    //            else
-    //            {
-    //                if (zoneMeshPointMap.found(e[1]))
-    //                {
-    //                    scalar curDelta = e.mag(mesh.points());
-    //                    avgDelta += curDelta;
-    //                    nDelta++;
-    //                    minDelta = min(minDelta, curDelta);
-    //                    maxDelta = max(maxDelta, curDelta);
-    //                }
-    //            }
-    //        }
-    //    }
-    //
-    //    avgDelta /= nDelta;
-    //}
+    avgDelta /= nDelta;
 
     if (debug)
     {
@@ -286,7 +299,6 @@ bool Foam::layerAdditionRemoval::changeTopology() const
         }
 
         // No topological changes allowed before first mesh motion
-        //
         oldLayerThickness_ = avgDelta;
 
         topologicalChange = false;
@@ -314,7 +326,7 @@ bool Foam::layerAdditionRemoval::changeTopology() const
                             << "Triggering layer removal" << endl;
                     }
 
-                    triggerRemoval_ = topoChanger().mesh().time().timeIndex();
+                    triggerRemoval_ = mesh.time().timeIndex();
 
                     // Old thickness looses meaning.
                     // Set it up to indicate layer removal
@@ -346,7 +358,7 @@ bool Foam::layerAdditionRemoval::changeTopology() const
                     << "Triggering layer addition" << endl;
             }
 
-            triggerAddition_ = topoChanger().mesh().time().timeIndex();
+            triggerAddition_ = mesh.time().timeIndex();
 
             // Old thickness looses meaning.
             // Set it up to indicate layer removal
@@ -377,8 +389,8 @@ void Foam::layerAdditionRemoval::setRefinement(polyTopoChange& ref) const
         if (debug)
         {
             Pout<< "layerAdditionRemoval::setRefinement(polyTopoChange&) "
-                << " for object " << name() << " : "
-                << "Clearing addressing after layer removal. " << endl;
+                << "for object " << name() << " : "
+                << "Clearing addressing after layer removal" << endl;
         }
 
         triggerRemoval_ = -1;
@@ -393,8 +405,8 @@ void Foam::layerAdditionRemoval::setRefinement(polyTopoChange& ref) const
         if (debug)
         {
             Pout<< "layerAdditionRemoval::setRefinement(polyTopoChange&) "
-                << " for object " << name() << " : "
-                << "Clearing addressing after layer addition. " << endl;
+                << "for object " << name() << " : "
+                << "Clearing addressing after layer addition" << endl;
         }
 
         triggerAddition_ = -1;
@@ -408,8 +420,8 @@ void Foam::layerAdditionRemoval::updateMesh(const mapPolyMesh&)
     if (debug)
     {
         Pout<< "layerAdditionRemoval::updateMesh(const mapPolyMesh&) "
-            << " for object " << name() << " : "
-            << "Clearing addressing on external request. ";
+            << "for object " << name() << " : "
+            << "Clearing addressing on external request";
 
         if (pointsPairingPtr_ || facesPairingPtr_)
         {
@@ -471,7 +483,8 @@ void Foam::layerAdditionRemoval::write(Ostream& os) const
         << faceZoneID_ << nl
         << minLayerThickness_ << nl
         << oldLayerThickness_ << nl
-        << maxLayerThickness_ << endl;
+        << maxLayerThickness_ << nl
+        << thicknessFromVolume_ << endl;
 }
 
 
@@ -486,6 +499,8 @@ void Foam::layerAdditionRemoval::writeDict(Ostream& os) const
         << token::END_STATEMENT << nl
         << "    maxLayerThickness " << maxLayerThickness_
         << token::END_STATEMENT << nl
+        << "    thicknessFromVolume " << thicknessFromVolume_
+        << token::END_STATEMENT << nl
         << "    oldLayerThickness " << oldLayerThickness_
         << token::END_STATEMENT << nl
         << "    active " << active()
diff --git a/src/dynamicMesh/layerAdditionRemoval/layerAdditionRemoval.H b/src/dynamicMesh/layerAdditionRemoval/layerAdditionRemoval.H
index daef765097..edb3218eae 100644
--- a/src/dynamicMesh/layerAdditionRemoval/layerAdditionRemoval.H
+++ b/src/dynamicMesh/layerAdditionRemoval/layerAdditionRemoval.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
@@ -65,6 +65,10 @@ class layerAdditionRemoval
         //- Max thickness of extrusion layer.  Triggers layer addition
         mutable scalar maxLayerThickness_;
 
+        //- Switch to calculate thickness as volume/area
+        //  If false, thickness calculated from edges
+        const bool thicknessFromVolume_;
+
         //- Layer thickness from previous step
         //  Used to decide whether to add or remove layers
         mutable scalar oldLayerThickness_;
@@ -78,7 +82,7 @@ class layerAdditionRemoval
         //- Layer removal trigger time index
         mutable label triggerRemoval_;
 
-        //- Layer addition trigger  time index
+        //- Layer addition trigger time index
         mutable label triggerAddition_;
 
 
@@ -120,6 +124,7 @@ class layerAdditionRemoval
             //- Clear addressing
             void clearAddressing() const;
 
+
         // Helpers
 
             //- Optionally read old thickness
@@ -149,10 +154,11 @@ public:
         (
             const word& name,
             const label index,
-            const polyTopoChanger& mme,
+            const polyTopoChanger& ptc,
             const word& zoneName,
             const scalar minThickness,
-            const scalar maxThickness
+            const scalar maxThickness,
+            const Switch thicknessFromVolume = true
         );
 
         //- Construct from dictionary
@@ -161,7 +167,7 @@ public:
             const word& name,
             const dictionary& dict,
             const label index,
-            const polyTopoChanger& mme
+            const polyTopoChanger& ptc
         );
 
 

From b66b9637ae2a82a8b2e46189312bb1998299567c Mon Sep 17 00:00:00 2001
From: Sergio Ferraris <sergio@alex.opencfd.co.uk>
Date: Tue, 3 Sep 2013 11:33:05 +0100
Subject: [PATCH 20/20] STY: Deleting psiCombustionModels.C

---
 .../psiCombustionModel/psiCombustionModels.C  | 44 -------------------
 1 file changed, 44 deletions(-)
 delete mode 100644 src/combustionModels/psiCombustionModel/psiCombustionModel/psiCombustionModels.C

diff --git a/src/combustionModels/psiCombustionModel/psiCombustionModel/psiCombustionModels.C b/src/combustionModels/psiCombustionModel/psiCombustionModel/psiCombustionModels.C
deleted file mode 100644
index e8d96e9b63..0000000000
--- a/src/combustionModels/psiCombustionModel/psiCombustionModel/psiCombustionModels.C
+++ /dev/null
@@ -1,44 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 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 "makeCombustionTypes.H"
-#include "psiCombustionModel.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace combustionModels
-{
-    makeCombustionTypes
-    (
-        infinitelyFastChemistry,
-        psiCombustionModel
-    );
-}
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //