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

src/OSspecific/POSIX/printStack.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -167,6 +167,26 @@ void getSymbolForRaw
 }
 
 
+void error::safePrintStack(std::ostream& os)
+{
+    // Get raw stack symbols
+    void *array[100];
+    size_t size = backtrace(array, 100);
+    char **strings = backtrace_symbols(array, size);
+
+    // See if they contain function between () e.g. "(__libc_start_main+0xd0)"
+    // and see if cplus_demangle can make sense of part before +
+    for (size_t i = 0; i < size; i++)
+    {
+        string msg(strings[i]);
+        fileName programFile;
+        word address;
+
+        os  << '#' << label(i) << '\t' << msg << std::endl;
+    }
+}
+
+
 void error::printStack(Ostream& os)
 {
     // Reads the starting addresses for the dynamically linked libraries

src/OpenFOAM/db/IOstreams/Sstreams/ISstream.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -29,8 +29,7 @@ Description
 
 SourceFiles
     ISstreamI.H
-    ISread.C
-    ISreadToken.C
+    ISstream.C
 
 \*---------------------------------------------------------------------------*/
 

src/OpenFOAM/db/Time/TimeIO.C

@@ -85,8 +85,8 @@ void Foam::Time::readDict()
         {
             if
             (
-                secondaryWriteControl_
-             == wcTimeStep && label(secondaryWriteInterval_) < 1
+                secondaryWriteControl_ == wcTimeStep
+             && label(secondaryWriteInterval_) < 1
             )
             {
                 FatalIOErrorIn("Time::readDict()", controlDict_)

src/OpenFOAM/db/dictionary/dictionary.C

@@ -189,7 +189,7 @@ Foam::dictionary::dictionary
     parent_(parentDict)
 {
     transfer(dict());
-    name() = parentDict.name() + "::" + name();
+    name() = parentDict.name() + '.' + name();
 }
 
 
@@ -622,7 +622,7 @@ Foam::dictionary Foam::dictionary::subOrEmptyDict
         }
         else
         {
-            return dictionary(*this, dictionary(name() + "::" + keyword));
+            return dictionary(*this, dictionary(name() + '.' + keyword));
         }
     }
     else
@@ -690,7 +690,7 @@ bool Foam::dictionary::add(entry* entryPtr, bool mergeEntry)
 
             if (hashedEntries_.insert(entryPtr->keyword(), entryPtr))
             {
-                entryPtr->name() = name() + "::" + entryPtr->keyword();
+                entryPtr->name() = name() + '.' + entryPtr->keyword();
 
                 if (entryPtr->keyword().isPattern())
                 {
@@ -718,7 +718,7 @@ bool Foam::dictionary::add(entry* entryPtr, bool mergeEntry)
 
     if (hashedEntries_.insert(entryPtr->keyword(), entryPtr))
     {
-        entryPtr->name() = name() + "::" + entryPtr->keyword();
+        entryPtr->name() = name() + '.' + entryPtr->keyword();
         IDLList<entry>::append(entryPtr);
 
         if (entryPtr->keyword().isPattern())
@@ -925,7 +925,7 @@ bool Foam::dictionary::changeKeyword
 
     // change name and HashTable, but leave DL-List untouched
     iter()->keyword() = newKeyword;
-    iter()->name() = name() + "::" + newKeyword;
+    iter()->name() = name() + '.' + newKeyword;
     hashedEntries_.erase(oldKeyword);
     hashedEntries_.insert(newKeyword, iter());
 

src/OpenFOAM/db/dictionary/dictionary.H

@@ -117,7 +117,7 @@ public:
         {
             const word scopedName = name_.name();
 
-            string::size_type i = scopedName.rfind(':');
+            string::size_type i = scopedName.rfind('.');
 
             if (i == scopedName.npos)
             {

src/OpenFOAM/db/dictionary/dictionaryIO.C

@@ -37,7 +37,7 @@ Foam::dictionary::dictionary
     Istream& is
 )
 :
-    dictionaryName(parentDict.name() + "::" + name),
+    dictionaryName(parentDict.name() + '.' + name),
     parent_(parentDict)
 {
     read(is);

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

@@ -27,6 +27,7 @@ License
 #include "dictionary.H"
 #include "IFstream.H"
 #include "addToMemberFunctionSelectionTable.H"
+#include "stringOps.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -68,10 +69,14 @@ namespace functionEntries
 
 Foam::fileName Foam::functionEntries::includeEntry::includeFileName
 (
-    Istream& is
+    Istream& is,
+    const dictionary& dict
 )
 {
     fileName fName(is);
+    // Substitute dictionary entries
+    stringOps::inplaceExpand(fName, dict);
+    // Substitute remaining environment variables
     fName.expand();
 
     if (fName.empty() || fName.isAbsolute())
@@ -94,7 +99,7 @@ bool Foam::functionEntries::includeEntry::execute
     Istream& is
 )
 {
-    const fileName fName(includeFileName(is));
+    const fileName fName(includeFileName(is, parentDict));
     IFstream ifs(fName);
 
     if (ifs)
@@ -129,7 +134,7 @@ bool Foam::functionEntries::includeEntry::execute
     Istream& is
 )
 {
-    const fileName fName(includeFileName(is));
+    const fileName fName(includeFileName(is, parentDict));
     IFstream ifs(fName);
 
     if (ifs)

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -77,7 +77,7 @@ protected:
     // Protected Member Functions
 
         //- Read the include fileName from Istream, expand and return
-        static fileName includeFileName(Istream&);
+        static fileName includeFileName(Istream&, const dictionary&);
 
 
 public:

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -69,7 +69,7 @@ bool Foam::functionEntries::includeIfPresentEntry::execute
     Istream& is
 )
 {
-    const fileName fName(includeFileName(is));
+    const fileName fName(includeFileName(is, parentDict));
     IFstream ifs(fName);
 
     if (ifs)
@@ -92,7 +92,7 @@ bool Foam::functionEntries::includeIfPresentEntry::execute
     Istream& is
 )
 {
-    const fileName fName(includeFileName(is));
+    const fileName fName(includeFileName(is, parentDict));
     IFstream ifs(fName);
 
     if (ifs)

src/OpenFOAM/db/dictionary/primitiveEntry/primitiveEntry.C

@@ -81,7 +81,7 @@ Foam::primitiveEntry::primitiveEntry(const keyType& key, const ITstream& is)
     entry(key),
     ITstream(is)
 {
-    name() += "::" + keyword();
+    name() += '.' + keyword();
 }
 
 

src/OpenFOAM/db/dictionary/primitiveEntry/primitiveEntryIO.C

@@ -186,7 +186,7 @@ Foam::primitiveEntry::primitiveEntry
     entry(key),
     ITstream
     (
-        is.name() + "::" + key,
+        is.name() + '.' + key,
         tokenList(10),
         is.format(),
         is.version()
@@ -201,7 +201,7 @@ Foam::primitiveEntry::primitiveEntry(const keyType& key, Istream& is)
     entry(key),
     ITstream
     (
-        is.name() + "::" + key,
+        is.name() + '.' + key,
         tokenList(10),
         is.format(),
         is.version()

src/OpenFOAM/db/error/error.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -162,6 +162,10 @@ public:
         operator dictionary() const;
 
 
+        //- Helper function to print a stack (if OpenFOAM IO not yet
+        //  initialised)
+        static void safePrintStack(std::ostream&);
+
         //- Helper function to print a stack
         static void printStack(Ostream&);
 

src/OpenFOAM/db/runTimeSelection/construction/runTimeSelectionTables.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -82,7 +82,13 @@ Description
         )                                                                     \
         {                                                                     \
             construct##argNames##ConstructorTables();                         \
-            argNames##ConstructorTablePtr_->insert(lookup, New);              \
+            if (!argNames##ConstructorTablePtr_->insert(lookup, New))         \
+            {                                                                 \
+                std::cerr<< "Duplicate entry " << lookup                      \
+                    << " in runtime selection table " << #baseType            \
+                    << std::endl;                                             \
+                error::safePrintStack(std::cerr);                             \
+            }                                                                 \
         }                                                                     \
                                                                               \
         ~add##argNames##ConstructorToTable()                                  \
@@ -167,11 +173,20 @@ Description
         )                                                                     \
         {                                                                     \
             construct##argNames##ConstructorTables();                         \
-            argNames##ConstructorTablePtr_->insert                            \
+            if                                                                \
             (                                                                 \
-                lookup,                                                       \
-                New##baseType                                                 \
-            );                                                                \
+               !argNames##ConstructorTablePtr_->insert                        \
+                (                                                             \
+                    lookup,                                                   \
+                    New##baseType                                             \
+                )                                                             \
+            )                                                                 \
+            {                                                                 \
+                std::cerr<< "Duplicate entry " << lookup                      \
+                    << " in runtime selection table " << #baseType            \
+                    << std::endl;                                             \
+                error::safePrintStack(std::cerr);                             \
+            }                                                                 \
         }                                                                     \
                                                                               \
         ~add##argNames##ConstructorToTable()                                  \

src/OpenFOAM/dimensionSet/dimensionSet.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -105,83 +105,6 @@ void Foam::dimensionSet::reset(const dimensionSet& ds)
 }
 
 
-Foam::string Foam::dimensionSet::asText() const
-{
-    OStringStream buf;
-
-    bool Dimensionless = true;
-
-    for (int Dimension=0; Dimension < dimensionSet::nDimensions-1; ++Dimension)
-    {
-        const scalar& expt = exponents_[Dimension];
-
-        if (expt < smallExponent && expt > -smallExponent)
-        {
-            continue;
-        }
-
-        if (Dimensionless)
-        {
-            Dimensionless = false;
-        }
-        else
-        {
-            buf << ' ';
-        }
-
-        // note: currently only handle SI
-        switch (Dimension)
-        {
-            case MASS:
-                buf << "kg";
-                break;
-
-            case LENGTH:
-                buf << "m";
-                break;
-
-            case TIME:
-                buf << "s";
-                break;
-
-            case TEMPERATURE:
-                buf << "K";
-                break;
-
-            case MOLES:
-                buf << "mol";
-                break;
-
-            case CURRENT:
-                buf << "A";
-                break;
-
-            case LUMINOUS_INTENSITY:
-                buf << "Cd";
-                break;
-
-            default:
-                buf << "??";  // this shouldn't be - flag as being weird
-                break;
-        }
-
-        if (expt != 1)
-        {
-            buf << '^' << expt;
-        }
-    }
-
-    if (Dimensionless)
-    {
-        return "none";
-    }
-    else
-    {
-        return buf.str();
-    }
-}
-
-
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 Foam::scalar Foam::dimensionSet::operator[](const dimensionType type) const
@@ -196,6 +119,18 @@ Foam::scalar& Foam::dimensionSet::operator[](const dimensionType type)
 }
 
 
+Foam::scalar Foam::dimensionSet::operator[](const label type) const
+{
+    return exponents_[type];
+}
+
+
+Foam::scalar& Foam::dimensionSet::operator[](const label type)
+{
+    return exponents_[type];
+}
+
+
 bool Foam::dimensionSet::operator==(const dimensionSet& ds) const
 {
     for (int Dimension=0; Dimension < nDimensions; ++Dimension)

src/OpenFOAM/dimensionSet/dimensionSet.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -40,10 +40,12 @@ SourceFiles
 #ifndef dimensionSet_H
 #define dimensionSet_H
 
-#include "scalar.H"
 #include "bool.H"
 #include "dimensionedScalarFwd.H"
 #include "className.H"
+#include "scalarField.H"
+#include "PtrList.H"
+#include "HashTable.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -53,6 +55,7 @@ namespace Foam
 // Forward declaration of friend functions and operators
 
 class dimensionSet;
+class dimensionSets;
 
 // Friend functions
 
@@ -140,6 +143,64 @@ public:
 
 private:
 
+    // Private classes
+
+        class tokeniser
+        {
+            // Private data
+
+                Istream& is_;
+
+                List<token> tokens_;
+
+                label start_;
+
+                label size_;
+
+            // Private Member Functions
+
+                void push(const token&);
+
+                token pop();
+
+                void unpop(const token&);
+
+        public:
+
+            // Constructors
+
+                 tokeniser(Istream&);
+
+            // Member Functions
+
+                Istream& stream()
+                {
+                    return is_;
+                }
+
+                bool hasToken() const;
+
+                token nextToken();
+
+                void putBack(const token&);
+
+                void splitWord(const word&);
+
+                static bool valid(char c);
+
+                static label priority(const token& t);
+
+        };
+
+
+        dimensionedScalar parse
+        (
+            const label lastPrior,
+            tokeniser& tis,
+            const HashTable<dimensionedScalar>&
+        ) const;
+
+
     // private data
 
         // dimensionSet stored as an array of dimension exponents
@@ -189,14 +250,56 @@ public:
 
         void reset(const dimensionSet&);
 
-        //- Return a text representation for added readability
-        string asText() const;
+
+    // I/O
+
+        //- Read using provided units. Used only in initial parsing
+        Istream& read
+        (
+            Istream& is,
+            scalar& multiplier,
+            const dictionary&
+        );
+
+        //- Read using provided units
+        Istream& read
+        (
+            Istream& is,
+            scalar& multiplier,
+            const HashTable<dimensionedScalar>&
+        );
+
+        //- Read using system units
+        Istream& read
+        (
+            Istream& is,
+            scalar& multiplier
+        );
+
+        //- Write using provided units
+        Ostream& write
+        (
+            Ostream& os,
+            scalar& multiplier,
+            const dimensionSets&
+        ) const;
+
+        //- Write using system units
+        Ostream& write
+        (
+            Ostream& os,
+            scalar& multiplier
+        ) const;
 
 
     // Member operators
 
         scalar operator[](const dimensionType) const;
         scalar& operator[](const dimensionType);
+
+        scalar operator[](const label) const;
+        scalar& operator[](const label);
+
         bool operator==(const dimensionSet&) const;
         bool operator!=(const dimensionSet&) const;
 

src/OpenFOAM/dimensionSet/dimensionSetIO.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,6 +25,7 @@ License
 
 #include "dimensionSet.H"
 #include "IOstreams.H"
+#include "dimensionedScalar.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -34,46 +35,670 @@ Foam::dimensionSet::dimensionSet(Istream& is)
 }
 
 
-// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+Foam::dimensionSet::tokeniser::tokeniser(Istream& is)
+:
+    is_(is),
+    tokens_(100),
+    start_(0),
+    size_(0)
+{}
 
-Foam::Istream& Foam::operator>>(Istream& is, dimensionSet& dset)
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::dimensionSet::tokeniser::push(const token& t)
 {
-    // Read begining of dimensionSet
-    if (token(is) != token::BEGIN_SQR)
+    label end = (start_+size_)%tokens_.size();
+    tokens_[end] = t;
+    if (size_ == tokens_.size())
     {
-        Info<< "expected a " << token::BEGIN_SQR << " in dimensionSet"
-            << endl << "in stream " << is.info() << endl;
+        start_ = tokens_.fcIndex(start_);
+    }
+    else
+    {
+        size_++;
+    }
+}
+
+
+Foam::token Foam::dimensionSet::tokeniser::pop()
+{
+    token t = tokens_[start_];
+    start_ = tokens_.fcIndex(start_);
+    --size_;
+    return t;
+}
+
+
+void Foam::dimensionSet::tokeniser::unpop(const token& t)
+{
+    ++size_;
+    start_ = tokens_.rcIndex(start_);
+    tokens_[start_] = t;
+}
+
+
+bool Foam::dimensionSet::tokeniser::hasToken() const
+{
+    return size_ || is_.good();
+}
+
+
+bool Foam::dimensionSet::tokeniser::valid(char c)
+{
+    return
+    (
+        !isspace(c)
+     && c != '"'   // string quote
+     && c != '\''  // string quote
+     && c != '/'   // div
+     && c != ';'   // end statement
+     && c != '{'   // beg subdict
+     && c != '}'   // end subdict
+     && c != '('   // beg expr
+     && c != ')'   // end expr
+     && c != '['   // beg dim
+     && c != ']'   // end dim
+     && c != '^'   // power
+     && c != '*'   // mult
+    );
+}
+
+
+Foam::label Foam::dimensionSet::tokeniser::priority(const token& t)
+{
+    if (!t.isPunctuation())
+    {
+        return 0;
+    }
+    else if
+    (
+        t.pToken() == token::MULTIPLY
+     || t.pToken() == token::DIVIDE
+    )
+    {
+        return 2;
+    }
+    else if (t.pToken() == '^')
+    {
+        return 3;
+    }
+    else
+    {
+        return 0;
+    }
+}
+
+
+void Foam::dimensionSet::tokeniser::splitWord(const word& w)
+{
+    size_t start = 0;
+    for (size_t i=0; i<w.size(); ++i)
+    {
+        if (!valid(w[i]))
+        {
+            if (i > start)
+            {
+                word subWord = w(start, i-start);
+                if (isdigit(subWord[0]) || subWord[0] == token::SUBTRACT)
+                {
+                    push(token(readLabel(IStringStream(subWord)())));
+                }
+                else
+                {
+                    push(token(subWord));
+                }
+            }
+            if (w[i] != token::SPACE)
+            {
+                if (isdigit(w[i]))
+                {
+                    push(token(readLabel(IStringStream(w[i])())));
+                }
+                else
+                {
+                    push(token::punctuationToken(w[i]));
+                }
+            }
+            start = i+1;
+        }
+    }
+    if (start < w.size())
+    {
+        word subWord = w(start, w.size()-start);
+        if (isdigit(subWord[0]) || subWord[0] == token::SUBTRACT)
+        {
+            push(token(readLabel(IStringStream(subWord)())));
+        }
+        else
+        {
+            push(token(subWord));
+        }
+    }
+}
+
+
+Foam::token Foam::dimensionSet::tokeniser::nextToken()
+{
+    if (size_ == 0)
+    {
+        token t(is_);
+        if (t.isWord())
+        {
+            splitWord(t.wordToken());
+            return pop();
+        }
+        else
+        {
+            return t;
+        }
+    }
+    else
+    {
+        return pop();
+    }
+}
+
+
+void Foam::dimensionSet::tokeniser::putBack(const token& t)
+{
+    if (size_ == 0)
+    {
+        push(t);
+    }
+    else
+    {
+        unpop(t);
+    }
+}
+
+
+Foam::dimensionedScalar Foam::dimensionSet::parse
+(
+    const label lastPrior,
+    tokeniser& tis,
+    const HashTable<dimensionedScalar>& readSet
+) const
+{
+    dimensionedScalar ds("", dimless, 1.0);
+
+    // Get initial token
+    token nextToken(tis.nextToken());
+
+    // Store type of last token read. Used to detect two consecutive
+    // symbols and assume multiplication
+    bool haveReadSymbol = false;
+
+
+    while (true)
+    {
+        if (nextToken.isWord())
+        {
+            const word& unitName = nextToken.wordToken();
+            const dimensionedScalar& unitDim = readSet[unitName];
+            ds.dimensions() *= unitDim.dimensions();
+            ds.value() *= unitDim.value();
+            haveReadSymbol = true;
+        }
+        else if (nextToken.isNumber())
+        {
+            // no dimensions, just value
+            ds.value() *= nextToken.number();
+            haveReadSymbol = true;
+        }
+        else if (nextToken.isPunctuation())
+        {
+            label nextPrior = tokeniser::priority(nextToken);
+
+            if (nextToken.pToken() == token::BEGIN_SQR)
+            {
+                // No idea when this will happen
+                tis.putBack(nextToken);
+                return ds;
+            }
+            else if (nextToken.pToken() == token::END_SQR)
+            {
+                tis.putBack(nextToken);
+                return ds;
+            }
+            else if (nextToken.pToken() == token::BEGIN_LIST)
+            {
+                dimensionedScalar sub(parse(nextPrior, tis, readSet));
+
+                token t = tis.nextToken();
+                if (!t.isPunctuation()  || t.pToken() !=  token::END_LIST)
+                {
+                    FatalIOErrorIn
+                    (
+                        "dimensionSet::parse"
+                        "(const label, tokeniser&"
+                        ", const HashTable<dimensionedScalar>&)",
+                        tis.stream()
+                    )   << "Illegal token " << t << exit(FatalIOError);
+                }
+
+                ds.dimensions() *= sub.dimensions();
+                ds.value() *= sub.value();
+
+                haveReadSymbol = true;
+            }
+            else if (nextToken.pToken() == token::END_LIST)
+            {
+                tis.putBack(nextToken);
+                return ds;
+            }
+            else if (nextToken.pToken() == token::MULTIPLY)
+            {
+                if (nextPrior > lastPrior)
+                {
+                    dimensionedScalar sub(parse(nextPrior, tis, readSet));
+
+                    ds.dimensions() *= sub.dimensions();
+                    ds.value() *= sub.value();
+                }
+                else
+                {
+                    // Restore token
+                    tis.putBack(nextToken);
+                    return ds;
+                }
+                haveReadSymbol = false;
+            }
+            else if (nextToken.pToken() == token::DIVIDE)
+            {
+                if (nextPrior > lastPrior)
+                {
+                    dimensionedScalar sub(parse(nextPrior, tis, readSet));
+
+                    ds.dimensions() /= sub.dimensions();
+                    ds.value() /= sub.value();
+                }
+                else
+                {
+                    tis.putBack(nextToken);
+                    return ds;
+                }
+                haveReadSymbol = false;
+            }
+            else if (nextToken.pToken() == '^')
+            {
+                if (nextPrior > lastPrior)
+                {
+                    dimensionedScalar exp(parse(nextPrior, tis, readSet));
+
+                    ds.dimensions().reset(pow(ds.dimensions(), exp.value()));
+                    ds.value() = Foam::pow(ds.value(), exp.value());
+                }
+                else
+                {
+                    tis.putBack(nextToken);
+                    return ds;
+                }
+                haveReadSymbol = false;
+            }
+            else
+            {
+                FatalIOErrorIn
+                (
+                    "dimensionSet::parse"
+                    "(const label, tokeniser&"
+                    ", const HashTable<dimensionedScalar>&)",
+                    tis.stream()
+                )   << "Illegal token " << nextToken << exit(FatalIOError);
+            }
+        }
+        else
+        {
+            FatalIOErrorIn
+            (
+                "dimensionSet::parse"
+                "(const label, tokeniser&"
+                ", const HashTable<dimensionedScalar>&)",
+                tis.stream()
+            )   << "Illegal token " << nextToken << exit(FatalIOError);
+        }
+
+
+        if (!tis.hasToken())
+        {
+            break;
+        }
+
+        nextToken = tis.nextToken();
+        if (nextToken.error())
+        {
+            break;
+        }
+
+        if (haveReadSymbol && (nextToken.isWord() || nextToken.isNumber()))
+        {
+            // Two consecutive symbols. Assume multiplication
+            tis.putBack(nextToken);
+            nextToken = token(token::MULTIPLY);
+        }
     }
 
-    // Read first five dimensions
-    for (int Dimension=0; Dimension<dimensionSet::CURRENT; Dimension++)
+    return ds;
+}
+
+
+Foam::Istream& Foam::dimensionSet::read
+(
+    Istream& is,
+    scalar& multiplier,
+    const HashTable<dimensionedScalar>& readSet
+)
+{
+    multiplier = 1.0;
+
+    // Read begining of dimensionSet
+    token startToken(is);
+
+    if (startToken != token::BEGIN_SQR)
     {
-        is >> dset.exponents_[Dimension];
+        FatalIOErrorIn
+        (
+            "dimensionSet::read"
+            "(Istream&, scalar&, const HashTable<dimensionedScalar>&)",
+            is
+        )   << "expected a " << token::BEGIN_SQR << " in dimensionSet"
+            << endl << "in stream " << is.info()
+            << exit(FatalIOError);
     }
 
     // Read next token
     token nextToken(is);
 
-    // If next token is another number
-    // read last two dimensions
-    // and then read another token for the end of the dimensionSet
-    if (nextToken.isNumber())
+    if (!nextToken.isNumber())
     {
-        dset.exponents_[dimensionSet::CURRENT] = nextToken.number();
-        is >> dset.exponents_[dimensionSet::LUMINOUS_INTENSITY];
-        is >> nextToken;
+        is.putBack(nextToken);
+
+        tokeniser tis(is);
+
+        dimensionedScalar ds(parse(0, tis, readSet));
+
+        multiplier = ds.value();
+        for (int i=0; i < dimensionSet::nDimensions; ++i)
+        {
+            exponents_[i] = ds.dimensions()[i];
+        }
     }
     else
     {
-        dset.exponents_[dimensionSet::CURRENT] = 0;
-        dset.exponents_[dimensionSet::LUMINOUS_INTENSITY] = 0;
+        // Read first five dimensions
+        exponents_[dimensionSet::MASS] = nextToken.number();
+        for (int Dimension=1; Dimension<dimensionSet::CURRENT; Dimension++)
+        {
+            is >> exponents_[Dimension];
+        }
+
+        // Read next token
+        token nextToken(is);
+
+        // If next token is another number
+        // read last two dimensions
+        // and then read another token for the end of the dimensionSet
+        if (nextToken.isNumber())
+        {
+            exponents_[dimensionSet::CURRENT] = nextToken.number();
+            is >> nextToken;
+            exponents_[dimensionSet::LUMINOUS_INTENSITY] = nextToken.number();
+            is >> nextToken;
+        }
+        else
+        {
+            exponents_[dimensionSet::CURRENT] = 0;
+            exponents_[dimensionSet::LUMINOUS_INTENSITY] = 0;
+        }
+
+        // Check end of dimensionSet
+        if (nextToken != token::END_SQR)
+        {
+            FatalIOErrorIn
+            (
+                "dimensionSet::read"
+                "(Istream&, scalar&, const HashTable<dimensionedScalar>&)",
+                is
+            )   << "expected a " << token::END_SQR << " in dimensionSet "
+                << endl << "in stream " << is.info()
+                << exit(FatalIOError);
+        }
+    }
+    // Check state of Istream
+    is.check("Istream& operator>>(Istream&, dimensionSet&)");
+
+    return is;
+}
+
+
+Foam::Istream& Foam::dimensionSet::read
+(
+    Istream& is,
+    scalar& multiplier
+)
+{
+    return read(is, multiplier, unitSet());
+}
+
+
+Foam::Istream& Foam::dimensionSet::read
+(
+    Istream& is,
+    scalar& multiplier,
+    const dictionary& readSet
+)
+{
+    multiplier = 1.0;
+
+    // Read begining of dimensionSet
+    token startToken(is);
+
+    if (startToken != token::BEGIN_SQR)
+    {
+        FatalIOErrorIn
+        (
+            "dimensionSet::read"
+            "(Istream&, scalar&, const dictionary&)",
+            is
+        )   << "expected a " << token::BEGIN_SQR << " in dimensionSet"
+            << endl << "in stream " << is.info()
+            << exit(FatalIOError);
     }
 
-    // Check end of dimensionSet
-    if (nextToken != token::END_SQR)
+    // Read next token
+    token nextToken(is);
+
+    if (nextToken.isWord())
     {
-        Info<< "expected a " << token::END_SQR << " in dimensionSet"
-            << endl << "in stream " << is.info() << endl;
+        bool continueParsing = true;
+        do
+        {
+            word symbolPow = nextToken.wordToken();
+            if (symbolPow[symbolPow.size()-1] == token::END_SQR)
+            {
+                symbolPow = symbolPow(0, symbolPow.size()-1);
+                continueParsing = false;
+            }
+
+
+            // Parse unit
+            dimensionSet symbolSet(dimless);
+
+            size_t index = symbolPow.find('^');
+            if (index != string::npos)
+            {
+                word symbol = symbolPow(0, index);
+                word exp = symbolPow(index+1, symbolPow.size()-index+1);
+                scalar exponent = readScalar(IStringStream(exp)());
+
+                dimensionedScalar s;
+                s.read(readSet[symbol], readSet);
+
+                symbolSet.reset(pow(s.dimensions(), exponent));
+                multiplier *= Foam::pow(s.value(), exponent);
+            }
+            else
+            {
+                dimensionedScalar s;
+                s.read(readSet[symbolPow], readSet);
+
+                symbolSet.reset(s.dimensions());
+                multiplier *= s.value();
+            }
+
+            // Add dimensions without checking
+            for (int i=0; i < dimensionSet::nDimensions; ++i)
+            {
+                exponents_[i] += symbolSet[i];
+            }
+
+            if (continueParsing)
+            {
+                nextToken = token(is);
+
+                if (!nextToken.isWord() || nextToken == token::END_SQR)
+                {
+                    continueParsing = false;
+                }
+            }
+        }
+        while (continueParsing);
+    }
+    else
+    {
+        // Read first five dimensions
+        exponents_[dimensionSet::MASS] = nextToken.number();
+        for (int Dimension=1; Dimension<dimensionSet::CURRENT; Dimension++)
+        {
+            is >> exponents_[Dimension];
+        }
+
+        // Read next token
+        token nextToken(is);
+
+        // If next token is another number
+        // read last two dimensions
+        // and then read another token for the end of the dimensionSet
+        if (nextToken.isNumber())
+        {
+            exponents_[dimensionSet::CURRENT] = nextToken.number();
+            is >> nextToken;
+            exponents_[dimensionSet::LUMINOUS_INTENSITY] = nextToken.number();
+            is >> nextToken;
+        }
+        else
+        {
+            exponents_[dimensionSet::CURRENT] = 0;
+            exponents_[dimensionSet::LUMINOUS_INTENSITY] = 0;
+        }
+
+        // Check end of dimensionSet
+        if (nextToken != token::END_SQR)
+        {
+            FatalIOErrorIn
+            (
+                "dimensionSet::read"
+                "(Istream&, scalar&, const dictionary&)",
+                is
+            )   << "expected a " << token::END_SQR << " in dimensionSet "
+                << endl << "in stream " << is.info()
+                << exit(FatalIOError);
+        }
+    }
+
+    // Check state of Istream
+    is.check("Istream& operator>>(Istream&, dimensionSet&)");
+
+    return is;
+}
+
+
+Foam::Ostream& Foam::dimensionSet::write
+(
+    Ostream& os,
+    scalar& multiplier,
+    const dimensionSets& writeUnits
+) const
+{
+    multiplier = 1.0;
+
+    os << token::BEGIN_SQR;
+
+    if (writeUnits.valid())
+    {
+        scalarField exponents(dimensionSet::nDimensions);
+        for (int d=0; d<dimensionSet::nDimensions; d++)
+        {
+            exponents[d] = exponents_[d];
+        }
+        writeUnits.coefficients(exponents);
+
+        bool hasPrinted = false;
+        forAll(exponents, i)
+        {
+            if (mag(exponents[i]) > smallExponent)
+            {
+                const dimensionedScalar& ds = writeUnits.units()[i];
+
+                if (hasPrinted)
+                {
+                    os  << token::SPACE;
+                }
+                hasPrinted = true;
+                os  << ds.name();
+                if (mag(exponents[i]-1) > smallExponent)
+                {
+                    os  << '^' << exponents[i];
+
+                    multiplier *= Foam::pow(ds.value(), exponents[i]);
+                }
+                else
+                {
+                    multiplier *= ds.value();
+                }
+            }
+        }
+    }
+    else
+    {
+        for (int d=0; d<dimensionSet::nDimensions-1; d++)
+        {
+            os << exponents_[d] << token::SPACE;
+        }
+        os << exponents_[dimensionSet::nDimensions-1];
+    }
+
+    os  << token::END_SQR;
+
+    // Check state of Ostream
+    os.check("Ostream& operator<<(Ostream&, const dimensionSet&)");
+
+    return os;
+}
+
+
+Foam::Ostream& Foam::dimensionSet::write
+(
+    Ostream& os,
+    scalar& multiplier
+) const
+{
+    return write(os, multiplier, writeUnitSet());
+}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+Foam::Istream& Foam::operator>>(Istream& is, dimensionSet& dset)
+{
+    scalar multiplier;
+    dset.read(is, multiplier);
+
+    if (mag(multiplier-1.0) > dimensionSet::smallExponent)
+    {
+        FatalIOErrorIn("Foam::operator>>(Istream&, dimensionSet&)", is)
+            << "Cannot use scaled units in dimensionSet"
+            << exit(FatalIOError);
     }
 
     // Check state of Istream
@@ -85,13 +710,8 @@ Foam::Istream& Foam::operator>>(Istream& is, dimensionSet& dset)
 
 Foam::Ostream& Foam::operator<<(Ostream& os, const dimensionSet& dset)
 {
-    os << token::BEGIN_SQR;
-
-    for (int Dimension=0; Dimension<dimensionSet::nDimensions-1; Dimension++)
-    {
-        os << dset.exponents_[Dimension] << token::SPACE;
-    }
-    os << dset.exponents_[dimensionSet::nDimensions-1] << token::END_SQR;
+    scalar multiplier;
+    dset.write(os, multiplier);
 
     // Check state of Ostream
     os.check("Ostream& operator<<(Ostream&, const dimensionSet&)");

src/OpenFOAM/dimensionSet/dimensionSets.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,34 +24,210 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "dimensionSet.H"
+#include "dimensionedScalar.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-const Foam::dimensionSet Foam::dimless(0, 0, 0, 0, 0, 0, 0);
+namespace Foam
+{
 
-const Foam::dimensionSet Foam::dimMass(1, 0, 0, 0, 0, 0, 0);
-const Foam::dimensionSet Foam::dimLength(0, 1, 0, 0, 0, 0, 0);
-const Foam::dimensionSet Foam::dimTime(0, 0, 1, 0, 0, 0, 0);
-const Foam::dimensionSet Foam::dimTemperature(0, 0, 0, 1, 0, 0, 0);
-const Foam::dimensionSet Foam::dimMoles(0, 0, 0, 0, 1, 0, 0);
-const Foam::dimensionSet Foam::dimCurrent(0, 0, 0, 0, 0, 1, 0);
-const Foam::dimensionSet Foam::dimLuminousIntensity(0, 0, 0, 0, 0, 0, 1);
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-const Foam::dimensionSet Foam::dimArea(sqr(dimLength));
-const Foam::dimensionSet Foam::dimVolume(pow3(dimLength));
-const Foam::dimensionSet Foam::dimVol(dimVolume);
+dictionary* dimensionSystemsPtr_(NULL);
 
-const Foam::dimensionSet Foam::dimVelocity(dimLength/dimTime);
-const Foam::dimensionSet Foam::dimAcceleration(dimVelocity/dimTime);
+dictionary& dimensionSystems()
+{
+    return debug::switchSet
+    (
+        "DimensionSets",
+        dimensionSystemsPtr_
+    );
+}
 
-const Foam::dimensionSet Foam::dimDensity(dimMass/dimVolume);
-const Foam::dimensionSet Foam::dimForce(dimMass*dimAcceleration);
-const Foam::dimensionSet Foam::dimEnergy(dimForce*dimLength);
-const Foam::dimensionSet Foam::dimPower(dimEnergy/dimTime);
 
-const Foam::dimensionSet Foam::dimPressure(dimForce/dimArea);
-const Foam::dimensionSet Foam::dimGasConstant(dimEnergy/dimMass/dimTemperature);
-const Foam::dimensionSet Foam::dimSpecificHeatCapacity(dimGasConstant);
+/* * * * * * * * * * * * * * * Static Member Data  * * * * * * * * * * * * * */
 
+autoPtr<HashTable<dimensionedScalar> > unitSetPtr_;
+autoPtr<dimensionSets> writeUnitSetPtr_;
+
+
+const HashTable<dimensionedScalar>& unitSet()
+{
+    if (!unitSetPtr_.valid())
+    {
+        const dictionary& dict = dimensionSystems();
+
+        if (!dict.found("unitSet"))
+        {
+            FatalIOErrorIn("unitSet()", dict)
+                << "Cannot find unitSet in dictionary " << dict.name()
+                << exit(FatalIOError);
+        }
+
+        const word unitSetCoeffs(word(dict.lookup("unitSet")) + "Coeffs");
+
+        if (!dict.found(unitSetCoeffs))
+        {
+            FatalIOErrorIn("unitSet()", dict)
+                << "Cannot find " << unitSetCoeffs << " in dictionary "
+                << dict.name() << exit(FatalIOError);
+        }
+
+        const dictionary& unitDict = dict.subDict(unitSetCoeffs);
+
+        unitSetPtr_.reset
+        (
+            new HashTable<dimensionedScalar>(unitDict.size())
+        );
+
+        forAllConstIter(dictionary, unitDict, iter)
+        {
+            if (iter().keyword() != "writeUnits")
+            {
+                dimensionedScalar dt;
+                dt.read(iter().stream(), unitDict);
+                bool ok = unitSetPtr_->insert(iter().keyword(), dt);
+                if (!ok)
+                {
+                    FatalIOErrorIn("unitSet()", dict)
+                        << "Duplicate unit " << iter().keyword()
+                        << " in DimensionSets dictionary"
+                        << exit(FatalIOError);
+                }
+            }
+        }
+
+        wordList writeUnitNames
+        (
+            unitDict.lookupOrDefault<wordList>
+            (
+                "writeUnits",
+                wordList(0)
+            )
+        );
+
+        writeUnitSetPtr_.reset
+        (
+            new dimensionSets
+            (
+                unitSetPtr_(),
+                writeUnitNames
+            )
+        );
+
+        if (writeUnitNames.size() != 0 && writeUnitNames.size() != 7)
+        {
+            FatalIOErrorIn("unitSet()", dict)
+                << "Cannot find entry \"writeUnits\" in " << unitDict.name()
+                << " or it is not a wordList of size 7"
+                << exit(FatalIOError);
+        }
+    }
+    return unitSetPtr_();
+}
+
+
+const dimensionSets& writeUnitSet()
+{
+    if (!writeUnitSetPtr_.valid())
+    {
+        (void)unitSet();
+    }
+    return writeUnitSetPtr_();
+}
+
+
+const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0);
+
+const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0);
+const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0);
+const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0);
+const dimensionSet dimTemperature(0, 0, 0, 1, 0, 0, 0);
+const dimensionSet dimMoles(0, 0, 0, 0, 1, 0, 0);
+const dimensionSet dimCurrent(0, 0, 0, 0, 0, 1, 0);
+const dimensionSet dimLuminousIntensity(0, 0, 0, 0, 0, 0, 1);
+
+const dimensionSet dimArea(sqr(dimLength));
+const dimensionSet dimVolume(pow3(dimLength));
+const dimensionSet dimVol(dimVolume);
+
+const dimensionSet dimVelocity(dimLength/dimTime);
+const dimensionSet dimAcceleration(dimVelocity/dimTime);
+
+const dimensionSet dimDensity(dimMass/dimVolume);
+const dimensionSet dimForce(dimMass*dimAcceleration);
+const dimensionSet dimEnergy(dimForce*dimLength);
+const dimensionSet dimPower(dimEnergy/dimTime);
+
+const dimensionSet dimPressure(dimForce/dimArea);
+const dimensionSet dimGasConstant(dimEnergy/dimMass/dimTemperature);
+const dimensionSet dimSpecificHeatCapacity(dimGasConstant);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::dimensionSets::dimensionSets
+(
+    const HashTable<dimensionedScalar>& units,
+    const wordList& unitNames
+)
+:
+    units_(unitNames.size()),
+    conversion_(unitNames.size()),
+    conversionPivots_(unitNames.size()),
+    valid_(false)
+{
+    forAll(unitNames, i)
+    {
+        units_.set
+        (
+            i,
+            new dimensionedScalar
+            (
+                units[unitNames[i]]
+            )
+        );
+    }
+
+    if (unitNames.size() == 7)
+    {
+        valid_ = true;
+
+        // Determine conversion from basic units to write units
+        for (label rowI = 0; rowI < conversion_.n(); rowI++)
+        {
+            scalar* row = conversion_[rowI];
+
+            for (label columnI = 0; columnI < conversion_.m(); columnI++)
+            {
+                const dimensionedScalar& dSet = units_[columnI];
+                row[columnI] = dSet.dimensions()[rowI];
+            }
+        }
+        conversionPivots_.setSize(conversion_.n());
+        LUDecompose(conversion_, conversionPivots_);
+
+        //- possibly some optimisation here to detect identity matri
+        //if
+        //(
+        //    conversionPivots_ == identity(conversionPivots_.size())
+        // && conversion_ == I)
+        //)
+        //{
+        //    identity_ = true;
+        //}
+    }
+}
+
+
+void Foam::dimensionSets::coefficients(scalarField& exponents) const
+{
+    LUBacksubstitute(conversion_, conversionPivots_, exponents);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
 
 // ************************************************************************* //

src/OpenFOAM/dimensionSet/dimensionSets.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -36,6 +36,10 @@ SourceFiles
 #ifndef dimensionSets_H
 #define dimensionSets_H
 
+#include "scalarMatrices.H"
+#include "dimensionedScalarFwd.H"
+#include "PtrList.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -69,6 +73,65 @@ extern const dimensionSet dimGasConstant;
 extern const dimensionSet dimSpecificHeatCapacity;
 
 
+
+class dimensionSets
+{
+    // Private data
+
+        //- Set of dimensions
+        PtrList<dimensionedScalar> units_;
+
+        //- LU decomposition of dimensions
+        scalarSquareMatrix conversion_;
+
+        //- See above
+        labelList conversionPivots_;
+
+        //- Is LU decomposition valid
+        bool valid_;
+
+public:
+
+    // Constructors
+
+        //- Construct from all units and set of units to use for inversion
+        //  (writing)
+        dimensionSets
+        (
+            const HashTable<dimensionedScalar>&,
+            const wordList& unitNames
+        );
+
+    // Member functions
+
+        //- Return the units
+        const PtrList<dimensionedScalar>& units() const
+        {
+            return units_;
+        }
+
+        //- Is there a valid inverse of the selected unit
+        bool valid() const
+        {
+            return valid_;
+        }
+
+        //- (if valid) obtain set of coefficients of unitNames
+        void coefficients(scalarField&) const;
+
+};
+
+
+//- Top level dictionary
+dictionary& dimensionSystems();
+
+//- Set of all dimensions
+const HashTable<dimensionedScalar>& unitSet();
+
+//- Set of units
+const dimensionSets& writeUnitSet();
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/dimensionedTypes/dimensionedType/dimensionedType.C

@@ -97,10 +97,10 @@ dimensioned<Type>::dimensioned
     Istream& is
 )
 :
-    name_(is),
-    dimensions_(is),
-    value_(pTraits<Type>(is))
-{}
+    dimensions_(dimless)
+{
+    read(is);
+}
 
 
 template <class Type>
@@ -111,9 +111,13 @@ dimensioned<Type>::dimensioned
 )
 :
     name_(name),
-    dimensions_(is),
-    value_(pTraits<Type>(is))
-{}
+    dimensions_(dimless)
+{
+    scalar multiplier;
+    dimensions_.read(is, multiplier);
+    is >> value_;
+    value_ *= multiplier;
+}
 
 
 template <class Type>
@@ -141,24 +145,29 @@ dimensioned<Type>::dimensioned
     }
 
     // If the dimensions are provided compare with the argument
+    scalar multiplier = 1.0;
+
     if (nextToken == token::BEGIN_SQR)
     {
-        dimensionSet dims(is);
+        dimensionSet dims(dimless);
+        dims.read(is, multiplier);
 
         if (dims != dimensions_)
         {
-            FatalErrorIn
+            FatalIOErrorIn
             (
                 "dimensioned<Type>::dimensioned"
-                "(const word&, const dimensionSet&, Istream&)"
+                "(const word&, const dimensionSet&, Istream&)",
+                is
             ) << "The dimensions " << dims
               << " provided do not match the required dimensions "
               << dimensions_
-              << abort(FatalError);
+              << abort(FatalIOError);
         }
     }
 
     is >> value_;
+    value_ *= multiplier;
 }
 
 
@@ -247,6 +256,83 @@ bool dimensioned<Type>::readIfPresent(const dictionary& dict)
 }
 
 
+template <class Type>
+Foam::Istream& dimensioned<Type>::read(Istream& is, const dictionary& readSet)
+{
+    // Read name
+    is >> name_;
+
+    // Read dimensionSet + multiplier
+    scalar mult;
+    dimensions_.read(is, mult, readSet);
+
+    // Read value
+    is >> value_;
+    value_ *= mult;
+
+    // Check state of Istream
+    is.check
+    (
+        "Istream& dimensioned<Type>::read(Istream& is, const dictionary&)"
+    );
+
+    return is;
+}
+
+
+template <class Type>
+Foam::Istream& dimensioned<Type>::read
+(
+    Istream& is,
+    const HashTable<dimensionedScalar>& readSet
+)
+{
+    // Read name
+    is >> name_;
+
+    // Read dimensionSet + multiplier
+    scalar mult;
+    dimensions_.read(is, mult, readSet);
+
+    // Read value
+    is >> value_;
+    value_ *= mult;
+
+    // Check state of Istream
+    is.check
+    (
+        "Istream& dimensioned<Type>::read"
+        "(Istream& is, const HashTable<dimensionedScalar>&)"
+    );
+
+    return is;
+}
+
+
+template <class Type>
+Foam::Istream& dimensioned<Type>::read(Istream& is)
+{
+    // Read name
+    is >> name_;
+
+    // Read dimensionSet + multiplier
+    scalar mult;
+    dimensions_.read(is, mult);
+
+    // Read value
+    is >> value_;
+    value_ *= mult;
+
+    // Check state of Istream
+    is.check
+    (
+        "Istream& dimensioned<Type>::read(Istream& is)"
+    );
+
+    return is;
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template <class Type>
@@ -414,13 +500,15 @@ Istream& operator>>(Istream& is, dimensioned<Type>& dt)
     }
 
     // If the dimensions are provided reset the dimensions to those read
+    scalar multiplier = 1.0;
     if (nextToken == token::BEGIN_SQR)
     {
-        is >> dt.dimensions_;
+        dt.dimensions_.read(is, multiplier);
     }
 
     // Read the value
     is >> dt.value_;
+    dt.value_ *= multiplier;
 
     // Check state of Istream
     is.check("Istream& operator>>(Istream&, dimensioned<Type>&)");
@@ -432,10 +520,17 @@ Istream& operator>>(Istream& is, dimensioned<Type>& dt)
 template <class Type>
 Ostream& operator<<(Ostream& os, const dimensioned<Type>& dt)
 {
-    // do a stream write op for a dimensions()et
-    os  << dt.name() << token::SPACE
-        << dt.dimensions() << token::SPACE
-        << dt.value();
+    // Write the name
+    os << dt.name() << token::SPACE;
+
+    // Write the dimensions
+    scalar mult;
+    dt.dimensions().write(os, mult);
+
+    os << token::SPACE;
+
+    // Write the value
+    os << dt.value()/mult;
 
     // Check state of Ostream
     os.check("Ostream& operator<<(Ostream&, const dimensioned<Type>&)");

src/OpenFOAM/dimensionedTypes/dimensionedType/dimensionedType.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -163,6 +163,17 @@ public:
         bool readIfPresent(const dictionary&);
 
 
+    // I/O
+
+        //- Read using provided units. Only used during startup.
+        Istream& read(Istream& is, const dictionary&);
+
+        //- Read using provided units
+        Istream& read(Istream& is, const HashTable<dimensionedScalar>&);
+
+        //- Read using system units
+        Istream& read(Istream& is);
+
     // Member operators
 
         //- Return a component as a dimensioned<cmptType>

src/OpenFOAM/fields/UniformDimensionedFields/UniformDimensionedField.C

@@ -48,8 +48,10 @@ Foam::UniformDimensionedField<Type>::UniformDimensionedField
     )
     {
         dictionary dict(readStream(typeName));
-        this->dimensions().reset(dict.lookup("dimensions"));
+        scalar multiplier;
+        this->dimensions().read(dict.lookup("dimensions"), multiplier);
         dict.lookup("value") >> this->value();
+        this->value() *= multiplier;
     }
 }
 
@@ -75,8 +77,10 @@ Foam::UniformDimensionedField<Type>::UniformDimensionedField
     dimensioned<Type>(regIOobject::name(), dimless, pTraits<Type>::zero)
 {
     dictionary dict(readStream(typeName));
-    this->dimensions().reset(dict.lookup("dimensions"));
+    scalar multiplier;
+    this->dimensions().read(dict.lookup("dimensions"), multiplier);
     dict.lookup("value") >> this->value();
+    this->value() *= multiplier;
 }
 
 
@@ -92,9 +96,11 @@ Foam::UniformDimensionedField<Type>::~UniformDimensionedField()
 template<class Type>
 bool Foam::UniformDimensionedField<Type>::writeData(Ostream& os) const
 {
-    os.writeKeyword("dimensions") << this->dimensions() << token::END_STATEMENT
+    scalar multiplier;
+    os.writeKeyword("dimensions");
+    this->dimensions().write(os, multiplier) << token::END_STATEMENT
         << nl;
-    os.writeKeyword("value") << this->value() << token::END_STATEMENT
+    os.writeKeyword("value") << this->value()/multiplier << token::END_STATEMENT
         << nl << nl;
 
     return (os.good());

src/OpenFOAM/fields/pointPatchFields/basic/fixedValue/fixedValuePointPatchField.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -121,6 +121,17 @@ public:
         }
 
 
+    // Member functions
+
+        // Access
+
+            //- Return true if this patch field fixes a value
+            virtual bool fixesValue() const
+            {
+                return true;
+            }
+
+
     // Member operators
 
         // Disable assignment operators

src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchField.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -278,6 +278,12 @@ public:
                 return internalField_;
             }
 
+            //- Return true if this patch field fixes a value
+            virtual bool fixesValue() const
+            {
+                return false;
+            }
+
             //- Return true if this patch field is coupled
             virtual bool coupled() const
             {

src/OpenFOAM/memory/tmp/tmp.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -79,6 +79,9 @@ public:
         //- Construct copy and increment reference count
         inline tmp(const tmp<T>&);
 
+        //- Construct copy transferring content of temporary if required
+        inline tmp(const tmp<T>&, bool allowTransfer);
+
 
     //- Destructor, delete object when reference count == 0
     inline ~tmp();

src/OpenFOAM/memory/tmp/tmpI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -68,6 +68,38 @@ inline Foam::tmp<T>::tmp(const tmp<T>& t)
 }
 
 
+template<class T>
+inline Foam::tmp<T>::tmp(const tmp<T>& t, bool allowTransfer)
+:
+    isTmp_(t.isTmp_),
+    ptr_(t.ptr_),
+    ref_(t.ref_)
+{
+    if (isTmp_)
+    {
+        if (allowTransfer)
+        {
+            const_cast<tmp<T>&>(t).ptr_ = 0;
+        }
+        else
+        {
+            if (ptr_)
+            {
+                ptr_->operator++();
+            }
+            else
+            {
+                FatalErrorIn
+                (
+                    "Foam::tmp<T>::tmp(const tmp<T>&, bool allowTransfer)"
+                )   << "attempted copy of a deallocated temporary"
+                    << abort(FatalError);
+            }
+        }
+    }
+}
+
+
 template<class T>
 inline Foam::tmp<T>::~tmp()
 {

src/OpenFOAM/primitives/strings/stringOps/stringOps.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -322,7 +322,12 @@ Foam::string& Foam::stringOps::inplaceExpand
 
 
                 // lookup in the dictionary
-                const entry* ePtr = dict.lookupEntryPtr(varName, true, true);
+                const entry* ePtr = dict.lookupScopedEntryPtr
+                (
+                    varName,
+                    true,
+                    false   // wildcards disabled. See primitiveEntry
+                );
 
                 // if defined - copy its entries
                 if (ePtr)

src/fieldSources/basicSource/pressureGradientExplicitSource/pressureGradientExplicitSource.C

@@ -210,7 +210,7 @@ void Foam::pressureGradientExplicitSource::setValue
             (
                 IOobject
                 (
-                    name_ + "::invA",
+                    name_ + ".invA",
                     mesh_.time().timeName(),
                     mesh_,
                     IOobject::NO_READ,

src/finiteVolume/Make/files

@@ -84,7 +84,6 @@ fvPatchFields = fields/fvPatchFields
 $(fvPatchFields)/fvPatchField/fvPatchFields.C
 
 basicFvPatchFields = $(fvPatchFields)/basic
-$(basicFvPatchFields)/basicSymmetry/basicSymmetryFvPatchFields.C
 $(basicFvPatchFields)/basicSymmetry/basicSymmetryFvPatchScalarField.C
 $(basicFvPatchFields)/calculated/calculatedFvPatchFields.C
 $(basicFvPatchFields)/coupled/coupledFvPatchFields.C

src/finiteVolume/fields/fvPatchFields/basic/basicSymmetry/basicSymmetryFvPatchField.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -58,10 +58,6 @@ class basicSymmetryFvPatchField
 
 public:
 
-    //- Runtime type information
-    TypeName(symmetryFvPatch::typeName_());
-
-
     // Constructors
 
         //- Construct from patch and internal field

src/finiteVolume/finiteVolume/fvSchemes/fvSchemes.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -269,98 +269,98 @@ Foam::fvSchemes::fvSchemes(const objectRegistry& obr)
     (
         ITstream
         (
-            objectPath() + "::ddtSchemes",
+            objectPath() + ".ddtSchemes",
             tokenList()
         )()
     ),
     defaultDdtScheme_
     (
-        ddtSchemes_.name() + "::default",
+        ddtSchemes_.name() + ".default",
         tokenList()
     ),
     d2dt2Schemes_
     (
         ITstream
         (
-            objectPath() + "::d2dt2Schemes",
+            objectPath() + ".d2dt2Schemes",
             tokenList()
         )()
     ),
     defaultD2dt2Scheme_
     (
-        d2dt2Schemes_.name() + "::default",
+        d2dt2Schemes_.name() + ".default",
         tokenList()
     ),
     interpolationSchemes_
     (
         ITstream
         (
-            objectPath() + "::interpolationSchemes",
+            objectPath() + ".interpolationSchemes",
             tokenList()
         )()
     ),
     defaultInterpolationScheme_
     (
-        interpolationSchemes_.name() + "::default",
+        interpolationSchemes_.name() + ".default",
         tokenList()
     ),
     divSchemes_
     (
         ITstream
         (
-            objectPath() + "::divSchemes",
+            objectPath() + ".divSchemes",
             tokenList()
         )()
     ),
     defaultDivScheme_
     (
-        divSchemes_.name() + "::default",
+        divSchemes_.name() + ".default",
         tokenList()
     ),
     gradSchemes_
     (
         ITstream
         (
-            objectPath() + "::gradSchemes",
+            objectPath() + ".gradSchemes",
             tokenList()
         )()
     ),
     defaultGradScheme_
     (
-        gradSchemes_.name() + "::default",
+        gradSchemes_.name() + ".default",
         tokenList()
     ),
     snGradSchemes_
     (
         ITstream
         (
-            objectPath() + "::snGradSchemes",
+            objectPath() + ".snGradSchemes",
             tokenList()
         )()
     ),
     defaultSnGradScheme_
     (
-        snGradSchemes_.name() + "::default",
+        snGradSchemes_.name() + ".default",
         tokenList()
     ),
     laplacianSchemes_
     (
         ITstream
         (
-            objectPath() + "::laplacianSchemes",
+            objectPath() + ".laplacianSchemes",
             tokenList()
         )()
     ),
     defaultLaplacianScheme_
     (
-        laplacianSchemes_.name() + "::default",
+        laplacianSchemes_.name() + ".default",
         tokenList()
     ),
     fluxRequired_
     (
         ITstream
         (
-            objectPath() + "::fluxRequired",
+            objectPath() + ".fluxRequired",
             tokenList()
         )()
     ),

src/finiteVolume/finiteVolume/fvc/fvcMeshPhi.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -148,4 +148,39 @@ void Foam::fvc::makeAbsolute
 }
 
 
+Foam::tmp<Foam::surfaceScalarField> Foam::fvc::absolute
+(
+    const tmp<surfaceScalarField>& tphi,
+    const volVectorField& U
+)
+{
+    if (tphi().mesh().moving())
+    {
+        return tphi + fvc::meshPhi(U);
+    }
+    else
+    {
+        return tmp<surfaceScalarField>(tphi, true);
+    }
+}
+
+
+Foam::tmp<Foam::surfaceScalarField> Foam::fvc::absolute
+(
+    const tmp<surfaceScalarField>& tphi,
+    const volScalarField& rho,
+    const volVectorField& U
+)
+{
+    if (tphi().mesh().moving())
+    {
+        return tphi + fvc::meshPhi(rho, U);
+    }
+    else
+    {
+        return tmp<surfaceScalarField>(tphi, true);
+    }
+}
+
+
 // ************************************************************************* //

src/finiteVolume/finiteVolume/fvc/fvcMeshPhi.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -116,6 +116,22 @@ namespace fvc
         const volScalarField& rho,
         const volVectorField& U
     );
+
+
+    //- Return the given relative flux in absolute form
+    tmp<surfaceScalarField> absolute
+    (
+        const tmp<surfaceScalarField>& tphi,
+        const volVectorField& U
+    );
+
+    //- Return the given relative flux in absolute form
+    tmp<surfaceScalarField> absolute
+    (
+        const tmp<surfaceScalarField>& tphi,
+        const volScalarField& rho,
+        const volVectorField& U
+    );
 }
 
 

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -54,6 +54,7 @@ namespace Foam
         scalar
     )
 
+/*
     makeLimitedSurfaceInterpolationTypeScheme
     (
         MUSCL,
@@ -71,6 +72,7 @@ namespace Foam
         rhoMagSqr,
         vector
     )
+*/
 }
 
 // ************************************************************************* //

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

@@ -346,7 +346,7 @@ Foam::KinematicCloud<CloudType>::KinematicCloud
         (
             IOobject
             (
-                this->name() + "::UTrans",
+                this->name() + ".UTrans",
                 this->db().time().timeName(),
                 this->db(),
                 IOobject::READ_IF_PRESENT,
@@ -362,7 +362,7 @@ Foam::KinematicCloud<CloudType>::KinematicCloud
         (
             IOobject
             (
-                this->name() + "::UCoeff",
+                this->name() + ".UCoeff",
                 this->db().time().timeName(),
                 this->db(),
                 IOobject::READ_IF_PRESENT,
@@ -426,7 +426,7 @@ Foam::KinematicCloud<CloudType>::KinematicCloud
         (
             IOobject
             (
-                this->name() + "::UTrans",
+                this->name() + ".UTrans",
                 this->db().time().timeName(),
                 this->db(),
                 IOobject::NO_READ,
@@ -442,7 +442,7 @@ Foam::KinematicCloud<CloudType>::KinematicCloud
         (
             IOobject
             (
-                name + "::UCoeff",
+                name + ".UCoeff",
                 this->db().time().timeName(),
                 this->db(),
                 IOobject::NO_READ,

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

@@ -564,7 +564,7 @@ Foam::KinematicCloud<CloudType>::theta() const
         (
             IOobject
             (
-                this->name() + "::theta",
+                this->name() + ".theta",
                 this->db().time().timeName(),
                 this->db(),
                 IOobject::NO_READ,
@@ -603,7 +603,7 @@ Foam::KinematicCloud<CloudType>::alpha() const
         (
             IOobject
             (
-                this->name() + "::alpha",
+                this->name() + ".alpha",
                 this->db().time().timeName(),
                 this->db(),
                 IOobject::NO_READ,
@@ -640,7 +640,7 @@ Foam::KinematicCloud<CloudType>::rhoEff() const
         (
             IOobject
             (
-                this->name() + "::rhoEff",
+                this->name() + ".rhoEff",
                 this->db().time().timeName(),
                 this->db(),
                 IOobject::NO_READ,

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

@@ -107,7 +107,7 @@ inline Foam::tmp<Foam::fvScalarMatrix> Foam::ReactingCloud<CloudType>::SYi
                 (
                     IOobject
                     (
-                        this->name() + "::rhoTrans",
+                        this->name() + ".rhoTrans",
                         this->db().time().timeName(),
                         this->db(),
                         IOobject::NO_READ,
@@ -155,7 +155,7 @@ Foam::ReactingCloud<CloudType>::Srho(const label i) const
         (
             IOobject
             (
-                this->name() + "::rhoTrans",
+                this->name() + ".rhoTrans",
                 this->db().time().timeName(),
                 this->db(),
                 IOobject::NO_READ,
@@ -192,7 +192,7 @@ Foam::ReactingCloud<CloudType>::Srho() const
         (
             IOobject
             (
-                this->name() + "::rhoTrans",
+                this->name() + ".rhoTrans",
                 this->db().time().timeName(),
                 this->db(),
                 IOobject::NO_READ,
@@ -236,7 +236,7 @@ Foam::ReactingCloud<CloudType>::Srho(volScalarField& rho) const
             (
                 IOobject
                 (
-                    this->name() + "::rhoTrans",
+                    this->name() + ".rhoTrans",
                     this->db().time().timeName(),
                     this->db(),
                     IOobject::NO_READ,

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

@@ -64,7 +64,7 @@ void Foam::ThermoCloud<CloudType>::setModels()
             (
                 IOobject
                 (
-                    this->name() + "::radAreaP",
+                    this->name() + ".radAreaP",
                     this->db().time().timeName(),
                     this->db(),
                     IOobject::READ_IF_PRESENT,
@@ -81,7 +81,7 @@ void Foam::ThermoCloud<CloudType>::setModels()
             (
                 IOobject
                 (
-                    this->name() + "::radT4",
+                    this->name() + ".radT4",
                     this->db().time().timeName(),
                     this->db(),
                     IOobject::READ_IF_PRESENT,
@@ -98,7 +98,7 @@ void Foam::ThermoCloud<CloudType>::setModels()
             (
                 IOobject
                 (
-                    this->name() + "::radAreaPT4",
+                    this->name() + ".radAreaPT4",
                     this->db().time().timeName(),
                     this->db(),
                     IOobject::READ_IF_PRESENT,
@@ -169,7 +169,7 @@ Foam::ThermoCloud<CloudType>::ThermoCloud
         (
             IOobject
             (
-                this->name() + "::hsTrans",
+                this->name() + ".hsTrans",
                 this->db().time().timeName(),
                 this->db(),
                 IOobject::READ_IF_PRESENT,
@@ -185,7 +185,7 @@ Foam::ThermoCloud<CloudType>::ThermoCloud
         (
             IOobject
             (
-                this->name() + "::hsCoeff",
+                this->name() + ".hsCoeff",
                 this->db().time().timeName(),
                 this->db(),
                 IOobject::READ_IF_PRESENT,
@@ -239,7 +239,7 @@ Foam::ThermoCloud<CloudType>::ThermoCloud
         (
             IOobject
             (
-                this->name() + "::hsTrans",
+                this->name() + ".hsTrans",
                 this->db().time().timeName(),
                 this->db(),
                 IOobject::NO_READ,
@@ -255,7 +255,7 @@ Foam::ThermoCloud<CloudType>::ThermoCloud
         (
             IOobject
             (
-                this->name() + "::hsCoeff",
+                this->name() + ".hsCoeff",
                 this->db().time().timeName(),
                 this->db(),
                 IOobject::NO_READ,
@@ -274,7 +274,7 @@ Foam::ThermoCloud<CloudType>::ThermoCloud
             (
                 IOobject
                 (
-                    this->name() + "::radAreaP",
+                    this->name() + ".radAreaP",
                     this->db().time().timeName(),
                     this->db(),
                     IOobject::NO_READ,
@@ -291,7 +291,7 @@ Foam::ThermoCloud<CloudType>::ThermoCloud
             (
                 IOobject
                 (
-                    this->name() + "::radT4",
+                    this->name() + ".radT4",
                     this->db().time().timeName(),
                     this->db(),
                     IOobject::NO_READ,
@@ -308,7 +308,7 @@ Foam::ThermoCloud<CloudType>::ThermoCloud
             (
                 IOobject
                 (
-                    this->name() + "::radAreaPT4",
+                    this->name() + ".radAreaPT4",
                     this->db().time().timeName(),
                     this->db(),
                     IOobject::NO_READ,

src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/LocalInteraction/LocalInteraction.C

@@ -46,8 +46,8 @@ Foam::LocalInteraction<CloudType>::LocalInteraction
 {
     if (writeFields_)
     {
-        word massEscapeName(this->owner().name() + "::massEscape");
-        word massStickName(this->owner().name() + "::massStick");
+        word massEscapeName(this->owner().name() + ".massEscape");
+        word massStickName(this->owner().name() + ".massStick");
         Info<< "    Interaction fields will be written to " << massEscapeName
             << " and " << massStickName << endl;
 
@@ -121,7 +121,7 @@ Foam::volScalarField& Foam::LocalInteraction<CloudType>::massEscape()
             (
                 IOobject
                 (
-                    this->owner().name() + "::massEscape",
+                    this->owner().name() + ".massEscape",
                     mesh.time().timeName(),
                     mesh,
                     IOobject::READ_IF_PRESENT,
@@ -150,7 +150,7 @@ Foam::volScalarField& Foam::LocalInteraction<CloudType>::massStick()
             (
                 IOobject
                 (
-                    this->owner().name() + "::massStick",
+                    this->owner().name() + ".massStick",
                     mesh.time().timeName(),
                     mesh,
                     IOobject::READ_IF_PRESENT,

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

@@ -710,7 +710,8 @@ void Foam::autoLayerDriver::setNumLayers
     const indirectPrimitivePatch& pp,
     pointField& patchDisp,
     labelList& patchNLayers,
-    List<extrudeMode>& extrudeStatus
+    List<extrudeMode>& extrudeStatus,
+    label& nAddedCells
 ) const
 {
     const fvMesh& mesh = meshRefiner_.mesh();
@@ -797,6 +798,24 @@ void Foam::autoLayerDriver::setNumLayers
         }
     }
 
+
+    // Calculate number of cells to create
+    nAddedCells = 0;
+    forAll(pp.localFaces(), faceI)
+    {
+        const face& f = pp.localFaces()[faceI];
+
+        // Get max of extrusion per point
+        label nCells = 0;
+        forAll(f, fp)
+        {
+            nCells = max(nCells, patchNLayers[f[fp]]);
+        }
+
+        nAddedCells += nCells;
+    }
+    reduce(nAddedCells, sumOp<label>());
+
     //reduce(nConflicts, sumOp<label>());
     //
     //Info<< "Set displacement to zero for " << nConflicts
@@ -2491,9 +2510,13 @@ void Foam::autoLayerDriver::addLayers
     // extrudeStatus = EXTRUDE.
     labelList patchNLayers(pp().nPoints(), 0);
 
+    // Ideal number of cells added
+    label nIdealAddedCells = 0;
+
     // Whether to add edge for all pp.localPoints.
     List<extrudeMode> extrudeStatus(pp().nPoints(), EXTRUDE);
 
+
     {
         // Get number of layer per point from number of layers per patch
         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -2506,7 +2529,8 @@ void Foam::autoLayerDriver::addLayers
 
             patchDisp,
             patchNLayers,
-            extrudeStatus
+            extrudeStatus,
+            nIdealAddedCells
         );
 
         // Precalculate mesh edge labels for patch edges
@@ -2776,6 +2800,7 @@ void Foam::autoLayerDriver::addLayers
         layerParams.nSmoothNormals(),
         layerParams.nSmoothSurfaceNormals(),
         layerParams.minMedianAxisAngleCos(),
+        layerParams.slipFeatureAngle(),
 
         dispVec,
         medialRatio,
@@ -3101,10 +3126,24 @@ void Foam::autoLayerDriver::addLayers
 
         label nExtruded = countExtrusion(pp, extrudeStatus);
         label nTotFaces = returnReduce(pp().size(), sumOp<label>());
+        label nAddedCells = 0;
+        {
+            forAll(flaggedCells, cellI)
+            {
+                if (flaggedCells[cellI])
+                {
+                    nAddedCells++;
+                }
+            }
+            reduce(nAddedCells, sumOp<label>());
+        }
         Info<< "Extruding " << nExtruded
             << " out of " << nTotFaces
             << " faces (" << 100.0*nExtruded/nTotFaces << "%)."
             << " Removed extrusion at " << nTotChanged << " faces."
+            << endl
+            << "Added " << nAddedCells << " out of " << nIdealAddedCells
+            << " cells (" << 100.0*nAddedCells/nIdealAddedCells << "%)."
             << endl;
 
         if (nTotChanged == 0)

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

@@ -192,7 +192,8 @@ class autoLayerDriver
                     const indirectPrimitivePatch& pp,
                     pointField& patchDisp,
                     labelList& patchNLayers,
-                    List<extrudeMode>& extrudeStatus
+                    List<extrudeMode>& extrudeStatus,
+                    label& nIdealAddedCells
                 ) const;
 
                 //- Helper function to make a pointVectorField with correct
@@ -469,6 +470,7 @@ class autoLayerDriver
                     const label nSmoothNormals,
                     const label nSmoothSurfaceNormals,
                     const scalar minMedianAxisAngleCos,
+                    const scalar featureAngle,
 
                     pointVectorField& dispVec,
                     pointScalarField& medialRatio,

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

@@ -33,9 +33,10 @@ Description
 #include "motionSmoother.H"
 #include "pointData.H"
 #include "PointEdgeWave.H"
-#include "OFstream.H"
+#include "OBJstream.H"
 #include "meshTools.H"
 #include "PatchTools.H"
+#include "unitConversion.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -775,6 +776,7 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo
     const label nSmoothNormals,
     const label nSmoothSurfaceNormals,
     const scalar minMedianAxisAngleCos,
+    const scalar featureAngle,
 
     pointVectorField& dispVec,
     pointScalarField& medialRatio,
@@ -904,7 +906,7 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo
                 const edge& e = edges[edgeI];
                 // Approximate medial axis location on edge.
                 //const point medialAxisPt = e.centre(points);
-                vector eVec = e.vec(mesh.points());
+                vector eVec = e.vec(points);
                 scalar eMag = mag(eVec);
                 if (eMag > VSMALL)
                 {
@@ -961,42 +963,111 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo
 
         labelHashSet adaptPatches(meshMover.adaptPatchIDs());
 
+
         forAll(patches, patchI)
         {
             const polyPatch& pp = patches[patchI];
-            //const pointPatchVectorField& pvf =
-            //    meshMover.displacement().boundaryField()[patchI];
+            const pointPatchVectorField& pvf =
+                meshMover.displacement().boundaryField()[patchI];
 
             if
             (
                 !pp.coupled()
              && !isA<emptyPolyPatch>(pp)
-             //&&  pvf.constraintType() != word::null //exclude fixedValue
              && !adaptPatches.found(patchI)
             )
             {
-                Info<< "Inserting points on patch " << pp.name() << endl;
-
                 const labelList& meshPoints = pp.meshPoints();
 
-                forAll(meshPoints, i)
+                // Check the type of the patchField. The types are
+                //  - fixedValue (0 or more layers) but the >0 layers have
+                //    already been handled in the adaptPatches loop
+                //  - constraint (but not coupled) types, e.g. symmetryPlane,
+                //    slip.
+                if (pvf.fixesValue())
                 {
-                    label pointI = meshPoints[i];
+                    // Disable all movement on fixedValue patchFields
+                    Info<< "Inserting all points on patch " << pp.name()
+                        << endl;
 
-                    if (!pointMedialDist[pointI].valid(dummyTrackData))
+                    forAll(meshPoints, i)
                     {
-                        maxPoints.append(pointI);
-                        maxInfo.append
-                        (
-                            pointData
+                        label pointI = meshPoints[i];
+                        if (!pointMedialDist[pointI].valid(dummyTrackData))
+                        {
+                            maxPoints.append(pointI);
+                            maxInfo.append
                             (
-                                points[pointI],
-                                0.0,
-                                pointI,         // passive data
-                                vector::zero    // passive data
-                            )
-                        );
-                        pointMedialDist[pointI] = maxInfo.last();
+                                pointData
+                                (
+                                    points[pointI],
+                                    0.0,
+                                    pointI,         // passive data
+                                    vector::zero    // passive data
+                                )
+                            );
+                            pointMedialDist[pointI] = maxInfo.last();
+                        }
+                    }
+                }
+                else
+                {
+                    // Based on geometry: analyse angle w.r.t. nearest moving
+                    // point. In the pointWallDist we transported the
+                    // normal as the passive vector. Note that this points
+                    // out of the originating wall so inside of the domain
+                    // on this patch.
+                    Info<< "Inserting points on patch " << pp.name()
+                        << " if angle to nearest layer patch > "
+                        << featureAngle << " degrees." << endl;
+
+                    scalar featureAngleCos = Foam::cos(degToRad(featureAngle));
+                    pointField pointNormals
+                    (
+                        PatchTools::pointNormals
+                        (
+                            mesh,
+                            pp,
+                            identity(pp.size())+pp.start()
+                        )
+                    );
+
+                    forAll(meshPoints, i)
+                    {
+                        label pointI = meshPoints[i];
+
+                        if (!pointMedialDist[pointI].valid(dummyTrackData))
+                        {
+                            // Check if angle not too large.
+                            scalar cosAngle =
+                            (
+                               -pointWallDist[pointI].v()
+                              & pointNormals[i]
+                            );
+                            if (cosAngle > featureAngleCos)
+                            {
+                                // Extrusion direction practically perpendicular
+                                // to the patch. Disable movement at the patch.
+
+                                maxPoints.append(pointI);
+                                maxInfo.append
+                                (
+                                    pointData
+                                    (
+                                        points[pointI],
+                                        0.0,
+                                        pointI,         // passive data
+                                        vector::zero    // passive data
+                                    )
+                                );
+                                pointMedialDist[pointI] = maxInfo.last();
+                            }
+                            else
+                            {
+                                // Extrusion direction makes angle with patch
+                                // so allow sliding - don't insert zero points
+                            }
+                        }
                     }
                 }
             }
@@ -1141,13 +1212,12 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance
     // reduce thickness where thickness/medial axis distance large
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    autoPtr<OFstream> str;
-    label vertI = 0;
+    autoPtr<OBJstream> str;
     if (debug)
     {
         str.reset
         (
-            new OFstream
+            new OBJstream
             (
                 mesh.time().path()
               / "thicknessRatioExcludePoints_"
@@ -1159,13 +1229,12 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance
             << str().name() << endl;
     }
 
-    autoPtr<OFstream> medialVecStr;
-    label medialVertI = 0;
+    autoPtr<OBJstream> medialVecStr;
     if (debug)
     {
         medialVecStr.reset
         (
-            new OFstream
+            new OBJstream
             (
                 mesh.time().path()
               / "thicknessRatioExcludeMedialVec_"
@@ -1227,21 +1296,19 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance
                 if (str.valid())
                 {
                     const point& pt = mesh.points()[pointI];
-                    meshTools::writeOBJ(str(), pt);
-                    vertI++;
-                    meshTools::writeOBJ(str(), pt+patchDisp[patchPointI]);
-                    vertI++;
-                    str()<< "l " << vertI-1 << ' ' << vertI << nl;
+                    str().write(linePointRef(pt, pt+patchDisp[patchPointI]));
                 }
                 if (medialVecStr.valid())
                 {
                     const point& pt = mesh.points()[pointI];
-                    meshTools::writeOBJ(medialVecStr(), pt);
-                    medialVertI++;
-                    meshTools::writeOBJ(medialVecStr(), medialVec[pointI]);
-                    medialVertI++;
-                    medialVecStr()<< "l " << medialVertI-1
-                        << ' ' << medialVertI << nl;
+                    medialVecStr().write
+                    (
+                        linePointRef
+                        (
+                            pt,
+                            medialVec[pointI]
+                        )
+                    );
                 }
             }
         }

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

@@ -178,6 +178,24 @@ class autoSnapDriver
                     vector& edgeOffset  // offset from pt to point on edge
                 ) const;
 
+
+                //- For any reverse (so from feature back to mesh) attraction:
+                //  add attraction if diagonal points on face attracted
+                void stringFeatureEdges
+                (
+                    const label iter,
+                    const scalar featureCos,
+
+                    const indirectPrimitivePatch& pp,
+                    const scalarField& snapDist,
+
+                    const vectorField& rawPatchAttraction,
+                    const List<pointConstraint>& rawPatchConstraints,
+
+                    vectorField& patchAttraction,
+                    List<pointConstraint>& patchConstraints
+                ) const;
+
                 //- Return hit if on multiple points
                 pointIndexHit findMultiPatchPoint
                 (

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

@@ -848,17 +848,17 @@ void Foam::autoSnapDriver::featureAttractionUsingReconstruction
         vector d = r.refPoint()-pt;
         d -= (d&n)*n;
 
-        // Correct for attraction to non-dominant face
-        correctAttraction
-        (
-            surfacePoints,
-            surfaceCount,
-            r.refPoint(),
-            n,                  // normalised normal
-            pt,
-
-            d                   // perpendicular offset vector
-        );
+        //// Correct for attraction to non-dominant face
+        //correctAttraction
+        //(
+        //    surfacePoints,
+        //    surfaceCount,
+        //    r.refPoint(),
+        //    n,                  // normalised normal
+        //    pt,
+        //
+        //    d                   // perpendicular offset vector
+        //);
 
         // Trim to snap distance
         if (magSqr(d) > sqr(snapDist[pointI]))
@@ -893,6 +893,15 @@ void Foam::autoSnapDriver::featureAttractionUsingReconstruction
         patchConstraint.applyConstraint(surfaceNormals[0]);
         patchConstraint.applyConstraint(surfaceNormals[1]);
         patchConstraint.applyConstraint(surfaceNormals[2]);
+
+        //Pout<< "# Feature point " << pt << nl;
+        //meshTools::writeOBJ(Pout, pt);
+        //meshTools::writeOBJ(Pout, surfacePoints[0]);
+        //meshTools::writeOBJ(Pout, surfacePoints[1]);
+        //meshTools::writeOBJ(Pout, surfacePoints[2]);
+        //Pout<< "l 1 2" << nl
+        //    << "l 1 3" << nl
+        //    << "l 1 4" << nl;
     }
 }
 
@@ -1001,6 +1010,195 @@ void Foam::autoSnapDriver::featureAttractionUsingReconstruction
 }
 
 
+void Foam::autoSnapDriver::stringFeatureEdges
+(
+    const label iter,
+    const scalar featureCos,
+
+    const indirectPrimitivePatch& pp,
+    const scalarField& snapDist,
+
+    const vectorField& rawPatchAttraction,
+    const List<pointConstraint>& rawPatchConstraints,
+
+    vectorField& patchAttraction,
+    List<pointConstraint>& patchConstraints
+) const
+{
+    // Snap edges to feature edges
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Walk existing edges and snap remaining ones (that are marked as
+    // feature edges in rawPatchConstraints)
+
+    // What this does is fill in any faces where not all points
+    // on the face are being attracted:
+    /*
+           +
+          / \
+         /   \
+      ---+    +---
+         \   /
+          \ /
+           +
+    */
+    // so the top and bottom will never get attracted since the nearest
+    // back from the feature edge will always be one of the left or right
+    // points since the face is diamond like. So here we walk the feature edges
+    // and add any non-attracted points.
+
+
+    while (true)
+    {
+        label nChanged = 0;
+
+        const labelListList& pointEdges = pp.pointEdges();
+        forAll(pointEdges, pointI)
+        {
+            if (patchConstraints[pointI].first() == 2)
+            {
+                const point& pt = pp.localPoints()[pointI];
+                const labelList& pEdges = pointEdges[pointI];
+                const vector& featVec = patchConstraints[pointI].second();
+
+                // Detect whether there are edges in both directions.
+                // (direction along the feature edge that is)
+                bool hasPos = false;
+                bool hasNeg = false;
+
+                forAll(pEdges, pEdgeI)
+                {
+                    const edge& e = pp.edges()[pEdges[pEdgeI]];
+                    label nbrPointI = e.otherVertex(pointI);
+
+                    if (patchConstraints[nbrPointI].first() > 1)
+                    {
+                        const point& nbrPt = pp.localPoints()[nbrPointI];
+                        const point featPt =
+                            nbrPt + patchAttraction[nbrPointI];
+                        const scalar cosAngle = (featVec & (featPt-pt));
+
+                        if (cosAngle > 0)
+                        {
+                            hasPos = true;
+                        }
+                        else
+                        {
+                            hasNeg = true;
+                        }
+                    }
+                }
+
+                if (!hasPos || !hasNeg)
+                {
+                    //Pout<< "**Detected feature string end at  "
+                    //    << pp.localPoints()[pointI] << endl;
+
+                    // No string. Assign best choice on either side
+                    label bestPosPointI = -1;
+                    scalar minPosDistSqr = GREAT;
+                    label bestNegPointI = -1;
+                    scalar minNegDistSqr = GREAT;
+
+                    forAll(pEdges, pEdgeI)
+                    {
+                        const edge& e = pp.edges()[pEdges[pEdgeI]];
+                        label nbrPointI = e.otherVertex(pointI);
+
+                        if
+                        (
+                            patchConstraints[nbrPointI].first() <= 1
+                         && rawPatchConstraints[nbrPointI].first() > 1
+                        )
+                        {
+                            const vector& nbrFeatVec =
+                                rawPatchConstraints[pointI].second();
+
+                            if (mag(featVec&nbrFeatVec) > featureCos)
+                            {
+                                // nbrPointI attracted to sameish feature
+                                // Note: also check on position.
+
+                                scalar d2 = magSqr
+                                (
+                                    rawPatchAttraction[nbrPointI]
+                                );
+
+                                const point featPt =
+                                    pp.localPoints()[nbrPointI]
+                                  + rawPatchAttraction[nbrPointI];
+                                const scalar cosAngle =
+                                    (featVec & (featPt-pt));
+
+                                if (cosAngle > 0)
+                                {
+                                    if (!hasPos && d2 < minPosDistSqr)
+                                    {
+                                        minPosDistSqr = d2;
+                                        bestPosPointI = nbrPointI;
+                                    }
+                                }
+                                else
+                                {
+                                    if (!hasNeg && d2 < minNegDistSqr)
+                                    {
+                                        minNegDistSqr = d2;
+                                        bestNegPointI = nbrPointI;
+                                    }
+                                }
+                            }
+                        }
+                    }
+
+                    if (bestPosPointI != -1)
+                    {
+                        // Use reconstructed-feature attraction. Use only
+                        // part of it since not sure...
+                        //const point& bestPt =
+                        //    pp.localPoints()[bestPosPointI];
+                        //Pout<< "**Overriding point " << bestPt
+                        //    << " on reconstructed feature edge at "
+                        //    << rawPatchAttraction[bestPosPointI]+bestPt
+                        //    << " to attracted-to-feature-edge." << endl;
+                        patchAttraction[bestPosPointI] =
+                            0.5*rawPatchAttraction[bestPosPointI];
+                        patchConstraints[bestPosPointI] =
+                            rawPatchConstraints[bestPosPointI];
+
+                        nChanged++;
+                    }
+                    if (bestNegPointI != -1)
+                    {
+                        // Use reconstructed-feature attraction. Use only
+                        // part of it since not sure...
+                        //const point& bestPt =
+                        //    pp.localPoints()[bestNegPointI];
+                        //Pout<< "**Overriding point " << bestPt
+                        //    << " on reconstructed feature edge at "
+                        //    << rawPatchAttraction[bestNegPointI]+bestPt
+                        //    << " to attracted-to-feature-edge." << endl;
+                        patchAttraction[bestNegPointI] =
+                            0.5*rawPatchAttraction[bestNegPointI];
+                        patchConstraints[bestNegPointI] =
+                            rawPatchConstraints[bestNegPointI];
+
+                        nChanged++;
+                    }
+                }
+            }
+        }
+
+
+        reduce(nChanged, sumOp<label>());
+        Info<< "Stringing feature edges : changed " << nChanged << " points"
+            << endl;
+        if (nChanged == 0)
+        {
+            break;
+        }
+    }
+}
+
+
 Foam::pointIndexHit Foam::autoSnapDriver::findNearFeatureEdge
 (
     const indirectPrimitivePatch& pp,
@@ -1125,7 +1323,7 @@ Foam::labelPair Foam::autoSnapDriver::findNearFeaturePoint
                 patchConstraints[oldPointI] = pointConstraint();
 
                 label edgeFeatI;
-                const pointIndexHit nearInfo = findNearFeatureEdge
+                findNearFeatureEdge
                 (
                     pp,
                     snapDist,
@@ -1721,7 +1919,7 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
                             // the old point to attract to nearest edge
                             // instead.
                             label edgeFeatI;
-                            const pointIndexHit nearInfo = findNearFeatureEdge
+                            findNearFeatureEdge
                             (
                                 pp,
                                 snapDist,
@@ -1752,7 +1950,7 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
                     //    << endl;
 
                     label featI;
-                    const pointIndexHit nearInfo = findNearFeatureEdge
+                    findNearFeatureEdge
                     (
                         pp,
                         snapDist,
@@ -2113,156 +2311,20 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
     // Walk existing edges and snap remaining ones (that are marked as
     // feature edges in allPatchConstraints)
 
-    while (true)
-    {
-        label nChanged = 0;
+    stringFeatureEdges
+    (
+        iter,
+        featureCos,
 
-        const labelListList& pointEdges = pp.pointEdges();
-        forAll(pointEdges, pointI)
-        {
-            if (patchConstraints[pointI].first() == 2)
-            {
-                const point& pt = pp.localPoints()[pointI];
-                const labelList& pEdges = pointEdges[pointI];
-                const vector& featVec = patchConstraints[pointI].second();
+        pp,
+        snapDist,
 
-                // Detect whether there are edges in both directions.
-                // (direction along the feature edge that is)
-                bool hasPos = false;
-                bool hasNeg = false;
-
-                forAll(pEdges, pEdgeI)
-                {
-                    const edge& e = pp.edges()[pEdges[pEdgeI]];
-                    label nbrPointI = e.otherVertex(pointI);
-
-                    if (patchConstraints[nbrPointI].first() > 1)
-                    {
-                        const point& nbrPt = pp.localPoints()[nbrPointI];
-                        const point featPt =
-                            nbrPt + patchAttraction[nbrPointI];
-                        const scalar cosAngle = (featVec & (featPt-pt));
-
-                        if (cosAngle > 0)
-                        {
-                            hasPos = true;
-                        }
-                        else
-                        {
-                            hasNeg = true;
-                        }
-                    }
-                }
-
-                if (!hasPos || !hasNeg)
-                {
-                    //Pout<< "**Detected feature string end at  "
-                    //    << pp.localPoints()[pointI] << endl;
-
-                    // No string. Assign best choice on either side
-                    label bestPosPointI = -1;
-                    scalar minPosDistSqr = GREAT;
-                    label bestNegPointI = -1;
-                    scalar minNegDistSqr = GREAT;
-
-                    forAll(pEdges, pEdgeI)
-                    {
-                        const edge& e = pp.edges()[pEdges[pEdgeI]];
-                        label nbrPointI = e.otherVertex(pointI);
-
-                        if
-                        (
-                            patchConstraints[nbrPointI].first() <= 1
-                         && allPatchConstraints[nbrPointI].first() > 1
-                        )
-                        {
-                            const vector& nbrFeatVec =
-                                allPatchConstraints[pointI].second();
-
-                            if (mag(featVec&nbrFeatVec) > featureCos)
-                            {
-                                // nbrPointI attracted to sameish feature
-                                // Note: also check on position.
-
-                                scalar d2 = magSqr
-                                (
-                                    allPatchAttraction[nbrPointI]
-                                );
-
-                                const point featPt =
-                                    pp.localPoints()[nbrPointI]
-                                  + allPatchAttraction[nbrPointI];
-                                const scalar cosAngle =
-                                    (featVec & (featPt-pt));
-
-                                if (cosAngle > 0)
-                                {
-                                    if (!hasPos && d2 < minPosDistSqr)
-                                    {
-                                        minPosDistSqr = d2;
-                                        bestPosPointI = nbrPointI;
-                                    }
-                                }
-                                else
-                                {
-                                    if (!hasNeg && d2 < minNegDistSqr)
-                                    {
-                                        minNegDistSqr = d2;
-                                        bestNegPointI = nbrPointI;
-                                    }
-                                }
-                            }
-                        }
-                    }
-
-                    if (bestPosPointI != -1)
-                    {
-                        // Use reconstructed-feature attraction. Use only
-                        // part of it since not sure...
-                        //const point& bestPt =
-                        //    pp.localPoints()[bestPosPointI];
-                        //Pout<< "**Overriding point " << bestPt
-                        //    << " on reconstructed feature edge at "
-                        //    << allPatchAttraction[bestPosPointI]+bestPt
-                        //    << " to attracted-to-feature-edge." << endl;
-                        patchAttraction[bestPosPointI] =
-                            0.5*allPatchAttraction[bestPosPointI];
-                        patchConstraints[bestPosPointI] =
-                            allPatchConstraints[bestPosPointI];
-
-                        nChanged++;
-                    }
-                    if (bestNegPointI != -1)
-                    {
-                        // Use reconstructed-feature attraction. Use only
-                        // part of it since not sure...
-                        //const point& bestPt =
-                        //    pp.localPoints()[bestNegPointI];
-                        //Pout<< "**Overriding point " << bestPt
-                        //    << " on reconstructed feature edge at "
-                        //    << allPatchAttraction[bestNegPointI]+bestPt
-                        //    << " to attracted-to-feature-edge." << endl;
-                        patchAttraction[bestNegPointI] =
-                            0.5*allPatchAttraction[bestNegPointI];
-                        patchConstraints[bestNegPointI] =
-                            allPatchConstraints[bestNegPointI];
-
-                        nChanged++;
-                    }
-                }
-            }
-        }
-
-
-        reduce(nChanged, sumOp<label>());
-        Info<< "Stringing feature edges : changed " << nChanged << " points"
-            << endl;
-        if (nChanged == 0)
-        {
-            break;
-        }
-    }
+        allPatchAttraction,
+        allPatchConstraints,
 
+        patchAttraction,
+        patchConstraints
+    );
 
 
     // Avoid diagonal attraction

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

@@ -61,6 +61,12 @@ Foam::layerParameters::layerParameters
         readScalar(dict.lookup("minThickness"))
     ),
     featureAngle_(readScalar(dict.lookup("featureAngle"))),
+    slipFeatureAngle_
+    (
+        dict.found("slipFeatureAngle")
+      ? readScalar(dict.lookup("slipFeatureAngle"))
+      : 0.5*featureAngle_
+    ),
     concaveAngle_
     (
         dict.lookupOrDefault("concaveAngle", defaultConcaveAngle)

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

@@ -79,6 +79,8 @@ class layerParameters
 
         scalar featureAngle_;
 
+        scalar slipFeatureAngle_;
+
         scalar concaveAngle_;
 
         label nGrow_;
@@ -178,6 +180,13 @@ public:
                 return featureAngle_;
             }
 
+            //- At non-patched sides allow mesh to slip if extrusion
+            //  direction makes angle larger than slipFeatureAngle.
+            scalar slipFeatureAngle() const
+            {
+                return slipFeatureAngle_;
+            }
+
             scalar concaveAngle() const
             {
                 return concaveAngle_;

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

@@ -312,7 +312,8 @@ Foam::label Foam::meshRefinement::markFeatureRefinement
                     {
                         Pout<< "Adding particle from point:" << pointI
                             << " coord:" << featureMesh.points()[pointI]
-                            << " pEdges:" << pointEdges[pointI]
+                            << " since number of emanating edges:"
+                            << pointEdges[pointI].size()
                             << endl;
                     }
 

src/meshTools/regionSplit/regionSplit.C

@@ -69,7 +69,7 @@ void Foam::regionSplit::transferCoupledFaceRegion
                   "regionSplit::transferCoupledFaceRegion"
                   "(const label, const label, labelList&, labelList&) const"
               )   << "Problem : coupled face " << faceI
-                  << " on patch " << mesh_.boundaryMesh().whichPatch(faceI)
+                  << " on patch " << mesh().boundaryMesh().whichPatch(faceI)
                   << " has region " << faceRegion[faceI]
                   << " but coupled face " << otherFaceI
                   << " has region " << faceRegion[otherFaceI]
@@ -109,7 +109,7 @@ void Foam::regionSplit::fillSeedMask
 
 
     // Collect faces on seed cell
-    const cell& cFaces = mesh_.cells()[seedCellID];
+    const cell& cFaces = mesh().cells()[seedCellID];
 
     label nFaces = 0;
 
@@ -144,7 +144,7 @@ void Foam::regionSplit::fillSeedMask
         {
             label faceI = changedFaces[i];
 
-            label own = mesh_.faceOwner()[faceI];
+            label own = mesh().faceOwner()[faceI];
 
             if (cellRegion[own] == -1)
             {
@@ -152,9 +152,9 @@ void Foam::regionSplit::fillSeedMask
                 changedCells.append(own);
             }
 
-            if (mesh_.isInternalFace(faceI))
+            if (mesh().isInternalFace(faceI))
             {
-                label nei = mesh_.faceNeighbour()[faceI];
+                label nei = mesh().faceNeighbour()[faceI];
 
                 if (cellRegion[nei] == -1)
                 {
@@ -178,7 +178,7 @@ void Foam::regionSplit::fillSeedMask
         {
             label cellI = changedCells[i];
 
-            const cell& cFaces = mesh_.cells()[cellI];
+            const cell& cFaces = mesh().cells()[cellI];
 
             forAll(cFaces, cFaceI)
             {
@@ -203,7 +203,7 @@ void Foam::regionSplit::fillSeedMask
         // Check for changes to any locally coupled face.
         // Global connections are done later.
 
-        const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+        const polyBoundaryMesh& patches = mesh().boundaryMesh();
 
         forAll(patches, patchI)
         {
@@ -274,7 +274,7 @@ Foam::label Foam::regionSplit::calcLocalRegionSplit
         {
             // Check that blockedFace is synced.
             boolList syncBlockedFace(blockedFace);
-            syncTools::swapFaceList(mesh_, syncBlockedFace);
+            syncTools::swapFaceList(mesh(), syncBlockedFace);
 
             forAll(syncBlockedFace, faceI)
             {
@@ -295,7 +295,7 @@ Foam::label Foam::regionSplit::calcLocalRegionSplit
     // Region per face.
     // -1 unassigned
     // -2 blocked
-    labelList faceRegion(mesh_.nFaces(), -1);
+    labelList faceRegion(mesh().nFaces(), -1);
 
     if (blockedFace.size())
     {
@@ -321,7 +321,7 @@ Foam::label Foam::regionSplit::calcLocalRegionSplit
     {
         // Find first unset cell
 
-        for (; unsetCellI < mesh_.nCells(); unsetCellI++)
+        for (; unsetCellI < mesh().nCells(); unsetCellI++)
         {
             if (cellRegion[unsetCellI] == -1)
             {
@@ -329,7 +329,7 @@ Foam::label Foam::regionSplit::calcLocalRegionSplit
             }
         }
 
-        if (unsetCellI >= mesh_.nCells())
+        if (unsetCellI >= mesh().nCells())
         {
             break;
         }
@@ -430,9 +430,9 @@ Foam::autoPtr<Foam::globalIndex> Foam::regionSplit::calcRegionSplit
         }
 
 
-        const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+        const polyBoundaryMesh& patches = mesh().boundaryMesh();
 
-        labelList nbrRegion(mesh_.nFaces()-mesh_.nInternalFaces(), -1);
+        labelList nbrRegion(mesh().nFaces()-mesh().nInternalFaces(), -1);
         forAll(patches, patchI)
         {
             const polyPatch& pp = patches[patchI];
@@ -444,7 +444,7 @@ Foam::autoPtr<Foam::globalIndex> Foam::regionSplit::calcRegionSplit
                 (
                     nbrRegion,
                     pp.size(),
-                    pp.start()-mesh_.nInternalFaces()
+                    pp.start()-mesh().nInternalFaces()
                 );
 
                 forAll(patchCells, i)
@@ -457,9 +457,9 @@ Foam::autoPtr<Foam::globalIndex> Foam::regionSplit::calcRegionSplit
                 }
             }
         }
-        syncTools::swapBoundaryFaceList(mesh_, nbrRegion);
+        syncTools::swapBoundaryFaceList(mesh(), nbrRegion);
 
-        Map<label> globalToMerged(mesh_.nFaces()-mesh_.nInternalFaces());
+        Map<label> globalToMerged(mesh().nFaces()-mesh().nInternalFaces());
 
         forAll(patches, patchI)
         {
@@ -472,7 +472,7 @@ Foam::autoPtr<Foam::globalIndex> Foam::regionSplit::calcRegionSplit
                 (
                     nbrRegion,
                     pp.size(),
-                    pp.start()-mesh_.nInternalFaces()
+                    pp.start()-mesh().nInternalFaces()
                 );
 
                 forAll(patchCells, i)
@@ -486,7 +486,7 @@ Foam::autoPtr<Foam::globalIndex> Foam::regionSplit::calcRegionSplit
                             //Pout<< "on patch:" << pp.name()
                             //    << " cell:" << patchCells[i]
                             //    << " at:"
-                            // << mesh_.cellCentres()[patchCells[i]]
+                            // << mesh().cellCentres()[patchCells[i]]
                             //    << " was:" << cellRegion[patchCells[i]]
                             //    << " nbr:" << patchNbrRegion[i]
                             //    << endl;
@@ -534,7 +534,7 @@ Foam::autoPtr<Foam::globalIndex> Foam::regionSplit::calcRegionSplit
     // 4a: count. Use a labelHashSet to count regions only once.
     label nCompact = 0;
     {
-        labelHashSet localRegion(mesh_.nFaces()-mesh_.nInternalFaces());
+        labelHashSet localRegion(mesh().nFaces()-mesh().nInternalFaces());
         forAll(cellRegion, cellI)
         {
             if
@@ -692,8 +692,8 @@ Foam::autoPtr<Foam::globalIndex> Foam::regionSplit::calcRegionSplit
 
 Foam::regionSplit::regionSplit(const polyMesh& mesh)
 :
-    labelList(mesh.nCells(), -1),
-    mesh_(mesh)
+    MeshObject<polyMesh, regionSplit>(mesh),
+    labelList(mesh.nCells(), -1)
 {
     globalNumberingPtr_ = calcRegionSplit
     (
@@ -710,8 +710,8 @@ Foam::regionSplit::regionSplit
     const boolList& blockedFace
 )
 :
-    labelList(mesh.nCells(), -1),
-    mesh_(mesh)
+    MeshObject<polyMesh, regionSplit>(mesh),
+    labelList(mesh.nCells(), -1)
 {
     globalNumberingPtr_ = calcRegionSplit
     (
@@ -729,8 +729,8 @@ Foam::regionSplit::regionSplit
     const List<labelPair>& explicitConnections
 )
 :
-    labelList(mesh.nCells(), -1),
-    mesh_(mesh)
+    MeshObject<polyMesh, regionSplit>(mesh),
+    labelList(mesh.nCells(), -1)
 {
     globalNumberingPtr_ = calcRegionSplit
     (

src/meshTools/regionSplit/regionSplit.H

@@ -98,6 +98,7 @@ SourceFiles
 #include "globalIndex.H"
 #include "labelPair.H"
 #include "boolList.H"
+#include "MeshObject.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -112,13 +113,11 @@ class polyMesh;
 
 class regionSplit
 :
+    public MeshObject<polyMesh, regionSplit>,
     public labelList
 {
     // Private data
 
-        //- Reference to mesh
-        const polyMesh& mesh_;
-
         mutable autoPtr<globalIndex> globalNumberingPtr_;
 
 
@@ -161,11 +160,13 @@ class regionSplit
             labelList& cellRegion
         ) const;
 
+
 public:
 
     //- Runtime type information
     ClassName("regionSplit");
 
+
     // Constructors
 
         //- Construct from mesh
@@ -185,6 +186,7 @@ public:
             const List<labelPair>&
         );
 
+
     // Member Functions
 
         //- Return global region numbering
@@ -198,7 +200,6 @@ public:
         {
             return globalNumbering().size();
         }
-
 };
 
 

src/postProcessing/functionObjects/field/turbulenceFields/turbulenceFields.C

@@ -141,7 +141,7 @@ void Foam::turbulenceFields::read(const dictionary& dict)
             Info<< "storing fields:" << nl;
             forAllConstIter(wordHashSet, fieldSet_, iter)
             {
-                Info<< "    " << modelName << "::" << iter.key() << nl;
+                Info<< "    " << modelName << '.' << iter.key() << nl;
             }
             Info<< endl;
         }

src/postProcessing/functionObjects/field/turbulenceFields/turbulenceFieldsTemplates.C

@@ -36,7 +36,7 @@ void Foam::turbulenceFields::processField
 {
     typedef GeometricField<Type, fvPatchField, volMesh> FieldType;
 
-    const word scopedName = modelName + "::" + fieldName;
+    const word scopedName = modelName + '.' + fieldName;
 
     if (obr_.foundObject<FieldType>(scopedName))
     {

src/randomProcesses/Kmesh/Kmesh.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -30,34 +30,27 @@ License
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-namespace Foam
+inline Foam::label Foam::Kmesh::index
+(
+    const label i,
+    const label j,
+    const label k,
+    const labelList& nn
+)
 {
-   //! \cond fileScope
-   inline label rep
-   (
-       const label i,
-       const label j,
-       const label k,
-       const labelList& nn
-   )
-   {
-       return (k + j*nn[2] + i*nn[1]*nn[2]);
-   }
-   //! \endcond
-
-} // End namespace Foam
+    return (k + j*nn[2] + i*nn[1]*nn[2]);
+}
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// from fvMesh
 Foam::Kmesh::Kmesh(const fvMesh& mesh)
 :
     vectorField(mesh.V().size()),
-    NN(vector::dim)
+    nn_(vector::dim)
 {
     boundBox box = mesh.bounds();
-    L = box.span();
+    l_ = box.span();
 
     vector cornerCellCentre = ::Foam::max(mesh.C().internalField());
     vector cellL = 2*(box.max() - cornerCellCentre);
@@ -65,18 +58,17 @@ Foam::Kmesh::Kmesh(const fvMesh& mesh)
     vector rdeltaByL;
     label nTot = 1;
 
-    label i;
-    forAll(NN, i)
+    forAll(nn_, i)
     {
-        NN[i] = label(L[i]/cellL[i] + 0.5);
-        nTot *= NN[i];
+        nn_[i] = label(l_[i]/cellL[i] + 0.5);
+        nTot *= nn_[i];
 
-        if (NN[i] > 1)
+        if (nn_[i] > 1)
         {
-            L[i] -= cellL[i];
+            l_[i] -= cellL[i];
         }
 
-        rdeltaByL[i] = NN[i]/(L[i]*L[i]);
+        rdeltaByL[i] = nn_[i]/(l_[i]*l_[i]);
     }
 
     if (nTot != mesh.nCells())
@@ -86,24 +78,31 @@ Foam::Kmesh::Kmesh(const fvMesh& mesh)
             << abort(FatalError);
     }
 
-    for (i=0; i<NN[0]; i++)
+    for (label i=0; i<nn_[0]; i++)
     {
-        scalar k1 = (i - NN[0]/2)*constant::mathematical::twoPi/L[0];
+        scalar k1 = (i - nn_[0]/2)*constant::mathematical::twoPi/l_[0];
 
-        for (label j=0; j<NN[1]; j++)
+        for (label j=0; j<nn_[1]; j++)
         {
-            scalar k2 = (j - NN[1]/2)*constant::mathematical::twoPi/L[1];
+            scalar k2 = (j - nn_[1]/2)*constant::mathematical::twoPi/l_[1];
 
-            for (label k=0; k<NN[2]; k++)
+            for (label k=0; k<nn_[2]; k++)
             {
-                scalar k3 = (k - NN[2]/2)*constant::mathematical::twoPi/L[2];
+                scalar k3 = (k - nn_[2]/2)*constant::mathematical::twoPi/l_[2];
 
-                (*this)[rep(i, j, k, NN)] = vector(k1, k2, k3);
+                (*this)[index(i, j, k, nn_)] = vector(k1, k2, k3);
             }
         }
     }
 
-    Kmax = mag((*this)[rep(NN[0]-1, NN[1]-1, NN[2]-1, NN)]);
+    kmax_ = mag
+    (
+        Foam::max
+        (
+            cmptMag((*this)[index(nn_[0]-1, nn_[1]-1, nn_[2]-1, nn_)]),
+            cmptMag((*this)[index(0, 0, 0, nn_)])
+        )
+    );
 }
 
 

src/randomProcesses/Kmesh/Kmesh.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -51,17 +51,29 @@ class Kmesh
 :
     public vectorField
 {
-
     // Private data
 
         //- Dimensions of box
-        vector L;
+        vector l_;
 
         //- Multi-dimensional addressing array
-        labelList NN;
+        labelList nn_;
 
         //- Maximum wavenumber
-        scalar Kmax;
+        scalar kmax_;
+
+
+    // Private functions
+
+        //- Return the linear index from the i-j-k indices
+        static inline label index
+        (
+            const label i,
+            const label j,
+            const label k,
+            const labelList& nn
+        );
+
 
 public:
 
@@ -77,17 +89,17 @@ public:
 
         const vector& sizeOfBox() const
         {
-            return L;
+            return l_;
         }
 
         const labelList& nn() const
         {
-            return NN;
+            return nn_;
         }
 
         scalar max() const
         {
-            return Kmax;
+            return kmax_;
         }
 };
 

src/regionModels/pyrolysisModels/Make/options

@@ -7,6 +7,7 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/solidThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/solidChemistryModel/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/radiationModels/lnInclude \
     -I$(LIB_SRC)/turbulenceModels \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/RAS/lnInclude \
@@ -26,4 +27,5 @@ LIB_LIBS = \
     -lcompressibleRASModels \
     -lcompressibleLESModels \
     -lLESdeltas \
-    -lregionModels
+    -lregionModels \
+    -lradiationModels

src/regionModels/pyrolysisModels/noPyrolysis/noPyrolysis.C

@@ -51,7 +51,8 @@ void noPyrolysis::constructThermoChemistry()
         solidChemistryModel::New(regionMesh()).ptr()
     );
 
-    solidThermo_.reset(&solidChemistry_->solid());
+    solidThermo_.reset(&solidChemistry_->solidThermo());
+    radiation_.reset(radiation::radiationModel::New(solidThermo_->T()).ptr());
 }
 
 bool noPyrolysis::read()
@@ -88,7 +89,8 @@ noPyrolysis::noPyrolysis(const word& modelType, const fvMesh& mesh)
 :
     pyrolysisModel(mesh),
     solidChemistry_(NULL),
-    solidThermo_(NULL)
+    solidThermo_(NULL),
+    radiation_(NULL)
 {
     if (active())
     {
@@ -105,7 +107,8 @@ noPyrolysis::noPyrolysis
 ):
     pyrolysisModel(mesh),
     solidChemistry_(NULL),
-    solidThermo_(NULL)
+    solidThermo_(NULL),
+    radiation_(NULL)
 {
     if (active())
     {
@@ -135,31 +138,31 @@ void noPyrolysis::evolveRegion()
 
 const volScalarField& noPyrolysis::rho() const
 {
-    return (solidThermo_->rho());
+    return solidThermo_->rho();
 }
 
 
 const volScalarField& noPyrolysis::T() const
 {
-    return (solidThermo_->T());
+    return solidThermo_->T();
 }
 
 
 const tmp<volScalarField> noPyrolysis::Cp() const
 {
-    return (solidThermo_->Cp());
+    return solidThermo_->Cp();
 }
 
 
 tmp<volScalarField> noPyrolysis::kappaRad() const
 {
-    return (solidThermo_->kappaRad());
+    return radiation_->absorptionEmission().a();
 }
 
 
 tmp<volScalarField> noPyrolysis::kappa() const
 {
-     return (solidThermo_->kappa());
+     return solidThermo_->kappa();
 }
 
 

src/regionModels/pyrolysisModels/noPyrolysis/noPyrolysis.H

@@ -38,6 +38,7 @@ SourceFiles
 #include "pyrolysisModel.H"
 #include "volFieldsFwd.H"
 #include "solidChemistryModel.H"
+#include "radiationModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -86,6 +87,9 @@ protected:
         //- Reference to solid thermo
         autoPtr<solidReactionThermo> solidThermo_;
 
+         //- Pointer to radiation model
+        autoPtr<radiation::radiationModel> radiation_;
+
 
 public:
 

src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C

@@ -363,7 +363,8 @@ reactingOneDim::reactingOneDim(const word& modelType, const fvMesh& mesh)
 :
     pyrolysisModel(modelType, mesh),
     solidChemistry_(solidChemistryModel::New(regionMesh())),
-    solidThermo_(solidChemistry_->solid()),
+    solidThermo_(solidChemistry_->solidThermo()),
+    radiation_(radiation::radiationModel::New(solidThermo_.T())),
     rho_
     (
         IOobject
@@ -474,7 +475,8 @@ reactingOneDim::reactingOneDim
 :
     pyrolysisModel(modelType, mesh, dict),
     solidChemistry_(solidChemistryModel::New(regionMesh())),
-    solidThermo_(solidChemistry_->solid()),
+    solidThermo_(solidChemistry_->solidThermo()),
+    radiation_(radiation::radiationModel::New(solidThermo_.T())),
     rho_
     (
         IOobject
@@ -656,7 +658,7 @@ const tmp<volScalarField> reactingOneDim::Cp() const
 
 tmp<volScalarField> reactingOneDim::kappaRad() const
 {
-    return solidThermo_.kappaRad();
+    return radiation_->absorptionEmission().a();
 }
 
 

src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.H

@@ -37,7 +37,7 @@ SourceFiles
 
 #include "pyrolysisModel.H"
 #include "solidChemistryModel.H"
-
+#include "radiationModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -81,6 +81,9 @@ protected:
         //- Reference to solid thermo
         solidReactionThermo& solidThermo_;
 
+        //- Pointer to radiation model
+        autoPtr<radiation::radiationModel> radiation_;
+
 
         // Reference to solid thermo properties
 
@@ -90,7 +93,7 @@ protected:
             //- List of solid components
             PtrList<volScalarField>& Ys_;
 
-            // Non-const access to temperature
+            // Non-const access to enthalpy
             volScalarField& h_;
 
 

src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/noRadiation/noRadiation.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -82,7 +82,7 @@ tmp<volScalarField> noRadiation::Shs()
         (
             IOobject
             (
-                typeName + "::Shs",
+                typeName + ".Shs",
                 owner().time().timeName(),
                 owner().regionMesh(),
                 IOobject::NO_READ,

src/regionModels/surfaceFilmModels/submodels/thermo/filmRadiationModel/standardRadiation/standardRadiation.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -115,7 +115,7 @@ tmp<volScalarField> standardRadiation::Shs()
         (
             IOobject
             (
-                typeName + "::Shs",
+                typeName + ".Shs",
                 owner().time().timeName(),
                 owner().regionMesh(),
                 IOobject::NO_READ,

src/regionModels/thermoBaffleModels/Make/options

@@ -8,6 +8,7 @@ EXE_INC = \
     -I$(LIB_SRC)/AMIInterpolation/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/radiationModels/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel/lnInclude
 
 LIB_LIBS = \
@@ -17,4 +18,5 @@ LIB_LIBS = \
     -lfiniteVolume \
     -lmeshTools \
     -lOpenFOAM \
-    -lsolidSpecie
+    -lsolidSpecie \
+    -lradiationModels

src/regionModels/thermoBaffleModels/derivedFvPatchFields/temperatureThermoBaffle/temperatureThermoBaffleFvPatchScalarField.C

@@ -46,7 +46,7 @@ temperatureThermoBaffleFvPatchScalarField
     turbulentTemperatureCoupledBaffleMixedFvPatchScalarField(p, iF),
     owner_(false),
     baffle_(),
-    solidThermoType_("undefined"),
+    solidThermoType_(new primitiveEntry("thermoType", "undefined")),
     dict_(dictionary::null)
 {}
 
@@ -69,7 +69,7 @@ temperatureThermoBaffleFvPatchScalarField
     ),
     owner_(ptf.owner_),
     baffle_(ptf.baffle_),
-    solidThermoType_(ptf.solidThermoType_),
+    solidThermoType_(ptf.solidThermoType_().clone()),
     dict_(ptf.dict_)
 {}
 
@@ -85,7 +85,7 @@ temperatureThermoBaffleFvPatchScalarField
     turbulentTemperatureCoupledBaffleMixedFvPatchScalarField(p, iF, dict),
     owner_(false),
     baffle_(),
-    solidThermoType_(),
+    solidThermoType_(new primitiveEntry("thermoType", "undefined")),
     dict_(dict)
 {
     if (!isA<mappedPatchBase>(patch().patch()))
@@ -126,7 +126,7 @@ temperatureThermoBaffleFvPatchScalarField
         Info << "Creating thermal baffle" <<  nbrMesh << endl;
         baffle_.reset(baffle::New(thisMesh, dict).ptr());
         owner_ = true;
-        dict.lookup("thermoType") >> solidThermoType_;
+        solidThermoType_ = dict.lookupEntry("thermoType", false, false).clone();
         baffle_->rename(nbrMesh);
     }
 }
@@ -142,7 +142,7 @@ temperatureThermoBaffleFvPatchScalarField
     turbulentTemperatureCoupledBaffleMixedFvPatchScalarField(ptf, iF),
     owner_(ptf.owner_),
     baffle_(ptf.baffle_),
-    solidThermoType_(ptf.solidThermoType_),
+    solidThermoType_(ptf.solidThermoType_().clone()),
     dict_(ptf.dict_)
 {}
 
@@ -219,8 +219,7 @@ void temperatureThermoBaffleFvPatchScalarField::write(Ostream& os) const
         os.writeKeyword(word(thermoModel + "Coeffs"));
         os << dict_.subDict(thermoModel + "Coeffs") << nl;
 
-        os.writeKeyword("thermoType") << solidThermoType_
-            << token::END_STATEMENT << nl;
+        os << solidThermoType_() << nl;
 
         os.writeKeyword("mixture");
         os << dict_.subDict("mixture") << nl;

src/regionModels/thermoBaffleModels/derivedFvPatchFields/temperatureThermoBaffle/temperatureThermoBaffleFvPatchScalarField.H

@@ -140,7 +140,7 @@ class temperatureThermoBaffleFvPatchScalarField
         autoPtr<regionModels::thermoBaffleModels::thermoBaffleModel> baffle_;
 
         //- Solid thermo type
-        word solidThermoType_;
+        autoPtr<entry> solidThermoType_;
 
         //- Dictionary
         dictionary dict_;

src/regionModels/thermoBaffleModels/noThermo/noThermo.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -88,7 +88,8 @@ void noThermo::evolveRegion()
 const tmp<volScalarField> noThermo::Cp() const
 {
     FatalErrorIn("const tmp<volScalarField>& noThermo::Cp() const")
-        << "Cp field not available for " << type() << abort(FatalError);
+        << "Cp field not available for " << type()
+        << abort(FatalError);
 
     return tmp<volScalarField>
     (
@@ -112,7 +113,8 @@ const tmp<volScalarField> noThermo::Cp() const
 const volScalarField& noThermo::kappaRad() const
 {
     FatalErrorIn("const volScalarField& noThermo::kappaRad() const")
-        << "kappa field not available for " << type() << abort(FatalError);
+        << "kappa field not available for " << type()
+        << abort(FatalError);
     return volScalarField::null();
 }
 
@@ -120,7 +122,8 @@ const volScalarField& noThermo::kappaRad() const
 const volScalarField& noThermo::rho() const
 {
     FatalErrorIn("const volScalarField& noThermo::rho() const")
-        << "rho field not available for " << type() << abort(FatalError);
+        << "rho field not available for " << type()
+        << abort(FatalError);
     return volScalarField::null();
 }
 
@@ -128,7 +131,8 @@ const volScalarField& noThermo::rho() const
 const volScalarField& noThermo::kappa() const
 {
    FatalErrorIn("const volScalarField& noThermo::kappa() const")
-        << "K field not available for " << type() << abort(FatalError);
+        << "K field not available for " << type()
+        << abort(FatalError);
     return volScalarField::null();
 }
 
@@ -136,7 +140,8 @@ const volScalarField& noThermo::kappa() const
 const volScalarField& noThermo::T() const
 {
     FatalErrorIn("const volScalarField& noThermo::T() const")
-        << "T field not available for " << type() << abort(FatalError);
+        << "T field not available for " << type()
+        << abort(FatalError);
     return volScalarField::null();
 }
 
@@ -144,7 +149,8 @@ const volScalarField& noThermo::T() const
 const solidThermo& noThermo::thermo() const
 {
     FatalErrorIn("const volScalarField& noThermo::T() const")
-        << "T field not available for " << type() << abort(FatalError);
+        << "T field not available for " << type()
+        << abort(FatalError);
     return reinterpret_cast<const solidThermo&>(null);
 }
 

src/regionModels/thermoBaffleModels/thermoBaffle2D/thermoBaffle2D.C

@@ -318,7 +318,7 @@ const tmp<volScalarField> thermoBaffle2D::Cp() const
 
 const volScalarField& thermoBaffle2D::kappaRad() const
 {
-    return thermo_->kappaRad();
+    return radiation_->absorptionEmission().a();
 }
 
 

src/regionModels/thermoBaffleModels/thermoBaffleModel/thermoBaffleModel.C

@@ -192,6 +192,9 @@ void thermoBaffleModel::init()
                 }
             }
         }
+
+        // Create radiation model
+        radiation_.reset(radiation::radiationModel::New(T()).ptr());
     }
 }
 
@@ -204,7 +207,8 @@ thermoBaffleModel::thermoBaffleModel(const fvMesh& mesh)
     thickness_(),
     delta_("delta", dimLength, 0.0),
     oneD_(false),
-    constantThickness_(true)
+    constantThickness_(true),
+    radiation_(NULL)
 {}
 
 
@@ -220,7 +224,8 @@ thermoBaffleModel::thermoBaffleModel
     thickness_(),
     delta_("delta", dimLength, 0.0),
     oneD_(false),
-    constantThickness_(dict.lookupOrDefault<bool>("constantThickness", true))
+    constantThickness_(dict.lookupOrDefault<bool>("constantThickness", true)),
+    radiation_(NULL)
 {
     init();
 }
@@ -232,7 +237,8 @@ thermoBaffleModel::thermoBaffleModel(const word& modelType, const fvMesh& mesh)
     thickness_(),
     delta_("delta", dimLength, 0.0),
     oneD_(false),
-    constantThickness_(lookupOrDefault<bool>("constantThickness", true))
+    constantThickness_(lookupOrDefault<bool>("constantThickness", true)),
+    radiation_(NULL)
 {
     init();
 }

src/regionModels/thermoBaffleModels/thermoBaffleModel/thermoBaffleModel.H

@@ -40,6 +40,8 @@ SourceFiles
 #include "volFieldsFwd.H"
 #include "solidThermo.H"
 #include "regionModel1D.H"
+#include "radiationModel.H"
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -88,6 +90,9 @@ protected:
         //- Is thickness constant
         bool constantThickness_;
 
+        //- Pointer to radiation model
+        autoPtr<radiation::radiationModel> radiation_;
+
 
     // Protected Member Functions
 

src/thermophysicalModels/basic/Make/files

@@ -1,11 +1,11 @@
 basicThermo/basicThermo.C
 fluidThermo/fluidThermo.C
 
-psiThermo/psiThermo/psiThermo.C
-psiThermo/hePsiThermo/hePsiThermos.C
+psiThermo/psiThermo.C
+psiThermo/psiThermos.C
 
-rhoThermo/rhoThermo/rhoThermo.C
-rhoThermo/heRhoThermo/heRhoThermos.C
+rhoThermo/rhoThermo.C
+rhoThermo/rhoThermos.C
 
 derivedFvPatchFields/fixedEnergy/fixedEnergyFvPatchScalarField.C
 derivedFvPatchFields/gradientEnergy/gradientEnergyFvPatchScalarField.C

src/thermophysicalModels/basic/basicThermo/basicThermo.H

@@ -105,10 +105,22 @@ public:
 
     // Selectors
 
+        //- Generic lookup for each of the related thermodynamics packages
+        template<class Thermo, class Table>
+        static typename Table::iterator lookupThermo
+        (
+            const dictionary& thermoDict,
+            Table* tablePtr
+        );
+
         //- Generic New for each of the related thermodynamics packages
         template<class Thermo>
         static autoPtr<Thermo> New(const fvMesh&);
 
+        //- Generic New for each of the related thermodynamics packages
+        template<class Thermo>
+        static autoPtr<Thermo> New(const fvMesh&, const dictionary&);
+
         //- Specialisation of the Generic New for basicThermo
         static autoPtr<basicThermo> New(const fvMesh&);
 

src/thermophysicalModels/basic/basicThermo/basicThermoTemplates.C

@@ -27,25 +27,13 @@ License
 
 // * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
 
-template<class Thermo>
-Foam::autoPtr<Thermo> Foam::basicThermo::New
+template<class Thermo, class Table>
+typename Table::iterator Foam::basicThermo::lookupThermo
 (
-    const fvMesh& mesh
+    const dictionary& thermoDict,
+    Table* tablePtr
 )
 {
-    IOdictionary thermoDict
-    (
-        IOobject
-        (
-            "thermophysicalProperties",
-            mesh.time().constant(),
-            mesh,
-            IOobject::MUST_READ_IF_MODIFIED,
-            IOobject::NO_WRITE,
-            false
-        )
-    );
-
     word thermoTypeName;
 
     if (thermoDict.isDict("thermoType"))
@@ -76,16 +64,13 @@ Foam::autoPtr<Thermo> Foam::basicThermo::New
           + word(thermoTypeDict.lookup("specie")) + ">>,"
           + word(thermoTypeDict.lookup("energy")) + ">>>";
 
-        Info<< thermoTypeName << endl;
-
         // Lookup the thermo package
-        typename Thermo::fvMeshConstructorTable::iterator cstrIter =
-            Thermo::fvMeshConstructorTablePtr_->find(thermoTypeName);
+        typename Table::iterator cstrIter = tablePtr->find(thermoTypeName);
 
         // Print error message if package not found in the table
-        if (cstrIter == Thermo::fvMeshConstructorTablePtr_->end())
+        if (cstrIter == tablePtr->end())
         {
-            FatalErrorIn(Thermo::typeName + "::New(const fvMesh&)")
+            FatalErrorIn(Thermo::typeName + "::New")
                 << "Unknown " << Thermo::typeName << " type " << nl
                 << "thermoType" << thermoTypeDict << nl << nl
                 << "Valid " << Thermo::typeName << " types are:" << nl << nl;
@@ -93,7 +78,7 @@ Foam::autoPtr<Thermo> Foam::basicThermo::New
             // Get the list of all the suitable thermo packages available
             wordList validThermoTypeNames
             (
-                Thermo::fvMeshConstructorTablePtr_->sortedToc()
+                tablePtr->sortedToc()
             );
 
             // Build a table of the thermo packages constituent parts
@@ -123,7 +108,7 @@ Foam::autoPtr<Thermo> Foam::basicThermo::New
             FatalError<< exit(FatalError);
         }
 
-        return autoPtr<Thermo>(cstrIter()(mesh));
+        return cstrIter;
     }
     else
     {
@@ -131,22 +116,69 @@ Foam::autoPtr<Thermo> Foam::basicThermo::New
 
         Info<< "Selecting thermodynamics package " << thermoTypeName << endl;
 
-        typename Thermo::fvMeshConstructorTable::iterator cstrIter =
-            Thermo::fvMeshConstructorTablePtr_->find(thermoTypeName);
+        typename Table::iterator cstrIter = tablePtr->find(thermoTypeName);
 
-        if (cstrIter == Thermo::fvMeshConstructorTablePtr_->end())
+        if (cstrIter == tablePtr->end())
         {
-            FatalErrorIn(Thermo::typeName + "::New(const fvMesh&)")
+            FatalErrorIn(Thermo::typeName + "::New")
                 << "Unknown " << Thermo::typeName << " type "
                 << thermoTypeName << nl << nl
                 << "Valid " << Thermo::typeName << " types are:" << nl
-                << Thermo::fvMeshConstructorTablePtr_->sortedToc() << nl
+                << tablePtr->sortedToc() << nl
                 << exit(FatalError);
         }
 
-        return autoPtr<Thermo>(cstrIter()(mesh));
+        return cstrIter;
     }
 }
 
 
+template<class Thermo>
+Foam::autoPtr<Thermo> Foam::basicThermo::New
+(
+    const fvMesh& mesh
+)
+{
+    IOdictionary thermoDict
+    (
+        IOobject
+        (
+            "thermophysicalProperties",
+            mesh.time().constant(),
+            mesh,
+            IOobject::MUST_READ_IF_MODIFIED,
+            IOobject::NO_WRITE,
+            false
+        )
+    );
+
+    typename Thermo::fvMeshConstructorTable::iterator cstrIter =
+        lookupThermo<Thermo, typename Thermo::fvMeshConstructorTable>
+        (
+            thermoDict,
+            Thermo::fvMeshConstructorTablePtr_
+        );
+
+    return autoPtr<Thermo>(cstrIter()(mesh));
+}
+
+
+template<class Thermo>
+Foam::autoPtr<Thermo> Foam::basicThermo::New
+(
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+{
+    typename Thermo::dictionaryConstructorTable::iterator cstrIter =
+        lookupThermo<Thermo, typename Thermo::dictionaryConstructorTable>
+        (
+            dict,
+            Thermo::dictionaryConstructorTablePtr_
+        );
+
+    return autoPtr<Thermo>(cstrIter()(mesh, dict));
+}
+
+
 // ************************************************************************* //

src/thermophysicalModels/basic/fluidThermo/makeThermo.H

@@ -37,7 +37,7 @@ Description
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#define makeThermo(BaseThermo,Cthermo,Mixture,Transport,Type,Thermo,EqnOfState,Specie)\
+#define makeThermoTypedefs(BaseThermo,Cthermo,Mixture,Transport,Type,Thermo,EqnOfState,Specie)\
                                                                               \
 typedef                                                                       \
     Transport                                                                 \
@@ -58,10 +58,8 @@ typedef                                                                       \
 typedef                                                                       \
     Cthermo                                                                   \
     <                                                                         \
-        Mixture                                                               \
-        <                                                                     \
-            Transport##Type##Thermo##EqnOfState##Specie                       \
-        >                                                                     \
+        BaseThermo,                                                           \
+        Mixture<Transport##Type##Thermo##EqnOfState##Specie>                  \
     > Cthermo##Mixture##Transport##Type##Thermo##EqnOfState##Specie;          \
                                                                               \
 defineTemplateTypeNameAndDebugWithName                                        \
@@ -73,7 +71,22 @@ defineTemplateTypeNameAndDebugWithName                                        \
       + ">>"                                                                  \
     ).c_str(),                                                                \
     0                                                                         \
-);                                                                            \
+);
+
+
+#define makeThermo(BaseThermo,Cthermo,Mixture,Transport,Type,Thermo,EqnOfState,Specie)\
+                                                                              \
+makeThermoTypedefs                                                            \
+(                                                                             \
+    BaseThermo,                                                               \
+    Cthermo,                                                                  \
+    Mixture,                                                                  \
+    Transport,                                                                \
+    Type,                                                                     \
+    Thermo,                                                                   \
+    EqnOfState,                                                               \
+    Specie                                                                    \
+)                                                                             \
                                                                               \
 addToRunTimeSelectionTable                                                    \
 (                                                                             \

src/thermophysicalModels/basic/psiThermo/hePsiThermo.C

@@ -27,8 +27,8 @@ License
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-template<class MixtureType>
-void Foam::hePsiThermo<MixtureType>::calculate()
+template<class BasicPsiThermo, class MixtureType>
+void Foam::hePsiThermo<BasicPsiThermo, MixtureType>::calculate()
 {
     const scalarField& hCells = this->he_.internalField();
     const scalarField& pCells = this->p_.internalField();
@@ -101,10 +101,10 @@ void Foam::hePsiThermo<MixtureType>::calculate()
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-template<class MixtureType>
-Foam::hePsiThermo<MixtureType>::hePsiThermo(const fvMesh& mesh)
+template<class BasicPsiThermo, class MixtureType>
+Foam::hePsiThermo<BasicPsiThermo, MixtureType>::hePsiThermo(const fvMesh& mesh)
 :
-    heThermo<psiThermo, MixtureType>(mesh)
+    heThermo<BasicPsiThermo, MixtureType>(mesh)
 {
     calculate();
 
@@ -115,19 +115,20 @@ Foam::hePsiThermo<MixtureType>::hePsiThermo(const fvMesh& mesh)
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-template<class MixtureType>
-Foam::hePsiThermo<MixtureType>::~hePsiThermo()
+template<class BasicPsiThermo, class MixtureType>
+Foam::hePsiThermo<BasicPsiThermo, MixtureType>::~hePsiThermo()
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class MixtureType>
-void Foam::hePsiThermo<MixtureType>::correct()
+template<class BasicPsiThermo, class MixtureType>
+void Foam::hePsiThermo<BasicPsiThermo, MixtureType>::correct()
 {
     if (debug)
     {
-        Info<< "entering hePsiThermo<MixtureType>::correct()" << endl;
+        Info<< "entering hePsiThermo<BasicPsiThermo, MixtureType>::correct()"
+            << endl;
     }
 
     // force the saving of the old-time values
@@ -137,7 +138,8 @@ void Foam::hePsiThermo<MixtureType>::correct()
 
     if (debug)
     {
-        Info<< "exiting hePsiThermo<MixtureType>::correct()" << endl;
+        Info<< "exiting hePsiThermo<BasicPsiThermo, MixtureType>::correct()"
+            << endl;
     }
 }
 

src/thermophysicalModels/basic/psiThermo/hePsiThermo.H

@@ -25,7 +25,7 @@ Class
     Foam::hePsiThermo
 
 Description
-    Enthalpy for a mixture based on compressibility
+    Energy for a mixture based on compressibility
 
 SourceFiles
     hePsiThermo.C
@@ -47,10 +47,10 @@ namespace Foam
                          Class hePsiThermo Declaration
 \*---------------------------------------------------------------------------*/
 
-template<class MixtureType>
+template<class BasicPsiThermo, class MixtureType>
 class hePsiThermo
 :
-    public heThermo<psiThermo, MixtureType>
+    public heThermo<BasicPsiThermo, MixtureType>
 {
     // Private Member Functions
 
@@ -58,7 +58,7 @@ class hePsiThermo
         void calculate();
 
         //- Construct as copy (not implemented)
-        hePsiThermo(const hePsiThermo<MixtureType>&);
+        hePsiThermo(const hePsiThermo<BasicPsiThermo, MixtureType>&);
 
 public:
 

src/thermophysicalModels/basic/rhoThermo/heRhoThermo.C

@@ -27,8 +27,8 @@ License
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-template<class MixtureType>
-void Foam::heRhoThermo<MixtureType>::calculate()
+template<class BasicPsiThermo, class MixtureType>
+void Foam::heRhoThermo<BasicPsiThermo, MixtureType>::calculate()
 {
     const scalarField& hCells = this->he().internalField();
     const scalarField& pCells = this->p_.internalField();
@@ -106,10 +106,10 @@ void Foam::heRhoThermo<MixtureType>::calculate()
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-template<class MixtureType>
-Foam::heRhoThermo<MixtureType>::heRhoThermo(const fvMesh& mesh)
+template<class BasicPsiThermo, class MixtureType>
+Foam::heRhoThermo<BasicPsiThermo, MixtureType>::heRhoThermo(const fvMesh& mesh)
 :
-    heThermo<rhoThermo, MixtureType>(mesh)
+    heThermo<BasicPsiThermo, MixtureType>(mesh)
 {
     calculate();
 }
@@ -117,15 +117,15 @@ Foam::heRhoThermo<MixtureType>::heRhoThermo(const fvMesh& mesh)
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-template<class MixtureType>
-Foam::heRhoThermo<MixtureType>::~heRhoThermo()
+template<class BasicPsiThermo, class MixtureType>
+Foam::heRhoThermo<BasicPsiThermo, MixtureType>::~heRhoThermo()
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class MixtureType>
-void Foam::heRhoThermo<MixtureType>::correct()
+template<class BasicPsiThermo, class MixtureType>
+void Foam::heRhoThermo<BasicPsiThermo, MixtureType>::correct()
 {
     if (debug)
     {

src/thermophysicalModels/basic/rhoThermo/heRhoThermo.H

@@ -25,7 +25,7 @@ Class
     Foam::heRhoThermo
 
 Description
-    Enthalpy for a mixture based on density
+    Energy for a mixture based on density
 
 SourceFiles
     heRhoThermo.C
@@ -47,10 +47,10 @@ namespace Foam
                          Class heRhoThermo Declaration
 \*---------------------------------------------------------------------------*/
 
-template<class MixtureType>
+template<class BasicPsiThermo, class MixtureType>
 class heRhoThermo
 :
-    public heThermo<rhoThermo, MixtureType>
+    public heThermo<BasicPsiThermo, MixtureType>
 {
     // Private Member Functions
 
@@ -58,7 +58,7 @@ class heRhoThermo
         void calculate();
 
         //- Construct as copy (not implemented)
-        heRhoThermo(const heRhoThermo<MixtureType>&);
+        heRhoThermo(const heRhoThermo<BasicPsiThermo, MixtureType>&);
 
 
 public:

src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModelTemplates.C

@@ -58,12 +58,11 @@ Foam::autoPtr<ChemistryModel> Foam::basicChemistryModel::New
 
         Info<< "Selecting chemistry type " << chemistryTypeDict << endl;
 
-        const int nCmpt = 8;
+        const int nCmpt = 7;
         const char* cmptNames[nCmpt] =
         {
             "chemistrySolver",
-            "chemistryModel",
-            "???ChemistryModel",
+            "chemistryThermo",
             "transport",
             "thermo",
             "equationOfState",
@@ -109,9 +108,8 @@ Foam::autoPtr<ChemistryModel> Foam::basicChemistryModel::New
         // Construct the name of the chemistry type from the components
         chemistryTypeName =
             word(chemistryTypeDict.lookup("chemistrySolver")) + '<'
-          + word(chemistryTypeDict.lookup("chemistryModel")) + '<'
-          + ChemistryModel::typeName + ','
-          + thermoTypeName + ">>";
+          + word(chemistryTypeDict.lookup("chemistryThermo")) + ','
+          + thermoTypeName + ">";
 
         typename ChemistryModel::fvMeshConstructorTable::iterator cstrIter =
             ChemistryModel::fvMeshConstructorTablePtr_->find(chemistryTypeName);
@@ -160,8 +158,7 @@ Foam::autoPtr<ChemistryModel> Foam::basicChemistryModel::New
             FatalError<< exit(FatalError);
         }
 
-        return autoPtr<ChemistryModel>
-            (cstrIter()(mesh, typeName, chemistryTypeName));
+        return autoPtr<ChemistryModel>(cstrIter()(mesh));
     }
     else
     {
@@ -183,8 +180,7 @@ Foam::autoPtr<ChemistryModel> Foam::basicChemistryModel::New
                 << exit(FatalError);
         }
 
-        return autoPtr<ChemistryModel>
-            (cstrIter()(mesh, typeName, chemistryTypeName));
+        return autoPtr<ChemistryModel>(cstrIter()(mesh));
     }
 }
 

src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C

@@ -23,26 +23,20 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "ODEChemistryModel.H"
-#include "chemistrySolver.H"
+#include "chemistryModel.H"
 #include "reactingMixture.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class CompType, class ThermoType>
-Foam::ODEChemistryModel<CompType, ThermoType>::ODEChemistryModel
+Foam::chemistryModel<CompType, ThermoType>::chemistryModel
 (
-    const fvMesh& mesh,
-    const word& ODEModelName,
-    const word& thermoTypeName
+    const fvMesh& mesh
 )
 :
-    CompType(mesh, thermoTypeName),
-
+    CompType(mesh),
     ODE(),
-
     Y_(this->thermo().composition().Y()),
-
     reactions_
     (
         dynamic_cast<const reactingMixture<ThermoType>&>(this->thermo())
@@ -68,7 +62,7 @@ Foam::ODEChemistryModel<CompType, ThermoType>::ODEChemistryModel
             (
                 IOobject
                 (
-                    "RR::" + Y_[fieldI].name(),
+                    "RR." + Y_[fieldI].name(),
                     mesh.time().timeName(),
                     mesh,
                     IOobject::NO_READ,
@@ -80,7 +74,7 @@ Foam::ODEChemistryModel<CompType, ThermoType>::ODEChemistryModel
         );
     }
 
-    Info<< "ODEChemistryModel: Number of species = " << nSpecie_
+    Info<< "chemistryModel: Number of species = " << nSpecie_
         << " and reactions = " << nReaction_ << endl;
 }
 
@@ -88,7 +82,7 @@ Foam::ODEChemistryModel<CompType, ThermoType>::ODEChemistryModel
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 template<class CompType, class ThermoType>
-Foam::ODEChemistryModel<CompType, ThermoType>::~ODEChemistryModel()
+Foam::chemistryModel<CompType, ThermoType>::~chemistryModel()
 {}
 
 
@@ -96,7 +90,7 @@ Foam::ODEChemistryModel<CompType, ThermoType>::~ODEChemistryModel()
 
 template<class CompType, class ThermoType>
 Foam::tmp<Foam::scalarField>
-Foam::ODEChemistryModel<CompType, ThermoType>::omega
+Foam::chemistryModel<CompType, ThermoType>::omega
 (
     const scalarField& c,
     const scalar T,
@@ -138,7 +132,7 @@ Foam::ODEChemistryModel<CompType, ThermoType>::omega
 
 
 template<class CompType, class ThermoType>
-Foam::scalar Foam::ODEChemistryModel<CompType, ThermoType>::omegaI
+Foam::scalar Foam::chemistryModel<CompType, ThermoType>::omegaI
 (
     const label index,
     const scalarField& c,
@@ -160,7 +154,7 @@ Foam::scalar Foam::ODEChemistryModel<CompType, ThermoType>::omegaI
 
 
 template<class CompType, class ThermoType>
-void Foam::ODEChemistryModel<CompType, ThermoType>::updateConcsInReactionI
+void Foam::chemistryModel<CompType, ThermoType>::updateConcsInReactionI
 (
     const label index,
     const scalar dt,
@@ -189,7 +183,7 @@ void Foam::ODEChemistryModel<CompType, ThermoType>::updateConcsInReactionI
 
 
 template<class CompType, class ThermoType>
-void Foam::ODEChemistryModel<CompType, ThermoType>::updateRRInReactionI
+void Foam::chemistryModel<CompType, ThermoType>::updateRRInReactionI
 (
     const label index,
     const scalar pr,
@@ -220,7 +214,7 @@ void Foam::ODEChemistryModel<CompType, ThermoType>::updateRRInReactionI
 
 
 template<class CompType, class ThermoType>
-Foam::scalar Foam::ODEChemistryModel<CompType, ThermoType>::omega
+Foam::scalar Foam::chemistryModel<CompType, ThermoType>::omega
 (
     const Reaction<ThermoType>& R,
     const scalarField& c,
@@ -338,7 +332,7 @@ Foam::scalar Foam::ODEChemistryModel<CompType, ThermoType>::omega
 
 
 template<class CompType, class ThermoType>
-void Foam::ODEChemistryModel<CompType, ThermoType>::derivatives
+void Foam::chemistryModel<CompType, ThermoType>::derivatives
 (
     const scalar time,
     const scalarField &c,
@@ -389,7 +383,7 @@ void Foam::ODEChemistryModel<CompType, ThermoType>::derivatives
 
 
 template<class CompType, class ThermoType>
-void Foam::ODEChemistryModel<CompType, ThermoType>::jacobian
+void Foam::chemistryModel<CompType, ThermoType>::jacobian
 (
     const scalar t,
     const scalarField& c,
@@ -532,7 +526,7 @@ void Foam::ODEChemistryModel<CompType, ThermoType>::jacobian
 
 template<class CompType, class ThermoType>
 Foam::tmp<Foam::volScalarField>
-Foam::ODEChemistryModel<CompType, ThermoType>::tc() const
+Foam::chemistryModel<CompType, ThermoType>::tc() const
 {
     scalar pf, cf, pr, cr;
     label lRef, rRef;
@@ -570,6 +564,8 @@ Foam::ODEChemistryModel<CompType, ThermoType>::tc() const
     );
 
     scalarField& tc = ttc();
+    const scalarField& T = this->thermo().T();
+    const scalarField& p = this->thermo().p();
 
     const label nReaction = reactions_.size();
 
@@ -578,8 +574,8 @@ Foam::ODEChemistryModel<CompType, ThermoType>::tc() const
         forAll(rho, celli)
         {
             scalar rhoi = rho[celli];
-            scalar Ti = this->thermo().T()[celli];
-            scalar pi = this->thermo().p()[celli];
+            scalar Ti = T[celli];
+            scalar pi = p[celli];
             scalarField c(nSpecie_);
             scalar cSum = 0.0;
 
@@ -615,7 +611,7 @@ Foam::ODEChemistryModel<CompType, ThermoType>::tc() const
 
 template<class CompType, class ThermoType>
 Foam::tmp<Foam::volScalarField>
-Foam::ODEChemistryModel<CompType, ThermoType>::Sh() const
+Foam::chemistryModel<CompType, ThermoType>::Sh() const
 {
     tmp<volScalarField> tSh
     (
@@ -657,7 +653,7 @@ Foam::ODEChemistryModel<CompType, ThermoType>::Sh() const
 
 template<class CompType, class ThermoType>
 Foam::tmp<Foam::volScalarField>
-Foam::ODEChemistryModel<CompType, ThermoType>::dQ() const
+Foam::chemistryModel<CompType, ThermoType>::dQ() const
 {
     tmp<volScalarField> tdQ
     (
@@ -689,7 +685,7 @@ Foam::ODEChemistryModel<CompType, ThermoType>::dQ() const
 
 
 template<class CompType, class ThermoType>
-Foam::label Foam::ODEChemistryModel<CompType, ThermoType>::nEqns() const
+Foam::label Foam::chemistryModel<CompType, ThermoType>::nEqns() const
 {
     // nEqns = number of species + temperature + pressure
     return nSpecie_ + 2;
@@ -697,7 +693,7 @@ Foam::label Foam::ODEChemistryModel<CompType, ThermoType>::nEqns() const
 
 
 template<class CompType, class ThermoType>
-void Foam::ODEChemistryModel<CompType, ThermoType>::calculate()
+void Foam::chemistryModel<CompType, ThermoType>::calculate()
 {
     if (!this->chemistry_)
     {
@@ -718,11 +714,14 @@ void Foam::ODEChemistryModel<CompType, ThermoType>::calculate()
         this->thermo().rho()
     );
 
+    const scalarField& T = this->thermo().T();
+    const scalarField& p = this->thermo().p();
+
     forAll(rho, celli)
     {
         const scalar rhoi = rho[celli];
-        const scalar Ti = this->thermo().T()[celli];
-        const scalar pi = this->thermo().p()[celli];
+        const scalar Ti = T[celli];
+        const scalar pi = p[celli];
 
         scalarField c(nSpecie_, 0.0);
         for (label i=0; i<nSpecie_; i++)
@@ -742,7 +741,7 @@ void Foam::ODEChemistryModel<CompType, ThermoType>::calculate()
 
 
 template<class CompType, class ThermoType>
-Foam::scalar Foam::ODEChemistryModel<CompType, ThermoType>::solve
+Foam::scalar Foam::chemistryModel<CompType, ThermoType>::solve
 (
     const scalar t0,
     const scalar deltaT
@@ -773,13 +772,16 @@ Foam::scalar Foam::ODEChemistryModel<CompType, ThermoType>::solve
 
     tmp<volScalarField> thc = this->thermo().hc();
     const scalarField& hc = thc();
+    const scalarField& he = this->thermo().he();
+    const scalarField& T = this->thermo().T();
+    const scalarField& p = this->thermo().p();
 
     forAll(rho, celli)
     {
         const scalar rhoi = rho[celli];
-        const scalar hi = this->thermo().he()[celli] + hc[celli];
-        const scalar pi = this->thermo().p()[celli];
-        scalar Ti = this->thermo().T()[celli];
+        const scalar hi = he[celli] + hc[celli];
+        const scalar pi = p[celli];
+        scalar Ti = T[celli];
 
         scalarField c(nSpecie_, 0.0);
         scalarField c0(nSpecie_, 0.0);
@@ -833,7 +835,7 @@ Foam::scalar Foam::ODEChemistryModel<CompType, ThermoType>::solve
 
 
 template<class CompType, class ThermoType>
-Foam::scalar Foam::ODEChemistryModel<CompType, ThermoType>::solve
+Foam::scalar Foam::chemistryModel<CompType, ThermoType>::solve
 (
     scalarField &c,
     const scalar T,
@@ -844,7 +846,7 @@ Foam::scalar Foam::ODEChemistryModel<CompType, ThermoType>::solve
 {
     notImplemented
     (
-        "ODEChemistryModel::solve"
+        "chemistryModel::solve"
         "("
             "scalarField&, "
             "const scalar, "

src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H

@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::ODEChemistryModel
+    Foam::chemistryModel
 
 Description
     Extends base chemistry model by adding a thermo package, and ODE functions.
@@ -30,13 +30,13 @@ Description
     terms.
 
 SourceFiles
-    ODEChemistryModelI.H
-    ODEChemistryModel.C
+    chemistryModelI.H
+    chemistryModel.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef ODEChemistryModel_H
-#define ODEChemistryModel_H
+#ifndef chemistryModel_H
+#define chemistryModel_H
 
 #include "Reaction.H"
 #include "ODE.H"
@@ -53,11 +53,11 @@ namespace Foam
 class fvMesh;
 
 /*---------------------------------------------------------------------------*\
-                      Class ODEChemistryModel Declaration
+                      Class chemistryModel Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class CompType, class ThermoType>
-class ODEChemistryModel
+class chemistryModel
 :
     public CompType,
     public ODE
@@ -65,10 +65,10 @@ class ODEChemistryModel
     // Private Member Functions
 
         //- Disallow copy constructor
-        ODEChemistryModel(const ODEChemistryModel&);
+        chemistryModel(const chemistryModel&);
 
         //- Disallow default bitwise assignment
-        void operator=(const ODEChemistryModel&);
+        void operator=(const chemistryModel&);
 
 
 protected:
@@ -104,22 +104,17 @@ protected:
 public:
 
     //- Runtime type information
-    TypeName("ODEChemistryModel");
+    TypeName("chemistryModel");
 
 
     // Constructors
 
-        //- Construct from components
-        ODEChemistryModel
-        (
-            const fvMesh& mesh,
-            const word& ODEModelName,
-            const word& thermoTypeName
-        );
+        //- Construct from mesh
+        chemistryModel(const fvMesh& mesh);
 
 
     //- Destructor
-    virtual ~ODEChemistryModel();
+    virtual ~chemistryModel();
 
 
     // Member Functions
@@ -263,12 +258,12 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#include "ODEChemistryModelI.H"
+#include "chemistryModelI.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #ifdef NoRepository
-#   include "ODEChemistryModel.C"
+#   include "chemistryModel.C"
 #endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModelI.H

@@ -30,7 +30,7 @@ License
 
 template<class CompType, class ThermoType>
 inline Foam::PtrList<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >&
-Foam::ODEChemistryModel<CompType, ThermoType>::RR()
+Foam::chemistryModel<CompType, ThermoType>::RR()
 {
     return RR_;
 }
@@ -38,7 +38,7 @@ Foam::ODEChemistryModel<CompType, ThermoType>::RR()
 
 template<class CompType, class ThermoType>
 inline const Foam::PtrList<Foam::Reaction<ThermoType> >&
-Foam::ODEChemistryModel<CompType, ThermoType>::reactions() const
+Foam::chemistryModel<CompType, ThermoType>::reactions() const
 {
     return reactions_;
 }
@@ -46,7 +46,7 @@ Foam::ODEChemistryModel<CompType, ThermoType>::reactions() const
 
 template<class CompType, class ThermoType>
 inline const Foam::PtrList<ThermoType>&
-Foam::ODEChemistryModel<CompType, ThermoType>::specieThermo() const
+Foam::chemistryModel<CompType, ThermoType>::specieThermo() const
 {
     return specieThermo_;
 }
@@ -54,7 +54,7 @@ Foam::ODEChemistryModel<CompType, ThermoType>::specieThermo() const
 
 template<class CompType, class ThermoType>
 inline Foam::label
-Foam::ODEChemistryModel<CompType, ThermoType>::nSpecie() const
+Foam::chemistryModel<CompType, ThermoType>::nSpecie() const
 {
     return nSpecie_;
 }
@@ -62,7 +62,7 @@ Foam::ODEChemistryModel<CompType, ThermoType>::nSpecie() const
 
 template<class CompType, class ThermoType>
 inline Foam::label
-Foam::ODEChemistryModel<CompType, ThermoType>::nReaction() const
+Foam::chemistryModel<CompType, ThermoType>::nReaction() const
 {
     return nReaction_;
 }
@@ -70,7 +70,7 @@ Foam::ODEChemistryModel<CompType, ThermoType>::nReaction() const
 
 template<class CompType, class ThermoType>
 inline const Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
-Foam::ODEChemistryModel<CompType, ThermoType>::RR
+Foam::chemistryModel<CompType, ThermoType>::RR
 (
     const label i
 ) const

src/thermophysicalModels/chemistryModel/chemistryModel/psiChemistryModel/psiChemistryModel.C

@@ -25,7 +25,6 @@ License
 
 #include "psiChemistryModel.H"
 #include "fvMesh.H"
-#include "Time.H"
 
 /* * * * * * * * * * * * * * * private static data * * * * * * * * * * * * * */
 
@@ -35,12 +34,12 @@ namespace Foam
     defineRunTimeSelectionTable(psiChemistryModel, fvMesh);
 }
 
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::psiChemistryModel::psiChemistryModel
 (
-    const fvMesh& mesh,
-    const word& thermoTypeName
+    const fvMesh& mesh
 )
 :
     basicChemistryModel(mesh),

src/thermophysicalModels/chemistryModel/chemistryModel/psiChemistryModel/psiChemistryModel.H

@@ -78,7 +78,7 @@ protected:
 public:
 
     //- Runtime type information
-    TypeName("psiChemistryModel");
+    TypeName("psi");
 
 
     //- Declare run-time constructor selection tables
@@ -87,19 +87,15 @@ public:
         autoPtr,
         psiChemistryModel,
         fvMesh,
-        (
-            const fvMesh& mesh,
-            const word& compTypeName,
-            const word& thermoTypeName
-        ),
-        (mesh, compTypeName, thermoTypeName)
+        (const fvMesh& mesh),
+        (mesh)
     );
 
 
     // Constructors
 
         //- Construct from mesh
-        psiChemistryModel(const fvMesh& mesh, const word& thermoTypeName);
+        psiChemistryModel(const fvMesh& mesh);
 
 
     //- Selector

src/thermophysicalModels/chemistryModel/chemistryModel/psiChemistryModel/psiChemistryModels.C

@@ -32,7 +32,7 @@ Description
 #include "makeChemistryModel.H"
 
 #include "psiChemistryModel.H"
-#include "ODEChemistryModel.H"
+#include "chemistryModel.H"
 #include "thermoPhysicsTypes.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -41,35 +41,35 @@ namespace Foam
 {
     makeChemistryModel
     (
-        ODEChemistryModel,
+        chemistryModel,
         psiChemistryModel,
         constGasThermoPhysics
     );
 
     makeChemistryModel
     (
-        ODEChemistryModel,
+        chemistryModel,
         psiChemistryModel,
         gasThermoPhysics
     );
 
     makeChemistryModel
     (
-        ODEChemistryModel,
+        chemistryModel,
         psiChemistryModel,
         constIncompressibleGasThermoPhysics
     );
 
     makeChemistryModel
     (
-        ODEChemistryModel,
+        chemistryModel,
         psiChemistryModel,
         incompressibleGasThermoPhysics
     );
 
     makeChemistryModel
     (
-        ODEChemistryModel,
+        chemistryModel,
         psiChemistryModel,
         icoPoly8ThermoPhysics
     );

src/thermophysicalModels/chemistryModel/chemistryModel/rhoChemistryModel/rhoChemistryModel.C

@@ -25,7 +25,6 @@ License
 
 #include "rhoChemistryModel.H"
 #include "fvMesh.H"
-#include "Time.H"
 
 /* * * * * * * * * * * * * * * private static data * * * * * * * * * * * * * */
 
@@ -35,12 +34,12 @@ namespace Foam
     defineRunTimeSelectionTable(rhoChemistryModel, fvMesh);
 }
 
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::rhoChemistryModel::rhoChemistryModel
 (
-    const fvMesh& mesh,
-    const word& thermoTypeName
+    const fvMesh& mesh
 )
 :
     basicChemistryModel(mesh),

src/thermophysicalModels/chemistryModel/chemistryModel/rhoChemistryModel/rhoChemistryModel.H

@@ -78,7 +78,7 @@ protected:
 public:
 
     //- Runtime type information
-    TypeName("rhoChemistryModel");
+    TypeName("rho");
 
 
     //- Declare run-time constructor selection tables
@@ -88,18 +88,16 @@ public:
         rhoChemistryModel,
         fvMesh,
         (
-            const fvMesh& mesh,
-            const word& compTypeName,
-            const word& thermoTypeName
+            const fvMesh& mesh
         ),
-        (mesh, compTypeName, thermoTypeName)
+        (mesh)
     );
 
 
     // Constructors
 
         //- Construct from mesh
-        rhoChemistryModel(const fvMesh& mesh, const word& thermoTypeName);
+        rhoChemistryModel(const fvMesh& mesh);
 
 
     //- Selector

src/thermophysicalModels/chemistryModel/chemistryModel/rhoChemistryModel/rhoChemistryModels.C

@@ -32,7 +32,7 @@ Description
 #include "makeChemistryModel.H"
 
 #include "rhoChemistryModel.H"
-#include "ODEChemistryModel.H"
+#include "chemistryModel.H"
 #include "thermoPhysicsTypes.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -41,35 +41,35 @@ namespace Foam
 {
     makeChemistryModel
     (
-        ODEChemistryModel,
+        chemistryModel,
         rhoChemistryModel,
         constGasThermoPhysics
     );
 
     makeChemistryModel
     (
-        ODEChemistryModel,
+        chemistryModel,
         rhoChemistryModel,
         gasThermoPhysics
     );
 
     makeChemistryModel
     (
-        ODEChemistryModel,
+        chemistryModel,
         rhoChemistryModel,
         constIncompressibleGasThermoPhysics
     );
 
     makeChemistryModel
     (
-        ODEChemistryModel,
+        chemistryModel,
         rhoChemistryModel,
         incompressibleGasThermoPhysics
     );
 
     makeChemistryModel
     (
-        ODEChemistryModel,
+        chemistryModel,
         rhoChemistryModel,
         icoPoly8ThermoPhysics
     );

src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.C

@@ -29,15 +29,13 @@ License
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-template<class ODEChemistryType>
-Foam::EulerImplicit<ODEChemistryType>::EulerImplicit
+template<class ChemistryModel>
+Foam::EulerImplicit<ChemistryModel>::EulerImplicit
 (
-    const fvMesh& mesh,
-    const word& ODEModelName,
-    const word& thermoType
+    const fvMesh& mesh
 )
 :
-    chemistrySolver<ODEChemistryType>(mesh, ODEModelName, thermoType),
+    chemistrySolver<ChemistryModel>(mesh),
     coeffsDict_(this->subDict("EulerImplicitCoeffs")),
     cTauChem_(readScalar(coeffsDict_.lookup("cTauChem"))),
     eqRateLimiter_(coeffsDict_.lookup("equilibriumRateLimiter"))
@@ -46,15 +44,15 @@ Foam::EulerImplicit<ODEChemistryType>::EulerImplicit
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-template<class ODEChemistryType>
-Foam::EulerImplicit<ODEChemistryType>::~EulerImplicit()
+template<class ChemistryModel>
+Foam::EulerImplicit<ChemistryModel>::~EulerImplicit()
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class ODEChemistryType>
-Foam::scalar Foam::EulerImplicit<ODEChemistryType>::solve
+template<class ChemistryModel>
+Foam::scalar Foam::EulerImplicit<ChemistryModel>::solve
 (
     scalarField &c,
     const scalar T,

src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.H

@@ -46,10 +46,10 @@ namespace Foam
                        Class EulerImplicit Declaration
 \*---------------------------------------------------------------------------*/
 
-template<class ODEChemistryType>
+template<class ChemistryModel>
 class EulerImplicit
 :
-    public chemistrySolver<ODEChemistryType>
+    public chemistrySolver<ChemistryModel>
 {
     // Private data
 
@@ -74,13 +74,8 @@ public:
 
     // Constructors
 
-        //- Construct from components
-        EulerImplicit
-        (
-            const fvMesh& mesh,
-            const word& ODEmodelName,
-            const word& thermoType
-        );
+        //- Construct from mesh
+        EulerImplicit(const fvMesh& mesh);
 
 
     //- Destructor

src/thermophysicalModels/chemistryModel/chemistrySolver/chemistrySolver/chemistrySolver.C

@@ -27,23 +27,20 @@ License
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-template<class ODEChemistryType>
-Foam::chemistrySolver<ODEChemistryType>::chemistrySolver
+template<class ChemistryModel>
+Foam::chemistrySolver<ChemistryModel>::chemistrySolver
 (
-    const fvMesh& mesh,
-    const word& ODEModelName,
-    const word& thermoType
+    const fvMesh& mesh
 )
 :
-    ODEChemistryType(mesh, ODEModelName, thermoType),
-    name_(ODEModelName)
+    ChemistryModel(mesh)
 {}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-template<class ODEChemistryType>
-Foam::chemistrySolver<ODEChemistryType>::~chemistrySolver()
+template<class ChemistryModel>
+Foam::chemistrySolver<ChemistryModel>::~chemistrySolver()
 {}
 
 

src/thermophysicalModels/chemistryModel/chemistrySolver/chemistrySolver/chemistrySolver.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,7 +35,7 @@ SourceFiles
 #ifndef chemistrySolver_H
 #define chemistrySolver_H
 
-#include "ODEChemistryModel.H"
+#include "chemistryModel.H"
 #include "IOdictionary.H"
 #include "scalarField.H"
 
@@ -48,30 +48,18 @@ namespace Foam
                       Class chemistrySolver Declaration
 \*---------------------------------------------------------------------------*/
 
-template<class ODEChemistryType>
+template<class ChemistryModel>
 class chemistrySolver
 :
-    public ODEChemistryType
+    public ChemistryModel
 {
-protected:
-
-    // Protected data
-
-        //- Name of the chemistry solver
-        const word name_;
-
 
 public:
 
     // Constructors
 
         //- Construct from components
-        chemistrySolver
-        (
-            const fvMesh& mesh,
-            const word& ODEModelName,
-            const word& thermoType
-        );
+        chemistrySolver(const fvMesh& mesh);
 
 
     //- Destructor

src/thermophysicalModels/chemistryModel/chemistrySolver/chemistrySolver/makeChemistrySolverTypes.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -28,7 +28,7 @@ License
 
 #include "chemistrySolver.H"
 
-#include "ODEChemistryModel.H"
+#include "chemistryModel.H"
 
 #include "noChemistrySolver.H"
 #include "EulerImplicit.H"
@@ -37,21 +37,22 @@ License
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#define makeChemistrySolverType(SS, ODEChem, Comp, Thermo)                    \
+#define makeChemistrySolverType(SS, Comp, Thermo)                             \
                                                                               \
-    typedef SS<ODEChem<Comp, Thermo> > SS##ODEChem##Comp##Thermo;             \
+    typedef SS<chemistryModel<Comp, Thermo> > SS##Comp##Thermo;               \
                                                                               \
     defineTemplateTypeNameAndDebugWithName                                    \
     (                                                                         \
-        SS##ODEChem##Comp##Thermo,                                            \
-        (#SS"<"#ODEChem"<"#Comp"," + Thermo::typeName() + ">>").c_str(),      \
+        SS##Comp##Thermo,                                                     \
+        (#SS"<" + word(Comp::typeName_())                                     \
+      + "," + Thermo::typeName() + ">").c_str(),                              \
         0                                                                     \
     );                                                                        \
                                                                               \
     addToRunTimeSelectionTable                                                \
     (                                                                         \
         Comp,                                                                 \
-        SS##ODEChem##Comp##Thermo,                                            \
+        SS##Comp##Thermo,                                                     \
         fvMesh                                                                \
     );
 
@@ -61,7 +62,6 @@ License
     makeChemistrySolverType                                                   \
     (                                                                         \
         noChemistrySolver,                                                    \
-        ODEChemistryModel,                                                    \
         CompChemModel,                                                        \
         Thermo                                                                \
     );                                                                        \
@@ -69,7 +69,6 @@ License
     makeChemistrySolverType                                                   \
     (                                                                         \
         EulerImplicit,                                                        \
-        ODEChemistryModel,                                                    \
         CompChemModel,                                                        \
         Thermo                                                                \
     );                                                                        \
@@ -77,7 +76,6 @@ License
     makeChemistrySolverType                                                   \
     (                                                                         \
         ode,                                                                  \
-        ODEChemistryModel,                                                    \
         CompChemModel,                                                        \
         Thermo                                                                \
     );                                                                        \
@@ -85,7 +83,6 @@ License
     makeChemistrySolverType                                                   \
     (                                                                         \
         sequential,                                                           \
-        ODEChemistryModel,                                                    \
         CompChemModel,                                                        \
         Thermo                                                                \
     );

src/thermophysicalModels/chemistryModel/chemistrySolver/noChemistrySolver/noChemistrySolver.C

@@ -28,29 +28,27 @@ License
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-template<class ODEChemistryType>
-Foam::noChemistrySolver<ODEChemistryType>::noChemistrySolver
+template<class ChemistryModel>
+Foam::noChemistrySolver<ChemistryModel>::noChemistrySolver
 (
-    const fvMesh& mesh,
-    const word& ODEModelName,
-    const word& thermoType
+    const fvMesh& mesh
 )
 :
-    chemistrySolver<ODEChemistryType>(mesh, ODEModelName, thermoType)
+    chemistrySolver<ChemistryModel>(mesh)
 {}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-template<class ODEChemistryType>
-Foam::noChemistrySolver<ODEChemistryType>::~noChemistrySolver()
+template<class ChemistryModel>
+Foam::noChemistrySolver<ChemistryModel>::~noChemistrySolver()
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class ODEChemistryType>
-Foam::scalar Foam::noChemistrySolver<ODEChemistryType>::solve
+template<class ChemistryModel>
+Foam::scalar Foam::noChemistrySolver<ChemistryModel>::solve
 (
     scalarField&,
     const scalar,

src/thermophysicalModels/chemistryModel/chemistrySolver/noChemistrySolver/noChemistrySolver.H

@@ -47,10 +47,10 @@ namespace Foam
                      Class noChemistrySolver Declaration
 \*---------------------------------------------------------------------------*/
 
-template<class ODEChemistryType>
+template<class ChemistryModel>
 class noChemistrySolver
 :
-    public chemistrySolver<ODEChemistryType>
+    public chemistrySolver<ChemistryModel>
 {
 
 public:
@@ -61,14 +61,8 @@ public:
 
     // Constructors
 
-
         //- Construct from components
-        noChemistrySolver
-        (
-            const fvMesh& mesh,
-            const word& ODEModelName,
-            const word& thermoType
-        );
+        noChemistrySolver(const fvMesh& mesh);
 
 
     //- Destructor

src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,19 +24,17 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "ode.H"
-#include "ODEChemistryModel.H"
+#include "chemistryModel.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-template<class ODEChemistryType>
-Foam::ode<ODEChemistryType>::ode
+template<class ChemistryModel>
+Foam::ode<ChemistryModel>::ode
 (
-    const fvMesh& mesh,
-    const word& ODEModelName,
-    const word& thermoType
+    const fvMesh& mesh
 )
 :
-    chemistrySolver<ODEChemistryType>(mesh, ODEModelName, thermoType),
+    chemistrySolver<ChemistryModel>(mesh),
     coeffsDict_(this->subDict("odeCoeffs")),
     solverName_(coeffsDict_.lookup("solver")),
     odeSolver_(ODESolver::New(solverName_, *this)),
@@ -46,15 +44,15 @@ Foam::ode<ODEChemistryType>::ode
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-template<class ODEChemistryType>
-Foam::ode<ODEChemistryType>::~ode()
+template<class ChemistryModel>
+Foam::ode<ChemistryModel>::~ode()
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class ODEChemistryType>
-Foam::scalar Foam::ode<ODEChemistryType>::solve
+template<class ChemistryModel>
+Foam::scalar Foam::ode<ChemistryModel>::solve
 (
     scalarField& c,
     const scalar T,

src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.H

@@ -47,10 +47,10 @@ namespace Foam
                             Class ode Declaration
 \*---------------------------------------------------------------------------*/
 
-template<class ODEChemistryType>
+template<class ChemistryModel>
 class ode
 :
-    public chemistrySolver<ODEChemistryType>
+    public chemistrySolver<ChemistryModel>
 {
     // Private data
 
@@ -66,18 +66,13 @@ class ode
 public:
 
     //- Runtime type information
-    TypeName("ODE");
+    TypeName("ode");
 
 
     // Constructors
 
-        //- Construct from components
-        ode
-        (
-            const fvMesh& mesh,
-            const word& ODEModeNewlName,
-            const word& thermoType
-        );
+        //- Construct from mesh
+        ode(const fvMesh& mesh);
 
 
     //- Destructor

src/thermophysicalModels/chemistryModel/chemistrySolver/sequential/sequential.C

@@ -28,15 +28,13 @@ License
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-template<class ODEChemistryType>
-Foam::sequential<ODEChemistryType>::sequential
+template<class ChemistryModel>
+Foam::sequential<ChemistryModel>::sequential
 (
-    const fvMesh& mesh,
-    const word& ODEModelName,
-    const word& thermoType
+    const fvMesh& mesh
 )
 :
-    chemistrySolver<ODEChemistryType>(mesh, ODEModelName, thermoType),
+    chemistrySolver<ChemistryModel>(mesh),
     coeffsDict_(this->subDict("sequentialCoeffs")),
     cTauChem_(readScalar(coeffsDict_.lookup("cTauChem"))),
     eqRateLimiter_(coeffsDict_.lookup("equilibriumRateLimiter"))
@@ -45,15 +43,15 @@ Foam::sequential<ODEChemistryType>::sequential
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-template<class ODEChemistryType>
-Foam::sequential<ODEChemistryType>::~sequential()
+template<class ChemistryModel>
+Foam::sequential<ChemistryModel>::~sequential()
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class ODEChemistryType>
-Foam::scalar Foam::sequential<ODEChemistryType>::solve
+template<class ChemistryModel>
+Foam::scalar Foam::sequential<ChemistryModel>::solve
 (
     scalarField &c,
     const scalar T,