Merge branch 'nextRelease' of ssh://git.code.sf.net/p/foam-extend/foam-extend-4.0 into nextRelease

applications/solvers/immersedBoundary/pimpleDyMIbFoam/Make/options

@@ -19,6 +19,6 @@ EXE_LIBS = \
     -lincompressibleTurbulenceModel \
     -lincompressibleRASModels \
     -lincompressibleLESModels \
-    -limmersedBoundaryTurbulence \
+    -lincompressibleImmersedBoundaryTurbulence \
     -llduSolvers \
     -L$(MESQUITE_LIB_DIR) -lmesquite

applications/utilities/immersedBoundary/writeIbMasks/writeIbMasks.C

@@ -33,6 +33,7 @@ Description
 #include "fvc.H"
 #include "fvMatrices.H"
 #include "immersedBoundaryFvPatch.H"
+#include "cellSet.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -59,6 +60,9 @@ void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
     label minLiveCell = -1;
     const scalarField& gammaIn = gamma.internalField();
 
+    // Collect dead cells
+    labelHashSet deadCellsHash;
+    
     forAll (mesh.boundary(), patchI)
     {
         if (isA<immersedBoundaryFvPatch>(mesh.boundary()[patchI]))
@@ -79,6 +83,9 @@ void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
                     minLiveCell = ibCells[dcI];
                 }
             }
+
+            // Collect dead cells
+            deadCellsHash.insert(ibPatch.ibPolyPatch().deadCells());
         }
     }
 
@@ -124,6 +131,17 @@ void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
     sGamma.write();
     gamma.write();
 
+
+    // Create dead cells set
+    {
+        cellSet
+        (
+            mesh,
+            "deadCells",
+            deadCellsHash
+        ).write();
+    }
+    
     // Check consistency of face area vectors
 
     Info<< nl << "Calculating divSf" << endl;
@@ -240,8 +258,6 @@ void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
             << "Sum normal areas: " << sum(openFaceAreas) << nl
             << "Sum iB areas: " << sum(ibVectors) << nl
             << endl;
-
-
     }
 
     Info<< endl;

applications/utilities/mesh/generation/extrudeMesh/extrudeMesh.C

@@ -25,7 +25,8 @@ Description
     Extrude mesh from existing patch (by default outwards facing normals;
     optional flips faces) or from patch read from file.
 
-    Note: Merges close points so be careful.
+    Note: Merges close points so be careful. This can be controlled with the
+    optional mergeTolerance parameter (1e-4) by default.
 
     Type of extrusion prescribed by run-time selectable model.
 
@@ -179,7 +180,11 @@ int main(int argc, char *argv[])
 
     const boundBox& bb = mesh.globalData().bb();
     const vector span = bb.span();
-    const scalar mergeDim = 1e-4 * bb.minDim();
+
+    // Read merge tolerance
+    const scalar mergeTolerance =
+        dict.lookupOrDefault<scalar>("mergeTolerance", 1e-4);
+    const scalar mergeDim = mergeTolerance*bb.minDim();
 
     Info<< "Mesh bounding box : " << bb << nl
         << "        with span : " << span << nl
@@ -200,7 +205,7 @@ int main(int argc, char *argv[])
 
     // Collapse edges
     // ~~~~~~~~~~~~~~
-
+    if (mergeTolerance > SMALL)
     {
         Info<< "Collapsing edges < " << mergeDim << " ..." << nl << endl;
 

applications/utilities/mesh/generation/extrudeMesh/extrudeProperties

@@ -70,5 +70,7 @@ sigmaRadialCoeffs
     pStrat          1;
 }
 
+mergeTolerance      1e-4;
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/Allwmake

@@ -80,7 +80,8 @@ wmake libso solidModels
 wmake libso dbns
 
 wmake libso immersedBoundary/immersedBoundary
-wmake libso immersedBoundary/immersedBoundaryTurbulence
+wmake libso immersedBoundary/immersedBoundaryTurbulence/incompressible
+wmake libso immersedBoundary/immersedBoundaryTurbulence/compressible
 wmake libso immersedBoundary/immersedBoundaryDynamicMesh
 
 wmake libso overset/oversetMesh

src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/decomposeMesh.C

@@ -459,7 +459,7 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
                 // cell with multiple faces, we need to collect neighbour cell
                 // indices on the other side and possibly swap the order of
                 // adding faces. The same "swapping" of insertion order needs to
-                // happend on the slave side, but now we sort on the remote
+                // happen on the slave side, but now we sort on the remote
                 // (master) data and not on the local (slave) data as we do for
                 // master processor. VV, 16/Feb/2019.
 
@@ -588,7 +588,7 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
                         // above, owner and neighbour proc are the same. In
                         // order to avoid using cellToProc list for slave
                         // processor (where ownCellI is actually found on the
-                        // other side), use curNbrPtc with patchFace indes
+                        // other side), use curNbrPtc with patchFace index
                         const label& ownerProc = curNbrPtc[patchFaceI];
 
                         // Add the patch face into the list in the correct

src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/fvFieldDecomposerDecomposeFields.C

@@ -158,9 +158,17 @@ fvFieldDecomposer::decomposeField
 
         const label patchStart = field.mesh().boundaryMesh()[patchI].start();
 
-        forAll (p, i)
+        // Bugfix
+        // Special patches soch as overset and immersed boundary have
+        // their faces and values in a separate list; the list of
+        // values does not correspond to faces.  Skip such patches.
+        // HJ, 4/Jul/2019
+        if (p.size() == field.mesh().boundaryMesh()[patchI].size())
         {
-            allFaceField[patchStart + i] = p[i];
+            forAll (p, i)
+            {
+                allFaceField[patchStart + i] = p[i];
+            }
         }
     }
 

src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesRebuild.C

@@ -1728,8 +1728,337 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
     // due to the presence of old/new patches
     globalMesh.addFvPatches(reconPatches, false);
 
-    // TODO: point, face and cell zones
+    // Recombine cell, face and point zones.
+    // Note 1: all zones have to be present on same processors in the same
+    // order. This is the result of the decomposition. See
+    // domainDecomposition::processorMesh member function
+    // Note 2: the code below could be written in a generic way by using a
+    // template helper member function, but it's not straightforward since we
+    // don't have a list of ZoneMeshes for all processors
 
+    // Get index for the first valid mesh
+    const label fvmID = firstValidMesh();
+
+    // First pass: count maximum number of cells, faces and points in zones
+    labelList nCellsPerZone(meshes_[fvmID].cellZones().size(), 0);
+    labelList nFacesPerZone(meshes_[fvmID].faceZones().size(), 0);
+    labelList nPointsPerZone(meshes_[fvmID].pointZones().size(), 0);
+
+    forAll (meshes_, procI)
+    {
+        if (meshes_.set(procI))
+        {
+            // Grab the current mesh
+            const polyMesh& curMesh = meshes_[procI];
+
+            // PART 1: Cell zones. Scope for clarity and safety
+            {
+                const cellZoneMesh& cz = curMesh.cellZones();
+
+                forAll (cz, zoneI)
+                {
+                    // Count number of cells in the zone
+                    nCellsPerZone[zoneI] += cz[zoneI].size();
+                }
+            }
+
+            // PART 2: Face zones. Scope for clarity and safety
+            {
+                const faceZoneMesh& fz = curMesh.faceZones();
+
+                forAll (fz, zoneI)
+                {
+                    // Count number of faces in the zone
+                    nFacesPerZone[zoneI] += fz[zoneI].size();
+                }
+            }
+
+            // PART 3: Point zones. Scope for clarity and safety
+            {
+                const pointZoneMesh& pz = curMesh.pointZones();
+
+                forAll (pz, zoneI)
+                {
+                    // Count number of points in the zone
+                    nPointsPerZone[zoneI] += pz[zoneI].size();
+                }
+            }
+        } // End if processor mesh set
+    } // End for all processor meshes
+
+    // Second pass: redistribute cells, faces and points in zones
+
+    // Create lists that contain labels of all cells/faces/points in a given
+    // zone coming from different processor meshes. Index is the zoneID, which
+    // is the same for all processor bits, while the other list collects all the
+    // cells/faces/points in a given zone.
+    labelListList reconCellZones(nCellsPerZone.size());
+    labelListList reconFaceZones(nFacesPerZone.size());
+    List<boolList> reconFaceZoneFlips(nFacesPerZone.size());
+    labelListList reconPointZones(nPointsPerZone.size());
+
+    // Size the lists appropriatelly for each zone
+    forAll (reconCellZones, zoneI)
+    {
+        reconCellZones[zoneI].setSize(nCellsPerZone[zoneI], -1);
+    }
+    forAll (reconFaceZones, zoneI)
+    {
+        reconFaceZones[zoneI].setSize(nFacesPerZone[zoneI], -1);
+        reconFaceZoneFlips[zoneI].setSize(nFacesPerZone[zoneI], false);
+    }
+    forAll (reconPointZones, zoneI)
+    {
+        reconPointZones[zoneI].setSize(nPointsPerZone[zoneI], -1);
+    }
+
+    // Reset counting lists for indexing during list item assignement and for
+    // collecting the final number of faces/points in face/pointZones (since
+    // these can be actually fewer if e.g. a processor face becomes an internal
+    // face).
+    nCellsPerZone = 0;
+    nFacesPerZone = 0;
+    nPointsPerZone = 0;
+
+    // Loop through all the meshes and collect cells/faces/points in the zones
+    forAll (meshes_, procI)
+    {
+        if (meshes_.set(procI))
+        {
+            // Grab the current mesh
+            const polyMesh& curMesh = meshes_[procI];
+
+            // PART 1: Cell zones. Scope for clarity and safety
+            {
+                const cellZoneMesh& cz = curMesh.cellZones();
+
+                // Get old-to-new cell addressing for this mesh
+                const labelList& curCellProcAddr = cellProcAddressing_[procI];
+
+                forAll (cz, zoneI)
+                {
+                    // Get "new" zone cell index
+                    label& nCells = nCellsPerZone[zoneI];
+
+                    // Reference to the new recon zone
+                    labelList& zoneReconCells = reconCellZones[zoneI];
+
+                    // Get all the cells in this zone
+                    const labelList& zoneCells = cz[zoneI];
+
+                    // Loop through the cells
+                    forAll (zoneCells, i)
+                    {
+                        // Get cell index in the processor mesh
+                        const label& oldCellI = zoneCells[i];
+
+                        // Get cell index in the new mesh
+                        const label& newCellI = curCellProcAddr[oldCellI];
+
+                        // Redundancy check: check if the the newCellI is
+                        // -1. This should not happen because cells have perfect
+                        // 1-to-1 mapping
+                        if (newCellI != -1)
+                        {
+                            // Insert the cell in the new recon zone and
+                            // increment the counter
+                            zoneReconCells[nCells++] = newCellI;
+                        }
+                        else
+                        {
+                            FatalErrorIn
+                            (
+                                "autoPtr<fvMesh>"
+                                "\n processorMeshesReconstructor::"
+                                "reconstructMesh(const Time& db)"
+                            )   << "Found unmapped cell while reconstructing"
+                                << " cell zones."
+                                << nl
+                                << "Cell from processor: " << procI << nl
+                                << "Cell zone name: " << cz[zoneI].name() << nl
+                                << "Index in the cell zone: " << i << nl
+                                << "Old cell index: " << oldCellI << nl
+                                << "New cell index: " << newCellI
+                                << abort(FatalError);
+                        }
+
+                    } // End for all cells in the zone
+                } // End for all cell zones
+            } // End scope for cell zone handling
+
+            // PART 2: Face zones. Scope for clarity and safety
+            {
+                const faceZoneMesh& fz = curMesh.faceZones();
+
+                // Get old-to-new face addressing for this mesh
+                const labelList& curFaceProcAddr = faceProcAddressing_[procI];
+
+                forAll (fz, zoneI)
+                {
+                    // Get "new" zone face index
+                    label& nFaces = nFacesPerZone[zoneI];
+
+                    // Reference to the new recon zone and flips in the zone
+                    labelList& zoneReconFaces = reconFaceZones[zoneI];
+                    boolList& zoneReconFaceFlips = reconFaceZoneFlips[zoneI];
+
+                    // Get all the faces in this zone and also their flips
+                    const labelList& zoneFaces = fz[zoneI];
+                    const boolList& zoneFlipMap = fz[zoneI].flipMap();
+
+                    // Loop through the faces
+                    forAll (zoneFaces, i)
+                    {
+                        // Get the face index in the processor mesh
+                        const label& oldFaceI = zoneFaces[i];
+
+                        // Get the face index in the new, reconstructed mesh
+                        const label& newFaceI = curFaceProcAddr[oldFaceI];
+
+                        // Check if the face is mapped.
+                        // Note:
+                        // 1. Need to decrement by 1 because of the face turning
+                        // 2. No need to handle negative new indices coming from
+                        //    slave processor because we'd end up with
+                        //    duplicate entries (two faces on two processors
+                        //    merged into a single one)
+                        if (newFaceI > 0)
+                        {
+                            // This is a face that's been correctly
+                            // mapped, insert the face in the new zone
+                            zoneReconFaces[nFaces] = newFaceI - 1;
+
+                            // Also store the flip map of the face. We don't
+                            // need to check whether the flip map has been
+                            // preserved because we only get the combined faces
+                            // that are inserted from master side.
+                            zoneReconFaceFlips[nFaces] = zoneFlipMap[i];
+
+                            // Increment the number of faces for this zone
+                            ++nFaces;
+                        }
+                    } // End for all faces in the zone
+                } // End for all face zones
+            } // End scope for face zone handling
+
+            // PART 3: Point zones. Scope for clarity and safety
+            {
+                const pointZoneMesh& fz = curMesh.pointZones();
+
+                // Get old-to-new point addressing for this mesh
+                const labelList& curPointProcAddr = pointProcAddressing_[procI];
+
+                forAll (fz, zoneI)
+                {
+                    // Get "new" zone point index
+                    label& nPoints = nPointsPerZone[zoneI];
+
+                    // Reference to the new recon zone
+                    labelList& zoneReconPoints = reconPointZones[zoneI];
+
+                    // Get all the points in this zone
+                    const labelList& zonePoints = fz[zoneI];
+
+                    // Loop through the points
+                    forAll (zonePoints, i)
+                    {
+                        // Get point index in the processor mesh
+                        const label& oldPointI = zonePoints[i];
+
+                        // Get point index in the new mesh
+                        const label& newPointI = curPointProcAddr[oldPointI];
+
+                        // Check if the point is mapped
+                        if (newPointI != -1)
+                        {
+                            // Insert the point in the new recon zone and
+                            // increment the counter
+                            zoneReconPoints[nPoints++] = newPointI;
+                        }
+
+                    } // End for all points in the zone
+                } // End for all point zones
+            } // End scope for point zone handling
+        } // End if the processor mesh is set
+    } // End for all processor meshes
+
+    // We need to resize the face and point zones to number of inserted
+    // faces/points because not all faces and points need to be
+    // inserted. There's nothing to do for cell zones because these are always
+    // mapped uniquely one-to-one
+    forAll (reconFaceZones, zoneI)
+    {
+        reconFaceZones[zoneI].setSize(nFacesPerZone[zoneI]);
+        reconFaceZoneFlips[zoneI].setSize(nFacesPerZone[zoneI]);
+    }
+    forAll (reconPointZones, zoneI)
+    {
+        reconPointZones[zoneI].setSize(nPointsPerZone[zoneI]);
+    }
+
+    // Now we have all the zones as ordinary lists without possible duplicate
+    // faces and points due to merging of processor boundaries. Create zone
+    // meshes
+
+    // PART 1: Cell zones
+    List<cellZone*> reconCz(reconCellZones.size());
+
+    // Loop through all the cell zones and create them
+    forAll (reconCz, zoneI)
+    {
+        // Notes:
+        // 1. Grab the name from the respective zone in the first valid mesh
+        // 2. Transfer the list of cell IDs, invalidating reconCellZones[zoneI]
+        reconCz[zoneI] = new cellZone
+        (
+            meshes_[fvmID].cellZones()[zoneI].name(),
+            reconCellZones[zoneI].xfer(),
+            zoneI,
+            globalMesh.cellZones()
+        );
+    }
+
+    // PART 2: Face zones
+    List<faceZone*> reconFz(reconFaceZones.size());
+
+    // Loop through all the face zones and create them
+    forAll (reconFz, zoneI)
+    {
+        // Notes:
+        // 1. Grab the name from the respective zone in the first valid mesh
+        // 2. Transfer the list of face IDs, invalidating reconFaceZones[zoneI]
+        reconFz[zoneI] = new faceZone
+        (
+            meshes_[fvmID].faceZones()[zoneI].name(),
+            reconFaceZones[zoneI].xfer(),
+            reconFaceZoneFlips[zoneI].xfer(),
+            zoneI,
+            globalMesh.faceZones()
+        );
+    }
+
+    // PART 3: Point zones
+    List<pointZone*> reconPz(reconPointZones.size());
+
+    // Loop through all the point zones and create them
+    forAll (reconPz, zoneI)
+    {
+        // Notes:
+        // 1. Grab the name from the respective zone in the first valid mesh
+        // 2. Transfer the list of point IDs, invalidating reconPointZones[zoneI]
+        reconPz[zoneI] = new pointZone
+        (
+            meshes_[fvmID].pointZones()[zoneI].name(),
+            reconPointZones[zoneI].xfer(),
+            zoneI,
+            globalMesh.pointZones()
+        );
+    }
+
+    // Add the zones into the mesh
+    globalMesh.addZones(reconPz, reconFz, reconCz);
+
+    // All done, return the global mesh pointer
     return globalMeshPtr;
 }
 

src/dynamicMesh/dynamicMesh/directTopoChange/directTopoChange/directActions/hexRef8.C

@@ -1698,7 +1698,7 @@ Foam::hexRef8::hexRef8(const polyMesh& mesh)
             IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
-        labelList(mesh_.nCells(), 0)
+        labelField(mesh_.nCells(), 0)
     ),
     pointLevel_
     (
@@ -1711,7 +1711,7 @@ Foam::hexRef8::hexRef8(const polyMesh& mesh)
             IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
-        labelList(mesh_.nPoints(), 0)
+        labelField(mesh_.nPoints(), 0)
     ),
     level0Edge_(getLevel0EdgeLength()),
     history_
@@ -1779,8 +1779,8 @@ Foam::hexRef8::hexRef8(const polyMesh& mesh)
 Foam::hexRef8::hexRef8
 (
     const polyMesh& mesh,
-    const labelList& cellLevel,
-    const labelList& pointLevel,
+    const labelField& cellLevel,
+    const labelField& pointLevel,
     const refinementHistory& history
 )
 :
@@ -1833,8 +1833,8 @@ Foam::hexRef8::hexRef8
     {
         FatalErrorIn
         (
-            "hexRef8::hexRef8(const polyMesh&, const labelList&"
-            ", const labelList&, const refinementHistory&)"
+            "hexRef8::hexRef8(const polyMesh&, const labelField&"
+            ", const labelField&, const refinementHistory&)"
         )   << "History enabled but number of visible cells in it "
             << history_.visibleCells().size()
             << " is not equal to the number of cells in the mesh "
@@ -1849,8 +1849,8 @@ Foam::hexRef8::hexRef8
     {
         FatalErrorIn
         (
-            "hexRef8::hexRef8(const polyMesh&, const labelList&"
-            ", const labelList&, const refinementHistory&)"
+            "hexRef8::hexRef8(const polyMesh&, const labelField&"
+            ", const labelField&, const refinementHistory&)"
         )   << "Incorrect cellLevel or pointLevel size." << endl
             << "Number of cells in mesh:" << mesh_.nCells()
             << " does not equal size of cellLevel:" << cellLevel_.size() << endl
@@ -1876,8 +1876,8 @@ Foam::hexRef8::hexRef8
 Foam::hexRef8::hexRef8
 (
     const polyMesh& mesh,
-    const labelList& cellLevel,
-    const labelList& pointLevel
+    const labelField& cellLevel,
+    const labelField& pointLevel
 )
 :
     mesh_(mesh),

src/dynamicMesh/dynamicMesh/directTopoChange/directTopoChange/directActions/hexRef8.H

@@ -43,6 +43,7 @@ SourceFiles
 #include "removeFaces.H"
 #include "refinementHistory.H"
 #include "PackedList.H"
+#include "labelIOField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -68,10 +69,10 @@ class hexRef8
         const polyMesh& mesh_;
 
         //- Per cell the refinement level
-        labelIOList cellLevel_;
+        labelIOField cellLevel_;
 
         //- Per point the refinement level
-        labelIOList pointLevel_;
+        labelIOField pointLevel_;
 
         //- Typical edge length between unrefined points
         const scalar level0Edge_;
@@ -326,8 +327,8 @@ public:
         hexRef8
         (
             const polyMesh& mesh,
-            const labelList& cellLevel,
-            const labelList& pointLevel,
+            const labelField& cellLevel,
+            const labelField& pointLevel,
             const refinementHistory& history
         );
 
@@ -335,8 +336,8 @@ public:
         hexRef8
         (
             const polyMesh& mesh,
-            const labelList& cellLevel,
-            const labelList& pointLevel
+            const labelField& cellLevel,
+            const labelField& pointLevel
         );
 
 
@@ -344,12 +345,12 @@ public:
 
         // Access
 
-            const labelIOList& cellLevel() const
+            const labelIOField& cellLevel() const
             {
                 return cellLevel_;
             }
 
-            const labelIOList& pointLevel() const
+            const labelIOField& pointLevel() const
             {
                 return pointLevel_;
             }

src/dynamicMesh/dynamicMesh/directTopoChange/meshCut/modifiers/multiDirRefinement/multiDirRefinement.C

@@ -264,8 +264,8 @@ void Foam::multiDirRefinement::refineHex8
     hexRef8 hexRefiner
     (
         mesh,
-        labelList(mesh.nCells(), 0),    // cellLevel
-        labelList(mesh.nPoints(), 0),   // pointLevel
+        labelField(mesh.nCells(), 0),    // cellLevel
+        labelField(mesh.nPoints(), 0),   // pointLevel
         refinementHistory
         (
             IOobject

src/dynamicMesh/dynamicMesh/polyMeshModifiers/refinement/polyhedralRefinement/polyhedralRefinement.C

@@ -2055,15 +2055,6 @@ void Foam::polyhedralRefinement::setCellsToRefine
         }
     }
 
-    // Remove all cells that would exceed the maximum refinement level
-    forAll (refineCell, cellI)
-    {
-        if (refineCell[cellI] && (cellLevel_[cellI] + 1 > maxRefinementLevel_))
-        {
-            refineCell[cellI] = false;
-        }
-    }
-
     // Make sure that the refinement is face consistent (2:1 consistency) and
     // point consistent (4:1 consistency) if necessary
 

src/dynamicMesh/dynamicMesh/polyMeshModifiers/refinement/prismatic2DRefinement/prismatic2DRefinement.C

@@ -2614,15 +2614,6 @@ void Foam::prismatic2DRefinement::setCellsToRefine
         }
     }
 
-    // Remove all cells that exceed the maximum refinement level
-    forAll (refineCell, cellI)
-    {
-        if (refineCell[cellI] && (cellLevel_[cellI] + 1 > maxRefinementLevel_))
-        {
-            refineCell[cellI] = false;
-        }
-    }
-
     // Make sure that the refinement is face consistent (2:1 consistency) and
     // point consistent (4:1 consistency) if necessary
 

src/dynamicMesh/dynamicMesh/polyMeshModifiers/refinement/refinement/refinement.C

@@ -1333,20 +1333,6 @@ bool Foam::refinement::changeTopology() const
 
 void Foam::refinement::setRefinement(polyTopoChange& ref) const
 {
-    // Make sure that the point levels are updated across coupled patches before
-    // setting refinement and unrefinement. Note: not sure why the sync is not
-    // performed correctly if I do it in updateMesh. This is a temporary
-    // solution, need to investigate in detail, but I assume something is not
-    // updated yet in that case. VV, 31/Jan/2018.
-    syncTools::syncPointList
-    (
-        mesh_,
-        pointLevel_,
-        maxEqOp<label>(),
-        label(0),   // Null value
-        true        // Apply separation for parallel cyclics
-    );
-
     // Set refinement and unrefinement
     this->setRefinementInstruction(ref);
     this->setUnrefinementInstruction(ref);

src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicPolyRefinementFvMesh.C

@@ -117,89 +117,105 @@ Foam::dynamicPolyRefinementFvMesh::dynamicPolyRefinementFvMesh
     // protectedInitialRefinement)
     refinementSelectionPtr_()
 {
-    // Only create topo modifiers if they haven't been read in
-    // HJ, 16/Oct/2018
-    if (topoChanger_.empty())
+    // Check whether we read polyMeshModifiers from
+    // constant/polyMesh/meshModifiers file in the base class
+    if (!topoChanger_.empty())
     {
-        // Only one topo changer engine
-        topoChanger_.setSize(1);
+        // Already initialized, warn the user that we'll neglect it
+        WarningIn
+        (
+            "dynamicPolyRefinementFvMesh::dynamicPolyRefinementFvMesh"
+            "\n("
+            "\n    const IOobject& io"
+            "\n    const word subDictName"
+            "\n)"
+        )  << "Using controls from constant/dynamicMeshDict instead of"
+           << " constant/polyMesh/meshModifiers."
+           << nl
+           << "To supress this warning, delete meshModifiers file."
+           << endl;
 
-        // Get number of valid geometric dimensions
-        const label nGeometricDirs = this->nGeometricD();
-
-        switch(nGeometricDirs)
-        {
-            case 3:
-                // Add the polyhedralRefinement engine for
-                // 3D isotropic refinement
-                Info<< "3D case detected. "
-                    << "Adding polyhedralRefinement topology modifier" << endl;
-                topoChanger_.set
-                (
-                    0,
-                    new polyhedralRefinement
-                    (
-                        "polyhedralRefinement",
-                        refinementDict_,
-                        0,
-                        topoChanger_
-                    )
-                );
-                break;
-
-            case 2:
-                // Add the prismatic2DRefinement engine for
-                // 2D isotropic refinement
-                Info<< "2D case detected. "
-                    << "Adding prismatic2DRefinement topology modifier" << endl;
-                topoChanger_.set
-                (
-                    0,
-                    new prismatic2DRefinement
-                    (
-                        "prismatic2DRefinement",
-                        refinementDict_,
-                        0,
-                        topoChanger_
-                    )
-                );
-                break;
-
-            case 1:
-                FatalErrorIn
-                (
-                    "dynamicPolyRefinementFvMesh::dynamicPolyRefinementFvMesh"
-                    "\n("
-                    "\n    const IOobject& io,"
-                    "\n    const word subDictName"
-                    "\n)"
-                )   << "1D case detected. No valid refinement strategy is"
-                    <<  " available for 1D cases."
-                    << abort(FatalError);
-                break;
-
-            default:
-                FatalErrorIn
-                (
-                    "dynamicPolyRefinementFvMesh::dynamicPolyRefinementFvMesh"
-                    "\n("
-                    "\n    const IOobject& io,"
-                    "\n    const word subDictName"
-                    "\n)"
-                )   << "Invalid number of geometric meshes detected: "
-                    << nGeometricDirs
-                    << nl << "It appears that this mesh is neither 1D, 2D or 3D."
-                    << nl << " the mesh."
-                    << abort(FatalError);
-
-        }
-
-        // Write mesh and modifiers
-        topoChanger_.writeOpt() = IOobject::AUTO_WRITE;
-        topoChanger_.write();
-        write();
+       // Clear the list
+       topoChanger_.clear();
     }
 
+    // Only one topo changer engine
+    topoChanger_.setSize(1);
+
+    // Get number of valid geometric dimensions
+    const label nGeometricDirs = this->nGeometricD();
+
+    switch(nGeometricDirs)
+    {
+        case 3:
+            // Add the polyhedralRefinement engine for
+            // 3D isotropic refinement
+            Info<< "3D case detected. "
+                << "Adding polyhedralRefinement topology modifier" << endl;
+            topoChanger_.set
+            (
+                0,
+                new polyhedralRefinement
+                (
+                    "polyhedralRefinement",
+                    refinementDict_,
+                    0,
+                    topoChanger_
+                )
+            );
+            break;
+
+        case 2:
+            // Add the prismatic2DRefinement engine for
+            // 2D isotropic refinement
+            Info<< "2D case detected. "
+                << "Adding prismatic2DRefinement topology modifier" << endl;
+            topoChanger_.set
+            (
+                0,
+                new prismatic2DRefinement
+                (
+                    "prismatic2DRefinement",
+                    refinementDict_,
+                    0,
+                    topoChanger_
+                )
+            );
+            break;
+
+        case 1:
+            FatalErrorIn
+            (
+                "dynamicPolyRefinementFvMesh::dynamicPolyRefinementFvMesh"
+                "\n("
+                "\n    const IOobject& io,"
+                "\n    const word subDictName"
+                "\n)"
+            )   << "1D case detected. No valid refinement strategy is"
+                <<  " available for 1D cases."
+                << abort(FatalError);
+            break;
+
+        default:
+            FatalErrorIn
+            (
+                "dynamicPolyRefinementFvMesh::dynamicPolyRefinementFvMesh"
+                "\n("
+                "\n    const IOobject& io,"
+                "\n    const word subDictName"
+                "\n)"
+            )   << "Invalid number of geometric meshes detected: "
+                << nGeometricDirs
+                << nl << "It appears that this mesh is neither 1D, 2D or 3D."
+                << abort(FatalError);
+
+    }
+
+    // Write mesh and modifiers
+    topoChanger_.writeOpt() = IOobject::AUTO_WRITE;
+    topoChanger_.write();
+    write();
+
     // Initialize refinement selection algorithm after modifiers
     refinementSelectionPtr_ = refinementSelection::New(*this, refinementDict_);
 }
@@ -271,13 +287,13 @@ bool Foam::dynamicPolyRefinementFvMesh::update()
             (
                 refinementSelectionPtr_->refinementCellCandidates()()
             );
-            Pout<< "Selected " << refCandidates.size()
+            Info<< "Selected " << refCandidates.size()
                 << " refinement candidates."
                 << endl;
         }
         else
         {
-            Pout<< "Skipping refinement for this time-step..." << endl;
+            Info<< "Skipping refinement for this time-step..." << endl;
         }
 
         // Set cells to refine. Note: refinement needs to make sure that face
@@ -297,13 +313,13 @@ bool Foam::dynamicPolyRefinementFvMesh::update()
             (
                 refinementSelectionPtr_->unrefinementPointCandidates()()
             );
-            Pout<< "Selected " << unrefCandidates.size()
+            Info<< "Selected " << unrefCandidates.size()
                 << " unrefinement candidates."
                 << endl;
         }
         else
         {
-            Pout<< "Skipping unrefinement for this time-step..." << endl;
+            Info<< "Skipping unrefinement for this time-step..." << endl;
         }
 
         // Set split points to unrefine around.
@@ -344,13 +360,13 @@ bool Foam::dynamicPolyRefinementFvMesh::update()
         // some topo changes
         if (sizeCellMap)
         {
-            Pout<< "Successfully performed polyhedral refinement. "
+            Info<< "Successfully performed polyhedral refinement. "
                 << "Changed from " << nOldCells << " to " << sizeCellMap
                 << " cells." << endl;
         }
         else
         {
-            Pout<< "Refinement/unrefinement not performed in this time step "
+            Info<< "Refinement/unrefinement not performed in this time step "
                 << "since no cells were selected." << endl;
         }
 
@@ -362,7 +378,9 @@ bool Foam::dynamicPolyRefinementFvMesh::update()
         // per time step
         curTimeIndex_ = time().timeIndex();
     }
-    Pout<< "No refinement/unrefinement" << endl;
+
+    Info<< "No refinement/unrefinement" << endl;
+
     // No refinement/unrefinement at this time step. Return false
     return false;
 }

src/dynamicMesh/topoChangerFvMesh/topoChangerFvMesh/topoChangerFvMeshLoadBalance.C

@@ -291,7 +291,7 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
                 true    // Create passive processor patches
             );
             fvMesh& procMesh = procMeshPtr();
-            procMesh.write();
+
             // Create a field decomposer
             fvFieldDecomposer fieldDecomposer
             (
@@ -304,7 +304,11 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
 
             if (procI != Pstream::myProcNo())
             {
-                Pout<< "Send mesh and fields to processor " << procI << endl;
+                if (debug)
+                {
+                    Pout<< "Send mesh and fields to processor " << procI
+                        << endl;
+                }
 
                 OPstream toProc
                 (
@@ -449,8 +453,6 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
         }
     }
 
-    sleep(2);
-    
     // Collect pieces of mesh and fields from other processors
     for (label procI = 0; procI < meshDecomp.nProcs(); procI++)
     {
@@ -459,7 +461,10 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
             // Check if there is a mesh to send
             if (migratedCells[procI][Pstream::myProcNo()] > 0)
             {
-                Pout<< "Receive mesh and fields from " << procI << endl;
+                if (debug)
+                {
+                    Pout<< "Receive mesh and fields from " << procI << endl;
+                }
 
                 // Note: communication can be optimised.  HJ, 27/Feb/2018
                 IPstream fromProc
@@ -839,6 +844,9 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
         oldPatchNMeshPoints         // oldPatchNMeshPoints
     );
 
+    // Remove point, face and cell zones from the original mesh
+    removeZones();
+
     // Reset fvMesh and patches
     resetFvPrimitives
     (
@@ -852,6 +860,55 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
         true                       // Valid boundary
     );
 
+    // Get pointers to cell/face/point zones from reconstructed mesh and "link"
+    // them to the new mesh
+
+    // Cell zones
+    const cellZoneMesh& reconCellZones = reconMesh.cellZones();
+    List<cellZone*> czs(reconCellZones.size());
+    forAll (czs, zoneI)
+    {
+        czs[zoneI] = new cellZone
+        (
+            reconCellZones[zoneI].name(),
+            reconCellZones[zoneI],
+            zoneI,
+            this->cellZones()
+        );
+    }
+
+    // Face zones
+    const faceZoneMesh& reconFaceZones = reconMesh.faceZones();
+    List<faceZone*> fzs(reconFaceZones.size());
+    forAll (fzs, zoneI)
+    {
+        fzs[zoneI] = new faceZone
+        (
+            reconFaceZones[zoneI].name(),
+            reconFaceZones[zoneI],
+            reconFaceZones[zoneI].flipMap(),
+            zoneI,
+            this->faceZones()
+        );
+    }
+
+    // Point zones
+    const pointZoneMesh& reconPointZones = reconMesh.pointZones();
+    List<pointZone*> pzs(reconPointZones.size());
+    forAll (pzs, zoneI)
+    {
+        pzs[zoneI] = new pointZone
+        (
+            reconPointZones[zoneI].name(),
+            reconPointZones[zoneI],
+            zoneI,
+            this->pointZones()
+        );
+    }
+
+    // Add the zones to the mesh
+    addZones(pzs, fzs, czs);
+
     // Create field reconstructor
     fvFieldReconstructor fieldReconstructor
     (

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

@@ -35,6 +35,9 @@ bool adjustTimeStep =
 scalar maxCo =
     runTime.controlDict().lookupOrDefault<scalar>("maxCo", 1.0);
 
+scalar minDeltaT =
+    runTime.controlDict().lookupOrDefault<scalar>("minDeltaT", SMALL);
+
 scalar maxDeltaT =
     runTime.controlDict().lookupOrDefault<scalar>("maxDeltaT", GREAT);
 

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

@@ -35,6 +35,9 @@ adjustTimeStep =
 maxCo =
     runTime.controlDict().lookupOrDefault<scalar>("maxCo", 1.0);
 
+minDeltaT =
+    runTime.controlDict().lookupOrDefault<scalar>("minDeltaT", SMALL);
+
 maxDeltaT =
     runTime.controlDict().lookupOrDefault<scalar>("maxDeltaT", GREAT);
 

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

@@ -38,10 +38,14 @@ if (adjustTimeStep)
 
     runTime.setDeltaT
     (
-        min
+        Foam::max
         (
-            deltaTFact*runTime.deltaT().value(),
-            maxDeltaT
+            minDeltaT,
+            Foam::min
+            (
+                deltaTFact*runTime.deltaT().value(),
+                maxDeltaT
+            )
         )
     );
 

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

@@ -36,10 +36,14 @@ if (adjustTimeStep)
     {
         runTime.setDeltaT
         (
-            min
+            Foam::min
             (
-                maxCo*runTime.deltaT().value()/CoNum,
-                maxDeltaT
+                runTime.deltaT().value(),
+                Foam::min
+                (
+                    maxCo*runTime.deltaT().value()/CoNum,
+                    maxDeltaT
+                )
             )
         );
     }

src/finiteVolume/cfdTools/general/solutionControl/solutionControl/solutionControl.C

@@ -1019,11 +1019,18 @@ void Foam::solutionControl::reconstructTransientVelocity
     const volScalarField& p
 ) const
 {
+    // If the mesh is moving, flux needs to be relative before boundary
+    // conditions for velocity are corrected. VV and IG, 4/Jan/2016.
+    fvc::makeRelative(phi, U);
+
     // Reconstruct the velocity using all the components from original equation
     U = 1.0/(1.0/rAU + ddtUEqn.A())*
         (
             U/rAU + ddtUEqn.H() - fvc::grad(p)
         );
+    
+    // Correct boundary conditions with relative flux
+    U.correctBoundaryConditions();
 
     // Get name and the corresponding index
     const word& UName = U.name();
@@ -1067,13 +1074,6 @@ void Foam::solutionControl::reconstructTransientVelocity
     // Now that the normal component is zero, add the normal component from
     // conservative flux
     faceU += phi*rSf;
-
-    // If the mesh is moving, flux needs to be relative before boundary
-    // conditions for velocity are corrected. VV and IG, 4/Jan/2016.
-    fvc::makeRelative(phi, U);
-
-    // Correct boundary conditions with relative flux
-    U.correctBoundaryConditions();
 }
 
 

src/finiteVolume/fields/fvPatchFields/derived/pressureInletVelocity/pressureInletVelocityFvPatchVectorField.C

@@ -113,7 +113,7 @@ void Foam::pressureInletVelocityFvPatchVectorField::updateCoeffs()
         // Flux not available, do not update
         InfoIn
         (
-            "void ::pressureInletVelocityFvPatchVectorField::updateCoeffs()"
+            "void pressureInletVelocityFvPatchVectorField::updateCoeffs()"
         )   << "Flux field " << phiName_ << " not found.  "
             << "Performing fixed value update" << endl;
 

src/finiteVolume/finiteVolume/fvc/fvcReconstruct.C

@@ -76,32 +76,66 @@ reconstruct
     GeometricField<GradType, fvPatchField, volMesh>& reconField =
         treconField();
 
-    // Note:
-    // 1) Reconstruction is only available in cell centres: there is no need
+    // Notes regarding boundary:
+    // 1. Reconstruction is only available in cell centres: there is no need
     //    to invert the tensor on the boundary
-    // 2) For boundaries, the only reconstructed data is the flux times
-    //    normal.  Based on this guess, boundary conditions can adjust
-    //    patch values
-    // HJ, 12/Aug/2011
+    // 2. For boundaries, the only reconstructed data is the flux times
+    //    normal. Based on this guess, boundary conditions can adjust
+    //    patch values. HJ, 12/Aug/2011
 
-    GeometricField<GradType, fvPatchField, volMesh> fluxTimesNormal =
-        surfaceSum((mesh.Sf()/mesh.magSf())*ssf);
+    // Notes regarding procedure:
+    // 1. hinv inverse must be used to stabilise the inverse on bad meshes
+    //    but it gives strange failures because it unnecessarily avoids
+    //    performing ordinary inverse for meshes with reasonably sized
+    //    determinant (e.g. if SfSf/magSf is small). HJ, 19/Aug/2015
+    // 2. hinv has been stabilised now. HJ, 22/Mar/2019
+    // 3. But we still need to make sure that the determinant is not extremely
+    //    small, which may happen for extremely small meshes. We avoid this
+    //    issue by dividing the reconstruction equation with magSf^2 (instead of
+    //    magSf), which basically makes the dyadic tensor that we need to invert
+    //    dimensionless. VV, 13/Jun/2019
 
-    // Note: hinv inverse must be used to stabilise the inverse on bad meshes
-    // but it gives strange failures.  Fixed hinv.  HJ, 22/Mar/2019
-    // HJ, 19/Aug/2015
+    // Here's a short derivation in a Latex--like notation, where:
+    // - Sf is the surface area vector
+    // - uf is the face velocity factor (or field to be reconstructed)
+    // - Ff is the face flux
+    // - magSf is the magnitude of the surface area vector
+    // - uP is the velocity field in the cell centre
+    // - G is the surface area dyadic tensor
+    // - Fn is the vector representing surface sum of directional fluxes
+    // - \dprod is a dot product
+
+    // 1. Sf \dprod uf = Ff
+    //    Multiply Eq (1) with Sf/magSf^2
+    // 2. \frac{Sf Sf}{magSf^2} \dprod uf = \frac{Sf Ff}{magSf^2}
+    //    Sum Eq (2) over all the faces
+    // 3. \sum_f(\frac{Sf Sf}{magSf^2} \dprod uf) = \sum_f(\frac{Sf Ff}{magSf^2})
+    //    Assume first order extrapolation of uf, e.g. uP = uf
+    // 4. \sum_f(\frac{Sf Sf}{magSf^2}) \dprod uP) = \sum_f(\frac{Sf Ff}{magSf^2})
+    //    Use shorthand notation
+    // 5. G \dprod uP = Fn
+    // 6. uP = G^-1 \dprod Fn
+
+    // Fn -> fluxTimesNormal
+    // G  -> G
+
+    // Calculate sum of the directional fluxes
+    const surfaceScalarField magSfSqr = sqr(mesh.magSf());
+    const GeometricField<GradType, fvPatchField, volMesh> fluxTimesNormal =
+        surfaceSum((mesh.Sf()/magSfSqr)*ssf);
+
+    // Calculate the G tensor
+    const volSymmTensorField G = surfaceSum(sqr(mesh.Sf())/magSfSqr);
+
+    // Finally calculate the reconstructed field using hinv for stabilisation on
+    // really bad fvMesh bits (uses ordinary inverse most of the time, see
+    // tensor.C)
     reconField.internalField() =
-    (
-        hinv
-        (
-            surfaceSum(sqr(mesh.Sf())/mesh.magSf())().internalField()
-        )
-      & fluxTimesNormal.internalField()
-    );
+        hinv(G.internalField()) & fluxTimesNormal.internalField();
 
+    // Boundary value update
     reconField.boundaryField() = fluxTimesNormal.boundaryField();
-
-    treconField().correctBoundaryConditions();
+    reconField.correctBoundaryConditions();
 
     return treconField;
 }

src/foam/meshes/polyMesh/polyMesh.C

@@ -43,6 +43,20 @@ defineTypeNameAndDebug(Foam::polyMesh, 0);
 Foam::word Foam::polyMesh::defaultRegion = "region0";
 Foam::word Foam::polyMesh::meshSubDir = "polyMesh";
 
+const Foam::debug::tolerancesSwitch
+Foam::polyMesh::emptyDirTol_
+(
+    "emptyDirectionTolerance",
+    1e-10
+);
+
+const Foam::debug::tolerancesSwitch
+Foam::polyMesh::wedgeDirTol_
+(
+    "wedgeDirectionTolerance",
+    1e-10
+);
+
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -93,20 +107,27 @@ void Foam::polyMesh::calcDirections() const
         globalEmptyDirVec /= mag(globalEmptyDirVec);
 
         // Check if all processors see the same 2-D from empty patches
-        if (mag(globalEmptyDirVec - emptyDirVec) > SMALL)
+        if (mag(globalEmptyDirVec - emptyDirVec) > emptyDirTol_())
         {
             FatalErrorIn
             (
                 "void polyMesh::calcDirections() const"
             )   << "Some processors detect different empty (2-D) "
                 << "directions.  Probably using empty patches on a "
-                << "bad parallel decomposition." << nl
-                << "Please check your geometry and empty patches"
+                << "bad parallel decomposition." << nl << nl
+                << "Global empty direction vector: " << globalEmptyDirVec << nl
+                << "Local empty direction vector: " << emptyDirVec << nl
+                << "If global and local empty direction vectors are different"
+                << " to a round-off error, try increasing the"
+                << " emptyDirectionTolerance in controlDict::Tolerances" << nl
+                << "The emptyDirectionTolerance for this run is: "
+                << emptyDirTol_() << nl << nl
+                << "Otherwise please check your geometry and empty patches."
                 << abort(FatalError);
         }
     }
 
-    if (mag(emptyDirVec) > SMALL)
+    if (mag(emptyDirVec) > emptyDirTol_())
     {
         for (direction cmpt = 0; cmpt < vector::nComponents; cmpt++)
         {
@@ -139,7 +160,7 @@ void Foam::polyMesh::calcDirections() const
         globalWedgeDirVec /= mag(globalWedgeDirVec);
 
         // Check if all processors see the same 2-D from wedge patches
-        if (mag(globalWedgeDirVec - wedgeDirVec) > SMALL)
+        if (mag(globalWedgeDirVec - wedgeDirVec) > wedgeDirTol_())
         {
             FatalErrorIn
             (
@@ -147,12 +168,20 @@ void Foam::polyMesh::calcDirections() const
             )   << "Some processors detect different wedge (2-D) "
                 << "directions.  Probably using wedge patches on a "
                 << "bad parallel decomposition." << nl
+                << "Global wedge direction vector: " << globalWedgeDirVec << nl
+                << "Local wedge direction vector: " << wedgeDirVec << nl
+                << "If global and local wedge direction vectors are different"
+                << " to a round-off error, try increasing the"
+                << " wedgeDirectionTolerance in controlDict::Tolerances" << nl
+                << "The wedgeDirectionTolerance for this run is: "
+                << wedgeDirTol_() << nl << nl
+                << "Otherwise please check your geometry and empty patches."
                 << "Please check your geometry and wedge patches"
                 << abort(FatalError);
         }
     }
 
-    if (mag(wedgeDirVec) > SMALL)
+    if (mag(wedgeDirVec) > wedgeDirTol_())
     {
         for (direction cmpt = 0; cmpt < vector::nComponents; cmpt++)
         {

src/foam/meshes/polyMesh/polyMesh.H

@@ -235,11 +235,23 @@ public:
     //- Runtime type information
     TypeName("polyMesh");
 
-    //- Return the default region name
-    static word defaultRegion;
 
-    //- Return the mesh sub-directory name (usually "polyMesh")
-    static word meshSubDir;
+    //- Static mesh data
+
+        //- Return the default region name
+        static word defaultRegion;
+
+        //- Return the mesh sub-directory name (usually "polyMesh")
+        static word meshSubDir;
+
+
+    //- Static data to control empty and wedge directions
+
+        //- Empty direction tolerance
+        static const debug::tolerancesSwitch emptyDirTol_;
+
+        //- Wedge direction tolerance
+        static const debug::tolerancesSwitch wedgeDirTol_;
 
 
     // Constructors

src/foam/meshes/primitiveMesh/PrimitivePatchTemplate/PrimitivePatchMeshData.C

@@ -73,53 +73,27 @@ calcMeshData() const
     // number of faces in the patch
     Map<label> markedPoints(4*this->size());
 
-
-    // Important:
-    // ~~~~~~~~~~
-    // In <= 1.5 the meshPoints would be in increasing order but this gives
-    // problems in processor point synchronisation where we have to find out
-    // how the opposite side would have allocated points.
-
-    // Note:
-    // ~~~~~
-    // This is all garbage.  All -ext versions will preserve strong ordering
-    // HJ, 17/Aug/2010
-
-    //- 1.5 code:
-    // If the point is used, set the mark to 1
+    // Note: foam-extend was relying on strong ordering of meshPoints on
+    // processor patches which causes point ordering errors when doing Dynamic
+    // Load Balancing. Because of this ordering error on two sides of processor
+    // boundaries, point fields ended up having wrong values near processor
+    // boundaries. Revert to unsorted version. VV, 11/July/2019.
+    dynamicLabelList meshPoints(2*this->size());
     forAll(*this, facei)
     {
         const Face& curPoints = this->operator[](facei);
 
         forAll(curPoints, pointi)
         {
-            markedPoints.insert(curPoints[pointi], -1);
+            if (markedPoints.insert(curPoints[pointi], meshPoints.size()))
+            {
+                meshPoints.append(curPoints[pointi]);
+            }
         }
     }
 
-    // Create the storage and store the meshPoints.  Mesh points are
-    // the ones marked by the usage loop above
-    meshPointsPtr_ = new labelList(markedPoints.toc());
-    labelList& pointPatch = *meshPointsPtr_;
-
-    // Sort the list to preserve compatibility with the old ordering
-    sort(pointPatch);
-
-    // For every point in map give it its label in mesh points
-    forAll(pointPatch, pointi)
-    {
-        markedPoints.find(pointPatch[pointi])() = pointi;
-    }
-
-    forAll(*this, faceI)
-    {
-        const Face& curPoints = this->operator[](faceI);
-
-        forAll (curPoints, pointI)
-        {
-            markedPoints.insert(curPoints[pointI], -1);
-        }
-    }
+    // Transfer to straight list (reuses storage)
+    meshPointsPtr_ = new labelList(meshPoints, true);
 
     // Create local faces. Note that we start off from copy of original face
     // list (even though vertices are overwritten below). This is done so

src/immersedBoundary/CMakeLists.txt

@@ -32,3 +32,4 @@
 
 add_subdirectory(immersedBoundary)
 add_subdirectory(immersedBoundaryDynamicMesh)
+add_subdirectory(immersedBoundaryTurbulence)

src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatch/immersedBoundaryFvPatch.C

@@ -124,7 +124,7 @@ void Foam::immersedBoundaryFvPatch::updatePhi
     scalar rDeltaT = 1.0/deltaT;
 
     // Multiply the raw mesh motion flux with the masking function
-    scalarField sGamma =
+    scalarField ibAreaRatio =
         mag(ibPolyPatch_.correctedFaceAreas())/mag(mesh.faceAreas());
 
     // Scaling of internal mesh flux field should be done only for the current
@@ -140,7 +140,7 @@ void Foam::immersedBoundaryFvPatch::updatePhi
         const label faceI = deadFaces[dfI];
         if (mesh.isInternalFace(faceI))
         {
-            phiIn[faceI] *= sGamma[faceI];
+            phiIn[faceI] = scalar(0);
         }
     }
  
@@ -150,7 +150,7 @@ void Foam::immersedBoundaryFvPatch::updatePhi
         const label faceI = cutFaces[cfI];
         if (mesh.isInternalFace(faceI))
         {
-            phiIn[faceI] *= sGamma[faceI];
+            phiIn[faceI] *= ibAreaRatio[faceI];
         }
     }
 
@@ -160,7 +160,7 @@ void Foam::immersedBoundaryFvPatch::updatePhi
         if (!isA<immersedBoundaryFvPatch>(mesh.boundary()[patchI]))
         {
             phi.boundaryField()[patchI] *=
-                mesh.boundary()[patchI].patchSlice(sGamma);
+                mesh.boundary()[patchI].patchSlice(ibAreaRatio);
         }
     }
 
@@ -231,12 +231,11 @@ void Foam::immersedBoundaryFvPatch::updatePhi
             // Attempt to correct via old volume
             scalar corrOldVol = newVols[cellI] - divPhi[cellI]*deltaT;
 
-            // Info<< "Flux maneouvre for cell " << cellI << ": "
+            // Pout<< "Flux maneouvre for cell " << cellI << ": "
             //     << " error: " << magDivPhi[cellI]
             //     << " V: " << newVols[cellI]
             //     << " V0: " << oldVols[cellI]
-            //     << " divPhi: " << divPhi[cellI]
-            //     << endl;
+            //     << " divPhi: " << divPhi[cellI];
 
             if (corrOldVol < SMALL)
             {
@@ -248,6 +247,10 @@ void Foam::immersedBoundaryFvPatch::updatePhi
             {
                 oldVols[cellI] = corrOldVol;
             }
+
+            // scalar corrDivMeshPhi =
+            //     mag((newVols[cellI] - oldVols[cellI]) - divPhi[cellI]*deltaT);
+            // Pout<< " Corrected: " << corrDivMeshPhi << endl;
         }
     }
 }

src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatchFields/basic/fixedValueIbFvPatchField/fixedValueIbFvPatchField.C

@@ -81,17 +81,8 @@ fixedValueIbFvPatchField<Type>::fixedValueIbFvPatchField
 
     if (!isType<immersedBoundaryFvPatch>(p))
     {
-        FatalIOErrorIn
-        (
-            "fixedValueIbFvPatchField<Type>::"
-            "fixedValueIbFvPatchField\n"
-            "(\n"
-            "    const fvPatch& p,\n"
-            "    const Field<Type>& field,\n"
-            "    const dictionary& dict\n"
-            ")\n",
-            dict
-        )   << "\n    patch type '" << p.type()
+        FatalIOErrorInFunction(dict)
+            << "\n    patch type '" << p.type()
             << "' not constraint type '" << typeName << "'"
             << "\n    for patch " << p.name()
             << " of field " << this->dimensionedInternalField().name()
@@ -128,17 +119,8 @@ fixedValueIbFvPatchField<Type>::fixedValueIbFvPatchField
     // HJ, 12/Apr/2012
     if (!isType<immersedBoundaryFvPatch>(p))
     {
-        FatalErrorIn
-        (
-            "fixedValueIbFvPatchField<Type>::"
-            "fixedValueIbFvPatchField\n"
-            "(\n"
-            "    const fixedValueIbFvPatchField<Type>&,\n"
-            "    const fvPatch& p,\n"
-            "    const DimensionedField<Type, volMesh>& iF,\n"
-            "    const fvPatchFieldMapper& mapper\n"
-            ")\n"
-        )   << "\n    patch type '" << p.type()
+        FatalErrorInFunction
+            << "\n    patch type '" << p.type()
             << "' not constraint type '" << typeName << "'"
             << "\n    for patch " << p.name()
             << " of field " << this->dimensionedInternalField().name()
@@ -257,6 +239,16 @@ void fixedValueIbFvPatchField<Type>::evaluate
 }
 
 
+template<class Type>
+void Foam::fixedValueIbFvPatchField<Type>::manipulateMatrix
+(
+    fvMatrix<Type>& matrix
+)
+{
+    this->setDeadValues(matrix);
+}
+
+
 template<class Type>
 void fixedValueIbFvPatchField<Type>::write(Ostream& os) const
 {

src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatchFields/basic/fixedValueIbFvPatchField/fixedValueIbFvPatchField.H

@@ -191,6 +191,9 @@ public:
                 const Pstream::commsTypes commsType = Pstream::blocking
             );
 
+            //- Manipulate a matrix
+            virtual void manipulateMatrix(fvMatrix<Type>& matrix);
+
 
         // I-O
 

src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatchFields/basic/mixedIbFvPatchField/mixedIbFvPatchField.C

@@ -87,17 +87,8 @@ mixedIbFvPatchField<Type>::mixedIbFvPatchField
 
     if (!isType<immersedBoundaryFvPatch>(p))
     {
-        FatalIOErrorIn
-        (
-            "mixedIbFvPatchField<Type>::"
-            "mixedIbFvPatchField\n"
-            "(\n"
-            "    const fvPatch& p,\n"
-            "    const Field<Type>& field,\n"
-            "    const dictionary& dict\n"
-            ")\n",
-            dict
-        )   << "\n    patch type '" << p.type()
+        FatalIOErrorInFunction(dict)
+            << "\n    patch type '" << p.type()
             << "' not constraint type '" << typeName << "'"
             << "\n    for patch " << p.name()
             << " of field " << this->dimensionedInternalField().name()
@@ -109,13 +100,11 @@ mixedIbFvPatchField<Type>::mixedIbFvPatchField
     // the patch is active
     // Initialise the value to avoid errors
     // HJ, 1/Dec/2017
-//    if (this->ibPatch().ibPolyPatch().active())
-    {
-        // Re-interpolate the data related to immersed boundary
-        this->updateIbValues();
 
-        mixedFvPatchField<Type>::evaluate();
-    }
+    // Re-interpolate the data related to immersed boundary
+    this->updateIbValues();
+
+    mixedFvPatchField<Type>::evaluate();
 }
 
 
@@ -143,17 +132,8 @@ mixedIbFvPatchField<Type>::mixedIbFvPatchField
     // HJ, 12/Apr/2012
     if (!isType<immersedBoundaryFvPatch>(p))
     {
-        FatalErrorIn
-        (
-            "mixedIbFvPatchField<Type>::"
-            "mixedIbFvPatchField\n"
-            "(\n"
-            "    const mixedIbFvPatchField<Type>&,\n"
-            "    const fvPatch& p,\n"
-            "    const DimensionedField<Type, volMesh>& iF,\n"
-            "    const fvPatchFieldMapper& mapper\n"
-            ")\n"
-        )   << "\n    patch type '" << p.type()
+        FatalErrorInFunction
+            << "\n    patch type '" << p.type()
             << "' not constraint type '" << typeName << "'"
             << "\n    for patch " << p.name()
             << " of field " << this->dimensionedInternalField().name()
@@ -277,6 +257,16 @@ void mixedIbFvPatchField<Type>::evaluate
 }
 
 
+template<class Type>
+void Foam::mixedIbFvPatchField<Type>::manipulateMatrix
+(
+    fvMatrix<Type>& matrix
+)
+{
+    this->setDeadValues(matrix);
+}
+
+
 template<class Type>
 void mixedIbFvPatchField<Type>::write(Ostream& os) const
 {

src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatchFields/basic/mixedIbFvPatchField/mixedIbFvPatchField.H

@@ -221,6 +221,9 @@ public:
                 const Pstream::commsTypes commsType = Pstream::blocking
             );
 
+            //- Manipulate a matrix
+            virtual void manipulateMatrix(fvMatrix<Type>& matrix);
+
 
         // I-O
 

src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatchFields/derived/movingImmersedBoundaryVelocity/movingImmersedBoundaryVelocityFvPatchVectorField.C

@@ -220,6 +220,15 @@ void Foam::movingImmersedBoundaryVelocityFvPatchVectorField::evaluate
 }
 
 
+void Foam::movingImmersedBoundaryVelocityFvPatchVectorField::manipulateMatrix
+(
+    fvVectorMatrix& matrix
+)
+{
+    setDeadValues(matrix);
+}
+
+
 void Foam::movingImmersedBoundaryVelocityFvPatchVectorField::write
 (
     Ostream& os

src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatchFields/derived/movingImmersedBoundaryVelocity/movingImmersedBoundaryVelocityFvPatchVectorField.H

@@ -163,6 +163,9 @@ public:
                 const Pstream::commsTypes commsType = Pstream::blocking
             );
 
+            //- Manipulate a matrix
+            virtual void manipulateMatrix(fvVectorMatrix& matrix);
+
 
         // I-O
 

src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatchFields/immersedBoundaryFieldBase/immersedBoundaryFieldBase.C

@@ -63,6 +63,39 @@ void Foam::immersedBoundaryFieldBase<Type>::setDeadValues
 }
 
 
+template<class Type>
+void Foam::immersedBoundaryFieldBase<Type>::setDeadValues
+(
+    fvMatrix<Type>& matrix
+) const
+{
+    // Fix the value in dead cells
+    if (setDeadValue_)
+    {
+        const labelList& dc = ibPatch_.ibPolyPatch().deadCells();
+
+        // Boost the diagonal of dead cells by the volume ratio
+        // Volume ratio is set to SMALL; revert for diagonal
+        // This should also guarantee strong diagonal dominance.
+        // HJ, 19/Jun/2019
+
+        scalarField& diag = matrix.diag();
+
+        forAll (dc, dcI)
+        {
+            diag[dc[dcI]] *= GREAT;
+        }
+
+        // Set values
+        matrix.setValues
+        (
+            dc,
+            Field<Type>(dc.size(), deadValue_)
+        );
+    }
+}
+
+
 template<class Type>
 void Foam::immersedBoundaryFieldBase<Type>::writeDeadData(Ostream& os) const
 {
@@ -107,30 +140,65 @@ void Foam::immersedBoundaryFieldBase<Type>::writeField
         if (Pstream::master())
         {
             // Assemble unique lists to correspond to a single surface
-            pointField allPoints(0);
-            Field<Type> completeField(0);
-            faceList allFaces(0);
-            label prevProcPatchSize = 0;
 
-            forAll(procPoints, procI)
+            // Count points and faces; currently unmerged
+            label nAllPoints = 0;
+            label nAllFaces = 0;
+
+            forAll (procPoints, procI)
             {
-                allPoints.append(procPoints[procI]);
-                completeField.append(procFields[procI]);
+                nAllPoints += procPoints[procI].size();
+                nAllFaces += procFaces[procI].size();
+            }
 
-                // Point labels in faces need to be incremented with respect to
-                // the size of the size of the previous processore patch
-                forAll(procFaces[procI], faceI)
+            pointField allPoints(nAllPoints);
+            faceList allFaces(nAllFaces);
+            Field<Type> completeField(nAllFaces);
+
+            // Reset counters
+            nAllPoints = 0;
+            nAllFaces = 0;
+
+            forAll (procPoints, procI)
+            {
+                const pointField& curPoints = procPoints[procI];
+                labelList renumberPoints(curPoints.size());
+
+                forAll (curPoints, cpI)
                 {
-                    face curFace = procFaces[procI][faceI];
-                    forAll(curFace, pointI)
-                    {
-                        curFace[pointI] += prevProcPatchSize;
-                    }
-                    allFaces.append(curFace);
+                    allPoints[nAllPoints] = curPoints[cpI];
+                    renumberPoints[cpI] = nAllPoints;
+                    nAllPoints++;
                 }
 
-                // Increment the total number of points
-                prevProcPatchSize += procPoints[procI].size();
+                const faceList& curFaces = procFaces[procI];
+                const Field<Type>& curField = procFields[procI];
+
+                forAll (curFaces, cfI)
+                {
+                    // Point labels in faces need to be renumbered with respect
+                    // to the size of the size of the previous processore patch
+                    const face& oldFace = curFaces[cfI];
+
+                    // Make a copy of face to renumber
+                    face renumberedFace(oldFace.size());
+
+                    forAll (oldFace, fpI)
+                    {
+                        renumberedFace[fpI] = renumberPoints[oldFace[fpI]];
+                    }
+
+                    allFaces[nAllFaces] = renumberedFace;
+                    completeField[nAllFaces] = curField[cfI];
+                    nAllFaces++;
+                }
+            }
+
+            if (nAllPoints != allPoints.size() || nAllFaces != allFaces.size())
+            {
+                FatalErrorInFunction
+                    << "Problem with merge of immersed boundary patch data"
+                    << abort(FatalError);
             }
 
             writerPtr->write

src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatchFields/immersedBoundaryFieldBase/immersedBoundaryFieldBase.H

@@ -123,6 +123,9 @@ public:
             //- Set values in dead cells
             void setDeadValues(Field<Type>& psiI) const;
 
+            //- Set matrix constraints in dead cells
+            void setDeadValues(fvMatrix<Type>& matrix) const;
+
 
         // I/O
 

src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatchFields/immersedBoundaryFvPatchField/immersedBoundaryFvPatchField.C

@@ -64,17 +64,8 @@ immersedBoundaryFvPatchField<Type>::immersedBoundaryFvPatchField
 {
     if (!isType<immersedBoundaryFvPatch>(p))
     {
-        FatalIOErrorIn
-        (
-            "immersedBoundaryFvPatchField<Type>::"
-            "immersedBoundaryFvPatchField\n"
-            "(\n"
-            "    const fvPatch& p,\n"
-            "    const Field<Type>& field,\n"
-            "    const dictionary& dict\n"
-            ")\n",
-            dict
-        )   << "\n    patch type '" << p.type()
+        FatalIOErrorInFunction(dict)
+            << "\n    patch type '" << p.type()
             << "' not constraint type '" << typeName << "'"
             << "\n    for patch " << p.name()
             << " of field " << this->dimensionedInternalField().name()
@@ -107,17 +98,8 @@ immersedBoundaryFvPatchField<Type>::immersedBoundaryFvPatchField
     // HJ, 12/Apr/2012
     if (!isType<immersedBoundaryFvPatch>(p))
     {
-        FatalErrorIn
-        (
-            "immersedBoundaryFvPatchField<Type>::"
-            "immersedBoundaryFvPatchField\n"
-            "(\n"
-            "    const immersedBoundaryFvPatchField<Type>&,\n"
-            "    const fvPatch& p,\n"
-            "    const DimensionedField<Type, volMesh>& iF,\n"
-            "    const fvPatchFieldMapper& mapper\n"
-            ")\n"
-        )   << "\n    patch type '" << p.type()
+        FatalErrorInFunction
+            << "\n    patch type '" << p.type()
             << "' not constraint type '" << typeName << "'"
             << "\n    for patch " << p.name()
             << " of field " << this->dimensionedInternalField().name()

src/immersedBoundary/immersedBoundary/immersedBoundaryPolyPatch/immersedBoundaryPolyPatch.C

@@ -1262,8 +1262,8 @@ void Foam::immersedBoundaryPolyPatch::calcCorrectedGeometry() const
 
     forAll (dc, dcI)
     {
-        // Set dead volume to small
-        V[dc[dcI]] = SMALL;
+        // Scale dead volume to small
+        V[dc[dcI]] *= SMALL;
     }
 
     // Correct for all cut faces
@@ -1286,7 +1286,7 @@ void Foam::immersedBoundaryPolyPatch::calcCorrectedGeometry() const
 
     forAll (df, dfI)
     {
-        // Set dead area to small
+        // Scale dead area to small
         Sf[df[dfI]] *= SMALL;
     }
 

src/immersedBoundary/immersedBoundary/immersedPoly/ImmersedFace.C

@@ -257,10 +257,11 @@ void Foam::ImmersedFace<Distance>::createSubfaces
         // than 3 points
         if (nDry < 3)
         {
-            FatalErrorInFunction
-                << "There are fewer than three points"
-                << " on the wet part of the face."
-                << abort(FatalError);
+            // The face is wet
+            isAllWet_ = true;
+            isAllDry_ = false;
+
+            drySubface_.clear();
         }
         else
         {
@@ -272,10 +273,11 @@ void Foam::ImmersedFace<Distance>::createSubfaces
 
         if (nWet < 3)
         {
-            FatalErrorInFunction
-                << "There are fewer than three points"
-                << " on the dry part of the face."
-                << abort(FatalError);
+            // The face is dry
+            isAllWet_ = false;
+            isAllDry_ = true;
+
+            wetSubface_.clear();
         }
         else
         {
@@ -342,7 +344,6 @@ Foam::ImmersedFace<Distance>::ImmersedFace
         absTol = minEdgeLength*immersedPoly::tolerance_();
     }
 
-
     // Check if all points are wet or dry, using absolute tolerance
     if (max(depth) < absTol)
     {

src/immersedBoundary/immersedBoundary/include/immersedBoundaryCourantNo.H

@@ -36,7 +36,7 @@ scalar velMag = 0.0;
 
 if (mesh.nInternalFaces())
 {
-    surfaceScalarField magPhi = pos(sGamma - 0.5)*mag(phi);
+    surfaceScalarField magPhi = sGamma*mag(phi);
 
     surfaceScalarField SfUfbyDelta =
         mesh.surfaceInterpolation::deltaCoeffs()*magPhi;
@@ -50,7 +50,7 @@ if (mesh.nInternalFaces())
     velMag = max(magPhi/mesh.magSf()).value();
 }
 
-Info<< "Courant Number mean: " << meanCoNum
+Info<< "Immersed boundary Courant Number mean: " << meanCoNum
     << " max: " << CoNum
     << " velocity magnitude: " << velMag
     << endl;

src/immersedBoundary/immersedBoundary/include/updateIbMasks.H

@@ -8,51 +8,14 @@
     gamma.correctBoundaryConditions();
     sGamma.correctBoundaryConditions();
 
+
+    // Reset gamma and sGamma
+    gamma == dimensionedScalar("one", dimless, 1);
+    sGamma == dimensionedScalar("one", dimless, 1);
+
     scalarField& gammaIn = gamma.internalField();
-
-    gammaIn = mesh.V()/mesh.cellVolumes();
-
     scalarField& sGammaIn = sGamma.internalField();
 
-    const surfaceScalarField& magSf = mesh.magSf();
-    const scalarField magFaceAreas = mag(mesh.faceAreas());
-
-    sGammaIn =  magSf.internalField()/
-        scalarField::subField(magFaceAreas, mesh.nInternalFaces());
-
-    const scalar maxGamma = max(gammaIn);
-    const scalar maxSGamma = min(sGammaIn);
-    if (maxGamma > 1 || maxSGamma > 1)
-    {
-        WarningIn
-        (
-            "updateIbMasks.H"
-        )
-            << "Unbounded gamma: max(gamma) = "
-            << maxGamma << ", max(sGamma) = " << maxSGamma << endl;
-    }
-
-    forAll (mesh.boundary(), patchI)
-    {
-        if (!isA<immersedBoundaryFvPatch>(mesh.boundary()[patchI]))
-        {
-            sGamma.boundaryField()[patchI] =
-                magSf.boundaryField()[patchI]/
-                mesh.boundary()[patchI].patchSlice(magFaceAreas);
-
-            gamma.boundaryField()[patchI] =
-                sGamma.boundaryField()[patchI];
-        }
-        else
-        {
-            sGamma.boundaryField()[patchI] =
-                scalarField(mesh.boundary()[patchI].size(), 1);
-
-            gamma.boundaryField()[patchI] =
-                scalarField(mesh.boundary()[patchI].size(), 1);
-        }
-    }
-
     // Adjust the mask in dead faces and dead cells
     forAll (mesh.boundary(), patchI)
     {

src/immersedBoundary/immersedBoundaryDynamicMesh/refineImmersedBoundaryMesh/refineImmersedBoundaryMesh.C

@@ -321,6 +321,7 @@ Foam::refineImmersedBoundaryMesh::refinementCells
         case IB_CELL_CELLS:
         {
             addIbCellCells(refCellSet);
+            [[fallthrough]];
         }
 
         case IB_CELLS:

src/immersedBoundary/immersedBoundaryTurbulence/CMakeLists.txt

@@ -0,0 +1,34 @@
+# --------------------------------------------------------------------------
+#   ========                 |
+#   \      /  F ield         | foam-extend: Open Source CFD
+#    \    /   O peration     | Version:     4.1
+#     \  /    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/>.
+#
+# Description
+#     CMakeLists.txt file for libraries and applications
+#
+# Author
+#     Hrvoje Jasak, Wikki Ltd, 2019. All rights reserved
+#
+#
+# --------------------------------------------------------------------------
+
+add_subdirectory(incompressible)
+add_subdirectory(compressible)

src/immersedBoundary/immersedBoundaryTurbulence/compressible/Make/files

@@ -0,0 +1,8 @@
+wallFunctions/immersedBoundaryKqRWallFunction/immersedBoundaryKqRWallFunctionFvPatchFields.C
+wallFunctions/immersedBoundaryEpsilonWallFunction/immersedBoundaryEpsilonWallFunctionFvPatchScalarField.C
+wallFunctions/immersedBoundaryOmegaWallFunctions/immersedBoundaryOmegaWallFunctionFvPatchScalarField.C
+
+wallFunctions/immersedBoundaryMutWallFunction/immersedBoundaryMutWallFunctionFvPatchScalarField.C
+wallFunctions/immersedBoundaryAlphatWallFunction/immersedBoundaryAlphatWallFunctionFvPatchScalarField.C
+
+LIB = $(FOAM_LIBBIN)/libcompressibleImmersedBoundaryTurbulence

src/immersedBoundary/immersedBoundaryTurbulence/compressible/Make/options

@@ -0,0 +1,17 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/immersedBoundary/immersedBoundary/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels \
+    -I$(LIB_SRC)/turbulenceModels/compressible/RAS/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude
+
+LIB_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lsampling \
+    -limmersedBoundary \
+    -lcompressibleTurbulenceModel \
+    -lbasicThermophysicalModels \
+    -lspecie

src/immersedBoundary/immersedBoundaryTurbulence/compressible/wallFunctions/immersedBoundaryAlphatWallFunction/immersedBoundaryAlphatWallFunctionFvPatchScalarField.C

@@ -0,0 +1,229 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "immersedBoundaryAlphatWallFunctionFvPatchScalarField.H"
+#include "RASModel.H"
+#include "fvPatchFieldMapper.H"
+#include "standAlonePatch.H"
+#include "surfaceWriter.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace compressible
+{
+namespace RASModels
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+immersedBoundaryAlphatWallFunctionFvPatchScalarField::
+immersedBoundaryAlphatWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    alphatWallFunctionFvPatchScalarField(p, iF),
+    immersedBoundaryFieldBase<scalar>(p, true, 1e-6)
+{}
+
+
+immersedBoundaryAlphatWallFunctionFvPatchScalarField::
+immersedBoundaryAlphatWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    alphatWallFunctionFvPatchScalarField(p, iF),   // Do not read mixed data
+    immersedBoundaryFieldBase<scalar>(p, true, 1e-6)
+{
+    this->readPatchType(dict);
+
+    if (!isType<immersedBoundaryFvPatch>(p))
+    {
+        FatalIOErrorInFunction(dict)
+            << "\n    patch type '" << p.type()
+            << "' not constraint type '" << typeName << "'"
+            << "\n    for patch " << p.name()
+            << " of field " << this->dimensionedInternalField().name()
+            << " in file " << this->dimensionedInternalField().objectPath()
+            << exit(FatalIOError);
+    }
+
+    scalarField::operator=(this->patchInternalField());
+}
+
+
+immersedBoundaryAlphatWallFunctionFvPatchScalarField::
+immersedBoundaryAlphatWallFunctionFvPatchScalarField
+(
+    const immersedBoundaryAlphatWallFunctionFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    alphatWallFunctionFvPatchScalarField(p, iF), // Do not map mixed data.  Set patchType later
+    immersedBoundaryFieldBase<scalar>(p, true, 1e-6)
+{
+    // Note: NO MAPPING.  Fields are created on the immersed boundary
+    // HJ, 12/Apr/2012
+    if (!isType<immersedBoundaryFvPatch>(p))
+    {
+        FatalErrorInFunction
+            << "\n    patch type '" << p.type()
+            << "' not constraint type '" << typeName << "'"
+            << "\n    for patch " << p.name()
+            << " of field " << this->dimensionedInternalField().name()
+            << " in file " << this->dimensionedInternalField().objectPath()
+            << exit(FatalIOError);
+    }
+
+    this->setPatchType(ptf);
+
+    // On creation of the mapped field, the internal field is dummy and
+    // cannot be used.  Initialise the value to avoid errors
+    // HJ, 1/Dec/2017
+    scalarField::operator=(scalar(0));
+}
+
+
+immersedBoundaryAlphatWallFunctionFvPatchScalarField::
+immersedBoundaryAlphatWallFunctionFvPatchScalarField
+(
+    const immersedBoundaryAlphatWallFunctionFvPatchScalarField& ptf
+)
+:
+    alphatWallFunctionFvPatchScalarField(ptf),
+    immersedBoundaryFieldBase<scalar>(ptf.ibPatch(), true, 1e-6)
+{
+    this->setPatchType(ptf);
+}
+
+
+immersedBoundaryAlphatWallFunctionFvPatchScalarField::
+immersedBoundaryAlphatWallFunctionFvPatchScalarField
+(
+    const immersedBoundaryAlphatWallFunctionFvPatchScalarField& ptf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    alphatWallFunctionFvPatchScalarField(ptf, iF),
+    immersedBoundaryFieldBase<scalar>(ptf.ibPatch(), true, 1e-6)
+{
+    this->setPatchType(ptf);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void immersedBoundaryAlphatWallFunctionFvPatchScalarField::autoMap
+(
+    const fvPatchFieldMapper&
+)
+{
+    scalarField::operator=(this->patchInternalField());
+}
+
+
+void immersedBoundaryAlphatWallFunctionFvPatchScalarField::rmap
+(
+    const fvPatchScalarField& ptf,
+    const labelList&
+)
+{}
+
+
+void immersedBoundaryAlphatWallFunctionFvPatchScalarField::updateOnMotion()
+{
+    if (size() != ibPatch().size())
+    {
+        // Use internal values, resizing the file if needed
+        scalarField::operator=(this->patchInternalField());
+    }
+}
+
+
+void immersedBoundaryAlphatWallFunctionFvPatchScalarField::evaluate
+(
+    const Pstream::commsTypes commsType
+)
+{
+    // Resize fields
+    if (size() != patch().size())
+    {
+        Info<< "Resizing immersedBoundaryAlphatWallFunction in evaluate"
+            << endl;
+
+        scalarField::operator=(patchInternalField());
+    }
+
+    // Get non-constant reference to internal field
+    scalarField& intField = const_cast<scalarField&>(this->internalField());
+
+    // Set dead value
+    this->setDeadValues(intField);
+
+    alphatWallFunctionFvPatchScalarField::evaluate(commsType);
+}
+
+
+void immersedBoundaryAlphatWallFunctionFvPatchScalarField::write
+(
+    Ostream& os
+) const
+{
+    fvPatchScalarField::write(os);
+    writeLocalEntries(os);
+
+    // The value entry needs to be written with zero size
+    scalarField::null().writeEntry("value", os);
+    // this->writeEntry("value", os);
+
+    writeField(*this);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeField
+(
+    fvPatchScalarField,
+    immersedBoundaryAlphatWallFunctionFvPatchScalarField
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace compressible
+} // End namespace Foam
+
+// ************************************************************************* //

src/immersedBoundary/immersedBoundaryTurbulence/compressible/wallFunctions/immersedBoundaryAlphatWallFunction/immersedBoundaryAlphatWallFunctionFvPatchScalarField.H

@@ -0,0 +1,186 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::compressible::RASModels::immersedBoundaryAlphatWallFunctionFvPatchScalarField
+
+Description
+    Boundary condition for turbulent (kinematic) viscosity when using wall
+    functions on immersed boundary patches
+    - replicates OpenFOAM v1.5 (and earlier) behaviour
+
+SourceFiles
+    immersedBoundaryAlphatWallFunctionFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef immersedBoundaryAlphatWallFunctionFvPatchScalarField_H
+#define immersedBoundaryAlphatWallFunctionFvPatchScalarField_H
+
+#include "alphatWallFunctionFvPatchScalarField.H"
+#include "immersedBoundaryFieldBase.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace compressible
+{
+namespace RASModels
+{
+
+/*---------------------------------------------------------------------------*\
+     Class immersedBoundaryAlphatWallFunctionFvPatchScalarField Declaration
+\*---------------------------------------------------------------------------*/
+
+class immersedBoundaryAlphatWallFunctionFvPatchScalarField
+:
+    public alphatWallFunctionFvPatchScalarField,
+    public immersedBoundaryFieldBase<scalar>
+{
+public:
+
+    //- Runtime type information
+    TypeName("immersedBoundaryAlphatWallFunction");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        immersedBoundaryAlphatWallFunctionFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        immersedBoundaryAlphatWallFunctionFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        //  immersedBoundaryAlphatWallFunctionFvPatchScalarField
+        //  onto a new patch
+        immersedBoundaryAlphatWallFunctionFvPatchScalarField
+        (
+            const immersedBoundaryAlphatWallFunctionFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        immersedBoundaryAlphatWallFunctionFvPatchScalarField
+        (
+            const immersedBoundaryAlphatWallFunctionFvPatchScalarField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new immersedBoundaryAlphatWallFunctionFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        immersedBoundaryAlphatWallFunctionFvPatchScalarField
+        (
+            const immersedBoundaryAlphatWallFunctionFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new immersedBoundaryAlphatWallFunctionFvPatchScalarField
+                (
+                    *this,
+                    iF
+                )
+            );
+        }
+
+
+    //- Destructor
+    virtual ~immersedBoundaryAlphatWallFunctionFvPatchScalarField()
+    {}
+
+
+    // Member functions
+
+        // Mapping functions
+
+            //- Map (and resize as needed) from self given a mapping object
+            virtual void autoMap
+            (
+                const fvPatchFieldMapper&
+            );
+
+            //- Reverse map the given fvPatchField onto this fvPatchField
+            virtual void rmap
+            (
+                const fvPatchScalarField&,
+                const labelList&
+            );
+
+            //- Update on mesh motion
+            virtual void updateOnMotion();
+
+
+        // Evaluation functions
+
+            //- Evaluate the patchField
+            virtual void evaluate
+            (
+                const Pstream::commsTypes  = Pstream::blocking
+            );
+
+
+        // I-O
+
+            //- Write
+            virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace compressible
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/immersedBoundary/immersedBoundaryTurbulence/compressible/wallFunctions/immersedBoundaryEpsilonWallFunction/immersedBoundaryEpsilonWallFunctionFvPatchScalarField.C

@@ -0,0 +1,296 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.1
+    \\  /    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 "immersedBoundaryEpsilonWallFunctionFvPatchScalarField.H"
+#include "RASModel.H"
+#include "fvPatchFieldMapper.H"
+#include "standAlonePatch.H"
+#include "surfaceWriter.H"
+#include "volFields.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace compressible
+{
+namespace RASModels
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+immersedBoundaryEpsilonWallFunctionFvPatchScalarField::
+immersedBoundaryEpsilonWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    epsilonWallFunctionFvPatchScalarField(p, iF),
+    immersedBoundaryFieldBase<scalar>(p, true, 0.09)
+{}
+
+
+immersedBoundaryEpsilonWallFunctionFvPatchScalarField::
+immersedBoundaryEpsilonWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    epsilonWallFunctionFvPatchScalarField(p, iF),   // Do not read mixed data
+    immersedBoundaryFieldBase<scalar>
+    (
+        p,
+        Switch(dict.lookup("setDeadValue")),
+        readScalar(dict.lookup("deadValue"))
+    )
+{
+    // Since patch does not read a dictionary, the patch type needs to be read
+    // manually.  HJ, 6/Sep/2018
+    this->readPatchType(dict);
+
+    if (!isType<immersedBoundaryFvPatch>(p))
+    {
+        FatalIOErrorInFunction(dict)
+            << "\n    patch type '" << p.type()
+            << "' not constraint type '" << typeName << "'"
+            << "\n    for patch " << p.name()
+            << " of field " << this->dimensionedInternalField().name()
+            << " in file " << this->dimensionedInternalField().objectPath()
+            << exit(FatalIOError);
+    }
+
+    scalarField::operator=(this->patchInternalField());
+}
+
+
+immersedBoundaryEpsilonWallFunctionFvPatchScalarField::
+immersedBoundaryEpsilonWallFunctionFvPatchScalarField
+(
+    const immersedBoundaryEpsilonWallFunctionFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    epsilonWallFunctionFvPatchScalarField(p, iF), // Do not map mixed data.  Set patchType later
+    immersedBoundaryFieldBase<scalar>
+    (
+        p,
+        ptf.setDeadValue(),
+        ptf.deadValue()
+    )
+{
+    // Note: NO MAPPING.  Fields are created on the immersed boundary
+    // HJ, 12/Apr/2012
+    if (!isType<immersedBoundaryFvPatch>(p))
+    {
+        FatalErrorInFunction
+            << "\n    patch type '" << p.type()
+            << "' not constraint type '" << typeName << "'"
+            << "\n    for patch " << p.name()
+            << " of field " << this->dimensionedInternalField().name()
+            << " in file " << this->dimensionedInternalField().objectPath()
+            << exit(FatalIOError);
+    }
+
+    // Copy the patch type since mixed data was not mapped
+    this->setPatchType(ptf);
+
+    // On creation of the mapped field, the internal field is dummy and
+    // cannot be used.  Initialise the value to avoid errors
+    // HJ, 1/Dec/2017
+    scalarField::operator=(SMALL);
+}
+
+
+immersedBoundaryEpsilonWallFunctionFvPatchScalarField::
+immersedBoundaryEpsilonWallFunctionFvPatchScalarField
+(
+    const immersedBoundaryEpsilonWallFunctionFvPatchScalarField& ptf
+)
+:
+    epsilonWallFunctionFvPatchScalarField(ptf),
+    immersedBoundaryFieldBase<scalar>
+    (
+        ptf.ibPatch(),
+        ptf.setDeadValue(),
+        ptf.deadValue()
+    )
+{}
+
+
+immersedBoundaryEpsilonWallFunctionFvPatchScalarField::
+immersedBoundaryEpsilonWallFunctionFvPatchScalarField
+(
+    const immersedBoundaryEpsilonWallFunctionFvPatchScalarField& ptf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    epsilonWallFunctionFvPatchScalarField(ptf, iF),
+    immersedBoundaryFieldBase<scalar>
+    (
+        ptf.ibPatch(),
+        ptf.setDeadValue(),
+        ptf.deadValue()
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void immersedBoundaryEpsilonWallFunctionFvPatchScalarField::autoMap
+(
+    const fvPatchFieldMapper&
+)
+{
+    scalarField::operator=(this->patchInternalField());
+
+    // Resize refValue as well.  HJ, 10/Jul/2018
+    refValue() = this->patchInternalField();
+}
+
+
+void immersedBoundaryEpsilonWallFunctionFvPatchScalarField::rmap
+(
+    const fvPatchScalarField& ptf,
+    const labelList&
+)
+{}
+
+
+void immersedBoundaryEpsilonWallFunctionFvPatchScalarField::updateOnMotion()
+{
+    if (size() != ibPatch().size())
+    {
+        // Use internal values, resizing the file if needed
+        scalarField::operator=(this->patchInternalField());
+
+        // Resize refValue as well.  HJ, 10/Jul/2018
+        refValue() = this->patchInternalField();
+    }
+}
+
+
+void immersedBoundaryEpsilonWallFunctionFvPatchScalarField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    // Resize fields
+    if (size() != patch().size())
+    {
+        Info<< "Resizing immersedBoundaryEpsilonWallFunction in evaluate: "
+            << "patch: " << patch().size() << " field: " << size()
+            << endl;
+
+        *this == patchInternalField();
+        refValue() = patchInternalField();
+    }
+
+    // If G field is present, execute evaluation
+    // Remove the warning from the IB patch
+    // HJ, 20/May/2018
+    if (db().foundObject<volScalarField>(GName()))
+    {
+        epsilonWallFunctionFvPatchScalarField::updateCoeffs();
+    }
+}
+
+
+void immersedBoundaryEpsilonWallFunctionFvPatchScalarField::evaluate
+(
+    const Pstream::commsTypes commsType
+)
+{
+    // Resize fields
+    if (size() != patch().size())
+    {
+        Info<< "Resizing immersedBoundaryEpsilonWallFunction in evaluate"
+            << endl;
+
+        *this == patchInternalField();
+        refValue() = patchInternalField();
+    }
+
+    // Get non-constant reference to internal field
+    scalarField& intField = const_cast<scalarField&>(this->internalField());
+
+    // Set dead value
+    this->setDeadValues(intField);
+
+    epsilonWallFunctionFvPatchScalarField::evaluate(commsType);
+}
+
+
+void immersedBoundaryEpsilonWallFunctionFvPatchScalarField::manipulateMatrix
+(
+    fvScalarMatrix& matrix
+)
+{
+    setDeadValues(matrix);
+
+    epsilonWallFunctionFvPatchScalarField::manipulateMatrix(matrix);
+}
+
+
+void immersedBoundaryEpsilonWallFunctionFvPatchScalarField::write
+(
+    Ostream& os
+) const
+{
+    fvPatchScalarField::write(os);
+    writeLocalEntries(os);
+    this->writeDeadData(os);
+
+    // The value entry needs to be written with zero size
+    scalarField::null().writeEntry("value", os);
+    // this->writeEntry("value", os);
+
+    writeField(*this);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeField
+(
+    fvPatchScalarField,
+    immersedBoundaryEpsilonWallFunctionFvPatchScalarField
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace compressible
+} // End namespace Foam
+
+// ************************************************************************* //

src/immersedBoundary/immersedBoundaryTurbulence/compressible/wallFunctions/immersedBoundaryEpsilonWallFunction/immersedBoundaryEpsilonWallFunctionFvPatchScalarField.H

@@ -0,0 +1,194 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.1
+    \\  /    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::compressible::RASModels::immersedBoundaryEpsilonWallFunctionFvPatchScalarField
+
+Description
+    Boundary condition for epsilon when using wall functions
+    on immersed boundary patches
+
+Author
+    Hrvoje Jasak, Wikki Ltd.  All rights reserved
+
+SourceFiles
+    immersedBoundaryEpsilonWallFunctionFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef immersedBoundaryEpsilonWallFunctionFvPatchScalarField_H
+#define immersedBoundaryEpsilonWallFunctionFvPatchScalarField_H
+
+#include "epsilonWallFunctionFvPatchScalarField.H"
+#include "immersedBoundaryFieldBase.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace compressible
+{
+namespace RASModels
+{
+
+/*---------------------------------------------------------------------------*\
+    Class immersedBoundaryEpsilonWallFunctionFvPatchScalarField Declaration
+\*---------------------------------------------------------------------------*/
+
+class immersedBoundaryEpsilonWallFunctionFvPatchScalarField
+:
+    public epsilonWallFunctionFvPatchScalarField,
+    public immersedBoundaryFieldBase<scalar>
+{
+public:
+
+    //- Runtime type information
+    TypeName("compressible::immersedBoundaryEpsilonWallFunction");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        immersedBoundaryEpsilonWallFunctionFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        immersedBoundaryEpsilonWallFunctionFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        //  immersedBoundaryEpsilonWallFunctionFvPatchScalarField
+        //  onto a new patch
+        immersedBoundaryEpsilonWallFunctionFvPatchScalarField
+        (
+            const immersedBoundaryEpsilonWallFunctionFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        immersedBoundaryEpsilonWallFunctionFvPatchScalarField
+        (
+            const immersedBoundaryEpsilonWallFunctionFvPatchScalarField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new immersedBoundaryEpsilonWallFunctionFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        immersedBoundaryEpsilonWallFunctionFvPatchScalarField
+        (
+            const immersedBoundaryEpsilonWallFunctionFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new immersedBoundaryEpsilonWallFunctionFvPatchScalarField
+                (
+                    *this,
+                    iF
+                )
+            );
+        }
+
+
+    //- Destructor
+    virtual ~immersedBoundaryEpsilonWallFunctionFvPatchScalarField()
+    {}
+
+
+    // Member functions
+
+        // Mapping functions
+
+            //- Map (and resize as needed) from self given a mapping object
+            virtual void autoMap
+            (
+                const fvPatchFieldMapper&
+            );
+
+            //- Reverse map the given fvPatchField onto this fvPatchField
+            virtual void rmap
+            (
+                const fvPatchScalarField&,
+                const labelList&
+            );
+
+            //- Update on mesh motion
+            virtual void updateOnMotion();
+
+
+        // Evaluation functions
+
+            //- Update the coefficients associated with the patch field
+            virtual void updateCoeffs();
+
+            //- Evaluate the patchField
+            virtual void evaluate
+            (
+                const Pstream::commsTypes  = Pstream::blocking
+            );
+
+            //- Manipulate a matrix
+            virtual void manipulateMatrix(fvScalarMatrix& matrix);
+
+
+        // I-O
+
+            //- Write
+            void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace compressible
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/immersedBoundary/immersedBoundaryTurbulence/compressible/wallFunctions/immersedBoundaryKqRWallFunction/immersedBoundaryKqRWallFunctionFvPatchField.C

@@ -0,0 +1,253 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "immersedBoundaryKqRWallFunctionFvPatchField.H"
+#include "fvPatchFieldMapper.H"
+#include "standAlonePatch.H"
+#include "surfaceWriter.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace compressible
+{
+namespace RASModels
+{
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+void immersedBoundaryKqRWallFunctionFvPatchField<Type>::checkType()
+{
+    if (!this->patch().isWall())
+    {
+        FatalErrorIn("immersedBoundaryKqRWallFunctionFvPatchField::checkType()")
+            << "Invalid wall function specification" << nl
+            << "    Patch type for patch " << this->patch().name()
+            << " must be wall" << nl
+            << "    Current patch type is " << this->patch().type()
+            << nl << endl
+            << abort(FatalError);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+immersedBoundaryKqRWallFunctionFvPatchField<Type>::
+immersedBoundaryKqRWallFunctionFvPatchField
+(
+    const fvPatch& p,
+    const DimensionedField<Type, volMesh>& iF
+)
+:
+    kqRWallFunctionFvPatchField<Type>(p, iF),
+    immersedBoundaryFieldBase<Type>(p, true, pTraits<Type>::one*SMALL)
+{
+    this->checkType();
+}
+
+
+template<class Type>
+immersedBoundaryKqRWallFunctionFvPatchField<Type>::
+immersedBoundaryKqRWallFunctionFvPatchField
+(
+    const fvPatch& p,
+    const DimensionedField<Type, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    kqRWallFunctionFvPatchField<Type>(p, iF),   // Do not read mixed data
+    immersedBoundaryFieldBase<Type>
+    (
+        p,
+        Switch(dict.lookup("setDeadValue")),
+        pTraits<Type>(dict.lookup("deadValue"))
+    )
+{
+    this->readPatchType(dict);
+    this->checkType();
+
+    Field<Type>::operator=(this->patchInternalField());
+}
+
+
+template<class Type>
+immersedBoundaryKqRWallFunctionFvPatchField<Type>::
+immersedBoundaryKqRWallFunctionFvPatchField
+(
+    const immersedBoundaryKqRWallFunctionFvPatchField& ptf,
+    const fvPatch& p,
+    const DimensionedField<Type, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    kqRWallFunctionFvPatchField<Type>(p, iF), // Do not map mixed data.  Set patchType later
+    immersedBoundaryFieldBase<Type>
+    (
+        p,
+        ptf.setDeadValue(),
+        ptf.deadValue()
+    )
+{
+    this->setPatchType(ptf);
+    this->checkType();
+
+    // cannot be used.  Initialise the value to avoid errors
+    // HJ, 1/Dec/2017
+    Field<Type>::operator=(pTraits<Type>::zero);
+}
+
+
+template<class Type>
+immersedBoundaryKqRWallFunctionFvPatchField<Type>::
+immersedBoundaryKqRWallFunctionFvPatchField
+(
+    const immersedBoundaryKqRWallFunctionFvPatchField& ptf
+)
+:
+    kqRWallFunctionFvPatchField<Type>(ptf),
+    immersedBoundaryFieldBase<Type>
+    (
+        ptf.ibPatch(),
+        ptf.setDeadValue(),
+        ptf.deadValue()
+    )
+{
+    this->setPatchType(ptf);
+    this->checkType();
+}
+
+
+template<class Type>
+immersedBoundaryKqRWallFunctionFvPatchField<Type>::
+immersedBoundaryKqRWallFunctionFvPatchField
+(
+    const immersedBoundaryKqRWallFunctionFvPatchField& ptf,
+    const DimensionedField<Type, volMesh>& iF
+)
+:
+    kqRWallFunctionFvPatchField<Type>(ptf, iF),
+    immersedBoundaryFieldBase<Type>
+    (
+        ptf.ibPatch(),
+        ptf.setDeadValue(),
+        ptf.deadValue()
+    )
+{
+    this->setPatchType(ptf);
+    this->checkType();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+void immersedBoundaryKqRWallFunctionFvPatchField<Type>::autoMap
+(
+    const fvPatchFieldMapper&
+)
+{
+    Field<Type>::operator=(this->patchInternalField());
+}
+
+
+template<class Type>
+void immersedBoundaryKqRWallFunctionFvPatchField<Type>::rmap
+(
+    const fvPatchField<Type>& ptf,
+    const labelList&
+)
+{}
+
+
+template<class Type>
+void immersedBoundaryKqRWallFunctionFvPatchField<Type>::updateOnMotion()
+{
+    if (this->size() != this->ibPatch().size())
+    {
+        // Use internal values, resizing the file if needed
+        Field<Type>::operator=(this->patchInternalField());
+    }
+}
+
+
+template<class Type>
+void immersedBoundaryKqRWallFunctionFvPatchField<Type>::evaluate
+(
+    const Pstream::commsTypes commsType
+)
+{
+    // Resize on evaluation if needed
+    if (this->size() != this->ibPatch().size())
+    {
+        // Use internal values, resizing the file if needed
+        Field<Type>::operator=(this->patchInternalField());
+    }
+
+    // Get non-constant reference to internal field
+    Field<Type>& intField = const_cast<Field<Type>&>(this->internalField());
+
+    // Set dead value
+    this->setDeadValues(intField);
+
+    kqRWallFunctionFvPatchField<Type>::evaluate(commsType);
+}
+
+
+template<class Type>
+void immersedBoundaryKqRWallFunctionFvPatchField<Type>::manipulateMatrix
+(
+    fvMatrix<Type>& matrix
+)
+{
+    this->setDeadValues(matrix);
+}
+
+
+template<class Type>
+void immersedBoundaryKqRWallFunctionFvPatchField<Type>::write(Ostream& os) const
+{
+    fvPatchField<Type>::write(os);
+    this->writeDeadData(os);
+
+    // The value entry needs to be written with zero size
+    Field<Type>::null().writeEntry("value", os);
+    // this->writeEntry("value", os);
+
+    this->writeField(*this);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace compressible
+} // End namespace Foam
+
+// ************************************************************************* //

src/immersedBoundary/immersedBoundaryTurbulence/compressible/wallFunctions/immersedBoundaryKqRWallFunction/immersedBoundaryKqRWallFunctionFvPatchField.H

@@ -0,0 +1,197 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::compressible::RASModels::immersedBoundaryKqRWallFunctionFvPatchField
+
+Description
+    Boundary condition for turbulence k, Q, and R when using wall functions.
+    Simply acts as a zero gradient condition.
+
+SourceFiles
+    immersedBoundaryKqRWallFunctionFvPatchField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef immersedBoundaryKqRWallFunctionFvPatchField_H
+#define immersedBoundaryKqRWallFunctionFvPatchField_H
+
+#include "kqRWallFunctionFvPatchFields.H"
+#include "immersedBoundaryFieldBase.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace compressible
+{
+namespace RASModels
+{
+
+/*---------------------------------------------------------------------------*\
+         Class immersedBoundaryKqRWallFunctionFvPatchField Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class immersedBoundaryKqRWallFunctionFvPatchField
+:
+    public kqRWallFunctionFvPatchField<Type>,
+    public immersedBoundaryFieldBase<Type>
+{
+    // Private member functions
+
+        //- Check the type of the patch
+        void checkType();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("compressible::immersedBoundaryKqRWallFunction");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        immersedBoundaryKqRWallFunctionFvPatchField
+        (
+            const fvPatch&,
+            const DimensionedField<Type, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        immersedBoundaryKqRWallFunctionFvPatchField
+        (
+            const fvPatch&,
+            const DimensionedField<Type, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        // immersedBoundaryKqRWallFunctionFvPatchField
+        //  onto a new patch
+        immersedBoundaryKqRWallFunctionFvPatchField
+        (
+            const immersedBoundaryKqRWallFunctionFvPatchField&,
+            const fvPatch&,
+            const DimensionedField<Type, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        immersedBoundaryKqRWallFunctionFvPatchField
+        (
+            const immersedBoundaryKqRWallFunctionFvPatchField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchField<Type> > clone() const
+        {
+            return tmp<fvPatchField<Type> >
+            (
+                new immersedBoundaryKqRWallFunctionFvPatchField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        immersedBoundaryKqRWallFunctionFvPatchField
+        (
+            const immersedBoundaryKqRWallFunctionFvPatchField&,
+            const DimensionedField<Type, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchField<Type> > clone
+        (
+            const DimensionedField<Type, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchField<Type> >
+            (
+                new immersedBoundaryKqRWallFunctionFvPatchField(*this, iF)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~immersedBoundaryKqRWallFunctionFvPatchField()
+    {}
+
+
+    // Member functions
+
+        // Mapping functions
+
+            //- Map (and resize as needed) from self given a mapping object
+            virtual void autoMap
+            (
+                const fvPatchFieldMapper&
+            );
+
+            //- Reverse map the given fvPatchField onto this fvPatchField
+            virtual void rmap
+            (
+                const fvPatchField<Type>&,
+                const labelList&
+            );
+
+            //- Update on mesh motion
+            virtual void updateOnMotion();
+
+
+        // Evaluation functions
+
+            //- Evaluate the patchField
+            virtual void evaluate
+            (
+                const Pstream::commsTypes commsType = Pstream::blocking
+            );
+
+            //- Manipulate a matrix
+            virtual void manipulateMatrix(fvMatrix<Type>& matrix);
+
+
+        // I-O
+
+            //- Write
+            void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace compressible
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "immersedBoundaryKqRWallFunctionFvPatchField.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/immersedBoundary/immersedBoundaryTurbulence/compressible/wallFunctions/immersedBoundaryKqRWallFunction/immersedBoundaryKqRWallFunctionFvPatchFields.C

@@ -0,0 +1,49 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "immersedBoundaryKqRWallFunctionFvPatchFields.H"
+#include "volFields.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace compressible
+{
+namespace RASModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+makePatchFields(immersedBoundaryKqRWallFunction);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace compressible
+} // End namespace Foam
+
+// ************************************************************************* //

src/immersedBoundary/immersedBoundaryTurbulence/compressible/wallFunctions/immersedBoundaryKqRWallFunction/immersedBoundaryKqRWallFunctionFvPatchFields.H

@@ -0,0 +1,55 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef immersedBoundaryKqRWallFunctionFvPatchFields_H
+#define immersedBoundaryKqRWallFunctionFvPatchFields_H
+
+#include "immersedBoundaryKqRWallFunctionFvPatchField.H"
+#include "fieldTypes.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace compressible
+{
+namespace RASModels
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeFieldTypedefs(immersedBoundaryKqRWallFunction)
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace compressible
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/immersedBoundary/immersedBoundaryTurbulence/compressible/wallFunctions/immersedBoundaryMutWallFunction/immersedBoundaryMutWallFunctionFvPatchScalarField.C

@@ -0,0 +1,226 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "immersedBoundaryMutWallFunctionFvPatchScalarField.H"
+#include "RASModel.H"
+#include "fvPatchFieldMapper.H"
+#include "standAlonePatch.H"
+#include "surfaceWriter.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace compressible
+{
+namespace RASModels
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+immersedBoundaryMutWallFunctionFvPatchScalarField::
+immersedBoundaryMutWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    mutkWallFunctionFvPatchScalarField(p, iF),
+    immersedBoundaryFieldBase<scalar>(p, true, 1e-6)
+{}
+
+
+immersedBoundaryMutWallFunctionFvPatchScalarField::
+immersedBoundaryMutWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    mutkWallFunctionFvPatchScalarField(p, iF),   // Do not read mixed data
+    immersedBoundaryFieldBase<scalar>(p, true, 1e-6)
+{
+    this->readPatchType(dict);
+
+    if (!isType<immersedBoundaryFvPatch>(p))
+    {
+        FatalIOErrorInFunction(dict)
+            << "\n    patch type '" << p.type()
+            << "' not constraint type '" << typeName << "'"
+            << "\n    for patch " << p.name()
+            << " of field " << this->dimensionedInternalField().name()
+            << " in file " << this->dimensionedInternalField().objectPath()
+            << exit(FatalIOError);
+    }
+
+    scalarField::operator=(this->patchInternalField());
+}
+
+
+immersedBoundaryMutWallFunctionFvPatchScalarField::
+immersedBoundaryMutWallFunctionFvPatchScalarField
+(
+    const immersedBoundaryMutWallFunctionFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    mutkWallFunctionFvPatchScalarField(p, iF), // Do not map mixed data.  Set patchType later
+    immersedBoundaryFieldBase<scalar>(p, true, 1e-6)
+{
+    // Note: NO MAPPING.  Fields are created on the immersed boundary
+    // HJ, 12/Apr/2012
+    if (!isType<immersedBoundaryFvPatch>(p))
+    {
+        FatalErrorInFunction
+            << "\n    patch type '" << p.type()
+            << "' not constraint type '" << typeName << "'"
+            << "\n    for patch " << p.name()
+            << " of field " << this->dimensionedInternalField().name()
+            << " in file " << this->dimensionedInternalField().objectPath()
+            << exit(FatalIOError);
+    }
+
+    this->setPatchType(ptf);
+
+    // On creation of the mapped field, the internal field is dummy and
+    // cannot be used.  Initialise the value to avoid errors
+    // HJ, 1/Dec/2017
+    scalarField::operator=(scalar(0));
+}
+
+
+immersedBoundaryMutWallFunctionFvPatchScalarField::
+immersedBoundaryMutWallFunctionFvPatchScalarField
+(
+    const immersedBoundaryMutWallFunctionFvPatchScalarField& ptf
+)
+:
+    mutkWallFunctionFvPatchScalarField(ptf),
+    immersedBoundaryFieldBase<scalar>(ptf.ibPatch(), true, 1e-6)
+{
+    this->setPatchType(ptf);
+}
+
+
+immersedBoundaryMutWallFunctionFvPatchScalarField::
+immersedBoundaryMutWallFunctionFvPatchScalarField
+(
+    const immersedBoundaryMutWallFunctionFvPatchScalarField& ptf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    mutkWallFunctionFvPatchScalarField(ptf, iF),
+    immersedBoundaryFieldBase<scalar>(ptf.ibPatch(), true, 1e-6)
+{
+    this->setPatchType(ptf);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void immersedBoundaryMutWallFunctionFvPatchScalarField::autoMap
+(
+    const fvPatchFieldMapper&
+)
+{
+    scalarField::operator=(this->patchInternalField());
+}
+
+
+void immersedBoundaryMutWallFunctionFvPatchScalarField::rmap
+(
+    const fvPatchScalarField& ptf,
+    const labelList&
+)
+{}
+
+
+void immersedBoundaryMutWallFunctionFvPatchScalarField::updateOnMotion()
+{
+    if (size() != ibPatch().size())
+    {
+        // Use internal values, resizing the file if needed
+        scalarField::operator=(this->patchInternalField());
+    }
+}
+
+
+void immersedBoundaryMutWallFunctionFvPatchScalarField::evaluate
+(
+    const Pstream::commsTypes commsType
+)
+{
+    // Resize fields
+    if (size() != patch().size())
+    {
+        Info<< "Resizing immersedBoundaryMutWallFunction in evaluate"
+            << endl;
+
+        scalarField::operator=(patchInternalField());
+    }
+
+    // Get non-constant reference to internal field
+    scalarField& intField = const_cast<scalarField&>(this->internalField());
+
+    // Set dead value
+    this->setDeadValues(intField);
+
+    mutkWallFunctionFvPatchScalarField::evaluate(commsType);
+}
+
+
+void immersedBoundaryMutWallFunctionFvPatchScalarField::write(Ostream& os) const
+{
+    fvPatchScalarField::write(os);
+    writeLocalEntries(os);
+
+    // The value entry needs to be written with zero size
+    scalarField::null().writeEntry("value", os);
+    // this->writeEntry("value", os);
+
+    writeField(*this);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeField
+(
+    fvPatchScalarField,
+    immersedBoundaryMutWallFunctionFvPatchScalarField
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace compressible
+} // End namespace Foam
+
+// ************************************************************************* //

src/immersedBoundary/immersedBoundaryTurbulence/compressible/wallFunctions/immersedBoundaryMutWallFunction/immersedBoundaryMutWallFunctionFvPatchScalarField.H

@@ -0,0 +1,182 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::compressible::RASModels::immersedBoundaryMutWallFunctionFvPatchScalarField
+
+Description
+    Boundary condition for turbulent (kinematic) viscosity when using wall
+    functions on immersed boundary patches
+    - replicates OpenFOAM v1.5 (and earlier) behaviour
+
+SourceFiles
+    immersedBoundaryMutWallFunctionFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef immersedBoundaryMutWallFunctionFvPatchScalarField_H
+#define immersedBoundaryMutWallFunctionFvPatchScalarField_H
+
+#include "mutkWallFunctionFvPatchScalarField.H"
+#include "immersedBoundaryFieldBase.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace compressible
+{
+namespace RASModels
+{
+
+/*---------------------------------------------------------------------------*\
+      Class immersedBoundaryMutWallFunctionFvPatchScalarField Declaration
+\*---------------------------------------------------------------------------*/
+
+class immersedBoundaryMutWallFunctionFvPatchScalarField
+:
+    public mutkWallFunctionFvPatchScalarField,
+    public immersedBoundaryFieldBase<scalar>
+{
+public:
+
+    //- Runtime type information
+    TypeName("immersedBoundaryMutWallFunction");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        immersedBoundaryMutWallFunctionFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        immersedBoundaryMutWallFunctionFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        //  immersedBoundaryMutWallFunctionFvPatchScalarField
+        //  onto a new patch
+        immersedBoundaryMutWallFunctionFvPatchScalarField
+        (
+            const immersedBoundaryMutWallFunctionFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        immersedBoundaryMutWallFunctionFvPatchScalarField
+        (
+            const immersedBoundaryMutWallFunctionFvPatchScalarField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new immersedBoundaryMutWallFunctionFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        immersedBoundaryMutWallFunctionFvPatchScalarField
+        (
+            const immersedBoundaryMutWallFunctionFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new immersedBoundaryMutWallFunctionFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~immersedBoundaryMutWallFunctionFvPatchScalarField()
+    {}
+
+
+    // Member functions
+
+        // Mapping functions
+
+            //- Map (and resize as needed) from self given a mapping object
+            virtual void autoMap
+            (
+                const fvPatchFieldMapper&
+            );
+
+            //- Reverse map the given fvPatchField onto this fvPatchField
+            virtual void rmap
+            (
+                const fvPatchScalarField&,
+                const labelList&
+            );
+
+            //- Update on mesh motion
+            virtual void updateOnMotion();
+
+
+        // Evaluation functions
+
+            //- Evaluate the patchField
+            virtual void evaluate
+            (
+                const Pstream::commsTypes  = Pstream::blocking
+            );
+
+
+        // I-O
+
+            //- Write
+            virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace compressible
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/immersedBoundary/immersedBoundaryTurbulence/compressible/wallFunctions/immersedBoundaryOmegaWallFunctions/immersedBoundaryOmegaWallFunctionFvPatchScalarField.C

@@ -0,0 +1,294 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "immersedBoundaryOmegaWallFunctionFvPatchScalarField.H"
+#include "RASModel.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace compressible
+{
+namespace RASModels
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+immersedBoundaryOmegaWallFunctionFvPatchScalarField::
+immersedBoundaryOmegaWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    omegaWallFunctionFvPatchScalarField(p, iF),
+    immersedBoundaryFieldBase<scalar>(p, true, 90)
+{}
+
+
+immersedBoundaryOmegaWallFunctionFvPatchScalarField::
+immersedBoundaryOmegaWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    omegaWallFunctionFvPatchScalarField(p, iF),    // Do not read mixed data
+    immersedBoundaryFieldBase<scalar>
+    (
+        p,
+        Switch(dict.lookup("setDeadValue")),
+        readScalar(dict.lookup("deadValue"))
+    )
+{
+    this->readPatchType(dict);
+
+    if (!isType<immersedBoundaryFvPatch>(p))
+    {
+        FatalIOErrorInFunction(dict)
+            << "\n    patch type '" << p.type()
+            << "' not constraint type '" << typeName << "'"
+            << "\n    for patch " << p.name()
+            << " of field " << this->dimensionedInternalField().name()
+            << " in file " << this->dimensionedInternalField().objectPath()
+            << exit(FatalIOError);
+    }
+
+    scalarField::operator=(this->patchInternalField());
+}
+
+
+immersedBoundaryOmegaWallFunctionFvPatchScalarField::
+immersedBoundaryOmegaWallFunctionFvPatchScalarField
+(
+    const immersedBoundaryOmegaWallFunctionFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    omegaWallFunctionFvPatchScalarField(p, iF), // Do not map mixed data.  Set patchType later
+    immersedBoundaryFieldBase<scalar>
+    (
+        p,
+        ptf.setDeadValue(),
+        ptf.deadValue()
+    )
+{
+    // Note: NO MAPPING.  Fields are created on the immersed boundary
+    // HJ, 12/Apr/2012
+    if (!isType<immersedBoundaryFvPatch>(p))
+    {
+        FatalErrorInFunction
+            << "\n    patch type '" << p.type()
+            << "' not constraint type '" << typeName << "'"
+            << "\n    for patch " << p.name()
+            << " of field " << this->dimensionedInternalField().name()
+            << " in file " << this->dimensionedInternalField().objectPath()
+            << exit(FatalIOError);
+    }
+
+    this->setPatchType(ptf);
+
+    // On creation of the mapped field, the internal field is dummy and
+    // cannot be used.  Initialise the value to avoid errors
+    // HJ, 1/Dec/2017
+    scalarField::operator=(SMALL);
+}
+
+
+immersedBoundaryOmegaWallFunctionFvPatchScalarField::
+immersedBoundaryOmegaWallFunctionFvPatchScalarField
+(
+    const immersedBoundaryOmegaWallFunctionFvPatchScalarField& ptf
+)
+:
+    omegaWallFunctionFvPatchScalarField(ptf),
+    immersedBoundaryFieldBase<scalar>
+    (
+        ptf.ibPatch(),
+        ptf.setDeadValue(),
+        ptf.deadValue()
+    )
+{
+    this->setPatchType(ptf);
+}
+
+
+immersedBoundaryOmegaWallFunctionFvPatchScalarField::
+immersedBoundaryOmegaWallFunctionFvPatchScalarField
+(
+    const immersedBoundaryOmegaWallFunctionFvPatchScalarField& ptf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    omegaWallFunctionFvPatchScalarField(ptf, iF),
+    immersedBoundaryFieldBase<scalar>
+    (
+        ptf.ibPatch(),
+        ptf.setDeadValue(),
+        ptf.deadValue()
+    )
+{
+    this->setPatchType(ptf);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void immersedBoundaryOmegaWallFunctionFvPatchScalarField::autoMap
+(
+    const fvPatchFieldMapper&
+)
+{
+    scalarField::operator=(this->patchInternalField());
+
+    // Resize refValue as well.  HJ, 10/Jul/2018
+    refValue() = this->patchInternalField();
+}
+
+
+void immersedBoundaryOmegaWallFunctionFvPatchScalarField::rmap
+(
+    const fvPatchScalarField& ptf,
+    const labelList&
+)
+{}
+
+
+void immersedBoundaryOmegaWallFunctionFvPatchScalarField::updateOnMotion()
+{
+    if (size() != ibPatch().size())
+    {
+        // Use internal values, resizing the file if needed
+        scalarField::operator=(this->patchInternalField());
+
+        // Resize refValue as well.  HJ, 10/Jul/2018
+        refValue() = this->patchInternalField();
+    }
+}
+
+
+void immersedBoundaryOmegaWallFunctionFvPatchScalarField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    // Resize fields
+    if (size() != this->ibPatch().size())
+    {
+        Info<< "Resizing immersedBoundaryOmegaWallFunction in evaluate: "
+            << "patch: " << patch().size() << " field: " << size()
+            << endl;
+
+        *this == patchInternalField();
+        refValue() = patchInternalField();
+    }
+
+    // If G field is present, execute evaluation
+    // Remove the warning from the IB patch
+    // HJ, 20/May/2018
+    if (db().foundObject<volScalarField>(GName()))
+    {
+        omegaWallFunctionFvPatchScalarField::updateCoeffs();
+    }
+}
+
+
+void immersedBoundaryOmegaWallFunctionFvPatchScalarField::evaluate
+(
+    const Pstream::commsTypes commsType
+)
+{
+    // Resize fields
+    if (size() != this->ibPatch().size())
+    {
+        Info<< "Resizing immersedBoundaryOmegaWallFunction in evaluate"
+            << endl;
+
+        *this == patchInternalField();
+        refValue() = patchInternalField();
+    }
+
+    // Get non-constant reference to internal field
+    scalarField& intField = const_cast<scalarField&>(this->internalField());
+
+    // Set dead value
+    this->setDeadValues(intField);
+
+    omegaWallFunctionFvPatchScalarField::evaluate(commsType);
+}
+
+
+void immersedBoundaryOmegaWallFunctionFvPatchScalarField::manipulateMatrix
+(
+    fvScalarMatrix& matrix
+)
+{
+    setDeadValues(matrix);
+
+    omegaWallFunctionFvPatchScalarField::manipulateMatrix(matrix);
+}
+
+
+void immersedBoundaryOmegaWallFunctionFvPatchScalarField::write
+(
+    Ostream& os
+) const
+{
+    fvPatchScalarField::write(os);
+    writeLocalEntries(os);
+    this->writeDeadData(os);
+
+    // The value entry needs to be written with zero size
+    scalarField::null().writeEntry("value", os);
+    // this->writeEntry("value", os);
+
+    writeField(*this);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeField
+(
+    fvPatchScalarField,
+    immersedBoundaryOmegaWallFunctionFvPatchScalarField
+);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace compressible
+} // End namespace Foam
+
+// ************************************************************************* //

src/immersedBoundary/immersedBoundaryTurbulence/compressible/wallFunctions/immersedBoundaryOmegaWallFunctions/immersedBoundaryOmegaWallFunctionFvPatchScalarField.H

@@ -0,0 +1,208 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::compressible::RASModels::omegaWallFunctionFvPatchScalarField
+
+Description
+    Provides a wall function boundary condition/constraint on omega
+
+    Computed value is:
+
+        omega = sqrt(omega_vis^2 + omega_log^2)
+
+    where
+        omega_vis = omega in viscous region
+        omega_log = omega in logarithmic region
+
+    Model described by Eq.(15) of:
+    @verbatim
+        Menter, F., Esch, T.
+        "Elements of Industrial Heat Transfer Prediction"
+        16th Brazilian Congress of Mechanical Engineering (COBEM),
+        Nov. 2001
+    @endverbatim
+
+Description
+    Boundary condition for omega when using wall functions
+    on immersed boundary patches
+
+Author
+    Hrvoje Jasak, Wikki Ltd.  All rights reserved
+
+SourceFiles
+    immersedBoundaryOmegaWallFunctionFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef immersedBoundaryOmegaWallFunctionFvPatchScalarField_H
+#define immersedBoundaryOmegaWallFunctionFvPatchScalarField_H
+
+#include "omegaWallFunctionFvPatchScalarField.H"
+#include "immersedBoundaryFieldBase.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace compressible
+{
+namespace RASModels
+{
+
+/*---------------------------------------------------------------------------*\
+     Class immersedBoundaryOmegaWallFunctionFvPatchScalarField Declaration
+\*---------------------------------------------------------------------------*/
+
+class immersedBoundaryOmegaWallFunctionFvPatchScalarField
+:
+    public omegaWallFunctionFvPatchScalarField,
+    public immersedBoundaryFieldBase<scalar>
+{
+public:
+
+    //- Runtime type information
+    TypeName("compressible::immersedBoundaryOmegaWallFunction");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        immersedBoundaryOmegaWallFunctionFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        immersedBoundaryOmegaWallFunctionFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        // immersedBoundaryOmegaWallFunctionFvPatchScalarField
+        //  onto a new patch
+        immersedBoundaryOmegaWallFunctionFvPatchScalarField
+        (
+            const immersedBoundaryOmegaWallFunctionFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        immersedBoundaryOmegaWallFunctionFvPatchScalarField
+        (
+            const immersedBoundaryOmegaWallFunctionFvPatchScalarField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new immersedBoundaryOmegaWallFunctionFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        immersedBoundaryOmegaWallFunctionFvPatchScalarField
+        (
+            const immersedBoundaryOmegaWallFunctionFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new immersedBoundaryOmegaWallFunctionFvPatchScalarField
+                (
+                    *this,
+                    iF
+                )
+            );
+        }
+
+
+    // Member functions
+
+        // Mapping functions
+
+            //- Map (and resize as needed) from self given a mapping object
+            virtual void autoMap
+            (
+                const fvPatchFieldMapper&
+            );
+
+            //- Reverse map the given fvPatchField onto this fvPatchField
+            virtual void rmap
+            (
+                const fvPatchScalarField&,
+                const labelList&
+            );
+
+            //- Update on mesh motion
+            virtual void updateOnMotion();
+
+
+        // Evaluation functions
+
+            //- Update the coefficients associated with the patch field
+            virtual void updateCoeffs();
+
+            //- Evaluate the patchField
+            virtual void evaluate
+            (
+                const Pstream::commsTypes  = Pstream::blocking
+            );
+
+            //- Manipulate a matrix
+            virtual void manipulateMatrix(fvScalarMatrix& matrix);
+
+
+        // I-O
+
+            //- Write
+            void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace compressible
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/immersedBoundary/immersedBoundaryTurbulence/incompressible/Make/files

@@ -3,7 +3,4 @@ wallFunctions/immersedBoundaryNutWallFunction/immersedBoundaryNutWallFunctionFvP
 wallFunctions/immersedBoundaryEpsilonWallFunction/immersedBoundaryEpsilonWallFunctionFvPatchScalarField.C
 wallFunctions/immersedBoundaryOmegaWallFunctions/immersedBoundaryOmegaWallFunctionFvPatchScalarField.C
 
-turbulenceModels/incompressible/RAS/backwardsCompability/wallFunctions/backwardsCompatibilityIbWallFunctions.C
-turbulenceModels/incompressible/RAS/immersedBoundaryKOmegaSST/immersedBoundaryKOmegaSST.C
-
-LIB = $(FOAM_LIBBIN)/libimmersedBoundaryTurbulence
+LIB = $(FOAM_LIBBIN)/libincompressibleImmersedBoundaryTurbulence

src/immersedBoundary/immersedBoundaryTurbulence/incompressible/wallFunctions/immersedBoundaryEpsilonWallFunction/immersedBoundaryEpsilonWallFunctionFvPatchScalarField.C

@@ -62,14 +62,31 @@ immersedBoundaryEpsilonWallFunctionFvPatchScalarField
     const dictionary& dict
 )
 :
-    epsilonWallFunctionFvPatchScalarField(p, iF, dict),
+    epsilonWallFunctionFvPatchScalarField(p, iF),   // Do not read mixed data
     immersedBoundaryFieldBase<scalar>
     (
         p,
         Switch(dict.lookup("setDeadValue")),
         readScalar(dict.lookup("deadValue"))
     )
-{}
+{
+    // Since patch does not read a dictionary, the patch type needs to be read
+    // manually.  HJ, 6/Sep/2018
+    this->readPatchType(dict);
+
+    if (!isType<immersedBoundaryFvPatch>(p))
+    {
+        FatalIOErrorInFunction(dict)
+            << "\n    patch type '" << p.type()
+            << "' not constraint type '" << typeName << "'"
+            << "\n    for patch " << p.name()
+            << " of field " << this->dimensionedInternalField().name()
+            << " in file " << this->dimensionedInternalField().objectPath()
+            << exit(FatalIOError);
+    }
+
+    scalarField::operator=(this->patchInternalField());
+}
 
 
 immersedBoundaryEpsilonWallFunctionFvPatchScalarField::
@@ -81,14 +98,35 @@ immersedBoundaryEpsilonWallFunctionFvPatchScalarField
     const fvPatchFieldMapper& mapper
 )
 :
-    epsilonWallFunctionFvPatchScalarField(ptf, p, iF, mapper),
+    epsilonWallFunctionFvPatchScalarField(p, iF), // Do not map mixed data.  Set patchType later
     immersedBoundaryFieldBase<scalar>
     (
         p,
         ptf.setDeadValue(),
         ptf.deadValue()
     )
-{}
+{
+    // Note: NO MAPPING.  Fields are created on the immersed boundary
+    // HJ, 12/Apr/2012
+    if (!isType<immersedBoundaryFvPatch>(p))
+    {
+        FatalErrorInFunction
+            << "\n    patch type '" << p.type()
+            << "' not constraint type '" << typeName << "'"
+            << "\n    for patch " << p.name()
+            << " of field " << this->dimensionedInternalField().name()
+            << " in file " << this->dimensionedInternalField().objectPath()
+            << exit(FatalIOError);
+    }
+
+    // Copy the patch type since mixed data was not mapped
+    this->setPatchType(ptf);
+
+    // On creation of the mapped field, the internal field is dummy and
+    // cannot be used.  Initialise the value to avoid errors
+    // HJ, 1/Dec/2017
+    scalarField::operator=(SMALL);
+}
 
 
 immersedBoundaryEpsilonWallFunctionFvPatchScalarField::
@@ -182,7 +220,6 @@ void immersedBoundaryEpsilonWallFunctionFvPatchScalarField::updateCoeffs()
     // HJ, 20/May/2018
     if (db().foundObject<volScalarField>(GName()))
     {
-        Info<< "Updating epsilon" << endl;
         epsilonWallFunctionFvPatchScalarField::updateCoeffs();
     }
 }
@@ -213,6 +250,17 @@ void immersedBoundaryEpsilonWallFunctionFvPatchScalarField::evaluate
 }
 
 
+void immersedBoundaryEpsilonWallFunctionFvPatchScalarField::manipulateMatrix
+(
+    fvScalarMatrix& matrix
+)
+{
+    setDeadValues(matrix);
+
+    epsilonWallFunctionFvPatchScalarField::manipulateMatrix(matrix);
+}
+
+
 void immersedBoundaryEpsilonWallFunctionFvPatchScalarField::write
 (
     Ostream& os

src/immersedBoundary/immersedBoundaryTurbulence/incompressible/wallFunctions/immersedBoundaryEpsilonWallFunction/immersedBoundaryEpsilonWallFunctionFvPatchScalarField.H

@@ -170,6 +170,9 @@ public:
                 const Pstream::commsTypes  = Pstream::blocking
             );
 
+            //- Manipulate a matrix
+            virtual void manipulateMatrix(fvScalarMatrix& matrix);
+
 
         // I-O
 

src/immersedBoundary/immersedBoundaryTurbulence/incompressible/wallFunctions/immersedBoundaryKqRWallFunction/immersedBoundaryKqRWallFunctionFvPatchField.C

@@ -82,7 +82,7 @@ immersedBoundaryKqRWallFunctionFvPatchField
     const dictionary& dict
 )
 :
-    kqRWallFunctionFvPatchField<Type>(p, iF, dict),
+    kqRWallFunctionFvPatchField<Type>(p, iF),   // Do not read mixed data
     immersedBoundaryFieldBase<Type>
     (
         p,
@@ -90,10 +90,10 @@ immersedBoundaryKqRWallFunctionFvPatchField
         pTraits<Type>(dict.lookup("deadValue"))
     )
 {
-    this->checkType();
     this->readPatchType(dict);
+    this->checkType();
 
-    *this == this->patchInternalField();
+    Field<Type>::operator=(this->patchInternalField());
 }
 
 
@@ -107,7 +107,7 @@ immersedBoundaryKqRWallFunctionFvPatchField
     const fvPatchFieldMapper& mapper
 )
 :
-    kqRWallFunctionFvPatchField<Type>(ptf, p, iF, mapper),
+    kqRWallFunctionFvPatchField<Type>(p, iF), // Do not map mixed data.  Set patchType later
     immersedBoundaryFieldBase<Type>
     (
         p,
@@ -117,6 +117,10 @@ immersedBoundaryKqRWallFunctionFvPatchField
 {
     this->setPatchType(ptf);
     this->checkType();
+
+    // cannot be used.  Initialise the value to avoid errors
+    // HJ, 1/Dec/2017
+    Field<Type>::operator=(pTraits<Type>::zero);
 }
 
 
@@ -216,6 +220,16 @@ void immersedBoundaryKqRWallFunctionFvPatchField<Type>::evaluate
 }
 
 
+template<class Type>
+void immersedBoundaryKqRWallFunctionFvPatchField<Type>::manipulateMatrix
+(
+    fvMatrix<Type>& matrix
+)
+{
+    this->setDeadValues(matrix);
+}
+
+
 template<class Type>
 void immersedBoundaryKqRWallFunctionFvPatchField<Type>::write(Ostream& os) const
 {

src/immersedBoundary/immersedBoundaryTurbulence/incompressible/wallFunctions/immersedBoundaryKqRWallFunction/immersedBoundaryKqRWallFunctionFvPatchField.H

@@ -167,6 +167,9 @@ public:
                 const Pstream::commsTypes commsType = Pstream::blocking
             );
 
+            //- Manipulate a matrix
+            virtual void manipulateMatrix(fvMatrix<Type>& matrix);
+
 
         // I-O
 

src/immersedBoundary/immersedBoundaryTurbulence/incompressible/wallFunctions/immersedBoundaryNutWallFunction/immersedBoundaryNutWallFunctionFvPatchScalarField.C

@@ -61,10 +61,23 @@ immersedBoundaryNutWallFunctionFvPatchScalarField
     const dictionary& dict
 )
 :
-    nutkWallFunctionFvPatchScalarField(p, iF, dict),
+    nutkWallFunctionFvPatchScalarField(p, iF),   // Do not read mixed data
     immersedBoundaryFieldBase<scalar>(p, true, 1e-6)
 {
     this->readPatchType(dict);
+
+    if (!isType<immersedBoundaryFvPatch>(p))
+    {
+        FatalIOErrorInFunction(dict)
+            << "\n    patch type '" << p.type()
+            << "' not constraint type '" << typeName << "'"
+            << "\n    for patch " << p.name()
+            << " of field " << this->dimensionedInternalField().name()
+            << " in file " << this->dimensionedInternalField().objectPath()
+            << exit(FatalIOError);
+    }
+
+    scalarField::operator=(this->patchInternalField());
 }
 
 
@@ -77,10 +90,28 @@ immersedBoundaryNutWallFunctionFvPatchScalarField
     const fvPatchFieldMapper& mapper
 )
 :
-    nutkWallFunctionFvPatchScalarField(ptf, p, iF, mapper),
+    nutkWallFunctionFvPatchScalarField(p, iF), // Do not map mixed data.  Set patchType later
     immersedBoundaryFieldBase<scalar>(p, true, 1e-6)
 {
+    // Note: NO MAPPING.  Fields are created on the immersed boundary
+    // HJ, 12/Apr/2012
+    if (!isType<immersedBoundaryFvPatch>(p))
+    {
+        FatalErrorInFunction
+            << "\n    patch type '" << p.type()
+            << "' not constraint type '" << typeName << "'"
+            << "\n    for patch " << p.name()
+            << " of field " << this->dimensionedInternalField().name()
+            << " in file " << this->dimensionedInternalField().objectPath()
+            << exit(FatalIOError);
+    }
+
     this->setPatchType(ptf);
+
+    // On creation of the mapped field, the internal field is dummy and
+    // cannot be used.  Initialise the value to avoid errors
+    // HJ, 1/Dec/2017
+    scalarField::operator=(scalar(0));
 }
 
 
@@ -151,7 +182,7 @@ void immersedBoundaryNutWallFunctionFvPatchScalarField::evaluate
         Info<< "Resizing immersedBoundaryNutWallFunction in evaluate"
             << endl;
 
-        *this == patchInternalField();
+        scalarField::operator=(patchInternalField());
     }
 
     // Get non-constant reference to internal field

src/immersedBoundary/immersedBoundaryTurbulence/incompressible/wallFunctions/immersedBoundaryOmegaWallFunctions/immersedBoundaryOmegaWallFunctionFvPatchScalarField.C

@@ -60,17 +60,28 @@ immersedBoundaryOmegaWallFunctionFvPatchScalarField
     const dictionary& dict
 )
 :
-    omegaWallFunctionFvPatchScalarField(p, iF, dict),
+    omegaWallFunctionFvPatchScalarField(p, iF),    // Do not read mixed data
     immersedBoundaryFieldBase<scalar>
     (
         p,
         Switch(dict.lookup("setDeadValue")),
         readScalar(dict.lookup("deadValue"))
-    )    
+    )
 {
     this->readPatchType(dict);
 
-    *this == this->patchInternalField();
+    if (!isType<immersedBoundaryFvPatch>(p))
+    {
+        FatalIOErrorInFunction(dict)
+            << "\n    patch type '" << p.type()
+            << "' not constraint type '" << typeName << "'"
+            << "\n    for patch " << p.name()
+            << " of field " << this->dimensionedInternalField().name()
+            << " in file " << this->dimensionedInternalField().objectPath()
+            << exit(FatalIOError);
+    }
+
+    scalarField::operator=(this->patchInternalField());
 }
 
 
@@ -83,7 +94,7 @@ immersedBoundaryOmegaWallFunctionFvPatchScalarField
     const fvPatchFieldMapper& mapper
 )
 :
-    omegaWallFunctionFvPatchScalarField(ptf, p, iF, mapper),
+    omegaWallFunctionFvPatchScalarField(p, iF), // Do not map mixed data.  Set patchType later
     immersedBoundaryFieldBase<scalar>
     (
         p,
@@ -91,7 +102,25 @@ immersedBoundaryOmegaWallFunctionFvPatchScalarField
         ptf.deadValue()
     )
 {
+    // Note: NO MAPPING.  Fields are created on the immersed boundary
+    // HJ, 12/Apr/2012
+    if (!isType<immersedBoundaryFvPatch>(p))
+    {
+        FatalErrorInFunction
+            << "\n    patch type '" << p.type()
+            << "' not constraint type '" << typeName << "'"
+            << "\n    for patch " << p.name()
+            << " of field " << this->dimensionedInternalField().name()
+            << " in file " << this->dimensionedInternalField().objectPath()
+            << exit(FatalIOError);
+    }
+
     this->setPatchType(ptf);
+
+    // On creation of the mapped field, the internal field is dummy and
+    // cannot be used.  Initialise the value to avoid errors
+    // HJ, 1/Dec/2017
+    scalarField::operator=(SMALL);
 }
 
 
@@ -220,6 +249,17 @@ void immersedBoundaryOmegaWallFunctionFvPatchScalarField::evaluate
 }
 
 
+void immersedBoundaryOmegaWallFunctionFvPatchScalarField::manipulateMatrix
+(
+    fvScalarMatrix& matrix
+)
+{
+    setDeadValues(matrix);
+
+    omegaWallFunctionFvPatchScalarField::manipulateMatrix(matrix);
+}
+
+
 void immersedBoundaryOmegaWallFunctionFvPatchScalarField::write
 (
     Ostream& os

src/immersedBoundary/immersedBoundaryTurbulence/incompressible/wallFunctions/immersedBoundaryOmegaWallFunctions/immersedBoundaryOmegaWallFunctionFvPatchScalarField.H

@@ -184,6 +184,9 @@ public:
                 const Pstream::commsTypes  = Pstream::blocking
             );
 
+            //- Manipulate a matrix
+            virtual void manipulateMatrix(fvScalarMatrix& matrix);
+
 
         // I-O
 

src/immersedBoundary/immersedBoundaryTurbulence/turbulenceModels/incompressible/RAS/backwardsCompability/wallFunctions/backwardsCompatibilityIbWallFunctions.C

@@ -1,336 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | foam-extend: Open Source CFD
-   \\    /   O peration     | Version:     4.1
-    \\  /    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 "backwardsCompatibilityIbWallFunctions.H"
-
-#include "calculatedFvPatchField.H"
-#include "nutkWallFunctionFvPatchScalarField.H"
-#include "nutLowReWallFunctionFvPatchScalarField.H"
-#include "epsilonWallFunctionFvPatchScalarField.H"
-#include "kqRWallFunctionFvPatchField.H"
-#include "RWallFunctionFvPatchSymmTensorField.H"
-#include "omegaWallFunctionFvPatchScalarField.H"
-#include "immersedBoundaryEpsilonWallFunctionFvPatchScalarField.H"
-#include "immersedBoundaryKqRWallFunctionFvPatchField.H"
-#include "immersedBoundaryOmegaWallFunctionFvPatchScalarField.H"
-#include "immersedBoundaryNutWallFunctionFvPatchScalarField.H"
-#include "immersedBoundaryFvPatch.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace incompressible
-{
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-tmp<volScalarField> autoCreateIbNut
-(
-    const word& fieldName,
-    const fvMesh& mesh,
-    const objectRegistry& obj
-)
-{
-    IOobject nutHeader
-    (
-        fieldName,
-        mesh.time().timeName(),
-        obj,
-        IOobject::MUST_READ,
-        IOobject::NO_WRITE,
-        false
-    );
-
-    if (nutHeader.headerOk())
-    {
-        return tmp<volScalarField>(new volScalarField(nutHeader, mesh));
-    }
-    else
-    {
-        Info<< "--> Creating " << fieldName
-            << " to employ run-time selectable wall functions" << endl;
-
-        const fvBoundaryMesh& bm = mesh.boundary();
-
-        wordList nutBoundaryTypes(bm.size());
-
-        forAll(bm, patchI)
-        {
-            if (isA<immersedBoundaryFvPatch>(bm[patchI]) && bm[patchI].isWall())
-            {
-                nutBoundaryTypes[patchI] =
-                    RASModels::immersedBoundaryNutWallFunctionFvPatchScalarField::typeName;
-            }
-            else if (bm[patchI].isWall())
-            {
-                nutBoundaryTypes[patchI] =
-                    RASModels::nutkWallFunctionFvPatchScalarField::typeName;
-            }
-            else
-            {
-                nutBoundaryTypes[patchI] =
-                    calculatedFvPatchField<scalar>::typeName;
-            }
-        }
-
-        tmp<volScalarField> nut
-        (
-            new volScalarField
-            (
-                IOobject
-                (
-                    fieldName,
-                    mesh.time().timeName(),
-                    mesh,
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE,
-                    false
-                ),
-                mesh,
-                dimensionedScalar("zero", dimArea/dimTime, 0.0),
-                nutBoundaryTypes
-            )
-        );
-
-        Info<< "    Writing new " << fieldName << endl;
-        nut().write();
-
-        return nut;
-    }
-}
-
-
-tmp<volScalarField> autoCreateIbNut
-(
-    const word& fieldName,
-    const fvMesh& mesh
-)
-{
-    return autoCreateIbNut(fieldName, mesh, mesh);
-}
-
-
-tmp<volScalarField> autoCreateIbEpsilon
-(
-    const word& fieldName,
-    const fvMesh& mesh,
-    const objectRegistry& obj
-)
-{
-    return
-        autoCreateWallFunctionField
-        <
-            scalar,
-            RASModels::epsilonWallFunctionFvPatchScalarField,
-            RASModels::immersedBoundaryEpsilonWallFunctionFvPatchScalarField
-        >
-        (
-            fieldName,
-            mesh,
-            obj
-        );
-}
-
-
-tmp<volScalarField> autoCreateIbEpsilon
-(
-    const word& fieldName,
-    const fvMesh& mesh
-)
-{
-    return
-        autoCreateWallFunctionField
-        <
-            scalar,
-            RASModels::epsilonWallFunctionFvPatchScalarField,
-            RASModels::immersedBoundaryEpsilonWallFunctionFvPatchScalarField
-        >
-        (
-            fieldName,
-            mesh,
-            mesh
-        );
-}
-
-
-tmp<volScalarField> autoCreateIbOmega
-(
-    const word& fieldName,
-    const fvMesh& mesh,
-    const objectRegistry& obj
-)
-{
-    return
-        autoCreateWallFunctionField
-        <
-            scalar,
-            RASModels::omegaWallFunctionFvPatchScalarField,
-            RASModels::immersedBoundaryOmegaWallFunctionFvPatchScalarField
-        >
-        (
-            fieldName,
-            mesh,
-            obj
-        );
-}
-
-
-tmp<volScalarField> autoCreateIbOmega
-(
-    const word& fieldName,
-    const fvMesh& mesh
-)
-{
-    return
-        autoCreateWallFunctionField
-        <
-            scalar,
-            RASModels::omegaWallFunctionFvPatchScalarField,
-            RASModels::immersedBoundaryOmegaWallFunctionFvPatchScalarField
-        >
-        (
-            fieldName,
-            mesh,
-            mesh
-        );
-}
-
-
-tmp<volScalarField> autoCreateIbK
-(
-    const word& fieldName,
-    const fvMesh& mesh,
-    const objectRegistry& obj
-)
-{
-    return
-        autoCreateWallFunctionField
-        <
-            scalar,
-            RASModels::kqRWallFunctionFvPatchField<scalar>,
-            RASModels::immersedBoundaryKqRWallFunctionFvPatchField<scalar>
-        >
-        (
-            fieldName,
-            mesh,
-            obj
-        );
-}
-
-
-tmp<volScalarField> autoCreateIbK
-(
-    const word& fieldName,
-    const fvMesh& mesh
-)
-{
-    return
-        autoCreateWallFunctionField
-        <
-            scalar,
-            RASModels::kqRWallFunctionFvPatchField<scalar>,
-            RASModels::immersedBoundaryKqRWallFunctionFvPatchField<scalar>
-        >
-        (
-            fieldName,
-            mesh,
-            mesh
-        );
-}
-
-
-tmp<volScalarField> autoCreateIbQ
-(
-    const word& fieldName,
-    const fvMesh& mesh,
-    const objectRegistry& obj
-)
-{
-    return
-        autoCreateWallFunctionField
-        <
-            scalar,
-            RASModels::kqRWallFunctionFvPatchField<scalar>,
-            RASModels::immersedBoundaryKqRWallFunctionFvPatchField<scalar>
-        >
-        (
-            fieldName,
-            mesh,
-            obj
-        );
-}
-
-
-tmp<volScalarField> autoCreateIbQ
-(
-    const word& fieldName,
-    const fvMesh& mesh
-)
-{
-    return
-        autoCreateWallFunctionField
-        <
-            scalar,
-            RASModels::kqRWallFunctionFvPatchField<scalar>,
-            RASModels::immersedBoundaryKqRWallFunctionFvPatchField<scalar>
-        >
-        (
-            fieldName,
-            mesh,
-            mesh
-        );
-}
-
-
-tmp<volSymmTensorField> autoCreateIbR
-(
-    const word& fieldName,
-    const fvMesh& mesh,
-    const objectRegistry& obj
-)
-{
-    return
-        autoCreateWallFunctionField
-        <
-            symmTensor,
-            RASModels::kqRWallFunctionFvPatchField<symmTensor>,
-            RASModels::immersedBoundaryKqRWallFunctionFvPatchField<symmTensor>
-        >
-        (
-            fieldName,
-            mesh,
-            obj
-        );
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace incompressible
-} // End namespace Foam
-
-// ************************************************************************* //
-

src/immersedBoundary/immersedBoundaryTurbulence/turbulenceModels/incompressible/RAS/backwardsCompability/wallFunctions/backwardsCompatibilityIbWallFunctions.H

@@ -1,158 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | foam-extend: Open Source CFD
-   \\    /   O peration     | Version:     4.1
-    \\  /    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::incompressible
-
-Description
-    Auto creation of fields to provide backwards compatibility with
-    runtime selectable wall functions
-
-SourceFiles
-    backwardsCompatibilityIbWallFunctions.C
-    backwardsCompatibilityIbWallFunctionsTemplates.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef backwardsCompatibilityIbWallFunctions_H
-#define backwardsCompatibilityIbWallFunctions_H
-
-#include "fvMesh.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace incompressible
-{
-    //- nut
-    tmp<volScalarField> autoCreateIbNut
-    (
-        const word& fieldName,
-        const fvMesh& mesh,
-        const objectRegistry& obj
-    );
-
-    tmp<volScalarField> autoCreateIbNut
-    (
-        const word& fieldName,
-        const fvMesh& mesh
-    );
-
-    //- epsilon
-    tmp<volScalarField> autoCreateIbEpsilon
-    (
-        const word& fieldName,
-        const fvMesh& mesh,
-        const objectRegistry& obj
-    );
-
-    tmp<volScalarField> autoCreateIbEpsilon
-    (
-        const word& fieldName,
-        const fvMesh& mesh
-    );
-
-    //- omega
-    tmp<volScalarField> autoCreateIbOmega
-    (
-        const word& fieldName,
-        const fvMesh& mesh,
-        const objectRegistry& obj
-    );
-
-    tmp<volScalarField> autoCreateIbOmega
-    (
-        const word& fieldName,
-        const fvMesh& mesh
-    );
-
-    //- k
-    tmp<volScalarField> autoCreateIbK
-    (
-        const word& fieldName,
-        const fvMesh& mesh,
-        const objectRegistry& obj
-    );
-
-    tmp<volScalarField> autoCreateIbK
-    (
-        const word& fieldName,
-        const fvMesh& mesh
-    );
-
-    //- Q
-    tmp<volScalarField> autoCreateIbQ
-    (
-        const word& fieldName,
-        const fvMesh& mesh,
-        const objectRegistry& obj
-    );
-
-    tmp<volScalarField> autoCreateIbQ
-    (
-        const word& fieldName,
-        const fvMesh& mesh
-    );
-
-    //- R
-    tmp<volSymmTensorField> autoCreateIbR
-    (
-        const word& fieldName,
-        const fvMesh& mesh,
-        const objectRegistry& obj
-    );
-
-    tmp<volSymmTensorField> autoCreateIbR
-    (
-        const word& fieldName,
-        const fvMesh& mesh
-    );
-
-    //- Helper function to create the new field
-    template<class Type, class PatchType, class ibPatchName >
-    tmp<GeometricField<Type, fvPatchField, volMesh> >
-    autoCreateIbWallFunctionField
-    (
-        const word& fieldName,
-        const fvMesh& mesh,
-        const objectRegistry& obj
-    );
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace incompressible
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#ifdef NoRepository
-#   include "backwardsCompatibilityIbWallFunctionsTemplates.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/immersedBoundary/immersedBoundaryTurbulence/turbulenceModels/incompressible/RAS/backwardsCompability/wallFunctions/backwardsCompatibilityIbWallFunctionsTemplates.C

@@ -1,184 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | foam-extend: Open Source CFD
-   \\    /   O peration     | Version:     4.1
-    \\  /    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 "backwardsCompatibilityIbWallFunctions.H"
-#include "foamTime.H"
-#include "OSspecific.H"
-#include "immersedBoundaryFvPatch.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace incompressible
-{
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-template<class Type, class PatchType, class ibPatchType>
-tmp<GeometricField<Type, fvPatchField, volMesh> >
-autoCreateWallFunctionField
-(
-    const word& fieldName,
-    const fvMesh& mesh,
-    const objectRegistry& obj
-)
-{
-    IOobject nutHeader
-    (
-        "nut",
-        mesh.time().timeName(),
-        obj,
-        IOobject::MUST_READ
-    );
-
-    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
-
-    if (nutHeader.headerOk())
-    {
-        return tmp<fieldType>
-        (
-            new fieldType
-            (
-                IOobject
-                (
-                    fieldName,
-                    mesh.time().timeName(),
-                    obj,
-                    IOobject::MUST_READ,
-                    IOobject::NO_WRITE,
-                    false
-                ),
-                mesh
-            )
-        );
-    }
-    else
-    {
-        Info<< "--> Upgrading " << fieldName
-            << " to employ run-time selectable wall functions" << endl;
-
-        // Read existing field
-        IOobject ioObj
-        (
-            fieldName,
-            mesh.time().timeName(),
-            obj,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE,
-            false
-        );
-
-        tmp<fieldType> fieldOrig
-        (
-            new fieldType
-            (
-                ioObj,
-                mesh
-            )
-        );
-
-        // rename file
-        Info<< "    Backup original " << fieldName << " to "
-            << fieldName << ".old" << endl;
-        mvBak(ioObj.objectPath(), "old");
-
-
-        PtrList<fvPatchField<Type> > newPatchFields(mesh.boundary().size());
-
-        const fvBoundaryMesh& bm = mesh.boundary();
-
-        forAll(newPatchFields, patchI)
-        {
-            if (isA<immersedBoundaryFvPatch>(bm[patchI]) && bm[patchI].isWall())
-            {
-                newPatchFields.set
-                (
-                    patchI,
-                    new ibPatchType
-                    (
-                        mesh.boundary()[patchI],
-                        fieldOrig().dimensionedInternalField()
-                    )
-                );
-                newPatchFields[patchI] == fieldOrig().boundaryField()[patchI];
-            }
-            else if (bm[patchI].isWall())
-            {
-                newPatchFields.set
-                (
-                    patchI,
-                    new PatchType
-                    (
-                        mesh.boundary()[patchI],
-                        fieldOrig().dimensionedInternalField()
-                    )
-                );
-                newPatchFields[patchI] == fieldOrig().boundaryField()[patchI];
-            }
-            else
-            {
-                newPatchFields.set
-                (
-                    patchI,
-                    fieldOrig().boundaryField()[patchI].clone()
-                );
-            }
-        }
-
-        tmp<fieldType> fieldNew
-        (
-            new fieldType
-            (
-                IOobject
-                (
-                    fieldName,
-                    mesh.time().timeName(),
-                    obj,
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE,
-                    false
-                ),
-                mesh,
-                fieldOrig().dimensions(),
-                fieldOrig().internalField(),
-                newPatchFields
-            )
-        );
-
-        Info<< "    Writing updated " << fieldName << endl;
-        fieldNew().write();
-
-        return fieldNew;
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace incompressible
-} // End namespace Foam
-
-// ************************************************************************* //

src/immersedBoundary/immersedBoundaryTurbulence/turbulenceModels/incompressible/RAS/immersedBoundaryKOmegaSST/immersedBoundaryKOmegaSST.C

@@ -1,483 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | foam-extend: Open Source CFD
-   \\    /   O peration     | Version:     4.1
-    \\  /    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 "immersedBoundaryKOmegaSST.H"
-#include "addToRunTimeSelectionTable.H"
-
-#include "backwardsCompatibilityIbWallFunctions.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace incompressible
-{
-namespace RASModels
-{
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-defineTypeNameAndDebug(immersedBoundaryKOmegaSST, 0);
-addToRunTimeSelectionTable(RASModel, immersedBoundaryKOmegaSST, dictionary);
-
-// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
-
-tmp<volScalarField> immersedBoundaryKOmegaSST::F1
-(
-const volScalarField& CDkOmega
-) const
-{
-    volScalarField CDkOmegaPlus = max
-    (
-        CDkOmega,
-        dimensionedScalar("1.0e-10", dimless/sqr(dimTime), 1.0e-10)
-    );
-
-    volScalarField arg1 = min
-    (
-        min
-        (
-            max
-            (
-                (scalar(1)/betaStar_)*sqrt(k_)/(omega_*y_),
-                scalar(500)*nu()/(sqr(y_)*omega_)
-            ),
-            (4*alphaOmega2_)*k_/(CDkOmegaPlus*sqr(y_))
-        ),
-        scalar(10)
-    );
-
-    return tanh(pow4(arg1));
-}
-
-
-tmp<volScalarField> immersedBoundaryKOmegaSST::F2() const
-{
-    volScalarField arg2 = min
-    (
-        max
-        (
-            (scalar(2)/betaStar_)*sqrt(k_)/(omega_*y_),
-            scalar(500)*nu()/(sqr(y_)*omega_)
-        ),
-        scalar(100)
-    );
-
-    return tanh(sqr(arg2));
-}
-
-
-tmp<volScalarField> immersedBoundaryKOmegaSST::F3() const
-{
-    tmp<volScalarField> arg3 = min
-    (
-        150*nu()/(omega_*sqr(y_)),
-        scalar(10)
-    );
-
-    return 1 - tanh(pow4(arg3));
-}
-
-
-tmp<volScalarField> immersedBoundaryKOmegaSST::F23() const
-{
-    tmp<volScalarField> f23(F2());
-
-    if (F3_)
-    {
-        f23() *= F3();
-    }
-
-    return f23;
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-immersedBoundaryKOmegaSST::immersedBoundaryKOmegaSST
-(
-    const volVectorField& U,
-    const surfaceScalarField& phi,
-    transportModel& transport,
-    const word& turbulenceModelName,
-    const word& modelName
-)
-:
-    RASModel(modelName, U, phi, transport, turbulenceModelName),
-
-    alphaK1_
-    (
-        dimensionedScalar::lookupOrAddToDict
-        (
-            "alphaK1",
-            coeffDict_,
-            0.85
-        )
-    ),
-    alphaK2_
-    (
-        dimensionedScalar::lookupOrAddToDict
-        (
-            "alphaK2",
-            coeffDict_,
-            1.0
-        )
-    ),
-    alphaOmega1_
-    (
-        dimensionedScalar::lookupOrAddToDict
-        (
-            "alphaOmega1",
-            coeffDict_,
-            0.5
-        )
-    ),
-    alphaOmega2_
-    (
-        dimensionedScalar::lookupOrAddToDict
-        (
-            "alphaOmega2",
-            coeffDict_,
-            0.856
-        )
-    ),
-    gamma1_
-    (
-        dimensionedScalar::lookupOrAddToDict
-        (
-            "gamma1",
-            coeffDict_,
-            5.0/9.0
-        )
-    ),
-    gamma2_
-    (
-        dimensionedScalar::lookupOrAddToDict
-        (
-            "gamma2",
-            coeffDict_,
-            0.44
-        )
-    ),
-    beta1_
-    (
-        dimensionedScalar::lookupOrAddToDict
-        (
-            "beta1",
-            coeffDict_,
-            0.075
-        )
-    ),
-    beta2_
-    (
-        dimensionedScalar::lookupOrAddToDict
-        (
-            "beta2",
-            coeffDict_,
-            0.0828
-        )
-    ),
-    betaStar_
-    (
-        dimensionedScalar::lookupOrAddToDict
-        (
-            "betaStar",
-            coeffDict_,
-            0.09
-        )
-    ),
-    a1_
-    (
-        dimensionedScalar::lookupOrAddToDict
-        (
-            "a1",
-            coeffDict_,
-            0.31
-        )
-    ),
-    b1_
-    (
-        dimensionedScalar::lookupOrAddToDict
-        (
-            "b1",
-            coeffDict_,
-            1.0
-        )
-    ),
-    c1_
-    (
-        dimensionedScalar::lookupOrAddToDict
-        (
-            "c1",
-            coeffDict_,
-            10.0
-        )
-    ),
-    F3_
-    (
-        Switch::lookupOrAddToDict
-        (
-            "F3",
-            coeffDict_,
-            false
-        )
-    ),
-
-    y_(mesh_),
-
-    k_
-    (
-        IOobject
-        (
-            "k",
-            runTime_.timeName(),
-            U_.db(),
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        autoCreateIbK("k", mesh_, U_.db())
-    ),
-    omega_
-    (
-        IOobject
-        (
-            "omega",
-            runTime_.timeName(),
-            U_.db(),
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        autoCreateIbOmega("omega", mesh_, U_.db())
-    ),
-    nut_
-    (
-        IOobject
-        (
-            "nut",
-            runTime_.timeName(),
-            U_.db(),
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        autoCreateIbNut("nut", mesh_, U_.db())
-    )
-{
-    bound(k_, k0_);
-    bound(omega_, omega0_);
-
-    nut_ =
-    (
-        a1_*k_/
-        max
-        (
-            a1_*omega_,
-            b1_*F23()*sqrt(2.0)*mag(symm(fvc::grad(U_)))
-        )
-    );
-    nut_.correctBoundaryConditions();
-
-    printCoeffs();
-}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-tmp<volSymmTensorField> immersedBoundaryKOmegaSST::R() const
-{
-    return tmp<volSymmTensorField>
-    (
-        new volSymmTensorField
-        (
-            IOobject
-            (
-                "R",
-                runTime_.timeName(),
-                U_.db(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            ((2.0/3.0)*I)*k_ - nut_*twoSymm(fvc::grad(U_)),
-            k_.boundaryField().types()
-        )
-    );
-}
-
-
-tmp<volSymmTensorField> immersedBoundaryKOmegaSST::devReff() const
-{
-    return tmp<volSymmTensorField>
-    (
-        new volSymmTensorField
-        (
-            IOobject
-            (
-                "devRhoReff",
-                runTime_.timeName(),
-                U_.db(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            -nuEff()*dev(twoSymm(fvc::grad(U_)))
-        )
-    );
-}
-
-
-tmp<fvVectorMatrix> immersedBoundaryKOmegaSST::divDevReff() const
-{
-    const volScalarField nuEffective = nuEff();
-
-    return
-    (
-      - fvm::laplacian(nuEffective, U_)
-      - (fvc::grad(U_) & fvc::grad(nuEffective))
-    );
-}
-
-
-bool immersedBoundaryKOmegaSST::read()
-{
-    if (RASModel::read())
-    {
-        alphaK1_.readIfPresent(coeffDict());
-        alphaK2_.readIfPresent(coeffDict());
-        alphaOmega1_.readIfPresent(coeffDict());
-        alphaOmega2_.readIfPresent(coeffDict());
-        gamma1_.readIfPresent(coeffDict());
-        gamma2_.readIfPresent(coeffDict());
-        beta1_.readIfPresent(coeffDict());
-        beta2_.readIfPresent(coeffDict());
-        betaStar_.readIfPresent(coeffDict());
-        a1_.readIfPresent(coeffDict());
-        b1_.readIfPresent(coeffDict());
-        c1_.readIfPresent(coeffDict());
-        F3_.readIfPresent("F3", coeffDict());
-
-        return true;
-    }
-    else
-    {
-        return false;
-    }
-}
-
-
-void immersedBoundaryKOmegaSST::correct()
-{
-    // Bound in case of topological change
-    // HJ, 22/Aug/2007
-    if (mesh_.changing())
-    {
-        bound(k_, k0_);
-        bound(omega_, omega0_);
-    }
-
-    RASModel::correct();
-
-    if (!turbulence_)
-    {
-        return;
-    }
-
-    if (mesh_.changing())
-    {
-        y_.correct();
-    }
-
-    const volScalarField S2(2*magSqr(symm(fvc::grad(U_))));
-    volScalarField G("RASModel::G", nut_*S2);
-
-    // Update omega and G at the wall
-    omega_.boundaryField().updateCoeffs();
-
-    const volScalarField CDkOmega
-    (
-        (2*alphaOmega2_)*(fvc::grad(k_) & fvc::grad(omega_))/omega_
-    );
-
-    const volScalarField F1(this->F1(CDkOmega));
-
-    // Turbulent frequency equation
-    fvScalarMatrix omegaEqn
-    (
-        fvm::ddt(omega_)
-      + fvm::div(phi_, omega_)
-      + fvm::SuSp(-fvc::div(phi_), omega_)
-      - fvm::laplacian(DomegaEff(F1), omega_)
-     ==
-        gamma(F1)
-       *min
-        (
-            S2,
-            (c1_/a1_)*betaStar_*omega_*max(a1_*omega_, b1_*F23()*sqrt(S2))
-        )
-      - fvm::Sp(beta(F1)*omega_, omega_)
-      - fvm::SuSp
-        (
-            (F1 - scalar(1))*CDkOmega/omega_,
-            omega_
-        )
-    );
-
-    omegaEqn.relax();
-
-    // No longer needed: matrix completes at the point of solution
-    // HJ, 17/Apr/2012
-//     omegaEqn.completeAssembly();
-
-    solve(omegaEqn);
-    bound(omega_, omega0_);
-
-    // Turbulent kinetic energy equation
-    fvScalarMatrix kEqn
-    (
-        fvm::ddt(k_)
-      + fvm::div(phi_, k_)
-      + fvm::SuSp(-fvc::div(phi_), k_)
-      - fvm::laplacian(DkEff(F1), k_)
-     ==
-        min(G, c1_*betaStar_*k_*omega_)
-      - fvm::Sp(betaStar_*omega_, k_)
-    );
-
-    kEqn.relax();
-    solve(kEqn);
-    bound(k_, k0_);
-
-    // Re-calculate viscosity
-    // Fixed sqrt(2) error.  HJ, 10/Jun/2015
-    nut_ = a1_*k_/max(a1_*omega_, b1_*F23()*sqrt(S2));
-    nut_ = min(nut_, nuRatio()*nu());
-    nut_.correctBoundaryConditions();
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace RASModels
-} // End namespace incompressible
-} // End namespace Foam
-
-// ************************************************************************* //

src/immersedBoundary/immersedBoundaryTurbulence/turbulenceModels/incompressible/RAS/immersedBoundaryKOmegaSST/immersedBoundaryKOmegaSST.H

@@ -1,313 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | foam-extend: Open Source CFD
-   \\    /   O peration     | Version:     4.1
-    \\  /    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::incompressible::RASModels::immersedBoundaryKOmegaSST
-
-Description
-    Implementation of the k-omega-SST turbulence model for incompressible
-    flows.
-
-    Turbulence model described in:
-    @verbatim
-        Menter, F., Esch, T.
-        "Elements of Industrial Heat Transfer Prediction"
-        16th Brazilian Congress of Mechanical Engineering (COBEM),
-        Nov. 2001.
-    @endverbatim
-
-    with updated coefficients from
-    @verbatim
-        Menter, F. R., Kuntz, M., and Langtry, R.,
-        "Ten Years of Industrial Experience with the SST Turbulence Model",
-        Turbulence, Heat and Mass Transfer 4, 2003,
-        pp. 625 - 632.
-    @endverbatim
-
-    but with the consistent production terms from the 2001 paper as form in the
-    2003 paper is a typo, see
-    @verbatim
-        http://turbmodels.larc.nasa.gov/sst.html
-    @endverbatim
-
-    and the addition of the optional F3 term for rough walls from
-    \verbatim
-        Hellsten, A.
-        "Some Improvements in Menter’s k-omega-SST turbulence model"
-        29th AIAA Fluid Dynamics Conference,
-        AIAA-98-2554,
-        June 1998.
-    \endverbatim
-
-    Note that this implementation is written in terms of alpha diffusion
-    coefficients rather than the more traditional sigma (alpha = 1/sigma) so
-    that the blending can be applied to all coefficuients in a consistent
-    manner.  The paper suggests that sigma is blended but this would not be
-    consistent with the blending of the k-epsilon and k-omega models.
-
-    Also note that the error in the last term of equation (2) relating to
-    sigma has been corrected.
-
-    The default model coefficients correspond to the following:
-    @verbatim
-        immersedBoundaryKOmegaSSTCoeffs
-        {
-            alphaK1     0.85;
-            alphaK2     1.0;
-            alphaOmega1 0.5;
-            alphaOmega2 0.856;
-            beta1       0.075;
-            beta2       0.0828;
-            betaStar    0.09;
-            gamma1      5/9;
-            gamma2      0.44;
-            a1          0.31;
-            b1          1.0;
-            c1          10.0;
-            F3          no;
-        }
-    @endverbatim
-
-SourceFiles
-    immersedBoundaryKOmegaSST.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef immersedBoundaryKOmegaSST_H
-#define immersedBoundaryKOmegaSST_H
-
-#include "RASModel.H"
-#include "wallDist.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace incompressible
-{
-namespace RASModels
-{
-
-/*---------------------------------------------------------------------------*\
-                           Class kOmega Declaration
-\*---------------------------------------------------------------------------*/
-
-class immersedBoundaryKOmegaSST
-:
-    public RASModel
-{
-    // Private data
-
-        // Model coefficients
-            dimensionedScalar alphaK1_;
-            dimensionedScalar alphaK2_;
-
-            dimensionedScalar alphaOmega1_;
-            dimensionedScalar alphaOmega2_;
-
-            dimensionedScalar gamma1_;
-            dimensionedScalar gamma2_;
-
-            dimensionedScalar beta1_;
-            dimensionedScalar beta2_;
-
-            dimensionedScalar betaStar_;
-
-            dimensionedScalar a1_;
-            dimensionedScalar b1_;
-            dimensionedScalar c1_;
-
-            Switch F3_;
-
-
-        //- Wall distance field
-        //  Note: different to wall distance in parent RASModel
-        wallDist y_;
-
-
-        // Fields
-
-            volScalarField k_;
-            volScalarField omega_;
-            volScalarField nut_;
-
-
-    // Private member functions
-
-        tmp<volScalarField> F1(const volScalarField& CDkOmega) const;
-        tmp<volScalarField> F2() const;
-        tmp<volScalarField> F3() const;
-        tmp<volScalarField> F23() const;
-
-        tmp<volScalarField> blend
-        (
-            const volScalarField& F1,
-            const dimensionedScalar& psi1,
-            const dimensionedScalar& psi2
-        ) const
-        {
-            return F1*(psi1 - psi2) + psi2;
-        }
-
-        tmp<volScalarField> alphaK
-        (
-            const volScalarField& F1
-        ) const
-        {
-            return blend(F1, alphaK1_, alphaK2_);
-        }
-
-        tmp<volScalarField> alphaOmega
-        (
-            const volScalarField& F1
-        ) const
-        {
-            return blend(F1, alphaOmega1_, alphaOmega2_);
-        }
-
-        tmp<volScalarField> beta
-        (
-            const volScalarField& F1
-        ) const
-        {
-            return blend(F1, beta1_, beta2_);
-        }
-
-        tmp<volScalarField> gamma
-        (
-            const volScalarField& F1
-        ) const
-        {
-            return blend(F1, gamma1_, gamma2_);
-        }
-
-
-public:
-
-    //- Runtime type information
-    TypeName("immersedBoundaryKOmegaSST");
-
-
-    // Constructors
-
-        //- Construct from components
-        immersedBoundaryKOmegaSST
-        (
-            const volVectorField& U,
-            const surfaceScalarField& phi,
-            transportModel& transport,
-            const word& turbulenceModelName = turbulenceModel::typeName,
-            const word& modelName = typeName
-        );
-
-
-    //- Destructor
-    virtual ~immersedBoundaryKOmegaSST()
-    {}
-
-
-    // Member Functions
-
-        //- Return the turbulence viscosity
-        virtual tmp<volScalarField> nut() const
-        {
-            return nut_;
-        }
-
-        //- Return the effective diffusivity for k
-        tmp<volScalarField> DkEff(const volScalarField& F1) const
-        {
-            return tmp<volScalarField>
-            (
-                new volScalarField("DkEff", alphaK(F1)*nut_ + nu())
-            );
-        }
-
-        //- Return the effective diffusivity for omega
-        tmp<volScalarField> DomegaEff(const volScalarField& F1) const
-        {
-            return tmp<volScalarField>
-            (
-                new volScalarField("DomegaEff", alphaOmega(F1)*nut_ + nu())
-            );
-        }
-
-        //- Return the turbulence kinetic energy
-        virtual tmp<volScalarField> k() const
-        {
-            return k_;
-        }
-
-        //- Return the turbulence specific dissipation rate
-        virtual tmp<volScalarField> omega() const
-        {
-            return omega_;
-        }
-
-        //- Return the turbulence kinetic energy dissipation rate
-        virtual tmp<volScalarField> epsilon() const
-        {
-            return tmp<volScalarField>
-            (
-                new volScalarField
-                (
-                    IOobject
-                    (
-                        "epsilon",
-                        mesh_.time().timeName(),
-                        mesh_
-                    ),
-                    betaStar_*k_*omega_,
-                    omega_.boundaryField().types()
-                )
-            );
-        }
-
-        //- Return the Reynolds stress tensor
-        virtual tmp<volSymmTensorField> R() const;
-
-        //- Return the effective stress tensor including the laminar stress
-        virtual tmp<volSymmTensorField> devReff() const;
-
-        //- Return the source term for the momentum equation
-        virtual tmp<fvVectorMatrix> divDevReff() const;
-
-        //- Solve the turbulence equations and correct the turbulence viscosity
-        virtual void correct();
-
-        //- Read RASProperties dictionary
-        virtual bool read();
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace RASModels
-} // namespace incompressible
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C

@@ -48,6 +48,7 @@ License
 #include "processorPointPatch.H"
 #include "globalIndex.H"
 #include "meshTools.H"
+#include "labelIOField.H"
 #include "OFstream.H"
 #include "geomDecomp.H"
 #include "Random.H"
@@ -860,7 +861,7 @@ Foam::meshRefinement::meshRefinement
     meshCutter_
     (
         mesh,
-        labelIOList
+        labelIOField
         (
             IOobject
             (
@@ -872,9 +873,9 @@ Foam::meshRefinement::meshRefinement
                 IOobject::NO_WRITE,
                 false
             ),
-            labelList(mesh_.nCells(), 0)
+            labelField(mesh_.nCells(), 0)
         ),
-        labelIOList
+        labelIOField
         (
             IOobject
             (
@@ -886,7 +887,7 @@ Foam::meshRefinement::meshRefinement
                 IOobject::NO_WRITE,
                 false
             ),
-            labelList(mesh_.nPoints(), 0)
+            labelField(mesh_.nPoints(), 0)
         ),
         refinementHistory
         (

src/overset/oversetMesh/Make/files

@@ -4,6 +4,8 @@ oversetFringe/oversetFringe/oversetFringe.C
 oversetFringe/oversetFringe/newOversetFringe.C
 oversetFringe/manualFringe/manualFringe.C
 oversetFringe/faceCellsFringe/faceCellsFringe.C
+oversetFringe/donorBasedLayeredOverlapFringe/donorBasedLayeredOverlapFringe.C
+oversetFringe/cuttingPatchFringe/cuttingPatchFringe.C
 oversetFringe/overlapFringe/overlapFringe/overlapFringe.C
 oversetFringe/overlapFringe/layeredOverlapFringe/layeredOverlapFringe.C
 oversetFringe/overlapFringe/adaptiveOverlapFringe/adaptiveOverlapFringe.C

src/overset/oversetMesh/Make/options

@@ -2,12 +2,10 @@ EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/surfMesh/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude \
-    -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude
+    -I$(LIB_SRC)/sampling/lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \
     -lmeshTools \
     -lsurfMesh \
-    -lsampling \
-    -ldynamicMesh
+    -lsampling

src/overset/oversetMesh/include/oversetAlphaCourantNo.H

@@ -38,21 +38,28 @@ scalar maxAlphaCo
 scalar alphaCoNum = 0.0;
 scalar meanAlphaCoNum = 0.0;
 
-surfaceScalarField alpha1f =
+const surfaceScalarField alpha1f =
     fvc::interpolate(min(max(alpha1, scalar(0)), scalar(1)));
 
-const dimensionedScalar alphaOffset("alphaOffset", dimless, dAlpha);
+const dimensionedScalar alphaOffset
+(
+    "alphaOffset",
+    dimless,
+    runTime.controlDict().lookupOrDefault("dAlpha", 0.01)
+);
 
 if (mesh.nInternalFaces())
 {
-    surfaceScalarField magAlphaPhi
+    const oversetMesh& om = oversetMesh::New(mesh);
+
+    const surfaceScalarField magAlphaPhi
     (
         pos(alpha1f - alphaOffset)*
         pos(scalar(1) - alphaOffset - alpha1f)*
-        mag(faceOversetMask*phi)
+        mag(om.sGamma()*phi)
     );
 
-    surfaceScalarField SfUfbyDelta =
+    const surfaceScalarField SfUfbyDelta =
         mesh.surfaceInterpolation::deltaCoeffs()*magAlphaPhi;
 
     const scalar deltaT = runTime.deltaT().value();

src/overset/oversetMesh/include/oversetCourantNo.H

@@ -36,9 +36,11 @@ scalar velMag = 0.0;
 
 if (mesh.nInternalFaces())
 {
-    surfaceScalarField magPhi = mag(faceOversetMask*phi);
+    const oversetMesh& om = oversetMesh::New(mesh);
 
-    surfaceScalarField SfUfbyDelta =
+    const surfaceScalarField magPhi = mag(om.sGamma()*phi);
+
+    const surfaceScalarField SfUfbyDelta =
         mesh.surfaceInterpolation::deltaCoeffs()*magPhi;
 
     CoNum = max(SfUfbyDelta/mesh.magSf())

src/overset/oversetMesh/oversetAdjustPhi/oversetAdjustPhi.C

@@ -59,7 +59,7 @@ void Foam::oversetAdjustPhi
     const unallocLabelList& neighbour = mesh.neighbour();
 
     // Get region split to identify separate mesh components
-    const labelList& regionID = om.regionID();
+    const scalarField& regionID = om.regionID().internalField();
 
     // Sum up incoming and outgoing flux
     scalar fringeIn = 0;
@@ -80,8 +80,12 @@ void Foam::oversetAdjustPhi
         {
             // Internal face
 
-            // Debug check
-            if (regionID[owner[curFace]] != regionID[neighbour[curFace]])
+            // Debug check. Note: use notEqual since we are comparing two
+            // scalars
+            if
+            (
+                notEqual(regionID[owner[curFace]], regionID[neighbour[curFace]])
+            )
             {
                 FatalErrorIn
                 (

src/overset/oversetMesh/oversetAdjustPhi/regionWiseOversetAdjustPhi.C

@@ -60,7 +60,7 @@ void Foam::regionWiseOversetAdjustPhi
     const unallocLabelList& neighbour = mesh.neighbour();
 
     // Get region split to identify separate mesh components
-    const labelList& regionID = om.regionID();
+    const scalarField& regionID = om.regionID().internalField();
 
     // Incoming and outgoing region fluxes
     scalarField regionIn(om.regions().size(), 0);
@@ -94,7 +94,7 @@ void Foam::regionWiseOversetAdjustPhi
             // scaled
             forAll (phip, i)
             {
-                // Get current region index
+                // Get current region index. Note: conversion to label
                 const label curRegion = regionID[fc[i]];
 
                 if (phip[i] < 0.0)
@@ -124,11 +124,12 @@ void Foam::regionWiseOversetAdjustPhi
         {
             // Internal face
 
-            // Get region index
+            // Get region index. Note: conversion to label
             const label curRegion = regionID[owner[curFace]];
 
-            // Check whether owner and neighbour belong to the same region
-            if (curRegion != regionID[neighbour[curFace]])
+            // Check whether owner and neighbour belong to the same region.
+            // Note: use notEqual function since we are comparing two scalars
+            if (notEqual(curRegion, regionID[neighbour[curFace]]))
             {
                 FatalErrorIn
                 (
@@ -210,7 +211,7 @@ void Foam::regionWiseOversetAdjustPhi
                 {
                     // Processor patch, master side
 
-                    // Get region index
+                    // Get region index. Note: conversion to label
                     const label curRegion =
                         regionID[mesh.boundary()[patchI].faceCells()[faceI]];
 
@@ -321,7 +322,7 @@ void Foam::regionWiseOversetAdjustPhi
         {
             // Internal face
 
-            // Get region index
+            // Get region index. Note: conversion to label
             const label curRegion = regionID[owner[curFace]];
 
             // Get reference to the flux for scaling
@@ -358,7 +359,7 @@ void Foam::regionWiseOversetAdjustPhi
 
             if (procPatch.master()) // Owner side
             {
-                // Get region index
+                // Get region index. Note: conversion to label
                 const label curRegion =
                     regionID[mesh.boundary()[patchI].faceCells()[faceI]];
 
@@ -374,7 +375,7 @@ void Foam::regionWiseOversetAdjustPhi
             }
             else // Neighbouring processor side
             {
-                // Get region index
+                // Get region index. Note: conversion to label
                 const label curRegion =
                     regionID[mesh.boundary()[patchI].faceCells()[faceI]];
 

src/overset/oversetMesh/oversetFringe/cuttingPatchFringe/cuttingPatchFringe.C

@@ -0,0 +1,649 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.1
+    \\  /    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 "cuttingPatchFringe.H"
+#include "oversetMesh.H"
+#include "oversetRegion.H"
+#include "polyPatchID.H"
+#include "addToRunTimeSelectionTable.H"
+#include "syncTools.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(cuttingPatchFringe, 0);
+    addToRunTimeSelectionTable
+    (
+        oversetFringe,
+        cuttingPatchFringe,
+        dictionary
+    );
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::cuttingPatchFringe::calcAddressing() const
+{
+    // Make sure that either acceptorsPtr is unnalocated or if it is allocated,
+    // that it is empty
+    if (acceptorsPtr_ && !acceptorsPtr_->empty())
+    {
+        FatalErrorIn
+        (
+            "void Foam::cuttingPatchFringe::calcAddressing() const"
+        )   << "Addressing already calculated"
+            << abort(FatalError);
+    }
+
+    // Get polyMesh
+    const polyMesh& mesh = this->mesh();
+
+    // Collect all cutting patches
+    labelHashSet patchIDs(cuttingPatchNames_.size());
+
+    forAll (cuttingPatchNames_, nameI)
+    {
+        // Get polyPatchID and check if valid
+        const polyPatchID cutPatch
+        (
+            cuttingPatchNames_[nameI],
+            mesh.boundaryMesh()
+        );
+
+        if (!cutPatch.active())
+        {
+            FatalErrorIn
+            (
+                "void cuttingPatchFringe::calcAddressing const"
+            )   << "Cutting patch " << cuttingPatchNames_[nameI]
+                << " cannot be found."
+                << abort(FatalError);
+        }
+
+        // Store patch ID in the set
+        patchIDs.insert(cutPatch.index());
+    }
+
+    if (debug)
+    {
+        Info<< "Starting cutting patch fringe assembly..." << endl;
+    }
+
+    // Note: similar code as in oversetRegion::calcHoleTriMesh. Consider
+    // refactoring. VV, 20/May/2019
+
+    // Make and invert local triSurface
+    triFaceList triFaces;
+    pointField triPoints;
+
+    // Memory management
+    {
+        triSurface ts = triSurfaceTools::triangulate
+        (
+            mesh.boundaryMesh(),
+            patchIDs
+        );
+
+        // Clean mutiple points and zero-sized triangles
+        ts.cleanup(false);
+
+        triFaces.setSize(ts.size());
+        triPoints = ts.points();
+
+        forAll (ts, tsI)
+        {
+            // Bugfix: no need to reverse the face because the normals point in
+            // the correct direction already. VV, 20/May/2019.
+            triFaces[tsI] = ts[tsI];
+        }
+    }
+
+    if (Pstream::parRun())
+    {
+        // Combine all faces and points into a single list
+
+        List<triFaceList> allTriFaces(Pstream::nProcs());
+        List<pointField> allTriPoints(Pstream::nProcs());
+
+        allTriFaces[Pstream::myProcNo()] = triFaces;
+        allTriPoints[Pstream::myProcNo()] = triPoints;
+
+        Pstream::gatherList(allTriFaces);
+        Pstream::scatterList(allTriFaces);
+
+        Pstream::gatherList(allTriPoints);
+        Pstream::scatterList(allTriPoints);
+
+        // Re-pack points and faces
+
+        label nTris = 0;
+        label nPoints = 0;
+
+        forAll (allTriFaces, procI)
+        {
+            nTris += allTriFaces[procI].size();
+            nPoints += allTriPoints[procI].size();
+        }
+
+        // Pack points
+        triPoints.setSize(nPoints);
+
+        // Prepare point renumbering array
+        labelListList renumberPoints(Pstream::nProcs());
+
+        nPoints = 0;
+
+        forAll (allTriPoints, procI)
+        {
+            const pointField& ptp = allTriPoints[procI];
+
+            renumberPoints[procI].setSize(ptp.size());
+
+            labelList& procRenumberPoints = renumberPoints[procI];
+
+            forAll (ptp, ptpI)
+            {
+                triPoints[nPoints] = ptp[ptpI];
+                procRenumberPoints[ptpI] = nPoints;
+
+                nPoints++;
+            }
+        }
+
+        // Pack triangles and renumber into complete points on the fly
+        triFaces.setSize(nTris);
+
+        nTris = 0;
+
+        forAll (allTriFaces, procI)
+        {
+            const triFaceList& ptf = allTriFaces[procI];
+
+            const labelList& procRenumberPoints = renumberPoints[procI];
+
+            forAll (ptf, ptfI)
+            {
+                const triFace& procFace = ptf[ptfI];
+
+                triFace& renumberFace = triFaces[nTris];
+
+                forAll (renumberFace, rfI)
+                {
+                    renumberFace[rfI] =
+                        procRenumberPoints[procFace[rfI]];
+                }
+
+                nTris++;
+            }
+        }
+    }
+
+    // Make a complete triSurface from local data
+    triSurface patchTriMesh(triFaces, triPoints);
+
+    // Clean up duplicate points and zero sized triangles
+    patchTriMesh.cleanup(false);
+
+    // Get this region
+    const oversetRegion& myRegion = this->region();
+
+    // Debug: write down the tri mesh
+    if (Pstream::master())
+    {
+        patchTriMesh.write
+        (
+            word
+            (
+                "patchTriSurface_region" + myRegion.name() + ".vtk"
+            )
+        );
+    }
+
+    // Create the tri surface search object
+    triSurfaceSearch patchTriSearch(patchTriMesh);
+
+    // Get cells in this region
+    const labelList& myRegionCells = myRegion.regionCells();
+
+    // Get cell centres for inside-outside search using search object
+    vectorField myCC(myRegionCells.size());
+
+    // Cell centres from polyMesh
+    const vectorField& cc = mesh.cellCentres();
+
+    forAll (myCC, i)
+    {
+        myCC[i] = cc[myRegionCells[i]];
+    }
+
+    // Inside mask: all cells within search object will be marked
+    boolList insideMask(mesh.nCells(), false);
+
+    // Get inside cells for cells in my region only
+    boolList myRegionInsideMask = patchTriSearch.calcInside(myCC);
+
+    // Note: insideMask has the size of all mesh cells and
+    // myRegionInsideMask has the size of cells in this region
+    forAll (myRegionInsideMask, i)
+    {
+        insideMask[myRegionCells[i]] = myRegionInsideMask[i];
+    }
+
+    // Make sure that the cut holes for this region are also properly marked as
+    // "inside". This may not be the case automatically for e.g. simulations
+    // with appendages
+    const labelList& cutRegionHoles = myRegion.cutHoles();
+    forAll (cutRegionHoles, i)
+    {
+        insideMask[cutRegionHoles[i]] = true;
+    }
+
+    // Get necessary mesh data (from polyMesh/primitiveMesh)
+    const cellList& meshCells = mesh.cells();
+    const unallocLabelList& owner = mesh.faceOwner();
+    const unallocLabelList& neighbour = mesh.faceNeighbour();
+
+    // Bool list for collecting faces with at least one unmarked
+    // cell (to determine the acceptors for the first iteration)
+    boolList hasUnmarkedCell(mesh.nFaces(), false);
+
+    // Loop through all cells
+    forAll (insideMask, cellI)
+    {
+        if (!insideMask[cellI])
+        {
+            // This cell is not inside (it is unmarked). Loop through
+            // its faces and set the flag
+            const cell& cFaces = meshCells[cellI];
+
+            forAll (cFaces, i)
+            {
+                // Set the mark for this global face
+                hasUnmarkedCell[cFaces[i]] = true;
+            }
+        }
+    }
+
+    // Sync the face list across processor boundaries
+    syncTools::syncFaceList
+    (
+        mesh,
+        hasUnmarkedCell,
+        orEqOp<bool>(),
+        true
+    );
+
+    // Mark-up for all inside faces
+    boolList insideFaceMask(mesh.nFaces(), false);
+
+    // Collect all acceptors for the first iteration (the cells that
+    // have at least one neighbour cell that is not marked)
+    labelHashSet acceptors(myRegionCells.size()/10);
+
+    // Loop again through all cells and collect marked ones into
+    // acceptors or holes, depending on whether they have unmarked cell
+    // as a neighbour (indicating an acceptor)
+    forAll (insideMask, cellI)
+    {
+        if (insideMask[cellI])
+        {
+            // This cell is inside the covered region
+            const cell& cFaces = meshCells[cellI];
+
+            forAll (cFaces, i)
+            {
+                // Get global face index
+                const label& faceI = cFaces[i];
+
+                // Check whether this neighbour is unmarked
+                if (hasUnmarkedCell[faceI])
+                {
+                    // This cell has unmarked neighbour, collect it into
+                    // the acceptor list
+                    acceptors.insert(cellI);
+
+                    // This cell is no longer "inside cell"
+                    insideMask[cellI] = false;
+                }
+                else
+                {
+                    // This is an "inside" face, mark it
+                    insideFaceMask[faceI] = true;
+                }
+            } // End for all faces
+        } // End if cell is inside
+    } // End for all cells
+
+    // Note: insideFaceMask already synced across processors because it relies
+    // on hasUnmarkedCell list, which has been synced just above
+
+    // Hash set containing new acceptors (for successive iterations)
+    labelHashSet newAcceptors(acceptors.size());
+
+    // Now that we have the initial set of acceptors (and holes), loop
+    // nLayers away from initial donors
+    for (label i = 0; i < nLayers_; ++i)
+    {
+        // Face markup for propagation
+        boolList propagateFace(mesh.nFaces(), false);
+
+        // Loop through all acceptors and mark faces that are around hole
+        // cells. This way, we make sure that we go towards the correct,
+        // inside direction
+        forAllConstIter (labelHashSet, acceptors, iter)
+        {
+            // Get the cell index and the cell
+            const label& cellI = iter.key();
+            const cell& cFaces = meshCells[cellI];
+
+            // Loop through all faces of the cell
+            forAll (cFaces, i)
+            {
+                // Get face index (global)
+                const label& faceI = cFaces[i];
+
+                if (insideFaceMask[faceI])
+                {
+                    // This is a hole face, we are moving in the right
+                    // direction. Mark the face for propagation
+                    propagateFace[faceI] = true;
+                }
+            } // End for all faces of the cell
+        } // End for all donor cells
+
+        // Sync the face list across processor boundaries
+        syncTools::syncFaceList
+        (
+            mesh,
+            propagateFace,
+            orEqOp<bool>(),
+            false
+        );
+
+        // Loop through all internal faces and append acceptors
+        for (label faceI = 0; faceI < mesh.nInternalFaces(); ++faceI)
+        {
+            if (propagateFace[faceI])
+            {
+                // Face is marked, select owner or neighbour
+                const label& own = owner[faceI];
+                const label& nei = neighbour[faceI];
+
+                // Either owner or neighbour may be hole, not both
+                if (insideMask[own])
+                {
+                    // Owner cell is a hole, insert it
+                    newAcceptors.insert(own);
+
+                    // Update hole mask
+                    insideMask[own] = false;
+                }
+                else if (insideMask[nei])
+                {
+                    // Neighbour cell is a hole, insert it
+                    newAcceptors.insert(nei);
+
+                    // Update hole mask
+                    insideMask[nei] = false;
+                }
+
+                // Also update hole face mask for next iteration
+                insideFaceMask[faceI] = false;
+            }
+        }
+
+        // Loop through boundary faces
+        for
+        (
+            label faceI = mesh.nInternalFaces();
+            faceI < mesh.nFaces();
+          ++faceI
+        )
+        {
+            if (propagateFace[faceI])
+            {
+                // Face is marked, select owner if this is the right
+                // side. Neighbour handled on the other side
+                const label& own = owner[faceI];
+
+                if (insideMask[own])
+                {
+                    // Face cell is a hole, insert it
+                    newAcceptors.insert(own);
+
+                    // Update hole mask
+                    insideMask[own] = false;
+                }
+
+                // Also update hole face mask for next iteration
+                insideFaceMask[faceI] = false;
+            }
+        }
+
+        // Transfer newAcceptors into acceptors for next iteration or
+        // for final assembly. Resize newAcceptors accordingly
+        acceptors.transfer(newAcceptors);
+        newAcceptors.resize(acceptors.size());
+
+    } // End for specified number of layers
+
+    // At this point, we have the final set of acceptors and we marked
+    // all cells that should be holes. Collect them into the list
+    dynamicLabelList fringeHoles(myRegionCells.size()/10);
+
+    forAll (insideMask, cellI)
+    {
+        if (insideMask[cellI])
+        {
+            fringeHoles.append(cellI);
+        }
+    }
+
+    // Set acceptors and holes from the data for all regions
+    acceptorsPtr_ = new labelList(acceptors.sortedToc());
+    fringeHolesPtr_ = new labelList(fringeHoles.xfer());
+
+    if (debug)
+    {
+        Info<< "In cuttingPatchFringe::calcAddressing() const" << nl
+            << "Found " << acceptorsPtr_->size() << " acceptors." << nl
+            << "Found " << fringeHolesPtr_->size() << " fringe holes." << endl;
+    }
+}
+
+
+void Foam::cuttingPatchFringe::clearAddressing() const
+{
+    deleteDemandDrivenData(fringeHolesPtr_);
+    deleteDemandDrivenData(acceptorsPtr_);
+    deleteDemandDrivenData(finalDonorAcceptorsPtr_);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from dictionary
+Foam::cuttingPatchFringe::cuttingPatchFringe
+(
+    const fvMesh& mesh,
+    const oversetRegion& region,
+    const dictionary& dict
+)
+:
+    oversetFringe(mesh, region, dict),
+    cuttingPatchNames_(dict.lookup("cuttingPatches")),
+    nLayers_(readLabel(dict.lookup("nLayers"))),
+    fringeHolesPtr_(nullptr),
+    acceptorsPtr_(nullptr),
+    finalDonorAcceptorsPtr_(nullptr)
+{
+    // Sanity check number of layers: must be greater than 0
+    if (nLayers_ < 1)
+    {
+        FatalIOErrorIn
+        (
+            "cuttingPatchFringe::"
+            "cuttingPatchFringe\n"
+            "(\n"
+            "    const fvMesh& mesh,\n"
+            "    const oversetRegion& region,\n"
+            "    const dictionary& dict\n"
+            ")",
+            dict
+        )   << "Invalid number of layers specified, nLayers = " << nLayers_
+            << nl
+            << "The number should be greater than 0."
+            << exit(FatalError);
+    }
+
+    // Preferably, the number of layers should be at least 2
+    if (nLayers_ == 1)
+    {
+        WarningIn
+        (
+            "cuttingPatchFringe::"
+            "cuttingPatchFringe\n"
+            "(\n"
+            "    const fvMesh& mesh,\n"
+            "    const oversetRegion& region,\n"
+            "    const dictionary& dict\n"
+            ")"
+        )   << "You have specified nLayers = " << nLayers_
+            << nl
+            << "We strongly advise to use at least 2 layers to avoid" << nl
+            << "possibility of having acceptors that cannot find decent" << nl
+            << "donors on the other side."
+            << endl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::cuttingPatchFringe::~cuttingPatchFringe()
+{
+    clearAddressing();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::cuttingPatchFringe::updateIteration
+(
+    donorAcceptorList& donorAcceptorRegionData
+) const
+{
+    // If the donorAcceptor list has been allocated, something went wrong with
+    // the iteration procedure (not-updated flag): this function has been called
+    // more than once, which should not happen for cuttingPatchFringe
+    if (finalDonorAcceptorsPtr_)
+    {
+        FatalErrorIn
+        (
+            "cuttingPatchFringe::"
+            "updateIteration(donorAcceptorList&) const"
+        )   << "finalDonorAcceptorPtr_ already allocated. Something went "
+            << "wrong with the iteration procedure (flag was not updated)."
+            << nl
+            << "This should not happen for cuttingPatchFringe."
+            << abort(FatalError);
+    }
+
+    // Allocate the list by reusing the argument list
+    finalDonorAcceptorsPtr_ = new donorAcceptorList
+    (
+        donorAcceptorRegionData,
+        true
+    );
+
+    // Set the flag to true and return
+    updateSuitableOverlapFlag(true);
+
+    return foundSuitableOverlap();
+}
+
+
+const Foam::labelList& Foam::cuttingPatchFringe::fringeHoles() const
+{
+    if (!fringeHolesPtr_)
+    {
+        calcAddressing();
+    }
+
+    return *fringeHolesPtr_;
+}
+
+
+const Foam::labelList& Foam::cuttingPatchFringe::candidateAcceptors() const
+{
+    if (!acceptorsPtr_)
+    {
+        calcAddressing();
+    }
+
+    return *acceptorsPtr_;
+}
+
+
+Foam::donorAcceptorList&
+Foam::cuttingPatchFringe::finalDonorAcceptors() const
+{
+    if (!finalDonorAcceptorsPtr_)
+    {
+        FatalErrorIn("cuttingPatchFringe::finalDonorAcceptors()")
+            << "finalDonorAcceptorPtr_ not allocated. Make sure you have"
+            << " called cuttingPatchFringe::updateIteration() before"
+            << " asking for final set of donor/acceptor pairs."
+            << abort(FatalError);
+    }
+
+    if (!foundSuitableOverlap())
+    {
+        FatalErrorIn("cuttingPatchFringe::finalDonorAcceptors()")
+            << "Attemted to access finalDonorAcceptors but suitable overlap "
+            << "has not been found. This is not allowed. "
+            << abort(FatalError);
+    }
+
+    return *finalDonorAcceptorsPtr_;
+}
+
+
+void Foam::cuttingPatchFringe::update() const
+{
+    Info<< "cuttingPatchFringe::update() const" << endl;
+
+    // Clear out
+    clearAddressing();
+
+    // Set flag to false and clear final donor/acceptors only
+    deleteDemandDrivenData(finalDonorAcceptorsPtr_);
+    updateSuitableOverlapFlag(false);
+}
+
+
+// ************************************************************************* //

src/overset/oversetMesh/oversetFringe/cuttingPatchFringe/cuttingPatchFringe.H

@@ -0,0 +1,151 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.1
+    \\  /    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
+    cuttingPatchFringe
+
+Description
+    Fringe algorithm based on collecting all the cells cut by its donor region
+    patch (assuming faceCells algorithm is used). Inside cells are marked and
+    moved nLayers towards the interior to define the acceptors as robustly as
+    possible.
+
+Author
+    Vuko Vukcevic, Wikki Ltd.  All rights reserved.
+
+SourceFiles
+    cuttingPatchFringe.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef cuttingPatchFringe_H
+#define cuttingPatchFringe_H
+
+#include "oversetFringe.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                   Class cuttingPatchFringe Declaration
+\*---------------------------------------------------------------------------*/
+
+class cuttingPatchFringe
+:
+    public oversetFringe
+{
+    // Private data
+
+        //- Names of the cuttingPatches
+        wordList cuttingPatchNames_;
+
+        //- How many layers to move away from connected region donors to define
+        //  acceptor (and holes)
+        label nLayers_;
+
+        //- Fringe hole cells
+        mutable labelList* fringeHolesPtr_;
+
+        //- Acceptor cells
+        mutable labelList* acceptorsPtr_;
+
+        //- Final donor/acceptor pairs for this region (fringe)
+        mutable donorAcceptorList* finalDonorAcceptorsPtr_;
+
+
+    // Private Member Functions
+
+        //- Calculate hole and acceptor addressing
+        void calcAddressing() const;
+
+        //- Clear addressing
+        void clearAddressing() const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("cuttingPatch");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        cuttingPatchFringe
+        (
+            const fvMesh& mesh,
+            const oversetRegion& region,
+            const dictionary& dict
+        );
+
+        //- Disallow default bitwise copy construct
+        cuttingPatchFringe
+        (
+            const cuttingPatchFringe&
+        ) = delete;
+
+        //- Disallow default bitwise assignment
+        void operator=(const cuttingPatchFringe&) = delete;
+
+
+    // Destructor
+
+        virtual ~cuttingPatchFringe();
+
+
+    // Member Functions
+
+        //- Update iteration. Note: invalidates parameter
+        virtual bool updateIteration
+        (
+            donorAcceptorList& donorAcceptorRegionData
+        ) const;
+
+        //- Return list of deactivated (hole) cells
+        //  Fringe hole cells are collected in addition to geometric hole
+        //  cells, which fall outside of all donor regions
+        virtual const labelList& fringeHoles() const;
+
+        //- Return list of acceptor cells
+        virtual const labelList& candidateAcceptors() const;
+
+        //- Return list of final donor acceptor pairs. Note: caller may
+        //  invalidate finalDonorAcceptorsPtr_ for optimisation purposes
+        virtual donorAcceptorList& finalDonorAcceptors() const;
+
+        //- Update the fringe
+        virtual void update() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/overset/oversetMesh/oversetFringe/donorBasedLayeredOverlapFringe/donorBasedLayeredOverlapFringe.C

@@ -0,0 +1,887 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.1
+    \\  /    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 "donorBasedLayeredOverlapFringe.H"
+#include "oversetMesh.H"
+#include "oversetRegion.H"
+#include "faceCellsFringe.H"
+#include "oversetRegion.H"
+#include "addToRunTimeSelectionTable.H"
+#include "syncTools.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(donorBasedLayeredOverlapFringe, 0);
+    addToRunTimeSelectionTable
+    (
+        oversetFringe,
+        donorBasedLayeredOverlapFringe,
+        dictionary
+    );
+}
+
+
+const Foam::debug::tolerancesSwitch
+Foam::donorBasedLayeredOverlapFringe::distTol_
+(
+    "donorBasedLayeredOverlapDistanceTolerance",
+    0.0
+);
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::donorBasedLayeredOverlapFringe::init() const
+{
+    // Set size of the list containing IDs
+    connectedRegionIDs_.setSize(connectedRegionNames_.size());
+
+    // Get list of all overset regions
+    const PtrList<oversetRegion>& allRegions =
+        this->region().overset().regions();
+
+    // Create list of all region names for easy lookup
+    wordList allRegionNames(allRegions.size());
+    forAll (allRegionNames, arI)
+    {
+        allRegionNames[arI] = allRegions[arI].name();
+    }
+
+    // Loop through all regions and check whether the overlap has been found
+    forAll (connectedRegionNames_, crI)
+    {
+        // Get name of this connected region
+        const word& crName = connectedRegionNames_[crI];
+
+        // Find this region in the list of all regions
+        const label regionID = findIndex(allRegionNames, crName);
+
+        if (regionID == -1)
+        {
+            FatalErrorIn("void donorBasedLayeredOverlapFringe::init() const")
+                << "Region " << crName << " not found in list of regions."
+                << "List of overset regions: " << allRegionNames
+                << abort(FatalError);
+        }
+
+        // Collect the region index in the list
+        connectedRegionIDs_[crI] = regionID;
+
+        // Sanity check: if the specified connected donor region has more than 1
+        // donor regions, this fringe algorithm is attempted to be used for
+        // something that's not intended. Issue an error
+        if (allRegions[regionID].donorRegions().size() != 1)
+        {
+            FatalErrorIn("void donorBasedLayeredOverlapFringe::init() const")
+                << "Region " << crName << " specified as connected region, but"
+                << " that region has "
+                << allRegions[regionID].donorRegions().size()
+                << " donor regions."
+                << abort(FatalError);
+        }
+
+        // Sanity check whether the donor region of connected region is actually
+        // this region
+        if (allRegions[regionID].donorRegions()[0] != this->region().index())
+        {
+            FatalErrorIn("void donorBasedLayeredOverlapFringe::init() const")
+                << "The donor region of region " << crName
+                << " should be only region " << this->region().name()
+                << abort(FatalError);
+        }
+    }
+
+    // Sanity check: number of (optionally) specified centre points must be
+    // equal to the number of connected regions
+    if
+    (
+        !regionCentrePoints_.empty()
+     && (regionCentrePoints_.size() != connectedRegionIDs_.size())
+    )
+    {
+        FatalErrorIn("void donorBasedLayeredOverlapFringe::init() const")
+            << "You have specified "
+            << regionCentrePoints_.size()
+            << " regionCentrePoints, while specifying "
+            << connectedRegionIDs_.size()
+            << " connectedRegions."
+            << nl
+            << "If you'd like to avoid using automatic centre point detection,"
+            << " make sure to specify centre points for all connected regions."
+            << abort(FatalError);
+    }
+}
+
+
+void Foam::donorBasedLayeredOverlapFringe::calcAddressing() const
+{
+    // Make sure that either acceptorsPtr is unnalocated or if it is allocated,
+    // that it is empty
+    if (acceptorsPtr_ && !acceptorsPtr_->empty())
+    {
+        FatalErrorIn
+        (
+            "void Foam::donorBasedLayeredOverlapFringe::calcAddressing() const"
+        )   << "Addressing already calculated"
+            << abort(FatalError);
+    }
+
+    if (!isInitialized_)
+    {
+        // This is the first call, initialize the data and set flag to true
+        init();
+        isInitialized_ = true;
+    }
+
+    // Get list of all overset regions
+    const PtrList<oversetRegion>& allRegions =
+        this->region().overset().regions();
+
+    // Loop through all connected regions and check whether the fringe overlap
+    // has been found for all of them
+    bool allFringesReady = true;
+    forAll (connectedRegionIDs_, crI)
+    {
+        // Get ID of this region
+        const label& regionID = connectedRegionIDs_[crI];
+
+        // Get the overset region
+        const oversetRegion& region = allRegions[regionID];
+
+        // Get the fringe algorithm from the region
+        const oversetFringe& fringe = region.fringe();
+
+        // If this is not faceCells fringe, issue a Warning. This fringe
+        // selection algorithm is intended to work only with faceCells fringe on
+        // the other side. VV, 9/Apr/2019
+        if (!isA<faceCellsFringe>(fringe))
+        {
+            WarningIn
+            (
+                "void Foam::donorBasedLayeredOverlapFringe::"
+                "updateIteration(donorAcceptorList&) const"
+            )   << "donorBasedLayeredOverlap fringe is designed to work"
+                << " with faceCells fringe as a connected region fringe."
+                << nl
+                << "Connected overset region " << region.name()
+                << " has " << fringe.type() << " fringe type. "
+                << nl
+                << "Proceed with care!"
+                << endl;
+        }
+
+        // Update flag collecting whether all connected regions found the
+        // overlap
+        allFringesReady &= fringe.foundSuitableOverlap();
+    }
+
+    // Sets containing all acceptors and all holes for all connected regions
+    const polyMesh& mesh = this->mesh();
+    labelHashSet allAcceptors(0.02*mesh.nCells());
+    labelHashSet allFringeHoles(0.02*mesh.nCells());
+
+    // Communicate state across processors
+    reduce(allFringesReady, andOp<bool>());
+
+    if (allFringesReady)
+    {
+        if (debug)
+        {
+            Info<< "All dependent fringes are ready."
+                << " Starting donor based layered overlap assembly..." << endl;
+        }
+
+        // Loop through connected regions
+        forAll (connectedRegionIDs_, crI)
+        {
+            // Get ID of this region
+            const label& regionID = connectedRegionIDs_[crI];
+
+            // Get fringe of the connected region
+            const oversetFringe& fringe = allRegions[regionID].fringe();
+
+            // The fringe should be finalized, which means we may take a const
+            // reference to its final donor acceptors
+            const donorAcceptorList& crDonorAcceptorPairs =
+                fringe.finalDonorAcceptors();
+
+            // Need to gather/scatter the donor-acceptor pairs across all
+            // processors because these pairs only represent acceptors found on
+            // my processor. It would be possible to optimize this a bit using
+            // the mapDistribute tool, but I don't think it will represent a big
+            // overhead, especially since it's done once.
+            List<donorAcceptorList> allDonorAcceptorPairs(Pstream::nProcs());
+
+            // Fill in my part and communicate
+            allDonorAcceptorPairs[Pstream::myProcNo()] = crDonorAcceptorPairs;
+            Pstream::gatherList(allDonorAcceptorPairs);
+            Pstream::scatterList(allDonorAcceptorPairs);
+
+            // Count approximate number of acceptors to guess the size for the
+            // hash set containing donors
+            label nAllAcceptors = 0;
+            forAll (allDonorAcceptorPairs, procI)
+            {
+                nAllAcceptors += allDonorAcceptorPairs[procI].size();
+            }
+
+            // Hash set containing donors
+            labelHashSet donors(6*nAllAcceptors);
+
+            // Initialize centre of the donors of this connected region in order
+            // to search in a given direction
+            vector centrePoint(vector::zero);
+
+            // Loop through all processors
+            forAll (allDonorAcceptorPairs, procI)
+            {
+                // Get all donor/acceptor pairs found on this processor
+                const donorAcceptorList& procDonorAcceptorPairs =
+                    allDonorAcceptorPairs[procI];
+
+                // Loop through all donor/acceptors
+                forAll (procDonorAcceptorPairs, daI)
+                {
+                    // Get this donor/acceptor pair
+                    const donorAcceptor& daPair = procDonorAcceptorPairs[daI];
+
+                    // Check whether all donors have been found
+                    if (!daPair.donorFound())
+                    {
+                        FatalErrorIn
+                        (
+                            "donorBasedLayeredOverlapFringe::"
+                            "updateIteration(donorAcceptorList&) const"
+                        )   << "Donor not found for donor/acceptor pair " << daI
+                            << nl
+                            << "Donor/acceptor data: " << daPair
+                            << nl
+                            << "In connected region: "
+                            << allRegions[regionID].name()
+                            << abort(FatalError);
+                    }
+
+                    // Mark donors on my processor from this connected region.
+                    // Note that the check has been made in constructor to make
+                    // sure that this region is the only donor region for the
+                    // connected region
+                    if (daPair.donorProcNo() == Pstream::myProcNo())
+                    {
+                        // Get donor index
+                        const label& dI = daPair.donorCell();
+
+                        // Insert donor into the hash set
+                        if (donors.insert(dI))
+                        {
+                            // Donor has been inserted (not previously found in
+                            // the hash set), add donor point to centre point
+                            // (the centre point will be calculated later on as
+                            // arithmetic mean)
+                            centrePoint += daPair.donorPoint();
+                        }
+
+                        // Loop through extended donor cells
+                        const donorAcceptor::DynamicLabelList& extDonors =
+                            daPair.extendedDonorCells();
+                        const donorAcceptor::DynamicPointList& extDonorPoints =
+                            daPair.extendedDonorPoints();
+
+                        forAll (extDonors, i)
+                        {
+                            // Get donor index
+                            const label& edI = extDonors[i];
+
+                            // Inser extended donor into the hash set
+                            if (donors.insert(edI))
+                            {
+                                // Donor has been inserted (not previously found
+                                // in the hash set), add extended donor point as
+                                // well
+                                centrePoint += extDonorPoints[i];
+                            }
+                        } // End for all extended donors
+                    } // End if this donor is on my processor
+                } // End for all (master) donor cells
+            } // End for all processors
+
+            // Use the centre point as specified by the user if it was specified
+            // (if the regionCentrePoints_ list is not empty). This avoids
+            // parallel communication as well.
+            if (!regionCentrePoints_.empty())
+            {
+                // Use specified centre point, discarding the data we calculated
+                // above
+                centrePoint = regionCentrePoints_[crI];
+            }
+            else
+            {
+                // User did not specify centre points and the centre point holds
+                // the sum of all the points. Reduce centre point and divide it
+                // with global number of unique donors
+                reduce(centrePoint, sumOp<vector>());
+                centrePoint /= returnReduce(donors.size(), sumOp<label>());
+            }
+
+            if (debug)
+            {
+                Info<< "Centre point for donors for region "
+                    << allRegions[regionID].name() << " is : " << centrePoint
+                    << endl;
+            }
+
+            // We now have a collection of all donors for this connected region
+            // and the centre point to move to. Let's collect the acceptors
+            labelHashSet acceptors(donors.size()); // Reasonable size estimate
+
+            // Get necessary mesh data (from polyMesh/primitiveMesh)
+            const vectorField& cc = mesh.cellCentres();
+            const vectorField& fc = mesh.faceCentres();
+            const cellList& meshCells = mesh.cells();
+            const unallocLabelList& owner = mesh.faceOwner();
+            const unallocLabelList& neighbour = mesh.faceNeighbour();
+
+            // Get bounding box of this region for additional check when marking
+            // acceptors
+            const boundBox& bb = allRegions[regionID].globalBounds();
+
+            // Mark cells that are eligible to be acceptors (not donors)
+            boolList eligibleCells(mesh.nCells(), true);
+            forAllConstIter (labelHashSet, donors, iter)
+            {
+                eligibleCells[iter.key()] = false;
+            }
+
+            // Loop nLayers away from initial donors
+            for (label i = 0; i < nLayers_; ++i)
+            {
+                // Face markup for propagation
+                boolList propagateFace(mesh.nFaces(), false);
+
+                // Loop through all donors and mark faces that are pointing
+                // towards the centre point and have an eligible neighbour
+                forAllConstIter (labelHashSet, donors, iter)
+                {
+                    // Get the cell index and the cell
+                    const label& cellI = iter.key();
+                    const cell& cFaces = meshCells[cellI];
+
+                    // Get cell centre of this donor and calculate distance to
+                    // centre point
+                    const vector& donorCentre = cc[cellI];
+                    const scalar donorCentreToRegionCentreDist =
+                        mag(donorCentre - centrePoint);
+
+                    // Loop through all faces of the cell
+                    forAll (cFaces, i)
+                    {
+                        // Get face index (global)
+                        const label& faceI = cFaces[i];
+
+                        // Get face centre and calculate distance to centre
+                        // point
+                        const vector& faceCentre = fc[faceI];
+                        const scalar faceCentreToRegionCentreDist =
+                            mag(faceCentre - centrePoint);
+
+                        if
+                        (
+                            // First, distance based criterium
+                            (
+                                faceCentreToRegionCentreDist
+                              - donorCentreToRegionCentreDist
+                              < distTol_
+                            )
+                            // Second, bounding box based criterium
+                         && (bb.contains(faceCentre))
+                        )
+                        {
+                            // Face is closer to the centre point than cell: we
+                            // are moving in the right direction. Mark the face
+                            propagateFace[faceI] = true;
+                        }
+                    } // End for all faces of the cell
+                } // End for all donor cells
+
+                // Sync the face list across processor boundaries
+                syncTools::syncFaceList
+                (
+                    mesh,
+                    propagateFace,
+                    orOp<bool>(),
+                    false
+                );
+
+                // Loop through all faces and append acceptors
+                for (label faceI = 0; faceI < mesh.nInternalFaces(); ++faceI)
+                {
+                    if (propagateFace[faceI])
+                    {
+                        // Face is marked, select owner or neighbour
+                        const label& own = owner[faceI];
+                        const label& nei = neighbour[faceI];
+
+                        // Either owner or neighbour may be eligible, not both
+                        if (eligibleCells[own])
+                        {
+                            // Owner cell is not a donor, insert it
+                            acceptors.insert(own);
+
+                            // Mark as ineligible
+                            eligibleCells[own] = false;
+                        }
+                        else if (eligibleCells[nei])
+                        {
+                            // Neighbour cell is not a donor, insert it
+                            acceptors.insert(nei);
+
+                            // Mark as ineligible
+                            eligibleCells[nei] = false;
+                        }
+                    }
+                }
+
+                // Loop through boundary faces
+                for
+                (
+                    label faceI = mesh.nInternalFaces();
+                    faceI < mesh.nFaces();
+                  ++faceI
+                )
+                {
+                    if (propagateFace[faceI])
+                    {
+                        // Face is marked, select owner if this is the right
+                        // side. Neighbour handled on the other side
+                        const label& own = owner[faceI];
+
+                        if (eligibleCells[own])
+                        {
+                            // Face cell is not a donor, insert it
+                            acceptors.insert(own);
+
+                            // Mark as ineligible
+                            eligibleCells[own] = false;
+                        }
+                    }
+                }
+
+                // Special treatment for all non-final iterations
+                if (i < nLayers_ - 1)
+                {
+                    // This is not the last iteration, transfer acceptors into
+                    // donors
+                    donors.transfer(acceptors);
+
+                    // Resize acceptors list
+                    acceptors.resize(donors.size());
+                }
+            } // End for specified number of layers
+
+            // At this point, we have the final set of acceptors and we marked
+            // all cells that are ineligible (either donor or acceptor). The
+            // remaining thing to do is to mark the interior holes
+
+            // Create a hash set that will contain fringe holes
+            labelHashSet fringeHoles(10*acceptors.size());
+
+            // Collect holes until there are no holes to collect. Note: for the
+            // first iteration, we will start with acceptors, otherwise, we will
+            // start from all fringe holes
+            label nAddedHoles;
+            bool firstIteration = true;
+            do
+            {
+                // Reset number of newly added holes
+                nAddedHoles = 0;
+
+                // Face markup for propagation
+                boolList propagateFace(mesh.nFaces(), false);
+
+                // Get collection of cells from which to start the search:
+                // acceptors for the first iteration and fringe holes otherwise
+                labelHashSet* curSetPtr = nullptr;
+                if (firstIteration)
+                {
+                    // Point current set to acceptors and switch off the flag
+                    curSetPtr = &acceptors;
+                    firstIteration = false;
+                }
+                else
+                {
+                    // Point current set to fringe holes
+                    curSetPtr = &fringeHoles;
+                }
+                const labelHashSet& curSet = *curSetPtr;
+
+                // Loop through all cells in the set and mark faces that are
+                // pointing towards the centre point and have eligible neighbour
+                // (not an acceptor or donor).
+                forAllConstIter (labelHashSet, curSet, iter)
+                {
+                    // Get the cell index and the cell
+                    const label& cellI = iter.key();
+                    const cell& cFaces = meshCells[cellI];
+
+                    // Note: there's no need to check for the distance here
+                    // because there's always at least one "buffer" layer
+                    // towards the outer side that consists of donors, which are
+                    // marked as ineligible at the beginning
+
+                    // Loop through all faces of the cell and mark all of them
+                    // for propagation
+                    forAll (cFaces, i)
+                    {
+                        propagateFace[cFaces[i]] = true;
+                    }
+                }
+
+                // Sync the face list across processor boundaries
+                syncTools::syncFaceList
+                (
+                    mesh,
+                    propagateFace,
+                    orOp<bool>(),
+                    false
+                );
+
+                // Loop through all faces and append interior holes
+                for (label faceI = 0; faceI < mesh.nInternalFaces(); ++faceI)
+                {
+                    if (propagateFace[faceI])
+                    {
+                        // Face is marked, select owner or neighbour
+                        const label& own = owner[faceI];
+                        const label& nei = neighbour[faceI];
+
+                        // Either owner or neighbour may be eligible, not both
+                        if (eligibleCells[own])
+                        {
+                            // Owner cell is not a hole (or an acceptor in the
+                            // first iteration), insert it
+                            if (fringeHoles.insert(own))
+                            {
+                                // Count number of added holes in this iteration
+                                ++nAddedHoles;
+                            }
+
+                            // Mark as ineligible
+                            eligibleCells[own] = false;
+                        }
+                        else if (eligibleCells[nei])
+                        {
+                            // Neighbour cell is not a hole (or an acceptor in
+                            // the first iteration), insert it
+                            if (fringeHoles.insert(nei))
+                            {
+                                // Count number of added holes in this iteration
+                                ++nAddedHoles;
+                            }
+
+                            // Mark as ineligible
+                            eligibleCells[nei] = false;
+                        }
+                    }
+                }
+
+                // Loop through boundary faces
+                for
+                (
+                    label faceI = mesh.nInternalFaces();
+                    faceI < mesh.nFaces();
+                  ++faceI
+                )
+                {
+                    if (propagateFace[faceI])
+                    {
+                        // Face is marked, select owner if this is the right
+                        // side. Neighbour handled on the other side
+                        const label& own = owner[faceI];
+
+                        if (eligibleCells[own])
+                        {
+                            // Face cell is not a hole (or an acceptor in the
+                            // first iteration), inser it
+                            if (fringeHoles.insert(own))
+                            {
+                                ++nAddedHoles;
+                            }
+
+                            // Mark as ineligible
+                            eligibleCells[own] = false;
+                        }
+                    }
+                }
+
+                // Need to reduce nAddedHoles in order to have synced loops in
+                // parallel
+                reduce(nAddedHoles, sumOp<label>());
+
+                // We moved one layer "inside" the fringe. Keep going until
+                // there are no more holes to add
+
+            } while (nAddedHoles != 0);
+
+            // Finally, we have collected all the fringe holes and acceptors
+            // from this connected region. Append them to the global sets
+            allAcceptors += acceptors;
+            allFringeHoles += fringeHoles;
+
+        } // End for all connected regions
+
+        // Set acceptors and holes from the data for all regions
+        acceptorsPtr_ = new labelList(allAcceptors.sortedToc());
+        fringeHolesPtr_ = new labelList(allFringeHoles.sortedToc());
+    }
+    else
+    {
+        // Connected fringes are not ready, allocate empty lists for acceptors
+        // and holes, which will be deleted when asked for again from the
+        // iterative procedure (see candidateAcceptors() and fringeHoles())
+        acceptorsPtr_ = new labelList;
+        fringeHolesPtr_ = new labelList;
+    }
+
+    if (debug)
+    {
+        Info<< "In donorBasedLayeredOverlapFringe::calcAddressing() const" << nl
+            << "Found " << acceptorsPtr_->size() << " acceptors." << nl
+            << "Found " << fringeHolesPtr_->size() << " fringe holes." << endl;
+    }
+}
+
+
+void Foam::donorBasedLayeredOverlapFringe::clearAddressing() const
+{
+    deleteDemandDrivenData(fringeHolesPtr_);
+    deleteDemandDrivenData(acceptorsPtr_);
+    deleteDemandDrivenData(finalDonorAcceptorsPtr_);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from dictionary
+Foam::donorBasedLayeredOverlapFringe::donorBasedLayeredOverlapFringe
+(
+    const fvMesh& mesh,
+    const oversetRegion& region,
+    const dictionary& dict
+)
+:
+    oversetFringe(mesh, region, dict),
+    connectedRegionNames_(dict.lookup("connectedRegions")),
+    connectedRegionIDs_(),
+    regionCentrePoints_
+    (
+        dict.lookupOrDefault<List<point> >
+        (
+            "regionCentrePoints",
+            List<point>(0)
+        )
+    ),
+    nLayers_(readLabel(dict.lookup("nLayers"))),
+    fringeHolesPtr_(nullptr),
+    acceptorsPtr_(nullptr),
+    finalDonorAcceptorsPtr_(nullptr),
+    isInitialized_(false)
+{
+    // Sanity check number of layers: must be greater than 0
+    if (nLayers_ < 1)
+    {
+        FatalIOErrorIn
+        (
+            "donorBasedLayeredOverlapFringe::"
+            "donorBasedLayeredOverlapFringe\n"
+            "(\n"
+            "    const fvMesh& mesh,\n"
+            "    const oversetRegion& region,\n"
+            "    const dictionary& dict\n"
+            ")",
+            dict
+        )   << "Invalid number of layers specified, nLayers = " << nLayers_
+            << nl
+            << "The number should be greater than 0."
+            << abort(FatalError);
+    }
+
+    // Recommendation: use at least two layers in order to avoid going defining
+    // acceptors on the wrong side and filling in the whole region with holes
+    if (nLayers_ == 1)
+    {
+        WarningIn
+        (
+            "donorBasedLayeredOverlapFringe::"
+            "donorBasedLayeredOverlapFringe\n"
+            "(\n"
+            "    const fvMesh& mesh,\n"
+            "    const oversetRegion& region,\n"
+            "    const dictionary& dict\n"
+            ")"
+        )   << "It is recommended to use at least 2 layers (nLayers = 2) in"
+            << " order to avoid defining acceptor cells on the wrong side"
+            << " of the region."
+            << nl
+            << "Be sure to check the fringe layer for this region."
+            << endl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::donorBasedLayeredOverlapFringe::~donorBasedLayeredOverlapFringe()
+{
+    clearAddressing();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::donorBasedLayeredOverlapFringe::updateIteration
+(
+    donorAcceptorList& donorAcceptorRegionData
+) const
+{
+    // If the donorAcceptor list has been allocated, something went wrong with
+    // the iteration procedure (not-updated flag): this function has been called
+    // more than once, which should not happen for
+    // donorBasedLayeredOverlapFringe
+    if (finalDonorAcceptorsPtr_)
+    {
+        FatalErrorIn
+        (
+            "donorBasedLayeredOverlapFringe::"
+            "updateIteration(donorAcceptorList&) const"
+        )   << "finalDonorAcceptorPtr_ already allocated. Something went "
+            << "wrong with the iteration procedure (flag was not updated)."
+            << nl
+            << "This should not happen for donorBasedLayeredOverlapFringe."
+            << abort(FatalError);
+    }
+
+    if
+    (
+        fringeHolesPtr_ && acceptorsPtr_
+     && returnReduce(!acceptorsPtr_->empty(), orOp<bool>())
+    )
+    {
+        // Note: we first check whether fringeHoles and acceptors pointers are
+        // allocated. If they are, there must be at least one processor with
+        // more than 0 acceptors in order for this fringe to be valid.
+
+        // Allocate the list by reusing the argument list
+        finalDonorAcceptorsPtr_ = new donorAcceptorList
+        (
+            donorAcceptorRegionData,
+            true
+        );
+
+        // Set the flag to true (for all processors due to reduction in the if
+        // statement)
+        updateSuitableOverlapFlag(true);
+    }
+    else
+    {
+        // Delete fringeHolesPtr and acceptorsPtr to trigger calculation of
+        // addressing, this time with other fringes up-to-date
+        deleteDemandDrivenData(fringeHolesPtr_);
+        deleteDemandDrivenData(acceptorsPtr_);
+    }
+
+    return foundSuitableOverlap();
+}
+
+
+const Foam::labelList& Foam::donorBasedLayeredOverlapFringe::fringeHoles() const
+{
+    // Note: updateIteration deletes fringeHolesPtr if the suitable overlap is
+    // not found, thus preparing it for the next iteration when the other
+    // fringes should be ready
+
+    if (!fringeHolesPtr_)
+    {
+        calcAddressing();
+    }
+
+    return *fringeHolesPtr_;
+}
+
+
+const Foam::labelList&
+Foam::donorBasedLayeredOverlapFringe::candidateAcceptors() const
+{
+    // Note: updateIteration deletes acceptorsPtr if the suitable overlap is
+    // not found, thus preparing it for the next iteration when the other
+    // fringes should be ready
+
+    if (!acceptorsPtr_)
+    {
+        calcAddressing();
+    }
+
+    return *acceptorsPtr_;
+}
+
+
+Foam::donorAcceptorList&
+Foam::donorBasedLayeredOverlapFringe::finalDonorAcceptors() const
+{
+    if (!finalDonorAcceptorsPtr_)
+    {
+        FatalErrorIn("donorBasedLayeredOverlapFringe::finalDonorAcceptors()")
+            << "finalDonorAcceptorPtr_ not allocated. Make sure you have"
+            << " called donorBasedLayeredOverlapFringe::updateIteration() before"
+            << " asking for final set of donor/acceptor pairs."
+            << abort(FatalError);
+    }
+
+    if (!foundSuitableOverlap())
+    {
+        FatalErrorIn("donorBasedLayeredOverlapFringe::finalDonorAcceptors()")
+            << "Attemted to access finalDonorAcceptors but suitable overlap "
+            << "has not been found. This is not allowed. "
+            << abort(FatalError);
+    }
+
+    return *finalDonorAcceptorsPtr_;
+}
+
+
+void Foam::donorBasedLayeredOverlapFringe::update() const
+{
+    Info<< "donorBasedLayeredOverlapFringe::update() const" << endl;
+
+    // Clear out
+    clearAddressing();
+
+    // Set flag to false and clear final donor/acceptors only
+    deleteDemandDrivenData(finalDonorAcceptorsPtr_);
+    updateSuitableOverlapFlag(false);
+}
+
+
+// ************************************************************************* //

src/overset/oversetMesh/oversetFringe/donorBasedLayeredOverlapFringe/donorBasedLayeredOverlapFringe.H

@@ -0,0 +1,190 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.1
+    \\  /    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
+    donorBasedLayeredOverlapFringe
+
+Description
+    Fringe algorithm based on donors from different regions. The algorithm waits
+    until the final donor/acceptor assembly has been performed for all regions
+    whose donor region is this region. Then, all donors are collected and the
+    acceptors are cells neighbouring the donors nLayers towards the interior.
+    Interior is defined either by user-specified points for each region or as a
+    centre of volume of the donor cells in particular region.
+
+    This fringe algorithm is intended to be used along with the faceCells fringe
+    on the other side, where the cell sizes are not significantly different from
+    each other (e.g. 10:1 cell ratio, where the finer cells are found on the
+    background mesh will probably be problematic to correctly set-up).
+
+    Note: based on distance tolerance (see distTol_ member), specified number of
+    layers and the bounding box of the connected donor region, it is possible
+    that acceptors end also on the wrong side (from the region centre as opposed
+    to towards the region centre). I do all I can to prevent this by looking at
+    the bounding box and hinting that at least 2 layers should be used. Be sure
+    to check the overset assembly for this region using calcOverset utility.
+
+Author
+    Vuko Vukcevic, Wikki Ltd.  All rights reserved.
+
+SourceFiles
+    donorBasedLayeredOverlapFringe.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef donorBasedLayeredOverlapFringe_H
+#define donorBasedLayeredOverlapFringe_H
+
+#include "oversetFringe.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                Class donorBasedLayeredOverlapFringe Declaration
+\*---------------------------------------------------------------------------*/
+
+class donorBasedLayeredOverlapFringe
+:
+    public oversetFringe
+{
+    // Private data
+
+        //- Names of connected regions. Looked up on construction
+        wordList connectedRegionNames_;
+
+        //- Regions IDs from which the donors will be collected as a starting
+        //  point. Note: initialized in init private member function because we
+        //  cannot initialize it in constructor. This is because certain overset
+        //  regions (and their fringes) may not be initialized at this point.
+        mutable labelList connectedRegionIDs_;
+
+        //- Optional list of points representing a rough estimate of the centre
+        //  for each underlying connected region. If these are not provided, the
+        //  centre is calculated as the centre of all donors for a given
+        //  connected region
+        List<point> regionCentrePoints_;
+
+        //- How many layers to move away from connected region donors to define
+        //  acceptor (and holes)
+        label nLayers_;
+
+        //- Fringe hole cells
+        mutable labelList* fringeHolesPtr_;
+
+        //- Acceptor cells
+        mutable labelList* acceptorsPtr_;
+
+        //- Final donor/acceptor pairs for this region (fringe)
+        mutable donorAcceptorList* finalDonorAcceptorsPtr_;
+
+        //- Initialization helper
+        mutable bool isInitialized_;
+
+
+    // Private static data
+
+        //- Distance tolerance to determine propagation direction. Note:
+        //  absolute value, default = 0. Might be useful is some strange cases
+        static const debug::tolerancesSwitch distTol_;
+
+
+    // Private Member Functions
+
+        //- Initialization
+        void init() const;
+
+        //- Calculate hole and acceptor addressing
+        void calcAddressing() const;
+
+        //- Clear addressing
+        void clearAddressing() const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("donorBasedLayeredOverlap");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        donorBasedLayeredOverlapFringe
+        (
+            const fvMesh& mesh,
+            const oversetRegion& region,
+            const dictionary& dict
+        );
+
+        //- Disallow default bitwise copy construct
+        donorBasedLayeredOverlapFringe
+        (
+            const donorBasedLayeredOverlapFringe&
+        ) = delete;
+
+        //- Disallow default bitwise assignment
+        void operator=(const donorBasedLayeredOverlapFringe&) = delete;
+
+
+    // Destructor
+
+        virtual ~donorBasedLayeredOverlapFringe();
+
+
+    // Member Functions
+
+        //- Update iteration. Note: invalidates parameter
+        virtual bool updateIteration
+        (
+            donorAcceptorList& donorAcceptorRegionData
+        ) const;
+
+        //- Return list of deactivated (hole) cells
+        //  Fringe hole cells are collected in addition to geometric hole
+        //  cells, which fall outside of all donor regions
+        virtual const labelList& fringeHoles() const;
+
+        //- Return list of acceptor cells
+        virtual const labelList& candidateAcceptors() const;
+
+        //- Return list of final donor acceptor pairs. Note: caller may
+        //  invalidate finalDonorAcceptorsPtr_ for optimisation purposes
+        virtual donorAcceptorList& finalDonorAcceptors() const;
+
+        //- Update the fringe
+        virtual void update() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/overset/oversetMesh/oversetFringe/faceCellsFringe/faceCellsFringe.C

@@ -72,7 +72,7 @@ void Foam::faceCellsFringe::calcAddressing() const
             (
                 "void faceCellsFringe::calcAddressing() const"
             )   << "Fringe patch " << patchNames_[nameI]
-                << " cannot be found"
+                << " cannot be found."
                 << abort(FatalError);
         }
 

src/overset/oversetMesh/oversetFringe/faceCellsFringe/faceCellsFringe.H

@@ -115,6 +115,12 @@ public:
 
     // Member Functions
 
+        //- Access to patch names
+        const wordList& patchNames() const
+        {
+            return patchNames_;
+        }
+
         //- Update iteration. Note: invalidates parameter
         virtual bool updateIteration
         (

src/overset/oversetMesh/oversetMesh/oversetMesh.H

@@ -121,8 +121,9 @@ private:
             //- Return overset type indicator field
             mutable volScalarField* oversetTypesPtr_;
 
-            //- Region ID: region index for each cell
-            mutable labelList* regionIDPtr_;
+            //- Region ID: region index for each cell as a volScalarField for
+            //  visualization. VV, 15/Apr/2019
+            mutable volScalarField* regionIDPtr_;
 
 
         // Overset discretisation support
@@ -286,7 +287,7 @@ public:
             const volScalarField& oversetTypes() const;
 
             //- Return region indicator
-            const labelList& regionID() const;
+            const volScalarField& regionID() const;
 
 
         // Overset discretisation support

src/overset/oversetMesh/oversetMesh/oversetMeshAddressing.C

@@ -27,6 +27,7 @@ License
 #include "surfaceFields.H"
 #include "volFields.H"
 #include "polyPatchID.H"
+#include "oversetFvPatchFields.H"
 #include "demandDrivenData.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -190,24 +191,52 @@ void Foam::oversetMesh::calcDomainMarkup() const
         }
     }
 
-    // Region ID
-    regionIDPtr_ = new labelList(mesh().nCells(), -1);
-    labelList& rID = *regionIDPtr_;
-
-    // Mark regions
+    // Region ID, initialized with -1 for sanity check later on
+    regionIDPtr_ = new volScalarField
+    (
+        IOobject
+        (
+            "regionIndex",
+            mesh().time().timeName(),
+            mesh(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh(),
+        dimensionedScalar("minusOne", dimless, -1.0),
+        "zeroGradient"
+    );
+    volScalarField& regionID = *regionIDPtr_;
+    scalarField& regionIDIn = regionID.internalField();
 
+    // Mark regions (internal field)
     forAll (regions_, regionI)
     {
         const labelList& curCells = regions_[regionI].zone();
 
         forAll (curCells, curCellI)
         {
-            rID[curCells[curCellI]] = regionI;
+            regionIDIn[curCells[curCellI]] = regionI;
+        }
+    }
+
+    // Update boundary values, making sure that we skip the overset patch
+    volScalarField::GeometricBoundaryField& regionIDb =
+        regionID.boundaryField();
+
+    forAll (regionIDb, patchI)
+    {
+        // Get the patch field
+        fvPatchScalarField& ripf = regionIDb[patchI];
+
+        if (!isA<oversetFvPatchScalarField>(ripf))
+        {
+            ripf = ripf.patchInternalField();
         }
     }
 
     // Check regions
-    if (min(rID) < 0)
+    if (min(regionID).value() < 0)
     {
         FatalErrorIn("void oversetMesh::calcDomainMarkup() const")
             << "Found cells without region ID.  Please check overset setup"
@@ -1166,7 +1195,7 @@ const Foam::volScalarField& Foam::oversetMesh::oversetTypes() const
 }
 
 
-const Foam::labelList& Foam::oversetMesh::regionID() const
+const Foam::volScalarField& Foam::oversetMesh::regionID() const
 {
     if (!regionIDPtr_)
     {

src/overset/oversetMesh/oversetRegion/oversetRegion.C

@@ -1730,6 +1730,20 @@ Foam::oversetRegion::~oversetRegion()
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
+const Foam::oversetFringe& Foam::oversetRegion::fringe() const
+{
+    if (fringePtr_.empty())
+    {
+        FatalErrorIn("const oversetFringe& oversetRegion::fringe() const")
+            << "Fringe pointer not allocated. It should have been initialized"
+            << " properly at construction. Something went wrong..."
+            << abort(FatalError);
+    }
+
+    return fringePtr_();
+}
+
+
 const Foam::labelList& Foam::oversetRegion::donorRegions() const
 {
     if (!donorRegionsPtr_)

src/overset/oversetMesh/oversetRegion/oversetRegion.H

@@ -327,6 +327,9 @@ public:
                 return zone();
             }
 
+            //- Return overset fringe algorithm for this region
+            const oversetFringe& fringe() const;
+
             //- Return list of donor region indices
             const labelList& donorRegions() const;
 

src/thermophysicalModels/basic/derivedFvPatchFields/gradientEnthalpy/gradientEnthalpyFvPatchScalarField.C

@@ -192,7 +192,6 @@ void Foam::gradientEnthalpyFvPatchScalarField::updateCoeffs()
 void Foam::gradientEnthalpyFvPatchScalarField::write(Ostream& os) const
 {
     fixedGradientFvPatchScalarField::write(os);
-    writeEntry("value", os);
 }
 
 

src/thermophysicalModels/basic/derivedFvPatchFields/gradientInternalEnergy/gradientInternalEnergyFvPatchScalarField.C

@@ -128,7 +128,6 @@ void gradientInternalEnergyFvPatchScalarField::updateCoeffs()
 void Foam::gradientInternalEnergyFvPatchScalarField::write(Ostream& os) const
 {
     fixedGradientFvPatchScalarField::write(os);
-    writeEntry("value", os);
 }
 
 

src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatWallFunction/alphatWallFunctionFvPatchScalarField.C

@@ -37,6 +37,15 @@ namespace compressible
 namespace RASModels
 {
 
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+void alphatWallFunctionFvPatchScalarField::writeLocalEntries(Ostream& os) const
+{
+    writeEntryIfDifferent<word>(os, "mut", "mut", mutName_);
+    os.writeKeyword("Prt") << Prt_ << token::END_STATEMENT << nl;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 alphatWallFunctionFvPatchScalarField::
@@ -127,8 +136,7 @@ void alphatWallFunctionFvPatchScalarField::updateCoeffs()
 void alphatWallFunctionFvPatchScalarField::write(Ostream& os) const
 {
     fvPatchField<scalar>::write(os);
-    writeEntryIfDifferent<word>(os, "mut", "mut", mutName_);
-    os.writeKeyword("Prt") << Prt_ << token::END_STATEMENT << nl;
+    writeLocalEntries(os);
     writeEntry("value", os);
 }
 

src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatWallFunction/alphatWallFunctionFvPatchScalarField.H

@@ -66,6 +66,14 @@ class alphatWallFunctionFvPatchScalarField
         scalar Prt_;
 
 
+protected:
+
+    // Protected Member Functions
+
+        //- Write local wall function variables
+        virtual void writeLocalEntries(Ostream&) const;
+
+
 public:
 
     //- Runtime type information
@@ -146,7 +154,7 @@ public:
         // I-O
 
             //- Write
-            void write(Ostream&) const;
+            virtual void write(Ostream&) const;
 };
 
 

src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C

@@ -54,6 +54,21 @@ void epsilonWallFunctionFvPatchScalarField::checkType()
 }
 
 
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+void epsilonWallFunctionFvPatchScalarField::writeLocalEntries(Ostream& os) const
+{
+    writeEntryIfDifferent<word>(os, "U", "U", UName_);
+    writeEntryIfDifferent<word>(os, "k", "k", kName_);
+    writeEntryIfDifferent<word>(os, "G", "RASModel::G", GName_);
+    writeEntryIfDifferent<word>(os, "mu", "mu", muName_);
+    writeEntryIfDifferent<word>(os, "mut", "mut", mutName_);
+    os.writeKeyword("Cmu") << Cmu_ << token::END_STATEMENT << nl;
+    os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
+    os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 epsilonWallFunctionFvPatchScalarField::epsilonWallFunctionFvPatchScalarField
@@ -251,7 +266,7 @@ void epsilonWallFunctionFvPatchScalarField::updateCoeffs()
 
     // TODO: perform averaging for cells sharing more than one boundary face
 
-    fixedInternalValueFvPatchField<scalar>::updateCoeffs();
+    fixedInternalValueFvPatchScalarField::updateCoeffs();
 }
 
 
@@ -260,22 +275,14 @@ void epsilonWallFunctionFvPatchScalarField::evaluate
     const Pstream::commsTypes commsType
 )
 {
-    fixedInternalValueFvPatchField<scalar>::evaluate(commsType);
+    fixedInternalValueFvPatchScalarField::evaluate(commsType);
 }
 
 
 void epsilonWallFunctionFvPatchScalarField::write(Ostream& os) const
 {
-    fixedInternalValueFvPatchField<scalar>::write(os);
-    writeEntryIfDifferent<word>(os, "U", "U", UName_);
-    writeEntryIfDifferent<word>(os, "k", "k", kName_);
-    writeEntryIfDifferent<word>(os, "G", "RASModel::G", GName_);
-    writeEntryIfDifferent<word>(os, "rho", "rho", rhoName_);
-    writeEntryIfDifferent<word>(os, "mu", "mu", muName_);
-    writeEntryIfDifferent<word>(os, "mut", "mut", mutName_);
-    os.writeKeyword("Cmu") << Cmu_ << token::END_STATEMENT << nl;
-    os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
-    os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
+    fixedInternalValueFvPatchScalarField::write(os);
+    writeLocalEntries(os);
 }
 
 

src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.H

@@ -37,7 +37,7 @@ SourceFiles
 #ifndef compressibleEpsilonWallFunctionFvPatchScalarField_H
 #define compressibleEpsilonWallFunctionFvPatchScalarField_H
 
-#include "fixedInternalValueFvPatchField.H"
+#include "fixedInternalValueFvPatchFields.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -92,6 +92,20 @@ class epsilonWallFunctionFvPatchScalarField
         void checkType();
 
 
+protected:
+
+    // Protected Member Functions
+
+        //- Write local wall function variables
+        virtual void writeLocalEntries(Ostream&) const;
+
+        //- Return name of turbulence generation field
+        inline const word& GName() const
+        {
+            return GName_;
+        }
+
+
 public:
 
     //- Runtime type information
@@ -169,7 +183,10 @@ public:
             virtual void updateCoeffs();
 
             //- Evaluate the patchField
-            virtual void evaluate(const Pstream::commsTypes);
+            virtual void evaluate
+            (
+                const Pstream::commsTypes  = Pstream::blocking
+            );
 
 
         // I-O

src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/kqRWallFunctions/kqRWallFunction/kqRWallFunctionFvPatchField.C

@@ -47,7 +47,8 @@ void kqRWallFunctionFvPatchField<Type>::checkType()
             << "Invalid wall function specification" << nl
             << "    Patch type for patch " << this->patch().name()
             << " must be wall" << nl
-            << "    Current patch type is " << this->patch().type() << nl << endl
+            << "    Current patch type is " << this->patch().type()
+            << nl << endl
             << abort(FatalError);
     }
 }
@@ -124,16 +125,6 @@ kqRWallFunctionFvPatchField<Type>::kqRWallFunctionFvPatchField
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class Type>
-void kqRWallFunctionFvPatchField<Type>::evaluate
-(
-    const Pstream::commsTypes commsType
-)
-{
-    zeroGradientFvPatchField<Type>::evaluate(commsType);
-}
-
-
 template<class Type>
 void kqRWallFunctionFvPatchField<Type>::write(Ostream& os) const
 {

src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/kqRWallFunctions/kqRWallFunction/kqRWallFunctionFvPatchField.H

@@ -134,15 +134,6 @@ public:
 
     // Member functions
 
-        // Evaluation functions
-
-            //- Evaluate the patchField
-            virtual void evaluate
-            (
-                const Pstream::commsTypes commsType=Pstream::Pstream::blocking
-            );
-
-
         // I-O
 
             //- Write

src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C

@@ -54,6 +54,22 @@ void omegaWallFunctionFvPatchScalarField::checkType()
 }
 
 
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+void omegaWallFunctionFvPatchScalarField::writeLocalEntries(Ostream& os) const
+{
+    writeEntryIfDifferent<word>(os, "U", "U", UName_);
+    writeEntryIfDifferent<word>(os, "k", "k", kName_);
+    writeEntryIfDifferent<word>(os, "G", "RASModel::G", GName_);
+    writeEntryIfDifferent<word>(os, "mu", "mu", muName_);
+    writeEntryIfDifferent<word>(os, "mut", "mut", mutName_);
+    os.writeKeyword("Cmu") << Cmu_ << token::END_STATEMENT << nl;
+    os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
+    os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
+    os.writeKeyword("beta1") << beta1_ << token::END_STATEMENT << nl;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 omegaWallFunctionFvPatchScalarField::omegaWallFunctionFvPatchScalarField
@@ -264,16 +280,7 @@ void omegaWallFunctionFvPatchScalarField::updateCoeffs()
 void omegaWallFunctionFvPatchScalarField::write(Ostream& os) const
 {
     fixedInternalValueFvPatchField<scalar>::write(os);
-    writeEntryIfDifferent<word>(os, "U", "U", UName_);
-    writeEntryIfDifferent<word>(os, "rho", "rho", rhoName_);
-    writeEntryIfDifferent<word>(os, "k", "k", kName_);
-    writeEntryIfDifferent<word>(os, "G", "RASModel::G", GName_);
-    writeEntryIfDifferent<word>(os, "mu", "mu", muName_);
-    writeEntryIfDifferent<word>(os, "mut", "mut", mutName_);
-    os.writeKeyword("Cmu") << Cmu_ << token::END_STATEMENT << nl;
-    os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
-    os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
-    os.writeKeyword("beta1") << beta1_ << token::END_STATEMENT << nl;
+    writeLocalEntries(os);
 }