From 72ec717d28ec82725bf6152d6853a87abd2afb91 Mon Sep 17 00:00:00 2001
From: Vuko Vukcevic <v.vukcevic@wikki.co.uk>
Date: Tue, 20 Mar 2018 14:28:51 +0100
Subject: [PATCH 01/15] Bugfix in pointPatchInterpolation.C

Certain demand driven data not properly cleaned-up on mesh motion (without topo
changes).
---
 .../pointPatchInterpolation/pointPatchInterpolation.C  | 10 ++--------
 1 file changed, 2 insertions(+), 8 deletions(-)

diff --git a/src/finiteVolume/interpolation/volPointInterpolation/pointPatchInterpolation/pointPatchInterpolation.C b/src/finiteVolume/interpolation/volPointInterpolation/pointPatchInterpolation/pointPatchInterpolation.C
index 8a346cb72..9a6ca03d6 100644
--- a/src/finiteVolume/interpolation/volPointInterpolation/pointPatchInterpolation/pointPatchInterpolation.C
+++ b/src/finiteVolume/interpolation/volPointInterpolation/pointPatchInterpolation/pointPatchInterpolation.C
@@ -407,14 +407,8 @@ const scalarListList& pointPatchInterpolation::patchPatchPointWeights() const
 
 bool pointPatchInterpolation::movePoints()
 {
-    if (!patchInterpolators_.empty())
-    {
-        // Reset patch interpolators (also lazy evaluated)
-        forAll (patchInterpolators_, patchi)
-        {
-            patchInterpolators_[patchi].movePoints();
-        }
-    }
+    // Clear out patch interpolators to force recalculation
+    patchInterpolators_.clear();
 
     // Clear out demand driven data
     clearOut();

From 049bc0ed951b64bd4dad07815c4a281ab1a1b50f Mon Sep 17 00:00:00 2001
From: Vuko Vukcevic <v.vukcevic@wikki.co.uk>
Date: Tue, 20 Mar 2018 17:27:01 +0100
Subject: [PATCH 02/15] Improvements to polyhedralRefinement

1. Extended protection for split points to nBufferLayers + 2 buffer layers in
   order to protect cells from unrefinement when we do refinement and
   unrefinement in a same go
2. Added separate entry for unrefineInterval along refineInterval to control
   refinement and unrefinement separately
---
 .../polyhedralRefinement.C                    |  12 +-
 .../dynamicPolyRefinementFvMesh.C             | 108 +++++++++++++-----
 .../dynamicPolyRefinementFvMesh.H             |   3 +
 3 files changed, 89 insertions(+), 34 deletions(-)

diff --git a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/polyhedralRefinement/polyhedralRefinement.C b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/polyhedralRefinement/polyhedralRefinement.C
index 2d75882b6..9bd53fae3 100644
--- a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/polyhedralRefinement/polyhedralRefinement.C
+++ b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/polyhedralRefinement/polyhedralRefinement.C
@@ -3555,9 +3555,7 @@ void Foam::polyhedralRefinement::setCellsToRefine
     // Transfer the contents into the data member (ordinary list)
     cellsToRefine_.transfer(cellsToRefineDynamic);
 
-    Info<< "polyhedralRefinement::setCellsToRefine"
-        << "(const labelList& refinementCellCandidates)" << nl
-        << "Selected " << returnReduce(cellsToRefine_.size(), sumOp<label>())
+    Info<< "Selected " << returnReduce(cellsToRefine_.size(), sumOp<label>())
         << " cells to refine." << endl;
 }
 
@@ -3709,9 +3707,9 @@ void Foam::polyhedralRefinement::setSplitPointsToUnrefine
     }
 
     // Extend protected cells across points using a specified number of buffer
-    // layers + 1 in order to stay far away from cells that are going to be
+    // layers + 2 in order to stay far away from cells that are going to be
     // refined
-    for (label i = 0; i < nBufferLayers_ + 1; ++i)
+    for (label i = 0; i < nBufferLayers_ + 2; ++i)
     {
         extendMarkedCellsAcrossPoints(protectedCell);
     }
@@ -3843,9 +3841,7 @@ void Foam::polyhedralRefinement::setSplitPointsToUnrefine
     // Transfer the contents into the data member (ordinary list)
     splitPointsToUnrefine_.transfer(splitPointsToUnrefineDynamic);
 
-    Info<< "polyhedralRefinement::setSplitPointsToUnrefine"
-        << "(const labelList& unrefinementPointCandidates)" << nl
-        << "Selected "
+    Info<< "Selected "
         << returnReduce(splitPointsToUnrefine_.size(), sumOp<label>())
         << " split points to unrefine." << endl;
 }
diff --git a/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicPolyRefinementFvMesh.C b/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicPolyRefinementFvMesh.C
index d3165ecef..62ebc8d3b 100644
--- a/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicPolyRefinementFvMesh.C
+++ b/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicPolyRefinementFvMesh.C
@@ -49,13 +49,23 @@ addToRunTimeSelectionTable
 
 void Foam::dynamicPolyRefinementFvMesh::readDict()
 {
-    // Read and check refinement interval
+    // Read and check refinement and unrefinement intervals
     refineInterval_ = readLabel(refinementDict_.lookup("refineInterval"));
     if (refineInterval_ < 1)
     {
         FatalErrorIn("dynamicPolyRefinementFvMesh::readDict()")
             << "Illegal refineInterval found: " << refineInterval_ << nl
-            << "The refineInterval controls the refinement/unrefinement"
+            << "The refineInterval controls the refinement"
+            << " trigerring within a certain time step and should be > 0"
+            << exit(FatalError);
+    }
+    
+    unrefineInterval_ = readLabel(refinementDict_.lookup("unrefineInterval"));
+    if (refineInterval_ < 1)
+    {
+        FatalErrorIn("dynamicPolyRefinementFvMesh::readDict()")
+            << "Illegal unrefineInterval found: " << refineInterval_ << nl
+            << "The unrefineInterval controls the unrefinement"
             << " trigerring within a certain time step and should be > 0"
             << exit(FatalError);
     }
@@ -86,6 +96,7 @@ Foam::dynamicPolyRefinementFvMesh::dynamicPolyRefinementFvMesh
         ).subDict(typeName + "Coeffs")
     ),
     refineInterval_(readLabel(refinementDict_.lookup("refineInterval"))),
+    unrefineInterval_(readLabel(refinementDict_.lookup("unrefineInterval"))),
     curTimeIndex_(-1),
 
     refinementSelectionPtr_(refinementSelection::New(*this, refinementDict_))
@@ -128,13 +139,22 @@ bool Foam::dynamicPolyRefinementFvMesh::update()
 
     // Performing refinement/unrefinement when:
     // 1. We are at the first time step
-    // 2. Time step is a multiplier of specified refineInterval
-    // 3. Only once per time step
+    // 2. Only once per time step
+    // 3. Time step is a multiplier of specified refineInterval or
+    //    unrefineInterval
+
+    // Get time index
+    const label timeID = time().timeIndex();
+
+    // Check whether to perform refinement and/or unrefinement
+    const bool performRefinement = timeID % refineInterval_ == 0;
+    const bool performUnrefinement = timeID % unrefineInterval_ == 0;
+
     if
     (
-        time().timeIndex() > 0
-     && time().timeIndex() % refineInterval_ == 0
-     && curTimeIndex_ < time().timeIndex()
+        timeID > 0
+     && curTimeIndex_ < timeID
+     && (performRefinement || performUnrefinement)
     )
     {
         // Update current time index to skip multiple topo changes per single
@@ -145,23 +165,47 @@ bool Foam::dynamicPolyRefinementFvMesh::update()
         polyhedralRefinement& polyRefModifier =
             refCast<polyhedralRefinement>(topoChanger_[0]);
 
-        // Get refinement candidates from refinement selection algorithm. Note:
-        // return type is Xfer<labelList> so there's no copying
-        const labelList refCandidates
-        (
-            refinementSelectionPtr_->refinementCellCandidates()
-        );
+        // Create empty list for refinement candidates
+        labelList refCandidates;
+
+        // Collect refinement candidates from refinement selection algorithm in
+        // case the refinement should be performed in this time step
+        if (performRefinement)
+        {
+            // Note: return type is Xfer<labelList> so there's no copying (two
+            // transfers are occuring)
+            refCandidates.transfer
+            (
+                refinementSelectionPtr_->refinementCellCandidates()()
+            );
+        }
+        else
+        {
+            Info<< "Skipping refinement for this time-step..." << endl;
+        }
 
         // Set cells to refine. Note: polyhedralRefinement ensures that face and
         // point consistent refinement is performed
         polyRefModifier.setCellsToRefine(refCandidates);
 
-        // Get unrefinement point candidates from refinement selection
-        // algorithm. Note: return type is Xfer<labelList> so there's no copying
-        const labelList unrefCandidates
-        (
-            refinementSelectionPtr_->unrefinementPointCandidates()
-        );
+        // Create empty list for unrefinement candidates
+        labelList unrefCandidates;
+
+        // Collect unrefinement candidates from refinement selection algorithm
+        // in case the unrefinement should be performed in this time step
+        if (performUnrefinement)
+        {
+            // Note: return type is Xfer<labelList> so there's no copying (two
+            // transfers are occuring)
+            unrefCandidates.transfer
+            (
+                refinementSelectionPtr_->unrefinementPointCandidates()()
+            );
+        }
+        else
+        {
+            Info<< "Skipping unrefinement for this time-step..." << endl;
+        }
 
         // Set split points to unrefine around.
         // Notes:
@@ -191,13 +235,25 @@ bool Foam::dynamicPolyRefinementFvMesh::update()
         // Perform refinement and unrefinement in one go
         autoPtr<mapPolyMesh> topoChangeMap = topoChanger_.changeMesh();
 
-        // Output cell balance
-        Info<< "Successfully performed polyhedral refinement. "
-            << "Changed from "
-            << returnReduce(topoChangeMap->nOldCells(), sumOp<label>())
-            << " to "
-            << returnReduce(topoChangeMap->cellMap().size(), sumOp<label>())
-            << " cells." << endl;
+        // Output cell balance if the topo change has been performed
+        const label nOldCells =
+            returnReduce(topoChangeMap->nOldCells(), sumOp<label>());
+        const label sizeCellMap =
+            returnReduce(topoChangeMap->cellMap().size(), sumOp<label>());
+
+        // If the size of cell map is different than zero, we actually performed
+        // some topo changes
+        if (sizeCellMap)
+        {
+            Info<< "Successfully performed polyhedral refinement. "
+                << "Changed from " << nOldCells << " to " << sizeCellMap
+                << " cells." << endl;
+        }
+        else
+        {
+            Info<< "Refinement/unrefinement not performed in this time step "
+                << "since no cells were selected." << endl;
+        }
 
         return topoChangeMap->morphing();
     }
diff --git a/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicPolyRefinementFvMesh.H b/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicPolyRefinementFvMesh.H
index 8ff1d53fd..ad6570683 100644
--- a/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicPolyRefinementFvMesh.H
+++ b/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicPolyRefinementFvMesh.H
@@ -65,6 +65,9 @@ class dynamicPolyRefinementFvMesh
         //- Refinement interval
         label refineInterval_;
 
+        //- Unrefinement interval
+        label unrefineInterval_;
+
         //- Current time index (helper variable to skip multiple topo changes in
         //  a single time step)
         label curTimeIndex_;

From 26d346cf7117a44441a707def851460cb2fa7605 Mon Sep 17 00:00:00 2001
From: Vuko Vukcevic <v.vukcevic@wikki.co.uk>
Date: Thu, 22 Mar 2018 09:49:51 +0100
Subject: [PATCH 03/15] Additional refinement selection algorithm for
 polyhedral AMR

minCellSizeRefinement selects all cells with cell size defined as cube root of
cell volume is larger than specified cell size. Therefore, it assumes that the
cells are cubes or close to cubes.
---
 src/dynamicMesh/topoChangerFvMesh/Make/files  |   1 +
 .../minCellSizeRefinement.C                   | 128 ++++++++++++++++++
 .../minCellSizeRefinement.H                   | 111 +++++++++++++++
 3 files changed, 240 insertions(+)
 create mode 100644 src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/refinementSelection/minCellSizeRefinement/minCellSizeRefinement.C
 create mode 100644 src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/refinementSelection/minCellSizeRefinement/minCellSizeRefinement.H

diff --git a/src/dynamicMesh/topoChangerFvMesh/Make/files b/src/dynamicMesh/topoChangerFvMesh/Make/files
index 8b031c30f..1fd494240 100644
--- a/src/dynamicMesh/topoChangerFvMesh/Make/files
+++ b/src/dynamicMesh/topoChangerFvMesh/Make/files
@@ -15,6 +15,7 @@ dynamicPolyRefinementFvMesh/dynamicPolyRefinementFvMesh.C
 dynamicPolyRefinementFvMesh/refinementSelection/refinementSelection/refinementSelection.C
 dynamicPolyRefinementFvMesh/refinementSelection/fieldBoundsRefinement/fieldBoundsRefinement.C
 dynamicPolyRefinementFvMesh/refinementSelection/minCellVolumeRefinement/minCellVolumeRefinement.C
+dynamicPolyRefinementFvMesh/refinementSelection/minCellSizeRefinement/minCellSizeRefinement.C
 dynamicPolyRefinementFvMesh/refinementSelection/minPatchDistanceRefinement/minPatchDistanceRefinement.C
 dynamicPolyRefinementFvMesh/refinementSelection/compositeRefinementSelection/compositeRefinementSelection.C
 
diff --git a/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/refinementSelection/minCellSizeRefinement/minCellSizeRefinement.C b/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/refinementSelection/minCellSizeRefinement/minCellSizeRefinement.C
new file mode 100644
index 000000000..b4164bec3
--- /dev/null
+++ b/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/refinementSelection/minCellSizeRefinement/minCellSizeRefinement.C
@@ -0,0 +1,128 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Author
+    Vuko Vukcevic, Wikki Ltd.  All rights reserved.
+
+\*---------------------------------------------------------------------------*/
+
+#include "minCellSizeRefinement.H"
+#include "addToRunTimeSelectionTable.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(minCellSizeRefinement, 0);
+addToRunTimeSelectionTable
+(
+    refinementSelection,
+    minCellSizeRefinement,
+    dictionary
+);
+
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::minCellSizeRefinement::minCellSizeRefinement
+(
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+:
+    refinementSelection(mesh, dict),
+    minDelta_(readScalar(coeffDict().lookup("minCellSize")))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor* * * * * * * * * * * * * * * * //
+
+Foam::minCellSizeRefinement::~minCellSizeRefinement()
+{}
+
+
+// * * * * * * * * * * * * Public Member Functions * * * * * * * * * * * * * //
+
+Foam::Xfer<Foam::labelList>
+Foam::minCellSizeRefinement::refinementCellCandidates() const
+{
+    // Get cell sizes: cube root of cell volume (assuming cube cells)
+    const scalarField cellSize = pow(mesh().V().field(), 1.0/3.0);
+
+    // Create storage for collection of cells. Assume that almost all of the
+    // cells will be marked to prevent excessive resizing.
+    dynamicLabelList refinementCandidates(mesh().nCells());
+
+    // Loop through cells and collect refinement candidates
+    forAll (cellSize, cellI)
+    {
+        if (cellSize[cellI] > minDelta_)
+        {
+            // Cell is larger than the specified minimum, append cell for
+            // potential refinement
+            refinementCandidates.append(cellI);
+        }
+    }
+
+    // Print out some information
+    Info<< "Selection algorithm " << type() << " selected "
+        << returnReduce(refinementCandidates.size(), sumOp<label>())
+        << " cells as refinement candidates."
+        << endl;
+
+    // Return the list in the Xfer container to prevent copying
+    return refinementCandidates.xfer();
+}
+
+
+Foam::Xfer<Foam::labelList>
+Foam::minCellSizeRefinement::unrefinementPointCandidates() const
+{
+    // Mark all points as unrefinement candidates since only split points may be
+    // considered for actual unrefinement and since this refinement criterion
+    // will be usually used in combination with others. VV, 15/Mar/2018.
+
+    // All points are unrefinement candidates
+    labelList unrefinementCandidates(mesh().nPoints());
+
+    forAll (unrefinementCandidates, pointI)
+    {
+        unrefinementCandidates[pointI] = pointI;
+    }
+
+    // Print out some information
+    Info<< "Selection algorithm " << type() << " selected "
+        << returnReduce(unrefinementCandidates.size(), sumOp<label>())
+        << " points as unrefinement candidates."
+        << endl;
+
+    // Return the list in the Xfer container to prevent copying
+    return unrefinementCandidates.xfer();
+}
+
+
+// ************************************************************************* //
diff --git a/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/refinementSelection/minCellSizeRefinement/minCellSizeRefinement.H b/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/refinementSelection/minCellSizeRefinement/minCellSizeRefinement.H
new file mode 100644
index 000000000..d030d8dc6
--- /dev/null
+++ b/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/refinementSelection/minCellSizeRefinement/minCellSizeRefinement.H
@@ -0,0 +1,111 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::minCellSizeRefinement
+
+Description
+    Selection of refinement cells based on a minimum cell size. Assumes that
+    cells are mostly cubic and compares the minimum cell size with V^1/3, where
+    V is the cell volume. If the cell size is larger than the specified one,
+    cell gets selected for refinement.
+
+SourceFiles
+    minCellSizeRefinement.C
+
+Author
+    Vuko Vukcevic, Wikki Ltd.  All rights reserved.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef minCellSizeRefinement_H
+#define minCellSizeRefinement_H
+
+#include "refinementSelection.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                   Class minCellSizeRefinement Declaration
+\*---------------------------------------------------------------------------*/
+
+class minCellSizeRefinement
+:
+    public refinementSelection
+{
+    // Private data
+
+        //- Minimum cell size
+        scalar minDelta_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        minCellSizeRefinement(const minCellSizeRefinement&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const minCellSizeRefinement&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("minCellSizeRefinement");
+
+
+    // Constructors
+
+        //- Construct from components
+        minCellSizeRefinement(const fvMesh& mesh, const dictionary& dict);
+
+
+    //- Destructor
+    virtual ~minCellSizeRefinement();
+
+
+    // Member Functions
+
+        // Selection of refinement/unrefinement candidates
+
+            //- Return transferable list of cells to refine
+            virtual Xfer<labelList> refinementCellCandidates() const;
+
+            //- Return transferable list of split points to unrefine
+            virtual Xfer<labelList> unrefinementPointCandidates() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

From 4a4056f9f572aff39521afc181df00a7f22239ad Mon Sep 17 00:00:00 2001
From: Vuko Vukcevic <v.vukcevic@wikki.co.uk>
Date: Thu, 22 Mar 2018 14:24:44 +0100
Subject: [PATCH 04/15] Improvements to polyhedral AMR

1. Removed redundant (and wrong) check when doing pointConsistentUnrefinement
2. Added separate entry for nRefinementBufferLayers and
   nUnrefinementBufferLayers
---
 .../polyhedralRefinement.C                    | 136 +++++++++++++-----
 .../polyhedralRefinement.H                    |   8 +-
 .../dynamicMeshDict                           |  12 +-
 3 files changed, 121 insertions(+), 35 deletions(-)

diff --git a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/polyhedralRefinement/polyhedralRefinement.C b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/polyhedralRefinement/polyhedralRefinement.C
index 9bd53fae3..42e9d83dc 100644
--- a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/polyhedralRefinement/polyhedralRefinement.C
+++ b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/polyhedralRefinement/polyhedralRefinement.C
@@ -3018,8 +3018,8 @@ Foam::label Foam::polyhedralRefinement::pointConsistentRefinement
             {
                 if (!willBeRefined)
                 {
-                    // Cell has not been marked for refinement, mark the cell for
-                    // refinement and increment the counter
+                    // Cell has not been marked for refinement, mark the cell
+                    // for refinement and increment the counter
                     willBeRefined = true;
                     ++nAddCells;
                 }
@@ -3029,10 +3029,12 @@ Foam::label Foam::polyhedralRefinement::pointConsistentRefinement
                     (
                         "label polyhedralRefinement::pointConsistentRefinement"
                         "(boolList cellsToRefine) const"
-                    )   << "Cell is marked for refinement, but the 4:1 point"
-                        << " consistency cannot be ensured." << nl
-                        << "Something went wrong before this step."
-                        << endl;
+                    )   << "Cell: " << curCellI << " is marked for refinement,"
+                        << "but the 4:1 consistency cannot be ensured." << nl
+                        << "This is probably because the refinement and "
+                        << "unrefinement regions are very close." << nl
+                        << "Try increasing nUnrefinementBufferLayers. "
+                        << abort(FatalError);
                 }
             }
         }
@@ -3090,7 +3092,10 @@ Foam::label Foam::polyhedralRefinement::faceConsistentUnrefinement
                     << "Owner: " << own << ", neighbour: " << nei
                     << nl
                     << "Owner level: " << ownLevel
-                    << ", neighbour level: " << neiLevel
+                    << ", neighbour level: " << neiLevel << nl
+                    << "This is probably because the refinement and "
+                    << "unrefinement regions are very close." << nl
+                    << "Try increasing nUnrefinementBufferLayers. "
                     << abort(FatalError);
             }
 
@@ -3115,7 +3120,10 @@ Foam::label Foam::polyhedralRefinement::faceConsistentUnrefinement
                     << "Owner: " << own << ", neighbour: " << nei
                     << nl
                     << "Owner level: " << ownLevel
-                    << ", neighbour level: " << neiLevel
+                    << ", neighbour level: " << neiLevel << nl
+                    << "This is probably because the refinement and "
+                    << "unrefinement regions are very close." << nl
+                    << "Try increasing nUnrefinementBufferLayers. "
                     << abort(FatalError);
             }
 
@@ -3168,7 +3176,10 @@ Foam::label Foam::polyhedralRefinement::faceConsistentUnrefinement
                     << "Owner: " << own
                     << nl
                     << "Owner level: " << curOwnLevel
-                    << ", neighbour level: " << neiLevel[i]
+                    << ", neighbour level: " << neiLevel[i] << nl
+                    << "This is probably because the refinement and "
+                    << "unrefinement regions are very close." << nl
+                    << "Try increasing nUnrefinementBufferLayers. "
                     << abort(FatalError);
             }
 
@@ -3263,18 +3274,25 @@ Foam::label Foam::polyhedralRefinement::pointConsistentUnrefinement
                     willBeUnrefined = false;
                     ++nRemCells;
                 }
-                else
-                {
-                    FatalErrorIn
-                    (
-                        "label polyhedralRefinement::"
-                        "pointConsistentUnrefinement"
-                        "(boolList cellsToRefine) const"
-                    )   << "Cell is not marked for unrefinement, but the 4:1"
-                        << " point consistency cannot be ensured." << nl
-                        << "Something went wrong before this step."
-                        << endl;
-                }
+                // Note: I am pretty sure that this check is redundant and
+                // causes unnecessary issues when more than 2 maximum refinement
+                // levels are used. In case of problems, take a detailed look
+                // here. VV, 22/Mar/2018.
+//                else
+//                {
+//                    FatalErrorIn
+//                    (
+//                        "label polyhedralRefinement::"
+//                        "pointConsistentUnrefinement"
+//                        "(boolList cellsToRefine) const"
+//                    )   << "Cell: " << curCellI << " is not marked for "
+//                        << "unrefinement, but the 4:1"
+//                        << " point consistency cannot be ensured." << nl
+//                        << "This is probably because the refinement and "
+//                        << "unrefinement regions are very close." << nl
+//                        << "Try increasing nUnrefinementBufferLayers. "
+//                        << abort(FatalError);
+//                }
             }
         }
     }
@@ -3332,7 +3350,14 @@ Foam::polyhedralRefinement::polyhedralRefinement
     (
         dict.lookupOrDefault<Switch>("pointBasedConsistency", true)
     ),
-    nBufferLayers_(readScalar(dict.lookup("nBufferLayers")))
+    nRefinementBufferLayers_
+    (
+        readScalar(dict.lookup("nRefinementBufferLayers"))
+    ),
+    nUnrefinementBufferLayers_
+    (
+        readScalar(dict.lookup("nUnrefinementBufferLayers"))
+    )
 {
     // Calculate level 0 edge length
     calcLevel0EdgeLength();
@@ -3424,8 +3449,8 @@ Foam::polyhedralRefinement::polyhedralRefinement
             << endl;
     }
 
-    // Check number of buffer layers
-    if (nBufferLayers_ < 0)
+    // Check number of refinement buffer layers
+    if (nRefinementBufferLayers_ < 0)
     {
         FatalErrorIn
         (
@@ -3436,11 +3461,53 @@ Foam::polyhedralRefinement::polyhedralRefinement
             "\n    const label index,"
             "\n    const polyTopoChanger& mme"
             "\n)"
-        )   << "Negative nBufferLayers specified."
+        )   << "Negative nRefinementBufferLayers specified."
             << nl
             << "This is not allowed."
             << abort(FatalError);
     }
+
+    // Check number of unrefinement buffer layers
+    if (nUnrefinementBufferLayers_ < 0)
+    {
+        FatalErrorIn
+        (
+            "polyhedralRefinement::polyhedralRefinement"
+            "\n("
+            "\n    const word& name,"
+            "\n    const dictionary& dict,"
+            "\n    const label index,"
+            "\n    const polyTopoChanger& mme"
+            "\n)"
+        )   << "Negative nUnrefinementBufferLayers specified."
+            << nl
+            << "This is not allowed."
+            << abort(FatalError);
+    }
+
+    // Check whether the number of unrefinement buffer layers is smaller than
+    // number of refinement buffer layers + 2
+    if (nUnrefinementBufferLayers_ < nRefinementBufferLayers_ + 2)
+    {
+        WarningIn
+        (
+            "polyhedralRefinement::polyhedralRefinement"
+            "\n("
+            "\n    const word& name,"
+            "\n    const dictionary& dict,"
+            "\n    const label index,"
+            "\n    const polyTopoChanger& mme"
+            "\n)"
+        )   << "Using " << nUnrefinementBufferLayers_
+            << " unrefinement buffer layers and " << nRefinementBufferLayers_
+            << " refinement buffer layers."
+            << nl
+            << "Make sure that the number of unrefinement buffer layers is "
+            << "at least nRefinementBufferLayers + 2" << nl
+            << "in order to avoid problems with point level inconsistency when "
+            << "refinement and unrefinement are performed in same iteration."
+            << endl;
+    }
 }
 
 
@@ -3495,8 +3562,9 @@ void Foam::polyhedralRefinement::setCellsToRefine
         }
     }
 
-    // Extend cells across faces using a specified number of buffer layers
-    for (label i = 0; i < nBufferLayers_; ++i)
+    // Extend cells across faces using a specified number of refinement buffer
+    // layers
+    for (label i = 0; i < nRefinementBufferLayers_; ++i)
     {
         extendMarkedCellsAcrossFaces(refineCell);
     }
@@ -3706,10 +3774,9 @@ void Foam::polyhedralRefinement::setSplitPointsToUnrefine
         protectedCell[cellsToRefine_[i]] = true;
     }
 
-    // Extend protected cells across points using a specified number of buffer
-    // layers + 2 in order to stay far away from cells that are going to be
-    // refined
-    for (label i = 0; i < nBufferLayers_ + 2; ++i)
+    // Extend protected cells across points using a specified number of
+    // unrefinement buffer layers
+    for (label i = 0; i < nUnrefinementBufferLayers_ + 2; ++i)
     {
         extendMarkedCellsAcrossPoints(protectedCell);
     }
@@ -4060,7 +4127,8 @@ void Foam::polyhedralRefinement::write(Ostream& os) const
         << maxCells_ << nl
         << maxRefinementLevel_ << nl
         << pointBasedConsistency_ << nl
-        << nBufferLayers_ << endl;
+        << nRefinementBufferLayers_ << nl
+        << nUnrefinementBufferLayers_ << endl;
 }
 
 
@@ -4079,7 +4147,9 @@ void Foam::polyhedralRefinement::writeDict(Ostream& os) const
         << token::END_STATEMENT << nl
         << "    pointBasedConsistency " << pointBasedConsistency_
         << token::END_STATEMENT << nl
-        << "    nBufferLayers " << nBufferLayers_
+        << "    nRefinementBufferLayers " << nRefinementBufferLayers_
+        << token::END_STATEMENT << nl
+        << "    nUnrefinementBufferLayers " << nUnrefinementBufferLayers_
         << token::END_STATEMENT << nl
         << "    active " << active()
         << token::END_STATEMENT << nl
diff --git a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/polyhedralRefinement/polyhedralRefinement.H b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/polyhedralRefinement/polyhedralRefinement.H
index 19251e98b..d4735760d 100644
--- a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/polyhedralRefinement/polyhedralRefinement.H
+++ b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/polyhedralRefinement/polyhedralRefinement.H
@@ -128,7 +128,13 @@ private:
             Switch pointBasedConsistency_;
 
             //- Number of buffer layers for refinement
-            label nBufferLayers_;
+            label nRefinementBufferLayers_;
+
+            //- Number of buffer layers for unrefinement, controlling how far
+            //  the unrefinement region needs to be from current refinement
+            //  region. Therefore, this should always be at least
+            //  nRefinementBufferLayers + 2 to avoid level inconsistencies
+            label nUnrefinementBufferLayers_;
 
 
     // Private Member Functions
diff --git a/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicMeshDict b/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicMeshDict
index 240fb9361..e86bd90d2 100644
--- a/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicMeshDict
+++ b/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicMeshDict
@@ -20,8 +20,13 @@ dynamicFvMesh   dynamicPolyRefinementFvMesh;
 dynamicPolyRefinementFvMeshCoeffs
 {
     // Dynamic mesh procedure controls
+
+    // Refine every refineInterval step
     refineInterval   1;
 
+    // Unrefine every unrefineInterval step
+    unrefineInterval 1;
+
     // Refinement selection criteria
     refinementSelection
     {
@@ -47,7 +52,12 @@ dynamicPolyRefinementFvMeshCoeffs
     maxRefinementLevel   3;
 
     // Number of buffer layers between refinement levels
-    nBufferLayers        1;
+    nRefinementBufferLayers        1;
+
+    // Number of buffer layers for unrefinement in order to run away from the
+    // region that is getting refined at the same time in order to avoid point
+    // level inconsistencies
+    nUnrefinementBufferLayers      4;
 
     // Whether to use point based consistency check. Needed when one allows more
     // than 2 refinement levels (automatically switched on in that case)

From 9ee5d4957523d7cd0db2e9fedb19cafe5e65a02d Mon Sep 17 00:00:00 2001
From: Vuko Vukcevic <v.vukcevic@wikki.co.uk>
Date: Tue, 17 Apr 2018 12:24:52 +0200
Subject: [PATCH 05/15] Updates to polyhedralRefinement and
 dynamicPolyRefinementFvMesh

1. Switched from pointBasedConsistency to edgeBasedConsistency in order to allow
   more than 2 refinement levels in a correct way (hopefully),
2. Note on parallelisation: getting edge cell addressing across processor
   boundaries is tricky, but if we use more stringent 1:1 face refinement across
   coupled boundaries, then the edge based consistency becomes a local operation
   (at least I'm 99% sure it does: hard to think of all variants with arbitrary
   polyhedra),
3. Added an option to skip unrefinement if both refinement/unrefinement are
   triggered for a single step. This is switched off by default since performing
   refinement/unrefinement in a single go should work without issues.
---
 .../polyhedralRefinement.C                    | 354 +++++++++---------
 .../polyhedralRefinement.H                    |  12 +-
 .../dynamicMeshDict                           |  12 +-
 .../dynamicPolyRefinementFvMesh.C             |  19 +-
 .../dynamicPolyRefinementFvMesh.H             |   4 +
 5 files changed, 219 insertions(+), 182 deletions(-)

diff --git a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/polyhedralRefinement/polyhedralRefinement.C b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/polyhedralRefinement/polyhedralRefinement.C
index 42e9d83dc..fd68a8353 100644
--- a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/polyhedralRefinement/polyhedralRefinement.C
+++ b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/polyhedralRefinement/polyhedralRefinement.C
@@ -2930,15 +2930,19 @@ Foam::label Foam::polyhedralRefinement::faceConsistentRefinement
         const label curOwnLevel =
             cellsToRefine[own] ? cellLevel_[own] + 1 : cellLevel_[own];
 
-        if (neiLevel[i] > (curOwnLevel + 1))
+        // Note: we are using more stringent 1:1 consistency across coupled
+        // boundaries in order to simplify handling of edge based consistency
+        // checks for parallel runs
+        if (neiLevel[i] > curOwnLevel)
         {
-            // Neighbour level is higher than owner level + 1, owner must be
+            // Neighbour level is higher than owner level, owner must be
             // marked for refinement
             cellsToRefine[own] = true;
             ++nAddCells;
         }
+
         // Note: other possibility (that owner level is higher than neighbour
-        // level + 1) is taken into account on the other side automatically
+        // level) is taken into account on the other side automatically
     }
 
     // Return number of added cells
@@ -2946,7 +2950,7 @@ Foam::label Foam::polyhedralRefinement::faceConsistentRefinement
 }
 
 
-Foam::label Foam::polyhedralRefinement::pointConsistentRefinement
+Foam::label Foam::polyhedralRefinement::edgeConsistentRefinement
 (
     boolList& cellsToRefine
 ) const
@@ -2954,92 +2958,75 @@ Foam::label Foam::polyhedralRefinement::pointConsistentRefinement
     // Count number of cells that will be added
     label nAddCells = 0;
 
-    // Maximum surrounding cell refinement level for each point
-    labelList maxRefLevel(mesh_.nPoints(), 0);
+    // Algorithm: loop over all edges and visit all unique cell pairs sharing
+    // this particular edge. Then, ensure 2:1 edge consistency by marking
+    // cell with lower level for refinement
 
-    // Get point cells
-    const labelListList& meshPointCells = mesh_.pointCells();
+    // Get edge cells
+    const labelListList& meshEdgeCells = mesh_.edgeCells();
 
-    // Loop through all points and collect maximum point level for each point
-    forAll (maxRefLevel, pointI)
+    // Loop through all mesh edges
+    forAll (meshEdgeCells, edgeI)
     {
-        // Get the cells for this point
-        const labelList& curCells = meshPointCells[pointI];
+        // Get current edge cells
+        const labelList& curEdgeCells = meshEdgeCells[edgeI];
 
-        // Get reference to maximum pooint level for this point
-        label& curMaxPointLevel = maxRefLevel[pointI];
-
-        // Find maximum refinement level for this point
-        forAll (curCells, i)
+        // Loop through all edge cells
+        forAll (curEdgeCells, i)
         {
-            // Get cell index and "future" cell level
-            const label& curCellI = curCells[i];
-            const label curCellLevel =
-                cellsToRefine[curCellI]
-              ? cellLevel_[curCellI] + 1
-              : cellLevel_[curCellI];
+            // Get first cell index
+            const label& cellI = curEdgeCells[i];
 
-            // Update maximum point level if the curCellLevel is higher
-            curMaxPointLevel = max(curMaxPointLevel, curCellLevel);
-        }
-    }
-
-    // Sync maximum refinement level across coupled boundaries
-    syncTools::syncPointList
-    (
-        mesh_,
-        maxRefLevel,
-        maxEqOp<label>(),
-        0,   // Null value
-        true // Apply separation for parallel cyclics
-    );
-
-    // Now that the levels are synced, go through all points and add cells to
-    // refine
-    forAll (maxRefLevel, pointI)
-    {
-        // Get the cells for this point
-        const labelList& curCells = meshPointCells[pointI];
-
-        // Loop through these point cells and set cells for refinement which
-        // would end up having refinement level smaller than maximum level - 1
-        forAll (curCells, i)
-        {
-            // Get cell index, reference to refinement flag and "future" cell
-            // level
-            const label& curCellI = curCells[i];
-            bool& willBeRefined = cellsToRefine[curCellI];
-            const label curCellLevel =
-                willBeRefined
-              ? cellLevel_[curCellI] + 1
-              : cellLevel_[curCellI];
-
-            if (curCellLevel < maxRefLevel[pointI] - 1)
+            // Loop through remaining edge cells
+            for (label j = i + 1; j < curEdgeCells.size(); ++j)
             {
-                if (!willBeRefined)
+                // Get second cell index
+                const label& cellJ = curEdgeCells[j];
+
+                // Get levels of the two cells. If the cell is marked for
+                // refinement, the level is current level + 1, otherwise it is
+                // equal to the current level
+
+                // Note: cellsToRefine flag for both cellI and cellJ might
+                // change, this is why we need to recalculate cellI level here
+                const label cellILevel =
+                    cellsToRefine[cellI]
+                  ? cellLevel_[cellI] + 1
+                  : cellLevel_[cellI];
+
+                const label cellJLevel =
+                    cellsToRefine[cellJ]
+                  ? cellLevel_[cellJ] + 1
+                  : cellLevel_[cellJ];
+
+                if (cellILevel > cellJLevel + 1)
                 {
-                    // Cell has not been marked for refinement, mark the cell
-                    // for refinement and increment the counter
-                    willBeRefined = true;
+                    // Level of cellI is higher than level of cellJ + 1, cellJ
+                    // must be marked for refinement
+                    cellsToRefine[cellJ] = true;
                     ++nAddCells;
                 }
-                else
+                else if (cellJLevel > cellILevel + 1)
                 {
-                    FatalErrorIn
-                    (
-                        "label polyhedralRefinement::pointConsistentRefinement"
-                        "(boolList cellsToRefine) const"
-                    )   << "Cell: " << curCellI << " is marked for refinement,"
-                        << "but the 4:1 consistency cannot be ensured." << nl
-                        << "This is probably because the refinement and "
-                        << "unrefinement regions are very close." << nl
-                        << "Try increasing nUnrefinementBufferLayers. "
-                        << abort(FatalError);
+                    // Level of cellJ is higher than level of cellI + 1, cellI
+                    // must be marked for refinement
+                    cellsToRefine[cellI] = true;
+                    ++nAddCells;
                 }
             }
         }
     }
 
+    // Note: in order to avoid very difficult and time-consuming parallelisation
+    // of edge cell connectivity and edge cell values, we enforce a more
+    // stringent face-based consistency across processor boundaries. Basically,
+    // if a face-based consistency of 1:1 (not 2:1 as for ordinary faces) is
+    // ensured, the edge-based consistency becomes a local operation (I'm not
+    // 100% sure to be honest since there are countless variants when dealing
+    // with arbitrary polyhedral cells).
+    // See faceConsistentRefinement for details. VV, 17/Apr/2018.
+
+    // Return number of added cells
     return nAddCells;
 }
 
@@ -3074,7 +3061,7 @@ Foam::label Foam::polyhedralRefinement::faceConsistentUnrefinement
         const label neiLevel =
             cellsToUnrefine[nei] ? cellLevel_[nei] - 1 : cellLevel_[nei];
 
-        if (ownLevel < (neiLevel -1))
+        if (ownLevel < (neiLevel - 1))
         {
             // Owner level is smaller than neighbour level - 1, we must not
             // unrefine owner
@@ -3158,9 +3145,12 @@ Foam::label Foam::polyhedralRefinement::faceConsistentUnrefinement
         const label curOwnLevel =
             cellsToUnrefine[own] ? cellLevel_[own] - 1 : cellLevel_[own];
 
-        if (curOwnLevel < (neiLevel[i] - 1))
+        // Note: we are using more stringent 1:1 consistency across coupled
+        // boundaries in order to simplify handling of edge based consistency
+        // checkes for parallel runs
+        if (curOwnLevel < neiLevel[i])
         {
-            // Owner level is smaller than neighbour level - 1, we must not
+            // Owner level is smaller than neighbour level, we must not
             // unrefine owner
 
             // Check whether the cell has not been marked for unrefinement
@@ -3188,7 +3178,7 @@ Foam::label Foam::polyhedralRefinement::faceConsistentUnrefinement
         }
 
         // Note: other possibility (that neighbour level is smaller than owner
-        // level - 1) is taken into account on the other side automatically
+        // level) is taken into account on the other side automatically
     }
 
     // Return number of local cells removed from unrefinement
@@ -3196,107 +3186,129 @@ Foam::label Foam::polyhedralRefinement::faceConsistentUnrefinement
 }
 
 
-Foam::label Foam::polyhedralRefinement::pointConsistentUnrefinement
+Foam::label Foam::polyhedralRefinement::edgeConsistentUnrefinement
 (
     boolList& cellsToUnrefine
 ) const
 {
-    // Count number of cells removed from unrefinement
+    // Count number of cells that will be removed
     label nRemCells = 0;
 
-    // Minimum cell refinement level for each point. Note: initialise with
-    // labelMax
-    labelList minRefLevel(mesh_.nPoints(), labelMax);
+    // Algorithm: loop over all edges and visit all unique cell pairs sharing
+    // this particular edge. Then, ensure 2:1 edge consistency by protecting the
+    // cell with lower level from unrefinement
 
-    // Get point cells
-    const labelListList& meshPointCells = mesh_.pointCells();
+    // Get edge cells
+    const labelListList& meshEdgeCells = mesh_.edgeCells();
 
-    // Loop through all points and collect minimum point level for each point
-    forAll (minRefLevel, pointI)
+    // Loop through all mesh edges
+    forAll (meshEdgeCells, edgeI)
     {
-        // Get the cell for this point
-        const labelList& curCells = meshPointCells[pointI];
+        // Get current edge cells
+        const labelList& curEdgeCells = meshEdgeCells[edgeI];
 
-        // Get reference to minimum point level for this point
-        label& curMinPointLevel = minRefLevel[pointI];
-
-        // Find minimum refinement level for this point
-        forAll (curCells, i)
+        // Loop through all edge cells
+        forAll (curEdgeCells, i)
         {
-            // Get cell index and "future" cell level
-            const label& curCellI = curCells[i];
-            const label curCellLevel =
-                cellsToUnrefine[curCellI]
-              ? cellLevel_[curCellI] - 1
-              : cellLevel_[curCellI];
+            // Get first cell index
+            const label& cellI = curEdgeCells[i];
 
-            // Update minimum point level if the curCellLevel is smaller
-            curMinPointLevel = min(curMinPointLevel, curCellLevel);
-        }
-    }
-
-    // Sync minimum refinement level across coupled boundaries
-    syncTools::syncPointList
-    (
-        mesh_,
-        minRefLevel,
-        minEqOp<label>(),
-        0,   // Null value
-        true // Apply separation for parallel cyclics
-    );
-
-    // Now that the levels are synced, go through all points and protect some
-    // cells from unrefinement
-    forAll (minRefLevel, pointI)
-    {
-        // Get the cells for this point
-        const labelList& curCells = meshPointCells[pointI];
-
-        // Loop through these point cells and protected cells from unrefinement
-        // which would end up having refinement level greater than level + 1
-        forAll (curCells, i)
-        {
-            // Get cell index, reference to unrefinement flag and "future" cell
-            // level
-            const label& curCellI = curCells[i];
-            bool& willBeUnrefined = cellsToUnrefine[curCellI];
-            const label curCellLevel =
-                willBeUnrefined
-              ? cellLevel_[curCellI] - 1
-              : cellLevel_[curCellI];
-
-            if (curCellLevel > minRefLevel[pointI] + 1)
+            // Loop through remaining edge cells
+            for (label j = i + 1; j < curEdgeCells.size(); ++j)
             {
-                if (willBeUnrefined)
+                // Get second cell index
+                const label& cellJ = curEdgeCells[j];
+
+                // Get levels of the two cells. If the cell is marked for
+                // unrefinement, the level is current level - 1, otherwise it is
+                // equal to the current level
+
+                // Note: cellsToUnrefine flag for both cellI and cellJ might
+                // change, this is why we need to recalculate cellI level here
+                const label cellILevel =
+                    cellsToUnrefine[cellI]
+                  ? cellLevel_[cellI] - 1
+                  : cellLevel_[cellI];
+
+                const label cellJLevel =
+                    cellsToUnrefine[cellJ]
+                  ? cellLevel_[cellJ] - 1
+                  : cellLevel_[cellJ];
+
+                if (cellILevel < cellJLevel - 1)
                 {
-                    // Cell has been marked for unrefinement, protect the cell
-                    // from unrefinement and increment the counter
-                    willBeUnrefined = false;
+                    // Level of cellI is smaller than level of cellJ - 1, cellI
+                    // must be protected from unrefinement
+
+                    // Check whether the cell has not been marked for
+                    // unrefinement
+                    if (!cellsToUnrefine[cellI])
+                    {
+                        FatalErrorIn
+                        (
+                            "label polyhedralRefinement::"
+                            "edgeConsistentUnrefinement"
+                            "(boolList& cellsToUnrefine)"
+                        )   << "Cell not marked for unrefinement, indicating a"
+                            << " previous unnoticed problem with unrefinement."
+                            << nl
+                            << "cellI: " << cellI << ", cellJ: " << cellJ
+                            << nl
+                            << "Level of cellI: " << cellILevel
+                            << ", level of cellJ: " << cellJLevel << nl
+                            << "This is probably because the refinement and "
+                            << "unrefinement regions are very close." << nl
+                            << "Try increasing nUnrefinementBufferLayers. "
+                            << abort(FatalError);
+                    }
+
+                    cellsToUnrefine[cellI] = false;
+                    ++nRemCells;
+                }
+                else if (cellJLevel < cellILevel - 1)
+                {
+                    // Level of cellJ is smaller than level of cellI - 1, cellJ
+                    // must be protected from unrefinement
+
+                    // Check whether the cell has not been marked for
+                    // unrefinement
+                    if (!cellsToUnrefine[cellJ])
+                    {
+                        FatalErrorIn
+                        (
+                            "label polyhedralRefinement::"
+                            "edgeConsistentUnrefinement"
+                            "(boolList& cellsToUnrefine)"
+                        )   << "Cell not marked for unrefinement, indicating a"
+                            << " previous unnoticed problem with unrefinement."
+                            << nl
+                            << "cellI: " << cellI << ", cellJ: " << cellJ
+                            << nl
+                            << "Level of cellI: " << cellILevel
+                            << ", level of cellJ: " << cellJLevel << nl
+                            << "This is probably because the refinement and "
+                            << "unrefinement regions are very close." << nl
+                            << "Try increasing nUnrefinementBufferLayers. "
+                            << abort(FatalError);
+                    }
+
+                    cellsToUnrefine[cellJ] = false;
                     ++nRemCells;
                 }
-                // Note: I am pretty sure that this check is redundant and
-                // causes unnecessary issues when more than 2 maximum refinement
-                // levels are used. In case of problems, take a detailed look
-                // here. VV, 22/Mar/2018.
-//                else
-//                {
-//                    FatalErrorIn
-//                    (
-//                        "label polyhedralRefinement::"
-//                        "pointConsistentUnrefinement"
-//                        "(boolList cellsToRefine) const"
-//                    )   << "Cell: " << curCellI << " is not marked for "
-//                        << "unrefinement, but the 4:1"
-//                        << " point consistency cannot be ensured." << nl
-//                        << "This is probably because the refinement and "
-//                        << "unrefinement regions are very close." << nl
-//                        << "Try increasing nUnrefinementBufferLayers. "
-//                        << abort(FatalError);
-//                }
             }
         }
     }
 
+    // Note: in order to avoid very difficult and time-consuming parallelisation
+    // of edge cell connectivity and edge cell values, we enforce a more
+    // stringent face-based consistency across processor boundaries. Basically,
+    // if a face-based consistency of 1:1 (not 2:1 as for ordinary faces) is
+    // ensured, the edge-based consistency becomes a local operation (I'm not
+    // 100% sure to be honest whether this is true all the time since there are
+    // countless variants when dealing with arbitrary polyhedral cells).
+    // See faceConsistentRefinement for details. VV, 3/Apr/2018.
+
+    // Return number of removed cells
     return nRemCells;
 }
 
@@ -3346,9 +3358,9 @@ Foam::polyhedralRefinement::polyhedralRefinement
     faceRemover_(mesh_, GREAT), // Merge boundary faces wherever possible
     maxCells_(readLabel(dict.lookup("maxCells"))),
     maxRefinementLevel_(readLabel(dict.lookup("maxRefinementLevel"))),
-    pointBasedConsistency_
+    edgeBasedConsistency_
     (
-        dict.lookupOrDefault<Switch>("pointBasedConsistency", true)
+        dict.lookupOrDefault<Switch>("edgeBasedConsistency", true)
     ),
     nRefinementBufferLayers_
     (
@@ -3426,7 +3438,7 @@ Foam::polyhedralRefinement::polyhedralRefinement
 
     // If the maximum refinementLevel is greater than 2 and the user insists on
     // not using point based refinement strategy, issue a warning
-    if (!pointBasedConsistency_ && maxRefinementLevel_ > 2)
+    if (!edgeBasedConsistency_ && maxRefinementLevel_ > 2)
     {
         WarningIn
         (
@@ -3445,7 +3457,7 @@ Foam::polyhedralRefinement::polyhedralRefinement
             << " 8:1 point conflicts."
             << nl
             << "In order to supress this message and use point based"
-            << " consistency checks, set pointBasedConsistency to true."
+            << " consistency checks, set edgeBasedConsistency to true."
             << endl;
     }
 
@@ -3583,11 +3595,11 @@ void Foam::polyhedralRefinement::setCellsToRefine
         // Reset counter at the beginning of each iteration
         nAddCells = 0;
 
-        if (pointBasedConsistency_)
+        if (edgeBasedConsistency_)
         {
-            // Check for 4:1 point based consistent refinement. Updates
+            // Check for 4:1 edge based consistent refinement. Updates
             // cellsToRefine and returns number of cells added in this iteration
-            nAddCells += pointConsistentRefinement(refineCell);
+            nAddCells += edgeConsistentRefinement(refineCell);
         }
 
         // Check for 2:1 face based consistent refinement. Updates cellsToRefine
@@ -3836,12 +3848,12 @@ void Foam::polyhedralRefinement::setSplitPointsToUnrefine
         // Reset number of removed cells from unrefinement for this iteration
         nRemCells = 0;
 
-        if (pointBasedConsistency_)
+        if (edgeBasedConsistency_)
         {
-            // Check for 4:1 point based consistent unrefinement. Updates
+            // Check for 4:1 edge based consistent unrefinement. Updates
             // cellsToUnrefine and returns number of removed cells from
             // unrefinement in this iteration
-            nRemCells += pointConsistentUnrefinement(cellsToUnrefine);
+            nRemCells += edgeConsistentUnrefinement(cellsToUnrefine);
         }
 
         // Check for 2:1 face based consistent unrefinement. Updates
@@ -4126,7 +4138,7 @@ void Foam::polyhedralRefinement::write(Ostream& os) const
         << name() << nl
         << maxCells_ << nl
         << maxRefinementLevel_ << nl
-        << pointBasedConsistency_ << nl
+        << edgeBasedConsistency_ << nl
         << nRefinementBufferLayers_ << nl
         << nUnrefinementBufferLayers_ << endl;
 }
@@ -4145,7 +4157,7 @@ void Foam::polyhedralRefinement::writeDict(Ostream& os) const
         << token::END_STATEMENT << nl
         << "    maxRefinementLevel " << maxRefinementLevel_
         << token::END_STATEMENT << nl
-        << "    pointBasedConsistency " << pointBasedConsistency_
+        << "    edgeBasedConsistency " << edgeBasedConsistency_
         << token::END_STATEMENT << nl
         << "    nRefinementBufferLayers " << nRefinementBufferLayers_
         << token::END_STATEMENT << nl
diff --git a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/polyhedralRefinement/polyhedralRefinement.H b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/polyhedralRefinement/polyhedralRefinement.H
index d4735760d..3e40eeb7d 100644
--- a/src/dynamicMesh/dynamicMesh/polyMeshModifiers/polyhedralRefinement/polyhedralRefinement.H
+++ b/src/dynamicMesh/dynamicMesh/polyMeshModifiers/polyhedralRefinement/polyhedralRefinement.H
@@ -124,8 +124,8 @@ private:
             //- Maximum number of refinement levels for a given cell
             label maxRefinementLevel_;
 
-            //- Switch whether to use point based consistency on refinement
-            Switch pointBasedConsistency_;
+            //- Switch whether to use edge based consistency on refinement
+            Switch edgeBasedConsistency_;
 
             //- Number of buffer layers for refinement
             label nRefinementBufferLayers_;
@@ -363,17 +363,17 @@ private:
             //  is obtained. Returns local number of cells changed
             label faceConsistentRefinement(boolList& cellsToRefine) const;
 
-            //- Updates cellsToRefine such that a point consistent 4:1 refinement
+            //- Updates cellsToRefine such that an edge consistent 4:1 refinement
             //  is obtained. Returns local number of cells changed
-            label pointConsistentRefinement(boolList& cellsToRefine) const;
+            label edgeConsistentRefinement(boolList& cellsToRefine) const;
 
             //- Updates cellsToUnrefine such that a face consistent 2:1
             //  unrefinement is obtained. Returns local number of cells changed
             label faceConsistentUnrefinement(boolList& cellsToUnrefine) const;
 
-            //- Updates cellsToUnrefine such that a point consistent 4:1
+            //- Updates cellsToUnrefine such that an edge consistent 4:1
             //  unrefinement is obtained. Returns local number of cells changed
-            label pointConsistentUnrefinement(boolList& cellsToUnrefine) const;
+            label edgeConsistentUnrefinement(boolList& cellsToUnrefine) const;
 
 
         // Copy control
diff --git a/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicMeshDict b/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicMeshDict
index e86bd90d2..21d60e808 100644
--- a/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicMeshDict
+++ b/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicMeshDict
@@ -27,6 +27,12 @@ dynamicPolyRefinementFvMeshCoeffs
     // Unrefine every unrefineInterval step
     unrefineInterval 1;
 
+    // Separate refinement/unrefinement steps. In case this is switched on,
+    // if both refinement and unrefinement should have been performed in a
+    // single step, unrefinement is skipped. Switched off by default, meaning
+    // that it should be safe to perform both at the same time
+    separateUpdates false;
+
     // Refinement selection criteria
     refinementSelection
     {
@@ -59,9 +65,9 @@ dynamicPolyRefinementFvMeshCoeffs
     // level inconsistencies
     nUnrefinementBufferLayers      4;
 
-    // Whether to use point based consistency check. Needed when one allows more
-    // than 2 refinement levels (automatically switched on in that case)
-    pointBasedRefinement yes;
+    // Whether to use edge based consistency check. Needed when one allows more
+    // than 2 refinement levels (automatically switched on)
+    edgeBasedConsistency           yes;
 }
 
 
diff --git a/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicPolyRefinementFvMesh.C b/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicPolyRefinementFvMesh.C
index 62ebc8d3b..12f360abd 100644
--- a/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicPolyRefinementFvMesh.C
+++ b/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicPolyRefinementFvMesh.C
@@ -59,7 +59,7 @@ void Foam::dynamicPolyRefinementFvMesh::readDict()
             << " trigerring within a certain time step and should be > 0"
             << exit(FatalError);
     }
-    
+
     unrefineInterval_ = readLabel(refinementDict_.lookup("unrefineInterval"));
     if (refineInterval_ < 1)
     {
@@ -69,6 +69,10 @@ void Foam::dynamicPolyRefinementFvMesh::readDict()
             << " trigerring within a certain time step and should be > 0"
             << exit(FatalError);
     }
+
+    // Read separate updates switch
+    separateUpdates_ =
+        refinementDict_.lookupOrDefault<Switch>("separateUpdates", false);
 }
 
 
@@ -97,6 +101,10 @@ Foam::dynamicPolyRefinementFvMesh::dynamicPolyRefinementFvMesh
     ),
     refineInterval_(readLabel(refinementDict_.lookup("refineInterval"))),
     unrefineInterval_(readLabel(refinementDict_.lookup("unrefineInterval"))),
+    separateUpdates_
+    (
+        refinementDict_.lookupOrDefault<Switch>("separateUpdates", false)
+    ),
     curTimeIndex_(-1),
 
     refinementSelectionPtr_(refinementSelection::New(*this, refinementDict_))
@@ -148,7 +156,14 @@ bool Foam::dynamicPolyRefinementFvMesh::update()
 
     // Check whether to perform refinement and/or unrefinement
     const bool performRefinement = timeID % refineInterval_ == 0;
-    const bool performUnrefinement = timeID % unrefineInterval_ == 0;
+
+    // Skip performing refinement/unrefinement in the same step if
+    // separateUpdates flag is switched on
+    bool performUnrefinement = timeID % unrefineInterval_ == 0;
+    if (performRefinement && separateUpdates_)
+    {
+        performUnrefinement = false;
+    }
 
     if
     (
diff --git a/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicPolyRefinementFvMesh.H b/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicPolyRefinementFvMesh.H
index ad6570683..25518740b 100644
--- a/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicPolyRefinementFvMesh.H
+++ b/src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicPolyRefinementFvMesh.H
@@ -68,6 +68,10 @@ class dynamicPolyRefinementFvMesh
         //- Unrefinement interval
         label unrefineInterval_;
 
+        //- Separate refinement/unrefinement: off by default, meaning that
+        //  refinement and unrefinement can be performed in the same step
+        Switch separateUpdates_;
+
         //- Current time index (helper variable to skip multiple topo changes in
         //  a single time step)
         label curTimeIndex_;

From d7a357b64697179e4daae49ee4280689e19d2c12 Mon Sep 17 00:00:00 2001
From: Vuko Vukcevic <v.vukcevic@wikki.co.uk>
Date: Tue, 24 Apr 2018 16:02:27 +0200
Subject: [PATCH 06/15] Bugfix in polyMesh::updateMesh

Need to set changing_ flag to true automatically when performin topo changes or
otherwise calling updateMesh
---
 src/foam/meshes/polyMesh/polyMeshUpdate.C | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/src/foam/meshes/polyMesh/polyMeshUpdate.C b/src/foam/meshes/polyMesh/polyMeshUpdate.C
index 5477b6ee2..b96511262 100644
--- a/src/foam/meshes/polyMesh/polyMeshUpdate.C
+++ b/src/foam/meshes/polyMesh/polyMeshUpdate.C
@@ -36,6 +36,9 @@ Description
 
 void Foam::polyMesh::updateMesh(const mapPolyMesh& mpm)
 {
+    // Mesh is changing
+    changing(true);
+
     // Update zones.  Since boundary depends on zones, they need to be
     // updated first.  HJ, 20/May/2014
     pointZones_.updateMesh();

From 8a14f3ce8ee774688666f81c3bc96ca675177c50 Mon Sep 17 00:00:00 2001
From: Vuko Vukcevic <v.vukcevic@wikki.co.uk>
Date: Tue, 24 Apr 2018 16:03:09 +0200
Subject: [PATCH 07/15] Bugfix in samplesCuttingPlane

Need to trigger the re-calculation of cutting plane on topo changes
---
 .../sampledSurface/isoSurface/isoSurface.H    |  3 +-
 .../isoSurface/sampledIsoSurface.C            |  4 +-
 .../sampledCuttingPlane/sampledCuttingPlane.C | 11 +++++
 .../sampledCuttingPlane/sampledCuttingPlane.H | 41 +++++++++++--------
 4 files changed, 38 insertions(+), 21 deletions(-)

diff --git a/src/sampling/sampledSurface/isoSurface/isoSurface.H b/src/sampling/sampledSurface/isoSurface/isoSurface.H
index fbd5c9ad7..f55c7d1c2 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurface.H
+++ b/src/sampling/sampledSurface/isoSurface/isoSurface.H
@@ -105,6 +105,7 @@ class isoSurface
         //- Reference to mesh
         const fvMesh& mesh_;
 
+        //- Point values
         const scalarField& pVals_;
 
         //- Input volScalarField with separated coupled patches rewritten
@@ -119,7 +120,6 @@ class isoSurface
         //- When to merge points
         const scalar mergeDistance_;
 
-
         //- Whether face might be cut
         List<cellCutType> faceCutType_;
 
@@ -446,7 +446,6 @@ public:
             const GeometricField<Type, fvPatchField, volMesh>& cCoords,
             const Field<Type>& pCoords
         ) const;
-
 };
 
 
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
index 28d73578e..7f889c52c 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
@@ -572,13 +572,13 @@ bool Foam::sampledIsoSurface::expire()
     facesPtr_.clear();
     subMeshPtr_.clear();
 
-    // already marked as expired
+    // Already marked as expired
     if (prevTimeIndex_ == -1)
     {
         return false;
     }
 
-    // force update
+    // Force update
     prevTimeIndex_ = -1;
     return true;
 }
diff --git a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
index 58e12ecb9..c8757fa1c 100644
--- a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
+++ b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
@@ -61,6 +61,8 @@ void Foam::sampledCuttingPlane::createGeometry()
     pointDistance_.clear();
     cellDistancePtr_.clear();
 
+    // Clear derived data
+    clearGeom();
 
     // Get any subMesh
     if (zoneID_.index() != -1 && !subMeshPtr_.valid())
@@ -311,6 +313,12 @@ Foam::sampledCuttingPlane::~sampledCuttingPlane()
 
 bool Foam::sampledCuttingPlane::needsUpdate() const
 {
+    // Update for changing mesh
+    if (mesh().changing())
+    {
+        needsUpdate_ = true;
+    }
+
     return needsUpdate_;
 }
 
@@ -327,6 +335,9 @@ bool Foam::sampledCuttingPlane::expire()
     // Clear any stored topologies
     facesPtr_.clear();
 
+    // Clear derived data
+    clearGeom();
+
     // already marked as expired
     if (needsUpdate_)
     {
diff --git a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H
index db0e340fc..52a21b9f1 100644
--- a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H
+++ b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H
@@ -68,10 +68,10 @@ class sampledCuttingPlane
         //- Whether to recalculate cell values as average of point values
         const Switch average_;
 
-        //- zone name/index (if restricted to zones)
+        //- Zone name/index (if restricted to zones)
         mutable cellZoneID zoneID_;
 
-        //- for zones: patch to put exposed faces into
+        //- For zones: patch to put exposed faces into
         mutable word exposedPatchName_;
 
         //- Track if the surface needs an update
@@ -90,7 +90,7 @@ class sampledCuttingPlane
         //- Constructed iso surface
         autoPtr<isoSurface> isoSurfPtr_;
 
-        //- triangles converted to faceList
+        //- Triangles converted to faceList
         mutable autoPtr<faceList> facesPtr_;
 
 
@@ -99,7 +99,7 @@ class sampledCuttingPlane
         //- Create iso surface
         void createGeometry();
 
-        //- sample field on faces
+        //- Sample field on faces
         template <class Type>
         tmp<Field<Type> > sampleField
         (
@@ -129,9 +129,8 @@ public:
         );
 
 
-    // Destructor
-
-        virtual ~sampledCuttingPlane();
+    //- Destructor
+    virtual ~sampledCuttingPlane();
 
 
     // Member Functions
@@ -177,67 +176,75 @@ public:
             return isoSurfPtr_();
         }
 
-        //- sample field on surface
+
+    // Sample
+
+        //- Sample field on surface
         virtual tmp<scalarField> sample
         (
             const volScalarField&
         ) const;
 
-        //- sample field on surface
+        //- Sample field on surface
         virtual tmp<vectorField> sample
         (
             const volVectorField&
         ) const;
 
-        //- sample field on surface
+        //- Sample field on surface
         virtual tmp<sphericalTensorField> sample
         (
             const volSphericalTensorField&
         ) const;
 
-        //- sample field on surface
+        //- Sample field on surface
         virtual tmp<symmTensorField> sample
         (
             const volSymmTensorField&
         ) const;
 
-        //- sample field on surface
+        //- Sample field on surface
         virtual tmp<tensorField> sample
         (
             const volTensorField&
         ) const;
 
 
-        //- interpolate field on surface
+    // Interpolate
+
+        //- Interpolate field on surface
         virtual tmp<scalarField> interpolate
         (
             const interpolation<scalar>&
         ) const;
 
-        //- interpolate field on surface
+        //- Interpolate field on surface
         virtual tmp<vectorField> interpolate
         (
             const interpolation<vector>&
         ) const;
 
-        //- interpolate field on surface
+        //- Interpolate field on surface
         virtual tmp<sphericalTensorField> interpolate
         (
             const interpolation<sphericalTensor>&
         ) const;
 
-        //- interpolate field on surface
+        //- Interpolate field on surface
         virtual tmp<symmTensorField> interpolate
         (
             const interpolation<symmTensor>&
         ) const;
 
-        //- interpolate field on surface
+        //- Interpolate field on surface
         virtual tmp<tensorField> interpolate
         (
             const interpolation<tensor>&
         ) const;
 
+
+    // Output
+
         //- Write
         virtual void print(Ostream&) const;
 };

From 047bc41b2730065b239cc2e6ad4d59935465158a Mon Sep 17 00:00:00 2001
From: Vuko Vukcevic <v.vukcevic@wikki.co.uk>
Date: Wed, 25 Apr 2018 10:00:27 +0200
Subject: [PATCH 08/15] Updates to polyMesh::changing_ flag

It's better to update it in polyTopoChanger::changeMesh() than in
polyMesh::updateMesh(). This way, it is consistent with directTopoChange engine
---
 .../polyTopoChanger/polyTopoChangerChangeMesh.C             | 6 ++++++
 src/foam/meshes/polyMesh/polyMeshUpdate.C                   | 3 ---
 2 files changed, 6 insertions(+), 3 deletions(-)

diff --git a/src/dynamicMesh/dynamicMesh/polyTopoChange/polyTopoChanger/polyTopoChangerChangeMesh.C b/src/dynamicMesh/dynamicMesh/polyTopoChange/polyTopoChanger/polyTopoChangerChangeMesh.C
index 2a7111c0e..933478830 100644
--- a/src/dynamicMesh/dynamicMesh/polyTopoChange/polyTopoChanger/polyTopoChangerChangeMesh.C
+++ b/src/dynamicMesh/dynamicMesh/polyTopoChange/polyTopoChanger/polyTopoChangerChangeMesh.C
@@ -2491,6 +2491,9 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChanger::changeMesh()
         // Increment the morph index
         morphIndex_++;
 
+        // Mark the mesh as changing
+        mesh_.changing(true);
+        
         return topoChangeMap;
     }
     else
@@ -2504,6 +2507,9 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChanger::changeMesh()
 
             // Sync mesh update
             mesh_.syncUpdateMesh();
+
+            // Mark the mesh as changing
+            mesh_.changing(true);
         }
 
         return autoPtr<mapPolyMesh>(new mapPolyMesh(mesh_));
diff --git a/src/foam/meshes/polyMesh/polyMeshUpdate.C b/src/foam/meshes/polyMesh/polyMeshUpdate.C
index b96511262..5477b6ee2 100644
--- a/src/foam/meshes/polyMesh/polyMeshUpdate.C
+++ b/src/foam/meshes/polyMesh/polyMeshUpdate.C
@@ -36,9 +36,6 @@ Description
 
 void Foam::polyMesh::updateMesh(const mapPolyMesh& mpm)
 {
-    // Mesh is changing
-    changing(true);
-
     // Update zones.  Since boundary depends on zones, they need to be
     // updated first.  HJ, 20/May/2014
     pointZones_.updateMesh();

From 6c22d372bdc252b80acdf1ba2973ae14eeb17d40 Mon Sep 17 00:00:00 2001
From: Vuko Vukcevic <v.vukcevic@wikki.co.uk>
Date: Wed, 25 Apr 2018 10:14:52 +0200
Subject: [PATCH 09/15] Updates to donorAcceptor's default constructor

Initializing values acceptor and donor data to invalid (-1) values
---
 src/overset/oversetMesh/donorAcceptor/donorAcceptor.H | 9 +++++++++
 1 file changed, 9 insertions(+)

diff --git a/src/overset/oversetMesh/donorAcceptor/donorAcceptor.H b/src/overset/oversetMesh/donorAcceptor/donorAcceptor.H
index 37692f350..d5f0a69b9 100644
--- a/src/overset/oversetMesh/donorAcceptor/donorAcceptor.H
+++ b/src/overset/oversetMesh/donorAcceptor/donorAcceptor.H
@@ -122,6 +122,15 @@ public:
 
         //- Construct null
         inline donorAcceptor()
+        :
+            acceptorCell_(-1),
+            acceptorProcNo_(-1),
+            acceptorPoint_(vector::zero),
+            donorCell_(-1),
+            donorProcNo_(-1),
+            donorPoint_(vector::zero),
+            extendedDonorCells_(),
+            extendedDonorPoints_()
         {}
 
         //- Construct from acceptor data

From 40530a7921a4c8523c1da9eb83f765334c29a39c Mon Sep 17 00:00:00 2001
From: Vuko Vukcevic <v.vukcevic@wikki.co.uk>
Date: Wed, 25 Apr 2018 12:37:23 +0200
Subject: [PATCH 10/15] Conservative bridged overlap GGI verification case

Channel test case (7 cells) with sudden contraction using bridged overlap with
potentialFoam.
---
 .../basicMesh-potentialFoam/0/U               |  49 ++++++++
 .../basicMesh-potentialFoam/0/p               |  40 +++++++
 .../basicMesh-potentialFoam/Allclean          |  11 ++
 .../basicMesh-potentialFoam/Allrun            |  13 +++
 .../basicMesh-potentialFoam/boundary          |  61 ++++++++++
 .../constant/polyMesh/blockMeshDict           | 110 ++++++++++++++++++
 .../constant/transportProperties              |  22 ++++
 .../basicMesh-potentialFoam/save/U            |  49 ++++++++
 .../basicMesh-potentialFoam/save/p            |  40 +++++++
 .../basicMesh-potentialFoam/setBatch.batch    |   3 +
 .../system/controlDict                        |  67 +++++++++++
 .../basicMesh-potentialFoam/system/fvSchemes  |  49 ++++++++
 .../basicMesh-potentialFoam/system/fvSolution |  52 +++++++++
 .../verification/GGI/bridgedOverlap/note      |   4 +
 14 files changed, 570 insertions(+)
 create mode 100644 validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/0/U
 create mode 100644 validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/0/p
 create mode 100755 validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/Allclean
 create mode 100755 validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/Allrun
 create mode 100644 validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/boundary
 create mode 100644 validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/constant/polyMesh/blockMeshDict
 create mode 100644 validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/constant/transportProperties
 create mode 100644 validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/save/U
 create mode 100644 validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/save/p
 create mode 100644 validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/setBatch.batch
 create mode 100644 validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/system/controlDict
 create mode 100644 validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/system/fvSchemes
 create mode 100644 validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/system/fvSolution
 create mode 100644 validationAndVerificationSuite/verification/GGI/bridgedOverlap/note

diff --git a/validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/0/U b/validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/0/U
new file mode 100644
index 000000000..0ba3887d1
--- /dev/null
+++ b/validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/0/U
@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | foam-extend: Open Source CFD                    |
+|  \\    /   O peration     | Version:     4.0                                |
+|   \\  /    A nd           | Web:         http://www.foam-extend.org         |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    location    "0";
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 1);
+
+boundaryField
+{
+    walls 
+    {
+        type            fixedValue;
+        value           uniform (0 0 1);
+    }
+    inlet
+    {
+        type            fixedValue;
+        value           uniform (0 0 1);
+    }
+    outlet
+    {
+        type            zeroGradient;
+    }
+    movingSlider
+    {
+        type            ggi;
+    }
+    staticSlider
+    {
+        type            ggi;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/0/p b/validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/0/p
new file mode 100644
index 000000000..93e85eb4e
--- /dev/null
+++ b/validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/0/p
@@ -0,0 +1,40 @@
+// -*- C++ -*-
+// File generated by PyFoam - sorry for the ugliness
+
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class volScalarField;
+ object p;
+}
+
+dimensions [ 0 2 -2 0 0 0 0 ];
+
+internalField uniform 0;
+
+boundaryField
+{
+  walls 
+  {
+      type  zeroGradient;
+  }
+  staticSlider
+  {
+    type ggi;
+  }
+  movingSlider
+  {
+    type ggi;
+  }
+  inlet
+  {
+      type  zeroGradient;
+  }
+  outlet
+  {
+    type fixedValue;
+      value uniform 0;
+  }
+} 	// ************************************************************************* //
+
diff --git a/validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/Allclean b/validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/Allclean
new file mode 100755
index 000000000..e8ea64fef
--- /dev/null
+++ b/validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/Allclean
@@ -0,0 +1,11 @@
+#!/bin/sh
+
+# Source tutorial clean functions
+. $WM_PROJECT_DIR/bin/tools/CleanFunctions
+
+cleanCase
+\rm -rf VTK
+\rm -rf 0
+\cp -r save 0
+\rm -f constant/polyMesh/boundary
+\rm -rf constant/polyMesh/sets
diff --git a/validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/Allrun b/validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/Allrun
new file mode 100755
index 000000000..96107a4e0
--- /dev/null
+++ b/validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/Allrun
@@ -0,0 +1,13 @@
+#!/bin/sh
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+application="potentialFoam"
+
+runApplication blockMesh
+\cp -r boundary constant/polyMesh/.
+
+runApplication setSet -batch setBatch.batch
+runApplication setsToZones -noFlipMap
+
+runApplication $application
diff --git a/validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/boundary b/validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/boundary
new file mode 100644
index 000000000..67f6acd06
--- /dev/null
+++ b/validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/boundary
@@ -0,0 +1,61 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | foam-extend: Open Source CFD                    |
+|  \\    /   O peration     | Version:     4.1                                |
+|   \\  /    A nd           | Web:         http://www.foam-extend.org         |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       polyBoundaryMesh;
+    location    "constant/polyMesh";
+    object      boundary;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+5
+(
+    walls
+    {
+        type            wall;
+        nFaces          24;
+        startFace       8;
+    }
+    inlet
+    {
+        type            patch;
+        nFaces          1;
+        startFace       32;
+    }
+    outlet
+    {
+        type            patch;
+        nFaces          3;
+        startFace       33;
+    }
+    movingSlider
+    {
+        type            ggi;
+        nFaces          3;
+        startFace       36;
+
+        shadowPatch     staticSlider;
+        zone            movingSliderZone;
+        bridgeOverlap   true;
+    }
+    staticSlider
+    {
+        type            ggi;
+        nFaces          1;
+        startFace       39;
+        
+        shadowPatch     movingSlider;
+        zone            staticSliderZone;
+        bridgeOverlap   true;
+    }
+)
+
+
+// ************************************************************************* //
diff --git a/validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/constant/polyMesh/blockMeshDict b/validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/constant/polyMesh/blockMeshDict
new file mode 100644
index 000000000..17e95b2e9
--- /dev/null
+++ b/validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/constant/polyMesh/blockMeshDict
@@ -0,0 +1,110 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | foam-extend: Open Source CFD                    |
+|  \\    /   O peration     | Version:     3.2                                |
+|   \\  /    A nd           | Web:         http://www.foam-extend.org         |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 1;
+
+vertices
+(
+    (0 0 0)
+    (1 0 0)
+    (1 1 0)
+    (0 1 0)
+    (0 0 5)
+    (1 0 5)
+    (1 1 5)
+    (0 1 5)
+    
+    (0 0 5)
+    (1 0 5)
+    (1 0.56 5)
+    (0 0.56 5)
+    (0 0 10)
+    (1 0 10)
+    (1 0.56 10)
+    (0 0.56 10)
+);
+
+blocks
+(
+    hex (0 1 2 3 4 5 6 7) stationary (1 1 2) simpleGrading (1 1 1)
+    hex (8 9 10 11 12 13 14 15) moving (1 3 2) simpleGrading (1 1 1)   
+);
+
+edges
+(
+    
+);
+
+boundary
+(
+    walls 
+    {
+        type wall;
+        faces
+        (            
+            (0 1 5 4)
+            (1 2 6 5)
+            (2 3 7 6)
+            (0 4 7 3)
+            (9 10 14 13)
+            (10 11 15 14)
+            (8 12 15 11)
+            (8 9 13 12)
+        );
+    }
+    
+    inlet
+    {
+        type patch;
+        faces
+        (
+            (0 3 2 1)            
+        );
+    }
+    
+    outlet
+    {
+        type patch;
+        faces
+        (
+            (12 13 14 15)            
+        );
+    }
+
+    movingSlider
+    {
+        type patch;
+        faces
+        (
+            (9 8 11 10)
+        );
+    }
+
+    staticSlider
+    {
+        type patch;
+        faces
+        (
+            (4 5 6 7)
+        );
+    }
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //
diff --git a/validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/constant/transportProperties b/validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/constant/transportProperties
new file mode 100644
index 000000000..02ce1d3df
--- /dev/null
+++ b/validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/constant/transportProperties
@@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | foam-extend: Open Source CFD                    |
+|  \\    /   O peration     | Version:     4.0                                |
+|   \\  /    A nd           | Web:         http://www.foam-extend.org         |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+transportModel  Newtonian;
+
+//nu              nu [0 2 -1 0 0 0 0] 1e-5;
+nu              nu [0 2 -1 0 0 0 0] 0;
+
+// ************************************************************************* //
diff --git a/validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/save/U b/validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/save/U
new file mode 100644
index 000000000..0ba3887d1
--- /dev/null
+++ b/validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/save/U
@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | foam-extend: Open Source CFD                    |
+|  \\    /   O peration     | Version:     4.0                                |
+|   \\  /    A nd           | Web:         http://www.foam-extend.org         |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    location    "0";
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 1);
+
+boundaryField
+{
+    walls 
+    {
+        type            fixedValue;
+        value           uniform (0 0 1);
+    }
+    inlet
+    {
+        type            fixedValue;
+        value           uniform (0 0 1);
+    }
+    outlet
+    {
+        type            zeroGradient;
+    }
+    movingSlider
+    {
+        type            ggi;
+    }
+    staticSlider
+    {
+        type            ggi;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/save/p b/validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/save/p
new file mode 100644
index 000000000..93e85eb4e
--- /dev/null
+++ b/validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/save/p
@@ -0,0 +1,40 @@
+// -*- C++ -*-
+// File generated by PyFoam - sorry for the ugliness
+
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class volScalarField;
+ object p;
+}
+
+dimensions [ 0 2 -2 0 0 0 0 ];
+
+internalField uniform 0;
+
+boundaryField
+{
+  walls 
+  {
+      type  zeroGradient;
+  }
+  staticSlider
+  {
+    type ggi;
+  }
+  movingSlider
+  {
+    type ggi;
+  }
+  inlet
+  {
+      type  zeroGradient;
+  }
+  outlet
+  {
+    type fixedValue;
+      value uniform 0;
+  }
+} 	// ************************************************************************* //
+
diff --git a/validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/setBatch.batch b/validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/setBatch.batch
new file mode 100644
index 000000000..6161b21e7
--- /dev/null
+++ b/validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/setBatch.batch
@@ -0,0 +1,3 @@
+faceSet staticSliderZone new patchToFace staticSlider
+faceSet movingSliderZone new patchToFace movingSlider
+quit
diff --git a/validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/system/controlDict b/validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/system/controlDict
new file mode 100644
index 000000000..640b64b63
--- /dev/null
+++ b/validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/system/controlDict
@@ -0,0 +1,67 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | foam-extend: Open Source CFD                    |
+|  \\    /   O peration     | Version:     3.1                                |
+|   \\  /    A nd           | Web:         http://www.extend-project.de       |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     icoFoam;
+
+startFrom       latestTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         1;
+
+deltaT          0.001;
+
+writeControl    timeStep;
+
+writeInterval   100;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  10;
+
+writeCompression compressed;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable yes;
+
+adjustTimeStep  no;
+
+maxCo           0.2;
+
+functions
+(
+    ggiCheck
+    {
+        type ggiCheck;
+        phi  phi;
+        functionObjectLibs ("libcheckFunctionObjects.so");
+    }
+);
+
+DebugSwitches
+{
+    ggi 1;
+    GGIInterpolation 1;
+}
+
+// ************************************************************************* //
diff --git a/validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/system/fvSchemes b/validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/system/fvSchemes
new file mode 100644
index 000000000..34ea0e6ee
--- /dev/null
+++ b/validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/system/fvSchemes
@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | foam-extend: Open Source CFD                    |
+|  \\    /   O peration     | Version:     3.1                                |
+|   \\  /    A nd           | Web:         http://www.extend-project.de       |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+}
+
+divSchemes
+{
+    div(phi,U)      Gauss upwind;
+
+    div((muEff*dev2(T(grad(U)))))  Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear limited 0.5;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         limited 0.5;
+}
+
+// ************************************************************************* //
diff --git a/validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/system/fvSolution b/validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/system/fvSolution
new file mode 100644
index 000000000..3d849e4cb
--- /dev/null
+++ b/validationAndVerificationSuite/verification/GGI/bridgedOverlap/basicMesh-potentialFoam/system/fvSolution
@@ -0,0 +1,52 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | foam-extend: Open Source CFD                    |
+|  \\    /   O peration     | Version:     3.1                                |
+|   \\  /    A nd           | Web:         http://www.extend-project.de       |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    p
+    {
+        solver           BiCGStab;
+        preconditioner   Cholesky;
+
+        minIter          0;
+        maxIter          1000;
+        tolerance        1e-10;
+        relTol           0.0;
+    }
+    U
+    {
+        solver           BiCGStab;
+        preconditioner   DILU;
+
+        minIter          0;
+        maxIter          1000;
+        tolerance        1e-07;
+        relTol           0;
+    }
+}
+
+SIMPLE
+{
+    nNonOrthogonalCorrectors 0;
+}
+
+PISO
+{
+    nCorrectors              2;
+    nNonOrthogonalCorrectors 0;
+}
+
+// ************************************************************************* //
diff --git a/validationAndVerificationSuite/verification/GGI/bridgedOverlap/note b/validationAndVerificationSuite/verification/GGI/bridgedOverlap/note
new file mode 100644
index 000000000..ba5adcf54
--- /dev/null
+++ b/validationAndVerificationSuite/verification/GGI/bridgedOverlap/note
@@ -0,0 +1,4 @@
+The conservative overlapping GGI needs certain top level changes for solvers
+other than potentialFoam. Basically, all parts that make up the flux on the RHS
+of pressure equation need to be interpolated together. This is not the case with
+the latest version of time and under-relaxation consistent solvers.

From ef69c12106d5bce13e2e9aa2e52387e11a1e1cf7 Mon Sep 17 00:00:00 2001
From: Vuko Vukcevic <v.vukcevic@wikki.co.uk>
Date: Wed, 25 Apr 2018 13:29:14 +0200
Subject: [PATCH 11/15] Graded normal extrude model

---
 src/mesh/extrudeModel/Make/files              |   1 +
 src/mesh/extrudeModel/Make/options            |   6 +-
 .../extrudeModel/gradedNormal/gradedNormal.C  | 143 ++++++++++++++++
 .../extrudeModel/gradedNormal/gradedNormal.H  | 158 ++++++++++++++++++
 4 files changed, 306 insertions(+), 2 deletions(-)
 create mode 100644 src/mesh/extrudeModel/gradedNormal/gradedNormal.C
 create mode 100644 src/mesh/extrudeModel/gradedNormal/gradedNormal.H

diff --git a/src/mesh/extrudeModel/Make/files b/src/mesh/extrudeModel/Make/files
index c8d0cc571..8941cb004 100644
--- a/src/mesh/extrudeModel/Make/files
+++ b/src/mesh/extrudeModel/Make/files
@@ -1,6 +1,7 @@
 extrudeModel/extrudeModel.C
 extrudeModel/newExtrudeModel.C
 linearNormal/linearNormal.C
+gradedNormal/gradedNormal.C
 linearRadial/linearRadial.C
 sigmaRadial/sigmaRadial.C
 wedge/wedge.C
diff --git a/src/mesh/extrudeModel/Make/options b/src/mesh/extrudeModel/Make/options
index c0517d317..d402ac311 100644
--- a/src/mesh/extrudeModel/Make/options
+++ b/src/mesh/extrudeModel/Make/options
@@ -1,8 +1,10 @@
 EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/dynamicMesh/lnInclude
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/ODE/lnInclude
 
 LIB_LIBS = \
     -lmeshTools \
-    -ldynamicMesh
+    -ldynamicMesh \
+    -lODE
 
diff --git a/src/mesh/extrudeModel/gradedNormal/gradedNormal.C b/src/mesh/extrudeModel/gradedNormal/gradedNormal.C
new file mode 100644
index 000000000..d2cb5c152
--- /dev/null
+++ b/src/mesh/extrudeModel/gradedNormal/gradedNormal.C
@@ -0,0 +1,143 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation, either version 3 of the License, or (at your
+    option) any later version.
+
+    foam-extend is distributed in the hope that it will be useful, but
+    WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+    General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "gradedNormal.H"
+#include "addToRunTimeSelectionTable.H"
+#include "BisectionRoot.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace extrudeModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(gradedNormal, 0);
+
+addToRunTimeSelectionTable(extrudeModel, gradedNormal, dictionary);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+gradedNormal::gradedNormal(const dictionary& dict)
+:
+    extrudeModel(typeName, dict),
+    thickness_(readScalar(coeffDict_.lookup("thickness"))),
+    delta0_(readScalar(coeffDict_.lookup("initialCellLength"))),
+    expansionFactor_(1.0)
+{
+    // Sanity checks
+    if (thickness_ <= SMALL)
+    {
+        FatalErrorIn("gradedNormal(const dictionary&)")
+            << "thickness should be positive: " << thickness_
+            << exit(FatalError);
+    }
+
+    if (delta0_ <= SMALL)
+    {
+        FatalErrorIn("gradedNormal(const dictionary&)")
+            << "initialCellLength should be positive: " << delta0_
+            << exit(FatalError);
+    }
+
+    const scalar maxExpFactor =
+        coeffDict_.lookupOrDefault<scalar>("maxExpansionFactor", 3.0);
+
+    if (maxExpFactor <= SMALL)
+    {
+        FatalErrorIn("gradedNormal(const dictionary&)")
+            << "maxExpansionFactor should be positive: " << maxExpFactor
+            << exit(FatalError);
+    }
+
+    const scalar bisectionTol =
+        coeffDict_.lookupOrDefault<scalar>("bisectionTol", 1e-5);
+
+    if (bisectionTol <= SMALL)
+    {
+        FatalErrorIn("gradedNormal(const dictionary&)")
+            << "bisectionTolerance should be positive: " << bisectionTol
+            << exit(FatalError);
+    }
+
+    // Create expansion factor equation represented as a function object
+    expansionFactorEqn eqn(*this);
+
+    // Calculate the expansionFactor using the bisection algorithm with the
+    // default tolerance of 1e-5
+    BisectionRoot<expansionFactorEqn> rootFinder
+    (
+        eqn,
+        bisectionTol
+    );
+
+    // Search for the root in [0, 3], where upper bound 3 is default
+    expansionFactor_ = rootFinder.root
+    (
+        0.0,
+        maxExpFactor
+    );
+
+    // Report the result
+    Info<< "Calculated expansion factor: " << expansionFactor_ << endl;
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+gradedNormal::~gradedNormal()
+{}
+
+
+// * * * * * * * * * * * * * * * * Operators * * * * * * * * * * * * * * * * //
+
+point gradedNormal::operator()
+(
+    const point& surfacePoint,
+    const vector& surfaceNormal,
+    const label layer
+) const
+{
+    scalar d = 0.0;
+
+    for (label i = 0; i < layer; ++i)
+    {
+        d += delta0_*pow(expansionFactor_, i + 1);
+    }
+
+    return surfacePoint + d*surfaceNormal;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace extrudeModels
+} // End namespace Foam
+
+// ************************************************************************* //
+
diff --git a/src/mesh/extrudeModel/gradedNormal/gradedNormal.H b/src/mesh/extrudeModel/gradedNormal/gradedNormal.H
new file mode 100644
index 000000000..bd2127bf0
--- /dev/null
+++ b/src/mesh/extrudeModel/gradedNormal/gradedNormal.H
@@ -0,0 +1,158 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation, either version 3 of the License, or (at your
+    option) any later version.
+
+    foam-extend is distributed in the hope that it will be useful, but
+    WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+    General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::extrudeModels::gradedNormal
+
+Description
+    Extrudes by transforming points normal to the surface. Input parameters are:
+    1. Extrusion thickness (in meters),
+    2. Initial delta (in meters),
+    3. Number of layers.
+
+    Expansion factor is calculated numerically using bisection algorithm.
+
+Author
+    Vuko Vukcevic. FMENA Zagreb.  All rights reserved.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef gradedNormal_H
+#define gradedNormal_H
+
+#include "point.H"
+#include "extrudeModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace extrudeModels
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class gradedNormal Declaration
+\*---------------------------------------------------------------------------*/
+
+class gradedNormal
+:
+    public extrudeModel
+{
+    // Private data
+
+        //- Layer thickness
+        scalar thickness_;
+
+        //- Initial delta (cell size at the extruded patch)
+        scalar delta0_;
+
+        //- Expansion factor
+        scalar expansionFactor_;
+
+
+    // Expansion factor equation (functor class used by BisectionRoot)
+    class expansionFactorEqn
+    {
+        // Private data
+
+            //- Const reference to underlying gradedNormal model
+            const gradedNormal& gnm_;
+
+
+    public:
+
+        //- Construct given gradedNormal model
+        explicit expansionFactorEqn(const gradedNormal& gnm)
+        :
+            gnm_(gnm)
+        {}
+
+
+        //- Function call operator
+        scalar operator()(const scalar& x) const
+        {
+            scalar result = gnm_.thickness();
+
+            for (label i = 0; i <= gnm_.nLayers(); ++i)
+            {
+                result -= gnm_.delta0()*pow(x, i);
+            }
+
+            return result;
+        }
+    };
+
+
+public:
+
+    //- Runtime type information
+    TypeName("gradedNormal");
+
+    // Constructors
+
+        //- Construct from components
+        gradedNormal(const dictionary& dict);
+
+
+    //- Destructor
+    ~gradedNormal();
+
+
+    // Access functions
+
+        //- Return const reference to thickness
+        const scalar& thickness() const
+        {
+            return thickness_;
+        }
+
+
+        //- Return const reference to initial delta
+        const scalar& delta0() const
+        {
+            return delta0_;
+        }
+
+
+    // Member Operators
+
+        point operator()
+        (
+            const point& surfacePoint,
+            const vector& surfaceNormal,
+            const label layer
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace extrudeModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
+

From b0811b398d361431018dabcf1103eb32f27ea4a2 Mon Sep 17 00:00:00 2001
From: Vuko Vukcevic <v.vukcevic@wikki.co.uk>
Date: Wed, 2 May 2018 15:30:17 +0200
Subject: [PATCH 12/15] Bugfix in argList.C

Correct initialization and override of global controlDict switches
---
 src/foam/global/argList/argList.C | 36 ++++++++++++++++++++++++++++++-
 1 file changed, 35 insertions(+), 1 deletion(-)

diff --git a/src/foam/global/argList/argList.C b/src/foam/global/argList/argList.C
index ab384ebc9..98ec14a3a 100644
--- a/src/foam/global/argList/argList.C
+++ b/src/foam/global/argList/argList.C
@@ -780,6 +780,40 @@ void Foam::argList::parse
     }
 
 
+    // Managing the overrides for the global control switches:
+    //
+    // Here is the order of precedence for the definition/overriding of the
+    // control switches, from lowest to highest:
+    //  - source code definitions from the various libraries/solvers
+    //  - file specified by the env. variable FOAM_GLOBAL_CONTROLDICT
+    //  - case's system/controlDict file
+    //  - command-line parameters
+    //
+    // First, we allow the users to specify the location of a centralized
+    // global controlDict dictionary using the environment variable
+    // FOAM_GLOBAL_CONTROLDICT.
+    fileName optionalGlobControlDictFileName =
+    getEnv("FOAM_GLOBAL_CONTROLDICT");
+
+    if (optionalGlobControlDictFileName.size() )
+    {
+        debug::updateCentralDictVars
+        (
+            optionalGlobControlDictFileName,
+            Pstream::master() && bannerEnabled
+        );
+    }
+
+    // Now that the rootPath_/globalCase_ directory is known (following the
+    // call to getRootCase()), we grab any global control switches overrides
+    // from the current case's controlDict.
+
+    debug::updateCentralDictVars
+    (
+        rootPath_/globalCase_/"system/controlDict",
+        Pstream::master() && bannerEnabled
+    );
+
     // Finally, a command-line override for central controlDict's variables.
     // This is the ultimate override for the global control switches.
 
@@ -793,7 +827,7 @@ void Foam::argList::parse
 
     forAll (globalControlDictSwitchSetNames, gI)
     {
-        word switchSetName = globalControlDictSwitchSetNames.names[gI];
+        const word switchSetName = globalControlDictSwitchSetNames.names[gI];
 
         if (optionFound(switchSetName))
         {

From 7ad4da9e8927ab7031bb955567b2aa13369e14c3 Mon Sep 17 00:00:00 2001
From: Vuko Vukcevic <v.vukcevic@wikki.co.uk>
Date: Thu, 3 May 2018 10:14:35 +0200
Subject: [PATCH 13/15] Update to vtkSurfaceWriter

Writing "double" instead of "float" more reliable for most paraview versions.

Author: Karl Meredith
---
 .../sampledSurface/writers/vtk/vtkSurfaceWriter.C  | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/src/sampling/sampledSurface/writers/vtk/vtkSurfaceWriter.C b/src/sampling/sampledSurface/writers/vtk/vtkSurfaceWriter.C
index 0a6b4d6cd..a36fcce2e 100644
--- a/src/sampling/sampledSurface/writers/vtk/vtkSurfaceWriter.C
+++ b/src/sampling/sampledSurface/writers/vtk/vtkSurfaceWriter.C
@@ -55,7 +55,7 @@ void Foam::vtkSurfaceWriter::writeGeometry
         << "DATASET POLYDATA" << nl;
 
     // Write vertex coords
-    os  << "POINTS " << points.size() << " float" << nl;
+    os  << "POINTS " << points.size() << " double" << nl;
     forAll(points, pointI)
     {
         const point& pt = points[pointI];
@@ -100,7 +100,7 @@ namespace Foam
         const Field<scalar>& values
     )
     {
-        os  << "1 " << values.size() << " float" << nl;
+        os  << "1 " << values.size() << " double" << nl;
 
         forAll(values, elemI)
         {
@@ -129,7 +129,7 @@ namespace Foam
         const Field<vector>& values
     )
     {
-        os  << "3 " << values.size() << " float" << nl;
+        os  << "3 " << values.size() << " double" << nl;
 
         forAll(values, elemI)
         {
@@ -147,7 +147,7 @@ namespace Foam
         const Field<sphericalTensor>& values
     )
     {
-        os  << "1 " << values.size() << " float" << nl;
+        os  << "1 " << values.size() << " double" << nl;
 
         forAll(values, elemI)
         {
@@ -164,7 +164,7 @@ namespace Foam
         const Field<symmTensor>& values
     )
     {
-        os  << "6 " << values.size() << " float" << nl;
+        os  << "6 " << values.size() << " double" << nl;
 
         forAll(values, elemI)
         {
@@ -185,7 +185,7 @@ namespace Foam
         const Field<tensor>& values
     )
     {
-        os  << "9 " << values.size() << " float" << nl;
+        os  << "9 " << values.size() << " double" << nl;
 
         forAll(values, elemI)
         {
@@ -210,7 +210,7 @@ void Foam::vtkSurfaceWriter::writeData
     const Field<Type>& values
 )
 {
-    os  << "1 " << values.size() << " float" << nl;
+    os  << "1 " << values.size() << " double" << nl;
 
     forAll(values, elemI)
     {

From 4d7859d8d27cfcbbbb2b0a4c455bf4683662e217 Mon Sep 17 00:00:00 2001
From: Vuko Vukcevic <v.vukcevic@wikki.co.uk>
Date: Wed, 9 May 2018 14:28:48 +0200
Subject: [PATCH 14/15] Merged thermalModels

There were two sets of thermalModels in conjugateHeatTransfer library and
solidModels library. They are now merged into the conjugateHeatTransferLibrary
---
 .../elasticThermalSolidFoam/Make/options      |   6 +-
 .../constantThermal/constantThermal.C         |  57 +++-
 .../constantThermal/constantThermal.H         |  12 +
 .../multiMaterialThermal.C                    |  75 +++++
 .../multiMaterialThermal.H                    |   9 +
 .../multiMaterialZonesThermal.C               |   4 +-
 .../multiMaterialZonesThermal.H               |   1 +
 .../thermalLaws/thermalLaw/newThermalLaw.C    |   4 +-
 .../thermalLaws/thermalLaw/thermalLaw.H       |   7 +
 .../thermalModel/thermalModel.C               |  10 +-
 .../thermalModel/thermalModel.H               |  13 +
 src/solidModels/Make/files                    |   7 -
 src/solidModels/Make/options                  |   2 +
 .../constantThermal/constantThermal.C         | 169 ----------
 .../constantThermal/constantThermal.H         | 131 --------
 .../multiMaterialThermal.C                    | 293 ------------------
 .../multiMaterialThermal.H                    | 125 --------
 .../thermalLaws/thermalLaw/newThermalLaw.C    |  79 -----
 .../thermalLaws/thermalLaw/thermalLaw.C       |  62 ----
 .../thermalLaws/thermalLaw/thermalLaw.H       | 170 ----------
 src/solidModels/thermalModel/thermalModel.C   |  90 ------
 src/solidModels/thermalModel/thermalModel.H   | 148 ---------
 22 files changed, 190 insertions(+), 1284 deletions(-)
 delete mode 100644 src/solidModels/thermalModel/thermalLaws/constantThermal/constantThermal.C
 delete mode 100644 src/solidModels/thermalModel/thermalLaws/constantThermal/constantThermal.H
 delete mode 100644 src/solidModels/thermalModel/thermalLaws/multiMaterialThermal/multiMaterialThermal.C
 delete mode 100644 src/solidModels/thermalModel/thermalLaws/multiMaterialThermal/multiMaterialThermal.H
 delete mode 100644 src/solidModels/thermalModel/thermalLaws/thermalLaw/newThermalLaw.C
 delete mode 100644 src/solidModels/thermalModel/thermalLaws/thermalLaw/thermalLaw.C
 delete mode 100644 src/solidModels/thermalModel/thermalLaws/thermalLaw/thermalLaw.H
 delete mode 100644 src/solidModels/thermalModel/thermalModel.C
 delete mode 100644 src/solidModels/thermalModel/thermalModel.H

diff --git a/applications/solvers/solidMechanics/elasticThermalSolidFoam/Make/options b/applications/solvers/solidMechanics/elasticThermalSolidFoam/Make/options
index 0aa531afe..c1dbfdb14 100644
--- a/applications/solvers/solidMechanics/elasticThermalSolidFoam/Make/options
+++ b/applications/solvers/solidMechanics/elasticThermalSolidFoam/Make/options
@@ -5,7 +5,8 @@ EXE_INC = \
     -I$(LIB_SRC)/dynamicMesh/topoChangerFvMesh/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/solidModels/lnInclude
+    -I$(LIB_SRC)/solidModels/lnInclude \
+    -I$(LIB_SRC)/conjugateHeatTransfer/lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \
@@ -13,4 +14,5 @@ EXE_LIBS = \
     -ldynamicFvMesh \
     -ldynamicMesh \
     -ltopoChangerFvMesh \
-    -lsolidModels
+    -lsolidModels \
+    -lconjugateHeatTransfer
diff --git a/src/conjugateHeatTransfer/thermalModel/thermalLaws/constantThermal/constantThermal.C b/src/conjugateHeatTransfer/thermalModel/thermalLaws/constantThermal/constantThermal.C
index f2833dbc7..346bd5e04 100644
--- a/src/conjugateHeatTransfer/thermalModel/thermalLaws/constantThermal/constantThermal.C
+++ b/src/conjugateHeatTransfer/thermalModel/thermalLaws/constantThermal/constantThermal.C
@@ -49,7 +49,9 @@ Foam::constantThermal::constantThermal
     thermalLaw(name, T, dict),
     rho_(dict.lookup("rho")),
     C_(dict.lookup("C")),
-    k_(dict.lookup("k"))
+    k_(dict.lookup("k")),
+    alpha_(dict.lookup("alpha")),
+    T0_(dict.lookup("T0"))
 {}
 
 
@@ -86,6 +88,7 @@ Foam::tmp<Foam::volScalarField> Foam::constantThermal::rho() const
     return tresult;
 }
 
+
 Foam::tmp<Foam::volScalarField> Foam::constantThermal::C() const
 {
     tmp<volScalarField> tresult
@@ -138,4 +141,56 @@ Foam::tmp<Foam::volScalarField> Foam::constantThermal::k() const
 }
 
 
+Foam::tmp<Foam::volScalarField> Foam::constantThermal::alpha() const
+{
+    tmp<volScalarField> tresult
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "alpha",
+                mesh().time().timeName(),
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            alpha_,
+            zeroGradientFvPatchScalarField::typeName
+        )
+    );
+
+    tresult().correctBoundaryConditions();
+
+    return tresult;
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::constantThermal::T0() const
+{
+    tmp<volScalarField> tresult
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "T0",
+                mesh().time().timeName(),
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            T0_,
+            zeroGradientFvPatchScalarField::typeName
+        )
+    );
+
+    tresult().correctBoundaryConditions();
+
+    return tresult;
+}
+
+
 // ************************************************************************* //
diff --git a/src/conjugateHeatTransfer/thermalModel/thermalLaws/constantThermal/constantThermal.H b/src/conjugateHeatTransfer/thermalModel/thermalLaws/constantThermal/constantThermal.H
index 0a9aa5b8f..c603f1bc7 100644
--- a/src/conjugateHeatTransfer/thermalModel/thermalLaws/constantThermal/constantThermal.H
+++ b/src/conjugateHeatTransfer/thermalModel/thermalLaws/constantThermal/constantThermal.H
@@ -64,6 +64,12 @@ class constantThermal
         //- Thermal conductivity
         dimensionedScalar k_;
 
+        //- Thermal expansion coefficient
+        dimensionedScalar alpha_;
+
+        //- Reference temperature
+        dimensionedScalar T0_;
+
 
     // Private Member Functions
 
@@ -109,6 +115,12 @@ public:
         //- Return thermal conductivity
         virtual tmp<volScalarField> k() const;
 
+        //- Return thermal expansion coefficient
+        virtual tmp<volScalarField> alpha() const;
+
+        //- Return reference temperature
+        virtual tmp<volScalarField> T0() const;
+
         virtual void correct()
         {}
 };
diff --git a/src/conjugateHeatTransfer/thermalModel/thermalLaws/multiMaterialThermal/multiMaterialThermal.C b/src/conjugateHeatTransfer/thermalModel/thermalLaws/multiMaterialThermal/multiMaterialThermal.C
index 84172e1c6..bb31c7288 100644
--- a/src/conjugateHeatTransfer/thermalModel/thermalLaws/multiMaterialThermal/multiMaterialThermal.C
+++ b/src/conjugateHeatTransfer/thermalModel/thermalLaws/multiMaterialThermal/multiMaterialThermal.C
@@ -109,6 +109,7 @@ void Foam::multiMaterialThermal::readLaws
     }
 }
 
+
 void Foam::multiMaterialThermal::checkLaws() const
 {
     const PtrList<thermalLaw>& laws = *this;
@@ -173,6 +174,7 @@ Foam::multiMaterialThermal::multiMaterialThermal
     checkLaws();
 }
 
+
 // Construct from dictionary and create default material field
 Foam::multiMaterialThermal::multiMaterialThermal
 (
@@ -244,6 +246,7 @@ Foam::tmp<Foam::volScalarField> Foam::multiMaterialThermal::rho() const
     return tresult;
 }
 
+
 Foam::tmp<Foam::volScalarField> Foam::multiMaterialThermal::C() const
 {
     tmp<volScalarField> tresult
@@ -310,6 +313,78 @@ Foam::tmp<Foam::volScalarField> Foam::multiMaterialThermal::k() const
 }
 
 
+Foam::tmp<Foam::volScalarField> Foam::multiMaterialThermal::alpha() const
+{
+    tmp<volScalarField> tresult
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "alpha",
+                mesh().time().timeName(),
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            dimensionedScalar("zeroE", dimless/dimTemperature, 0),
+            zeroGradientFvPatchScalarField::typeName
+        )
+    );
+    volScalarField& result = tresult();
+
+    // Accumulate data for all fields
+    const PtrList<thermalLaw>& laws = *this;
+
+    forAll (laws, lawI)
+    {
+        result.internalField() +=
+            indicator(lawI)*laws[lawI].alpha()().internalField();
+    }
+
+    result.correctBoundaryConditions();
+
+    return tresult;
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::multiMaterialThermal::T0() const
+{
+    tmp<volScalarField> tresult
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "T0",
+                mesh().time().timeName(),
+                mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh(),
+            dimensionedScalar("zeroT0", dimTemperature, 0),
+            zeroGradientFvPatchScalarField::typeName
+        )
+    );
+    volScalarField& result = tresult();
+
+    // Accumulate data for all fields
+    const PtrList<thermalLaw>& laws = *this;
+
+    forAll (laws, lawI)
+    {
+        result.internalField() +=
+            indicator(lawI)*laws[lawI].T0()().internalField();
+    }
+
+    result.correctBoundaryConditions();
+
+    return tresult;
+}
+
+
 void Foam::multiMaterialThermal::correct()
 {
     PtrList<thermalLaw>& laws = *this;
diff --git a/src/conjugateHeatTransfer/thermalModel/thermalLaws/multiMaterialThermal/multiMaterialThermal.H b/src/conjugateHeatTransfer/thermalModel/thermalLaws/multiMaterialThermal/multiMaterialThermal.H
index 67a2a4b9e..cdea16dae 100644
--- a/src/conjugateHeatTransfer/thermalModel/thermalLaws/multiMaterialThermal/multiMaterialThermal.H
+++ b/src/conjugateHeatTransfer/thermalModel/thermalLaws/multiMaterialThermal/multiMaterialThermal.H
@@ -66,6 +66,7 @@ class multiMaterialThermal
         //- Calculate indicator field given index
         tmp<volScalarField> indicator(const label index) const;
 
+
 protected:
 
     // Protected data
@@ -73,6 +74,7 @@ protected:
         //- Material indicator field
         volScalarField materials_;
 
+
     // Protected member functions
 
         //- Check that the laws and the material field are okay
@@ -113,6 +115,7 @@ public:
             const scalar defaultMaterial
         );
 
+
     // Destructor
 
         virtual ~multiMaterialThermal();
@@ -129,6 +132,12 @@ public:
         //- Return thermal conductivity
         virtual tmp<volScalarField> k() const;
 
+        //- Return thermal expansion coefficient
+        virtual tmp<volScalarField> alpha() const;
+
+        //- Return reference temperature
+        virtual tmp<volScalarField> T0() const;
+
         //- Correct the rheological model
         virtual void correct();
 };
diff --git a/src/conjugateHeatTransfer/thermalModel/thermalLaws/multiMaterialZonesThermal/multiMaterialZonesThermal.C b/src/conjugateHeatTransfer/thermalModel/thermalLaws/multiMaterialZonesThermal/multiMaterialZonesThermal.C
index c1f3607cb..936aa8094 100644
--- a/src/conjugateHeatTransfer/thermalModel/thermalLaws/multiMaterialZonesThermal/multiMaterialZonesThermal.C
+++ b/src/conjugateHeatTransfer/thermalModel/thermalLaws/multiMaterialZonesThermal/multiMaterialZonesThermal.C
@@ -54,13 +54,13 @@ Foam::multiMaterialZonesThermal::multiMaterialZonesThermal
     PtrList<thermalLaw>& laws = *this;
     forAll (laws, lawI)
     {
-        wordList zones (lawEntries[lawI].dict().lookup("zones"));
+        const wordList zones(lawEntries[lawI].dict().lookup("zones"));
 
         forAll(zones, zoneI)
         {
             const label zoneID = mesh().cellZones().findZoneID(zones[zoneI]);
 
-            if ( zoneID < 0 )
+            if (zoneID < 0)
             {
                 FatalErrorIn
                 (
diff --git a/src/conjugateHeatTransfer/thermalModel/thermalLaws/multiMaterialZonesThermal/multiMaterialZonesThermal.H b/src/conjugateHeatTransfer/thermalModel/thermalLaws/multiMaterialZonesThermal/multiMaterialZonesThermal.H
index 7db5fc697..1a8b95691 100644
--- a/src/conjugateHeatTransfer/thermalModel/thermalLaws/multiMaterialZonesThermal/multiMaterialZonesThermal.H
+++ b/src/conjugateHeatTransfer/thermalModel/thermalLaws/multiMaterialZonesThermal/multiMaterialZonesThermal.H
@@ -61,6 +61,7 @@ class multiMaterialZonesThermal
         //- Disallow default bitwise assignment
         void operator=(const multiMaterialZonesThermal&);
 
+
 public:
 
     //- Runtime type information
diff --git a/src/conjugateHeatTransfer/thermalModel/thermalLaws/thermalLaw/newThermalLaw.C b/src/conjugateHeatTransfer/thermalModel/thermalLaws/thermalLaw/newThermalLaw.C
index a68810a4d..d1755e342 100644
--- a/src/conjugateHeatTransfer/thermalModel/thermalLaws/thermalLaw/newThermalLaw.C
+++ b/src/conjugateHeatTransfer/thermalModel/thermalLaws/thermalLaw/newThermalLaw.C
@@ -44,7 +44,7 @@ autoPtr<thermalLaw> thermalLaw::New
     const dictionary& dict
 )
 {
-    word rheoTypeName = dict.lookup("type");
+    const word rheoTypeName = dict.lookup("type");
 
     Info<< "Selecting thermal model " << rheoTypeName << endl;
 
@@ -64,7 +64,7 @@ autoPtr<thermalLaw> thermalLaw::New
         )   << "Unknown thermalLaw type "
             << rheoTypeName << endl << endl
             << "Valid  thermalLaws are : " << endl
-            << dictionaryConstructorTablePtr_->toc()
+            << dictionaryConstructorTablePtr_->sortedToc()
             << exit(FatalIOError);
     }
 
diff --git a/src/conjugateHeatTransfer/thermalModel/thermalLaws/thermalLaw/thermalLaw.H b/src/conjugateHeatTransfer/thermalModel/thermalLaws/thermalLaw/thermalLaw.H
index 3f2c1f60d..0207cb0c5 100644
--- a/src/conjugateHeatTransfer/thermalModel/thermalLaws/thermalLaw/thermalLaw.H
+++ b/src/conjugateHeatTransfer/thermalModel/thermalLaws/thermalLaw/thermalLaw.H
@@ -83,6 +83,7 @@ protected:
             return T_.mesh();
         }
 
+
 public:
 
     //- Runtime type information
@@ -156,6 +157,12 @@ public:
         //- Return thermal conductivity
         virtual tmp<volScalarField> k() const = 0;
 
+        //- Return thermal expansion coefficient
+        virtual tmp<volScalarField> alpha() const = 0;
+
+        //- Return reference temperature
+        virtual tmp<volScalarField> T0() const = 0;
+
         //- Correct the thermal model
         virtual void correct() = 0;
 };
diff --git a/src/conjugateHeatTransfer/thermalModel/thermalModel.C b/src/conjugateHeatTransfer/thermalModel/thermalModel.C
index e8c866685..1a50c0358 100644
--- a/src/conjugateHeatTransfer/thermalModel/thermalModel.C
+++ b/src/conjugateHeatTransfer/thermalModel/thermalModel.C
@@ -118,9 +118,10 @@ void thermalModel::modifyResistance
     }
 }
 
+
 tmp<fvScalarMatrix> thermalModel::laplacian(volScalarField& T)
 {
-    word kScheme ("laplacian(k,T)");
+    const word kScheme("laplacian(k,T)");
 
     surfaceScalarField kf = fvc::interpolate(lawPtr_->k());
 
@@ -128,10 +129,11 @@ tmp<fvScalarMatrix> thermalModel::laplacian(volScalarField& T)
 
     return tmp<fvScalarMatrix>
     (
-        new fvScalarMatrix( fvm::laplacian(kf, T, kScheme) )
+        new fvScalarMatrix(fvm::laplacian(kf, T, kScheme))
     );
 }
 
+
 tmp<volScalarField> thermalModel::S() const
 {
     tmp<volScalarField> tsource
@@ -155,15 +157,17 @@ tmp<volScalarField> thermalModel::S() const
             )
         )
     );
+    volScalarField& source = tsource();
 
     forAll(sourcePtr_, sourceI)
     {
-        sourcePtr_[sourceI].addSource(tsource());
+        sourcePtr_[sourceI].addSource(source);
     }
 
     return tsource;
 }
 
+
 bool thermalModel::read()
 {
     if (regIOobject::read())
diff --git a/src/conjugateHeatTransfer/thermalModel/thermalModel.H b/src/conjugateHeatTransfer/thermalModel/thermalModel.H
index b248b62bf..48368ea48 100644
--- a/src/conjugateHeatTransfer/thermalModel/thermalModel.H
+++ b/src/conjugateHeatTransfer/thermalModel/thermalModel.H
@@ -76,6 +76,7 @@ class thermalModel
         //- Thermal Source
         PtrList<thermalSource> sourcePtr_;
 
+
     // Private Member Functions
 
         //- Disallow copy construct
@@ -128,6 +129,18 @@ public:
             return lawPtr_->k();
         }
 
+        //- Return thermal expansion coefficient
+        tmp<volScalarField> alpha() const
+        {
+            return lawPtr_->alpha();
+        }
+
+        //- Return reference temperature
+        tmp<volScalarField> T0() const
+        {
+            return lawPtr_->T0();
+        }
+
         //- Add resistance for gaps
         void modifyResistance(surfaceScalarField&) const;
 
diff --git a/src/solidModels/Make/files b/src/solidModels/Make/files
index c0b293a29..7bd7e1242 100644
--- a/src/solidModels/Make/files
+++ b/src/solidModels/Make/files
@@ -79,12 +79,5 @@ $(rheologyLaws)/multiMaterial/multiMaterial.C
 $(rheologyLaws)/orthotropicLinearElastic/orthotropicLinearElastic.C
 $(rheologyLaws)/PronyViscoelastic/PronyViscoelastic.C
 
-thermalModel/thermalModel.C
-thermalLaws = thermalModel/thermalLaws
-$(thermalLaws)/thermalLaw/thermalLaw.C
-$(thermalLaws)/thermalLaw/newThermalLaw.C
-$(thermalLaws)/constantThermal/constantThermal.C
-$(thermalLaws)/multiMaterialThermal/multiMaterialThermal.C
-
 
 LIB = $(FOAM_LIBBIN)/libsolidModels
diff --git a/src/solidModels/Make/options b/src/solidModels/Make/options
index 95481d075..86d7896f8 100644
--- a/src/solidModels/Make/options
+++ b/src/solidModels/Make/options
@@ -3,6 +3,7 @@ EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/finiteArea/lnInclude \
     -I$(LIB_SRC)/lagrangian/lnInclude \
+    -I$(LIB_SRC)/conjugateHeatTransfer/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \
@@ -12,6 +13,7 @@ LIB_LIBS = \
     -lfiniteVolume \
     -lmeshTools \
     -lfiniteArea \
+    -lconjugateHeatTransfer \
     -ldynamicFvMesh \
     -ldynamicMesh \
     -ltopoChangerFvMesh \
diff --git a/src/solidModels/thermalModel/thermalLaws/constantThermal/constantThermal.C b/src/solidModels/thermalModel/thermalLaws/constantThermal/constantThermal.C
deleted file mode 100644
index aa421d8a8..000000000
--- a/src/solidModels/thermalModel/thermalLaws/constantThermal/constantThermal.C
+++ /dev/null
@@ -1,169 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | foam-extend: Open Source CFD
-   \\    /   O peration     | Version:     4.0
-    \\  /    A nd           | Web:         http://www.foam-extend.org
-     \\/     M anipulation  | For copyright notice see file Copyright
--------------------------------------------------------------------------------
-License
-    This file is part of foam-extend.
-
-    foam-extend is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation, either version 3 of the License, or (at your
-    option) any later version.
-
-    foam-extend is distributed in the hope that it will be useful, but
-    WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-    General Public License for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "constantThermal.H"
-#include "addToRunTimeSelectionTable.H"
-#include "zeroGradientFvPatchFields.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    defineTypeNameAndDebug(constantThermal, 0);
-    addToRunTimeSelectionTable(thermalLaw, constantThermal, dictionary);
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-// Construct from dictionary
-Foam::constantThermal::constantThermal
-(
-    const word& name,
-    const volScalarField& T,
-    const dictionary& dict
-)
-:
-    thermalLaw(name, T, dict),
-    C_(dict.lookup("C")),
-    k_(dict.lookup("k")),
-    alpha_(dict.lookup("alpha")),
-    T0_(dict.lookup("T0"))
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::constantThermal::~constantThermal()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-Foam::tmp<Foam::volScalarField> Foam::constantThermal::C() const
-{
-    tmp<volScalarField> tresult
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "C",
-                mesh().time().timeName(),
-                mesh(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh(),
-            C_,
-            zeroGradientFvPatchScalarField::typeName
-        )
-    );
-
-    tresult().correctBoundaryConditions();
-
-    return tresult;
-}
-
-
-Foam::tmp<Foam::volScalarField> Foam::constantThermal::k() const
-{
-    tmp<volScalarField> tresult
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "k",
-                mesh().time().timeName(),
-                mesh(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh(),
-            k_,
-            zeroGradientFvPatchScalarField::typeName
-        )
-    );
-
-    tresult().correctBoundaryConditions();
-
-    return tresult;
-}
-
-
-Foam::tmp<Foam::volScalarField> Foam::constantThermal::alpha() const
-{
-    tmp<volScalarField> tresult
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "alpha",
-                mesh().time().timeName(),
-                mesh(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh(),
-            alpha_,
-            zeroGradientFvPatchScalarField::typeName
-        )
-    );
-
-    tresult().correctBoundaryConditions();
-
-    return tresult;
-}
-
-
-Foam::tmp<Foam::volScalarField> Foam::constantThermal::T0() const
-{
-    tmp<volScalarField> tresult
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "T0",
-                mesh().time().timeName(),
-                mesh(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh(),
-            T0_,
-            zeroGradientFvPatchScalarField::typeName
-        )
-    );
-
-    tresult().correctBoundaryConditions();
-
-    return tresult;
-}
-
-
-// ************************************************************************* //
diff --git a/src/solidModels/thermalModel/thermalLaws/constantThermal/constantThermal.H b/src/solidModels/thermalModel/thermalLaws/constantThermal/constantThermal.H
deleted file mode 100644
index 12e833239..000000000
--- a/src/solidModels/thermalModel/thermalLaws/constantThermal/constantThermal.H
+++ /dev/null
@@ -1,131 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | foam-extend: Open Source CFD
-   \\    /   O peration     | Version:     4.0
-    \\  /    A nd           | Web:         http://www.foam-extend.org
-     \\/     M anipulation  | For copyright notice see file Copyright
--------------------------------------------------------------------------------
-License
-    This file is part of foam-extend.
-
-    foam-extend is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation, either version 3 of the License, or (at your
-    option) any later version.
-
-    foam-extend is distributed in the hope that it will be useful, but
-    WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-    General Public License for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    constantThermal
-
-Description
-    Constant thermal properties
-
-Author
-    Hrvoje Jasak, Wikki Ltd.  All rights reserved.
-
-SourceFiles
-    constantThermal.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef constantThermal_H
-#define constantThermal_H
-
-#include "thermalLaw.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                         Class constantThermal Declaration
-\*---------------------------------------------------------------------------*/
-
-class constantThermal
-:
-    public thermalLaw
-{
-    // Private data
-
-        //- Specific heat capacity
-        dimensionedScalar C_;
-
-        //- Thermal conductivity
-        dimensionedScalar k_;
-
-        //- Thermal expansion coefficient
-        dimensionedScalar alpha_;
-
-        //- Reference temperature
-        dimensionedScalar T0_;
-
-
-    // Private Member Functions
-
-        //- Disallow default bitwise copy construct
-        constantThermal(const constantThermal&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const constantThermal&);
-
-
-public:
-
-    //- Runtime type information
-    TypeName("constant");
-
-    // Static data members
-
-
-    // Constructors
-
-        //- Construct from dictionary
-        constantThermal
-        (
-            const word& name,
-            const volScalarField& T,
-            const dictionary& dict
-        );
-
-
-    // Destructor
-
-        virtual ~constantThermal();
-
-
-    // Member Functions
-
-        //- Return specific heat capacity
-        virtual tmp<volScalarField> C() const;
-
-        //- Return thermal conductivity
-        virtual tmp<volScalarField> k() const;
-
-        //- Return thermal expansion coefficient
-        virtual tmp<volScalarField> alpha() const;
-
-        //- Return reference temperature
-        virtual tmp<volScalarField> T0() const;
-
-        virtual void correct()
-        {}
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/solidModels/thermalModel/thermalLaws/multiMaterialThermal/multiMaterialThermal.C b/src/solidModels/thermalModel/thermalLaws/multiMaterialThermal/multiMaterialThermal.C
deleted file mode 100644
index 39ed73d99..000000000
--- a/src/solidModels/thermalModel/thermalLaws/multiMaterialThermal/multiMaterialThermal.C
+++ /dev/null
@@ -1,293 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | foam-extend: Open Source CFD
-   \\    /   O peration     | Version:     4.0
-    \\  /    A nd           | Web:         http://www.foam-extend.org
-     \\/     M anipulation  | For copyright notice see file Copyright
--------------------------------------------------------------------------------
-License
-    This file is part of foam-extend.
-
-    foam-extend is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation, either version 3 of the License, or (at your
-    option) any later version.
-
-    foam-extend is distributed in the hope that it will be useful, but
-    WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-    General Public License for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "multiMaterialThermal.H"
-#include "addToRunTimeSelectionTable.H"
-#include "zeroGradientFvPatchFields.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    defineTypeNameAndDebug(multiMaterialThermal, 0);
-    addToRunTimeSelectionTable(thermalLaw, multiMaterialThermal, dictionary);
-}
-
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-Foam::tmp<Foam::scalarField> Foam::multiMaterialThermal::indicator
-(
-    const label i
-) const
-{
-    const scalarField& mat = materials_.internalField();
-
-    tmp<scalarField> tresult(new scalarField(mat.size(), 0.0));
-    scalarField& result = tresult();
-
-    forAll (mat, matI)
-    {
-        if (mat[matI] > i - SMALL && mat[matI] < i + 1 - SMALL)
-        {
-            result[matI] = 1.0;
-        }
-    }
-
-    return tresult;
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-// Construct from dictionary
-Foam::multiMaterialThermal::multiMaterialThermal
-(
-    const word& name,
-    const volScalarField& T,
-    const dictionary& dict
-)
-:
-    thermalLaw(name, T, dict),
-    PtrList<thermalLaw>(),
-    materials_
-    (
-        IOobject
-        (
-            "materials",
-            mesh().time().timeName(),
-            mesh(),
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh()
-    )
-{
-    PtrList<thermalLaw>& laws = *this;
-
-    PtrList<entry> lawEntries(dict.lookup("laws"));
-    laws.setSize(lawEntries.size());
-
-    forAll (laws, lawI)
-    {
-        laws.set
-        (
-            lawI,
-            thermalLaw::New
-            (
-                lawEntries[lawI].keyword(),
-                T,
-                lawEntries[lawI].dict()
-            )
-        );
-    }
-
-    if
-    (
-        min(materials_).value() < 0
-     || max(materials_).value() > laws.size() + SMALL
-    )
-    {
-        FatalErrorIn
-        (
-            "multiMaterialThermal::multiMaterialThermal\n"
-            "(\n"
-            "    const word& name,\n"
-            "    const volScalarField& T,\n"
-            "    const dictionary& dict\n"
-            ")"
-        )   << "Invalid definition of material indicator field.  "
-            << "Number of materials: " << laws.size()
-            << " max index: " << max(materials_)
-            << ".  Should be " << laws.size() - 1
-            << abort(FatalError);
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::multiMaterialThermal::~multiMaterialThermal()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-Foam::tmp<Foam::volScalarField> Foam::multiMaterialThermal::C() const
-{
-    tmp<volScalarField> tresult
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "C",
-                mesh().time().timeName(),
-                mesh(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh(),
-            dimensionedScalar("zeroC", dimSpecificHeatCapacity, 0),
-            zeroGradientFvPatchScalarField::typeName
-        )
-    );
-    volScalarField& result = tresult();
-
-    // Accumulate data for all fields
-    const PtrList<thermalLaw>& laws = *this;
-
-    forAll (laws, lawI)
-    {
-        result.internalField() +=
-            indicator(lawI)*laws[lawI].C()().internalField();
-    }
-
-    result.correctBoundaryConditions();
-
-    return tresult;
-}
-
-
-Foam::tmp<Foam::volScalarField> Foam::multiMaterialThermal::k() const
-{
-    tmp<volScalarField> tresult
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "k",
-                mesh().time().timeName(),
-                mesh(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh(),
-            dimensionedScalar("zerok", dimThermalConductivity, 0),
-            zeroGradientFvPatchScalarField::typeName
-        )
-    );
-    volScalarField& result = tresult();
-
-    // Accumulate data for all fields
-    const PtrList<thermalLaw>& laws = *this;
-
-    forAll (laws, lawI)
-    {
-        result.internalField() +=
-            indicator(lawI)*laws[lawI].k()().internalField();
-    }
-
-    result.correctBoundaryConditions();
-
-    return tresult;
-}
-
-
-Foam::tmp<Foam::volScalarField> Foam::multiMaterialThermal::alpha() const
-{
-    tmp<volScalarField> tresult
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "alpha",
-                mesh().time().timeName(),
-                mesh(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh(),
-            dimensionedScalar("zeroE", dimless/dimTemperature, 0),
-            zeroGradientFvPatchScalarField::typeName
-        )
-    );
-    volScalarField& result = tresult();
-
-    // Accumulate data for all fields
-    const PtrList<thermalLaw>& laws = *this;
-
-    forAll (laws, lawI)
-    {
-        result.internalField() +=
-            indicator(lawI)*laws[lawI].alpha()().internalField();
-    }
-
-    result.correctBoundaryConditions();
-
-    return tresult;
-}
-
-
-Foam::tmp<Foam::volScalarField> Foam::multiMaterialThermal::T0() const
-{
-    tmp<volScalarField> tresult
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "T0",
-                mesh().time().timeName(),
-                mesh(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh(),
-            dimensionedScalar("zeroT0", dimTemperature, 0),
-            zeroGradientFvPatchScalarField::typeName
-        )
-    );
-    volScalarField& result = tresult();
-
-    // Accumulate data for all fields
-    const PtrList<thermalLaw>& laws = *this;
-
-    forAll (laws, lawI)
-    {
-        result.internalField() +=
-            indicator(lawI)*laws[lawI].T0()().internalField();
-    }
-
-    result.correctBoundaryConditions();
-
-    return tresult;
-}
-
-
-void Foam::multiMaterialThermal::correct()
-{
-    PtrList<thermalLaw>& laws = *this;
-
-    forAll (laws, lawI)
-    {
-        laws[lawI].correct();
-    }
-}
-
-
-// ************************************************************************* //
diff --git a/src/solidModels/thermalModel/thermalLaws/multiMaterialThermal/multiMaterialThermal.H b/src/solidModels/thermalModel/thermalLaws/multiMaterialThermal/multiMaterialThermal.H
deleted file mode 100644
index 44e0ef02a..000000000
--- a/src/solidModels/thermalModel/thermalLaws/multiMaterialThermal/multiMaterialThermal.H
+++ /dev/null
@@ -1,125 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | foam-extend: Open Source CFD
-   \\    /   O peration     | Version:     4.0
-    \\  /    A nd           | Web:         http://www.foam-extend.org
-     \\/     M anipulation  | For copyright notice see file Copyright
--------------------------------------------------------------------------------
-License
-    This file is part of foam-extend.
-
-    foam-extend is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation, either version 3 of the License, or (at your
-    option) any later version.
-
-    foam-extend is distributed in the hope that it will be useful, but
-    WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-    General Public License for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    multiMaterialThermal
-
-Description
-    Zoned multi-material thermal model controlled by an indicator field.
-
-Author
-    Hrvoje Jasak, Wikki Ltd.  All rights reserved.
-
-SourceFiles
-    multiMaterialThermal.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef multiMaterialThermal_H
-#define multiMaterialThermal_H
-
-#include "thermalLaw.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                    Class multiMaterialThermal Declaration
-\*---------------------------------------------------------------------------*/
-
-class multiMaterialThermal
-:
-    public thermalLaw,
-    public PtrList<thermalLaw>
-{
-    // Private data
-
-        //- Material indicator field
-        volScalarField materials_;
-
-    // Private Member Functions
-
-        //- Disallow default bitwise copy construct
-        multiMaterialThermal(const multiMaterialThermal&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const multiMaterialThermal&);
-
-
-        //- Calculate indicator field given index
-        tmp<scalarField> indicator(const label index) const;
-
-public:
-
-    //- Runtime type information
-    TypeName("multiMaterial");
-
-    // Static data members
-
-
-    // Constructors
-
-        //- Construct from dictionary
-        multiMaterialThermal
-        (
-            const word& name,
-            const volScalarField& T,
-            const dictionary& dict
-        );
-
-
-    // Destructor
-
-        virtual ~multiMaterialThermal();
-
-
-    // Member Functions
-
-        //- Return specific heat capacity
-        virtual tmp<volScalarField> C() const;
-
-        //- Return thermal conductivity
-        virtual tmp<volScalarField> k() const;
-
-        //- Return thermal expansion coefficient
-        virtual tmp<volScalarField> alpha() const;
-
-        //- Return reference temperature
-        virtual tmp<volScalarField> T0() const;
-
-        //- Correct the rheological model
-        virtual void correct();
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/solidModels/thermalModel/thermalLaws/thermalLaw/newThermalLaw.C b/src/solidModels/thermalModel/thermalLaws/thermalLaw/newThermalLaw.C
deleted file mode 100644
index baa5ad0fc..000000000
--- a/src/solidModels/thermalModel/thermalLaws/thermalLaw/newThermalLaw.C
+++ /dev/null
@@ -1,79 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | foam-extend: Open Source CFD
-   \\    /   O peration     | Version:     4.0
-    \\  /    A nd           | Web:         http://www.foam-extend.org
-     \\/     M anipulation  | For copyright notice see file Copyright
--------------------------------------------------------------------------------
-License
-    This file is part of foam-extend.
-
-    foam-extend is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation, either version 3 of the License, or (at your
-    option) any later version.
-
-    foam-extend is distributed in the hope that it will be useful, but
-    WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-    General Public License for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    thermalLaw
-
-\*---------------------------------------------------------------------------*/
-
-#include "thermalLaw.H"
-#include "volFields.H"
-#include "surfaceFields.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-autoPtr<thermalLaw> thermalLaw::New
-(
-    const word& name,
-    const volScalarField& T,
-    const dictionary& dict
-)
-{
-    word rheoTypeName = dict.lookup("type");
-
-    Info<< "Selecting thermal model " << rheoTypeName << endl;
-
-    dictionaryConstructorTable::iterator cstrIter =
-        dictionaryConstructorTablePtr_->find(rheoTypeName);
-
-    if (cstrIter == dictionaryConstructorTablePtr_->end())
-    {
-        FatalIOErrorIn
-        (
-            "thermalLaw::New(\n"
-            "    const word& name,\n"
-            "    const volScalarField& T,\n"
-            "    const dictionary& dict\n"
-            ")",
-            dict
-        )   << "Unknown thermalLaw type "
-            << rheoTypeName << endl << endl
-            << "Valid  thermalLaws are : " << endl
-            << dictionaryConstructorTablePtr_->sortedToc()
-            << exit(FatalIOError);
-    }
-
-    return autoPtr<thermalLaw>(cstrIter()(name, T, dict));
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// ************************************************************************* //
diff --git a/src/solidModels/thermalModel/thermalLaws/thermalLaw/thermalLaw.C b/src/solidModels/thermalModel/thermalLaws/thermalLaw/thermalLaw.C
deleted file mode 100644
index b04f5c2d9..000000000
--- a/src/solidModels/thermalModel/thermalLaws/thermalLaw/thermalLaw.C
+++ /dev/null
@@ -1,62 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | foam-extend: Open Source CFD
-   \\    /   O peration     | Version:     4.0
-    \\  /    A nd           | Web:         http://www.foam-extend.org
-     \\/     M anipulation  | For copyright notice see file Copyright
--------------------------------------------------------------------------------
-License
-    This file is part of foam-extend.
-
-    foam-extend is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation, either version 3 of the License, or (at your
-    option) any later version.
-
-    foam-extend is distributed in the hope that it will be useful, but
-    WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-    General Public License for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    thermalLaw
-
-\*---------------------------------------------------------------------------*/
-
-#include "thermalLaw.H"
-#include "volFields.H"
-#include "fvc.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-defineTypeNameAndDebug(thermalLaw, 0);
-defineRunTimeSelectionTable(thermalLaw, dictionary);
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-thermalLaw::thermalLaw
-(
-    const word& name,
-    const volScalarField& T,
-    const dictionary& dict
-)
-:
-    name_(name),
-    T_(T)
-{}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// ************************************************************************* //
diff --git a/src/solidModels/thermalModel/thermalLaws/thermalLaw/thermalLaw.H b/src/solidModels/thermalModel/thermalLaws/thermalLaw/thermalLaw.H
deleted file mode 100644
index adc96a615..000000000
--- a/src/solidModels/thermalModel/thermalLaws/thermalLaw/thermalLaw.H
+++ /dev/null
@@ -1,170 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | foam-extend: Open Source CFD
-   \\    /   O peration     | Version:     4.0
-    \\  /    A nd           | Web:         http://www.foam-extend.org
-     \\/     M anipulation  | For copyright notice see file Copyright
--------------------------------------------------------------------------------
-License
-    This file is part of foam-extend.
-
-    foam-extend is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation, either version 3 of the License, or (at your
-    option) any later version.
-
-    foam-extend is distributed in the hope that it will be useful, but
-    WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-    General Public License for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    thermalLaw
-
-Description
-    Thermal material properties for solids.
-
-Author
-    Hrvoje Jasak, Wikki Ltd.  All rights reserved.
-
-SourceFiles
-    thermalLaw.C
-    newThermalLaw.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef thermalLaw_H
-#define thermalLaw_H
-
-#include "IOdictionary.H"
-#include "typeInfo.H"
-#include "runTimeSelectionTables.H"
-#include "volFields.H"
-#include "tmp.H"
-#include "autoPtr.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                         Class thermalLaw Declaration
-\*---------------------------------------------------------------------------*/
-
-class thermalLaw
-{
-    // Private data
-
-        //- Name
-        const word name_;
-
-        //- Reference to stress field
-        const volScalarField& T_;
-
-
-    // Private Member Functions
-
-        //- Disallow copy construct
-        thermalLaw(const thermalLaw&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const thermalLaw&);
-
-
-protected:
-
-        //- Return reference to mesh
-        const fvMesh& mesh() const
-        {
-            return T_.mesh();
-        }
-
-public:
-
-    //- Runtime type information
-    TypeName("thermalLaw");
-
-
-    // Declare run-time constructor selection table
-
-        declareRunTimeSelectionTable
-        (
-            autoPtr,
-            thermalLaw,
-            dictionary,
-            (
-                const word name,
-                const volScalarField& T,
-                const dictionary& dict
-            ),
-            (name, T, dict)
-        );
-
-
-    // Selectors
-
-        //- Return a reference to the selected thermal model
-        static autoPtr<thermalLaw> New
-        (
-            const word& name,
-            const volScalarField& T,
-            const dictionary& dict
-        );
-
-
-    // Constructors
-
-        //- Construct from dictionary
-        thermalLaw
-        (
-            const word& name,
-            const volScalarField& T,
-            const dictionary& dict
-        );
-
-
-    // Destructor
-
-        virtual ~thermalLaw()
-        {}
-
-
-    // Member Functions
-
-        //- Return name
-        const word& name() const
-        {
-            return name_;
-        }
-
-
-        //- Return specific heat capacity
-        virtual tmp<volScalarField> C() const = 0;
-
-        //- Return thermal conductivity
-        virtual tmp<volScalarField> k() const = 0;
-
-        //- Return thermal expansion coefficient
-        virtual tmp<volScalarField> alpha() const = 0;
-
-        //- Return reference temperature
-        virtual tmp<volScalarField> T0() const = 0;
-
-        //- Correct the rheological model
-        virtual void correct() = 0;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/solidModels/thermalModel/thermalModel.C b/src/solidModels/thermalModel/thermalModel.C
deleted file mode 100644
index dbc7d6dc3..000000000
--- a/src/solidModels/thermalModel/thermalModel.C
+++ /dev/null
@@ -1,90 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | foam-extend: Open Source CFD
-   \\    /   O peration     | Version:     4.0
-    \\  /    A nd           | Web:         http://www.foam-extend.org
-     \\/     M anipulation  | For copyright notice see file Copyright
--------------------------------------------------------------------------------
-License
-    This file is part of foam-extend.
-
-    foam-extend is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation, either version 3 of the License, or (at your
-    option) any later version.
-
-    foam-extend is distributed in the hope that it will be useful, but
-    WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-    General Public License for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    thermalModel
-
-Description
-    Thermal  material properties for solids.
-
-Author
-    Hrvoje Jasak, Wikki Ltd.  All rights reserved.
-
-\*---------------------------------------------------------------------------*/
-
-#include "thermalModel.H"
-#include "volFields.H"
-#include "fvc.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-defineTypeNameAndDebug(thermalModel, 0);
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-thermalModel::thermalModel(const volScalarField& T)
-:
-    IOdictionary
-    (
-        IOobject
-        (
-            "thermalProperties",
-            T.time().constant(),
-            T.db(),
-            IOobject::MUST_READ_IF_MODIFIED,
-            IOobject::NO_WRITE
-        )
-    ),
-    T_(T),
-    lawPtr_(thermalLaw::New("law", T_, subDict("thermal")))
-{}
-
-
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
-
-bool thermalModel::read()
-{
-    if (regIOobject::read())
-    {
-        lawPtr_ = thermalLaw::New("law", T_, subDict("thermal"));
-
-        return true;
-    }
-    else
-    {
-        return false;
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// ************************************************************************* //
diff --git a/src/solidModels/thermalModel/thermalModel.H b/src/solidModels/thermalModel/thermalModel.H
deleted file mode 100644
index bd0358f1f..000000000
--- a/src/solidModels/thermalModel/thermalModel.H
+++ /dev/null
@@ -1,148 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | foam-extend: Open Source CFD
-   \\    /   O peration     | Version:     4.0
-    \\  /    A nd           | Web:         http://www.foam-extend.org
-     \\/     M anipulation  | For copyright notice see file Copyright
--------------------------------------------------------------------------------
-License
-    This file is part of foam-extend.
-
-    foam-extend is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation, either version 3 of the License, or (at your
-    option) any later version.
-
-    foam-extend is distributed in the hope that it will be useful, but
-    WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-    General Public License for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    thermalModel
-
-Description
-    Thermal material properties for solids.
-
-Author
-    Hrvoje Jasak, Wikki Ltd.  All rights reserved.
-
-SourceFiles
-    thermalModel.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef thermalModel_H
-#define thermalModel_H
-
-#include "IOdictionary.H"
-#include "typeInfo.H"
-#include "runTimeSelectionTables.H"
-#include "volFields.H"
-#include "tmp.H"
-#include "thermalLaw.H"
-#include "Switch.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                         Class thermalModel Declaration
-\*---------------------------------------------------------------------------*/
-
-class thermalModel
-:
-    public IOdictionary
-{
-    // Private data
-
-        //- Reference to temperature field
-        const volScalarField& T_;
-
-        //- Thermal law
-        autoPtr<thermalLaw> lawPtr_;
-
-
-    // Private Member Functions
-
-        //- Disallow copy construct
-        thermalModel(const thermalModel&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const thermalModel&);
-
-
-public:
-
-    //- Runtime type information
-    TypeName("thermalModel");
-
-    // Constructors
-
-        //- Construct from dictionary
-        thermalModel(const volScalarField& T);
-
-
-    // Destructor
-
-        virtual ~thermalModel()
-        {}
-
-
-    // Member Functions
-
-        //- Return thermal law
-        const thermalLaw& law() const
-        {
-            return lawPtr_();
-        }
-
-        //- Return specific heat capacity
-        tmp<volScalarField> C() const
-        {
-            return lawPtr_->C();
-        }
-
-        //- Return thermal conductivity
-        tmp<volScalarField> k() const
-        {
-            return lawPtr_->k();
-        }
-
-        //- Return thermal expansion coefficient
-        tmp<volScalarField> alpha() const
-        {
-            return lawPtr_->alpha();
-        }
-
-        //- Return reference temperature
-        tmp<volScalarField> T0() const
-        {
-            return lawPtr_->T0();
-        }
-
-        //- Correct the thermal model
-        void correct()
-        {
-            lawPtr_->correct();
-        }
-
-        //- Read thermalProperties dictionary
-        virtual bool read();
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

From 579268ae2d76f5a2e94e3b4e36ed613447fad4fe Mon Sep 17 00:00:00 2001
From: Vuko Vukcevic <v.vukcevic@wikki.co.uk>
Date: Wed, 9 May 2018 15:08:01 +0200
Subject: [PATCH 15/15] Updated tutorials for using global thermalModels

---
 .../conjugateHeatFoam/conjugateCavity/0/solid/materials      | 2 +-
 .../conjugateCavity/constant/solid/thermalProperties         | 4 ++++
 .../conjugateCavity/0/solid/materials                        | 2 +-
 .../conjugateCavity/constant/solid/thermalProperties         | 4 ++++
 .../conjugateHeatedSolid/constant/solid/thermalProperties    | 4 ++++
 .../hotCylinder/constant/polyMesh/boundary                   | 2 +-
 .../hotCylinder/constant/thermalProperties                   | 5 +++++
 7 files changed, 20 insertions(+), 3 deletions(-)

diff --git a/tutorials/coupled/conjugateHeatFoam/conjugateCavity/0/solid/materials b/tutorials/coupled/conjugateHeatFoam/conjugateCavity/0/solid/materials
index 3867d9034..f7e452ef5 100644
--- a/tutorials/coupled/conjugateHeatFoam/conjugateCavity/0/solid/materials
+++ b/tutorials/coupled/conjugateHeatFoam/conjugateCavity/0/solid/materials
@@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | foam-extend: Open Source CFD                    |
-|  \\    /   O peration     | Version:     4.0                                |
+|  \\    /   O peration     | Version:     4.1                                |
 |   \\  /    A nd           | Web:         http://www.foam-extend.org         |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
diff --git a/tutorials/coupled/conjugateHeatFoam/conjugateCavity/constant/solid/thermalProperties b/tutorials/coupled/conjugateHeatFoam/conjugateCavity/constant/solid/thermalProperties
index 1e109f27a..453abfdac 100644
--- a/tutorials/coupled/conjugateHeatFoam/conjugateCavity/constant/solid/thermalProperties
+++ b/tutorials/coupled/conjugateHeatFoam/conjugateCavity/constant/solid/thermalProperties
@@ -27,6 +27,10 @@ thermal
             C               C [0 2 -2 -1 0] 100;
             k               k [1 1 -3 -1  0] 10;
 
+            // Values not used
+            alpha           alpha [0 0 0 -1 0 0 0] 0.0;
+            T0              T0 [0 0 0 1 0 0 0] 0.0;
+
             zones           ( solidBlock );
         }
     );
diff --git a/tutorials/coupled/conjugateHeatSimpleFoam/conjugateCavity/0/solid/materials b/tutorials/coupled/conjugateHeatSimpleFoam/conjugateCavity/0/solid/materials
index 3867d9034..f7e452ef5 100644
--- a/tutorials/coupled/conjugateHeatSimpleFoam/conjugateCavity/0/solid/materials
+++ b/tutorials/coupled/conjugateHeatSimpleFoam/conjugateCavity/0/solid/materials
@@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | foam-extend: Open Source CFD                    |
-|  \\    /   O peration     | Version:     4.0                                |
+|  \\    /   O peration     | Version:     4.1                                |
 |   \\  /    A nd           | Web:         http://www.foam-extend.org         |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
diff --git a/tutorials/coupled/conjugateHeatSimpleFoam/conjugateCavity/constant/solid/thermalProperties b/tutorials/coupled/conjugateHeatSimpleFoam/conjugateCavity/constant/solid/thermalProperties
index 1e109f27a..453abfdac 100644
--- a/tutorials/coupled/conjugateHeatSimpleFoam/conjugateCavity/constant/solid/thermalProperties
+++ b/tutorials/coupled/conjugateHeatSimpleFoam/conjugateCavity/constant/solid/thermalProperties
@@ -27,6 +27,10 @@ thermal
             C               C [0 2 -2 -1 0] 100;
             k               k [1 1 -3 -1  0] 10;
 
+            // Values not used
+            alpha           alpha [0 0 0 -1 0 0 0] 0.0;
+            T0              T0 [0 0 0 1 0 0 0] 0.0;
+
             zones           ( solidBlock );
         }
     );
diff --git a/tutorials/coupled/conjugateHeatSimpleFoam/conjugateHeatedSolid/constant/solid/thermalProperties b/tutorials/coupled/conjugateHeatSimpleFoam/conjugateHeatedSolid/constant/solid/thermalProperties
index a76ecd825..aad9d08ae 100644
--- a/tutorials/coupled/conjugateHeatSimpleFoam/conjugateHeatedSolid/constant/solid/thermalProperties
+++ b/tutorials/coupled/conjugateHeatSimpleFoam/conjugateHeatedSolid/constant/solid/thermalProperties
@@ -20,6 +20,10 @@ thermal
     rho             rho [1 -3 0 0 0] 10;
     C               C [0 2 -2 -1 0] 100;
     k               k [1 1 -3 -1  0] 10;
+            
+    // Values not used
+    alpha           alpha [0 0 0 -1 0 0 0] 0.0;
+    T0              T0 [0 0 0 1 0 0 0] 0.0;
 
     gaps
     (
diff --git a/tutorials/solidMechanics/elasticThermalSolidFoam/hotCylinder/constant/polyMesh/boundary b/tutorials/solidMechanics/elasticThermalSolidFoam/hotCylinder/constant/polyMesh/boundary
index 87bebf7b8..1bbb32d19 100644
--- a/tutorials/solidMechanics/elasticThermalSolidFoam/hotCylinder/constant/polyMesh/boundary
+++ b/tutorials/solidMechanics/elasticThermalSolidFoam/hotCylinder/constant/polyMesh/boundary
@@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | foam-extend: Open Source CFD                    |
-|  \\    /   O peration     | Version:     4.0                                |
+|  \\    /   O peration     | Version:     4.1                                |
 |   \\  /    A nd           | Web:         http://www.foam-extend.org         |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
diff --git a/tutorials/solidMechanics/elasticThermalSolidFoam/hotCylinder/constant/thermalProperties b/tutorials/solidMechanics/elasticThermalSolidFoam/hotCylinder/constant/thermalProperties
index d7f384729..4365bad3c 100644
--- a/tutorials/solidMechanics/elasticThermalSolidFoam/hotCylinder/constant/thermalProperties
+++ b/tutorials/solidMechanics/elasticThermalSolidFoam/hotCylinder/constant/thermalProperties
@@ -21,6 +21,11 @@ thermal
     k        k [1 1 -3 -1 0 0 0] 250;
     alpha    alpha [0 0 0 -1 0 0 0] 2.3e-05;
     T0       T0 [0 0 0 1 0 0 0] 0;
+
+    // Data not used
+    rho      rho [1 -3 0 0 0 0 0] 1;
+    gaps     ();
+    sources  ();
 }