Merge branch 'master' into cvm

Conflicts:
	tutorials/lagrangian/reactingParcelFilmFoam/evaporationTest/system/extrudeToRegionMeshDict
	tutorials/lagrangian/reactingParcelFilmFoam/panel/system/extrudeToRegionMeshDict

src/OSspecific/POSIX/regExp.H

@@ -62,6 +62,7 @@ class regExp
         //- Precompiled regular expression
         mutable regex_t* preg_;
 
+
     // Private Member Functions
 
         //- Disallow default bitwise copy construct
@@ -72,6 +73,8 @@ class regExp
 
 public:
 
+    // Static Member Functions
+
         //- Is character a regular expression meta-character?
         //  any character: '.' \n
         //  quantifiers: '*', '+', '?' \n
@@ -102,66 +105,69 @@ public:
         //- Construct from std::string (or string), optionally ignoring case
         regExp(const std::string&, const bool ignoreCase=false);
 
+
     //- Destructor
     ~regExp();
 
 
     // Member functions
 
-    //- Access
+        // Access
 
-        //- Return true if a precompiled expression does not exist
-        inline bool empty() const
-        {
-            return !preg_;
-        }
+            //- Return true if a precompiled expression does not exist
+            inline bool empty() const
+            {
+                return !preg_;
+            }
 
-        //- Does a precompiled expression exist?
-        inline bool exists() const
-        {
-            return preg_ ? true : false;
-        }
+            //- Does a precompiled expression exist?
+            inline bool exists() const
+            {
+                return preg_ ? true : false;
+            }
 
-        //- Return the number of (groups)
-        inline int ngroups() const
-        {
-            return preg_ ? preg_->re_nsub : 0;
-        }
+            //- Return the number of (groups)
+            inline int ngroups() const
+            {
+                return preg_ ? preg_->re_nsub : 0;
+            }
 
 
-    //- Editing
+        // Editing
 
-        //- Compile pattern into a regular expression, optionally ignoring case
-        void set(const char*, const bool ignoreCase=false) const;
+            //- Compile pattern into a regular expression, optionally ignoring
+            //  case
+            void set(const char*, const bool ignoreCase=false) const;
 
-        //- Compile pattern into a regular expression, optionally ignoring case
-        void set(const std::string&, const bool ignoreCase=false) const;
+            //- Compile pattern into a regular expression, optionally ignoring
+            //  case
+            void set(const std::string&, const bool ignoreCase=false) const;
 
 
-        //- Release precompiled expression.
-        //  Returns true if precompiled expression existed before clear
-        bool clear() const;
+            //- Release precompiled expression.
+            //  Returns true if precompiled expression existed before clear
+            bool clear() const;
 
 
-    //- Searching
+        // Searching
 
-        //- Find position within string.
-        //  Returns the index where it begins or string::npos if not found
-        std::string::size_type find(const std::string& str) const;
+            //- Find position within string.
+            //  Returns the index where it begins or string::npos if not found
+            std::string::size_type find(const std::string& str) const;
 
-        //- Return true if it matches the entire string
-        //  The begin-of-line (^) and end-of-line ($) anchors are implicit
-        bool match(const std::string&) const;
+            //- Return true if it matches the entire string
+            //  The begin-of-line (^) and end-of-line ($) anchors are implicit
+            bool match(const std::string&) const;
 
-        //- Return true if it matches and sets the sub-groups matched
-        //  The begin-of-line (^) and end-of-line ($) anchors are implicit
-        bool match(const string&, List<string>& groups) const;
+            //- Return true if it matches and sets the sub-groups matched
+            //  The begin-of-line (^) and end-of-line ($) anchors are implicit
+            bool match(const string&, List<string>& groups) const;
 
-        //- Return true if the regex was found within string
-        bool search(const std::string& str) const
-        {
-            return std::string::npos != find(str);
-        }
+            //- Return true if the regex was found within string
+            bool search(const std::string& str) const
+            {
+                return std::string::npos != find(str);
+            }
 
 
     // Member Operators
@@ -173,7 +179,6 @@ public:
         //- Assign and compile pattern from string
         //  Always case sensitive
         void operator=(const std::string&);
-
 };
 
 

src/OpenFOAM/db/IOstreams/token/token.H

@@ -435,11 +435,11 @@ Ostream& operator<<(Ostream&, const token::compound&);
 
 
 #define defineCompoundTypeName(Type, Name)                                    \
-    typedef token::Compound<Type > tokenCompound##Name##_;                    \
+    typedef token::Compound<Type> tokenCompound##Name##_;                     \
     defineTemplateTypeNameAndDebugWithName(tokenCompound##Name##_, #Type, 0);
 
 #define addCompoundToRunTimeSelectionTable(Type, Name)                        \
-    token::compound::addIstreamConstructorToTable<token::Compound<Type > >    \
+    token::compound::addIstreamConstructorToTable<token::Compound<Type> >     \
         add##Name##IstreamConstructorToTable_;
 
 

src/OpenFOAM/db/Time/Time.C

@@ -228,6 +228,8 @@ Foam::Time::Time
 
     objectRegistry(*this),
 
+    libs_(),
+
     controlDict_
     (
         IOobject
@@ -257,9 +259,10 @@ Foam::Time::Time
     graphFormat_("raw"),
     runTimeModifiable_(true),
 
-    readLibs_(controlDict_, "libs"),
     functionObjects_(*this)
 {
+    libs_.open(controlDict_, "libs");
+
     // Explicitly set read flags on objectRegistry so anything constructed
     // from it reads as well (e.g. fvSolution).
     readOpt() = IOobject::MUST_READ_IF_MODIFIED;
@@ -313,6 +316,8 @@ Foam::Time::Time
 
     objectRegistry(*this),
 
+    libs_(),
+
     controlDict_
     (
         IOobject
@@ -343,9 +348,11 @@ Foam::Time::Time
     graphFormat_("raw"),
     runTimeModifiable_(true),
 
-    readLibs_(controlDict_, "libs"),
     functionObjects_(*this)
 {
+    libs_.open(controlDict_, "libs");
+
+
     // Explicitly set read flags on objectRegistry so anything constructed
     // from it reads as well (e.g. fvSolution).
     readOpt() = IOobject::MUST_READ_IF_MODIFIED;
@@ -401,6 +408,8 @@ Foam::Time::Time
 
     objectRegistry(*this),
 
+    libs_(),
+
     controlDict_
     (
         IOobject
@@ -430,9 +439,10 @@ Foam::Time::Time
     graphFormat_("raw"),
     runTimeModifiable_(true),
 
-    readLibs_(controlDict_, "libs"),
     functionObjects_(*this)
-{}
+{
+    libs_.open(controlDict_, "libs");
+}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

src/OpenFOAM/db/Time/Time.H

@@ -75,6 +75,10 @@ class Time
         //- file-change monitor for all registered files
         mutable autoPtr<fileMonitor> monitorPtr_;
 
+        //- Any loaded dynamic libraries. Make sure to construct before
+        //  reading controlDict.
+        dlLibraryTable libs_;
+
         //- The controlDict
         IOdictionary controlDict_;
 
@@ -166,9 +170,6 @@ private:
         //- Is runtime modification of dictionaries allowed?
         Switch runTimeModifiable_;
 
-        //- Instantiate a dummy class to cause the reading of dynamic libraries
-        dlLibraryTable::readDlLibrary readLibs_;
-
         //- Function objects executed at start and on ++, +=
         mutable functionObjectList functionObjects_;
 
@@ -375,6 +376,12 @@ public:
                 return functionObjects_;
             }
 
+            //- External access to the loaded libraries
+            dlLibraryTable& libs()
+            {
+                return libs_;
+            }
+
             //- Return true if time currently being sub-cycled, otherwise false
             bool subCycling() const
             {

src/OpenFOAM/db/dictionary/functionEntries/codeStream/codeStream.C

@@ -36,6 +36,7 @@ License
 #include "Time.H"
 #include "PstreamReduceOps.H"
 
+#include "long.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -61,6 +62,36 @@ const Foam::word Foam::functionEntries::codeStream::codeTemplateC
     = "codeStreamTemplate.C";
 
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+const Foam::dictionary& Foam::functionEntries::codeStream::topDict
+(
+    const dictionary& dict
+)
+{
+    const dictionary& p = dict.parent();
+
+    if (&p != &dict && !p.name().empty())
+    {
+        return topDict(p);
+    }
+    else
+    {
+        return dict;
+    }
+}
+
+
+Foam::dlLibraryTable& Foam::functionEntries::codeStream::libs
+(
+    const dictionary& dict
+)
+{
+    const IOdictionary& d = static_cast<const IOdictionary&>(topDict(dict));
+    return const_cast<Time&>(d.time()).libs();
+}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 bool Foam::functionEntries::codeStream::execute
@@ -99,8 +130,11 @@ bool Foam::functionEntries::codeStream::execute
     const fileName libPath = dynCode.libPath();
 
     // see if library is loaded
-    void* lib = dlLibraryTable::findLibrary(libPath);
-
+    void* lib = NULL;
+    if (isA<IOdictionary>(topDict(parentDict)))
+    {
+        lib = libs(parentDict).findLibrary(libPath);
+    }
 
     if (!lib)
     {
@@ -110,9 +144,23 @@ bool Foam::functionEntries::codeStream::execute
 
     // nothing loaded
     // avoid compilation if possible by loading an existing library
-    if (!lib && dlLibraryTable::open(libPath, false))
+    if (!lib)
     {
-        lib = dlLibraryTable::findLibrary(libPath);
+        if (isA<IOdictionary>(topDict(parentDict)))
+        {
+            // Cached access to dl libs. Guarantees clean up upon destruction
+            // of Time.
+            dlLibraryTable& dlLibs = libs(parentDict);
+            if (dlLibs.open(libPath, false))
+            {
+                lib = dlLibs.findLibrary(libPath);
+            }
+        }
+        else
+        {
+            // Uncached opening of libPath
+            lib = dlOpen(libPath);
+        }
     }
 
 
@@ -167,19 +215,30 @@ bool Foam::functionEntries::codeStream::execute
         // all processes must wait for compile to finish
         reduce(create, orOp<bool>());
 
-        if (!dlLibraryTable::open(libPath, false))
+        if (isA<IOdictionary>(topDict(parentDict)))
         {
-            FatalIOErrorIn
-            (
-                "functionEntries::codeStream::execute(..)",
-                parentDict
-            )   << "Failed loading library " << libPath << nl
-                << "Did you add all libraries to the 'libs' entry"
-                << " in system/controlDict?"
-                << exit(FatalIOError);
-        }
+            // Cached access to dl libs. Guarantees clean up upon destruction
+            // of Time.
+            dlLibraryTable& dlLibs = libs(parentDict);
+            if (!dlLibs.open(libPath, false))
+            {
+                FatalIOErrorIn
+                (
+                    "functionEntries::codeStream::execute(..)",
+                    parentDict
+                )   << "Failed loading library " << libPath << nl
+                    << "Did you add all libraries to the 'libs' entry"
+                    << " in system/controlDict?"
+                    << exit(FatalIOError);
+            }
 
-        lib = dlLibraryTable::findLibrary(libPath);
+            lib = dlLibs.findLibrary(libPath);
+        }
+        else
+        {
+            // Uncached opening of libPath
+            lib = dlOpen(libPath);
+        }
     }
 
 

src/OpenFOAM/db/dictionary/functionEntries/codeStream/codeStream.H

@@ -98,6 +98,8 @@ SourceFiles
 
 namespace Foam
 {
+class dlLibraryTable;
+
 namespace functionEntries
 {
 
@@ -113,9 +115,14 @@ class codeStream
         //- Interpreter function type
         typedef void (*streamingFunctionType)(Ostream&, const dictionary&);
 
-
     // Private Member Functions
 
+        //- Helper function: parent (of parent etc.) of dictionary up to the top
+        static const dictionary& topDict(const dictionary&);
+
+        //- Helper function: access to dlLibraryTable of Time
+        static dlLibraryTable& libs(const dictionary& dict);
+
         //- Disallow default bitwise copy construct
         codeStream(const codeStream&);
 

src/OpenFOAM/db/dynamicLibrary/dlLibraryTable/dlLibraryTable.C

@@ -25,10 +25,11 @@ License
 
 #include "dlLibraryTable.H"
 #include "OSspecific.H"
+#include "long.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-Foam::dlLibraryTable Foam::dlLibraryTable::loadedLibraries;
+defineTypeNameAndDebug(Foam::dlLibraryTable, 0);
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -39,11 +40,13 @@ Foam::dlLibraryTable::dlLibraryTable()
 {}
 
 
-Foam::dlLibraryTable::readDlLibrary::readDlLibrary
+Foam::dlLibraryTable::dlLibraryTable
 (
     const dictionary& dict,
     const word& libsEntry
 )
+:
+    HashTable<fileName, void*, Hash<void*> >()
 {
     open(dict, libsEntry);
 }
@@ -58,7 +61,11 @@ Foam::dlLibraryTable::~dlLibraryTable()
         // bug in dlclose - does not call static destructors of
         // loaded library when actually unloading the library.
         // See https://bugzilla.novell.com/show_bug.cgi?id=680125 and 657627.
-        // Seems related to using a non-system compiler!
+        if (debug)
+        {
+            Info<< "dlLibraryTable::~dlLibraryTable() : closing " << iter()
+                << " with handle " << long(iter.key()) << endl;
+        }
         dlClose(iter.key());
     }
 }
@@ -76,6 +83,12 @@ bool Foam::dlLibraryTable::open
     {
         void* functionLibPtr = dlOpen(functionLibName);
 
+        if (debug)
+        {
+            Info<< "dlLibraryTable::open : opened " << functionLibName
+                << " resulting in handle " << long(functionLibPtr) << endl;
+        }
+
         if (!functionLibPtr)
         {
             if (verbose)
@@ -91,14 +104,7 @@ bool Foam::dlLibraryTable::open
         }
         else
         {
-            if (loadedLibraries.insert(functionLibPtr, functionLibName))
-            {
-                return true;
-            }
-            else
-            {
-                return false;
-            }
+            return insert(functionLibPtr, functionLibName);
         }
     }
     else
@@ -117,7 +123,13 @@ bool Foam::dlLibraryTable::close
     void* libPtr = findLibrary(functionLibName);
     if (libPtr)
     {
-        loadedLibraries.erase(libPtr);
+        if (debug)
+        {
+            Info<< "dlLibraryTable::close : closing " << functionLibName
+                << " with handle " << long(libPtr) << endl;
+        }
+
+        erase(libPtr);
 
         if (!dlClose(libPtr))
         {
@@ -141,7 +153,7 @@ bool Foam::dlLibraryTable::close
 
 void* Foam::dlLibraryTable::findLibrary(const fileName& functionLibName)
 {
-    forAllConstIter(dlLibraryTable, loadedLibraries, iter)
+    forAllConstIter(dlLibraryTable, *this, iter)
     {
         if (iter() == functionLibName)
         {

src/OpenFOAM/db/dynamicLibrary/dlLibraryTable/dlLibraryTable.H

@@ -63,30 +63,18 @@ class dlLibraryTable
 
 public:
 
-    // Static data members
-
-        //- Static data someStaticData
-        static dlLibraryTable loadedLibraries;
-
-
-    // Public classes
-
-        //- Class whose construction causes the reading of dynamic libraries
-        class readDlLibrary
-        {
-        public:
-
-            //- Read all the libraries listed in the 'libsEntry' entry in the
-            //  given dictionary if present
-            readDlLibrary(const dictionary&, const word& libsEntry);
-        };
-
+    // Declare name of the class and its debug switch
+    ClassName("dlLibraryTable");
 
     // Constructors
 
         //- Construct null
         dlLibraryTable();
 
+        //- Construct from dictionary and name of 'libs' entry giving
+        //  the libraries to load
+        dlLibraryTable(const dictionary&, const word&);
+
 
     //- Destructor
     ~dlLibraryTable();
@@ -95,23 +83,23 @@ public:
     // Member Functions
 
         //- Open the named library, optionally with warnings if problems occur
-        static bool open(const fileName& name, const bool verbose = true);
+        bool open(const fileName& name, const bool verbose = true);
 
         //- Close the named library, optionally with warnings if problems occur
-        static bool close(const fileName& name, const bool verbose = true);
+        bool close(const fileName& name, const bool verbose = true);
 
         //- Find the handle of the named library
-        static void* findLibrary(const fileName& name);
+        void* findLibrary(const fileName& name);
 
         //- Open all the libraries listed in the 'libsEntry' entry in the
         //  given dictionary if present
-        static bool open(const dictionary&, const word& libsEntry);
+        bool open(const dictionary&, const word& libsEntry);
 
         //- Open all the libraries listed in the 'libsEntry' entry in the
         //  given dictionary if present and check the additions
         //  to the given constructor table
         template<class TablePtr>
-        static bool open
+        bool open
         (
             const dictionary&,
             const word& libsEntry,

src/OpenFOAM/db/dynamicLibrary/dynamicCode/dynamicCodeContext.C

@@ -46,7 +46,6 @@ Foam::dynamicCodeContext::dynamicCodeContext(const dictionary& dict)
         const entry& codeEntry = dict.lookupEntry("code", false, false);
         code_ = stringOps::trim(codeEntry.stream());
         stringOps::inplaceExpand(code_, dict);
-        addLineDirective(code_, codeEntry.startLineNumber(), dict.name());
     }
 
     // note: removes any leading/trailing whitespace
@@ -64,7 +63,6 @@ Foam::dynamicCodeContext::dynamicCodeContext(const dictionary& dict)
     {
         include_ = stringOps::trim(includePtr->stream());
         stringOps::inplaceExpand(include_, dict);
-        addLineDirective(include_, includePtr->startLineNumber(), dict.name());
     }
 
     // optional
@@ -92,6 +90,28 @@ Foam::dynamicCodeContext::dynamicCodeContext(const dictionary& dict)
     OSHA1stream os;
     os  << include_ << options_ << libs_ << localCode_ << code_;
     sha1_ = os.digest();
+
+
+
+    // Add line number after calculating sha1 since includes processorDDD
+    // in path which differs between processors.
+
+    {
+        const entry& codeEntry = dict.lookupEntry("code", false, false);
+        addLineDirective(code_, codeEntry.startLineNumber(), dict.name());
+    }
+    if (includePtr)
+    {
+        addLineDirective(include_, includePtr->startLineNumber(), dict.name());
+    }
+    if (optionsPtr)
+    {
+        addLineDirective(options_, optionsPtr->startLineNumber(), dict.name());
+    }
+    if (libsPtr)
+    {
+        addLineDirective(libs_, libsPtr->startLineNumber(), dict.name());
+    }
 }
 
 

src/OpenFOAM/db/functionObjects/functionObject/functionObject.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,6 +26,7 @@ License
 #include "functionObject.H"
 #include "dictionary.H"
 #include "dlLibraryTable.H"
+#include "Time.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -57,7 +58,7 @@ Foam::autoPtr<Foam::functionObject> Foam::functionObject::New
         Info<< "Selecting function " << functionType << endl;
     }
 
-    dlLibraryTable::open
+    const_cast<Time&>(t).libs().open
     (
         functionDict,
         "functionObjectLibs",

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -104,7 +104,7 @@ const Foam::GAMGAgglomeration& Foam::GAMGAgglomeration::New
     {
         const word agglomeratorType(controlDict.lookup("agglomerator"));
 
-        dlLibraryTable::open
+        const_cast<Time&>(mesh.thisDb().time()).libs().open
         (
             controlDict,
             "geometricGAMGAgglomerationLibs",
@@ -159,7 +159,7 @@ const Foam::GAMGAgglomeration& Foam::GAMGAgglomeration::New
     {
         const word agglomeratorType(controlDict.lookup("agglomerator"));
 
-        dlLibraryTable::open
+        const_cast<Time&>(mesh.thisDb().time()).libs().open
         (
             controlDict,
             "algebraicGAMGAgglomerationLibs",

src/OpenFOAM/meshes/polyMesh/polyMesh.C

@@ -36,6 +36,7 @@ License
 #include "polyMeshTetDecomposition.H"
 #include "indexedOctree.H"
 #include "treeDataCell.H"
+#include "SubField.H"
 
 #include "pointMesh.H"
 

src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C

@@ -55,6 +55,8 @@ Foam::label Foam::polyMeshTetDecomposition::findSharedBasePoint
 
     const point& oCc = pC[oCI];
 
+    List<scalar> tetQualities(2, 0.0);
+
     forAll(f, faceBasePtI)
     {
         scalar thisBaseMinTetQuality = VGREAT;
@@ -66,8 +68,6 @@ Foam::label Foam::polyMeshTetDecomposition::findSharedBasePoint
             label facePtI = (tetPtI + faceBasePtI) % f.size();
             label otherFacePtI = f.fcIndex(facePtI);
 
-            List<scalar> tetQualities(2, 0.0);
-
             {
                 // owner cell tet
                 label ptAI = f[facePtI];

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C

@@ -35,6 +35,7 @@ License
 #include "Time.H"
 #include "diagTensor.H"
 #include "transformField.H"
+#include "SubField.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 

src/OpenFOAM/meshes/polyMesh/syncTools/syncToolsTemplates.C

@@ -31,6 +31,7 @@ License
 #include "contiguous.H"
 #include "transform.H"
 #include "transformList.H"
+#include "SubField.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 

src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H

@@ -56,8 +56,6 @@ SourceFiles
 #include "DynamicList.H"
 #include "edgeList.H"
 #include "pointField.H"
-#include "SubField.H"
-#include "SubList.H"
 #include "faceList.H"
 #include "cellList.H"
 #include "cellShapeList.H"
@@ -65,8 +63,6 @@ SourceFiles
 #include "boolList.H"
 #include "HashSet.H"
 #include "Map.H"
-#include "EdgeMap.H"
-#include "boundBox.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -29,7 +29,7 @@ License
 #include "ListOps.H"
 #include "unitConversion.H"
 #include "SortableList.H"
-
+#include "EdgeMap.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshClear.C

@@ -27,6 +27,8 @@ License
 #include "indexedOctree.H"
 #include "treeDataCell.H"
 #include "demandDrivenData.H"
+#include "indexedOctree.H"
+#include "treeDataCell.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFindCell.C

@@ -25,6 +25,7 @@ License
 
 #include "primitiveMesh.H"
 #include "cell.H"
+#include "boundBox.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 

src/OpenFOAM/meshes/primitiveShapes/plane/plane.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -96,7 +96,8 @@ void Foam::plane::calcPntAndVec
             "    const point&,\n"
             "    const point&\n"
             ")\n"
-        ) << "Bad points." << abort(FatalError);
+        )   << "Bad points:" << point1 << ' ' << point2 << ' ' << point3
+            << abort(FatalError);
     }
 
     unitVector_ = line12 ^ line23;
@@ -112,7 +113,8 @@ void Foam::plane::calcPntAndVec
             "    const point&,\n"
             "    const point&\n"
             ")\n"
-        )   << "Plane normal defined with zero length"
+        )   << "Plane normal defined with zero length" << nl
+            << "Bad points:" << point1 << ' ' << point2 << ' ' << point3
             << abort(FatalError);
     }
 
@@ -137,7 +139,7 @@ Foam::plane::plane(const vector& normalVector)
     else
     {
         FatalErrorIn("plane::plane(const vector&)")
-            << "plane normal has zero length"
+            << "plane normal has zero length. basePoint:" << basePoint_
             << abort(FatalError);
     }
 }
@@ -158,7 +160,7 @@ Foam::plane::plane(const point& basePoint, const vector& normalVector)
     else
     {
         FatalErrorIn("plane::plane(const point&, const vector&)")
-            << "plane normal has zero length"
+            << "plane normal has zero length. basePoint:" << basePoint_
             << abort(FatalError);
     }
 }
@@ -228,8 +230,7 @@ Foam::plane::plane(const dictionary& dict)
         (
             "plane::plane(const dictionary&)",
             dict
-        )
-            << "Invalid plane type: " << planeType
+        )   << "Invalid plane type: " << planeType
             << abort(FatalIOError);
     }
 }
@@ -250,7 +251,7 @@ Foam::plane::plane(Istream& is)
     else
     {
         FatalErrorIn("plane::plane(Istream& is)")
-            << "plane normal has zero length"
+            << "plane normal has zero length. basePoint:" << basePoint_
             << abort(FatalError);
     }
 }

src/OpenFOAM/primitives/strings/wordRe/wordRe.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2009-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2009-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -52,6 +52,7 @@ SourceFiles
 
 #include "word.H"
 #include "regExp.H"
+#include "keyType.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -60,7 +61,6 @@ namespace Foam
 
 // Forward declaration of friend functions and operators
 class wordRe;
-
 class Istream;
 class Ostream;
 
@@ -110,6 +110,7 @@ public:
         //- Test string for regular expression meta characters
         static inline bool isPattern(const string&);
 
+
     // Constructors
 
         //- Construct null
@@ -118,71 +119,83 @@ public:
         //- Construct as copy
         inline wordRe(const wordRe&);
 
+        //- Construct from keyType
+        inline wordRe(const keyType&, const compOption=LITERAL);
+
         //- Construct as copy of word
         inline wordRe(const word&);
 
         //- Construct as copy of character array
         //  Optionally specify how it should be treated.
-        inline wordRe(const char*, const compOption=LITERAL);
+        inline wordRe(const char*, const compOption = LITERAL);
 
         //- Construct as copy of string.
         //  Optionally specify how it should be treated.
-        inline wordRe(const string&, const compOption=LITERAL);
+        inline wordRe(const string&, const compOption = LITERAL);
 
         //- Construct as copy of std::string
         //  Optionally specify how it should be treated.
-        inline wordRe(const std::string&, const compOption=LITERAL);
+        inline wordRe(const std::string&, const compOption = LITERAL);
 
         //- Construct from Istream
         //  Words are treated as literals, strings with an auto-test
         wordRe(Istream&);
 
+
     // Member functions
 
-    //- Access
+        // Access
 
-        //- Should be treated as a match rather than a literal string?
-        inline bool isPattern() const;
+            //- Should be treated as a match rather than a literal string?
+            inline bool isPattern() const;
 
-    //- Infrastructure
 
-        //- Compile the regular expression
-        inline bool compile() const;
+        // Infrastructure
 
-        //- Possibly compile the regular expression, with greater control
-        inline bool compile(const compOption) const;
+            //- Compile the regular expression
+            inline bool compile() const;
 
-        //- Recompile an existing regular expression
-        inline bool recompile() const;
+            //- Possibly compile the regular expression, with greater control
+            inline bool compile(const compOption) const;
 
-        //- Frees precompiled regular expression, making wordRe a literal.
-        //  Optionally strips invalid word characters
-        inline void uncompile(const bool doStripInvalid=false) const;
+            //- Recompile an existing regular expression
+            inline bool recompile() const;
 
-    //- Editing
+            //- Frees precompiled regular expression, making wordRe a literal.
+            //  Optionally strips invalid word characters
+            inline void uncompile(const bool doStripInvalid = false) const;
 
-        //- Copy string, auto-test for regular expression or other options
-        inline void set(const std::string&, const compOption=DETECT);
 
-        //- Copy string, auto-test for regular expression or other options
-        inline void set(const char*, const compOption=DETECT);
+        // Editing
 
-        //- Clear string and precompiled regular expression
-        inline void clear();
+            //- Copy string, auto-test for regular expression or other options
+            inline void set(const std::string&, const compOption = DETECT);
 
-    //- Searching
+            //- Copy string, auto-test for regular expression or other options
+            inline void set(const char*, const compOption = DETECT);
 
-        //- Smart match as regular expression or as a string
-        //  Optionally force a literal match only
-        inline bool match(const std::string&, bool literalMatch=false) const;
+            //- Clear string and precompiled regular expression
+            inline void clear();
 
-    //- Miscellaneous
 
-        //- Return a string with quoted meta-characters
-        inline string quotemeta() const;
+        // Searching
 
-        //- Output some basic info
-        Ostream& info(Ostream&) const;
+            //- Smart match as regular expression or as a string
+            //  Optionally force a literal match only
+            inline bool match
+            (
+                const std::string&,
+                bool literalMatch = false
+            ) const;
+
+
+        // Miscellaneous
+
+            //- Return a string with quoted meta-characters
+            inline string quotemeta() const;
+
+            //- Output some basic info
+            Ostream& info(Ostream&) const;
 
 
     // Member operators
@@ -196,6 +209,10 @@ public:
             //- Copy word, never a regular expression
             inline const wordRe& operator=(const word&);
 
+            //- Copy keyType, auto-test for regular expression
+            //  Always case sensitive
+            inline const wordRe& operator=(const keyType&);
+
             //- Copy string, auto-test for regular expression
             //  Always case sensitive
             inline const wordRe& operator=(const string&);

src/OpenFOAM/primitives/strings/wordRe/wordReI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2009-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2009-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -65,6 +65,18 @@ inline Foam::wordRe::wordRe(const word& str)
 {}
 
 
+inline Foam::wordRe::wordRe(const keyType& str, const compOption opt)
+:
+    word(str, false),
+    re_()
+{
+    if (str.isPattern())
+    {
+        compile(opt);
+    }
+}
+
+
 inline Foam::wordRe::wordRe(const char* str, const compOption opt)
 :
     word(str, false),
@@ -236,6 +248,17 @@ inline const Foam::wordRe& Foam::wordRe::operator=(const word& str)
 }
 
 
+inline const Foam::wordRe& Foam::wordRe::operator=(const keyType& str)
+{
+    string::operator=(str);
+    if (str.isPattern())
+    {
+        compile();
+    }
+    return *this;
+}
+
+
 inline const Foam::wordRe& Foam::wordRe::operator=(const string& str)
 {
     string::operator=(str);

src/dynamicFvMesh/dynamicFvMesh/dynamicFvMeshNew.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -49,7 +49,7 @@ Foam::autoPtr<Foam::dynamicFvMesh> Foam::dynamicFvMesh::New(const IOobject& io)
 
     Info<< "Selecting dynamicFvMesh " << dynamicFvMeshTypeName << endl;
 
-    dlLibraryTable::open
+    const_cast<Time&>(io.time()).libs().open
     (
         dict,
         "dynamicFvMeshLibs",

src/dynamicMesh/motionSolver/motionSolver.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -79,7 +79,7 @@ Foam::autoPtr<Foam::motionSolver> Foam::motionSolver::New(const polyMesh& mesh)
 
     Info<< "Selecting motion solver: " << solverTypeName << endl;
 
-    dlLibraryTable::open
+    const_cast<Time&>(mesh.time()).libs().open
     (
         solverDict,
         "motionSolverLibs",

src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C

@@ -81,30 +81,6 @@ void Foam::addPatchCellLayer::addVertex
             f[fp++] = pointI;
         }
     }
-
-    // Check for duplicates.
-    if (debug)
-    {
-        label n = 0;
-        for (label i = 0; i < fp; i++)
-        {
-            if (f[i] == pointI)
-            {
-                n++;
-
-                if (n == 2)
-                {
-                    f.setSize(fp);
-                    FatalErrorIn
-                    (
-                        "addPatchCellLayer::addVertex(const label, face&"
-                        ", label&)"
-                    )   << "Point " << pointI << " already present in face "
-                        << f << abort(FatalError);
-                }
-            }
-        }
-    }
 }
 
 
@@ -1641,11 +1617,11 @@ void Foam::addPatchCellLayer::setRefinement
                         {
                             label offset =
                                 addedPoints_[vStart].size() - numEdgeSideFaces;
-                            for (label ioff = offset; ioff > 0; ioff--)
+                            for (label ioff = offset-1; ioff >= 0; ioff--)
                             {
                                 addVertex
                                 (
-                                    addedPoints_[vStart][ioff-1],
+                                    addedPoints_[vStart][ioff],
                                     newFace,
                                     newFp
                                 );
@@ -1660,6 +1636,51 @@ void Foam::addPatchCellLayer::setRefinement
 
                         newFace.setSize(newFp);
 
+                        if (debug)
+                        {
+                            labelHashSet verts(2*newFace.size());
+                            forAll(newFace, fp)
+                            {
+                                if (!verts.insert(newFace[fp]))
+                                {
+                                    FatalErrorIn
+                                    (
+                                        "addPatchCellLayer::setRefinement(..)"
+                                    )   << "Duplicate vertex in face"
+                                        << " to be added." << nl
+                                        << "newFace:" << newFace << nl
+                                        << "points:"
+                                        <<  UIndirectList<point>
+                                            (
+                                                meshMod.points(),
+                                                newFace
+                                            ) << nl
+                                        << "Layer:" << i
+                                        << " out of:" << numEdgeSideFaces << nl
+                                        << "ExtrudeEdge:" << meshEdgeI
+                                        << " at:"
+                                        <<  mesh_.edges()[meshEdgeI].line
+                                            (
+                                                mesh_.points()
+                                            ) << nl
+                                        << "string:" << stringedVerts
+                                        << "stringpoints:"
+                                        << UIndirectList<point>
+                                            (
+                                                pp.localPoints(),
+                                                stringedVerts
+                                            ) << nl
+                                        << "stringNLayers:"
+                                        <<  UIndirectList<label>
+                                            (
+                                                nPointLayers,
+                                                stringedVerts
+                                            ) << nl
+                                        << abort(FatalError);
+                                }
+                            }
+                        }
+
                         label nbrFaceI = nbrFace
                         (
                             pp.edgeFaces(),

src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -109,6 +109,12 @@ bool Foam::combineFaces::validFace
         return false;
     }
 
+    bool isNonManifold = bigFace.checkPointManifold(false, NULL);
+    if (isNonManifold)
+    {
+        return false;
+    }
+
     // Check for convexness
     face f(getOutsideFace(bigFace));
 
@@ -984,6 +990,7 @@ void Foam::combineFaces::setUnrefinement
                 zoneFlip                        // face flip in zone
             )
         );
+        restoredFaces.insert(masterFaceI, masterFaceI);
 
         // Add the previously removed faces
         for (label i = 1; i < faces.size(); i++)
@@ -991,7 +998,7 @@ void Foam::combineFaces::setUnrefinement
             //Pout<< "Restoring removed face with vertices " << faces[i]
             //    << endl;
 
-            meshMod.setAction
+            label faceI = meshMod.setAction
             (
                 polyAddFace
                 (
@@ -1007,6 +1014,7 @@ void Foam::combineFaces::setUnrefinement
                     zoneFlip                // zoneFlip
                 )
             );
+            restoredFaces.insert(faceI, masterFaceI);
         }
 
         // Clear out restored set

src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -302,6 +302,35 @@ void Foam::polyTopoChange::getMergeSets
 }
 
 
+bool Foam::polyTopoChange::hasValidPoints(const face& f) const
+{
+    forAll(f, fp)
+    {
+        if (f[fp] < 0 || f[fp] >= points_.size())
+        {
+            return false;
+        }
+    }
+    return true;
+}
+
+
+Foam::pointField Foam::polyTopoChange::facePoints(const face& f) const
+{
+    pointField points(f.size());
+    forAll(f, fp)
+    {
+        if (f[fp] < 0 && f[fp] >= points_.size())
+        {
+            FatalErrorIn("polyTopoChange::facePoints(const face&) const")
+                << "Problem." << abort(FatalError);
+        }
+        points[fp] = points_[f[fp]];
+    }
+    return points;
+}
+
+
 void Foam::polyTopoChange::checkFace
 (
     const face& f,
@@ -329,7 +358,14 @@ void Foam::polyTopoChange::checkFace
                 << "f:" << f
                 << " faceI(-1 if added face):" << faceI
                 << " own:" << own << " nei:" << nei
-                << " patchI:" << patchI << abort(FatalError);
+                << " patchI:" << patchI << nl;
+            if (hasValidPoints(f))
+            {
+                FatalError
+                        << "points (removed points marked with "
+                        << vector::max << ") " << facePoints(f);
+            }
+            FatalError << abort(FatalError);
         }
     }
     else
@@ -344,7 +380,14 @@ void Foam::polyTopoChange::checkFace
                 << "f:" << f
                 << " faceI(-1 if added face):" << faceI
                 << " own:" << own << " nei:" << nei
-                << " patchI:" << patchI << abort(FatalError);
+                << " patchI:" << patchI << nl;
+            if (hasValidPoints(f))
+            {
+                FatalError
+                        << "points (removed points marked with "
+                        << vector::max << ") : " << facePoints(f);
+            }
+            FatalError << abort(FatalError);
         }
 
         if (nei <= own)
@@ -358,7 +401,14 @@ void Foam::polyTopoChange::checkFace
                 << "f:" << f
                 << " faceI(-1 if added face):" << faceI
                 << " own:" << own << " nei:" << nei
-                << " patchI:" << patchI << abort(FatalError);
+                << " patchI:" << patchI << nl;
+            if (hasValidPoints(f))
+            {
+                FatalError
+                        << "points (removed points marked with "
+                        << vector::max << ") : " << facePoints(f);
+            }
+            FatalError << abort(FatalError);
         }
     }
 
@@ -373,7 +423,14 @@ void Foam::polyTopoChange::checkFace
             << "f:" << f
             << " faceI(-1 if added face):" << faceI
             << " own:" << own << " nei:" << nei
-            << " patchI:" << patchI << abort(FatalError);
+            << " patchI:" << patchI << nl;
+            if (hasValidPoints(f))
+            {
+                FatalError
+                        << "points (removed points marked with "
+                        << vector::max << ") : " << facePoints(f);
+            }
+            FatalError << abort(FatalError);
     }
     if (faceI >= 0 && faceI < faces_.size() && faceRemoved(faceI))
     {
@@ -386,7 +443,14 @@ void Foam::polyTopoChange::checkFace
             << "f:" << f
             << " faceI(-1 if added face):" << faceI
             << " own:" << own << " nei:" << nei
-            << " patchI:" << patchI << abort(FatalError);
+            << " patchI:" << patchI << nl;
+            if (hasValidPoints(f))
+            {
+                FatalError
+                        << "points (removed points marked with "
+                        << vector::max << ") : " << facePoints(f);
+            }
+            FatalError << abort(FatalError);
     }
     forAll(f, fp)
     {
@@ -401,7 +465,14 @@ void Foam::polyTopoChange::checkFace
                 << "f:" << f
                 << " faceI(-1 if added face):" << faceI
                 << " own:" << own << " nei:" << nei
-                << " patchI:" << patchI << abort(FatalError);
+                << " patchI:" << patchI << nl;
+            if (hasValidPoints(f))
+            {
+                FatalError
+                        << "points (removed points marked with "
+                        << vector::max << ") : " << facePoints(f);
+            }
+            FatalError << abort(FatalError);
         }
     }
 }
@@ -729,8 +800,16 @@ void Foam::polyTopoChange::getFaceOrder
                 << " neighbour " << faceNeighbour_[faceI]
                 << " region " << region_[faceI] << endl
                 << "This is usually caused by not specifying a patch for"
-                << " a boundary face."
-                << abort(FatalError);
+                << " a boundary face." << nl
+                << "Switch on the polyTopoChange::debug flag to catch"
+                << " this error earlier." << nl;
+            if (hasValidPoints(faces_[faceI]))
+            {
+                FatalError
+                        << "points (removed points marked with "
+                        << vector::max << ") " << facePoints(faces_[faceI]);
+            }
+            FatalError << abort(FatalError);
         }
     }
 }

src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -239,6 +239,12 @@ class polyTopoChange
             List<objectMap>& cellsFromCells
         );
 
+        //- Are all face vertices valid
+        bool hasValidPoints(const face&) const;
+
+        //- Return face points
+        pointField facePoints(const face& f) const;
+
         //- Check inputs to modFace or addFace
         void checkFace
         (

src/finiteVolume/Make/files

@@ -208,9 +208,6 @@ $(interpolation)/interpolationPoint/pointMVCWeight.C
 $(interpolation)/interpolationPoint/makeInterpolationPoint.C
 
 volPointInterpolation = interpolation/volPointInterpolation
-/*
-$(volPointInterpolation)/pointPatchInterpolation/pointPatchInterpolation.C
-*/
 $(volPointInterpolation)/volPointInterpolation.C
 
 surfaceInterpolation = interpolation/surfaceInterpolation
@@ -353,6 +350,11 @@ $(general)/findRefCell/findRefCell.C
 $(general)/adjustPhi/adjustPhi.C
 $(general)/bound/bound.C
 
+solutionControl = $(general)/solutionControl
+$(solutionControl)/solutionControl/solutionControl.C
+$(solutionControl)/simpleControl/simpleControl.C
+$(solutionControl)/pimpleControl/pimpleControl.C
+
 porousMedia = $(general)/porousMedia
 $(porousMedia)/porousZone.C
 $(porousMedia)/porousZones.C

src/finiteVolume/cfdTools/general/include/readPIMPLEControls.H

@@ -1,17 +0,0 @@
-    const dictionary& pimple = mesh.solutionDict().subDict("PIMPLE");
-
-    const int nOuterCorr =
-        pimple.lookupOrDefault<int>("nOuterCorrectors", 1);
-
-    const int nCorr =
-        pimple.lookupOrDefault<int>("nCorrectors", 1);
-
-    const int nNonOrthCorr =
-        pimple.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);
-
-    const bool momentumPredictor =
-        pimple.lookupOrDefault("momentumPredictor", true);
-
-    const bool transonic =
-        pimple.lookupOrDefault("transonic", false);
-

src/finiteVolume/cfdTools/general/include/readPISOControls.H

@@ -1,17 +1,17 @@
-    const dictionary& piso = mesh.solutionDict().subDict("PISO");
+    const dictionary& pisoDict = mesh.solutionDict().subDict("PISO");
 
     const int nOuterCorr =
-        piso.lookupOrDefault<int>("nOuterCorrectors", 1);
+        pisoDict.lookupOrDefault<int>("nOuterCorrectors", 1);
 
     const int nCorr =
-        piso.lookupOrDefault<int>("nCorrectors", 1);
+        pisoDict.lookupOrDefault<int>("nCorrectors", 1);
 
     const int nNonOrthCorr =
-        piso.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);
+        pisoDict.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);
 
     const bool momentumPredictor =
-        piso.lookupOrDefault("momentumPredictor", true);
+        pisoDict.lookupOrDefault("momentumPredictor", true);
 
     const bool transonic =
-        piso.lookupOrDefault("transonic", false);
+        pisoDict.lookupOrDefault("transonic", false);
 

src/finiteVolume/cfdTools/general/include/readSIMPLEControls.H

@@ -1,11 +0,0 @@
-    const dictionary& simple = mesh.solutionDict().subDict("SIMPLE");
-
-    const int nNonOrthCorr =
-        simple.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);
-
-    const bool momentumPredictor =
-        simple.lookupOrDefault("momentumPredictor", true);
-
-    const bool transonic =
-        simple.lookupOrDefault("transonic", false);
-

src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C

@@ -0,0 +1,166 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "pimpleControl.H"
+#include "Switch.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(pimpleControl, 0);
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+void Foam::pimpleControl::read()
+{
+    solutionControl::read(false);
+
+    // Read solution controls
+    const dictionary& pimpleDict = dict();
+    nOuterCorr_ = pimpleDict.lookupOrDefault<label>("nOuterCorrectors", 1);
+    nCorr_ = pimpleDict.lookupOrDefault<label>("nCorrectors", 1);
+    turbOnFinalIterOnly_ =
+        pimpleDict.lookupOrDefault<Switch>("turbOnFinalIterOnly", true);
+}
+
+
+bool Foam::pimpleControl::criteriaSatisfied()
+{
+    if ((corr_ == 0) || residualControl_.empty() || finalIter())
+    {
+        return false;
+    }
+
+    bool firstIter = corr_ == 1;
+
+    bool achieved = true;
+    const dictionary& solverDict = mesh_.solverPerformanceDict();
+
+    forAllConstIter(dictionary, solverDict, iter)
+    {
+        const word& variableName = iter().keyword();
+        label fieldI = applyToField(variableName);
+        if (fieldI != -1)
+        {
+            const List<lduMatrix::solverPerformance> sp(iter().stream());
+            const scalar residual = sp.last().initialResidual();
+
+            if (firstIter)
+            {
+                residualControl_[fieldI].initialResidual =
+                    sp.first().initialResidual();
+            }
+
+            bool absCheck = residual < residualControl_[fieldI].absTol;
+
+            bool relCheck = false;
+
+            scalar relative = 0.0;
+            if (!firstIter)
+            {
+                scalar iniRes =
+                    residualControl_[fieldI].initialResidual
+                  + ROOTVSMALL;
+
+                relative = residual/iniRes;
+
+                relCheck = relative < residualControl_[fieldI].relTol;
+            }
+
+            achieved = achieved && (absCheck || relCheck);
+
+            if (debug)
+            {
+                Info<< algorithmName_ << "loop statistics:" << endl;
+
+                Info<< "    " << variableName << " iter " << corr_
+                    << ": ini res = "
+                    << residualControl_[fieldI].initialResidual
+                    << ", abs tol = " << residual
+                    << " (" << residualControl_[fieldI].absTol << ")"
+                    << ", rel tol = " << relative
+                    << " (" << residualControl_[fieldI].relTol << ")"
+                    << endl;
+            }
+        }
+    }
+
+    return achieved;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::pimpleControl::pimpleControl(fvMesh& mesh)
+:
+    solutionControl(mesh, "PIMPLE"),
+    nOuterCorr_(0),
+    nCorr_(0),
+    corr_(0),
+    turbOnFinalIterOnly_(true)
+{
+    read();
+
+    if (nOuterCorr_ > 1)
+    {
+        Info<< nl;
+        if (!residualControl_.empty())
+        {
+            Info<< algorithmName_ << ": max iterations = " << nOuterCorr_
+                << endl;
+            forAll(residualControl_, i)
+            {
+                Info<< "    field " << residualControl_[i].name << token::TAB
+                    << ": relTol " << residualControl_[i].relTol
+                    << ", tolerance " << residualControl_[i].absTol
+                    << nl;
+            }
+            Info<< endl;
+        }
+        else
+        {
+            Info<< algorithmName_ << ": no residual control data found. " << nl
+                << "Calculations will employ " << nOuterCorr_
+                << " corrector loops" << nl << endl;
+        }
+    }
+    else
+    {
+        Info<< nl << algorithmName_ << ": Operating solver in PISO mode" << nl
+            << endl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::pimpleControl::~pimpleControl()
+{}
+
+
+// ************************************************************************* //

src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.H

@@ -0,0 +1,158 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::pimpleControl
+
+Description
+    PIMPLE control class to supply convergence information/checks for
+    the PIMPLE loop.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef pimpleControl_H
+#define pimpleControl_H
+
+#include "solutionControl.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class pimpleControl Declaration
+\*---------------------------------------------------------------------------*/
+
+class pimpleControl
+:
+    public solutionControl
+{
+protected:
+
+    // Protected data
+
+        // Solution controls
+
+            //- Maximum number of PIMPLE correctors
+            label nOuterCorr_;
+
+            //- Maximum number of PISO correctors
+            label nCorr_;
+
+            //- Current PIMPLE corrector
+            label corr_;
+
+            //- Flag to indicate whether to only solve turbulence on final iter
+            bool turbOnFinalIterOnly_;
+
+
+    // Protected Member Functions
+
+        //- Read constrols from fvSolution dictionary
+        virtual void read();
+
+        //- Return true if all convergence checks are satified
+        virtual bool criteriaSatisfied();
+
+        //- Disallow default bitwise copy construct
+        pimpleControl(const pimpleControl&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const pimpleControl&);
+
+
+public:
+
+
+    // Static Data Members
+
+        //- Run-time type information
+        TypeName("pimpleControl");
+
+
+    // Constructors
+
+        //- Construct from mesh
+        pimpleControl(fvMesh& mesh);
+
+
+    //- Destructor
+    virtual ~pimpleControl();
+
+
+    // Member Functions
+
+        // Access
+
+            //- Current corrector index
+            inline label corr() const;
+
+            //- Maximum number of PIMPLE correctors
+            inline label nOuterCorr() const;
+
+            //- Maximum number of PISO correctors
+            inline label nCorr() const;
+
+
+        // Solution control
+
+            //- Loop start
+            inline bool start();
+
+            //- Loop loop
+            inline bool loop();
+
+            //- Helper function to identify final PIMPLE (outer) iteration
+            inline bool finalIter() const;
+
+            //- Helper function to identify final inner iteration
+            inline bool finalInnerIter
+            (
+                const label corr,
+                const label nonOrth
+            ) const;
+
+            //- Helper function to identify whether to solve for turbulence
+            inline bool turbCorr() const;
+
+
+    // Member Operators
+
+        void operator++(int);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "pimpleControlI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H

@@ -0,0 +1,130 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline Foam::label Foam::pimpleControl::corr() const
+{
+    return corr_;
+}
+
+
+inline Foam::label Foam::pimpleControl::nOuterCorr() const
+{
+    return nOuterCorr_;
+}
+
+
+inline Foam::label Foam::pimpleControl::nCorr() const
+{
+    return nCorr_;
+}
+
+
+inline bool Foam::pimpleControl::start()
+{
+    corr_ = 0;
+
+    return true;
+}
+
+
+inline bool Foam::pimpleControl::loop()
+{
+    read();
+
+    if (criteriaSatisfied())
+    {
+        Info<< algorithmName_ << ": converged in " << corr_ << " iterations"
+            << endl;
+        return false;
+    }
+    else
+    {
+        if (finalIter())
+        {
+            mesh_.data::add("finalIteration", true);
+        }
+
+        if (corr_ < nOuterCorr_)
+        {
+            if (nOuterCorr_ != 1)
+            {
+                Info<< algorithmName_ << ": iteration " << corr_ + 1 << endl;
+            }
+
+            return true;
+        }
+        else
+        {
+            if ((!residualControl_.empty()) && (nOuterCorr_ != 1))
+            {
+                Info<< algorithmName_ << ": not converged within "
+                    << nOuterCorr_ << " iterations" << endl;
+            }
+
+            return false;
+        }
+    }
+}
+
+
+inline bool Foam::pimpleControl::finalIter() const
+{
+    return corr_ == nOuterCorr_-1;
+}
+
+
+inline bool Foam::pimpleControl::finalInnerIter
+(
+    const label corr,
+    const label nonOrth
+) const
+{
+    return
+        corr_ == nOuterCorr_-1
+     && corr == nCorr_-1
+     && nonOrth == nNonOrthCorr_;
+}
+
+
+inline bool Foam::pimpleControl::turbCorr() const
+{
+    return !turbOnFinalIterOnly_ || finalIter();
+}
+
+
+inline void Foam::pimpleControl::operator++(int)
+{
+    if (finalIter())
+    {
+        mesh_.data::remove("finalIteration");
+    }
+
+    corr_++;
+}
+
+
+// ************************************************************************* //

src/finiteVolume/cfdTools/general/solutionControl/simpleControl/simpleControl.C

@@ -0,0 +1,119 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "simpleControl.H"
+#include "Time.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(simpleControl, 0);
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+void Foam::simpleControl::read()
+{
+    solutionControl::read(true);
+}
+
+
+bool Foam::simpleControl::criteriaSatisfied()
+{
+    if (residualControl_.empty())
+    {
+        return false;
+    }
+
+    bool achieved = true;
+    const dictionary& solverDict = mesh_.solverPerformanceDict();
+
+    forAllConstIter(dictionary, solverDict, iter)
+    {
+        const word& variableName = iter().keyword();
+        label fieldI = applyToField(variableName);
+        if (fieldI != -1)
+        {
+            const List<lduMatrix::solverPerformance> sp(iter().stream());
+            const scalar residual = sp.first().initialResidual();
+
+            bool absCheck = residual < residualControl_[fieldI].absTol;
+            achieved = achieved && absCheck;
+
+            if (debug)
+            {
+                Info<< algorithmName_ << " solution statistics:" << endl;
+
+                Info<< "    " << variableName << ": tolerance = " << residual
+                    << " (" << residualControl_[fieldI].absTol << ")"
+                    << endl;
+            }
+        }
+    }
+
+    return achieved;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::simpleControl::simpleControl(fvMesh& mesh)
+:
+    solutionControl(mesh, "SIMPLE"),
+    initialised_(false)
+{
+    read();
+
+    Info<< nl;
+
+    if (residualControl_.size() > 0)
+    {
+        Info<< algorithmName_ << ": convergence criteria" << nl;
+        forAll(residualControl_, i)
+        {
+            Info<< "    field " << residualControl_[i].name << token::TAB
+                << " tolerance " << residualControl_[i].absTol
+                << nl;
+        }
+        Info<< endl;
+    }
+    else
+    {
+        Info<< algorithmName_ << ": no convergence criteria found. "
+            << "Calculations will run for " << mesh_.time().endTime().value()
+            << " steps." << nl << endl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::simpleControl::~simpleControl()
+{}
+
+
+// ************************************************************************* //

src/finiteVolume/cfdTools/general/solutionControl/simpleControl/simpleControl.H

@@ -22,18 +22,18 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::pimpleLoop
+    Foam::simpleControl
 
 Description
-    PIMPLE loop class to formalise the iteration and automate the handling
-    of the "finalIteration" mesh data entry.
+    SIMPLE control class to supply convergence information/checks for
+    the SIMPLE loop.
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef pimpleLoop_H
-#define pimpleLoop_H
+#ifndef simpleControl_H
+#define simpleControl_H
 
-#include "fvMesh.H"
+#include "solutionControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -41,79 +41,62 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                         Class pimpleLoop Declaration
+                         Class simpleControl Declaration
 \*---------------------------------------------------------------------------*/
 
-class pimpleLoop
+class simpleControl
+:
+    public solutionControl
 {
-    // Private data
 
-        //- Reference to the mesh
-        fvMesh& mesh_;
+protected:
 
-        //- Number of PIMPLE correctors
-        const int nCorr_;
+    // Protected Data
 
-        //- Current PIMPLE corrector
-        int corr_;
+        //- Initialised flag
+        bool initialised_;
 
 
-    // Private Member Functions
+    // Protected Member Functions
+
+        //- Read constrols from fvSolution dictionary
+        void read();
+
+        //- Return true if all convergence checks are satified
+        bool criteriaSatisfied();
 
         //- Disallow default bitwise copy construct
-        pimpleLoop(const pimpleLoop&);
+        simpleControl(const simpleControl&);
 
         //- Disallow default bitwise assignment
-        void operator=(const pimpleLoop&);
+        void operator=(const simpleControl&);
 
 
 public:
 
+
+    // Static Data Members
+
+        //- Run-time type information
+        TypeName("simpleControl");
+
+
     // Constructors
 
-        //- Construct from components
-        pimpleLoop(fvMesh& mesh, const int nCorr)
-        :
-            mesh_(mesh),
-            nCorr_(nCorr),
-            corr_(0)
-        {}
+        //- Construct from mesh
+        simpleControl(fvMesh& mesh);
 
 
     //- Destructor
-    ~pimpleLoop()
-    {}
+    virtual ~simpleControl();
 
 
     // Member Functions
 
-        bool loop()
-        {
-            if (finalIter())
-            {
-                mesh_.data::add("finalIteration", true);
-            }
+        // Solution control
 
-            return corr_ < nCorr_;
-        }
-
-        bool finalIter() const
-        {
-            return corr_ == nCorr_-1;
-        }
-
-
-    // Member Operators
-
-        void operator++(int)
-        {
-            if (finalIter())
-            {
-                mesh_.data::remove("finalIteration");
-            }
-
-            corr_++;
-        }
+            //- Loop loop
+            inline bool loop();
 };
 
 
@@ -123,6 +106,10 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#include "simpleControlI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //

src/finiteVolume/cfdTools/general/solutionControl/simpleControl/simpleControlI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,20 +23,34 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "residualControlFunctionObject.H"
+#include "Time.H"
 
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-namespace Foam
+inline bool Foam::simpleControl::loop()
 {
-    defineNamedTemplateTypeNameAndDebug(residualControlFunctionObject, 0);
+    read();
 
-    addToRunTimeSelectionTable
-    (
-        functionObject,
-        residualControlFunctionObject,
-        dictionary
-    );
+    Time& time = const_cast<Time&>(mesh_.time());
+
+    if (initialised_)
+    {
+        if (criteriaSatisfied())
+        {
+            Info<< nl << algorithmName_ << " solution converged in "
+                << time.timeName() << " iterations" << nl << endl;
+
+            // Set to finalise calculation
+            time.writeAndEnd();
+        }
+    }
+    else
+    {
+        initialised_ = true;
+    }
+
+    return time.loop();
 }
 
+
 // ************************************************************************* //

src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.C

@@ -0,0 +1,124 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "solutionControl.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(solutionControl, 0);
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+void Foam::solutionControl::read(const bool absTolOnly)
+{
+    const dictionary& solnDict = this->dict();
+
+    // Read solution controls
+    nNonOrthCorr_ =
+        solnDict.lookupOrDefault<label>("nNonOrthogonalCorrectors", 0);
+    momentumPredictor_ = solnDict.lookupOrDefault("momentumPredictor", true);
+    transonic_ = solnDict.lookupOrDefault("transonic", false);
+
+    // Read residual information
+    const dictionary residualDict(solnDict.subOrEmptyDict("residualControl"));
+    DynamicList<fieldData> data(residualDict.toc().size());
+    wordHashSet fieldNames;
+
+    forAllConstIter(dictionary, residualDict, iter)
+    {
+        if (fieldNames.insert(iter().keyword()))
+        {
+            fieldData fd;
+            fd.name = iter().keyword().c_str();
+
+            if (absTolOnly)
+            {
+                fd.absTol = readScalar(residualDict.lookup(iter().keyword()));
+                fd.relTol = -1;
+                fd.initialResidual = -1;
+            }
+            else
+            {
+                if (iter().isDict())
+                {
+                    const dictionary& fieldDict(iter().dict());
+                    fd.absTol = readScalar(fieldDict.lookup("tolerance"));
+                    fd.relTol = readScalar(fieldDict.lookup("relTol"));
+                    fd.initialResidual = 0.0;
+                }
+                else
+                {
+                    FatalErrorIn("bool Foam::solutionControl::read()")
+                        << "Residual data for " << iter().keyword()
+                        << " must be specified as a dictionary";
+                }
+            }
+
+            data.append(fd);
+        }
+    }
+
+    residualControl_.transfer(data);
+}
+
+
+Foam::label Foam::solutionControl::applyToField(const word& fieldName) const
+{
+    forAll(residualControl_, i)
+    {
+        if (residualControl_[i].name.match(fieldName))
+        {
+            return i;
+        }
+    }
+
+    return -1;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::solutionControl::solutionControl(fvMesh& mesh, const word& algorithmName)
+:
+    mesh_(mesh),
+    residualControl_(),
+    algorithmName_(algorithmName),
+    nNonOrthCorr_(0),
+    momentumPredictor_(true),
+    transonic_(false)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::solutionControl::~solutionControl()
+{}
+
+
+// ************************************************************************* //

src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.H

@@ -0,0 +1,155 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::solutionControl
+
+Description
+    Base class for solution control classes
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef solutionControl_H
+#define solutionControl_H
+
+#include "fvMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                  Class solutionControl Declaration
+\*---------------------------------------------------------------------------*/
+
+class solutionControl
+{
+public:
+
+    struct fieldData
+    {
+        wordRe name;
+        scalar absTol;
+        scalar relTol;
+        scalar initialResidual;
+    };
+
+
+protected:
+
+    // Protected data
+
+        //- Reference to the mesh database
+        fvMesh& mesh_;
+
+        //- List of residual data per field
+        List<fieldData> residualControl_;
+
+        //- The dictionary name, e.g. SIMPLE, PIMPLE
+        const word algorithmName_;
+
+
+        // Solution controls
+
+            //- Maximum number of non-orthogonal correctors
+            label nNonOrthCorr_;
+
+            //- Flag to indicate to solve for momentum
+            bool momentumPredictor_;
+
+            //- Flag to indictae to solve using transonic algorithm
+            bool transonic_;
+
+
+    // Protected Member Functions
+
+        //- Read constrols from fvSolution dictionary
+        virtual void read(const bool absTolOnly);
+
+        //- Return index of field in residualControl_ if present
+        virtual label applyToField(const word& fieldName) const;
+
+        //- Return true if all convergence checks are satified
+        virtual bool criteriaSatisfied() = 0;
+
+        //- Disallow default bitwise copy construct
+        solutionControl(const solutionControl&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const solutionControl&);
+
+
+public:
+
+
+    // Static Data Members
+
+        //- Run-time type information
+        TypeName("solutionControl");
+
+
+    // Constructors
+
+        //- Construct from mesh
+        solutionControl(fvMesh& mesh, const word& algorithmName);
+
+
+    //- Destructor
+    virtual ~solutionControl();
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return the solution dictionary
+            inline const dictionary& dict() const;
+
+
+        // Solution control
+
+            //- Maximum number of non-orthogonal correctors
+            inline label nNonOrthCorr() const;
+
+            //- Flag to indicate to solve for momentum
+            inline bool momentumPredictor() const;
+
+            //- Flag to indictae to solve using transonic algorithm
+            inline bool transonic() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "solutionControlI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControlI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -21,29 +21,32 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
-Typedef
-    Foam::IOresidualControl
-
-Description
-    Instance of the generic IOOutputFilter for residualControl.
-
 \*---------------------------------------------------------------------------*/
 
-#ifndef IOresidualControl_H
-#define IOresidualControl_H
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-#include "residualControl.H"
-#include "IOOutputFilter.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
+inline const Foam::dictionary& Foam::solutionControl::dict() const
 {
-    typedef IOOutputFilter<residualControl> IOresidualControl;
+    return mesh_.solutionDict().subDict(algorithmName_);
 }
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#endif
+inline Foam::label Foam::solutionControl::nNonOrthCorr() const
+{
+    return nNonOrthCorr_;
+}
+
+
+inline bool Foam::solutionControl::momentumPredictor() const
+{
+    return momentumPredictor_;
+}
+
+
+inline bool Foam::solutionControl::transonic() const
+{
+    return transonic_;
+}
+
 
 // ************************************************************************* //

src/finiteVolume/fields/fvPatchFields/constraint/empty/emptyFvPatchField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -138,18 +138,21 @@ emptyFvPatchField<Type>::emptyFvPatchField
 template<class Type>
 void emptyFvPatchField<Type>::updateCoeffs()
 {
-    if
-    (
-        this->patch().patch().size()
-      % this->dimensionedInternalField().mesh().nCells()
-    )
-    {
-        FatalErrorIn("emptyFvPatchField<Type>::updateCoeffs()")
-            << "This mesh contains patches of type empty but is not 1D or 2D\n"
-               "    by virtue of the fact that the number of faces of this\n"
-               "    empty patch is not divisible by the number of cells."
-            << exit(FatalError);
-    }
+    //- Check moved to checkMesh. Test here breaks down if multiple empty
+    //  patches.
+    //if
+    //(
+    //    this->patch().patch().size()
+    //  % this->dimensionedInternalField().mesh().nCells()
+    //)
+    //{
+    //    FatalErrorIn("emptyFvPatchField<Type>::updateCoeffs()")
+    //        << "This mesh contains patches of type empty but is not"
+    //        << "1D or 2D\n"
+    //           "    by virtue of the fact that the number of faces of this\n"
+    //           "    empty patch is not divisible by the number of cells."
+    //        << exit(FatalError);
+    //}
 }
 
 

src/finiteVolume/fields/fvPatchFields/derived/codedFixedValue/codedFixedValueFvPatchField.C

@@ -59,60 +59,66 @@ void* Foam::codedFixedValueFvPatchField<Type>::loadLibrary
     const fileName& libPath,
     const string& globalFuncName,
     const dictionary& contextDict
-)
+) const
 {
     void* lib = 0;
 
     // avoid compilation by loading an existing library
-    if (!libPath.empty() && dlLibraryTable::open(libPath, false))
+    if (!libPath.empty())
     {
-        lib = dlLibraryTable::findLibrary(libPath);
+        dlLibraryTable& libs = const_cast<Time&>(this->db().time()).libs();
 
-        // verify the loaded version and unload if needed
-        if (lib)
+        if (libs.open(libPath, false))
         {
-            // provision for manual execution of code after loading
-            if (dlSymFound(lib, globalFuncName))
-            {
-                loaderFunctionType function =
-                    reinterpret_cast<loaderFunctionType>
-                    (
-                        dlSym(lib, globalFuncName)
-                    );
+            lib = libs.findLibrary(libPath);
 
-                if (function)
+            // verify the loaded version and unload if needed
+            if (lib)
+            {
+                // provision for manual execution of code after loading
+                if (dlSymFound(lib, globalFuncName))
                 {
-                    (*function)(true);    // force load
+                    loaderFunctionType function =
+                        reinterpret_cast<loaderFunctionType>
+                        (
+                            dlSym(lib, globalFuncName)
+                        );
+
+                    if (function)
+                    {
+                        (*function)(true);    // force load
+                    }
+                    else
+                    {
+                        FatalIOErrorIn
+                        (
+                            "codedFixedValueFvPatchField<Type>::"
+                            "updateLibrary()",
+                            contextDict
+                        )   << "Failed looking up symbol " << globalFuncName
+                            << nl << "from " << libPath << exit(FatalIOError);
+                    }
                 }
                 else
                 {
                     FatalIOErrorIn
                     (
-                        "codedFixedValueFvPatchField<Type>::updateLibrary()",
+                        "codedFixedValueFvPatchField<Type>::loadLibrary()",
                         contextDict
                     )   << "Failed looking up symbol " << globalFuncName << nl
                         << "from " << libPath << exit(FatalIOError);
-                }
-            }
-            else
-            {
-                FatalIOErrorIn
-                (
-                    "codedFixedValueFvPatchField<Type>::loadLibrary()",
-                    contextDict
-                )   << "Failed looking up symbol " << globalFuncName << nl
-                    << "from " << libPath << exit(FatalIOError);
 
-                lib = 0;
-                if (!dlLibraryTable::close(libPath, false))
-                {
-                    FatalIOErrorIn
-                    (
-                        "codedFixedValueFvPatchField<Type>::loadLibrary()",
-                        contextDict
-                    )   << "Failed unloading library "
-                        << libPath
-                        << exit(FatalIOError);
+                    lib = 0;
+                    if (!libs.close(libPath, false))
+                    {
+                        FatalIOErrorIn
+                        (
+                            "codedFixedValueFvPatchField<Type>::loadLibrary()",
+                            contextDict
+                        )   << "Failed unloading library "
+                            << libPath
+                            << exit(FatalIOError);
+                    }
                 }
             }
         }
@@ -128,15 +134,19 @@ void Foam::codedFixedValueFvPatchField<Type>::unloadLibrary
     const fileName& libPath,
     const string& globalFuncName,
     const dictionary& contextDict
-)
+) const
 {
     void* lib = 0;
 
-    if (!libPath.empty())
+    if (libPath.empty())
     {
-        lib = dlLibraryTable::findLibrary(libPath);
+        return;
     }
 
+    dlLibraryTable& libs = const_cast<Time&>(this->db().time()).libs();
+
+    lib = libs.findLibrary(libPath);
+
     if (!lib)
     {
         return;
@@ -166,7 +176,7 @@ void Foam::codedFixedValueFvPatchField<Type>::unloadLibrary
         }
     }
 
-    if (!dlLibraryTable::close(libPath, false))
+    if (!libs.close(libPath, false))
     {
         FatalIOErrorIn
         (
@@ -334,7 +344,7 @@ void Foam::codedFixedValueFvPatchField<Type>::updateLibrary() const
 
 
     // the correct library was already loaded => we are done
-    if (dlLibraryTable::findLibrary(libPath))
+    if (const_cast<Time&>(this->db().time()).libs().findLibrary(libPath))
     {
         return;
     }

src/finiteVolume/fields/fvPatchFields/derived/codedFixedValue/codedFixedValueFvPatchField.H

@@ -130,20 +130,20 @@ class codedFixedValueFvPatchField
         );
 
         //- Load specified library and execute globalFuncName(true)
-        static void* loadLibrary
+        void* loadLibrary
         (
             const fileName& libPath,
             const string& globalFuncName,
             const dictionary& contextDict
-        );
+        ) const;
 
         //- Execute globalFuncName(false) and unload specified library
-        static void unloadLibrary
+        void unloadLibrary
         (
             const fileName& libPath,
             const string& globalFuncName,
             const dictionary& contextDict
-        );
+        ) const;
 
         //- Set the rewrite vars controlling the Type
         static void setFieldTemplates(dynamicCode& dynCode);

src/finiteVolume/finiteVolume/snGradSchemes/correctedSnGrad/correctedSnGrad.C

@@ -56,7 +56,7 @@ Foam::fv::correctedSnGrad<Type>::fullGradCorrection
             gradScheme<Type>::New
             (
                 mesh,
-                mesh.gradScheme(vf.name())
+                mesh.gradScheme("grad(" + vf.name() + ')')
             )().grad(vf, "grad(" + vf.name() + ')')
         );
     tssf().rename("snGradCorr(" + vf.name() + ')');
@@ -108,7 +108,7 @@ Foam::fv::correctedSnGrad<Type>::correction
                 gradScheme<typename pTraits<Type>::cmptType>::New
                 (
                     mesh,
-                    mesh.gradScheme(ssf.name())
+                    mesh.gradScheme("grad(" + ssf.name() + ')')
                 )()
                 //gaussGrad<typename pTraits<Type>::cmptType>(mesh)
                .grad(vf.component(cmpt))

src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearUpwind/linearUpwind.H

@@ -111,7 +111,16 @@ public:
                     mesh.gradScheme(gradSchemeName_)
                 )
             )
-        {}
+        {
+            if (!schemeData.eof())
+            {
+                IOWarningIn("linearUpwind(const fvMesh&, Istream&)", schemeData)
+                    << "unexpected additional entries in stream." << nl
+                    << "    Only the name of the gradient scheme in the"
+                       " 'gradSchemes' dictionary should be specified."
+                    << endl;
+            }
+        }
 
         //- Construct from faceFlux and Istream
         linearUpwind
@@ -131,7 +140,16 @@ public:
                     mesh.gradScheme(gradSchemeName_)
                 )
             )
-        {}
+        {
+            if (!schemeData.eof())
+            {
+                IOWarningIn("linearUpwind(const fvMesh&, Istream&)", schemeData)
+                    << "unexpected additional entries in stream." << nl
+                    << "    Only the name of the gradient scheme in the"
+                       " 'gradSchemes' dictionary should be specified."
+                    << endl;
+            }
+        }
 
 
     // Member Functions

src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearUpwind/linearUpwindV.H

@@ -110,7 +110,13 @@ public:
                     mesh.gradScheme(gradSchemeName_)
                 )
             )
-        {}
+        {
+            IOWarningIn("linearUpwindV(const fvMesh&, Istream&)", schemeData)
+                << "unexpected additional entries in stream." << nl
+                << "    Only the name of the gradient scheme in the"
+                   " 'gradSchemes' dictionary should be specified."
+                << endl;
+        }
 
         //- Construct from faceFlux and Istream
         linearUpwindV
@@ -130,7 +136,13 @@ public:
                     mesh.gradScheme(gradSchemeName_)
                 )
             )
-        {}
+        {
+            IOWarningIn("linearUpwindV(const fvMesh&, Istream&)", schemeData)
+                << "unexpected additional entries in stream." << nl
+                << "    Only the name of the gradient scheme in the"
+                   " 'gradSchemes' dictionary should be specified."
+                << endl;
+        }
 
 
     // Member Functions

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

@@ -100,13 +100,7 @@ template<class CloudType>
 template<class TrackData>
 void Foam::KinematicCloud<CloudType>::solve(TrackData& td)
 {
-    if (solution_.transient())
-    {
-        td.cloud().preEvolve();
-
-        evolveCloud(td);
-    }
-    else
+    if (solution_.steadyState())
     {
         td.cloud().storeState();
 
@@ -114,7 +108,21 @@ void Foam::KinematicCloud<CloudType>::solve(TrackData& td)
 
         evolveCloud(td);
 
-        td.cloud().relaxSources(td.cloud().cloudCopy());
+        if (solution_.coupled())
+        {
+            td.cloud().relaxSources(td.cloud().cloudCopy());
+        }
+    }
+    else
+    {
+        td.cloud().preEvolve();
+
+        evolveCloud(td);
+
+        if (solution_.coupled())
+        {
+            td.cloud().scaleSources();
+        }
     }
 
     td.cloud().info();
@@ -258,6 +266,7 @@ void Foam::KinematicCloud<CloudType>::cloudReset(KinematicCloud<CloudType>& c)
     injectionModel_.reset(c.injectionModel_.ptr());
     patchInteractionModel_.reset(c.patchInteractionModel_.ptr());
     postProcessingModel_.reset(c.postProcessingModel_.ptr());
+    surfaceFilmModel_.reset(c.surfaceFilmModel_.ptr());
 
     UIntegrator_.reset(c.UIntegrator_.ptr());
 }
@@ -556,11 +565,23 @@ void Foam::KinematicCloud<CloudType>::relax
 ) const
 {
     const scalar coeff = solution_.relaxCoeff(name);
-
     field = field0 + coeff*(field - field0);
 }
 
 
+template<class CloudType>
+template<class Type>
+void Foam::KinematicCloud<CloudType>::scale
+(
+    DimensionedField<Type, volMesh>& field,
+    const word& name
+) const
+{
+    const scalar coeff = solution_.relaxCoeff(name);
+    field *= coeff;
+}
+
+
 template<class CloudType>
 void Foam::KinematicCloud<CloudType>::relaxSources
 (
@@ -568,6 +589,15 @@ void Foam::KinematicCloud<CloudType>::relaxSources
 )
 {
     this->relax(UTrans_(), cloudOldTime.UTrans(), "U");
+    this->relax(UCoeff_(), cloudOldTime.UCoeff(), "U");
+}
+
+
+template<class CloudType>
+void Foam::KinematicCloud<CloudType>::scaleSources()
+{
+    this->scale(UTrans_(), "U");
+    this->scale(UCoeff_(), "U");
 }
 
 

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H

@@ -498,9 +498,20 @@ public:
                 const word& name
             ) const;
 
+            //- Scale field
+            template<class Type>
+            void scale
+            (
+                DimensionedField<Type, volMesh>& field,
+                const word& name
+            ) const;
+
             //- Apply relaxation to (steady state) cloud sources
             void relaxSources(const KinematicCloud<CloudType>& cloudOldTime);
 
+            //- Apply scaling to (transient) cloud sources
+            void scaleSources();
+
             //- Evolve the cloud
             void evolve();
 

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H

@@ -370,42 +370,22 @@ Foam::KinematicCloud<CloudType>::theta() const
                 false
             ),
             mesh_,
-            dimensionedScalar("zero", dimless, 0.0)
+            dimensionedScalar("zero", dimless, 0.0),
+            zeroGradientFvPatchScalarField::typeName
         )
     );
 
     volScalarField& theta = ttheta();
-    theta.boundaryField() == 0;
-
     forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
     {
         const parcelType& p = iter();
         const label cellI = p.cell();
 
-        if ((p.face() != -1))
-        {
-            const label patchI = p.patch(p.face());
-            if (patchI != -1)
-            {
-                scalarField& thetap = theta.boundaryField()[patchI];
-                const label faceI = p.patchFace(patchI, p.face());
-                thetap[faceI] += p.nParticle()*p.areaP();
-            }
-        }
-
         theta[cellI] += p.nParticle()*p.volume();
     }
 
     theta.internalField() /= mesh_.V();
-
-    forAll(theta.boundaryField(), patchI)
-    {
-        scalarField& thetap = theta.boundaryField()[patchI];
-        if (thetap.size() > 0)
-        {
-            thetap /= mesh_.magSf().boundaryField()[patchI];
-        }
-    }
+    theta.correctBoundaryConditions();
 
     return ttheta;
 }

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

@@ -116,40 +116,47 @@ void Foam::cloudSolution::read()
         dict_.lookup("calcFrequency") >> calcFrequency_;
         dict_.lookup("maxCo") >> maxCo_;
         dict_.lookup("maxTrackTime") >> maxTrackTime_;
-        dict_.subDict("sourceTerms").lookup("resetOnStartup")
-            >> resetSourcesOnStartup_;
+
+        if (coupled_)
+        {
+            dict_.subDict("sourceTerms").lookup("resetOnStartup")
+                >> resetSourcesOnStartup_;
+        }
     }
 
-    const dictionary&
-        schemesDict(dict_.subDict("sourceTerms").subDict("schemes"));
-
-    wordList vars(schemesDict.toc());
-    schemes_.setSize(vars.size());
-    forAll(vars, i)
+    if (coupled_)
     {
-        // read solution variable name
-        schemes_[i].first() = vars[i];
+        const dictionary&
+            schemesDict(dict_.subDict("sourceTerms").subDict("schemes"));
 
-        // set semi-implicit (1) explicit (0) flag
-        Istream& is = schemesDict.lookup(vars[i]);
-        const word scheme(is);
-        if (scheme == "semiImplicit")
+        wordList vars(schemesDict.toc());
+        schemes_.setSize(vars.size());
+        forAll(vars, i)
         {
-            schemes_[i].second().first() = true;
-        }
-        else if (scheme == "explicit")
-        {
-            schemes_[i].second().first() = false;
-        }
-        else
-        {
-            FatalErrorIn("void cloudSolution::read()")
-                << "Invalid scheme " << scheme << ". Valid schemes are "
-                << "explicit and semiImplicit" << exit(FatalError);
-        }
+            // read solution variable name
+            schemes_[i].first() = vars[i];
 
-        // read under-relaxation factor
-        is  >> schemes_[i].second().second();
+            // set semi-implicit (1) explicit (0) flag
+            Istream& is = schemesDict.lookup(vars[i]);
+            const word scheme(is);
+            if (scheme == "semiImplicit")
+            {
+                schemes_[i].second().first() = true;
+            }
+            else if (scheme == "explicit")
+            {
+                schemes_[i].second().first() = false;
+            }
+            else
+            {
+                FatalErrorIn("void cloudSolution::read()")
+                    << "Invalid scheme " << scheme << ". Valid schemes are "
+                    << "explicit and semiImplicit" << exit(FatalError);
+            }
+
+            // read under-relaxation factor
+            is  >> schemes_[i].second().second();
+        }
     }
 }
 

src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C

@@ -295,10 +295,10 @@ void Foam::ReactingCloud<CloudType>::relaxSources
     const ReactingCloud<CloudType>& cloudOldTime
 )
 {
-    typedef DimensionedField<scalar, volMesh> dsfType;
-
     CloudType::relaxSources(cloudOldTime);
 
+    typedef DimensionedField<scalar, volMesh> dsfType;
+
     forAll(rhoTrans_, fieldI)
     {
         dsfType& rhoT = rhoTrans_[fieldI];
@@ -308,6 +308,21 @@ void Foam::ReactingCloud<CloudType>::relaxSources
 }
 
 
+template<class CloudType>
+void Foam::ReactingCloud<CloudType>::scaleSources()
+{
+    CloudType::scaleSources();
+
+    typedef DimensionedField<scalar, volMesh> dsfType;
+
+    forAll(rhoTrans_, fieldI)
+    {
+        dsfType& rhoT = rhoTrans_[fieldI];
+        this->scale(rhoT, "rho");
+    }
+}
+
+
 template<class CloudType>
 void Foam::ReactingCloud<CloudType>::evolve()
 {

src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.H

@@ -287,6 +287,9 @@ public:
             //- Apply relaxation to (steady state) cloud sources
             void relaxSources(const ReactingCloud<CloudType>& cloudOldTime);
 
+            //- Apply scaling to (transient) cloud sources
+            void scaleSources();
+
             //- Evolve the cloud
             void evolve();
 

src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C

@@ -290,6 +290,17 @@ void Foam::ThermoCloud<CloudType>::relaxSources
     CloudType::relaxSources(cloudOldTime);
 
     this->relax(hsTrans_(), cloudOldTime.hsTrans(), "hs");
+    this->relax(hsCoeff_(), cloudOldTime.hsCoeff(), "hs");
+}
+
+
+template<class CloudType>
+void Foam::ThermoCloud<CloudType>::scaleSources()
+{
+    CloudType::scaleSources();
+
+    this->scale(hsTrans_(), "hs");
+    this->scale(hsCoeff_(), "hs");
 }
 
 

src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H

@@ -301,6 +301,9 @@ public:
             //- Apply relaxation to (steady state) cloud sources
             void relaxSources(const ThermoCloud<CloudType>& cloudOldTime);
 
+            //- Apply scaling to (transient) cloud sources
+            void scaleSources();
+
             //- Evolve the cloud
             void evolve();
 

src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H

@@ -518,7 +518,7 @@ inline Foam::scalar Foam::KinematicParcel<ParcelType>::Re
     const scalar muc
 ) const
 {
-    return rhoc*mag(U - Uc_)*d/muc;
+    return rhoc*mag(U - Uc_)*d/(muc + ROOTVSMALL);
 }
 
 

src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C

@@ -63,7 +63,7 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::CpEff
 
 template<class ParcelType>
 template<class TrackData>
-Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::HEff
+Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::HsEff
 (
     TrackData& td,
     const scalar p,
@@ -74,9 +74,9 @@ Foam::scalar Foam::ReactingMultiphaseParcel<ParcelType>::HEff
 ) const
 {
     return
-        this->Y_[GAS]*td.cloud().composition().H(idG, YGas_, p, T)
-      + this->Y_[LIQ]*td.cloud().composition().H(idL, YLiquid_, p, T)
-      + this->Y_[SLD]*td.cloud().composition().H(idS, YSolid_, p, T);
+        this->Y_[GAS]*td.cloud().composition().Hs(idG, YGas_, p, T)
+      + this->Y_[LIQ]*td.cloud().composition().Hs(idL, YLiquid_, p, T)
+      + this->Y_[SLD]*td.cloud().composition().Hs(idS, YSolid_, p, T);
 }
 
 
@@ -326,8 +326,6 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         updateMassFractions(mass0, dMassGas, dMassLiquid, dMassSolid);
 
 
-
-
     // Heat transfer
     // ~~~~~~~~~~~~~
 
@@ -384,33 +382,41 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         // Transfer mass lost from particle to carrier mass source
         forAll(YGas_, i)
         {
+            scalar dm = np0*dMassGas[i];
             label gid = composition.localToGlobalCarrierId(GAS, i);
-            td.cloud().rhoTrans(gid)[cellI] += np0*dMassGas[i];
-//            td.cloud().hsTrans()[cellI] +=
-//                np0*dMassGas[i]*composition.carrier().Hs(gid, T0);
+            scalar hs = composition.carrier().Hs(gid, T0);
+            td.cloud().rhoTrans(gid)[cellI] += dm;
+            td.cloud().UTrans()[cellI] += dm*U0;
+            td.cloud().hsTrans()[cellI] += dm*hs;
         }
         forAll(YLiquid_, i)
         {
+            scalar dm = np0*dMassLiquid[i];
             label gid = composition.localToGlobalCarrierId(LIQ, i);
-            td.cloud().rhoTrans(gid)[cellI] += np0*dMassLiquid[i];
-//            td.cloud().hsTrans()[cellI] +=
-//                np0*dMassLiquid[i]*composition.carrier().Hs(gid, T0);
+            scalar hs = composition.carrier().Hs(gid, T0);
+            td.cloud().rhoTrans(gid)[cellI] += dm;
+            td.cloud().UTrans()[cellI] += dm*U0;
+            td.cloud().hsTrans()[cellI] += dm*hs;
         }
 /*
         // No mapping between solid components and carrier phase
         forAll(YSolid_, i)
         {
+            scalar dm = np0*dMassSolid[i];
             label gid = composition.localToGlobalCarrierId(SLD, i);
-            td.cloud().rhoTrans(gid)[cellI] += np0*dMassSolid[i];
-//            td.cloud().hsTrans()[cellI] +=
-//                np0*dMassSolid[i]*composition.carrier().Hs(gid, T0);
+            scalar hs = composition.carrier().Hs(gid, T0);
+            td.cloud().rhoTrans(gid)[cellI] += dm;
+            td.cloud().UTrans()[cellI] += dm*U0;
+            td.cloud().hsTrans()[cellI] += dm*hs;
         }
 */
         forAll(dMassSRCarrier, i)
         {
-            td.cloud().rhoTrans(i)[cellI] += np0*dMassSRCarrier[i];
-//            td.cloud().hsTrans()[cellI] +=
-//                np0*dMassSRCarrier[i]*composition.carrier().Hs(i, T0);
+            scalar dm = np0*dMassSRCarrier[i];
+            scalar hs = composition.carrier().Hs(i, T0);
+            td.cloud().rhoTrans(i)[cellI] += dm;
+            td.cloud().UTrans()[cellI] += dm*U0;
+            td.cloud().hsTrans()[cellI] += dm*hs;
         }
 
         // Update momentum transfer
@@ -430,36 +436,38 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     // Remove the particle when mass falls below minimum threshold
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    if (mass1 < td.cloud().constProps().minParticleMass())
+    if (np0*mass1 < td.cloud().constProps().minParticleMass())
     {
         td.keepParticle = false;
 
         if (td.cloud().solution().coupled())
         {
+            scalar dm = np0*mass1;
+
             // Absorb parcel into carrier phase
             forAll(YGas_, i)
             {
                 label gid = composition.localToGlobalCarrierId(GAS, i);
-                td.cloud().rhoTrans(gid)[cellI] += np0*mass1*YMix[GAS]*YGas_[i];
+                td.cloud().rhoTrans(gid)[cellI] += dm*YMix[GAS]*YGas_[i];
             }
             forAll(YLiquid_, i)
             {
                 label gid = composition.localToGlobalCarrierId(LIQ, i);
-                td.cloud().rhoTrans(gid)[cellI] +=
-                    np0*mass1*YMix[LIQ]*YLiquid_[i];
+                td.cloud().rhoTrans(gid)[cellI] += dm*YMix[LIQ]*YLiquid_[i];
             }
 /*
             // No mapping between solid components and carrier phase
             forAll(YSolid_, i)
             {
                 label gid = composition.localToGlobalCarrierId(SLD, i);
-                td.cloud().rhoTrans(gid)[cellI] +=
-                    np0*mass1*YMix[SLD]*YSolid_[i];
+                td.cloud().rhoTrans(gid)[cellI] += dm*YMix[SLD]*YSolid_[i];
             }
 */
-            td.cloud().UTrans()[cellI] += np0*mass1*U1;
-            td.cloud().hsTrans()[cellI] +=
-                np0*mass1*HEff(td, pc, T1, idG, idL, idS); // using total h
+            td.cloud().UTrans()[cellI] += dm*U1;
+
+            td.cloud().hsTrans()[cellI] += dm*HsEff(td, pc, T1, idG, idL, idS);
+
+            td.cloud().addToMassPhaseChange(dm);
         }
     }
 
@@ -540,28 +548,35 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
 
     Sh -= dMassTot*td.cloud().constProps().LDevol()/dt;
 
-    // Molar average molecular weight of carrier mix
-    const scalar Wc = this->rhoc_*specie::RR*this->Tc_/this->pc_;
-
     // Update molar emissions
-    forAll(dMassDV, i)
+    if (td.cloud().heatTransfer().BirdCorrection())
     {
+        // Molar average molecular weight of carrier mix
+        const scalar Wc =
+            max(SMALL, this->rhoc_*specie::RR*this->Tc_/this->pc_);
+
         // Note: hardcoded gaseous diffusivities for now
         // TODO: add to carrier thermo
         const scalar beta = sqr(cbrt(15.0) + cbrt(15.0));
-        const label id = composition.localToGlobalCarrierId(GAS, i);
-        const scalar Cp = composition.carrier().Cp(id, Ts);
-        const scalar W = composition.carrier().W(id);
-        const scalar Ni = dMassDV[i]/(this->areaS(d)*dt*W);
 
-        // Dab calc'd using API vapour mass diffusivity function
-        const scalar Dab =
-            3.6059e-3*(pow(1.8*Ts, 1.75))*sqrt(1.0/W + 1.0/Wc)/(this->pc_*beta);
+        forAll(dMassDV, i)
+        {
+            const label id = composition.localToGlobalCarrierId(GAS, i);
+            const scalar Cp = composition.carrier().Cp(id, Ts);
+            const scalar W = composition.carrier().W(id);
+            const scalar Ni = dMassDV[i]/(this->areaS(d)*dt*W);
 
-        N += Ni;
-        NCpW += Ni*Cp*W;
-        Cs[id] += Ni*d/(2.0*Dab);
-     }
+            // Dab calc'd using API vapour mass diffusivity function
+            const scalar Dab =
+                3.6059e-3*(pow(1.8*Ts, 1.75))
+               *sqrt(1.0/W + 1.0/Wc)
+               /(this->pc_*beta);
+
+            N += Ni;
+            NCpW += Ni*Cp*W;
+            Cs[id] += Ni*d/(2.0*Dab);
+        }
+    }
 }
 
 

src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.H

@@ -134,9 +134,9 @@ private:
             const label idS
         ) const;
 
-        //- Return the mixture effective enthalpy
+        //- Return the mixture effective sensible enthalpy
         template<class TrackData>
-        scalar HEff
+        scalar HsEff
         (
             TrackData& td,
             const scalar p,

src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C

@@ -213,9 +213,17 @@ void Foam::ReactingParcel<ParcelType>::correctSurfaceValues
         sumYiCbrtW += Ys[i]*cbrtW;
     }
 
+    Cps = max(Cps, ROOTVSMALL);
+
     rhos *= pc_/(specie::RR*T);
+    rhos = max(rhos, ROOTVSMALL);
+
     mus /= sumYiSqrtW;
+    mus = max(mus, ROOTVSMALL);
+
     kappas /= sumYiCbrtW;
+    kappas = max(kappas, ROOTVSMALL);
+
     Prs = Cps*mus/kappas;
 }
 
@@ -335,7 +343,9 @@ void Foam::ReactingParcel<ParcelType>::calc
     Res = this->Re(U0, d0, rhos, mus);
 
     // Update particle component mass and mass fractions
-    scalar mass1 = updateMassFraction(mass0, dMassPC, Y_);
+    scalarField dMass(dMassPC);
+
+    scalar mass1 = updateMassFraction(mass0, dMass, Y_);
 
 
     // Heat transfer
@@ -385,19 +395,21 @@ void Foam::ReactingParcel<ParcelType>::calc
             Spu
         );
 
-    dUTrans += 0.5*(mass0 - mass1)*(U0 + U1);
 
     // Accumulate carrier phase source terms
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     if (td.cloud().solution().coupled())
     {
-        // Transfer mass lost from particle to carrier mass source
-        forAll(dMassPC, i)
+        // Transfer mass lost to carrier mass and enthalpy sources
+        forAll(dMass, i)
         {
+            scalar dm = np0*dMass[i];
             label gid = composition.localToGlobalCarrierId(0, i);
-            td.cloud().rhoTrans(gid)[cellI] += np0*dMassPC[i];
-//            td.cloud().hsTrans()[cellI] +=
-//                np0*dMassPC[i]*composition.carrier().Hs(gid, T0);
+            scalar hs = composition.carrier().Hs(gid, T0);
+
+            td.cloud().rhoTrans(gid)[cellI] += dm;
+            td.cloud().UTrans()[cellI] += dm*U0;
+            td.cloud().hsTrans()[cellI] += dm*hs;
         }
 
         // Update momentum transfer
@@ -416,21 +428,27 @@ void Foam::ReactingParcel<ParcelType>::calc
 
     // Remove the particle when mass falls below minimum threshold
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    if (mass1 < td.cloud().constProps().minParticleMass())
+    if (np0*mass1 < td.cloud().constProps().minParticleMass())
     {
         td.keepParticle = false;
 
         if (td.cloud().solution().coupled())
         {
+            scalar dm = np0*mass1;
+
             // Absorb parcel into carrier phase
             forAll(Y_, i)
             {
+                scalar dmi = dm*Y_[i];
                 label gid = composition.localToGlobalCarrierId(0, i);
-                td.cloud().rhoTrans(gid)[cellI] += np0*mass1*Y_[i];
+                scalar hs = composition.carrier().Hs(gid, T1);
+
+                td.cloud().rhoTrans(gid)[cellI] += dmi;
+                td.cloud().hsTrans()[cellI] += dmi*hs;
             }
-            td.cloud().UTrans()[cellI] += np0*mass1*U1;
-            td.cloud().hsTrans()[cellI] +=
-                np0*mass1*composition.H(0, Y_, pc_, T1);
+            td.cloud().UTrans()[cellI] += dm*U1;
+
+            td.cloud().addToMassPhaseChange(dm);
         }
     }
 
@@ -517,33 +535,44 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
     // Add to cumulative phase change mass
     td.cloud().addToMassPhaseChange(this->nParticle_*dMassTot);
 
-    // Average molecular weight of carrier mix - assumes perfect gas
-    const scalar Wc = this->rhoc_*specie::RR*this->Tc_/this->pc_;
-
-    forAll(YComponents, i)
+    forAll(dMassPC, i)
     {
         const label idc = composition.localToGlobalCarrierId(idPhase, i);
         const label idl = composition.globalIds(idPhase)[i];
 
         const scalar dh = td.cloud().phaseChange().dh(idc, idl, pc_, T);
         Sh -= dMassPC[i]*dh/dt;
+    }
 
-        // Update particle surface thermo properties
-        const scalar Dab =
-            composition.liquids().properties()[idl].D(pc_, Ts, Wc);
 
-        const scalar Cp = composition.carrier().Cp(idc, Ts);
-        const scalar W = composition.carrier().W(idc);
-        const scalar Ni = dMassPC[i]/(this->areaS(d)*dt*W);
+    // Update molar emissions
+    if (td.cloud().heatTransfer().BirdCorrection())
+    {
+        // Average molecular weight of carrier mix - assumes perfect gas
+        const scalar Wc = this->rhoc_*specie::RR*this->Tc_/this->pc_;
 
-        // Molar flux of species coming from the particle (kmol/m^2/s)
-        N += Ni;
 
-        // Sum of Ni*Cpi*Wi of emission species
-        NCpW += Ni*Cp*W;
+        forAll(dMassPC, i)
+        {
+            const label idc = composition.localToGlobalCarrierId(idPhase, i);
+            const label idl = composition.globalIds(idPhase)[i];
 
-        // Concentrations of emission species
-        Cs[idc] += Ni*d/(2.0*Dab);
+            const scalar Cp = composition.carrier().Cp(idc, Ts);
+            const scalar W = composition.carrier().W(idc);
+            const scalar Ni = dMassPC[i]/(this->areaS(d)*dt*W);
+
+            const scalar Dab =
+                composition.liquids().properties()[idl].D(pc_, Ts, Wc);
+
+            // Molar flux of species coming from the particle (kmol/m^2/s)
+            N += Ni;
+
+            // Sum of Ni*Cpi*Wi of emission species
+            NCpW += Ni*Cp*W;
+
+            // Concentrations of emission species
+            Cs[idc] += Ni*d/(2.0*Dab);
+        }
     }
 }
 

src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.C

@@ -153,7 +153,10 @@ void Foam::ThermoParcel<ParcelType>::calcSurfaceValues
     mus = td.muInterp().interpolate(this->position(), tetIs)/TRatio;
 
     Pr = td.cloud().constProps().Pr();
+    Pr = max(ROOTVSMALL, Pr);
+
     kappas = Cpc_*mus/Pr;
+    kappas = max(ROOTVSMALL, kappas);
 }
 
 

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

@@ -148,7 +148,8 @@ void Foam::SurfaceFilmModel<CloudType>::inject(TrackData& td)
     const labelList& filmPatches = filmModel.intCoupledPatchIDs();
     const labelList& primaryPatches = filmModel.primaryPatchIDs();
 
-    const polyBoundaryMesh& pbm = this->owner().mesh().boundaryMesh();
+    const fvMesh& mesh = this->owner().mesh();
+    const polyBoundaryMesh& pbm = mesh.boundaryMesh();
 
     forAll(filmPatches, i)
     {
@@ -163,6 +164,10 @@ void Foam::SurfaceFilmModel<CloudType>::inject(TrackData& td)
 
         cacheFilmFields(filmPatchI, primaryPatchI, distMap, filmModel);
 
+        const vectorField& Cf = mesh.C().boundaryField()[primaryPatchI];
+        const vectorField& Sf = mesh.Sf().boundaryField()[primaryPatchI];
+        const scalarField& magSf = mesh.magSf().boundaryField()[primaryPatchI];
+
         forAll(injectorCellsPatch, j)
         {
             if (diameterParcelPatch_[j] > 0)
@@ -179,7 +184,15 @@ void Foam::SurfaceFilmModel<CloudType>::inject(TrackData& td)
                 const label tetFaceI = this->owner().mesh().cells()[cellI][0];
                 const label tetPtI = 1;
 
-                const point& pos = this->owner().mesh().C()[cellI];
+//                const point& pos = this->owner().mesh().C()[cellI];
+
+                const scalar offset =
+                    max
+                    (
+                        diameterParcelPatch_[j],
+                        deltaFilmPatch_[primaryPatchI][j]
+                    );
+                const point pos = Cf[j] - 1.1*offset*Sf[j]/magSf[j];
 
                 // Create a new parcel
                 parcelType* pPtr =
@@ -217,11 +230,11 @@ void Foam::SurfaceFilmModel<CloudType>::cacheFilmFields
     const regionModels::surfaceFilmModels::surfaceFilmModel& filmModel
 )
 {
-    massParcelPatch_ = filmModel.massForPrimary().boundaryField()[filmPatchI];
+    massParcelPatch_ = filmModel.cloudMassTrans().boundaryField()[filmPatchI];
     distMap.distribute(massParcelPatch_);
 
     diameterParcelPatch_ =
-        filmModel.diametersForPrimary().boundaryField()[filmPatchI];
+        filmModel.cloudDiameterTrans().boundaryField()[filmPatchI];
     distMap.distribute(diameterParcelPatch_);
 
     UFilmPatch_ = filmModel.Us().boundaryField()[filmPatchI];

src/lagrangian/intermediate/submodels/Thermodynamic/SurfaceFilmModel/ThermoSurfaceFilm/ThermoSurfaceFilm.C

@@ -391,8 +391,8 @@ void Foam::ThermoSurfaceFilm<CloudType>::splashInteraction
     const scalar dBarSplash = 1/cbrt(6.0)*cbrt(mRatio/Ns)*d + ROOTVSMALL;
 
     // cumulative diameter splash distribution
-    const scalar dMax = cbrt(mRatio)*d;
-    const scalar dMin = 0.001*dMax;
+    const scalar dMax = 0.9*cbrt(mRatio)*d;
+    const scalar dMin = 0.1*dMax;
     const scalar K = exp(-dMin/dBarSplash) - exp(-dMax/dBarSplash);
 
     // surface energy of secondary parcels [J]
@@ -437,7 +437,7 @@ void Foam::ThermoSurfaceFilm<CloudType>::splashInteraction
 
     // magnitude of the normal velocity of the first splashed parcel
     const scalar magUns0 =
-        sqrt(2.0*parcelsPerSplash_*EKs/mSplash/(1 + coeff1/sqr(coeff2)));
+        sqrt(2.0*parcelsPerSplash_*EKs/mSplash/(1.0 + coeff1/sqr(coeff2)));
 
     // Set splashed parcel properties
     forAll(dNew, i)
@@ -467,7 +467,7 @@ void Foam::ThermoSurfaceFilm<CloudType>::splashInteraction
 
         // Apply correction to velocity for 2-D cases
         meshTools::constrainDirection(mesh, mesh.solutionD(), pPtr->U());
-
+Info<< "NEW PARTICLE: " << *pPtr << endl;
         // Add the new parcel
         this->owner().addParticle(pPtr);
 

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

@@ -2615,9 +2615,12 @@ void Foam::autoLayerDriver::addLayers
             extrudeStatus
         );
 
-        Info<< "Extruding " << countExtrusion(pp, extrudeStatus)
-            << " out of " << returnReduce(pp().size(), sumOp<label>())
-            << " faces. Removed extrusion at " << nTotChanged << " faces."
+        label nExtruded = countExtrusion(pp, extrudeStatus);
+        label nTotFaces = returnReduce(pp().size(), sumOp<label>());
+        Info<< "Extruding " << nExtruded
+            << " out of " << nTotFaces
+            << " faces (" << 100.0*nExtruded/nTotFaces << "%)."
+            << " Removed extrusion at " << nTotChanged << " faces."
             << endl;
 
         if (nTotChanged == 0)

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

@@ -722,10 +722,15 @@ void Foam::autoSnapDriver::preSmoothPatch
     if (debug)
     {
         const_cast<Time&>(mesh.time())++;
-        Pout<< "Writing patch smoothed mesh to time " << meshRefiner_.timeName()
-            << endl;
-
-        mesh.write();
+        Pout<< "Writing patch smoothed mesh to time "
+            << meshRefiner_.timeName() << '.' << endl;
+        meshRefiner_.write
+        (
+            debug,
+            mesh.time().path()/meshRefiner_.timeName()
+        );
+        Pout<< "Dumped mesh in = "
+            << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
     }
 
     Info<< "Patch points smoothed in = "
@@ -956,25 +961,6 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
         vector(GREAT, GREAT, GREAT)     // null value (note: cannot use VGREAT)
     );
 
-
-    // Check for displacement being outwards.
-    outwardsDisplacement(pp, patchDisp);
-
-    // Set initial distribution of displacement field (on patches) from
-    // patchDisp and make displacement consistent with b.c. on displacement
-    // pointVectorField.
-    meshMover.setDisplacement(patchDisp);
-
-    if (debug)
-    {
-        dumpMove
-        (
-            mesh.time().path()/"patchDisplacement.obj",
-            pp.localPoints(),
-            pp.localPoints() + patchDisp
-        );
-    }
-
     return patchDisp;
 }
 
@@ -1128,8 +1114,8 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::autoSnapDriver::repatchToSurface
     indirectPrimitivePatch& pp = ppPtr();
 
     // Divide surfaces into zoned and unzoned
-    labelList zonedSurfaces = meshRefiner_.surfaces().getNamedSurfaces();
-    labelList unzonedSurfaces = meshRefiner_.surfaces().getUnnamedSurfaces();
+    labelList zonedSurfaces = surfaces.getNamedSurfaces();
+    labelList unzonedSurfaces = surfaces.getUnnamedSurfaces();
 
 
     // Faces that do not move
@@ -1344,19 +1330,6 @@ void Foam::autoSnapDriver::doSnap
         // Pre-smooth patch vertices (so before determining nearest)
         preSmoothPatch(snapParams, nInitErrors, baffles, meshMover);
 
-        if (debug)
-        {
-            Pout<< "Writing patch smoothed mesh to time "
-                << meshRefiner_.timeName() << '.' << endl;
-            meshRefiner_.write
-            (
-                debug,
-                mesh.time().path()/meshRefiner_.timeName()
-            );
-            Pout<< "Dumped mesh in = "
-                << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
-        }
-
 
         for (label iter = 0; iter < nFeatIter; iter++)
         {
@@ -1382,6 +1355,15 @@ void Foam::autoSnapDriver::doSnap
                 );
             }
 
+            // Check for displacement being outwards.
+            outwardsDisplacement(ppPtr(), disp);
+
+            // Set initial distribution of displacement field (on patches)
+            // from patchDisp and make displacement consistent with b.c.
+            // on displacement pointVectorField.
+            meshMover.setDisplacement(disp);
+
+
             if (debug&meshRefinement::OBJINTERSECTIONS)
             {
                 dumpMove
@@ -1467,7 +1449,7 @@ void Foam::autoSnapDriver::doSnap
         meshRefiner_.write
         (
             debug,
-            mesh.time().path()/meshRefiner_.timeName()
+            meshRefiner_.timeName()
         );
     }
 }

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

@@ -117,12 +117,12 @@ class autoSnapDriver
                     const List<pointConstraint>& constraints,
                     vectorField& disp
                 ) const;
-                void calcNearest
-                (
-                    const pointField& points,
-                    vectorField& disp,
-                    vectorField& surfaceNormal
-                ) const;
+                //void calcNearest
+                //(
+                //    const pointField& points,
+                //    vectorField& disp,
+                //    vectorField& surfaceNormal
+                //) const;
                 void calcNearestFace
                 (
                     const label iter,

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

@@ -135,63 +135,6 @@ void Foam::autoSnapDriver::smoothAndConstrain
 }
 
 
-void Foam::autoSnapDriver::calcNearest
-(
-    const pointField& points,
-    vectorField& disp,
-    vectorField& surfaceNormal
-) const
-{
-    const refinementSurfaces& surfaces = meshRefiner_.surfaces();
-
-    // Displacement and orientation per pp face.
-    disp.setSize(points.size());
-    disp = vector::zero;
-    surfaceNormal.setSize(points.size());
-    surfaceNormal = vector::zero;
-
-    {
-        // Divide surfaces into zoned and unzoned
-        labelList zonedSurfaces =
-            meshRefiner_.surfaces().getNamedSurfaces();
-        labelList unzonedSurfaces =
-            meshRefiner_.surfaces().getUnnamedSurfaces();
-
-        {
-            List<pointIndexHit> hitInfo;
-            labelList hitSurface;
-            labelList hitRegion;
-            surfaces.findNearestRegion
-            (
-                unzonedSurfaces,
-                points,
-                sqr(scalarField(points.size(), GREAT)),// sqr of attract dist
-                hitSurface,
-                hitInfo,
-                hitRegion,
-                surfaceNormal
-            );
-
-            forAll(hitInfo, i)
-            {
-                if (hitInfo[i].hit())
-                {
-                    disp[i] =
-                        hitInfo[i].hitPoint()
-                      - points[i];
-                }
-                else
-                {
-                    WarningIn("calcNearest(..)")
-                        << "Did not hit anything from face:" << i
-                        << " at:" << points[i] << endl;
-                }
-            }
-        }
-    }
-}
-
-
 void Foam::autoSnapDriver::calcNearestFace
 (
     const label iter,
@@ -202,6 +145,7 @@ void Foam::autoSnapDriver::calcNearestFace
 ) const
 {
     const fvMesh& mesh = meshRefiner_.mesh();
+    const refinementSurfaces& surfaces = meshRefiner_.surfaces();
 
     // Displacement and orientation per pp face.
     faceDisp.setSize(pp.size());
@@ -209,7 +153,148 @@ void Foam::autoSnapDriver::calcNearestFace
     faceSurfaceNormal.setSize(pp.size());
     faceSurfaceNormal = vector::zero;
 
-    calcNearest(pp.faceCentres(), faceDisp, faceSurfaceNormal);
+    // Divide surfaces into zoned and unzoned
+    labelList zonedSurfaces = surfaces.getNamedSurfaces();
+    labelList unzonedSurfaces = surfaces.getUnnamedSurfaces();
+
+    // Per pp face the current surface snapped to
+    labelList snapSurf(pp.size(), -1);
+
+
+    // Do zoned surfaces
+    // ~~~~~~~~~~~~~~~~~
+    // Zoned faces only attract to corresponding surface
+
+    // Extract faces per zone
+    const wordList& faceZoneNames = surfaces.faceZoneNames();
+
+    forAll(zonedSurfaces, i)
+    {
+        label zoneSurfI = zonedSurfaces[i];
+
+        // Get indices of faces on pp that are also in zone
+        label zoneI = mesh.faceZones().findZoneID(faceZoneNames[zoneSurfI]);
+        if (zoneI == -1)
+        {
+            FatalErrorIn
+            (
+                "autoSnapDriver::calcNearestFace(..)"
+            )   << "Problem. Cannot find zone " << faceZoneNames[zoneSurfI]
+                << exit(FatalError);
+        }
+        const faceZone& fZone = mesh.faceZones()[zoneI];
+        PackedBoolList isZonedFace(mesh.nFaces());
+        forAll(fZone, i)
+        {
+            isZonedFace[fZone[i]] = 1;
+        }
+
+        DynamicList<label> ppFaces(fZone.size());
+        DynamicList<label> meshFaces(fZone.size());
+        forAll(pp.addressing(), i)
+        {
+            if (isZonedFace[pp.addressing()[i]])
+            {
+                snapSurf[i] = zoneSurfI;
+                ppFaces.append(i);
+                meshFaces.append(pp.addressing()[i]);
+            }
+        }
+
+        //Pout<< "For faceZone " << fZone.name()
+        //    << " found " << ppFaces.size() << " out of " << pp.size()
+        //    << endl;
+
+        pointField fc
+        (
+            indirectPrimitivePatch
+            (
+                IndirectList<face>(mesh.faces(), meshFaces),
+                mesh.points()
+            ).faceCentres()
+        );
+
+        List<pointIndexHit> hitInfo;
+        labelList hitSurface;
+        labelList hitRegion;
+        vectorField hitNormal;
+        surfaces.findNearestRegion
+        (
+            labelList(1, zoneSurfI),
+            fc,
+            sqr(scalarField(fc.size(), GREAT)),// sqr of attract dist
+            hitSurface,
+            hitInfo,
+            hitRegion,
+            hitNormal
+        );
+
+        forAll(hitInfo, hitI)
+        {
+            if (hitInfo[hitI].hit())
+            {
+                label faceI = ppFaces[hitI];
+                faceDisp[faceI] = hitInfo[hitI].hitPoint() - fc[hitI];
+                faceSurfaceNormal[faceI] = hitNormal[hitI];
+            }
+        }
+    }
+
+
+    // Do unzoned surfaces
+    // ~~~~~~~~~~~~~~~~~~~
+    // Unzoned faces attract to any unzoned surface
+
+    DynamicList<label> ppFaces(pp.size());
+    DynamicList<label> meshFaces(pp.size());
+    forAll(pp.addressing(), i)
+    {
+        if (snapSurf[i] == -1)
+        {
+            ppFaces.append(i);
+            meshFaces.append(pp.addressing()[i]);
+        }
+    }
+    //Pout<< "Found " << ppFaces.size() << " unzoned faces out of " << pp.size()
+    //    << endl;
+
+    pointField fc
+    (
+        indirectPrimitivePatch
+        (
+            IndirectList<face>(mesh.faces(), meshFaces),
+            mesh.points()
+        ).faceCentres()
+    );
+
+    List<pointIndexHit> hitInfo;
+    labelList hitSurface;
+    labelList hitRegion;
+    vectorField hitNormal;
+    surfaces.findNearestRegion
+    (
+        unzonedSurfaces,
+        fc,
+        sqr(scalarField(fc.size(), GREAT)),// sqr of attract dist
+        hitSurface,
+        hitInfo,
+        hitRegion,
+        hitNormal
+    );
+
+    forAll(hitInfo, hitI)
+    {
+        if (hitInfo[hitI].hit())
+        {
+            label faceI = ppFaces[hitI];
+            faceDisp[faceI] = hitInfo[hitI].hitPoint() - fc[hitI];
+            faceSurfaceNormal[faceI] = hitNormal[hitI];
+        }
+    }
+
+
+    // Determine rotation
+    // ~~~~~~~~~~~~~~~~~~
 
     // Determine rotation axis
     faceRotation.setSize(pp.size());
@@ -243,104 +328,6 @@ void Foam::autoSnapDriver::calcNearestFace
 }
 
 
-void Foam::autoSnapDriver::interpolateFaceToPoint
-(
-    const label iter,
-    const indirectPrimitivePatch& pp,
-    const vectorField& faceSurfaceNormal,
-    const vectorField& faceDisp,
-    const vectorField& faceRotation,
-
-    vectorField& patchDisp,
-    vectorField& patchRotationDisp
-) const
-{
-    const fvMesh& mesh = meshRefiner_.mesh();
-
-    // Displacement due to normal snapping
-    patchDisp.setSize(pp.nPoints());
-    patchDisp = vector::zero;
-    // Displacement due to rotation
-    patchRotationDisp.setSize(pp.nPoints());
-    patchRotationDisp = vector::zero;
-
-
-    // Unweighted interpolation
-    // ~~~~~~~~~~~~~~~~~~~~~~~~
-
-    labelList nPatchDisp(pp.nPoints(), 0.0);
-
-    forAll(pp, faceI)
-    {
-        const face& f = pp.localFaces()[faceI];
-        const point& fc = pp.faceCentres()[faceI];
-
-        forAll(f, fp)
-        {
-            label pointI = f[fp];
-            const point& pt = pp.localPoints()[pointI];
-
-            patchDisp[pointI] += faceDisp[faceI];
-            patchRotationDisp[pointI] += (faceRotation[faceI]^(pt-fc));
-            nPatchDisp[pointI]++;
-        }
-    }
-
-    syncTools::syncPointList
-    (
-        mesh,
-        pp.meshPoints(),
-        patchDisp,
-        plusEqOp<vector>(),
-        vector::zero,
-        mapDistribute::transform()
-    );
-    syncTools::syncPointList
-    (
-        mesh,
-        pp.meshPoints(),
-        patchRotationDisp,
-        plusEqOp<vector>(),
-        vector::zero,
-        mapDistribute::transform()
-    );
-    syncTools::syncPointList
-    (
-        mesh,
-        pp.meshPoints(),
-        nPatchDisp,
-        plusEqOp<label>(),
-        0,
-        mapDistribute::transform()
-    );
-
-    forAll(patchDisp, pointI)
-    {
-        patchDisp[pointI] /= nPatchDisp[pointI];
-        patchRotationDisp[pointI] /= nPatchDisp[pointI];
-    }
-
-
-    if (debug&meshRefinement::OBJINTERSECTIONS)
-    {
-        dumpMove
-        (
-            mesh.time().path()
-          / "patchDisp_" + name(iter) + ".obj",
-            pp.localPoints(),
-            pp.localPoints() + patchDisp
-        );
-        dumpMove
-        (
-            mesh.time().path()
-          / "patchRotationDisp_" + name(iter) + ".obj",
-            pp.localPoints(),
-            pp.localPoints() + patchRotationDisp
-        );
-    }
-}
-
-
 // Gets passed in offset to nearest point on feature edge. Calculates
 // if the point has a different number of faces on either side of the feature
 // and if so attracts the point to that non-dominant plane.
@@ -500,8 +487,6 @@ void Foam::autoSnapDriver::binFeatureFaces
 
     const label pointI,
 
-//    const vectorField& faceSurfaceNormal,
-//    const vectorField& faceDisp,
     const List<List<point> >& pointFaceNormals,
     const List<List<point> >& pointFaceDisp,
     const List<List<point> >& pointFaceCentres,
@@ -547,9 +532,6 @@ void Foam::autoSnapDriver::featureAttractionUsingReconstruction
     const indirectPrimitivePatch& pp,
     const scalarField& snapDist,
 
-    //const vectorField& faceSurfaceNormal,
-    //const vectorField& faceDisp,
-    //const vectorField& faceRotation,
     const List<List<point> >& pointFaceNormals,
     const List<List<point> >& pointFaceDisp,
     const List<List<point> >& pointFaceCentres,
@@ -609,8 +591,6 @@ void Foam::autoSnapDriver::featureAttractionUsingReconstruction
             snapDist,
             pointI,
 
-            //faceSurfaceNormal,
-            //faceDisp,
             pointFaceNormals,
             pointFaceDisp,
             pointFaceCentres,
@@ -741,9 +721,6 @@ void Foam::autoSnapDriver::determineAllFeatures
     const indirectPrimitivePatch& pp,
     const scalarField& snapDist,
 
-    //const vectorField& faceSurfaceNormal,
-    //const vectorField& faceDisp,
-    //const vectorField& faceRotation,
     const List<List<point> >& pointFaceNormals,
     const List<List<point> >& pointFaceDisp,
     const List<List<point> >& pointFaceCentres,
@@ -1778,33 +1755,8 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
         faceRotation
     );
 
-    //
-    // Interpolate face-wise data to points
-    //
-
-
-    // Displacement due to snapping of face centres
-    pointField patchDisp(localPoints.size(), vector::zero);
-    // Displacement due to rotation
-    pointField patchRotationDisp(localPoints.size(), vector::zero);
-
-    interpolateFaceToPoint
-    (
-        iter,
-        pp,
-        faceSurfaceNormal,
-        faceDisp,
-        faceRotation,
-
-        patchDisp,
-        patchRotationDisp
-    );
-
-    // Override patchDisp with exact nearest
-    {
-        pointField pointSurfaceNormal;
-        calcNearest(pp.localPoints(), patchDisp, pointSurfaceNormal);
-    }
+    // Start off with nearest point on surface.
+    vectorField patchDisp = nearestDisp;
 
 
     // Collect (possibly remote) face-wise data on coupled points.
@@ -1818,7 +1770,6 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
 
     List<List<point> > pointFaceNormals(pp.nPoints());
     List<List<point> > pointFaceDisp(pp.nPoints());
-    List<List<point> > pointFaceRotation(pp.nPoints());
     List<List<point> > pointFaceCentres(pp.nPoints());
 
     // Fill local data
@@ -1829,15 +1780,12 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
         pNormals.setSize(pFaces.size());
         List<point>& pDisp = pointFaceDisp[pointI];
         pDisp.setSize(pFaces.size());
-        List<point>& pRot = pointFaceRotation[pointI];
-        pRot.setSize(pFaces.size());
         List<point>& pFc = pointFaceCentres[pointI];
         pFc.setSize(pFaces.size());
         forAll(pFaces, i)
         {
             pNormals[i] = faceSurfaceNormal[pFaces[i]];
             pDisp[i] = faceDisp[pFaces[i]];
-            pRot[i] = faceRotation[pFaces[i]];
             pFc[i] = pp.faceCentres()[pFaces[i]];
         }
     }
@@ -1861,15 +1809,6 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
         listTransform()
     );
     syncTools::syncPointList
-    (
-        mesh,
-        pp.meshPoints(),
-        pointFaceRotation,
-        listPlusEqOp(),
-        List<point>(),
-        listTransform()
-    );
-    syncTools::syncPointList
     (
         mesh,
         pp.meshPoints(),
@@ -1881,60 +1820,6 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
 
 
 
-//XXXXXXX
-//    // Indices of pp points that are coupled
-//    labelList coupledPointLabels(pp.nPoints());
-//    // Corresponding coupled point index
-//    labelList coupledPointAddr(pp.nPoints());
-//    {
-//        label coupledI = 0;
-//        forAll(pp.meshPoints(), pointI)
-//        {
-//            label meshPointI = pp.meshPoints()[pointI];
-//            Map<label>::const_iterator fnd = coupledPatchMP.find(meshPointI);
-//            if (fnd != coupledPatchMP.end())
-//            {
-//                coupledPointLabels[coupledI] = pointI;
-//                coupledPointAddr[coupledI] = fnd();
-//                coupledI++;
-//            }
-//        }
-//        coupledPointLabels.setSize(coupledI);
-//        coupledPointAddr.setSize(coupledI);
-//    }
-//
-//    List<List<point> > pointFaceNormals(map.constructSize());
-//    List<List<point> > pointFaceDisp(map.constructSize());
-//    List<List<point> > pointFaceRotation(map.constructSize());
-//    List<List<point> > pointFaceCentres(map.constructSize());
-//
-//    // Collect local data.
-//    forAll(coupledPointLabels, coupledI)
-//    {
-//        label pointI = coupledPointLabels[coupledI];
-//        label coupledPointI = coupledPointAddr[coupledI];
-//        const labelList& pFaces = pp.pointFaces()[pointI];
-//        List<point>& pNormals = pointFaceNormals[coupledPointI];
-//        pNormals.setSize(pFaces.size());
-//        List<point>& pDisp = pointFaceDisp[coupledPointI];
-//        pDisp.setSize(pFaces.size());
-//        List<point>& pRot = pointFaceRotation[coupledPointI];
-//        pRot.setSize(pFaces.size());
-//        List<point>& pFc = pointFaceCentres[coupledPointI];
-//        pFc.setSize(pFaces.size());
-//        forAll(pFaces, i)
-//        {
-//            pNormals[i] = faceSurfaceNormal[pFaces[i]];
-//            pDisp[i] = pointFaceDisp[pFaces[i]];
-//            pRot[i] = pointFaceRotation[pFaces[i]];
-//            pFc[i] = pp.faceCentres()[pFaces[i]];
-//        }
-//    }
-//
-//    // Pull remote data into local slots
-//XXXXXXX
-
-
     // Nearest feature
     vectorField patchAttraction(localPoints.size(), vector::zero);
     // Constraints at feature
@@ -1973,9 +1858,6 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
         << "     linear   : max:" << gMax(patchDisp)
         << " avg:" << gAverage(patchDisp)
         << endl
-        << "     rotation : max:" << gMax(patchRotationDisp)
-        << " avg:" << gAverage(patchRotationDisp)
-        << endl
         << "     feature  : max:" << gMax(patchAttraction)
         << " avg:" << gAverage(patchAttraction)
         << endl;
@@ -1985,20 +1867,13 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
     // - patchDisp          : point movement to go to nearest point on surface
     //                       (either direct or through interpolation of
     //                        face nearest)
-    // - patchRotationDisp  : point movement to align face normal with surface.
-    //
     // - patchAttraction    : direct attraction to features
     // - patchConstraints   : type of features
 
-    // Use any combination of patchDisp, patchRotationDisp and direct feature
+    // Use any combination of patchDisp and direct feature
     // attraction.
 
 
-    // Bit of patchRotation
-    //patchDisp = 0.3*patchRotationDisp+0.7*patchDisp;
-    // Disable patchRotationDisp
-    patchRotationDisp = vector::zero;
-
     // Mix with direct feature attraction
     forAll(patchConstraints, pointI)
     {
@@ -2017,13 +1892,6 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
     //    pp.localPoints(),
     //    pp.localPoints() + patchDisp
     //);
-    //dumpMove
-    //(
-    //    mesh.time().path()
-    //  / "rotationPatchDisp_" + name(iter) + ".obj",
-    //    pp.localPoints(),
-    //    pp.localPoints() + patchRotationDisp
-    //);
 
 
 
@@ -2130,15 +1998,6 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
         vector(GREAT, GREAT, GREAT)     // null value (note: cannot use VGREAT)
     );
 
-
-    // Check for displacement being outwards.
-    outwardsDisplacement(pp, patchDisp);
-
-    // Set initial distribution of displacement field (on patches) from
-    // patchDisp and make displacement consistent with b.c. on displacement
-    // pointVectorField.
-    meshMover.setDisplacement(patchDisp);
-
     return patchDisp;
 }
 

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementMerge.C

@@ -277,7 +277,11 @@ Foam::label Foam::meshRefinement::mergePatchFacesUndo
         }
     }
 
-    // Get all sets of faces that can be merged
+    // Get all sets of faces that can be merged. Since only faces on the same
+    // patch get merged there is no risk of e.g. patchID faces getting merged
+    // with original patches (or even processor patches). There is a risk
+    // though of original patchfaces getting merged if they make a small
+    // angle. Would be pretty weird starting mesh though.
     labelListList allFaceSets
     (
         faceCombiner.getMergeSets
@@ -349,10 +353,25 @@ Foam::label Foam::meshRefinement::mergePatchFacesUndo
 
         faceCombiner.updateMesh(map);
 
-        updateMesh(map, labelList(0));
+        // Get the kept faces that need to be recalculated.
+        // Merging two boundary faces might shift the cell centre
+        // (unless the faces are absolutely planar)
+        labelHashSet retestFaces(2*allFaceSets.size());
+
+        forAll(allFaceSets, setI)
+        {
+            label oldMasterI = allFaceSets[setI][0];
+            retestFaces.insert(map().reverseFaceMap()[oldMasterI]);
+        }
+        updateMesh(map, growFaceCellFace(retestFaces));
 
         if (debug)
         {
+            // Check sync
+            Pout<< "Checking sync after initial merging " << nFaceSets
+                << " faces." << endl;
+            checkData();
+
             Pout<< "Writing initial merged-faces mesh to time "
                 << timeName() << nl << endl;
             write();
@@ -529,11 +548,22 @@ Foam::label Foam::meshRefinement::mergePatchFacesUndo
             inplaceMapKey(map().reverseFaceMap(), restoredFaces);
             inplaceMapKey(map().reverseCellMap(), restoredCells);
 
+
+            // Get the kept faces that need to be recalculated.
+            // Merging two boundary faces might shift the cell centre
+            // (unless the faces are absolutely planar)
+            labelHashSet retestFaces(2*restoredFaces.size());
+
+            forAllConstIter(Map<label>, restoredFaces, iter)
+            {
+                retestFaces.insert(iter.key());
+            }
+
             // Experimental:restore all points/face/cells in maps
             updateMesh
             (
                 map,
-                labelList(0),       // changedFaces
+                growFaceCellFace(retestFaces),
                 restoredPoints,
                 restoredFaces,
                 restoredCells
@@ -541,6 +571,11 @@ Foam::label Foam::meshRefinement::mergePatchFacesUndo
 
             if (debug)
             {
+                // Check sync
+                Pout<< "Checking sync after restoring " << retestFaces.size()
+                    << " faces." << endl;
+                checkData();
+
                 Pout<< "Writing merged-faces mesh to time "
                     << timeName() << nl << endl;
                 write();
@@ -591,7 +626,26 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::doRemovePoints
     mesh_.setInstance(timeName());
 
     pointRemover.updateMesh(map);
-    updateMesh(map, labelList(0));
+
+
+    // Retest all affected faces and all the cells using them
+    labelHashSet retestFaces(pointRemover.savedFaceLabels().size());
+    forAll(pointRemover.savedFaceLabels(), i)
+    {
+        label faceI = pointRemover.savedFaceLabels()[i];
+        if (faceI >= 0)
+        {
+            retestFaces.insert(faceI);
+        }
+    }
+    updateMesh(map, growFaceCellFace(retestFaces));
+
+    if (debug)
+    {
+        // Check sync
+        Pout<< "Checking sync after removing points." << endl;
+        checkData();
+    }
 
     return map;
 }
@@ -645,7 +699,25 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::doRestorePoints
     mesh_.setInstance(timeName());
 
     pointRemover.updateMesh(map);
-    updateMesh(map, labelList(0));
+
+    labelHashSet retestFaces(2*facesToRestore.size());
+    forAll(facesToRestore, i)
+    {
+        label faceI = map().reverseFaceMap()[facesToRestore[i]];
+        if (faceI >= 0)
+        {
+            retestFaces.insert(faceI);
+        }
+    }
+    updateMesh(map, growFaceCellFace(retestFaces));
+
+    if (debug)
+    {
+        // Check sync
+        Pout<< "Checking sync after restoring points on "
+            << facesToRestore.size() << " faces." << endl;
+        checkData();
+    }
 
     return map;
 }

src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C

@@ -86,7 +86,16 @@ Foam::refinementFeatures::refinementFeatures
         const edgeList& edges = eMesh.edges();
 
         // Calculate bb of all points
-        const treeBoundBox bb(points);
+        treeBoundBox bb(points);
+
+        // 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);
 
         edgeTrees_.set
         (
@@ -120,7 +129,8 @@ Foam::refinementFeatures::refinementFeatures
         }
 
         Info<< "Detected " << featurePoints.size()
-            << " featurePoints out of " << points.size() << endl;
+            << " featurePoints out of " << points.size()
+            << " on feature " << eMesh.name() << endl;
 
         pointTrees_.set
         (
@@ -155,26 +165,30 @@ void Foam::refinementFeatures::findNearestEdge
     forAll(edgeTrees_, featI)
     {
         const indexedOctree<treeDataEdge>& tree = edgeTrees_[featI];
-        forAll(samples, sampleI)
+
+        if (tree.shapes().size() > 0)
         {
-            const point& sample = samples[sampleI];
-
-            scalar distSqr;
-            if (nearInfo[sampleI].hit())
+            forAll(samples, sampleI)
             {
-                distSqr = magSqr(nearInfo[sampleI].hitPoint()-sample);
-            }
-            else
-            {
-                distSqr = nearestDistSqr[sampleI];
-            }
+                const point& sample = samples[sampleI];
 
-            pointIndexHit info = tree.findNearest(sample, distSqr);
+                scalar distSqr;
+                if (nearInfo[sampleI].hit())
+                {
+                    distSqr = magSqr(nearInfo[sampleI].hitPoint()-sample);
+                }
+                else
+                {
+                    distSqr = nearestDistSqr[sampleI];
+                }
 
-            if (info.hit())
-            {
-                nearInfo[sampleI] = info;
-                nearFeature[sampleI] = featI;
+                pointIndexHit info = tree.findNearest(sample, distSqr);
+
+                if (info.hit())
+                {
+                    nearInfo[sampleI] = info;
+                    nearFeature[sampleI] = featI;
+                }
             }
         }
     }
@@ -197,33 +211,37 @@ void Foam::refinementFeatures::findNearestPoint
     forAll(pointTrees_, featI)
     {
         const indexedOctree<treeDataPoint>& tree = pointTrees_[featI];
-        forAll(samples, sampleI)
+
+        if (tree.shapes().pointLabels().size() > 0)
         {
-            const point& sample = samples[sampleI];
-
-            scalar distSqr;
-            if (nearFeature[sampleI] != -1)
+            forAll(samples, sampleI)
             {
-                label nearFeatI = nearFeature[sampleI];
-                const indexedOctree<treeDataPoint>& nearTree =
-                    pointTrees_[nearFeatI];
-                label featPointI =
-                    nearTree.shapes().pointLabels()[nearIndex[sampleI]];
-                const point& featPt =
-                    operator[](nearFeatI).points()[featPointI];
-                distSqr = magSqr(featPt-sample);
-            }
-            else
-            {
-                distSqr = nearestDistSqr[sampleI];
-            }
+                const point& sample = samples[sampleI];
 
-            pointIndexHit info = tree.findNearest(sample, distSqr);
+                scalar distSqr;
+                if (nearFeature[sampleI] != -1)
+                {
+                    label nearFeatI = nearFeature[sampleI];
+                    const indexedOctree<treeDataPoint>& nearTree =
+                        pointTrees_[nearFeatI];
+                    label featPointI =
+                        nearTree.shapes().pointLabels()[nearIndex[sampleI]];
+                    const point& featPt =
+                        operator[](nearFeatI).points()[featPointI];
+                    distSqr = magSqr(featPt-sample);
+                }
+                else
+                {
+                    distSqr = nearestDistSqr[sampleI];
+                }
 
-            if (info.hit())
-            {
-                nearFeature[sampleI] = featI;
-                nearIndex[sampleI] = info.index();
+                pointIndexHit info = tree.findNearest(sample, distSqr);
+
+                if (info.hit())
+                {
+                    nearFeature[sampleI] = featI;
+                    nearIndex[sampleI] = info.index();
+                }
             }
         }
     }

src/meshTools/cellDist/patchWave/patchWave.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -27,6 +27,7 @@ License
 #include "polyMesh.H"
 #include "wallPoint.H"
 #include "globalMeshData.H"
+#include "SubField.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 

src/meshTools/cellQuality/cellQuality.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,6 +26,7 @@ License
 
 #include "cellQuality.H"
 #include "unitConversion.H"
+#include "SubField.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 

src/meshTools/directMapped/directMappedPolyPatch/directMappedPatchBase.C

@@ -36,6 +36,7 @@ License
 #include "polyPatch.H"
 #include "Time.H"
 #include "mapDistribute.H"
+#include "SubField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/meshTools/indexedOctree/treeDataPoint.C

@@ -37,7 +37,8 @@ defineTypeNameAndDebug(Foam::treeDataPoint, 0);
 
 Foam::treeDataPoint::treeDataPoint(const pointField& points)
 :
-    points_(points)
+    points_(points),
+    useSubset_(false)
 {}
 
 
@@ -48,7 +49,8 @@ Foam::treeDataPoint::treeDataPoint
 )
 :
     points_(points),
-    pointLabels_(pointLabels)
+    pointLabels_(pointLabels),
+    useSubset_(true)
 {}
 
 
@@ -56,7 +58,7 @@ Foam::treeDataPoint::treeDataPoint
 
 Foam::pointField Foam::treeDataPoint::shapePoints() const
 {
-    if (pointLabels_.size())
+    if (useSubset_)
     {
         return pointField(points_, pointLabels_);
     }
@@ -86,7 +88,7 @@ bool Foam::treeDataPoint::overlaps
     const treeBoundBox& cubeBb
 ) const
 {
-    label pointI = (pointLabels_.size() ? pointLabels_[index] : index);
+    label pointI = (useSubset_ ? pointLabels_[index] : index);
     return cubeBb.contains(points_[pointI]);
 }
 
@@ -106,7 +108,7 @@ void Foam::treeDataPoint::findNearest
     forAll(indices, i)
     {
         const label index = indices[i];
-        label pointI = (pointLabels_.size() ? pointLabels_[index] : index);
+        label pointI = (useSubset_ ? pointLabels_[index] : index);
 
         const point& pt = points_[pointI];
 
@@ -141,7 +143,7 @@ void Foam::treeDataPoint::findNearest
     forAll(indices, i)
     {
         const label index = indices[i];
-        label pointI = (pointLabels_.size() ? pointLabels_[index] : index);
+        label pointI = (useSubset_ ? pointLabels_[index] : index);
 
         const point& shapePt = points_[pointI];
 

src/meshTools/indexedOctree/treeDataPoint.H

@@ -65,6 +65,8 @@ class treeDataPoint
         //- Subset of points to work on (or empty)
         const labelList pointLabels_;
 
+        const bool useSubset_;
+
 public:
 
     // Declare name of the class and its debug switch
@@ -88,7 +90,7 @@ public:
             {
                 return
                 (
-                    pointLabels_.size()
+                    useSubset_
                   ? pointLabels_.size()
                   : points_.size()
                 );

src/meshTools/twoDPointCorrector/twoDPointCorrector.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -27,6 +27,7 @@ License
 #include "polyMesh.H"
 #include "wedgePolyPatch.H"
 #include "emptyPolyPatch.H"
+#include "SubField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/parallel/decompose/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.C

@@ -27,6 +27,7 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "SortableList.H"
 #include "globalIndex.H"
+#include "SubField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C

@@ -117,6 +117,7 @@ License
 #include "Time.H"
 #include "OFstream.H"
 #include "globalIndex.H"
+#include "SubField.H"
 
 extern "C"
 {

src/parallel/decompose/scotchDecomp/Make/options

@@ -1,4 +1,11 @@
+# Note including of mplib compilation rules. This
+# is purely to avoid scotch.h including mpicxx.h
+# which causes problems.
+sinclude $(GENERAL_RULES)/mplib$(WM_MPLIB)
+sinclude $(RULES)/mplib$(WM_MPLIB)
+
 EXE_INC = \
+    $(PFLAGS) $(PINC) \
     -I$(SCOTCH_ARCH_PATH)/include \
     -I/usr/include/scotch \
     -I../decompositionMethods/lnInclude

src/parallel/decompose/scotchDecomp/scotchDecomp.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2009-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2009-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -65,7 +65,7 @@ License
     The current default mapping strategy in Scotch can be seen by using the
     "-vs" option of program gmap. It is, to date:
 
-    b
+    r
     {
         job=t,
         map=t,
@@ -76,11 +76,17 @@ License
             {
                 asc=b
                 {
-                    bnd=d{pass=40,dif=1,rem=1}f{move=80,pass=-1,bal=0.005},
-                    org=f{move=80,pass=-1,bal=0.005},
+                    bnd=
+                    (
+                        d{pass=40,dif=1,rem=1}
+                     |
+                    )
+                    f{move=80,pass=-1,bal=0.002491},
+                    org=f{move=80,pass=-1,bal=0.002491},
                     width=3
                 },
-                low=h{pass=10}f{move=80,pass=-1,bal=0.0005},
+                low=h{pass=10}
+                f{move=80,pass=-1,bal=0.002491},
                 type=h,
                 vert=80,
                 rat=0.8
@@ -89,11 +95,17 @@ License
             {
                 asc=b
                 {
-                    bnd=d{pass=40,dif=1,rem=1}f{move=80,pass=-1,bal=0.005},
-                    org=f{move=80,pass=-1,bal=0.005},
+                    bnd=
+                    (
+                        d{pass=40,dif=1,rem=1}
+                      |
+                    )
+                    f{move=80,pass=-1,bal=0.002491},
+                    org=f{move=80,pass=-1,bal=0.002491},
                     width=3
                 },
-                low=h{pass=10}f{move=80,pass=-1,bal=0.0005},
+                low=h{pass=10}
+                f{move=80,pass=-1,bal=0.002491},
                 type=h,
                 vert=80,
                 rat=0.8
@@ -102,6 +114,9 @@ License
     }
 
 
+    Note: instead of gmap run gpart <nProcs> -vs <grfFile>
+    where <grfFile> can be obtained by running with 'writeGraph=true'
+
 \*---------------------------------------------------------------------------*/
 
 #include "scotchDecomp.H"
@@ -110,6 +125,7 @@ License
 #include "Time.H"
 #include "OFstream.H"
 #include "globalIndex.H"
+#include "SubField.H"
 
 extern "C"
 {

src/postProcessing/functionObjects/jobControl/Make/files

@@ -1,7 +1,4 @@
 abortCalculation/abortCalculation.C
 abortCalculation/abortCalculationFunctionObject.C
 
-residualControl/residualControl.C
-residualControl/residualControlFunctionObject.C
-
 LIB = $(FOAM_LIBBIN)/libjobControl

src/postProcessing/functionObjects/jobControl/residualControl/controlDict

@@ -1,68 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    location    "system";
-    object      controlDict;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-application     icoFoam;
-
-startFrom       startTime;
-
-startTime       0;
-
-stopAt          endTime;
-
-endTime         0.5;
-
-deltaT          0.005;
-
-writeControl    timeStep;
-
-writeInterval   20;
-
-purgeWrite      0;
-
-writeFormat     ascii;
-
-writePrecision  6;
-
-writeCompression off;
-
-timeFormat      general;
-
-timePrecision   6;
-
-runTimeModifiable true;
-
-functions
-{
-    convergenceChecks
-    {
-        type            residualControl;
-        functionObjectLibs ( "libjobControl.so" );
-        outputControl   timeStep;
-        outputInterval  1;
-
-        maxResiduals
-        {
-            p       5e-4;
-            U       1e-3;
-
-            // possibly check turbulence fields
-            "(k|epsilon|omega)" 1e-3;
-        }
-    }
-}
-
-// ************************************************************************* //

src/postProcessing/functionObjects/jobControl/residualControl/residualControl.C

@@ -1,173 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "residualControl.H"
-#include "fvMesh.H"
-#include "Time.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-defineTypeNameAndDebug(Foam::residualControl, 0);
-
-
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
-bool Foam::residualControl::checkCriteria(const bool verbose) const
-{
-    bool achieved = true;
-    const fvMesh& mesh = static_cast<const fvMesh&>(obr_);
-    const dictionary& solverDict = mesh.solverPerformanceDict();
-
-    forAllConstIter(dictionary, solverDict, iter)
-    {
-        const word& variableName = iter().keyword();
-        scalar maxResidual;
-
-        if (maxResiduals_.readIfPresent(variableName, maxResidual))
-        {
-            // use the residual from the first solution
-            const scalar eqnResidual =
-                List<lduMatrix::solverPerformance>
-                (
-                    iter().stream()
-                ).first().initialResidual();
-
-            achieved = achieved && (eqnResidual < maxResidual);
-
-            if (verbose)
-            {
-                Info<< "    " << variableName
-                    << ": requested max residual = " << maxResidual
-                    << ", eqn residual = " << eqnResidual << nl;
-            }
-        }
-    }
-
-    return achieved;
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::residualControl::residualControl
-(
-    const word& name,
-    const objectRegistry& obr,
-    const dictionary& dict,
-    const bool
-)
-:
-    name_(name),
-    obr_(obr),
-    active_(true),
-    maxResiduals_(),
-    criteriaSatisfied_(false)
-{
-    // Only active if a fvMesh is available
-    if (isA<fvMesh>(obr_))
-    {
-        read(dict);
-    }
-    else
-    {
-        active_ = false;
-        WarningIn
-        (
-            "residualControl::residualControl"
-            "("
-                "const word&, "
-                "const objectRegistry&, "
-                "const dictionary&, "
-                "const bool"
-            ")"
-        )   << "No fvMesh available, deactivating."
-            << nl << endl;
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::residualControl::~residualControl()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::residualControl::read(const dictionary& dict)
-{
-    if (active_)
-    {
-        maxResiduals_ = dict.subDict("maxResiduals");
-    }
-}
-
-
-void Foam::residualControl::execute()
-{
-    if (active_)
-    {
-        criteriaSatisfied_ = checkCriteria(false);
-
-        if (criteriaSatisfied_)
-        {
-            Info<< "Convergence criteria satisfied - finalising run" << nl
-                << endl;
-
-            checkCriteria(true);
-
-            Info<< endl;
-
-            const fvMesh& mesh = static_cast<const fvMesh&>(obr_);
-            Time& time = const_cast<Time&>(mesh.time());
-            time.writeAndEnd();
-        }
-    }
-}
-
-
-void Foam::residualControl::end()
-{
-    if (active_)
-    {
-        if (criteriaSatisfied_)
-        {
-            Info<< "Residual control criteria satisfied" << nl;
-        }
-        else
-        {
-            Info<< "Residual control criteria not satisfied" << nl;
-        }
-    }
-}
-
-
-void Foam::residualControl::write()
-{
-    // do nothing
-}
-
-
-// ************************************************************************* //

src/postProcessing/functionObjects/jobControl/residualControl/residualControl.H

@@ -1,152 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::residualControl
-
-Description
-    Function object that allows users to set target convergence criteria, and
-    stop the run if the conditions are satisfied.
-
-SourceFiles
-    residualControl.C
-    IOresidualControl.H
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef residualControl_H
-#define residualControl_H
-
-#include "dictionary.H"
-#include "pointFieldFwd.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// Forward declaration of classes
-class objectRegistry;
-class dictionary;
-class mapPolyMesh;
-
-/*---------------------------------------------------------------------------*\
-                       Class residualControl Declaration
-\*---------------------------------------------------------------------------*/
-
-class residualControl
-{
-protected:
-
-    // Private data
-
-        //- Name of this object
-        word name_;
-
-        //- Reference to object registry
-        const objectRegistry& obr_;
-
-        //- On/off switch - on if obr_ is an fvMesh object
-        bool active_;
-
-        //- Dictionary of variable names vs max residual
-        dictionary maxResiduals_;
-
-        //- Flag to indicate whether convergence criteria have been met
-        bool criteriaSatisfied_;
-
-
-    // Protected Member Functions
-
-        //- Perform residual control checks
-        bool checkCriteria(const bool verbose) const;
-
-        //- Disallow default bitwise copy construct
-        residualControl(const residualControl&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const residualControl&);
-
-
-public:
-
-    //- Runtime type information
-    TypeName("residualControl");
-
-
-    // Constructors
-
-        //- Construct for given objectRegistry and dictionary.
-        //  Allow the possibility to load fields from files
-        residualControl
-        (
-            const word& name,
-            const objectRegistry& obr,
-            const dictionary& dict,
-            const bool loadFromFilesUnused = false
-        );
-
-
-    //- Destructor
-    virtual ~residualControl();
-
-
-    // Member Functions
-
-        //- Return name of the residual criteria check name
-        virtual const word& name() const
-        {
-            return name_;
-        }
-
-        //- Read the system calls
-        virtual void read(const dictionary&);
-
-        //- Check the residual criteria at each time-step
-        virtual void execute();
-
-        //- Report residual criteria check at the final time-loop
-        virtual void end();
-
-        //- Write, not used
-        virtual void write();
-
-        //- Update for changes of mesh
-        virtual void updateMesh(const mapPolyMesh&)
-        {}
-
-        //- Update for changes of mesh
-        virtual void movePoints(const pointField&)
-        {}
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/postProcessing/functionObjects/jobControl/residualControl/residualControlFunctionObject.H

@@ -1,54 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Typedef
-    Foam::residualControlFunctionObject
-
-Description
-    FunctionObject wrapper around residualControl to allow them to be created
-    via the functions entry within controlDict.
-
-SourceFiles
-    residualControlFunctionObject.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef residualControlFunctionObject_H
-#define residualControlFunctionObject_H
-
-#include "residualControl.H"
-#include "OutputFilterFunctionObject.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    typedef OutputFilterFunctionObject<residualControl>
-        residualControlFunctionObject;
-}
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/postProcessing/functionObjects/utilities/codedFunctionObject/codedFunctionObject.C

@@ -54,60 +54,65 @@ void* Foam::codedFunctionObject::loadLibrary
     const fileName& libPath,
     const string& globalFuncName,
     const dictionary& contextDict
-)
+) const
 {
     void* lib = 0;
 
     // avoid compilation by loading an existing library
-    if (!libPath.empty() && dlLibraryTable::open(libPath, false))
+    if (!libPath.empty())
     {
-        lib = dlLibraryTable::findLibrary(libPath);
+        dlLibraryTable& libs = const_cast<Time&>(time_).libs();
 
-        // verify the loaded version and unload if needed
-        if (lib)
+        if (libs.open(libPath, false))
         {
-            // provision for manual execution of code after loading
-            if (dlSymFound(lib, globalFuncName))
-            {
-                loaderFunctionType function =
-                    reinterpret_cast<loaderFunctionType>
-                    (
-                        dlSym(lib, globalFuncName)
-                    );
+            lib = libs.findLibrary(libPath);
 
-                if (function)
+            // verify the loaded version and unload if needed
+            if (lib)
+            {
+                // provision for manual execution of code after loading
+                if (dlSymFound(lib, globalFuncName))
                 {
-                    (*function)(true);    // force load
+                    loaderFunctionType function =
+                        reinterpret_cast<loaderFunctionType>
+                        (
+                            dlSym(lib, globalFuncName)
+                        );
+
+                    if (function)
+                    {
+                        (*function)(true);    // force load
+                    }
+                    else
+                    {
+                        FatalIOErrorIn
+                        (
+                            "codedFunctionObject::updateLibrary()",
+                            contextDict
+                        )   << "Failed looking up symbol " << globalFuncName
+                            << nl << "from " << libPath << exit(FatalIOError);
+                    }
                 }
                 else
                 {
                     FatalIOErrorIn
                     (
-                        "codedFunctionObject::updateLibrary()",
+                        "codedFunctionObject::loadLibrary()",
                         contextDict
                     )   << "Failed looking up symbol " << globalFuncName << nl
                         << "from " << libPath << exit(FatalIOError);
-                }
-            }
-            else
-            {
-                FatalIOErrorIn
-                (
-                    "codedFunctionObject::loadLibrary()",
-                    contextDict
-                )   << "Failed looking up symbol " << globalFuncName << nl
-                    << "from " << libPath << exit(FatalIOError);
 
-                lib = 0;
-                if (!dlLibraryTable::close(libPath, false))
-                {
-                    FatalIOErrorIn
-                    (
-                        "codedFunctionObject::loadLibrary()",
-                        contextDict
-                    )   << "Failed unloading library "
-                        << libPath
-                        << exit(FatalIOError);
+                    lib = 0;
+                    if (!libs.close(libPath, false))
+                    {
+                        FatalIOErrorIn
+                        (
+                            "codedFunctionObject::loadLibrary()",
+                            contextDict
+                        )   << "Failed unloading library "
+                            << libPath
+                            << exit(FatalIOError);
+                    }
                 }
             }
         }
@@ -122,15 +127,19 @@ void Foam::codedFunctionObject::unloadLibrary
     const fileName& libPath,
     const string& globalFuncName,
     const dictionary& contextDict
-)
+) const
 {
     void* lib = 0;
 
-    if (!libPath.empty())
+    if (libPath.empty())
     {
-        lib = dlLibraryTable::findLibrary(libPath);
+        return;
     }
 
+    dlLibraryTable& libs = const_cast<Time&>(time_).libs();
+
+    lib = libs.findLibrary(libPath);
+
     if (!lib)
     {
         return;
@@ -160,7 +169,7 @@ void Foam::codedFunctionObject::unloadLibrary
         }
     }
 
-    if (!dlLibraryTable::close(libPath, false))
+    if (!libs.close(libPath, false))
     {
         FatalIOErrorIn
         (
@@ -274,7 +283,7 @@ void Foam::codedFunctionObject::updateLibrary() const
 
 
     // the correct library was already loaded => we are done
-    if (dlLibraryTable::findLibrary(libPath))
+    if (const_cast<Time&>(time_).libs().findLibrary(libPath))
     {
         return;
     }

src/postProcessing/functionObjects/utilities/codedFunctionObject/codedFunctionObject.H

@@ -88,20 +88,20 @@ protected:
         typedef void (*loaderFunctionType)(bool);
 
         //- Load specified library and execute globalFuncName(true)
-        static void* loadLibrary
+        void* loadLibrary
         (
             const fileName& libPath,
             const string& globalFuncName,
             const dictionary& contextDict
-        );
+        ) const;
 
         //- Execute globalFuncName(false) and unload specified library
-        static void unloadLibrary
+        void unloadLibrary
         (
             const fileName& libPath,
             const string& globalFuncName,
             const dictionary& contextDict
-        );
+        ) const;
 
 
         //- Create library based on the dynamicCodeContext

src/regionModels/regionModel/regionModel1D/regionModel1D.C

@@ -114,7 +114,6 @@ void Foam::regionModels::regionModel1D::initialise()
             boundaryFaceCells_[localPyrolysisFaceI].transfer(cellIDs);
 
             localPyrolysisFaceI++;
-
             nLayers_ = nCells;
         }
     }

src/regionModels/regionModel/regionModel1D/regionModel1D.H

@@ -155,14 +155,14 @@ public:
                 //- Return the global boundary face IDs oppossite coupled patch
                 inline const labelList& boundaryFaceOppositeFace() const;
 
-                //- Return the number of layers in the region
-                inline label nLayers() const;
-
 
             // Geometry
 
                 //- Return the face area magnitudes / [m2]
                 inline const surfaceScalarField& nMagSf() const;
+
+                //- Return the number of layers in the region
+                inline label nLayers() const;
 };
 
 

src/regionModels/regionModel/regionModel1D/regionModel1DI.H

@@ -49,12 +49,6 @@ Foam::regionModels::regionModel1D::boundaryFaceOppositeFace() const
 }
 
 
-inline Foam::label Foam::regionModels::regionModel1D::nLayers() const
-{
-    return nLayers_;
-}
-
-
 inline const Foam::surfaceScalarField&
 Foam::regionModels::regionModel1D::nMagSf() const
 {
@@ -71,4 +65,10 @@ Foam::regionModels::regionModel1D::nMagSf() const
 }
 
 
+inline Foam::label Foam::regionModels::regionModel1D::nLayers() const
+{
+    return nLayers_;
+}
+
+
 // ************************************************************************* //

src/regionModels/regionModel/singleLayerRegion/singleLayerRegion.C

@@ -111,14 +111,12 @@ void Foam::regionModels::singleLayerRegion::initialise()
 
     if (nBoundaryFaces != regionMesh().nCells())
     {
-        /*
         FatalErrorIn("singleLayerRegion::initialise()")
             << "Number of primary region coupled boundary faces not equal to "
             << "the number of cells in the local region" << nl << nl
             << "Number of cells = " << regionMesh().nCells() << nl
             << "Boundary faces  = " << nBoundaryFaces << nl
             << abort(FatalError);
-        */
     }
 
     scalarField passiveMagSf(magSf.size(), 0.0);
@@ -178,12 +176,11 @@ Foam::regionModels::singleLayerRegion::singleLayerRegion
     bool readFields
 )
 :
-    regionModel(mesh, regionType, modelName, readFields),
+    regionModel(mesh, regionType, modelName, false),
     nHatPtr_(NULL),
     magSfPtr_(NULL),
     passivePatchIDs_()
 {
-    Info << "singleLayerRegion" << endl;
     if (active_)
     {
         constructMeshObjects();

src/regionModels/surfaceFilmModels/Make/files

@@ -12,9 +12,10 @@ submodels/subModelBase.C
 KINEMATICMODELS=submodels/kinematic
 $(KINEMATICMODELS)/injectionModel/injectionModel/injectionModel.C
 $(KINEMATICMODELS)/injectionModel/injectionModel/injectionModelNew.C
-$(KINEMATICMODELS)/injectionModel/noInjection/noInjection.C
-$(KINEMATICMODELS)/injectionModel/cloudInjection/cloudInjection.C
+$(KINEMATICMODELS)/injectionModel/injectionModelList/injectionModelList.C
+$(KINEMATICMODELS)/injectionModel/drippingInjection/drippingInjection.C
 $(KINEMATICMODELS)/injectionModel/removeInjection/removeInjection.C
+$(KINEMATICMODELS)/injectionModel/curvatureSeparation/curvatureSeparation.C
 
 THERMOMODELS=submodels/thermo
 $(THERMOMODELS)/phaseChangeModel/phaseChangeModel/phaseChangeModel.C

src/regionModels/surfaceFilmModels/derivedFvPatchFields/filmHeightInletVelocity/filmHeightInletVelocityFvPatchVectorField.C

@@ -123,10 +123,10 @@ void Foam::filmHeightInletVelocityFvPatchVectorField::updateCoeffs()
     const fvPatchField<scalar>& deltafp =
         patch().lookupPatchField<volScalarField, scalar>(deltafName_);
 
-    const vectorField& n = patch().nf();
+    vectorField n(patch().nf());
     const scalarField& magSf = patch().magSf();
 
-    operator==(deltafp*n*phip/(rhop*magSf*sqr(deltafp) + ROOTVSMALL));
+    operator==(n*phip/(rhop*magSf*deltafp + ROOTVSMALL));
 
     fixedValueFvPatchVectorField::updateCoeffs();
 }

src/regionModels/surfaceFilmModels/derivedFvPatchFields/wallFunctions/alphatFilmWallFunction/alphatFilmWallFunctionFvPatchScalarField.C

@@ -156,9 +156,9 @@ void alphatFilmWallFunctionFvPatchScalarField::updateCoeffs()
 
     const mapDistribute& distMap = filmModel.mappedPatches()[filmPatchI].map();
 
-    scalarField mDotFilm =
-        filmModel.massPhaseChangeForPrimary().boundaryField()[filmPatchI];
-    distMap.distribute(mDotFilm);
+    tmp<volScalarField> mDotFilm(filmModel.primaryMassTrans());
+    scalarField mDotFilmp = mDotFilm().boundaryField()[filmPatchI];
+    distMap.distribute(mDotFilmp);
 
     // Retrieve RAS turbulence model
     const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
@@ -185,7 +185,7 @@ void alphatFilmWallFunctionFvPatchScalarField::updateCoeffs()
         scalar Pr = muw[faceI]/alphaw[faceI];
 
         scalar factor = 0.0;
-        scalar mStar = mDotFilm[faceI]/(y[faceI]*uTau);
+        scalar mStar = mDotFilmp[faceI]/(y[faceI]*uTau);
         if (yPlus > yPlusCrit_)
         {
             scalar expTerm = exp(min(50.0, yPlusCrit_*mStar*Pr));
@@ -209,6 +209,7 @@ void alphatFilmWallFunctionFvPatchScalarField::updateCoeffs()
 }
 
 
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 void alphatFilmWallFunctionFvPatchScalarField::write(Ostream& os) const

src/regionModels/surfaceFilmModels/derivedFvPatchFields/wallFunctions/mutkFilmWallFunction/mutkFilmWallFunctionFvPatchScalarField.C

@@ -72,9 +72,9 @@ tmp<scalarField> mutkFilmWallFunctionFvPatchScalarField::calcUTau
 
     const mapDistribute& distMap = filmModel.mappedPatches()[filmPatchI].map();
 
-    scalarField mDotFilm =
-        filmModel.massPhaseChangeForPrimary().boundaryField()[filmPatchI];
-    distMap.distribute(mDotFilm);
+    tmp<volScalarField> mDotFilm(filmModel.primaryMassTrans());
+    scalarField mDotFilmp = mDotFilm().boundaryField()[filmPatchI];
+    distMap.distribute(mDotFilmp);
 
 
     // Retrieve RAS turbulence model
@@ -95,7 +95,7 @@ tmp<scalarField> mutkFilmWallFunctionFvPatchScalarField::calcUTau
 
         scalar yPlus = y[faceI]*ut/(muw[faceI]/rhow[faceI]);
 
-        scalar mStar = mDotFilm[faceI]/(y[faceI]*ut);
+        scalar mStar = mDotFilmp[faceI]/(y[faceI]*ut);
 
         scalar factor = 0.0;
         if (yPlus > yPlusCrit_)

src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayerI.H

@@ -24,7 +24,8 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "kinematicSingleLayer.H"
+#include "surfaceInterpolate.H"
+#include "fvcSurfaceIntegrate.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -163,9 +164,24 @@ inline const volScalarField& kinematicSingleLayer::muPrimary() const
 }
 
 
-inline injectionModel& kinematicSingleLayer::injection()
+inline injectionModelList& kinematicSingleLayer::injection()
 {
-    return injection_();
+    return injection_;
+}
+
+
+inline tmp<volScalarField> kinematicSingleLayer::mass() const
+{
+    return rho_*delta_*magSf();
+}
+
+
+inline tmp<volScalarField> kinematicSingleLayer::netMass() const
+{
+    dimensionedScalar d0("SMALL", dimLength, ROOTVSMALL);
+    return
+        fvc::surfaceSum(phi_/(fvc::interpolate(delta_) + d0))*time().deltaT()
+      + rho_*delta_*magSf();
 }
 
 

src/regionModels/surfaceFilmModels/noFilm/noFilm.C

@@ -103,6 +103,15 @@ const volScalarField& noFilm::delta() const
 }
 
 
+const volScalarField& noFilm::sigma() const
+{
+    FatalErrorIn("const volScalarField& noFilm::sigma() const")
+        << "sigma field not available for " << type() << abort(FatalError);
+
+    return volScalarField::null();
+}
+
+
 const volVectorField& noFilm::U() const
 {
     FatalErrorIn("const volVectorField& noFilm::U() const")
@@ -184,32 +193,42 @@ const volScalarField& noFilm::kappa() const
 }
 
 
-const volScalarField& noFilm::massForPrimary() const
+tmp<volScalarField> noFilm::primaryMassTrans() const
 {
-    FatalErrorIn("const volScalarField& noFilm::massForPrimary() const")
-        << "massForPrimary field not available for " << type()
-        << abort(FatalError);
-
-    return volScalarField::null();
-}
-
-
-const volScalarField& noFilm::diametersForPrimary() const
-{
-    FatalErrorIn("const volScalarField& noFilm::diametersForPrimary() const")
-        << "diametersForPrimary field not available for " << type()
-        << abort(FatalError);
-
-    return volScalarField::null();
-}
-
-
-const volScalarField& noFilm::massPhaseChangeForPrimary() const
-{
-    FatalErrorIn
+    return tmp<volScalarField>
     (
-        "const volScalarField& noFilm::massPhaseChangeForPrimary() const"
-    )   << "massPhaseChange field not available for " << type()
+        new volScalarField
+        (
+            IOobject
+            (
+                "noFilm::primaryMassTrans",
+                time().timeName(),
+                primaryMesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            primaryMesh(),
+            dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
+        )
+    );
+}
+
+
+const volScalarField& noFilm::cloudMassTrans() const
+{
+    FatalErrorIn("const volScalarField& noFilm::cloudMassTrans() const")
+        << "cloudMassTrans field not available for " << type()
+        << abort(FatalError);
+
+    return volScalarField::null();
+}
+
+
+const volScalarField& noFilm::cloudDiameterTrans() const
+{
+    FatalErrorIn("const volScalarField& noFilm::cloudDiameterTrans() const")
+        << "cloudDiameterTrans field not available for " << type()
         << abort(FatalError);
 
     return volScalarField::null();
@@ -238,7 +257,7 @@ tmp<DimensionedField<scalar, volMesh> > noFilm::Srho() const
 }
 
 
-tmp<DimensionedField<scalar, volMesh> > noFilm::Srho(const label) const
+tmp<DimensionedField<scalar, volMesh> > noFilm::Srho(const label i) const
 {
     return tmp<DimensionedField<scalar, volMesh> >
     (
@@ -246,7 +265,7 @@ tmp<DimensionedField<scalar, volMesh> > noFilm::Srho(const label) const
         (
             IOobject
             (
-                "noFilm::Srho(i)",
+                "noFilm::Srho(" + Foam::name(i) + ")",
                 time().timeName(),
                 primaryMesh(),
                 IOobject::NO_READ,

src/regionModels/surfaceFilmModels/noFilm/noFilm.H

@@ -111,49 +111,52 @@ public:
             );
 
 
-       // Fields
+        // Fields
 
-           //- Return the film thickness [m]
-           virtual const volScalarField& delta() const;
+            //- Return the film thickness [m]
+            virtual const volScalarField& delta() const;
 
-           //- Return the film velocity [m/s]
-           virtual const volVectorField& U() const;
+            //- Return const access to the surface tension / [m/s2]
+            inline const volScalarField& sigma() const;
 
-           //- Return the film density [kg/m3]
-           virtual const volScalarField& rho() const;
+            //- Return the film velocity [m/s]
+            virtual const volVectorField& U() const;
 
-           //- Return the film surface velocity [m/s]
-           virtual const volVectorField& Us() const;
+            //- Return the film density [kg/m3]
+            virtual const volScalarField& rho() const;
 
-           //- Return the film wall velocity [m/s]
-           virtual const volVectorField& Uw() const;
+            //- Return the film surface velocity [m/s]
+            virtual const volVectorField& Us() const;
 
-           //- Return the film mean temperature [K]
-           virtual const volScalarField& T() const;
+            //- Return the film wall velocity [m/s]
+            virtual const volVectorField& Uw() const;
 
-           //- Return the film surface temperature [K]
-           virtual const volScalarField& Ts() const;
+            //- Return the film mean temperature [K]
+            virtual const volScalarField& T() const;
 
-           //- Return the film wall temperature [K]
-           virtual const volScalarField& Tw() const;
+            //- Return the film surface temperature [K]
+            virtual const volScalarField& Ts() const;
 
-           //- Return the film specific heat capacity [J/kg/K]
-           virtual const volScalarField& Cp() const;
+            //- Return the film wall temperature [K]
+            virtual const volScalarField& Tw() const;
 
-           //- Return the film thermal conductivity [W/m/K]
-           virtual const volScalarField& kappa() const;
+            //- Return the film specific heat capacity [J/kg/K]
+            virtual const volScalarField& Cp() const;
+
+            //- Return the film thermal conductivity [W/m/K]
+            virtual const volScalarField& kappa() const;
 
 
-           // Transfer fields - to the primary region
+            // Transfer fields - to the primary region
 
-               //- Return the film mass available for transfer
-               virtual const volScalarField& massForPrimary() const;
+                //- Return mass transfer source - Eulerian phase only
+                virtual tmp<volScalarField> primaryMassTrans() const;
 
-               //- Return the parcel diameters originating from film
-               virtual const volScalarField& diametersForPrimary() const;
+                //- Return the film mass available for transfer
+                virtual const volScalarField& cloudMassTrans() const;
 
-               //- Return the film mass evolved via phase change
-               virtual const volScalarField& massPhaseChangeForPrimary() const;
+                //- Return the parcel diameters originating from film
+                virtual const volScalarField& cloudDiameterTrans() const;
 
 
         // Source fields

src/regionModels/surfaceFilmModels/submodels/kinematic/injectionModel/curvatureSeparation/curvatureSeparation.C

@@ -0,0 +1,355 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "curvatureSeparation.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvMesh.H"
+#include "Time.H"
+#include "volFields.H"
+#include "kinematicSingleLayer.H"
+#include "surfaceInterpolate.H"
+#include "fvcDiv.H"
+#include "fvcGrad.H"
+#include "stringListOps.H"
+#include "cyclicPolyPatch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(curvatureSeparation, 0);
+addToRunTimeSelectionTable
+(
+    injectionModel,
+    curvatureSeparation,
+    dictionary
+);
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+tmp<volScalarField> curvatureSeparation::calcInvR1
+(
+    const volVectorField& U
+) const
+{
+    // method 1
+/*
+    tmp<volScalarField> tinvR1
+    (
+        new volScalarField("invR1", fvc::div(owner().nHat()))
+    );
+*/
+
+    // method 2
+    dimensionedScalar smallU("smallU", dimVelocity, ROOTVSMALL);
+    volVectorField UHat(U/(mag(U) + smallU));
+    tmp<volScalarField> tinvR1
+    (
+        new volScalarField("invR1", UHat & (UHat & gradNHat_))
+    );
+
+
+    scalarField& invR1 = tinvR1().internalField();
+
+    // apply defined patch radii
+    const scalar rMin = 1e-6;
+    const fvMesh& mesh = owner().regionMesh();
+    const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+    forAll(definedPatchRadii_, i)
+    {
+        label patchI = definedPatchRadii_[i].first();
+        scalar definedInvR1 = 1.0/max(rMin, definedPatchRadii_[i].second());
+        UIndirectList<scalar>(invR1, pbm[patchI].faceCells()) = definedInvR1;
+    }
+
+    // filter out large radii
+    const scalar rMax = 1e6;
+    forAll(invR1, i)
+    {
+        if (mag(invR1[i]) < 1/rMax)
+        {
+            invR1[i] = -1.0;
+        }
+    }
+
+    if (debug && mesh.time().outputTime())
+    {
+        tinvR1().write();
+    }
+
+    return tinvR1;
+}
+
+
+tmp<scalarField> curvatureSeparation::calcCosAngle
+(
+    const surfaceScalarField& phi
+) const
+{
+    const fvMesh& mesh = owner().regionMesh();
+    const vectorField nf(mesh.Sf()/mesh.magSf());
+    const unallocLabelList& own = mesh.owner();
+    const unallocLabelList& nbr = mesh.neighbour();
+
+    scalarField phiMax(mesh.nCells(), -GREAT);
+    scalarField cosAngle(mesh.nCells(), 0.0);
+    forAll(nbr, faceI)
+    {
+        label cellO = own[faceI];
+        label cellN = nbr[faceI];
+
+        if (phi[faceI] > phiMax[cellO])
+        {
+            phiMax[cellO] = phi[faceI];
+            cosAngle[cellO] = -gHat_ & nf[faceI];
+        }
+        if (-phi[faceI] > phiMax[cellN])
+        {
+            phiMax[cellN] = -phi[faceI];
+            cosAngle[cellN] = -gHat_ & -nf[faceI];
+        }
+    }
+
+    forAll(phi.boundaryField(), patchI)
+    {
+        const fvsPatchScalarField& phip = phi.boundaryField()[patchI];
+        const fvPatch& pp = phip.patch();
+        const labelList& faceCells = pp.faceCells();
+        const vectorField nf(pp.nf());
+        forAll(phip, i)
+        {
+            label cellI = faceCells[i];
+            if (phip[i] > phiMax[cellI])
+            {
+                phiMax[cellI] = phip[i];
+                cosAngle[cellI] = -gHat_ & nf[i];
+            }
+        }
+    }
+/*
+    // correction for cyclics - use cyclic pairs' face normal instead of
+    // local face normal
+    const fvBoundaryMesh& pbm = mesh.boundary();
+    forAll(phi.boundaryField(), patchI)
+    {
+        if (isA<cyclicPolyPatch>(pbm[patchI]))
+        {
+            const scalarField& phip = phi.boundaryField()[patchI];
+            const vectorField nf(pbm[patchI].nf());
+            const labelList& faceCells = pbm[patchI].faceCells();
+            const label sizeBy2 = pbm[patchI].size()/2;
+
+            for (label face0=0; face0<sizeBy2; face0++)
+            {
+                label face1 = face0 + sizeBy2;
+                label cell0 = faceCells[face0];
+                label cell1 = faceCells[face1];
+
+                // flux leaving half 0, entering half 1
+                if (phip[face0] > phiMax[cell0])
+                {
+                    phiMax[cell0] = phip[face0];
+                    cosAngle[cell0] = -gHat_ & -nf[face1];
+                }
+
+                // flux leaving half 1, entering half 0
+                if (-phip[face1] > phiMax[cell1])
+                {
+                    phiMax[cell1] = -phip[face1];
+                    cosAngle[cell1] = -gHat_ & nf[face0];
+                }
+            }
+        }
+    }
+*/
+    // checks
+    if (debug && mesh.time().outputTime())
+    {
+        volScalarField volCosAngle
+        (
+            IOobject
+            (
+                "cosAngle",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ
+            ),
+            mesh,
+            dimensionedScalar("zero", dimless, 0.0),
+            zeroGradientFvPatchScalarField::typeName
+        );
+        volCosAngle.internalField() = cosAngle;
+        volCosAngle.correctBoundaryConditions();
+        volCosAngle.write();
+    }
+
+    return max(min(cosAngle, scalar(1.0)), scalar(-1.0));
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+curvatureSeparation::curvatureSeparation
+(
+    const surfaceFilmModel& owner,
+    const dictionary& dict
+)
+:
+    injectionModel(type(), owner, dict),
+    gradNHat_(fvc::grad(owner.nHat())),
+    deltaByR1Min_(coeffs().lookupOrDefault<scalar>("deltaByR1Min", 0.0)),
+    definedPatchRadii_(),
+    magG_(mag(owner.g().value())),
+    gHat_(owner.g().value()/magG_)
+{
+    List<Tuple2<word, scalar> > prIn(coeffs().lookup("definedPatchRadii"));
+    const wordList& allPatchNames = owner.regionMesh().boundaryMesh().names();
+
+    DynamicList<Tuple2<label, scalar> > prData(allPatchNames.size());
+
+    labelHashSet uniquePatchIDs;
+
+    forAllReverse(prIn, i)
+    {
+        labelList patchIDs = findStrings(prIn[i].first(), allPatchNames);
+        forAll(patchIDs, j)
+        {
+            const label patchI = patchIDs[j];
+
+            if (!uniquePatchIDs.found(patchI))
+            {
+                const scalar radius = prIn[i].second();
+                prData.append(Tuple2<label, scalar>(patchI, radius));
+
+                uniquePatchIDs.insert(patchI);
+            }
+        }
+    }
+
+    definedPatchRadii_.transfer(prData);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+curvatureSeparation::~curvatureSeparation()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void curvatureSeparation::correct
+(
+    scalarField& availableMass,
+    scalarField& massToInject,
+    scalarField& diameterToInject
+)
+{
+    const kinematicSingleLayer& film =
+        refCast<const kinematicSingleLayer>(this->owner());
+    const fvMesh& mesh = film.regionMesh();
+
+    const volScalarField& delta = film.delta();
+    const volVectorField& U = film.U();
+    const surfaceScalarField& phi = film.phi();
+    const volScalarField& rho = film.rho();
+    const scalarField magSqrU(magSqr(film.U()));
+    const volScalarField& sigma = film.sigma();
+
+    const scalarField invR1(calcInvR1(U));
+    const scalarField cosAngle(calcCosAngle(phi));
+
+    // calculate force balance
+    const scalar Fthreshold = 1e-10;
+    scalarField Fnet(mesh.nCells(), 0.0);
+    scalarField separated(mesh.nCells(), 0.0);
+    forAll(invR1, i)
+    {
+        if ((invR1[i] > 0) && (delta[i]*invR1[i] > deltaByR1Min_))
+        {
+            scalar R1 = 1.0/(invR1[i] + ROOTVSMALL);
+            scalar R2 = R1 + delta[i];
+
+            // inertial force
+            scalar Fi = -delta[i]*rho[i]*magSqrU[i]*72.0/60.0*invR1[i];
+
+            // body force
+            scalar Fb =
+              - 0.5*rho[i]*magG_*invR1[i]*(sqr(R1) - sqr(R2))*cosAngle[i];
+
+            // surface force
+            scalar Fs = sigma[i]/R2;
+
+            Fnet[i] = Fi + Fb + Fs;
+
+            if (Fnet[i] + Fthreshold < 0)
+            {
+                separated[i] = 1.0;
+            }
+        }
+    }
+
+    // inject all available mass
+    massToInject = separated*availableMass;
+    diameterToInject = separated*delta;
+    availableMass -= separated*availableMass;
+
+    if (debug && mesh.time().outputTime())
+    {
+        volScalarField volFnet
+        (
+            IOobject
+            (
+                "Fnet",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ
+            ),
+            mesh,
+            dimensionedScalar("zero", dimForce, 0.0),
+            zeroGradientFvPatchScalarField::typeName
+        );
+        volFnet.internalField() = Fnet;
+        volFnet.correctBoundaryConditions();
+        volFnet.write();
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //

src/regionModels/surfaceFilmModels/submodels/kinematic/injectionModel/curvatureSeparation/curvatureSeparation.H

@@ -0,0 +1,158 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::curvatureSeparation
+
+Description
+    Curvature film separation model
+
+    Assesses film curvature via the mesh geometry and calculates a force
+    balance of the form:
+
+        F_sum = F_inertial + F_body + F_surface
+
+    If F_sum < 0, the film separates. Similarly, if F_sum > 0 the film will
+    remain attached.
+
+    Based on description given by
+        Owen and D. J. Ryley. The flow of thin liquid films around corners.
+        International Journal of Multiphase Flow, 11(1):51-62, 1985.
+
+
+SourceFiles
+    curvatureSeparation.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef curvatureSeparation_H
+#define curvatureSeparation_H
+
+#include "injectionModel.H"
+#include "surfaceFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+/*---------------------------------------------------------------------------*\
+                   Class curvatureSeparation Declaration
+\*---------------------------------------------------------------------------*/
+
+class curvatureSeparation
+:
+    public injectionModel
+{
+private:
+
+    // Private member functions
+
+        //- Disallow default bitwise copy construct
+        curvatureSeparation(const curvatureSeparation&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const curvatureSeparation&);
+
+
+protected:
+
+    // Protected data
+
+        //- Gradient of surface normals
+        volTensorField gradNHat_;
+
+        //- Minimum gravity driven film thickness (non-dimensionalised delta/R1)
+        scalar deltaByR1Min_;
+
+        //- List of radii for patches - if patch not defined, radius
+        // calculated based on mesh geometry
+        List<Tuple2<label, scalar> > definedPatchRadii_;
+
+        //- Magnitude of gravity vector
+        scalar magG_;
+
+        //- Direction of gravity vector
+        vector gHat_;
+
+
+    // Protected Member Functions
+
+        //- Calculate local (inverse) radius of curvature
+        tmp<volScalarField> calcInvR1(const volVectorField& U) const;
+
+        //- Calculate the cosine of the angle between gravity vector and
+        //  cell out flow direction
+        tmp<scalarField> calcCosAngle(const surfaceScalarField& phi) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("curvatureSeparation");
+
+
+    // Constructors
+
+        //- Construct from surface film model
+        curvatureSeparation
+        (
+            const surfaceFilmModel& owner,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~curvatureSeparation();
+
+
+    // Member Functions
+
+        // Evolution
+
+            //- Correct
+            virtual void correct
+            (
+                scalarField& availableMass,
+                scalarField& massToInject,
+                scalarField& diameterToInject
+            );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //