diff --git a/applications/solvers/basic/laplacianFoam/overLaplacianDyMFoam/createFields.H b/applications/solvers/basic/laplacianFoam/overLaplacianDyMFoam/createFields.H
index fc66a21b162a5710cc040d6118d2d3c78ba7010b..62036fcc76a46baa39344d2212790fe8b65076c0 100644
--- a/applications/solvers/basic/laplacianFoam/overLaplacianDyMFoam/createFields.H
+++ b/applications/solvers/basic/laplacianFoam/overLaplacianDyMFoam/createFields.H
@@ -13,23 +13,6 @@
         mesh
     );
 
-    // Add overset specific interpolations
-    {
-        dictionary oversetDict;
-        oversetDict.add("T", true);
-
-        const_cast<dictionary&>
-        (
-            mesh.schemesDict()
-        ).add
-        (
-            "oversetInterpolationRequired",
-            oversetDict,
-            true
-        );
-    }
-
-
     Info<< "Reading transportProperties\n" << endl;
 
     IOdictionary transportProperties
diff --git a/applications/solvers/basic/potentialFoam/overPotentialFoam/createFields.H b/applications/solvers/basic/potentialFoam/overPotentialFoam/createFields.H
index 80f04a8a68d47db1f17f151cce525365ab63b226..d39f3044b710ee9cf7c07aba1a8e0a99c0da20ff 100644
--- a/applications/solvers/basic/potentialFoam/overPotentialFoam/createFields.H
+++ b/applications/solvers/basic/potentialFoam/overPotentialFoam/createFields.H
@@ -121,21 +121,14 @@ mesh.setFluxRequired(Phi.name());
 
 #include "createMRF.H"
 
-// Add overset specific interpolations
+// Add solver-specific interpolations
 {
-    dictionary oversetDict;
-    oversetDict.add("Phi", true);
-    oversetDict.add("U", true);
+    wordHashSet& nonInt =
+        const_cast<wordHashSet&>(Stencil::New(mesh).nonInterpolatedFields());
+
+    nonInt.insert("cellMask");
+    nonInt.insert("interpolatedCells");
 
-    const_cast<dictionary&>
-    (
-        mesh.schemesDict()
-    ).add
-    (
-        "oversetInterpolationRequired",
-        oversetDict,
-        true
-    );
 }
 
 // Mask field for zeroing out contributions on hole cells
diff --git a/applications/test/complex/Test-complex.C b/applications/test/complex/Test-complex.C
index 0cf31f9c857b6627ff4b026f5638f0006059b5c9..bc4d3d05dc736941d6822e3f5c3a5565e309fd8f 100644
--- a/applications/test/complex/Test-complex.C
+++ b/applications/test/complex/Test-complex.C
@@ -25,12 +25,11 @@ Application
 
 Description
     Some tests for complex numbers
+
 \*---------------------------------------------------------------------------*/
 
 #include "argList.H"
-#include "complex.H"
-#include "complexVector.H"
-#include "Field.H"
+#include "complexFields.H"
 
 using namespace Foam;
 
@@ -52,15 +51,26 @@ int main(int argc, char *argv[])
         << "complex(scalar) : " << complex(3.14519) << nl
         << nl;
 
+    std::complex<scalar> c1(10, -3);
+    Info<< "std::complex : " << c1 << nl;
+    Info<< "sin: " << std::sin(c1) << nl;
+
+
     Info<< "complexVector::zero : " << complexVector::zero << nl
         << "complexVector::one  : " << complexVector::one << nl
         << nl;
 
-    // Comparison
-
-    for (complex c :  { complex{1, 0}, complex{1, 2}} )
+    for (complex c : { complex{1, 0}, complex{1, 2}} )
     {
+        Info<< nl;
         print1(c);
+
+        Info<< "sin: " << sin(c) << nl;
+        Info<< "pow(3): " << pow(c, 3) << nl;
+        Info<< "pow3: " << pow3(c) << nl;
+        Info<< "log: " << log(c) << nl;
+        Info<< "pow025: " << pow025(c) << nl;
+
         // TDB: allow implicit construct from scalar?
         //
         // if (c == 1.0)
@@ -69,22 +79,46 @@ int main(int argc, char *argv[])
         // }
     }
 
-    Field<complex> fld1(3, complex(2.0, 1.0));
+    complexField fld1(3, complex(2.0, 1.0));
+    complexField fld2(fld1);
+
+    for (complex& c : fld2)
+    {
+        c = ~c;
+    }
 
     Info<< "Field " << flatOutput(fld1) << nl;
+    Info<< "Conjugate: " << flatOutput(fld2) << nl;
+
+    // Some arbitrary change
+    for (complex& c : fld2)
+    {
+        c.Im() *= 5;
+    }
+
 
     Info<< "sum = " << sum(fld1) << nl;
     // Not yet Info<< "min = " << min(fld1) << nl;
 
     fld1 *= 10;
-    Info<< "Multiply: " << flatOutput(fld1) << nl;
+    Info<< "scalar multiply: " << flatOutput(fld1) << nl;
 
-    for (complex& c : fld1)
-    {
-        c = ~c;
-    }
+    fld1 /= 10;
+    Info<< "scalar divide: " << flatOutput(fld1) << nl;
+
+    Info<< "sin: " << sin(fld1) << nl;
+
+    Info<< "operator + : " << (fld1 + fld2) << nl;
+
+    // Some operators are still incomplete
+
+    // Info<< "operator * : " << (fld1 * fld2) << nl;
+    // Info<< "operator / : " << (fld1 / fld2) << nl;
+    Info<< "operator / : " << (fld1 / 2) << nl;
+    // Info<< "operator / : " << (fld1 / fld2) << nl;
+    Info<< "sqrt   : " << sqrt(fld1) << nl;
+    // Info<< "pow(2) : " << pow(fld1, 2) << nl;
 
-    Info<< "Conjugate: " << flatOutput(fld1) << nl;
 
     Info<< "\nEnd\n" << endl;
     return 0;
diff --git a/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C b/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C
index 3fa384024651a9534951f3f9270d37b2b41d7e83..aafa48c28d8157e9c48c2178ecf1add77e77f61e 100644
--- a/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C
+++ b/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C
@@ -217,6 +217,7 @@ int main(int argc, char *argv[])
     );
 
     #include "addRegionOption.H"
+    argList::addOption("dict", "file", "Use alternative extrudeMeshDict");
     #include "setRootCase.H"
     #include "createTimeExtruded.H"
 
@@ -236,15 +237,19 @@ int main(int argc, char *argv[])
             << runTimeExtruded.timeName() << nl << endl;
     }
 
-
-    IOdictionary dict
+    const IOdictionary dict
     (
-        IOobject
+        IOobject::selectIO
         (
-            "extrudeMeshDict",
-            runTimeExtruded.system(),
-            runTimeExtruded,
-            IOobject::MUST_READ_IF_MODIFIED
+            IOobject
+            (
+                "extrudeMeshDict",
+                runTimeExtruded.system(),
+                runTimeExtruded,
+                IOobject::MUST_READ_IF_MODIFIED,
+                IOobject::NO_WRITE
+            ),
+            args.opt<fileName>("dict", "")
         )
     );
 
diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
index 7b0eba46d66f08086f755bfead4f93710c75770b..d3040ffac3b7ad04499d8477745ff565fe150cb5 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2015-2018 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2015-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2016 OpenFOAM Foundation
@@ -1691,19 +1691,41 @@ int main(int argc, char *argv[])
     }
 
 
-    const bool mergePatchFaces
-    (
-        meshDict.lookupOrDefault("mergePatchFaces", true)
-    );
-
-    if (!mergePatchFaces)
+    // How to treat co-planar faces
+    meshRefinement::FaceMergeType mergeType =
+        meshRefinement::FaceMergeType::GEOMETRIC;
     {
-        Info<< "Not merging patch-faces of cell to preserve"
-            << " (split)hex cell shape."
-            << nl << endl;
+        const bool mergePatchFaces
+        (
+            meshDict.lookupOrDefault("mergePatchFaces", true)
+        );
+
+        if (!mergePatchFaces)
+        {
+            Info<< "Not merging patch-faces of cell to preserve"
+                << " (split)hex cell shape."
+                << nl << endl;
+            mergeType = meshRefinement::FaceMergeType::NONE;
+        }
+        else
+        {
+            const bool mergeAcrossPatches
+            (
+                meshDict.lookupOrDefault("mergeAcrossPatches", false)
+            );
+
+            if (mergeAcrossPatches)
+            {
+                Info<< "Merging co-planar patch-faces of cells"
+                    << ", regardless of patch assignment"
+                    << nl << endl;
+                mergeType = meshRefinement::FaceMergeType::IGNOREPATCH;
+            }
+        }
     }
 
 
+
     if (wantRefine)
     {
         cpuTime timer;
@@ -1732,7 +1754,7 @@ int main(int argc, char *argv[])
             refineParams,
             snapParams,
             refineParams.handleSnapProblems(),
-            mergePatchFaces,        // merge co-planar faces
+            mergeType,
             motionDict
         );
 
@@ -1784,7 +1806,7 @@ int main(int argc, char *argv[])
         (
             snapDict,
             motionDict,
-            mergePatchFaces,
+            mergeType,
             curvature,
             planarAngle,
             snapParams
@@ -1851,7 +1873,7 @@ int main(int argc, char *argv[])
             layerDict,
             motionDict,
             layerParams,
-            mergePatchFaces,
+            mergeType,
             preBalance,
             decomposer,
             distributor
diff --git a/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C b/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
index e3bfa82f355c976b9d647c9b749503355c4bad17..9f66e9c61989d54977275ef55a114d18770ec672 100644
--- a/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
+++ b/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
@@ -55,7 +55,7 @@ Description
 #include "faceSet.H"
 #include "pointSet.H"
 #include "processorMeshes.H"
-#include "hexRef8.H"
+#include "hexRef8Data.H"
 
 #ifdef HAVE_ZOLTAN
     #include "zoltanRenumber.H"
@@ -1327,8 +1327,24 @@ int main(int argc, char *argv[])
 
     // Remove old procAddressing files
     processorMeshes::removeFiles(mesh);
-    // Remove refinement data
-    hexRef8::removeFiles(mesh);
+
+    // Update refinement data
+    hexRef8Data refData
+    (
+        IOobject
+        (
+            "dummy",
+            mesh.facesInstance(),
+            polyMesh::meshSubDir,
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE,
+            false
+        )
+    );
+    refData.updateMesh(map());
+    refData.write();
+
     // Update sets
     topoSet::updateMesh(mesh.facesInstance(), map(), cellSets);
     topoSet::updateMesh(mesh.facesInstance(), map(), faceSets);
diff --git a/src/OpenFOAM/Make/files b/src/OpenFOAM/Make/files
index 54c3b2a19a9691c8614659b8416d7caa8784e2e4..b3c9ecfdcc38b88ff3b8b2f3a7fd4a92f5c510be 100644
--- a/src/OpenFOAM/Make/files
+++ b/src/OpenFOAM/Make/files
@@ -672,7 +672,8 @@ $(Fields)/quaternionField/quaternionField.C
 $(Fields)/quaternionField/quaternionIOField.C
 $(Fields)/triadField/triadField.C
 $(Fields)/triadField/triadIOField.C
-$(Fields)/complexFields/complexFields.C
+$(Fields)/complex/complexField.C
+$(Fields)/complex/complexVectorField.C
 $(Fields)/transformField/transformField.C
 $(Fields)/fieldTypes.C
 
diff --git a/src/OpenFOAM/containers/Bits/bitSet/bitSet.H b/src/OpenFOAM/containers/Bits/bitSet/bitSet.H
index 2c5c283136ce51a4bcbb12439a46bdd4d423a472..00d943059686761e69f56e7efc65f14de9bbdea6 100644
--- a/src/OpenFOAM/containers/Bits/bitSet/bitSet.H
+++ b/src/OpenFOAM/containers/Bits/bitSet/bitSet.H
@@ -513,10 +513,6 @@ public:
         //- Move assignment
         inline bitSet& operator=(bitSet&& bitset);
 
-        //- Complement operator.
-        //  Return a copy of the existing set with all its bits flipped.
-        inline bitSet operator~() const;
-
         //- Bitwise-AND all the bits in other with the bits in this bitset.
         //  The operands may have dissimilar sizes without affecting the size
         //  of the set.
@@ -555,14 +551,6 @@ public:
         {
             return *this;
         }
-
-
-    // Housekeeping
-
-        //- Deprecated(2018-04) compatibility name for PackedBoolList
-        //  \deprecated(2018-04) - use toc() method
-        inline labelList used() const { return toc(); }
-
 };
 
 
@@ -576,6 +564,9 @@ Ostream& operator<<(Ostream& os, const bitSet& bitset);
 Ostream& operator<<(Ostream& os, const InfoProxy<bitSet>& info);
 
 
+//- Bitset complement, returns a copy of the bitset with all its bits flipped
+inline bitSet operator~(const bitSet& bitset);
+
 //- Bitwise-AND of two bitsets.
 //  See bitSet::operator&= for more details.
 inline bitSet operator&(const bitSet& a, const bitSet& b);
diff --git a/src/OpenFOAM/containers/Bits/bitSet/bitSetI.H b/src/OpenFOAM/containers/Bits/bitSet/bitSetI.H
index afc7d93f9719743939df5ad61387773b0172e991..5312caa3d1b1c086890ab94a5423e34e9158e568 100644
--- a/src/OpenFOAM/containers/Bits/bitSet/bitSetI.H
+++ b/src/OpenFOAM/containers/Bits/bitSet/bitSetI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2018-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -663,14 +663,6 @@ inline Foam::bitSet& Foam::bitSet::operator=(bitSet&& bitset)
 }
 
 
-inline Foam::bitSet Foam::bitSet::operator~() const
-{
-    bitSet result(*this);
-    result.flip();
-    return result;
-}
-
-
 inline Foam::bitSet& Foam::bitSet::operator&=(const bitSet& other)
 {
     return andEq(other);
@@ -697,6 +689,14 @@ inline Foam::bitSet& Foam::bitSet::operator-=(const bitSet& other)
 
 // * * * * * * * * * * * * * * * Global Operators  * * * * * * * * * * * * * //
 
+inline Foam::bitSet Foam::operator~(const bitSet& bitset)
+{
+    bitSet result(bitset);
+    result.flip();
+    return result;
+}
+
+
 inline Foam::bitSet Foam::operator&(const bitSet& a, const bitSet& b)
 {
     bitSet result(a);
diff --git a/src/OpenFOAM/fields/Fields/complex/complexField.C b/src/OpenFOAM/fields/Fields/complex/complexField.C
new file mode 100644
index 0000000000000000000000000000000000000000..5662ab40763d26425db2041d6bf1eb9b5f9c857f
--- /dev/null
+++ b/src/OpenFOAM/fields/Fields/complex/complexField.C
@@ -0,0 +1,168 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2004-2010, 2019 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+                            | Copyright (C) 2011 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "complexField.H"
+#include "addToRunTimeSelectionTable.H"
+
+#define TEMPLATE
+#include "FieldFunctionsM.C"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineCompoundTypeName(List<complex>, complexList);
+    addCompoundToRunTimeSelectionTable(List<complex>, complexList);
+}
+
+
+// * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
+
+Foam::complexField Foam::ComplexField
+(
+    const UList<scalar>& re,
+    const UList<scalar>& im
+)
+{
+    complexField cf(re.size());
+
+    forAll(cf, i)
+    {
+        cf[i].Re() = re[i];
+        cf[i].Im() = im[i];
+    }
+
+    return cf;
+}
+
+
+Foam::complexField Foam::ReComplexField(const UList<scalar>& re)
+{
+    complexField cf(re.size());
+
+    forAll(cf, i)
+    {
+        cf[i].Re() = re[i];
+        cf[i].Im() = 0.0;
+    }
+
+    return cf;
+}
+
+
+Foam::complexField Foam::ImComplexField(const UList<scalar>& im)
+{
+    complexField cf(im.size());
+
+    forAll(cf, i)
+    {
+        cf[i].Re() = 0.0;
+        cf[i].Im() = im[i];
+    }
+
+    return cf;
+}
+
+
+Foam::scalarField Foam::ReImSum(const UList<complex>& cf)
+{
+    scalarField sf(cf.size());
+
+    forAll(sf, i)
+    {
+        sf[i] = cf[i].Re() + cf[i].Im();
+    }
+
+    return sf;
+}
+
+
+Foam::scalarField Foam::Re(const UList<complex>& cf)
+{
+    scalarField sf(cf.size());
+
+    forAll(sf, i)
+    {
+        sf[i] = cf[i].Re();
+    }
+
+    return sf;
+}
+
+
+Foam::scalarField Foam::Im(const UList<complex>& cf)
+{
+    scalarField sf(cf.size());
+
+    forAll(sf, i)
+    {
+        sf[i] = cf[i].Im();
+    }
+
+    return sf;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+UNARY_FUNCTION(complex, complex, pow3)
+UNARY_FUNCTION(complex, complex, pow4)
+UNARY_FUNCTION(complex, complex, pow5)
+UNARY_FUNCTION(complex, complex, pow6)
+UNARY_FUNCTION(complex, complex, pow025)
+UNARY_FUNCTION(complex, complex, sqrt)
+UNARY_FUNCTION(complex, complex, exp)
+UNARY_FUNCTION(complex, complex, log)
+UNARY_FUNCTION(complex, complex, log10)
+UNARY_FUNCTION(complex, complex, sin)
+UNARY_FUNCTION(complex, complex, cos)
+UNARY_FUNCTION(complex, complex, tan)
+UNARY_FUNCTION(complex, complex, asin)
+UNARY_FUNCTION(complex, complex, acos)
+UNARY_FUNCTION(complex, complex, atan)
+UNARY_FUNCTION(complex, complex, sinh)
+UNARY_FUNCTION(complex, complex, cosh)
+UNARY_FUNCTION(complex, complex, tanh)
+UNARY_FUNCTION(complex, complex, asinh)
+UNARY_FUNCTION(complex, complex, acosh)
+UNARY_FUNCTION(complex, complex, atanh)
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "undefFieldFunctionsM.H"
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/fields/Fields/complex/complexField.H b/src/OpenFOAM/fields/Fields/complex/complexField.H
new file mode 100644
index 0000000000000000000000000000000000000000..c9b9351364f672d9a133159619638b9a4a408978
--- /dev/null
+++ b/src/OpenFOAM/fields/Fields/complex/complexField.H
@@ -0,0 +1,114 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2004-2010, 2019 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+                            | Copyright (C) 2011 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Typedef
+    Foam::complexField
+
+Description
+    Specialisation of Field\<T\> for complex.
+
+SourceFiles
+    complexField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef complexField_H
+#define complexField_H
+
+#include "complex.H"
+#include "scalarField.H"
+
+#define TEMPLATE
+#include "FieldFunctionsM.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+typedef Field<complex> complexField;
+
+//- Zip up two lists of values into a list of complex
+complexField ComplexField
+(
+    const UList<scalar>& re,
+    const UList<scalar>& im
+);
+
+//- Create complex field from a list of real (using imag == 0)
+complexField ReComplexField(const UList<scalar>& re);
+
+//- Create complex field from a list of imag (using real == 0)
+complexField ImComplexField(const UList<scalar>& im);
+
+//- Extract real component
+scalarField Re(const UList<complex>& cf);
+
+//- Extract imag component
+scalarField Im(const UList<complex>& cf);
+
+//- Sum real and imag components
+scalarField ReImSum(const UList<complex>& cf);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+UNARY_FUNCTION(complex, complex, pow3)
+UNARY_FUNCTION(complex, complex, pow4)
+UNARY_FUNCTION(complex, complex, pow5)
+UNARY_FUNCTION(complex, complex, pow6)
+UNARY_FUNCTION(complex, complex, pow025)
+UNARY_FUNCTION(complex, complex, sqrt)
+UNARY_FUNCTION(complex, complex, exp)
+UNARY_FUNCTION(complex, complex, log)
+UNARY_FUNCTION(complex, complex, log10)
+UNARY_FUNCTION(complex, complex, sin)
+UNARY_FUNCTION(complex, complex, cos)
+UNARY_FUNCTION(complex, complex, tan)
+UNARY_FUNCTION(complex, complex, asin)
+UNARY_FUNCTION(complex, complex, acos)
+UNARY_FUNCTION(complex, complex, atan)
+UNARY_FUNCTION(complex, complex, sinh)
+UNARY_FUNCTION(complex, complex, cosh)
+UNARY_FUNCTION(complex, complex, tanh)
+UNARY_FUNCTION(complex, complex, asinh)
+UNARY_FUNCTION(complex, complex, acosh)
+UNARY_FUNCTION(complex, complex, atanh)
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "undefFieldFunctionsM.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/fields/Fields/complex/complexFields.H b/src/OpenFOAM/fields/Fields/complex/complexFields.H
new file mode 100644
index 0000000000000000000000000000000000000000..ef031db0e439415c81a91dcfc5676be0a7571ab6
--- /dev/null
+++ b/src/OpenFOAM/fields/Fields/complex/complexFields.H
@@ -0,0 +1,44 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+InClass
+    Foam::complexFields
+
+Description
+    Specialisations of Field\<T\> for complex and complexVector
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef complexFields_H
+#define complexFields_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "complexField.H"
+#include "complexVectorField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/fields/Fields/complex/complexVectorField.C b/src/OpenFOAM/fields/Fields/complex/complexVectorField.C
new file mode 100644
index 0000000000000000000000000000000000000000..3ac98022f5fc1c5b49cdf280f2bfd3f0f8145d2b
--- /dev/null
+++ b/src/OpenFOAM/fields/Fields/complex/complexVectorField.C
@@ -0,0 +1,156 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2004-2010, 2019 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+                            | Copyright (C) 2011 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "complexVectorField.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineCompoundTypeName(List<complexVector>, complexVectorList);
+    addCompoundToRunTimeSelectionTable(List<complexVector>, complexVectorList);
+}
+
+
+// * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
+
+Foam::complexVectorField Foam::ComplexField
+(
+    const UList<vector>& re,
+    const UList<vector>& im
+)
+{
+    complexVectorField cvf(re.size());
+
+    for (direction cmpt=0; cmpt<vector::nComponents; ++cmpt)
+    {
+        forAll(cvf, i)
+        {
+            cvf[i].component(cmpt).Re() = re[i].component(cmpt);
+            cvf[i].component(cmpt).Im() = im[i].component(cmpt);
+        }
+    }
+
+    return cvf;
+}
+
+
+Foam::complexVectorField Foam::ReComplexField(const UList<vector>& re)
+{
+    complexVectorField cvf(re.size());
+
+    for (direction cmpt=0; cmpt<vector::nComponents; ++cmpt)
+    {
+        forAll(cvf, i)
+        {
+            cvf[i].component(cmpt).Re() = re[i].component(cmpt);
+            cvf[i].component(cmpt).Im() = 0.0;
+        }
+    }
+
+    return cvf;
+}
+
+
+Foam::complexVectorField Foam::ImComplexField(const UList<vector>& im)
+{
+    complexVectorField cvf(im.size());
+
+    for (direction cmpt=0; cmpt<vector::nComponents; ++cmpt)
+    {
+        forAll(cvf, i)
+        {
+            cvf[i].component(cmpt).Re() = 0.0;
+            cvf[i].component(cmpt).Im() = im[i].component(cmpt);
+        }
+    }
+
+    return cvf;
+}
+
+
+Foam::vectorField Foam::ReImSum(const UList<complexVector>& cvf)
+{
+    vectorField vf(cvf.size());
+
+    for (direction cmpt=0; cmpt<vector::nComponents; ++cmpt)
+    {
+        forAll(cvf, i)
+        {
+            vf[i].component(cmpt) =
+                cvf[i].component(cmpt).Re() + cvf[i].component(cmpt).Im();
+        }
+    }
+
+    return vf;
+}
+
+
+Foam::vectorField Foam::Re(const UList<complexVector>& cvf)
+{
+    vectorField vf(cvf.size());
+
+    for (direction cmpt=0; cmpt<vector::nComponents; ++cmpt)
+    {
+        forAll(cvf, i)
+        {
+            vf[i].component(cmpt) = cvf[i].component(cmpt).Re();
+        }
+    }
+
+    return vf;
+}
+
+
+Foam::vectorField Foam::Im(const UList<complexVector>& cvf)
+{
+    vectorField vf(cvf.size());
+
+    for (direction cmpt=0; cmpt<vector::nComponents; ++cmpt)
+    {
+        forAll(cvf, i)
+        {
+            vf[i].component(cmpt) = cvf[i].component(cmpt).Im();
+        }
+    }
+
+    return vf;
+}
+
+
+Foam::complexVectorField Foam::operator^
+(
+    const UList<vector>& vf,
+    const UList<complexVector>& cvf
+)
+{
+    return ComplexField(vf^Re(cvf), vf^Im(cvf));
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/fields/Fields/complexFields/complexFields.H b/src/OpenFOAM/fields/Fields/complex/complexVectorField.H
similarity index 65%
rename from src/OpenFOAM/fields/Fields/complexFields/complexFields.H
rename to src/OpenFOAM/fields/Fields/complex/complexVectorField.H
index f937e656495b49349f91097a6480cdb524a5876c..80158aeb3a13e2ca42a05b9c55e84c9ba73778d5 100644
--- a/src/OpenFOAM/fields/Fields/complexFields/complexFields.H
+++ b/src/OpenFOAM/fields/Fields/complex/complexVectorField.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011 OpenFOAM Foundation
@@ -23,12 +23,6 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
-Typedef
-    Foam::complexField
-
-Description
-    Specialisation of Field\<T\> for complex.
-
 Typedef
     Foam::complexVectorField
 
@@ -36,16 +30,16 @@ Description
     Specialisation of Field\<T\> for complexVector.
 
 SourceFiles
-    complexFields.C
+    complexVectorField.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef complexFields_H
-#define complexFields_H
+#ifndef complexVectorField_H
+#define complexVectorField_H
 
 #include "complex.H"
 #include "complexVector.H"
-#include "primitiveFields.H"
+#include "vectorField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -54,31 +48,35 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-typedef Field<complex> complexField;
+typedef Field<complexVector> complexVectorField;
 
-complexField ComplexField(const UList<scalar>&, const UList<scalar>&);
-complexField ReComplexField(const UList<scalar>&);
-complexField ImComplexField(const UList<scalar>&);
-scalarField Re(const UList<complex>&);
-scalarField Im(const UList<complex>&);
-scalarField ReImSum(const UList<complex>&);
+//- Zip up two lists of values into a list of complex
+complexVectorField ComplexField
+(
+    const UList<vector>& re,
+    const UList<vector>& im
+);
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+//- Create complex field from a list of real (using imag == 0)
+complexVectorField ReComplexField(const UList<vector>& re);
 
-typedef Field<complexVector> complexVectorField;
+//- Create complex field from a list of imag (using real == 0)
+complexVectorField ImComplexField(const UList<vector>& im);
+
+//- Extract real component
+vectorField Re(const UList<complexVector>& cvf);
+
+//- Extract imag component
+vectorField Im(const UList<complexVector>& cvf);
 
-complexVectorField ComplexField(const UList<vector>&, const UList<vector>&);
-complexVectorField ReComplexField(const UList<vector>&);
-complexVectorField ImComplexField(const UList<vector>&);
-vectorField Re(const UList<complexVector>&);
-vectorField Im(const UList<complexVector>&);
-vectorField ReImSum(const UList<complexVector>&);
+//- Sum real and imag components
+vectorField ReImSum(const UList<complexVector>& cvf);
 
 complexVectorField operator^
 (
-    const UList<vector>&,
-    const UList<complexVector>&
+    const UList<vector>& vf,
+    const UList<complexVector>& cvf
 );
 
 
diff --git a/src/OpenFOAM/fields/Fields/complexFields/complexFields.C b/src/OpenFOAM/fields/Fields/complexFields/complexFields.C
deleted file mode 100644
index 604c8fe2fb4750894fe9d56764240679a77d3c2b..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/fields/Fields/complexFields/complexFields.C
+++ /dev/null
@@ -1,248 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           |
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-                            | Copyright (C) 2011 OpenFOAM Foundation
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Description
-    Specialisation of Field\<T\> for complex and complexVector.
-
-\*---------------------------------------------------------------------------*/
-
-#include "complexFields.H"
-#include "addToRunTimeSelectionTable.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-defineCompoundTypeName(List<complex>, complexList);
-addCompoundToRunTimeSelectionTable(List<complex>, complexList);
-
-complexField ComplexField(const UList<scalar>& re, const UList<scalar>& im)
-{
-    complexField cf(re.size());
-
-    forAll(cf, i)
-    {
-        cf[i].Re() = re[i];
-        cf[i].Im() = im[i];
-    }
-
-    return cf;
-}
-
-
-complexField ReComplexField(const UList<scalar>& sf)
-{
-    complexField cf(sf.size());
-
-    forAll(cf, i)
-    {
-        cf[i].Re() = sf[i];
-        cf[i].Im() = 0.0;
-    }
-
-    return cf;
-}
-
-
-complexField ImComplexField(const UList<scalar>& sf)
-{
-    complexField cf(sf.size());
-
-    forAll(cf, i)
-    {
-        cf[i].Re() = 0.0;
-        cf[i].Im() = sf[i];
-    }
-
-    return cf;
-}
-
-
-scalarField ReImSum(const UList<complex>& cf)
-{
-    scalarField sf(cf.size());
-
-    forAll(sf, i)
-    {
-        sf[i] = cf[i].Re() + cf[i].Im();
-    }
-
-    return sf;
-}
-
-
-scalarField Re(const UList<complex>& cf)
-{
-    scalarField sf(cf.size());
-
-    forAll(sf, i)
-    {
-        sf[i] = cf[i].Re();
-    }
-
-    return sf;
-}
-
-
-scalarField Im(const UList<complex>& cf)
-{
-    scalarField sf(cf.size());
-
-    forAll(sf, i)
-    {
-        sf[i] = cf[i].Im();
-    }
-
-    return sf;
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-defineCompoundTypeName(List<complexVector>, complexVectorList);
-addCompoundToRunTimeSelectionTable(List<complexVector>, complexVectorList);
-
-complexVectorField ComplexField
-(
-    const UList<vector>& re,
-    const UList<vector>& im
-)
-{
-    complexVectorField cvf(re.size());
-
-    for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
-    {
-        forAll(cvf, i)
-        {
-            cvf[i].component(cmpt).Re() = re[i].component(cmpt);
-            cvf[i].component(cmpt).Im() = im[i].component(cmpt);
-        }
-    }
-
-    return cvf;
-}
-
-
-complexVectorField ReComplexField(const UList<vector>& vf)
-{
-    complexVectorField cvf(vf.size());
-
-    for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
-    {
-        forAll(cvf, i)
-        {
-            cvf[i].component(cmpt).Re() = vf[i].component(cmpt);
-            cvf[i].component(cmpt).Im() = 0.0;
-        }
-    }
-
-    return cvf;
-}
-
-
-complexVectorField ImComplexField(const UList<vector>& vf)
-{
-    complexVectorField cvf(vf.size());
-
-    for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
-    {
-        forAll(cvf, i)
-        {
-            cvf[i].component(cmpt).Re() = 0.0;
-            cvf[i].component(cmpt).Im() = vf[i].component(cmpt);
-        }
-    }
-
-    return cvf;
-}
-
-
-vectorField ReImSum(const UList<complexVector>& cvf)
-{
-    vectorField vf(cvf.size());
-
-    for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
-    {
-        forAll(cvf, i)
-        {
-            vf[i].component(cmpt) =
-                cvf[i].component(cmpt).Re() + cvf[i].component(cmpt).Im();
-        }
-    }
-
-    return vf;
-}
-
-
-vectorField Re(const UList<complexVector>& cvf)
-{
-    vectorField vf(cvf.size());
-
-    for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
-    {
-        forAll(cvf, i)
-        {
-            vf[i].component(cmpt) = cvf[i].component(cmpt).Re();
-        }
-    }
-
-    return vf;
-}
-
-
-vectorField Im(const UList<complexVector>& cvf)
-{
-    vectorField vf(cvf.size());
-
-    for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
-    {
-        forAll(cvf, i)
-        {
-            vf[i].component(cmpt) = cvf[i].component(cmpt).Im();
-        }
-    }
-
-    return vf;
-}
-
-
-complexVectorField operator^
-(
-    const UList<vector>& vf,
-    const UList<complexVector>& cvf
-)
-{
-    return ComplexField(vf^Re(cvf), vf^Im(cvf));
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/PBiCGStab/PBiCGStab.C b/src/OpenFOAM/matrices/lduMatrix/solvers/PBiCGStab/PBiCGStab.C
index adf8b3054a49f2b900bc97127b642bcb0b0c998b..1e343565fb7ce81d4126db542fe2ee3b768e2389 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/PBiCGStab/PBiCGStab.C
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/PBiCGStab/PBiCGStab.C
@@ -208,7 +208,11 @@ Foam::solverPerformance Foam::PBiCGStab::solve
             solverPerf.finalResidual() =
                 gSumMag(sA, matrix().mesh().comm())/normFactor;
 
-            if (solverPerf.checkConvergence(tolerance_, relTol_))
+            if
+            (
+                solverPerf.nIterations() >= minIter_
+             && solverPerf.checkConvergence(tolerance_, relTol_)
+            )
             {
                 for (label cell=0; cell<nCells; cell++)
                 {
diff --git a/src/OpenFOAM/meshes/boundBox/boundBoxI.H b/src/OpenFOAM/meshes/boundBox/boundBoxI.H
index b762840bc2a27d465497724bf15ee4763dc64195..a71dba711a75241acbff01f88d6aab9b9daada8f 100644
--- a/src/OpenFOAM/meshes/boundBox/boundBoxI.H
+++ b/src/OpenFOAM/meshes/boundBox/boundBoxI.H
@@ -159,7 +159,7 @@ inline Foam::scalar Foam::boundBox::avgDim() const
 }
 
 
-Foam::label Foam::boundBox::nDim() const
+inline Foam::label Foam::boundBox::nDim() const
 {
     label ngood = 0;
 
diff --git a/src/OpenFOAM/primitives/Vector/complexVector/complexVector.H b/src/OpenFOAM/primitives/Vector/complexVector/complexVector.H
index cbcf718501d26f1a7fd1de2b5d0bd211a7bad4cb..488ea35c01dc5fd3c1d5ab87ecaeb3fe7b8a8617 100644
--- a/src/OpenFOAM/primitives/Vector/complexVector/complexVector.H
+++ b/src/OpenFOAM/primitives/Vector/complexVector/complexVector.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011 OpenFOAM Foundation
@@ -27,7 +27,7 @@ Typedef
     Foam::complexVector
 
 Description
-    complexVector obtained from generic Vector.
+    A Vector of complex values with 'scalar' precision.
 
 SourceFiles
     complexVectorI.H
@@ -48,6 +48,7 @@ namespace Foam
     typedef Vector<complex> complexVector;
 }
 
+// Functions
 #include "complexVectorI.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/OpenFOAM/primitives/Vector/complexVector/complexVectorI.H b/src/OpenFOAM/primitives/Vector/complexVector/complexVectorI.H
index 5a059aed07dab4d9c81bbc7f9b049f0b867c5140..4f257d65439b24b9b0c3da13abd942e3219accca 100644
--- a/src/OpenFOAM/primitives/Vector/complexVector/complexVectorI.H
+++ b/src/OpenFOAM/primitives/Vector/complexVector/complexVectorI.H
@@ -23,10 +23,6 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
-Description
-    complexVector specific part of 3D complexVector obtained from
-    generic Vector.
-
 \*---------------------------------------------------------------------------*/
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -34,7 +30,7 @@ Description
 namespace Foam
 {
 
-// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * Global Operators  * * * * * * * * * * * * * //
 
 inline complexVector operator*(const complex& v1, const complexVector& v2)
 {
@@ -80,8 +76,7 @@ inline complexVector operator/(const complex& v1, const complexVector& v2)
 }
 
 
-// complexVector dot product
-
+//- Dot product for complexVector
 inline complex operator&(const complexVector& v1, const complexVector& v2)
 {
     return complex
@@ -93,8 +88,7 @@ inline complex operator&(const complexVector& v1, const complexVector& v2)
 }
 
 
-// complexVector cross product
-
+//- Cross product for complexVector
 inline complexVector operator^(const complexVector& v1, const complexVector& v2)
 {
     return complexVector
@@ -106,8 +100,7 @@ inline complexVector operator^(const complexVector& v1, const complexVector& v2)
 }
 
 
-// complexVector cross product
-
+//- Cross product for complexVector
 inline complexVector operator^(const vector& v1, const complexVector& v2)
 {
     return complexVector
diff --git a/src/OpenFOAM/primitives/complex/complex.C b/src/OpenFOAM/primitives/complex/complex.C
index a85f596fdccd0b336d80ea8f35aaea85eba2b077..1368df2ae440e46ec5e7eab9c54f5b9f31ffe380 100644
--- a/src/OpenFOAM/primitives/complex/complex.C
+++ b/src/OpenFOAM/primitives/complex/complex.C
@@ -32,7 +32,7 @@ License
 
 const char* const Foam::complex::typeName = "complex";
 const Foam::complex Foam::complex::zero(0, 0);
-const Foam::complex Foam::complex::one(1, 1);
+const Foam::complex Foam::complex::one(1, 0);
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -55,12 +55,15 @@ Foam::word Foam::name(const complex& c)
 
 Foam::Istream& Foam::operator>>(Istream& is, complex& c)
 {
-    is.readBegin("complex");
-
-    is  >> c.re >> c.im;
+    scalar r, i;
 
+    is.readBegin("complex");
+    is >> r >> i;
     is.readEnd("complex");
 
+    c.real(r);
+    c.imag(i);
+
     is.check(FUNCTION_NAME);
     return is;
 }
@@ -69,7 +72,7 @@ Foam::Istream& Foam::operator>>(Istream& is, complex& c)
 Foam::Ostream& Foam::operator<<(Ostream& os, const complex& c)
 {
     os  << token::BEGIN_LIST
-        << c.re << token::SPACE << c.im
+        << c.real() << token::SPACE << c.imag()
         << token::END_LIST;
 
     return os;
diff --git a/src/OpenFOAM/primitives/complex/complex.H b/src/OpenFOAM/primitives/complex/complex.H
index c4b258a06ed42350ddaae2bb9d0372860c97cd82..48e3505d1935ec874d945127aa9b1014c55bf627 100644
--- a/src/OpenFOAM/primitives/complex/complex.H
+++ b/src/OpenFOAM/primitives/complex/complex.H
@@ -38,6 +38,7 @@ SourceFiles
 #ifndef complex_H
 #define complex_H
 
+#include <complex>
 #include "scalar.H"
 #include "word.H"
 #include "zero.H"
@@ -68,9 +69,6 @@ inline complex operator*(const complex&, const scalar);
 inline complex operator/(const complex&, const scalar);
 inline complex operator/(const scalar, const complex&);
 
-Istream& operator>>(Istream& is, complex& c);
-Ostream& operator<<(Ostream& os, const complex& c);
-
 
 /*---------------------------------------------------------------------------*\
                            Class complex Declaration
@@ -98,7 +96,7 @@ public:
         //- A complex zero (0,0)
         static const complex zero;
 
-        //- A complex one (1,1)
+        //- A complex one (1,0)
         static const complex one;
 
 
@@ -119,6 +117,12 @@ public:
         //- Construct from real and imaginary parts
         inline constexpr complex(const scalar r, const scalar i) noexcept;
 
+        //- Construct from std::complex
+        inline complex(const std::complex<float>& c);
+
+        //- Construct from std::complex
+        inline complex(const std::complex<double>& c);
+
         //- Construct from Istream
         explicit complex(Istream& is);
 
@@ -172,6 +176,13 @@ public:
 
     // Member Operators
 
+        //- Conversion to std::complex
+        inline operator std::complex<scalar>() const
+        {
+            return std::complex<scalar>(re, im);
+        }
+
+
         //- Copy assignment
         inline void operator=(const complex& c);
 
@@ -189,12 +200,6 @@ public:
         inline void operator*=(const scalar s);
         inline void operator/=(const scalar s);
 
-        //- Conjugate
-        inline complex operator~() const;
-
-        //- Conjugate
-        inline complex operator!() const;
-
         inline bool operator==(const complex& c) const;
         inline bool operator!=(const complex& c) const;
 
@@ -224,17 +229,19 @@ public:
         friend complex operator*(const complex& c, const scalar s);
         friend complex operator/(const complex& c, const scalar s);
         friend complex operator/(const scalar s, const complex& c);
+};
 
 
-    // IOstream Operators
+// * * * * * * * * * * * * * * * Global Operators  * * * * * * * * * * * * * //
 
-        friend Istream& operator>>(Istream& is, complex& c);
-        friend Ostream& operator<<(Ostream& os, const complex& c);
-};
+Istream& operator>>(Istream& is, complex& c);
+Ostream& operator<<(Ostream& os, const complex& c);
 
+//- Complex conjugate
+inline complex operator~(const complex& c);
 
-// * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * * //
 
+// * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * * //
 
 //- Return string representation of complex
 word name(const complex& c);
diff --git a/src/OpenFOAM/primitives/complex/complexI.H b/src/OpenFOAM/primitives/complex/complexI.H
index 7104ac84b7eb09fadae81d9cc1f4c0e9b8b938a0..25b6cc5a49236ab682683d152bdf00dcce4f18c6 100644
--- a/src/OpenFOAM/primitives/complex/complexI.H
+++ b/src/OpenFOAM/primitives/complex/complexI.H
@@ -55,6 +55,20 @@ inline constexpr Foam::complex::complex(const scalar r, const scalar i) noexcept
 {}
 
 
+inline Foam::complex::complex(const std::complex<float>& c)
+:
+    re(c.real()),
+    im(c.imag())
+{}
+
+
+inline Foam::complex::complex(const std::complex<double>& c)
+:
+    re(c.real()),
+    im(c.imag())
+{}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 inline void Foam::complex::real(scalar val)
@@ -174,33 +188,28 @@ inline void Foam::complex::operator/=(const scalar s)
 }
 
 
-inline Foam::complex Foam::complex::operator~() const
+inline bool Foam::complex::operator==(const complex& c) const
 {
-    return conjugate();
+    return (equal(re, c.re) && equal(im, c.im));
 }
 
 
-inline Foam::complex Foam::complex::operator!() const
+inline bool Foam::complex::operator!=(const complex& c) const
 {
-    return conjugate();
+    return !operator==(c);
 }
 
 
-inline bool Foam::complex::operator==(const complex& c) const
-{
-    return (equal(re, c.re) && equal(im, c.im));
-}
+// * * * * * * * * * * * * * * * Global Operators  * * * * * * * * * * * * * //
 
-
-inline bool Foam::complex::operator!=(const complex& c) const
+inline Foam::complex Foam::operator~(const complex& c)
 {
-    return !operator==(c);
+    return c.conjugate();
 }
 
 
 // * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
 
-
 namespace Foam
 {
 
@@ -341,9 +350,97 @@ inline complex operator/(const scalar s, const complex& c)
     return complex(s/c.re, s/c.im);
 }
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+// Complex transcendental functions
+namespace Foam
+{
+    #define transFunc(func)                                                   \
+    inline complex func(const Foam::complex& z)                               \
+    {                                                                         \
+        return std:: func (std::complex<scalar>(z));                          \
+    }
+
+transFunc(sqrt)
+transFunc(exp)
+transFunc(log)
+transFunc(log10)
+transFunc(sin)
+transFunc(cos)
+transFunc(tan)
+transFunc(asin)
+transFunc(acos)
+transFunc(atan)
+transFunc(sinh)
+transFunc(cosh)
+transFunc(tanh)
+transFunc(asinh)
+transFunc(acosh)
+transFunc(atanh)
+
+// Special treatment for pow()
+inline complex pow(const complex& x, const complex& y)
+{
+    return std::pow(std::complex<scalar>(x), std::complex<scalar>(y));
+}
+
+
+// Combinations of complex and real
+#define powFuncs(type2)                                                       \
+    inline complex pow(const complex& x, const type2& y)                      \
+    {                                                                         \
+        return std::pow(std::complex<scalar>(x), scalar(y));                  \
+    }                                                                         \
+                                                                              \
+    inline Foam::complex pow(const type2& x, const complex& y)                \
+    {                                                                         \
+        return std::pow(scalar(x), std::complex<scalar>(y));                  \
+    }
+
+
+powFuncs(float)
+powFuncs(double)
+powFuncs(int)
+powFuncs(long)
+
+
+inline complex pow3(const complex& c)
+{
+    return c*sqr(c);
+}
+
+
+inline complex pow4(const complex& c)
+{
+    return sqr(sqr(c));
+}
+
+
+inline complex pow5(const complex& c)
+{
+    return c*pow4(c);
+}
+
+
+inline complex pow6(const complex& c)
+{
+    return pow3(sqr(c));
+}
+
+
+inline complex pow025(const complex& c)
+{
+    return sqrt(sqrt(c));
+}
+
 } // End namespace Foam
 
+#undef transFunc
+#undef powFuncs
+
 // ************************************************************************* //
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.C
index ca9427f357a932e7bdc21bea027cac789fe539a4..67551474268068c8ed10a5dc432a9d5212371c32 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.C
@@ -129,6 +129,7 @@ bool Foam::combineFaces::validFace
 void Foam::combineFaces::regioniseFaces
 (
     const scalar minCos,
+    const bool mergeAcrossPatches,
     const label celli,
     const labelList& cEdges,
     Map<label>& faceRegion
@@ -143,16 +144,31 @@ void Foam::combineFaces::regioniseFaces
         label f0, f1;
         meshTools::getEdgeFaces(mesh_, celli, edgeI, f0, f1);
 
+        const vector& a0 = mesh_.faceAreas()[f0];
+        const vector& a1 = mesh_.faceAreas()[f1];
+
         const label p0 = patches.whichPatch(f0);
         const label p1 = patches.whichPatch(f1);
 
         // Face can be merged if
-        // - same non-coupled patch
         // - small angle
-        if (p0 != -1 && p0 == p1 && !patches[p0].coupled())
+        // - mergeAcrossPatches=false : same non-coupled patch
+        // - mergeAcrossPatches=true  : always
+        if
+        (
+            p0 != -1
+         && p1 != -1
+         && !patches[p0].coupled()
+         && !patches[p1].coupled()
+        )
         {
-            const vector f0Normal = normalised(mesh_.faceAreas()[f0]);
-            const vector f1Normal = normalised(mesh_.faceAreas()[f1]);
+            if (!mergeAcrossPatches && (p0 != p1))
+            {
+                continue;
+            }
+
+            const vector f0Normal = normalised(a0);
+            const vector f1Normal = normalised(a1);
 
             if ((f0Normal & f1Normal) > minCos)
             {
@@ -285,7 +301,8 @@ Foam::labelListList Foam::combineFaces::getMergeSets
 (
     const scalar featureCos,
     const scalar minConcaveCos,
-    const labelHashSet& boundaryCells
+    const labelHashSet& boundaryCells,
+    const bool mergeAcrossPatches
 ) const
 {
     // Lists of faces that can be merged.
@@ -303,7 +320,14 @@ Foam::labelListList Foam::combineFaces::getMergeSets
 
         // Region per face
         Map<label> faceRegion(cFaces.size());
-        regioniseFaces(featureCos, celli, cEdges, faceRegion);
+        regioniseFaces
+        (
+            featureCos,
+            mergeAcrossPatches,
+            celli,
+            cEdges,
+            faceRegion
+        );
 
         // Now we have in faceRegion for every face the region with planar
         // face sharing the same region. We now check whether the resulting
@@ -338,7 +362,7 @@ Foam::labelListList Foam::combineFaces::getMergeSets
 
             // For every set check if it forms a valid convex face
 
-            forAllConstIters(regionToFaces, iter)
+            forAllIters(regionToFaces, iter)
             {
                 // Make face out of setFaces
                 indirectPrimitivePatch bigFace
@@ -354,7 +378,33 @@ Foam::labelListList Foam::combineFaces::getMergeSets
                 // Only store if -only one outside loop -which forms convex face
                 if (validFace(minConcaveCos, bigFace))
                 {
-                    allFaceSets.append(iter.val());
+                    labelList& faceIDs = iter.val();
+
+                    // For cross-patch merging we want to make the
+                    // largest face the one to decide the final patch
+                    // (i.e. master face)
+                    if (mergeAcrossPatches)
+                    {
+                        const vectorField& areas = mesh_.faceAreas();
+
+                        label maxIndex = 0;
+                        scalar maxMagSqr = magSqr(areas[faceIDs[0]]);
+                        for (label i = 1; i < faceIDs.size(); ++i)
+                        {
+                            const scalar a2 = magSqr(areas[faceIDs[i]]);
+                            if (a2 > maxMagSqr)
+                            {
+                                maxMagSqr = a2;
+                                maxIndex = i;
+                            }
+                        }
+                        if (maxIndex != 0)
+                        {
+                            Swap(faceIDs[0], faceIDs[maxIndex]);
+                        }
+                    }
+
+                    allFaceSets.append(faceIDs);
                 }
             }
         }
@@ -367,7 +417,8 @@ Foam::labelListList Foam::combineFaces::getMergeSets
 Foam::labelListList Foam::combineFaces::getMergeSets
 (
     const scalar featureCos,
-    const scalar minConcaveCos
+    const scalar minConcaveCos,
+    const bool mergeAcrossPatches
 ) const
 {
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
@@ -388,7 +439,13 @@ Foam::labelListList Foam::combineFaces::getMergeSets
         }
     }
 
-    return getMergeSets(featureCos, minConcaveCos, boundaryCells);
+    return getMergeSets
+    (
+        featureCos,
+        minConcaveCos,
+        boundaryCells,
+        mergeAcrossPatches
+    );
 }
 
 
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.H b/src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.H
index 2eb86e87ad16e75b9000f0d33eff2e9009e376cc..851def914ea051a2376983d23072d42f6d5219de 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.H
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C)2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2016 OpenFOAM Foundation
@@ -103,6 +103,7 @@ class combineFaces
         void regioniseFaces
         (
             const scalar minCos,
+            const bool mergeAcrossPatches,
             const label celli,
             const labelList& cEdges,
             Map<label>& faceRegion
@@ -155,20 +156,27 @@ public:
         // Helper functions
 
             //- Extract lists of all (non-coupled) boundary faces on selected
-            //  cells that can be merged. Uses getFaceRegions.
+            //  cells that can be merged. Uses getFaceRegions. Optionally
+            //  allow faces-on-different-patches to be merged (into the largest
+            //  area face - could be improved). Note: causes a problem in
+            //  undoing - all restored faces get the patch/zone from the
+            //  master face.
             labelListList getMergeSets
             (
                 const scalar featureCos,
                 const scalar minConcaveCos,
-                const labelHashSet& boundaryCells
+                const labelHashSet& boundaryCells,
+                const bool mergeAcrossPatches = false
             ) const;
 
             //- Extract lists of all (non-coupled) boundary faces that can
-            //  be merged. Uses getFaceRegions.
+            //  be merged. Uses getFaceRegions. See note above about
+            //  mergeAcrossPatches.
             labelListList getMergeSets
             (
                 const scalar featureCos,
-                const scalar minConcaveCos
+                const scalar minConcaveCos,
+                const bool mergeAcrossPatches = false
             ) const;
 
             //- Gets outside of patch as a face (in mesh point labels)
@@ -197,7 +205,8 @@ public:
             //  Returns maps from added restored point to
             //  original point label (i.e. content of savedPointLabels_).
             //  (only restoredPoints are actually set; rest are just for
-            //   generalness)
+            //   generalness). See note above about restoring faces from
+            //  different patches (mergeAcrossPatches)
             void setUnrefinement
             (
                 const labelList& masterFaces,
diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/limitedSnGrad/limitedSnGrad.C b/src/finiteVolume/finiteVolume/snGradSchemes/limitedSnGrad/limitedSnGrad.C
index c4586a74f13e774d6bad94fa892a7f1076b6e803..8f435086592c60dab86496adb8630db61b41d426 100644
--- a/src/finiteVolume/finiteVolume/snGradSchemes/limitedSnGrad/limitedSnGrad.C
+++ b/src/finiteVolume/finiteVolume/snGradSchemes/limitedSnGrad/limitedSnGrad.C
@@ -30,6 +30,7 @@ License
 #include "volFields.H"
 #include "surfaceFields.H"
 #include "localMax.H"
+#include "fvcCellReduce.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -82,6 +83,40 @@ limitedSnGrad<Type>::correction
             << "limiter min: " << min(limiter.primitiveField())
             << " max: "<< max(limiter.primitiveField())
             << " avg: " << average(limiter.primitiveField()) << endl;
+
+
+        if (fv::debug & 2)
+        {
+            static scalar oldTime = -1;
+            static label subIter = 0;
+            if (vf.mesh().time().value() != oldTime)
+            {
+                oldTime = vf.mesh().time().value();
+                subIter = 0;
+            }
+            else
+            {
+                ++subIter;
+            }
+            word fieldName("limiter_" + Foam::name(subIter));
+
+            GeometricField<scalar, fvPatchField, volMesh> volLimiter
+            (
+                IOobject
+                (
+                    fieldName,
+                    vf.mesh().time().timeName(),
+                    vf.mesh(),
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                fvc::cellReduce(limiter, minEqOp<scalar>(), scalar(1.0))
+            );
+            Info<< "Writing limiter field to " << volLimiter.objectPath()
+                << endl;
+            volLimiter.write();
+        }
     }
 
     return limiter*corr;
diff --git a/src/mesh/blockMesh/PDRblockMesh/PDRblock.H b/src/mesh/blockMesh/PDRblockMesh/PDRblock.H
index 91e16b66e28195fd0bbdb7f8d085787f24cfd7c5..5f85f21ca9d57796cf33d13975e7370bc1dd9614 100644
--- a/src/mesh/blockMesh/PDRblockMesh/PDRblock.H
+++ b/src/mesh/blockMesh/PDRblockMesh/PDRblock.H
@@ -110,6 +110,12 @@ public:
             //- True if the location is within the range
             inline bool contains(const scalar p) const;
 
+            //- The first() value is considered the min value.
+            inline const scalar& min() const;
+
+            //- The last() value is considered the max value.
+            inline const scalar& max() const;
+
             //- Mid-point location, zero for an empty list.
             inline scalar centre() const;
 
@@ -305,6 +311,12 @@ public:
         //- Cell dimensions at i,j,k position.
         inline vector span(const labelVector& ijk) const;
 
+        //- Grid point at i,j,k position.
+        inline point grid(const label i, const label j, const label k) const;
+
+        //- Grid point at i,j,k position.
+        inline point grid(const labelVector& ijk) const;
+
         //- Cell centre at i,j,k position.
         inline point C(const label i, const label j, const label k) const;
 
diff --git a/src/mesh/blockMesh/PDRblockMesh/PDRblockI.H b/src/mesh/blockMesh/PDRblockMesh/PDRblockI.H
index a624ee68a650bde8e6ce62ce4e6f11da076bd25a..7a04ce09ea09dd3c4df179d6386ee5d1af5644ed 100644
--- a/src/mesh/blockMesh/PDRblockMesh/PDRblockI.H
+++ b/src/mesh/blockMesh/PDRblockMesh/PDRblockI.H
@@ -62,6 +62,18 @@ inline bool Foam::PDRblock::location::contains(const scalar p) const
 }
 
 
+inline const Foam::scalar& Foam::PDRblock::location::min() const
+{
+    return scalarList::empty() ? pTraits<scalar>::rootMax : first();
+}
+
+
+inline const Foam::scalar& Foam::PDRblock::location::max() const
+{
+    return scalarList::empty() ? pTraits<scalar>::rootMin : last();
+}
+
+
 inline Foam::scalar Foam::PDRblock::location::centre() const
 {
     return scalarList::empty() ? 0 : (0.5*first() + 0.5*last());
@@ -206,6 +218,29 @@ inline Foam::vector Foam::PDRblock::span(const labelVector& ijk) const
 }
 
 
+inline Foam::point Foam::PDRblock::grid
+(
+    const label i,
+    const label j,
+    const label k
+) const
+{
+    return point(grid_.x()[i], grid_.y()[j], grid_.z()[k]);
+}
+
+
+inline Foam::point Foam::PDRblock::grid(const labelVector& ijk) const
+{
+    return
+        point
+        (
+            grid_.x()[ijk.x()],
+            grid_.y()[ijk.y()],
+            grid_.z()[ijk.z()]
+        );
+}
+
+
 inline Foam::point Foam::PDRblock::C
 (
     const label i,
diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H
index 56067d2b9ff847d8b26c6d4f5f6f151816af53b2..98c8650c2c43c748dba8e1a1e06f31ecc54255bf 100644
--- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H
+++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2015-2019 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2015-2017 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2017 OpenFOAM Foundation
@@ -129,6 +129,16 @@ public:
             REMOVE = 4      //!< set value to -1 any face that was refined
         };
 
+        //- Enumeration for what to do with co-planar patch faces on a single
+        //  cell
+        enum FaceMergeType
+        {
+            NONE,           // no merging
+            GEOMETRIC,      // use feature angle
+            IGNOREPATCH     // use feature angle, allow merging of different
+                            // patches
+        };
+
 
 private:
 
@@ -1463,7 +1473,8 @@ public:
                     const scalar minCos,
                     const scalar concaveCos,
                     const label mergeSize,
-                    const labelList& patchIDs
+                    const labelList& patchIDs,
+                    const meshRefinement::FaceMergeType mergeType
                 );
 
                 //- Merge coplanar faces. preserveFaces is != -1 for faces
@@ -1474,7 +1485,8 @@ public:
                     const scalar concaveCos,
                     const labelList& patchIDs,
                     const dictionary& motionDict,
-                    const labelList& preserveFaces
+                    const labelList& preserveFaces,
+                    const meshRefinement::FaceMergeType mergeType
                 );
 
                 autoPtr<mapPolyMesh> doRemovePoints
diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementMerge.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementMerge.C
index 99c7faa7aad5f03b4784b6e54f789ea458000bbb..6e3479fad0119fa37f6c0471a60992773e9214c1 100644
--- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementMerge.C
+++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementMerge.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016-2018 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2016-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2014 OpenFOAM Foundation
@@ -42,7 +42,8 @@ Foam::label Foam::meshRefinement::mergePatchFaces
     const scalar minCos,
     const scalar concaveCos,
     const label mergeSize,
-    const labelList& patchIDs
+    const labelList& patchIDs,
+    const meshRefinement::FaceMergeType mergeType
 )
 {
     // Patch face merging engine
@@ -73,7 +74,8 @@ Foam::label Foam::meshRefinement::mergePatchFaces
         (
             minCos,
             concaveCos,
-            boundaryCells
+            boundaryCells,
+            (mergeType == FaceMergeType::IGNOREPATCH) // merge across patches?
         )
     );
 
@@ -249,7 +251,8 @@ Foam::label Foam::meshRefinement::mergePatchFacesUndo
     const scalar concaveCos,
     const labelList& patchIDs,
     const dictionary& motionDict,
-    const labelList& preserveFaces
+    const labelList& preserveFaces,
+    const meshRefinement::FaceMergeType mergeType
 )
 {
     // Patch face merging engine
@@ -286,7 +289,8 @@ Foam::label Foam::meshRefinement::mergePatchFacesUndo
         (
             minCos,
             concaveCos,
-            boundaryCells
+            boundaryCells,
+            (mergeType == FaceMergeType::IGNOREPATCH) // merge across patches?
         )
     );
 
diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.C b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.C
index e67f69fd58263a110b72e960ca372914d239129d..6518115a1331057dcc1d4e9af16138802e2ff24e 100644
--- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.C
+++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.C
@@ -3219,7 +3219,8 @@ Foam::snappyLayerDriver::snappyLayerDriver
 void Foam::snappyLayerDriver::mergePatchFacesUndo
 (
     const layerParameters& layerParams,
-    const dictionary& motionDict
+    const dictionary& motionDict,
+    const meshRefinement::FaceMergeType mergeType
 )
 {
     // Clip to 30 degrees. Not helpful!
@@ -3260,7 +3261,8 @@ void Foam::snappyLayerDriver::mergePatchFacesUndo
         concaveCos,
         meshRefiner_.meshedPatches(),
         motionDict,
-        duplicateFace
+        duplicateFace,
+        mergeType   // How to merge co-planar patch faces
     );
 
     nChanged += meshRefiner_.mergeEdgesUndo(minCos, motionDict);
@@ -4635,7 +4637,7 @@ void Foam::snappyLayerDriver::doLayers
     const dictionary& shrinkDict,
     const dictionary& motionDict,
     const layerParameters& layerParams,
-    const bool mergePatchFaces,
+    const meshRefinement::FaceMergeType mergeType,
     const bool preBalance,
     decompositionMethod& decomposer,
     fvMeshDistribute& distributor
@@ -4653,9 +4655,13 @@ void Foam::snappyLayerDriver::doLayers
     Info<< "Using mesh parameters " << motionDict << nl << endl;
 
     // Merge coplanar boundary faces
-    if (mergePatchFaces)
+    if
+    (
+        mergeType == meshRefinement::FaceMergeType::GEOMETRIC
+     || mergeType == meshRefinement::FaceMergeType::IGNOREPATCH
+    )
     {
-        mergePatchFacesUndo(layerParams, motionDict);
+        mergePatchFacesUndo(layerParams, motionDict, mergeType);
     }
 
 
diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.H b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.H
index db8af3985a16f589ced255453f0c7e7837c387e3..231837ab5f27a8bf60714e2b88743248dac9602b 100644
--- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.H
+++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.H
@@ -634,7 +634,8 @@ public:
             void mergePatchFacesUndo
             (
                 const layerParameters& layerParams,
-                const dictionary& motionDict
+                const dictionary& motionDict,
+                const meshRefinement::FaceMergeType mergeType
             );
 
             //- Add cell layers
@@ -654,12 +655,11 @@ public:
                 const dictionary& shrinkDict,
                 const dictionary& motionDict,
                 const layerParameters& layerParams,
-                const bool mergePatchFaces,         // merging patch faces
+                const meshRefinement::FaceMergeType mergeType,
                 const bool preBalance,              // balance before adding?
                 decompositionMethod& decomposer,
                 fvMeshDistribute& distributor
             );
-
 };
 
 
diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C
index 75ba0a0f155ce5ed6ac022ae3dead97705a88f7c..cde402e87224c8caffe0afa487246d37a5158e29 100644
--- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C
+++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C
@@ -2794,7 +2794,7 @@ void Foam::snappyRefineDriver::addFaceZones
 
 void Foam::snappyRefineDriver::mergePatchFaces
 (
-    const bool geometricMerge,
+    const meshRefinement::FaceMergeType mergeType,
     const refinementParameters& refineParams,
     const dictionary& motionDict
 )
@@ -2812,7 +2812,11 @@ void Foam::snappyRefineDriver::mergePatchFaces
 
     const fvMesh& mesh = meshRefiner_.mesh();
 
-    if (geometricMerge)
+    if
+    (
+        mergeType == meshRefinement::FaceMergeType::GEOMETRIC
+     || mergeType == meshRefinement::FaceMergeType::IGNOREPATCH
+    )
     {
         meshRefiner_.mergePatchFacesUndo
         (
@@ -2820,7 +2824,8 @@ void Foam::snappyRefineDriver::mergePatchFaces
             Foam::cos(degToRad(45.0)),
             meshRefiner_.meshedPatches(),
             motionDict,
-            labelList(mesh.nFaces(), -1)
+            labelList(mesh.nFaces(), -1),
+            mergeType
         );
     }
     else
@@ -2831,7 +2836,8 @@ void Foam::snappyRefineDriver::mergePatchFaces
             Foam::cos(degToRad(45.0)),
             Foam::cos(degToRad(45.0)),
             4,          // only merge faces split into 4
-            meshRefiner_.meshedPatches()
+            meshRefiner_.meshedPatches(),
+            meshRefinement::FaceMergeType::GEOMETRIC // no merge across patches
         );
     }
 
@@ -2855,7 +2861,7 @@ void Foam::snappyRefineDriver::doRefine
     const refinementParameters& refineParams,
     const snapParameters& snapParams,
     const bool prepareForSnapping,
-    const bool doMergePatchFaces,
+    const meshRefinement::FaceMergeType mergeType,
     const dictionary& motionDict
 )
 {
@@ -3063,7 +3069,7 @@ void Foam::snappyRefineDriver::doRefine
     // Do something about cells with refined faces on the boundary
     if (prepareForSnapping)
     {
-        mergePatchFaces(doMergePatchFaces, refineParams, motionDict);
+        mergePatchFaces(mergeType, refineParams, motionDict);
     }
 
 
diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.H b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.H
index eeec7ed89abfdef23d169e9260fbffc49987aff1..da43baf51f0b71ca9b0ced2a4c5edb58cd92cc52 100644
--- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.H
+++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.H
@@ -43,6 +43,7 @@ SourceFiles
 #include "writer.H"
 #include "DynamicList.H"
 #include "labelVector.H"
+#include "meshRefinement.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -53,7 +54,6 @@ namespace Foam
 class refinementParameters;
 class snapParameters;
 
-class meshRefinement;
 class decompositionMethod;
 class fvMeshDistribute;
 class fvMesh;
@@ -226,7 +226,7 @@ class snappyRefineDriver
         //- Merge refined boundary faces (from exposing coarser cell)
         void mergePatchFaces
         (
-            const bool geometricMerge,
+            const meshRefinement::FaceMergeType mergeType,
             const refinementParameters& refineParams,
             const dictionary& motionDict
         );
@@ -268,7 +268,7 @@ public:
             const refinementParameters& refineParams,
             const snapParameters& snapParams,
             const bool prepareForSnapping,
-            const bool mergePatchFaces,
+            const meshRefinement::FaceMergeType mergeType,
             const dictionary& motionDict
         );
 
diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.C b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.C
index a30549e7656da547a7aae1fb880ff414213ff373..7f2ff128e8bf2c56c0a7c208f8c915fded36a038 100644
--- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.C
+++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.C
@@ -2528,7 +2528,7 @@ void Foam::snappySnapDriver::doSnap
 (
     const dictionary& snapDict,
     const dictionary& motionDict,
-    const bool mergePatchFaces,
+    const meshRefinement::FaceMergeType mergeType,
     const scalar featureCos,
     const scalar planarAngle,
     const snapParameters& snapParams
@@ -3031,7 +3031,11 @@ void Foam::snappySnapDriver::doSnap
         repatchToSurface(snapParams, adaptPatchIDs, duplicateFace);
     }
 
-    if (mergePatchFaces)
+    if
+    (
+        mergeType == meshRefinement::FaceMergeType::GEOMETRIC
+     || mergeType == meshRefinement::FaceMergeType::IGNOREPATCH
+    )
     {
         labelList duplicateFace(getInternalOrBaffleDuplicateFace());
 
@@ -3044,7 +3048,8 @@ void Foam::snappySnapDriver::doSnap
             featureCos,     // concaveCos
             meshRefiner_.meshedPatches(),
             motionDict,
-            duplicateFace   // faces not to merge
+            duplicateFace,  // faces not to merge
+            mergeType
         );
 
         nChanged += meshRefiner_.mergeEdgesUndo(featureCos, motionDict);
diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.H b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.H
index c1175df7c05a52fe4fd61b4d8c7301aca65012c1..c459a2b47449da2aeef8fc2c866db7460f390cb6 100644
--- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.H
+++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.H
@@ -781,12 +781,11 @@ public:
             (
                 const dictionary& snapDict,
                 const dictionary& motionDict,
-                const bool mergePatchFaces,
+                const meshRefinement::FaceMergeType mergeType,
                 const scalar featureCos,
                 const scalar planarAngle,
                 const snapParameters& snapParams
             );
-
 };
 
 
diff --git a/tutorials/basic/overLaplacianDyMFoam/heatTransfer/system/fvSchemes b/tutorials/basic/overLaplacianDyMFoam/heatTransfer/system/fvSchemes
index aaf2392bd190561b549ceb29b25b81d56e70a03a..ac3a013d5e0cc8888c3acafbd627cdb2b7219b10 100644
--- a/tutorials/basic/overLaplacianDyMFoam/heatTransfer/system/fvSchemes
+++ b/tutorials/basic/overLaplacianDyMFoam/heatTransfer/system/fvSchemes
@@ -74,9 +74,13 @@ oversetInterpolation
     //searchBoxDivisions  (100 100 1);
 }
 
-oversetInterpolationRequired
+
+// Any additional fields that should not be interpolated
+oversetInterpolationSuppresed
 {
-    // Any additional fields that require overset interpolation
+    // For backwards compatibility; should not really be suppressed
+    grad(T);
 }
 
+
 // ************************************************************************* //
diff --git a/tutorials/basic/overPotentialFoam/cylinder/cylinderAndBackground/system/fvSchemes b/tutorials/basic/overPotentialFoam/cylinder/cylinderAndBackground/system/fvSchemes
index e25ea6f5fcdd60a3d96dd8700b59a37f0b6ee808..7a4160973b67ffa684785d943a6151c90b3e4f91 100644
--- a/tutorials/basic/overPotentialFoam/cylinder/cylinderAndBackground/system/fvSchemes
+++ b/tutorials/basic/overPotentialFoam/cylinder/cylinderAndBackground/system/fvSchemes
@@ -53,5 +53,4 @@ oversetInterpolation
 }
 
 
-
 // ************************************************************************* //
diff --git a/tutorials/incompressible/pisoFoam/LES/motorBike/lesFiles/controlDict b/tutorials/incompressible/pisoFoam/LES/motorBike/lesFiles/controlDict
index b8096e514db550d8b434919d21e63a6b6d19a639..283724bc2d63c640a67f1d991628970012e79e35 100644
--- a/tutorials/incompressible/pisoFoam/LES/motorBike/lesFiles/controlDict
+++ b/tutorials/incompressible/pisoFoam/LES/motorBike/lesFiles/controlDict
@@ -49,7 +49,7 @@ runTimeModifiable true;
 
 functions
 {
-    #include "cuttingPlane"
+    #include "samples"
     #include "streamLines"
     #include "forceCoeffs"
 }
diff --git a/wmake/makefiles/general b/wmake/makefiles/general
index 6b7fc6e4b6ce8d7d47b8d26d71d56c341aee902f..198e5c1e59a3eb6c7b56fca8a0be5d902fbd127a 100644
--- a/wmake/makefiles/general
+++ b/wmake/makefiles/general
@@ -31,7 +31,7 @@ endif
 
 
 #------------------------------------------------------------------------------
-# Unset suffices list (suffix rules are not used)
+# No default suffix rules used
 #------------------------------------------------------------------------------
 
 .SUFFIXES:
diff --git a/wmake/makefiles/info b/wmake/makefiles/info
new file mode 100644
index 0000000000000000000000000000000000000000..90fa2b550df45c70030d48dda837ca346a36ffc4
--- /dev/null
+++ b/wmake/makefiles/info
@@ -0,0 +1,74 @@
+#----------------------------*- makefile-gmake -*------------------------------
+# =========                 |
+# \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+#  \\    /   O peration     |
+#   \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
+#    \\/     M anipulation  |
+#------------------------------------------------------------------------------
+# License
+#     This file is part of OpenFOAM, licensed under GNU General Public License
+#     <http://www.gnu.org/licenses/>.
+#
+# File
+#     wmake/makefiles/info
+#
+# Description
+#     Makefile to generate information.
+#     Used by wmake -show-*
+#
+#------------------------------------------------------------------------------
+
+#------------------------------------------------------------------------------
+# Use POSIX shell. Default to POSIX for the OS.
+#------------------------------------------------------------------------------
+
+SHELL   = /bin/sh
+
+#------------------------------------------------------------------------------
+# No default suffix rules used
+#------------------------------------------------------------------------------
+
+.SUFFIXES:
+
+
+#------------------------------------------------------------------------------
+# Compilation rules
+#------------------------------------------------------------------------------
+
+GENERAL_RULES = $(WM_DIR)/rules/General
+include $(GENERAL_RULES)/general
+
+
+#------------------------------------------------------------------------------
+# Display information
+#------------------------------------------------------------------------------
+
+export WM_VERSION
+
+
+.PHONY: compile
+compile:
+	@echo "$(strip $(CC) $(c++FLAGS))"
+
+.PHONY: api
+api:
+	@echo "$${WM_VERSION#*=}"
+
+.PHONY: c
+c:
+	@echo "$(strip $(cc))"
+
+.PHONY: cflags
+cflags:
+	@echo "$(strip $(cFLAGS))"
+
+.PHONY: cxx
+cxx:
+	@echo "$(strip $(CC))"
+
+.PHONY: cxxflags
+cxxflags:
+	@echo "$(strip $(c++FLAGS))"
+
+
+#----------------------------- vim: set ft=make: ------------------------------
diff --git a/wmake/wmake b/wmake/wmake
index c05e25c3f1d73af7b718022433487d5c4ebdd562..7591fc94512939c7335ca0f287a49311cbcd3925 100755
--- a/wmake/wmake
+++ b/wmake/wmake
@@ -77,6 +77,11 @@ options:
   -pwd              Print root directory containing a Make/ directory and exit
   -update           Update lnInclude directories, dep files, remove deprecated
                     files and directories
+  -show             Identical to -show-compile
+  -show-api         Print api value and exit
+  -show-compile     Print C++ compiler value/flags and exit
+  -show-cxx         Print C++ compiler value and exit
+  -show-cxxflags    Print C++ compiler flags and exit
   -h | -help        Print the usage
 
 
@@ -121,7 +126,7 @@ allCores()
 #------------------------------------------------------------------------------
 
 # Default to compiling the local target only
-unset all update optPrintRootDir
+unset all optShow update optPrintRootDir
 
 while [ "$#" -gt 0 ]
 do
@@ -133,6 +138,14 @@ do
         -s | -silent)
             export WM_QUIET=true
             ;;
+        -show | -show-compile)
+            $make -f $WM_DIR/makefiles/info compile
+            optShow=true
+            ;;
+        -show-api | -show-cxx | -show-cxxflags | -show-c | -show-cflags)
+            $make -f $WM_DIR/makefiles/info "${1#-show-}"
+            optShow=true
+            ;;
         -a | -all | all)
             all=all
             ;;
@@ -147,12 +160,12 @@ do
 
             [ "$nCores" = 0 ] && allCores
             export WM_NCOMPPROCS=$nCores
-            echo "Compiling enabled on $WM_NCOMPPROCS cores"
+            echo "Compiling enabled on $WM_NCOMPPROCS cores" 1>&2
             ;;
         # Parallel compilation on specified number of cores
         -j[1-9]*)
             export WM_NCOMPPROCS=${1#-j}
-            echo "Compiling enabled on $WM_NCOMPPROCS cores"
+            echo "Compiling enabled on $WM_NCOMPPROCS cores" 1>&2
             ;;
         # Keep going, ignoring errors
         -k | -keep-going | -non-stop)
@@ -190,6 +203,11 @@ do
     shift
 done
 
+if [ "$optShow" = true ]
+then
+    exit 0
+fi
+
 
 #------------------------------------------------------------------------------
 # Check environment variables