From 18d3b81960a7b344e1687c1095ab0178e165582d Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Mon, 13 May 2013 09:38:08 +0100
Subject: [PATCH 01/33] ENH: pressurTools function object - simplified usage

---
 .../utilities/pressureTools/pressureTools.C                | 7 +++----
 .../utilities/pressureTools/pressureTools.H                | 6 +-----
 2 files changed, 4 insertions(+), 9 deletions(-)

diff --git a/src/postProcessing/functionObjects/utilities/pressureTools/pressureTools.C b/src/postProcessing/functionObjects/utilities/pressureTools/pressureTools.C
index 607d52e503..181099c007 100644
--- a/src/postProcessing/functionObjects/utilities/pressureTools/pressureTools.C
+++ b/src/postProcessing/functionObjects/utilities/pressureTools/pressureTools.C
@@ -70,7 +70,7 @@ Foam::dimensionedScalar Foam::pressureTools::rhoScale
     }
     else
     {
-        return dimensionedScalar("rhoRef", dimDensity, rhoRef_);
+        return dimensionedScalar("rhoRef", dimDensity, rhoInf_);
     }
 }
 
@@ -100,7 +100,7 @@ Foam::tmp<Foam::volScalarField> Foam::pressureTools::rho
                         IOobject::NO_WRITE
                     ),
                     p.mesh(),
-                    dimensionedScalar("zero", dimDensity, rhoRef_)
+                    dimensionedScalar("zero", dimDensity, rhoInf_)
                 )
             );
     }
@@ -193,7 +193,6 @@ Foam::pressureTools::pressureTools
     pName_("p"),
     UName_("U"),
     rhoName_("rho"),
-    rhoRef_(1.0),
     calcTotal_(false),
     pRef_(0.0),
     calcCoeff_(false),
@@ -242,7 +241,7 @@ void Foam::pressureTools::read(const dictionary& dict)
 
         if (p.dimensions() != dimPressure)
         {
-            dict.lookup("rhoRef") >> rhoRef_;
+            dict.lookup("rhoRef") >> rhoInf_;
         } 
 
         dict.lookup("calcTotal") >> calcTotal_;
diff --git a/src/postProcessing/functionObjects/utilities/pressureTools/pressureTools.H b/src/postProcessing/functionObjects/utilities/pressureTools/pressureTools.H
index 11f490429d..82a44921ed 100644
--- a/src/postProcessing/functionObjects/utilities/pressureTools/pressureTools.H
+++ b/src/postProcessing/functionObjects/utilities/pressureTools/pressureTools.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -92,7 +92,6 @@ Description
     \table
         Property     | Description             | Required    | Default value
         type         | type name: pressureTools| yes         |
-        rhoRef       | Reference density for incompressible cases | no | 1
         calcTotal    | Calculate total coefficient | yes     |
         pRef         | Reference pressure for total pressure | no | 0.0
         calcCoeff    | Calculate pressure coefficient | yes  |
@@ -150,9 +149,6 @@ class pressureTools
         //- Name of density field, default is "rho"
         word rhoName_;
 
-        //- Reference density employed for incompressible cases
-        scalar rhoRef_;
-
 
         // Total pressure calculation
 

From cf085c7f12b3ee0648abe61249a905fe995e8532 Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Mon, 13 May 2013 09:38:30 +0100
Subject: [PATCH 02/33] STYLE: Corrected fatal error message

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

diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ConeNozzleInjection/ConeNozzleInjection.C b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ConeNozzleInjection/ConeNozzleInjection.C
index cebed76580..97c3c819a0 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ConeNozzleInjection/ConeNozzleInjection.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ConeNozzleInjection/ConeNozzleInjection.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
@@ -167,7 +167,8 @@ Foam::ConeNozzleInjection<CloudType>::ConeNozzleInjection
             "Foam::ConeNozzleInjection<CloudType>::ConeNozzleInjection"
             "("
                 "const dictionary&, "
-                "CloudType&"
+                "CloudType&, "
+                "const word&"
             ")"
         )<< "innerNozzleDiameter >= outerNozzleDiameter" << nl
          << exit(FatalError);
@@ -369,6 +370,7 @@ void Foam::ConeNozzleInjection<CloudType>::setPositionAndCell
                     "const scalar, "
                     "vector&, "
                     "label&, "
+                    "label&, "
                     "label&"
                 ")"
             )<< "Unknown injectionMethod type" << nl

From c08917f1963ed53974a560d9cb3ce2c945c855d7 Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Mon, 13 May 2013 09:38:54 +0100
Subject: [PATCH 03/33] STYLE: Code formatting

---
 .../functions/DataEntry/DataEntry/DataEntry.H | 104 +++++++++---------
 1 file changed, 51 insertions(+), 53 deletions(-)

diff --git a/src/OpenFOAM/primitives/functions/DataEntry/DataEntry/DataEntry.H b/src/OpenFOAM/primitives/functions/DataEntry/DataEntry/DataEntry.H
index 9cfdb66e25..a9b1686d58 100644
--- a/src/OpenFOAM/primitives/functions/DataEntry/DataEntry/DataEntry.H
+++ b/src/OpenFOAM/primitives/functions/DataEntry/DataEntry/DataEntry.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -142,62 +142,60 @@ public:
             virtual void convertTimeBase(const Time& t);
 
 
-public:
+        // Evaluation
+
+            //- Return value as a function of (scalar) independent variable
+            virtual Type value(const scalar x) const;
+
+            //- Return value as a function of (scalar) independent variable
+            virtual tmp<Field<Type> > value(const scalarField& x) const;
+
+            //- Integrate between two (scalar) values
+            virtual Type integrate(const scalar x1, const scalar x2) const;
+
+            //- Integrate between two (scalar) values
+            virtual tmp<Field<Type> > integrate
+            (
+                const scalarField& x1,
+                const scalarField& x2
+            ) const;
+
+            //- Return dimensioned type
+            virtual dimensioned<Type> dimValue(const scalar x) const;
+
+            //- Return dimensioned type as a function of (scalar)
+            virtual tmp<Field<dimensioned<Type> > > dimValue
+            (
+                const scalarField& x
+            ) const;
+
+            //- Integrate between two scalars and return a dimensioned type
+            virtual dimensioned<Type> dimIntegrate
+            (
+                const scalar x1,
+                const scalar x2
+            ) const;
+
+            //- Integrate between two scalar fields and return a field of
+            //  dimensioned type
+            virtual tmp<Field<dimensioned<Type> > > dimIntegrate
+            (
+                const scalarField& x1,
+                const scalarField& x2
+            ) const;
 
 
-    // Evaluation
+        // I/O
 
-        //- Return value as a function of (scalar) independent variable
-        virtual Type value(const scalar x) const;
+            //- Ostream Operator
+            friend Ostream& operator<< <Type>
+            (
+                Ostream& os,
+                const DataEntry<Type>& de
+            );
 
-        //- Return value as a function of (scalar) independent variable
-        virtual tmp<Field<Type> > value(const scalarField& x) const;
-
-        //- Integrate between two (scalar) values
-        virtual Type integrate(const scalar x1, const scalar x2) const;
-
-        //- Integrate between two (scalar) values
-        virtual tmp<Field<Type> > integrate
-        (
-            const scalarField& x1,
-            const scalarField& x2
-        ) const;
-
-        //- Return dimensioned type
-        virtual dimensioned<Type> dimValue(const scalar x) const;
-
-        //- Return dimensioned type as a function of (scalar)
-        virtual tmp<Field<dimensioned<Type> > > dimValue
-        (
-            const scalarField& x
-        ) const;
-
-        //- Integrate between two scalars and returns a dimensioned type
-        virtual dimensioned<Type> dimIntegrate
-        (
-            const scalar x1,
-            const scalar x2
-        ) const;
-
-        //- Integrate between two scalars and returns list of dimensioned type
-        virtual tmp<Field<dimensioned<Type> > > dimIntegrate
-        (
-            const scalarField& x1,
-            const scalarField& x2
-        ) const;
-
-
-    // I/O
-
-        //- Ostream Operator
-        friend Ostream& operator<< <Type>
-        (
-            Ostream& os,
-            const DataEntry<Type>& de
-        );
-
-        //- Write in dictionary format
-        virtual void writeData(Ostream& os) const;
+            //- Write in dictionary format
+            virtual void writeData(Ostream& os) const;
 };
 
 

From c6d61b01799bd5672227eb28b2cea39d8d302b04 Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Mon, 13 May 2013 12:10:27 +0100
Subject: [PATCH 04/33] ENH: Forces function object - refactored divRhoEff
 function

---
 .../functionObjects/forces/forces/forces.C    | 49 +++++++------------
 1 file changed, 17 insertions(+), 32 deletions(-)

diff --git a/src/postProcessing/functionObjects/forces/forces/forces.C b/src/postProcessing/functionObjects/forces/forces/forces.C
index c159b652e8..dd164d99a4 100644
--- a/src/postProcessing/functionObjects/forces/forces/forces.C
+++ b/src/postProcessing/functionObjects/forces/forces/forces.C
@@ -28,16 +28,12 @@ License
 #include "dictionary.H"
 #include "Time.H"
 #include "wordReList.H"
-
-#include "incompressible/singlePhaseTransportModel/singlePhaseTransportModel.H"
-#include "incompressible/RAS/RASModel/RASModel.H"
-#include "incompressible/LES/LESModel/LESModel.H"
-
-#include "fluidThermo.H"
-#include "compressible/RAS/RASModel/RASModel.H"
-#include "compressible/LES/LESModel/LESModel.H"
-
+#include "fvcGrad.H"
 #include "porosityModel.H"
+#include "fluidThermo.H"
+#include "incompressible/turbulenceModel/turbulenceModel.H"
+#include "compressible/turbulenceModel/turbulenceModel.H"
+#include "incompressible/singlePhaseTransportModel/singlePhaseTransportModel.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -70,38 +66,27 @@ void Foam::forces::writeFileHeader(const label i)
 
 Foam::tmp<Foam::volSymmTensorField> Foam::forces::devRhoReff() const
 {
-    if (obr_.foundObject<compressible::RASModel>("RASProperties"))
-    {
-        const compressible::RASModel& ras
-            = obr_.lookupObject<compressible::RASModel>("RASProperties");
+    typedef compressible::turbulenceModel cmpTurbModel;
+    typedef incompressible::turbulenceModel icoTurbModel;
 
-        return ras.devRhoReff();
-    }
-    else if (obr_.foundObject<incompressible::RASModel>("RASProperties"))
+    if (obr_.foundObject<cmpTurbModel>(cmpTurbModel::typeName))
     {
-        const incompressible::RASModel& ras
-            = obr_.lookupObject<incompressible::RASModel>("RASProperties");
+        const cmpTurbModel& turb =
+            obr_.lookupObject<cmpTurbModel>(cmpTurbModel::typeName);
 
-        return rho()*ras.devReff();
+        return turb.devRhoReff();
     }
-    else if (obr_.foundObject<compressible::LESModel>("LESProperties"))
+    else if (obr_.foundObject<icoTurbModel>(icoTurbModel::typeName))
     {
-        const compressible::LESModel& les =
-        obr_.lookupObject<compressible::LESModel>("LESProperties");
+        const incompressible::turbulenceModel& turb =
+            obr_.lookupObject<icoTurbModel>(icoTurbModel::typeName);
 
-        return les.devRhoReff();
+        return rho()*turb.devReff();
     }
-    else if (obr_.foundObject<incompressible::LESModel>("LESProperties"))
-    {
-        const incompressible::LESModel& les
-            = obr_.lookupObject<incompressible::LESModel>("LESProperties");
-
-        return rho()*les.devReff();
-    }
-    else if (obr_.foundObject<fluidThermo>("thermophysicalProperties"))
+    else if (obr_.foundObject<fluidThermo>(fluidThermo::typeName))
     {
         const fluidThermo& thermo =
-             obr_.lookupObject<fluidThermo>("thermophysicalProperties");
+            obr_.lookupObject<fluidThermo>(fluidThermo::typeName);
 
         const volVectorField& U = obr_.lookupObject<volVectorField>(UName_);
 

From a41c30ff0435179d39dc7420d771330233a4d37b Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Mon, 13 May 2013 15:19:31 +0100
Subject: [PATCH 05/33] ENH: Added ranomise to lookup table injection models

---
 .../KinematicLookupTableInjection.C               | 15 +++++++++++++--
 .../KinematicLookupTableInjection.H               |  5 ++++-
 .../ReactingLookupTableInjection.C                | 15 +++++++++++++--
 .../ReactingLookupTableInjection.H                |  5 ++++-
 .../ReactingMultiphaseLookupTableInjection.C      | 15 +++++++++++++--
 .../ReactingMultiphaseLookupTableInjection.H      |  5 ++++-
 .../ThermoLookupTableInjection.C                  | 15 +++++++++++++--
 .../ThermoLookupTableInjection.H                  |  5 ++++-
 8 files changed, 68 insertions(+), 12 deletions(-)

diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/KinematicLookupTableInjection/KinematicLookupTableInjection.C b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/KinematicLookupTableInjection/KinematicLookupTableInjection.C
index c40a95fbac..293fd7d740 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/KinematicLookupTableInjection/KinematicLookupTableInjection.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/KinematicLookupTableInjection/KinematicLookupTableInjection.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
@@ -43,6 +43,7 @@ Foam::KinematicLookupTableInjection<CloudType>::KinematicLookupTableInjection
     (
         readScalar(this->coeffDict().lookup("parcelsPerSecond"))
     ),
+    randomise_(readBool(this->coeffDict().lookup("randomise"))),
     injectors_
     (
         IOobject
@@ -87,6 +88,7 @@ Foam::KinematicLookupTableInjection<CloudType>::KinematicLookupTableInjection
     inputFileName_(im.inputFileName_),
     duration_(im.duration_),
     parcelsPerSecond_(im.parcelsPerSecond_),
+    randomise_(im.randomise_),
     injectors_(im.injectors_),
     injectorCells_(im.injectorCells_),
     injectorTetFaces_(im.injectorTetFaces_),
@@ -177,7 +179,16 @@ void Foam::KinematicLookupTableInjection<CloudType>::setPositionAndCell
     label& tetPtI
 )
 {
-    label injectorI = parcelI*injectorCells_.size()/nParcels;
+    label injectorI = 0;
+    if (randomise_)
+    {
+        cachedRandom& rnd = this->owner().rndGen();
+        injectorI = rnd.position<label>(0, injectorCells_.size() - 1);
+    }
+    else
+    {
+        injectorI = parcelI*injectorCells_.size()/nParcels;
+    }
 
     position = injectors_[injectorI].x();
     cellOwner = injectorCells_[injectorI];
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/KinematicLookupTableInjection/KinematicLookupTableInjection.H b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/KinematicLookupTableInjection/KinematicLookupTableInjection.H
index d156f4793a..f17410bc1d 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/KinematicLookupTableInjection/KinematicLookupTableInjection.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/KinematicLookupTableInjection/KinematicLookupTableInjection.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -78,6 +78,9 @@ class KinematicLookupTableInjection
         //- Number of parcels per injector - common to all injection sources
         const scalar parcelsPerSecond_;
 
+        //- Flag to indicate to randomise injection positions
+        bool randomise_;
+
         //- List of injectors
         kinematicParcelInjectionDataIOList injectors_;
 
diff --git a/src/lagrangian/intermediate/submodels/Reacting/InjectionModel/ReactingLookupTableInjection/ReactingLookupTableInjection.C b/src/lagrangian/intermediate/submodels/Reacting/InjectionModel/ReactingLookupTableInjection/ReactingLookupTableInjection.C
index 711b3cec7a..1259da7e37 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/InjectionModel/ReactingLookupTableInjection/ReactingLookupTableInjection.C
+++ b/src/lagrangian/intermediate/submodels/Reacting/InjectionModel/ReactingLookupTableInjection/ReactingLookupTableInjection.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
@@ -42,6 +42,7 @@ Foam::ReactingLookupTableInjection<CloudType>::ReactingLookupTableInjection
     (
         readScalar(this->coeffDict().lookup("parcelsPerSecond"))
     ),
+    randomise_(readBool(this->coeffDict().lookup("randomise"))),
     injectors_
     (
         IOobject
@@ -86,6 +87,7 @@ Foam::ReactingLookupTableInjection<CloudType>::ReactingLookupTableInjection
     inputFileName_(im.inputFileName_),
     duration_(im.duration_),
     parcelsPerSecond_(im.parcelsPerSecond_),
+    randomise_(im.randomise_),
     injectors_(im.injectors_),
     injectorCells_(im.injectorCells_),
     injectorTetFaces_(im.injectorTetFaces_),
@@ -176,7 +178,16 @@ void Foam::ReactingLookupTableInjection<CloudType>::setPositionAndCell
     label& tetPtI
 )
 {
-    label injectorI = parcelI*injectorCells_.size()/nParcels;
+    label injectorI = 0;
+    if (randomise_)
+    {
+        cachedRandom& rnd = this->owner().rndGen();
+        injectorI = rnd.position<label>(0, injectorCells_.size() - 1);
+    }
+    else
+    {
+        injectorI = parcelI*injectorCells_.size()/nParcels;
+    }
 
     position = injectors_[injectorI].x();
     cellOwner = injectorCells_[injectorI];
diff --git a/src/lagrangian/intermediate/submodels/Reacting/InjectionModel/ReactingLookupTableInjection/ReactingLookupTableInjection.H b/src/lagrangian/intermediate/submodels/Reacting/InjectionModel/ReactingLookupTableInjection/ReactingLookupTableInjection.H
index ef062fd72b..654109c509 100644
--- a/src/lagrangian/intermediate/submodels/Reacting/InjectionModel/ReactingLookupTableInjection/ReactingLookupTableInjection.H
+++ b/src/lagrangian/intermediate/submodels/Reacting/InjectionModel/ReactingLookupTableInjection/ReactingLookupTableInjection.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -81,6 +81,9 @@ class ReactingLookupTableInjection
         //- Number of parcels per injector - common to all injection sources
         const scalar parcelsPerSecond_;
 
+        //- Flag to indicate to randomise injection positions
+        bool randomise_;
+
         //- List of injectors
         reactingParcelInjectionDataIOList injectors_;
 
diff --git a/src/lagrangian/intermediate/submodels/ReactingMultiphase/InjectionModel/ReactingMultiphaseLookupTableInjection/ReactingMultiphaseLookupTableInjection.C b/src/lagrangian/intermediate/submodels/ReactingMultiphase/InjectionModel/ReactingMultiphaseLookupTableInjection/ReactingMultiphaseLookupTableInjection.C
index bd630d879d..179b846caf 100644
--- a/src/lagrangian/intermediate/submodels/ReactingMultiphase/InjectionModel/ReactingMultiphaseLookupTableInjection/ReactingMultiphaseLookupTableInjection.C
+++ b/src/lagrangian/intermediate/submodels/ReactingMultiphase/InjectionModel/ReactingMultiphaseLookupTableInjection/ReactingMultiphaseLookupTableInjection.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
@@ -43,6 +43,7 @@ ReactingMultiphaseLookupTableInjection
     (
         readScalar(this->coeffDict().lookup("parcelsPerSecond"))
     ),
+    randomise_(readBool(this->coeffDict().lookup("randomise"))),
     injectors_
     (
         IOobject
@@ -88,6 +89,7 @@ ReactingMultiphaseLookupTableInjection
     inputFileName_(im.inputFileName_),
     duration_(im.duration_),
     parcelsPerSecond_(im.parcelsPerSecond_),
+    randomise_(im.randomise_),
     injectors_(im.injectors_),
     injectorCells_(im.injectorCells_),
     injectorTetFaces_(im.injectorTetFaces_),
@@ -182,7 +184,16 @@ void Foam::ReactingMultiphaseLookupTableInjection<CloudType>::setPositionAndCell
     label& tetPtI
 )
 {
-    label injectorI = parcelI*injectorCells_.size()/nParcels;
+    label injectorI = 0;
+    if (randomise_)
+    {
+        cachedRandom& rnd = this->owner().rndGen();
+        injectorI = rnd.position<label>(0, injectorCells_.size() - 1);
+    }
+    else
+    {
+        injectorI = parcelI*injectorCells_.size()/nParcels;
+    }
 
     position = injectors_[injectorI].x();
     cellOwner = injectorCells_[injectorI];
diff --git a/src/lagrangian/intermediate/submodels/ReactingMultiphase/InjectionModel/ReactingMultiphaseLookupTableInjection/ReactingMultiphaseLookupTableInjection.H b/src/lagrangian/intermediate/submodels/ReactingMultiphase/InjectionModel/ReactingMultiphaseLookupTableInjection/ReactingMultiphaseLookupTableInjection.H
index 7e6c328558..81f41d42de 100644
--- a/src/lagrangian/intermediate/submodels/ReactingMultiphase/InjectionModel/ReactingMultiphaseLookupTableInjection/ReactingMultiphaseLookupTableInjection.H
+++ b/src/lagrangian/intermediate/submodels/ReactingMultiphase/InjectionModel/ReactingMultiphaseLookupTableInjection/ReactingMultiphaseLookupTableInjection.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -84,6 +84,9 @@ class ReactingMultiphaseLookupTableInjection
         //- Number of parcels per injector - common to all injection sources
         const scalar parcelsPerSecond_;
 
+        //- Flag to indicate to randomise injection positions
+        bool randomise_;
+
         //- List of injectors
         reactingMultiphaseParcelInjectionDataIOList injectors_;
 
diff --git a/src/lagrangian/intermediate/submodels/Thermodynamic/InjectionModel/ThermoLookupTableInjection/ThermoLookupTableInjection.C b/src/lagrangian/intermediate/submodels/Thermodynamic/InjectionModel/ThermoLookupTableInjection/ThermoLookupTableInjection.C
index fdc0a3b4a3..682fc9c3ea 100644
--- a/src/lagrangian/intermediate/submodels/Thermodynamic/InjectionModel/ThermoLookupTableInjection/ThermoLookupTableInjection.C
+++ b/src/lagrangian/intermediate/submodels/Thermodynamic/InjectionModel/ThermoLookupTableInjection/ThermoLookupTableInjection.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
@@ -43,6 +43,7 @@ Foam::ThermoLookupTableInjection<CloudType>::ThermoLookupTableInjection
     (
         readScalar(this->coeffDict().lookup("parcelsPerSecond"))
     ),
+    randomise_(readBool(this->coeffDict().lookup("randomise"))),
     injectors_
     (
         IOobject
@@ -87,6 +88,7 @@ Foam::ThermoLookupTableInjection<CloudType>::ThermoLookupTableInjection
     inputFileName_(im.inputFileName_),
     duration_(im.duration_),
     parcelsPerSecond_(im.parcelsPerSecond_),
+    randomise_(im.randomise_),
     injectors_(im.injectors_),
     injectorCells_(im.injectorCells_),
     injectorTetFaces_(im.injectorTetFaces_),
@@ -177,7 +179,16 @@ void Foam::ThermoLookupTableInjection<CloudType>::setPositionAndCell
     label& tetPtI
 )
 {
-    label injectorI = parcelI*injectorCells_.size()/nParcels;
+    label injectorI = 0;
+    if (randomise_)
+    {
+        cachedRandom& rnd = this->owner().rndGen();
+        injectorI = rnd.position<label>(0, injectorCells_.size() - 1);
+    }
+    else
+    {
+        injectorI = parcelI*injectorCells_.size()/nParcels;
+    }
 
     position = injectors_[injectorI].x();
     cellOwner = injectorCells_[injectorI];
diff --git a/src/lagrangian/intermediate/submodels/Thermodynamic/InjectionModel/ThermoLookupTableInjection/ThermoLookupTableInjection.H b/src/lagrangian/intermediate/submodels/Thermodynamic/InjectionModel/ThermoLookupTableInjection/ThermoLookupTableInjection.H
index efb83c3614..42c48f290d 100644
--- a/src/lagrangian/intermediate/submodels/Thermodynamic/InjectionModel/ThermoLookupTableInjection/ThermoLookupTableInjection.H
+++ b/src/lagrangian/intermediate/submodels/Thermodynamic/InjectionModel/ThermoLookupTableInjection/ThermoLookupTableInjection.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -80,6 +80,9 @@ class ThermoLookupTableInjection
         //- Number of parcels per injector - common to all injection sources
         const scalar parcelsPerSecond_;
 
+        //- Flag to indicate to randomise injection positions
+        bool randomise_;
+
         //- List of injectors
         kinematicParcelInjectionDataIOList injectors_;
 

From 1ecf2b471a2ee77054bfed1933e105f493c6c6ff Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Mon, 13 May 2013 15:32:04 +0100
Subject: [PATCH 06/33] ENH: updated initialisation for fieldAverage function
 object - mantis #841

---
 .../fieldAverage/fieldAverage/fieldAverage.C      | 15 +++++++++++----
 .../fieldAverage/fieldAverage/fieldAverage.H      |  3 +++
 2 files changed, 14 insertions(+), 4 deletions(-)

diff --git a/src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverage.C b/src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverage.C
index d714747310..f19b7f69d1 100644
--- a/src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverage.C
+++ b/src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverage.C
@@ -311,6 +311,7 @@ Foam::fieldAverage::fieldAverage
     prevTimeIndex_(-1),
     resetOnRestart_(false),
     resetOnOutput_(false),
+    initialised_(false),
     faItems_(),
     meanScalarFields_(),
     meanVectorFields_(),
@@ -361,11 +362,7 @@ void Foam::fieldAverage::read(const dictionary& dict)
         dict.readIfPresent("resetOnOutput", resetOnOutput_);
         dict.lookup("fields") >> faItems_;
 
-        initialize();
         readAveragingProperties();
-
-        // ensure first averaging works unconditionally
-        prevTimeIndex_ = -1;
     }
 }
 
@@ -374,6 +371,16 @@ void Foam::fieldAverage::execute()
 {
     if (active_)
     {
+        if (!initialised_)
+        {
+            initialize();
+
+            // ensure first averaging works unconditionally
+            prevTimeIndex_ = -1;
+
+            initialised_ = true;
+        }
+
         calcAverages();
     }
 }
diff --git a/src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverage.H b/src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverage.H
index f8b689757d..4a6d75e06b 100644
--- a/src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverage.H
+++ b/src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverage.H
@@ -171,6 +171,9 @@ protected:
         //- Reset the averaging process on output flag
         Switch resetOnOutput_;
 
+        //- Initialised flag
+        bool initialised_;
+
         //- List of field average items, describing what averages to be
         //  calculated and output
         List<fieldAverageItem> faItems_;

From 78b9b93f77b4e71580ee13836673be3c59981334 Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Mon, 13 May 2013 17:29:38 +0100
Subject: [PATCH 07/33] BUG: Corrected particle lift force calculation - mantis
 #822

---
 .../Kinematic/ParticleForces/Lift/LiftForce/LiftForce.C       | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Lift/LiftForce/LiftForce.C b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Lift/LiftForce/LiftForce.C
index 290d2f3c3d..a4faceeec7 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Lift/LiftForce/LiftForce.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Lift/LiftForce/LiftForce.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -141,7 +141,7 @@ Foam::forceSuSp Foam::LiftForce<CloudType>::calcCoupled
 
     scalar Cl = this->Cl(p, curlUc, Re, muc);
 
-    value.Su() = mass/p.rho()*p.d()/2.0*p.rhoc()*Cl*((p.Uc() - p.U())^curlUc);
+    value.Su() = mass/p.rho()*p.rhoc()*Cl*((p.Uc() - p.U())^curlUc);
 
     return value;
 }

From fc73c4629ae9852a0dbd284f83fc0c4f4bfc58fb Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Mon, 13 May 2013 17:50:10 +0100
Subject: [PATCH 08/33] STYLE: minor code formatting

---
 .../cfdTools/general/adjustPhi/adjustPhi.C       | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/src/finiteVolume/cfdTools/general/adjustPhi/adjustPhi.C b/src/finiteVolume/cfdTools/general/adjustPhi/adjustPhi.C
index 206e7c4e95..27622aeeb7 100644
--- a/src/finiteVolume/cfdTools/general/adjustPhi/adjustPhi.C
+++ b/src/finiteVolume/cfdTools/general/adjustPhi/adjustPhi.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
@@ -56,11 +56,7 @@ bool Foam::adjustPhi
 
             if (!isA<processorFvsPatchScalarField>(phip))
             {
-                if
-                (
-                    Up.fixesValue()
-                 && !isA<inletOutletFvPatchVectorField>(Up)
-                )
+                if (Up.fixesValue() && !isA<inletOutletFvPatchVectorField>(Up))
                 {
                     forAll(phip, i)
                     {
@@ -113,8 +109,12 @@ bool Foam::adjustPhi
         {
             FatalErrorIn
             (
-                "adjustPhi(surfaceScalarField& phi, const volVectorField& U,"
-                "const volScalarField& p"
+                "adjustPhi"
+                "("
+                    "surfaceScalarField&, "
+                    "const volVectorField&,"
+                    "volScalarField&"
+                ")"
             )   << "Continuity error cannot be removed by adjusting the"
                    " outflow.\nPlease check the velocity boundary conditions"
                    " and/or run potentialFoam to initialise the outflow." << nl

From d99266c2e6ba6106e53fc95699fd6d59a0cac98c Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Mon, 13 May 2013 17:51:23 +0100
Subject: [PATCH 09/33] ENH: Refactored film shear stress model -> film
 turbulence model

---
 src/regionModels/surfaceFilmModels/Make/files |   4 +
 .../kinematicSingleLayer.C                    |  21 ++-
 .../kinematicSingleLayer.H                    |   9 +-
 .../kinematicSingleLayerI.H                   |   8 +-
 .../filmTurbulenceModel/filmTurbulenceModel.C |  73 ++++++++
 .../filmTurbulenceModel/filmTurbulenceModel.H | 148 ++++++++++++++++
 .../filmTurbulenceModelNew.C                  |  77 +++++++++
 .../filmTurbulenceModel/laminar/laminar.C     | 160 ++++++++++++++++++
 .../filmTurbulenceModel/laminar/laminar.H     | 118 +++++++++++++
 9 files changed, 607 insertions(+), 11 deletions(-)
 create mode 100644 src/regionModels/surfaceFilmModels/submodels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.C
 create mode 100644 src/regionModels/surfaceFilmModels/submodels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.H
 create mode 100644 src/regionModels/surfaceFilmModels/submodels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModelNew.C
 create mode 100644 src/regionModels/surfaceFilmModels/submodels/kinematic/filmTurbulenceModel/laminar/laminar.C
 create mode 100644 src/regionModels/surfaceFilmModels/submodels/kinematic/filmTurbulenceModel/laminar/laminar.H

diff --git a/src/regionModels/surfaceFilmModels/Make/files b/src/regionModels/surfaceFilmModels/Make/files
index 1151e2f7fd..f14a789338 100644
--- a/src/regionModels/surfaceFilmModels/Make/files
+++ b/src/regionModels/surfaceFilmModels/Make/files
@@ -24,6 +24,10 @@ $(KINEMATICMODELS)/injectionModel/drippingInjection/drippingInjection.C
 $(KINEMATICMODELS)/injectionModel/removeInjection/removeInjection.C
 $(KINEMATICMODELS)/injectionModel/curvatureSeparation/curvatureSeparation.C
 
+$(KINEMATICMODELS)/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.C
+$(KINEMATICMODELS)/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModelNew.C
+$(KINEMATICMODELS)/filmTurbulenceModel/laminar/laminar.C
+
 THERMOMODELS=submodels/thermo
 $(THERMOMODELS)/phaseChangeModel/phaseChangeModel/phaseChangeModel.C
 $(THERMOMODELS)/phaseChangeModel/phaseChangeModel/phaseChangeModelNew.C
diff --git a/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C b/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C
index fecefcd2e3..00a3ac39af 100644
--- a/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C
+++ b/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C
@@ -163,7 +163,7 @@ tmp<volScalarField> kinematicSingleLayer::pu()
         (
             IOobject
             (
-                "pu",
+                typeName + ":pu",
                 time_.timeName(),
                 regionMesh(),
                 IOobject::NO_READ,
@@ -185,7 +185,7 @@ tmp<volScalarField> kinematicSingleLayer::pp()
         (
             IOobject
             (
-                "pp",
+                typeName + ":pp",
                 time_.timeName(),
                 regionMesh(),
                 IOobject::NO_READ,
@@ -216,6 +216,8 @@ void kinematicSingleLayer::updateSubmodels()
     // Update source fields
     const dimensionedScalar deltaT = time().deltaT();
     rhoSp_ += cloudMassTrans_/magSf()/deltaT;
+
+    turbulence_->correct();
 }
 
 
@@ -282,9 +284,7 @@ void kinematicSingleLayer::updateSurfaceVelocities()
     Uw_ -= nHat()*(Uw_ & nHat());
     Uw_.correctBoundaryConditions();
 
-    // apply quadratic profile to surface velocity (scale by sqrt(2))
-    Us_ = 1.414*U_;
-    Us_.correctBoundaryConditions();
+    Us_ = turbulence_->Us();
 }
 
 
@@ -309,6 +309,7 @@ tmp<Foam::fvVectorMatrix> kinematicSingleLayer::solveMomentum
 //      - fvm::SuSp(rhoSp_, U_)
       - rhoSp_*U_
       + forces_.correct(U_)
+      + turbulence_->Su(U_)
     );
 
     fvVectorMatrix& UEqn = tUEqn();
@@ -779,6 +780,8 @@ kinematicSingleLayer::kinematicSingleLayer
 
     injection_(*this, coeffs_),
 
+    turbulence_(filmTurbulenceModel::New(*this, coeffs_)),
+
     forces_(*this, coeffs_),
 
     addedMassTotal_(0.0)
@@ -1020,7 +1023,7 @@ tmp<volScalarField> kinematicSingleLayer::primaryMassTrans() const
         (
             IOobject
             (
-                "kinematicSingleLayer::primaryMassTrans",
+                typeName + ":primaryMassTrans",
                 time().timeName(),
                 primaryMesh(),
                 IOobject::NO_READ,
@@ -1075,7 +1078,7 @@ tmp<DimensionedField<scalar, volMesh> > kinematicSingleLayer::Srho() const
         (
             IOobject
             (
-                "kinematicSingleLayer::Srho",
+                typeName + ":Srho",
                 time().timeName(),
                 primaryMesh(),
                 IOobject::NO_READ,
@@ -1100,7 +1103,7 @@ tmp<DimensionedField<scalar, volMesh> > kinematicSingleLayer::Srho
         (
             IOobject
             (
-                "kinematicSingleLayer::Srho(" + Foam::name(i) + ")",
+                typeName + ":Srho(" + Foam::name(i) + ")",
                 time().timeName(),
                 primaryMesh(),
                 IOobject::NO_READ,
@@ -1122,7 +1125,7 @@ tmp<DimensionedField<scalar, volMesh> > kinematicSingleLayer::Sh() const
         (
             IOobject
             (
-                "kinematicSingleLayer::Sh",
+                typeName + ":Sh",
                 time().timeName(),
                 primaryMesh(),
                 IOobject::NO_READ,
diff --git a/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.H b/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.H
index 079999d0cd..2c9e1020a3 100644
--- a/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.H
+++ b/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -43,6 +43,7 @@ SourceFiles
 
 #include "injectionModelList.H"
 #include "forceList.H"
+#include "filmTurbulenceModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -197,6 +198,9 @@ protected:
             //- Cloud injection
             injectionModelList injection_;
 
+            //- Turbulence model
+            autoPtr<filmTurbulenceModel> turbulence_;
+
             //- List of film forces
             forceList forces_;
 
@@ -444,6 +448,9 @@ public:
             //- Injection
             inline injectionModelList& injection();
 
+            //- Turbulence
+            inline const filmTurbulenceModel& turbulence() const;
+
 
         // Helper functions
 
diff --git a/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayerI.H b/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayerI.H
index 6a77b614bc..15f6398ad8 100644
--- a/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayerI.H
+++ b/src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayerI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -169,6 +169,12 @@ inline injectionModelList& kinematicSingleLayer::injection()
 }
 
 
+inline const filmTurbulenceModel& kinematicSingleLayer::turbulence() const
+{
+    return turbulence_();
+}
+
+
 inline tmp<volScalarField> kinematicSingleLayer::mass() const
 {
     return rho_*delta_*magSf();
diff --git a/src/regionModels/surfaceFilmModels/submodels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.C b/src/regionModels/surfaceFilmModels/submodels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.C
new file mode 100644
index 0000000000..a77307e881
--- /dev/null
+++ b/src/regionModels/surfaceFilmModels/submodels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.C
@@ -0,0 +1,73 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "filmTurbulenceModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(filmTurbulenceModel, 0);
+defineRunTimeSelectionTable(filmTurbulenceModel, dictionary);
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+filmTurbulenceModel::filmTurbulenceModel(const surfaceFilmModel& owner)
+:
+    subModelBase(owner)
+{}
+
+
+filmTurbulenceModel::filmTurbulenceModel
+(
+    const word& type,
+    const surfaceFilmModel& owner,
+    const dictionary& dict
+)
+:
+    subModelBase(type, owner, dict)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+filmTurbulenceModel::~filmTurbulenceModel()
+{}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/regionModels/surfaceFilmModels/submodels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.H b/src/regionModels/surfaceFilmModels/submodels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.H
new file mode 100644
index 0000000000..571bf83e8d
--- /dev/null
+++ b/src/regionModels/surfaceFilmModels/submodels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.H
@@ -0,0 +1,148 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::regionModels::surfaceFilmModels::filmTurbulenceModel
+
+Description
+    Base class for film turbulence models
+
+SourceFiles
+    filmTurbulenceModel.C
+    filmTurbulenceModelNew.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef filmTurbulenceModel_H
+#define filmTurbulenceModel_H
+
+#include "subModelBase.H"
+#include "runTimeSelectionTables.H"
+#include "fvMatricesFwd.H"
+#include "volFieldsFwd.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class filmTurbulenceModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class filmTurbulenceModel
+:
+    public subModelBase
+{
+private:
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        filmTurbulenceModel(const filmTurbulenceModel&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const filmTurbulenceModel&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("filmTurbulenceModel");
+
+
+    // Declare runtime constructor selection table
+
+         declareRunTimeSelectionTable
+         (
+             autoPtr,
+             filmTurbulenceModel,
+             dictionary,
+             (
+                const surfaceFilmModel& owner,
+                const dictionary& dict
+             ),
+             (owner, dict)
+         );
+
+    // Constructors
+
+        //- Construct null
+        filmTurbulenceModel(const surfaceFilmModel& owner);
+
+        //- Construct from type name, dictionary and surface film model
+        filmTurbulenceModel
+        (
+            const word& type,
+            const surfaceFilmModel& owner,
+            const dictionary& dict
+        );
+
+
+    // Selectors
+
+        //- Return a reference to the selected injection model
+        static autoPtr<filmTurbulenceModel> New
+        (
+            const surfaceFilmModel& owner,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~filmTurbulenceModel();
+
+
+    // Member Functions
+
+        // Evolution
+
+            //- Return the film surface velocity
+            virtual tmp<volVectorField> Us() const = 0;
+
+            //- Return the film turbulence viscosity
+            virtual tmp<volScalarField> mut() const = 0;
+
+            //- Correct/update the model
+            virtual void correct() = 0;
+
+            //- Return the source for the film momentum equation
+            virtual tmp<fvVectorMatrix> Su(volVectorField& U) const = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/regionModels/surfaceFilmModels/submodels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModelNew.C b/src/regionModels/surfaceFilmModels/submodels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModelNew.C
new file mode 100644
index 0000000000..ba939ef798
--- /dev/null
+++ b/src/regionModels/surfaceFilmModels/submodels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModelNew.C
@@ -0,0 +1,77 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "filmTurbulenceModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
+
+autoPtr<filmTurbulenceModel> filmTurbulenceModel::New
+(
+    const surfaceFilmModel& model,
+    const dictionary& dict
+)
+{
+    const word modelType(dict.lookup("turbulence"));
+
+    Info<< "        " << modelType << endl;
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(modelType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "filmTurbulenceModel::New"
+            "("
+                "const surfaceFilmModel&, "
+                "const dictionary&"
+            ")"
+        )   << "Unknown filmTurbulenceModel type " << modelType
+            << nl << nl << "Valid filmTurbulenceModel types are:" << nl
+            << dictionaryConstructorTablePtr_->toc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<filmTurbulenceModel>(cstrIter()(model, dict));
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/regionModels/surfaceFilmModels/submodels/kinematic/filmTurbulenceModel/laminar/laminar.C b/src/regionModels/surfaceFilmModels/submodels/kinematic/filmTurbulenceModel/laminar/laminar.C
new file mode 100644
index 0000000000..108109702b
--- /dev/null
+++ b/src/regionModels/surfaceFilmModels/submodels/kinematic/filmTurbulenceModel/laminar/laminar.C
@@ -0,0 +1,160 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "laminar.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvMesh.H"
+#include "fvMatrices.H"
+#include "Time.H"
+#include "volFields.H"
+#include "fvmSup.H"
+#include "kinematicSingleLayer.H"
+#include "zeroGradientFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(laminar, 0);
+addToRunTimeSelectionTable(filmTurbulenceModel, laminar, dictionary);
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+laminar::laminar
+(
+    const surfaceFilmModel& owner,
+    const dictionary& dict
+)
+:
+    filmTurbulenceModel(type(), owner, dict),
+    Cf_(readScalar(coeffs_.lookup("Cf")))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+laminar::~laminar()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+tmp<volVectorField> laminar::Us() const
+{
+    tmp<volVectorField> tUs
+    (
+        new volVectorField
+        (
+            IOobject
+            (
+                typeName + ".Us",
+                owner_.regionMesh().time().timeName(),
+                owner_.regionMesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            owner_.regionMesh(),
+            dimensionedVector("zero", dimVelocity, vector::zero),
+            zeroGradientFvPatchVectorField::typeName
+        )
+    );
+
+    // apply quadratic profile
+    tUs() = Foam::sqrt(2.0)*owner_.U();
+    tUs().correctBoundaryConditions();
+
+    return tUs;
+}
+
+
+tmp<volScalarField> laminar::mut() const
+{
+    return tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                typeName + ".mut",
+                owner_.regionMesh().time().timeName(),
+                owner_.regionMesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            owner_.regionMesh(),
+            dimensionedScalar("zero", dimMass/dimLength/dimTime, 0.0)
+        )
+    );
+}
+
+
+void laminar::correct()
+{
+    // do nothing
+}
+
+tmp<fvVectorMatrix> laminar::Su(volVectorField& U) const
+{
+    // local reference to film model
+    const kinematicSingleLayer& film =
+        static_cast<const kinematicSingleLayer&>(owner_);
+
+    // local references to film fields
+    const volScalarField& mu = film.mu();
+    const volVectorField& Uw = film.Uw();
+    const volVectorField& Us = film.Us();
+    const volScalarField& delta = film.delta();
+    const volVectorField& Up = film.UPrimary();
+    const volScalarField& rhop = film.rhoPrimary();
+
+    // employ simple coeff-based model
+    volScalarField Cs("Cs", Cf_*rhop*mag(Up - Us));
+
+    dimensionedScalar d0("SMALL", delta.dimensions(), SMALL);
+    volScalarField Cw("Cw", mu/(0.3333*(delta + d0)));
+    Cw.min(5000.0);
+
+    return
+    (
+       - fvm::Sp(Cs, U) + Cs*Us // surface contribution
+       - fvm::Sp(Cw, U) + Cw*Uw // wall contribution
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/regionModels/surfaceFilmModels/submodels/kinematic/filmTurbulenceModel/laminar/laminar.H b/src/regionModels/surfaceFilmModels/submodels/kinematic/filmTurbulenceModel/laminar/laminar.H
new file mode 100644
index 0000000000..397030efff
--- /dev/null
+++ b/src/regionModels/surfaceFilmModels/submodels/kinematic/filmTurbulenceModel/laminar/laminar.H
@@ -0,0 +1,118 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::regionModels::surfaceFilmModels::laminar
+
+Description
+    Film laminar turbulence model.
+
+SourceFiles
+    laminar.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef laminar_H
+#define laminar_H
+
+#include "filmTurbulenceModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class laminar Declaration
+\*---------------------------------------------------------------------------*/
+
+class laminar
+:
+    public filmTurbulenceModel
+{
+private:
+
+    // Private Data
+
+        //- Surface roughness coefficient
+        scalar Cf_;
+
+
+    // Private member functions
+
+        //- Disallow default bitwise copy construct
+        laminar(const laminar&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const laminar&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("laminar");
+
+
+    // Constructors
+
+        //- Construct from surface film model
+        laminar(const surfaceFilmModel& owner, const dictionary& dict);
+
+
+    //- Destructor
+    virtual ~laminar();
+
+
+    // Member Functions
+
+        // Evolution
+
+            //- Return the film surface velocity
+            virtual tmp<volVectorField> Us() const;
+
+            //- Return the film turbulence viscosity
+            virtual tmp<volScalarField> mut() const;
+
+            //- Correct/update the model
+            virtual void correct();
+
+            //- Return the source for the film momentum equation
+            virtual tmp<fvVectorMatrix> Su(volVectorField& U) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

From 14934a9bfc52cc3b918f77fab9ed2178bcfbf2e3 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 14 May 2013 16:45:59 +0100
Subject: [PATCH 10/33] BUG: faceAreaPairGAMGAgglomeration: corrected face area
 -> weight conversion

---
 .../faceAreaPairGAMGAgglomeration.C                          | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/src/finiteVolume/fvMatrices/solvers/GAMGSymSolver/GAMGAgglomerations/faceAreaPairGAMGAgglomeration/faceAreaPairGAMGAgglomeration.C b/src/finiteVolume/fvMatrices/solvers/GAMGSymSolver/GAMGAgglomerations/faceAreaPairGAMGAgglomeration/faceAreaPairGAMGAgglomeration.C
index 9bfcda020c..42a35f5c1b 100644
--- a/src/finiteVolume/fvMatrices/solvers/GAMGSymSolver/GAMGAgglomerations/faceAreaPairGAMGAgglomeration/faceAreaPairGAMGAgglomeration.C
+++ b/src/finiteVolume/fvMatrices/solvers/GAMGSymSolver/GAMGAgglomerations/faceAreaPairGAMGAgglomeration/faceAreaPairGAMGAgglomeration.C
@@ -90,8 +90,6 @@ Foam::faceAreaPairGAMGAgglomeration::faceAreaPairGAMGAgglomeration
 :
     pairGAMGAgglomeration(mesh, controlDict)
 {
-    vectorField n(faceAreas/mag(faceAreas));
-
     //agglomerate(mesh, sqrt(mag(faceAreas)));
     agglomerate
     (
@@ -100,7 +98,8 @@ Foam::faceAreaPairGAMGAgglomeration::faceAreaPairGAMGAgglomeration
         (
             cmptMultiply
             (
-                n,
+                faceAreas
+               /sqrt(mag(faceAreas)),
                 vector(1, 1.01, 1.02)
                 //vector::one
             )

From 6d3cfe9848aba235d36735bf848452354cc6a511 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 14 May 2013 20:42:54 +0100
Subject: [PATCH 11/33] ENH: GAMGAgglomeration: enable debug printing by
 default

---
 etc/controlDict | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/etc/controlDict b/etc/controlDict
index b9c41cb3a4..fba02069de 100644
--- a/etc/controlDict
+++ b/etc/controlDict
@@ -135,7 +135,7 @@ DebugSwitches
     FDIC                0;
     FaceCellWave        0;
     GAMG                0;
-    GAMGAgglomeration   0;
+    GAMGAgglomeration   1;
     GAMGInterface       0;
     GAMGInterfaceField  0;
     Gamma               0;

From 8009617b23ec4d96c5cc8652d73bef48a2ad9a66 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 14 May 2013 20:44:40 +0100
Subject: [PATCH 12/33] ENH: wordIOList: return columnsizes when printing table

---
 .../primitives/strings/word/wordIOList.C       | 18 +++++++++++++++---
 .../primitives/strings/word/wordIOList.H       |  3 ++-
 2 files changed, 17 insertions(+), 4 deletions(-)

diff --git a/src/OpenFOAM/primitives/strings/word/wordIOList.C b/src/OpenFOAM/primitives/strings/word/wordIOList.C
index ede987af53..58456cb0ac 100644
--- a/src/OpenFOAM/primitives/strings/word/wordIOList.C
+++ b/src/OpenFOAM/primitives/strings/word/wordIOList.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -38,12 +38,17 @@ namespace Foam
 }
 
 
-void Foam::printTable(const List<wordList>& wll, Ostream& os)
+void Foam::printTable
+(
+    const List<wordList>& wll,
+    List<string::size_type>& columnWidth,
+    Ostream& os
+)
 {
     if (!wll.size()) return;
 
     // Find the maximum word length for each column
-    List<string::size_type> columnWidth(wll[0].size(), string::size_type(0));
+    columnWidth.setSize(wll[0].size(), string::size_type(0));
     forAll(columnWidth, j)
     {
         forAll(wll, i)
@@ -75,4 +80,11 @@ void Foam::printTable(const List<wordList>& wll, Ostream& os)
 }
 
 
+void Foam::printTable(const List<wordList>& wll, Ostream& os)
+{
+    List<string::size_type> columnWidth;
+    printTable(wll, columnWidth, os);
+}
+
+
 // ************************************************************************* //
diff --git a/src/OpenFOAM/primitives/strings/word/wordIOList.H b/src/OpenFOAM/primitives/strings/word/wordIOList.H
index c560a33a95..da9d5aa426 100644
--- a/src/OpenFOAM/primitives/strings/word/wordIOList.H
+++ b/src/OpenFOAM/primitives/strings/word/wordIOList.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -43,6 +43,7 @@ namespace Foam
     typedef IOList<wordList> wordListIOList;
 
     // Print word list list as a table
+    void printTable(const List<wordList>&, List<string::size_type>&, Ostream&);
     void printTable(const List<wordList>&, Ostream&);
 }
 

From 1df4f9e8f6555b41091391b6ed26817f8ee07454 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 14 May 2013 20:48:38 +0100
Subject: [PATCH 13/33] ENH: lduPrimitiveMesh: make constructor easier

---
 src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.C | 6 ++++--
 src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.H | 2 +-
 2 files changed, 5 insertions(+), 3 deletions(-)

diff --git a/src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.C b/src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.C
index 070edd4b69..507fc425ef 100644
--- a/src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.C
+++ b/src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.C
@@ -254,7 +254,7 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
     const label nCells,
     labelList& l,
     labelList& u,
-    const Xfer<PtrList<const lduInterface> >& primitiveInterfaces,
+    PtrList<const lduInterface>& primitiveInterfaces,
     const lduSchedule& ps,
     const label comm
 )
@@ -262,10 +262,12 @@ Foam::lduPrimitiveMesh::lduPrimitiveMesh
     lduAddressing(nCells),
     lowerAddr_(l, true),
     upperAddr_(u, true),
-    primitiveInterfaces_(primitiveInterfaces),
+    primitiveInterfaces_(0),
     patchSchedule_(ps),
     comm_(comm)
 {
+    primitiveInterfaces_.transfer(primitiveInterfaces);
+
     // Create interfaces
     interfaces_.setSize(primitiveInterfaces_.size());
     forAll(primitiveInterfaces_, i)
diff --git a/src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.H b/src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.H
index 66c26dc37e..5bd7901c8b 100644
--- a/src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.H
+++ b/src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.H
@@ -133,7 +133,7 @@ public:
             const label nCells,
             labelList& l,
             labelList& u,
-            const Xfer<PtrList<const lduInterface> >& primitiveInterfaces,
+            PtrList<const lduInterface>& primitiveInterfaces,
             const lduSchedule& ps,
             const label comm
         );

From f3b7047a75e3dfcef4b859b2e29eba0e4a7ff55c Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 14 May 2013 20:50:33 +0100
Subject: [PATCH 14/33] ENH: procFaces: new GAMGProcAgglomeration method

---
 src/OpenFOAM/Make/files                       |   6 +-
 .../GAMGAgglomerateLduAddressing.C            |   2 +-
 .../GAMGAgglomeration/GAMGAgglomeration.C     |  99 +++++-
 .../pairGAMGAgglomeration.H                   |  17 +-
 .../procFacesGAMGProcAgglomeration.C          | 333 ++++++++++++++++++
 .../procFacesGAMGProcAgglomeration.H          | 136 +++++++
 6 files changed, 566 insertions(+), 27 deletions(-)
 create mode 100644 src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGProcAgglomerations/procFacesGAMGProcAgglomeration/procFacesGAMGProcAgglomeration.C
 create mode 100644 src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGProcAgglomerations/procFacesGAMGProcAgglomeration/procFacesGAMGProcAgglomeration.H

diff --git a/src/OpenFOAM/Make/files b/src/OpenFOAM/Make/files
index 19a8138fe5..6d346617c1 100644
--- a/src/OpenFOAM/Make/files
+++ b/src/OpenFOAM/Make/files
@@ -331,10 +331,8 @@ eagerGAMGProcAgglomeration = $(GAMGProcAgglomerations)/eagerGAMGProcAgglomeratio
 $(eagerGAMGProcAgglomeration)/eagerGAMGProcAgglomeration.C
 noneGAMGProcAgglomeration = $(GAMGProcAgglomerations)/noneGAMGProcAgglomeration
 $(noneGAMGProcAgglomeration)/noneGAMGProcAgglomeration.C
-/*
-cellFaceRatioGAMGProcAgglomeration = $(GAMGProcAgglomerations)/cellFaceRatioGAMGProcAgglomeration
-$(cellFaceRatioGAMGProcAgglomeration)/cellFaceRatioGAMGProcAgglomeration.C
-*/
+procFacesGAMGProcAgglomeration = $(GAMGProcAgglomerations)/procFacesGAMGProcAgglomeration
+$(procFacesGAMGProcAgglomeration)/procFacesGAMGProcAgglomeration.C
 
 
 meshes/lduMesh/lduMesh.C
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomerateLduAddressing.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomerateLduAddressing.C
index e6781ed23f..066aa39d3b 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomerateLduAddressing.C
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomerateLduAddressing.C
@@ -351,7 +351,7 @@ void Foam::GAMGAgglomeration::agglomerateLduAddressing
     );
 
 
-    if (debug)
+    if (debug & 2)
     {
         Pout<< "GAMGAgglomeration :"
             << " agglomerated level " << fineLevelIndex
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.C
index d1551eb399..a53a5fe9ab 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.C
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.C
@@ -27,9 +27,9 @@ License
 #include "lduMesh.H"
 #include "lduMatrix.H"
 #include "Time.H"
-#include "dlLibraryTable.H"
 #include "GAMGInterface.H"
 #include "GAMGProcAgglomeration.H"
+#include "IOmanip.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -67,23 +67,94 @@ void Foam::GAMGAgglomeration::compactLevels(const label nCreatedLevels)
         procBoundaryFaceMap_.setSize(nCreatedLevels);
 
         procAgglomeratorPtr_().agglomerate();
-    }
 
-    if (debug)
-    {
-        for (label levelI = 0; levelI <= size(); levelI++)
+
+        if (debug)
         {
-            if (hasMeshLevel(levelI))
+
+            Info<< "GAMGAgglomeration:" << nl
+                << "    local agglomerator     : " << type() << nl
+                << "    processor agglomerator : "
+                << procAgglomeratorPtr_().type() << nl
+                << nl;
+
+            Info<< setw(40) << "nCells"
+                << setw(24) << "nInterfaces"
+                << setw(24) << "Ratio" << nl
+                << setw(8) << "Level"
+                << setw(8) << "nProcs"
+                << "        "
+                << setw(8) << "avg"
+                << setw(8) << "max"
+                << "        "
+                << setw(8) << "avg"
+                << setw(8) << "max"
+                << "            " << setw(4) << "avg"
+                << "    " << setw(4) << "max"
+                << nl
+                << setw(8) << "-----"
+                << setw(8) << "------"
+                << "        "
+                << setw(8) << "---"
+                << setw(8) << "---"
+                << "        "
+                << setw(8) << "---"
+                << setw(8) << "---"
+                << "            " << setw(4) << "---"
+                << "    " << setw(4) << "---"
+                << nl;
+
+            for (label levelI = 0; levelI <= size(); levelI++)
             {
-                const lduMesh& fineMesh = meshLevel(levelI);
-                Pout<< "Level " << levelI << " fine mesh:"<< nl;
-                Pout<< fineMesh.info() << endl;
-            }
-            else
-            {
-                Pout<< "Level " << levelI << " has no fine mesh:" << nl
-                    << endl;
+                label nProcs = 0;
+                label nCells = 0;
+                label nInterfaces = 0;
+                label nIntFaces = 0;
+                scalar ratio = 0.0;
+
+                if (hasMeshLevel(levelI))
+                {
+                    nProcs = 1;
+
+                    const lduMesh& fineMesh = meshLevel(levelI);
+                    nCells = fineMesh.lduAddr().size();
+
+                    const lduInterfacePtrsList interfaces =
+                        fineMesh.interfaces();
+                    forAll(interfaces, i)
+                    {
+                        if (interfaces.set(i))
+                        {
+                            nInterfaces++;
+                            nIntFaces += interfaces[i].faceCells().size();
+                        }
+                    }
+                    ratio = scalar(nIntFaces)/nCells;
+                }
+
+                label totNprocs = returnReduce(nProcs, sumOp<label>());
+
+                label maxNCells = returnReduce(nCells, maxOp<label>());
+                label totNCells = returnReduce(nCells, sumOp<label>());
+
+                label maxNInt = returnReduce(nInterfaces, maxOp<label>());
+                label totNInt = returnReduce(nInterfaces, sumOp<label>());
+
+                scalar maxRatio = returnReduce(ratio, maxOp<scalar>());
+                scalar totRatio = returnReduce(ratio, sumOp<scalar>());
+
+                Info<< setw(8) << levelI
+                    << setw(8) << totNprocs
+                    << setw(16) << totNCells/totNprocs
+                    << setw(8) << maxNCells
+                    << setw(16) << totNInt/totNprocs
+                    << setw(8) << maxNInt
+                    << "        "
+                    << setw(8) << setprecision(4) << totRatio/totNprocs
+                    << setw(8) << setprecision(4) << maxRatio
+                    << nl;
             }
+            Info<< endl;
         }
     }
 }
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/pairGAMGAgglomeration/pairGAMGAgglomeration.H b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/pairGAMGAgglomeration/pairGAMGAgglomeration.H
index 2867cac2f6..59de4b1e4d 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/pairGAMGAgglomeration/pairGAMGAgglomeration.H
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/pairGAMGAgglomeration/pairGAMGAgglomeration.H
@@ -61,14 +61,6 @@ protected:
 
     // Protected Member Functions
 
-        //- Calculate and return agglomeration of given level
-        tmp<labelField> agglomerate
-        (
-            label& nCoarseCells,
-            const lduAddressing& fineMatrixAddressing,
-            const scalarField& faceWeights
-        );
-
         //- Agglomerate all levels starting from the given face weights
         void agglomerate
         (
@@ -97,6 +89,15 @@ public:
             const lduMesh& mesh,
             const dictionary& controlDict
         );
+
+        //- Calculate and return agglomeration
+        static tmp<labelField> agglomerate
+        (
+            label& nCoarseCells,
+            const lduAddressing& fineMatrixAddressing,
+            const scalarField& faceWeights
+        );
+
 };
 
 
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGProcAgglomerations/procFacesGAMGProcAgglomeration/procFacesGAMGProcAgglomeration.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGProcAgglomerations/procFacesGAMGProcAgglomeration/procFacesGAMGProcAgglomeration.C
new file mode 100644
index 0000000000..6e13d52c99
--- /dev/null
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGProcAgglomerations/procFacesGAMGProcAgglomeration/procFacesGAMGProcAgglomeration.C
@@ -0,0 +1,333 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "procFacesGAMGProcAgglomeration.H"
+#include "addToRunTimeSelectionTable.H"
+#include "GAMGAgglomeration.H"
+#include "Random.H"
+#include "lduMesh.H"
+#include "processorLduInterface.H"
+#include "processorGAMGInterface.H"
+#include "pairGAMGAgglomeration.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(procFacesGAMGProcAgglomeration, 0);
+
+    addToRunTimeSelectionTable
+    (
+        GAMGProcAgglomeration,
+        procFacesGAMGProcAgglomeration,
+        GAMGAgglomeration
+    );
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+// Create single cell mesh
+Foam::autoPtr<Foam::lduPrimitiveMesh>
+Foam::procFacesGAMGProcAgglomeration::singleCellMesh
+(
+    const label singleCellMeshComm,
+    const lduMesh& mesh,
+    scalarField& faceWeights
+) const
+{
+    // Count number of faces per processor
+    List<Map<label> > procFaces(UPstream::nProcs(mesh.comm()));
+    Map<label>& myNeighbours = procFaces[UPstream::myProcNo(mesh.comm())];
+
+    {
+        const lduInterfacePtrsList interfaces(mesh.interfaces());
+        forAll(interfaces, intI)
+        {
+            if (interfaces.set(intI))
+            {
+                const processorLduInterface& pp =
+                    refCast<const processorLduInterface>
+                    (
+                        interfaces[intI]
+                    );
+                label size = interfaces[intI].faceCells().size();
+                myNeighbours.insert(pp.neighbProcNo(), size);
+            }
+        }
+    }
+
+    Pstream::gatherList(procFaces, Pstream::msgType(), mesh.comm());
+    Pstream::scatterList(procFaces, Pstream::msgType(), mesh.comm());
+
+    autoPtr<lduPrimitiveMesh> singleCellMeshPtr;
+
+    if (Pstream::master(mesh.comm()))
+    {
+        // I am master
+        label nCells = Pstream::nProcs(mesh.comm());
+
+        DynamicList<label> l(3*nCells);
+        DynamicList<label> u(3*nCells);
+        DynamicList<scalar> weight(3*nCells);
+
+        DynamicList<label> nbrs;
+        DynamicList<scalar> weights;
+
+        forAll(procFaces, procI)
+        {
+            const Map<label>& neighbours = procFaces[procI];
+
+            // Add all the higher processors
+            nbrs.clear();
+            weights.clear();
+            forAllConstIter(Map<label>, neighbours, iter)
+            {
+                if (iter.key() > procI)
+                {
+                    nbrs.append(iter.key());
+                    weights.append(iter());
+                }
+                sort(nbrs);
+                forAll(nbrs, i)
+                {
+                    l.append(procI);
+                    u.append(nbrs[i]);
+                    weight.append(weights[i]);
+                }
+            }
+        }
+
+        faceWeights.transfer(weight);
+
+        PtrList<const lduInterface> primitiveInterfaces(0);
+        const lduSchedule ps(0);
+
+        singleCellMeshPtr.reset
+        (
+            new lduPrimitiveMesh
+            (
+                nCells,
+                l,
+                u,
+                primitiveInterfaces,
+                ps,
+                singleCellMeshComm
+            )
+        );
+    }
+    return singleCellMeshPtr;
+}
+
+
+Foam::tmp<Foam::labelField>
+Foam::procFacesGAMGProcAgglomeration::processorAgglomeration
+(
+    const lduMesh& mesh
+) const
+{
+    label singleCellMeshComm = UPstream::allocateCommunicator
+    (
+        mesh.comm(),
+        labelList(1, 0)            // only processor 0
+    );
+
+    scalarField faceWeights;
+    autoPtr<lduPrimitiveMesh> singleCellMeshPtr
+    (
+        singleCellMesh
+        (
+            singleCellMeshComm,
+            mesh,
+            faceWeights
+        )
+    );
+
+    tmp<labelField> tfineToCoarse(new labelField(0));
+    labelField& fineToCoarse = tfineToCoarse();
+
+    if (singleCellMeshPtr.valid())
+    {
+        // On master call the agglomerator
+        const lduPrimitiveMesh& singleCellMesh = singleCellMeshPtr();
+
+        label nCoarseProcs;
+        fineToCoarse = pairGAMGAgglomeration::agglomerate
+        (
+            nCoarseProcs,
+            singleCellMesh,
+            faceWeights
+        );
+
+        labelList coarseToMaster(nCoarseProcs, labelMax);
+        forAll(fineToCoarse, cellI)
+        {
+            label coarseI = fineToCoarse[cellI];
+            coarseToMaster[coarseI] = min(coarseToMaster[coarseI], cellI);
+        }
+
+        // Sort according to master and redo restriction
+        labelList newToOld;
+        sortedOrder(coarseToMaster, newToOld);
+        labelList oldToNew(invert(newToOld.size(), newToOld));
+
+        fineToCoarse = UIndirectList<label>(oldToNew, fineToCoarse)();
+    }
+
+    Pstream::scatter(fineToCoarse, Pstream::msgType(), mesh.comm());
+    UPstream::freeCommunicator(singleCellMeshComm);
+
+    return tfineToCoarse;
+}
+
+
+bool Foam::procFacesGAMGProcAgglomeration::doProcessorAgglomeration
+(
+    const lduMesh& mesh
+) const
+{
+    // Check the need for further agglomeration on all processors
+    bool doAgg = mesh.lduAddr().size() < nAgglomeratingCells_;
+    mesh.reduce(doAgg, orOp<bool>());
+    return doAgg;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::procFacesGAMGProcAgglomeration::procFacesGAMGProcAgglomeration
+(
+    GAMGAgglomeration& agglom,
+    const dictionary& controlDict
+)
+:
+    GAMGProcAgglomeration(agglom, controlDict),
+    nAgglomeratingCells_(readLabel(controlDict.lookup("nAgglomeratingCells")))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::procFacesGAMGProcAgglomeration::~procFacesGAMGProcAgglomeration()
+{
+    forAllReverse(comms_, i)
+    {
+        if (comms_[i] != -1)
+        {
+            UPstream::freeCommunicator(comms_[i]);
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::procFacesGAMGProcAgglomeration::agglomerate()
+{
+    if (debug)
+    {
+        Pout<< nl << "Starting mesh overview" << endl;
+        printStats(Pout, agglom_);
+    }
+
+    if (agglom_.size() >= 1)
+    {
+        Random rndGen(0);
+
+        for
+        (
+            label fineLevelIndex = 2;
+            fineLevelIndex < agglom_.size();
+            fineLevelIndex++
+        )
+        {
+            if (agglom_.hasMeshLevel(fineLevelIndex))
+            {
+                // Get the fine mesh
+                const lduMesh& levelMesh = agglom_.meshLevel(fineLevelIndex);
+
+                label levelComm = levelMesh.comm();
+                label nProcs = UPstream::nProcs(levelComm);
+
+                if (nProcs > 1 && doProcessorAgglomeration(levelMesh))
+                {
+                    tmp<labelField> tprocAgglomMap
+                    (
+                        processorAgglomeration(levelMesh)
+                    );
+                    const labelField& procAgglomMap = tprocAgglomMap();
+
+                    // Master processor
+                    labelList masterProcs;
+                    // Local processors that agglomerate. agglomProcIDs[0] is in
+                    // masterProc.
+                    List<int> agglomProcIDs;
+                    GAMGAgglomeration::calculateRegionMaster
+                    (
+                        levelComm,
+                        procAgglomMap,
+                        masterProcs,
+                        agglomProcIDs
+                    );
+
+                    // Allocate a communicator for the processor-agglomerated
+                    // matrix
+                    comms_.append
+                    (
+                        UPstream::allocateCommunicator
+                        (
+                            levelComm,
+                            masterProcs
+                        )
+                    );
+
+
+                    // Use procesor agglomeration maps to do the actual
+                    // collecting.
+                    GAMGProcAgglomeration::agglomerate
+                    (
+                        fineLevelIndex,
+                        procAgglomMap,
+                        masterProcs,
+                        agglomProcIDs,
+                        comms_.last()
+                    );
+                }
+            }
+        }
+    }
+
+    // Print a bit
+    if (debug)
+    {
+        Pout<< nl << "Agglomerated mesh overview" << endl;
+        printStats(Pout, agglom_);
+    }
+
+    return true;
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGProcAgglomerations/procFacesGAMGProcAgglomeration/procFacesGAMGProcAgglomeration.H b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGProcAgglomerations/procFacesGAMGProcAgglomeration/procFacesGAMGProcAgglomeration.H
new file mode 100644
index 0000000000..5c574c4158
--- /dev/null
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGProcAgglomerations/procFacesGAMGProcAgglomeration/procFacesGAMGProcAgglomeration.H
@@ -0,0 +1,136 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::procFacesGAMGProcAgglomeration
+
+Description
+    Processor agglomeration of GAMGAgglomerations. Needs nAgglomeratingCells
+    which is when to start agglomerating processors. Processors get agglomerated
+    by constructing a single cell mesh for each processor with each
+    processor interface a face. This then gets agglomerated using
+    the pairGAMGAgglomeration algorithm with the number of faces
+    on the original processor interface as face weight.
+
+SourceFiles
+    procFacesGAMGProcAgglomeration.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef procFacesGAMGProcAgglomeration_H
+#define procFacesGAMGProcAgglomeration_H
+
+#include "GAMGProcAgglomeration.H"
+#include "DynamicList.H"
+#include "labelField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class GAMGAgglomeration;
+class lduMesh;
+class lduPrimitiveMesh;
+
+/*---------------------------------------------------------------------------*\
+                    Class procFacesGAMGProcAgglomeration Declaration
+\*---------------------------------------------------------------------------*/
+
+class procFacesGAMGProcAgglomeration
+:
+    public GAMGProcAgglomeration
+{
+    // Private data
+
+        //- When to processor agglomerate
+        const label nAgglomeratingCells_;
+
+        //- Allocated communicators
+        DynamicList<label> comms_;
+
+    // Private Member Functions
+
+        //- Return (on master) all single-cell meshes collected. single-cell
+        //  meshes are just one cell with all proc faces intact.
+        autoPtr<lduPrimitiveMesh> singleCellMesh
+        (
+            const label singleCellMeshComm,
+            const lduMesh& mesh,
+            scalarField& faceWeights
+        ) const;
+
+        //- Construct processor agglomeration: for every processor the
+        //  coarse processor-cluster it agglomerates onto
+        tmp<labelField> processorAgglomeration(const lduMesh&) const;
+
+        //- Do we need to agglomerate across processors?
+        bool doProcessorAgglomeration(const lduMesh&) const;
+
+        //- Disallow default bitwise copy construct
+        procFacesGAMGProcAgglomeration
+        (
+            const procFacesGAMGProcAgglomeration&
+        );
+
+        //- Disallow default bitwise assignment
+        void operator=(const procFacesGAMGProcAgglomeration&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("procFaces");
+
+
+    // Constructors
+
+        //- Construct given agglomerator and controls
+        procFacesGAMGProcAgglomeration
+        (
+            GAMGAgglomeration& agglom,
+            const dictionary& controlDict
+        );
+
+
+    //- Destructor
+    virtual ~procFacesGAMGProcAgglomeration();
+
+
+    // Member Functions
+
+       //- Modify agglomeration. Return true if modified
+        virtual bool agglomerate();
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

From caebdcc8953471c3b00832475448ff31674a73aa Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 14 May 2013 20:56:13 +0100
Subject: [PATCH 15/33] STYLE: meshObject: fix start date

---
 src/OpenFOAM/meshes/MeshObject/meshObject.C | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/OpenFOAM/meshes/MeshObject/meshObject.C b/src/OpenFOAM/meshes/MeshObject/meshObject.C
index 7130f38b5d..a616f777b0 100644
--- a/src/OpenFOAM/meshes/MeshObject/meshObject.C
+++ b/src/OpenFOAM/meshes/MeshObject/meshObject.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

From fe32d73d5af22623afdc8422bfa038759d494bf4 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Thu, 16 May 2013 11:17:54 +0100
Subject: [PATCH 16/33] BUG: distributedTriSurfaceMesh: handle zero sized local
 surfaces

---
 .../searchableSurface/triSurfaceMesh.C        | 34 +++++++-----
 .../triSurfaceSearch/triSurfaceRegionSearch.C | 55 +++++++++++--------
 .../triSurfaceSearch/triSurfaceSearch.C       | 41 ++++++++------
 .../distributedTriSurfaceMesh.C               | 24 +++++++-
 .../distributedTriSurfaceMesh.H               |  6 +-
 5 files changed, 99 insertions(+), 61 deletions(-)

diff --git a/src/meshTools/searchableSurface/triSurfaceMesh.C b/src/meshTools/searchableSurface/triSurfaceMesh.C
index 2222981793..6009ec6670 100644
--- a/src/meshTools/searchableSurface/triSurfaceMesh.C
+++ b/src/meshTools/searchableSurface/triSurfaceMesh.C
@@ -433,24 +433,28 @@ Foam::triSurfaceMesh::edgeTree() const
           + nInternalEdges()
         );
 
-        treeBoundBox bb;
-        label nPoints;
-        PatchTools::calcBounds
-        (
-            static_cast<const triSurface&>(*this),
-            bb,
-            nPoints
-        );
+        treeBoundBox bb(vector::zero, vector::zero);
 
-        // Random number generator. Bit dodgy since not exactly random ;-)
-        Random rndGen(65431);
+        if (bEdges.size())
+        {
+            label nPoints;
+            PatchTools::calcBounds
+            (
+                static_cast<const triSurface&>(*this),
+                bb,
+                nPoints
+            );
 
-        // Slightly extended bb. Slightly off-centred just so on symmetric
-        // geometry there are less face/edge aligned items.
+            // Random number generator. Bit dodgy since not exactly random ;-)
+            Random rndGen(65431);
 
-        bb = bb.extend(rndGen, 1e-4);
-        bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
-        bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+            // Slightly extended bb. Slightly off-centred just so on symmetric
+            // geometry there are less face/edge aligned items.
+
+            bb = bb.extend(rndGen, 1e-4);
+            bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+            bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+        }
 
         scalar oldTol = indexedOctree<treeDataEdge>::perturbTol();
         indexedOctree<treeDataEdge>::perturbTol() = tolerance();
diff --git a/src/meshTools/triSurface/triSurfaceSearch/triSurfaceRegionSearch.C b/src/meshTools/triSurface/triSurfaceSearch/triSurfaceRegionSearch.C
index 738cb8dcdd..34dc89f920 100644
--- a/src/meshTools/triSurface/triSurfaceSearch/triSurfaceRegionSearch.C
+++ b/src/meshTools/triSurface/triSurfaceSearch/triSurfaceRegionSearch.C
@@ -121,33 +121,40 @@ Foam::triSurfaceRegionSearch::treeByRegion() const
             );
 
             // Calculate bb without constructing local point numbering.
-            treeBoundBox bb;
-            label nPoints;
-            PatchTools::calcBounds
-            (
-                indirectRegionPatches_[regionI],
-                bb,
-                nPoints
-            );
+            treeBoundBox bb(vector::zero, vector::zero);
 
-//            if (nPoints != surface().points().size())
-//            {
-//                WarningIn("triSurfaceRegionSearch::treeByRegion() const")
-//                    << "Surface does not have compact point numbering. "
-//                    << "Of " << surface().points().size()
-//                    << " only " << nPoints
-//                    << " are used. This might give problems in some routines."
-//                    << endl;
-//            }
+            if (indirectRegionPatches_[regionI].size())
+            {
+                label nPoints;
+                PatchTools::calcBounds
+                (
+                    indirectRegionPatches_[regionI],
+                    bb,
+                    nPoints
+                );
 
-            // Random number generator. Bit dodgy since not exactly random ;-)
-            Random rndGen(65431);
+    //            if (nPoints != surface().points().size())
+    //            {
+    //                WarningIn("triSurfaceRegionSearch::treeByRegion() const")
+    //                    << "Surface does not have compact point numbering. "
+    //                    << "Of " << surface().points().size()
+    //                    << " only " << nPoints
+    //                    << " are used."
+    //                    << " This might give problems in some routines."
+    //                    << endl;
+    //            }
 
-            // Slightly extended bb. Slightly off-centred just so on symmetric
-            // geometry there are fewer face/edge aligned items.
-            bb = bb.extend(rndGen, 1e-4);
-            bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
-            bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+                // Random number generator. Bit dodgy since not exactly
+                // random ;-)
+                Random rndGen(65431);
+
+                // Slightly extended bb. Slightly off-centred just so
+                // on symmetric geometry there are fewer face/edge
+                // aligned items.
+                bb = bb.extend(rndGen, 1e-4);
+                bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+                bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+            }
 
             treeByRegion_.set
             (
diff --git a/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.C b/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.C
index fbff791040..43877e3393 100644
--- a/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.C
+++ b/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.C
@@ -200,28 +200,33 @@ Foam::triSurfaceSearch::tree() const
     if (treePtr_.empty())
     {
         // Calculate bb without constructing local point numbering.
-        treeBoundBox bb;
-        label nPoints;
-        PatchTools::calcBounds(surface(), bb, nPoints);
+        treeBoundBox bb(vector::zero, vector::zero);
 
-        if (nPoints != surface().points().size())
+        if (surface().size())
         {
-            WarningIn("triSurfaceSearch::tree() const")
-                << "Surface does not have compact point numbering."
-                << " Of " << surface().points().size() << " only " << nPoints
-                << " are used. This might give problems in some routines."
-                << endl;
+            label nPoints;
+            PatchTools::calcBounds(surface(), bb, nPoints);
+
+            if (nPoints != surface().points().size())
+            {
+                WarningIn("triSurfaceSearch::tree() const")
+                    << "Surface does not have compact point numbering."
+                    << " Of " << surface().points().size()
+                    << " only " << nPoints
+                    << " are used. This might give problems in some routines."
+                    << endl;
+            }
+
+            // Random number generator. Bit dodgy since not exactly random ;-)
+            Random rndGen(65431);
+
+            // Slightly extended bb. Slightly off-centred just so on symmetric
+            // geometry there are less face/edge aligned items.
+            bb = bb.extend(rndGen, 1e-4);
+            bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+            bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
         }
 
-        // Random number generator. Bit dodgy since not exactly random ;-)
-        Random rndGen(65431);
-
-        // Slightly extended bb. Slightly off-centred just so on symmetric
-        // geometry there are less face/edge aligned items.
-        bb = bb.extend(rndGen, 1e-4);
-        bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
-        bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
-
         scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
         indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;
 
diff --git a/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C b/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C
index 80308fac5e..7c80bd1a56 100644
--- a/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C
+++ b/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C
@@ -2360,10 +2360,28 @@ bool Foam::distributedTriSurfaceMesh::writeObject
     // Make sure dictionary goes to same directory as surface
     const_cast<fileName&>(dict_.instance()) = searchableSurface::instance();
 
+    // Copy of triSurfaceMesh::writeObject except for the sorting of
+    // triangles by region. This is done so we preserve region names,
+    // even if locally we have zero triangles.
+    {
+        fileName fullPath(searchableSurface::objectPath());
+
+        if (!mkDir(fullPath.path()))
+        {
+            return false;
+        }
+
+        // Important: preserve any zero-sized patches
+        triSurface::write(fullPath, true);
+
+        if (!isFile(fullPath))
+        {
+            return false;
+        }
+    }
+
     // Dictionary needs to be written in ascii - binary output not supported.
-    return
-        triSurfaceMesh::writeObject(fmt, ver, cmp)
-     && dict_.writeObject(IOstream::ASCII, ver, cmp);
+    return dict_.writeObject(IOstream::ASCII, ver, cmp);
 }
 
 
diff --git a/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.H b/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.H
index 32a062a573..7e4d27ef0c 100644
--- a/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.H
+++ b/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.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
@@ -446,6 +446,10 @@ public:
         // regIOobject implementation
 
             //- Write using given format, version and compression
+            //  Do not use the triSurfaceMesh::writeObject since it
+            //  would filter out empty regions. These need to be preserved
+            //  in case we want to make decisions based on the number of
+            //  regions.
             virtual bool writeObject
             (
                 IOstream::streamFormat fmt,

From 5f68e587e6fec4ff2f7b16e2e47367ff7a9e5e87 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Thu, 16 May 2013 11:28:02 +0100
Subject: [PATCH 17/33] ENH: viewFactorsGen: exclude AMI

---
 .../utilities/preProcessing/viewFactorsGen/searchingEngine.H   | 3 ++-
 .../utilities/preProcessing/viewFactorsGen/viewFactorsGen.C    | 1 +
 2 files changed, 3 insertions(+), 1 deletion(-)

diff --git a/applications/utilities/preProcessing/viewFactorsGen/searchingEngine.H b/applications/utilities/preProcessing/viewFactorsGen/searchingEngine.H
index ed1d9031ec..ef5b264a54 100644
--- a/applications/utilities/preProcessing/viewFactorsGen/searchingEngine.H
+++ b/applications/utilities/preProcessing/viewFactorsGen/searchingEngine.H
@@ -26,7 +26,8 @@ dict.add("mergeDistance", SMALL);
 labelHashSet includePatches;
 forAll(patches, patchI)
 {
-    if (!isA<processorPolyPatch>(patches[patchI]))
+    const polyPatch& pp = patches[patchI];
+    if (!pp.coupled() && !isA<cyclicAMIPolyPatch>(pp))
     {
         includePatches.insert(patchI);
     }
diff --git a/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C b/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
index c9aa219d13..8276528945 100644
--- a/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
+++ b/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
@@ -44,6 +44,7 @@ Description
 #include "volFields.H"
 #include "surfaceFields.H"
 #include "distributedTriSurfaceMesh.H"
+#include "cyclicAMIPolyPatch.H"
 #include "triSurfaceTools.H"
 #include "mapDistribute.H"
 

From 25a373bd395d1e8335849cd159a1f29c7f7d216f Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Thu, 16 May 2013 12:49:16 +0100
Subject: [PATCH 18/33] ENH: snappyHexMesh: enable free-standing baffle merging

---
 .../snappyHexMesh/snappyHexMeshDict           |   6 +
 .../autoHexMeshDriver/autoRefineDriver.C      |   5 +-
 .../refinementParameters.C                    |  11 +-
 .../refinementParameters.H                    |  11 +-
 .../meshRefinement/meshRefinement.H           |  35 +-
 .../meshRefinement/meshRefinementBaffles.C    | 612 +++++++++++-------
 .../meshRefinement/meshRefinementRefine.C     |  25 +-
 7 files changed, 445 insertions(+), 260 deletions(-)

diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
index f448faa40c..ca4a9dc6b5 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
@@ -180,6 +180,12 @@ castellatedMeshControls
     // - used if feature snapping (see snapControls below) is used
     resolveFeatureAngle 30;
 
+    // Planar angle:
+    // - used to determine if neighbouring surface normals
+    //   are roughly the same so e.g. free-standing baffles can be merged
+    // If not specified same as resolveFeatureAngle
+    //planarAngle 30;
+
 
     // Region-wise refinement
     // ~~~~~~~~~~~~~~~~~~~~~~
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
index 484f851c0f..75ea461102 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
@@ -522,6 +522,7 @@ void Foam::autoRefineDriver::baffleAndSplitMesh
         false,                          // perpendicular edge connected cells
         scalarField(0),                 // per region perpendicular angle
         !handleSnapProblems,            // merge free standing baffles?
+        refineParams.planarAngle(),
         motionDict,
         const_cast<Time&>(mesh.time()),
         globalToMasterPatch_,
@@ -606,8 +607,8 @@ void Foam::autoRefineDriver::splitAndMergeBaffles
         handleSnapProblems,
         handleSnapProblems,                 // remove perp edge connected cells
         perpAngle,                          // perp angle
-        false,                              // merge free standing baffles?
-        //true,                               // merge free standing baffles?
+        true,                               // merge free standing baffles?
+        refineParams.planarAngle(),         // planar angle
         motionDict,
         const_cast<Time&>(mesh.time()),
         globalToMasterPatch_,
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C
index ab2140a4ec..7bf8d8cf6c 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.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
@@ -41,6 +41,7 @@ Foam::refinementParameters::refinementParameters
     maxLocalCells_(readLabel(dict.lookup("procCellLimit"))),
     minRefineCells_(readLabel(dict.lookup("minimumRefine"))),
     curvature_(readScalar(dict.lookup("curvature"))),
+    planarAngle_(dict.lookupOrDefault("planarAngle", curvature_)),
     nBufferLayers_(readLabel(dict.lookup("nBufferLayers"))),
     keepPoints_(dict.lookup("keepPoints")),
     allowFreeStandingZoneFaces_(dict.lookup("allowFreeStandingZoneFaces")),
@@ -53,6 +54,14 @@ Foam::refinementParameters::refinementParameters(const dictionary& dict)
     maxGlobalCells_(readLabel(dict.lookup("maxGlobalCells"))),
     maxLocalCells_(readLabel(dict.lookup("maxLocalCells"))),
     minRefineCells_(readLabel(dict.lookup("minRefinementCells"))),
+    planarAngle_
+    (
+        dict.lookupOrDefault
+        (
+            "planarAngle",
+            readScalar(dict.lookup("resolveFeatureAngle"))
+        )
+    ),
     nBufferLayers_(readLabel(dict.lookup("nCellsBetweenLevels"))),
     keepPoints_(pointField(1, dict.lookup("locationInMesh"))),
     allowFreeStandingZoneFaces_(dict.lookup("allowFreeStandingZoneFaces")),
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.H
index 66867e7f55..4fb05644e8 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.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
@@ -67,6 +67,9 @@ class refinementParameters
         //- Curvature
         scalar curvature_;
 
+        //- Planarity criterion
+        scalar planarAngle_;
+
         //- Number of layers between different refinement levels
         const label nBufferLayers_;
 
@@ -129,6 +132,12 @@ public:
                 return curvature_;
             }
 
+            //- Angle when two intersections are considered to be planar
+            scalar planarAngle() const
+            {
+                return planarAngle_;
+            }
+
             //- Number of layers between different refinement levels
             label nBufferLayers() const
             {
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
index 6ab861d8a4..60e6261ccc 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
@@ -253,6 +253,14 @@ private:
                 label& nRefine
             );
 
+            //- Mark every cell with level of feature passing through it
+            //  (or -1 if not passed through). Uses tracking.
+            void markFeatureCellLevel
+            (
+                const point& keepPoint,
+                labelList& maxFeatureLevel
+            ) const;
+
             //- Calculate list of cells to refine based on intersection of
             //  features.
             label markFeatureRefinement
@@ -337,16 +345,6 @@ private:
                 const labelList& globalToSlavePatch
             ) const;
 
-            //- Geometric test on see whether face needs to be baffled:
-            //  is face boundary face and perpendicular to surface normal?
-            bool validBaffleTopology
-            (
-                const label faceI,
-                const vector& n1,
-                const vector& n2,
-                const vector& testDir
-            ) const;
-
             //- Determine patches for baffles
             void getBafflePatches
             (
@@ -435,7 +433,11 @@ private:
 
             //- Extract those baffles (duplicate) faces that are on the edge
             //  of a baffle region. These are candidates for merging.
-            List<labelPair> freeStandingBaffles(const List<labelPair>&) const;
+            List<labelPair> freeStandingBaffles
+            (
+                const List<labelPair>&,
+                const scalar freeStandingAngle
+            ) const;
 
 
         // Zone handling
@@ -489,6 +491,16 @@ private:
                 labelList& namedSurfaceIndex
             ) const;
 
+            //- Remove any loose standing cells
+            void handleSnapProblems
+            (
+                const bool removeEdgeConnectedCells,
+                const scalarField& perpendicularAngle,
+                const dictionary& motionDict,
+                Time& runTime,
+                const labelList& globalToMasterPatch,
+                const labelList& globalToSlavePatch
+            );
 
         //- Disallow default bitwise copy construct
         meshRefinement(const meshRefinement&);
@@ -707,6 +719,7 @@ public:
                 const bool removeEdgeConnectedCells,
                 const scalarField& perpendicularAngle,
                 const bool mergeFreeStanding,
+                const scalar freeStandingAngle,
                 const dictionary& motionDict,
                 Time& runTime,
                 const labelList& globalToMasterPatch,
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
index 339e51cd7f..7a4f87753d 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
@@ -44,6 +44,7 @@ License
 #include "regionSplit.H"
 #include "removeCells.H"
 #include "unitConversion.H"
+#include "OBJstream.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -217,43 +218,43 @@ Foam::label Foam::meshRefinement::getBafflePatch
 }
 
 
-// Check if we are a boundary face and normal of surface does
-// not align with test vector. In this case there'd probably be
-// a freestanding 'baffle' so we might as well not create it.
-// Note that since it is not a proper baffle we cannot detect it
-// afterwards so this code cannot be merged with the
-// filterDuplicateFaces code.
-bool Foam::meshRefinement::validBaffleTopology
-(
-    const label faceI,
-    const vector& n1,
-    const vector& n2,
-    const vector& testDir
-) const
-{
-
-    label patchI = mesh_.boundaryMesh().whichPatch(faceI);
-    if (patchI == -1 || mesh_.boundaryMesh()[patchI].coupled())
-    {
-        return true;
-    }
-    else if (mag(n1&n2) > cos(degToRad(30)))
-    {
-        // Both normals aligned. Check that test vector perpendicularish to
-        // surface normal
-        scalar magTestDir = mag(testDir);
-        if (magTestDir > VSMALL)
-        {
-            if (mag(n1&(testDir/magTestDir)) < cos(degToRad(45)))
-            {
-                //Pout<< "** disabling baffling face "
-                //    << mesh_.faceCentres()[faceI] << endl;
-                return false;
-            }
-        }
-    }
-    return true;
-}
+//// Check if we are a boundary face and normal of surface does
+//// not align with test vector. In this case there'd probably be
+//// a freestanding 'baffle' so we might as well not create it.
+//// Note that since it is not a proper baffle we cannot detect it
+//// afterwards so this code cannot be merged with the
+//// filterDuplicateFaces code.
+//bool Foam::meshRefinement::validBaffleTopology
+//(
+//    const label faceI,
+//    const vector& n1,
+//    const vector& n2,
+//    const vector& testDir
+//) const
+//{
+//
+//    label patchI = mesh_.boundaryMesh().whichPatch(faceI);
+//    if (patchI == -1 || mesh_.boundaryMesh()[patchI].coupled())
+//    {
+//        return true;
+//    }
+//    else if (mag(n1&n2) > cos(degToRad(30)))
+//    {
+//        // Both normals aligned. Check that test vector perpendicularish to
+//        // surface normal
+//        scalar magTestDir = mag(testDir);
+//        if (magTestDir > VSMALL)
+//        {
+//            if (mag(n1&(testDir/magTestDir)) < cos(degToRad(45)))
+//            {
+//                //Pout<< "** disabling baffling face "
+//                //    << mesh_.faceCentres()[faceI] << endl;
+//                return false;
+//            }
+//        }
+//    }
+//    return true;
+//}
 
 
 // Determine patches for baffles on all intersected unnamed faces
@@ -800,7 +801,8 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createZoneBaffles
 // region.
 Foam::List<Foam::labelPair> Foam::meshRefinement::freeStandingBaffles
 (
-    const List<labelPair>& couples
+    const List<labelPair>& couples,
+    const scalar planarAngle
 ) const
 {
     // All duplicate faces on edge of the patch are to be merged.
@@ -809,6 +811,22 @@ Foam::List<Foam::labelPair> Foam::meshRefinement::freeStandingBaffles
     labelList nBafflesPerEdge(mesh_.nEdges(), 0);
 
 
+    // This algorithm is quite tricky. We don't want to use edgeFaces and
+    // also want it to run in parallel so it is now an algorithm over
+    // all (boundary) faces instead.
+    // We want to pick up any edges that are only used by the baffle
+    // or internal faces but not by any other boundary faces. So
+    // - increment count on an edge by 1 if it is used by any (uncoupled)
+    //   boundary face.
+    // - increment count on an edge by 1000000 if it is used by a baffle face
+    // - sum in parallel
+    //
+    // So now any edge that is used by baffle faces only will have the
+    // value 2*1000000+2*1.
+
+
+    const label baffleValue = 1000000;
+
 
     // Count number of boundary faces per edge
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -847,18 +865,22 @@ Foam::List<Foam::labelPair> Foam::meshRefinement::freeStandingBaffles
 
     forAll(couples, i)
     {
-        const labelList& fEdges0 = mesh_.faceEdges(couples[i].first(), fe0);
-
-        forAll(fEdges0, fEdgeI)
         {
-            nBafflesPerEdge[fEdges0[fEdgeI]] += 1000000;
+            label f0 = couples[i].first();
+            const labelList& fEdges0 = mesh_.faceEdges(f0, fe0);
+            forAll(fEdges0, fEdgeI)
+            {
+                nBafflesPerEdge[fEdges0[fEdgeI]] += baffleValue;
+            }
         }
 
-        const labelList& fEdges1 = mesh_.faceEdges(couples[i].second(), fe1);
-
-        forAll(fEdges1, fEdgeI)
         {
-            nBafflesPerEdge[fEdges1[fEdgeI]] += 1000000;
+            label f1 = couples[i].second();
+            const labelList& fEdges1 = mesh_.faceEdges(f1, fe1);
+            forAll(fEdges1, fEdgeI)
+            {
+                nBafflesPerEdge[fEdges1[fEdgeI]] += baffleValue;
+            }
         }
     }
 
@@ -873,8 +895,8 @@ Foam::List<Foam::labelPair> Foam::meshRefinement::freeStandingBaffles
 
 
     // Baffles which are not next to other boundaries and baffles will have
-    // nBafflesPerEdge value 2*1000000+2*1 (from 2 boundary faces which are
-    // both baffle faces)
+    // nBafflesPerEdge value 2*baffleValue+2*1 (from 2 boundary faces which
+    // are both baffle faces)
 
     List<labelPair> filteredCouples(couples.size());
     label filterI = 0;
@@ -889,15 +911,15 @@ Foam::List<Foam::labelPair> Foam::meshRefinement::freeStandingBaffles
          == patches.whichPatch(couple.second())
         )
         {
-            const labelList& fEdges = mesh_.faceEdges(couples[i].first());
+            const labelList& fEdges = mesh_.faceEdges(couple.first());
 
             forAll(fEdges, fEdgeI)
             {
                 label edgeI = fEdges[fEdgeI];
 
-                if (nBafflesPerEdge[edgeI] == 2*1000000+2*1)
+                if (nBafflesPerEdge[edgeI] == 2*baffleValue+2*1)
                 {
-                    filteredCouples[filterI++] = couples[i];
+                    filteredCouples[filterI++] = couple;
                     break;
                 }
             }
@@ -905,111 +927,127 @@ Foam::List<Foam::labelPair> Foam::meshRefinement::freeStandingBaffles
     }
     filteredCouples.setSize(filterI);
 
-    //Info<< "freeStandingBaffles : from "
-    //    << returnReduce(couples.size(), sumOp<label>())
-    //    << " down to "
-    //    << returnReduce(filteredCouples.size(), sumOp<label>())
-    //    << " baffles." << nl << endl;
+
+    label nFiltered = returnReduce(filteredCouples.size(), sumOp<label>());
+
+    Info<< "freeStandingBaffles : detected "
+        << nFiltered
+        << " free-standing baffles out of "
+        << returnReduce(couples.size(), sumOp<label>())
+        << nl << endl;
 
 
+    if (nFiltered > 0)
+    {
+        // Collect segments
+        // ~~~~~~~~~~~~~~~~
 
-//XXXXXX
-//    {
-//        // Collect segments
-//        // ~~~~~~~~~~~~~~~~
-//
-//        pointField start(filteredCouples.size());
-//        pointField end(filteredCouples.size());
-//
-//        const pointField& cellCentres = mesh_.cellCentres();
-//
-//        forAll(filteredCouples, i)
-//        {
-//            const labelPair& couple = couples[i];
-//            start[i] = cellCentres[mesh_.faceOwner()[couple.first()]];
-//            end[i] = cellCentres[mesh_.faceOwner()[couple.second()]];
-//        }
-//
-//        // Extend segments a bit
-//        {
-//            const vectorField smallVec(Foam::sqrt(SMALL)*(end-start));
-//            start -= smallVec;
-//            end += smallVec;
-//        }
-//
-//
-//        // Do test for intersections
-//        // ~~~~~~~~~~~~~~~~~~~~~~~~~
-//        labelList surface1;
-//        List<pointIndexHit> hit1;
-//        labelList region1;
-//        vectorField normal1;
-//
-//        labelList surface2;
-//        List<pointIndexHit> hit2;
-//        labelList region2;
-//        vectorField normal2;
-//
-//        surfaces_.findNearestIntersection
-//        (
-//            surfacesToBaffle,
-//            start,
-//            end,
-//
-//            surface1,
-//            hit1,
-//            region1,
-//            normal1,
-//
-//            surface2,
-//            hit2,
-//            region2,
-//            normal2
-//        );
-//
-//        forAll(testFaces, i)
-//        {
-//            if (hit1[i].hit() && hit2[i].hit())
-//            {
-//                bool createBaffle = true;
-//
-//                label faceI = couples[i].first();
-//                label patchI = mesh_.boundaryMesh().whichPatch(faceI);
-//                if (patchI != -1 && !mesh_.boundaryMesh()[patchI].coupled())
-//                {
-//                    // Check if we are a boundary face and normal of surface
-//                    // does
-//                    // not align with test vector. In this case there'd
-//                    // probably be
-//                    // a freestanding 'baffle' so we might as well not
-//                    // create it.
-//                    // Note that since it is not a proper baffle we cannot
-//                    // detect it
-//                    // afterwards so this code cannot be merged with the
-//                    // filterDuplicateFaces code.
-//                    if (mag(normal1[i]&normal2[i]) > cos(degToRad(30)))
-//                    {
-//                        // Both normals aligned
-//                        vector n = end[i]-start[i];
-//                        scalar magN = mag(n);
-//                        if (magN > VSMALL)
-//                        {
-//                            n /= magN;
-//
-//                            if (mag(normal1[i]&n) < cos(degToRad(45)))
-//                            {
-//                                Pout<< "** disabling baffling face "
-//                                    << mesh_.faceCentres()[faceI] << endl;
-//                                createBaffle = false;
-//                            }
-//                        }
-//                    }
-//                }
-//
-//
-//        }
-//XXXXXX
+        pointField start(filteredCouples.size());
+        pointField end(filteredCouples.size());
 
+        const pointField& cellCentres = mesh_.cellCentres();
+
+        forAll(filteredCouples, i)
+        {
+            const labelPair& couple = filteredCouples[i];
+            start[i] = cellCentres[mesh_.faceOwner()[couple.first()]];
+            end[i] = cellCentres[mesh_.faceOwner()[couple.second()]];
+        }
+
+        // Extend segments a bit
+        {
+            const vectorField smallVec(Foam::sqrt(SMALL)*(end-start));
+            start -= smallVec;
+            end += smallVec;
+        }
+
+
+        // Do test for intersections
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~
+        labelList surface1;
+        List<pointIndexHit> hit1;
+        labelList region1;
+        vectorField normal1;
+
+        labelList surface2;
+        List<pointIndexHit> hit2;
+        labelList region2;
+        vectorField normal2;
+
+        surfaces_.findNearestIntersection
+        (
+            identity(surfaces_.surfaces().size()),
+            start,
+            end,
+
+            surface1,
+            hit1,
+            region1,
+            normal1,
+
+            surface2,
+            hit2,
+            region2,
+            normal2
+        );
+
+        //mkDir(mesh_.time().path()/timeName());
+        //OBJstream str
+        //(
+        //    mesh_.time().path()/timeName()/"flatBaffles.obj"
+        //);
+
+        const scalar planarAngleCos = Foam::cos(degToRad(planarAngle));
+
+        label filterI = 0;
+        forAll(filteredCouples, i)
+        {
+            const labelPair& couple = filteredCouples[i];
+
+            if
+            (
+                hit1[i].hit()
+             && hit2[i].hit()
+             && (
+                    surface1[i] != surface2[i]
+                 || hit1[i].index() != hit2[i].index()
+                )
+            )
+            {
+                // Two different hits. Check angle.
+                //str.write
+                //(
+                //    linePointRef(hit1[i].hitPoint(), hit2[i].hitPoint()),
+                //    normal1[i],
+                //    normal2[i]
+                //);
+
+                if ((normal1[i]&normal2[i]) > planarAngleCos)
+                {
+                    // Both normals aligned
+                    vector n = end[i]-start[i];
+                    scalar magN = mag(n);
+                    if (magN > VSMALL)
+                    {
+                        filteredCouples[filterI++] = couple;
+                    }
+                }
+            }
+            else if (hit1[i].hit() || hit2[i].hit())
+            {
+                // Single hit. Do not include in freestanding baffles.
+            }
+        }
+
+        filteredCouples.setSize(filterI);
+
+        Info<< "freeStandingBaffles : detected "
+            << returnReduce(filterI, sumOp<label>())
+            << " planar (within " << planarAngle
+            << " degrees) free-standing baffles out of "
+            << nFiltered
+            << nl << endl;
+    }
 
     return filteredCouples;
 }
@@ -1832,15 +1870,102 @@ void Foam::meshRefinement::makeConsistentFaceIndex
 }
 
 
+void Foam::meshRefinement::handleSnapProblems
+(
+    const bool removeEdgeConnectedCells,
+    const scalarField& perpendicularAngle,
+    const dictionary& motionDict,
+    Time& runTime,
+    const labelList& globalToMasterPatch,
+    const labelList& globalToSlavePatch
+)
+{
+    Info<< nl
+        << "Introducing baffles to block off problem cells" << nl
+        << "----------------------------------------------" << nl
+        << endl;
+
+    labelList facePatch
+    (
+        markFacesOnProblemCells
+        (
+            motionDict,
+            removeEdgeConnectedCells,
+            perpendicularAngle,
+            globalToMasterPatch
+        )
+    );
+    Info<< "Analyzed problem cells in = "
+        << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
+
+    if (debug&meshRefinement::MESH)
+    {
+        faceSet problemTopo(mesh_, "problemFacesTopo", 100);
+
+        const labelList facePatchTopo
+        (
+            markFacesOnProblemCells
+            (
+                motionDict,
+                removeEdgeConnectedCells,
+                perpendicularAngle,
+                globalToMasterPatch
+            )
+        );
+        forAll(facePatchTopo, faceI)
+        {
+            if (facePatchTopo[faceI] != -1)
+            {
+                problemTopo.insert(faceI);
+            }
+        }
+        problemTopo.instance() = timeName();
+        Pout<< "Dumping " << problemTopo.size()
+            << " problem faces to " << problemTopo.objectPath() << endl;
+        problemTopo.write();
+    }
+
+    Info<< "Introducing baffles to delete problem cells." << nl << endl;
+
+    if (debug)
+    {
+        runTime++;
+    }
+
+    // Create baffles with same owner and neighbour for now.
+    createBaffles(facePatch, facePatch);
+
+    if (debug)
+    {
+        // Debug:test all is still synced across proc patches
+        checkData();
+    }
+    Info<< "Created baffles in = "
+        << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
+
+    printMeshInfo(debug, "After introducing baffles");
+
+    if (debug&meshRefinement::MESH)
+    {
+        Pout<< "Writing extra baffled mesh to time "
+            << timeName() << endl;
+        write(debug, runTime.path()/"extraBaffles");
+        Pout<< "Dumped debug data in = "
+            << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 // Split off unreachable areas of mesh.
 void Foam::meshRefinement::baffleAndSplitMesh
 (
-    const bool handleSnapProblems,
+    const bool doHandleSnapProblems,
     const bool removeEdgeConnectedCells,
     const scalarField& perpendicularAngle,
     const bool mergeFreeStanding,
+    const scalar freeStandingAngle,
     const dictionary& motionDict,
     Time& runTime,
     const labelList& globalToMasterPatch,
@@ -1902,82 +2027,92 @@ void Foam::meshRefinement::baffleAndSplitMesh
 
     // Create some additional baffles where we want surface cells removed.
 
-    if (handleSnapProblems)
+    if (doHandleSnapProblems)
     {
-        Info<< nl
-            << "Introducing baffles to block off problem cells" << nl
-            << "----------------------------------------------" << nl
-            << endl;
+        //Info<< nl
+        //    << "Introducing baffles to block off problem cells" << nl
+        //    << "----------------------------------------------" << nl
+        //    << endl;
+        //
+        //labelList facePatch
+        //(
+        //    markFacesOnProblemCells
+        //    (
+        //        motionDict,
+        //        removeEdgeConnectedCells,
+        //        perpendicularAngle,
+        //        globalToMasterPatch
+        //    )
+        //    //markFacesOnProblemCellsGeometric(motionDict)
+        //);
+        //Info<< "Analyzed problem cells in = "
+        //    << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
+        //
+        //if (debug&meshRefinement::MESH)
+        //{
+        //    faceSet problemTopo(mesh_, "problemFacesTopo", 100);
+        //
+        //    const labelList facePatchTopo
+        //    (
+        //        markFacesOnProblemCells
+        //        (
+        //            motionDict,
+        //            removeEdgeConnectedCells,
+        //            perpendicularAngle,
+        //            globalToMasterPatch
+        //        )
+        //    );
+        //    forAll(facePatchTopo, faceI)
+        //    {
+        //        if (facePatchTopo[faceI] != -1)
+        //        {
+        //            problemTopo.insert(faceI);
+        //        }
+        //    }
+        //    problemTopo.instance() = timeName();
+        //    Pout<< "Dumping " << problemTopo.size()
+        //        << " problem faces to " << problemTopo.objectPath() << endl;
+        //    problemTopo.write();
+        //}
+        //
+        //Info<< "Introducing baffles to delete problem cells." << nl << endl;
+        //
+        //if (debug)
+        //{
+        //    runTime++;
+        //}
+        //
+        //// Create baffles with same owner and neighbour for now.
+        //createBaffles(facePatch, facePatch);
+        //
+        //if (debug)
+        //{
+        //    // Debug:test all is still synced across proc patches
+        //    checkData();
+        //}
+        //Info<< "Created baffles in = "
+        //    << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
+        //
+        //printMeshInfo(debug, "After introducing baffles");
+        //
+        //if (debug&meshRefinement::MESH)
+        //{
+        //    Pout<< "Writing extra baffled mesh to time "
+        //        << timeName() << endl;
+        //    write(debug, runTime.path()/"extraBaffles");
+        //    Pout<< "Dumped debug data in = "
+        //        << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
+        //}
 
-        labelList facePatch
+        handleSnapProblems
         (
-            markFacesOnProblemCells
-            (
-                motionDict,
-                removeEdgeConnectedCells,
-                perpendicularAngle,
-                globalToMasterPatch
-            )
-            //markFacesOnProblemCellsGeometric(motionDict)
+            removeEdgeConnectedCells,
+            perpendicularAngle,
+            motionDict,
+            runTime,
+            globalToMasterPatch,
+            globalToSlavePatch
         );
-        Info<< "Analyzed problem cells in = "
-            << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
-
-        if (debug&meshRefinement::MESH)
-        {
-            faceSet problemTopo(mesh_, "problemFacesTopo", 100);
-
-            const labelList facePatchTopo
-            (
-                markFacesOnProblemCells
-                (
-                    motionDict,
-                    removeEdgeConnectedCells,
-                    perpendicularAngle,
-                    globalToMasterPatch
-                )
-            );
-            forAll(facePatchTopo, faceI)
-            {
-                if (facePatchTopo[faceI] != -1)
-                {
-                    problemTopo.insert(faceI);
-                }
-            }
-            problemTopo.instance() = timeName();
-            Pout<< "Dumping " << problemTopo.size()
-                << " problem faces to " << problemTopo.objectPath() << endl;
-            problemTopo.write();
-        }
-
-        Info<< "Introducing baffles to delete problem cells." << nl << endl;
-
-        if (debug)
-        {
-            runTime++;
-        }
-
-        // Create baffles with same owner and neighbour for now.
-        createBaffles(facePatch, facePatch);
-
-        if (debug)
-        {
-            // Debug:test all is still synced across proc patches
-            checkData();
-        }
-        Info<< "Created baffles in = "
-            << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
-
-        printMeshInfo(debug, "After introducing baffles");
-
-        if (debug&meshRefinement::MESH)
-        {
-            Pout<< "Writing extra baffled mesh to time "
-                << timeName() << endl;
-            write(debug, runTime.path()/"extraBaffles");
-            Pout<< "Dumped debug data in = "
-                << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
-        }
     }
 
 
@@ -2036,7 +2171,8 @@ void Foam::meshRefinement::baffleAndSplitMesh
                 (
                     identity(mesh_.nFaces()-mesh_.nInternalFaces())
                    +mesh_.nInternalFaces()
-                )
+                ),
+                freeStandingAngle
             )
         );
 
@@ -2052,6 +2188,18 @@ void Foam::meshRefinement::baffleAndSplitMesh
             // from them.
             mergeBaffles(couples);
 
+            // Detect any problem cells resulting from merging of baffles
+            // and delete them
+            handleSnapProblems
+            (
+                removeEdgeConnectedCells,
+                perpendicularAngle,
+                motionDict,
+                runTime,
+                globalToMasterPatch,
+                globalToSlavePatch
+            );
+
             if (debug)
             {
                 // Debug:test all is still synced across proc patches
@@ -2663,18 +2811,6 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
 
             if (surface1[i] != -1)
             {
-                //- Not allowed not to create baffle - is vital for regioning.
-                //  Have logic instead at erosion!
-                //bool createBaffle = validBaffleTopology
-                //(
-                //    faceI,
-                //    normal1[i],
-                //    normal2[i],
-                //    end[i]-start[i]
-                //);
-                //
-
-
                 // If both hit should probably choose 'nearest'
                 if
                 (
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
index 42e7983e5a..ee734bcebc 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
@@ -244,14 +244,10 @@ bool Foam::meshRefinement::markForRefine
 }
 
 
-// Calculates list of cells to refine based on intersection with feature edge.
-Foam::label Foam::meshRefinement::markFeatureRefinement
+void Foam::meshRefinement::markFeatureCellLevel
 (
     const point& keepPoint,
-    const label nAllowRefine,
-
-    labelList& refineCell,
-    label& nRefine
+    labelList& maxFeatureLevel
 ) const
 {
     // We want to refine all cells containing a feature edge.
@@ -384,7 +380,7 @@ Foam::label Foam::meshRefinement::markFeatureRefinement
 
 
     // Largest refinement level of any feature passed through
-    labelList maxFeatureLevel(mesh_.nCells(), -1);
+    maxFeatureLevel = labelList(mesh_.nCells(), -1);
 
     // Database to pass into trackedParticle::move
     trackedParticle::trackingData td(cloud, maxFeatureLevel);
@@ -467,8 +463,23 @@ Foam::label Foam::meshRefinement::markFeatureRefinement
         // Track all particles to their end position.
         cloud.move(td, GREAT);
     }
+}
 
 
+// Calculates list of cells to refine based on intersection with feature edge.
+Foam::label Foam::meshRefinement::markFeatureRefinement
+(
+    const point& keepPoint,
+    const label nAllowRefine,
+
+    labelList& refineCell,
+    label& nRefine
+) const
+{
+    // Largest refinement level of any feature passed through
+    labelList maxFeatureLevel;
+    markFeatureCellLevel(keepPoint, maxFeatureLevel);
+
     // See which cells to refine. maxFeatureLevel will hold highest level
     // of any feature edge that passed through.
 

From 08c1cac832ab94fa648c8741374119557b9421de Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Thu, 16 May 2013 15:58:15 +0100
Subject: [PATCH 19/33] ENH: Added option to (de)activate Ashford correction in
 turbulnce SA models

---
 .../LES/SpalartAllmaras/SpalartAllmaras.C     | 65 +++++++++++++++----
 .../LES/SpalartAllmaras/SpalartAllmaras.H     | 25 +++++--
 .../RAS/SpalartAllmaras/SpalartAllmaras.C     | 53 ++++++++++++---
 .../RAS/SpalartAllmaras/SpalartAllmaras.H     | 12 +++-
 .../LES/SpalartAllmaras/SpalartAllmaras.C     | 59 ++++++++++++++---
 .../LES/SpalartAllmaras/SpalartAllmaras.H     | 15 ++++-
 .../RAS/SpalartAllmaras/SpalartAllmaras.C     | 55 +++++++++++++---
 .../RAS/SpalartAllmaras/SpalartAllmaras.H     | 14 ++--
 8 files changed, 247 insertions(+), 51 deletions(-)

diff --git a/src/turbulenceModels/compressible/LES/SpalartAllmaras/SpalartAllmaras.C b/src/turbulenceModels/compressible/LES/SpalartAllmaras/SpalartAllmaras.C
index 6994f926e3..324d857a31 100644
--- a/src/turbulenceModels/compressible/LES/SpalartAllmaras/SpalartAllmaras.C
+++ b/src/turbulenceModels/compressible/LES/SpalartAllmaras/SpalartAllmaras.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
@@ -63,20 +63,51 @@ tmp<volScalarField> SpalartAllmaras::fv1() const
 
 tmp<volScalarField> SpalartAllmaras::fv2() const
 {
-    return 1.0/pow3(scalar(1) + rho()*nuTilda_/(mu()*Cv2_));
+    if (ashfordCorrection_)
+    {
+        return 1.0/pow3(scalar(1) + rho()*nuTilda_/(mu()*Cv2_));
+    }
+    else
+    {
+        const volScalarField chi("chi", rho()*nuTilda_/mu());
+        return 1.0 - chi/(1.0 + chi*fv1());
+    }
 }
 
 
 tmp<volScalarField> SpalartAllmaras::fv3() const
 {
-    volScalarField chi(rho()*nuTilda_/mu());
-    volScalarField chiByCv2((1/Cv2_)*chi);
+    if (ashfordCorrection_)
+    {
+        volScalarField chi(rho()*nuTilda_/mu());
+        volScalarField chiByCv2((1/Cv2_)*chi);
 
-    return
-        (scalar(1) + chi*fv1())
-       *(1/Cv2_)
-       *(3*(scalar(1) + chiByCv2) + sqr(chiByCv2))
-       /pow3(scalar(1) + chiByCv2);
+        return
+            (scalar(1) + chi*fv1())
+           *(1/Cv2_)
+           *(3*(scalar(1) + chiByCv2) + sqr(chiByCv2))
+           /pow3(scalar(1) + chiByCv2);
+    }
+    else
+    {
+        return tmp<volScalarField>
+        (
+            new volScalarField
+            (
+                IOobject
+                (
+                    "fv3",
+                    mesh_.time().timeName(),
+                    mesh_,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                mesh_,
+                dimensionedScalar("fv3", dimless, 1),
+                zeroGradientFvPatchScalarField::typeName
+            )
+        );
+    }
 }
 
 
@@ -222,6 +253,8 @@ SpalartAllmaras::SpalartAllmaras
         )
     ),
 
+    ashfordCorrection_(coeffDict_.lookupOrDefault("ashfordCorrection", true)),
+
     nuTilda_
     (
         IOobject
@@ -265,6 +298,11 @@ SpalartAllmaras::SpalartAllmaras
     updateSubGridScaleFields();
 
     printCoeffs();
+
+    if (ashfordCorrection_)
+    {
+        Info<< "    Employing Ashford correction" << endl;
+    }
 }
 
 
@@ -344,16 +382,19 @@ bool SpalartAllmaras::read()
     {
         sigmaNut_.readIfPresent(coeffDict());
         Prt_.readIfPresent(coeffDict());
+
         Cb1_.readIfPresent(coeffDict());
         Cb2_.readIfPresent(coeffDict());
-        Cw1_ = Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_;
-        Cw2_.readIfPresent(coeffDict());
-        Cw3_.readIfPresent(coeffDict());
         Cv1_.readIfPresent(coeffDict());
         Cv2_.readIfPresent(coeffDict());
         CDES_.readIfPresent(coeffDict());
         ck_.readIfPresent(coeffDict());
         kappa_.readIfPresent(*this);
+        Cw1_ = Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_;
+        Cw2_.readIfPresent(coeffDict());
+        Cw3_.readIfPresent(coeffDict());
+
+        ashfordCorrection_.readIfPresent("ashfordCorrection", coeffDict());
 
         return true;
     }
diff --git a/src/turbulenceModels/compressible/LES/SpalartAllmaras/SpalartAllmaras.H b/src/turbulenceModels/compressible/LES/SpalartAllmaras/SpalartAllmaras.H
index 1897578942..11d24d3388 100644
--- a/src/turbulenceModels/compressible/LES/SpalartAllmaras/SpalartAllmaras.H
+++ b/src/turbulenceModels/compressible/LES/SpalartAllmaras/SpalartAllmaras.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -30,6 +30,15 @@ Group
 Description
     SpalartAllmaras for compressible flows
 
+    Extended according to
+    \verbatim
+        "An Unstructured Grid Generation and Adaptive Solution Technique
+        for High Reynolds Number Compressible Flows"
+        G.A. Ashford,
+        Ph.D. thesis, University of Michigan, 1996.
+    \endverbatim
+    by using the optional flag \c ashfordCorrection
+
 SourceFiles
     SpalartAllmaras.C
 
@@ -77,12 +86,16 @@ class SpalartAllmaras
             dimensionedScalar Cw3_;
 
 
-    // Fields
+        //- Optional flag to activate the Ashford correction
+        Switch ashfordCorrection_;
 
-        volScalarField nuTilda_;
-        volScalarField dTilda_;
-        volScalarField muSgs_;
-        volScalarField alphaSgs_;
+
+        // Fields
+
+            volScalarField nuTilda_;
+            volScalarField dTilda_;
+            volScalarField muSgs_;
+            volScalarField alphaSgs_;
 
 
     // Private Member Functions
diff --git a/src/turbulenceModels/compressible/RAS/SpalartAllmaras/SpalartAllmaras.C b/src/turbulenceModels/compressible/RAS/SpalartAllmaras/SpalartAllmaras.C
index 6c5f3aae0c..9e5cf2d11f 100644
--- a/src/turbulenceModels/compressible/RAS/SpalartAllmaras/SpalartAllmaras.C
+++ b/src/turbulenceModels/compressible/RAS/SpalartAllmaras/SpalartAllmaras.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
@@ -63,7 +63,14 @@ tmp<volScalarField> SpalartAllmaras::fv2
     const volScalarField& fv1
 ) const
 {
-    return 1.0/pow3(scalar(1) + chi/Cv2_);
+    if (ashfordCorrection_)
+    {
+        return 1.0/pow3(scalar(1) + chi/Cv2_);
+    }
+    else
+    {
+        return 1.0 - chi/(1.0 + chi*fv1);
+    }
 }
 
 
@@ -73,13 +80,36 @@ tmp<volScalarField> SpalartAllmaras::fv3
     const volScalarField& fv1
 ) const
 {
-    const volScalarField chiByCv2((1/Cv2_)*chi);
+    if (ashfordCorrection_)
+    {
+        const volScalarField chiByCv2((1/Cv2_)*chi);
 
-    return
-        (scalar(1) + chi*fv1)
-       *(1/Cv2_)
-       *(3*(scalar(1) + chiByCv2) + sqr(chiByCv2))
-       /pow3(scalar(1) + chiByCv2);
+        return
+            (scalar(1) + chi*fv1)
+           *(1/Cv2_)
+           *(3*(scalar(1) + chiByCv2) + sqr(chiByCv2))
+           /pow3(scalar(1) + chiByCv2);
+    }
+    else
+    {
+        return tmp<volScalarField>
+        (
+            new volScalarField
+            (
+                IOobject
+                (
+                    "fv3",
+                    mesh_.time().timeName(),
+                    mesh_,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                mesh_,
+                dimensionedScalar("fv3", dimless, 1),
+                zeroGradientFvPatchScalarField::typeName
+            )
+        );
+    }
 }
 
 
@@ -252,6 +282,11 @@ SpalartAllmaras::SpalartAllmaras
     alphat_.correctBoundaryConditions();
 
     printCoeffs();
+
+    if (ashfordCorrection_)
+    {
+        Info<< "    Employing Ashford correction" << endl;
+    }
 }
 
 
@@ -372,6 +407,8 @@ bool SpalartAllmaras::read()
         Cv1_.readIfPresent(coeffDict());
         Cv2_.readIfPresent(coeffDict());
 
+        ashfordCorrection_.readIfPresent("ashfordCorrection", coeffDict());
+
         return true;
     }
     else
diff --git a/src/turbulenceModels/compressible/RAS/SpalartAllmaras/SpalartAllmaras.H b/src/turbulenceModels/compressible/RAS/SpalartAllmaras/SpalartAllmaras.H
index ea829fc88d..ce1054c482 100644
--- a/src/turbulenceModels/compressible/RAS/SpalartAllmaras/SpalartAllmaras.H
+++ b/src/turbulenceModels/compressible/RAS/SpalartAllmaras/SpalartAllmaras.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,14 +37,16 @@ Description
         P.R. Spalart,
         S.R. Allmaras,
         La Recherche Aerospatiale, No. 1, 1994, pp. 5-21.
+    \endverbatim
 
-        Extended according to:
-
+    Extended according to:
+    \verbatim
         "An Unstructured Grid Generation and Adaptive Solution Technique
         for High Reynolds Number Compressible Flows"
         G.A. Ashford,
         Ph.D. thesis, University of Michigan, 1996.
     \endverbatim
+    by using the optional flag \c ashfordCorrection
 
     The default model coefficients correspond to the following:
     \verbatim
@@ -110,6 +112,10 @@ protected:
             dimensionedScalar Cv2_;
 
 
+        //- Optional flag to activate the Ashford correction
+        Switch ashfordCorrection_;
+
+
         // Fields
 
             volScalarField nuTilda_;
diff --git a/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.C b/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.C
index 4bbf67702a..39afb02f30 100644
--- a/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.C
+++ b/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.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
@@ -60,20 +60,51 @@ tmp<volScalarField> SpalartAllmaras::fv1() const
 
 tmp<volScalarField> SpalartAllmaras::fv2() const
 {
-    return 1/pow3(scalar(1) + nuTilda_/(Cv2_*nu()));
+    if (ashfordCorrection_)
+    {
+        return 1/pow3(scalar(1) + nuTilda_/(Cv2_*nu()));
+    }
+    else
+    {
+        const volScalarField chi("chi", nuTilda_/nu());
+        return 1.0 - chi/(1.0 + chi*fv1());
+    }
 }
 
 
 tmp<volScalarField> SpalartAllmaras::fv3() const
 {
-    const volScalarField chi("chi", nuTilda_/nu());
-    const volScalarField chiByCv2(chi/Cv2_);
+    if (ashfordCorrection_)
+    {
+        const volScalarField chi("chi", nuTilda_/nu());
+        const volScalarField chiByCv2(chi/Cv2_);
 
-    return
-        (scalar(1) + chi*fv1())
-       *(1/Cv2_)
-       *(3*(scalar(1) + chiByCv2) + sqr(chiByCv2))
-       /pow3(scalar(1) + chiByCv2);
+        return
+            (scalar(1) + chi*fv1())
+           *(1/Cv2_)
+           *(3*(scalar(1) + chiByCv2) + sqr(chiByCv2))
+           /pow3(scalar(1) + chiByCv2);
+    }
+    else
+    {
+        return tmp<volScalarField>
+        (
+            new volScalarField
+            (
+                IOobject
+                (
+                    "fv3",
+                    mesh_.time().timeName(),
+                    mesh_,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                mesh_,
+                dimensionedScalar("fv3", dimless, 1),
+                zeroGradientFvPatchScalarField::typeName
+            )
+        );
+    }
 }
 
 
@@ -246,6 +277,8 @@ SpalartAllmaras::SpalartAllmaras
         )
     ),
 
+    ashfordCorrection_(coeffDict_.lookupOrDefault("ashfordCorrection", true)),
+
     y_(mesh_),
 
     nuTilda_
@@ -275,6 +308,11 @@ SpalartAllmaras::SpalartAllmaras
     )
 {
     updateSubGridScaleFields();
+
+    if (ashfordCorrection_)
+    {
+        Info<< "    Employing Ashford correction" << endl;
+    }
 }
 
 
@@ -376,6 +414,7 @@ bool SpalartAllmaras::read()
     {
         sigmaNut_.readIfPresent(coeffDict());
         kappa_.readIfPresent(*this);
+
         Cb1_.readIfPresent(coeffDict());
         Cb2_.readIfPresent(coeffDict());
         Cv1_.readIfPresent(coeffDict());
@@ -386,6 +425,8 @@ bool SpalartAllmaras::read()
         Cw2_.readIfPresent(coeffDict());
         Cw3_.readIfPresent(coeffDict());
 
+        ashfordCorrection_.readIfPresent("ashfordCorrection", coeffDict());
+
         return true;
     }
     else
diff --git a/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.H b/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.H
index fe4fe280f3..3f3b39e12b 100644
--- a/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.H
+++ b/src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -30,6 +30,15 @@ Group
 Description
     SpalartAllmaras DES (SA + LES) turbulence model for incompressible flows
 
+    Extended according to
+    \verbatim
+        "An Unstructured Grid Generation and Adaptive Solution Technique
+        for High Reynolds Number Compressible Flows"
+        G.A. Ashford,
+        Ph.D. thesis, University of Michigan, 1996.
+    \endverbatim
+    by using the optional flag \c ashfordCorrection
+
 SourceFiles
     SpalartAllmaras.C
 
@@ -90,6 +99,10 @@ protected:
             dimensionedScalar Cw3_;
 
 
+        //- Optional flag to activate the Ashford correction
+        Switch ashfordCorrection_;
+
+
         // Fields
 
             wallDist y_;
diff --git a/src/turbulenceModels/incompressible/RAS/SpalartAllmaras/SpalartAllmaras.C b/src/turbulenceModels/incompressible/RAS/SpalartAllmaras/SpalartAllmaras.C
index 76046a80e2..999978f3c9 100644
--- a/src/turbulenceModels/incompressible/RAS/SpalartAllmaras/SpalartAllmaras.C
+++ b/src/turbulenceModels/incompressible/RAS/SpalartAllmaras/SpalartAllmaras.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
@@ -61,7 +61,14 @@ tmp<volScalarField> SpalartAllmaras::fv2
     const volScalarField& fv1
 ) const
 {
-    return 1.0/pow3(scalar(1) + chi/Cv2_);
+    if (ashfordCorrection_)
+    {
+        return 1.0/pow3(scalar(1) + chi/Cv2_);
+    }
+    else
+    {
+        return 1.0 - chi/(1.0 + chi*fv1);
+    }
 }
 
 
@@ -71,13 +78,36 @@ tmp<volScalarField> SpalartAllmaras::fv3
     const volScalarField& fv1
 ) const
 {
-    const volScalarField chiByCv2((1/Cv2_)*chi);
+    if (ashfordCorrection_)
+    {
+        const volScalarField chiByCv2((1/Cv2_)*chi);
 
-    return
-        (scalar(1) + chi*fv1)
-       *(1/Cv2_)
-       *(3*(scalar(1) + chiByCv2) + sqr(chiByCv2))
-       /pow3(scalar(1) + chiByCv2);
+        return
+            (scalar(1) + chi*fv1)
+           *(1/Cv2_)
+           *(3*(scalar(1) + chiByCv2) + sqr(chiByCv2))
+           /pow3(scalar(1) + chiByCv2);
+    }
+    else
+    {
+        return tmp<volScalarField>
+        (
+            new volScalarField
+            (
+                IOobject
+                (
+                    "fv3",
+                    mesh_.time().timeName(),
+                    mesh_,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                mesh_,
+                dimensionedScalar("fv3", dimless, 1),
+                zeroGradientFvPatchScalarField::typeName
+            )
+        );
+    }
 }
 
 
@@ -195,6 +225,8 @@ SpalartAllmaras::SpalartAllmaras
         )
     ),
 
+    ashfordCorrection_(coeffDict_.lookupOrDefault("ashfordCorrection", true)),
+
     nuTilda_
     (
         IOobject
@@ -224,6 +256,11 @@ SpalartAllmaras::SpalartAllmaras
     d_(mesh_)
 {
     printCoeffs();
+
+    if (ashfordCorrection_)
+    {
+        Info<< "    Employing Ashford correction" << endl;
+    }
 }
 
 
@@ -368,6 +405,8 @@ bool SpalartAllmaras::read()
         Cv1_.readIfPresent(coeffDict());
         Cv2_.readIfPresent(coeffDict());
 
+        ashfordCorrection_.readIfPresent("ashfordCorrection", coeffDict());
+
         return true;
     }
     else
diff --git a/src/turbulenceModels/incompressible/RAS/SpalartAllmaras/SpalartAllmaras.H b/src/turbulenceModels/incompressible/RAS/SpalartAllmaras/SpalartAllmaras.H
index d7479c2dd5..597b0859c4 100644
--- a/src/turbulenceModels/incompressible/RAS/SpalartAllmaras/SpalartAllmaras.H
+++ b/src/turbulenceModels/incompressible/RAS/SpalartAllmaras/SpalartAllmaras.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,14 +37,16 @@ Description
         P.R. Spalart,
         S.R. Allmaras,
         La Recherche Aerospatiale, No. 1, 1994, pp. 5-21.
+    \endverbatim
 
-        Extended according to:
-
+    Extended according to
+    \verbatim
         "An Unstructured Grid Generation and Adaptive Solution Technique
         for High Reynolds Number Compressible Flows"
         G.A. Ashford,
         Ph.D. thesis, University of Michigan, 1996.
     \endverbatim
+    using the optional flag \c ashfordCorrection
 
     The default model coefficients correspond to the following:
     \verbatim
@@ -82,7 +84,7 @@ namespace RASModels
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class SpalartAllmaras Declaration
+                       Class SpalartAllmaras Declaration
 \*---------------------------------------------------------------------------*/
 
 class SpalartAllmaras
@@ -108,6 +110,10 @@ protected:
             dimensionedScalar Cv2_;
 
 
+        //- Optional flag to activate the Ashford correction
+        Switch ashfordCorrection_;
+
+
         // Fields
 
             volScalarField nuTilda_;

From 153e0c690b8bb1c2798b4d296264c721bded672d Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Thu, 16 May 2013 20:28:47 +0100
Subject: [PATCH 20/33] ENH: GAMGAgglomeration: better printing

---
 .../GAMGAgglomeration/GAMGAgglomeration.C     | 22 +++++++++++++++++--
 1 file changed, 20 insertions(+), 2 deletions(-)

diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.C
index a53a5fe9ab..0b950fd433 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.C
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.C
@@ -79,8 +79,9 @@ void Foam::GAMGAgglomeration::compactLevels(const label nCreatedLevels)
                 << nl;
 
             Info<< setw(40) << "nCells"
+                << setw(24) << "nFaces/nCells"
                 << setw(24) << "nInterfaces"
-                << setw(24) << "Ratio" << nl
+                << setw(24) << "nIntFaces/nCells" << nl
                 << setw(8) << "Level"
                 << setw(8) << "nProcs"
                 << "        "
@@ -89,6 +90,9 @@ void Foam::GAMGAgglomeration::compactLevels(const label nCreatedLevels)
                 << "        "
                 << setw(8) << "avg"
                 << setw(8) << "max"
+                << "        "
+                << setw(8) << "avg"
+                << setw(8) << "max"
                 << "            " << setw(4) << "avg"
                 << "    " << setw(4) << "max"
                 << nl
@@ -108,6 +112,7 @@ void Foam::GAMGAgglomeration::compactLevels(const label nCreatedLevels)
             {
                 label nProcs = 0;
                 label nCells = 0;
+                scalar faceCellRatio = 0;
                 label nInterfaces = 0;
                 label nIntFaces = 0;
                 scalar ratio = 0.0;
@@ -118,6 +123,8 @@ void Foam::GAMGAgglomeration::compactLevels(const label nCreatedLevels)
 
                     const lduMesh& fineMesh = meshLevel(levelI);
                     nCells = fineMesh.lduAddr().size();
+                    faceCellRatio =
+                        scalar(fineMesh.lduAddr().lowerAddr().size())/nCells;
 
                     const lduInterfacePtrsList interfaces =
                         fineMesh.interfaces();
@@ -137,6 +144,11 @@ void Foam::GAMGAgglomeration::compactLevels(const label nCreatedLevels)
                 label maxNCells = returnReduce(nCells, maxOp<label>());
                 label totNCells = returnReduce(nCells, sumOp<label>());
 
+                scalar maxFaceCellRatio =
+                    returnReduce(faceCellRatio, maxOp<scalar>());
+                scalar totFaceCellRatio =
+                    returnReduce(faceCellRatio, sumOp<scalar>());
+
                 label maxNInt = returnReduce(nInterfaces, maxOp<label>());
                 label totNInt = returnReduce(nInterfaces, sumOp<label>());
 
@@ -145,9 +157,15 @@ void Foam::GAMGAgglomeration::compactLevels(const label nCreatedLevels)
 
                 Info<< setw(8) << levelI
                     << setw(8) << totNprocs
+
                     << setw(16) << totNCells/totNprocs
                     << setw(8) << maxNCells
-                    << setw(16) << totNInt/totNprocs
+
+                    << "        "
+                    << setw(8) << setprecision(4) << totFaceCellRatio/totNprocs
+                    << setw(8) << setprecision(4) << maxFaceCellRatio
+                    << "        "
+                    << setw(8) << totNInt/totNprocs
                     << setw(8) << maxNInt
                     << "        "
                     << setw(8) << setprecision(4) << totRatio/totNprocs

From 941e1778dc581440e60df250359421308904af8e Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Thu, 16 May 2013 20:29:34 +0100
Subject: [PATCH 21/33] BUG: timeVaryingMappedFixedValueFvPatchField: clone
 only if existing

---
 .../timeVaryingMappedFixedValueFvPatchField.C | 21 ++++++++++++++++---
 1 file changed, 18 insertions(+), 3 deletions(-)

diff --git a/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C
index 6ea997f3fd..a09f526123 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C
@@ -79,7 +79,12 @@ timeVaryingMappedFixedValueFvPatchField
     endSampleTime_(-1),
     endSampledValues_(0),
     endAverage_(pTraits<Type>::zero),
-    offset_(ptf.offset_().clone().ptr())
+    offset_
+    (
+        ptf.offset_.valid()
+      ? ptf.offset_().clone().ptr()
+      : NULL
+    )
 {}
 
 
@@ -138,7 +143,12 @@ timeVaryingMappedFixedValueFvPatchField
     endSampleTime_(ptf.endSampleTime_),
     endSampledValues_(ptf.endSampledValues_),
     endAverage_(ptf.endAverage_),
-    offset_(ptf.offset_().clone().ptr())
+    offset_
+    (
+        ptf.offset_.valid()
+      ? ptf.offset_().clone().ptr()
+      : NULL
+    )
 {}
 
 
@@ -163,7 +173,12 @@ timeVaryingMappedFixedValueFvPatchField
     endSampleTime_(ptf.endSampleTime_),
     endSampledValues_(ptf.endSampledValues_),
     endAverage_(ptf.endAverage_),
-    offset_(ptf.offset_().clone().ptr())
+    offset_
+    (
+        ptf.offset_.valid()
+      ? ptf.offset_().clone().ptr()
+      : NULL
+    )
 {}
 
 

From 02f727c0eb4864a60c122747c4a0e474fd441081 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Fri, 17 May 2013 11:19:58 +0100
Subject: [PATCH 22/33] ENH: GAMGAgglomeration: corrected message

---
 .../GAMGAgglomeration/GAMGAgglomeration.C                  | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.C
index 0b950fd433..75af8eb20a 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.C
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.C
@@ -279,9 +279,9 @@ const Foam::GAMGAgglomeration& Foam::GAMGAgglomeration::New
                 "(const lduMesh& mesh, const dictionary& controlDict)"
             )   << "Unknown GAMGAgglomeration type "
                 << agglomeratorType << ".\n"
-                << "Valid algebraic GAMGAgglomeration types are :"
+                << "Valid matrix GAMGAgglomeration types are :"
                 << lduMatrixConstructorTablePtr_->sortedToc() << endl
-                << "Valid algebraic GAMGAgglomeration types are :"
+                << "Valid geometric GAMGAgglomeration types are :"
                 << lduMeshConstructorTablePtr_->sortedToc()
                 << exit(FatalError);
         }
@@ -374,7 +374,8 @@ Foam::autoPtr<Foam::GAMGAgglomeration> Foam::GAMGAgglomeration::New
         FatalErrorIn
         (
             "GAMGAgglomeration::New"
-            "(const lduMesh& mesh, const dictionary& controlDict)"
+            "(const lduMesh& mesh, const scalarField&"
+            ", const vectorField&, const dictionary& controlDict)"
         )   << "Unknown GAMGAgglomeration type "
             << agglomeratorType << ".\n"
             << "Valid geometric GAMGAgglomeration types are :"

From 6b656e63f787d1e928bc108ef5e08e467a47d478 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Fri, 17 May 2013 11:20:26 +0100
Subject: [PATCH 23/33] BUG: MGridGenAgglomerate: make sure face data is
 correct size

---
 .../MGridGenGamgAgglomeration/MGridGenGAMGAgglomeration.C   | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/src/fvAgglomerationMethods/MGridGenGamgAgglomeration/MGridGenGAMGAgglomeration.C b/src/fvAgglomerationMethods/MGridGenGamgAgglomeration/MGridGenGAMGAgglomeration.C
index 46bcbc8953..59c45bc11b 100644
--- a/src/fvAgglomerationMethods/MGridGenGamgAgglomeration/MGridGenGAMGAgglomeration.C
+++ b/src/fvAgglomerationMethods/MGridGenGamgAgglomeration/MGridGenGAMGAgglomeration.C
@@ -59,7 +59,11 @@ Foam::MGridGenGAMGAgglomeration::MGridGenGAMGAgglomeration
 
     // Start geometric agglomeration from the cell volumes and areas of the mesh
     scalarField* VPtr = const_cast<scalarField*>(&fvMesh_.cellVolumes());
-    vectorField* SfPtr = const_cast<vectorField*>(&fvMesh_.faceAreas());
+    SubField<vector> Sf(fvMesh_.faceAreas(), fvMesh_.nInternalFaces());
+    vectorField* SfPtr = const_cast<vectorField*>
+    (
+        &Sf.operator const vectorField&()
+    );
 
     // Create the boundary area cell field
     scalarField* SbPtr(new scalarField(fvMesh_.nCells(), 0));

From 099ad24c21141c2e2484aaf39c9e2b1dceb59511 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Mon, 20 May 2013 09:18:23 +0100
Subject: [PATCH 24/33] STYLE: structuredRenumber: indentation

---
 .../IO/writeRegisteredObject/writeRegisteredObject.H          | 4 ++--
 .../renumberMethods/structuredRenumber/structuredRenumber.C   | 2 +-
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/postProcessing/functionObjects/IO/writeRegisteredObject/writeRegisteredObject.H b/src/postProcessing/functionObjects/IO/writeRegisteredObject/writeRegisteredObject.H
index f4bd97bd36..014d4ce90c 100644
--- a/src/postProcessing/functionObjects/IO/writeRegisteredObject/writeRegisteredObject.H
+++ b/src/postProcessing/functionObjects/IO/writeRegisteredObject/writeRegisteredObject.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -38,7 +38,7 @@ Description
         type        writeRegisteredObject;
         functionObjectLibs ("libIOFunctionObjects.so");
         ...
-        objectNames (obj1 obj2);   
+        objectNames (obj1 obj2);
     }
     \endverbatim
 
diff --git a/src/renumber/renumberMethods/structuredRenumber/structuredRenumber.C b/src/renumber/renumberMethods/structuredRenumber/structuredRenumber.C
index b0a13a8399..289adba5c9 100644
--- a/src/renumber/renumberMethods/structuredRenumber/structuredRenumber.C
+++ b/src/renumber/renumberMethods/structuredRenumber/structuredRenumber.C
@@ -202,7 +202,7 @@ Foam::labelList Foam::structuredRenumber::renumber
             if (depthFirst_)
             {
                 orderedToOld[nLayers*cellData[cellI].data()+layerI] = cellI;
-                                }
+            }
             else
             {
                 orderedToOld[cellData[cellI].data()+nLayers*layerI] = cellI;

From 88ceb576544eaaffb552dd9082c434a9cfaa3464 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Mon, 20 May 2013 09:29:51 +0100
Subject: [PATCH 25/33] BUG: foamToEnsight: 4 digit files

---
 .../postProcessing/dataConversion/foamToEnsight/foamToEnsight.C | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/applications/utilities/postProcessing/dataConversion/foamToEnsight/foamToEnsight.C b/applications/utilities/postProcessing/dataConversion/foamToEnsight/foamToEnsight.C
index edb29d3998..d50bf96448 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToEnsight/foamToEnsight.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToEnsight/foamToEnsight.C
@@ -240,7 +240,7 @@ int main(int argc, char *argv[])
     wordHashSet allCloudNames;
     if (Pstream::master())
     {
-        word geomFileName = prepend + "000";
+        word geomFileName = prepend + "0000";
 
         // test pre check variable if there is a moving mesh
         if (meshMoving)

From 45955e958d4f16e788a42bcab2a9a2c0680896d7 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Mon, 20 May 2013 10:45:54 +0100
Subject: [PATCH 26/33] ENH: timeVaryingMapped: do not reread table

---
 .../timeVaryingMappedFixedValueFvPatchField.C             | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C
index a09f526123..df9c571397 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C
@@ -46,6 +46,7 @@ timeVaryingMappedFixedValueFvPatchField
     fieldTableName_(iF.name()),
     setAverage_(false),
     perturb_(0),
+    mapperPtr_(NULL),
     sampleTimes_(0),
     startSampleTime_(-1),
     startSampledValues_(0),
@@ -71,7 +72,7 @@ timeVaryingMappedFixedValueFvPatchField
     fieldTableName_(ptf.fieldTableName_),
     setAverage_(ptf.setAverage_),
     perturb_(ptf.perturb_),
-    mapperPtr_(ptf.mapperPtr_),
+    mapperPtr_(NULL),
     sampleTimes_(0),
     startSampleTime_(-1),
     startSampledValues_(0),
@@ -135,7 +136,7 @@ timeVaryingMappedFixedValueFvPatchField
     fieldTableName_(ptf.fieldTableName_),
     setAverage_(ptf.setAverage_),
     perturb_(ptf.perturb_),
-    mapperPtr_(ptf.mapperPtr_),
+    mapperPtr_(NULL),
     sampleTimes_(ptf.sampleTimes_),
     startSampleTime_(ptf.startSampleTime_),
     startSampledValues_(ptf.startSampledValues_),
@@ -152,7 +153,6 @@ timeVaryingMappedFixedValueFvPatchField
 {}
 
 
-
 template<class Type>
 timeVaryingMappedFixedValueFvPatchField<Type>::
 timeVaryingMappedFixedValueFvPatchField
@@ -165,7 +165,7 @@ timeVaryingMappedFixedValueFvPatchField
     fieldTableName_(ptf.fieldTableName_),
     setAverage_(ptf.setAverage_),
     perturb_(ptf.perturb_),
-    mapperPtr_(ptf.mapperPtr_),
+    mapperPtr_(NULL),
     sampleTimes_(ptf.sampleTimes_),
     startSampleTime_(ptf.startSampleTime_),
     startSampledValues_(ptf.startSampledValues_),

From 3cdc5d40854c3c8fdf087f663e0f5c498129c922 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 21 May 2013 14:20:41 +0100
Subject: [PATCH 27/33] STYLE: processorPatch: remove unused comm argument

---
 .../constraint/processor/processorPolyPatch.C      | 14 --------------
 .../constraint/processor/processorPolyPatch.H      |  1 -
 .../processorCyclic/processorCyclicPolyPatch.C     |  2 --
 .../processorCyclic/processorCyclicPolyPatch.H     |  1 -
 4 files changed, 18 deletions(-)

diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.C
index 515a78bbb8..8a29a7e06d 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.C
@@ -54,14 +54,12 @@ Foam::processorPolyPatch::processorPolyPatch
     const label start,
     const label index,
     const polyBoundaryMesh& bm,
-    const label comm,
     const int myProcNo,
     const int neighbProcNo,
     const transformType transform
 )
 :
     coupledPolyPatch(name, size, start, index, bm, typeName, transform),
-//    comm_(comm),
     myProcNo_(myProcNo),
     neighbProcNo_(neighbProcNo),
     neighbFaceCentres_(),
@@ -80,10 +78,6 @@ Foam::processorPolyPatch::processorPolyPatch
 )
 :
     coupledPolyPatch(name, dict, index, bm, patchType),
-//    comm_
-//    (
-//        dict.lookupOrDefault("communicator", UPstream::worldComm)
-//    ),
     myProcNo_(readLabel(dict.lookup("myProcNo"))),
     neighbProcNo_(readLabel(dict.lookup("neighbProcNo"))),
     neighbFaceCentres_(),
@@ -99,7 +93,6 @@ Foam::processorPolyPatch::processorPolyPatch
 )
 :
     coupledPolyPatch(pp, bm),
-//    comm_(pp.comm_),
     myProcNo_(pp.myProcNo_),
     neighbProcNo_(pp.neighbProcNo_),
     neighbFaceCentres_(),
@@ -118,7 +111,6 @@ Foam::processorPolyPatch::processorPolyPatch
 )
 :
     coupledPolyPatch(pp, bm, index, newSize, newStart),
-//    comm_(pp.comm_),
     myProcNo_(pp.myProcNo_),
     neighbProcNo_(pp.neighbProcNo_),
     neighbFaceCentres_(),
@@ -137,7 +129,6 @@ Foam::processorPolyPatch::processorPolyPatch
 )
 :
     coupledPolyPatch(pp, bm, index, mapAddressing, newStart),
-//    comm_(pp.comm_),
     myProcNo_(pp.myProcNo_),
     neighbProcNo_(pp.neighbProcNo_),
     neighbFaceCentres_(),
@@ -1091,11 +1082,6 @@ bool Foam::processorPolyPatch::order
 void Foam::processorPolyPatch::write(Ostream& os) const
 {
     coupledPolyPatch::write(os);
-//    if (comm_ != UPstream::worldComm)
-//    {
-//        os.writeKeyword("communicator") << comm_
-//            << token::END_STATEMENT << nl;
-//    }
     os.writeKeyword("myProcNo") << myProcNo_
         << token::END_STATEMENT << nl;
     os.writeKeyword("neighbProcNo") << neighbProcNo_
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.H b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.H
index c5a6226dc8..6a3d202a42 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.H
@@ -134,7 +134,6 @@ public:
             const label start,
             const label index,
             const polyBoundaryMesh& bm,
-            const label comm,
             const int myProcNo,
             const int neighbProcNo,
             const transformType transform = UNKNOWN  // transformation type
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processorCyclic/processorCyclicPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processorCyclic/processorCyclicPolyPatch.C
index fd8ed4d30f..137339e631 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processorCyclic/processorCyclicPolyPatch.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processorCyclic/processorCyclicPolyPatch.C
@@ -46,7 +46,6 @@ Foam::processorCyclicPolyPatch::processorCyclicPolyPatch
     const label start,
     const label index,
     const polyBoundaryMesh& bm,
-    const label comm,
     const int myProcNo,
     const int neighbProcNo,
     const word& referPatchName,
@@ -60,7 +59,6 @@ Foam::processorCyclicPolyPatch::processorCyclicPolyPatch
         start,
         index,
         bm,
-        comm,
         myProcNo,
         neighbProcNo,
         transform
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processorCyclic/processorCyclicPolyPatch.H b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processorCyclic/processorCyclicPolyPatch.H
index 7f749b3286..76a833ba61 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processorCyclic/processorCyclicPolyPatch.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processorCyclic/processorCyclicPolyPatch.H
@@ -119,7 +119,6 @@ public:
             const label start,
             const label index,
             const polyBoundaryMesh& bm,
-            const label comm,
             const int myProcNo,
             const int neighbProcNo,
             const word& referPatchName,

From f115890dbad2d2be701a4009cfbb730fc33eaa96 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 21 May 2013 14:22:43 +0100
Subject: [PATCH 28/33] ENH: autoPtr: added reUse flag, improved error message

---
 src/OpenFOAM/memory/autoPtr/autoPtr.H         |  6 +++-
 src/OpenFOAM/memory/autoPtr/autoPtrI.H        | 30 +++++++++++++++----
 src/OpenFOAM/memory/tmp/tmpI.H                | 15 ++++++----
 .../thermoSingleLayer/thermoSingleLayerI.H    |  3 +-
 4 files changed, 42 insertions(+), 12 deletions(-)

diff --git a/src/OpenFOAM/memory/autoPtr/autoPtr.H b/src/OpenFOAM/memory/autoPtr/autoPtr.H
index 7732ecc1cf..0d5944128e 100644
--- a/src/OpenFOAM/memory/autoPtr/autoPtr.H
+++ b/src/OpenFOAM/memory/autoPtr/autoPtr.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
@@ -67,6 +67,10 @@ public:
         //  setting the arguments pointer to NULL
         inline autoPtr(const autoPtr<T>&);
 
+        //- Construct either by transfering pointer or cloning. Should
+        //  only be called with type that supports cloning.
+        inline autoPtr(const autoPtr<T>&, const bool reUse);
+
 
     //- Destructor, delete object if pointer is not NULL
     inline ~autoPtr();
diff --git a/src/OpenFOAM/memory/autoPtr/autoPtrI.H b/src/OpenFOAM/memory/autoPtr/autoPtrI.H
index df7ea51d1a..55e032fcd4 100644
--- a/src/OpenFOAM/memory/autoPtr/autoPtrI.H
+++ b/src/OpenFOAM/memory/autoPtr/autoPtrI.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
@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "error.H"
+#include <typeinfo>
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -43,6 +44,21 @@ inline Foam::autoPtr<T>::autoPtr(const autoPtr<T>& ap)
 }
 
 
+template<class T>
+inline Foam::autoPtr<T>::autoPtr(const autoPtr<T>& ap, const bool reUse)
+{
+    if (reUse)
+    {
+        ptr_ = ap.ptr_;
+        ap.ptr_ = 0;
+    }
+    else
+    {
+        ptr_ = ap().clone().ptr();
+    }
+}
+
+
 template<class T>
 inline Foam::autoPtr<T>::~autoPtr()
 {
@@ -81,7 +97,8 @@ inline void Foam::autoPtr<T>::set(T* p)
     if (ptr_)
     {
         FatalErrorIn("void Foam::autoPtr<T>::set(T*)")
-            << "object already allocated"
+            << "object of type " << typeid(T).name()
+            << " already allocated"
             << abort(FatalError);
     }
 
@@ -116,7 +133,8 @@ inline T& Foam::autoPtr<T>::operator()()
     if (!ptr_)
     {
         FatalErrorIn("T& Foam::autoPtr<T>::operator()()")
-            << "object is not allocated"
+            << "object of type " << typeid(T).name()
+            << " is not allocated"
             << abort(FatalError);
     }
 
@@ -130,7 +148,8 @@ inline const T& Foam::autoPtr<T>::operator()() const
     if (!ptr_)
     {
         FatalErrorIn("const T& Foam::autoPtr<T>::operator()() const")
-            << "object is not allocated"
+            << "object of type " << typeid(T).name()
+            << " is not allocated"
             << abort(FatalError);
     }
 
@@ -151,7 +170,8 @@ inline T* Foam::autoPtr<T>::operator->()
     if (!ptr_)
     {
         FatalErrorIn("Foam::autoPtr<T>::operator->()")
-            << "object is not allocated"
+            << "object of type " << typeid(T).name()
+            << " is not allocated"
             << abort(FatalError);
     }
 
diff --git a/src/OpenFOAM/memory/tmp/tmpI.H b/src/OpenFOAM/memory/tmp/tmpI.H
index c3a51f3564..ad34bb4b17 100644
--- a/src/OpenFOAM/memory/tmp/tmpI.H
+++ b/src/OpenFOAM/memory/tmp/tmpI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "error.H"
+#include <typeinfo>
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -62,6 +63,7 @@ inline Foam::tmp<T>::tmp(const tmp<T>& t)
         {
             FatalErrorIn("Foam::tmp<T>::tmp(const tmp<T>&)")
                 << "attempted copy of a deallocated temporary"
+                << " of type " << typeid(T).name()
                 << abort(FatalError);
         }
     }
@@ -93,6 +95,7 @@ inline Foam::tmp<T>::tmp(const tmp<T>& t, bool allowTransfer)
                 (
                     "Foam::tmp<T>::tmp(const tmp<T>&, bool allowTransfer)"
                 )   << "attempted copy of a deallocated temporary"
+                    << " of type " << typeid(T).name()
                     << abort(FatalError);
             }
         }
@@ -149,7 +152,7 @@ inline T* Foam::tmp<T>::ptr() const
          if (!ptr_)
          {
              FatalErrorIn("Foam::tmp<T>::ptr() const")
-                 << "temporary deallocated"
+                 << "temporary of type " << typeid(T).name() << " deallocated"
                  << abort(FatalError);
          }
 
@@ -188,7 +191,7 @@ inline T& Foam::tmp<T>::operator()()
         if (!ptr_)
         {
             FatalErrorIn("T& Foam::tmp<T>::operator()()")
-                << "temporary deallocated"
+                << "temporary of type " << typeid(T).name() << " deallocated"
                 << abort(FatalError);
         }
 
@@ -217,7 +220,7 @@ inline const T& Foam::tmp<T>::operator()() const
         if (!ptr_)
         {
             FatalErrorIn("const T& Foam::tmp<T>::operator()() const")
-                << "temporary deallocated"
+                << "temporary of type " << typeid(T).name() << " deallocated"
                 << abort(FatalError);
         }
 
@@ -245,7 +248,7 @@ inline T* Foam::tmp<T>::operator->()
          if (!ptr_)
          {
              FatalErrorIn("Foam::tmp<T>::operator->()")
-                 << "temporary deallocated"
+                 << "temporary of type " << typeid(T).name() << " deallocated"
                  << abort(FatalError);
          }
 
@@ -294,6 +297,7 @@ inline void Foam::tmp<T>::operator=(const tmp<T>& t)
         {
             FatalErrorIn("Foam::tmp<T>::operator=(const tmp<T>&)")
                 << "attempted copy of a deallocated temporary"
+                << " of type " << typeid(T).name()
                 << abort(FatalError);
         }
     }
@@ -301,6 +305,7 @@ inline void Foam::tmp<T>::operator=(const tmp<T>& t)
     {
         FatalErrorIn("Foam::tmp<T>::operator=(const tmp<T>&)")
             << "attempted to assign to a const reference to constant object"
+            << " of type " << typeid(T).name()
             << abort(FatalError);
     }
 }
diff --git a/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayerI.H b/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayerI.H
index 659a315f49..0905176adf 100644
--- a/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayerI.H
+++ b/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayerI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,6 +25,7 @@ License
 
 #include "thermoSingleLayer.H"
 #include "heatTransferModel.H"
+#include "filmRadiationModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

From 5f034c923520749de525149046a73495bbd2643f Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 21 May 2013 14:24:44 +0100
Subject: [PATCH 29/33] ENH: processorPolyPatch: remove unused argument

---
 .../utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C | 1 -
 src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C             | 2 --
 2 files changed, 3 deletions(-)

diff --git a/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C b/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C
index f4ce2cc95d..f2c091e4ad 100644
--- a/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C
+++ b/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C
@@ -493,7 +493,6 @@ int main(int argc, char *argv[])
                         mesh.nFaces(),      // start
                         patchI,             // index
                         mesh.boundaryMesh(),// polyBoundaryMesh
-                        Pstream::worldComm, // communicator
                         Pstream::myProcNo(),// myProcNo
                         nbrProcI            // neighbProcNo
                     )
diff --git a/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C b/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C
index f6e8bc02dc..61cc9986a8 100644
--- a/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C
+++ b/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C
@@ -945,7 +945,6 @@ void Foam::fvMeshDistribute::addProcPatches
                         mesh_.nFaces(),
                         mesh_.boundaryMesh().size(),
                         mesh_.boundaryMesh(),
-                        Pstream::worldComm,
                         Pstream::myProcNo(),
                         nbrProc[bFaceI]
                     );
@@ -989,7 +988,6 @@ void Foam::fvMeshDistribute::addProcPatches
                         mesh_.nFaces(),
                         mesh_.boundaryMesh().size(),
                         mesh_.boundaryMesh(),
-                        Pstream::worldComm,
                         Pstream::myProcNo(),
                         nbrProc[bFaceI],
                         cycName,

From 975602a937f14909a2593f9c6ec23d85ef285f6f Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 21 May 2013 15:31:07 +0100
Subject: [PATCH 30/33] ENH: createBafflesDict: added comment

---
 .../manipulation/createBaffles/createBafflesDict   | 14 ++++++++++++++
 1 file changed, 14 insertions(+)

diff --git a/applications/utilities/mesh/manipulation/createBaffles/createBafflesDict b/applications/utilities/mesh/manipulation/createBaffles/createBafflesDict
index b392f07642..ef2eb96227 100644
--- a/applications/utilities/mesh/manipulation/createBaffles/createBafflesDict
+++ b/applications/utilities/mesh/manipulation/createBaffles/createBafflesDict
@@ -14,6 +14,20 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+// Sample for creating baffles:
+// - usually converting internal faces into two boundary faces
+// - or converting boundary faces into a boundary face
+//   (internalFacesOnly=false)(though should use really createPatch
+//    to do this)
+// - can also create duplicate (overlapping) sets of baffles:
+//      - internalFacesOnly = false
+//      - have 4 entries in patches:
+//          - master
+//          - slave
+//          - additional master
+//          - additional slave
+
+
 // Whether to convert internal faces only (so leave boundary faces intact).
 // This is only relevant if your face selection type can pick up boundary
 // faces.

From 2e5227a6f8bd8e2b7dbaa3f585fc32f386a1d94a Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 21 May 2013 15:33:17 +0100
Subject: [PATCH 31/33] ENH: Test-ExtendedStencil: more examples

---
 .../extendedStencil/Test-ExtendedStencil.C    | 60 ++++++++++++++++++-
 1 file changed, 59 insertions(+), 1 deletion(-)

diff --git a/applications/test/extendedStencil/Test-ExtendedStencil.C b/applications/test/extendedStencil/Test-ExtendedStencil.C
index 72c1afb2c1..0e50f8710f 100644
--- a/applications/test/extendedStencil/Test-ExtendedStencil.C
+++ b/applications/test/extendedStencil/Test-ExtendedStencil.C
@@ -56,6 +56,8 @@ Description
 //#include "centredCFCFaceToCellStencilObject.H"
 
 #include "centredCECCellToCellStencilObject.H"
+#include "centredCFCCellToCellStencilObject.H"
+#include "centredCPCCellToCellStencilObject.H"
 
 using namespace Foam;
 
@@ -460,7 +462,63 @@ int main(int argc, char *argv[])
         {
             writeStencilOBJ
             (
-                runTime.path()/"centredCell" + Foam::name(cellI) + ".obj",
+                runTime.path()/"centredCECCell" + Foam::name(cellI) + ".obj",
+                mesh.cellCentres()[cellI],
+                stencilPoints[cellI]
+            );
+        }
+    }
+    {
+        const extendedCentredCellToCellStencil& addressing =
+        centredCPCCellToCellStencilObject::New
+        (
+            mesh
+        );
+
+        Info<< "cellCellCell:" << endl;
+        writeStencilStats(addressing.stencil());
+
+        // Collect stencil cell centres
+        List<List<point> > stencilPoints(mesh.nCells());
+        addressing.collectData
+        (
+            mesh.C(),
+            stencilPoints
+        );
+
+        forAll(stencilPoints, cellI)
+        {
+            writeStencilOBJ
+            (
+                runTime.path()/"centredCPCCell" + Foam::name(cellI) + ".obj",
+                mesh.cellCentres()[cellI],
+                stencilPoints[cellI]
+            );
+        }
+    }
+    {
+        const extendedCentredCellToCellStencil& addressing =
+        centredCFCCellToCellStencilObject::New
+        (
+            mesh
+        );
+
+        Info<< "cellCellCell:" << endl;
+        writeStencilStats(addressing.stencil());
+
+        // Collect stencil cell centres
+        List<List<point> > stencilPoints(mesh.nCells());
+        addressing.collectData
+        (
+            mesh.C(),
+            stencilPoints
+        );
+
+        forAll(stencilPoints, cellI)
+        {
+            writeStencilOBJ
+            (
+                runTime.path()/"centredCFCCell" + Foam::name(cellI) + ".obj",
                 mesh.cellCentres()[cellI],
                 stencilPoints[cellI]
             );

From bfc3b12920bbdb65c3cf942a5168ba1e17a58720 Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Wed, 22 May 2013 09:23:24 +0100
Subject: [PATCH 32/33] ENH: faceSource|cellSource function objects: added
 warning if field not processed

---
 .../field/fieldValues/cellSource/cellSource.C | 23 ++++++++++++++-----
 .../field/fieldValues/faceSource/faceSource.C | 21 +++++++++++++----
 2 files changed, 33 insertions(+), 11 deletions(-)

diff --git a/src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSource.C b/src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSource.C
index 3d7977c488..595c16084d 100644
--- a/src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSource.C
+++ b/src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSource.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
@@ -222,11 +222,22 @@ void Foam::fieldValues::cellSource::write()
 
         forAll(fields_, i)
         {
-            writeValues<scalar>(fields_[i]);
-            writeValues<vector>(fields_[i]);
-            writeValues<sphericalTensor>(fields_[i]);
-            writeValues<symmTensor>(fields_[i]);
-            writeValues<tensor>(fields_[i]);
+            const word& fieldName = fields_[i];
+            bool processed = false;
+
+            processed = processed || writeValues<scalar>(fieldName);
+            processed = processed || writeValues<vector>(fieldName);
+            processed = processed || writeValues<sphericalTensor>(fieldName);
+            processed = processed || writeValues<symmTensor>(fieldName);
+            processed = processed || writeValues<tensor>(fieldName);
+
+            if (!processed)
+            {
+                WarningIn("void Foam::fieldValues::cellSource::write()")
+                    << "Requested field " << fieldName
+                    << " not found in database and not processed"
+                    << endl;
+            }
         }
 
         if (Pstream::master())
diff --git a/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.C b/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.C
index dd54de718a..ee3b626562 100644
--- a/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.C
+++ b/src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.C
@@ -630,11 +630,22 @@ void Foam::fieldValues::faceSource::write()
 
         forAll(fields_, i)
         {
-            writeValues<scalar>(fields_[i]);
-            writeValues<vector>(fields_[i]);
-            writeValues<sphericalTensor>(fields_[i]);
-            writeValues<symmTensor>(fields_[i]);
-            writeValues<tensor>(fields_[i]);
+            const word& fieldName = fields_[i];
+            bool processed = false;
+
+            processed = processed || writeValues<scalar>(fieldName);
+            processed = processed || writeValues<vector>(fieldName);
+            processed = processed || writeValues<sphericalTensor>(fieldName);
+            processed = processed || writeValues<symmTensor>(fieldName);
+            processed = processed || writeValues<tensor>(fieldName);
+
+            if (!processed)
+            {
+                WarningIn("void Foam::fieldValues::faceSource::write()")
+                    << "Requested field " << fieldName
+                    << " not found in database and not processed"
+                    << endl;
+            }
         }
 
         if (Pstream::master())

From 9ffb6c3f99e6a0621cecddaf5e5c157993b989f7 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Wed, 22 May 2013 17:26:05 +0100
Subject: [PATCH 33/33] STYLE: processor*patch: unused argument

---
 .../parallelProcessing/decomposePar/domainDecomposition.C       | 2 --
 1 file changed, 2 deletions(-)

diff --git a/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C b/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C
index f563978b60..15cbd85d11 100644
--- a/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C
+++ b/applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C
@@ -448,7 +448,6 @@ bool Foam::domainDecomposition::writeDecomposition()
                             curStart,
                             nPatches,
                             procMesh.boundaryMesh(),
-                            Pstream::worldComm,
                             procI,
                             curNeighbourProcessors[procPatchI]
                         );
@@ -476,7 +475,6 @@ bool Foam::domainDecomposition::writeDecomposition()
                             curStart,
                             nPatches,
                             procMesh.boundaryMesh(),
-                            Pstream::worldComm,
                             procI,
                             curNeighbourProcessors[procPatchI],
                             referPatch,