Merge master, fixing conflicts

src/OpenFOAM/Make/files

@@ -1,5 +1,6 @@
 global/global.Cver
 global/dimensionedConstants/dimensionedConstants.C
+global/dimensionedConstants/constants/constants.C
 global/argList/argList.C
 global/clock/clock.C
 

src/OpenFOAM/db/IOstreams/Pstreams/IPstream.H

@@ -117,12 +117,6 @@ public:
                 const std::streamsize bufSize
             );
 
-            //- Non-blocking receives: wait until all have finished.
-            static void waitRequests();
-
-            //- Non-blocking receives: has request i finished?
-            static bool finishedRequest(const label i);
-
             //- Return next token from stream
             Istream& read(token&);
 

src/OpenFOAM/db/IOstreams/Pstreams/OPstream.H

@@ -115,12 +115,6 @@ public:
                 const std::streamsize bufSize
             );
 
-            //- Non-blocking writes: wait until all have finished.
-            static void waitRequests();
-
-            //- Non-blocking writes: has request i finished?
-            static bool finishedRequest(const label i);
-
             //- Write next token to stream
             Ostream& write(const token&);
 

src/OpenFOAM/db/IOstreams/Pstreams/Pstream.H

@@ -264,6 +264,12 @@ public:
         //  Spawns slave processes and initialises inter-communication
         static bool init(int& argc, char**& argv);
 
+        //- Non-blocking comms: wait until all have finished.
+        static void waitRequests();
+
+        //- Non-blocking comms: has request i finished?
+        static bool finishedRequest(const label i);
+
         //- Is this a parallel run?
         static bool parRun()
         {

src/OpenFOAM/db/Time/Time.C

@@ -115,7 +115,14 @@ void Foam::Time::setControls()
         {
             if (timeDirs.size())
             {
-                startTime_ = timeDirs[0].value();
+                if (timeDirs[0].name() == constant() && timeDirs.size() >= 2)
+                {
+                    startTime_ = timeDirs[1].value();
+                }
+                else
+                {
+                    startTime_ = timeDirs[0].value();
+                }
             }
         }
         else if (startFrom == "latestTime")

src/OpenFOAM/db/functionObjects/functionObjectList/functionObjectList.C

@@ -106,6 +106,20 @@ void Foam::functionObjectList::clear()
 }
 
 
+Foam::label Foam::functionObjectList::findObjectID(const word& name) const
+{
+    forAll(*this, objectI)
+    {
+        if (operator[](objectI).name() == name)
+        {
+            return objectI;
+        }
+    }
+
+    return -1;
+}
+
+
 void Foam::functionObjectList::on()
 {
     execution_ = true;
@@ -142,14 +156,9 @@ bool Foam::functionObjectList::execute()
             read();
         }
 
-        forAllIter
-        (
-            PtrList<functionObject>,
-            static_cast<PtrList<functionObject>&>(*this),
-            iter
-        )
+        forAll(*this, objectI)
         {
-            ok = iter().execute() && ok;
+            ok = operator[](objectI).execute() && ok;
         }
     }
 
@@ -168,14 +177,9 @@ bool Foam::functionObjectList::end()
             read();
         }
 
-        forAllIter
-        (
-            PtrList<functionObject>,
-            static_cast<PtrList<functionObject>&>(*this),
-            iter
-        )
+        forAll(*this, objectI)
         {
-            ok = iter().end() && ok;
+            ok = operator[](objectI).end() && ok;
         }
     }
 

src/OpenFOAM/db/functionObjects/functionObjectList/functionObjectList.H

@@ -120,9 +120,8 @@ public:
         );
 
 
-    // Destructor
-
-        virtual ~functionObjectList();
+    //- Destructor
+    virtual ~functionObjectList();
 
 
     // Member Functions
@@ -139,6 +138,8 @@ public:
         //- Clear the list of function objects
         virtual void clear();
 
+        //- Find the ID of a given function object by name
+        virtual label findObjectID(const word& name) const;
 
         //- Switch the function objects on
         virtual void on();
@@ -161,7 +162,6 @@ public:
 
         //- Read and set the function objects if their data have changed
         virtual bool read();
-
 };
 
 

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

@@ -300,8 +300,7 @@ evaluate()
         // Block for any outstanding requests
         if (Pstream::defaultCommsType == Pstream::nonBlocking)
         {
-            IPstream::waitRequests();
-            OPstream::waitRequests();
+            Pstream::waitRequests();
         }
 
         forAll(*this, patchi)

src/OpenFOAM/global/dimensionedConstants/constants/atomic/atomicConstants.C

@@ -0,0 +1,132 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "mathConstants.H"
+#include "universalConstants.H"
+#include "electromagneticConstants.H"
+#include "atomicConstants.H"
+
+#include "dimensionedConstants.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+const char* Foam::constant::atomic::group = "atomic";
+
+
+const Foam::dimensionedScalar Foam::constant::atomic::alpha
+(
+    dimensionedConstant
+    (
+        group,
+        "alpha",
+        dimensionedScalar
+        (
+            "alpha",
+            sqr(constant::electromagnetic::e)
+           /(
+                dimensionedScalar("C", dimless, 2.0)
+               *constant::electromagnetic::epsilon0
+               *constant::universal::h
+               *constant::universal::c
+            )
+        )
+    )
+);
+
+
+const Foam::dimensionedScalar Foam::constant::atomic::Rinf
+(
+    dimensionedConstant
+    (
+        group,
+        "Rinf",
+        dimensionedScalar
+        (
+            "Rinf",
+            sqr(alpha)*me*constant::universal::c
+           /(dimensionedScalar("C", dimless, 2.0)*constant::universal::h)
+        )
+    )
+);
+
+
+const Foam::dimensionedScalar Foam::constant::atomic::a0
+(
+    dimensionedConstant
+    (
+        group,
+        "a0",
+        dimensionedScalar
+        (
+            "a0",
+            alpha
+           /(dimensionedScalar("C", dimless, 4.0*constant::math::pi)*Rinf)
+        )
+    )
+);
+
+
+const Foam::dimensionedScalar Foam::constant::atomic::re
+(
+    dimensionedConstant
+    (
+        group,
+        "re",
+        dimensionedScalar
+        (
+            "re",
+            sqr(constant::electromagnetic::e)
+           /(
+                dimensionedScalar("C", dimless, 4.0*constant::math::pi)
+               *constant::electromagnetic::epsilon0
+               *me
+               *sqr(constant::universal::c)
+            )
+        )
+    )
+);
+
+
+const Foam::dimensionedScalar Foam::constant::atomic::Eh
+(
+    dimensionedConstant
+    (
+        group,
+        "Eh",
+        dimensionedScalar
+        (
+            "Eh",
+            dimensionedScalar("C", dimless, 2.0)
+           *Rinf*constant::universal::h*constant::universal::c
+        )
+    )
+);
+
+
+// ************************************************************************* //
+
+
+

src/OpenFOAM/global/dimensionedConstants/constants/atomic/atomicConstants.H

@@ -0,0 +1,81 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Namespace
+    Foam::constant::atom
+
+Description
+    Atomic constants
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef atomicConstants_H
+#define atomicConstants_H
+
+#include "dimensionedScalar.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace constant
+{
+namespace atomic
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    //- Group name for atomic constants
+    extern const char* group;
+
+    //- Fine-structure constant: default SI units: []
+    extern const dimensionedScalar alpha;
+
+    //- Rydberg constant: default SI units: [1/m]
+    extern const dimensionedScalar Rinf;
+
+    //- Bohr radius: default SI units: [m]
+    extern const dimensionedScalar a0;
+
+    //- Classical electron radius: default SI units: [m]
+    extern const dimensionedScalar re;
+
+    //- Hartree energy: default SI units: [J]
+    extern const dimensionedScalar Eh;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace atomic
+} // end namespace constant
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
+
+
+

src/OpenFOAM/global/dimensionedConstants/constants/constants.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2009-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,17 +22,21 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-\*---------------------------------------------------------------------------*/
+Description
+    Collection of dimensioned constants
 
-#include "createReactingCloudTypes.H"
-#include "BasicReactingMultiphaseCloud.H"
+\*---------------------------------------------------------------------------*/
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-namespace Foam
-{
-    createReactingCloudType(BasicReactingMultiphaseCloud);
-};
+// Constants supplied in the main controlDict
+#include "fundamentalConstants.C"
+
+// Derived constants
+#include "universalConstants.C"
+#include "electromagneticConstants.C"
+#include "atomicConstants.C"
+#include "physicoChemicalConstants.C"
 
 
 // ************************************************************************* //

src/OpenFOAM/global/dimensionedConstants/constants/constants.H

@@ -0,0 +1,59 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Namespace
+    Foam::constant
+
+Description
+    Collection of constants
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef constants_H
+#define constants_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Dimensionless coefficents
+
+// Mathematical constants
+#include "mathConstants.H"
+
+
+// Dimensioned constants
+
+// Fundamental constants
+#include "fundamentalConstants.H"
+
+// Derived constants
+#include "universalConstants.H"
+#include "electromagneticConstants.H"
+#include "atomicConstants.H"
+#include "physicoChemicalConstants.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/global/dimensionedConstants/constants/electromagnetic/electromagneticConstants.C

@@ -0,0 +1,164 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "mathConstants.H"
+#include "universalConstants.H"
+#include "electromagneticConstants.H"
+#include "atomicConstants.H"
+
+#include "dimensionedConstants.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+const char* Foam::constant::electromagnetic::group = "electromagnetic";
+
+
+const Foam::dimensionedScalar Foam::constant::electromagnetic::mu0
+(
+    dimensionedConstant
+    (
+        group,
+        "mu0",
+        dimensionedScalar
+        (
+            "mu0",
+            dimless,
+            4.0*constant::math::pi*1e-07
+        )
+    )
+);
+
+
+const Foam::dimensionedScalar Foam::constant::electromagnetic::epsilon0
+(
+    dimensionedConstant
+    (
+        group,
+        "epsilon0",
+        dimensionedScalar
+        (
+            "epsilon0",
+            dimensionedScalar("C", dimless, 1.0)
+           /(mu0*sqr(constant::universal::c))
+        )
+    )
+);
+
+
+const Foam::dimensionedScalar Foam::constant::electromagnetic::Z0
+(
+    dimensionedConstant
+    (
+        group,
+        "Z0",
+        dimensionedScalar
+        (
+            "Z0",
+            mu0*constant::universal::c
+        )
+    )
+);
+
+
+const Foam::dimensionedScalar Foam::constant::electromagnetic::kappa
+(
+    dimensionedConstant
+    (
+        group,
+        "kappa",
+        dimensionedScalar
+        (
+            "kappa",
+            dimensionedScalar("C", dimless, 1.0/(4.0*constant::math::pi))
+           /epsilon0
+        )
+    )
+);
+
+
+const Foam::dimensionedScalar Foam::constant::electromagnetic::G0
+(
+    dimensionedConstant
+    (
+        group,
+        "G0",
+        dimensionedScalar
+        (
+            "G0",
+            dimensionedScalar("C", dimless, 2)*sqr(e)/constant::universal::h
+        )
+    )
+);
+
+
+const Foam::dimensionedScalar Foam::constant::electromagnetic::KJ
+(
+    dimensionedConstant
+    (
+        group,
+        "KJ",
+        dimensionedScalar
+        (
+            "KJ",
+            dimensionedScalar("C", dimless, 2)*e/constant::universal::h
+        )
+    )
+);
+
+
+const Foam::dimensionedScalar Foam::constant::electromagnetic::phi0
+(
+    dimensionedConstant
+    (
+        group,
+        "phi0",
+        dimensionedScalar
+        (
+            "phi0",
+            constant::universal::h/(dimensionedScalar("C", dimless, 2)*e)
+        )
+    )
+);
+
+
+const Foam::dimensionedScalar Foam::constant::electromagnetic::RK
+(
+    dimensionedConstant
+    (
+        group,
+        "RK",
+        dimensionedScalar
+        (
+            "RK",
+            constant::universal::h/sqr(e)
+        )
+    )
+);
+
+
+// ************************************************************************* //
+
+

src/OpenFOAM/global/dimensionedConstants/constants/electromagnetic/electromagneticConstants.H

@@ -0,0 +1,89 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Namespace
+    Foam::constant::em
+
+Description
+    Electromagnetic constants
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef electromagneticConstants_H
+#define electromagneticConstants_H
+
+#include "dimensionedScalar.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace constant
+{
+namespace electromagnetic
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    //- Group name for electromagnetic constants
+    extern const char* group;
+
+    //- Magnetic constant/permeability of free space: default SI units: [H/m]
+    extern const dimensionedScalar mu0;
+
+    //- Electric constant: default SI units: [F/m]
+    extern const dimensionedScalar epsilon0;
+
+    //- Characteristic impedance of a vacuum: default SI units: [ohm]
+    extern const dimensionedScalar Z0;
+
+    //- Coulomb constant: default SI units: [N.m2/C2]
+    extern const dimensionedScalar kappa;
+
+    //- Conductance quantum: default SI units: [S]
+    extern const dimensionedScalar G0;
+
+    //- Josephson constant: default SI units: [Hz/V]
+    extern const dimensionedScalar KJ;
+
+    //- Magnetic flux quantum: default SI units: [Wb]
+    extern const dimensionedScalar phi0;
+
+    //- von Klitzing constant: default SI units: [ohm]
+    extern const dimensionedScalar RK;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace electromagnetic
+} // end namespace constant
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
+
+

src/OpenFOAM/global/dimensionedConstants/constants/fundamental/fundamentalConstants.C

@@ -0,0 +1,128 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+    Fundamental dimensioned constants
+
+\*---------------------------------------------------------------------------*/
+
+#include "fundamentalConstants.H"
+
+#include "universalConstants.H"
+#include "electromagneticConstants.H"
+#include "atomicConstants.H"
+#include "physicoChemicalConstants.H"
+
+#include "dimensionedConstants.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Universal constants
+
+const Foam::dimensionedScalar Foam::constant::universal::c
+(
+    dimensionedConstant(universal::group, "c")
+);
+
+
+const Foam::dimensionedScalar Foam::constant::universal::G
+(
+    dimensionedConstant(universal::group, "G")
+);
+
+
+const Foam::dimensionedScalar Foam::constant::universal::h
+(
+    dimensionedConstant(universal::group, "h")
+);
+
+
+// Electromagnetic
+
+const Foam::dimensionedScalar Foam::constant::electromagnetic::e
+(
+    dimensionedConstant(electromagnetic::group, "e")
+);
+
+
+// Atomic
+
+const Foam::dimensionedScalar Foam::constant::atomic::me
+(
+    dimensionedConstant(atomic::group, "me")
+);
+
+
+const Foam::dimensionedScalar Foam::constant::atomic::mp
+(
+    dimensionedConstant(atomic::group, "mp")
+);
+
+
+// Physico-chemical
+
+const Foam::dimensionedScalar Foam::constant::physicoChemical::mu
+(
+    dimensionedConstant(physicoChemical::group, "mu")
+);
+
+
+const Foam::dimensionedScalar Foam::constant::physicoChemical::NA
+(
+//    dimensionedConstant(physicoChemical::group, "NA")
+    dimensionedConstant
+    (
+        physicoChemical::group,
+        "NA",
+        dimensionedScalar
+        (
+            "NA",
+            dimless/dimMoles,
+            6.0221417930e+23
+        )
+    )
+);
+
+
+const Foam::dimensionedScalar Foam::constant::physicoChemical::k
+(
+    dimensionedConstant(physicoChemical::group, "k")
+);
+
+
+// Standard
+
+const Foam::dimensionedScalar Foam::constant::standard::Pstd
+(
+    dimensionedConstant("standard", "Pstd")
+);
+
+
+const Foam::dimensionedScalar Foam::constant::standard::Tstd
+(
+    dimensionedConstant("standard", "Tstd")
+);
+
+
+// ************************************************************************* //

src/OpenFOAM/global/dimensionedConstants/constants/fundamental/fundamentalConstants.H

@@ -0,0 +1,101 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+    Fundamental dimensioned constants
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef fundamentalConstants_H
+#define fundamentalConstants_H
+
+#include "dimensionedScalar.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace constant
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace universal
+{
+    //- Speed of light in a vacuum
+    extern const dimensionedScalar c;
+
+    //- Newtonian constant of gravitation
+    extern const dimensionedScalar G;
+
+    //- Planck constant
+    extern const dimensionedScalar h;
+}
+
+namespace electromagnetic
+{
+    //- Elementary charge
+    extern const dimensionedScalar e;
+}
+
+namespace atomic
+{
+    //- Electron mass
+    extern const dimensionedScalar me;
+
+    //- Proton mass
+    extern const dimensionedScalar mp;
+}
+
+namespace physicoChemical
+{
+    //- Atomic mass unit
+    extern const dimensionedScalar mu;
+
+    //- Avagadro number
+    extern const dimensionedScalar NA;
+
+    //- Boltzmann constant
+    extern const dimensionedScalar k;
+}
+
+namespace standard
+{
+    //- Standard pressure
+    extern const dimensionedScalar Pstd;
+
+    //- Standard temperature
+    extern const dimensionedScalar Tstd;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace constant
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/global/dimensionedConstants/constants/math/mathConstants.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2009-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,15 +23,15 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Namespace
-    Foam::mathematicalConstant
+    Foam::constant::math
 
 Description
-    Mathematical constants such as pi, e.
+    mathematical constants
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef mathematicalConstants_H
-#define mathematicalConstants_H
+#ifndef mathConstants_H
+#define mathConstants_H
 
 #include "scalar.H"
 
@@ -39,19 +39,24 @@ Description
 
 namespace Foam
 {
+namespace constant
+{
+namespace math
+{
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-namespace mathematicalConstant
-{
+    static word group = "math";
+
     const scalar e(M_E);
     const scalar pi(M_PI);
     const scalar twoPi(2*pi);
     const scalar piByTwo(0.5*pi);
-}
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+} // End namespace math
+} // end namespace constant
 } // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -59,3 +64,6 @@ namespace mathematicalConstant
 #endif
 
 // ************************************************************************* //
+
+
+

src/OpenFOAM/global/dimensionedConstants/constants/physicoChemical/physicoChemicalConstants.C

@@ -0,0 +1,132 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "mathConstants.H"
+#include "universalConstants.H"
+#include "electromagneticConstants.H"
+#include "physicoChemicalConstants.H"
+
+#include "dimensionedConstants.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+const char* Foam::constant::physicoChemical::group = "physicoChemical";
+
+
+const Foam::dimensionedScalar Foam::constant::physicoChemical::R
+(
+    dimensionedConstant
+    (
+        group,
+        "R",
+        dimensionedScalar
+        (
+            "R",
+            NA*k
+        )
+    )
+);
+
+
+const Foam::dimensionedScalar Foam::constant::physicoChemical::F
+(
+    dimensionedConstant
+    (
+        group,
+        "F",
+        dimensionedScalar
+        (
+            "F",
+            NA*constant::electromagnetic::e
+        )
+    )
+);
+
+
+const Foam::dimensionedScalar Foam::constant::physicoChemical::sigma
+(
+    dimensionedConstant
+    (
+        group,
+        "sigma",
+        dimensionedScalar
+        (
+            "sigma",
+            dimensionedScalar("C", dimless, sqr(constant::math::pi)/60.0)
+           *pow4(k)/(pow3(constant::universal::hr)*sqr(constant::universal::c))
+        )
+    )
+);
+
+
+const Foam::dimensionedScalar Foam::constant::physicoChemical::b
+(
+    dimensionedConstant
+    (
+        group,
+        "b",
+        dimensionedScalar
+        (
+            "b",
+            (constant::universal::h*constant::universal::c/k)
+           /dimensionedScalar("C", dimless, 4.965114231)
+        )
+    )
+);
+
+
+const Foam::dimensionedScalar Foam::constant::physicoChemical::c1
+(
+    dimensionedConstant
+    (
+        group,
+        "c1",
+        dimensionedScalar
+        (
+            "c1",
+            dimensionedScalar("C", dimless, constant::math::twoPi)
+           *constant::universal::h*sqr(constant::universal::c)
+        )
+    )
+);
+
+
+const Foam::dimensionedScalar Foam::constant::physicoChemical::c2
+(
+    dimensionedConstant
+    (
+        group,
+        "c2",
+        dimensionedScalar
+        (
+            "c2",
+            constant::universal::h*constant::universal::c/k
+        )
+    )
+);
+
+
+// ************************************************************************* //

src/OpenFOAM/global/dimensionedConstants/constants/physicoChemical/physicoChemicalConstants.H

@@ -0,0 +1,85 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Namespace
+    Foam::constant::phys
+
+Description
+    Physico-chemical constants
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef physicoChemicalConstants_H
+#define physicoChemicalConstants_H
+
+#include "dimensionedScalar.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace constant
+{
+namespace physicoChemical
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    //- Group name for physico-chemical constants
+    extern const char* group;
+
+    //- Universal gas constant: default SI units: [J/mol/K]
+    extern const dimensionedScalar R;
+
+    //- Faraday constant: default SI units: [C/mol]
+    extern const dimensionedScalar F;
+
+    //- Stefan-Boltzmann constant: default SI units: [W/m2/K4]
+    extern const dimensionedScalar sigma;
+
+    //- Wien displacement law constant: default SI units: [m.K]
+    extern const dimensionedScalar b;
+
+    //- First radiation constant: default SI units: [W/m2]
+    extern const dimensionedScalar c1;
+
+    //- Second radiation constant: default SI units: [m.K]
+    extern const dimensionedScalar c2;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace physicoChemical
+} // end namespace constant
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
+
+
+
+

src/OpenFOAM/global/dimensionedConstants/constants/universal/universalConstants.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2009-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,25 +24,30 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "radiationConstants.H"
+#include "universalConstants.H"
+#include "mathConstants.H"
+
 #include "dimensionedConstants.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-//- Stefan-Boltzmann constant (default in [J/(K4 m2 s)])
-const Foam::dimensionedScalar Foam::radiation::sigmaSB
+const char* Foam::constant::universal::group = "universal";
+
+
+const Foam::dimensionedScalar Foam::constant::universal::hr
 (
-    Foam::dimensionedConstant
+    dimensionedConstant
     (
-        "sigmaSB",
+        group,
+        "hr",
         dimensionedScalar
         (
-            "sigmaSB",
-            dimensionSet(1, 0, -3, -4, 0, 0, 0),
-            5.670E-08
+            "hr",
+            h/(dimensionedScalar("C", dimless, constant::math::twoPi))
         )
     )
 );
 
 
 // ************************************************************************* //
+

src/OpenFOAM/global/dimensionedConstants/constants/universal/universalConstants.H

@@ -0,0 +1,66 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Namespace
+    Foam::constant::uni
+
+Description
+    Universal constants
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef universalConstants_H
+#define universalConstants_H
+
+#include "dimensionedScalar.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace constant
+{
+namespace universal
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    //- Group name for universal constants
+    extern const char* group;
+
+    //- Reduced Planck constant: default SI units: [J/s]
+    extern const dimensionedScalar hr;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace universal
+} // End namespace constant
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
+

src/OpenFOAM/global/dimensionedConstants/dimensionedConstants.C

@@ -44,6 +44,40 @@ dictionary& dimensionedConstants()
     );
 }
 
+
+dimensionedScalar dimensionedConstant
+(
+    const word& group,
+    const word& varName
+)
+{
+    dictionary& dict = dimensionedConstants();
+
+    // Check that the entries exist.
+    // Note: should make FatalError robust instead!
+
+    if (!dict.found("unitSet"))
+    {
+        std::cerr<< "Cannot find unitSet in dictionary " << dict.name()
+            << std::endl;
+    }
+
+    const word unitSetCoeffs(word(dict.lookup("unitSet")) + "Coeffs");
+
+    if (!dict.found(unitSetCoeffs))
+    {
+        std::cerr<< "Cannot find " << unitSetCoeffs << " in dictionary "
+            << dict.name() << std::endl;
+    }
+
+    dictionary& unitDict = dict.subDict(unitSetCoeffs);
+
+    dictionary& groupDict = unitDict.subDict(group);
+
+    return dimensionedScalar(groupDict.lookup(varName));
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/global/dimensionedConstants/dimensionedConstants.H

@@ -40,6 +40,7 @@ SourceFiles
 #define dimensionedConstants_H
 
 #include "dictionary.H"
+#include "dimensionedScalar.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -50,20 +51,42 @@ namespace Foam
 
 dictionary& dimensionedConstants();
 
+
+dimensionedScalar dimensionedConstant(const word& group, const word& varName);
+
+
 template<class T>
 T dimensionedConstant
 (
-    const char* switchName,
-    const T defaultValue
+    const word& group,
+    const word& varName,
+    const T& defaultValue
 )
 {
-    if (dimensionedConstants().found(switchName))
+    dictionary& dict = dimensionedConstants();
+
+    const word unitSet(dict.lookup("unitSet"));
+
+    dictionary& unitDict(dict.subDict(unitSet + "Coeffs"));
+
+    if (unitDict.found(group))
     {
-        return pTraits<T>(dimensionedConstants().lookup(switchName));
+        dictionary& groupDict = unitDict.subDict(group);
+        if (groupDict.found(varName))
+        {
+            return pTraits<T>(groupDict.lookup(varName));
+        }
+        else
+        {
+            groupDict.add(varName, defaultValue);
+            return defaultValue;
+        }
     }
     else
     {
-        dimensionedConstants().add(switchName, defaultValue);
+        unitDict.add(group, dictionary::null);
+        unitDict.subDict(group).add(varName, defaultValue);
+
         return defaultValue;
     }
 }

src/OpenFOAM/meshes/meshShapes/face/face.C

@@ -27,7 +27,7 @@ License
 #include "face.H"
 #include "triFace.H"
 #include "triPointRef.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -99,13 +99,13 @@ Foam::label Foam::face::mostConcaveAngle
         if ((edgeNormal & n) > 0)
         {
             // Concave angle.
-            angle = mathematicalConstant::pi + edgeAngle;
+            angle = constant::math::pi + edgeAngle;
         }
         else
         {
             // Convex angle. Note '-' to take into account that rightEdge
             // and leftEdge are head-to-tail connected.
-            angle = mathematicalConstant::pi - edgeAngle;
+            angle = constant::math::pi - edgeAngle;
         }
 
         if (angle > maxAngle)
@@ -214,7 +214,7 @@ Foam::label Foam::face::split
         label index = fcIndex(fcIndex(startIndex));
 
         label minIndex = index;
-        scalar minDiff = Foam::mathematicalConstant::pi;
+        scalar minDiff = constant::math::pi;
 
         for(label i = 0; i < size() - 3; i++)
         {

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

@@ -339,8 +339,7 @@ void Foam::mapDistribute::compact(const boolList& elemIsUsed)
 
         // Wait for all to finish
 
-        OPstream::waitRequests();
-        IPstream::waitRequests();
+        Pstream::waitRequests();
 
 
         // Compact out all submap entries that are referring to unused elements

src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C

@@ -287,8 +287,7 @@ void Foam::mapDistribute::distribute
 
 
             // Wait till all finished
-            IPstream::waitRequests();
-            OPstream::waitRequests();
+            Pstream::waitRequests();
 
             // Consume
             for (label domain = 0; domain < Pstream::nProcs(); domain++)
@@ -413,8 +412,7 @@ void Foam::mapDistribute::distribute
 
             // Wait for all to finish
 
-            OPstream::waitRequests();
-            IPstream::waitRequests();
+            Pstream::waitRequests();
 
             // Collect neighbour fields
 
@@ -717,8 +715,7 @@ void Foam::mapDistribute::distribute
 
 
             // Wait till all finished
-            IPstream::waitRequests();
-            OPstream::waitRequests();
+            Pstream::waitRequests();
 
             // Consume
             for (label domain = 0; domain < Pstream::nProcs(); domain++)
@@ -842,8 +839,7 @@ void Foam::mapDistribute::distribute
 
             // Wait for all to finish
 
-            OPstream::waitRequests();
-            IPstream::waitRequests();
+            Pstream::waitRequests();
 
             // Collect neighbour fields
 

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

@@ -27,7 +27,7 @@ License
 #include "primitiveMesh.H"
 #include "pyramidPointFaceRef.H"
 #include "ListOps.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 #include "SortableList.H"
 
 
@@ -410,7 +410,7 @@ bool Foam::primitiveMesh::checkFaceOrthogonality
 
     // Severe nonorthogonality threshold
     const scalar severeNonorthogonalityThreshold =
-        ::cos(nonOrthThreshold_/180.0*mathematicalConstant::pi);
+        ::cos(nonOrthThreshold_/180.0*constant::math::pi);
 
     scalar minDDotS = GREAT;
 
@@ -472,9 +472,9 @@ bool Foam::primitiveMesh::checkFaceOrthogonality
             if (debug || report)
             {
                 Info<< "    Mesh non-orthogonality Max: "
-                    << ::acos(minDDotS)/mathematicalConstant::pi*180.0
+                    << ::acos(minDDotS)/constant::math::pi*180.0
                     << " average: " <<
-                    ::acos(sumDDotS/neiSize)/mathematicalConstant::pi*180.0
+                    ::acos(sumDDotS/neiSize)/constant::math::pi*180.0
                     << endl;
             }
         }
@@ -839,7 +839,7 @@ bool Foam::primitiveMesh::checkFaceAngles
             << exit(FatalError);
     }
 
-    const scalar maxSin = Foam::sin(maxDeg/180.0*mathematicalConstant::pi);
+    const scalar maxSin = Foam::sin(maxDeg/180.0*constant::math::pi);
 
     const pointField& p = points();
     const faceList& fcs = faces();
@@ -916,7 +916,7 @@ bool Foam::primitiveMesh::checkFaceAngles
     {
         scalar maxConcaveDegr =
             Foam::asin(Foam::min(1.0, maxEdgeSin))
-           *180.0/mathematicalConstant::pi;
+           *180.0/constant::math::pi;
 
         if (debug || report)
         {

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

@@ -31,7 +31,7 @@ Description
 #include "primitiveMesh.H"
 #include "pyramidPointFaceRef.H"
 #include "cell.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -200,7 +200,7 @@ bool Foam::primitiveMesh::checkMeshMotion
                 )   << "Severe non-orthogonality in mesh motion for face "
                     << faceI
                     << " between cells " << own[faceI] << " and " << nei[faceI]
-                    << ": Angle = " << ::acos(dDotS)/mathematicalConstant::pi*180.0
+                    << ": Angle = " << ::acos(dDotS)/constant::math::pi*180.0
                     << " deg." << endl;
 
                 nDotProductErrors++;

src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H

@@ -26,7 +26,7 @@ License
 
 #include "IOstreams.H"
 #include "pointHit.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -307,9 +307,9 @@ inline scalar triangle<Point, PointRef>::quality() const
     return
         mag()
       / (
-            mathematicalConstant::pi
-          * Foam::sqr(circumRadius())
-          * 0.413497
+            constant::math::pi
+           *Foam::sqr(circumRadius())
+           *0.413497
           + VSMALL
         );
 }

src/OpenFOAM/primitives/SymmTensor/labelSymmTensor/labelSymmTensor.C

@@ -25,7 +25,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "labelSymmTensor.H"
-#include "mathematicalConstants.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/OpenFOAM/primitives/Tensor/tensor/tensor.C

@@ -25,7 +25,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "tensor.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -150,10 +150,8 @@ vector eigenValues(const tensor& t)
                 scalar aBy3 = a/3;
 
                 i = m2SqrtQ*cos(theta/3) - aBy3;
-                ii = m2SqrtQ*cos((theta + mathematicalConstant::twoPi)/3)
-                    - aBy3;
-                iii = m2SqrtQ*cos((theta - mathematicalConstant::twoPi)/3)
-                    - aBy3;
+                ii = m2SqrtQ*cos((theta + constant::math::twoPi)/3) - aBy3;
+                iii = m2SqrtQ*cos((theta - constant::math::twoPi)/3) - aBy3;
             }
             else
             {
@@ -345,10 +343,8 @@ vector eigenValues(const symmTensor& t)
                 scalar aBy3 = a/3;
 
                 i = m2SqrtQ*cos(theta/3) - aBy3;
-                ii = m2SqrtQ*cos((theta + mathematicalConstant::twoPi)/3)
-                    - aBy3;
-                iii = m2SqrtQ*cos((theta - mathematicalConstant::twoPi)/3)
-                    - aBy3;
+                ii = m2SqrtQ*cos((theta + constant::math::twoPi)/3) - aBy3;
+                iii = m2SqrtQ*cos((theta - constant::math::twoPi)/3) - aBy3;
             }
             else
             {

src/OpenFOAM/primitives/Tensor2D/tensor2D/tensor2D.C

@@ -25,7 +25,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "tensor2D.H"
-#include "mathematicalConstants.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/OpenFOAM/primitives/transform/transform.H

@@ -34,7 +34,7 @@ Description
 #define transform_H
 
 #include "tensor.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -104,7 +104,7 @@ inline Tensor<Cmpt> transform(const tensor& tt, const Tensor<Cmpt>& t)
         (tt.yx()*t.xx() + tt.yy()*t.yx() + tt.yz()*t.zx())*tt.zx()
       + (tt.yx()*t.xy() + tt.yy()*t.yy() + tt.yz()*t.zy())*tt.zy()
       + (tt.yx()*t.xz() + tt.yy()*t.yz() + tt.yz()*t.zz())*tt.zz(),
-        
+
         (tt.zx()*t.xx() + tt.zy()*t.yx() + tt.zz()*t.zx())*tt.xx()
       + (tt.zx()*t.xy() + tt.zy()*t.yy() + tt.zz()*t.zy())*tt.xy()
       + (tt.zx()*t.xz() + tt.zy()*t.yz() + tt.zz()*t.zz())*tt.xz(),
@@ -155,7 +155,7 @@ inline SymmTensor<Cmpt> transform(const tensor& tt, const SymmTensor<Cmpt>& st)
         (tt.yx()*st.xx() + tt.yy()*st.xy() + tt.yz()*st.xz())*tt.zx()
       + (tt.yx()*st.xy() + tt.yy()*st.yy() + tt.yz()*st.yz())*tt.zy()
       + (tt.yx()*st.xz() + tt.yy()*st.yz() + tt.yz()*st.zz())*tt.zz(),
-        
+
         (tt.zx()*st.xx() + tt.zy()*st.xy() + tt.zz()*st.xz())*tt.zx()
       + (tt.zx()*st.xy() + tt.zy()*st.yy() + tt.zz()*st.yz())*tt.zy()
       + (tt.zx()*st.xz() + tt.zy()*st.yz() + tt.zz()*st.zz())*tt.zz()
@@ -205,11 +205,11 @@ inline scalar pseudoAngle
 
     if (sin < -SMALL)
     {
-        return (3.0 + cos)*mathematicalConstant::piByTwo;
+        return (3.0 + cos)*constant::math::piByTwo;
     }
     else
     {
-        return (1.0 - cos)*mathematicalConstant::piByTwo;
+        return (1.0 - cos)*constant::math::piByTwo;
     }
 }
 

src/Pstream/dummy/IPread.C

@@ -86,17 +86,4 @@ int Foam::IPstream::read
 }
 
 
-void Foam::IPstream::waitRequests()
-{}
-
-
-bool Foam::IPstream::finishedRequest(const label)
-{
-    notImplemented("IPstream::finishedRequest()");
-    return false;
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 // ************************************************************************* //

src/Pstream/dummy/OPwrite.C

@@ -65,17 +65,4 @@ bool Foam::OPstream::write
 }
 
 
-void Foam::OPstream::waitRequests()
-{}
-
-
-bool Foam::OPstream::finishedRequest(const label)
-{
-    notImplemented("OPstream::finishedRequest()");
-    return false;
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 // ************************************************************************* //

src/Pstream/dummy/Pstream.C

@@ -61,6 +61,17 @@ void Foam::Pstream::abort()
 void Foam::reduce(scalar&, const sumOp<scalar>&)
 {}
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+void Foam::Pstream::waitRequests()
+{}
+
+
+bool Foam::Pstream::finishedRequest(const label i)
+{
+    notImplemented("Pstream::finishedRequest()");
+    return false;
+}
+
 
 // ************************************************************************* //

src/Pstream/mpi/IPread.C

@@ -30,6 +30,7 @@ Description
 #include "mpi.h"
 
 #include "IPstream.H"
+#include "PstreamGlobals.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -37,7 +38,7 @@ Description
 
 // Outstanding non-blocking operations.
 //! @cond fileScope
-Foam::DynamicList<MPI_Request> IPstream_outstandingRequests_;
+//Foam::DynamicList<MPI_Request> IPstream_outstandingRequests_;
 //! @endcond fileScope
 
 // * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * * //
@@ -185,7 +186,7 @@ Foam::label Foam::IPstream::read
             return 0;
         }
 
-        IPstream_outstandingRequests_.append(request);
+        PstreamGlobals::outstandingRequests_.append(request);
 
         // Assume the message is completely received.
         return bufSize;
@@ -204,52 +205,6 @@ Foam::label Foam::IPstream::read
 }
 
 
-void Foam::IPstream::waitRequests()
-{
-    if (IPstream_outstandingRequests_.size())
-    {
-        if
-        (
-            MPI_Waitall
-            (
-                IPstream_outstandingRequests_.size(),
-                IPstream_outstandingRequests_.begin(),
-                MPI_STATUSES_IGNORE
-            )
-        )
-        {
-            FatalErrorIn
-            (
-                "IPstream::waitRequests()"
-            )   << "MPI_Waitall returned with error" << endl;
-        }
-
-        IPstream_outstandingRequests_.clear();
-    }
-}
-
-
-bool Foam::IPstream::finishedRequest(const label i)
-{
-    if (i >= IPstream_outstandingRequests_.size())
-    {
-        FatalErrorIn
-        (
-            "IPstream::finishedRequest(const label)"
-        )   << "There are " << IPstream_outstandingRequests_.size()
-            << " outstanding send requests and you are asking for i=" << i
-            << nl
-            << "Maybe you are mixing blocking/non-blocking comms?"
-            << Foam::abort(FatalError);
-    }
-
-    int flag;
-    MPI_Test(&IPstream_outstandingRequests_[i], &flag, MPI_STATUS_IGNORE);
-
-    return flag != 0;
-}
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 // ************************************************************************* //

src/Pstream/mpi/Make/files

@@ -1,5 +1,6 @@
 OPwrite.C
 IPread.C
 Pstream.C
+PstreamGlobals.C
 
 LIB = $(FOAM_MPI_LIBBIN)/libPstream

src/Pstream/mpi/OPwrite.C

@@ -30,13 +30,7 @@ Description
 #include "mpi.h"
 
 #include "OPstream.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-// Outstanding non-blocking operations.
-//! @cond fileScope
-Foam::DynamicList<MPI_Request> OPstream_outstandingRequests_;
-//! @endcond fileScope
+#include "PstreamGlobals.H"
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
@@ -126,7 +120,7 @@ bool Foam::OPstream::write
             &request
         );
 
-        OPstream_outstandingRequests_.append(request);
+        PstreamGlobals::outstandingRequests_.append(request);
     }
     else
     {
@@ -142,52 +136,6 @@ bool Foam::OPstream::write
 }
 
 
-void Foam::OPstream::waitRequests()
-{
-    if (OPstream_outstandingRequests_.size())
-    {
-        if
-        (
-            MPI_Waitall
-            (
-                OPstream_outstandingRequests_.size(),
-                OPstream_outstandingRequests_.begin(),
-                MPI_STATUSES_IGNORE
-            )
-        )
-        {
-            FatalErrorIn
-            (
-                "OPstream::waitRequests()"
-            )   << "MPI_Waitall returned with error" << Foam::endl;
-        }
-
-        OPstream_outstandingRequests_.clear();
-    }
-}
-
-
-bool Foam::OPstream::finishedRequest(const label i)
-{
-    if (i >= OPstream_outstandingRequests_.size())
-    {
-        FatalErrorIn
-        (
-            "OPstream::finishedRequest(const label)"
-        )   << "There are " << OPstream_outstandingRequests_.size()
-            << " outstanding send requests and you are asking for i=" << i
-            << nl
-            << "Maybe you are mixing blocking/non-blocking comms?"
-            << Foam::abort(FatalError);
-    }
-
-    int flag;
-    MPI_Test(&OPstream_outstandingRequests_[i], &flag, MPI_STATUS_IGNORE);
-
-    return flag != 0;
-}
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 // ************************************************************************* //

src/Pstream/mpi/Pstream.C

@@ -29,6 +29,7 @@ License
 #include "Pstream.H"
 #include "PstreamReduceOps.H"
 #include "OSspecific.H"
+#include "PstreamGlobals.H"
 
 #include <cstring>
 #include <cstdlib>
@@ -130,6 +131,19 @@ void Foam::Pstream::exit(int errnum)
     delete[] buff;
 #   endif
 
+    if (PstreamGlobals::outstandingRequests_.size())
+    {
+        label n = PstreamGlobals::outstandingRequests_.size();
+        PstreamGlobals::outstandingRequests_.clear();
+
+        WarningIn("Pstream::exit(int)")
+            << "There are still " << n << " outstanding MPI_Requests." << endl
+            << "This means that your code exited before doing a"
+            << " Pstream::waitRequests()." << endl
+            << "This should not happen for a normal code exit."
+            << endl;
+    }
+
     if (errnum == 0)
     {
         MPI_Finalize();
@@ -422,6 +436,57 @@ void Foam::reduce(scalar& Value, const sumOp<scalar>& bop)
 }
 
 
+void Foam::Pstream::waitRequests()
+{
+    if (PstreamGlobals::outstandingRequests_.size())
+    {
+        if
+        (
+            MPI_Waitall
+            (
+                PstreamGlobals::outstandingRequests_.size(),
+                PstreamGlobals::outstandingRequests_.begin(),
+                MPI_STATUSES_IGNORE
+            )
+        )
+        {
+            FatalErrorIn
+            (
+                "Pstream::waitRequests()"
+            )   << "MPI_Waitall returned with error" << Foam::endl;
+        }
+
+        PstreamGlobals::outstandingRequests_.clear();
+    }
+}
+
+
+bool Foam::Pstream::finishedRequest(const label i)
+{
+    if (i >= PstreamGlobals::outstandingRequests_.size())
+    {
+        FatalErrorIn
+        (
+            "Pstream::finishedRequest(const label)"
+        )   << "There are " << PstreamGlobals::outstandingRequests_.size()
+            << " outstanding send requests and you are asking for i=" << i
+            << nl
+            << "Maybe you are mixing blocking/non-blocking comms?"
+            << Foam::abort(FatalError);
+    }
+
+    int flag;
+    MPI_Test
+    (
+       &PstreamGlobals::outstandingRequests_[i],
+       &flag,
+        MPI_STATUS_IGNORE
+    );
+
+    return flag != 0;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 // ************************************************************************* //

src/Pstream/mpi/PstreamGlobals.C

@@ -22,32 +22,24 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-InNamespace
-    Foam::radiation
-
-Description
-    Constants used in radiation modelling
-
 \*---------------------------------------------------------------------------*/
 
-#ifndef radiationConstants_H
-#define radiationConstants_H
-
-#include "dimensionedScalar.H"
+#include "PstreamGlobals.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
-namespace radiation
-{
-    //- Stefan-Boltzmann constant [J/(K4 m2 s)]
-    extern const dimensionedScalar sigmaSB;
-}
-}
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+// Outstanding non-blocking operations.
+//! @cond fileScope
+DynamicList<MPI_Request> PstreamGlobals::outstandingRequests_;
+//! @endcond fileScope
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#endif
+} // End namespace Foam
 
 // ************************************************************************* //

src/Pstream/mpi/PstreamGlobals.H

@@ -0,0 +1,69 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Namespace
+    Foam::PstreamGlobals
+
+Description
+    Global functions and variables for working with parallel streams,
+    but principally for gamma/mpi
+
+SourceFiles
+    PstreamGlobals.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PstreamGlobals_H
+#define PstreamGlobals_H
+
+#include "mpi.h"
+
+#include "DynamicList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class PstreamGlobals Declaration
+\*---------------------------------------------------------------------------*/
+
+namespace PstreamGlobals
+{
+
+extern DynamicList<MPI_Request> outstandingRequests_;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C

@@ -34,7 +34,7 @@ Description
 #include "removePoints.H"
 #include "pointFields.H"
 #include "motionSmoother.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 #include "pointSet.H"
 #include "faceSet.H"
 #include "cellSet.H"
@@ -2479,17 +2479,11 @@ void Foam::autoLayerDriver::mergePatchFacesUndo
     const dictionary& motionDict
 )
 {
-    scalar minCos = Foam::cos
-    (
-        layerParams.featureAngle()
-      * mathematicalConstant::pi/180.0
-    );
+    scalar minCos =
+        Foam::cos(layerParams.featureAngle()*constant::math::pi/180.0);
 
-    scalar concaveCos = Foam::cos
-    (
-        layerParams.concaveAngle()
-      * mathematicalConstant::pi/180.0
-    );
+    scalar concaveCos =
+        Foam::cos(layerParams.concaveAngle()*constant::math::pi/180.0);
 
     Info<< nl
         << "Merging all faces of a cell" << nl
@@ -2588,7 +2582,7 @@ void Foam::autoLayerDriver::addLayers
     (
         pp,
         meshEdges,
-        layerParams.featureAngle()*mathematicalConstant::pi/180.0,
+        layerParams.featureAngle()*constant::math::pi/180.0,
 
         patchDisp,
         patchNLayers,
@@ -2987,6 +2981,7 @@ void Foam::autoLayerDriver::addLayers
         (
             invExpansionRatio,
             pp,
+            labelList(0),       // exposed patchIDs, not used for adding layers
             nPatchFaceLayers,   // layers per face
             nPatchPointLayers,  // layers per point
             firstDisp,          // thickness of layer nearest internal mesh

src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C

@@ -35,6 +35,7 @@ License
 #include "refinementSurfaces.H"
 #include "shellSurfaces.H"
 #include "mapDistributePolyMesh.H"
+#include "mathConstants.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -678,8 +679,8 @@ void Foam::autoRefineDriver::mergePatchFaces
 
     meshRefiner_.mergePatchFaces
     (
-        Foam::cos(45*mathematicalConstant::pi/180.0),
-        Foam::cos(45*mathematicalConstant::pi/180.0),
+        Foam::cos(45*constant::math::pi/180.0),
+        Foam::cos(45*constant::math::pi/180.0),
         meshRefiner_.meshedPatches()
     );
 
@@ -688,7 +689,7 @@ void Foam::autoRefineDriver::mergePatchFaces
         meshRefiner_.checkData();
     }
 
-    meshRefiner_.mergeEdges(Foam::cos(45*mathematicalConstant::pi/180.0));
+    meshRefiner_.mergeEdges(Foam::cos(45*constant::math::pi/180.0));
 
     if (debug)
     {

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

@@ -26,7 +26,7 @@ License
 
 #include "layerParameters.H"
 #include "polyBoundaryMesh.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 #include "refinementSurfaces.H"
 #include "searchableSurfaces.H"
 #include "regExp.H"
@@ -192,12 +192,7 @@ Foam::layerParameters::layerParameters
     ),
     layerTerminationCos_
     (
-        Foam::cos
-        (
-            0.5
-          * featureAngle_
-          * mathematicalConstant::pi/180.
-        )
+        Foam::cos(0.5*featureAngle_*constant::math::pi/180.0)
     ),
     maxThicknessToMedialRatio_
     (
@@ -206,7 +201,7 @@ Foam::layerParameters::layerParameters
     minMedianAxisAngleCos_
     (
         Foam::cos(readScalar(dict.lookup("minMedianAxisAngle")))
-      * mathematicalConstant::pi/180.
+       *constant::math::pi/180.0
     ),
     nBufferCellsNoExtrude_
     (
@@ -274,12 +269,7 @@ Foam::layerParameters::layerParameters
     ),
     layerTerminationCos_
     (
-        Foam::cos
-        (
-            0.5
-          * featureAngle_
-          * mathematicalConstant::pi/180.
-        )
+        Foam::cos(0.5*featureAngle_*constant::math::pi/180.0)
     ),
     maxThicknessToMedialRatio_
     (
@@ -288,7 +278,7 @@ Foam::layerParameters::layerParameters
     minMedianAxisAngleCos_
     (
         Foam::cos(readScalar(dict.lookup("minMedianAxisAngle")))
-      * mathematicalConstant::pi/180.
+       *constant::math::pi/180.0
     ),
     nBufferCellsNoExtrude_
     (

src/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C

@@ -25,7 +25,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "refinementParameters.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 #include "polyMesh.H"
 #include "globalIndex.H"
 
@@ -63,7 +63,7 @@ Foam::refinementParameters::refinementParameters(const dictionary& dict)
     }
     else
     {
-        curvature_ = Foam::cos(featAngle*mathematicalConstant::pi/180.0);
+        curvature_ = Foam::cos(featAngle*constant::math::pi/180.0);
     }
 }
 

src/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C

@@ -37,6 +37,7 @@ License
 #include "searchableSurfaces.H"
 #include "polyMeshGeometry.H"
 #include "IOmanip.H"
+#include "mathConstants.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -246,10 +247,7 @@ Foam::Map<Foam::label> Foam::meshRefinement::findEdgeConnectedProblemCells
             nearestRegion[i]
         );
 
-        scalar angle =
-            perpendicularAngle[region]
-          / 180.0
-          * mathematicalConstant::pi;
+        scalar angle = perpendicularAngle[region]/180.0*constant::math::pi;
 
         if (angle >= 0)
         {
@@ -306,7 +304,7 @@ bool Foam::meshRefinement::isCollapsedFace
         vector d = ownCc - mesh_.cellCentres()[nei];
 
         scalar dDotS = (d & s)/(mag(d)*magS + VSMALL);
-        
+
         if (dDotS < maxNonOrtho)
         {
             return true;

src/decompositionMethods/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.C

@@ -198,9 +198,7 @@ void Foam::hierarchGeomDecomp::findBinary
     label high = values.size();
 
     // Safeguards to avoid infinite loop.
-    label lowPrev = -1;
-    label midPrev = -1;
-    label highPrev = -1;
+    scalar midValuePrev = VGREAT;
 
     while (true)
     {
@@ -208,10 +206,10 @@ void Foam::hierarchGeomDecomp::findBinary
 
         if (debug)
         {
-            Pout<< "low:" << low << " lowValue:" << lowValue
+            Pout<< "    low:" << low << " lowValue:" << lowValue
                 << " high:" << high << " highValue:" << highValue
-                << " mid:" << mid << " midValue:" << midValue << nl
-                << "globalSize:" << size << " wantedSize:" << wantedSize
+                << " mid:" << mid << " midValue:" << midValue << endl
+                << "    globalSize:" << size << " wantedSize:" << wantedSize
                 << " sizeTol:" << sizeTol << endl;
         }
 
@@ -235,10 +233,7 @@ void Foam::hierarchGeomDecomp::findBinary
         mid = findLower(values, midValue, low, high);
 
         // Safeguard if same as previous.
-        bool hasNotChanged =
-            (mid == midPrev)
-         && (low == lowPrev)
-         && (high == highPrev);
+        bool hasNotChanged = (mag(midValue-midValuePrev) < SMALL);
 
         if (returnReduce(hasNotChanged, andOp<bool>()))
         {
@@ -248,9 +243,7 @@ void Foam::hierarchGeomDecomp::findBinary
             break;
         }
 
-        midPrev = mid;
-        lowPrev = low;
-        highPrev = high;
+        midValuePrev = midValue;
     }
 }
 
@@ -280,9 +273,7 @@ void Foam::hierarchGeomDecomp::findBinary
     label high = values.size();
 
     // Safeguards to avoid infinite loop.
-    label lowPrev = -1;
-    label midPrev = -1;
-    label highPrev = -1;
+    scalar midValuePrev = VGREAT;
 
     while (true)
     {
@@ -294,10 +285,10 @@ void Foam::hierarchGeomDecomp::findBinary
 
         if (debug)
         {
-            Pout<< "low:" << low << " lowValue:" << lowValue
+            Pout<< "    low:" << low << " lowValue:" << lowValue
                 << " high:" << high << " highValue:" << highValue
-                << " mid:" << mid << " midValue:" << midValue << nl
-                << "globalSize:" << weightedSize 
+                << " mid:" << mid << " midValue:" << midValue << endl
+                << "    globalSize:" << weightedSize 
                 << " wantedSize:" << wantedSize
                 << " sizeTol:" << sizeTol << endl;
         }
@@ -322,10 +313,7 @@ void Foam::hierarchGeomDecomp::findBinary
         mid = findLower(values, midValue, low, high);
 
         // Safeguard if same as previous.
-        bool hasNotChanged =
-            (mid == midPrev)
-         && (low == lowPrev)
-         && (high == highPrev);
+        bool hasNotChanged = (mag(midValue-midValuePrev) < SMALL);
 
         if (returnReduce(hasNotChanged, andOp<bool>()))
         {
@@ -335,9 +323,7 @@ void Foam::hierarchGeomDecomp::findBinary
             break;
         }
 
-        midPrev = mid;
-        lowPrev = low;
-        highPrev = high;
+        midValuePrev = midValue;
     }
 }
 
@@ -460,11 +446,11 @@ void Foam::hierarchGeomDecomp::sortComponent
         if (debug)
         {
             Pout<< "For component " << compI << ", bin " << bin
-                << " copying" << nl
+                << " copying" << endl
                 << "from " << leftCoord << " at local index "
-                << leftIndex << nl
+                << leftIndex << endl
                 << "to " << rightCoord << " localSize:"
-                << localSize << nl
+                << localSize << endl
                 << endl;
         }
 
@@ -643,11 +629,11 @@ void Foam::hierarchGeomDecomp::sortComponent
         if (debug)
         {
             Pout<< "For component " << compI << ", bin " << bin
-                << " copying" << nl
+                << " copying" << endl
                 << "from " << leftCoord << " at local index "
-                << leftIndex << nl
+                << leftIndex << endl
                 << "to " << rightCoord << " localSize:"
-                << localSize << nl
+                << localSize << endl
                 << endl;
         }
 
@@ -788,7 +774,6 @@ Foam::labelList Foam::hierarchGeomDecomp::decompose
 
     pointField rotatedPoints = rotDelta_ & points;
 
-
     // Calculate tolerance of cell distribution. For large cases finding
     // distibution to the cell exact would cause too many iterations so allow
     // some slack.

src/dynamicFvMesh/dynamicInkJetFvMesh/dynamicInkJetFvMesh.C

@@ -27,7 +27,7 @@ License
 #include "dynamicInkJetFvMesh.H"
 #include "addToRunTimeSelectionTable.H"
 #include "volFields.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -90,7 +90,7 @@ Foam::dynamicInkJetFvMesh::~dynamicInkJetFvMesh()
 bool Foam::dynamicInkJetFvMesh::update()
 {
     scalar scalingFunction =
-        0.5*(::cos(2*mathematicalConstant::pi*frequency_*time().value()) - 1.0);
+        0.5*(::cos(constant::math::twoPi*frequency_*time().value()) - 1.0);
 
     Info<< "Mesh scaling. Time = " << time().value() << " scaling: "
         << scalingFunction << endl;
@@ -113,7 +113,7 @@ bool Foam::dynamicInkJetFvMesh::update()
 
     fvMesh::movePoints(newPoints);
 
-    volVectorField& U = 
+    volVectorField& U =
         const_cast<volVectorField&>(lookupObject<volVectorField>("U"));
     U.correctBoundaryConditions();
 

src/dynamicFvMesh/solidBodyMotionFvMesh/solidBodyMotionFunctions/SDA/SDA.C

@@ -26,9 +26,9 @@ License
 
 #include "SDA.H"
 #include "addToRunTimeSelectionTable.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 
-using namespace Foam::mathematicalConstant;
+using namespace Foam::constant::math;
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -70,18 +70,18 @@ Foam::septernion Foam::solidBodyMotionFunctions::SDA::transformation() const
     scalar time = time_.value();
 
     scalar Tpi = Tp_ + dTp_*(time/dTi_);   // Current roll period [sec]
-    scalar wr = 2*pi/Tpi; // Current Freq [/sec]
+    scalar wr = twoPi/Tpi; // Current Freq [/sec]
 
     // Current Phase for roll [rad]
     scalar r = dTp_/dTi_;
     scalar u = Tp_ + r*time;
-    scalar phr = 2*pi*((Tp_/u - 1) + log(mag(u)) - log(Tp_))/r;
+    scalar phr = twoPi*((Tp_/u - 1) + log(mag(u)) - log(Tp_))/r;
 
     // Current Phase for Sway [rad]
     scalar phs = phr + pi;
 
     // Current Phase for Heave [rad]
-    scalar phh = phr + pi/2;
+    scalar phh = phr + piByTwo;
 
     scalar rollA = max(rollAmax_*exp(-sqr(Tpi - Tpn_)/(2*Q_)), rollAmin_);
 

src/dynamicFvMesh/solidBodyMotionFvMesh/solidBodyMotionFunctions/SKA/SKA.C

@@ -29,9 +29,9 @@ License
 #include "Tuple2.H"
 #include "IFstream.H"
 #include "interpolateXY.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 
-using namespace Foam::mathematicalConstant;
+using namespace Foam::constant::math;
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 

src/dynamicFvMesh/solidBodyMotionFvMesh/solidBodyMotionFvMesh.C

@@ -27,11 +27,8 @@ License
 #include "solidBodyMotionFvMesh.H"
 #include "addToRunTimeSelectionTable.H"
 #include "volFields.H"
-#include "mathematicalConstants.H"
 #include "transformField.H"
 
-using namespace Foam::mathematicalConstant;
-
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam

src/dynamicMesh/meshCut/cellLooper/hexCellLooper.C

@@ -28,7 +28,6 @@ License
 #include "cellFeatures.H"
 #include "polyMesh.H"
 #include "cellModeller.H"
-#include "mathematicalConstants.H"
 #include "plane.H"
 #include "ListOps.H"
 #include "meshTools.H"
@@ -36,18 +35,12 @@ License
 
 #include "addToRunTimeSelectionTable.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
 defineTypeNameAndDebug(hexCellLooper, 0);
-
 addToRunTimeSelectionTable(cellLooper, hexCellLooper, word);
-
-
 }
 
 

src/dynamicMesh/meshCut/cellLooper/topoCellLooper.C

@@ -27,7 +27,7 @@ License
 #include "topoCellLooper.H"
 #include "cellFeatures.H"
 #include "polyMesh.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 #include "DynamicList.H"
 #include "ListOps.H"
 #include "meshTools.H"
@@ -35,9 +35,8 @@ License
 
 #include "addToRunTimeSelectionTable.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
 namespace Foam
 {
    defineTypeNameAndDebug(topoCellLooper, 0);
@@ -46,7 +45,7 @@ namespace Foam
 
 // Angle for polys to be considered splitHexes.
 const Foam::scalar Foam::topoCellLooper::featureCos =
-    Foam::cos(10.0 * mathematicalConstant::pi/180.0);
+    Foam::cos(10.0*constant::math::pi/180.0);
 
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

src/dynamicMesh/meshCut/meshModifiers/undoableMeshCutter/undoableMeshCutter.C

@@ -32,16 +32,14 @@ License
 #include "cellCuts.H"
 #include "splitCell.H"
 #include "mapPolyMesh.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 #include "meshTools.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-
 defineTypeNameAndDebug(undoableMeshCutter, 0);
-
 }
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
@@ -193,8 +191,8 @@ Foam::undoableMeshCutter::undoableMeshCutter
     liveSplitCells_(mesh.nCells()/100 + 100),
     faceRemover_
     (
-        mesh,   
-        Foam::cos(30./180. * mathematicalConstant::pi)
+        mesh,
+        Foam::cos(30.0/180.0*constant::math::pi)
     )
 {}
 

src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C

@@ -27,6 +27,7 @@ License
 #include "polyMeshGeometry.H"
 #include "pyramidPointFaceRef.H"
 #include "syncTools.H"
+#include "mathConstants.H"
 
 namespace Foam
 {
@@ -247,7 +248,7 @@ Foam::scalar Foam::polyMeshGeometry::checkNonOrtho
                     << " between cells " << mesh.faceOwner()[faceI]
                     << " and " << nei
                     << ": Angle = "
-                    << ::acos(dDotS)/mathematicalConstant::pi*180.0
+                    << ::acos(dDotS)/constant::math::pi*180.0
                     << " deg." << endl;
             }
 
@@ -268,7 +269,7 @@ Foam::scalar Foam::polyMeshGeometry::checkNonOrtho
                     << " between cells " << mesh.faceOwner()[faceI]
                     << " and " << nei
                     << ": Angle = "
-                    << ::acos(dDotS)/mathematicalConstant::pi*180.0
+                    << ::acos(dDotS)/constant::math::pi*180.0
                     << " deg." << endl;
             }
 
@@ -368,7 +369,7 @@ bool Foam::polyMeshGeometry::checkFaceDotProduct
 
     // Severe nonorthogonality threshold
     const scalar severeNonorthogonalityThreshold =
-        ::cos(orthWarn/180.0*mathematicalConstant::pi);
+        ::cos(orthWarn/180.0*constant::math::pi);
 
 
     // Calculate coupled cell centre
@@ -503,9 +504,9 @@ bool Foam::polyMeshGeometry::checkFaceDotProduct
         if (nDDotS > 0)
         {
             Info<< "Mesh non-orthogonality Max: "
-                << ::acos(minDDotS)/mathematicalConstant::pi*180.0
+                << ::acos(minDDotS)/constant::math::pi*180.0
                 << " average: " <<
-                   ::acos(sumDDotS/nDDotS)/mathematicalConstant::pi*180.0
+                   ::acos(sumDDotS/nDDotS)/constant::math::pi*180.0
                 << endl;
         }
     }
@@ -1170,7 +1171,7 @@ bool Foam::polyMeshGeometry::checkVolRatio
     {
         label face0 = baffles[i].first();
         label face1 = baffles[i].second();
-        
+
         scalar ownVol = mag(cellVolumes[own[face0]]);
 
         scalar neiVol = mag(cellVolumes[own[face1]]);
@@ -1257,7 +1258,7 @@ bool Foam::polyMeshGeometry::checkFaceAngles
             << abort(FatalError);
     }
 
-    const scalar maxSin = Foam::sin(maxDeg/180.0*mathematicalConstant::pi);
+    const scalar maxSin = Foam::sin(maxDeg/180.0*constant::math::pi);
 
     const faceList& fcs = mesh.faces();
 
@@ -1338,7 +1339,7 @@ bool Foam::polyMeshGeometry::checkFaceAngles
         {
             scalar maxConcaveDegr =
                 Foam::asin(Foam::min(1.0, maxEdgeSin))
-             * 180.0/mathematicalConstant::pi;
+              *180.0/constant::math::pi;
 
             Info<< "There are " << nConcave
                 << " faces with concave angles between consecutive"
@@ -1823,7 +1824,7 @@ bool Foam::polyMeshGeometry::checkCellDeterminant
         forAll(cFaces, cFaceI)
         {
             label faceI = cFaces[cFaceI];
-    
+
             scalar magArea = mag(faceAreas[faceI]);
 
             magAreaSum += magArea;
@@ -1850,7 +1851,7 @@ bool Foam::polyMeshGeometry::checkCellDeterminant
             nWarnDet++;
         }
     }
-    
+
     reduce(minDet, minOp<scalar>());
     reduce(sumDet, sumOp<scalar>());
     reduce(nSumDet, sumOp<label>());

src/dynamicMesh/perfectInterface/perfectInterface.C

@@ -38,6 +38,7 @@ Description
 #include "polyModifyFace.H"
 #include "polyRemovePoint.H"
 #include "polyRemoveFace.H"
+#include "indirectPrimitivePatch.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -61,7 +62,7 @@ const Foam::scalar Foam::perfectInterface::tol_ = 1E-3;
 
 Foam::pointField Foam::perfectInterface::calcFaceCentres
 (
-    const primitivePatch& pp
+    const indirectPrimitivePatch& pp
 )
 {
     const pointField& points = pp.points();
@@ -155,6 +156,295 @@ bool Foam::perfectInterface::changeTopology() const
 }
 
 
+void Foam::perfectInterface::setRefinement
+(
+    const indirectPrimitivePatch& pp0,
+    const indirectPrimitivePatch& pp1,
+    polyTopoChange& ref
+) const
+{
+    const polyMesh& mesh = topoChanger().mesh();
+
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    // Some aliases
+    const edgeList& edges0 = pp0.edges();
+    const pointField& pts0 = pp0.localPoints();
+    const pointField& pts1 = pp1.localPoints();    
+    const labelList& meshPts0 = pp0.meshPoints();
+    const labelList& meshPts1 = pp1.meshPoints();
+
+
+    // Get local dimension as fraction of minimum edge length
+
+    scalar minLen = GREAT;
+
+    forAll(edges0, edgeI)
+    {
+        minLen = min(minLen, edges0[edgeI].mag(pts0));
+    }
+    scalar typDim = tol_*minLen;
+
+    if (debug)
+    {
+        Pout<< "typDim:" << typDim << " edges0:" << edges0.size()
+            << " pts0:" << pts0.size() << " pts1:" << pts1.size()
+            << " pp0:" << pp0.size() << " pp1:" << pp1.size() << endl;
+    }
+
+
+    // Determine pointMapping in mesh point labels. Uses geometric
+    // comparison to find correspondence between patch points.
+
+    labelList renumberPoints(mesh.points().size());
+    forAll(renumberPoints, i)
+    {
+        renumberPoints[i] = i;
+    }
+    {
+        labelList from1To0Points(pts1.size());
+
+        bool matchOk = matchPoints
+        (
+            pts1,
+            pts0,
+            scalarField(pts1.size(), typDim),   // tolerance
+            true,                               // verbose
+            from1To0Points
+        );
+
+        if (!matchOk)
+        {
+            FatalErrorIn
+            (
+                "perfectInterface::setRefinement(polyTopoChange& ref) const"
+            )   << "Points on patch sides do not match to within tolerance "
+                << typDim << exit(FatalError);
+        }
+
+        forAll(pts1, i)
+        {
+            renumberPoints[meshPts1[i]] = meshPts0[from1To0Points[i]];
+        }
+    }
+
+
+
+    // Calculate correspondence between patch faces
+
+    labelList from0To1Faces(pp1.size());
+
+    bool matchOk = matchPoints
+    (
+        calcFaceCentres(pp0),
+        calcFaceCentres(pp1),
+        scalarField(pp0.size(), typDim),    // tolerance
+        true,                               // verbose
+        from0To1Faces
+    );
+
+    if (!matchOk)
+    {
+        FatalErrorIn
+        (
+            "perfectInterface::setRefinement(polyTopoChange& ref) const"
+        )   << "Face centres of patch sides do not match to within tolerance "
+            << typDim << exit(FatalError);
+    }
+
+
+
+    // Now
+    // - renumber faces using pts1 (except patch1 faces)
+    // - remove patch1 faces. Remember cell label on owner side.
+    // - modify patch0 faces to be internal.
+
+    // 1. Get faces to be renumbered
+    labelHashSet affectedFaces(2*pp1.size());
+    forAll(meshPts1, i)
+    {
+        label meshPointI = meshPts1[i];
+
+        if (meshPointI != renumberPoints[meshPointI])
+        {
+            const labelList& pFaces = mesh.pointFaces()[meshPointI];
+
+            forAll(pFaces, pFaceI)
+            {
+                affectedFaces.insert(pFaces[pFaceI]);
+            }
+        }
+    }
+    forAll(pp1, i)
+    {
+        affectedFaces.erase(pp1.addressing()[i]);
+    }
+    // Remove patch0 from renumbered faces. Should not be nessecary since
+    // patch0 and 1 should not share any point (if created by mergeMeshing)
+    // so affectedFaces should not contain any patch0 faces but you can
+    // never be sure what the user is doing.
+    forAll(pp0, i)
+    {
+        label faceI = pp0.addressing()[i];
+
+        if (affectedFaces.erase(faceI))
+        {
+            WarningIn
+            (
+                "perfectInterface::setRefinement(polyTopoChange&) const"
+            )   << "Found face " << faceI << " vertices "
+                << mesh.faces()[faceI] << " whose points are"
+                << " used both by master patch and slave patch" << endl;
+        }
+    }
+
+
+    // 2. Renumber (non patch0/1) faces.
+    for
+    (
+        labelHashSet::const_iterator iter = affectedFaces.begin();
+        iter != affectedFaces.end();
+        ++iter
+    )
+    {
+        label faceI = iter.key();
+
+        const face& f = mesh.faces()[faceI];
+
+        face newFace(f.size());
+
+        forAll(newFace, fp)
+        {
+            newFace[fp] = renumberPoints[f[fp]];
+        }
+
+        label nbr = -1;
+
+        label patchI = -1;
+
+        if (mesh.isInternalFace(faceI))
+        {
+            nbr = mesh.faceNeighbour()[faceI];
+        }
+        else
+        {
+            patchI = patches.whichPatch(faceI);
+        }
+
+        label zoneID = mesh.faceZones().whichZone(faceI);
+
+        bool zoneFlip = false;
+
+        if (zoneID >= 0)
+        {
+            const faceZone& fZone = mesh.faceZones()[zoneID];
+
+            zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
+        }
+
+        ref.setAction
+        (
+            polyModifyFace
+            (
+                newFace,                    // modified face
+                faceI,                      // label of face being modified
+                mesh.faceOwner()[faceI],    // owner
+                nbr,                        // neighbour
+                false,                      // face flip
+                patchI,                     // patch for face
+                false,                      // remove from zone
+                zoneID,                     // zone for face
+                zoneFlip                    // face flip in zone
+            )
+        );
+    }
+
+
+    // 3. Remove patch1 points
+    forAll(meshPts1, i)
+    {
+        label meshPointI = meshPts1[i];
+
+        if (meshPointI != renumberPoints[meshPointI])
+        {
+            ref.setAction(polyRemovePoint(meshPointI));
+        }
+    }
+
+
+    // 4. Remove patch1 faces
+    forAll(pp1, i)
+    {
+        label faceI = pp1.addressing()[i];
+        ref.setAction(polyRemoveFace(faceI));
+    }
+
+
+    // 5. Modify patch0 faces for new points (not really nessecary; see
+    // comment above about patch1 and patch0 never sharing points) and
+    // becoming internal.
+    const boolList& mfFlip =
+        mesh.faceZones()[faceZoneID_.index()].flipMap();
+
+    forAll(pp0, i)
+    {
+        label faceI = pp0.addressing()[i];
+
+        const face& f = mesh.faces()[faceI];
+
+        face newFace(f.size());
+
+        forAll(newFace, fp)
+        {
+            newFace[fp] = renumberPoints[f[fp]];
+        }
+
+        label own = mesh.faceOwner()[faceI];
+
+        label pp1FaceI = pp1.addressing()[from0To1Faces[i]];
+
+        label nbr = mesh.faceOwner()[pp1FaceI];
+
+        if (own < nbr)
+        {
+            ref.setAction
+            (
+                polyModifyFace
+                (
+                    newFace,                // modified face
+                    faceI,                  // label of face being modified
+                    own,                    // owner
+                    nbr,                    // neighbour
+                    false,                  // face flip
+                    -1,                     // patch for face
+                    false,                  // remove from zone
+                    faceZoneID_.index(),    // zone for face
+                    mfFlip[i]               // face flip in zone
+                )
+            );
+        }
+        else
+        {
+            ref.setAction
+            (
+                polyModifyFace
+                (
+                    newFace.reverseFace(),  // modified face
+                    faceI,                  // label of face being modified
+                    nbr,                    // owner
+                    own,                    // neighbour
+                    true,                   // face flip
+                    -1,                     // patch for face
+                    false,                  // remove from zone
+                    faceZoneID_.index(),    // zone for face
+                    !mfFlip[i]              // face flip in zone
+                )
+            );
+        }
+    }
+}
+
+
 void Foam::perfectInterface::setRefinement(polyTopoChange& ref) const
 {
     if (debug)
@@ -176,286 +466,25 @@ void Foam::perfectInterface::setRefinement(polyTopoChange& ref) const
         const polyMesh& mesh = topoChanger().mesh();
 
         const polyBoundaryMesh& patches = mesh.boundaryMesh();
-
-        const polyPatch& pp0 = patches[masterPatchID_.index()];
-        const polyPatch& pp1 = patches[slavePatchID_.index()];
-
-        // Some aliases
-        const edgeList& edges0 = pp0.edges();
-        const pointField& pts0 = pp0.localPoints();
-        const pointField& pts1 = pp1.localPoints();    
-        const labelList& meshPts0 = pp0.meshPoints();
-        const labelList& meshPts1 = pp1.meshPoints();
-                        
-
-        // Get local dimension as fraction of minimum edge length
-
-        scalar minLen = GREAT;
-
-        forAll(edges0, edgeI)
-        {
-            minLen = min(minLen, edges0[edgeI].mag(pts0));
-        }
-        scalar typDim = tol_*minLen;
-
-        if (debug)
-        {
-            Pout<< "typDim:" << typDim << " edges0:" << edges0.size()
-                << " pts0:" << pts0.size() << " pts1:" << pts1.size()
-                << " pp0:" << pp0.size() << " pp1:" << pp1.size() << endl;
-        }
+        const polyPatch& patch0 = patches[masterPatchID_.index()];
+        const polyPatch& patch1 = patches[slavePatchID_.index()];
 
 
-        // Determine pointMapping in mesh point labels. Uses geometric
-        // comparison to find correspondence between patch points.
-
-        labelList renumberPoints(mesh.points().size());
-        forAll(renumberPoints, i)
-        {
-            renumberPoints[i] = i;
-        }
-        {
-            labelList from1To0Points(pts1.size());
-
-            bool matchOk = matchPoints
-            (
-                pts1,
-                pts0,
-                scalarField(pts1.size(), typDim),   // tolerance
-                true,                               // verbose
-                from1To0Points
-            );
-
-            if (!matchOk)
-            {
-                FatalErrorIn
-                (
-                    "perfectInterface::setRefinement(polyTopoChange& ref) const"
-                )   << "Points on patches " << pp0.name() << " and "
-                    << pp1.name() << " do not match to within tolerance "
-                    << typDim << exit(FatalError);
-            }
-
-            forAll(pts1, i)
-            {
-                renumberPoints[meshPts1[i]] = meshPts0[from1To0Points[i]];
-            }
-        }
-
-
-
-        // Calculate correspondence between patch faces
-
-        labelList from0To1Faces(pp1.size());
-
-        bool matchOk = matchPoints
+        labelList pp0Labels(identity(patch0.size())+patch0.start());
+        indirectPrimitivePatch pp0
         (
-            calcFaceCentres(pp0),
-            calcFaceCentres(pp1),
-            scalarField(pp0.size(), typDim),    // tolerance
-            true,                               // verbose
-            from0To1Faces
+            IndirectList<face>(mesh.faces(), pp0Labels),
+            mesh.points()
         );
 
-        if (!matchOk)
-        {
-            FatalErrorIn
-            (
-                "perfectInterface::setRefinement(polyTopoChange& ref) const"
-            )   << "Face centres of patches " << pp0.name() << " and "
-                << pp1.name() << " do not match to within tolerance " << typDim
-                << exit(FatalError);
-        }
-
-
-
-        // Now
-        // - renumber faces using pts1 (except patch1 faces)
-        // - remove patch1 faces. Remember cell label on owner side.
-        // - modify patch0 faces to be internal.
-
-        // 1. Get faces to be renumbered
-        labelHashSet affectedFaces(2*pp1.size());
-        forAll(meshPts1, i)
-        {
-            label meshPointI = meshPts1[i];
-
-            if (meshPointI != renumberPoints[meshPointI])
-            {
-                const labelList& pFaces = mesh.pointFaces()[meshPointI];
-
-                forAll(pFaces, pFaceI)
-                {
-                    affectedFaces.insert(pFaces[pFaceI]);
-                }
-            }
-        }
-        forAll(pp1, i)
-        {
-            affectedFaces.erase(pp1.start() + i);
-        }
-        // Remove patch0 from renumbered faces. Should not be nessecary since
-        // patch0 and 1 should not share any point (if created by mergeMeshing)
-        // so affectedFaces should not contain any patch0 faces but you can
-        // never be sure what the user is doing.
-        forAll(pp0, i)
-        {
-            if (affectedFaces.erase(pp0.start() + i))
-            {
-                WarningIn
-                (
-                    "perfectInterface::setRefinement(polyTopoChange&) const"
-                )   << "Found face " << pp0.start() + i << " vertices "
-                    << mesh.faces()[pp0.start() + i] << " whose points are"
-                    << " used both by master patch " << pp0.name()
-                    << " and slave patch " << pp1.name()
-                    << endl;
-            }
-        }
-
-
-        // 2. Renumber (non patch0/1) faces.
-        for
+        labelList pp1Labels(identity(patch1.size())+patch1.start());
+        indirectPrimitivePatch pp1
         (
-            labelHashSet::const_iterator iter = affectedFaces.begin();
-            iter != affectedFaces.end();
-            ++iter
-        )
-        {
-            label faceI = iter.key();
+            IndirectList<face>(mesh.faces(), pp1Labels),
+            mesh.points()
+        );
 
-            const face& f = mesh.faces()[faceI];
-
-            face newFace(f.size());
-
-            forAll(newFace, fp)
-            {
-                newFace[fp] = renumberPoints[f[fp]];
-            }
-
-            label nbr = -1;
-
-            label patchI = -1;
-
-            if (mesh.isInternalFace(faceI))
-            {
-                nbr = mesh.faceNeighbour()[faceI];
-            }
-            else
-            {
-                patchI = patches.whichPatch(faceI);
-            }
-
-            label zoneID = mesh.faceZones().whichZone(faceI);
-
-            bool zoneFlip = false;
-
-            if (zoneID >= 0)
-            {
-                const faceZone& fZone = mesh.faceZones()[zoneID];
-
-                zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
-            }
-
-            ref.setAction
-            (
-                polyModifyFace
-                (
-                    newFace,                    // modified face
-                    faceI,                      // label of face being modified
-                    mesh.faceOwner()[faceI],    // owner
-                    nbr,                        // neighbour
-                    false,                      // face flip
-                    patchI,                     // patch for face
-                    false,                      // remove from zone
-                    zoneID,                     // zone for face
-                    zoneFlip                    // face flip in zone
-                )
-            );
-        }
-
-
-        // 3. Remove patch1 points
-        forAll(meshPts1, i)
-        {
-            label meshPointI = meshPts1[i];
-
-            if (meshPointI != renumberPoints[meshPointI])
-            {
-                ref.setAction(polyRemovePoint(meshPointI));
-            }
-        }
-
-
-        // 4. Remove patch1 faces
-        forAll(pp1, i)
-        {
-            ref.setAction(polyRemoveFace(pp1.start() + i));
-        }
-
-
-        // 5. Modify patch0 faces for new points (not really nessecary; see
-        // comment above about patch1 and patch0 never sharing points) and
-        // becoming internal.
-        const boolList& mfFlip =
-            mesh.faceZones()[faceZoneID_.index()].flipMap();
-
-        forAll(pp0, i)
-        {
-            label faceI = pp0.start() + i;
-
-            const face& f = mesh.faces()[faceI];
-
-            face newFace(f.size());
-
-            forAll(newFace, fp)
-            {
-                newFace[fp] = renumberPoints[f[fp]];
-            }
-
-            label own = mesh.faceOwner()[faceI];
-
-            label pp1FaceI = pp1.start() + from0To1Faces[i];
-
-            label nbr = mesh.faceOwner()[pp1FaceI];
-
-            if (own < nbr)
-            {
-                ref.setAction
-                (
-                    polyModifyFace
-                    (
-                        newFace,                // modified face
-                        faceI,                  // label of face being modified
-                        own,                    // owner
-                        nbr,                    // neighbour
-                        false,                  // face flip
-                        -1,                     // patch for face
-                        false,                  // remove from zone
-                        faceZoneID_.index(),    // zone for face
-                        mfFlip[i]               // face flip in zone
-                    )
-                );
-            }
-            else
-            {
-                ref.setAction
-                (
-                    polyModifyFace
-                    (
-                        newFace.reverseFace(),  // modified face
-                        faceI,                  // label of face being modified
-                        nbr,                    // owner
-                        own,                    // neighbour
-                        false,                  // face flip
-                        -1,                     // patch for face
-                        false,                  // remove from zone
-                        faceZoneID_.index(),    // zone for face
-                        !mfFlip[i]              // face flip in zone
-                    )
-                );
-            }
-        }
+        setRefinement(pp0, pp1, ref);
     }
 }
 

src/dynamicMesh/perfectInterface/perfectInterface.H

@@ -40,6 +40,7 @@ SourceFiles
 #include "polyMeshModifier.H"
 #include "polyPatchID.H"
 #include "ZoneIDs.H"
+#include "indirectPrimitivePatch.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -75,7 +76,7 @@ class perfectInterface
     // Private Member Functions
 
         //- Calculate face centres on patch
-        static pointField calcFaceCentres(const primitivePatch&);
+        static pointField calcFaceCentres(const indirectPrimitivePatch&);
 
 
         //- Disallow default bitwise copy construct
@@ -128,6 +129,17 @@ public:
         //  into the topological change
         virtual void setRefinement(polyTopoChange&) const;
 
+        //- Insert the layer addition/removal instructions
+        //  into the topological change. Uses only mesh, not any of the
+        //  patch and zone indices. Bit of a workaround - used in extruding
+        //  a mesh.
+        virtual void setRefinement
+        (
+            const indirectPrimitivePatch& pp0,
+            const indirectPrimitivePatch& pp1,
+            polyTopoChange&
+        ) const;
+
         //- Modify motion points to comply with the topological change
         virtual void modifyMotionPoints(pointField& motionPoints) const;
 

src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C

@@ -335,17 +335,19 @@ Foam::label Foam::addPatchCellLayer::addSideFace
     label inflateEdgeI = -1;
 
     // Check mesh faces using edge
-    forAll(meshFaces, i)
+    if (addToMesh_)
     {
-        if (mesh_.isInternalFace(meshFaces[i]))
+        forAll(meshFaces, i)
         {
-            // meshEdge uses internal faces so ok to inflate from it
-            inflateEdgeI = meshEdgeI;
-            break;
+            if (mesh_.isInternalFace(meshFaces[i]))
+            {
+                // meshEdge uses internal faces so ok to inflate from it
+                inflateEdgeI = meshEdgeI;
+                break;
+            }
         }
     }
 
-
     // Get my mesh face and its zone.
     label meshFaceI = pp.addressing()[ownFaceI];
     label zoneI = mesh_.faceZones().whichZone(meshFaceI);
@@ -504,9 +506,14 @@ Foam::label Foam::addPatchCellLayer::addSideFace
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct from mesh
-Foam::addPatchCellLayer::addPatchCellLayer(const polyMesh& mesh)
+Foam::addPatchCellLayer::addPatchCellLayer
+(
+    const polyMesh& mesh,
+    const bool addToMesh
+)
 :
     mesh_(mesh),
+    addToMesh_(addToMesh),
     addedPoints_(0),
     layerFaces_(0)
 {}
@@ -551,6 +558,7 @@ void Foam::addPatchCellLayer::setRefinement
 (
     const scalarField& expansionRatio,
     const indirectPrimitivePatch& pp,
+    const labelList& exposedPatchID,
     const labelList& nFaceLayers,
     const labelList& nPointLayers,
     const vectorField& firstLayerDisp,
@@ -620,7 +628,7 @@ void Foam::addPatchCellLayer::setRefinement
     DynamicList<label> ef;
 
     // Precalculate mesh edges for pp.edges.
-    labelList meshEdges(calcMeshEdges(mesh_, pp));
+    const labelList meshEdges(pp.meshEdges(mesh_.edges(), mesh_.pointEdges()));
 
     if (debug)
     {
@@ -827,6 +835,19 @@ void Foam::addPatchCellLayer::setRefinement
     }
 
 
+    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+
+    // Precalculated patchID for each patch face
+    labelList patchID(pp.size());
+
+    forAll(pp, patchFaceI)
+    {
+        label meshFaceI = pp.addressing()[patchFaceI];
+
+        patchID[patchFaceI] = patches.whichPatch(meshFaceI);
+    }
+
+
     // From master point (in patch point label) to added points (in mesh point
     // label)
     addedPoints_.setSize(pp.nPoints());
@@ -857,6 +878,33 @@ void Foam::addPatchCellLayer::setRefinement
     // Create new points
     //
 
+    // If creating new mesh: copy existing patch points
+    labelList copiedPatchPoints;
+    if (!addToMesh_)
+    {
+        copiedPatchPoints.setSize(firstLayerDisp.size());
+        forAll(firstLayerDisp, patchPointI)
+        {
+            if (addedPoints_[patchPointI].size())
+            {
+                label meshPointI = meshPoints[patchPointI];
+                label zoneI = mesh_.pointZones().whichZone(meshPointI);
+                copiedPatchPoints[patchPointI] = meshMod.setAction
+                (
+                    polyAddPoint
+                    (
+                        mesh_.points()[meshPointI],         // point
+                        -1,         // master point
+                        zoneI,      // zone for point
+                        true        // supports a cell
+                    )
+                );
+            }
+        }
+    }
+
+
+    // Create points for additional layers
     forAll(firstLayerDisp, patchPointI)
     {
         if (addedPoints_[patchPointI].size())
@@ -878,7 +926,7 @@ void Foam::addPatchCellLayer::setRefinement
                     polyAddPoint
                     (
                         pt,         // point
-                        meshPointI, // master point
+                        (addToMesh_ ? meshPointI : -1), // master point
                         zoneI,      // zone for point
                         true        // supports a cell
                     )
@@ -922,34 +970,15 @@ void Foam::addPatchCellLayer::setRefinement
                         -1,             // master point
                         -1,             // master edge
                         -1,             // master face
-                        mesh_.faceOwner()[meshFaceI],   // master cell id
+                        (addToMesh_ ? mesh_.faceOwner()[meshFaceI] : -1),//master
                         ownZoneI        // zone for cell
                     )
                 );
-
-                //Pout<< "For patchFace:" << patchFaceI
-                //    << " meshFace:" << pp.addressing()[patchFaceI]
-                //    << " layer:" << i << " added cell:"
-                //    << addedCells[patchFaceI][i]
-                //    << endl;
             }
         }
     }
 
 
-    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-
-    // Precalculated patchID for each patch face
-    labelList patchID(pp.size());
-
-    forAll(pp, patchFaceI)
-    {
-        label meshFaceI = pp.addressing()[patchFaceI];
-
-        patchID[patchFaceI] = patches.whichPatch(meshFaceI);
-    }
-
-
 
     // Create faces on top of the original patch faces.
     // These faces are created from original patch faces outwards so in order
@@ -965,7 +994,6 @@ void Foam::addPatchCellLayer::setRefinement
         if (addedCells[patchFaceI].size())
         {
             layerFaces_[patchFaceI].setSize(addedCells[patchFaceI].size() + 1);
-            layerFaces_[patchFaceI][0] = meshFaceI;
 
             label zoneI = mesh_.faceZones().whichZone(meshFaceI);
 
@@ -981,7 +1009,12 @@ void Foam::addPatchCellLayer::setRefinement
                     if (addedPoints_[f[fp]].empty())
                     {
                         // Keep original point
-                        newFace[fp] = meshPoints[f[fp]];
+                        newFace[fp] =
+                        (
+                            addToMesh_
+                          ? meshPoints[f[fp]]
+                          : copiedPatchPoints[f[fp]]
+                        );
                     }
                     else
                     {
@@ -1021,18 +1054,13 @@ void Foam::addPatchCellLayer::setRefinement
                         nei,                        // neighbour
                         -1,                         // master point
                         -1,                         // master edge
-                        meshFaceI,                  // master face for addition
+                        (addToMesh_ ? meshFaceI : -1), // master face
                         false,                      // flux flip
                         patchI,                     // patch for face
                         zoneI,                      // zone for face
                         false                       // face zone flip
                     )
                 );
-                //Pout<< "Added inbetween face " << newFace
-                //    << " own:" << addedCells[patchFaceI][i]
-                //    << " nei:" << nei
-                //    << " patch:" << patchI
-                //    << endl;
             }
         }
     }
@@ -1040,35 +1068,81 @@ void Foam::addPatchCellLayer::setRefinement
     //
     // Modify old patch faces to be on the inside
     //
-    forAll(pp, patchFaceI)
+
+    if (addToMesh_)
     {
-        if (addedCells[patchFaceI].size())
+        forAll(pp, patchFaceI)
         {
-            label meshFaceI = pp.addressing()[patchFaceI];
+            if (addedCells[patchFaceI].size())
+            {
+                label meshFaceI = pp.addressing()[patchFaceI];
 
-            label zoneI = mesh_.faceZones().whichZone(meshFaceI);
+                layerFaces_[patchFaceI][0] = meshFaceI;
 
-            meshMod.setAction
-            (
-                polyModifyFace
+                label zoneI = mesh_.faceZones().whichZone(meshFaceI);
+
+                meshMod.setAction
                 (
-                    pp[patchFaceI],                 // modified face
-                    meshFaceI,                      // label of face
-                    mesh_.faceOwner()[meshFaceI],   // owner
-                    addedCells[patchFaceI][0],      // neighbour
-                    false,                          // face flip
-                    -1,                             // patch for face
-                    false,                          // remove from zone
-                    zoneI,                          // zone for face
-                    false                           // face flip in zone
-                )
-            );
-            //Pout<< "Modified old patch face " << meshFaceI
-            //    << " own:" << mesh_.faceOwner()[meshFaceI]
-            //    << " nei:" << addedCells[patchFaceI][0]
-            //    << endl;
+                    polyModifyFace
+                    (
+                        pp[patchFaceI],                 // modified face
+                        meshFaceI,                      // label of face
+                        mesh_.faceOwner()[meshFaceI],   // owner
+                        addedCells[patchFaceI][0],      // neighbour
+                        false,                          // face flip
+                        -1,                             // patch for face
+                        false,                          // remove from zone
+                        zoneI,                          // zone for face
+                        false                           // face flip in zone
+                    )
+                );
+            }
         }
     }
+    else
+    {
+        // If creating new mesh: reverse original faces and put them
+        // in the exposed patch ID.
+        forAll(pp, patchFaceI)
+        {
+            if (nFaceLayers[patchFaceI] > 0)
+            {
+                label meshFaceI = pp.addressing()[patchFaceI];
+                label zoneI = mesh_.faceZones().whichZone(meshFaceI);
+                bool zoneFlip = false;
+                if (zoneI != -1)
+                {
+                    const faceZone& fz = mesh_.faceZones()[zoneI];
+                    zoneFlip = !fz.flipMap()[fz.whichFace(meshFaceI)];
+                }
+
+                // Reverse and renumber old patch face.
+                face f(pp.localFaces()[patchFaceI].reverseFace());
+                forAll(f, fp)
+                {
+                    f[fp] = copiedPatchPoints[f[fp]];
+                }
+
+                layerFaces_[patchFaceI][0] = meshMod.setAction
+                (
+                    polyAddFace
+                    (
+                        f,                          // modified face
+                        addedCells[patchFaceI][0],  // owner
+                        -1,                         // neighbour
+                        -1,                         // masterPoint
+                        -1,                         // masterEdge
+                        -1,                         // masterFace
+                        true,                       // face flip
+                        exposedPatchID[patchFaceI], // patch for face
+                        zoneI,                      // zone for face
+                        zoneFlip                    // face flip in zone
+                    )
+                );
+            }
+        }
+    }
+
 
 
     //
@@ -1255,7 +1329,16 @@ void Foam::addPatchCellLayer::setRefinement
                         forAll(stringedVerts, stringedI)
                         {
                             label v = stringedVerts[stringedI];
-                            addVertex(meshPoints[v], newFace, newFp);
+                            addVertex
+                            (
+                                (
+                                    addToMesh_
+                                  ? meshPoints[v]
+                                  : copiedPatchPoints[v]
+                                ),
+                                newFace,
+                                newFp
+                            );
                         }
                     }
                     else
@@ -1276,7 +1359,16 @@ void Foam::addPatchCellLayer::setRefinement
                             }
                             else
                             {
-                                addVertex(meshPoints[v], newFace, newFp);
+                                addVertex
+                                (
+                                    (
+                                        addToMesh_
+                                      ? meshPoints[v]
+                                      : copiedPatchPoints[v]
+                                    ),
+                                    newFace,
+                                    newFp
+                                );
                             }
                         }
                     }
@@ -1316,7 +1408,16 @@ void Foam::addPatchCellLayer::setRefinement
                         }
                         else
                         {
-                            addVertex(meshPoints[v], newFace, newFp);
+                            addVertex
+                            (
+                                (
+                                    addToMesh_
+                                  ? meshPoints[v]
+                                  : copiedPatchPoints[v]
+                                ),
+                                newFace,
+                                newFp
+                            );
                         }
                     }
 
@@ -1446,30 +1547,4 @@ void Foam::addPatchCellLayer::updateMesh
 }
 
 
-Foam::labelList Foam::addPatchCellLayer::calcMeshEdges
-(
-    const primitiveMesh& mesh,
-    const indirectPrimitivePatch& pp
-)
-{
-    labelList meshEdges(pp.nEdges());
-
-    forAll(meshEdges, patchEdgeI)
-    {
-        const edge& e = pp.edges()[patchEdgeI];
-
-        label v0 = pp.meshPoints()[e[0]];
-        label v1 = pp.meshPoints()[e[1]];
-        meshEdges[patchEdgeI] = meshTools::findEdge
-        (
-            mesh.edges(),
-            mesh.pointEdges()[v0],
-            v0,
-            v1
-        );
-    }
-    return meshEdges;
-}
-
-
 // ************************************************************************* //

src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.H

@@ -26,7 +26,8 @@ Class
     Foam::addPatchCellLayer
 
 Description
-    Adds layers of cells to outside of polyPatch.
+    Adds layers of cells to outside of polyPatch. Can optionally create
+    stand-alone extruded mesh (addToMesh=false).
 
     Call setRefinement with offset vector for every patch point and number
     of layers per patch face and number of layers per patch point.
@@ -164,6 +165,9 @@ class addPatchCellLayer
         //- Reference to mesh
         const polyMesh& mesh_;
 
+        //- Add layers to existing mesh or create new mesh
+        const bool addToMesh_;
+
         //- For all patchpoints: list of added points (size 0 or nLayers)
         //  First point in list is one nearest to original point in patch,
         //  last one is the new point on the surface.
@@ -253,7 +257,7 @@ public:
     // Constructors
 
         //- Construct from mesh.
-        addPatchCellLayer(const polyMesh& mesh);
+        addPatchCellLayer(const polyMesh& mesh, const bool addToMesh = true);
 
 
     // Member Functions
@@ -291,8 +295,10 @@ public:
             //- Play commands into polyTopoChange to create layers on top
             //  of indirectPrimitivePatch (have to be outside faces).
             //  Gets displacement per patch point.
-            //  - nPointLayers : number of layers per (patch)point
-            //  - nFaceLayers : number of layers per (patch) face
+            //  - exposedPatchID : only used if creating a new mesh (addToMesh=false)
+            //    gives per pp face the patch the exposed face should get.
+            //  - nPointLayers : number of layers per (patch)point.
+            //  - nFaceLayers : number of layers per (patch) face.
             //  - firstDisplacement : displacement per point for first
             //    layer of points (i.e. nearest to original mesh). If zero
             //    do not add point.
@@ -305,14 +311,15 @@ public:
             //        get a cell should firstDisplacement be <> 0
             //  Note: cells get added from owner cells of patch faces
             //        (instead of e.g. from patch faces)
-           void setRefinement
-           (
-               const scalarField& expansionRatio,
-               const indirectPrimitivePatch& pp,
-               const labelList& nFaceLayers,
-               const labelList& nPointLayers,
-               const vectorField& firstLayerDisp,
-               polyTopoChange& meshMod
+            void setRefinement
+            (
+                const scalarField& expansionRatio,
+                const indirectPrimitivePatch& pp,
+                const labelList& exposedPatchID,
+                const labelList& nFaceLayers,
+                const labelList& nPointLayers,
+                const vectorField& firstLayerDisp,
+                polyTopoChange& meshMod
             );
 
 
@@ -329,6 +336,7 @@ public:
                 (
                     scalarField(pp.nPoints(), 1.0), // expansion ration
                     pp,
+                    labelList(0),
                     labelList(pp.size(), nLayers),
                     labelList(pp.nPoints(), nLayers),
                     overallDisplacement / nLayers,
@@ -346,15 +354,6 @@ public:
                 const labelList& faceMap,   // new to old patch faces
                 const labelList& pointMap   // new to old patch points
             );
-
-        // Helper
-
-            //- Per patch edge the corresponding mesh edge
-            static labelList calcMeshEdges
-            (
-                const primitiveMesh& mesh,
-                const indirectPrimitivePatch&
-            );
 };
 
 

src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C

@@ -749,6 +749,11 @@ void Foam::polyTopoChange::getFaceOrder
                 " const"
             )   << "Did not determine new position"
                 << " for face " << faceI
+                << " owner " << faceOwner_[faceI]
+                << " neighbour " << faceNeighbour_[faceI]
+                << " region " << region_[faceI] << endl
+                << "This is usually caused by not specifying a patch for"
+                << " a boundary face."
                 << abort(FatalError);
         }
     }
@@ -3176,7 +3181,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChange::makeMesh
 (
     autoPtr<fvMesh>& newMeshPtr,
     const IOobject& io,
-    const fvMesh& mesh,
+    const polyMesh& mesh,
     const bool syncParallel,
     const bool orderCells,
     const bool orderPoints

src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.H

@@ -585,7 +585,7 @@ public:
             (
                 autoPtr<fvMesh>& newMesh,
                 const IOobject& io,
-                const fvMesh& mesh,
+                const polyMesh& mesh,
                 const bool syncParallel = true,
                 const bool orderCells = false,
                 const bool orderPoints = false

src/engine/engineTime/engineTime.C

@@ -25,7 +25,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "engineTime.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -125,7 +125,7 @@ bool Foam::engineTime::read()
 
 Foam::scalar Foam::engineTime::degToRad(const scalar deg) const
 {
-    return mathematicalConstant::pi*deg/180.0;
+    return constant::math::pi*deg/180.0;
 }
 
 
@@ -239,6 +239,4 @@ Foam::scalar Foam::engineTime::timeToUserTime(const scalar t) const
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 // ************************************************************************* //

src/engine/include/StCorr.H

@@ -27,16 +27,16 @@
                         )
                     );
 
-                    Ak = sphereFraction*4.0*mathematicalConstant::pi
+                    Ak = sphereFraction*4.0*constant::math::pi
                        *pow
                         (
                             3.0*Vk
-                           /(sphereFraction*4.0*mathematicalConstant::pi),
+                           /(sphereFraction*4.0*constant::math::pi),
                             2.0/3.0
                         );
                 }
                 break;
-            
+
                 case 2:
                 {
                     // Assume it is part-circular
@@ -56,11 +56,11 @@
                         )
                     );
 
-                    Ak = circleFraction*mathematicalConstant::pi*thickness
+                    Ak = circleFraction*constant::math::pi*thickness
                        *sqrt
                         (
                             4.0*Vk
-                           /(circleFraction*thickness*mathematicalConstant::pi)
+                           /(circleFraction*thickness*constant::math::pi)
                         );
                 }
                 break;

src/finiteVolume/cfdTools/general/SRF/SRFModel/rpm/rpm.C

@@ -22,13 +22,11 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Description
-
 \*---------------------------------------------------------------------------*/
 
 #include "rpm.H"
 #include "addToRunTimeSelectionTable.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -58,7 +56,7 @@ Foam::SRF::rpm::rpm
     rpm_(readScalar(SRFModelCoeffs_.lookup("rpm")))
 {
     // Initialise the angular velocity
-    omega_.value() = axis_*rpm_*2.0*mathematicalConstant::pi/60.0;
+    omega_.value() = axis_*rpm_*constant::math::twoPi/60.0;
 }
 
 
@@ -78,7 +76,7 @@ bool Foam::SRF::rpm::read()
         SRFModelCoeffs_.lookup("rpm") >> rpm_;
 
         // Update angular velocity
-        omega_.value() = axis_*rpm_*(2.0*mathematicalConstant::pi/60.0);
+        omega_.value() = axis_*rpm_*(constant::math::twoPi/60.0);
 
         return true;
     }

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

@@ -14,7 +14,7 @@
 #include "fixedValueFvPatchFields.H"
 #include "adjustPhi.H"
 #include "findRefCell.H"
-#include "mathematicalConstants.H"
+#include "constants.H"
 
 #include "OSspecific.H"
 #include "argList.H"

src/finiteVolume/fields/fvPatchFields/derived/oscillatingFixedValue/oscillatingFixedValueFvPatchField.C

@@ -25,7 +25,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "oscillatingFixedValueFvPatchField.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -39,8 +39,8 @@ scalar oscillatingFixedValueFvPatchField<Type>::currentScale() const
 {
     return
         1.0
-      + amplitude_*
-        sin(2*mathematicalConstant::pi*frequency_*this->db().time().value());
+      + amplitude_
+       *sin(constant::math::twoPi*frequency_*this->db().time().value());
 }
 
 

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

@@ -45,7 +45,7 @@ Foam::linearUpwind<Type>::correction
         (
             IOobject
             (
-                "linearUpwindCorrection(" + vf.name() + ')',
+                "linearUpwind::correction(" + vf.name() + ')',
                 mesh.time().timeName(),
                 mesh,
                 IOobject::NO_READ,

src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/linearUpwind/linearUpwindV.C

@@ -38,93 +38,96 @@ Foam::linearUpwindV<Type>::correction
     const GeometricField<Type, fvPatchField, volMesh>& vf
 ) const
 {
-   const fvMesh& mesh = this->mesh();
+    const fvMesh& mesh = this->mesh();
 
-   tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > tsfCorr
-   (
-       new GeometricField<Type, fvsPatchField, surfaceMesh>
-       (
-           IOobject
-           (
-               vf.name(),
-               mesh.time().timeName(),
-               mesh
-           ),
-           mesh,
-           dimensioned<Type>
-           (
-               vf.name(),
-               vf.dimensions(),
-               pTraits<Type>::zero
-           )
-       )
-   );
+    tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > tsfCorr
+    (
+        new GeometricField<Type, fvsPatchField, surfaceMesh>
+        (
+            IOobject
+            (
+                "linearUpwindV::correction(" + vf.name() + ')',
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh,
+            dimensioned<Type>
+            (
+                vf.name(),
+                vf.dimensions(),
+                pTraits<Type>::zero
+            )
+        )
+    );
 
-   GeometricField<Type, fvsPatchField, surfaceMesh>& sfCorr = tsfCorr();
+    GeometricField<Type, fvsPatchField, surfaceMesh>& sfCorr = tsfCorr();
 
-   const surfaceScalarField& faceFlux = this->faceFlux_;
-   const surfaceScalarField& w = mesh.weights();
+    const surfaceScalarField& faceFlux = this->faceFlux_;
+    const surfaceScalarField& w = mesh.weights();
 
-   const labelList& own = mesh.owner();
-   const labelList& nei = mesh.neighbour();
+    const labelList& own = mesh.owner();
+    const labelList& nei = mesh.neighbour();
 
-   const vectorField& C = mesh.C();
-   const vectorField& Cf = mesh.Cf();
+    const vectorField& C = mesh.C();
+    const vectorField& Cf = mesh.Cf();
 
-   GeometricField
-       <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>
-       gradVf = gradScheme_().grad(vf);
+    GeometricField
+        <typename outerProduct<vector, Type>::type, fvPatchField, volMesh>
+        gradVf = gradScheme_().grad(vf);
 
-   forAll(faceFlux, facei)
-   {
-       vector maxCorr;
+    forAll(faceFlux, facei)
+    {
+        vector maxCorr;
 
-       if (faceFlux[facei] > 0.0)
-       {
-           maxCorr =
-               (1.0 - w[facei])
+        if (faceFlux[facei] > 0.0)
+        {
+            maxCorr =
+                (1.0 - w[facei])
                *(vf[nei[facei]] - vf[own[facei]]);
 
-           sfCorr[facei] =
-               (Cf[facei] - C[own[facei]]) & gradVf[own[facei]];
-       }
-       else
-       {
-           maxCorr =
-               w[facei]*(vf[own[facei]] - vf[nei[facei]]);
+            sfCorr[facei] =
+                (Cf[facei] - C[own[facei]]) & gradVf[own[facei]];
+        }
+        else
+        {
+            maxCorr =
+                w[facei]*(vf[own[facei]] - vf[nei[facei]]);
 
-           sfCorr[facei] =
+            sfCorr[facei] =
                (Cf[facei] - C[nei[facei]]) & gradVf[nei[facei]];
-       }
+        }
 
-       scalar sfCorrs = magSqr(sfCorr[facei]);
-       scalar maxCorrs = sfCorr[facei] & maxCorr;
+        scalar sfCorrs = magSqr(sfCorr[facei]);
+        scalar maxCorrs = sfCorr[facei] & maxCorr;
 
-       if (sfCorrs > 0)
-       {
-           if (maxCorrs < 0)
-           {
-               sfCorr[facei] = vector::zero;
-           }
-           else if (sfCorrs > maxCorrs)
-           {
-               sfCorr[facei] *= maxCorrs/(sfCorrs + VSMALL);
-           }
-       }
-       else if (sfCorrs < 0)
-       {
-           if (maxCorrs > 0)
-           {
-               sfCorr[facei] = vector::zero;
-           }
-           else if (sfCorrs < maxCorrs)
-           {
-               sfCorr[facei] *= maxCorrs/(sfCorrs - VSMALL);
-           }
-       }
-   }
+        if (sfCorrs > 0)
+        {
+            if (maxCorrs < 0)
+            {
+                sfCorr[facei] = vector::zero;
+            }
+            else if (sfCorrs > maxCorrs)
+            {
+                sfCorr[facei] *= maxCorrs/(sfCorrs + VSMALL);
+            }
+        }
+        else if (sfCorrs < 0)
+        {
+            if (maxCorrs > 0)
+            {
+                sfCorr[facei] = vector::zero;
+            }
+            else if (sfCorrs < maxCorrs)
+            {
+                sfCorr[facei] *= maxCorrs/(sfCorrs - VSMALL);
+            }
+        }
+    }
 
-   return tsfCorr;
+    return tsfCorr;
 }
 
 

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

@@ -131,9 +131,12 @@ public:
                 (
                     IOobject
                     (
-                        vf.name(),
+                        "cubic::correction(" + vf.name() +')',
                         mesh.time().timeName(),
-                        mesh
+                        mesh,
+                        IOobject::NO_READ,
+                        IOobject::NO_WRITE,
+                        false
                     ),
                     surfaceInterpolationScheme<Type>::interpolate(vf, kSc, -kSc)
                 )

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

@@ -133,7 +133,7 @@ public:
                 (
                     IOobject
                     (
-                        vf.name(),
+                        "localMax::interpolate(" + vf.name() + ')',
                         mesh.time().timeName(),
                         mesh
                     ),

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

@@ -133,7 +133,7 @@ public:
                 (
                     IOobject
                     (
-                        vf.name(),
+                        "localMin::interpolate(" + vf.name() + ')',
                         mesh.time().timeName(),
                         mesh
                     ),

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

@@ -144,7 +144,7 @@ public:
                 (
                     IOobject
                     (
-                        vf.name(),
+                        "skewCorrected::skewCorrection(" + vf.name() + ')',
                         mesh.time().timeName(),
                         mesh
                     ),

src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.C

@@ -32,7 +32,7 @@ Description
 #include "surfaceFields.H"
 #include "demandDrivenData.H"
 #include "coupledFvPatch.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -353,7 +353,7 @@ void surfaceInterpolation::makeCorrectionVectors() const
                     (sum(magSf*mag(corrVecs))/sum(magSf)).value(),
                     1.0
                 )
-            )*180.0/mathematicalConstant::pi;
+            )*180.0/constant::math::pi;
     }
 
     if (debug)

src/lagrangian/basic/Cloud/CloudIO.C

@@ -139,7 +139,13 @@ void Foam::Cloud<ParticleType>::readFields()
 
 template<class ParticleType>
 void Foam::Cloud<ParticleType>::writeFields() const
-{}
+{
+    if (this->size())
+    {
+        const ParticleType& p = *this->first();
+        ParticleType::writeFields(p.cloud());
+    }
+}
 
 
 template<class ParticleType>

src/lagrangian/basic/indexedParticle/indexedParticle.H

@@ -47,7 +47,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class indexedParticle Declaration
+                      Class indexedParticle Declaration
 \*---------------------------------------------------------------------------*/
 
 class indexedParticle

src/lagrangian/basic/indexedParticle/indexedParticleCloud.C

@@ -45,20 +45,17 @@ defineTemplateTypeNameAndDebug(Cloud<indexedParticle>, 0);
 Foam::indexedParticleCloud::indexedParticleCloud
 (
     const polyMesh& mesh,
-    const word& cloudName
+    const word& cloudName,
+    bool readFields
 )
 :
     Cloud<indexedParticle>(mesh, cloudName, false)
 {
-    indexedParticle::readFields(*this);
+    if (readFields)
+    {
+        indexedParticle::readFields(*this);
+    }
 }
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::indexedParticleCloud::writeFields() const
-{
-    indexedParticle::writeFields(*this);
-}
-
 // ************************************************************************* //

src/lagrangian/basic/indexedParticle/indexedParticleCloud.H

@@ -45,7 +45,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class indexedParticleCloud Declaration
+                    Class indexedParticleCloud Declaration
 \*---------------------------------------------------------------------------*/
 
 class indexedParticleCloud
@@ -69,14 +69,9 @@ public:
         indexedParticleCloud
         (
             const polyMesh&,
-            const word& cloudName = "defaultCloud"
+            const word& cloudName = "defaultCloud",
+            bool readFields = true
         );
-
-
-    // Member Functions
-
-        //- Write fields
-        virtual void writeFields() const;
 };
 
 

src/lagrangian/basic/passiveParticle/passiveParticleCloud.C

@@ -43,12 +43,16 @@ defineTemplateTypeNameAndDebug(Cloud<passiveParticle>, 0);
 Foam::passiveParticleCloud::passiveParticleCloud
 (
     const polyMesh& mesh,
-    const word& cloudName
+    const word& cloudName,
+    bool readFields
 )
 :
     Cloud<passiveParticle>(mesh, cloudName, false)
 {
-    readFields();
+    if (readFields)
+    {
+        passiveParticle::readFields(*this);
+    }
 }
 
 
@@ -63,18 +67,4 @@ Foam::passiveParticleCloud::passiveParticleCloud
 {}
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::passiveParticleCloud::readFields()
-{
-    passiveParticle::readFields(*this);
-}
-
-
-void Foam::passiveParticleCloud::writeFields() const
-{
-    passiveParticle::writeFields(*this);
-}
-
-
 // ************************************************************************* //

src/lagrangian/basic/passiveParticle/passiveParticleCloud.H

@@ -45,7 +45,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class passiveParticleCloud Declaration
+                   Class passiveParticleCloud Declaration
 \*---------------------------------------------------------------------------*/
 
 class passiveParticleCloud
@@ -72,7 +72,8 @@ public:
         passiveParticleCloud
         (
             const polyMesh&,
-            const word& cloudName = "defaultCloud"
+            const word& cloudName = "defaultCloud",
+            bool readFields = true
         );
 
         //- Construct from mesh, cloud name, and a list of particles
@@ -82,15 +83,6 @@ public:
             const word& cloudName,
             const IDLList<passiveParticle>& particles
         );
-
-
-    // Member Functions
-
-        //- Read fields
-        virtual void readFields();
-
-        //- Write fields
-        virtual void writeFields() const;
 };
 
 

src/lagrangian/coalCombustion/CoalCloud/CoalCloud.C

@@ -1,70 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation; either version 2 of the License, or (at your
-    option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "CoalCloud.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-template<class ThermoType>
-Foam::CoalCloud<ThermoType>::CoalCloud
-(
-    const word& cloudName,
-    const volScalarField& rho,
-    const volVectorField& U,
-    const dimensionedVector& g,
-    basicThermo& thermo
-)
-:
-    ReactingMultiphaseCloud<CoalParcel<ThermoType> >
-    (
-        cloudName,
-        rho,
-        U,
-        g,
-        thermo
-    )
-{
-    CoalParcel<ThermoType>::readFields(*this);
-}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class ThermoType>
-Foam::CoalCloud<ThermoType>::~CoalCloud()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class ThermoType>
-void Foam::CoalCloud<ThermoType>::writeFields() const
-{
-    CoalParcel<ThermoType>::writeFields(*this);
-}
-
-
-// ************************************************************************* //

src/lagrangian/coalCombustion/CoalCloud/CoalCloud.H

@@ -26,9 +26,7 @@ Class
     CoalCloud
 
 Description
-
-SourceFiles
-    CoalCloud.C
+    Coal cloud templated on the type of carrier phase thermodynamics
 
 \*---------------------------------------------------------------------------*/
 
@@ -37,69 +35,21 @@ SourceFiles
 
 #include "ReactingMultiphaseCloud.H"
 #include "CoalParcel.H"
+#include "thermoPhysicsTypes.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
+    typedef ReactingMultiphaseCloud<CoalParcel<constGasThermoPhysics> >
+        constThermoCoalCloud;
 
-/*---------------------------------------------------------------------------*\
-                          Class CoalCloud Declaration
-\*---------------------------------------------------------------------------*/
+    typedef ReactingMultiphaseCloud<CoalParcel<gasThermoPhysics> >
+        thermoCoalCloud;
 
-template<class ThermoType>
-class CoalCloud
-:
-    public ReactingMultiphaseCloud<CoalParcel<ThermoType> >
-{
-    // Private Member Functions
-
-        //- Disallow default bitwise copy construct
-        CoalCloud(const CoalCloud&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const CoalCloud&);
-
-
-public:
-
-    //-Runtime type information
-    TypeName("CoalCloud");
-
-
-    // Constructors
-
-        //- Construct given carrier gas fields
-        CoalCloud
-        (
-            const word& cloudName,
-            const volScalarField& rho,
-            const volVectorField& U,
-            const dimensionedVector& g,
-            basicThermo& thermo
-        );
-
-
-    //- Destructor
-    ~CoalCloud();
-
-
-    // Member Functions
-
-        //- Write fields
-        virtual void writeFields() const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#ifdef NoRepository
-    #include "CoalCloud.C"
-#endif
+    typedef ReactingMultiphaseCloud<CoalParcel<icoPoly8ThermoPhysics> >
+        icoPoly8ThermoCoalCloud;
+}
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/lagrangian/coalCombustion/Make/files

@@ -1,6 +1,3 @@
-/* Coal cloud */
-CoalCloud/defineCoalCloud.C
-
 /* Coal parcel and sub-models */
 CoalParcel/defineCoalParcel.C
 CoalParcel/makeCoalParcelSubmodels.C

src/lagrangian/coalCombustion/include/createCoalParcelTypes.H

@@ -92,8 +92,7 @@ License
     (                                                                         \
         ReactingMultiphaseCloud<ParcelType##ThermoType>,                      \
         0                                                                     \
-    );                                                                        \
-    defineParcelTypeNameAndDebug(CoalCloud<ParcelType##ThermoType>, 0);
+    );
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationDiffusionLimitedRate/COxidationDiffusionLimitedRate.C

@@ -25,6 +25,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "COxidationDiffusionLimitedRate.H"
+#include "mathConstants.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -102,7 +103,7 @@ Foam::scalar Foam::COxidationDiffusionLimitedRate<CloudType>::calculate
     const scalarField& YLiquid,
     const scalarField& YSolid,
     const scalarField& YMixture,
-    const scalarField& dMassVolatile,
+    const scalar N,
     scalarField& dMassGas,
     scalarField& dMassLiquid,
     scalarField& dMassSolid,
@@ -123,8 +124,7 @@ Foam::scalar Foam::COxidationDiffusionLimitedRate<CloudType>::calculate
     const scalar YO2 = this->owner().mcCarrierThermo().Y(O2GlobalId_)[cellI];
 
     // Change in C mass [kg]
-    scalar dmC =
-        4.0*mathematicalConstant::pi*d*D_*YO2*Tc*rhoc/(Sb_*(T + Tc))*dt;
+    scalar dmC = 4.0*constant::math::pi*d*D_*YO2*Tc*rhoc/(Sb_*(T + Tc))*dt;
 
     // Limit mass transfer by availability of C
     dmC = min(mass*fComb, dmC);

src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationDiffusionLimitedRate/COxidationDiffusionLimitedRate.H

@@ -131,7 +131,7 @@ public:
             const scalarField& YLiquid,
             const scalarField& YSolid,
             const scalarField& YMixture,
-            const scalarField& dMassVolatile,
+            const scalar N,
             scalarField& dMassGas,
             scalarField& dMassLiquid,
             scalarField& dMassSolid,

src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationKineticDiffusionLimitedRate/COxidationKineticDiffusionLimitedRate.C

@@ -25,6 +25,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "COxidationKineticDiffusionLimitedRate.H"
+#include "mathConstants.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -110,7 +111,7 @@ Foam::scalar Foam::COxidationKineticDiffusionLimitedRate<CloudType>::calculate
     const scalarField& YLiquid,
     const scalarField& YSolid,
     const scalarField& YMixture,
-    const scalarField& dMassVolatile,
+    const scalar N,
     scalarField& dMassGas,
     scalarField& dMassLiquid,
     scalarField& dMassSolid,
@@ -137,7 +138,7 @@ Foam::scalar Foam::COxidationKineticDiffusionLimitedRate<CloudType>::calculate
     const scalar Rk = C2_*exp(-E_/(specie::RR*Tc));
 
     // Particle surface area
-    const scalar Ap = mathematicalConstant::pi*sqr(d);
+    const scalar Ap = constant::math::pi*sqr(d);
 
     // Change in C mass [kg]
     scalar dmC = Ap*rhoc*specie::RR*Tc*YO2/WO2_*D0*Rk/(D0 + Rk);

src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationKineticDiffusionLimitedRate/COxidationKineticDiffusionLimitedRate.H

@@ -139,7 +139,7 @@ public:
             const scalarField& YLiquid,
             const scalarField& YSolid,
             const scalarField& YMixture,
-            const scalarField& dMassVolatile,
+            const scalar N,
             scalarField& dMassGas,
             scalarField& dMassLiquid,
             scalarField& dMassSolid,

src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationMurphyShaddix/COxidationMurphyShaddix.C

@@ -25,6 +25,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "COxidationMurphyShaddix.H"
+#include "mathConstants.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -106,7 +107,7 @@ Foam::scalar Foam::COxidationMurphyShaddix<CloudType>::calculate
     const scalarField& YLiquid,
     const scalarField& YSolid,
     const scalarField& YMixture,
-    const scalarField& dMassVolatile,
+    const scalar N,
     scalarField& dMassGas,
     scalarField& dMassLiquid,
     scalarField& dMassSolid,
@@ -133,7 +134,7 @@ Foam::scalar Foam::COxidationMurphyShaddix<CloudType>::calculate
     }
 
     // Particle surface area [m^2]
-    const scalar Ap = mathematicalConstant::pi*sqr(d);
+    const scalar Ap = constant::math::pi*sqr(d);
 
     // Calculate diffision constant at continuous phase temperature
     // and density [m^2/s]
@@ -142,9 +143,6 @@ Foam::scalar Foam::COxidationMurphyShaddix<CloudType>::calculate
     // Far field partial pressure O2 [Pa]
     const scalar ppO2 = rhoO2/WO2_*specie::RR*Tc;
 
-    // Molar emission rate of volatiles per unit surface area
-    const scalar qVol = sum(dMassVolatile)/(WVol_*Ap);
-
     // Total molar concentration of the carrier phase [kmol/m^3]
     const scalar C = pc/(specie::RR*Tc);
 
@@ -173,7 +171,7 @@ Foam::scalar Foam::COxidationMurphyShaddix<CloudType>::calculate
     while ((mag(qCs - qCsOld)/qCs > tolerance_) && (iter <= maxIters_))
     {
         qCsOld = qCs;
-        const scalar PO2Surface = ppO2*exp(-(qCs + qVol)*d/(2*C*D));
+        const scalar PO2Surface = ppO2*exp(-(qCs + N)*d/(2*C*D));
         qCs = A_*exp(-E_/(specie::RR*T))*pow(PO2Surface, n_);
         qCs = (100.0*qCs + iter*qCsOld)/(100.0 + iter);
         qCs = min(qCs, qCsLim);

src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationMurphyShaddix/COxidationMurphyShaddix.H

@@ -152,7 +152,7 @@ public:
             const scalarField& YLiquid,
             const scalarField& YSolid,
             const scalarField& YMixture,
-            const scalarField& dMassVolatile,
+            const scalar N,
             scalarField& dMassGas,
             scalarField& dMassLiquid,
             scalarField& dMassSolid,

src/lagrangian/dieselSpray/injector/commonRailInjector/commonRailInjector.C

@@ -27,22 +27,20 @@ License
 #include "commonRailInjector.H"
 #include "addToRunTimeSelectionTable.H"
 #include "Random.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 namespace Foam
 {
+    defineTypeNameAndDebug(commonRailInjector, 0);
 
-defineTypeNameAndDebug(commonRailInjector, 0);
-
-addToRunTimeSelectionTable
-(
-    injectorType,
-    commonRailInjector,
-    dictionary
-);
+    addToRunTimeSelectionTable
+    (
+        injectorType,
+        commonRailInjector,
+        dictionary
+    );
 }
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -71,38 +69,58 @@ Foam::commonRailInjector::commonRailInjector
     averageParcelMass_(mass_/nParcels_),
     pressureIndependentVelocity_(false)
 {
-
     // convert CA to real time
     forAll(massFlowRateProfile_, i)
     {
-        massFlowRateProfile_[i][0] = t.userTimeToTime(massFlowRateProfile_[i][0]);
+        massFlowRateProfile_[i][0] =
+            t.userTimeToTime(massFlowRateProfile_[i][0]);
         velocityProfile_[i][0] = t.userTimeToTime(massFlowRateProfile_[i][0]);
     }
 
     forAll(injectionPressureProfile_, i)
     {
-        injectionPressureProfile_[i][0] = t.userTimeToTime(injectionPressureProfile_[i][0]);
+        injectionPressureProfile_[i][0] =
+            t.userTimeToTime(injectionPressureProfile_[i][0]);
     }
 
-    if (mag(injectionPressureProfile_[0][0]-massFlowRateProfile_[0][0]) > SMALL)
+    if
+    (
+        mag(injectionPressureProfile_[0][0] - massFlowRateProfile_[0][0])
+      > SMALL
+    )
     {
-        FatalError << "commonRailInjector::commonRailInjector(const time& t, const dictionary dict) " << endl
-                << " start-time entries for injectionPressureProfile and massFlowRateProfile do no match"
-                << abort(FatalError);
+        FatalErrorIn
+        (
+            "commonRailInjector::commonRailInjector"
+            "(const time& t, const dictionary dict)"
+        )   << " start-time entries for injectionPressureProfile and "
+            << "massFlowRateProfile do no match"
+            << abort(FatalError);
     }
-    Info << "injectionPressureProfile_.size() = " << injectionPressureProfile_.size()
+    Info<< "injectionPressureProfile_.size() = "
+        << injectionPressureProfile_.size()
         << ", massFlowRateProfile_.size() = " << massFlowRateProfile_.size()
         << endl;
 
-    if (mag(injectionPressureProfile_[injectionPressureProfile_.size()-1][0]-massFlowRateProfile_[massFlowRateProfile_.size()-1][0]) > SMALL)
+    if
+    (
+        mag(injectionPressureProfile_[injectionPressureProfile_.size()-1][0]
+      - massFlowRateProfile_[massFlowRateProfile_.size()-1][0])
+      > SMALL
+    )
     {
-        FatalError << "commonRailInjector::commonRailInjector(const time& t, const dictionary dict) " << endl
-                << " end-time entries for injectionPressureProfile and massFlowRateProfile do no match"
-                << abort(FatalError);
+        FatalErrorIn
+        (
+            "commonRailInjector::commonRailInjector"
+            "(const time& t, const dictionary dict)"
+        )   << "End-time entries for injectionPressureProfile and "
+            << "massFlowRateProfile do no match"
+            << abort(FatalError);
     }
 
     scalar integratedMFR = integrateTable(massFlowRateProfile_);
-    scalar integratedP = integrateTable(injectionPressureProfile_)/(teoi()-tsoi());
+    scalar integratedP =
+        integrateTable(injectionPressureProfile_)/(teoi() - tsoi());
 
     forAll(massFlowRateProfile_, i)
     {
@@ -113,7 +131,6 @@ Foam::commonRailInjector::commonRailInjector
         TProfile_[i][1] = T_;
 
         CdProfile_[i][0] = massFlowRateProfile_[i][0];
-
     }
 
     forAll(injectionPressureProfile_, i)
@@ -122,7 +139,7 @@ Foam::commonRailInjector::commonRailInjector
     }
     // Normalize the direction vector
     direction_ /= mag(direction_);
-    
+
     setTangentialVectors();
 
     // check molar fractions
@@ -134,18 +151,22 @@ Foam::commonRailInjector::commonRailInjector
 
     if (mag(Xsum - 1.0) > SMALL)
     {
-        Info << "Warning!!!\n commonRailInjector::commonRailInjector(const time& t, Istream& is)"
-            << "X does not add up to 1.0, correcting molar fractions."
+        WarningIn
+        (
+            "commonRailInjector::commonRailInjector"
+            "(const time& t, const dictionary dict)"
+        )   << "X does not add up to 1.0, correcting molar fractions."
             << endl;
         forAll(X_, i)
         {
             X_[i] /= Xsum;
         }
     }
-    Info << "end constructor. in commonRail" << endl;
 
+    Info << "end constructor. in commonRail" << endl;
 }
 
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::commonRailInjector::~commonRailInjector()
@@ -181,17 +202,20 @@ Foam::label Foam::commonRailInjector::nParcelsToInject
 ) const
 {
 
-    scalar mInj = mass_*(fractionOfInjection(time1)-fractionOfInjection(time0));
+    scalar mInj =
+        mass_*(fractionOfInjection(time1) - fractionOfInjection(time0));
     label nParcels = label(mInj/averageParcelMass_ + 0.49);
-    
+
     return nParcels;
 }
 
+
 const Foam::vector Foam::commonRailInjector::position(const label n) const
 {
     return position_;
 }
 
+
 Foam::vector Foam::commonRailInjector::position
 (
     const label n,
@@ -220,33 +244,35 @@ Foam::vector Foam::commonRailInjector::position
     {
         // otherwise, disc injection
         scalar iRadius = d_*rndGen.scalar01();
-        scalar iAngle = 2.0*mathematicalConstant::pi*rndGen.scalar01();
+        scalar iAngle = constant::math::twoPi*rndGen.scalar01();
 
         return
-        ( 
+        (
             position_
           + iRadius
-          * (
-              tangentialInjectionVector1_*cos(iAngle)
-            + tangentialInjectionVector2_*sin(iAngle)
-          )
+           *(
+                tangentialInjectionVector1_*cos(iAngle)
+              + tangentialInjectionVector2_*sin(iAngle)
+            )
         );
-        
     }
 
     return position_;
 }
 
+
 Foam::label Foam::commonRailInjector::nHoles() const
 {
     return 1;
 }
 
+
 Foam::scalar Foam::commonRailInjector::d() const
 {
     return d_;
 }
 
+
 const Foam::vector& Foam::commonRailInjector::direction
 (
     const label i,
@@ -256,6 +282,7 @@ const Foam::vector& Foam::commonRailInjector::direction
     return direction_;
 }
 
+
 Foam::scalar Foam::commonRailInjector::mass
 (
     const scalar time0,
@@ -264,47 +291,55 @@ Foam::scalar Foam::commonRailInjector::mass
     const scalar angleOfWedge
 ) const
 {
-    scalar mInj = mass_*(fractionOfInjection(time1)-fractionOfInjection(time0));
+    scalar mInj =
+        mass_*(fractionOfInjection(time1) - fractionOfInjection(time0));
 
-    // correct mass if calculation is 2D 
+    // correct mass if calculation is 2D
     if (twoD)
     {
-        mInj *= 0.5*angleOfWedge/mathematicalConstant::pi;
+        mInj *= 0.5*angleOfWedge/constant::math::pi;
     }
 
     return mInj;
 }
 
+
 Foam::scalar Foam::commonRailInjector::mass() const
 {
     return mass_;
 }
 
+
 const Foam::scalarField& Foam::commonRailInjector::X() const
 {
     return X_;
 }
 
+
 Foam::List<Foam::commonRailInjector::pair> Foam::commonRailInjector::T() const
 {
     return TProfile_;
 }
 
+
 Foam::scalar Foam::commonRailInjector::T(const scalar time) const
 {
     return T_;
 }
 
+
 Foam::scalar Foam::commonRailInjector::tsoi() const
 {
     return massFlowRateProfile_[0][0];
 }
 
+
 Foam::scalar Foam::commonRailInjector::teoi() const
 {
     return massFlowRateProfile_[massFlowRateProfile_.size()-1][0];
 }
 
+
 Foam::scalar Foam::commonRailInjector::massFlowRate
 (
     const scalar time
@@ -313,6 +348,7 @@ Foam::scalar Foam::commonRailInjector::massFlowRate
     return getTableValue(massFlowRateProfile_, time);
 }
 
+
 Foam::scalar Foam::commonRailInjector::injectionPressure
 (
     const scalar time
@@ -321,6 +357,7 @@ Foam::scalar Foam::commonRailInjector::injectionPressure
     return getTableValue(injectionPressureProfile_, time);
 }
 
+
 Foam::scalar Foam::commonRailInjector::velocity
 (
     const scalar time
@@ -329,11 +366,14 @@ Foam::scalar Foam::commonRailInjector::velocity
     return getTableValue(velocityProfile_, time);
 }
 
-Foam::List<Foam::commonRailInjector::pair> Foam::commonRailInjector::CdProfile() const
+
+Foam::List<Foam::commonRailInjector::pair> Foam::commonRailInjector::CdProfile()
+const
 {
     return CdProfile_;
 }
 
+
 Foam::scalar Foam::commonRailInjector::Cd
 (
     const scalar time
@@ -342,11 +382,13 @@ Foam::scalar Foam::commonRailInjector::Cd
     return getTableValue(CdProfile_, time);
 }
 
+
 Foam::scalar Foam::commonRailInjector::fractionOfInjection(const scalar time) const
 {
     return integrateTable(massFlowRateProfile_, time)/mass_;
 }
 
+
 Foam::scalar Foam::commonRailInjector::injectedMass
 (
     const scalar t
@@ -362,13 +404,14 @@ void Foam::commonRailInjector::correctProfiles
     const scalar referencePressure
 )
 {
-    scalar A = 0.25*mathematicalConstant::pi*pow(d_, 2.0);
+    scalar A = 0.25*constant::math::pi*sqr(d_);
     scalar pDummy = 1.0e+5;
     scalar rho = fuel.rho(pDummy, T_, X_);
 
     forAll(velocityProfile_, i)
     {
-        scalar Pinj = getTableValue(injectionPressureProfile_, velocityProfile_[i][0]);
+        scalar Pinj =
+            getTableValue(injectionPressureProfile_, velocityProfile_[i][0]);
         scalar Vinj = sqrt(2.0*(Pinj - referencePressure)/rho);
         scalar mfr = massFlowRateProfile_[i][1]/(rho*A);
         scalar Cd = mfr/Vinj;
@@ -377,14 +420,17 @@ void Foam::commonRailInjector::correctProfiles
     }
 }
 
+
 Foam::vector Foam::commonRailInjector::tan1(const label n) const
 {
     return tangentialInjectionVector1_;
 }
 
+
 Foam::vector Foam::commonRailInjector::tan2(const label n) const
 {
     return tangentialInjectionVector2_;
 }
 
+
 // ************************************************************************* //

src/lagrangian/dieselSpray/injector/commonRailInjector/commonRailInjector.H

@@ -47,7 +47,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class commonRailInjector Declaration
+                     Class commonRailInjector Declaration
 \*---------------------------------------------------------------------------*/
 
 class commonRailInjector
@@ -110,16 +110,11 @@ public:
     // Constructors
 
         //- Construct from components
-        commonRailInjector
-        (
-            const Time& t,
-            const dictionary& dict
-        );
+        commonRailInjector(const Time& t, const dictionary& dict);
 
 
-    // Destructor
-
-        ~commonRailInjector();
+    //- Destructor
+    virtual ~commonRailInjector();
 
 
     // Member Functions
@@ -146,7 +141,7 @@ public:
             const vector& axisOfWedgeNormal,
             Random& rndGen
         ) const;
-    
+
         //- Return the number of holes
         label nHoles() const;
 
@@ -201,14 +196,14 @@ public:
         {
             return injectionPressureProfile_;
         }
-    
+
         scalar injectionPressure(const scalar time) const;
 
         List<pair> velocityProfile() const
         {
             return velocityProfile_;
         }
-        
+
         scalar velocity(const scalar time) const;
 
         List<pair> CdProfile() const;
@@ -227,7 +222,6 @@ public:
         {
             return pressureIndependentVelocity_;
         }
-
 };
 
 

src/lagrangian/dieselSpray/injector/definedInjector/definedInjector.C

@@ -27,21 +27,19 @@ License
 #include "definedInjector.H"
 #include "addToRunTimeSelectionTable.H"
 #include "Random.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 namespace Foam
 {
+    defineTypeNameAndDebug(definedInjector, 0);
 
-defineTypeNameAndDebug(definedInjector, 0);
-
-addToRunTimeSelectionTable
-(
-    injectorType,
-    definedInjector,
-    dictionary
-);
-
+    addToRunTimeSelectionTable
+    (
+        injectorType,
+        definedInjector,
+        dictionary
+    );
 }
 
 
@@ -73,12 +71,13 @@ Foam::definedInjector::definedInjector
     // convert CA to real time - mass flow rate profile
     forAll(massFlowRateProfile_, i)
     {
-        massFlowRateProfile_[i][0] = t.userTimeToTime(massFlowRateProfile_[i][0]);
-	    // dummy
-	    injectionPressureProfile_[i][0] = massFlowRateProfile_[i][0];
-	    injectionPressureProfile_[i][1] = 0.0;
-	    CdProfile_[i][0] = massFlowRateProfile_[i][0];
-	    CdProfile_[i][1] = 1.0;
+        massFlowRateProfile_[i][0] =
+            t.userTimeToTime(massFlowRateProfile_[i][0]);
+        // dummy
+        injectionPressureProfile_[i][0] = massFlowRateProfile_[i][0];
+        injectionPressureProfile_[i][1] = 0.0;
+        CdProfile_[i][0] = massFlowRateProfile_[i][0];
+        CdProfile_[i][1] = 1.0;
     }
 
     forAll(velocityProfile_, i)
@@ -89,15 +88,28 @@ Foam::definedInjector::definedInjector
     // check if time entries match
     if (mag(massFlowRateProfile_[0][0]-velocityProfile_[0][0]) > SMALL)
     {
-        FatalError << "definedInjector::definedInjector(const time& t, const dictionary dict) " << endl
-            << " start-times do not match for velocityProfile and massFlowRateProfile."
+        FatalErrorIn
+        (
+            "definedInjector::definedInjector"
+            "(const time& t, const dictionary dict)"
+        )   << "Start-times do not match for velocityProfile and "
+            << "massFlowRateProfile." << nl
             << abort(FatalError);
     }
 
-    if (mag(massFlowRateProfile_[massFlowRateProfile_.size()-1][0]-velocityProfile_[velocityProfile_.size()-1][0]) > SMALL)
+    if
+    (
+        mag(massFlowRateProfile_[massFlowRateProfile_.size()-1][0]
+      - velocityProfile_[velocityProfile_.size()-1][0])
+      > SMALL
+    )
     {
-        FatalError << "definedInjector::definedInjector(const time& t, const dictionary dict) " << endl
-            << " end-times do not match for velocityProfile and massFlowRateProfile."
+        FatalErrorIn
+        (
+            "definedInjector::definedInjector"
+            "(const time& t, const dictionary dict)"
+        )   << "End-times do not match for velocityProfile and "
+            << "massFlowRateProfile."
             << abort(FatalError);
     }
 
@@ -112,7 +124,7 @@ Foam::definedInjector::definedInjector
 
     // Normalize the direction vector
     direction_ /= mag(direction_);
-    
+
     setTangentialVectors();
 
     // check molar fractions
@@ -126,7 +138,8 @@ Foam::definedInjector::definedInjector
     {
         WarningIn
         (
-            "definedInjector::definedInjector(const time& t, Istream& is)"
+            "definedInjector::definedInjector"
+            "(const time& t, const dictionary dict)"
         )   << "X does not add up to 1.0, correcting molar fractions."
             << endl;
 
@@ -174,15 +187,17 @@ Foam::label Foam::definedInjector::nParcelsToInject
 
     scalar mInj = mass_*(fractionOfInjection(time1)-fractionOfInjection(time0));
     label nParcels = label(mInj/averageParcelMass_ + 0.49);
-    
+
     return nParcels;
 }
 
+
 const Foam::vector Foam::definedInjector::position(const label n) const
 {
     return position_;
 }
 
+
 Foam::vector Foam::definedInjector::position
 (
     const label n,
@@ -211,10 +226,10 @@ Foam::vector Foam::definedInjector::position
     {
         // otherwise, disc injection
         scalar iRadius = d_*rndGen.scalar01();
-        scalar iAngle = 2.0*mathematicalConstant::pi*rndGen.scalar01();
+        scalar iAngle = constant::math::twoPi*rndGen.scalar01();
 
         return
-        ( 
+        (
             position_
           + iRadius
           * (
@@ -222,22 +237,25 @@ Foam::vector Foam::definedInjector::position
             + tangentialInjectionVector2_*sin(iAngle)
           )
         );
-        
+
     }
 
     return position_;
 }
 
+
 Foam::label Foam::definedInjector::nHoles() const
 {
     return 1;
 }
 
+
 Foam::scalar Foam::definedInjector::d() const
 {
     return d_;
 }
 
+
 const Foam::vector& Foam::definedInjector::direction
 (
     const label i,
@@ -247,6 +265,7 @@ const Foam::vector& Foam::definedInjector::direction
     return direction_;
 }
 
+
 Foam::scalar Foam::definedInjector::mass
 (
     const scalar time0,
@@ -255,62 +274,73 @@ Foam::scalar Foam::definedInjector::mass
     const scalar angleOfWedge
 ) const
 {
-    scalar mInj = mass_*(fractionOfInjection(time1)-fractionOfInjection(time0));
+    scalar mInj =
+        mass_*(fractionOfInjection(time1) - fractionOfInjection(time0));
 
-    // correct mass if calculation is 2D 
+    // correct mass if calculation is 2D
     if (twoD)
     {
-        mInj *= 0.5*angleOfWedge/mathematicalConstant::pi;
+        mInj *= 0.5*angleOfWedge/constant::math::pi;
     }
 
     return mInj;
 }
 
+
 Foam::scalar Foam::definedInjector::mass() const
 {
     return mass_;
 }
 
+
 Foam::scalar Foam::definedInjector::massFlowRate(const scalar time) const
 {
     return getTableValue(massFlowRateProfile_, time);
 }
 
+
 Foam::scalar Foam::definedInjector::injectionPressure(const scalar time) const
 {
     return getTableValue(injectionPressureProfile_, time);
 }
 
+
 Foam::scalar Foam::definedInjector::Cd(const scalar time) const
 {
     return getTableValue(CdProfile_, time);
 }
 
+
 const Foam::scalarField& Foam::definedInjector::X() const
 {
     return X_;
 }
 
+
 Foam::List<Foam::definedInjector::pair> Foam::definedInjector::T() const
 {
     return TProfile_;
 }
 
+
 Foam::scalar Foam::definedInjector::T(const scalar time) const
 {
     return T_;
 }
 
+
 Foam::scalar Foam::definedInjector::tsoi() const
 {
     return massFlowRateProfile_[0][0];
 }
 
+
 Foam::scalar Foam::definedInjector::teoi() const
 {
     return massFlowRateProfile_[massFlowRateProfile_.size()-1][0];
 }
 
+
 Foam::scalar Foam::definedInjector::fractionOfInjection
 (
     const scalar time
@@ -319,6 +349,7 @@ Foam::scalar Foam::definedInjector::fractionOfInjection
     return integrateTable(massFlowRateProfile_, time)/mass_;
 }
 
+
 Foam::scalar Foam::definedInjector::velocity
 (
     const scalar time
@@ -327,6 +358,7 @@ Foam::scalar Foam::definedInjector::velocity
     return getTableValue(velocityProfile_, time);
 }
 
+
 Foam::scalar Foam::definedInjector::injectedMass
 (
     const scalar t
@@ -335,13 +367,14 @@ Foam::scalar Foam::definedInjector::injectedMass
     return mass_*fractionOfInjection(t);
 }
 
+
 void Foam::definedInjector::correctProfiles
 (
     const liquidMixture& fuel,
     const scalar referencePressure
 )
 {
-    scalar A = 0.25*mathematicalConstant::pi*pow(d_, 2.0);
+    scalar A = 0.25*constant::math::pi*sqr(d_);
     scalar pDummy = 1.0e+5;
     scalar rho = fuel.rho(pDummy, T_, X_);
 
@@ -354,14 +387,17 @@ void Foam::definedInjector::correctProfiles
     }
 }
 
+
 Foam::vector Foam::definedInjector::tan1(const label n) const
 {
     return tangentialInjectionVector1_;
 }
 
+
 Foam::vector Foam::definedInjector::tan2(const label n) const
 {
     return tangentialInjectionVector2_;
 }
 
+
 // ************************************************************************* //

src/lagrangian/dieselSpray/injector/definedInjector/definedInjector.H

@@ -46,7 +46,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class definedInjector Declaration
+                      Class definedInjector Declaration
 \*---------------------------------------------------------------------------*/
 
 class definedInjector
@@ -69,8 +69,8 @@ private:
         scalar T_;
         label nParcels_;
         scalarField X_;
-        List<pair> massFlowRateProfile_;                         // aj
-        List<pair> velocityProfile_;                    // aj
+        List<pair> massFlowRateProfile_;
+        List<pair> velocityProfile_;
         List<pair> injectionPressureProfile_;
         List<pair> CdProfile_;
         List<pair> TProfile_;
@@ -101,6 +101,7 @@ private:
         //- Return the instantaneous injection velocity
         scalar injectionVelocity(const scalar) const;
 
+
 public:
 
     //- Runtime type information
@@ -110,16 +111,11 @@ public:
     // Constructors
 
         //- Construct from components
-        definedInjector
-        (
-            const Time& t,
-            const dictionary& dict
-        );
+        definedInjector(const Time& t, const dictionary& dict);
 
 
-    // Destructor
-
-        ~definedInjector();
+    //- Destructor
+    virtual ~definedInjector();
 
 
     // Member Functions
@@ -146,7 +142,7 @@ public:
             const vector& axisOfWedgeNormal,
             Random& rndGen
         ) const;
-    
+
         //- Return the number of holes
         label nHoles() const;
 
@@ -201,14 +197,14 @@ public:
         {
             return injectionPressureProfile_;
         }
-    
+
         scalar injectionPressure(const scalar time) const;
 
         List<pair> velocityProfile() const
         {
             return velocityProfile_;
         }
-        
+
         scalar velocity(const scalar time) const;
 
         List<pair> CdProfile() const
@@ -231,7 +227,6 @@ public:
         {
             return pressureIndependentVelocity_;
         }
-
 };
 
 

src/lagrangian/dieselSpray/injector/injector/injector.C

@@ -31,10 +31,9 @@ License
 namespace Foam
 {
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+// No code
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/lagrangian/dieselSpray/injector/injector/injector.H

@@ -47,7 +47,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class injector Declaration
+                          Class injector Declaration
 \*---------------------------------------------------------------------------*/
 
 class injector
@@ -62,7 +62,7 @@ class injector
 public:
 
     // Constructors
-    
+
         //- Constructor from Istream
         injector(const Time& t, Istream& is);
 

src/lagrangian/dieselSpray/injector/injectorType/injectorType.C

@@ -30,10 +30,8 @@ License
 
 namespace Foam
 {
-
-defineTypeNameAndDebug(injectorType, 0);
-defineRunTimeSelectionTable(injectorType, dictionary);
-
+    defineTypeNameAndDebug(injectorType, 0);
+    defineRunTimeSelectionTable(injectorType, dictionary);
 }
 
 
@@ -47,6 +45,7 @@ Foam::injectorType::injectorType
 )
 {}
 
+
 // * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
 
 Foam::autoPtr<Foam::injectorType> Foam::injectorType::New
@@ -55,30 +54,25 @@ Foam::autoPtr<Foam::injectorType> Foam::injectorType::New
     const dictionary& dict
 )
 {
-    word injectorTypeName
-    (
-        dict.lookup("injectorType")
-    );
+    word injectorTypeName(dict.lookup("injectorType"));
 
-    Info<< "Selecting injectorType "
-         << injectorTypeName << endl;
+    Info<< "Selecting injectorType " << injectorTypeName << endl;
 
     dictionaryConstructorTable::iterator cstrIter =
         dictionaryConstructorTablePtr_->find(injectorTypeName);
 
     if (cstrIter == dictionaryConstructorTablePtr_->end())
     {
-        FatalError
-            << "injectorType::New(const dictionary&) : " << endl
-            << "    unknown injectorType type "
+        FatalErrorIn("injectorType::New(const dictionary&)")
+            << "Unknown injectorType type "
             << injectorTypeName
-            << ", constructor not in hash table" << endl << endl
-            << "    Valid injector types are :" << endl;
-        Info<< dictionaryConstructorTablePtr_->sortedToc() << abort(FatalError);
+            << ", constructor not in hash table" << nl << nl
+            << "    Valid injector types are:" << nl
+            << dictionaryConstructorTablePtr_->sortedToc()
+            << abort(FatalError);
     }
 
     return autoPtr<injectorType>(cstrIter()(t, dict));
-
 }
 
 
@@ -87,6 +81,7 @@ Foam::autoPtr<Foam::injectorType> Foam::injectorType::New
 Foam::injectorType::~injectorType()
 {}
 
+
 Foam::scalar Foam::injectorType::getTableValue
 (
     const List<pair>& table,
@@ -120,12 +115,13 @@ Foam::scalar Foam::injectorType::getTableValue
             i++;
         }
         // value sits bewteen table[i][0] and table[i+1][0]
-        return table[i][1] 
+        return table[i][1]
                + (value-table[i][0])/(table[i+1][0]-table[i][0])
                * (table[i+1][1]-table[i][1]);
     }
 }
 
+
 Foam::scalar Foam::injectorType::integrateTable
 (
     const List<pair>& table,
@@ -162,6 +158,7 @@ Foam::scalar Foam::injectorType::integrateTable
     return sum;
 }
 
+
 Foam::scalar Foam::injectorType::integrateTable
 (
     const List<pair>& table
@@ -178,4 +175,5 @@ Foam::scalar Foam::injectorType::integrateTable
     return integratedTable;
 }
 
+
 // ************************************************************************* //

src/lagrangian/dieselSpray/injector/injectorType/injectorType.H

@@ -53,12 +53,12 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class injectorType Declaration
+                       Class injectorType Declaration
 \*---------------------------------------------------------------------------*/
 
 class injectorType
 {
-        typedef VectorSpace<Vector<scalar>, scalar, 2> pair;
+    typedef VectorSpace<Vector<scalar>, scalar, 2> pair;
 
 public:
 
@@ -79,15 +79,11 @@ public:
             (t, dict)
         );
 
-    
+
     // Constructors
 
         //- Construct from components
-        injectorType
-        (
-            const Time& t,
-            const dictionary& dict
-        );
+        injectorType(const Time& t, const dictionary& dict);
 
 
     // Selectors
@@ -99,13 +95,12 @@ public:
         );
 
 
-    // Destructor
-
-        virtual ~injectorType();
+    //- Destructor
+    virtual ~injectorType();
 
 
     // Member Functions
-    
+
         //- Return number of particles to inject
         virtual label nParcelsToInject
         (
@@ -131,14 +126,14 @@ public:
 
         //- Return the number of holes
         virtual label nHoles() const = 0;
-    
+
         //- Return the injector diameter
         virtual scalar d() const = 0;
 
         //- Return the injection direction for hole i
         virtual const vector& direction
         (
-            const label i, 
+            const label i,
             const scalar time
         ) const = 0;
 
@@ -172,7 +167,7 @@ public:
         //- Return the instantaneous velocity
         virtual scalar velocity(const scalar time) const = 0;
 
-        //- Return the discharge coefficient        
+        //- Return the discharge coefficient
         virtual List<pair> CdProfile() const = 0;
 
         //- Return the instantaneous discharge coefficient
@@ -197,10 +192,12 @@ public:
 
         virtual bool pressureIndependentVelocity() const = 0;
 
-        //- Return a vector perpendicular to the injection direction and tan2 for hole n
+        //- Return a vector perpendicular to the injection direction and tan2
+        //  for hole n
         virtual vector tan1(const label n) const = 0;
 
-        //- Return a vector perpendicular to the injection direction and tan1 for hole n
+        //- Return a vector perpendicular to the injection direction and tan1
+        //  for hole n
         virtual vector tan2(const label n) const = 0;
 
         scalar getTableValue
@@ -225,7 +222,6 @@ public:
             const liquidMixture& fuel,
             const scalar referencePressure
         ) = 0;
-
 };
 
 

src/lagrangian/dieselSpray/injector/multiHoleInjector/multiHoleInjector.C

@@ -27,22 +27,20 @@ License
 #include "multiHoleInjector.H"
 #include "addToRunTimeSelectionTable.H"
 #include "Random.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 namespace Foam
 {
+    defineTypeNameAndDebug(multiHoleInjector, 0);
 
-defineTypeNameAndDebug(multiHoleInjector, 0);
-
-addToRunTimeSelectionTable
-(
-    injectorType,
-    multiHoleInjector,
-    dictionary
-);
+    addToRunTimeSelectionTable
+    (
+        injectorType,
+        multiHoleInjector,
+        dictionary
+    );
 }
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -79,31 +77,43 @@ Foam::multiHoleInjector::multiHoleInjector
     tangentialInjectionVector1_(nHoles_),
     tangentialInjectionVector2_(nHoles_)
 {
-
-
     // check if time entries for soi and eoi match
-    if (mag(massFlowRateProfile_[0][0]-TProfile_[0][0]) > SMALL)
+    if (mag(massFlowRateProfile_[0][0] - TProfile_[0][0]) > SMALL)
     {
-        FatalError << "multiHoleInjector::multiHoleInjector(const time& t, const dictionary dict) " << endl
-            << " start-times do not match for TemperatureProfile and massFlowRateProfile."
+        FatalErrorIn
+        (
+            "multiHoleInjector::multiHoleInjector"
+            "(const time& t, const dictionary dict)"
+        )   << "Start-times do not match for TemperatureProfile and "
+            << "massFlowRateProfile."
             << abort(FatalError);
     }
 
-    if (mag(massFlowRateProfile_[massFlowRateProfile_.size()-1][0]-TProfile_[TProfile_.size()-1][0]) > SMALL)
+    if
+    (
+        mag(massFlowRateProfile_[massFlowRateProfile_.size()-1][0]
+      - TProfile_[TProfile_.size()-1][0])
+      > SMALL
+    )
     {
-        FatalError << "multiHoleInjector::multiHoleInjector(const time& t, const dictionary dict) " << endl
-            << " end-times do not match for TemperatureProfile and massFlowRateProfile."
+        FatalErrorIn
+        (
+            "multiHoleInjector::multiHoleInjector"
+            "(const time& t, const dictionary dict)"
+        )   << "End-times do not match for TemperatureProfile and "
+            << "massFlowRateProfile."
             << abort(FatalError);
     }
 
     // convert CA to real time
     forAll(massFlowRateProfile_, i)
     {
-        massFlowRateProfile_[i][0] = t.userTimeToTime(massFlowRateProfile_[i][0]);
+        massFlowRateProfile_[i][0] =
+            t.userTimeToTime(massFlowRateProfile_[i][0]);
         velocityProfile_[i][0] = massFlowRateProfile_[i][0];
         injectionPressureProfile_[i][0] = massFlowRateProfile_[i][0];
     }
-    
+
     forAll(TProfile_, i)
     {
         TProfile_[i][0] = t.userTimeToTime(TProfile_[i][0]);
@@ -115,7 +125,7 @@ Foam::multiHoleInjector::multiHoleInjector
     {
         // correct the massFlowRateProfile to match the injected mass
         massFlowRateProfile_[i][1] *= mass_/integratedMFR;
-        
+
         CdProfile_[i][0] = massFlowRateProfile_[i][0];
         CdProfile_[i][1] = Cd_;
     }
@@ -131,17 +141,20 @@ Foam::multiHoleInjector::multiHoleInjector
 
     if (mag(Xsum - 1.0) > SMALL)
     {
-        Info << "Warning!!!\n multiHoleInjector::multiHoleInjector(const time& t, Istream& is)"
-            << "X does not add up to 1.0, correcting molar fractions."
+        WarningIn
+        (
+            "multiHoleInjector::multiHoleInjector"
+            "(const time& t, const dictionary dict)"
+        )   << "X does not add up to 1.0, correcting molar fractions."
             << endl;
         forAll(X_, i)
         {
             X_[i] /= Xsum;
         }
     }
-
 }
 
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::multiHoleInjector::~multiHoleInjector()
@@ -152,7 +165,7 @@ Foam::multiHoleInjector::~multiHoleInjector()
 
 void Foam::multiHoleInjector::setTangentialVectors()
 {
-    scalar pi180 = mathematicalConstant::pi/180.0;
+    scalar pi180 = constant::math::pi/180.0;
     scalar alpha = xyAngle_*pi180;
     scalar phi = zAngle_*pi180;
 
@@ -172,7 +185,7 @@ void Foam::multiHoleInjector::setTangentialVectors()
 
     scalar angle = 0.0;
     scalar u = umbrellaAngle_*pi180/2.0;
-    for(label i=0; i<nHoles_; i++)
+    for (label i=0; i<nHoles_; i++)
     {
         angle += angleSpacing_[i];
         scalar v = angle*pi180;
@@ -183,27 +196,26 @@ void Foam::multiHoleInjector::setTangentialVectors()
             dp /= mag(dp);
         }
         position_[i] = centerPosition_ + 0.5*nozzleTipDiameter_*dp;
-//        Info << "i = " << i << ", dir = " << direction_[i] << ", pos = " << position_[i] << endl;
     }
 
     Random rndGen(label(0));
 
-    for(label i=0; i<nHoles_; i++)
+    for (label i=0; i<nHoles_; i++)
     {
         vector tangent(vector::zero);
         scalar magV = 0;
         while (magV < SMALL)
         {
             vector testThis = rndGen.vector01();
-            
+
             tangent = testThis - (testThis & direction_[i])*direction_[i];
             magV = mag(tangent);
         }
 
         tangentialInjectionVector1_[i] = tangent/magV;
-        tangentialInjectionVector2_[i] = direction_[i] ^ tangentialInjectionVector1_[i];
-
-    }   
+        tangentialInjectionVector2_[i] =
+            direction_[i] ^ tangentialInjectionVector1_[i];
+    }
 }
 
 
@@ -213,18 +225,20 @@ Foam::label Foam::multiHoleInjector::nParcelsToInject
     const scalar time1
 ) const
 {
-
-    scalar mInj = mass_*(fractionOfInjection(time1)-fractionOfInjection(time0));
+    scalar mInj =
+        mass_*(fractionOfInjection(time1) - fractionOfInjection(time0));
     label nParcels = label(mInj/averageParcelMass_ + 0.49);
-    
+
     return nParcels;
 }
 
+
 const Foam::vector Foam::multiHoleInjector::position(const label n) const
 {
     return position_[n];
 }
 
+
 Foam::vector Foam::multiHoleInjector::position
 (
     const label n,
@@ -253,10 +267,10 @@ Foam::vector Foam::multiHoleInjector::position
     {
         // otherwise, disc injection
         scalar iRadius = d_*rndGen.scalar01();
-        scalar iAngle = 2.0*mathematicalConstant::pi*rndGen.scalar01();
+        scalar iAngle = constant::math::twoPi*rndGen.scalar01();
 
         return
-        ( 
+        (
             position_[n]
           + iRadius
           * (
@@ -264,21 +278,23 @@ Foam::vector Foam::multiHoleInjector::position
             + tangentialInjectionVector2_[n]*sin(iAngle)
           )
         );
-        
     }
     return position_[0];
 }
 
+
 Foam::label Foam::multiHoleInjector::nHoles() const
 {
     return nHoles_;
 }
 
+
 Foam::scalar Foam::multiHoleInjector::d() const
 {
     return d_;
 }
 
+
 const Foam::vector& Foam::multiHoleInjector::direction
 (
     const label i,
@@ -288,6 +304,7 @@ const Foam::vector& Foam::multiHoleInjector::direction
     return direction_[i];
 }
 
+
 Foam::scalar Foam::multiHoleInjector::mass
 (
     const scalar time0,
@@ -296,47 +313,55 @@ Foam::scalar Foam::multiHoleInjector::mass
     const scalar angleOfWedge
 ) const
 {
-    scalar mInj = mass_*(fractionOfInjection(time1)-fractionOfInjection(time0));
+    scalar mInj =
+        mass_*(fractionOfInjection(time1) - fractionOfInjection(time0));
 
-    // correct mass if calculation is 2D 
+    // correct mass if calculation is 2D
     if (twoD)
     {
-        mInj *= 0.5*angleOfWedge/mathematicalConstant::pi;
+        mInj *= 0.5*angleOfWedge/constant::math::pi;
     }
 
     return mInj;
 }
 
+
 Foam::scalar Foam::multiHoleInjector::mass() const
 {
     return mass_;
 }
 
+
 const Foam::scalarField& Foam::multiHoleInjector::X() const
 {
     return X_;
 }
 
+
 Foam::List<Foam::multiHoleInjector::pair> Foam::multiHoleInjector::T() const
 {
     return TProfile_;
 }
 
+
 Foam::scalar Foam::multiHoleInjector::T(const scalar time) const
 {
     return getTableValue(TProfile_, time);
 }
 
+
 Foam::scalar Foam::multiHoleInjector::tsoi() const
 {
     return massFlowRateProfile_[0][0];
 }
 
+
 Foam::scalar Foam::multiHoleInjector::teoi() const
 {
     return massFlowRateProfile_[massFlowRateProfile_.size()-1][0];
 }
 
+
 Foam::scalar Foam::multiHoleInjector::massFlowRate
 (
     const scalar time
@@ -345,6 +370,7 @@ Foam::scalar Foam::multiHoleInjector::massFlowRate
     return getTableValue(massFlowRateProfile_, time);
 }
 
+
 Foam::scalar Foam::multiHoleInjector::injectionPressure
 (
     const scalar time
@@ -353,6 +379,7 @@ Foam::scalar Foam::multiHoleInjector::injectionPressure
     return getTableValue(injectionPressureProfile_, time);
 }
 
+
 Foam::scalar Foam::multiHoleInjector::velocity
 (
     const scalar time
@@ -361,11 +388,14 @@ Foam::scalar Foam::multiHoleInjector::velocity
     return getTableValue(velocityProfile_, time);
 }
 
-Foam::List<Foam::multiHoleInjector::pair> Foam::multiHoleInjector::CdProfile() const
+
+Foam::List<Foam::multiHoleInjector::pair> Foam::multiHoleInjector::CdProfile()
+const
 {
     return CdProfile_;
 }
 
+
 Foam::scalar Foam::multiHoleInjector::Cd
 (
     const scalar time
@@ -374,11 +404,16 @@ Foam::scalar Foam::multiHoleInjector::Cd
     return Cd_;
 }
 
-Foam::scalar Foam::multiHoleInjector::fractionOfInjection(const scalar time) const
+
+Foam::scalar Foam::multiHoleInjector::fractionOfInjection
+(
+    const scalar time
+) const
 {
     return integrateTable(massFlowRateProfile_, time)/mass_;
 }
 
+
 Foam::scalar Foam::multiHoleInjector::injectedMass
 (
     const scalar t
@@ -394,8 +429,7 @@ void Foam::multiHoleInjector::correctProfiles
     const scalar referencePressure
 )
 {
-
-    scalar A = nHoles_*0.25*mathematicalConstant::pi*pow(d_, 2.0);
+    scalar A = nHoles_*0.25*constant::math::pi*sqr(d_);
 
     forAll(velocityProfile_, i)
     {
@@ -407,14 +441,17 @@ void Foam::multiHoleInjector::correctProfiles
     }
 }
 
+
 Foam::vector Foam::multiHoleInjector::tan1(const label n) const
 {
     return tangentialInjectionVector1_[n];
 }
 
+
 Foam::vector Foam::multiHoleInjector::tan2(const label n) const
 {
     return tangentialInjectionVector2_[n];
 }
 
+
 // ************************************************************************* //

src/lagrangian/dieselSpray/injector/multiHoleInjector/multiHoleInjector.H

@@ -47,7 +47,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class multiHoleInjector Declaration
+                     Class multiHoleInjector Declaration
 \*---------------------------------------------------------------------------*/
 
 class multiHoleInjector
@@ -123,9 +123,8 @@ public:
         );
 
 
-    // Destructor
-
-        ~multiHoleInjector();
+    //- Destructor
+    virtual ~multiHoleInjector();
 
 
     // Member Functions
@@ -152,7 +151,7 @@ public:
             const vector& axisOfWedgeNormal,
             Random& rndGen
         ) const;
-    
+
         //- Return the number of holes
         label nHoles() const;
 
@@ -207,14 +206,14 @@ public:
         {
             return injectionPressureProfile_;
         }
-    
+
         scalar injectionPressure(const scalar time) const;
 
         List<pair> velocityProfile() const
         {
             return velocityProfile_;
         }
-        
+
         scalar velocity(const scalar time) const;
 
         List<pair> CdProfile() const;
@@ -233,7 +232,6 @@ public:
         {
             return pressureIndependentVelocity_;
         }
-
 };