Merge branch 'develop-pre-release' of develop.openfoam.com:Development/OpenFOAM-plus into develop-pre-release

applications/solvers/compressible/rhoPimpleFoam/Make/options

@@ -22,4 +22,5 @@ EXE_LIBS = \
     -lfvOptions \
     -ldynamicFvMesh \
     -ltopoChangerFvMesh \
-    -ldynamicMesh
+    -ldynamicMesh \
+    -latmosphericModels

applications/solvers/compressible/rhoSimpleFoam/Make/options

@@ -17,4 +17,5 @@ EXE_LIBS = \
     -lfiniteVolume \
     -lsampling \
     -lmeshTools \
-    -lfvOptions
+    -lfvOptions \
+    -latmosphericModels

applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/Make/options

@@ -18,4 +18,5 @@ EXE_LIBS = \
     -lfiniteVolume \
     -lfvOptions \
     -lsampling \
-    -lmeshTools
+    -lmeshTools \
+    -latmosphericModels

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/Make/options

@@ -16,4 +16,5 @@ EXE_LIBS = \
     -lfiniteVolume \
     -lsampling \
     -lmeshTools \
-    -lfvOptions
+    -lfvOptions \
+    -latmosphericModels

applications/solvers/heatTransfer/buoyantPimpleFoam/Make/options

@@ -18,4 +18,5 @@ EXE_LIBS = \
     -lradiationModels \
     -lspecie \
     -lturbulenceModels \
-    -lcompressibleTurbulenceModels
+    -lcompressibleTurbulenceModels \
+    -latmosphericModels

applications/solvers/heatTransfer/buoyantSimpleFoam/Make/options

@@ -19,4 +19,5 @@ EXE_LIBS = \
     -lradiationModels \
     -lturbulenceModels \
     -lcompressibleTurbulenceModels \
-    -lmeshTools
+    -lmeshTools \
+    -latmosphericModels

applications/solvers/incompressible/pimpleFoam/Make/options

@@ -19,4 +19,5 @@ EXE_LIBS = \
     -ldynamicFvMesh \
     -ltopoChangerFvMesh \
     -ldynamicMesh \
-    -lmeshTools
+    -lmeshTools \
+    -latmosphericModels

applications/solvers/incompressible/simpleFoam/Make/options

@@ -15,4 +15,5 @@ EXE_LIBS = \
     -lfiniteVolume \
     -lmeshTools \
     -lfvOptions \
-    -lsampling
+    -lsampling \
+    -latmosphericModels

applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C

@@ -426,7 +426,7 @@ Foam::multiphaseSystem::multiphaseSystem
                 iter(),
                 *phases_.lookup(iter.key().first()),
                 *phases_.lookup(iter.key().second())
-            ).ptr()
+            )
         );
     }
 
@@ -664,7 +664,7 @@ Foam::tmp<Foam::volVectorField> Foam::multiphaseSystem::Svm
 Foam::autoPtr<Foam::multiphaseSystem::dragCoeffFields>
 Foam::multiphaseSystem::dragCoeffs() const
 {
-    autoPtr<dragCoeffFields> dragCoeffsPtr(new dragCoeffFields);
+    auto dragCoeffsPtr = autoPtr<dragCoeffFields>::New();
 
     forAllConstIter(dragModelTable, dragModels_, iter)
     {
@@ -706,7 +706,7 @@ Foam::multiphaseSystem::dragCoeffs() const
             }
         }
 
-        dragCoeffsPtr().insert(iter.key(), Kptr);
+        dragCoeffsPtr().set(iter.key(), Kptr);
     }
 
     return dragCoeffsPtr;

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C

@@ -70,7 +70,7 @@ HeatAndMassTransferPhaseSystem
 
         // Initially assume no mass transfer
 
-        dmdt_.insert
+        dmdt_.set
         (
             pair,
             new volScalarField
@@ -88,7 +88,7 @@ HeatAndMassTransferPhaseSystem
             )
         );
 
-        dmdtExplicit_.insert
+        dmdtExplicit_.set
         (
             pair,
             new volScalarField
@@ -107,7 +107,7 @@ HeatAndMassTransferPhaseSystem
         volScalarField H1(heatTransferModels_[pair][pair.first()]->K());
         volScalarField H2(heatTransferModels_[pair][pair.second()]->K());
 
-        Tf_.insert
+        Tf_.set
         (
             pair,
             new volScalarField
@@ -257,18 +257,12 @@ template<class BasePhaseSystem>
 Foam::autoPtr<Foam::phaseSystem::heatTransferTable>
 Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::heatTransfer() const
 {
-    autoPtr<phaseSystem::heatTransferTable> eqnsPtr
+    auto eqnsPtr = autoPtr<phaseSystem::heatTransferTable>::New();
-    (
+    auto& eqns = *eqnsPtr;
-        new phaseSystem::heatTransferTable()
-    );
 
-    phaseSystem::heatTransferTable& eqns = eqnsPtr();
+    for (const phaseModel& phase : this->phaseModels_)
-
-    forAll(this->phaseModels_, phasei)
     {
-        const phaseModel& phase = this->phaseModels_[phasei];
+        eqns.set
-
-        eqns.insert
         (
             phase.name(),
             new fvScalarMatrix(phase.thermo().he(), dimEnergy/dimTime)

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.C

@@ -107,18 +107,12 @@ template<class BasePhaseSystem>
 Foam::autoPtr<Foam::phaseSystem::heatTransferTable>
 Foam::HeatTransferPhaseSystem<BasePhaseSystem>::heatTransfer() const
 {
-    autoPtr<phaseSystem::heatTransferTable> eqnsPtr
+    auto eqnsPtr = autoPtr<phaseSystem::heatTransferTable>::New();
-    (
+    auto& eqns = *eqnsPtr;
-        new phaseSystem::heatTransferTable()
-    );
 
-    phaseSystem::heatTransferTable& eqns = eqnsPtr();
+    for (const phaseModel& phase : this->phaseModels_)
-
-    forAll(this->phaseModels_, phasei)
     {
-        const phaseModel& phase = this->phaseModels_[phasei];
+        eqns.set
-
-        eqns.insert
         (
             phase.name(),
             new fvScalarMatrix(phase.thermo().he(), dimEnergy/dimTime)

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.C

@@ -62,22 +62,16 @@ Foam::InterfaceCompositionPhaseChangePhaseSystem<BasePhaseSystem>::
 massTransfer() const
 {
     // Create a mass transfer matrix for each species of each phase
-    autoPtr<phaseSystem::massTransferTable> eqnsPtr
+    auto eqnsPtr = autoPtr<phaseSystem::massTransferTable>::New();
-    (
+    auto& eqns = *eqnsPtr;
-        new phaseSystem::massTransferTable()
-    );
 
-    phaseSystem::massTransferTable& eqns = eqnsPtr();
+    for (const phaseModel& phase : this->phaseModels_)
-
-    forAll(this->phaseModels_, phasei)
     {
-        const phaseModel& phase = this->phaseModels_[phasei];
-
         const PtrList<volScalarField>& Yi = phase.Y();
 
         forAll(Yi, i)
         {
-            eqns.insert
+            eqns.set
             (
                 Yi[i].name(),
                 new fvScalarMatrix(Yi[i], dimMass/dimTime)

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C

@@ -87,7 +87,7 @@ MomentumTransferPhaseSystem
         const phasePair& pair =
             *(this->phasePairs_[dragModelIter.key()]);
 
-        Kds_.insert
+        Kds_.set
         (
             pair,
             new volScalarField
@@ -103,7 +103,7 @@ MomentumTransferPhaseSystem
         const phasePair& pair =
             *(this->phasePairs_[virtualMassModelIter.key()]);
 
-        Vms_.insert
+        Vms_.set
         (
             pair,
             new volScalarField
@@ -373,18 +373,12 @@ Foam::autoPtr<Foam::phaseSystem::momentumTransferTable>
 Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransfer() const
 {
     // Create a momentum transfer matrix for each phase
-    autoPtr<phaseSystem::momentumTransferTable> eqnsPtr
+    auto eqnsPtr = autoPtr<phaseSystem::momentumTransferTable>::New();
-    (
+    auto& eqns = *eqnsPtr;
-        new phaseSystem::momentumTransferTable()
-    );
 
-    phaseSystem::momentumTransferTable& eqns = eqnsPtr();
+    for (const phaseModel& phase : this->phaseModels_)
-
-    forAll(this->phaseModels_, phasei)
     {
-        const phaseModel& phase = this->phaseModels_[phasei];
+        eqns.set
-
-        eqns.insert
         (
             phase.name(),
             new fvVectorMatrix(phase.U(), dimMass*dimVelocity/dimTime)
@@ -492,11 +486,8 @@ template<class BasePhaseSystem>
 Foam::autoPtr<Foam::PtrList<Foam::volVectorField>>
 Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Fs() const
 {
-    autoPtr<PtrList<volVectorField>> tFs
+    auto tFs = autoPtr<PtrList<volVectorField>>::New(this->phases().size());
-    (
+    auto& Fs = *tFs;
-        new PtrList<volVectorField>(this->phases().size())
-    );
-    PtrList<volVectorField>& Fs = tFs();
 
     // Add the lift force
     forAllConstIters(liftModels_, modelIter)
@@ -570,11 +561,9 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::phiDs
     const PtrList<volScalarField>& rAUs
 ) const
 {
-    autoPtr<PtrList<surfaceScalarField>> tphiDs
+    auto tphiDs =
-    (
+        autoPtr<PtrList<surfaceScalarField>>::New(this->phases().size());
-        new PtrList<surfaceScalarField>(this->phases().size())
+    auto& phiDs = *tphiDs;
-    );
-    PtrList<surfaceScalarField>& phiDs = tphiDs();
 
     // Add the turbulent dispersion force
     forAllConstIters(turbulentDispersionModels_, turbulentDispersionModelIter)

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C

@@ -52,7 +52,7 @@ ThermalPhaseChangePhaseSystem
         }
 
         // Initially assume no mass transfer
-        iDmdt_.insert
+        iDmdt_.set
         (
             pair,
             new volScalarField
@@ -186,22 +186,16 @@ Foam::autoPtr<Foam::phaseSystem::massTransferTable>
 Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::massTransfer() const
 {
     // Create a mass transfer matrix for each species of each phase
-    autoPtr<phaseSystem::massTransferTable> eqnsPtr
+    auto eqnsPtr = autoPtr<phaseSystem::massTransferTable>::New();
-    (
+    auto& eqns = *eqnsPtr;
-        new phaseSystem::massTransferTable()
-    );
 
-    phaseSystem::massTransferTable& eqns = eqnsPtr();
+    for (const phaseModel& phase : this->phaseModels_)
-
-    forAll(this->phaseModels_, phasei)
     {
-        const phaseModel& phase = this->phaseModels_[phasei];
-
         const PtrList<volScalarField>& Yi = phase.Y();
 
         forAll(Yi, i)
         {
-            eqns.insert
+            eqns.set
             (
                 Yi[i].name(),
                 new fvScalarMatrix(Yi[i], dimMass/dimTime)

applications/test/HashPtrTable/Make/files

@@ -1,3 +1,3 @@
-Test-hashPtrTable.C
+Test-HashPtrTable.C
 
-EXE = $(FOAM_USER_APPBIN)/Test-hashPtrTable
+EXE = $(FOAM_USER_APPBIN)/Test-HashPtrTable

applications/test/HashPtrTable/Test-HashPtrTable.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2017-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,6 +27,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include <iostream>
+#include "autoPtr.H"
 #include "HashPtrTable.H"
 
 using namespace Foam;
@@ -42,7 +43,7 @@ void printTable(const HashPtrTable<T>& table)
         Info<< iter.key() << " = ";
         if (ptr)
         {
-            Info<< *ptr;
+            Info<< *ptr << " (" << long(ptr) << ")";
         }
         else
         {
@@ -77,11 +78,11 @@ void printTable(const HashPtrTable<T>& table)
 int main()
 {
     HashPtrTable<double> myTable;
-    myTable.insert("abc", new double(42.1));
+    myTable.set("abc", new double(42.1));
-    myTable.insert("def", nullptr);
+    myTable.set("def", nullptr);
-    myTable.insert("pi", new double(3.14159));
+    myTable.set("pi", new double(3.14159));
-    myTable.insert("natlog", new double(2.718282));
+    myTable.set("natlog", new double(2.718282));
-    myTable.insert("sqrt2", new double(1.414214));
+    myTable.insert("sqrt2", autoPtr<double>::New(1.414214));
 
     // Info<< myTable << endl;
     printTable(myTable);
@@ -105,6 +106,25 @@ int main()
 
     printTable(myTable);
 
+    HashPtrTable<double> moved(std::move(copy));
+
+    Info<< nl << "test movable" << nl;
+    Info<<"input:" << nl;
+    printTable(copy);
+
+    Info<<"output:" << nl;
+    printTable(moved);
+
+    HashPtrTable<double> other;
+
+    Info<<"move assign" << nl;
+
+    other = std::move(moved);
+    printTable(other);
+
+    Info<<"old" << nl;
+    printTable(moved);
+
     return 0;
 }
 

applications/test/PtrMap/Make/files

@@ -0,0 +1,3 @@
+Test-PtrMap.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-PtrMap

applications/test/PtrMap/Test-PtrMap.C

@@ -0,0 +1,130 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Description
+
+
+\*---------------------------------------------------------------------------*/
+
+#include <iostream>
+#include "PtrMap.H"
+
+using namespace Foam;
+
+template<class T>
+void printTable(const PtrMap<T>& table)
+{
+    Info<< table.size() << nl << "(" << nl;
+
+    forAllConstIters(table, iter)
+    {
+        const T* ptr = iter.object();
+        Info<< iter.key() << " = ";
+        if (ptr)
+        {
+            Info<< *ptr << " (" << long(ptr) << ")";
+        }
+        else
+        {
+            Info<< "nullptr";
+        }
+        Info<< nl;
+    }
+
+    Info<< ")" << endl;
+
+    // Values only, with for-range
+    Info<< "values (";
+    for (auto val : table)
+    {
+        Info<< ' ';
+        if (val)
+        {
+            Info<< *val;
+        }
+        else
+        {
+            Info<< "nullptr";
+        }
+    }
+    Info<< " )" << nl;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+//  Main program:
+
+int main()
+{
+    PtrMap<double> myTable;
+    myTable.set(1, new double(42.1));
+    myTable.set(2, nullptr);
+    myTable.set(3, new double(3.14159));
+    myTable.set(4, new double(2.718282));
+    myTable.set(4, new double(1.414214));
+
+    // Info<< myTable << endl;
+    printTable(myTable);
+
+    PtrMap<double> copy(myTable);
+
+    // Info<< copy << endl;
+    printTable(copy);
+    Info<< copy << endl;
+
+    Info<<"\nerase some existing and non-existing entries" << nl;
+
+    auto iter = myTable.find(3);
+    myTable.erase(iter);
+
+    iter = myTable.find(1000);  // unknown key
+    myTable.erase(iter);
+
+    myTable.erase(1);
+    iter = myTable.find(100000);  // unknown key
+
+    printTable(myTable);
+
+    PtrMap<double> moved(std::move(copy));
+
+    Info<< nl << "test movable" << nl;
+    Info<<"input:" << nl;
+    printTable(copy);
+
+    Info<<"output:" << nl;
+    printTable(moved);
+
+    PtrMap<double> other;
+
+    Info<<"move assign" << nl;
+
+    other = std::move(moved);
+    printTable(other);
+
+    Info<<"old" << nl;
+    printTable(moved);
+
+    return 0;
+}
+
+// ************************************************************************* //

applications/utilities/mesh/advanced/PDRMesh/PDRMesh.C

@@ -813,28 +813,33 @@ int main(int argc, char *argv[])
 
     // Read all fields in time and constant directories
     IOobjectList objects(mesh, runTime.timeName());
-    IOobjectList timeObjects(IOobjectList(mesh, mesh.facesInstance()));
-    forAllConstIter(IOobjectList, timeObjects, iter)
     {
+        IOobjectList timeObjects(mesh, mesh.facesInstance());
+
+        // Transfer specific types
+        forAllIters(timeObjects, iter)
+        {
+            autoPtr<IOobject> objPtr(timeObjects.remove(iter));
+            const auto& obj = *objPtr;
+
             if
             (
-            iter()->headerClassName() == volScalarField::typeName
+                obj.headerClassName() == volScalarField::typeName
-         || iter()->headerClassName() == volVectorField::typeName
+             || obj.headerClassName() == volVectorField::typeName
-         || iter()->headerClassName() == volSphericalTensorField::typeName
+             || obj.headerClassName() == volSphericalTensorField::typeName
-         || iter()->headerClassName() == volTensorField::typeName
+             || obj.headerClassName() == volTensorField::typeName
-         || iter()->headerClassName() == volSymmTensorField::typeName
+             || obj.headerClassName() == volSymmTensorField::typeName
-         || iter()->headerClassName() == surfaceScalarField::typeName
+             || obj.headerClassName() == surfaceScalarField::typeName
-         || iter()->headerClassName() == surfaceVectorField::typeName
+             || obj.headerClassName() == surfaceVectorField::typeName
-         || iter()->headerClassName()
+             || obj.headerClassName() == surfaceSphericalTensorField::typeName
-            == surfaceSphericalTensorField::typeName
+             || obj.headerClassName() == surfaceSymmTensorField::typeName
-         || iter()->headerClassName() == surfaceSymmTensorField::typeName
+             || obj.headerClassName() == surfaceTensorField::typeName
-         || iter()->headerClassName() == surfaceTensorField::typeName
             )
             {
-            objects.add(*iter());
+                objects.add(objPtr);
+            }
         }
     }
-
     // Read vol fields and subset.
 
     wordList scalarNames(objects.names(volScalarField::typeName));

applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C

@@ -1082,7 +1082,7 @@ int main(int argc, char *argv[])
 
     if (!limitDict.empty())
     {
-        Info<< "Read refinement shells in = "
+        Info<< "Read limit shells in = "
             << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
     }
 

applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict

@@ -49,6 +49,14 @@ geometry
         max (3.5 2 0.5);
     }
 
+    // Shell for directional refinement
+    wakeBox
+    {
+        type searchableBox;
+        min (1.5 1 -0.5);
+        max (3.5 2 0.5);
+    }
+
     sphere.stl
     {
         type triSurfaceMesh;
@@ -268,10 +276,40 @@ castellatedMeshControls
         //    // All three settings can be overridden on a surface by
         //    // surface basis in the refinementSurfaces section.
         //    gapLevel (<numGapCells> <minLevel> <maxlevel>);
+        //    // Optional: when doing the gapLevel refinement directly remove
+        //    //           based on orientation w.r.t. gap. This limits the
+        //    //           amount of cells before doing the 'locationInMesh'
+        //    //           cell selection. Default is 'mixed' i.e. keep cells
+        //    //           whilst doing the gap-level refinement.
+        //    //gapMode inside;  // inside/outside/mixed
+        //}
+
+        //wakeBox
+        //{
+        //    mode        inside;
+        //    // Dummy base level
+        //    levels      ((10000 0));
+        //
+        //    // Optional directional refinement (after all other refinement)
+        //    // Directional refinement
+        //    // for all cells according to 'mode' ('inside' or 'outside';
+        //    // 'distance' not supported) and within certain range. E.g.
+        //    //  - for all cells with level 2-5
+        //    //  - do one split in x direction
+        //    levelIncrement  (2 5 (1 0 0));
+        //
+        //    // Note
+        //    // - ignores 'levels' and gap* settings.
+        //    // - the cellLevel/pointLevels files are no longer consistent
+        //    //   with the mesh, the resulting mesh is no longer compatible
+        //    //   with e.g. dynamic refinement/unrefinement.
+        //    // - cellLevel will include any directional refinement
+        //    //   (i.e. it will be the maximum of all three directions)
         //}
     }
 
 
+
     // Optionally limit refinement in geometric region. This limits all
     // refinement (from features, refinementSurfaces, refinementRegions)
     // in a given geometric region. The syntax is exactly the same as for the

etc/caseDicts/postProcessing/fields/processorField

@@ -0,0 +1,19 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Web:      www.OpenFOAM.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+Description
+    Writes a volScalarField with processor number. Useful for postprocessing.
+
+\*---------------------------------------------------------------------------*/
+
+type            processorField;
+libs            ("libfieldFunctionObjects.so");
+
+executeControl  writeTime;
+writeControl    writeTime;
+
+// ************************************************************************* //

etc/codeTemplates/dynamicCode/codedPoints0MotionSolverTemplate.C

@@ -0,0 +1,115 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) YEAR AUTHOR,AFFILIATION
+     \\/     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 "codedPoints0MotionSolverTemplate.H"
+#include "fvCFD.H"
+#include "unitConversion.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(${typeName}Points0MotionSolver, 0);
+
+addRemovableToRunTimeSelectionTable
+(
+    motionSolver,
+    ${typeName}Points0MotionSolver,
+    dictionary
+);
+
+
+// * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
+
+extern "C"
+{
+    // dynamicCode:
+    // SHA1 = ${SHA1sum}
+    //
+    // unique function name that can be checked if the correct library version
+    // has been loaded
+    void ${typeName}_${SHA1sum}(bool load)
+    {
+        if (load)
+        {
+            // code that can be explicitly executed after loading
+        }
+        else
+        {
+            // code that can be explicitly executed before unloading
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+//{{{ begin localCode
+${localCode}
+//}}} end localCode
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+${typeName}Points0MotionSolver::${typeName}Points0MotionSolver
+(
+    const polyMesh& mesh,
+    const IOdictionary& dict
+)
+:
+    points0MotionSolver(mesh, dict, "${typeName}")
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+${typeName}Points0MotionSolver::~${typeName}Points0MotionSolver()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+tmp<pointField> ${typeName}Points0MotionSolver::curPoints() const
+{
+    if (${verbose:-false})
+    {
+        Info<<"curPoints ${typeName} sha1: ${SHA1sum}\n";
+    }
+
+//{{{ begin code
+    ${code}
+//}}} end code
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

etc/codeTemplates/dynamicCode/codedPoints0MotionSolverTemplate.H

@@ -0,0 +1,103 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) YEAR AUTHOR,AFFILIATION
+     \\/     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/>.
+
+Description
+    Template for use with dynamic code generation of a
+    points0MotionSolver.
+
+SourceFiles
+    codedPoints0MotionSolverTemplate.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef codedPoints0MotionSolverTemplate_H
+#define codedPoints0MotionSolverTemplate_H
+
+#include "points0MotionSolver.H"
+
+//{{{ begin codeInclude
+${codeInclude}
+//}}} end codeInclude
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                     A templated codedPoints0MotionSolver
+\*---------------------------------------------------------------------------*/
+
+class ${typeName}Points0MotionSolver
+:
+    public points0MotionSolver
+{
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        ${typeName}Points0MotionSolver(const ${typeName}Points0MotionSolver&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const ${typeName}Points0MotionSolver&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("${typeName}");
+
+
+    // Constructors
+
+        //- Construct from mesh and dictionary
+        ${typeName}Points0MotionSolver
+        (
+            const polyMesh& mesh,
+            const IOdictionary&
+        );
+
+
+    //- Destructor
+    virtual ~${typeName}Points0MotionSolver();
+
+
+    // Member Functions
+
+        //- Provide current points for motion. Uses current motion field
+        virtual tmp<pointField> curPoints() const;
+
+        //- Solve for motion
+        virtual void solve()
+        {};
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/Allwmake

@@ -82,6 +82,7 @@ wmake $targetType sixDoFRigidBodyState
 wmake $targetType rigidBodyDynamics
 wmake $targetType rigidBodyMeshMotion
 wmake $targetType semiPermeableBaffle
+wmake $targetType atmosphericModels
 
 # Needs access to Turbulence
 

src/OpenFOAM/Make/files

@@ -628,36 +628,36 @@ fields/cloud/cloud.C
 
 Fields = fields/Fields
 
+$(Fields)/Field/FieldBase.C
 $(Fields)/labelField/labelField.C
-$(Fields)/scalarField/scalarField.C
-$(Fields)/vectorField/vectorField.C
-$(Fields)/vector2DField/vector2DField.C
-$(Fields)/sphericalTensorField/sphericalTensorField.C
-$(Fields)/diagTensorField/diagTensorField.C
-$(Fields)/symmTensorField/symmTensorField.C
-$(Fields)/tensorField/tensorField.C
-$(Fields)/quaternionField/quaternionField.C
-$(Fields)/triadField/triadField.C
-$(Fields)/complexFields/complexFields.C
-
 $(Fields)/labelField/labelIOField.C
 $(Fields)/labelField/labelFieldIOField.C
+$(Fields)/scalarField/scalarField.C
 $(Fields)/scalarField/scalarIOField.C
 $(Fields)/scalarField/scalarFieldIOField.C
+$(Fields)/vectorField/vectorField.C
 $(Fields)/vectorField/vectorIOField.C
 $(Fields)/vectorField/vectorFieldIOField.C
+$(Fields)/vector2DField/vector2DField.C
 $(Fields)/vector2DField/vector2DIOField.C
 $(Fields)/vector2DField/vector2DFieldIOField.C
+$(Fields)/sphericalTensorField/sphericalTensorField.C
 $(Fields)/sphericalTensorField/sphericalTensorIOField.C
 $(Fields)/sphericalTensorField/sphericalTensorFieldIOField.C
+$(Fields)/diagTensorField/diagTensorField.C
 $(Fields)/diagTensorField/diagTensorIOField.C
 $(Fields)/diagTensorField/diagTensorFieldIOField.C
+$(Fields)/symmTensorField/symmTensorField.C
 $(Fields)/symmTensorField/symmTensorIOField.C
 $(Fields)/symmTensorField/symmTensorFieldIOField.C
+$(Fields)/tensorField/tensorField.C
 $(Fields)/tensorField/tensorIOField.C
 $(Fields)/tensorField/tensorFieldIOField.C
+$(Fields)/quaternionField/quaternionField.C
 $(Fields)/quaternionField/quaternionIOField.C
+$(Fields)/triadField/triadField.C
 $(Fields)/triadField/triadIOField.C
+$(Fields)/complexFields/complexFields.C
 $(Fields)/transformField/transformField.C
 
 

src/OpenFOAM/containers/HashTables/HashPtrTable/HashPtrTable.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2017-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,20 +28,6 @@ License
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-template<class T, class Key, class Hash>
-Foam::HashPtrTable<T, Key, Hash>::HashPtrTable()
-:
-    parent_type()
-{}
-
-
-template<class T, class Key, class Hash>
-Foam::HashPtrTable<T, Key, Hash>::HashPtrTable(const label size)
-:
-    parent_type(size)
-{}
-
-
 template<class T, class Key, class Hash>
 Foam::HashPtrTable<T, Key, Hash>::HashPtrTable
 (
@@ -55,16 +41,26 @@ Foam::HashPtrTable<T, Key, Hash>::HashPtrTable
         const T* ptr = iter.object();
         if (ptr)
         {
-            this->insert(iter.key(), new T(*ptr));
+            this->set(iter.key(), new T(*ptr));
         }
         else
         {
-            this->insert(iter.key(), nullptr);
+            this->set(iter.key(), nullptr);
         }
     }
 }
 
 
+template<class T, class Key, class Hash>
+Foam::HashPtrTable<T, Key, Hash>::HashPtrTable
+(
+    HashPtrTable<T, Key, Hash>&& ht
+)
+:
+    parent_type(std::move(ht))
+{}
+
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 template<class T, class Key, class Hash>
@@ -77,13 +73,13 @@ Foam::HashPtrTable<T, Key, Hash>::~HashPtrTable()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class T, class Key, class Hash>
-T* Foam::HashPtrTable<T, Key, Hash>::remove(iterator& iter)
+Foam::autoPtr<T> Foam::HashPtrTable<T, Key, Hash>::remove(iterator& iter)
 {
     if (iter.found())
     {
-        T* ptr = iter.object();
+        autoPtr<T> aptr(iter.object());
         this->parent_type::erase(iter);
-        return ptr;
+        return aptr;
     }
 
     return nullptr;
@@ -140,11 +136,10 @@ void Foam::HashPtrTable<T, Key, Hash>::operator=
     const HashPtrTable<T, Key, Hash>& rhs
 )
 {
-    // Check for assignment to self
     if (this == &rhs)
     {
         FatalErrorInFunction
-            << "attempted assignment to self"
+            << "attempted copy assignment to self"
             << abort(FatalError);
     }
 
@@ -155,16 +150,34 @@ void Foam::HashPtrTable<T, Key, Hash>::operator=
         const T* ptr = iter.object();
         if (ptr)
         {
-            this->insert(iter.key(), new T(*ptr));
+            this->set(iter.key(), new T(*ptr));
         }
         else
         {
-            this->insert(iter.key(), nullptr);
+            this->set(iter.key(), nullptr);
         }
     }
 }
 
 
+template<class T, class Key, class Hash>
+void Foam::HashPtrTable<T, Key, Hash>::operator=
+(
+    HashPtrTable<T, Key, Hash>&& rhs
+)
+{
+    if (this == &rhs)
+    {
+        FatalErrorInFunction
+            << "attempted move assignment to self"
+            << abort(FatalError);
+    }
+
+    this->clear();
+    this->transfer(rhs);
+}
+
+
 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 
 #include "HashPtrTableIO.C"

src/OpenFOAM/containers/HashTables/HashPtrTable/HashPtrTable.H

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2017-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -43,11 +43,12 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declarations
+
 class Istream;
 class Ostream;
 
-// Forward declaration of friend functions and operators
+template<class T> class autoPtr;
-
 template<class T, class Key, class Hash> class HashPtrTable;
 
 template<class T, class Key, class Hash>
@@ -92,10 +93,10 @@ public:
     // Constructors
 
         //- Construct null with default table capacity
-        HashPtrTable();
+        inline HashPtrTable();
 
         //- Construct given initial table capacity
-        explicit HashPtrTable(const label size);
+        explicit inline HashPtrTable(const label size);
 
         //- Construct from Istream using given Istream constructor class
         template<class INew>
@@ -107,9 +108,12 @@ public:
         //- Construct from dictionary with default dictionary constructor class
         explicit HashPtrTable(const dictionary& dict);
 
-        //- Construct as copy
+        //- Copy construct
         HashPtrTable(const this_type& ht);
 
+        //- Move construct
+        HashPtrTable(this_type&& ht);
+
 
     //- Destructor
     ~HashPtrTable();
@@ -121,16 +125,18 @@ public:
 
         //- Remove and return the pointer specified by given iterator.
         //  Includes a safeguard against the end-iterator.
-        T* remove(iterator& iter);
+        autoPtr<T> remove(iterator& iter);
 
-        //- Erase an entry specified by given iterator.
+        //- Erase an entry specified by given iterator and delete the
+        //- allocated pointer.
         //  Includes a safeguard against the end-iterator.
         bool erase(iterator& iter);
 
-        //- Erase an entry specified by the given key
+        //- Erase an entry specified by the given key and delete the
+        //- allocated pointer.
         bool erase(const Key& key);
 
-        //- Clear all entries from table and deleting any allocated pointers
+        //- Clear all entries from table and delete any allocated pointers
         void clear();
 
         //- Write
@@ -142,6 +148,9 @@ public:
         //- Copy assignment
         void operator=(const this_type& rhs);
 
+        //- Move assignment
+        void operator=(this_type&& rhs);
+
 
     // IOstream Operators
 
@@ -156,6 +165,31 @@ public:
             Ostream& os,
             const HashPtrTable<T, Key, Hash>& tbl
         );
+
+
+    // Housekeeping
+
+        //- No insert() with raw pointers (potential memory leaks)
+        //  Use insert() with autoPtr or set()
+        inline bool insert(const Key&, T*) = delete;
+
+        //- Insert a new entry, not overwriting existing entries.
+        //  \return True if the entry inserted, which means that it did
+        //  not previously exist in the table.
+        inline bool insert(const Key& key, autoPtr<T>& aptr);
+
+        //- Insert a new entry, not overwriting existing entries.
+        inline bool insert(const Key& key, autoPtr<T>&& aptr);
+
+        //- Use set(), not insert() to avoid potential memory leaks
+        //- Assign a new entry, overwriting existing entries.
+        inline bool set(const Key& key, T* ptr);
+
+        //- Assign a new entry, overwriting existing entries.
+        inline bool set(const Key& key, autoPtr<T>& aptr);
+
+        //- Assign a new entry, overwriting existing entries.
+        inline bool set(const Key& key, autoPtr<T>&& aptr);
 };
 
 
@@ -165,6 +199,10 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#include "HashPtrTableI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #ifdef NoRepository
     #include "HashPtrTable.C"
 #endif

src/OpenFOAM/containers/HashTables/HashPtrTable/HashPtrTableI.H

@@ -0,0 +1,113 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "autoPtr.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class T, class Key, class Hash>
+inline Foam::HashPtrTable<T, Key, Hash>::HashPtrTable()
+:
+    parent_type()
+{}
+
+
+template<class T, class Key, class Hash>
+inline Foam::HashPtrTable<T, Key, Hash>::HashPtrTable(const label size)
+:
+    parent_type(size)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+template<class T, class Key, class Hash>
+inline bool Foam::HashPtrTable<T, Key, Hash>::insert
+(
+    const Key& key,
+    autoPtr<T>& aptr
+)
+{
+    if (parent_type::insert(key, aptr.get()))
+    {
+        aptr.release(); // Now owned by HashPtrTable
+        return true;
+    }
+
+    return false;
+}
+
+
+template<class T, class Key, class Hash>
+inline bool Foam::HashPtrTable<T, Key, Hash>::insert
+(
+    const Key& key,
+    autoPtr<T>&& aptr
+)
+{
+    if (parent_type::insert(key, aptr.get()))
+    {
+        aptr.release(); // Now owned by HashPtrTable
+        return true;
+    }
+
+    return false;
+}
+
+
+template<class T, class Key, class Hash>
+inline bool Foam::HashPtrTable<T, Key, Hash>::set
+(
+    const Key& key,
+    T* ptr
+)
+{
+    return this->parent_type::set(key, ptr);
+}
+
+
+template<class T, class Key, class Hash>
+inline bool Foam::HashPtrTable<T, Key, Hash>::set
+(
+    const Key& key,
+    autoPtr<T>& aptr
+)
+{
+    return this->set(key, aptr.release());
+}
+
+
+template<class T, class Key, class Hash>
+inline bool Foam::HashPtrTable<T, Key, Hash>::set
+(
+    const Key& key,
+    autoPtr<T>&& aptr
+)
+{
+    return this->set(key, aptr.release());
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/containers/HashTables/HashPtrTable/HashPtrTableIO.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2017-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -65,7 +65,7 @@ void Foam::HashPtrTable<T, Key, Hash>::read(Istream& is, const INew& inew)
                 {
                     Key key;
                     is >> key;
-                    this->insert(key, inew(key, is).ptr());
+                    this->set(key, inew(key, is).ptr());
 
                     is.fatalCheck
                     (
@@ -110,7 +110,7 @@ void Foam::HashPtrTable<T, Key, Hash>::read(Istream& is, const INew& inew)
             is.putBack(lastToken);
             Key key;
             is >> key;
-            this->insert(key, inew(key, is).ptr());
+            this->set(key, inew(key, is).ptr());
 
             is.fatalCheck
             (
@@ -147,7 +147,7 @@ void Foam::HashPtrTable<T, Key, Hash>::read
     {
         const word& k = iter().keyword();
 
-        this->insert(k, inew(dict.subDict(k)).ptr());
+        this->set(k, inew(dict.subDict(k)).ptr());
     }
 }
 

src/OpenFOAM/containers/HashTables/HashTable/HashTable.C

@@ -298,7 +298,7 @@ Foam::label Foam::HashTable<T, Key, Hash>::countEntries
 
 
 template<class T, class Key, class Hash>
-bool Foam::HashTable<T, Key, Hash>::set
+bool Foam::HashTable<T, Key, Hash>::setEntry
 (
     const Key& key,
     const T& obj,

src/OpenFOAM/containers/HashTables/HashTable/HashTable.H

@@ -165,8 +165,10 @@ class HashTable
             }
 
         private:
-            //- Disallow default bitwise copy construct / assignment
+            //- No copy construct
             pair_entry(const pair_entry&) = delete;
+
+            //- No copy assignment
             void operator=(const pair_entry&) = delete;
         };
 
@@ -214,8 +216,10 @@ class HashTable
 
 
         private:
-            //- Disallow default bitwise copy construct / assignment
+            //- No copy construct
             unary_entry(const unary_entry&) = delete;
+
+            //- No copy assignment
             void operator=(const unary_entry&) = delete;
         };
 
@@ -249,7 +253,7 @@ class HashTable
 
         //- Assign a new hash-entry to a possibly already existing key.
         //  \return True if the new entry was set.
-        bool set(const Key& key, const T& obj, const bool overwrite);
+        bool setEntry(const Key& key, const T& obj, const bool overwrite);
 
 
 public:
@@ -606,10 +610,10 @@ public:
         //- Return existing entry or insert a new entry.
         inline T& operator()(const Key& key, const T& deflt);
 
-        //- Copy assignment
+        //- Copy assign
         void operator=(const HashTable<T, Key, Hash>& rhs);
 
-        //- Copy assignment from an initializer list
+        //- Copy assign from an initializer list
         void operator=(std::initializer_list<std::pair<Key, T>> rhs);
 
         //- Move assign

src/OpenFOAM/containers/HashTables/HashTable/HashTableI.H

@@ -121,7 +121,7 @@ inline bool Foam::HashTable<T, Key, Hash>::insert
     const T& obj
 )
 {
-    return this->set(key, obj, false);  // No overwrite
+    return this->setEntry(key, obj, false);  // No overwrite
 }
 
 
@@ -132,7 +132,7 @@ inline bool Foam::HashTable<T, Key, Hash>::set
     const T& obj
 )
 {
-    return this->set(key, obj, true);   // Overwrite
+    return this->setEntry(key, obj, true);   // Overwrite
 }
 
 

src/OpenFOAM/containers/HashTables/PtrMap/PtrMap.H

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2017-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -80,11 +80,32 @@ public:
             parent_type(is)
         {}
 
-        //- Construct as copy
+        //- Copy construct
         PtrMap(const this_type& map)
         :
             parent_type(map)
         {}
+
+        //- Move construct
+        PtrMap(this_type&& map)
+        :
+            parent_type(std::move(map))
+        {}
+
+
+    // Member Operators
+
+        //- Copy assignment
+        void operator=(const this_type& rhs)
+        {
+            parent_type::operator=(rhs);
+        }
+
+        //- Move assignment
+        void operator=(this_type&& rhs)
+        {
+            parent_type::operator=(std::move(rhs));
+        }
 };
 
 

src/OpenFOAM/db/IOobjectList/IOobjectList.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2016-2017 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2016-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -52,7 +52,7 @@ namespace Foam
                     InfoInFunction << "Found " << iter.key() << endl;
                 }
 
-                results.insert
+                results.set
                 (
                     iter.key(),
                     new IOobject(*(iter.object()))
@@ -123,9 +123,27 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::IOobjectList::IOobjectList(const label nIoObjects)
+Foam::IOobjectList::IOobjectList()
 :
-    HashPtrTable<IOobject>(nIoObjects)
+    HashPtrTable<IOobject>()
+{}
+
+
+Foam::IOobjectList::IOobjectList(const label nObjects)
+:
+    HashPtrTable<IOobject>(nObjects)  // Could also use 2*nObjects instead
+{}
+
+
+Foam::IOobjectList::IOobjectList(const IOobjectList& list)
+:
+    HashPtrTable<IOobject>(list)
+{}
+
+
+Foam::IOobjectList::IOobjectList(IOobjectList&& list)
+:
+    HashPtrTable<IOobject>(std::move(list))
 {}
 
 
@@ -142,7 +160,7 @@ Foam::IOobjectList::IOobjectList
     HashPtrTable<IOobject>()
 {
     word newInstance;
-    fileNameList ObjectNames = fileHandler().readObjects
+    fileNameList objNames = fileHandler().readObjects
     (
         db,
         instance,
@@ -150,9 +168,9 @@ Foam::IOobjectList::IOobjectList
         newInstance
     );
 
-    for (const auto& objName : ObjectNames)
+    for (const auto& objName : objNames)
     {
-        IOobject* objectPtr = new IOobject
+        auto objectPtr = autoPtr<IOobject>::New
         (
             objName,
             newInstance,
@@ -183,35 +201,35 @@ Foam::IOobjectList::IOobjectList
         {
             insert(objectPtr->name(), objectPtr);
         }
-        else
-        {
-            delete objectPtr;
     }
 }
-}
-
-
-Foam::IOobjectList::IOobjectList(const IOobjectList& iolist)
-:
-    HashPtrTable<IOobject>(iolist)
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::IOobjectList::~IOobjectList()
-{}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-bool Foam::IOobjectList::add(IOobject& io)
+bool Foam::IOobjectList::add(autoPtr<IOobject>& objectPtr)
 {
-    return insert(io.name(), &io);
+    if (objectPtr.valid())
+    {
+        return insert(objectPtr->name(), objectPtr);
+    }
+
+    return false;
 }
 
 
-bool Foam::IOobjectList::remove(IOobject& io)
+bool Foam::IOobjectList::add(autoPtr<IOobject>&& objectPtr)
+{
+    if (objectPtr.valid())
+    {
+        return insert(objectPtr->name(), objectPtr);
+    }
+
+    return false;
+}
+
+
+bool Foam::IOobjectList::remove(const IOobject& io)
 {
     return erase(io.name());
 }
@@ -219,7 +237,7 @@ bool Foam::IOobjectList::remove(IOobject& io)
 
 Foam::IOobject* Foam::IOobjectList::lookup(const word& name) const
 {
-    const_iterator iter = find(name);
+    const_iterator iter = cfind(name);
 
     if (iter.found())
     {
@@ -267,7 +285,7 @@ Foam::IOobjectList Foam::IOobjectList::lookupClass(const word& clsName) const
                 InfoInFunction << "Found " << iter.key() << endl;
             }
 
-            results.insert
+            results.set
             (
                 iter.key(),
                 new IOobject(*(iter.object()))
@@ -369,6 +387,14 @@ Foam::wordList Foam::IOobjectList::sortedNames
 }
 
 
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+void Foam::IOobjectList::operator=(IOobjectList&& list)
+{
+    transfer(list);
+}
+
+
 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
 
 Foam::Ostream& Foam::operator<<(Ostream& os, const IOobjectList& list)

src/OpenFOAM/db/IOobjectList/IOobjectList.H

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2016-2017 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2016-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -45,7 +45,7 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declaration of friend functions and operators
+// Forward declarations
 class IOobjectList;
 Ostream& operator<<(Ostream& os, const IOobjectList& list);
 
@@ -58,18 +58,21 @@ class IOobjectList
 :
     public HashPtrTable<IOobject>
 {
-    // Private Member Functions
-
-        //- No copy assignment
-        void operator=(const IOobjectList&) = delete;
-
-
 public:
 
     // Constructors
 
-        //- Construct given an initial estimate for the number of entries
+        //- Construct null with default (128) table capacity
-        explicit IOobjectList(const label nIoObjects = 128);
+        IOobjectList();
+
+        //- Construct given initial table capacity
+        explicit IOobjectList(const label nObjects);
+
+        //- Copy construct
+        IOobjectList(const IOobjectList& list);
+
+        //- Move construct
+        IOobjectList(IOobjectList&& list);
 
         //- Construct from objectRegistry and instance path
         IOobjectList
@@ -82,12 +85,9 @@ public:
             bool registerObject = true
         );
 
-        //- Construct as copy
-        IOobjectList(const IOobjectList& iolist);
-
 
     //- Destructor
-    ~IOobjectList();
+    ~IOobjectList() = default;
 
 
     // Member functions
@@ -95,10 +95,16 @@ public:
     // Basic methods
 
         //- Add an IOobject to the list
-        bool add(IOobject& io);
+        bool add(autoPtr<IOobject>& objectPtr);
+
+        //- Add an IOobject to the list
+        bool add(autoPtr<IOobject>&& objectPtr);
+
+        //- Remove an IOobject from the list, by iterator
+        using HashPtrTable<IOobject>::remove;
 
         //- Remove an IOobject from the list
-        bool remove(IOobject& io);
+        bool remove(const IOobject& io);
 
 
     // Lookup
@@ -230,6 +236,15 @@ public:
         wordList sortedNames(const word& clsName, const wordRes& matcher) const;
 
 
+    // Member Operators
+
+        //- No copy assignment
+        void operator=(const IOobjectList&) = delete;
+
+        //- Move assignment
+        void operator=(IOobjectList&& list);
+
+
     // Ostream Operator
 
         friend Ostream& operator<<(Ostream& os, const IOobjectList& list);

src/OpenFOAM/fields/DimensionedFields/DimensionedField/SubDimensionedFieldI.H

@@ -63,7 +63,7 @@ inline Foam::SubDimensionedField<Type, GeoMesh>::SubDimensionedField
     const SubDimensionedField<Type, GeoMesh>& sfield
 )
 :
-    tmp<SubDimensionedField<Type, GeoMesh>>::refCount(),
+    refCount(),
     SubField<Type>(sfield)
 {}
 

src/OpenFOAM/fields/FieldFields/FieldField/FieldField.C

@@ -21,9 +21,6 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
-Description
-    Generic fieldField type.
-
 \*---------------------------------------------------------------------------*/
 
 #include "FieldField.H"
@@ -137,7 +134,7 @@ FieldField<Field, Type>::FieldField
 template<template<class> class Field, class Type>
 FieldField<Field, Type>::FieldField(const FieldField<Field, Type>& f)
 :
-    tmp<FieldField<Field, Type>>::refCount(),
+    refCount(),
     PtrList<Field<Type>>(f)
 {}
 

src/OpenFOAM/fields/FieldFields/FieldField/FieldField.H

@@ -73,7 +73,7 @@ Ostream& operator<<
 template<template<class> class Field, class Type>
 class FieldField
 :
-    public tmp<FieldField<Field, Type>>::refCount,
+    public refCount,
     public PtrList<Field<Type>>
 {
 

src/OpenFOAM/fields/Fields/Field/Field.C

@@ -29,12 +29,6 @@ License
 #include "contiguous.H"
 #include "mapDistributeBase.H"
 
-// * * * * * * * * * * * * * * * Static Members  * * * * * * * * * * * * * * //
-
-template<class Type>
-const char* const Foam::Field<Type>::typeName("Field");
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class Type>

src/OpenFOAM/fields/Fields/Field/Field.H

@@ -47,6 +47,7 @@ SourceFiles
 #include "VectorSpace.H"
 #include "scalarList.H"
 #include "labelList.H"
+#include "FieldBase.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -73,7 +74,7 @@ Ostream& operator<<(Ostream&, const tmp<Field<Type>>&);
 template<class Type>
 class Field
 :
-    public tmp<Field<Type>>::refCount,
+    public FieldBase,
     public List<Type>
 {
 
@@ -86,11 +87,6 @@ public:
     typedef SubField<Type> subField;
 
 
-    // Static data members
-
-        static const char* const typeName;
-
-
     // Static Member Functions
 
         //- Return nullObject reference field

src/OpenFOAM/fields/Fields/Field/FieldBase.C

@@ -0,0 +1,33 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "FieldBase.H"
+
+// * * * * * * * * * * * * * * * Static Members  * * * * * * * * * * * * * * //
+
+const char* const Foam::FieldBase::typeName("Field");
+
+
+// ************************************************************************* //

src/OpenFOAM/fields/Fields/Field/FieldBase.H

@@ -0,0 +1,77 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::FieldBase
+
+Description
+    Template invariant parts for Field
+
+SourceFiles
+    FieldBase.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef FieldBase_H
+#define FieldBase_H
+
+#include "refCount.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class FieldBase Declaration
+\*---------------------------------------------------------------------------*/
+
+class FieldBase
+:
+    public refCount
+{
+public:
+
+    // Static data members
+
+        //- Typename for Field
+        static const char* const typeName;
+
+
+    // Constructors
+
+        //- Construct null, with refCount zero
+        constexpr FieldBase() noexcept
+        :
+            refCount()
+        {}
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/fields/Fields/Field/SubField.H

@@ -59,7 +59,7 @@ template<class Type> class SubField;
 template<class Type>
 class SubField
 :
-    public tmp<SubField<Type>>::refCount,
+    public refCount,
     public SubList<Type>
 {
 

src/OpenFOAM/fields/Fields/Field/SubFieldI.H

@@ -74,7 +74,7 @@ inline Foam::SubField<Type>::SubField
     const SubField<Type>& sfield
 )
 :
-    tmp<SubField<Type>>::refCount(),
+    refCount(),
     SubList<Type>(sfield)
 {}
 

src/OpenFOAM/global/fileOperations/masterUncollatedFileOperation/masterUncollatedFileOperation.C

@@ -2319,7 +2319,7 @@ Foam::instantList Foam::fileOperations::masterUncollatedFileOperation::findTimes
 
         instantList* tPtr = new instantList(std::move(times));
 
-        times_.insert(directory, tPtr);
+        times_.set(directory, tPtr);
 
         if (debug)
         {

src/OpenFOAM/graph/graph.C

@@ -66,7 +66,7 @@ void Foam::graph::readCurves(Istream& is)
         y[i] = xyData[i].y_;
     }
 
-    insert
+    set
     (
         wordify(yName_),
         new curve(wordify(yName_), curve::curveStyle::CONTINUOUS, y)
@@ -105,7 +105,11 @@ Foam::graph::graph
     yName_(yName),
     x_(x)
 {
-    insert(wordify(yName), new curve(yName, curve::curveStyle::CONTINUOUS, y));
+    set
+    (
+        wordify(yName),
+        new curve(yName, curve::curveStyle::CONTINUOUS, y)
+    );
 }
 
 

src/OpenFOAM/meshes/polyMesh/polyMeshIO.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2015-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -207,6 +207,7 @@ Foam::polyMesh::readUpdateState Foam::polyMesh::readUpdate()
         if (boundaryChanged)
         {
             WarningInFunction
+                << "Number of patches has changed.  This may have "
                 << "unexpected consequences.  Proceed with care." << endl;
 
             boundary_.clear();

src/OpenFOAM/primitives/functions/Function1/Function1/Function1.C

@@ -38,18 +38,11 @@ Foam::Function1<Type>::Function1(const word& entryName)
 template<class Type>
 Foam::Function1<Type>::Function1(const Function1<Type>& de)
 :
-    tmp<Function1<Type>>::refCount(),
+    refCount(),
     name_(de.name_)
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class Type>
-Foam::Function1<Type>::~Function1()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type>
@@ -111,8 +104,8 @@ Foam::FieldFunction1<Function1Type>::value
     const scalarField& x
 ) const
 {
-    tmp<Field<Type>> tfld(new Field<Type>(x.size()));
+    auto tfld = tmp<Field<Type>>::New(x.size());
-    Field<Type>& fld = tfld.ref();
+    auto& fld = tfld.ref();
 
     forAll(x, i)
     {
@@ -152,8 +145,8 @@ Foam::FieldFunction1<Function1Type>::integrate
     const scalarField& x2
 ) const
 {
-    tmp<Field<Type>> tfld(new Field<Type>(x1.size()));
+    auto tfld = tmp<Field<Type>>::New(x1.size());
-    Field<Type>& fld = tfld.ref();
+    auto& fld = tfld.ref();
 
     forAll(x1, i)
     {

src/OpenFOAM/primitives/functions/Function1/Function1/Function1.H

@@ -61,11 +61,11 @@ template<class Type> Ostream& operator<<(Ostream&, const Function1<Type>&);
 template<class Type>
 class Function1
 :
-    public tmp<Function1<Type>>::refCount
+    public refCount
 {
     // Private Member Functions
 
-        //- Disallow default bitwise assignment
+        //- No copy assignment
         void operator=(const Function1<Type>&) = delete;
 
 
@@ -119,7 +119,7 @@ public:
 
 
     //- Destructor
-    virtual ~Function1();
+    virtual ~Function1() = default;
 
 
     // Member Functions
@@ -194,8 +194,7 @@ public:
 
 
     //- Destructor
-    virtual ~FieldFunction1()
+    virtual ~FieldFunction1() = default;
-    {}
 
 
     // Member Functions

src/TurbulenceModels/turbulenceModels/Make/files

@@ -38,7 +38,6 @@ $(nutWallFunctions)/nutWallFunction/nutWallFunctionFvPatchScalarField.C
 
 $(nutWallFunctions)/nutkWallFunction/nutkWallFunctionFvPatchScalarField.C
 $(nutWallFunctions)/nutkRoughWallFunction/nutkRoughWallFunctionFvPatchScalarField.C
-$(nutWallFunctions)/nutkAtmRoughWallFunction/nutkAtmRoughWallFunctionFvPatchScalarField.C
 
 $(nutWallFunctions)/nutUBlendedWallFunction/nutUBlendedWallFunctionFvPatchScalarField.C
 $(nutWallFunctions)/nutUWallFunction/nutUWallFunctionFvPatchScalarField.C
@@ -70,11 +69,5 @@ RASBCs = RAS/derivedFvPatchFields
 $(RASBCs)/turbulentMixingLengthDissipationRateInlet/turbulentMixingLengthDissipationRateInletFvPatchScalarField.C
 $(RASBCs)/turbulentMixingLengthFrequencyInlet/turbulentMixingLengthFrequencyInletFvPatchScalarField.C
 
-/* Atmospheric boundary layer BCs */
-$(RASBCs)/atmBoundaryLayer/atmBoundaryLayer.C
-$(RASBCs)/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.C
-$(RASBCs)/atmBoundaryLayerInletK/atmBoundaryLayerInletKFvPatchScalarField.C
-$(RASBCs)/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.C
-
 
 LIB = $(FOAM_LIBBIN)/libturbulenceModels

src/atmosphericModels/Make/files

@@ -0,0 +1,10 @@
+derivedFvPatchFields/atmBoundaryLayer/atmBoundaryLayer.C
+derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.C
+derivedFvPatchFields/atmBoundaryLayerInletK/atmBoundaryLayerInletKFvPatchScalarField.C
+derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.C
+derivedFvPatchFields/nutkAtmRoughWallFunction/nutkAtmRoughWallFunctionFvPatchScalarField.C
+
+atmosphericTurbulentTransportModels.C
+porosityModels/powerLawLopesdaCosta/powerLawLopesdaCosta.C
+
+LIB = $(FOAM_LIBBIN)/libatmosphericModels

src/atmosphericModels/Make/options

@@ -0,0 +1,18 @@
+EXE_INC = \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/surfMesh/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude
+
+LIB_LIBS = \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
+    -lincompressibleTransportModels \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lsurfMesh \
+    -lfvOptions

src/atmosphericModels/atmosphericTurbulentTransportModels.C

@@ -0,0 +1,35 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 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 "turbulentTransportModels.H"
+
+// -------------------------------------------------------------------------- //
+// RAS models
+// -------------------------------------------------------------------------- //
+
+#include "kEpsilonLopesdaCosta.H"
+makeRASModel(kEpsilonLopesdaCosta);
+
+// ************************************************************************* //

src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayer/atmBoundaryLayer.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2014-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2018 OpenFOAM Foundation
      \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -168,8 +168,7 @@ public:
         //- Construct from the coordinates field and dictionary
         atmBoundaryLayer(const vectorField& p, const dictionary&);
 
-        //- Construct by mapping given
+        //- Construct by mapping given atmBoundaryLayer onto a new patch
-        // atmBoundaryLayer onto a new patch
         atmBoundaryLayer
         (
             const atmBoundaryLayer&,

src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -43,7 +43,7 @@ atmBoundaryLayerInletEpsilonFvPatchScalarField
     const DimensionedField<scalar, volMesh>& iF
 )
 :
-    fixedValueFvPatchScalarField(p, iF),
+    inletOutletFvPatchScalarField(p, iF),
     atmBoundaryLayer()
 {}
 
@@ -56,10 +56,23 @@ atmBoundaryLayerInletEpsilonFvPatchScalarField
     const dictionary& dict
 )
 :
-    fixedValueFvPatchScalarField(p, iF, dict, false),
+    inletOutletFvPatchScalarField(p, iF),
     atmBoundaryLayer(patch().Cf(), dict)
 {
-    scalarField::operator=(epsilon(patch().Cf()));
+    phiName_ = dict.lookupOrDefault<word>("phi", "phi");
+
+    refValue() = epsilon(patch().Cf());
+    refGrad() = 0;
+    valueFraction() = 1;
+
+    if (dict.found("value"))
+    {
+        scalarField::operator=(scalarField("value", dict, p.size()));
+    }
+    else
+    {
+        scalarField::operator=(refValue());
+    }
 }
 
 
@@ -72,7 +85,7 @@ atmBoundaryLayerInletEpsilonFvPatchScalarField
     const fvPatchFieldMapper& mapper
 )
 :
-    fixedValueFvPatchScalarField(psf, p, iF, mapper),
+    inletOutletFvPatchScalarField(psf, p, iF, mapper),
     atmBoundaryLayer(psf, mapper)
 {}
 
@@ -84,7 +97,7 @@ atmBoundaryLayerInletEpsilonFvPatchScalarField
     const DimensionedField<scalar, volMesh>& iF
 )
 :
-    fixedValueFvPatchScalarField(psf, iF),
+    inletOutletFvPatchScalarField(psf, iF),
     atmBoundaryLayer(psf)
 {}
 
@@ -96,7 +109,7 @@ void atmBoundaryLayerInletEpsilonFvPatchScalarField::autoMap
     const fvPatchFieldMapper& m
 )
 {
-    fixedValueFvPatchScalarField::autoMap(m);
+    inletOutletFvPatchScalarField::autoMap(m);
     atmBoundaryLayer::autoMap(m);
 }
 
@@ -107,7 +120,7 @@ void atmBoundaryLayerInletEpsilonFvPatchScalarField::rmap
     const labelList& addr
 )
 {
-    fixedValueFvPatchScalarField::rmap(psf, addr);
+    inletOutletFvPatchScalarField::rmap(psf, addr);
 
     const atmBoundaryLayerInletEpsilonFvPatchScalarField& blpsf =
         refCast<const atmBoundaryLayerInletEpsilonFvPatchScalarField>(psf);
@@ -120,6 +133,7 @@ void atmBoundaryLayerInletEpsilonFvPatchScalarField::write(Ostream& os) const
 {
     fvPatchScalarField::write(os);
     atmBoundaryLayer::write(os);
+    os.writeEntryIfDifferent<word>("phi", "phi", phiName_);
     writeEntry("value", os);
 }
 

src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -60,7 +60,7 @@ SourceFiles
 #define atmBoundaryLayerInletEpsilonFvPatchScalarField_H
 
 #include "fvPatchFields.H"
-#include "fixedValueFvPatchFields.H"
+#include "inletOutletFvPatchFields.H"
 #include "atmBoundaryLayer.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -74,7 +74,7 @@ namespace Foam
 
 class atmBoundaryLayerInletEpsilonFvPatchScalarField
 :
-    public fixedValueFvPatchScalarField,
+    public inletOutletFvPatchScalarField,
     public atmBoundaryLayer
 {
 

src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletK/atmBoundaryLayerInletKFvPatchScalarField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -43,7 +43,7 @@ atmBoundaryLayerInletKFvPatchScalarField
     const DimensionedField<scalar, volMesh>& iF
 )
 :
-    fixedValueFvPatchScalarField(p, iF),
+    inletOutletFvPatchScalarField(p, iF),
     atmBoundaryLayer()
 {}
 
@@ -56,10 +56,23 @@ atmBoundaryLayerInletKFvPatchScalarField
     const dictionary& dict
 )
 :
-    fixedValueFvPatchScalarField(p, iF, dict, false),
+    inletOutletFvPatchScalarField(p, iF),
     atmBoundaryLayer(patch().Cf(), dict)
 {
-    scalarField::operator=(k(patch().Cf()));
+    phiName_ = dict.lookupOrDefault<word>("phi", "phi");
+
+    refValue() = k(patch().Cf());
+    refGrad() = 0;
+    valueFraction() = 1;
+
+    if (dict.found("value"))
+    {
+        scalarField::operator=(scalarField("value", dict, p.size()));
+    }
+    else
+    {
+        scalarField::operator=(refValue());
+    }
 }
 
 
@@ -72,7 +85,7 @@ atmBoundaryLayerInletKFvPatchScalarField
     const fvPatchFieldMapper& mapper
 )
 :
-    fixedValueFvPatchScalarField(psf, p, iF, mapper),
+    inletOutletFvPatchScalarField(psf, p, iF, mapper),
     atmBoundaryLayer(psf, mapper)
 {}
 
@@ -84,7 +97,7 @@ atmBoundaryLayerInletKFvPatchScalarField
     const DimensionedField<scalar, volMesh>& iF
 )
 :
-    fixedValueFvPatchScalarField(psf, iF),
+    inletOutletFvPatchScalarField(psf, iF),
     atmBoundaryLayer(psf)
 {}
 
@@ -96,7 +109,7 @@ void atmBoundaryLayerInletKFvPatchScalarField::autoMap
     const fvPatchFieldMapper& m
 )
 {
-    fixedValueFvPatchScalarField::autoMap(m);
+    inletOutletFvPatchScalarField::autoMap(m);
     atmBoundaryLayer::autoMap(m);
 }
 
@@ -107,7 +120,7 @@ void atmBoundaryLayerInletKFvPatchScalarField::rmap
     const labelList& addr
 )
 {
-    fixedValueFvPatchScalarField::rmap(psf, addr);
+    inletOutletFvPatchScalarField::rmap(psf, addr);
 
     const atmBoundaryLayerInletKFvPatchScalarField& blpsf =
         refCast<const atmBoundaryLayerInletKFvPatchScalarField>(psf);
@@ -120,6 +133,7 @@ void atmBoundaryLayerInletKFvPatchScalarField::write(Ostream& os) const
 {
     fvPatchScalarField::write(os);
     atmBoundaryLayer::write(os);
+    os.writeEntryIfDifferent<word>("phi", "phi", phiName_);
     writeEntry("value", os);
 }
 

src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletK/atmBoundaryLayerInletKFvPatchScalarField.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2014-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -60,7 +60,7 @@ SourceFiles
 #define atmBoundaryLayerInletKFvPatchScalarField_H
 
 #include "fvPatchFields.H"
-#include "fixedValueFvPatchFields.H"
+#include "inletOutletFvPatchFields.H"
 #include "atmBoundaryLayer.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -74,7 +74,7 @@ namespace Foam
 
 class atmBoundaryLayerInletKFvPatchScalarField
 :
-    public fixedValueFvPatchScalarField,
+    public inletOutletFvPatchScalarField,
     public atmBoundaryLayer
 {
 

src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -43,7 +43,7 @@ atmBoundaryLayerInletVelocityFvPatchVectorField
     const DimensionedField<vector, volMesh>& iF
 )
 :
-    fixedValueFvPatchVectorField(p, iF),
+    inletOutletFvPatchVectorField(p, iF),
     atmBoundaryLayer()
 {}
 
@@ -56,10 +56,23 @@ atmBoundaryLayerInletVelocityFvPatchVectorField
     const dictionary& dict
 )
 :
-    fixedValueFvPatchVectorField(p, iF, dict, false),
+    inletOutletFvPatchVectorField(p, iF),
     atmBoundaryLayer(patch().Cf(), dict)
 {
-    vectorField::operator=(U(patch().Cf()));
+    phiName_ = dict.lookupOrDefault<word>("phi", "phi");
+
+    refValue() = U(patch().Cf());
+    refGrad() = Zero;
+    valueFraction() = 1;
+
+    if (dict.found("value"))
+    {
+        vectorField::operator=(vectorField("value", dict, p.size()));
+    }
+    else
+    {
+        vectorField::operator=(refValue());
+    }
 }
 
 
@@ -72,7 +85,7 @@ atmBoundaryLayerInletVelocityFvPatchVectorField
     const fvPatchFieldMapper& mapper
 )
 :
-    fixedValueFvPatchVectorField(pvf, p, iF, mapper),
+    inletOutletFvPatchVectorField(pvf, p, iF, mapper),
     atmBoundaryLayer(pvf, mapper)
 {}
 
@@ -84,7 +97,7 @@ atmBoundaryLayerInletVelocityFvPatchVectorField
     const DimensionedField<vector, volMesh>& iF
 )
 :
-    fixedValueFvPatchVectorField(pvf, iF),
+    inletOutletFvPatchVectorField(pvf, iF),
     atmBoundaryLayer(pvf)
 {}
 
@@ -96,7 +109,7 @@ void atmBoundaryLayerInletVelocityFvPatchVectorField::autoMap
     const fvPatchFieldMapper& m
 )
 {
-    fixedValueFvPatchVectorField::autoMap(m);
+    inletOutletFvPatchVectorField::autoMap(m);
     atmBoundaryLayer::autoMap(m);
 }
 
@@ -107,7 +120,7 @@ void atmBoundaryLayerInletVelocityFvPatchVectorField::rmap
     const labelList& addr
 )
 {
-    fixedValueFvPatchVectorField::rmap(pvf, addr);
+    inletOutletFvPatchVectorField::rmap(pvf, addr);
 
     const atmBoundaryLayerInletVelocityFvPatchVectorField& blpvf =
         refCast<const atmBoundaryLayerInletVelocityFvPatchVectorField>(pvf);
@@ -120,6 +133,7 @@ void atmBoundaryLayerInletVelocityFvPatchVectorField::write(Ostream& os) const
 {
     fvPatchVectorField::write(os);
     atmBoundaryLayer::write(os);
+    os.writeEntryIfDifferent<word>("phi", "phi", phiName_);
     writeEntry("value", os);
 }
 

src/atmosphericModels/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -61,7 +61,7 @@ SourceFiles
 #define atmBoundaryLayerInletVelocityFvPatchVectorField_H
 
 #include "fvPatchFields.H"
-#include "fixedValueFvPatchFields.H"
+#include "inletOutletFvPatchFields.H"
 #include "atmBoundaryLayer.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -75,7 +75,7 @@ namespace Foam
 
 class atmBoundaryLayerInletVelocityFvPatchVectorField
 :
-    public fixedValueFvPatchVectorField,
+    public inletOutletFvPatchVectorField,
     public atmBoundaryLayer
 {
 

src/atmosphericModels/derivedFvPatchFields/nutkAtmRoughWallFunction/nutkAtmRoughWallFunctionFvPatchScalarField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

src/atmosphericModels/kEpsilonLopesdaCosta/kEpsilonLopesdaCosta.C

@@ -0,0 +1,477 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 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 "kEpsilonLopesdaCosta.H"
+#include "fvOptions.H"
+#include "explicitPorositySource.H"
+#include "bound.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace RASModels
+{
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+void kEpsilonLopesdaCosta<BasicTurbulenceModel>::setPorosityCoefficient
+(
+    volScalarField::Internal& C,
+    const porosityModels::powerLawLopesdaCosta& pm
+)
+{
+    if (pm.dict().found(C.name()))
+    {
+        const labelList& cellZoneIDs = pm.cellZoneIDs();
+
+        const scalar Cpm = readScalar(pm.dict().lookup(C.name()));
+
+        forAll(cellZoneIDs, zonei)
+        {
+            const labelList& cells =
+                this->mesh_.cellZones()[cellZoneIDs[zonei]];
+
+            forAll(cells, i)
+            {
+                const label celli = cells[i];
+                C[celli] = Cpm;
+            }
+        }
+    }
+}
+
+
+template<class BasicTurbulenceModel>
+void kEpsilonLopesdaCosta<BasicTurbulenceModel>::setCdSigma
+(
+    volScalarField::Internal& C,
+    const porosityModels::powerLawLopesdaCosta& pm
+)
+{
+    if (pm.dict().found(C.name()))
+    {
+        const labelList& cellZoneIDs = pm.cellZoneIDs();
+        const scalarField& Sigma = pm.Sigma();
+
+        const scalar Cpm = readScalar(pm.dict().lookup(C.name()));
+
+        forAll(cellZoneIDs, zonei)
+        {
+            const labelList& cells =
+                this->mesh_.cellZones()[cellZoneIDs[zonei]];
+
+            forAll(cells, i)
+            {
+                const label celli = cells[i];
+                C[celli] = Cpm*Sigma[celli];
+            }
+        }
+    }
+}
+
+
+template<class BasicTurbulenceModel>
+void kEpsilonLopesdaCosta<BasicTurbulenceModel>::setPorosityCoefficients()
+{
+    fv::options::optionList& fvOptions(fv::options::New(this->mesh_));
+
+    forAll(fvOptions, i)
+    {
+        if (isA<fv::explicitPorositySource>(fvOptions[i]))
+        {
+            const fv::explicitPorositySource& eps =
+                refCast<const fv::explicitPorositySource>(fvOptions[i]);
+
+            if (isA<porosityModels::powerLawLopesdaCosta>(eps.model()))
+            {
+                const porosityModels::powerLawLopesdaCosta& pm =
+                    refCast<const porosityModels::powerLawLopesdaCosta>
+                    (
+                        eps.model()
+                    );
+
+                setPorosityCoefficient(Cmu_, pm);
+                setPorosityCoefficient(C1_, pm);
+                setPorosityCoefficient(C2_, pm);
+                setPorosityCoefficient(sigmak_, pm);
+                setPorosityCoefficient(sigmaEps_, pm);
+
+                setCdSigma(CdSigma_, pm);
+                setPorosityCoefficient(betap_, pm);
+                setPorosityCoefficient(betad_, pm);
+                setPorosityCoefficient(C4_, pm);
+                setPorosityCoefficient(C5_, pm);
+            }
+        }
+    }
+}
+
+
+template<class BasicTurbulenceModel>
+void kEpsilonLopesdaCosta<BasicTurbulenceModel>::correctNut()
+{
+    this->nut_ = Cmu_*sqr(k_)/epsilon_;
+    this->nut_.correctBoundaryConditions();
+    fv::options::New(this->mesh_).correct(this->nut_);
+
+    BasicTurbulenceModel::correctNut();
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<fvScalarMatrix> kEpsilonLopesdaCosta<BasicTurbulenceModel>::kSource
+(
+    const volScalarField::Internal& magU,
+    const volScalarField::Internal& magU3
+) const
+{
+    return fvm::Su(CdSigma_*(betap_*magU3 - betad_*magU*k_()), k_);
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<fvScalarMatrix>
+kEpsilonLopesdaCosta<BasicTurbulenceModel>::epsilonSource
+(
+    const volScalarField::Internal& magU,
+    const volScalarField::Internal& magU3
+) const
+{
+    return fvm::Su
+    (
+        CdSigma_
+       *(C4_*betap_*epsilon_()/k_()*magU3 - C5_*betad_*magU*epsilon_()),
+        epsilon_
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+kEpsilonLopesdaCosta<BasicTurbulenceModel>::kEpsilonLopesdaCosta
+(
+    const alphaField& alpha,
+    const rhoField& rho,
+    const volVectorField& U,
+    const surfaceScalarField& alphaRhoPhi,
+    const surfaceScalarField& phi,
+    const transportModel& transport,
+    const word& propertiesName,
+    const word& type
+)
+:
+    eddyViscosity<RASModel<BasicTurbulenceModel>>
+    (
+        type,
+        alpha,
+        rho,
+        U,
+        alphaRhoPhi,
+        phi,
+        transport,
+        propertiesName
+    ),
+
+    Cmu_
+    (
+        IOobject
+        (
+            "Cmu",
+            this->runTime_.timeName(),
+            this->mesh_
+        ),
+        this->mesh_,
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Cmu",
+            this->coeffDict_,
+            0.09
+        )
+    ),
+    C1_
+    (
+        IOobject
+        (
+            "C1",
+            this->runTime_.timeName(),
+            this->mesh_
+        ),
+        this->mesh_,
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "C1",
+            this->coeffDict_,
+            1.44
+        )
+    ),
+    C2_
+    (
+        IOobject
+        (
+            "C2",
+            this->runTime_.timeName(),
+            this->mesh_
+        ),
+        this->mesh_,
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "C2",
+            this->coeffDict_,
+            1.92
+        )
+    ),
+    sigmak_
+    (
+        IOobject
+        (
+            "sigmak",
+            this->runTime_.timeName(),
+            this->mesh_
+        ),
+        this->mesh_,
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "sigmak",
+            this->coeffDict_,
+            1.0
+        )
+    ),
+    sigmaEps_
+    (
+        IOobject
+        (
+            "sigmaEps",
+            this->runTime_.timeName(),
+            this->mesh_
+        ),
+        this->mesh_,
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "sigmaEps",
+            this->coeffDict_,
+            1.3
+        )
+    ),
+
+    CdSigma_
+    (
+        IOobject
+        (
+            "CdSigma",
+            this->runTime_.timeName(),
+            this->mesh_
+        ),
+        this->mesh_,
+        dimensionedScalar("CdSigma", dimless/dimLength, 0)
+    ),
+    betap_
+    (
+        IOobject
+        (
+            "betap",
+            this->runTime_.timeName(),
+            this->mesh_
+        ),
+        this->mesh_,
+        dimensionedScalar("betap", dimless, 0)
+    ),
+    betad_
+    (
+        IOobject
+        (
+            "betad",
+            this->runTime_.timeName(),
+            this->mesh_
+        ),
+        this->mesh_,
+        dimensionedScalar("betad", dimless, 0)
+    ),
+    C4_
+    (
+        IOobject
+        (
+            "C4",
+            this->runTime_.timeName(),
+            this->mesh_
+        ),
+        this->mesh_,
+        dimensionedScalar("C4", dimless, 0)
+    ),
+    C5_
+    (
+        IOobject
+        (
+            "C5",
+            this->runTime_.timeName(),
+            this->mesh_
+        ),
+        this->mesh_,
+        dimensionedScalar("C5", dimless, 0)
+    ),
+
+    k_
+    (
+        IOobject
+        (
+            IOobject::groupName("k", alphaRhoPhi.group()),
+            this->runTime_.timeName(),
+            this->mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        this->mesh_
+    ),
+    epsilon_
+    (
+        IOobject
+        (
+            IOobject::groupName("epsilon", alphaRhoPhi.group()),
+            this->runTime_.timeName(),
+            this->mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        this->mesh_
+    )
+{
+    bound(k_, this->kMin_);
+    bound(epsilon_, this->epsilonMin_);
+
+    if (type == typeName)
+    {
+        this->printCoeffs(type);
+    }
+
+    setPorosityCoefficients();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+bool kEpsilonLopesdaCosta<BasicTurbulenceModel>::read()
+{
+    if (eddyViscosity<RASModel<BasicTurbulenceModel>>::read())
+    {
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+template<class BasicTurbulenceModel>
+void kEpsilonLopesdaCosta<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& nut = this->nut_;
+    fv::options& fvOptions(fv::options::New(this->mesh_));
+
+    eddyViscosity<RASModel<BasicTurbulenceModel>>::correct();
+
+    volScalarField::Internal divU
+    (
+        fvc::div(fvc::absolute(this->phi(), U))().v()
+    );
+
+    tmp<volTensorField> tgradU = fvc::grad(U);
+    volScalarField::Internal G
+    (
+        this->GName(),
+        nut.v()*(dev(twoSymm(tgradU().v())) && tgradU().v())
+    );
+    tgradU.clear();
+
+    // Update epsilon and G at the wall
+    epsilon_.boundaryFieldRef().updateCoeffs();
+
+    volScalarField::Internal magU(mag(U));
+    volScalarField::Internal magU3(pow3(magU));
+
+    // Dissipation equation
+    tmp<fvScalarMatrix> epsEqn
+    (
+        fvm::ddt(alpha, rho, epsilon_)
+      + fvm::div(alphaRhoPhi, epsilon_)
+      - fvm::laplacian(alpha*rho*DepsilonEff(), epsilon_)
+     ==
+        C1_*alpha()*rho()*G*epsilon_()/k_()
+      - fvm::SuSp(((2.0/3.0)*C1_)*alpha()*rho()*divU, epsilon_)
+      - fvm::Sp(C2_*alpha()*rho()*epsilon_()/k_(), epsilon_)
+      + epsilonSource(magU, magU3)
+      + fvOptions(alpha, rho, epsilon_)
+    );
+
+    epsEqn.ref().relax();
+    fvOptions.constrain(epsEqn.ref());
+    epsEqn.ref().boundaryManipulate(epsilon_.boundaryFieldRef());
+    solve(epsEqn);
+    fvOptions.correct(epsilon_);
+    bound(epsilon_, this->epsilonMin_);
+
+    // Turbulent kinetic energy equation
+    tmp<fvScalarMatrix> kEqn
+    (
+        fvm::ddt(alpha, rho, k_)
+      + fvm::div(alphaRhoPhi, k_)
+      - fvm::laplacian(alpha*rho*DkEff(), k_)
+     ==
+        alpha()*rho()*G
+      - fvm::SuSp((2.0/3.0)*alpha()*rho()*divU, k_)
+      - fvm::Sp(alpha()*rho()*epsilon_()/k_(), k_)
+      + kSource(magU, magU3)
+      + fvOptions(alpha, rho, k_)
+    );
+
+    kEqn.ref().relax();
+    fvOptions.constrain(kEqn.ref());
+    solve(kEqn);
+    fvOptions.correct(k_);
+    bound(k_, this->kMin_);
+
+    correctNut();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace Foam
+
+// ************************************************************************* //

src/atmosphericModels/kEpsilonLopesdaCosta/kEpsilonLopesdaCosta.H

@@ -0,0 +1,245 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 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::RASModels::kEpsilonLopesdaCosta
+
+Group
+    grpRASTurbulence
+
+Description
+    Variant of the standard k-epsilon turbulence model with additional source
+    terms to handle the changes in turbulence in porous regions represented by
+    the powerLawLopesdaCosta porosity model.
+
+    Reference:
+    \verbatim
+        Costa, J. C. P. L. D. (2007).
+        Atmospheric flow over forested and non-forested complex terrain.
+    \endverbatim
+
+    The default model coefficients are
+    \verbatim
+        kEpsilonLopesdaCostaCoeffs
+        {
+            Cmu         0.09;
+            C1          1.44;
+            C2          1.92;
+            sigmak      1.0;
+            sigmaEps    1.3;
+        }
+    \endverbatim
+
+See also
+    Foam::RASModels::kEpsilon
+    Foam::porosityModels::powerLawLopesdaCosta
+
+SourceFiles
+    kEpsilonLopesdaCosta.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef kEpsilonLopesdaCosta_H
+#define kEpsilonLopesdaCosta_H
+
+#include "RASModel.H"
+#include "eddyViscosity.H"
+#include "powerLawLopesdaCosta.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace RASModels
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class kEpsilonLopesdaCosta Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class BasicTurbulenceModel>
+class kEpsilonLopesdaCosta
+:
+    public eddyViscosity<RASModel<BasicTurbulenceModel>>
+{
+    // Private Member Functions
+
+        // Disallow default bitwise copy construct and assignment
+        kEpsilonLopesdaCosta(const kEpsilonLopesdaCosta&);
+        void operator=(const kEpsilonLopesdaCosta&);
+
+
+protected:
+
+    // Protected data
+
+        // Standard model coefficients
+
+            volScalarField Cmu_;
+            volScalarField::Internal C1_;
+            volScalarField::Internal C2_;
+            volScalarField sigmak_;
+            volScalarField sigmaEps_;
+
+        // Lopes da Costa porosity coefficients
+
+            volScalarField::Internal CdSigma_;
+            volScalarField::Internal betap_;
+            volScalarField::Internal betad_;
+            volScalarField::Internal C4_;
+            volScalarField::Internal C5_;
+
+
+        // Fields
+
+            volScalarField k_;
+            volScalarField epsilon_;
+
+
+    // Protected Member Functions
+
+        void setPorosityCoefficient
+        (
+            volScalarField::Internal& C,
+            const porosityModels::powerLawLopesdaCosta& pm
+        );
+
+        void setCdSigma
+        (
+            volScalarField::Internal& C,
+            const porosityModels::powerLawLopesdaCosta& pm
+        );
+
+        void setPorosityCoefficients();
+
+        virtual void correctNut();
+
+        virtual tmp<fvScalarMatrix> kSource
+        (
+            const volScalarField::Internal& magU,
+            const volScalarField::Internal& magU3
+        ) const;
+
+        virtual tmp<fvScalarMatrix> epsilonSource
+        (
+            const volScalarField::Internal& magU,
+            const volScalarField::Internal& magU3
+        ) const;
+
+
+public:
+
+    typedef typename BasicTurbulenceModel::alphaField alphaField;
+    typedef typename BasicTurbulenceModel::rhoField rhoField;
+    typedef typename BasicTurbulenceModel::transportModel transportModel;
+
+
+    //- Runtime type information
+    TypeName("kEpsilonLopesdaCosta");
+
+
+    // Constructors
+
+        //- Construct from components
+        kEpsilonLopesdaCosta
+        (
+            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 ~kEpsilonLopesdaCosta()
+    {}
+
+
+    // Member Functions
+
+        //- Re-read model coefficients if they have changed
+        virtual bool read();
+
+        //- Return the effective diffusivity for k
+        tmp<volScalarField> DkEff() const
+        {
+            return tmp<volScalarField>
+            (
+                new volScalarField
+                (
+                    "DkEff",
+                    (this->nut_/sigmak_ + this->nu())
+                )
+            );
+        }
+
+        //- Return the effective diffusivity for epsilon
+        tmp<volScalarField> DepsilonEff() const
+        {
+            return tmp<volScalarField>
+            (
+                new volScalarField
+                (
+                    "DepsilonEff",
+                    (this->nut_/sigmaEps_ + this->nu())
+                )
+            );
+        }
+
+        //- Return the turbulence kinetic energy
+        virtual tmp<volScalarField> k() const
+        {
+            return k_;
+        }
+
+        //- Return the turbulence kinetic energy dissipation rate
+        virtual tmp<volScalarField> epsilon() const
+        {
+            return epsilon_;
+        }
+
+        //- Solve the turbulence equations and correct the turbulence viscosity
+        virtual void correct();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "kEpsilonLopesdaCosta.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/atmosphericModels/porosityModels/powerLawLopesdaCosta/powerLawLopesdaCosta.C

@@ -0,0 +1,419 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 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 "addToRunTimeSelectionTable.H"
+#include "powerLawLopesdaCosta.H"
+#include "geometricOneField.H"
+#include "fvMatrices.H"
+#include "triSurfaceMesh.H"
+#include "triSurfaceTools.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    namespace porosityModels
+    {
+        defineTypeNameAndDebug(powerLawLopesdaCosta, 0);
+        addToRunTimeSelectionTable(porosityModel, powerLawLopesdaCosta, mesh);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::porosityModels::powerLawLopesdaCostaZone::powerLawLopesdaCostaZone
+(
+    const word& name,
+    const word& modelType,
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+:
+    zoneName_(name + ":porousZone")
+{
+    dictionary coeffs(dict.optionalSubDict(modelType + "Coeffs"));
+
+    // Vertical direction within porous region
+    vector zDir(coeffs.lookup("zDir"));
+
+    // Span of the search for the cells in the porous region
+    scalar searchSpan(readScalar(coeffs.lookup("searchSpan")));
+
+    // Top surface file name defining the extent of the porous region
+    word topSurfaceFileName(coeffs.lookup("topSurface"));
+
+    // List of the ground patches defining the lower surface
+    // of the porous region
+    List<wordRe> groundPatches(coeffs.lookup("groundPatches"));
+
+    // Functional form of the porosity surface area per unit volume as a
+    // function of the normalized vertical position
+    autoPtr<Function1<scalar >> SigmaFunc
+    (
+        Function1<scalar>::New("Sigma", dict)
+    );
+
+    // Searchable triSurface for the top of the porous region
+    triSurfaceMesh searchSurf
+    (
+        IOobject
+        (
+            topSurfaceFileName,
+            mesh.time().constant(),
+            "triSurface",
+            mesh.time()
+        )
+    );
+
+    // Limit the porous cell search to those within the lateral and vertical
+    // extent of the top surface
+
+    boundBox surfBounds(searchSurf.points());
+    boundBox searchBounds
+    (
+        surfBounds.min() - searchSpan*zDir, surfBounds.max()
+    );
+
+    const pointField& C = mesh.cellCentres();
+
+    // List of cells within the porous region
+    labelList porousCells(C.size());
+    label porousCelli = 0;
+
+    forAll(C, celli)
+    {
+        if (searchBounds.contains(C[celli]))
+        {
+            porousCells[porousCelli++] = celli;
+        }
+    }
+
+    porousCells.setSize(porousCelli);
+
+    pointField start(porousCelli);
+    pointField end(porousCelli);
+
+    forAll(porousCells, porousCelli)
+    {
+        start[porousCelli] = C[porousCells[porousCelli]];
+        end[porousCelli] = start[porousCelli] + searchSpan*zDir;
+    }
+
+    // Field of distance between the cell centre and the corresponding point
+    // on the porous region top surface
+    scalarField zTop(porousCelli);
+
+    // Search the porous cells for those with a corresponding point on the
+    // porous region top surface
+    List<pointIndexHit> hitInfo;
+    searchSurf.findLine(start, end, hitInfo);
+
+    porousCelli = 0;
+
+    forAll(porousCells, celli)
+    {
+        const pointIndexHit& hit = hitInfo[celli];
+
+        if (hit.hit())
+        {
+            porousCells[porousCelli] = porousCells[celli];
+
+            zTop[porousCelli] =
+                (hit.hitPoint() - C[porousCells[porousCelli]]) & zDir;
+
+            porousCelli++;
+        }
+    }
+
+    // Resize the porous cells list and height field
+    porousCells.setSize(porousCelli);
+    zTop.setSize(porousCelli);
+
+    // Create a triSurface for the ground patch(s)
+    triSurface groundSurface
+    (
+        triSurfaceTools::triangulate
+        (
+            mesh.boundaryMesh(),
+            mesh.boundaryMesh().patchSet(groundPatches),
+            searchBounds
+        )
+    );
+
+    // Combine the ground triSurfaces across all processors
+    if (Pstream::parRun())
+    {
+        List<List<labelledTri>> groundSurfaceProcTris(Pstream::nProcs());
+        List<pointField> groundSurfaceProcPoints(Pstream::nProcs());
+
+        groundSurfaceProcTris[Pstream::myProcNo()] = groundSurface;
+        groundSurfaceProcPoints[Pstream::myProcNo()] = groundSurface.points();
+
+        Pstream::gatherList(groundSurfaceProcTris);
+        Pstream::scatterList(groundSurfaceProcTris);
+        Pstream::gatherList(groundSurfaceProcPoints);
+        Pstream::scatterList(groundSurfaceProcPoints);
+
+        label nTris = 0;
+        forAll(groundSurfaceProcTris, i)
+        {
+            nTris += groundSurfaceProcTris[i].size();
+        }
+
+        List<labelledTri> groundSurfaceTris(nTris);
+        label trii = 0;
+        label offset = 0;
+        forAll(groundSurfaceProcTris, i)
+        {
+            forAll(groundSurfaceProcTris[i], j)
+            {
+                groundSurfaceTris[trii] = groundSurfaceProcTris[i][j];
+                groundSurfaceTris[trii][0] += offset;
+                groundSurfaceTris[trii][1] += offset;
+                groundSurfaceTris[trii][2] += offset;
+                trii++;
+            }
+            offset += groundSurfaceProcPoints[i].size();
+        }
+
+        label nPoints = 0;
+        forAll(groundSurfaceProcPoints, i)
+        {
+            nPoints += groundSurfaceProcPoints[i].size();
+        }
+
+        pointField groundSurfacePoints(nPoints);
+
+        label pointi = 0;
+        forAll(groundSurfaceProcPoints, i)
+        {
+            forAll(groundSurfaceProcPoints[i], j)
+            {
+                groundSurfacePoints[pointi++] = groundSurfaceProcPoints[i][j];
+            }
+        }
+
+        groundSurface = triSurface(groundSurfaceTris, groundSurfacePoints);
+    }
+
+    // Create a searchable triSurface for the ground
+    triSurfaceSearch groundSearch(groundSurface);
+
+    // Search the porous cells for the corresponding point on the ground
+
+    start.setSize(porousCelli);
+    end.setSize(porousCelli);
+
+    forAll(porousCells, porousCelli)
+    {
+        start[porousCelli] = C[porousCells[porousCelli]];
+        end[porousCelli] = start[porousCelli] - searchSpan*zDir;
+    }
+
+    groundSearch.findLine(start, end, hitInfo);
+
+    scalarField zBottom(porousCelli);
+
+    forAll(porousCells, porousCelli)
+    {
+        const pointIndexHit& hit = hitInfo[porousCelli];
+
+        if (hit.hit())
+        {
+            zBottom[porousCelli] =
+                (C[porousCells[porousCelli]] - hit.hitPoint()) & zDir;
+        }
+    }
+
+    // Create the normalized height field
+    scalarField zNorm(zBottom/(zBottom + zTop));
+
+    // Create the porosity surface area per unit volume zone field
+    Sigma_ = SigmaFunc->value(zNorm);
+
+    // Create the porous region cellZone and add to the mesh cellZones
+
+    cellZoneMesh& cellZones = const_cast<cellZoneMesh&>(mesh.cellZones());
+
+    label zoneID = cellZones.findZoneID(zoneName_);
+
+    if (zoneID == -1)
+    {
+        zoneID = cellZones.size();
+        cellZones.setSize(zoneID + 1);
+
+        cellZones.set
+        (
+            zoneID,
+            new cellZone
+            (
+                zoneName_,
+                porousCells,
+                zoneID,
+                cellZones
+            )
+        );
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "Unable to create porous cellZone " << zoneName_
+            << ": zone already exists"
+            << abort(FatalError);
+    }
+}
+
+
+Foam::porosityModels::powerLawLopesdaCosta::powerLawLopesdaCosta
+(
+    const word& name,
+    const word& modelType,
+    const fvMesh& mesh,
+    const dictionary& dict,
+    const word& dummyCellZoneName
+)
+:
+    powerLawLopesdaCostaZone(name, modelType, mesh, dict),
+    porosityModel
+    (
+        name,
+        modelType,
+        mesh,
+        dict,
+        powerLawLopesdaCostaZone::zoneName_
+    ),
+    Cd_(readScalar(coeffs_.lookup("Cd"))),
+    C1_(readScalar(coeffs_.lookup("C1"))),
+    rhoName_(coeffs_.lookupOrDefault<word>("rho", "rho"))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::porosityModels::powerLawLopesdaCosta::~powerLawLopesdaCosta()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::scalarField&
+Foam::porosityModels::powerLawLopesdaCostaZone::Sigma() const
+{
+    return Sigma_;
+}
+
+
+void Foam::porosityModels::powerLawLopesdaCosta::calcTransformModelData()
+{}
+
+
+void Foam::porosityModels::powerLawLopesdaCosta::calcForce
+(
+    const volVectorField& U,
+    const volScalarField& rho,
+    const volScalarField& mu,
+    vectorField& force
+) const
+{
+    scalarField Udiag(U.size(), 0.0);
+    const scalarField& V = mesh_.V();
+
+    apply(Udiag, V, rho, U);
+
+    force = Udiag*U;
+}
+
+
+void Foam::porosityModels::powerLawLopesdaCosta::correct
+(
+    fvVectorMatrix& UEqn
+) const
+{
+    const vectorField& U = UEqn.psi();
+    const scalarField& V = mesh_.V();
+    scalarField& Udiag = UEqn.diag();
+
+    if (UEqn.dimensions() == dimForce)
+    {
+        const volScalarField& rho =
+            mesh_.lookupObject<volScalarField>(rhoName_);
+
+        apply(Udiag, V, rho, U);
+    }
+    else
+    {
+        apply(Udiag, V, geometricOneField(), U);
+    }
+}
+
+
+void Foam::porosityModels::powerLawLopesdaCosta::correct
+(
+    fvVectorMatrix& UEqn,
+    const volScalarField& rho,
+    const volScalarField& mu
+) const
+{
+    const vectorField& U = UEqn.psi();
+    const scalarField& V = mesh_.V();
+    scalarField& Udiag = UEqn.diag();
+
+    apply(Udiag, V, rho, U);
+}
+
+
+void Foam::porosityModels::powerLawLopesdaCosta::correct
+(
+    const fvVectorMatrix& UEqn,
+    volTensorField& AU
+) const
+{
+    const vectorField& U = UEqn.psi();
+
+    if (UEqn.dimensions() == dimForce)
+    {
+        const volScalarField& rho =
+            mesh_.lookupObject<volScalarField>(rhoName_);
+
+        apply(AU, rho, U);
+    }
+    else
+    {
+        apply(AU, geometricOneField(), U);
+    }
+}
+
+
+bool Foam::porosityModels::powerLawLopesdaCosta::writeData(Ostream& os) const
+{
+    os  << indent << name_ << endl;
+    dict_.write(os);
+
+    return true;
+}
+
+
+// ************************************************************************* //

src/atmosphericModels/porosityModels/powerLawLopesdaCosta/powerLawLopesdaCosta.H

@@ -0,0 +1,229 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 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::powerLawLopesdaCosta
+
+Description
+    Variant of the power law porosity model with spatially varying
+    drag coefficient
+
+    given by:
+
+        \f[
+            S = -\rho C_d \Sigma |U|^{(C_1 - 1)} U
+        \f]
+
+    where
+    \vartable
+        \Sigma | Porosity surface area per unit volume
+        C_d    | Model linear coefficient
+        C_1    | Model exponent coefficient
+    \endvartable
+
+    Reference:
+    \verbatim
+        Costa, J. C. P. L. D. (2007).
+        Atmospheric flow over forested and non-forested complex terrain.
+    \endverbatim
+
+See also
+    Foam::RASModels::kEpsilonLopesdaCosta
+
+SourceFiles
+    powerLawLopesdaCosta.C
+    powerLawLopesdaCostaTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef powerLawLopesdaCosta_H
+#define powerLawLopesdaCosta_H
+
+#include "porosityModel.H"
+#include "Function1.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace porosityModels
+{
+
+/*---------------------------------------------------------------------------*\
+                  Class powerLawLopesdaCostaZone Declaration
+\*---------------------------------------------------------------------------*/
+
+class powerLawLopesdaCostaZone
+{
+protected:
+
+    // Protected data
+
+        //- Automatically generated zone name for this porous zone
+        const word zoneName_;
+
+        //- Porosity surface area per unit volume zone field
+        scalarField Sigma_;
+
+
+public:
+
+    //- Constructor
+    powerLawLopesdaCostaZone
+    (
+        const word& name,
+        const word& modelType,
+        const fvMesh& mesh,
+        const dictionary& dict
+    );
+
+    // Member Functions
+
+        //- Return the porosity surface area per unit volume zone field
+        const scalarField& Sigma() const;
+};
+
+
+/*---------------------------------------------------------------------------*\
+                      Class powerLawLopesdaCosta Declaration
+\*---------------------------------------------------------------------------*/
+
+class powerLawLopesdaCosta
+:
+    public powerLawLopesdaCostaZone,
+    public porosityModel
+{
+    // Private data
+
+        //- Cd coefficient
+        scalar Cd_;
+
+        //- C1 coefficient
+        scalar C1_;
+
+        //- Name of density field
+        word rhoName_;
+
+
+    // Private Member Functions
+
+        //- Apply resistance
+        template<class RhoFieldType>
+        void apply
+        (
+            scalarField& Udiag,
+            const scalarField& V,
+            const RhoFieldType& rho,
+            const vectorField& U
+        ) const;
+
+        //- Apply resistance
+        template<class RhoFieldType>
+        void apply
+        (
+            tensorField& AU,
+            const RhoFieldType& rho,
+            const vectorField& U
+        ) const;
+
+        //- Disallow default bitwise copy construct
+        powerLawLopesdaCosta(const powerLawLopesdaCosta&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const powerLawLopesdaCosta&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("powerLawLopesdaCosta");
+
+    //- Constructor
+    powerLawLopesdaCosta
+    (
+        const word& name,
+        const word& modelType,
+        const fvMesh& mesh,
+        const dictionary& dict,
+        const word& cellZoneName
+    );
+
+    //- Destructor
+    virtual ~powerLawLopesdaCosta();
+
+
+    // Member Functions
+
+        //- Transform the model data wrt mesh changes
+        virtual void calcTransformModelData();
+
+        //- Calculate the porosity force
+        virtual void calcForce
+        (
+            const volVectorField& U,
+            const volScalarField& rho,
+            const volScalarField& mu,
+            vectorField& force
+        ) const;
+
+        //- Add resistance
+        virtual void correct(fvVectorMatrix& UEqn) const;
+
+        //- Add resistance
+        virtual void correct
+        (
+            fvVectorMatrix& UEqn,
+            const volScalarField& rho,
+            const volScalarField& mu
+        ) const;
+
+        //- Add resistance
+        virtual void correct
+        (
+            const fvVectorMatrix& UEqn,
+            volTensorField& AU
+        ) const;
+
+
+    // I-O
+
+        //- Write
+        bool writeData(Ostream& os) const;
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace porosityModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "powerLawLopesdaCostaTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/atmosphericModels/porosityModels/powerLawLopesdaCosta/powerLawLopesdaCostaTemplates.C

@@ -0,0 +1,87 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 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 "volFields.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class RhoFieldType>
+void Foam::porosityModels::powerLawLopesdaCosta::apply
+(
+    scalarField& Udiag,
+    const scalarField& V,
+    const RhoFieldType& rho,
+    const vectorField& U
+) const
+{
+    const scalar C1m1b2 = (C1_ - 1.0)/2.0;
+
+    forAll(cellZoneIDs_, zonei)
+    {
+        const labelList& cells = mesh_.cellZones()[cellZoneIDs_[zonei]];
+
+        forAll(cells, i)
+        {
+            const label celli = cells[i];
+
+            Udiag[celli] +=
+                V[celli]*rho[celli]
+               *Cd_*Sigma_[i]*pow(magSqr(U[celli]), C1m1b2);
+        }
+    }
+}
+
+
+template<class RhoFieldType>
+void Foam::porosityModels::powerLawLopesdaCosta::apply
+(
+    tensorField& AU,
+    const RhoFieldType& rho,
+    const vectorField& U
+) const
+{
+    const scalar C1m1b2 = (C1_ - 1.0)/2.0;
+
+    forAll(cellZoneIDs_, zonei)
+    {
+        const labelList& cells = mesh_.cellZones()[cellZoneIDs_[zonei]];
+
+        forAll(cells, i)
+        {
+            const label celli = cells[i];
+
+            AU[celli] =
+                AU[celli]
+              + I
+               *(
+                   0.5*rho[celli]*Cd_*Sigma_[i]
+                  *pow(magSqr(U[celli]), C1m1b2)
+                );
+        }
+    }
+}
+
+
+// ************************************************************************* //

src/dynamicMesh/Make/files

@@ -112,6 +112,7 @@ motionSolvers/velocity/velocityDisplacement/velocityDisplacementMotionSolver.C
 motionSolvers/componentVelocity/componentVelocityMotionSolver.C
 motionSolvers/displacement/solidBody/solidBodyMotionSolver.C
 motionSolvers/displacement/solidBody/multiSolidBodyMotionSolver.C
+motionSolvers/displacement/codedPoints0/codedPoints0MotionSolver.C
 
 solidBodyMotionFunctions = motionSolvers/displacement/solidBody/solidBodyMotionFunctions
 $(solidBodyMotionFunctions)/solidBodyMotionFunction/solidBodyMotionFunction.C

src/dynamicMesh/motionSolvers/displacement/codedPoints0/codedPoints0MotionSolver.C

@@ -0,0 +1,199 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "codedPoints0MotionSolver.H"
+#include "dictionary.H"
+#include "Time.H"
+#include "SHA1Digest.H"
+#include "dynamicCode.H"
+#include "dynamicCodeContext.H"
+#include "stringOps.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(codedPoints0MotionSolver, 0);
+
+    addToRunTimeSelectionTable
+    (
+        motionSolver,
+        codedPoints0MotionSolver,
+        dictionary
+    );
+}
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+void Foam::codedPoints0MotionSolver::prepare
+(
+    dynamicCode& dynCode,
+    const dynamicCodeContext& context
+) const
+{
+    // Set additional rewrite rules
+    dynCode.setFilterVariable("typeName", name_);
+
+    // Compile filtered C template
+    dynCode.addCompileFile("codedPoints0MotionSolverTemplate.C");
+
+    // Copy filtered H template
+    dynCode.addCopyFile("codedPoints0MotionSolverTemplate.H");
+
+    // Debugging: make BC verbose
+    // dynCode.setFilterVariable("verbose", "true");
+    // Info<<"compile " << name_ << " sha1: "
+    //     << context.sha1() << endl;
+
+    // Define Make/options
+    dynCode.setMakeOptions
+    (
+        "EXE_INC = -g \\\n"
+        "-I$(LIB_SRC)/finiteVolume/lnInclude \\\n"
+        "-I$(LIB_SRC)/fvMotionSolvers/lnInclude \\\n"
+        "-I$(LIB_SRC)/dynamicMesh/lnInclude \\\n"
+        "-I$(LIB_SRC)/meshTools/lnInclude \\\n"
+      + context.options()
+      + "\n\nLIB_LIBS = \\\n"
+      + "    -lfvMotionSolvers \\\n"
+      + context.libs()
+    );
+}
+
+
+Foam::dlLibraryTable& Foam::codedPoints0MotionSolver::libs() const
+{
+    return const_cast<Time&>(mesh().time()).libs();
+}
+
+
+Foam::string Foam::codedPoints0MotionSolver::description() const
+{
+    return "points0MotionSolver " + name();
+}
+
+
+void Foam::codedPoints0MotionSolver::clearRedirect() const
+{
+    redirectMotionSolverPtr_.clear();
+}
+
+
+const Foam::dictionary&
+Foam::codedPoints0MotionSolver::codeDict() const
+{
+    return coeffDict();
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::codedPoints0MotionSolver::codedPoints0MotionSolver
+(
+    const polyMesh& mesh,
+    const IOdictionary& dict
+)
+:
+    motionSolver(mesh, dict, typeName),
+    codedBase()
+{
+    // Backward compatibility
+    if (dict.found("redirectType"))
+    {
+        dict.lookup("redirectType") >> name_;
+    }
+    else
+    {
+        dict.lookup("name") >> name_;
+    }
+
+    updateLibrary(name_);
+    redirectMotionSolver();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::codedPoints0MotionSolver::~codedPoints0MotionSolver()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::motionSolver&
+Foam::codedPoints0MotionSolver::redirectMotionSolver() const
+{
+    if (!redirectMotionSolverPtr_.valid())
+    {
+        // Get the dictionary for the solver and override the
+        // solver name (in case it is not a subdictionary and contains
+        // the 'coded' as the motionSolver)
+        dictionary constructDict(coeffDict());
+        constructDict.set("solver", name_);
+        constructDict.set("motionSolver", name_);
+
+        IOobject io(*this);
+        io.readOpt() = IOobject::NO_READ;
+
+        redirectMotionSolverPtr_ = motionSolver::New
+        (
+            mesh(),
+            IOdictionary(io, constructDict)
+        );
+    }
+    return redirectMotionSolverPtr_();
+}
+
+
+Foam::tmp<Foam::pointField> Foam::codedPoints0MotionSolver::curPoints() const
+{
+    updateLibrary(name_);
+    return redirectMotionSolver().curPoints();
+}
+
+
+void Foam::codedPoints0MotionSolver::solve()
+{
+    updateLibrary(name_);
+    redirectMotionSolver().solve();
+}
+
+
+void Foam::codedPoints0MotionSolver::movePoints(const pointField& fld)
+{
+    updateLibrary(name_);
+    return redirectMotionSolver().movePoints(fld);
+}
+
+
+void Foam::codedPoints0MotionSolver::updateMesh(const mapPolyMesh& mpm)
+{
+    updateLibrary(name_);
+    return redirectMotionSolver().updateMesh(mpm);
+}
+
+
+// ************************************************************************* //

src/dynamicMesh/motionSolvers/displacement/codedPoints0/codedPoints0MotionSolver.H

@@ -0,0 +1,176 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::codedPoints0MotionSolver
+
+Group
+
+Description
+    Provides a general interface to enable dynamic code compilation of
+    mesh motion solvers.
+
+    The entries are:
+    \plaintable
+       codeInclude | include files
+       codeOptions | include paths; inserted into EXE_INC in Make/options
+       codeLibs    | link line; inserted into LIB_LIBS in Make/options
+       localCode   | c++; local static functions;
+       code        | c++; upon motionSolver::curPoints();
+    \endplaintable
+
+    Note that the dynamically generated motionSolver is an points0MotionSolver,
+    i.e. it has access to the points0 (points from constant/polyMesh).
+
+Usage
+    Example of simplistic specification:
+    \verbatim
+    difference
+    {
+        type coded;
+        // Name of on-the-fly generated motion solver
+        name myMotion;
+        code
+        #{
+            const pointField& p0 = points0();
+            const tensor rotT =
+                rotationTensor(vector(1, 0, 0), vector(0, 1, 0));
+            return transform(rotT, p0);
+        #};
+    }
+    \endverbatim
+
+See also
+    Foam::motionSolver
+    Foam::codedBase
+
+SourceFiles
+    codedPoints0MotionSolver.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef codedPoints0MotionSolver_H
+#define codedPoints0MotionSolver_H
+
+#include "motionSolver.H"
+#include "codedBase.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+/*---------------------------------------------------------------------------*\
+                  Class codedPoints0MotionSolver Declaration
+\*---------------------------------------------------------------------------*/
+
+class codedPoints0MotionSolver
+:
+    public motionSolver,
+    public codedBase
+{
+protected:
+
+    // Protected data
+
+        //- Name of redirected motion solver
+        word name_;
+
+        //- Underlying motionSolver
+        mutable autoPtr<motionSolver> redirectMotionSolverPtr_;
+
+
+    // Protected Member Functions
+
+        //- Get the loaded dynamic libraries
+        virtual dlLibraryTable& libs() const;
+
+        //- Adapt the context for the current object
+        virtual void prepare(dynamicCode&, const dynamicCodeContext&) const;
+
+        // Return a description (type + name) for the output
+        virtual string description() const;
+
+        // Clear any redirected objects
+        virtual void clearRedirect() const;
+
+        // Get the dictionary to initialize the codeContext
+        virtual const dictionary& codeDict() const;
+
+
+private:
+
+        //- Disallow default bitwise copy construct
+        codedPoints0MotionSolver(const codedPoints0MotionSolver&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const codedPoints0MotionSolver&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("coded");
+
+
+    // Constructors
+
+        //- Construct from mesh and dictionary
+        codedPoints0MotionSolver
+        (
+            const polyMesh& mesh,
+            const IOdictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~codedPoints0MotionSolver();
+
+
+    // Member Functions
+
+        //- Dynamically compiled motionSolver
+        motionSolver& redirectMotionSolver() const;
+
+        //- Return point location obtained from the current motion field
+        virtual tmp<pointField> curPoints() const;
+
+        //- Solve for motion
+        virtual void solve();
+
+        //- Update local data for geometry changes
+        virtual void movePoints(const pointField&);
+
+        //-  Update local data for topology changes
+        virtual void updateMesh(const mapPolyMesh&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.C

@@ -1564,6 +1564,7 @@ void Foam::hexRef8::walkFaceFromMid
 Foam::label Foam::hexRef8::faceConsistentRefinement
 (
     const bool maxSet,
+    const labelUList& cellLevel,
     bitSet& refineCell
 ) const
 {
@@ -1573,10 +1574,10 @@ Foam::label Foam::hexRef8::faceConsistentRefinement
     for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
     {
         label own = mesh_.faceOwner()[facei];
-        label nei = mesh_.faceNeighbour()[facei];
+        label ownLevel = cellLevel[own] + refineCell.get(own);
 
-        label ownLevel = cellLevel_[own] + refineCell.get(own);
+        label nei = mesh_.faceNeighbour()[facei];
-        label neiLevel = cellLevel_[nei] + refineCell.get(nei);
+        label neiLevel = cellLevel[nei] + refineCell.get(nei);
 
         if (ownLevel > (neiLevel+1))
         {
@@ -1613,7 +1614,7 @@ Foam::label Foam::hexRef8::faceConsistentRefinement
     {
         label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
 
-        neiLevel[i] = cellLevel_[own] + refineCell.get(own);
+        neiLevel[i] = cellLevel[own] + refineCell.get(own);
     }
 
     // Swap to neighbour
@@ -1623,7 +1624,7 @@ Foam::label Foam::hexRef8::faceConsistentRefinement
     forAll(neiLevel, i)
     {
         label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
-        label ownLevel = cellLevel_[own] + refineCell.get(own);
+        label ownLevel = cellLevel[own] + refineCell.get(own);
 
         if (ownLevel > (neiLevel[i]+1))
         {
@@ -1650,6 +1651,7 @@ Foam::label Foam::hexRef8::faceConsistentRefinement
 // Debug: check if wanted refinement is compatible with 2:1
 void Foam::hexRef8::checkWantedRefinementLevels
 (
+    const labelUList& cellLevel,
     const labelList& cellsToRefine
 ) const
 {
@@ -1658,10 +1660,10 @@ void Foam::hexRef8::checkWantedRefinementLevels
     for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
     {
         label own = mesh_.faceOwner()[facei];
-        label nei = mesh_.faceNeighbour()[facei];
+        label ownLevel = cellLevel[own] + refineCell.get(own);
 
-        label ownLevel = cellLevel_[own] + refineCell.get(own);
+        label nei = mesh_.faceNeighbour()[facei];
-        label neiLevel = cellLevel_[nei] + refineCell.get(nei);
+        label neiLevel = cellLevel[nei] + refineCell.get(nei);
 
         if (mag(ownLevel-neiLevel) > 1)
         {
@@ -1669,11 +1671,11 @@ void Foam::hexRef8::checkWantedRefinementLevels
             dumpCell(nei);
             FatalErrorInFunction
                 << "cell:" << own
-                << " current level:" << cellLevel_[own]
+                << " current level:" << cellLevel[own]
                 << " level after refinement:" << ownLevel
                 << nl
                 << "neighbour cell:" << nei
-                << " current level:" << cellLevel_[nei]
+                << " current level:" << cellLevel[nei]
                 << " level after refinement:" << neiLevel
                 << nl
                 << "which does not satisfy 2:1 constraints anymore."
@@ -1689,7 +1691,7 @@ void Foam::hexRef8::checkWantedRefinementLevels
     {
         label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
 
-        neiLevel[i] = cellLevel_[own] + refineCell.get(own);
+        neiLevel[i] = cellLevel[own] + refineCell.get(own);
     }
 
     // Swap to neighbour
@@ -1701,7 +1703,7 @@ void Foam::hexRef8::checkWantedRefinementLevels
         label facei = i + mesh_.nInternalFaces();
 
         label own = mesh_.faceOwner()[facei];
-        label ownLevel = cellLevel_[own] + refineCell.get(own);
+        label ownLevel = cellLevel[own] + refineCell.get(own);
 
         if (mag(ownLevel - neiLevel[i]) > 1)
         {
@@ -1715,7 +1717,7 @@ void Foam::hexRef8::checkWantedRefinementLevels
                 << " on patch " << patchi << " "
                 << mesh_.boundaryMesh()[patchi].name()
                 << " owner cell " << own
-                << " current level:" << cellLevel_[own]
+                << " current level:" << cellLevel[own]
                 << " level after refinement:" << ownLevel
                 << nl
                 << " (coupled) neighbour cell will get refinement "
@@ -2251,6 +2253,7 @@ Foam::hexRef8::hexRef8
 
 Foam::labelList Foam::hexRef8::consistentRefinement
 (
+    const labelUList& cellLevel,
     const labelList& cellsToRefine,
     const bool maxSet
 ) const
@@ -2264,7 +2267,12 @@ Foam::labelList Foam::hexRef8::consistentRefinement
 
     while (true)
     {
-        label nChanged = faceConsistentRefinement(maxSet, refineCell);
+        label nChanged = faceConsistentRefinement
+        (
+            maxSet,
+            cellLevel,
+            refineCell
+        );
 
         reduce(nChanged, sumOp<label>());
 
@@ -2286,7 +2294,7 @@ Foam::labelList Foam::hexRef8::consistentRefinement
 
     if (debug)
     {
-        checkWantedRefinementLevels(newCellsToRefine);
+        checkWantedRefinementLevels(cellLevel, newCellsToRefine);
     }
 
     return newCellsToRefine;
@@ -3089,11 +3097,11 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement2
             refineCell.set(celli);
         }
     }
-    faceConsistentRefinement(true, refineCell);
+    faceConsistentRefinement(true, cellLevel_, refineCell);
 
     while (true)
     {
-        label nChanged = faceConsistentRefinement(true, refineCell);
+        label nChanged = faceConsistentRefinement(true, cellLevel_, refineCell);
 
         reduce(nChanged, sumOp<label>());
 
@@ -3141,7 +3149,7 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement2
 
         const bitSet savedRefineCell(refineCell);
 
-        label nChanged = faceConsistentRefinement(true, refineCell);
+        label nChanged = faceConsistentRefinement(true, cellLevel_, refineCell);
 
         {
             cellSet cellsOut2

src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.H

@@ -300,11 +300,16 @@ class hexRef8
         label faceConsistentRefinement
         (
             const bool maxSet,
+            const labelUList& cellLevel,
             bitSet& refineCell
         ) const;
 
         //- Check wanted refinement for 2:1 consistency
-        void checkWantedRefinementLevels(const labelList&) const;
+        void checkWantedRefinementLevels
+        (
+            const labelUList& cellLevel,
+            const labelList&
+        ) const;
 
 
         // Cellshape recognition
@@ -420,10 +425,25 @@ public:
             //  removes cells to refine (maxSet = false)
             labelList consistentRefinement
             (
+                const labelUList& cellLevel,
                 const labelList& cellsToRefine,
                 const bool maxSet
             ) const;
 
+            //- Given valid mesh and current cell level and proposed
+            //  cells to refine calculate any clashes (due to 2:1) and return
+            //  ok list of cells to refine.
+            //  Either adds cells to refine to set (maxSet = true) or
+            //  removes cells to refine (maxSet = false)
+            labelList consistentRefinement
+            (
+                const labelList& cellsToRefine,
+                const bool maxSet
+            ) const
+            {
+                return consistentRefinement(cellLevel_, cellsToRefine, maxSet);
+            }
+
             //- Like consistentRefinement but slower:
             //
             //  - specify number of cells between consecutive refinement levels

src/finiteArea/faMatrices/faMatrix/faMatrix.H

@@ -71,7 +71,7 @@ Ostream& operator<<(Ostream&, const faMatrix<Type>&);
 template<class Type>
 class faMatrix
 :
-    public tmp<faMatrix<Type>>::refCount,
+    public refCount,
     public lduMatrix
 {
     // Private data

src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -223,6 +223,9 @@ public:
         //- Return const access to the cell zone IDs
         inline const labelList& cellZoneIDs() const;
 
+        //- Return dictionary used for model construction
+        const dictionary& dict() const;
+
         //- Transform the model data wrt mesh changes
         virtual void transformModelData();
 

src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModelI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,6 +35,12 @@ inline bool Foam::porosityModel::active() const
 }
 
 
+inline const Foam::dictionary& Foam::porosityModel::dict() const
+{
+    return dict_;
+}
+
+
 inline const Foam::labelList& Foam::porosityModel::cellZoneIDs() const
 {
     return cellZoneIDs_;

src/finiteVolume/fields/fvPatchFields/derived/outletMappedUniformInlet/outletMappedUniformInletFvPatchField.C

@@ -43,6 +43,21 @@ outletMappedUniformInletFvPatchField
 {}
 
 
+template<class Type>
+Foam::outletMappedUniformInletFvPatchField<Type>::
+outletMappedUniformInletFvPatchField
+(
+    const fvPatch& p,
+    const DimensionedField<Type, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFvPatchField<Type>(p, iF, dict),
+    outletPatchName_(dict.lookup("outletPatch")),
+    phiName_(dict.lookupOrDefault<word>("phi", "phi"))
+{}
+
+
 template<class Type>
 Foam::outletMappedUniformInletFvPatchField<Type>::
 outletMappedUniformInletFvPatchField
@@ -59,21 +74,6 @@ outletMappedUniformInletFvPatchField
 {}
 
 
-template<class Type>
-Foam::outletMappedUniformInletFvPatchField<Type>::
-outletMappedUniformInletFvPatchField
-(
-    const fvPatch& p,
-    const DimensionedField<Type, volMesh>& iF,
-    const dictionary& dict
-)
-:
-    fixedValueFvPatchField<Type>(p, iF, dict),
-    outletPatchName_(dict.lookup("outletPatch")),
-    phiName_(dict.lookupOrDefault<word>("phi", "phi"))
-{}
-
-
 template<class Type>
 Foam::outletMappedUniformInletFvPatchField<Type>::
 outletMappedUniformInletFvPatchField

src/finiteVolume/fields/fvPatchFields/derived/outletMappedUniformInlet/outletMappedUniformInletFvPatchField.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -29,14 +29,14 @@ Group
 
 Description
     This boundary conditon averages the field over the "outlet" patch specified
-    by name "outletPatch" and applies this as the uniform value of the
+    by name "outletPatch" and applies this as the uniform value of the field
-    field over this patch.
+    over this patch.
 
 Usage
     \table
         Property        | Description             | Required    | Default value
-        outletPatch     | name of outlet patch    | yes         |
+        outletPatch     | Name of outlet patch    | yes         |
-        phi             | flux field name         | no          | phi
+        phi             | Flux field name         | no          | phi
     \endtable
 
     Example of the boundary condition specification:

src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.H

@@ -51,7 +51,7 @@ Usage
         setAverage   | Switch to activate setting of average value | no | false
         perturb      | Perturb points for regular geometries    | no | 1e-5
         points       | Name of points file        | no          | points
-        fieldTableName | Alternative field name to sample | no | this field name
+        fieldTable   | Alternative field name to sample | no | this field name
         mapMethod    | Type of mapping            | no | planarInterpolation
         offset       | Offset to mapped values    | no | Zero
     \endtable

src/finiteVolume/finiteVolume/convectionSchemes/convectionScheme/convectionScheme.C

@@ -46,7 +46,7 @@ namespace fv
 template<class Type>
 convectionScheme<Type>::convectionScheme(const convectionScheme& cs)
 :
-    tmp<convectionScheme<Type>>::refCount(),
+    refCount(),
     mesh_(cs.mesh_)
 {}
 

src/finiteVolume/finiteVolume/convectionSchemes/convectionScheme/convectionScheme.H

@@ -67,7 +67,7 @@ namespace fv
 template<class Type>
 class convectionScheme
 :
-    public tmp<convectionScheme<Type>>::refCount
+    public refCount
 {
     // Private data
 

src/finiteVolume/finiteVolume/d2dt2Schemes/d2dt2Scheme/d2dt2Scheme.C

@@ -21,9 +21,6 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
-Description
-    Abstract base class for finite volume calculus d2dt2 schemes.
-
 \*---------------------------------------------------------------------------*/
 
 #include "fv.H"
@@ -84,13 +81,6 @@ tmp<d2dt2Scheme<Type>> d2dt2Scheme<Type>::New
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class Type>
-d2dt2Scheme<Type>::~d2dt2Scheme()
-{}
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace fv

src/finiteVolume/finiteVolume/d2dt2Schemes/d2dt2Scheme/d2dt2Scheme.H

@@ -25,7 +25,7 @@ Class
     Foam::fv::d2dt2Scheme
 
 Description
-    Abstract base class for d2dt2 schemes.
+    Abstract base class for finite volume d2dt2 schemes.
 
 SourceFiles
     d2dt2Scheme.C
@@ -64,7 +64,7 @@ namespace fv
 template<class Type>
 class d2dt2Scheme
 :
-    public tmp<d2dt2Scheme<Type>>::refCount
+    public refCount
 {
 
 protected:
@@ -76,11 +76,11 @@ protected:
 
     // Private Member Functions
 
-        //- Disallow copy construct
+        //- No copy construct
-        d2dt2Scheme(const d2dt2Scheme&);
+        d2dt2Scheme(const d2dt2Scheme&) = delete;
 
-        //- Disallow default bitwise assignment
+        //- No copy assignment
-        void operator=(const d2dt2Scheme&);
+        void operator=(const d2dt2Scheme&) = delete;
 
 
 public:
@@ -127,7 +127,7 @@ public:
 
 
     //- Destructor
-    virtual ~d2dt2Scheme();
+    virtual ~d2dt2Scheme() = default;
 
 
     // Member Functions

src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.C

@@ -84,13 +84,6 @@ tmp<ddtScheme<Type>> ddtScheme<Type>::New
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class Type>
-ddtScheme<Type>::~ddtScheme()
-{}
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 template<class Type>
@@ -120,14 +113,11 @@ tmp<fvMatrix<Type>> ddtScheme<Type>::fvmDdt
 {
     NotImplemented;
 
-    return tmp<fvMatrix<Type>>
+    return tmp<fvMatrix<Type>>::New
-    (
-        new fvMatrix<Type>
     (
         vf,
         alpha.dimensions()*rho.dimensions()
         *vf.dimensions()*dimVol/dimTime
-        )
     );
 }
 

src/finiteVolume/finiteVolume/ddtSchemes/ddtScheme/ddtScheme.H

@@ -30,7 +30,6 @@ Group
 Description
     Abstract base class for ddt schemes.
 
-
 SourceFiles
     ddtScheme.C
 
@@ -68,7 +67,7 @@ namespace fv
 template<class Type>
 class ddtScheme
 :
-    public tmp<ddtScheme<Type>>::refCount
+    public refCount
 {
 
 protected:
@@ -83,11 +82,11 @@ protected:
 
     // Private Member Functions
 
-        //- Disallow copy construct
+        //- No copy construct
-        ddtScheme(const ddtScheme&);
+        ddtScheme(const ddtScheme&) = delete;
 
-        //- Disallow default bitwise assignment
+        //- No copy assignment
-        void operator=(const ddtScheme&);
+        void operator=(const ddtScheme&) = delete;
 
 
 public:
@@ -136,7 +135,7 @@ public:
 
 
     //- Destructor
-    virtual ~ddtScheme();
+    virtual ~ddtScheme() = default;
 
 
     // Member Functions

src/finiteVolume/finiteVolume/divSchemes/divScheme/divScheme.C

@@ -85,16 +85,6 @@ tmp<divScheme<Type>> divScheme<Type>::New
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class Type>
-divScheme<Type>::~divScheme()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace fv

src/finiteVolume/finiteVolume/divSchemes/divScheme/divScheme.H

@@ -67,7 +67,7 @@ namespace fv
 template<class Type>
 class divScheme
 :
-    public tmp<divScheme<Type>>::refCount
+    public refCount
 {
 
 protected:
@@ -80,11 +80,11 @@ protected:
 
     // Private Member Functions
 
-        //- Disallow copy construct
+        //- No copy construct
-        divScheme(const divScheme&);
+        divScheme(const divScheme&) = delete;
 
-        //- Disallow default bitwise assignment
+        //- No copy assignment
-        void operator=(const divScheme&);
+        void operator=(const divScheme&) = delete;
 
 
 public:
@@ -133,7 +133,7 @@ public:
 
 
     //- Destructor
-    virtual ~divScheme();
+    virtual ~divScheme() = default;
 
 
     // Member Functions

src/finiteVolume/finiteVolume/gradSchemes/gradScheme/gradScheme.C

@@ -72,12 +72,6 @@ Foam::tmp<Foam::fv::gradScheme<Type>> Foam::fv::gradScheme<Type>::New
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class Type>
-Foam::fv::gradScheme<Type>::~gradScheme()
-{}
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 template<class Type>

src/finiteVolume/finiteVolume/gradSchemes/gradScheme/gradScheme.H

@@ -60,7 +60,7 @@ namespace fv
 template<class Type>
 class gradScheme
 :
-    public tmp<gradScheme<Type>>::refCount
+    public refCount
 {
     // Private data
 
@@ -69,11 +69,11 @@ class gradScheme
 
     // Private Member Functions
 
-        //- Disallow copy construct
+        //- No copy construct
-        gradScheme(const gradScheme&);
+        gradScheme(const gradScheme&) = delete;
 
-        //- Disallow default bitwise assignment
+        //- No copy assignment
-        void operator=(const gradScheme&);
+        void operator=(const gradScheme&) = delete;
 
 
 public:
@@ -114,7 +114,7 @@ public:
 
 
     //- Destructor
-    virtual ~gradScheme();
+    virtual ~gradScheme() = default;
 
 
     // Member Functions

src/finiteVolume/finiteVolume/laplacianSchemes/laplacianScheme/laplacianScheme.C

@@ -83,13 +83,6 @@ tmp<laplacianScheme<Type, GType>> laplacianScheme<Type, GType>::New
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class Type, class GType>
-laplacianScheme<Type, GType>::~laplacianScheme()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type, class GType>

src/finiteVolume/finiteVolume/laplacianSchemes/laplacianScheme/laplacianScheme.H

@@ -68,7 +68,7 @@ namespace fv
 template<class Type, class GType>
 class laplacianScheme
 :
-    public tmp<laplacianScheme<Type, GType>>::refCount
+    public refCount
 {
 
 protected:
@@ -84,11 +84,11 @@ private:
 
     // Private Member Functions
 
-        //- Disallow copy construct
+        //- No copy construct
-        laplacianScheme(const laplacianScheme&);
+        laplacianScheme(const laplacianScheme&) = delete;
 
-        //- Disallow default bitwise assignment
+        //- No copy assignment
-        void operator=(const laplacianScheme&);
+        void operator=(const laplacianScheme&) = delete;
 
 
 public:
@@ -162,7 +162,7 @@ public:
 
 
     //- Destructor
-    virtual ~laplacianScheme();
+    virtual ~laplacianScheme() = default;
 
 
     // Member Functions

src/finiteVolume/finiteVolume/snGradSchemes/snGradScheme/snGradScheme.C

@@ -85,13 +85,6 @@ tmp<snGradScheme<Type>> snGradScheme<Type>::New
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class Type>
-snGradScheme<Type>::~snGradScheme()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type>

src/finiteVolume/finiteVolume/snGradSchemes/snGradScheme/snGradScheme.H

@@ -63,7 +63,7 @@ namespace fv
 template<class Type>
 class snGradScheme
 :
-    public tmp<snGradScheme<Type>>::refCount
+    public refCount
 {
     // Private data
 
@@ -73,11 +73,11 @@ class snGradScheme
 
     // Private Member Functions
 
-        //- Disallow copy construct
+        //- No copy construct
-        snGradScheme(const snGradScheme&);
+        snGradScheme(const snGradScheme&) = delete;
 
-        //- Disallow default bitwise assignment
+        //- No copy assignment
-        void operator=(const snGradScheme&);
+        void operator=(const snGradScheme&) = delete;
 
 
 public:
@@ -118,7 +118,7 @@ public:
 
 
     //- Destructor
-    virtual ~snGradScheme();
+    virtual ~snGradScheme() = default;
 
 
     // Member Functions

src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C

@@ -321,7 +321,7 @@ Foam::fvMatrix<Type>::fvMatrix
 template<class Type>
 Foam::fvMatrix<Type>::fvMatrix(const fvMatrix<Type>& fvm)
 :
-    tmp<fvMatrix<Type>>::refCount(),
+    refCount(),
     lduMatrix(fvm),
     psi_(fvm.psi_),
     dimensions_(fvm.dimensions_),

src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H

@@ -114,7 +114,7 @@ template<class T> class UIndirectList;
 template<class Type>
 class fvMatrix
 :
-    public tmp<fvMatrix<Type>>::refCount,
+    public refCount,
     public lduMatrix
 {
     // Private data

src/finiteVolume/interpolation/surfaceInterpolation/multivariateSchemes/multivariateSurfaceInterpolationScheme/multivariateSurfaceInterpolationScheme.C

@@ -21,9 +21,6 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
-Description
-    Abstract base class for surface interpolation schemes.
-
 \*---------------------------------------------------------------------------*/
 
 #include "fv.H"
@@ -85,12 +82,4 @@ Foam::multivariateSurfaceInterpolationScheme<Type>::New
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class Type>
-Foam::multivariateSurfaceInterpolationScheme<Type>::
-~multivariateSurfaceInterpolationScheme()
-{}
-
-
 // ************************************************************************* //

src/finiteVolume/interpolation/surfaceInterpolation/multivariateSchemes/multivariateSurfaceInterpolationScheme/multivariateSurfaceInterpolationScheme.H

@@ -50,7 +50,7 @@ namespace Foam
 template<class Type>
 class multivariateSurfaceInterpolationScheme
 :
-    public tmp<multivariateSurfaceInterpolationScheme<Type>>::refCount
+    public refCount
 {
 
 public:
@@ -85,14 +85,14 @@ private:
 
     // Private Member Functions
 
-        //- Disallow default bitwise copy construct
+        //- No copy construct
         multivariateSurfaceInterpolationScheme
         (
             const multivariateSurfaceInterpolationScheme&
-        );
+        ) = delete;
 
-        //- Disallow default bitwise assignment
+        //- No copy assignment
-        void operator=(const multivariateSurfaceInterpolationScheme&);
+        void operator=(const multivariateSurfaceInterpolationScheme&) = delete;
 
 
 public:
@@ -143,7 +143,7 @@ public:
 
 
     //- Destructor
-    virtual ~multivariateSurfaceInterpolationScheme();
+    virtual ~multivariateSurfaceInterpolationScheme() = default;
 
 
     // Member Functions