Merge remote-tracking branch 'origin/merge-foundation' into adjust-config-foundation-merge

applications/solvers/combustion/reactingFoam/createFields.H

@@ -1,3 +1,5 @@
+#include "createRDeltaT.H"
+
 Info<< "Creating reaction model\n" << endl;
 
 autoPtr<combustionModels::psiCombustionModel> reaction

applications/solvers/combustion/reactingFoam/reactingFoam.C

@@ -52,7 +52,6 @@ int main(int argc, char *argv[])
     #include "createMesh.H"
     #include "createControl.H"
     #include "createTimeControls.H"
-    #include "createRDeltaT.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
     #include "createFieldRefs.H"

applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/createFields.H

@@ -1,3 +1,5 @@
+#include "createRDeltaT.H"
+
 Info<< "Creating reaction model\n" << endl;
 
 autoPtr<combustionModels::rhoCombustionModel> reaction

applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/rhoReactingBuoyantFoam.C

@@ -53,7 +53,6 @@ int main(int argc, char *argv[])
     #include "createMesh.H"
     #include "createControl.H"
     #include "createTimeControls.H"
-    #include "createRDeltaT.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
     #include "createFieldRefs.H"

applications/solvers/combustion/reactingFoam/rhoReactingFoam/createFields.H

@@ -1,3 +1,5 @@
+#include "createRDeltaT.H"
+
 Info<< "Creating reaction model\n" << endl;
 
 autoPtr<combustionModels::rhoCombustionModel> reaction

applications/solvers/combustion/reactingFoam/rhoReactingFoam/rhoReactingFoam.C

@@ -53,7 +53,6 @@ int main(int argc, char *argv[])
     #include "createMesh.H"
     #include "createControl.H"
     #include "createTimeControls.H"
-    #include "createRDeltaT.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
     #include "createFieldRefs.H"

applications/solvers/compressible/rhoCentralFoam/createFields.H

@@ -1,3 +1,5 @@
+#include "createRDeltaT.H"
+
 Info<< "Reading thermophysical properties\n" << endl;
 
 autoPtr<psiThermo> pThermo

applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.C

@@ -39,6 +39,8 @@ Description
 #include "turbulentFluidThermoModel.H"
 #include "fixedRhoFvPatchScalarField.H"
 #include "directionInterpolate.H"
+#include "localEulerDdtScheme.H"
+#include "fvcSmooth.H"
 #include "motionSolver.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C

@@ -54,7 +54,6 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "createFieldRefs.H"
     #include "createTimeControls.H"
-    #include "createRDeltaT.H"
 
     turbulence->validate();
 

applications/solvers/compressible/rhoPimpleFoam/createFields.H

@@ -1,3 +1,5 @@
+#include "createRDeltaT.H"
+
 Info<< "Reading thermophysical properties\n" << endl;
 
 autoPtr<psiThermo> pThermo

applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/rhoPimpleDyMFoam.C

@@ -57,7 +57,6 @@ int main(int argc, char *argv[])
     #include "createTime.H"
     #include "createDynamicFvMesh.H"
     #include "createControl.H"
-    #include "createRDeltaT.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
     #include "createFieldRefs.H"

applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C

@@ -56,7 +56,6 @@ int main(int argc, char *argv[])
     #include "createMesh.H"
     #include "createControl.H"
     #include "createTimeControls.H"
-    #include "createRDeltaT.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
     #include "createFieldRefs.H"

applications/solvers/heatTransfer/chtMultiRegionFoam/include/setInitialMultiRegionDeltaT.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Global
-    setInitialDeltaT
+    setInitialMultiRegionDeltaT
 
 Description
     Set the initial timestep for the CHT MultiRegion solver.
@@ -48,7 +48,7 @@ if (adjustTimeStep)
             min
             (
                 min(maxCo/CoNum, maxDi/DiNum)*runTime.deltaT().value(),
-                maxDeltaT
+                min(runTime.deltaTValue(), maxDeltaT)
             )
         );
         Info<< "deltaT = " <<  runTime.deltaT().value() << endl;

applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointOutletPressure/adjointOutletPressureFvPatchScalarField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -63,13 +63,8 @@ adjointOutletPressureFvPatchScalarField
     const dictionary& dict
 )
 :
-    fixedValueFvPatchScalarField(p, iF)
+    fixedValueFvPatchScalarField(p, iF, dict)
-{
+{}
-    fvPatchField<scalar>::operator=
-    (
-        scalarField("value", dict, p.size())
-    );
-}
 
 
 Foam::adjointOutletPressureFvPatchScalarField::

applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointOutletVelocity/adjointOutletVelocityFvPatchVectorField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -50,10 +50,8 @@ adjointOutletVelocityFvPatchVectorField
     const dictionary& dict
 )
 :
-    fixedValueFvPatchVectorField(p, iF)
+    fixedValueFvPatchVectorField(p, iF, dict)
-{
+{}
-    fvPatchVectorField::operator=(vectorField("value", dict, p.size()));
-}
 
 
 Foam::adjointOutletVelocityFvPatchVectorField::

applications/solvers/lagrangian/DPMFoam/DPMTurbulenceModels/DPMTurbulenceModels.C

@@ -28,8 +28,8 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "makeTurbulenceModel.H"
 
-#include "laminar.H"
+#include "laminarModel.H"
-#include "turbulentTransportModel.H"
+#include "RASModel.H"
 #include "LESModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -52,6 +52,10 @@ makeBaseTurbulenceModel
     singlePhaseTransportModel
 );
 
+#define makeLaminarModel(Type)                                                 \
+    makeTemplatedTurbulenceModel                                               \
+    (singlePhaseTransportModelPhaseIncompressibleTurbulenceModel, laminar, Type)
+
 #define makeRASModel(Type)                                                     \
     makeTemplatedTurbulenceModel                                               \
     (singlePhaseTransportModelPhaseIncompressibleTurbulenceModel, RAS, Type)
@@ -60,6 +64,9 @@ makeBaseTurbulenceModel
     makeTemplatedTurbulenceModel                                               \
     (singlePhaseTransportModelPhaseIncompressibleTurbulenceModel, LES, Type)
 
+#include "Stokes.H"
+makeLaminarModel(Stokes);
+
 #include "kEpsilon.H"
 makeRASModel(kEpsilon);
 

applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C

@@ -56,7 +56,6 @@ int main(int argc, char *argv[])
     #include "createMesh.H"
     #include "createControl.H"
     #include "createTimeControls.H"
-    #include "createRDeltaT.H"
     #include "createFields.H"
     #include "createFieldRefs.H"
     #include "createFvOptions.H"

applications/solvers/lagrangian/coalChemistryFoam/createFields.H

@@ -1,3 +1,5 @@
+#include "createRDeltaT.H"
+
 #include "readGravitationalAcceleration.H"
 
 Info<< "Creating combustion model\n" << endl;

applications/solvers/lagrangian/reactingParcelFilmFoam/Make/options

@@ -1,4 +1,5 @@
 EXE_INC = \
+    -I. \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I${LIB_SRC}/sampling/lnInclude \
     -I${LIB_SRC}/meshTools/lnInclude \
@@ -25,7 +26,6 @@ EXE_INC = \
     -I$(LIB_SRC)/combustionModels/lnInclude \
     -I$(FOAM_SOLVERS)/combustion/reactingFoam
 
-
 EXE_LIBS = \
     -lfiniteVolume \
     -lfvOptions \

applications/solvers/lagrangian/reactingParcelFilmFoam/createFieldRefs.H

@@ -1 +1,5 @@
 const label inertIndex(composition.species()[inertSpecie]);
+
+const volScalarField& T = thermo.T();
+const volScalarField& psi = thermo.psi();
+filmModelType& surfaceFilm = tsurfaceFilm();

applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H

@@ -37,8 +37,6 @@ volScalarField rho
 );
 
 volScalarField& p = thermo.p();
-const volScalarField& T = thermo.T();
-const volScalarField& psi = thermo.psi();
 
 Info<< "\nReading field U\n" << endl;
 volVectorField U

applications/solvers/lagrangian/reactingParcelFilmFoam/createSurfaceFilmModel.H

@@ -3,5 +3,3 @@ Info<< "\nConstructing surface film model" << endl;
 typedef regionModels::surfaceFilmModels::surfaceFilmModel filmModelType;
 
 autoPtr<filmModelType> tsurfaceFilm(filmModelType::New(mesh, g));
-filmModelType& surfaceFilm = tsurfaceFilm();
-

applications/solvers/lagrangian/reactingParcelFoam/createFields.H

@@ -1,3 +1,5 @@
+#include "createRDeltaT.H"
+
 #include "readGravitationalAcceleration.H"
 
 Info<< "Creating combustion model\n" << endl;

applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C

@@ -55,7 +55,6 @@ int main(int argc, char *argv[])
     #include "createMesh.H"
     #include "createControl.H"
     #include "createTimeControls.H"
-    #include "createRDeltaT.H"
     #include "createFields.H"
     #include "createFieldRefs.H"
     #include "createFvOptions.H"

applications/solvers/multiphase/MPPICInterFoam/CompressibleTwoPhaseMixtureTurbulenceModels/CompressibleTwoPhaseMixtureTurbulenceModels.C

@@ -28,7 +28,6 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "makeTurbulenceModel.H"
 
-#include "laminar.H"
 #include "turbulentTransportModel.H"
 #include "LESModel.H"
 
@@ -50,13 +49,32 @@ makeBaseTurbulenceModel
     immiscibleIncompressibleTwoPhaseMixture
 );
 
+#define makeLaminarModel(Type)                                                 \
+    makeTemplatedTurbulenceModel                                               \
+    (                                                                          \
+       immiscibleIncompressibleTwoPhaseMixturePhaseCompressibleTurbulenceModel,\
+        laminar,                                                               \
+        Type                                                                   \
+    )
+
 #define makeRASModel(Type)                                                     \
     makeTemplatedTurbulenceModel                                               \
-    (immiscibleIncompressibleTwoPhaseMixturePhaseCompressibleTurbulenceModel, RAS, Type)
+    (                                                                          \
+       immiscibleIncompressibleTwoPhaseMixturePhaseCompressibleTurbulenceModel,\
+        RAS,                                                                   \
+        Type                                                                   \
+    )
 
 #define makeLESModel(Type)                                                     \
     makeTemplatedTurbulenceModel                                               \
-    (immiscibleIncompressibleTwoPhaseMixturePhaseCompressibleTurbulenceModel, LES, Type)
+    (                                                                          \
+       immiscibleIncompressibleTwoPhaseMixturePhaseCompressibleTurbulenceModel,\
+        LES,                                                                   \
+        Type                                                                   \
+    )
+
+#include "Stokes.H"
+makeLaminarModel(Stokes);
 
 #include "kEpsilon.H"
 makeRASModel(kEpsilon);
@@ -70,4 +88,5 @@ makeLESModel(kEqn);
 #include "kOmega.H"
 makeRASModel(kOmega);
 
+
 // ************************************************************************* //

applications/solvers/multiphase/driftFluxFoam/compressibleTurbulenceModels.C

@@ -28,7 +28,7 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "makeTurbulenceModel.H"
 
-#include "laminar.H"
+#include "laminarModel.H"
 #include "RASModel.H"
 #include "LESModel.H"
 
@@ -52,6 +52,14 @@ makeBaseTurbulenceModel
     incompressibleTwoPhaseInteractingMixture
 );
 
+#define makeLaminarModel(Type)                                                 \
+    makeTemplatedTurbulenceModel                                               \
+    (                                                                          \
+        incompressibleTwoPhaseInteractingMixtureCompressibleTurbulenceModel,   \
+        laminar,                                                               \
+        Type                                                                   \
+    )
+
 #define makeRASModel(Type)                                                     \
     makeTemplatedTurbulenceModel                                               \
     (                                                                          \
@@ -68,6 +76,9 @@ makeBaseTurbulenceModel
         Type                                                                   \
     )
 
+#include "Stokes.H"
+makeLaminarModel(Stokes);
+
 #include "kEpsilon.H"
 makeRASModel(kEpsilon);
 

applications/solvers/multiphase/interFoam/createFields.H

@@ -1,3 +1,5 @@
+#include "createRDeltaT.H"
+
 Info<< "Reading field p_rgh\n" << endl;
 volScalarField p_rgh
 (

applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C

@@ -63,7 +63,6 @@ int main(int argc, char *argv[])
     #include "createControl.H"
     #include "createTimeControls.H"
     #include "createDyMControls.H"
-    #include "createRDeltaT.H"
     #include "createFields.H"
     #include "createFvOptions.H"
 

applications/solvers/multiphase/interFoam/interFoam.C

@@ -65,7 +65,6 @@ int main(int argc, char *argv[])
     #include "createMesh.H"
     #include "createControl.H"
     #include "createTimeControls.H"
-    #include "createRDeltaT.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
     #include "createFvOptions.H"

applications/solvers/multiphase/interFoam/interMixingFoam/createFields.H

@@ -1,3 +1,5 @@
+#include "createRDeltaT.H"
+
 Info<< "Reading field p_rgh\n" << endl;
 volScalarField p_rgh
 (

applications/solvers/multiphase/interFoam/interMixingFoam/interMixingFoam.C

@@ -55,7 +55,6 @@ int main(int argc, char *argv[])
     #include "createMesh.H"
     #include "createControl.H"
     #include "createTimeControls.H"
-    #include "createRDeltaT.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
     #include "createFvOptions.H"

applications/solvers/multiphase/reactingEulerFoam/interfacialModels/dragModels/Tenneti/Tenneti.H

@@ -58,6 +58,7 @@ namespace dragModels
 {
 
 class SchillerNaumann;
+
 /*---------------------------------------------------------------------------*\
                            Class Tenneti Declaration
 \*---------------------------------------------------------------------------*/
@@ -68,7 +69,7 @@ class Tenneti
 {
     // Private data
 
-        //- Ergun drag model
+        //- SchillerNaumann drag model
         autoPtr<SchillerNaumann> SchillerNaumann_;
 
         //- Residual Reynolds Number

applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/createFields.H

@@ -1,3 +1,4 @@
+#include "createRDeltaT.H"
 #include "readGravitationalAcceleration.H"
 #include "readhRef.H"
 

applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseCompressibleTurbulenceModels/multiphaseCompressibleTurbulenceModels.C

@@ -27,7 +27,7 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "makeTurbulenceModel.H"
 
-#include "laminar.H"
+#include "laminarModel.H"
 #include "RASModel.H"
 #include "LESModel.H"
 
@@ -53,6 +53,10 @@ makeBaseTurbulenceModel
     phaseModel
 );
 
+#define makeLaminarModel(Type)                                                 \
+    makeTemplatedLaminarModel                                                  \
+    (phaseModelPhaseCompressibleTurbulenceModel, laminar, Type)
+
 #define makeRASModel(Type)                                                     \
     makeTemplatedTurbulenceModel                                               \
     (phaseModelPhaseCompressibleTurbulenceModel, RAS, Type)
@@ -61,6 +65,9 @@ makeBaseTurbulenceModel
     makeTemplatedTurbulenceModel                                               \
     (phaseModelPhaseCompressibleTurbulenceModel, LES, Type)
 
+#include "Stokes.H"
+makeLaminarModel(Stokes);
+
 #include "kEpsilon.H"
 makeRASModel(kEpsilon);
 

applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/reactingMultiphaseEulerFoam.C

@@ -53,7 +53,6 @@ int main(int argc, char *argv[])
     #include "createMesh.H"
     #include "createControl.H"
     #include "createTimeControls.H"
-    #include "createRDeltaT.H"
     #include "createFields.H"
     #include "createFieldRefs.H"
 

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/createFields.H

@@ -1,3 +1,4 @@
+#include "createRDeltaT.H"
 #include "readGravitationalAcceleration.H"
 #include "readhRef.H"
 

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/reactingTwoPhaseEulerFoam.C

@@ -54,7 +54,6 @@ int main(int argc, char *argv[])
     #include "createMesh.H"
     #include "createControl.H"
     #include "createTimeControls.H"
-    #include "createRDeltaT.H"
     #include "createFields.H"
     #include "createFieldRefs.H"
 

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/phaseCompressibleTurbulenceModels.C

@@ -27,7 +27,7 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "makeTurbulenceModel.H"
 
-#include "laminar.H"
+#include "laminarModel.H"
 #include "RASModel.H"
 #include "LESModel.H"
 
@@ -53,6 +53,10 @@ makeBaseTurbulenceModel
     phaseModel
 );
 
+#define makeLaminarModel(Type)                                                 \
+    makeTemplatedLaminarModel                                                  \
+    (phaseModelPhaseCompressibleTurbulenceModel, laminar, Type)
+
 #define makeRASModel(Type)                                                     \
     makeTemplatedTurbulenceModel                                               \
     (phaseModelPhaseCompressibleTurbulenceModel, RAS, Type)
@@ -61,6 +65,9 @@ makeBaseTurbulenceModel
     makeTemplatedTurbulenceModel                                               \
     (phaseModelPhaseCompressibleTurbulenceModel, LES, Type)
 
+#include "Stokes.H"
+makeLaminarModel(Stokes);
+
 #include "kEpsilon.H"
 makeRASModel(kEpsilon);
 

applications/solvers/multiphase/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/phaseCompressibleTurbulenceModels.C

@@ -32,7 +32,7 @@ License
 #include "ThermalDiffusivity.H"
 #include "EddyDiffusivity.H"
 
-#include "laminar.H"
+#include "laminarModel.H"
 #include "RASModel.H"
 #include "LESModel.H"
 
@@ -58,6 +58,10 @@ makeBaseTurbulenceModel
     phaseModel
 );
 
+#define makeLaminarModel(Type)                                                 \
+    makeTemplatedLaminarModel                                                  \
+    (phaseModelPhaseCompressibleTurbulenceModel, laminar, Type)
+
 #define makeRASModel(Type)                                                     \
     makeTemplatedTurbulenceModel                                               \
     (phaseModelPhaseCompressibleTurbulenceModel, RAS, Type)
@@ -66,6 +70,9 @@ makeBaseTurbulenceModel
     makeTemplatedTurbulenceModel                                               \
     (phaseModelPhaseCompressibleTurbulenceModel, LES, Type)
 
+#include "Stokes.H"
+makeLaminarModel(Stokes);
+
 #include "kEpsilon.H"
 makeRASModel(kEpsilon);
 

applications/utilities/mesh/generation/blockMesh/blockMesh.C

@@ -86,6 +86,27 @@ int main(int argc, char *argv[])
         "specify alternative dictionary for the blockMesh description"
     );
 
+    argList::addNote
+    (
+        "Block description\n"
+        "\n"
+        "  For a given block, the correspondence between the ordering of\n"
+        "  vertex labels and face labels is shown below.\n"
+        "  For vertex numbering in the sequence 0 to 7 (block, centre):\n"
+        "    faces 0 (f0) and 1 are left and right, respectively;\n"
+        "    faces 2 and 3 are bottom and top;\n"
+        "    and faces 4 and 5 are front the back:\n"
+        "\n"
+        "           4 ---- 5\n"
+        "      f3   |\\     |\\   f5\n"
+        "      |    | 7 ---- 6   \\\n"
+        "      |    0 |--- 1 |    \\\n"
+        "      |     \\|     \\|    f4\n"
+        "      f2     3 ---- 2\n"
+        "\n"
+        "            f0 ----- f1\n"
+    );
+
     #include "addRegionOption.H"
     #include "setRootCase.H"
     #include "createTime.H"

applications/utilities/mesh/generation/foamyMesh/foamyHexMeshBackgroundMesh/foamyHexMeshBackgroundMesh.C

@@ -3,7 +3,7 @@
  \\      /   F ield          | OpenFOAM: The Open Source CFD Toolbox
   \\    /    O peration      |
    \\  /     A nd            | Copyright (C) 2012-2016 OpenFOAM Foundation
-    \\/      M anipulation   |
+    \\/      M anipulation   | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -640,11 +640,14 @@ int main(int argc, char *argv[])
                 geometry,
                 surfaces
             );
+
             // Patch fields
+            volScalarField::Boundary& cellDistanceBf =
+                cellDistance.boundaryFieldRef();
             forAll(fvm.C().boundaryField(), patchi)
             {
                 const pointField& cc = fvm.C().boundaryField()[patchi];
-                fvPatchScalarField& fld = cellDistance.boundaryField()[patchi];
+                fvPatchScalarField& fld = cellDistanceBf[patchi];
                 scalarField patchDistSqr
                 (
                     fld.patch().patchInternalField(distSqr)

applications/utilities/miscellaneous/foamList/Make/options

@@ -6,7 +6,8 @@ EXE_INC = \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
 
 EXE_LIBS = \
     -lbarotropicCompressibilityModel \

applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposer.C

@@ -42,7 +42,7 @@ Foam::lagrangianFieldDecomposer::lagrangianFieldDecomposer
 )
 :
     procMesh_(procMesh),
-    positions_(procMesh, cloudName, false),
+    positions_(procMesh, cloudName, IDLList<passiveParticle>()),
     particleIndices_(lagrangianPositions.size())
 {
     label pi = 0;

applications/utilities/parallelProcessing/decomposePar/lagrangianFieldDecomposer.H

@@ -101,7 +101,6 @@ public:
             const label cloudI,
             const IOobjectList& lagrangianObjects,
             PtrList<PtrList<IOField<Type>>>& lagrangianFields
-//            PtrList<IOField<Type>>& lagrangianFields
         );
 
         template<class Type>
@@ -113,7 +112,6 @@ public:
             <
                 PtrList<CompactIOField<Field<Type>, Type>>
             >& lagrangianFields
-//            PtrList<CompactIOField<Field<Type>, Type >>& lagrangianFields
         );
 
 

applications/utilities/postProcessing/dataConversion/foamToGMV/gmvOutput.H

@@ -79,7 +79,7 @@ for (label i=0; i < nTypes; i++)
             )
         )
         {
-            Cloud<passiveParticle> particles(mesh);
+            Cloud<passiveParticle> particles(mesh, cloud::defaultName);
 
             IOobjectList objects(mesh, runTime.timeName(), cloud::prefix);
 

applications/utilities/postProcessing/graphics/ensightFoamReader/USERD_get_maxsize_info.H

@@ -73,7 +73,7 @@ int USERD_get_maxsize_info
         {
             // Get the maximum number of spray parcels
             // and store it
-            Cloud<passiveParticle> lagrangian(*meshPtr);
+            Cloud<passiveParticle> lagrangian(*meshPtr, cloud::defaultName);
 
             if (lagrangian.size() > nMaxParcels)
             {

applications/utilities/postProcessing/graphics/ensightFoamReader/USERD_set_filenames.H

@@ -151,7 +151,7 @@ int USERD_set_filenames
     {
         runTime.setTime(timeDirs[n+1], n+1);
 
-        Cloud<passiveParticle> lagrangian(*meshPtr);
+        Cloud<passiveParticle> lagrangian(*meshPtr, cloud::defaultName);
 
         n++;
         if (lagrangian.size())
@@ -177,7 +177,7 @@ int USERD_set_filenames
 
         delete sprayPtr;
 
-        sprayPtr = new Cloud<passiveParticle>(*meshPtr);
+        sprayPtr = new Cloud<passiveParticle>(*meshPtr, cloud::defaultName);
 
         IOobjectList objects(*meshPtr, runTime.timeName(), cloud::prefix);
 

applications/utilities/postProcessing/graphics/ensightFoamReader/USERD_set_time_set_and_step.H

@@ -49,9 +49,7 @@ void USERD_set_time_set_and_step
         if (Numparts_available > nPatches+1)
         {
             delete sprayPtr;
-            sprayPtr = new Cloud<passiveParticle>(*meshPtr);
+            sprayPtr = new Cloud<passiveParticle>(*meshPtr, cloud::defaultName);
         }
     }
 }
-
-

applications/utilities/preProcessing/mapFields/createTimes.H

@@ -1,11 +1,17 @@
     Info<< "\nCreate databases as time" << endl;
 
+    const auto caseDirOrig = getEnv("FOAM_CASE");
+    const auto caseNameOrig = getEnv("FOAM_CASE_NAME");
+    setEnv("FOAM_CASE", rootDirSource/caseDirSource, true);
+    setEnv("FOAM_CASE_NAME", caseDirSource, true);
     Time runTimeSource
     (
         Time::controlDictName,
         rootDirSource,
         caseDirSource
     );
+    setEnv("FOAM_CASE", caseDirOrig, true);
+    setEnv("FOAM_CASE_NAME", caseNameOrig, true);
 
     Time runTimeTarget
     (

applications/utilities/preProcessing/mapFields/mapFields.C

@@ -292,8 +292,8 @@ int main(int argc, char *argv[])
     fileName rootDirTarget(args.rootPath());
     fileName caseDirTarget(args.globalCaseName());
 
-    const fileName casePath = args[1];
+    fileName casePath = args[1];
-    const fileName rootDirSource = casePath.path();
+    const fileName rootDirSource = casePath.path().toAbsolute();
     const fileName caseDirSource = casePath.name();
 
     Info<< "Source: " << rootDirSource << " " << caseDirSource << endl;

bin/foamEtcFile

@@ -103,7 +103,7 @@ OpenFOAM-*)         # standard naming convention OpenFOAM-<VERSION>
     version="${projectDirName##OpenFOAM-}"
     ;;
 
-openfoam[0-9]* | openfoamdev)     # debian naming convention 'openfoam<VERSION>'
+openfoam[0-9]* | openfoam-dev)     # debian naming convention 'openfoam<VERSION>'
     versionNum="${projectDirName##openfoam}"
     case "$versionNum" in
     ??)         # convert 2 digit version number to decimal delineated

bin/tools/pre-commit-hook

@@ -3,8 +3,8 @@
 # =========                 |
 # \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
 #  \\    /   O peration     |
-#   \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+#   \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-#    \\/     M anipulation  |
+#    \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 #------------------------------------------------------------------------------
 # License
 #     This file is part of OpenFOAM.
@@ -48,11 +48,14 @@
 #     test the specified files/directories for standards conformance.
 #
 #------------------------------------------------------------------------------
+
 hookName="pre-commit"
+headerSeparator="-----------------------------------"
+
 die()
 {
     echo "$hookName hook failure" 1>&2
-    echo '-----------------------------------' 1>&2
+    echo $headerSeparator 1>&2
     echo '' 1>&2
     echo "$@" 1>&2
     echo '' 1>&2
@@ -105,7 +108,7 @@ dieOnBadFiles()
     if [ -n "$badFiles" ]
     then
         echo "$hookName hook failure" 1>&2
-        echo '-----------------------------------' 1>&2
+        echo $headerSeparator 1>&2
         echo "$@" 1>&2
         echo '' 1>&2
         echo "File(s):" 1>&2
@@ -267,6 +270,49 @@ checkLineLengthNonDirective()
 }
 
 
+#
+# check for non-standard code patterns
+#
+checkNonStandardCodePatterns()
+{
+    echo "$hookName: checking for non-standard code ..." 1>&2
+
+    scope=$(gitScope $@)
+
+    badFiles=$(
+    for f in $fileList
+    do
+        # limit to *.[CH] files
+        case "$f" in
+        (*.[CH])
+            # Directly report the incorrect markers
+            git grep -n --color \
+                -e '> >' -e '\bNULL\b' \
+                $scope"$f"
+        ;;
+        esac
+    done
+    )
+
+    dieOnBadFiles "$(cat<<MESSAGE
+Please revise the files reported below for the following non-standard code:
+
+  1. Spaced ending of multi-level template parameters are not allowed, such as:
+
+      List<List<scalar> >
+
+    which instead should be:
+
+      List<List<scalar>>
+
+  2. The use of the 'NULL' macro should be replaced by 'nullptr'
+
+$headerSeparator
+MESSAGE
+    )"
+}
+
+
 #
 # check that OpenFOAM Foundation copyright is current
 #
@@ -278,27 +324,27 @@ checkCopyright()
     badFiles=$(
     for f in $fileList
     do
-        startYear=`grep "Copyright.*OpenFOAM" $f | sed 's/[^0-9]*\([0-9]*\).*/\1/g'`
+        startYear=`grep "Copyright.*OpenCFD" $f | sed 's/[^0-9]*\([0-9]*\).*/\1/g'`
-        endYear=`grep "Copyright.*-.*OpenFOAM" $f | sed 's/[^-]*-\([0-9]*\).*/\1/g'`
+        endYear=`grep "Copyright.*-.*OpenCFD" $f | sed 's/[^-]*-\([0-9]*\).*/\1/g'`
         #echo "startYear=$startYear endYear=$endYear"
         if [ "$startYear" != "" ]
         then
             if [ "$endYear" != "" ]
             then
-                # Date is of type 2011-2012 OpenFOAM Foundation
+                # Date is of type 2011-2012 OpenCFD Ltd.
                 if [ "$year" != "$endYear" ]
                 then
                     echo "Updated copyright for: $f" 1>&2
                     echo "$f"
-                    sed -i -e "s/$startYear-$endYear OpenFOAM/$startYear-$year OpenFOAM/g" $f
+                    sed -i -e "s/$startYear-$endYear OpenCFD/$startYear-$year OpenCFD/g" $f
                 fi
             else
-                # Date is of type 2011 OpenFOAM Foundation
+                # Date is of type 2011 OpenCFD Ltd.
                 if [ "$year" != "$startYear" ]
                 then
                     echo "$f"
                     echo "Updated copyright for: $f" 1>&2
-                    sed -i -e "s/$startYear OpenFOAM/$startYear-$year OpenFOAM/g" $f
+                    sed -i -e "s/$startYear OpenCFD/$startYear-$year OpenCFD/g" $f
                 fi
             fi
         fi
@@ -322,6 +368,9 @@ checkIllegalCode
 # ensure code conforms to 80 columns max
 checkLineLengthNonDirective
 
+# check for non-standard code patterns
+checkNonStandardCodePatterns
+
 checkCopyright
 
 exit 0

etc/codeTemplates/dynamicCode/codedFvOptionTemplate.C

@@ -104,7 +104,7 @@ ${typeName}FvOption${SourceType}
     const fvMesh& mesh
 )
 :
-    option(name, modelType, dict, mesh)
+    cellSetOption(name, modelType, dict, mesh)
 {
     if (${verbose:-false})
     {

etc/codeTemplates/dynamicCode/codedFvOptionTemplate.H

@@ -105,7 +105,7 @@ SourceFiles
 #ifndef codedFvOptionTemplate_H
 #define codedFvOptionTemplate_H
 
-#include "fvOption.H"
+#include "cellSetOption.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -121,7 +121,7 @@ namespace fv
 
 class ${typeName}FvOption${SourceType}
 :
-    public option
+    public cellSetOption
 {
 public:
 

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

@@ -174,7 +174,8 @@ bool Foam::functionObjectList::readFunctionObject
 (
     const string& funcNameArgs,
     dictionary& functionsDict,
-    HashSet<word>& requiredFields
+    HashSet<word>& requiredFields,
+    const word& region
 )
 {
     // Parse the optional functionObject arguments:
@@ -314,6 +315,12 @@ bool Foam::functionObjectList::readFunctionObject
         funcDict.set(entry::New(entryStream).ptr());
     }
 
+    // Insert the region name if specified
+    if (region != word::null)
+    {
+        funcDict.set("region", region);
+    }
+
     // Merge this functionObject dictionary into functionsDict
     dictionary funcArgsDict;
     funcArgsDict.add(string::validate<word>(funcNameArgs), funcDict);
@@ -377,6 +384,14 @@ Foam::autoPtr<Foam::functionObjectList> Foam::functionObjectList::New
 
     dictionary& functionsDict = controlDict.subDict("functions");
 
+    word region = word::null;
+
+    // Set the region name if specified
+    if (args.optionFound("region"))
+    {
+        region = args["region"];
+    }
+
     if
     (
         args.optionFound("dict")
@@ -402,7 +417,13 @@ Foam::autoPtr<Foam::functionObjectList> Foam::functionObjectList::New
 
         if (args.optionFound("func"))
         {
-            readFunctionObject(args["func"], functionsDict, requiredFields);
+            readFunctionObject
+            (
+                args["func"],
+                functionsDict,
+                requiredFields,
+                region
+            );
         }
 
         if (args.optionFound("funcs"))
@@ -411,7 +432,13 @@ Foam::autoPtr<Foam::functionObjectList> Foam::functionObjectList::New
 
             forAll(funcs, i)
             {
-                readFunctionObject(funcs[i], functionsDict, requiredFields);
+                readFunctionObject
+                (
+                    funcs[i],
+                    functionsDict,
+                    requiredFields,
+                    region
+                );
             }
         }
 
@@ -701,12 +728,14 @@ bool Foam::functionObjectList::read()
                             "functionObject::" + objPtr->name() + "::read"
                         );
 
-                        ok = objPtr->read(dict) && ok;
+                        enabled = objPtr->read(dict);
+                        ok = enabled && ok;
                     }
                 }
-                else
+
+                if (!enabled)
                 {
-                    // Delete the disabled functionObject
+                    // Delete the disabled/invalid(read) functionObject
                     delete objPtr;
                     objPtr = nullptr;
                     continue;

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

@@ -234,7 +234,8 @@ public:
         (
             const string& funcNameArgs0,
             dictionary& functionsDict,
-            HashSet<word>& requiredFields
+            HashSet<word>& requiredFields,
+            const word& region = word::null
         );
 
         //- Read and set the function objects if their data have changed

src/OpenFOAM/fields/cloud/cloud.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -30,10 +30,10 @@ License
 
 namespace Foam
 {
-defineTypeNameAndDebug(cloud, 0);
+    defineTypeNameAndDebug(cloud, 0);
 
-const word cloud::prefix("lagrangian");
+    const word cloud::prefix("lagrangian");
-word cloud::defaultName("defaultCloud");
+    word cloud::defaultName("defaultCloud");
 }
 
 
@@ -45,7 +45,7 @@ Foam::cloud::cloud(const objectRegistry& obr, const word& cloudName)
     (
         IOobject
         (
-            (cloudName.size() ? cloudName : defaultName),
+            cloudName,
             obr.time().timeName(),
             prefix,
             obr,

src/OpenFOAM/fields/cloud/cloud.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -78,7 +78,7 @@ public:
     // Constructors
 
         //- Construct for the given objectRegistry and named cloud instance
-        cloud(const objectRegistry&, const word& cloudName = "");
+        cloud(const objectRegistry&, const word& cloudName = defaultName);
 
 
     //- Destructor

src/OpenFOAM/interpolations/patchToPatchInterpolation/PatchToPatchInterpolate.C

@@ -54,7 +54,7 @@ PatchToPatchInterpolation<FromPatch, ToPatch>::pointInterpolate
         new Field<Type>(toPatch_.nPoints(), Zero)
     );
 
-    Field<Type>& result = tresult();
+    Field<Type>& result = tresult.ref();
 
     const List<typename FromPatch::FaceType>& fromPatchLocalFaces =
         fromPatch_.localFaces();
@@ -118,7 +118,7 @@ PatchToPatchInterpolation<FromPatch, ToPatch>::faceInterpolate
         new Field<Type>(toPatch_.size(), Zero)
     );
 
-    Field<Type>& result = tresult();
+    Field<Type>& result = tresult.ref();
 
     const labelListList& fromPatchFaceFaces = fromPatch_.faceFaces();
 

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

@@ -211,9 +211,9 @@ bool Foam::GAMGAgglomeration::continueAgglomerating
 ) const
 {
     const label nTotalCoarseCells = returnReduce(nCoarseCells, sumOp<label>());
-    if (nTotalCoarseCells >= Pstream::nProcs()*nCellsInCoarsestLevel_)
+    if (nTotalCoarseCells < Pstream::nProcs()*nCellsInCoarsestLevel_)
     {
-        return true;
+        return false;
     }
     else
     {

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.H

@@ -42,7 +42,7 @@ Description
       - Coarse matrix scaling: performed by correction scaling, using steepest
         descent optimisation.
       - Type of cycle: V-cycle with optional pre-smoothing.
-      - Coarsest-level matrix solved using PCG or PBiCG.
+      - Coarsest-level matrix solved using PCG or PBiCGStab.
 
 SourceFiles
     GAMGSolver.C
@@ -303,9 +303,9 @@ class GAMGSolver
             const scalar relTol
         ) const;
 
-        //- Create and return the dictionary to specify the PBiCG solver
+        //- Create and return the dictionary to specify the PBiCGStab solver
         //  to solve the coarsest level
-        dictionary PBiCGsolverDict
+        dictionary PBiCGStabSolverDict
         (
             const scalar tol,
             const scalar relTol

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C

@@ -25,7 +25,7 @@ License
 
 #include "GAMGSolver.H"
 #include "PCG.H"
-#include "PBiCG.H"
+#include "PBiCGStab.H"
 #include "SubField.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -533,13 +533,13 @@ Foam::dictionary Foam::GAMGSolver::PCGsolverDict
 }
 
 
-Foam::dictionary Foam::GAMGSolver::PBiCGsolverDict
+Foam::dictionary Foam::GAMGSolver::PBiCGStabSolverDict
 (
     const scalar tol,
     const scalar relTol
 ) const
 {
-    dictionary dict(IStringStream("solver PBiCG; preconditioner DILU;")());
+    dictionary dict(IStringStream("solver PBiCGStab; preconditioner DILU;")());
     dict.add("tolerance", tol);
     dict.add("relTol", relTol);
 
@@ -612,14 +612,14 @@ void Foam::GAMGSolver::solveCoarsestLevel
     //
     //            if (allMatrix.asymmetric())
     //            {
-    //                coarseSolverPerf = PBiCG
+    //                coarseSolverPerf = PBiCGStab
     //                (
     //                    "coarsestLevelCorr",
     //                    allMatrix,
     //                    procInterfaceLevelsBouCoeffs_[coarsestLevel],
     //                    procInterfaceLevelsIntCoeffs_[coarsestLevel],
     //                    procInterfaceLevels_[coarsestLevel],
-    //                    PBiCGsolverDict(tolerance_, relTol_)
+    //                    PBiCGStabSolverDict(tolerance_, relTol_)
     //                ).solve
     //                (
     //                    coarsestCorrField,
@@ -673,14 +673,14 @@ void Foam::GAMGSolver::solveCoarsestLevel
 
         if (matrixLevels_[coarsestLevel].asymmetric())
         {
-            coarseSolverPerf = PBiCG
+            coarseSolverPerf = PBiCGStab
             (
                 "coarsestLevelCorr",
                 matrixLevels_[coarsestLevel],
                 interfaceLevelsBouCoeffs_[coarsestLevel],
                 interfaceLevelsIntCoeffs_[coarsestLevel],
                 interfaceLevels_[coarsestLevel],
-                PBiCGsolverDict(tolerance_, relTol_)
+                PBiCGStabSolverDict(tolerance_, relTol_)
             ).solve
             (
                 coarsestCorrField,

src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C

@@ -24,19 +24,14 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "globalMeshData.H"
-#include "Time.H"
 #include "Pstream.H"
 #include "PstreamCombineReduceOps.H"
 #include "processorPolyPatch.H"
-#include "demandDrivenData.H"
 #include "globalPoints.H"
 #include "polyMesh.H"
 #include "mapDistribute.H"
 #include "labelIOList.H"
-#include "PackedList.H"
 #include "mergePoints.H"
-#include "matchPoints.H"
-#include "OFstream.H"
 #include "globalIndexAndTransform.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -577,7 +572,7 @@ void Foam::globalMeshData::calcPointConnectivity
     labelPairList myData(globalPointSlavesMap().constructSize());
     forAll(slaves, pointi)
     {
-        myData[pointi] = globalIndexAndTransform::encode
+        myData[pointi] = transforms.encode
         (
             Pstream::myProcNo(),
             pointi,
@@ -624,9 +619,9 @@ void Foam::globalMeshData::calcPointConnectivity
                 );
                 // Add transform to connectivity
                 const labelPair& n = myData[pTransformSlaves[i]];
-                label proci = globalIndexAndTransform::processor(n);
+                label proci = transforms.processor(n);
-                label index = globalIndexAndTransform::index(n);
+                label index = transforms.index(n);
-                pConnectivity[connI++] = globalIndexAndTransform::encode
+                pConnectivity[connI++] = transforms.encode
                 (
                     proci,
                     index,
@@ -678,6 +673,8 @@ void Foam::globalMeshData::calcGlobalPointEdges
     const globalIndex& globalEdgeNumbers = globalEdgeNumbering();
     const labelListList& slaves = globalPointSlaves();
     const labelListList& transformedSlaves = globalPointTransformedSlaves();
+    const globalIndexAndTransform& transforms = globalTransforms();
+
 
     // Create local version
     globalPointEdges.setSize(globalPointSlavesMap().constructSize());
@@ -697,11 +694,11 @@ void Foam::globalMeshData::calcGlobalPointEdges
         forAll(pEdges, i)
         {
             label otherPointi = edges[pEdges[i]].otherVertex(pointi);
-            globalPPoints[i] = globalIndexAndTransform::encode
+            globalPPoints[i] = transforms.encode
             (
                 Pstream::myProcNo(),
                 otherPointi,
-                globalTransforms().nullTransformIndex()
+                transforms.nullTransformIndex()
             );
         }
     }
@@ -790,9 +787,9 @@ void Foam::globalMeshData::calcGlobalPointEdges
                 {
                     // Add transform to connectivity
                     const labelPair& n = otherData[j];
-                    label proci = globalIndexAndTransform::processor(n);
+                    label proci = transforms.processor(n);
-                    label index = globalIndexAndTransform::index(n);
+                    label index = transforms.index(n);
-                    globalPPoints[sz++] = globalIndexAndTransform::encode
+                    globalPPoints[sz++] = transforms.encode
                     (
                         proci,
                         index,
@@ -834,16 +831,18 @@ Foam::label Foam::globalMeshData::findTransform
     const label localPoint
 ) const
 {
-    const label remoteProci = globalIndexAndTransform::processor(remotePoint);
+    const globalIndexAndTransform& transforms = globalTransforms();
-    const label remoteIndex = globalIndexAndTransform::index(remotePoint);
+
+    const label remoteProci = transforms.processor(remotePoint);
+    const label remoteIndex = transforms.index(remotePoint);
 
     label remoteTransformI = -1;
     label localTransformI = -1;
     forAll(info, i)
     {
-        label proci = globalIndexAndTransform::processor(info[i]);
+        label proci = transforms.processor(info[i]);
-        label pointi = globalIndexAndTransform::index(info[i]);
+        label pointi = transforms.index(info[i]);
-        label transformI = globalIndexAndTransform::transformIndex(info[i]);
+        label transformI = transforms.transformIndex(info[i]);
 
         if (proci == Pstream::myProcNo() && pointi == localPoint)
         {
@@ -875,7 +874,7 @@ Foam::label Foam::globalMeshData::findTransform
             << abort(FatalError);
     }
 
-    return globalTransforms().subtractTransformIndex
+    return transforms.subtractTransformIndex
     (
         remoteTransformI,
         localTransformI
@@ -893,6 +892,7 @@ void Foam::globalMeshData::calcGlobalEdgeSlaves() const
 
     const edgeList& edges = coupledPatch().edges();
     const globalIndex& globalEdgeNumbers = globalEdgeNumbering();
+    const globalIndexAndTransform& transforms = globalTransforms();
 
 
     // The whole problem with deducting edge-connectivity from
@@ -941,11 +941,11 @@ void Foam::globalMeshData::calcGlobalEdgeSlaves() const
         // Append myself.
         eEdges.append
         (
-            globalIndexAndTransform::encode
+            transforms.encode
             (
                 Pstream::myProcNo(),
                 edgeI,
-                globalTransforms().nullTransformIndex()
+                transforms.nullTransformIndex()
             )
         );
 
@@ -986,7 +986,7 @@ void Foam::globalMeshData::calcGlobalEdgeSlaves() const
                         label proci = globalEdgeNumbers.whichProcID(pEdges0[i]);
                         eEdges.append
                         (
-                            globalIndexAndTransform::encode
+                            transforms.encode
                             (
                                 proci,
                                 globalEdgeNumbers.toLocal(proci, pEdges0[i]),
@@ -999,7 +999,11 @@ void Foam::globalMeshData::calcGlobalEdgeSlaves() const
         }
 
         allEdgeConnectivity[edgeI].transfer(eEdges);
-        sort(allEdgeConnectivity[edgeI], globalIndexAndTransform::less());
+        sort
+        (
+            allEdgeConnectivity[edgeI],
+            globalIndexAndTransform::less(transforms)
+        );
     }
 
     // We now have - in allEdgeConnectivity - a list of edges which are shared
@@ -1020,10 +1024,10 @@ void Foam::globalMeshData::calcGlobalEdgeSlaves() const
             if
             (
                 (
-                    globalIndexAndTransform::processor(masterInfo)
+                    transforms.processor(masterInfo)
                  == Pstream::myProcNo()
                 )
-             && (globalIndexAndTransform::index(masterInfo) == edgeI)
+             && (transforms.index(masterInfo) == edgeI)
             )
             {
                 // Sort into transformed and untransformed
@@ -1039,14 +1043,14 @@ void Foam::globalMeshData::calcGlobalEdgeSlaves() const
                 for (label i = 1; i < edgeInfo.size(); i++)
                 {
                     const labelPair& info = edgeInfo[i];
-                    label proci = globalIndexAndTransform::processor(info);
+                    label proci = transforms.processor(info);
-                    label index = globalIndexAndTransform::index(info);
+                    label index = transforms.index(info);
-                    label transform = globalIndexAndTransform::transformIndex
+                    label transform = transforms.transformIndex
                     (
                         info
                     );
 
-                    if (transform == globalTransforms().nullTransformIndex())
+                    if (transform == transforms.nullTransformIndex())
                     {
                         eEdges[nonTransformI++] = globalEdgeNumbers.toGlobal
                         (
@@ -1078,7 +1082,7 @@ void Foam::globalMeshData::calcGlobalEdgeSlaves() const
             globalEdgeNumbers,
             globalEdgeSlaves,
 
-            globalTransforms(),
+            transforms,
             transformedEdges,
             globalEdgeTransformedSlavesPtr_(),
 
@@ -1351,6 +1355,7 @@ void Foam::globalMeshData::calcGlobalPointBoundaryFaces() const
     const labelListList& pointSlaves = globalPointSlaves();
     const labelListList& pointTransformSlaves =
         globalPointTransformedSlaves();
+    const globalIndexAndTransform& transforms = globalTransforms();
 
 
     // Any faces coming in through transformation
@@ -1432,7 +1437,7 @@ void Foam::globalMeshData::calcGlobalPointBoundaryFaces() const
                         label proci = globalIndices.whichProcID(slave);
                         label facei = globalIndices.toLocal(proci, slave);
 
-                        myBFaces[n++] = globalIndexAndTransform::encode
+                        myBFaces[n++] = transforms.encode
                         (
                             proci,
                             facei,
@@ -1466,7 +1471,7 @@ void Foam::globalMeshData::calcGlobalPointBoundaryFaces() const
             globalIndices,
             globalPointBoundaryFaces,
 
-            globalTransforms(),
+            transforms,
             transformedFaces,
             globalPointTransformedBoundaryFacesPtr_(),
 
@@ -1581,6 +1586,7 @@ void Foam::globalMeshData::calcGlobalPointBoundaryCells() const
     const labelListList& pointSlaves = globalPointSlaves();
     const labelListList& pointTransformSlaves =
         globalPointTransformedSlaves();
+    const globalIndexAndTransform& transforms = globalTransforms();
 
     List<labelPairList> transformedCells(pointSlaves.size());
 
@@ -1660,7 +1666,7 @@ void Foam::globalMeshData::calcGlobalPointBoundaryCells() const
                     {
                         label proci = globalIndices.whichProcID(slave);
                         label celli = globalIndices.toLocal(proci, slave);
-                        myBCells[n++] = globalIndexAndTransform::encode
+                        myBCells[n++] = transforms.encode
                         (
                             proci,
                             celli,
@@ -1693,7 +1699,7 @@ void Foam::globalMeshData::calcGlobalPointBoundaryCells() const
             globalIndices,
             globalPointBoundaryCells,
 
-            globalTransforms(),
+            transforms,
             transformedCells,
             globalPointTransformedBoundaryCellsPtr_(),
 

src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H

@@ -72,7 +72,6 @@ See also
     mapDistribute
     globalIndexAndTransform
 
-
 SourceFiles
     globalMeshData.C
     globalMeshDataTemplates.C

src/OpenFOAM/meshes/polyMesh/globalMeshData/globalPoints.C

@@ -64,15 +64,15 @@ Foam::label Foam::globalPoints::findSamePoint
     const labelPair& info
 ) const
 {
-    const label proci = globalIndexAndTransform::processor(info);
+    const label proci = globalTransforms_.processor(info);
-    const label index = globalIndexAndTransform::index(info);
+    const label index = globalTransforms_.index(info);
 
     forAll(allInfo, i)
     {
         if
         (
-            globalIndexAndTransform::processor(allInfo[i]) == proci
+            globalTransforms_.processor(allInfo[i]) == proci
-         && globalIndexAndTransform::index(allInfo[i]) == index
+         && globalTransforms_.index(allInfo[i]) == index
         )
         {
             return i;
@@ -98,21 +98,21 @@ Foam::labelPairList Foam::globalPoints::addSendTransform
     forAll(info, i)
     {
         //Pout<< "    adding send transform to" << nl
-        //    << "    proc:" << globalIndexAndTransform::processor(info[i])
+        //    << "    proc:" << globalTransforms_.processor(info[i])
         //    << nl
-        //    << "    index:" << globalIndexAndTransform::index(info[i]) << nl
+        //    << "    index:" << globalTransforms_.index(info[i]) << nl
         //    << "    trafo:"
         //    <<  globalTransforms_.decodeTransformIndex
-        //        (globalIndexAndTransform::transformIndex(info[i]))
+        //        (globalTransforms_.transformIndex(info[i]))
         //    << endl;
 
-        sendInfo[i] = globalIndexAndTransform::encode
+        sendInfo[i] = globalTransforms_.encode
         (
-            globalIndexAndTransform::processor(info[i]),
+            globalTransforms_.processor(info[i]),
-            globalIndexAndTransform::index(info[i]),
+            globalTransforms_.index(info[i]),
             globalTransforms_.addToTransformIndex
             (
-                globalIndexAndTransform::transformIndex(info[i]),
+                globalTransforms_.transformIndex(info[i]),
                 patchi,
                 true,           // patchi is sending side
                 tol             // tolerance for comparison
@@ -204,11 +204,11 @@ bool Foam::globalPoints::mergeInfo
             }
             else
             {
-                label myTransform = globalIndexAndTransform::transformIndex
+                label myTransform = globalTransforms_.transformIndex
                 (
                     myInfo[index]
                 );
-                label nbrTransform = globalIndexAndTransform::transformIndex
+                label nbrTransform = globalTransforms_.transformIndex
                 (
                     nbrInfo[i]
                 );
@@ -294,7 +294,7 @@ bool Foam::globalPoints::mergeInfo
         labelPairList knownInfo
         (
             1,
-            globalIndexAndTransform::encode
+            globalTransforms_.encode
             (
                 Pstream::myProcNo(),
                 localPointi,
@@ -356,9 +356,9 @@ void Foam::globalPoints::printProcPoint
     const labelPair& pointInfo
 ) const
 {
-    label proci = globalIndexAndTransform::processor(pointInfo);
+    label proci = globalTransforms_.processor(pointInfo);
-    label index = globalIndexAndTransform::index(pointInfo);
+    label index = globalTransforms_.index(pointInfo);
-    label trafoI = globalIndexAndTransform::transformIndex(pointInfo);
+    label trafoI = globalTransforms_.transformIndex(pointInfo);
 
     Pout<< "    proc:" << proci;
     Pout<< " localpoint:";
@@ -421,7 +421,7 @@ void Foam::globalPoints::initOwnPoints
                     labelPairList knownInfo
                     (
                         1,
-                        globalIndexAndTransform::encode
+                        globalTransforms_.encode
                         (
                             Pstream::myProcNo(),
                             localPointi,
@@ -457,7 +457,7 @@ void Foam::globalPoints::initOwnPoints
                     labelPairList knownInfo
                     (
                         1,
-                        globalIndexAndTransform::encode
+                        globalTransforms_.encode
                         (
                             Pstream::myProcNo(),
                             localPointi,
@@ -750,8 +750,8 @@ void Foam::globalPoints::remove
             // is in it. This would be an ordinary connection and can be
             // handled by normal face-face connectivity.
 
-            label proc0 = globalIndexAndTransform::processor(pointInfo[0]);
+            label proc0 = globalTransforms_.processor(pointInfo[0]);
-            label proc1 = globalIndexAndTransform::processor(pointInfo[1]);
+            label proc1 = globalTransforms_.processor(pointInfo[1]);
 
             if
             (
@@ -759,14 +759,14 @@ void Foam::globalPoints::remove
                     proc0 == Pstream::myProcNo()
                  && directNeighbours.found
                     (
-                        globalIndexAndTransform::index(pointInfo[0])
+                        globalTransforms_.index(pointInfo[0])
                     )
                 )
              || (
                     proc1 == Pstream::myProcNo()
                  && directNeighbours.found
                     (
-                        globalIndexAndTransform::index(pointInfo[1])
+                        globalTransforms_.index(pointInfo[1])
                     )
                 )
             )
@@ -776,14 +776,14 @@ void Foam::globalPoints::remove
                 {
                     //Pout<< "Removing direct neighbour:"
                     //    << mesh_.points()
-                    //       [globalIndexAndTransform::index(pointInfo[0])]
+                    //       [globalTransforms_.index(pointInfo[0])]
                     //    << endl;
                 }
                 else if (proc1 == Pstream::myProcNo())
                 {
                     //Pout<< "Removing direct neighbour:"
                     //    << mesh_.points()
-                    //       [globalIndexAndTransform::index(pointInfo[1])]
+                    //       [globalTransforms_.index(pointInfo[1])]
                     //    << endl;
                 }
             }
@@ -812,11 +812,11 @@ void Foam::globalPoints::remove
             // So this meshPoint will have info of size one only.
             if
             (
-                globalIndexAndTransform::processor(pointInfo[0])
+                globalTransforms_.processor(pointInfo[0])
              != Pstream::myProcNo()
              || !directNeighbours.found
                 (
-                    globalIndexAndTransform::index(pointInfo[0])
+                    globalTransforms_.index(pointInfo[0])
                 )
             )
             {
@@ -995,7 +995,7 @@ void Foam::globalPoints::calculateSharedPoints
     forAllConstIter(Map<label>, meshToProcPoint_, iter)
     {
         labelPairList& pointInfo = procPoints_[iter()];
-        sort(pointInfo, globalIndexAndTransform::less());
+        sort(pointInfo, globalIndexAndTransform::less(globalTransforms_));
     }
 
 
@@ -1017,10 +1017,10 @@ void Foam::globalPoints::calculateSharedPoints
             if
             (
                 (
-                    globalIndexAndTransform::processor(masterInfo)
+                    globalTransforms_.processor(masterInfo)
                  == Pstream::myProcNo()
                 )
-             && (globalIndexAndTransform::index(masterInfo) == iter.key())
+             && (globalTransforms_.index(masterInfo) == iter.key())
             )
             {
                 labelList& pPoints = pointPoints_[iter.key()];
@@ -1035,9 +1035,9 @@ void Foam::globalPoints::calculateSharedPoints
                 for (label i = 1; i < pointInfo.size(); i++)
                 {
                     const labelPair& info = pointInfo[i];
-                    label proci = globalIndexAndTransform::processor(info);
+                    label proci = globalTransforms_.processor(info);
-                    label index = globalIndexAndTransform::index(info);
+                    label index = globalTransforms_.index(info);
-                    label transform = globalIndexAndTransform::transformIndex
+                    label transform = globalTransforms_.transformIndex
                     (
                         info
                     );

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

@@ -271,10 +271,10 @@ Foam::mapDistribute::mapDistribute
     forAll(transformedElements, i)
     {
         labelPair elem = transformedElements[i];
-        label proci = globalIndexAndTransform::processor(elem);
+        label proci = globalTransforms.processor(elem);
         if (proci != Pstream::myProcNo())
         {
-            label index = globalIndexAndTransform::index(elem);
+            label index = globalTransforms.index(elem);
             label nCompact = compactMap[proci].size();
             compactMap[proci].insert(index, nCompact);
         }
@@ -302,7 +302,7 @@ Foam::mapDistribute::mapDistribute
     forAll(transformedElements, i)
     {
         labelPair elem = transformedElements[i];
-        label trafoI = globalIndexAndTransform::transformIndex(elem);
+        label trafoI = globalTransforms.transformIndex(elem);
         nPerTransform[trafoI]++;
     }
     // Offset per transformIndex
@@ -322,9 +322,9 @@ Foam::mapDistribute::mapDistribute
     forAll(transformedElements, i)
     {
         labelPair elem = transformedElements[i];
-        label proci = globalIndexAndTransform::processor(elem);
+        label proci = globalTransforms.processor(elem);
-        label index = globalIndexAndTransform::index(elem);
+        label index = globalTransforms.index(elem);
-        label trafoI = globalIndexAndTransform::transformIndex(elem);
+        label trafoI = globalTransforms.transformIndex(elem);
 
         // Get compact index for untransformed element
         label rawElemI =
@@ -379,10 +379,10 @@ Foam::mapDistribute::mapDistribute
 
         forAll(elems, i)
         {
-            label proci = globalIndexAndTransform::processor(elems[i]);
+            label proci = globalTransforms.processor(elems[i]);
             if (proci != Pstream::myProcNo())
             {
-                label index = globalIndexAndTransform::index(elems[i]);
+                label index = globalTransforms.index(elems[i]);
                 label nCompact = compactMap[proci].size();
                 compactMap[proci].insert(index, nCompact);
             }
@@ -414,7 +414,7 @@ Foam::mapDistribute::mapDistribute
 
         forAll(elems, i)
         {
-            label trafoI = globalIndexAndTransform::transformIndex(elems[i]);
+            label trafoI = globalTransforms.transformIndex(elems[i]);
             nPerTransform[trafoI]++;
         }
     }
@@ -439,9 +439,9 @@ Foam::mapDistribute::mapDistribute
 
         forAll(elems, i)
         {
-            label proci = globalIndexAndTransform::processor(elems[i]);
+            label proci = globalTransforms.processor(elems[i]);
-            label index = globalIndexAndTransform::index(elems[i]);
+            label index = globalTransforms.index(elems[i]);
-            label trafoI = globalIndexAndTransform::transformIndex(elems[i]);
+            label trafoI = globalTransforms.transformIndex(elems[i]);
 
             // Get compact index for untransformed element
             label rawElemI =

src/OpenFOAM/primitives/globalIndexAndTransform/globalIndexAndTransform.C

@@ -25,13 +25,14 @@ License
 
 #include "globalIndexAndTransform.H"
 #include "cyclicPolyPatch.H"
+#include "DynamicField.H"
+#include "globalMeshData.H"
 
 // * * * * * * * * * * * * Private Static Data Members * * * * * * * * * * * //
 
 namespace Foam
 {
-defineTypeNameAndDebug(globalIndexAndTransform, 0);
+    defineTypeNameAndDebug(globalIndexAndTransform, 0);
-const label globalIndexAndTransform::base_ = 32;
 }
 
 
@@ -127,10 +128,8 @@ void Foam::globalIndexAndTransform::determineTransforms()
 {
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
-    transforms_ = List<vectorTensorTransform>(6);
+    DynamicList<vectorTensorTransform> localTransforms;
-    scalarField maxTol(6);
+    DynamicField<scalar> localTols;
-
-    label nextTrans = 0;
 
     label dummyMatch = -1;
 
@@ -170,7 +169,7 @@ void Foam::globalIndexAndTransform::determineTransforms()
                         (
                             matchTransform
                             (
-                                transforms_,
+                                localTransforms,
                                 dummyMatch,
                                 transform,
                                 cpp.matchTolerance(),
@@ -178,15 +177,8 @@ void Foam::globalIndexAndTransform::determineTransforms()
                             ) == 0
                         )
                         {
-                            if (nextTrans == 6)
+                            localTransforms.append(transform);
-                            {
+                            localTols.append(cpp.matchTolerance());
-                                FatalErrorInFunction
-                                    << "More than six unsigned transforms"
-                                    << " detected:" << nl << transforms_
-                                    << exit(FatalError);
-                            }
-                            transforms_[nextTrans] = transform;
-                            maxTol[nextTrans++] = cpp.matchTolerance();
                         }
                     }
                 }
@@ -207,7 +199,7 @@ void Foam::globalIndexAndTransform::determineTransforms()
                         (
                             matchTransform
                             (
-                                transforms_,
+                                localTransforms,
                                 dummyMatch,
                                 transform,
                                 cpp.matchTolerance(),
@@ -215,15 +207,8 @@ void Foam::globalIndexAndTransform::determineTransforms()
                             ) == 0
                         )
                         {
-                            if (nextTrans == 6)
+                            localTransforms.append(transform);
-                            {
+                            localTols.append(cpp.matchTolerance());
-                                FatalErrorInFunction
-                                    << "More than six unsigned transforms"
-                                    << " detected:" << nl << transforms_
-                                    << exit(FatalError);
-                            }
-                            transforms_[nextTrans] = transform;
-                            maxTol[nextTrans++] = cpp.matchTolerance();
                         }
                     }
                 }
@@ -233,21 +218,18 @@ void Foam::globalIndexAndTransform::determineTransforms()
 
 
     // Collect transforms on master
-
     List<List<vectorTensorTransform>> allTransforms(Pstream::nProcs());
-    allTransforms[Pstream::myProcNo()] = transforms_;
+    allTransforms[Pstream::myProcNo()] = localTransforms;
     Pstream::gatherList(allTransforms);
 
     // Collect matching tolerance on master
     List<scalarField> allTols(Pstream::nProcs());
-    allTols[Pstream::myProcNo()] = maxTol;
+    allTols[Pstream::myProcNo()] = localTols;
     Pstream::gatherList(allTols);
 
     if (Pstream::master())
     {
-        transforms_ = List<vectorTensorTransform>(3);
+        localTransforms.clear();
-
-        label nextTrans = 0;
 
         forAll(allTransforms, proci)
         {
@@ -264,45 +246,23 @@ void Foam::globalIndexAndTransform::determineTransforms()
                     (
                         matchTransform
                         (
-                            transforms_,
+                            localTransforms,
                             dummyMatch,
                             transform,
                             allTols[proci][pSVI],
                             true
-                        ) ==  0
+                        ) == 0
                     )
                     {
-                        transforms_[nextTrans++] = transform;
+                        localTransforms.append(transform);
-                    }
-
-                    if (nextTrans > 3)
-                    {
-                        FatalErrorInFunction
-                            << "More than three independent basic "
-                            << "transforms detected:" << nl
-                            << allTransforms
-                            << transforms_
-                            << exit(FatalError);
                     }
                 }
             }
         }
-
-        transforms_.setSize(nextTrans);
     }
 
+    transforms_.transfer(localTransforms);
     Pstream::scatter(transforms_);
-
-    if (transforms_.size() > 3)
-    {
-        WarningInFunction
-            << "More than three independent basic "
-            << "transforms detected:" << nl
-            << transforms_ << nl
-            << "This is not a space filling tiling and will probably"
-            << " give problems for e.g. lagrangian tracking or interpolation"
-            << endl;
-    }
 }
 
 
@@ -351,16 +311,12 @@ void Foam::globalIndexAndTransform::determinePatchTransformSign()
 {
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
-    patchTransformSign_.setSize(patches.size(), Pair<label>(-1, 0));
+    patchTransformSign_.setSize(patches.size(), labelPair(-1, 0));
-
-    label matchTransI = -1;
 
     forAll(patches, patchi)
     {
         const polyPatch& pp = patches[patchi];
 
-        // Pout<< nl << patchi << " " << pp.name() << endl;
-
         // Note: special check for unordered cyclics. These are in fact
         // transform bcs and should probably be split off.
         if
@@ -375,15 +331,12 @@ void Foam::globalIndexAndTransform::determinePatchTransformSign()
             )
         )
         {
-            const coupledPolyPatch& cpp =
+            const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>(pp);
-            refCast<const coupledPolyPatch>(pp);
 
             if (cpp.separated())
             {
                 const vectorField& sepVecs = cpp.separation();
 
-                // Pout<< "sepVecs " << sepVecs << endl;
-
                 // This loop is implicitly expecting only a single
                 // value for separation()
                 forAll(sepVecs, sVI)
@@ -394,6 +347,7 @@ void Foam::globalIndexAndTransform::determinePatchTransformSign()
                     {
                         vectorTensorTransform t(sepVec);
 
+                        label matchTransI;
                         label sign = matchTransform
                         (
                             transforms_,
@@ -402,22 +356,8 @@ void Foam::globalIndexAndTransform::determinePatchTransformSign()
                             cpp.matchTolerance(),
                             true
                         );
-
-                        // Pout<< sign << " " << matchTransI << endl;
-
-                        // List<label> permutation(transforms_.size(), 0);
-
-                        // permutation[matchTransI] = sign;
-
-                        // Pout<< encodeTransformIndex(permutation) << nl
-                        //     << transformPermutations_
-                        //        [
-                        //            encodeTransformIndex(permutation)
-                        //        ]
-                        //     << endl;
-
                         patchTransformSign_[patchi] =
-                            Pair<label>(matchTransI, sign);
+                            labelPair(matchTransI, sign);
                     }
                 }
 
@@ -426,8 +366,6 @@ void Foam::globalIndexAndTransform::determinePatchTransformSign()
             {
                 const tensorField& transTensors = cpp.reverseT();
 
-                // Pout<< "transTensors " << transTensors << endl;
-
                 // This loop is implicitly expecting only a single
                 // value for reverseT()
                 forAll(transTensors, tTI)
@@ -438,6 +376,7 @@ void Foam::globalIndexAndTransform::determinePatchTransformSign()
                     {
                         vectorTensorTransform t(transT);
 
+                        label matchTransI;
                         label sign = matchTransform
                         (
                             transforms_,
@@ -447,37 +386,65 @@ void Foam::globalIndexAndTransform::determinePatchTransformSign()
                             true
                         );
 
-                        // Pout<< sign << " " << matchTransI << endl;
-
-                        // List<label> permutation(transforms_.size(), 0);
-
-                        // permutation[matchTransI] = sign;
-
-                        // Pout<< encodeTransformIndex(permutation) << nl
-                        //     << transformPermutations_
-                        //        [
-                        //            encodeTransformIndex(permutation)
-                        //        ]
-                        //     << endl;
-
                         patchTransformSign_[patchi] =
-                            Pair<label>(matchTransI, sign);
+                            labelPair(matchTransI, sign);
                     }
                 }
             }
         }
     }
+}
 
-    // Pout<< patchTransformSign_ << endl;
+
+bool Foam::globalIndexAndTransform::uniqueTransform
+(
+    const point& pt,
+    labelPairList& trafos,
+    const label patchi,
+    const labelPair& patchTrafo
+) const
+{
+    if (findIndex(trafos, patchTrafo) == -1)
+    {
+        // New transform. Check if already have 3
+        if (trafos.size() == 3)
+        {
+            if (patchi > -1)
+            {
+                WarningInFunction
+                    << "Point " << pt
+                    << " is on patch " << mesh_.boundaryMesh()[patchi].name();
+            }
+            else
+            {
+                WarningInFunction
+                    << "Point " << pt << " is on a coupled patch";
+            }
+            Warning
+                << " with transformation " << patchTrafo
+                << " but also on 3 other patches with transforms "
+                << trafos << nl
+                << "This is not a space filling tiling and might"
+                << " indicate a setup problem and give problems"
+                << " for e.g. lagrangian tracking or interpolation" << endl;
+
+            // Already warned so no need to extend more
+            trafos.clear();
+            return false;
+        }
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
 }
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::globalIndexAndTransform::globalIndexAndTransform
+Foam::globalIndexAndTransform::globalIndexAndTransform(const polyMesh& mesh)
-(
-    const polyMesh& mesh
-)
 :
     mesh_(mesh),
     transforms_(),
@@ -546,13 +513,102 @@ Foam::globalIndexAndTransform::globalIndexAndTransform
         Info<< "nullTransformIndex:" << nullTransformIndex() << endl
             << endl;
     }
+
+
+    if (transforms_.size() > 0)
+    {
+        // Check that the transforms are space filling : any point
+        // can only use up to three transforms
+
+        const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+
+
+        // 1. Collect transform&sign per point and do local check
+
+        List<labelPairList> pointToTrafos(mesh_.nPoints());
+
+        forAll(patches, patchi)
+        {
+            const polyPatch& pp = patches[patchi];
+
+            const labelPair& transSign = patchTransformSign_[patchi];
+
+            if (transSign.first() > -1)
+            {
+                const labelList& mp = pp.meshPoints();
+                forAll(mp, i)
+                {
+                    labelPairList& trafos = pointToTrafos[mp[i]];
+
+                    bool newTransform = uniqueTransform
+                    (
+                        mesh_.points()[mp[i]],
+                        trafos,
+                        patchi,
+                        transSign
+                    );
+
+                    if (newTransform)
+                    {
+                        trafos.append(transSign);
+                    }
+                }
+            }
+        }
+
+
+        // Synchronise across collocated (= untransformed) points
+        // TBD: there is a big problem in that globalMeshData uses
+        //      globalIndexAndTransform. Triggers recursion.
+        if (false)
+        {
+            const globalMeshData& gmd = mesh_.globalData();
+            const indirectPrimitivePatch& cpp = gmd.coupledPatch();
+            const labelList& meshPoints = cpp.meshPoints();
+            const mapDistribute& slavesMap = gmd.globalCoPointSlavesMap();
+            const labelListList& slaves = gmd.globalCoPointSlaves();
+
+            List<labelPairList> elems(slavesMap.constructSize());
+            forAll(meshPoints, i)
+            {
+                elems[i] = pointToTrafos[meshPoints[i]];
+            }
+
+            // Pull slave data onto master. No need to update transformed slots.
+            slavesMap.distribute(elems, false);
+
+            // Combine master data with slave data
+            forAll(slaves, i)
+            {
+                labelPairList& trafos = elems[i];
+
+                const labelList& slavePoints = slaves[i];
+
+                // Combine master with untransformed slave data
+                forAll(slavePoints, j)
+                {
+                    const labelPairList& slaveTrafos = elems[slavePoints[j]];
+
+                    forAll(slaveTrafos, slaveI)
+                    {
+                        bool newTransform = uniqueTransform
+                        (
+                            mesh_.points()[meshPoints[i]],
+                            trafos,
+                            -1,
+                            slaveTrafos[slaveI]
+                        );
+
+                        if (newTransform)
+                        {
+                            trafos.append(slaveTrafos[slaveI]);
+                        }
+                    }
+                }
+            }
+        }
+    }
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::globalIndexAndTransform::~globalIndexAndTransform()
-{}
-
-
 // ************************************************************************* //

src/OpenFOAM/primitives/globalIndexAndTransform/globalIndexAndTransform.H

@@ -28,7 +28,9 @@ Description
     Determination and storage of the possible independent transforms
     introduced by coupledPolyPatches, as well as all of the possible
     permutations of these transforms generated by the presence of
-    multiple coupledPolyPatches, i.e. more than one cyclic boundary.
+    multiple coupledPolyPatches, i.e. more than one cyclic boundary. Note that
+    any given point can be on maximum 3 transforms only (and these transforms
+    have to be perpendicular)
 
     Also provides global index encoding and decoding for entity
     (i.e. cell) index, processor index and transform index (0 or
@@ -71,8 +73,15 @@ public:
         //  - transform
         class less
         {
+            const globalIndexAndTransform& gi_;
+
         public:
 
+            less(const globalIndexAndTransform& gi)
+            :
+                gi_(gi)
+            {}
+
             inline bool operator()(const labelPair&, const labelPair&) const;
         };
 
@@ -87,10 +96,8 @@ private:
         //- The possible independent (non-permuted) transforms of the
         //  geometry, i.e. for a geometry with two cyclics, this
         //  stores the two transforms, not the eight permutations.
-        //  There may not be more than three transforms in the range
+        //  Any point can not be on more than three transforms but overall
-        //  of coupledPolyPatch geometries (separated XOR
+        //  the mesh can have more than three.
-        //  non-parallel) and symmetries (cuboid periodicity only)
-        //  supported.
         List<vectorTensorTransform> transforms_;
 
         //- The permutations of the transforms, stored for lookup
@@ -105,13 +112,7 @@ private:
         //- Mapping from patch index to which transform it matches (or
         //  -1 for none) (.first()) and what sign to use for it,
         //  i.e. +/- 1 (.second()).
-        List<Pair<label>> patchTransformSign_;
+        List<labelPair> patchTransformSign_;
-
-
-    // Private static data
-
-        //- Number of spaces to reserve for transform encoding
-        static const label base_;
 
 
     // Private Member Functions
@@ -124,7 +125,7 @@ private:
         void determineTransformPermutations();
 
         //- Determine which patch uses which transform (if any) and which
-        //- Sign to use
+        //  sign to use
         void determinePatchTransformSign();
 
         //- Test a list of reference transforms to see if the test
@@ -139,16 +140,14 @@ private:
             bool checkBothSigns
         ) const;
 
-        //- Encode transform index. Hardcoded to 3 independent transforms max.
+        //- Return true if transform is not yet present in trafos. Issues
-        inline label encodeTransformIndex
+        //  warning if too many transforms
+        bool uniqueTransform
         (
-            const FixedList<Foam::label, 3>& permutationIndices
+            const point& pt,
-        ) const;
+            labelPairList& trafos,
-
+            const label patchi,
-        //- Decode transform index. Hardcoded to 3 independent transforms max.
+            const labelPair& patchTrafo
-        inline FixedList<label, 3> decodeTransformIndex
-        (
-            const label transformIndex
         ) const;
 
         //- Disallow default bitwise copy construct
@@ -160,10 +159,6 @@ private:
 
 public:
 
-        //- Declare friendship with the entry class for IO
-        friend class globalPoints;
-
-
     // Declare name of the class and its debug switch
     ClassName("globalIndexAndTransform");
 
@@ -174,10 +169,6 @@ public:
         globalIndexAndTransform(const polyMesh& mesh);
 
 
-    //- Destructor
-    ~globalIndexAndTransform();
-
-
     // Member Functions
 
         //- Generate a transform index from the permutation indices of
@@ -185,9 +176,12 @@ public:
         //  only be -1, 0 or +1.
         inline label encodeTransformIndex
         (
-            const List<label>& permutationIndices
+            const labelList& permutationIndices
         ) const;
 
+        //- Decode transform index
+        inline labelList decodeTransformIndex(const label transformIndex) const;
+
         //- Add patch transformation to transformIndex. Return new
         //  transformIndex. (by default the patch is the sending, not the
         //  receiving, patch)
@@ -221,74 +215,71 @@ public:
         ) const;
 
         //- Encode index and bare index as components on own processor
-        inline static labelPair encode
+        inline labelPair encode
         (
             const label index,
             const label transformIndex
-        );
+        ) const;
 
         //- Encode index and bare index as components on given processor
-        inline static labelPair encode
+        inline labelPair encode
         (
             const label proci,
             const label index,
             const label transformIndex
-        );
+        ) const;
 
         //- Index carried by the object
-        inline static label index(const labelPair& globalIAndTransform);
+        inline label index(const labelPair& globalIAndTransform) const;
 
         //- Which processor does this come from?
-        inline static label processor(const labelPair& globalIAndTransform);
+        inline label processor(const labelPair& globalIAndTransform) const;
 
         //- Transform carried by the object
-        inline static label transformIndex
+        inline label transformIndex(const labelPair& globalIAndTransform) const;
+
+        //- Return the number of independent transforms
+        inline label nIndependentTransforms() const;
+
+        //- Return access to the stored independent transforms
+        inline const List<vectorTensorTransform>& transforms() const;
+
+        //- Return access to the permuted transforms
+        inline const List<vectorTensorTransform>&
+        transformPermutations() const;
+
+        //- Return the transformIndex (index in transformPermutations)
+        //  of the identity transform
+        inline label nullTransformIndex() const;
+
+        //- Return access to the per-patch transform-sign pairs
+        inline const List<labelPair>& patchTransformSign() const;
+
+        //- Access the overall (permuted) transform corresponding
+        //  to the transformIndex
+        inline const vectorTensorTransform& transform
         (
-            const labelPair& globalIAndTransform
+            label transformIndex
-        );
+        ) const;
 
-        // Access
+        //- Access the all of the indices of the transform
+        //  permutations corresponding the transforms of the
+        //  listed patch indices. This only allows a maximum of three
+        //  transformations (since routine is used to transform points and
+        //  any given point can only be on 3 or less transforms)
+        inline labelList transformIndicesForPatches
+        (
+            const labelHashSet& patchIs
+        ) const;
 
-            //- Return the number of independent transforms
+        //- Apply all of the transform permutations
-            inline label nIndependentTransforms() const;
+        //  corresponding the transforms of the listed patch
-
+        //  indices to the supplied point
-            //- Return access to the stored independent transforms
+        inline pointField transformPatches
-            inline const List<vectorTensorTransform>& transforms() const;
+        (
-
+            const labelHashSet& patchIs,
-            //- Return access to the permuted transforms
+            const point& pt
-            inline const List<vectorTensorTransform>&
+        ) const;
-            transformPermutations() const;
-
-            //- Return the transformIndex (index in transformPermutations)
-            //  of the identity transform
-            inline label nullTransformIndex() const;
-
-            //- Return access to the per-patch transform-sign pairs
-            inline const List<Pair<label>>& patchTransformSign() const;
-
-            //- Access the overall (permuted) transform corresponding
-            //  to the transformIndex
-            inline const vectorTensorTransform& transform
-            (
-                label transformIndex
-            ) const;
-
-            //- Access the all of the indices of the transform
-            //  permutations corresponding the transforms of the
-            //  listed patch indices
-            inline labelList transformIndicesForPatches
-            (
-                const labelHashSet& patchIs
-            ) const;
-
-            //- Apply all of the transform permutations
-            //  corresponding the transforms of the listed patch
-            //  indices to the supplied point
-            inline pointField transformPatches
-            (
-                const labelHashSet& patchIs,
-                const point& pt
-            ) const;
 
 };
 

src/OpenFOAM/primitives/globalIndexAndTransform/globalIndexAndTransformI.H

@@ -33,8 +33,8 @@ bool Foam::globalIndexAndTransform::less::operator()
     const labelPair& b
 ) const
 {
-    label procA = globalIndexAndTransform::processor(a);
+    label procA = gi_.processor(a);
-    label procB = globalIndexAndTransform::processor(b);
+    label procB = gi_.processor(b);
 
     if (procA < procB)
     {
@@ -47,8 +47,8 @@ bool Foam::globalIndexAndTransform::less::operator()
     else
     {
         // Equal proc.
-        label indexA = globalIndexAndTransform::index(a);
+        label indexA = gi_.index(a);
-        label indexB = globalIndexAndTransform::index(b);
+        label indexB = gi_.index(b);
 
         if (indexA < indexB)
         {
@@ -61,8 +61,8 @@ bool Foam::globalIndexAndTransform::less::operator()
         else
         {
             // Equal index
-            label transformA = globalIndexAndTransform::transformIndex(a);
+            label transformA = gi_.transformIndex(a);
-            label transformB = globalIndexAndTransform::transformIndex(b);
+            label transformB = gi_.transformIndex(b);
 
             return transformA < transformB;
         }
@@ -72,7 +72,7 @@ bool Foam::globalIndexAndTransform::less::operator()
 
 Foam::label Foam::globalIndexAndTransform::encodeTransformIndex
 (
-    const List<label>& permutationIndices
+    const labelList& permutationIndices
 ) const
 {
     if (permutationIndices.size() != transforms_.size())
@@ -106,68 +106,20 @@ Foam::label Foam::globalIndexAndTransform::encodeTransformIndex
 }
 
 
-Foam::label Foam::globalIndexAndTransform::encodeTransformIndex
+Foam::labelList Foam::globalIndexAndTransform::decodeTransformIndex
-(
-    const FixedList<Foam::label, 3>& permutation
-) const
-{
-    if (nIndependentTransforms() == 0)
-    {
-        return 0;
-    }
-    if (nIndependentTransforms() == 1)
-    {
-        return permutation[0]+1;
-    }
-    else if (nIndependentTransforms() == 2)
-    {
-        return (permutation[1]+1)*3 + (permutation[0]+1);
-    }
-    else
-    {
-        return
-            (permutation[2]+1)*9
-          + (permutation[1]+1)*3
-          + (permutation[0]+1);
-    }
-}
-
-
-Foam::FixedList<Foam::label, 3>
-Foam::globalIndexAndTransform::decodeTransformIndex
 (
     const label transformIndex
 ) const
 {
-    FixedList<label, 3> permutation(label(0));
+    labelList permutation(transforms_.size(), 0);
 
     label t = transformIndex;
-    if (nIndependentTransforms() > 0)
+    forAll(permutation, i)
     {
-        permutation[0] = (t%3)-1;
+        permutation[i] = (t%3)-1;
-        if (nIndependentTransforms() > 1)
+        t /= 3;
-        {
-            t /= 3;
-            permutation[1] = (t%3)-1;
-            if (nIndependentTransforms() > 2)
-            {
-                t /= 3;
-                permutation[2] = (t%3)-1;
-            }
-        }
     }
 
-    #ifdef FULLDEBUG
-    t /= 3;
-    if (t != 0)
-    {
-        FatalErrorInFunction
-            << "transformIndex : " << transformIndex
-            << " has more than 3 fields."
-            << abort(FatalError);
-    }
-    #endif
-
     return permutation;
 }
 
@@ -180,15 +132,28 @@ Foam::label Foam::globalIndexAndTransform::addToTransformIndex
     const scalar tol
 ) const
 {
-    const Pair<label>& transSign = patchTransformSign_[patchi];
+    const labelPair& transSign = patchTransformSign_[patchi];
 
     label matchTransI = transSign.first();
 
-    // Hardcoded for max 3 transforms only!
+    if (matchTransI >= transforms_.size())
-
-    if (matchTransI > -1 && matchTransI < 3)
     {
-        FixedList<label, 3> permutation = decodeTransformIndex(transformIndex);
+        FatalErrorInFunction
+            << "patch:" << mesh_.boundaryMesh()[patchi].name()
+            << " transform:" << matchTransI
+            << " out of possible transforms:" << transforms_
+            << exit(FatalError);
+        return labelMin;
+    }
+    else if (matchTransI == -1)
+    {
+        // No additional transformation for this patch
+        return transformIndex;
+    }
+    else
+    {
+        // Decode current set of transforms
+        labelList permutation(decodeTransformIndex(transformIndex));
 
 
         // Add patch transform
@@ -267,10 +232,6 @@ Foam::label Foam::globalIndexAndTransform::addToTransformIndex
 
         return encodeTransformIndex(permutation);
     }
-    else
-    {
-        return transformIndex;
-    }
 }
 
 
@@ -287,7 +248,7 @@ Foam::label Foam::globalIndexAndTransform::minimumTransformIndex
 
 
     // Count number of transforms
-    FixedList<label, 3> permutation0 = decodeTransformIndex(transformIndex0);
+    labelList permutation0(decodeTransformIndex(transformIndex0));
     label n0 = 0;
     forAll(permutation0, i)
     {
@@ -297,7 +258,7 @@ Foam::label Foam::globalIndexAndTransform::minimumTransformIndex
         }
     }
 
-    FixedList<label, 3> permutation1 = decodeTransformIndex(transformIndex1);
+    labelList permutation1(decodeTransformIndex(transformIndex1));
     label n1 = 0;
     forAll(permutation1, i)
     {
@@ -324,8 +285,8 @@ Foam::label Foam::globalIndexAndTransform::subtractTransformIndex
     const label transformIndex1
 ) const
 {
-    FixedList<label, 3> permutation0 = decodeTransformIndex(transformIndex0);
+    labelList permutation0(decodeTransformIndex(transformIndex0));
-    FixedList<label, 3> permutation1 = decodeTransformIndex(transformIndex1);
+    labelList permutation1(decodeTransformIndex(transformIndex1));
 
     forAll(permutation0, i)
     {
@@ -340,7 +301,7 @@ Foam::labelPair Foam::globalIndexAndTransform::encode
 (
     const label index,
     const label transformIndex
-)
+) const
 {
     return encode(Pstream::myProcNo(), index, transformIndex);
 }
@@ -351,21 +312,22 @@ Foam::labelPair Foam::globalIndexAndTransform::encode
     const label proci,
     const label index,
     const label transformIndex
-)
+) const
 {
-    if (transformIndex < 0 || transformIndex >= base_)
+    if (transformIndex < 0 || transformIndex >= transformPermutations_.size())
     {
         FatalErrorInFunction
             << "TransformIndex " << transformIndex
             << " is outside allowed range of 0 to "
-            << base_ - 1
+            << transformPermutations_.size() - 1
             << abort(FatalError);
     }
 
-    if (proci > labelMax/base_)
+    if (proci > labelMax/transformPermutations_.size())
     {
         FatalErrorInFunction
-            << "Overflow : encoding processor " << proci << " in base " << base_
+            << "Overflow : encoding processor " << proci
+            << " in base " << transformPermutations_.size()
             << " exceeds capability of label (" << labelMax
             << "). Please recompile with larger datatype for label."
             << exit(FatalError);
@@ -374,7 +336,7 @@ Foam::labelPair Foam::globalIndexAndTransform::encode
     return labelPair
     (
         index,
-        transformIndex + proci*base_
+        transformIndex + proci*transformPermutations_.size()
     );
 }
 
@@ -382,7 +344,7 @@ Foam::labelPair Foam::globalIndexAndTransform::encode
 Foam::label Foam::globalIndexAndTransform::index
 (
     const labelPair& globalIAndTransform
-)
+) const
 {
     return globalIAndTransform.first();
 }
@@ -391,18 +353,18 @@ Foam::label Foam::globalIndexAndTransform::index
 Foam::label Foam::globalIndexAndTransform::processor
 (
     const labelPair& globalIAndTransform
-)
+) const
 {
-    return globalIAndTransform.second()/base_;
+    return globalIAndTransform.second()/transformPermutations_.size();
 }
 
 
 Foam::label Foam::globalIndexAndTransform::transformIndex
 (
     const labelPair& globalIAndTransform
-)
+) const
 {
-    return globalIAndTransform.second() % base_;
+    return globalIAndTransform.second()%transformPermutations_.size();
 }
 
 
@@ -432,7 +394,7 @@ Foam::label Foam::globalIndexAndTransform::nullTransformIndex() const
 }
 
 
-const Foam::List<Foam::Pair<Foam::label>>&
+const Foam::labelPairList&
 Foam::globalIndexAndTransform::patchTransformSign() const
 {
     return patchTransformSign_;
@@ -453,7 +415,7 @@ Foam::labelList Foam::globalIndexAndTransform::transformIndicesForPatches
     const labelHashSet& patchis
 ) const
 {
-    List<label> permutation(transforms_.size(), 0);
+    labelList permutation(transforms_.size(), 0);
 
     labelList selectedTransformIs(0);
 
@@ -466,7 +428,7 @@ Foam::labelList Foam::globalIndexAndTransform::transformIndicesForPatches
     {
         label patchi = iter.key();
 
-        const Pair<label>& transSign = patchTransformSign_[patchi];
+        const labelPair& transSign = patchTransformSign_[patchi];
 
         label matchTransI = transSign.first();
 
@@ -520,7 +482,7 @@ Foam::labelList Foam::globalIndexAndTransform::transformIndicesForPatches
         }
         case 2:
         {
-            List<label> tempPermutation = permutation;
+            labelList tempPermutation = permutation;
 
             label a = 0;
             label b = 1;
@@ -565,7 +527,7 @@ Foam::labelList Foam::globalIndexAndTransform::transformIndicesForPatches
         }
         case 3:
         {
-            List<label> tempPermutation = permutation;
+            labelList tempPermutation = permutation;
 
             tempPermutation[0] = 0;
             tempPermutation[1] = 0;

src/TurbulenceModels/compressible/turbulentFluidThermoModels/makeTurbulenceModel.H

@@ -34,7 +34,8 @@ License
     {                                                                          \
         typedef TDModel<BaseModel<Transport>>                                  \
             Transport##BaseModel;                                              \
-        typedef laminar<Transport##BaseModel> Laminar##Transport##BaseModel;   \
+        typedef laminarModel<Transport##BaseModel>                             \
+            laminar##Transport##BaseModel;                                     \
         typedef RASModel<EddyDiffusivity<Transport##BaseModel>>                \
             RAS##Transport##BaseModel;                                         \
         typedef LESModel<EddyDiffusivity<Transport##BaseModel>>                \
@@ -65,14 +66,18 @@ License
             Transport##BaseModel;                                              \
                                                                                \
                                                                                \
-        typedef laminar<Transport##BaseModel> Laminar##Transport##BaseModel;   \
+        typedef laminarModel<Transport##BaseModel>                             \
+            laminar##Transport##BaseModel;                                     \
                                                                                \
-        defineNamedTemplateTypeNameAndDebug(Laminar##Transport##BaseModel, 0); \
+        defineNamedTemplateTypeNameAndDebug(laminar##Transport##BaseModel, 0); \
+                                                                               \
+        defineTemplateRunTimeSelectionTable                                    \
+        (laminar##Transport##BaseModel, dictionary);                           \
                                                                                \
         addToRunTimeSelectionTable                                             \
         (                                                                      \
             Transport##baseModel,                                              \
-            Laminar##Transport##BaseModel,                                     \
+            laminar##Transport##BaseModel,                                     \
             dictionary                                                         \
         );                                                                     \
                                                                                \
@@ -110,6 +115,27 @@ License
     }
 
 
+#define makeTemplatedLaminarModel(BaseModel, SType, Type)                      \
+    typedef Foam::SType##Models::Type<Foam::BaseModel>                         \
+        Type##SType##BaseModel;                                                \
+    defineNamedTemplateTypeNameAndDebug(Type##SType##BaseModel, 0);            \
+                                                                               \
+    namespace Foam                                                             \
+    {                                                                          \
+        namespace SType##Models                                                \
+        {                                                                      \
+            typedef Type<BaseModel> Type##SType##BaseModel;                    \
+                                                                               \
+            addToRunTimeSelectionTable                                         \
+            (                                                                  \
+                SType##BaseModel,                                              \
+                Type##SType##BaseModel,                                        \
+                dictionary                                                     \
+            );                                                                 \
+        }                                                                      \
+    }
+
+
 #define makeTemplatedTurbulenceModel(BaseModel, SType, Type)                   \
     typedef Foam::SType##Models::Type<Foam::EddyDiffusivity<Foam::BaseModel>> \
         Type##SType##BaseModel;                                                \

src/TurbulenceModels/compressible/turbulentFluidThermoModels/turbulentFluidThermoModel.H

@@ -46,6 +46,7 @@ SourceFiles
 #include "CompressibleTurbulenceModel.H"
 #include "ThermalDiffusivity.H"
 #include "EddyDiffusivity.H"
+#include "laminarModel.H"
 #include "RASModel.H"
 #include "LESModel.H"
 #include "fluidThermo.H"
@@ -59,6 +60,7 @@ namespace Foam
         typedef ThermalDiffusivity<CompressibleTurbulenceModel<fluidThermo>>
             turbulenceModel;
 
+        typedef laminarModel<turbulenceModel> laminarModel;
         typedef RASModel<EddyDiffusivity<turbulenceModel>> RASModel;
         typedef LESModel<EddyDiffusivity<turbulenceModel>> LESModel;
 

src/TurbulenceModels/compressible/turbulentFluidThermoModels/turbulentFluidThermoModels.C

@@ -37,6 +37,18 @@ makeBaseTurbulenceModel
     fluidThermo
 );
 
+
+// -------------------------------------------------------------------------- //
+// Laminar models
+// -------------------------------------------------------------------------- //
+
+#include "Stokes.H"
+makeLaminarModel(Stokes);
+
+#include "Maxwell.H"
+makeLaminarModel(Maxwell);
+
+
 // -------------------------------------------------------------------------- //
 // RAS models
 // -------------------------------------------------------------------------- //

src/TurbulenceModels/compressible/turbulentFluidThermoModels/turbulentFluidThermoModels.H

@@ -32,7 +32,7 @@ License
 #include "ThermalDiffusivity.H"
 #include "EddyDiffusivity.H"
 
-#include "laminar.H"
+#include "laminarModel.H"
 #include "RASModel.H"
 #include "LESModel.H"
 
@@ -48,6 +48,10 @@ makeTurbulenceModelTypes
     fluidThermo
 );
 
+#define makeLaminarModel(Type)                                                 \
+    makeTemplatedLaminarModel                                                  \
+    (fluidThermoCompressibleTurbulenceModel, laminar, Type)
+
 #define makeRASModel(Type)                                                     \
     makeTemplatedTurbulenceModel                                               \
     (fluidThermoCompressibleTurbulenceModel, RAS, Type)

src/TurbulenceModels/incompressible/turbulentTransportModels/turbulentTransportModel.H

@@ -44,6 +44,7 @@ SourceFiles
 #define turbulentTransportModel_H
 
 #include "IncompressibleTurbulenceModel.H"
+#include "laminarModel.H"
 #include "RASModel.H"
 #include "LESModel.H"
 #include "incompressible/transportModel/transportModel.H"
@@ -56,6 +57,7 @@ namespace Foam
     {
         typedef IncompressibleTurbulenceModel<transportModel> turbulenceModel;
 
+        typedef laminarModel<turbulenceModel> laminarModel;
         typedef RASModel<turbulenceModel> RASModel;
         typedef LESModel<turbulenceModel> LESModel;
 

src/TurbulenceModels/incompressible/turbulentTransportModels/turbulentTransportModels.C

@@ -36,6 +36,18 @@ makeBaseTurbulenceModel
     transportModel
 );
 
+
+// -------------------------------------------------------------------------- //
+// Laminar models
+// -------------------------------------------------------------------------- //
+
+#include "Stokes.H"
+makeLaminarModel(Stokes);
+
+#include "Maxwell.H"
+makeLaminarModel(Maxwell);
+
+
 // -------------------------------------------------------------------------- //
 // RAS models
 // -------------------------------------------------------------------------- //

src/TurbulenceModels/incompressible/turbulentTransportModels/turbulentTransportModels.H

@@ -28,7 +28,7 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "makeTurbulenceModel.H"
 
-#include "laminar.H"
+#include "laminarModel.H"
 #include "RASModel.H"
 #include "LESModel.H"
 
@@ -43,6 +43,10 @@ makeTurbulenceModelTypes
     transportModel
 );
 
+#define makeLaminarModel(Type)                                                 \
+    makeTemplatedTurbulenceModel                                               \
+    (transportModelIncompressibleTurbulenceModel, laminar, Type)
+
 #define makeRASModel(Type)                                                     \
     makeTemplatedTurbulenceModel                                               \
     (transportModelIncompressibleTurbulenceModel, RAS, Type)

src/TurbulenceModels/schemes/Make/options

@@ -5,7 +5,8 @@ EXE_INC = \
     -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
 
 LIB_LIBS = \
     -lcompressibleTransportModels \
@@ -14,4 +15,5 @@ LIB_LIBS = \
     -lincompressibleTransportModels \
     -lcompressibleTurbulenceModels \
     -lfluidThermophysicalModels \
-    -lfiniteVolume
+    -lfiniteVolume \
+    -lmeshTools

src/TurbulenceModels/turbulenceModels/Base/kOmegaSST/kOmegaSSTBase.C

@@ -142,14 +142,31 @@ template<class BasicEddyViscosityModel>
 tmp<volScalarField::Internal>
 kOmegaSSTBase<BasicEddyViscosityModel>::epsilonByk
 (
-    const volScalarField::Internal& F1,
+    const volScalarField& F1,
-    const volScalarField::Internal& F2
+    const volTensorField& gradU
 ) const
 {
     return betaStar_*omega_();
 }
 
 
+template<class BasicEddyViscosityModel>
+tmp<volScalarField::Internal> kOmegaSSTBase<BasicEddyViscosityModel>::GbyNu
+(
+    const volScalarField::Internal& GbyNu0,
+    const volScalarField::Internal& F2,
+    const volScalarField::Internal& S2
+) const
+{
+    return min
+    (
+        GbyNu0,
+        (c1_/a1_)*betaStar_*omega_()
+       *max(a1_*omega_(), b1_*F2*sqrt(S2))
+    );
+}
+
+
 template<class BasicEddyViscosityModel>
 tmp<fvScalarMatrix> kOmegaSSTBase<BasicEddyViscosityModel>::kSource() const
 {
@@ -181,9 +198,9 @@ tmp<fvScalarMatrix> kOmegaSSTBase<BasicEddyViscosityModel>::omegaSource() const
 template<class BasicEddyViscosityModel>
 tmp<fvScalarMatrix> kOmegaSSTBase<BasicEddyViscosityModel>::Qsas
 (
-    const volScalarField& S2,
+    const volScalarField::Internal& S2,
-    const volScalarField& gamma,
+    const volScalarField::Internal& gamma,
-    const volScalarField& beta
+    const volScalarField::Internal& beta
 ) const
 {
     return tmp<fvScalarMatrix>
@@ -422,13 +439,12 @@ void kOmegaSSTBase<BasicEddyViscosityModel>::correct()
 
     BasicEddyViscosityModel::correct();
 
-    volScalarField divU(fvc::div(fvc::absolute(this->phi(), U)));
+    volScalarField::Internal divU(fvc::div(fvc::absolute(this->phi(), U)));
 
     tmp<volTensorField> tgradU = fvc::grad(U);
     volScalarField S2(2*magSqr(symm(tgradU())));
-    volScalarField GbyNu((tgradU() && dev(twoSymm(tgradU()))));
+    volScalarField::Internal GbyNu0((tgradU() && dev(twoSymm(tgradU()))));
-    volScalarField G(this->GName(), nut*GbyNu);
+    volScalarField::Internal G(this->GName(), nut*GbyNu0);
-    tgradU.clear();
 
     // Update omega and G at the wall
     omega_.boundaryFieldRef().updateCoeffs();
@@ -439,10 +455,11 @@ void kOmegaSSTBase<BasicEddyViscosityModel>::correct()
     );
 
     volScalarField F1(this->F1(CDkOmega));
+    volScalarField F23(this->F23());
 
     {
-        volScalarField gamma(this->gamma(F1));
+        volScalarField::Internal gamma(this->gamma(F1));
-        volScalarField beta(this->beta(F1));
+        volScalarField::Internal beta(this->beta(F1));
 
         // Turbulent frequency equation
         tmp<fvScalarMatrix> omegaEqn
@@ -451,20 +468,15 @@ void kOmegaSSTBase<BasicEddyViscosityModel>::correct()
           + fvm::div(alphaRhoPhi, omega_)
           - fvm::laplacian(alpha*rho*DomegaEff(F1), omega_)
          ==
-            alpha*rho*gamma
+            alpha()*rho()*gamma*GbyNu(GbyNu0, F23(), S2())
-           *min
+          - fvm::SuSp((2.0/3.0)*alpha()*rho()*gamma*divU, omega_)
-            (
+          - fvm::Sp(alpha()*rho()*beta*omega_(), omega_)
-                GbyNu,
-                (c1_/a1_)*betaStar_*omega_*max(a1_*omega_, b1_*F23()*sqrt(S2))
-            )
-          - fvm::SuSp((2.0/3.0)*alpha*rho*gamma*divU, omega_)
-          - fvm::Sp(alpha*rho*beta*omega_, omega_)
           - fvm::SuSp
             (
-                alpha*rho*(F1 - scalar(1))*CDkOmega/omega_,
+                alpha()*rho()*(F1() - scalar(1))*CDkOmega()/omega_(),
                 omega_
             )
-          + Qsas(S2, gamma, beta)
+          + Qsas(S2(), gamma, beta)
           + omegaSource()
           + fvOptions(alpha, rho, omega_)
         );
@@ -484,13 +496,15 @@ void kOmegaSSTBase<BasicEddyViscosityModel>::correct()
       + fvm::div(alphaRhoPhi, k_)
       - fvm::laplacian(alpha*rho*DkEff(F1), k_)
      ==
-        min(alpha*rho*G, (c1_*betaStar_)*alpha*rho*k_*omega_)
+        alpha()*rho()*Pk(G)
-      - fvm::SuSp((2.0/3.0)*alpha*rho*divU, k_)
+      - fvm::SuSp((2.0/3.0)*alpha()*rho()*divU, k_)
-      - fvm::Sp(alpha*rho*betaStar_*omega_, k_)
+      - fvm::Sp(alpha()*rho()*epsilonByk(F1, tgradU()), k_)
       + kSource()
       + fvOptions(alpha, rho, k_)
     );
 
+    tgradU.clear();
+
     kEqn.ref().relax();
     fvOptions.constrain(kEqn.ref());
     solve(kEqn);

src/TurbulenceModels/turbulenceModels/Base/kOmegaSST/kOmegaSSTBase.H

@@ -100,9 +100,6 @@ SourceFiles
 #ifndef kOmegaSSTBase_H
 #define kOmegaSSTBase_H
 
-#include "RASModel.H"
-#include "eddyViscosity.H"
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -121,7 +118,7 @@ class kOmegaSSTBase
 
         // Disallow default bitwise copy construct and assignment
         kOmegaSSTBase(const kOmegaSSTBase&);
-        kOmegaSSTBase& operator=(const kOmegaSSTBase&);
+        void operator=(const kOmegaSSTBase&);
 
 
 protected:
@@ -164,10 +161,10 @@ protected:
 
     // Protected Member Functions
 
-        tmp<volScalarField> F1(const volScalarField& CDkOmega) const;
+        virtual tmp<volScalarField> F1(const volScalarField& CDkOmega) const;
-        tmp<volScalarField> F2() const;
+        virtual tmp<volScalarField> F2() const;
-        tmp<volScalarField> F3() const;
+        virtual tmp<volScalarField> F3() const;
-        tmp<volScalarField> F23() const;
+        virtual tmp<volScalarField> F23() const;
 
         tmp<volScalarField> blend
         (
@@ -179,6 +176,16 @@ protected:
             return F1*(psi1 - psi2) + psi2;
         }
 
+        tmp<volScalarField::Internal> blend
+        (
+            const volScalarField::Internal& F1,
+            const dimensionedScalar& psi1,
+            const dimensionedScalar& psi2
+        ) const
+        {
+            return F1*(psi1 - psi2) + psi2;
+        }
+
         tmp<volScalarField> alphaK(const volScalarField& F1) const
         {
             return blend(F1, alphaK1_, alphaK2_);
@@ -189,12 +196,18 @@ protected:
             return blend(F1, alphaOmega1_, alphaOmega2_);
         }
 
-        tmp<volScalarField> beta(const volScalarField& F1) const
+        tmp<volScalarField::Internal> beta
+        (
+            const volScalarField::Internal& F1
+        ) const
         {
             return blend(F1, beta1_, beta2_);
         }
 
-        tmp<volScalarField> gamma(const volScalarField& F1) const
+        tmp<volScalarField::Internal> gamma
+        (
+            const volScalarField::Internal& F1
+        ) const
         {
             return blend(F1, gamma1_, gamma2_);
         }
@@ -212,8 +225,16 @@ protected:
         //- Return epsilon/k which for standard RAS is betaStar*omega
         virtual tmp<volScalarField::Internal> epsilonByk
         (
-            const volScalarField::Internal& F1,
+            const volScalarField& F1,
-            const volScalarField::Internal& F2
+            const volTensorField& gradU
+        ) const;
+
+        //- Return G/nu
+        virtual tmp<volScalarField::Internal> GbyNu
+        (
+            const volScalarField::Internal& GbyNu0,
+            const volScalarField::Internal& F2,
+            const volScalarField::Internal& S2
         ) const;
 
         virtual tmp<fvScalarMatrix> kSource() const;
@@ -222,9 +243,9 @@ protected:
 
         virtual tmp<fvScalarMatrix> Qsas
         (
-            const volScalarField& S2,
+            const volScalarField::Internal& S2,
-            const volScalarField& gamma,
+            const volScalarField::Internal& gamma,
-            const volScalarField& beta
+            const volScalarField::Internal& beta
         ) const;
 
 
@@ -326,7 +347,7 @@ public:
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #ifdef NoRepository
-#   include "kOmegaSSTBase.C"
+    #include "kOmegaSSTBase.C"
 #endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/TurbulenceModels/turbulenceModels/DES/kOmegaSSTDDES/kOmegaSSTDDES.H

@@ -72,7 +72,7 @@ class kOmegaSSTDDES
 
         // Disallow default bitwise copy construct and assignment
         kOmegaSSTDDES(const kOmegaSSTDDES&);
-        kOmegaSSTDDES& operator=(const kOmegaSSTDDES&);
+        void operator=(const kOmegaSSTDDES&);
 
 
 protected:

src/TurbulenceModels/turbulenceModels/DES/kOmegaSSTDES/kOmegaSSTDES.C

@@ -66,6 +66,30 @@ tmp<volScalarField> kOmegaSSTDES<BasicTurbulenceModel>::dTilda
 }
 
 
+template<class BasicTurbulenceModel>
+tmp<volScalarField::Internal> kOmegaSSTDES<BasicTurbulenceModel>::epsilonByk
+(
+    const volScalarField& F1,
+    const volTensorField& gradU
+) const
+{
+    volScalarField CDES(this->CDES(F1));
+    return sqrt(this->k_())/dTilda(mag(gradU), CDES)()();
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField::Internal> kOmegaSSTDES<BasicTurbulenceModel>::GbyNu
+(
+    const volScalarField::Internal& GbyNu0,
+    const volScalarField::Internal& F2,
+    const volScalarField::Internal& S2
+) const
+{
+    return GbyNu0; // Unlimited
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class BasicTurbulenceModel>
@@ -150,96 +174,6 @@ bool kOmegaSSTDES<BasicTurbulenceModel>::read()
 }
 
 
-template<class BasicTurbulenceModel>
-void kOmegaSSTDES<BasicTurbulenceModel>::correct()
-{
-    if (!this->turbulence_)
-    {
-        return;
-    }
-
-    // Local references
-    const alphaField& alpha = this->alpha_;
-    const rhoField& rho = this->rho_;
-    const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
-    const volVectorField& U = this->U_;
-    volScalarField& k = this->k_;
-    volScalarField& omega = this->omega_;
-    volScalarField& nut = this->nut_;
-
-    DESModel<BasicTurbulenceModel>::correct();
-
-    volScalarField divU(fvc::div(fvc::absolute(this->phi(), U)));
-
-    tmp<volTensorField> tgradU = fvc::grad(U);
-    volScalarField magGradU(mag(tgradU()));
-    volScalarField S2(2*magSqr(symm(tgradU())));
-    volScalarField GbyNu((tgradU() && dev(twoSymm(tgradU()))));
-    volScalarField G(this->GName(), nut*GbyNu);
-    tgradU.clear();
-
-    // Update omega and G at the wall
-    omega.boundaryFieldRef().updateCoeffs();
-
-    volScalarField CDkOmega
-    (
-        (2*this->alphaOmega2_)*(fvc::grad(k) & fvc::grad(omega))/omega
-    );
-
-    volScalarField F1(this->F1(CDkOmega));
-
-    {
-        volScalarField gamma(this->gamma(F1));
-        volScalarField beta(this->beta(F1));
-
-        // Turbulent frequency equation
-        tmp<fvScalarMatrix> omegaEqn
-        (
-            fvm::ddt(alpha, rho, omega)
-          + fvm::div(alphaRhoPhi, omega)
-          - fvm::laplacian(alpha*rho*this->DomegaEff(F1), omega)
-         ==
-            alpha*rho*gamma*GbyNu // Using unlimited GybNu
-          - fvm::SuSp((2.0/3.0)*alpha*rho*gamma*divU, omega)
-          - fvm::Sp(alpha*rho*beta*omega, omega)
-          - fvm::SuSp(alpha*rho*(F1 - scalar(1))*CDkOmega/omega, omega)
-          + this->omegaSource()
-        );
-
-        omegaEqn.ref().relax();
-
-        omegaEqn.ref().boundaryManipulate(omega.boundaryFieldRef());
-
-        solve(omegaEqn);
-        bound(omega, this->omegaMin_);
-    }
-
-    {
-        volScalarField CDES(this->CDES(F1));
-        volScalarField dTilda(this->dTilda(magGradU, CDES));
-
-        // Turbulent kinetic energy equation
-        tmp<fvScalarMatrix> kEqn
-        (
-            fvm::ddt(alpha, rho, k)
-          + fvm::div(alphaRhoPhi, k)
-          - fvm::laplacian(alpha*rho*this->DkEff(F1), k)
-         ==
-            min(alpha*rho*G, (this->c1_*this->betaStar_)*alpha*rho*k*omega)
-          - fvm::SuSp((2.0/3.0)*alpha*rho*divU, k)
-          - fvm::Sp(alpha*rho*sqrt(k)/dTilda, k)  // modified for DES
-          + this->kSource()
-        );
-
-        kEqn.ref().relax();
-        solve(kEqn);
-        bound(k, this->kMin_);
-    }
-
-    this->correctNut(S2);
-}
-
-
 template<class BasicTurbulenceModel>
 tmp<volScalarField> kOmegaSSTDES<BasicTurbulenceModel>::LESRegion() const
 {

src/TurbulenceModels/turbulenceModels/DES/kOmegaSSTDES/kOmegaSSTDES.H

@@ -74,7 +74,7 @@ class kOmegaSSTDES
 
         // Disallow default bitwise copy construct and assignment
         kOmegaSSTDES(const kOmegaSSTDES&);
-        kOmegaSSTDES& operator=(const kOmegaSSTDES&);
+        void operator=(const kOmegaSSTDES&);
 
 
 protected:
@@ -106,6 +106,21 @@ protected:
             const volScalarField& CDES
         ) const;
 
+        //- Return epsilon/k
+        virtual tmp<volScalarField::Internal> epsilonByk
+        (
+            const volScalarField& F1,
+            const volTensorField& gradU
+        ) const;
+
+        //- Return G/nu
+        virtual tmp<volScalarField::Internal> GbyNu
+        (
+            const volScalarField::Internal& GbyNu0,
+            const volScalarField::Internal& F2,
+            const volScalarField::Internal& S2
+        ) const;
+
 
 public:
 
@@ -144,9 +159,6 @@ public:
         //- Re-read model coefficients if they have changed
         virtual bool read();
 
-        //- Solve the turbulence equations and correct the turbulence viscosity
-        virtual void correct();
-
         //- Return the LES field indicator
         virtual tmp<volScalarField> LESRegion() const;
 };
@@ -159,7 +171,7 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 #ifdef NoRepository
-#   include "kOmegaSSTDES.C"
+    #include "kOmegaSSTDES.C"
 #endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/TurbulenceModels/turbulenceModels/DES/kOmegaSSTIDDES/kOmegaSSTIDDES.H

@@ -85,7 +85,7 @@ class kOmegaSSTIDDES
 
         // Disallow default bitwise copy construct and assignment
         kOmegaSSTIDDES(const kOmegaSSTIDDES&);
-        kOmegaSSTIDDES& operator=(const kOmegaSSTIDDES&);
+        void operator=(const kOmegaSSTIDDES&);
 
 
 protected:
@@ -160,7 +160,7 @@ public:
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #ifdef NoRepository
-#   include "kOmegaSSTIDDES.C"
+    #include "kOmegaSSTIDDES.C"
 #endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/TurbulenceModels/turbulenceModels/RAS/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -56,7 +56,7 @@ atmBoundaryLayerInletEpsilonFvPatchScalarField
     const dictionary& dict
 )
 :
-    fixedValueFvPatchScalarField(p, iF),
+    fixedValueFvPatchScalarField(p, iF, dict, false),
     atmBoundaryLayer(patch().Cf(), dict)
 {
     scalarField::operator=(epsilon(patch().Cf()));

src/TurbulenceModels/turbulenceModels/RAS/derivedFvPatchFields/atmBoundaryLayerInletK/atmBoundaryLayerInletKFvPatchScalarField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -56,7 +56,7 @@ atmBoundaryLayerInletKFvPatchScalarField
     const dictionary& dict
 )
 :
-    fixedValueFvPatchScalarField(p, iF),
+    fixedValueFvPatchScalarField(p, iF, dict, false),
     atmBoundaryLayer(patch().Cf(), dict)
 {
     scalarField::operator=(k(patch().Cf()));

src/TurbulenceModels/turbulenceModels/RAS/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -56,7 +56,7 @@ atmBoundaryLayerInletVelocityFvPatchVectorField
     const dictionary& dict
 )
 :
-    fixedValueFvPatchVectorField(p, iF),
+    fixedValueFvPatchVectorField(p, iF, dict, false),
     atmBoundaryLayer(patch().Cf(), dict)
 {
     vectorField::operator=(U(patch().Cf()));

src/TurbulenceModels/turbulenceModels/RAS/kOmegaSSTLM/kOmegaSSTLM.C

@@ -60,13 +60,13 @@ tmp<volScalarField::Internal> kOmegaSSTLM<BasicTurbulenceModel>::Pk
 template<class BasicTurbulenceModel>
 tmp<volScalarField::Internal> kOmegaSSTLM<BasicTurbulenceModel>::epsilonByk
 (
-    const volScalarField::Internal& F1,
+    const volScalarField& F1,
-    const volScalarField::Internal& F2
+    const volTensorField& gradU
 ) const
 {
     return
         min(max(gammaIntEff_, scalar(0.1)), scalar(1))
-       *kOmegaSST<BasicTurbulenceModel>::epsilonByk(F1, F2);
+       *kOmegaSST<BasicTurbulenceModel>::epsilonByk(F1, gradU);
 }
 
 

src/TurbulenceModels/turbulenceModels/RAS/kOmegaSSTLM/kOmegaSSTLM.H

@@ -164,8 +164,8 @@ protected:
         //- Modified form of the k-omega SST epsilon/k
         virtual tmp<volScalarField::Internal> epsilonByk
         (
-            const volScalarField::Internal& F1,
+            const volScalarField& F1,
-            const volScalarField::Internal& F2
+            const volTensorField& gradU
         ) const;
 
         //- Freestream blending-function

src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/fixedShearStress/fixedShearStressFvPatchVectorField.C

@@ -50,7 +50,7 @@ Foam::fixedShearStressFvPatchVectorField::fixedShearStressFvPatchVectorField
     const dictionary& dict
 )
 :
-    fixedValueFvPatchVectorField(p, iF),
+    fixedValueFvPatchVectorField(p, iF, dict, false),
     tau0_(dict.lookupOrDefault<vector>("tau", Zero))
 {
     fvPatchField<vector>::operator=(patchInternalField());

src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.H

@@ -28,7 +28,7 @@ Group
     grpWallFunctions
 
 Description
-    This boundary condition provides a wall constraint on turbulnce specific
+    This boundary condition provides a wall constraint on turbulence specific
     dissipation, omega for both low and high Reynolds number turbulence models.
 
     The near-wall omega may be either blended between the viscous region and

src/TurbulenceModels/turbulenceModels/laminar/Maxwell/Maxwell.C

@@ -0,0 +1,240 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "Maxwell.H"
+#include "fvOptions.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace laminarModels
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+Maxwell<BasicTurbulenceModel>::Maxwell
+(
+    const alphaField& alpha,
+    const rhoField& rho,
+    const volVectorField& U,
+    const surfaceScalarField& alphaRhoPhi,
+    const surfaceScalarField& phi,
+    const transportModel& transport,
+    const word& propertiesName,
+    const word& type
+)
+:
+    laminarModel<BasicTurbulenceModel>
+    (
+        type,
+        alpha,
+        rho,
+        U,
+        alphaRhoPhi,
+        phi,
+        transport,
+        propertiesName
+    ),
+
+    nuM_
+    (
+        dimensioned<scalar>
+        (
+            "nuM",
+            dimViscosity,
+            this->coeffDict_.lookup("nuM")
+        )
+    ),
+
+    lambda_
+    (
+        dimensioned<scalar>
+        (
+            "lambda",
+            dimTime,
+            this->coeffDict_.lookup("lambda")
+        )
+    ),
+
+    sigma_
+    (
+        IOobject
+        (
+            IOobject::groupName("sigma", U.group()),
+            this->runTime_.timeName(),
+            this->mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        this->mesh_
+    )
+{
+    if (type == typeName)
+    {
+        this->printCoeffs(type);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+bool Maxwell<BasicTurbulenceModel>::read()
+{
+    if (laminarModel<BasicTurbulenceModel>::read())
+    {
+        nuM_.readIfPresent(this->coeffDict());
+        lambda_.readIfPresent(this->coeffDict());
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+template<class BasicTurbulenceModel>
+tmp<Foam::volSymmTensorField>
+Maxwell<BasicTurbulenceModel>::R() const
+{
+    return sigma_;
+}
+
+template<class BasicTurbulenceModel>
+tmp<Foam::volSymmTensorField>
+Maxwell<BasicTurbulenceModel>::devRhoReff() const
+{
+    return tmp<volSymmTensorField>
+    (
+        new volSymmTensorField
+        (
+            IOobject
+            (
+                IOobject::groupName("devRhoReff", this->U_.group()),
+                this->runTime_.timeName(),
+                this->mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            this->alpha_*this->rho_*sigma_
+          - (this->alpha_*this->rho_*this->nu())
+           *dev(twoSymm(fvc::grad(this->U_)))
+        )
+    );
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<Foam::fvVectorMatrix>
+Maxwell<BasicTurbulenceModel>::divDevRhoReff
+(
+    volVectorField& U
+) const
+{
+    return
+    (
+        fvc::div
+        (
+            this->alpha_*this->rho_*this->nuM_*fvc::grad(U)
+        )
+      + fvc::div(this->alpha_*this->rho_*sigma_)
+      - fvc::div(this->alpha_*this->rho_*this->nu()*dev2(T(fvc::grad(U))))
+      - fvm::laplacian(this->alpha_*this->rho_*nu0(), U)
+    );
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<Foam::fvVectorMatrix>
+Maxwell<BasicTurbulenceModel>::divDevRhoReff
+(
+    const volScalarField& rho,
+    volVectorField& U
+) const
+{
+    return
+    (
+        fvc::div
+        (
+            this->alpha_*rho*this->nuM_*fvc::grad(U)
+        )
+      + fvc::div(this->alpha_*rho*sigma_)
+      - fvc::div(this->alpha_*rho*this->nu()*dev2(T(fvc::grad(U))))
+      - fvm::laplacian(this->alpha_*rho*nu0(), U)
+    );
+}
+
+
+template<class BasicTurbulenceModel>
+void Maxwell<BasicTurbulenceModel>::correct()
+{
+    // Local references
+    const alphaField& alpha = this->alpha_;
+    const rhoField& rho = this->rho_;
+    const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
+    const volVectorField& U = this->U_;
+    volSymmTensorField& sigma = this->sigma_;
+    fv::options& fvOptions(fv::options::New(this->mesh_));
+
+    laminarModel<BasicTurbulenceModel>::correct();
+
+    tmp<volTensorField> tgradU(fvc::grad(U));
+    const volTensorField& gradU = tgradU();
+    dimensionedScalar rLambda = 1.0/(lambda_);
+
+    // Note sigma is positive on lhs of momentum eqn
+    volSymmTensorField P
+    (
+        twoSymm(sigma & gradU)
+      - nuM_*rLambda*twoSymm(gradU)
+    );
+
+    // Viscoelastic stress equation
+    tmp<fvSymmTensorMatrix> sigmaEqn
+    (
+        fvm::ddt(alpha, rho, sigma)
+      + fvm::div(alphaRhoPhi, sigma)
+      + fvm::Sp(alpha*rho*rLambda, sigma)
+      ==
+        alpha*rho*P
+      + fvOptions(alpha, rho, sigma)
+    );
+
+    sigmaEqn.ref().relax();
+    fvOptions.constrain(sigmaEqn.ref());
+    solve(sigmaEqn);
+    fvOptions.correct(sigma_);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace laminarModels
+} // End namespace Foam
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/laminar/Maxwell/Maxwell.H

@@ -0,0 +1,174 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::laminarModels::Maxwell
+
+Group
+    grpLaminar
+
+Description
+    Maxwell model for viscoelasticity using the upper-convected time
+    derivative of the stress tensor.
+    See http://en.wikipedia.org/wiki/Upper-convected_Maxwell_model
+
+    The model includes an additional viscosity (nu) from the transport
+    model from which it is instantiated, which makes it equivalent to
+    the Oldroyd-B model for the case of an incompressible transport
+    model (where nu is non-zero).
+    See https://en.wikipedia.org/wiki/Oldroyd-B_model
+
+    Reference:
+    \verbatim
+        Amoreira, L. J., & Oliveira, P. J. (2010).
+        Comparison of different formulations for the numerical calculation
+        of unsteady incompressible viscoelastic fluid flow.
+        Adv. Appl. Math. Mech, 4, 483-502.
+    \endverbatim
+
+SourceFiles
+    Maxwell.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Maxwell_H
+#define Maxwell_H
+
+#include "laminarModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace laminarModels
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class Maxwell Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class BasicTurbulenceModel>
+class Maxwell
+:
+    public laminarModel<BasicTurbulenceModel>
+{
+
+protected:
+
+    // Protected data
+
+        // Model coefficients
+
+            dimensionedScalar nuM_;
+            dimensionedScalar lambda_;
+
+
+        // Fields
+
+            volSymmTensorField sigma_;
+
+
+    // Protected Member Functions
+
+        //- Return the turbulence viscosity
+        tmp<volScalarField> nu0() const
+        {
+            return this->nu() + nuM_;
+        }
+
+
+public:
+
+    typedef typename BasicTurbulenceModel::alphaField alphaField;
+    typedef typename BasicTurbulenceModel::rhoField rhoField;
+    typedef typename BasicTurbulenceModel::transportModel transportModel;
+
+
+    //- Runtime type information
+    TypeName("Maxwell");
+
+
+    // Constructors
+
+        //- Construct from components
+        Maxwell
+        (
+            const alphaField& alpha,
+            const rhoField& rho,
+            const volVectorField& U,
+            const surfaceScalarField& alphaRhoPhi,
+            const surfaceScalarField& phi,
+            const transportModel& transport,
+            const word& propertiesName = turbulenceModel::propertiesName,
+            const word& type = typeName
+        );
+
+
+    //- Destructor
+    virtual ~Maxwell()
+    {}
+
+
+    // Member Functions
+
+        //- Read model coefficients if they have changed
+        virtual bool read();
+
+        //- Return the Reynolds stress tensor
+        virtual tmp<volSymmTensorField> R() const;
+
+        //- Return the effective stress tensor
+        virtual tmp<volSymmTensorField> devRhoReff() const;
+
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff(volVectorField& U) const;
+
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff
+        (
+            const volScalarField& rho,
+            volVectorField& U
+        ) const;
+
+        //- Solve the turbulence equations and correct eddy-Viscosity and
+        //  related properties
+        virtual void correct();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace laminarModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "Maxwell.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/laminar/Stokes/Stokes.C

@@ -23,17 +23,24 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "laminar.H"
+#include "Stokes.H"
 #include "volFields.H"
 #include "surfaceFields.H"
 #include "fvcGrad.H"
 #include "fvcDiv.H"
 #include "fvmLaplacian.H"
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace laminarModels
+{
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class BasicTurbulenceModel>
-Foam::laminar<BasicTurbulenceModel>::laminar
+Stokes<BasicTurbulenceModel>::Stokes
 (
     const alphaField& alpha,
     const rhoField& rho,
@@ -44,7 +51,7 @@ Foam::laminar<BasicTurbulenceModel>::laminar
     const word& propertiesName
 )
 :
-    linearViscousStress<BasicTurbulenceModel>
+    linearViscousStress<laminarModel<BasicTurbulenceModel>>
     (
         typeName,
         alpha,
@@ -58,57 +65,26 @@ Foam::laminar<BasicTurbulenceModel>::laminar
 {}
 
 
-// * * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * //
-
-template<class BasicTurbulenceModel>
-Foam::autoPtr<Foam::laminar<BasicTurbulenceModel>>
-Foam::laminar<BasicTurbulenceModel>::New
-(
-    const alphaField& alpha,
-    const rhoField& rho,
-    const volVectorField& U,
-    const surfaceScalarField& alphaRhoPhi,
-    const surfaceScalarField& phi,
-    const transportModel& transport,
-    const word& propertiesName
-)
-{
-    return autoPtr<laminar>
-    (
-        new laminar
-        (
-            alpha,
-            rho,
-            U,
-            alphaRhoPhi,
-            phi,
-            transport,
-            propertiesName
-        )
-    );
-}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class BasicTurbulenceModel>
-const Foam::dictionary&
+const dictionary&
-Foam::laminar<BasicTurbulenceModel>::coeffDict() const
+Stokes<BasicTurbulenceModel>::coeffDict() const
 {
     return dictionary::null;
 }
 
 
 template<class BasicTurbulenceModel>
-bool Foam::laminar<BasicTurbulenceModel>::read()
+bool Stokes<BasicTurbulenceModel>::read()
 {
     return true;
 }
 
 
 template<class BasicTurbulenceModel>
-Foam::tmp<Foam::volScalarField>
+tmp<volScalarField>
-Foam::laminar<BasicTurbulenceModel>::nut() const
+Stokes<BasicTurbulenceModel>::nut() const
 {
     return tmp<volScalarField>
     (
@@ -131,8 +107,8 @@ Foam::laminar<BasicTurbulenceModel>::nut() const
 
 
 template<class BasicTurbulenceModel>
-Foam::tmp<Foam::scalarField>
+tmp<scalarField>
-Foam::laminar<BasicTurbulenceModel>::nut
+Stokes<BasicTurbulenceModel>::nut
 (
     const label patchi
 ) const
@@ -145,8 +121,8 @@ Foam::laminar<BasicTurbulenceModel>::nut
 
 
 template<class BasicTurbulenceModel>
-Foam::tmp<Foam::volScalarField>
+tmp<volScalarField>
-Foam::laminar<BasicTurbulenceModel>::nuEff() const
+Stokes<BasicTurbulenceModel>::nuEff() const
 {
     return tmp<volScalarField>
     (
@@ -159,8 +135,8 @@ Foam::laminar<BasicTurbulenceModel>::nuEff() const
 
 
 template<class BasicTurbulenceModel>
-Foam::tmp<Foam::scalarField>
+tmp<scalarField>
-Foam::laminar<BasicTurbulenceModel>::nuEff
+Stokes<BasicTurbulenceModel>::nuEff
 (
     const label patchi
 ) const
@@ -170,8 +146,8 @@ Foam::laminar<BasicTurbulenceModel>::nuEff
 
 
 template<class BasicTurbulenceModel>
-Foam::tmp<Foam::volScalarField>
+tmp<volScalarField>
-Foam::laminar<BasicTurbulenceModel>::k() const
+Stokes<BasicTurbulenceModel>::k() const
 {
     return tmp<volScalarField>
     (
@@ -194,8 +170,8 @@ Foam::laminar<BasicTurbulenceModel>::k() const
 
 
 template<class BasicTurbulenceModel>
-Foam::tmp<Foam::volScalarField>
+tmp<volScalarField>
-Foam::laminar<BasicTurbulenceModel>::epsilon() const
+Stokes<BasicTurbulenceModel>::epsilon() const
 {
     return tmp<volScalarField>
     (
@@ -221,8 +197,8 @@ Foam::laminar<BasicTurbulenceModel>::epsilon() const
 
 
 template<class BasicTurbulenceModel>
-Foam::tmp<Foam::volSymmTensorField>
+tmp<volSymmTensorField>
-Foam::laminar<BasicTurbulenceModel>::R() const
+Stokes<BasicTurbulenceModel>::R() const
 {
     return tmp<volSymmTensorField>
     (
@@ -248,10 +224,15 @@ Foam::laminar<BasicTurbulenceModel>::R() const
 
 
 template<class BasicTurbulenceModel>
-void Foam::laminar<BasicTurbulenceModel>::correct()
+void Stokes<BasicTurbulenceModel>::correct()
 {
-    BasicTurbulenceModel::correct();
+    laminarModel<BasicTurbulenceModel>::correct();
 }
 
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace laminarModels
+} // End namespace Foam
+
 // ************************************************************************* //

src/TurbulenceModels/turbulenceModels/laminar/Stokes/Stokes.H

@@ -22,34 +22,37 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::laminar
+    Foam::Stokes
 
 Description
-    Turbulence model for laminar flow.
+    Turbulence model for Stokes flow.
 
 SourceFiles
-    laminar.C
+    Stokes.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef laminar_H
+#ifndef Stokes_H
-#define laminar_H
+#define Stokes_H
 
+#include "laminarModel.H"
 #include "linearViscousStress.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
+namespace laminarModels
+{
 
 /*---------------------------------------------------------------------------* \
-                           Class laminar Declaration
+                           Class Stokes Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class BasicTurbulenceModel>
-class laminar
+class Stokes
 :
-    public linearViscousStress<BasicTurbulenceModel>
+    public linearViscousStress<laminarModel<BasicTurbulenceModel>>
 {
 
 public:
@@ -60,13 +63,13 @@ public:
 
 
     //- Runtime type information
-    TypeName("laminar");
+    TypeName("Stokes");
 
 
     // Constructors
 
         //- Construct from components
-        laminar
+        Stokes
         (
             const alphaField& alpha,
             const rhoField& rho,
@@ -81,7 +84,7 @@ public:
     // Selectors
 
         //- Return a reference to the selected turbulence model
-        static autoPtr<laminar> New
+        static autoPtr<Stokes> New
         (
             const alphaField& alpha,
             const rhoField& rho,
@@ -94,7 +97,7 @@ public:
 
 
     //- Destructor
-    virtual ~laminar()
+    virtual ~Stokes()
     {}
 
 
@@ -106,41 +109,42 @@ public:
         //- Read turbulenceProperties dictionary
         virtual bool read();
 
-        //- Return the turbulence viscosity, i.e. 0 for laminar flow
+        //- Return the turbulence viscosity, i.e. 0 for Stokes flow
         virtual tmp<volScalarField> nut() const;
 
         //- Return the turbulence viscosity on patch
         virtual tmp<scalarField> nut(const label patchi) const;
 
-        //- Return the effective viscosity, i.e. the laminar viscosity
+        //- Return the effective viscosity, i.e. the Stokes viscosity
         virtual tmp<volScalarField> nuEff() const;
 
         //- Return the effective viscosity on patch
         virtual tmp<scalarField> nuEff(const label patchi) const;
 
-        //- Return the turbulence kinetic energy, i.e. 0 for laminar flow
+        //- Return the turbulence kinetic energy, i.e. 0 for Stokes flow
         virtual tmp<volScalarField> k() const;
 
         //- Return the turbulence kinetic energy dissipation rate,
-        //  i.e. 0 for laminar flow
+        //  i.e. 0 for Stokes flow
         virtual tmp<volScalarField> epsilon() const;
 
-        //- Return the Reynolds stress tensor, i.e. 0 for laminar flow
+        //- Return the Reynolds stress tensor, i.e. 0 for Stokes flow
         virtual tmp<volSymmTensorField> R() const;
 
-        //- Correct the laminar viscosity
+        //- Correct the Stokes viscosity
         virtual void correct();
 };
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+} // End namespace laminarModels
 } // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #ifdef NoRepository
-    #include "laminar.C"
+    #include "Stokes.C"
 #endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/TurbulenceModels/turbulenceModels/laminar/laminarModel/laminarModel.C

@@ -0,0 +1,339 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "laminarModel.H"
+#include "Stokes.H"
+
+// * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+void Foam::laminarModel<BasicTurbulenceModel>::printCoeffs(const word& type)
+{
+    if (printCoeffs_)
+    {
+        Info<< type << "Coeffs" << coeffDict_ << endl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+Foam::laminarModel<BasicTurbulenceModel>::laminarModel
+(
+    const word& type,
+    const alphaField& alpha,
+    const rhoField& rho,
+    const volVectorField& U,
+    const surfaceScalarField& alphaRhoPhi,
+    const surfaceScalarField& phi,
+    const transportModel& transport,
+    const word& propertiesName
+)
+:
+    BasicTurbulenceModel
+    (
+        type,
+        alpha,
+        rho,
+        U,
+        alphaRhoPhi,
+        phi,
+        transport,
+        propertiesName
+    ),
+
+    laminarDict_(this->subOrEmptyDict("laminar")),
+    printCoeffs_(laminarDict_.lookupOrDefault<Switch>("printCoeffs", false)),
+    coeffDict_(laminarDict_.subOrEmptyDict(type + "Coeffs"))
+{
+    // Force the construction of the mesh deltaCoeffs which may be needed
+    // for the construction of the derived models and BCs
+    this->mesh_.deltaCoeffs();
+}
+
+
+// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+Foam::autoPtr<Foam::laminarModel<BasicTurbulenceModel>>
+Foam::laminarModel<BasicTurbulenceModel>::New
+(
+    const alphaField& alpha,
+    const rhoField& rho,
+    const volVectorField& U,
+    const surfaceScalarField& alphaRhoPhi,
+    const surfaceScalarField& phi,
+    const transportModel& transport,
+    const word& propertiesName
+)
+{
+    IOdictionary modelDict
+    (
+        IOobject
+        (
+            IOobject::groupName(propertiesName, U.group()),
+            U.time().constant(),
+            U.db(),
+            IOobject::MUST_READ_IF_MODIFIED,
+            IOobject::NO_WRITE,
+            false
+        )
+    );
+
+    if (modelDict.found("laminar"))
+    {
+        // get model name, but do not register the dictionary
+        // otherwise it is registered in the database twice
+        const word modelType
+        (
+            modelDict.subDict("laminar").lookup("laminarModel")
+        );
+
+        Info<< "Selecting laminar stress model " << modelType << endl;
+
+        typename dictionaryConstructorTable::iterator cstrIter =
+            dictionaryConstructorTablePtr_->find(modelType);
+
+        if (cstrIter == dictionaryConstructorTablePtr_->end())
+        {
+            FatalErrorInFunction
+                << "Unknown laminarModel type "
+                << modelType << nl << nl
+                << "Valid laminarModel types:" << endl
+                << dictionaryConstructorTablePtr_->sortedToc()
+                << exit(FatalError);
+        }
+
+        return autoPtr<laminarModel>
+        (
+            cstrIter()
+            (
+                alpha,
+                rho,
+                U,
+                alphaRhoPhi,
+                phi,
+                transport, propertiesName)
+        );
+    }
+    else
+    {
+        Info<< "Selecting laminar stress model "
+            << laminarModels::Stokes<BasicTurbulenceModel>::typeName << endl;
+
+        return autoPtr<laminarModel>
+        (
+            new laminarModels::Stokes<BasicTurbulenceModel>
+            (
+                alpha,
+                rho,
+                U,
+                alphaRhoPhi,
+                phi,
+                transport,
+                propertiesName
+            )
+        );
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+bool Foam::laminarModel<BasicTurbulenceModel>::read()
+{
+    if (BasicTurbulenceModel::read())
+    {
+        laminarDict_ <<= this->subDict("laminar");
+
+        if
+        (
+            const dictionary* dictPtr =
+                laminarDict_.subDictPtr(type() + "Coeffs")
+        )
+        {
+            coeffDict_ <<= *dictPtr;
+        }
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+template<class BasicTurbulenceModel>
+Foam::tmp<Foam::volScalarField>
+Foam::laminarModel<BasicTurbulenceModel>::nut() const
+{
+    return tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                IOobject::groupName("nut", this->U_.group()),
+                this->runTime_.timeName(),
+                this->mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            this->mesh_,
+            dimensionedScalar("nut", dimViscosity, 0.0)
+        )
+    );
+}
+
+
+template<class BasicTurbulenceModel>
+Foam::tmp<Foam::scalarField>
+Foam::laminarModel<BasicTurbulenceModel>::nut
+(
+    const label patchi
+) const
+{
+    return tmp<scalarField>
+    (
+        new scalarField(this->mesh_.boundary()[patchi].size(), 0.0)
+    );
+}
+
+
+template<class BasicTurbulenceModel>
+Foam::tmp<Foam::volScalarField>
+Foam::laminarModel<BasicTurbulenceModel>::nuEff() const
+{
+    return tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject::groupName("nuEff", this->U_.group()), this->nu()
+        )
+    );
+}
+
+
+template<class BasicTurbulenceModel>
+Foam::tmp<Foam::scalarField>
+Foam::laminarModel<BasicTurbulenceModel>::nuEff
+(
+    const label patchi
+) const
+{
+    return this->nu(patchi);
+}
+
+
+template<class BasicTurbulenceModel>
+Foam::tmp<Foam::volScalarField>
+Foam::laminarModel<BasicTurbulenceModel>::k() const
+{
+    return tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                IOobject::groupName("k", this->U_.group()),
+                this->runTime_.timeName(),
+                this->mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            this->mesh_,
+            dimensionedScalar("k", sqr(this->U_.dimensions()), 0.0)
+        )
+    );
+}
+
+
+template<class BasicTurbulenceModel>
+Foam::tmp<Foam::volScalarField>
+Foam::laminarModel<BasicTurbulenceModel>::epsilon() const
+{
+    return tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                IOobject::groupName("epsilon", this->U_.group()),
+                this->runTime_.timeName(),
+                this->mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            this->mesh_,
+            dimensionedScalar
+            (
+                "epsilon", sqr(this->U_.dimensions())/dimTime, 0.0
+            )
+        )
+    );
+}
+
+
+template<class BasicTurbulenceModel>
+Foam::tmp<Foam::volSymmTensorField>
+Foam::laminarModel<BasicTurbulenceModel>::R() const
+{
+    return tmp<volSymmTensorField>
+    (
+        new volSymmTensorField
+        (
+            IOobject
+            (
+                IOobject::groupName("R", this->U_.group()),
+                this->runTime_.timeName(),
+                this->mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            this->mesh_,
+            dimensionedSymmTensor
+            (
+                "R", sqr(this->U_.dimensions()), Zero
+            )
+        )
+    );
+}
+
+
+template<class BasicTurbulenceModel>
+void Foam::laminarModel<BasicTurbulenceModel>::correct()
+{
+    BasicTurbulenceModel::correct();
+}
+
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/laminar/laminarModel/laminarModel.H

@@ -0,0 +1,208 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::laminarModel
+
+Description
+    Templated abstract base class for laminar transport models
+
+SourceFiles
+    laminarModel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef laminarModel_H
+#define laminarModel_H
+
+#include "TurbulenceModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class laminarModel Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class BasicTurbulenceModel>
+class laminarModel
+:
+    public BasicTurbulenceModel
+{
+
+protected:
+
+    // Protected data
+
+        //- laminar coefficients dictionary
+        dictionary laminarDict_;
+
+        //- Flag to print the model coeffs at run-time
+        Switch printCoeffs_;
+
+        //- Model coefficients dictionary
+        dictionary coeffDict_;
+
+
+    // Protected Member Functions
+
+        //- Print model coefficients
+        virtual void printCoeffs(const word& type);
+
+
+private:
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        laminarModel(const laminarModel&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const laminarModel&);
+
+
+public:
+
+    typedef typename BasicTurbulenceModel::alphaField alphaField;
+    typedef typename BasicTurbulenceModel::rhoField rhoField;
+    typedef typename BasicTurbulenceModel::transportModel transportModel;
+
+
+    //- Runtime type information
+    TypeName("laminar");
+
+
+    // Declare run-time constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            laminarModel,
+            dictionary,
+            (
+                const alphaField& alpha,
+                const rhoField& rho,
+                const volVectorField& U,
+                const surfaceScalarField& alphaRhoPhi,
+                const surfaceScalarField& phi,
+                const transportModel& transport,
+                const word& propertiesName
+            ),
+            (alpha, rho, U, alphaRhoPhi, phi, transport, propertiesName)
+        );
+
+
+    // Constructors
+
+        //- Construct from components
+        laminarModel
+        (
+            const word& type,
+            const alphaField& alpha,
+            const rhoField& rho,
+            const volVectorField& U,
+            const surfaceScalarField& alphaRhoPhi,
+            const surfaceScalarField& phi,
+            const transportModel& transport,
+            const word& propertiesName
+        );
+
+
+    // Selectors
+
+        //- Return a reference to the selected laminar model
+        static autoPtr<laminarModel> New
+        (
+            const alphaField& alpha,
+            const rhoField& rho,
+            const volVectorField& U,
+            const surfaceScalarField& alphaRhoPhi,
+            const surfaceScalarField& phi,
+            const transportModel& transport,
+            const word& propertiesName = turbulenceModel::propertiesName
+        );
+
+
+    //- Destructor
+    virtual ~laminarModel()
+    {}
+
+
+    // Member Functions
+
+        //- Read model coefficients if they have changed
+        virtual bool read();
+
+
+        // Access
+
+            //- Const access to the coefficients dictionary
+            virtual const dictionary& coeffDict() const
+            {
+                return coeffDict_;
+            }
+
+            //- Return the turbulence viscosity, i.e. 0 for laminar flow
+            virtual tmp<volScalarField> nut() const;
+
+            //- Return the turbulence viscosity on patch
+            virtual tmp<scalarField> nut(const label patchi) const;
+
+            //- Return the effective viscosity, i.e. the laminar viscosity
+            virtual tmp<volScalarField> nuEff() const;
+
+            //- Return the effective viscosity on patch
+            virtual tmp<scalarField> nuEff(const label patchi) const;
+
+            //- Return the turbulence kinetic energy, i.e. 0 for laminar flow
+            virtual tmp<volScalarField> k() const;
+
+            //- Return the turbulence kinetic energy dissipation rate,
+            //  i.e. 0 for laminar flow
+            virtual tmp<volScalarField> epsilon() const;
+
+            //- Return the Reynolds stress tensor, i.e. 0 for laminar flow
+            virtual tmp<volSymmTensorField> R() const;
+
+            //- Correct the laminar transport
+            virtual void correct();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "laminarModel.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/laminar/laminarModel/laminarModelDoc.H

@@ -0,0 +1,32 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+\defgroup grpLaminar laminar transport model
+@{
+    \ingroup grpTurbulence
+    This group contains laminar models.
+@}
+
+\*---------------------------------------------------------------------------*/