From ee7a5233e77a6108a3babb2b939914e44b9d016b Mon Sep 17 00:00:00 2001
From: Pascal Beckstein <pascal@becksteins.de>
Date: Wed, 10 Feb 2016 18:15:52 +0100
Subject: [PATCH 01/11] Added ALL option for cellSource pointToCell of
 meshTools

---
 .../cellSources/pointToCell/pointToCell.C     | 62 +++++++++++++++++--
 .../cellSources/pointToCell/pointToCell.H     |  6 +-
 2 files changed, 59 insertions(+), 9 deletions(-)

diff --git a/src/meshTools/sets/cellSources/pointToCell/pointToCell.C b/src/meshTools/sets/cellSources/pointToCell/pointToCell.C
index 7c78a1db2..d522a722c 100644
--- a/src/meshTools/sets/cellSources/pointToCell/pointToCell.C
+++ b/src/meshTools/sets/cellSources/pointToCell/pointToCell.C
@@ -46,18 +46,21 @@ addToRunTimeSelectionTable(topoSetSource, pointToCell, istream);
 Foam::topoSetSource::addToUsageTable Foam::pointToCell::usage_
 (
     pointToCell::typeName,
-    "\n    Usage: pointToCell <pointSet> any\n\n"
-    "    Select all cells with any point in the pointSet\n\n"
+    "\n    Usage: pointToCell <pointSet> any|all\n\n"
+    "    Select cells with\n"
+    "    -any point in the pointSet\n"
+    "    -all points in the pointSet\n\n"
 );
 
 template<>
-const char* Foam::NamedEnum<Foam::pointToCell::pointAction, 1>::names[] =
+const char* Foam::NamedEnum<Foam::pointToCell::pointAction, 2>::names[] =
 {
-    "any"
+    "any",
+    "all"
 };
 
 
-const Foam::NamedEnum<Foam::pointToCell::pointAction, 1>
+const Foam::NamedEnum<Foam::pointToCell::pointAction, 2>
     Foam::pointToCell::pointActionNames_;
 
 
@@ -69,9 +72,9 @@ void Foam::pointToCell::combine(topoSet& set, const bool add) const
     pointSet loadedSet(mesh_, setName_);
 
 
-    // Handle any selection
     if (option_ == ANY)
     {
+        // Add cells with any point in loadedSet
         for
         (
             pointSet::const_iterator iter = loadedSet.begin();
@@ -89,6 +92,53 @@ void Foam::pointToCell::combine(topoSet& set, const bool add) const
             }
         }
     }
+    else if (option_ == ALL)
+    {
+        // Add all cells whose points are all in set.
+
+        Map<label> numPoints(loadedSet.size());
+
+        forAllConstIter(pointSet, loadedSet, iter)
+        {
+            label pointI = iter.key();
+
+            const labelList& pCells = mesh_.pointCells()[pointI];
+
+            forAll(pCells, pCellI)
+            {
+                label cellI = pCells[pCellI];
+
+                Map<label>::iterator fndCell = numPoints.find(cellI);
+
+                if (fndCell == numPoints.end())
+                {
+                    numPoints.insert(cellI, 1);
+                }
+                else
+                {
+                    fndCell()++;
+                }
+            }
+        }
+
+
+        // Include cells that are referenced as many times as there are points
+        // in cell -> all points of cell
+        for
+        (
+            Map<label>::const_iterator iter = numPoints.begin();
+            iter != numPoints.end();
+            ++iter
+        )
+        {
+            label cellI = iter.key();
+
+            if (iter() == mesh_.cellPoints()[cellI].size())
+            {
+                addOrDelete(set, cellI, add);
+            }
+        }
+    }
 }
 
 
diff --git a/src/meshTools/sets/cellSources/pointToCell/pointToCell.H b/src/meshTools/sets/cellSources/pointToCell/pointToCell.H
index 0f7a168e4..7150f8337 100644
--- a/src/meshTools/sets/cellSources/pointToCell/pointToCell.H
+++ b/src/meshTools/sets/cellSources/pointToCell/pointToCell.H
@@ -55,8 +55,8 @@ public:
         //- Enumeration defining the valid options
         enum pointAction
         {
-            ANY     // Cells using any point in set
-            //ALL   // Possible extension: cells whose all points are in set
+            ANY,    // Cells using any point in set
+            ALL     // Cells whose all points are in set
         };
 
 private:
@@ -64,7 +64,7 @@ private:
         //- Add usage string
         static addToUsageTable usage_;
 
-        static const NamedEnum<pointAction, 1> pointActionNames_;
+        static const NamedEnum<pointAction, 2> pointActionNames_;
 
         //- Name of set to use
         word setName_;

From 244a23ff0af0651989679f4a9f66244934033998 Mon Sep 17 00:00:00 2001
From: Pascal Beckstein <pascal@becksteins.de>
Date: Wed, 10 Feb 2016 18:08:07 +0100
Subject: [PATCH 02/11] Added tetPolyPatches based on directMapped

---
 src/tetFiniteElement/Make/files               |  3 +
 .../directMapped/directMappedTetPolyPatch.C   | 52 +++++++++++
 .../directMapped/directMappedTetPolyPatch.H   | 86 +++++++++++++++++++
 .../directMappedWallTetPolyPatch.C            | 52 +++++++++++
 .../directMappedWallTetPolyPatch.H            | 84 ++++++++++++++++++
 5 files changed, 277 insertions(+)
 create mode 100644 src/tetFiniteElement/tetPolyMesh/tetPolyPatches/derived/directMapped/directMappedTetPolyPatch.C
 create mode 100644 src/tetFiniteElement/tetPolyMesh/tetPolyPatches/derived/directMapped/directMappedTetPolyPatch.H
 create mode 100644 src/tetFiniteElement/tetPolyMesh/tetPolyPatches/derived/directMapped/directMappedWallTetPolyPatch.C
 create mode 100644 src/tetFiniteElement/tetPolyMesh/tetPolyPatches/derived/directMapped/directMappedWallTetPolyPatch.H

diff --git a/src/tetFiniteElement/Make/files b/src/tetFiniteElement/Make/files
index bcf6dd84a..dfda74e3f 100644
--- a/src/tetFiniteElement/Make/files
+++ b/src/tetFiniteElement/Make/files
@@ -17,6 +17,7 @@ cyclicGgiTetPolyPatch = $(constraintTetPolyPatches)/cyclicGgi
 mixingPlaneTetPolyPatch = $(constraintTetPolyPatches)/mixingPlane
 globalTetPolyPatch    = $(constraintTetPolyPatches)/global
 wallTetPolyPatch      = $(derivedTetPolyPatches)/wall
+directMappedTetPolyPatch = $(derivedTetPolyPatches)/directMapped
 
 MapTetFemFields       = $(tetPolyMesh)/MapTetFemFields
 
@@ -39,6 +40,8 @@ $(mixingPlaneTetPolyPatch)/mixingPlaneTetPolyPatch.C
 $(globalTetPolyPatch)/globalTetPolyPatch.C
 $(globalTetPolyPatch)/calcGlobalTetPolyPatchAddr.C
 $(wallTetPolyPatch)/wallTetPolyPatch.C
+$(directMappedTetPolyPatch)/directMappedTetPolyPatch.C
+$(directMappedTetPolyPatch)/directMappedWallTetPolyPatch.C
 $(tetPolyBoundaryMesh)/tetPolyBoundaryMesh.C
 $(tetPolyMesh)/tetPolyMeshLduAddressing.C
 $(tetPolyMesh)/tetPolyMesh.C
diff --git a/src/tetFiniteElement/tetPolyMesh/tetPolyPatches/derived/directMapped/directMappedTetPolyPatch.C b/src/tetFiniteElement/tetPolyMesh/tetPolyPatches/derived/directMapped/directMappedTetPolyPatch.C
new file mode 100644
index 000000000..d7372e945
--- /dev/null
+++ b/src/tetFiniteElement/tetPolyMesh/tetPolyPatches/derived/directMapped/directMappedTetPolyPatch.C
@@ -0,0 +1,52 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation, either version 3 of the License, or (at your
+    option) any later version.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "directMappedTetPolyPatch.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(directMappedTetPolyPatch, 0);
+
+// Add the patch constructor functions to the hash tables
+addToRunTimeSelectionTable
+(
+    faceTetPolyPatch,
+    directMappedTetPolyPatch,
+    polyPatch
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
+
diff --git a/src/tetFiniteElement/tetPolyMesh/tetPolyPatches/derived/directMapped/directMappedTetPolyPatch.H b/src/tetFiniteElement/tetPolyMesh/tetPolyPatches/derived/directMapped/directMappedTetPolyPatch.H
new file mode 100644
index 000000000..48995e57e
--- /dev/null
+++ b/src/tetFiniteElement/tetPolyMesh/tetPolyPatches/derived/directMapped/directMappedTetPolyPatch.H
@@ -0,0 +1,86 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation, either version 3 of the License, or (at your
+    option) any later version.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::directMappedTetPolyPatch
+
+Description
+    Direct mapped patch.
+
+SourceFiles
+    directMappedTetPolyPatch.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef directMappedTetPolyPatch_H
+#define directMappedTetPolyPatch_H
+
+#include "faceTetPolyPatch.H"
+#include "directMappedPolyPatch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                     Class directMappedTetPolyPatch Declaration
+\*---------------------------------------------------------------------------*/
+
+class directMappedTetPolyPatch
+:
+    public faceTetPolyPatch
+{
+
+
+
+public:
+
+    //- Runtime type information
+    TypeName(directMappedPolyPatch::typeName_());
+
+
+    // Constructors
+
+        //- Construct from polyPatch
+        directMappedTetPolyPatch
+        (
+            const polyPatch& patch,
+            const tetPolyBoundaryMesh& bm
+        )
+        :
+            faceTetPolyPatch(patch, bm)
+        {}
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
+
diff --git a/src/tetFiniteElement/tetPolyMesh/tetPolyPatches/derived/directMapped/directMappedWallTetPolyPatch.C b/src/tetFiniteElement/tetPolyMesh/tetPolyPatches/derived/directMapped/directMappedWallTetPolyPatch.C
new file mode 100644
index 000000000..609a17458
--- /dev/null
+++ b/src/tetFiniteElement/tetPolyMesh/tetPolyPatches/derived/directMapped/directMappedWallTetPolyPatch.C
@@ -0,0 +1,52 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation, either version 3 of the License, or (at your
+    option) any later version.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "directMappedWallTetPolyPatch.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(directMappedWallTetPolyPatch, 0);
+
+// Add the patch constructor functions to the hash tables
+addToRunTimeSelectionTable
+(
+    faceTetPolyPatch,
+    directMappedWallTetPolyPatch,
+    polyPatch
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
+
diff --git a/src/tetFiniteElement/tetPolyMesh/tetPolyPatches/derived/directMapped/directMappedWallTetPolyPatch.H b/src/tetFiniteElement/tetPolyMesh/tetPolyPatches/derived/directMapped/directMappedWallTetPolyPatch.H
new file mode 100644
index 000000000..f6db4a97d
--- /dev/null
+++ b/src/tetFiniteElement/tetPolyMesh/tetPolyPatches/derived/directMapped/directMappedWallTetPolyPatch.H
@@ -0,0 +1,84 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation, either version 3 of the License, or (at your
+    option) any later version.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::directMappedWallTetPolyPatch
+
+Description
+    Direct mapped wall patch.
+
+SourceFiles
+    directMappedWallTetPolyPatch.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef directMappedWallTetPolyPatch_H
+#define directMappedWallTetPolyPatch_H
+
+#include "wallTetPolyPatch.H"
+#include "directMappedWallPolyPatch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                   Class directMappedWallTetPolyPatch Declaration
+\*---------------------------------------------------------------------------*/
+
+class directMappedWallTetPolyPatch
+:
+    public wallTetPolyPatch
+{
+
+public:
+
+    //- Runtime type information
+    TypeName(directMappedWallPolyPatch::typeName_());
+
+
+    // Constructors
+
+        //- Construct from polyPatch
+        directMappedWallTetPolyPatch
+        (
+            const polyPatch& patch,
+            const tetPolyBoundaryMesh& bm
+        )
+        :
+            wallTetPolyPatch(patch, bm)
+        {}
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
+

From 51047b302d9b65e14c4ce6f825f70b09f2c5ec37 Mon Sep 17 00:00:00 2001
From: Pascal Beckstein <pascal@becksteins.de>
Date: Thu, 11 Feb 2016 11:54:10 +0100
Subject: [PATCH 03/11] Make faMesh decomposition work for regions other than
 defaultRegion

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

diff --git a/applications/utilities/parallelProcessing/decomposePar/faMeshDecomposition.C b/applications/utilities/parallelProcessing/decomposePar/faMeshDecomposition.C
index ec06c6705..effd4cd06 100644
--- a/applications/utilities/parallelProcessing/decomposePar/faMeshDecomposition.C
+++ b/applications/utilities/parallelProcessing/decomposePar/faMeshDecomposition.C
@@ -55,7 +55,7 @@ void faMeshDecomposition::distributeFaces()
         (
             IOobject
             (
-                fvMesh::defaultRegion,
+                GeoMesh<polyMesh>::mesh_.name(),
                 processorDb.timeName(),
                 processorDb
             )
@@ -252,7 +252,7 @@ void faMeshDecomposition::decomposeMesh(const bool filterEmptyPatches)
         (
             IOobject
             (
-                fvMesh::defaultRegion,
+                GeoMesh<polyMesh>::mesh_.name(),
                 processorDb.timeName(),
                 processorDb
             )
@@ -1073,7 +1073,7 @@ void faMeshDecomposition::decomposeMesh(const bool filterEmptyPatches)
         (
             IOobject
             (
-                fvMesh::defaultRegion,
+                GeoMesh<polyMesh>::mesh_.name(),
                 processorDb.timeName(),
                 processorDb
             )
@@ -1191,7 +1191,7 @@ bool faMeshDecomposition::writeDecomposition()
         (
             IOobject
             (
-                fvMesh::defaultRegion,
+                GeoMesh<polyMesh>::mesh_.name(),
                 processorDb.timeName(),
                 processorDb
             )

From b6a5ca25502a094358957a9d1678dc1f07ddab69 Mon Sep 17 00:00:00 2001
From: Pascal Beckstein <pascal@becksteins.de>
Date: Thu, 11 Feb 2016 15:44:59 +0100
Subject: [PATCH 04/11] Add region option for all finiteArea utilities

---
 .../utilities/finiteArea/checkFaMesh/checkFaMesh.C       | 4 +++-
 .../finiteArea/createFaMeshFromStl/createFaMeshFromStl.C | 2 ++
 .../finiteArea/createFaMeshFromStl/makePolyMesh.H        | 9 ++++++++-
 .../utilities/finiteArea/makeFaMesh/makeFaMesh.C         | 4 +++-
 4 files changed, 16 insertions(+), 3 deletions(-)

diff --git a/applications/utilities/finiteArea/checkFaMesh/checkFaMesh.C b/applications/utilities/finiteArea/checkFaMesh/checkFaMesh.C
index 5a155a806..66458407c 100644
--- a/applications/utilities/finiteArea/checkFaMesh/checkFaMesh.C
+++ b/applications/utilities/finiteArea/checkFaMesh/checkFaMesh.C
@@ -37,9 +37,11 @@ using namespace Foam;
 
 int main(int argc, char *argv[])
 {
+#   include "addRegionOption.H"
+
 #   include "setRootCase.H"
 #   include "createTime.H"
-#   include "createMesh.H"
+#   include "createNamedMesh.H"
 #   include "createFaMesh.H"
 
     Info<< "Time = " << runTime.timeName() << nl << endl;
diff --git a/applications/utilities/finiteArea/createFaMeshFromStl/createFaMeshFromStl.C b/applications/utilities/finiteArea/createFaMeshFromStl/createFaMeshFromStl.C
index 1ea27f407..e5253a028 100644
--- a/applications/utilities/finiteArea/createFaMeshFromStl/createFaMeshFromStl.C
+++ b/applications/utilities/finiteArea/createFaMeshFromStl/createFaMeshFromStl.C
@@ -50,6 +50,8 @@ using namespace Foam;
 
 int main(int argc, char *argv[])
 {
+#   include "addRegionOption.H"
+
     argList::validArgs.append("STL mesh file");
 
     argList args(argc, argv);
diff --git a/applications/utilities/finiteArea/createFaMeshFromStl/makePolyMesh.H b/applications/utilities/finiteArea/createFaMeshFromStl/makePolyMesh.H
index 13872316c..decbd99ac 100644
--- a/applications/utilities/finiteArea/createFaMeshFromStl/makePolyMesh.H
+++ b/applications/utilities/finiteArea/createFaMeshFromStl/makePolyMesh.H
@@ -23,11 +23,18 @@
         faceRegion[faceI] = patch[faceI].region();
     }
 
+    word polyRegionName;
+
+    if (!args.optionReadIfPresent("region", polyRegionName))
+    {
+        polyRegionName = polyMesh::defaultRegion;
+    }
+
     polyMesh mesh
     (
         IOobject
         (
-            polyMesh::defaultRegion,
+            polyRegionName,
             runTime.constant(),
             runTime
         ),
diff --git a/applications/utilities/finiteArea/makeFaMesh/makeFaMesh.C b/applications/utilities/finiteArea/makeFaMesh/makeFaMesh.C
index 725b1defb..372e122e4 100644
--- a/applications/utilities/finiteArea/makeFaMesh/makeFaMesh.C
+++ b/applications/utilities/finiteArea/makeFaMesh/makeFaMesh.C
@@ -63,9 +63,11 @@ public:
 
 int main(int argc, char *argv[])
 {
+#   include "addRegionOption.H"
+
 #   include "setRootCase.H"
 #   include "createTime.H"
-#   include "createMesh.H"
+#   include "createNamedMesh.H"
 
     // Reading faMeshDefinition dictionary
     IOdictionary faMeshDefinition

From 50574d8b1c8ebfcac3c4f17746d9a0916809fbf9 Mon Sep 17 00:00:00 2001
From: Pascal Beckstein <pascal@becksteins.de>
Date: Thu, 11 Feb 2016 17:08:17 +0100
Subject: [PATCH 05/11] Correctly ignore changes to the configuration and
 include files from KDevelop

---
 .gitignore | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/.gitignore b/.gitignore
index e4914e5a7..79eb72b26 100644
--- a/.gitignore
+++ b/.gitignore
@@ -13,6 +13,9 @@
 # file-browser settings - anywhere
 .directory
 
+# KDevelop include-path files - anywhere
+.kdev_include_paths
+
 # CVS recovered versions - anywhere
 .#*
 
@@ -104,8 +107,8 @@ etc/cshrc
 etc/bashrc
 etc/settings.csh
 etc/settings.sh
-etc/pref.csh
-etc/pref.sh
+etc/prefs.csh
+etc/prefs.sh
 
 # make sure that this settings file is not used
 etc/bashrc.preset

From 512b75bef50fe80deebb37b78767b41f52fee38f Mon Sep 17 00:00:00 2001
From: Pascal Beckstein <pascal@becksteins.de>
Date: Thu, 7 Apr 2016 13:54:24 +0200
Subject: [PATCH 06/11] Fix for assignment of value for
 UniformDimensionedScalarField if constructed from Istream via dictionary
 lookup

---
 .../fields/UniformDimensionedFields/UniformDimensionedField.C   | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/foam/fields/UniformDimensionedFields/UniformDimensionedField.C b/src/foam/fields/UniformDimensionedFields/UniformDimensionedField.C
index a2c19e36c..5be501332 100644
--- a/src/foam/fields/UniformDimensionedFields/UniformDimensionedField.C
+++ b/src/foam/fields/UniformDimensionedFields/UniformDimensionedField.C
@@ -61,7 +61,7 @@ Foam::UniformDimensionedField<Type>::UniformDimensionedField
 {
     dictionary dict(readStream(typeName));
     this->dimensions().reset(dict.lookup("dimensions"));
-    this->value() = dict.lookup("value");
+    this->value() = pTraits<Type>(dict.lookup("value"));
 }
 
 

From 80aa1d588f65fa4a1e0e2049d9a439e91a69112b Mon Sep 17 00:00:00 2001
From: Pascal Beckstein <pascal@becksteins.de>
Date: Wed, 10 Feb 2016 18:38:26 +0100
Subject: [PATCH 07/11] Fix: Prevent erroneous warning about failed library
 loading if new dynamicMesh is constructed multiple times

---
 .../dynamicFvMesh/newDynamicFvMesh.C          | 24 +++++++++++++++----
 1 file changed, 20 insertions(+), 4 deletions(-)

diff --git a/src/dynamicMesh/dynamicFvMesh/dynamicFvMesh/newDynamicFvMesh.C b/src/dynamicMesh/dynamicFvMesh/dynamicFvMesh/newDynamicFvMesh.C
index f163d165e..629690e58 100644
--- a/src/dynamicMesh/dynamicFvMesh/dynamicFvMesh/newDynamicFvMesh.C
+++ b/src/dynamicMesh/dynamicFvMesh/dynamicFvMesh/newDynamicFvMesh.C
@@ -37,13 +37,29 @@ Foam::autoPtr<Foam::dynamicFvMesh> Foam::dynamicFvMesh::New(const IOobject& io)
     forAll(libNames,i) {
         const word libName("lib"+libNames[i]+".so");
 
-        bool ok=dlLibraryTable::open(libName);
-        if(!ok) {
-            WarningIn("dynamicFvMesh::New(const IOobject& io)")
-                << "Loading of dynamic mesh library " << libName
+        bool libLoaded = false;
+
+        dlLibraryTable& ll = dlLibraryTable::loadedLibraries;
+
+        forAllConstIter(dlLibraryTable, ll, llI)
+        {
+            if (ll(llI.key()) == libName)
+            {
+                libLoaded = true;
+                break;
+            }
+        }
+
+        if (!libLoaded)
+        {
+            bool ok=dlLibraryTable::open(libName);
+            if(!ok) {
+                WarningIn("dynamicFvMesh::New(const IOobject& io)")
+                    << "Loading of dynamic mesh library " << libName
                     << " unsuccesful. Some dynamic mesh  methods may not be "
                     << " available"
                     << endl;
+            }
         }
     }
 

From a90d26789073c4687abbbe3afad6920167237cb2 Mon Sep 17 00:00:00 2001
From: Pascal Beckstein <pascal@becksteins.de>
Date: Sat, 23 Apr 2016 21:46:38 +0200
Subject: [PATCH 08/11] Fix public access to curMotionTimeIndex of polyMesh

---
 src/foam/meshes/polyMesh/polyMesh.H | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/foam/meshes/polyMesh/polyMesh.H b/src/foam/meshes/polyMesh/polyMesh.H
index ba7c7a396..01fe6e65e 100644
--- a/src/foam/meshes/polyMesh/polyMesh.H
+++ b/src/foam/meshes/polyMesh/polyMesh.H
@@ -425,7 +425,7 @@ public:
             //- Return current motion time index
             label curMotionTimeIndex() const
             {
-                return moving_;
+                return curMotionTimeIndex_;
             }
 
             //- Set the mesh to be moving

From 1344797239c6491f689c0d028c0c8c6eb3f116ea Mon Sep 17 00:00:00 2001
From: Pascal Beckstein <pascal@becksteins.de>
Date: Wed, 10 Feb 2016 18:50:02 +0100
Subject: [PATCH 09/11] Fix: Correct REMOVE option for setSet utility

---
 .../mesh/manipulation/setSet/setSet.C          | 18 +++++++++++-------
 1 file changed, 11 insertions(+), 7 deletions(-)

diff --git a/applications/utilities/mesh/manipulation/setSet/setSet.C b/applications/utilities/mesh/manipulation/setSet/setSet.C
index 86ace70ca..abdfe6ba1 100644
--- a/applications/utilities/mesh/manipulation/setSet/setSet.C
+++ b/applications/utilities/mesh/manipulation/setSet/setSet.C
@@ -325,7 +325,7 @@ void printAllSets(const polyMesh& mesh, Ostream& os)
 
 
 // Physically remove a set
-void removeSet
+bool removeSet
 (
     const polyMesh& mesh,
     const word& setType,
@@ -346,7 +346,11 @@ void removeSet
         fileName object = objects[setName]->objectPath();
         Info<< "Removing file " << object << endl;
         rm(object);
+
+        return true;
     }
+
+    return false;
 }
 
 
@@ -392,11 +396,7 @@ bool doCommand
 
         IOobject::readOption r;
 
-        if (action == topoSetSource::REMOVE)
-        {
-            removeSet(mesh, setType, setName);
-        }
-        else if
+        if
         (
             (action == topoSetSource::NEW)
          || (action == topoSetSource::CLEAR)
@@ -414,7 +414,11 @@ bool doCommand
             currentSet.resize(max(currentSet.size(), typSize));
         }
 
-        if (!currentSetPtr.valid())
+        if (action == topoSetSource::REMOVE)
+        {
+            ok = removeSet(mesh, setType, setName);
+        }
+        else if (!currentSetPtr.valid())
         {
             Info<< "    Cannot construct/load set "
                 << topoSet::localPath(mesh, setName) << endl;

From e655cbc255752da9a3b5c57b3588722dd6bf6cb3 Mon Sep 17 00:00:00 2001
From: Pascal Beckstein <pascal@becksteins.de>
Date: Wed, 10 Feb 2016 18:07:37 +0100
Subject: [PATCH 10/11] Added new faceSource rotatedBoxToFace for meshTools

---
 src/meshTools/Make/files                      |   1 +
 .../rotatedBoxToFace/rotatedBoxToFace.C       | 193 ++++++++++++++++++
 .../rotatedBoxToFace/rotatedBoxToFace.H       | 133 ++++++++++++
 3 files changed, 327 insertions(+)
 create mode 100644 src/meshTools/sets/faceSources/rotatedBoxToFace/rotatedBoxToFace.C
 create mode 100644 src/meshTools/sets/faceSources/rotatedBoxToFace/rotatedBoxToFace.H

diff --git a/src/meshTools/Make/files b/src/meshTools/Make/files
index fcd0d492c..8675cad0d 100644
--- a/src/meshTools/Make/files
+++ b/src/meshTools/Make/files
@@ -73,6 +73,7 @@ $(faceSources)/boundaryToFace/boundaryToFace.C
 $(faceSources)/zoneToFace/zoneToFace.C
 $(faceSources)/setToFace/setToFace.C
 $(faceSources)/boxToFace/boxToFace.C
+$(faceSources)/rotatedBoxToFace/rotatedBoxToFace.C
 
 pointSources = sets/pointSources
 $(pointSources)/labelToPoint/labelToPoint.C
diff --git a/src/meshTools/sets/faceSources/rotatedBoxToFace/rotatedBoxToFace.C b/src/meshTools/sets/faceSources/rotatedBoxToFace/rotatedBoxToFace.C
new file mode 100644
index 000000000..2544eed89
--- /dev/null
+++ b/src/meshTools/sets/faceSources/rotatedBoxToFace/rotatedBoxToFace.C
@@ -0,0 +1,193 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation, either version 3 of the License, or (at your
+    option) any later version.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "rotatedBoxToFace.H"
+#include "polyMesh.H"
+#include "cellModeller.H"
+
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(rotatedBoxToFace, 0);
+
+addToRunTimeSelectionTable(topoSetSource, rotatedBoxToFace, word);
+
+addToRunTimeSelectionTable(topoSetSource, rotatedBoxToFace, istream);
+
+}
+
+
+Foam::topoSetSource::addToUsageTable Foam::rotatedBoxToFace::usage_
+(
+    rotatedBoxToFace::typeName,
+    "\n    Usage: rotatedBoxToFace (originx originy originz)"
+    " (ix iy iz) (jx jy jz) (kx ky kz)\n\n"
+    "    Select all faces with faceCentre within parallelopiped\n\n"
+);
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::rotatedBoxToFace::combine(topoSet& set, const bool add) const
+{
+    // Define a cell for the box
+    pointField boxPoints(8);
+    boxPoints[0] = origin_;
+    boxPoints[1] = origin_ + i_;
+    boxPoints[2] = origin_ + i_ + j_;
+    boxPoints[3] = origin_ + j_;
+    boxPoints[4] = origin_ + k_;
+    boxPoints[5] = origin_ + k_ + i_;
+    boxPoints[6] = origin_ + k_ + i_ + j_;
+    boxPoints[7] = origin_ + k_ + j_;
+
+    labelList boxVerts(8);
+    forAll(boxVerts, i)
+    {
+        boxVerts[i] = i;
+    }
+
+    const cellModel& hex = *(cellModeller::lookup("hex"));
+
+    // Get outwards pointing faces.
+    faceList boxFaces(cellShape(hex, boxVerts).faces());
+
+    // Precalculate normals
+    vectorField boxFaceNormals(boxFaces.size());
+    forAll(boxFaces, i)
+    {
+        boxFaceNormals[i] = boxFaces[i].normal(boxPoints);
+
+        Pout<< "Face:" << i << " position:" << boxFaces[i].centre(boxPoints)
+            << " normal:" << boxFaceNormals[i] << endl;
+    }
+
+    // Check whether face centre is inside all faces of box.
+
+    const pointField& ctrs = mesh_.faceCentres();
+
+    forAll(ctrs, faceI)
+    {
+        bool inside = true;
+
+        forAll(boxFaces, i)
+        {
+            const face& f = boxFaces[i];
+
+            if (((ctrs[faceI] - boxPoints[f[0]]) & boxFaceNormals[i]) > 0)
+            {
+                inside = false;
+                break;
+            }
+        }
+
+        if (inside)
+        {
+            addOrDelete(set, faceI, add);
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+Foam::rotatedBoxToFace::rotatedBoxToFace
+(
+    const polyMesh& mesh,
+    const vector& origin,
+    const vector& i,
+    const vector& j,
+    const vector& k
+)
+:
+    topoSetSource(mesh),
+    origin_(origin),
+    i_(i),
+    j_(j),
+    k_(k)
+{}
+
+
+// Construct from dictionary
+Foam::rotatedBoxToFace::rotatedBoxToFace
+(
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+:
+    topoSetSource(mesh),
+    origin_(dict.lookup("origin")),
+    i_(dict.lookup("i")),
+    j_(dict.lookup("j")),
+    k_(dict.lookup("k"))
+{}
+
+
+// Construct from Istream
+Foam::rotatedBoxToFace::rotatedBoxToFace(const polyMesh& mesh, Istream& is)
+:
+    topoSetSource(mesh),
+    origin_(is),
+    i_(is),
+    j_(is),
+    k_(is)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::rotatedBoxToFace::~rotatedBoxToFace()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::rotatedBoxToFace::applyToSet
+(
+    const topoSetSource::setAction action,
+    topoSet& set
+) const
+{
+    if ((action == topoSetSource::NEW) || (action == topoSetSource::ADD))
+    {
+        Info<< "    Adding faces with center within rotated box " << endl;
+
+        combine(set, true);
+    }
+    else if (action == topoSetSource::DELETE)
+    {
+        Info<< "    Removing faces with center within rotated box " << endl;
+
+        combine(set, false);
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/meshTools/sets/faceSources/rotatedBoxToFace/rotatedBoxToFace.H b/src/meshTools/sets/faceSources/rotatedBoxToFace/rotatedBoxToFace.H
new file mode 100644
index 000000000..ac66b696b
--- /dev/null
+++ b/src/meshTools/sets/faceSources/rotatedBoxToFace/rotatedBoxToFace.H
@@ -0,0 +1,133 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation, either version 3 of the License, or (at your
+    option) any later version.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::rotatedBoxToFace
+
+Description
+    A topoSetSource to select faces based on face centres inside
+    rotated/skewed box (parallelopiped?).
+
+    Box defined as origin and i,j,k vectors.
+    E.g. box rotated 45 degrees around z-axis with sizes sqrt(0.2^2+0.2^2)
+    (and extra large, 200 in z direction):
+    @verbatim
+       origin   ( 0.4 0.4 -100);
+       i        ( 0.2 0.2    0);
+       j        (-0.2 0.2    0);
+       k        ( 0.0 0.0  100);
+    @endverbatim
+
+SourceFiles
+    rotatedBoxToFace.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef rotatedBoxToFace_H
+#define rotatedBoxToFace_H
+
+#include "topoSetSource.H"
+#include "treeBoundBox.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                     Class rotatedBoxToFace Declaration
+\*---------------------------------------------------------------------------*/
+
+class rotatedBoxToFace
+:
+    public topoSetSource
+{
+
+    // Private data
+
+        //- Add usage string
+        static addToUsageTable usage_;
+
+
+        //- skewed box
+        const vector origin_;
+        const vector i_;
+        const vector j_;
+        const vector k_;
+
+
+    // Private Member Functions
+
+        void combine(topoSet& set, const bool add) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("rotatedBoxToFace");
+
+    // Constructors
+
+        //- Construct from components
+        rotatedBoxToFace
+        (
+            const polyMesh& mesh,
+            const vector& origin,
+            const vector& i,
+            const vector& j,
+            const vector& k
+        );
+
+        //- Construct from dictionary
+        rotatedBoxToFace(const polyMesh& mesh, const dictionary& dict);
+
+        //- Construct from Istream
+        rotatedBoxToFace(const polyMesh& mesh, Istream&);
+
+
+    // Destructor
+
+        virtual ~rotatedBoxToFace();
+
+
+    // Member Functions
+
+        virtual void applyToSet
+        (
+            const topoSetSource::setAction action,
+            topoSet&
+        ) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
+

From 244e6aa02817e5141bb0dba9f91c8eadee6b17aa Mon Sep 17 00:00:00 2001
From: Henrik Rusche <h.rusche@wikki-gmbh.de>
Date: Wed, 6 Jul 2016 23:32:35 +0200
Subject: [PATCH 11/11] BUGFIX: Switching real gas tutorial back to Peng
 Robinson / Some clean-up

---
 .../aungierRedlichKwong/aungierRedlichKwong.C | 10 +--
 .../aungierRedlichKwong/aungierRedlichKwong.H |  6 +-
 .../aungierRedlichKwongI.H                    | 23 ++---
 .../pengRobinson/pengRobinson.C               | 10 +--
 .../pengRobinson/pengRobinson.H               | 10 +--
 .../pengRobinson/pengRobinsonI.H              | 84 ++++++++-----------
 .../redlichKwong/redlichKwong.C               |  2 -
 .../redlichKwong/redlichKwong.H               |  2 -
 .../redlichKwong/redlichKwongI.H              | 61 ++++++--------
 .../soaveRedlichKwong/soaveRedlichKwong.C     | 14 ++--
 .../soaveRedlichKwong/soaveRedlichKwong.H     |  8 +-
 .../soaveRedlichKwong/soaveRedlichKwongI.H    | 59 ++++++-------
 .../constant/thermophysicalProperties         |  8 +-
 .../ras/t-junction/system/controlDict         |  2 +-
 14 files changed, 125 insertions(+), 174 deletions(-)

diff --git a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwong.C b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwong.C
index d1164ec30..983244872 100755
--- a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwong.C
+++ b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwong.C
@@ -56,12 +56,12 @@ Foam::aungierRedlichKwong::aungierRedlichKwong(Istream& is)
     //CL: Only uses the default values
     rhoMin_(1e-3),
     rhoMax_(1500),
-    // Starting GUESS for the density by ideal gas law
-    rhostd_(this->rho(this->Pstd(), this->Tstd(), this->Pstd()/(this->Tstd()*this->R()))),
     aSave(0.0),
     daSave(0.0),
     d2aSave(0.0),
-    TSave(0.0)
+    TSave(0.0),
+    // Starting GUESS for the density by ideal gas law
+    rhostd_(this->rho(this->Pstd(), this->Tstd(), this->Pstd()/(this->Tstd()*this->R())))
 {
     is.check("aungierRedlichKwong::aungierRedlichKwong(Istream& is)");
 }
@@ -91,12 +91,12 @@ Foam::aungierRedlichKwong::aungierRedlichKwong(const dictionary& dict)
     //CL: therefore, rho can be larger than rhoMax and smaller than rhoMin
     rhoMin_(dict.subDict("equationOfState").lookupOrDefault("rhoMin",1e-3)),
     rhoMax_(dict.subDict("equationOfState").lookupOrDefault("rhoMax",1500)),
-    // Starting GUESS for the density by ideal gas law
-    rhostd_(this->rho(this->Pstd(), this->Tstd(), this->Pstd()/(this->Tstd()*this->R()))),
     aSave(0.0),
     daSave(0.0),
     d2aSave(0.0),
     TSave(0.0)
+    // Starting GUESS for the density by ideal gas law
+    rhostd_(this->rho(this->Pstd(), this->Tstd(), this->Pstd()/(this->Tstd()*this->R()))),
 {
     is.check("aungierRedlichKwong::aungierRedlichKwong(Istream& is)");
 }
diff --git a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwong.H b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwong.H
index 472ccf69c..264dae224 100755
--- a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwong.H
+++ b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwong.H
@@ -87,9 +87,6 @@ class aungierRedlichKwong
     scalar rhoMin_;
     scalar rhoMax_;
 
-    //- Density @STD, initialise after a0, b, c, n!
-    scalar rhostd_;
-
     //CL: Variables to save the values of a, dadT and d2adT2 of the mixture
     //CL: Variables must corrected for changing temperatures
     mutable scalar aSave;
@@ -99,6 +96,9 @@ class aungierRedlichKwong
     //CL: save the temperature for which the save coefficients (aSave, daSave, d2aSave) are correct
     mutable scalar TSave;
 
+    //- Density @STD, initialise after a0, b, c, n!
+    scalar rhostd_;
+
     //Protected functions
     //CL: function updates the coefficients (aSave, daSave, d2aSave)
     inline void updateModelCoefficients(const  scalar T) const;
diff --git a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwongI.H b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwongI.H
index 887717b8b..f4e37cfa6 100755
--- a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwongI.H
+++ b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/aungierRedlichKwong/aungierRedlichKwongI.H
@@ -68,11 +68,11 @@ inline aungierRedlichKwong::aungierRedlichKwong(const word& name, const aungierR
     c2_(pg.c2_),
     rhoMin_(pg.rhoMin_),
     rhoMax_(pg.rhoMax_),
-    rhostd_(pg.rhostd_),
     aSave(0.0),
     daSave(0.0),
     d2aSave(0.0),
-    TSave(0.0)
+    TSave(0.0),
+    rhostd_(pg.rhostd_)
 {}
 
 
@@ -237,13 +237,16 @@ inline scalar aungierRedlichKwong::dpdv(const scalar rho, const scalar T) const
 {
     scalar Vm = this->W()/rho;
     scalar Vm2 = Vm*Vm;
+    scalar bpVm = b_ + Vm;
+    scalar bpVm2 = bpVm*bpVm;
+    scalar cpVm = c_ + Vm;
 
     return
     (
-        a(T)*(b3_ - 2*b2_*c_ + b_*(c_ + Vm)*(c_ - 3*Vm) + 2*Vm*pow(c_ + Vm, 2))
-      - this->RR()*T*Vm2*(b2_ + 2*b_*Vm + Vm2)
+        a(T)*(b3_ - 2*b2_*c_ + b_*cpVm*(c_ - 3*Vm) + 2*Vm*cpVm*cpVm)
+      - this->RR()*T*Vm2*bpVm2
     )
-   /(Vm2*pow(b_ - c_ - Vm, 2)*pow(b_ + Vm, 2));
+   /(Vm2*pow(b_ - c_ - Vm, 2)*bpVm2);
 }
 
 
@@ -284,8 +287,7 @@ inline scalar aungierRedlichKwong::integral_p_dv
 {
     scalar Vm = this->W()/rho;
 
-    return this->RR()*T*log(-b_ + c_ + Vm) + a(T)/b_*(log(b_ + Vm) - log(Vm));
-    //return this->RR()*T*log(-b_ + c_ + Vm) + a(T)*log(b_ + Vm)/b_ - a(T)*log(Vm)/b_;
+    return this->RR()*T*log(Vm - b_ + c_) + a(T)/b_*log((b_ + Vm)/Vm);
 }
 
 
@@ -299,8 +301,7 @@ inline scalar aungierRedlichKwong::integral_dpdT_dv
 {
     scalar Vm = this->W()/rho;
 
-    return this->RR()*log(-b_ + c_ + Vm) + dadT(T)/b_*(log(b_ + Vm) - log(Vm));
-    //return this->RR()*log(-b_ + c_ + Vm) + dadT(T)*log(b_ + Vm)/b_ - dadT(T)*log(Vm)/b_;
+    return this->RR()*log(Vm - b_ + c_) + dadT(T)/b_*log((b_ + Vm)/Vm);
 }
 
 
@@ -309,7 +310,7 @@ inline scalar aungierRedlichKwong::d2pdT2(const scalar rho,const scalar T) const
 {
     scalar Vm = this->W()/rho;
 
-    return -d2adT2(T)/(Vm*(Vm + b_));
+    return -d2adT2(T)/(Vm*(b_ + Vm));
 }
 
 
@@ -400,7 +401,7 @@ inline scalar aungierRedlichKwong::integral_d2pdT2_dv
 {
     scalar Vm = this->W()/rho;
 
-    return d2adT2(T)*log(b_ + Vm)/b_ - d2adT2(T)*log(Vm)/b_;
+    return d2adT2(T)/b_*log((b_ + Vm)/Vm);
 }
 
 
diff --git a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinson.C b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinson.C
index 1beb1cf01..719920f50 100755
--- a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinson.C
+++ b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinson.C
@@ -47,19 +47,15 @@ Foam::pengRobinson::pengRobinson(Istream& is)
     b_(0.077796*this->RR()*Tcrit_/pcrit_),
     n_(0.37464 + 1.54226*azentricFactor_ - 0.26992*pow(azentricFactor_, 2)),
     b2_(b_*b_),
-    b3_(b2_*b_),
-    b4_(b3_*b_),
-    b5_(b4_*b_),
-    b6_(b5_*b_),
     //CL: Only uses the default values
     rhoMin_(1e-3),
     rhoMax_(1500),
-    // Starting GUESS for the density by ideal gas law
-    rhostd_(this->rho(this->Pstd(), this->Tstd(), this->Pstd()/(this->Tstd()*this->R()))),
     aSave(0.0),
     daSave(0.0),
     d2aSave(0.0),
-    TSave(0.0)
+    TSave(0.0),
+    // Starting GUESS for the density by ideal gas law
+    rhostd_(this->rho(this->Pstd(), this->Tstd(), this->Pstd()/(this->Tstd()*this->R())))
 {
     is.check("pengRobinson::pengRobinson(Istream& is)");
 }
diff --git a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinson.H b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinson.H
index 0829ae642..49c75b0d7 100755
--- a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinson.H
+++ b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinson.H
@@ -79,18 +79,11 @@ class pengRobinson
 
     //CL: pow of constants b_ used in the code e.g. b2_=b*b;
     scalar b2_;
-    scalar b3_;
-    scalar b4_;
-    scalar b5_;
-    scalar b6_;
 
     //CL: rhoMin and rhoMax are only used as boundaries for the bisection methode (see rho function)
     scalar rhoMin_;
     scalar rhoMax_;
 
-    //- Density @STD, initialise after a0, b!
-    scalar rhostd_;
-
     //CL: Variables to save the values of a, dadT and d2adT2 of the mixture
     //CL: Variables must corrected for changing temperatures
     mutable scalar aSave;
@@ -100,6 +93,9 @@ class pengRobinson
     //CL: save the temperature for which the save coefficients (aSave, daSave, d2aSave) are correct
     mutable scalar TSave;
 
+    //- Density @STD, initialise after a0, b!
+    scalar rhostd_;
+
     //Protected functions
     //CL: function updates the coefficients (aSave, daSave, d2aSave)
     inline void updateModelCoefficients(const  scalar T) const;
diff --git a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinsonI.H b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinsonI.H
index 39d74d245..26147d68e 100755
--- a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinsonI.H
+++ b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/pengRobinson/pengRobinsonI.H
@@ -60,17 +60,13 @@ inline pengRobinson::pengRobinson(const word& name, const pengRobinson& pr)
     b_(pr.b_),
     n_(pr.n_),
     b2_(pr.b2_),
-    b3_(pr.b3_),
-    b4_(pr.b4_),
-    b5_(pr.b5_),
-    b6_(pr.b6_),
     rhoMin_(pr.rhoMin_),
     rhoMax_(pr.rhoMax_),
-    rhostd_(pr.rhostd_),
     aSave(0.0),
     daSave(0.0),
     d2aSave(0.0),
-    TSave(0.0)
+    TSave(0.0),
+    rhostd_(pr.rhostd_)
 {}
 
 
@@ -123,7 +119,7 @@ inline void pengRobinson::updateModelCoefficients(const scalar T) const
     d2aSave = a0_*n_*(n_ + 1)*TTc/(2*T*T);
 
     //CL: saving the temperature at which the coefficients are valid
-    TSave=T;
+    TSave = T;
 }
 
 
@@ -131,7 +127,7 @@ inline void pengRobinson::updateModelCoefficients(const scalar T) const
 inline scalar pengRobinson::a(const scalar T) const
 {
     //CL: check if a has already been calculated for this temperature
-    if(TSave==T)
+    if(TSave == T)
     {
         return aSave;
     }
@@ -148,7 +144,7 @@ inline scalar pengRobinson::a(const scalar T) const
 inline scalar pengRobinson::dadT(const scalar T) const
 {
     // check if a has already been calculated for this temperature
-    if(TSave==T)
+    if(TSave == T)
     {
         return daSave;
     }
@@ -165,7 +161,7 @@ inline scalar pengRobinson::dadT(const scalar T) const
 inline scalar pengRobinson::d2adT2(const scalar T) const
 {
     // check if a has already been calculated for this temperature
-    if(TSave==T)
+    if(TSave == T)
     {
         return d2aSave;
     }
@@ -210,22 +206,15 @@ inline scalar pengRobinson::p(const scalar rho, const scalar T) const
 inline scalar pengRobinson::dpdv(const scalar rho, const scalar T) const
 {
     scalar Vm = this->W()/rho;
-    scalar Vm2 = Vm*Vm;
-    scalar Vm3 = Vm2*Vm;
-    scalar Vm4 = Vm3*Vm;
 
-    return
-    (
-        2*a(T)*
-        (
-            b3_ - b2_*Vm - b_*Vm2 + Vm3
-        )
-      - this->RR()*T*
-        (
-            b4_ - 4*b3_*Vm + 2*b2_*Vm2 + 4*b_*Vm3 + Vm4
-        )
-    )
-   /pow(b_ - Vm, 8);
+    scalar bmVm = b_ - Vm;
+    scalar bmVm2 = bmVm*bmVm;
+    scalar bpVm = b_ + Vm;
+
+    scalar t1 = bmVm2 - 2*Vm*Vm;
+    scalar t12 = t1*t1;
+
+    return (2*a(T)*bmVm2*bpVm - this->RR()*T*t12)/(bmVm2*t12);
 }
 
 
@@ -285,7 +274,7 @@ inline scalar pengRobinson::integral_dpdT_dv
 ) const
 {
     scalar Vm = this->W()/rho;
-    scalar root2=pow(2, 0.5);
+    scalar root2 = pow(2, 0.5);
 
     return
       - root2*dadT(T)*log(b_*(1 - root2) + Vm)/(4*b_)
@@ -308,23 +297,16 @@ inline scalar pengRobinson::d2pdv2(const scalar rho, const scalar T) const
 {
     scalar Vm = this->W()/rho;
     scalar Vm2 = Vm*Vm;
-    scalar Vm3 = Vm2*Vm;
-    scalar Vm4 = Vm3*Vm;
-    scalar Vm5 = Vm4*Vm;
-    scalar Vm6 = Vm5*Vm;
 
-    return 2*
-    (
-        a(T)*
-        (
-            5*b5_ - 9*b4_*Vm + 4*b2_*Vm3 + 3*b_*Vm4 - 3*Vm5
-        )
-      - this->RR()*T*
-        (
-            b6_ - 6*b5_*Vm + 9*b4_*Vm2 + 4*b3_*Vm3 - 9*b2_*Vm4 - 6*b_*Vm5 - Vm6
-        )
-    )
-    /pow(b_ - Vm, 9);
+    scalar bmVm = b_ - Vm;
+    scalar bmVm2 = bmVm*bmVm;
+    scalar bmVm3 = bmVm2*bmVm;
+
+    scalar t1 = bmVm2 - 2*Vm2;
+    scalar t12 = t1*t1;
+    scalar t13 = t12*t1;
+
+    return 2*(a(T)*bmVm3*(5*b_ + 6*b_*Vm + 3*Vm2) - this->RR()*T*t13)/(bmVm*t1);
 }
 
 
@@ -360,16 +342,16 @@ inline scalar pengRobinson::d2pdvdT
 ) const
 {
     scalar Vm = this->W()/rho;
-    scalar Vm2 = Vm*Vm;
-    scalar Vm3 = Vm2*Vm;
-    scalar Vm4 = Vm3*Vm;
 
-    return
-    (
-        2*dadT(T)*(b3_ - b2_*Vm - b_*Vm2 + Vm3)
-      - this->RR()*(b4_ - 4*b3_*Vm + 2*b2_*Vm2 + 4*b_*Vm3 + Vm4)
-    )
-   /pow(b_ - Vm, 6);
+    scalar bmVm = b_ - Vm;
+    scalar bmVm2 = bmVm*bmVm;
+    scalar bpVm = b_ + Vm;
+
+    scalar t1 = bmVm2 - 2*Vm*Vm;
+    scalar t12 = t1*t1;
+    scalar t14 = t12*t12;
+
+    return (2*dadT(T)*bmVm2*bpVm - this->RR()*t14)/(bmVm2*t12);
 }
 
 
diff --git a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/redlichKwong/redlichKwong.C b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/redlichKwong/redlichKwong.C
index 07f0ddc03..ea8af3415 100755
--- a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/redlichKwong/redlichKwong.C
+++ b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/redlichKwong/redlichKwong.C
@@ -45,8 +45,6 @@ Foam::redlichKwong::redlichKwong(Istream& is)
     a_(0.42748*pow(this->RR(), 2)*pow(Tcrit_, 2.5)/pcrit_),
     b_(0.08664*this->RR()*Tcrit_/pcrit_),
     b2_(pow(b_,2)),
-    b3_(pow(b_,3)),
-    b5_(pow(b_,5)),
     //CL: Only uses the default values
     rhoMin_(1e-3),
     rhoMax_(1500),
diff --git a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/redlichKwong/redlichKwong.H b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/redlichKwong/redlichKwong.H
index 300afd1a4..73473af84 100755
--- a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/redlichKwong/redlichKwong.H
+++ b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/redlichKwong/redlichKwong.H
@@ -69,8 +69,6 @@ class redlichKwong
 
     //CL: pow of constants b_ used in the code e.g. b2_=b*b;
     scalar b2_;
-    scalar b3_;
-    scalar b5_;
 
     //CL: rhoMin and rhoMax are only used as boundaries for the bisection methode (see rho function)
     scalar rhoMin_;
diff --git a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/redlichKwong/redlichKwongI.H b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/redlichKwong/redlichKwongI.H
index 9511a3733..d5d40c9a4 100755
--- a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/redlichKwong/redlichKwongI.H
+++ b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/redlichKwong/redlichKwongI.H
@@ -57,8 +57,6 @@ inline redlichKwong::redlichKwong(const word& name, const redlichKwong& rk)
     a_(rk.a_),
     b_(rk.b_),
     b2_(rk.b2_),
-    b3_(rk.b3_),
-    b5_(rk.b5_),
     rhoMin_(rk.rhoMin_),
     rhoMax_(rk.rhoMax_),
     rhostd_(rk.rhostd_)
@@ -132,15 +130,13 @@ inline scalar redlichKwong::dpdv(const scalar rho, const scalar T) const
 {
     scalar Vm = this->W()/rho;
     scalar Vm2 = Vm*Vm;
-    scalar Vm3 = Vm2*Vm;
     scalar T05 = pow(T, 0.5);
-    scalar bpVm2 = pow(b_ + Vm, 2);
+    scalar bmVm = b_ - Vm;
+    scalar bmVm2 = bmVm*bmVm;
+    scalar bpVm = b_ + Vm;
+    scalar bpVm2 = bpVm*bpVm;
 
-    return
-    (
-        a_*(b3_ - 3*b_*Vm2 + 2*Vm3) - this->RR()*T*T05*Vm2*bpVm2
-    )
-   /(T05*Vm2*bpVm2*pow(b_ - Vm, 2));
+    return (a_*bmVm2*bpVm - this->RR()*T*T05*Vm2*bpVm2)/(T05*Vm2*bpVm2*bmVm2);
 }
 
 
@@ -181,8 +177,7 @@ inline scalar redlichKwong::integral_p_dv
 {
     scalar Vm = this->W()/rho;
 
-    return this->RR()*T*log(Vm - b_) + a_/(b_*sqrt(T))*(log(b_ + Vm) - log(Vm));
-    //return this->RR()*T*log(Vm - b_) + (a_*log(b_ + Vm))/(b_*sqrt(T)) - (a_*log(Vm))/(b_*sqrt(T));
+    return this->RR()*T*log(Vm - b_) + a_/(b_*sqrt(T))*log((b_ + Vm)/Vm);
 }
 
 
@@ -196,8 +191,7 @@ inline scalar redlichKwong::integral_dpdT_dv
 {
     scalar Vm = this->W()/rho;
 
-    return this->RR()*log(Vm - b_) - a_/(2*b_*pow(T, 1.5))*(log(b_ + Vm) - log(Vm));
-    //return this->RR()*log(Vm - b_) - (a_*log(b_ + Vm))/(2*b_*T15_) + (a_*log(Vm))/(2*b_*T15_);
+    return this->RR()*log(Vm - b_) - a_/(2*b_*pow(T, 1.5))*log((b_ + Vm)/Vm);
 }
 
 
@@ -205,6 +199,7 @@ inline scalar redlichKwong::integral_dpdT_dv
 inline scalar redlichKwong::d2pdT2(const scalar rho, const scalar T) const
 {
     scalar Vm = this->W()/rho;
+
     return -0.75*a_/(pow(T, 2.5)*Vm*(b_ + Vm));
 }
 
@@ -215,19 +210,16 @@ inline scalar redlichKwong::d2pdv2(const scalar rho, const scalar T) const
     scalar Vm = this->W()/rho;
     scalar Vm2 = Vm*Vm;
     scalar Vm3 = Vm2*Vm;
-    scalar Vm4 = Vm3*Vm;
-    scalar Vm5 = Vm4*Vm;
     scalar T05 = pow(T, 0.5);
+    scalar bmVm = b_ - Vm;
+    scalar bmVm2 = bmVm*bmVm;
+    scalar bmVm3 = bmVm2*bmVm;
+    scalar bpVm = b_ + Vm;
+    scalar bpVm2 = bpVm*bpVm;
+    scalar bpVm3 = bpVm2*bpVm;
 
-    return
-    (
-        2*
-        (
-            a_*(b5_ - 3*b3_*Vm2 - b2_*Vm3 + 6*b_*Vm4 - 3*Vm5)
-          + this->RR()*T*T05*Vm3*(b3_ + 3*b2_*Vm + 3*b_*Vm2 + Vm3)
-        )
-        /(T05*Vm3*pow((b_ + Vm), 3)*pow(Vm - b_, 3))
-    );
+    return 2*(a_*bmVm3*(b2_ + 3*b_*Vm + 3*Vm2) + this->RR()*T*T05*Vm3*bpVm3)
+        /(T05*Vm3*bpVm3*bmVm3);
 }
 
 
@@ -256,18 +248,14 @@ inline scalar redlichKwong::d2pdvdT(const scalar rho, const scalar T) const
 {
     scalar Vm = this->W()/rho;
     scalar Vm2 = Vm*Vm;
-    scalar Vm3 = Vm2*Vm;
     scalar T15_ = pow(T, 1.5);
+    scalar bmVm = b_ - Vm;
+    scalar bmVm2 = bmVm*bmVm;
+    scalar bpVm = b_ + Vm;
+    scalar bpVm2 = bpVm*bpVm;
 
-    return
-    -(
-        0.5*
-        (
-            a_*(b3_ - 3*b_*Vm2 + 2*Vm3)
-          + 2*this->RR()*T15_*Vm2*(b2_ + 2*b_*Vm + Vm2)
-        )
-    )
-    /(T15_*Vm2*pow(b_ + Vm, 2)*pow(b_ - Vm, 2));
+    return -(0.5*a_*bmVm2*(b_ + 2*Vm) + this->RR()*T15_*Vm2*bpVm2)
+        /(T15_*Vm2*bpVm2*bmVm2);
 }
 
 
@@ -279,9 +267,10 @@ inline scalar redlichKwong::integral_d2pdT2_dv
     const scalar T
 ) const
 {
-    scalar T25_=pow(T,2.5);
     scalar Vm = this->W()/rho;
-    return 0.75*a_*log(b_  +  Vm)/(T25_*b_) - 0.75*a_*log(Vm)/(T25_*b_);
+    scalar T25 = pow(T, 2.5);
+
+    return 0.75*a_/(T25*b_)*log((b_ +  Vm)/Vm);
 }
 
 
diff --git a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwong.C b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwong.C
index 449744218..1b2a60056 100755
--- a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwong.C
+++ b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwong.C
@@ -47,17 +47,15 @@ Foam::soaveRedlichKwong::soaveRedlichKwong(Istream& is)
     b_(0.08664*this->RR()*Tcrit_/pcrit_),
     n_(0.48508 + 1.55171*azentricFactor_ - 0.15613*pow(azentricFactor_, 2)),
     b2_(b_*b_),
-    b3_(b2_*b_),
-    b5_(b2_*b3_),
     //CL: Only uses the default values
     rhoMin_(1e-3),
     rhoMax_(1500),
-    // Starting GUESS for the density by ideal gas law
-    rhostd_(this->rho(this->Pstd(), this->Tstd(), this->Pstd()/(this->Tstd()*this->R()))),
     aSave(0.0),
     daSave(0.0),
     d2aSave(0.0),
-    TSave(0.0)
+    TSave(0.0),
+    // Starting GUESS for the density by ideal gas law
+    rhostd_(this->rho(this->Pstd(), this->Tstd(), this->Pstd()/(this->Tstd()*this->R())))
 {
     is.check("soaveRedlichKwong::soaveRedlichKwong(Istream& is)");
 }
@@ -83,12 +81,12 @@ soaveRedlichKwong::soaveRedlichKwong(const dictionary& dict)
     //CL: therefore, rho can be larger than rhoMax and smaller than rhoMin
     rhoMin_(dict.subDict("equationOfState").lookupOrDefault("rhoMin",1e-3)),
     rhoMax_(dict.subDict("equationOfState").lookupOrDefault("rhoMax",1500)),
-    // Starting GUESS for the density by ideal gas law
-    rhostd_(this->rho(this->Pstd(), this->Tstd(), this->Pstd()/(this->Tstd()*this->R()))),
     aSave(0.0),
     daSave(0.0),
     d2aSave(0.0),
-    TSave(0.0)
+    TSave(0.0),
+    // Starting GUESS for the density by ideal gas law
+    rhostd_(this->rho(this->Pstd(), this->Tstd(), this->Pstd()/(this->Tstd()*this->R())))
 {
     is.check("soaveRedlichKwong::soaveRedlichKwong(Istream& is)");
 }
diff --git a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwong.H b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwong.H
index a33b8da83..adec662f1 100755
--- a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwong.H
+++ b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwong.H
@@ -76,16 +76,11 @@ class soaveRedlichKwong
 
     //CL: pow of constants b_ used in the code e.g. b2_=b*b;
     scalar b2_;
-    scalar b3_;
-    scalar b5_;
 
     //CL: rhoMin and rhoMax are only used as boundaries for the bisection methode (see rho function)
     scalar rhoMin_;
     scalar rhoMax_;
 
-    //- Density @STD, initialise after a0, b!
-    scalar rhostd_;
-
     //CL: Variables to save the values of a, dadT and d2adT2 of the mixture
     //CL: Variables must corrected for changing temperatures
     mutable scalar aSave;
@@ -95,6 +90,9 @@ class soaveRedlichKwong
     //CL: save the temperature for which the save coefficients (aSave, daSave, d2aSave) are correct
     mutable scalar TSave;
 
+    //- Density @STD, initialise after a0, b!
+    scalar rhostd_;
+
     //Protected functions
     //CL: function updates the coefficients (aSave, daSave, d2aSave)
     inline void updateModelCoefficients(const  scalar T) const;
diff --git a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwongI.H b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwongI.H
index 9a478d085..106fd1f6f 100755
--- a/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwongI.H
+++ b/src/thermophysicalModels/specie/equationOfState/cubicEquationOfState/soaveRedlichKwong/soaveRedlichKwongI.H
@@ -60,15 +60,13 @@ inline soaveRedlichKwong::soaveRedlichKwong(const word& name, const soaveRedlich
     b_(srk.b_),
     n_(srk.n_),
     b2_(srk.b2_),
-    b3_(srk.b3_),
-    b5_(srk.b5_),
     rhoMin_(srk.rhoMin_),
     rhoMax_(srk.rhoMax_),
-    rhostd_(srk.rhostd_),
     aSave(0.0),
     daSave(0.0),
     d2aSave(0.0),
-    TSave(0.0)
+    TSave(0.0),
+    rhostd_(srk.rhostd_)
 {}
 
 
@@ -129,7 +127,7 @@ inline void soaveRedlichKwong::updateModelCoefficients(const scalar T) const
 inline scalar soaveRedlichKwong::a(const scalar T) const
 {
     //CL: check if a has already been calculated for this temperature
-    if(TSave==T)
+    if(TSave == T)
     {
         return aSave;
     }
@@ -146,7 +144,7 @@ inline scalar soaveRedlichKwong::a(const scalar T) const
 inline scalar soaveRedlichKwong::dadT(const scalar T) const
 {
     // check if a has already been calculated for this temperature
-    if(TSave==T)
+    if(TSave == T)
     {
         return daSave;
     }
@@ -163,7 +161,7 @@ inline scalar soaveRedlichKwong::dadT(const scalar T) const
 inline scalar soaveRedlichKwong::d2adT2(const scalar T) const
 {
     // check if a has already been calculated for this temperature
-    if(TSave==T)
+    if(TSave == T)
     {
         return d2aSave;
     }
@@ -208,13 +206,12 @@ inline scalar soaveRedlichKwong::dpdv(const scalar rho, const scalar T) const
 {
     scalar Vm = this->W()/rho;
     scalar Vm2 = Vm*Vm;
-    scalar Vm3 = Vm2*Vm;
+    scalar bmVm = b_ - Vm;
+    scalar bmVm2 = bmVm*bmVm;
+    scalar bpVm = b_ + Vm;
+    scalar bpVm2 = bpVm*bpVm;
 
-    return
-    (
-        a(T)*(b3_ - 3*b_*Vm2 + 2*Vm3) - this->RR()*T*Vm2*(b2_ + 2*b_*Vm + Vm2)
-    )
-   /(Vm2*pow(b_ + Vm, 2)*pow(b_ - Vm, 2));
+    return(a(T)*bmVm2*(b_ + 2*Vm) - this->RR()*T*Vm2*bpVm2)/(Vm2*bpVm2*bmVm2);
 }
 
 
@@ -255,8 +252,7 @@ inline scalar soaveRedlichKwong::integral_p_dv
 {
     scalar Vm = this->W()/rho;
 
-    return this->RR()*T*log(Vm - b_) + a(T)/b_*(log(b_ + Vm) - log(Vm));
-    //return this->RR()*T*log(Vm - b_) + a(T)*log(b_ + Vm)/b_ - a(T)*log(Vm)/b_;
+    return this->RR()*T*log(Vm - b_) + a(T)/b_*log((b_ + Vm)/Vm);
 }
 
 
@@ -271,8 +267,7 @@ inline scalar soaveRedlichKwong::integral_dpdT_dv
 {
     scalar Vm = this->W()/rho;
 
-    return this->RR()*log(Vm - b_) + dadT(T)/b_*(log(b_ + Vm) - log(Vm));
-    //return this->RR()*log(Vm - b_) + dadT(T)*log(b_ + Vm)/b_ - dadT(T)*log(Vm)/b_;
+    return this->RR()*log(Vm - b_) + dadT(T)/b_*log((b_ + Vm)/Vm);
 }
 
 
@@ -291,15 +286,15 @@ inline scalar soaveRedlichKwong::d2pdv2(const scalar rho, const scalar T) const
     scalar Vm = this->W()/rho;
     scalar Vm2 = Vm*Vm;
     scalar Vm3 = Vm2*Vm;
-    scalar Vm4 = Vm3*Vm;
-    scalar Vm5 = Vm4*Vm;
+    scalar bmVm = b_ - Vm;
+    scalar bmVm2 = bmVm*bmVm;
+    scalar bmVm3 = bmVm2*bmVm;
+    scalar bpVm = b_ + Vm;
+    scalar bpVm2 = bpVm*bpVm;
+    scalar bpVm3 = bpVm2*bpVm;
 
-    return 2*
-    (
-        a(T)*(b5_ - 3*b3_*Vm2 - b2_*Vm3 + 6*b_*Vm4 - 3*Vm5)
-      + this->RR()*T*Vm3*(b3_ + 3*b2_*Vm + 3*b_*Vm2 + Vm3)
-    )
-   /(Vm3*pow(b_ + Vm, 3)*pow(Vm - b_, 3));
+    return 2*(a(T)*bmVm3*(b2_ + 3*b_*Vm + 3*Vm2) + this->RR()*T*Vm3*bpVm3)
+        /(Vm3*bpVm3*bmVm3);
 }
 
 
@@ -336,13 +331,13 @@ inline scalar soaveRedlichKwong::d2pdvdT
 {
     scalar Vm = this->W()/rho;
     scalar Vm2 = Vm*Vm;
-    scalar Vm3 = Vm2*Vm;
+    scalar bmVm = b_ - Vm;
+    scalar bmVm2 = bmVm*bmVm;
+    scalar bpVm = b_ + Vm;
+    scalar bpVm2 = bpVm*bpVm;
 
-    return
-    (
-        dadT(T)*(b3_ - 3*b_*Vm2 + 2*Vm3) - this->RR()*Vm2*(b2_ + 2*b_*Vm + Vm2)
-    )
-    /(Vm2*pow(b_ + Vm, 2)*pow(b_ - Vm, 2));
+    return (dadT(T)*bmVm2*(b_ + 2*Vm) - this->RR()*Vm2*bmVm2)
+        /(Vm2*bpVm2*bmVm2);
 }
 
 
@@ -356,7 +351,7 @@ inline scalar soaveRedlichKwong::integral_d2pdT2_dv
 {
     scalar Vm = this->W()/rho;
 
-    return d2adT2(T)*log(b_ + Vm)/b_ - d2adT2(T)*log(Vm)/b_;
+    return d2adT2(T)/b_*log((b_ +  Vm)/Vm);
 }
 
 
diff --git a/tutorials/compressible/realFluidPisoFoam/ras/t-junction/constant/thermophysicalProperties b/tutorials/compressible/realFluidPisoFoam/ras/t-junction/constant/thermophysicalProperties
index 6df25c688..74073cbc4 100755
--- a/tutorials/compressible/realFluidPisoFoam/ras/t-junction/constant/thermophysicalProperties
+++ b/tutorials/compressible/realFluidPisoFoam/ras/t-junction/constant/thermophysicalProperties
@@ -18,14 +18,14 @@ FoamFile
 
 //CL: List of possible real gas models
 
-thermoType  realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<redlichKwong>>>>>;
-mixture         CO2 1 44.01 73.773e5 304.13 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801e-13 1.4792e-06 116;
+//thermoType  realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<redlichKwong>>>>>;
+//mixture         CO2 1 44.01 73.773e5 304.13 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801e-13 1.4792e-06 116;
 
 //thermoType  realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<aungierRedlichKwong>>>>>;
 //mixture CO2 1 44.01 73.773e5 304.13 467.6 0.22394  49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801E-13 1.4792e-06 116;
 
-//thermoType  realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<pengRobinson>>>>>;
-//mixture        CO2 1 44.01 73.773e5 304.13 0.22394 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801e-13 1.4792e-06 116;
+thermoType  realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<pengRobinson>>>>>;
+mixture        CO2 1 44.01 73.773e5 304.13 0.22394 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801e-13 1.4792e-06 116;
 
 //thermoType  realGasHThermo<pureMixture<sutherlandTransport<realGasSpecieThermo<nasaHeatCapacityPolynomial<soaveRedlichKwong>>>>>;
 //mixture         CO2 1 44.01 73.773e5 304.13 0.22394 49436.5054 -626.411601 5.30172524 0.002503813816 -0.0000002127308728 -0.000000000768998878 2.849677801e-13 1.4792e-06 116;
diff --git a/tutorials/compressible/realFluidPisoFoam/ras/t-junction/system/controlDict b/tutorials/compressible/realFluidPisoFoam/ras/t-junction/system/controlDict
index 9a3b44200..f3e29a18e 100644
--- a/tutorials/compressible/realFluidPisoFoam/ras/t-junction/system/controlDict
+++ b/tutorials/compressible/realFluidPisoFoam/ras/t-junction/system/controlDict
@@ -45,7 +45,7 @@ timePrecision   10;
 
 adjustTimeStep  yes;
 
-maxCo           0.25;
+maxCo           0.1;
 
 maxDeltaT       1e-2;