ENH: Added relaxed nonOrthogonal laplacian corrector scheme

- had calculated boundaries (default) for the evaluated variables,
  which meant they retained their initialized values (usually Zero).

  This normally goes unnoticed, since the boundary values are largely
  irrelevant in the volField expressions. However, when applying
  functions that balk at a zero value - eg, log() - this raises a
  floating point exception.

  These boundary should be zeroGradient, since the evaluated variables
  correspond to the internalField only. Could continue to use
  calculated, but then need to set the calculated boundary values from
  the patch internal field manually.

Allwmake

@@ -80,18 +80,6 @@ then
     echo
 elif [ -d "$WM_PROJECT_DIR/modules" ]
 then
-    echo "========================================"
-    echo "Compile OpenFOAM modules"
-    echo "    ignoring possible compilation errors"
-    echo "    make certain to check the output file"
-    echo
-    set +e
-    export WM_CONTINUE_ON_ERROR=true
-
-    # Default build into OpenFOAM project locations
-    : "${FOAM_MODULE_PREFIX:=${FOAM_LIBBIN%/*}}"
-    export FOAM_MODULE_PREFIX
-
     (cd "$WM_PROJECT_DIR/modules" 2>/dev/null && wmake -all)
 fi
 

META-INFO/api-info

@@ -1,2 +1,2 @@
 api=2006
-patch=0
+patch=201012

applications/utilities/mesh/generation/foamyMesh/foamyHexMesh/Make/options

@@ -26,6 +26,8 @@ EXE_INC = \
 EXE_LIBS = \
     ${CGAL_LIBS} \
     -lconformalVoronoiMesh \
+    -lfileFormats \
+    -lsurfMesh \
     -lmeshTools \
     -ldecompositionMethods \
     -ldecompose \

applications/utilities/mesh/manipulation/checkMesh/writeFields.C

@@ -616,6 +616,17 @@ void Foam::writeFields
     }
     if (selectedFields.found("faceZone"))
     {
+        // Determine for each face the zone index (scalar for ease of
+        // manipulation)
+        scalarField zoneID(mesh.nFaces(), -1);
+        const faceZoneMesh& czs = mesh.faceZones();
+        for (const auto& zone : czs)
+        {
+            UIndirectList<scalar>(zoneID, zone) = zone.index();
+        }
+
+
+        // Split into internal and boundary values
         surfaceScalarField faceZone
         (
             IOobject
@@ -632,10 +643,13 @@ void Foam::writeFields
             calculatedFvsPatchScalarField::typeName
         );
 
-        const faceZoneMesh& czs = mesh.faceZones();
-        for (const auto& zone : czs)
+        faceZone.primitiveFieldRef() =
+            SubField<scalar>(zoneID, mesh.nInternalFaces());
+        surfaceScalarField::Boundary& bfld = faceZone.boundaryFieldRef();
+        for (auto& pfld : bfld)
         {
-            UIndirectList<scalar>(faceZone, zone) = zone.index();
+            const fvPatch& fvp = pfld.patch();
+            pfld == SubField<scalar>(zoneID, fvp.size(), fvp.start());
         }
 
         //faceZone.correctBoundaryConditions();

applications/utilities/mesh/manipulation/mergeOrSplitBaffles/mergeOrSplitBaffles.C

@@ -286,7 +286,7 @@ int main(int argc, char *argv[])
 
     #include "addOverwriteOption.H"
     #include "addRegionOption.H"
-    argList::addBoolOption
+    argList::addOption
     (
         "dict",
         "file",

applications/utilities/postProcessing/dataConversion/foamToEnsight/convertLagrangian.H

@@ -67,10 +67,10 @@ if (doLagrangian)
             }
         }
 
-        forAllConstIters(theseCloudFields, fieldIter)
+        // Field order may differ on individual processors, so sort by name
+        for (const word& fieldName : theseCloudFields.sortedToc())
         {
-            const word& fieldName = fieldIter.key();
-            const word& fieldType = fieldIter.val();
+            const word& fieldType = theseCloudFields[fieldName];
 
             IOobject fieldObject
             (
@@ -88,9 +88,13 @@ if (doLagrangian)
                 // but that combination does not work.
                 // So check the header and sync globally
 
+                const bool parRun = Pstream::parRun();
+                Pstream::parRun() = false;
+
                 fieldExists =
                     fieldObject.typeHeaderOk<IOField<scalar>>(false);
 
+                Pstream::parRun() = parRun;
                 reduce(fieldExists, orOp<bool>());
             }
 

applications/utilities/postProcessing/lumped/lumpedPointForces/Make/options

@@ -6,5 +6,6 @@ EXE_INC = \
 
 EXE_LIBS = \
     -lfiniteVolume \
+    -lfileFormats \
     -lmeshTools \
     -llumpedPointMotion

applications/utilities/postProcessing/lumped/lumpedPointMovement/Make/options

@@ -6,5 +6,6 @@ EXE_INC = \
 
 EXE_LIBS = \
     -lfiniteVolume \
+    -lfileFormats \
     -lmeshTools \
     -llumpedPointMotion

applications/utilities/preProcessing/PDRsetFields/Allwmake

@@ -5,6 +5,6 @@ cd "${0%/*}" || exit                                # Run from this directory
 #------------------------------------------------------------------------------
 
 wmake $targetType pdrFields
-wmake PDRsetFields
+wmake $targetType PDRsetFields
 
 #------------------------------------------------------------------------------

applications/utilities/preProcessing/setAlphaField/alphaFieldFunctions/implicitFunctions/cylinder/cylinderImplicitFunction.C

@@ -69,13 +69,12 @@ Foam::implicitFunctions::cylinderImplicitFunction::cylinderImplicitFunction
     const dictionary& dict
 )
 :
-    cylinderImplicitFunction
-    (
-        dict.get<point>("origin"),
-        dict.get<scalar>("radius"),
-        dict.getOrDefault<scalar>("scale", 1),
-        dict.get<vector>("direction")
-    )
+    // __INTEL_COMPILER bug with inheriting constructors?? (issue #1821)
+    origin_(dict.get<point>("origin")),
+    radius_(dict.get<scalar>("radius")),
+    scale_(dict.getOrDefault<scalar>("scale", 1)),
+    direction_(dict.get<vector>("direction").normalise()),
+    project_(tensor::I - direction_*direction_) // outer product
 {}
 
 

applications/utilities/preProcessing/setAlphaField/alphaFieldFunctions/implicitFunctions/sin/sinImplicitFunction.C

@@ -72,15 +72,13 @@ Foam::implicitFunctions::sinImplicitFunction::sinImplicitFunction
     const dictionary& dict
 )
 :
-    sinImplicitFunction
-    (
-        dict.get<scalar>("period"),
-        dict.getOrDefault<scalar>("phase", 0),
-        dict.get<scalar>("amplitude"),
-        dict.get<vector>("direction"),
-        dict.get<vector>("up"),
-        dict.get<vector>("origin")
-    )
+    // __INTEL_COMPILER bug with inheriting constructors?? (issue #1821)
+    period_(dict.get<scalar>("period")),
+    phase_(dict.getOrDefault<scalar>("phase", 0)),
+    amplitude_(dict.get<scalar>("amplitude")),
+    up_(dict.get<vector>("up").normalise()),
+    direction_(dict.get<vector>("direction").normalise()),
+    origin_(dict.get<vector>("origin"))
 {}
 
 

applications/utilities/preProcessing/setExprFields/setExprFields.C

@@ -898,7 +898,7 @@ int main(int argc, char *argv[])
                 );
             }
         }
-        else if (exprDictPtr.valid())
+        else
         {
             FatalErrorInFunction
                 << "No command-line or dictionary??" << nl << endl

modules/Allwmake

@@ -7,9 +7,22 @@ targetType=libso
 # Default build into OpenFOAM project locations unless specified with
 # -prefix or FOAM_MODULE_PREFIX env varable
 
-: "${FOAM_MODULE_PREFIX:=${FOAM_LIBBIN%/*}}"
+# Long form to avoid dash 0.5.8 error (issue #1757)
+[ -n "$FOAM_MODULE_PREFIX" ] || FOAM_MODULE_PREFIX="${FOAM_LIBBIN%/*}"
 export FOAM_MODULE_PREFIX
 
+echo "========================================"
+echo "Compile OpenFOAM modules"
+echo "prefix = $FOAM_MODULE_PREFIX"
+echo
+echo "    ignoring possible compilation errors"
+echo "    make certain to check the output file"
+echo
+set +e
+export WM_CONTINUE_ON_ERROR=true
+
+#------------------------------------------------------------------------------
+
 for moduleName in $(./list-modules)
 do
     if [ -d "$moduleName" ]

src/ODE/ODESolvers/ODESolver/ODESolverNew.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -49,7 +49,7 @@ Foam::autoPtr<Foam::ODESolver> Foam::ODESolver::New
             "ODESolver",
             solverType,
             *dictionaryConstructorTablePtr_
-        ) << exit(FatalError);
+        ) << exit(FatalIOError);
     }
 
     return autoPtr<ODESolver>(cstrIter()(odes, dict));

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

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018 OpenCFD Ltd.
+    Copyright (C) 2018-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -86,7 +86,19 @@ inline bool Foam::HashPtrTable<T, Key, Hash>::set
     T* ptr
 )
 {
-    return this->parent_type::set(key, ptr);
+    // Newer: const T* old = this->get(key);
+
+    iterator iter(this->find(key));
+    const T* old = (iter.good() ? iter.val() : nullptr);
+
+    const bool ok = this->parent_type::set(key, ptr);
+
+    if (ok && old != ptr)
+    {
+        delete const_cast<T*>(old);
+    }
+
+    return ok;
 }
 
 

src/OpenFOAM/containers/PtrLists/PtrListOps/PtrListOpsTemplates.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -79,7 +79,7 @@ void Foam::sortedOrder
 template<class T>
 void Foam::sort(UPtrList<T>& list)
 {
-    labelList order(input.size());
+    labelList order(list.size());
     sortedOrder(list, order);
     list.sortOrder(order, false);  // false = allow nullptr
 }
@@ -88,7 +88,7 @@ void Foam::sort(UPtrList<T>& list)
 template<class T, class Compare>
 void Foam::sort(UPtrList<T>& list, const Compare& comp)
 {
-    labelList order(input.size());
+    labelList order(list.size());
     sortedOrder(list, order, comp);
     list.sortOrder(order, false);  // false = allow nullptr
 }

src/OpenFOAM/interpolations/interpolationTable/interpolationTable.H

@@ -43,14 +43,13 @@ Description
 
     Read csv format:
     \verbatim
-        readerType      csv;
-        file            "<constant>/p0vsTime.csv";
-        hasHeaderLine   true;   // skip first line
-        timeColumn      0;      // time is in column 0
-        valueColumns    (1);    // value starts in column 1
+        readerType          csv;
+        file                "<constant>/p0vsTime.csv";
+        hasHeaderLine       true;   // skip first line
+        refColumn           0;      // reference (eg, time) is in column 0
+        componentColumns    (1);    // component values starts in column 1
     \endverbatim
 
-
 Note
     - Accessing an empty list results in an error.
     - Accessing a list with a single element always returns the same value.

src/OpenFOAM/interpolations/interpolationTable/tableReaders/csv/csvTableReader.C

@@ -37,6 +37,7 @@ template<class Type>
 Foam::labelList Foam::csvTableReader<Type>::getComponentColumns
 (
     const word& name,
+    std::initializer_list<std::pair<const char*,int>> compat,
     const dictionary& dict
 )
 {
@@ -45,7 +46,7 @@ Foam::labelList Foam::csvTableReader<Type>::getComponentColumns
 
     labelList cols;
 
-    ITstream& is = dict.lookup(name);
+    ITstream& is = dict.lookupCompat(name, compat);
     is.format(IOstream::ASCII);
     is >> cols;
     dict.checkITstream(is, name);
@@ -113,8 +114,11 @@ Foam::csvTableReader<Type>::csvTableReader(const dictionary& dict)
 :
     tableReader<Type>(dict),
     headerLine_(dict.get<bool>("hasHeaderLine")),
-    refColumn_(dict.get<label>("timeColumn")),
-    componentColumns_(getComponentColumns("valueColumns", dict)),
+    refColumn_(dict.getCompat<label>("refColumn", {{"timeColumn", 1912}})),
+    componentColumns_
+    (
+        getComponentColumns("componentColumns", {{"valueColumns", 1912}}, dict)
+    ),
     separator_(dict.getOrDefault<string>("separator", ",")[0])
 {}
 
@@ -218,9 +222,9 @@ void Foam::csvTableReader<Type>::write(Ostream& os) const
     tableReader<Type>::write(os);
 
     os.writeEntry("hasHeaderLine", headerLine_);
-    os.writeEntry("timeColumn", refColumn_);
+    os.writeEntry("refColumn", refColumn_);
 
-    // Force writing labelList in ascii
+    // Force writing labelList in ASCII
     const enum IOstream::streamFormat fmt = os.format();
     os.format(IOstream::ASCII);
     os.writeEntry("componentColumns", componentColumns_);

src/OpenFOAM/interpolations/interpolationTable/tableReaders/csv/csvTableReader.H

@@ -76,6 +76,7 @@ class csvTableReader
         static labelList getComponentColumns
         (
             const word& name,
+            std::initializer_list<std::pair<const char*,int>> compat,
             const dictionary& dict
         );
 

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

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -86,7 +86,7 @@ Foam::mapPolyMesh::mapPolyMesh(const polyMesh& mesh)
 
     forAll(faceZoneFaceMap_, zonei)
     {
-        pointZoneMap_[zonei] = identity(mesh.faceZones()[zonei].size());
+        faceZoneFaceMap_[zonei] = identity(mesh.faceZones()[zonei].size());
     }
 
     forAll(cellZoneMap_, zonei)

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

@@ -264,7 +264,7 @@ void Foam::Function1Types::CSV<Type>::writeData(Ostream& os) const
     os.writeEntry("nHeaderLine", nHeaderLine_);
     os.writeEntry("refColumn", refColumn_);
 
-    // Force writing labelList in ascii
+    // Force writing labelList in ASCII
     const enum IOstream::streamFormat fmt = os.format();
     os.format(IOstream::ASCII);
     os.writeEntry("componentColumns", componentColumns_);

src/Pstream/Allwmake-mpi

@@ -1,6 +1,6 @@
 #!/bin/sh
 cd "${0%/*}" || exit                                # Run from this directory
-. ${WM_PROJECT_DIR:?}/wmake/scripts/AllwmakeParseArguments -no-recursion
+. ${WM_PROJECT_DIR:?}/wmake/scripts/AllwmakeParseArguments -no-recursion "$@"
 . ${WM_PROJECT_DIR:?}/wmake/scripts/wmakeFunctions  # Require wmake functions
 
 #------------------------------------------------------------------------------

src/TurbulenceModels/turbulenceModels/RAS/realizableKE/realizableKE.C

@@ -190,6 +190,15 @@ realizableKE<BasicTurbulenceModel>::realizableKE
             1.2
         )
     ),
+    implicitKEpsilonSource_
+    (
+        Switch::lookupOrAddToDict
+        (
+            "implicitKEpsilonSource",
+            this->coeffDict_,
+            false
+        )
+    ),
 
     k_
     (
@@ -237,6 +246,11 @@ bool realizableKE<BasicTurbulenceModel>::read()
         C2_.readIfPresent(this->coeffDict());
         sigmak_.readIfPresent(this->coeffDict());
         sigmaEps_.readIfPresent(this->coeffDict());
+        implicitKEpsilonSource_.readIfPresent
+        (
+            "implicitKEpsilonSource",
+            this->coeffDict()
+        );
 
         return true;
     }
@@ -288,6 +302,41 @@ void realizableKE<BasicTurbulenceModel>::correct()
         )
     );
 
+    tmp<fvScalarMatrix> tEpsilonSourceEq
+    (
+        new fvScalarMatrix
+        (
+            epsilon_,
+            dimVolume*rho.dimensions()*epsilon_.dimensions()/dimTime
+        )
+    );
+
+
+    tmp<volScalarField> Cmu(rCmu(tgradU(), S2, magS));
+
+    if (implicitKEpsilonSource_)
+    {
+        tmp<volScalarField> tEpsSp
+        (
+           - C1*alpha*rho*eta/sqrt(this->nut_/Cmu())
+           + alpha*rho*C2_/(sqrt(this->nut_/Cmu()) + sqrt(this->nu()))
+        );
+
+        tEpsilonSourceEq.ref() +=
+            0.5*tEpsSp()*pow(epsilon_, 1.5)
+          - fvm::Sp(1.5*tEpsSp()*sqrt(epsilon_), epsilon_);
+    }
+    else
+    {
+        tEpsilonSourceEq.ref() +=
+            C1*alpha*rho*magS*epsilon_
+          - fvm::Sp
+            (
+                C2_*alpha*rho*epsilon_/(k_ + sqrt(nuLimited*epsilon_)),
+                epsilon_
+            );
+    }
+
     // Dissipation equation
     tmp<fvScalarMatrix> epsEqn
     (
@@ -295,12 +344,7 @@ void realizableKE<BasicTurbulenceModel>::correct()
       + fvm::div(alphaRhoPhi, epsilon_)
       - fvm::laplacian(alpha*rho*DepsilonEff(), epsilon_)
      ==
-        C1*alpha*rho*magS*epsilon_
-      - fvm::Sp
-        (
-            C2_*alpha*rho*epsilon_/(k_ + sqrt(nuLimited*epsilon_)),
-            epsilon_
-        )
+        tEpsilonSourceEq()
       + epsilonSource()
       + fvOptions(alpha, rho, epsilon_)
     );
@@ -314,6 +358,26 @@ void realizableKE<BasicTurbulenceModel>::correct()
 
 
     // Turbulent kinetic energy equation
+    tmp<fvScalarMatrix> tkSourceEq
+    (
+        new fvScalarMatrix
+        (
+            k_,
+            dimVolume*rho.dimensions()*k_.dimensions()/dimTime
+        )
+    );
+
+    if (implicitKEpsilonSource_)
+    {
+        tkSourceEq.ref() +=
+            alpha*rho*G + alpha*alpha*rho*rho*Cmu()/this->nut_*sqr(k_)
+          - fvm::Sp(2.0*alpha*alpha*rho*rho*Cmu()*k_/this->nut_, k_);
+    }
+    else
+    {
+        tkSourceEq.ref() += alpha*rho*G - fvm::Sp(alpha*rho*epsilon_/k_, k_);
+    }
+
 
     tmp<fvScalarMatrix> kEqn
     (
@@ -321,9 +385,8 @@ void realizableKE<BasicTurbulenceModel>::correct()
       + fvm::div(alphaRhoPhi, k_)
       - fvm::laplacian(alpha*rho*DkEff(), k_)
      ==
-        alpha*rho*G
       - fvm::SuSp(2.0/3.0*alpha*rho*divU, k_)
-      - fvm::Sp(alpha*rho*epsilon_/k_, k_)
+      + tkSourceEq()
       + kSource()
       + fvOptions(alpha, rho, k_)
     );

src/TurbulenceModels/turbulenceModels/RAS/realizableKE/realizableKE.H

@@ -97,6 +97,9 @@ protected:
             dimensionedScalar sigmak_;
             dimensionedScalar sigmaEps_;
 
+            //- Use implicit k-Epilon source
+            Switch implicitKEpsilonSource_;
+
 
         // Fields
 

src/conversion/ccm/Allwmake

@@ -12,7 +12,7 @@ warning="==> skip optional libccm adapter"
 # Link with static libccmio only (fewer issues)
 if have_ccmio
 then
-    wmake libso || echo "$warning (build issues detected)"
+    wmake $targetType || echo "$warning (build issues detected)"
 else
     echo $warning
 fi

src/dynamicFaMesh/interfaceTrackingFvMesh/fvPatchFields/freeSurfacePressure/freeSurfacePressureFvPatchScalarField.H

@@ -34,7 +34,7 @@ Description
     calculated as:
 
         \f[
-            p = pa - \vec{g} & \vec{r}
+            p = pa - \vec{g} \cdot \vec{r}
         \f]
 
     where

src/finiteVolume/expressions/base/fvExprDriverTemplates.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2010-2018 Bernhard Gschaider <bgschaid@hfd-research.com>
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -261,7 +261,14 @@ Foam::tmp<GeomField> Foam::expressions::fvExprDriver::getOrReadFieldImpl
 
         tfield.reset
         (
-            GeomField::New(name, meshRef, dimensioned<Type>(Zero))
+            GeomField::New
+            (
+                name,
+                meshRef,
+                dimensioned<Type>(Zero),
+                // Patch is zeroGradient (volFields) or calculated (other)
+                defaultBoundaryType(GeomField::null())
+            )
         );
 
         GeomField& fld = tfield.ref();

src/finiteVolume/expressions/fields/base/patchExprFieldBase.C

@@ -31,6 +31,7 @@ License
 #include "fvMesh.H"
 #include "fvPatch.H"
 #include "pointMesh.H"
+#include "stringOps.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -97,7 +98,6 @@ Foam::expressions::patchExprFieldBase::patchExprFieldBase
         else
         {
             // No gradient expression - same as Zero
-
             if (debug_)
             {
                 Info<< "No gradientExpr" << nl;
@@ -106,7 +106,14 @@ Foam::expressions::patchExprFieldBase::patchExprFieldBase
 
         if (dict.readIfPresent("fractionExpr", expr))
         {
-            if (!expr.empty() && expr != "0")
+            stringOps::inplaceTrim(expr);
+
+            if (expr == "0" || expr == "1")
+            {
+                // Special cases, handled with more efficiency
+                fracExpr_ = expr;
+            }
+            else if (!expr.empty())
             {
                 if (isPointVal)
                 {
@@ -116,11 +123,9 @@ Foam::expressions::patchExprFieldBase::patchExprFieldBase
                 fracExpr_ = expressions::exprString(expr, dict);
             }
         }
-        else
-        {
-            // No fraction expression - same as 1 (one)
-            // Mixed BC may elect to simply ignore gradient expression
-        }
+
+        // No fraction expression? - defer treatment to inherited BC
+        // Mixed BC may elect to simply ignore gradient expression
     }
 }
 

src/finiteVolume/expressions/fields/fvPatchFields/exprFixedValueFvPatchField.C

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Original code Copyright (C) 2009-2018 Bernhard Gschaider
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2009-2018 Bernhard Gschaider
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -89,7 +89,7 @@ Foam::exprFixedValueFvPatchField<Type>::exprFixedValueFvPatchField
     setDebug();
     DebugInFunction << nl;
 
-    // Basic sanity
+    // Require valueExpr
     if (this->valueExpr_.empty())
     {
         FatalIOErrorInFunction(dict)
@@ -97,6 +97,7 @@ Foam::exprFixedValueFvPatchField<Type>::exprFixedValueFvPatchField
             << exit(FatalIOError);
     }
 
+
     driver_.readDict(dict);
 
     if (dict.found("value"))
@@ -110,11 +111,12 @@ Foam::exprFixedValueFvPatchField<Type>::exprFixedValueFvPatchField
     {
         (*this) == this->patchInternalField();
 
+        #ifdef FULLDEBUG
         WarningInFunction
             << "No value defined for "
             << this->internalField().name() << " on "
-            << this->patch().name() << " - setting to internalField value "
-            << nl;
+            << this->patch().name() << " - using patch internal field" << endl;
+        #endif
     }
 
     if (this->evalOnConstruct_)
@@ -161,31 +163,29 @@ Foam::exprFixedValueFvPatchField<Type>::exprFixedValueFvPatchField
 template<class Type>
 void Foam::exprFixedValueFvPatchField<Type>::updateCoeffs()
 {
-    if (debug)
-    {
-        InfoInFunction
-            << "Value: " << this->valueExpr_ << nl
-            << "Variables: ";
-        driver_.writeVariableStrings(Info) << endl;
-    }
-
     if (this->updated())
     {
         return;
     }
 
-    DebugInFunction
-        << "updating" << nl;
+    if (debug)
+    {
+        InfoInFunction
+            << "Value: " << this->valueExpr_ << nl
+            << "Variables: ";
+        driver_.writeVariableStrings(Info) << nl;
+        Info<< "... updating" << endl;
+    }
+
 
     // Expression evaluation
     {
+        bool evalValue = (!this->valueExpr_.empty() && this->valueExpr_ != "0");
+
+
         driver_.clearVariables();
 
-        if (this->valueExpr_.empty())
-        {
-            (*this) == Zero;
-        }
-        else
+        if (evalValue)
         {
             tmp<Field<Type>> tresult(driver_.evaluate<Type>(this->valueExpr_));
 
@@ -196,6 +196,10 @@ void Foam::exprFixedValueFvPatchField<Type>::updateCoeffs()
 
             (*this) == tresult;
         }
+        else
+        {
+            (*this) == Zero;
+        }
     }
 
     fixedValueFvPatchField<Type>::updateCoeffs();

src/finiteVolume/expressions/fields/fvPatchFields/exprMixedFvPatchField.C

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Original code Copyright (C) 2009-2018 Bernhard Gschaider
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2009-2018 Bernhard Gschaider
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -91,20 +91,51 @@ Foam::exprMixedFvPatchField<Type>::exprMixedFvPatchField
     setDebug();
     DebugInFunction << nl;
 
-    // Basic sanity checks
+    // Require one or both of valueExpr, gradientExpr
     if (this->valueExpr_.empty() && this->gradExpr_.empty())
     {
-        if (this->valueExpr_.empty())
+        FatalIOErrorInFunction(dict)
+            << "For " << this->internalField().name() << " on "
+            << this->patch().name() << nl
+            << "Require either or both: valueExpr and gradientExpr" << nl
+            << exit(FatalIOError);
+    }
+
+    if (this->fracExpr_.empty())
+    {
+        // No fractionExpr. Expect only one of valueExpr or gradientExpr
+        if (!this->valueExpr_.empty() && !this->gradExpr_.empty())
         {
-            FatalIOErrorInFunction(dict)
-                << "The valueExpr was not defined!" << nl
-                << exit(FatalIOError);
+            IOWarningInFunction(dict)
+                << "For " << this->internalField().name() << " on "
+                << this->patch().name() << nl
+                << "Recommend using fractionExpr when specifying both"
+                << " valueExpr and gradientExpr. Assuming a value of 1."
+                << nl << endl;
         }
+    }
+    else if (this->fracExpr_ == "0")
+    {
+        // Gradient only. Expect gradientExpr
         if (this->gradExpr_.empty())
         {
-            FatalIOErrorInFunction(dict)
-                << "The gradientExpr was not defined!" << nl
-                << exit(FatalIOError);
+            IOWarningInFunction(dict)
+                << "For " << this->internalField().name() << " on "
+                << this->patch().name() << nl
+                << "Gradient only, but did not specify gradientExpr."
+                << nl << endl;
+        }
+    }
+    else if (this->fracExpr_ == "1")
+    {
+        // Value only. Expect valueExpr
+        if (this->valueExpr_.empty())
+        {
+            IOWarningInFunction(dict)
+                << "For " << this->internalField().name() << " on "
+                << this->patch().name() << nl
+                << "Value only, but did not specify valueExpr."
+                << nl << endl;
         }
     }
 
@@ -114,14 +145,12 @@ Foam::exprMixedFvPatchField<Type>::exprMixedFvPatchField
     // Similar to fvPatchField constructor, which we have bypassed
     dict.readIfPresent("patchType", this->patchType());
 
+    bool needsRefValue = true;
     if (dict.found("refValue"))
     {
+        needsRefValue = false;
         this->refValue() = Field<Type>("refValue", dict, p.size());
     }
-    else
-    {
-        this->refValue() = this->patchInternalField();
-    }
 
     if (dict.found("value"))
     {
@@ -130,22 +159,27 @@ Foam::exprMixedFvPatchField<Type>::exprMixedFvPatchField
             Field<Type>("value", dict, p.size())
         );
 
-        if (!dict.found("refValue"))
+        if (needsRefValue)
         {
             // Ensure refValue has a sensible value for the "update" below
-            this->refValue() = Field<Type>("value", dict, p.size());
+            this->refValue() = static_cast<const Field<Type>&>(*this);
         }
     }
     else
     {
+        if (needsRefValue)
+        {
+            this->refValue() = this->patchInternalField();
+        }
+
         fvPatchField<Type>::operator=(this->refValue());
 
+        #ifdef FULLDEBUG
         WarningInFunction
             << "No value defined for "
-            << this->internalField().name()
-            << " on " << this->patch().name() << " therefore using "
-            << "the internal field next to the patch"
-            << endl;
+            << this->internalField().name() << " on "
+            << this->patch().name() << " - using patch internal field" << endl;
+        #endif
     }
 
 
@@ -164,7 +198,7 @@ Foam::exprMixedFvPatchField<Type>::exprMixedFvPatchField
     }
     else
     {
-        this->valueFraction() = 1;
+        this->valueFraction() = scalar(1);
     }
 
 
@@ -234,6 +268,11 @@ Foam::exprMixedFvPatchField<Type>::exprMixedFvPatchField
 template<class Type>
 void Foam::exprMixedFvPatchField<Type>::updateCoeffs()
 {
+    if (this->updated())
+    {
+        return;
+    }
+
     if (debug)
     {
         InfoInFunction
@@ -241,45 +280,67 @@ void Foam::exprMixedFvPatchField<Type>::updateCoeffs()
             << "Gradient: " << this->gradExpr_ << nl
             << "Fraction: " << this->fracExpr_ << nl
             << "Variables: ";
-
-        driver_.writeVariableStrings(Info) << endl;
+        driver_.writeVariableStrings(Info) << nl;
+        Info<< "... updating" << endl;
     }
 
-    if (this->updated())
-    {
-        return;
-    }
-
-    DebugInFunction << " - updating" << nl;
-
 
     // Expression evaluation
     {
-        driver_.clearVariables();
+        bool evalValue = (!this->valueExpr_.empty() && this->valueExpr_ != "0");
+        bool evalGrad = (!this->gradExpr_.empty() && this->gradExpr_ != "0");
+        bool evalFrac = (!this->fracExpr_.empty());
+        scalar fraction = 1;
+
+        // Have one or both of valueExpr, gradientExpr (checked in constructor)
 
         if (this->valueExpr_.empty())
         {
-            this->refValue() = Zero;
+            // No value expression -> gradient only
+            fraction = 0;
+            evalValue = false;
+            evalFrac = false;
         }
-        else
-        {
-            this->refValue() = driver_.evaluate<Type>(this->valueExpr_);
-        }
-
-        bool evalGrad = !this->gradExpr_.empty();
-
-        if (this->fracExpr_.empty() || this->fracExpr_ == "1")
+        else if (this->gradExpr_.empty())
         {
+            // No gradient expression -> value only
+            fraction = 1;
             evalGrad = false;
-            this->valueFraction() = scalar(1);
+            evalFrac = false;
+        }
+        else if (this->fracExpr_.empty())
+        {
+            // No fractionExpr, but has both valueExpr and gradientExpr
+            // -> treat as value only (warning in constructor)
+            fraction = 1;
+            evalGrad = false;
+            evalFrac = false;
         }
         else if (this->fracExpr_ == "0")
         {
-            this->valueFraction() = Zero;
+            // Gradient only
+            fraction = 0;
+            evalValue = false;
+            evalFrac = false;
+        }
+        else if (this->fracExpr_ == "1")
+        {
+            // Value only
+            fraction = 1;
+            evalGrad = false;
+            evalFrac = false;
+        }
+
+
+        driver_.clearVariables();
+
+        if (evalValue)
+        {
+            this->refValue() = driver_.evaluate<Type>(this->valueExpr_);
         }
         else
         {
-            this->valueFraction() = driver_.evaluate<scalar>(this->fracExpr_);
+            this->refValue() = Zero;
         }
 
         if (evalGrad)
@@ -290,6 +351,15 @@ void Foam::exprMixedFvPatchField<Type>::updateCoeffs()
         {
             this->refGrad() = Zero;
         }
+
+        if (evalFrac)
+        {
+            this->valueFraction() = driver_.evaluate<scalar>(this->fracExpr_);
+        }
+        else
+        {
+            this->valueFraction() = fraction;
+        }
     }
 
     mixedFvPatchField<Type>::updateCoeffs();

src/finiteVolume/fields/fvPatchFields/derived/pressureDirectedInletOutletVelocity/pressureDirectedInletOutletVelocityFvPatchVectorField.H

@@ -30,7 +30,7 @@ Group
     grpInletBoundaryConditions grpOutletBoundaryConditions
 
 Description
-    This velocity inlet/outlet boundary condition is applied to pressure
+    This velocity inlet/outlet boundary condition is applied to velocity
     boundaries where the pressure is specified.  A zero-gradient condition is
     applied for outflow (as defined by the flux); for inflow, the velocity
     is obtained from the flux with the specified inlet direction.

src/finiteVolume/fields/fvPatchFields/derived/pressureInletOutletVelocity/pressureInletOutletVelocityFvPatchVectorField.H

@@ -30,7 +30,7 @@ Group
     grpInletBoundaryConditions grpOutletBoundaryConditions
 
 Description
-    This velocity inlet/outlet boundary condition is applied to pressure
+    This velocity inlet/outlet boundary condition is applied to velocity
     boundaries where the pressure is specified.  A zero-gradient condition is
     applied for outflow (as defined by the flux); for inflow, the velocity is
     obtained from the patch-face normal component of the internal-cell value.

src/finiteVolume/finiteVolume/laplacianSchemes/relaxedNonOrthoGaussLaplacianScheme/relaxedNonOrthoGaussLaplacianScheme.C

@@ -0,0 +1,276 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "relaxedNonOrthoGaussLaplacianScheme.H"
+#include "surfaceInterpolate.H"
+#include "fvcDiv.H"
+#include "fvcGrad.H"
+#include "fvMatrices.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace fv
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type, class GType>
+tmp<fvMatrix<Type>>
+relaxedNonOrthoGaussLaplacianScheme<Type, GType>::fvmLaplacianUncorrected
+(
+    const surfaceScalarField& gammaMagSf,
+    const surfaceScalarField& deltaCoeffs,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    tmp<fvMatrix<Type>> tfvm
+    (
+        new fvMatrix<Type>
+        (
+            vf,
+            deltaCoeffs.dimensions()*gammaMagSf.dimensions()*vf.dimensions()
+        )
+    );
+    fvMatrix<Type>& fvm = tfvm.ref();
+
+    fvm.upper() = deltaCoeffs.primitiveField()*gammaMagSf.primitiveField();
+    fvm.negSumDiag();
+
+    forAll(vf.boundaryField(), patchi)
+    {
+        const fvPatchField<Type>& pvf = vf.boundaryField()[patchi];
+        const fvsPatchScalarField& pGamma = gammaMagSf.boundaryField()[patchi];
+        const fvsPatchScalarField& pDeltaCoeffs =
+            deltaCoeffs.boundaryField()[patchi];
+
+        if (pvf.coupled())
+        {
+            fvm.internalCoeffs()[patchi] =
+                pGamma*pvf.gradientInternalCoeffs(pDeltaCoeffs);
+            fvm.boundaryCoeffs()[patchi] =
+               -pGamma*pvf.gradientBoundaryCoeffs(pDeltaCoeffs);
+        }
+        else
+        {
+            fvm.internalCoeffs()[patchi] = pGamma*pvf.gradientInternalCoeffs();
+            fvm.boundaryCoeffs()[patchi] = -pGamma*pvf.gradientBoundaryCoeffs();
+        }
+    }
+
+    return tfvm;
+}
+
+
+template<class Type, class GType>
+tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
+relaxedNonOrthoGaussLaplacianScheme<Type, GType>::gammaSnGradCorr
+(
+    const surfaceVectorField& SfGammaCorr,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    const fvMesh& mesh = this->mesh();
+
+    tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> tgammaSnGradCorr
+    (
+        new GeometricField<Type, fvsPatchField, surfaceMesh>
+        (
+            IOobject
+            (
+                "gammaSnGradCorr("+vf.name()+')',
+                vf.instance(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            SfGammaCorr.dimensions()
+           *vf.dimensions()*mesh.deltaCoeffs().dimensions()
+        )
+    );
+
+    for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
+    {
+        tgammaSnGradCorr.ref().replace
+        (
+            cmpt,
+            fvc::dotInterpolate(SfGammaCorr, fvc::grad(vf.component(cmpt)))
+        );
+    }
+
+    return tgammaSnGradCorr;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type, class GType>
+tmp<GeometricField<Type, fvPatchField, volMesh>>
+relaxedNonOrthoGaussLaplacianScheme<Type, GType>::fvcLaplacian
+(
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    const fvMesh& mesh = this->mesh();
+
+    tmp<GeometricField<Type, fvPatchField, volMesh>> tLaplacian
+    (
+        fvc::div(this->tsnGradScheme_().snGrad(vf)*mesh.magSf())
+    );
+
+    tLaplacian.ref().rename("laplacian(" + vf.name() + ')');
+
+    return tLaplacian;
+}
+
+
+template<class Type, class GType>
+tmp<fvMatrix<Type>>
+relaxedNonOrthoGaussLaplacianScheme<Type, GType>::fvmLaplacian
+(
+    const GeometricField<GType, fvsPatchField, surfaceMesh>& gamma,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    const fvMesh& mesh = this->mesh();
+
+    typedef GeometricField<Type, fvsPatchField, surfaceMesh> SType;
+
+    const surfaceVectorField Sn(mesh.Sf()/mesh.magSf());
+
+    const surfaceVectorField SfGamma(mesh.Sf() & gamma);
+    const GeometricField<scalar, fvsPatchField, surfaceMesh> SfGammaSn
+    (
+        SfGamma & Sn
+    );
+    const surfaceVectorField SfGammaCorr(SfGamma - SfGammaSn*Sn);
+
+    tmp<fvMatrix<Type>> tfvm = fvmLaplacianUncorrected
+    (
+        SfGammaSn,
+        this->tsnGradScheme_().deltaCoeffs(vf),
+        vf
+    );
+    fvMatrix<Type>& fvm = tfvm.ref();
+
+    tmp<SType> tfaceFluxCorrection = gammaSnGradCorr(SfGammaCorr, vf);
+
+    if (this->tsnGradScheme_().corrected())
+    {
+        tfaceFluxCorrection.ref() +=
+            SfGammaSn*this->tsnGradScheme_().correction(vf);
+    }
+
+    const word corrName(tfaceFluxCorrection().name());
+
+    tmp<SType> trelaxedCorrection(new SType(tfaceFluxCorrection()));
+
+    const word oldName(corrName + "_0");
+    const scalar relax(vf.mesh().equationRelaxationFactor(oldName));
+
+    const objectRegistry& obr = vf.db();
+    if (obr.foundObject<SType>(oldName))
+    {
+        SType& oldCorrection = obr.lookupObjectRef<SType>(oldName);
+
+        trelaxedCorrection.ref() *= relax;
+        trelaxedCorrection.ref() += (1.0-relax)*oldCorrection;
+
+        oldCorrection = tfaceFluxCorrection;
+    }
+    else
+    {
+        SType* s = new SType(oldName, tfaceFluxCorrection);
+        s->store();
+    }
+
+    fvm.source() -=
+        mesh.V()
+       *fvc::div
+        (
+            trelaxedCorrection()
+        )().primitiveField();
+
+    if (mesh.fluxRequired(vf.name()))
+    {
+        fvm.faceFluxCorrectionPtr() = trelaxedCorrection.ptr();
+    }
+
+    return tfvm;
+}
+
+
+template<class Type, class GType>
+tmp<GeometricField<Type, fvPatchField, volMesh>>
+relaxedNonOrthoGaussLaplacianScheme<Type, GType>::fvcLaplacian
+(
+    const GeometricField<GType, fvsPatchField, surfaceMesh>& gamma,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    const fvMesh& mesh = this->mesh();
+
+    const surfaceVectorField Sn(mesh.Sf()/mesh.magSf());
+    const surfaceVectorField SfGamma(mesh.Sf() & gamma);
+    const GeometricField<scalar, fvsPatchField, surfaceMesh> SfGammaSn
+    (
+        SfGamma & Sn
+    );
+    const surfaceVectorField SfGammaCorr(SfGamma - SfGammaSn*Sn);
+
+    tmp<GeometricField<Type, fvPatchField, volMesh>> tLaplacian
+    (
+        fvc::div
+        (
+            SfGammaSn*this->tsnGradScheme_().snGrad(vf)
+          + gammaSnGradCorr(SfGammaCorr, vf)
+        )
+    );
+
+    tLaplacian.ref().rename
+    (
+        "laplacian(" + gamma.name() + ',' + vf.name() + ')'
+    );
+
+    return tLaplacian;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/finiteVolume/finiteVolume/laplacianSchemes/relaxedNonOrthoGaussLaplacianScheme/relaxedNonOrthoGaussLaplacianScheme.H

@@ -0,0 +1,190 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::fv::relaxedNonOrthoGaussLaplacianScheme
+
+Description
+    Basic second-order laplacian using face-gradients and Gauss' theorem.
+
+SourceFiles
+    relaxedNonOrthoGaussLaplacianScheme.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef relaxedNonOrthoGaussLaplacianScheme_H
+#define relaxedNonOrthoGaussLaplacianScheme_H
+
+#include "laplacianScheme.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace fv
+{
+
+/*---------------------------------------------------------------------------*\
+             Class relaxedNonOrthoGaussLaplacianScheme Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type, class GType>
+class relaxedNonOrthoGaussLaplacianScheme
+:
+    public fv::laplacianScheme<Type, GType>
+{
+    // Private Member Functions
+
+        tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> gammaSnGradCorr
+        (
+            const surfaceVectorField& SfGammaCorr,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        );
+
+        //- No copy construct
+        relaxedNonOrthoGaussLaplacianScheme
+        (
+            const relaxedNonOrthoGaussLaplacianScheme&
+        ) = delete;
+
+        //- No copy assignment
+        void operator=(const relaxedNonOrthoGaussLaplacianScheme&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("relaxedNonOrthoGauss");
+
+
+    // Constructors
+
+        //- Construct null
+        relaxedNonOrthoGaussLaplacianScheme(const fvMesh& mesh)
+        :
+            laplacianScheme<Type, GType>(mesh)
+        {}
+
+        //- Construct from Istream
+        relaxedNonOrthoGaussLaplacianScheme(const fvMesh& mesh, Istream& is)
+        :
+            laplacianScheme<Type, GType>(mesh, is)
+        {}
+
+        //- Construct from mesh, interpolation and snGradScheme schemes
+        relaxedNonOrthoGaussLaplacianScheme
+        (
+            const fvMesh& mesh,
+            const tmp<surfaceInterpolationScheme<GType>>& igs,
+            const tmp<snGradScheme<Type>>& sngs
+        )
+        :
+            laplacianScheme<Type, GType>(mesh, igs, sngs)
+        {}
+
+
+    //- Destructor
+    virtual ~relaxedNonOrthoGaussLaplacianScheme() = default;
+
+
+    // Member Functions
+
+        static tmp<fvMatrix<Type>> fvmLaplacianUncorrected
+        (
+            const surfaceScalarField& gammaMagSf,
+            const surfaceScalarField& deltaCoeffs,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        );
+
+        tmp<GeometricField<Type, fvPatchField, volMesh>> fvcLaplacian
+        (
+            const GeometricField<Type, fvPatchField, volMesh>&
+        );
+
+        tmp<fvMatrix<Type>> fvmLaplacian
+        (
+            const GeometricField<GType, fvsPatchField, surfaceMesh>&,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        );
+
+        tmp<GeometricField<Type, fvPatchField, volMesh>> fvcLaplacian
+        (
+            const GeometricField<GType, fvsPatchField, surfaceMesh>&,
+            const GeometricField<Type, fvPatchField, volMesh>&
+        );
+};
+
+
+// Use macros to emulate partial-specialisation of the Laplacian functions
+// for scalar diffusivity gamma
+
+#define defineFvmLaplacianScalarGamma(Type)                                    \
+                                                                               \
+template<>                                                                     \
+tmp<fvMatrix<Type>>                                                            \
+relaxedNonOrthoGaussLaplacianScheme<Type, scalar>::fvmLaplacian                \
+(                                                                              \
+    const GeometricField<scalar, fvsPatchField, surfaceMesh>&,                 \
+    const GeometricField<Type, fvPatchField, volMesh>&                         \
+);                                                                             \
+                                                                               \
+template<>                                                                     \
+tmp<GeometricField<Type, fvPatchField, volMesh>>                               \
+relaxedNonOrthoGaussLaplacianScheme<Type, scalar>::fvcLaplacian                \
+(                                                                              \
+    const GeometricField<scalar, fvsPatchField, surfaceMesh>&,                 \
+    const GeometricField<Type, fvPatchField, volMesh>&                         \
+);
+
+
+defineFvmLaplacianScalarGamma(scalar);
+defineFvmLaplacianScalarGamma(vector);
+defineFvmLaplacianScalarGamma(sphericalTensor);
+defineFvmLaplacianScalarGamma(symmTensor);
+defineFvmLaplacianScalarGamma(tensor);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "relaxedNonOrthoGaussLaplacianScheme.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/finiteVolume/laplacianSchemes/relaxedNonOrthoGaussLaplacianScheme/relaxedNonOrthoGaussLaplacianSchemes.C

@@ -0,0 +1,140 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "relaxedNonOrthoGaussLaplacianScheme.H"
+#include "fvMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makeFvLaplacianScheme(relaxedNonOrthoGaussLaplacianScheme)
+
+#define declareFvmLaplacianScalarGamma(Type)                                   \
+                                                                               \
+template<>                                                                     \
+Foam::tmp<Foam::fvMatrix<Foam::Type>>                                          \
+Foam::fv::relaxedNonOrthoGaussLaplacianScheme<Foam::Type, Foam::scalar>::      \
+fvmLaplacian                                                                   \
+(                                                                              \
+    const GeometricField<scalar, fvsPatchField, surfaceMesh>& gamma,           \
+    const GeometricField<Type, fvPatchField, volMesh>& vf                      \
+)                                                                              \
+{                                                                              \
+    const fvMesh& mesh = this->mesh();                                         \
+                                                                               \
+    typedef GeometricField<Type, fvsPatchField, surfaceMesh> SType;            \
+                                                                               \
+    GeometricField<scalar, fvsPatchField, surfaceMesh> gammaMagSf              \
+    (                                                                          \
+        gamma*mesh.magSf()                                                     \
+    );                                                                         \
+                                                                               \
+    tmp<fvMatrix<Type>> tfvm = fvmLaplacianUncorrected                         \
+    (                                                                          \
+        gammaMagSf,                                                            \
+        this->tsnGradScheme_().deltaCoeffs(vf),                                \
+        vf                                                                     \
+    );                                                                         \
+    fvMatrix<Type>& fvm = tfvm.ref();                                          \
+                                                                               \
+    if (this->tsnGradScheme_().corrected())                                    \
+    {                                                                          \
+        tmp<SType> tCorr(this->tsnGradScheme_().correction(vf));               \
+        const word corrName(tCorr().name());                                   \
+        tmp<SType> tfaceFluxCorrection(gammaMagSf*tCorr);                      \
+                                                                               \
+        tmp<SType> trelaxedCorrection(new SType(tfaceFluxCorrection()));       \
+                                                                               \
+        const word oldName(corrName + "_0");                                   \
+        const scalar relax(vf.mesh().equationRelaxationFactor(corrName));      \
+        const objectRegistry& obr = vf.db();                                   \
+        if (obr.foundObject<SType>(oldName))                                   \
+        {                                                                      \
+            SType& oldCorrection = obr.lookupObjectRef<SType>(oldName);        \
+            trelaxedCorrection.ref() *= relax;                                 \
+            trelaxedCorrection.ref() += (1.0-relax)*oldCorrection;             \
+                                                                               \
+            oldCorrection = trelaxedCorrection();                              \
+        }                                                                      \
+        else                                                                   \
+        {                                                                      \
+            SType* s = new SType(oldName, tfaceFluxCorrection);                \
+            s->store();                                                        \
+        }                                                                      \
+                                                                               \
+        tmp<Field<Type>> tcorr                                                 \
+        (                                                                      \
+            mesh.V()                                                           \
+           *fvc::div                                                           \
+            (                                                                  \
+                trelaxedCorrection()                                           \
+            )().primitiveField()                                               \
+        );                                                                     \
+                                                                               \
+        fvm.source() -= tcorr();                                               \
+                                                                               \
+        if (mesh.fluxRequired(vf.name()))                                      \
+        {                                                                      \
+            fvm.faceFluxCorrectionPtr() = trelaxedCorrection.ptr();            \
+        }                                                                      \
+    }                                                                          \
+                                                                               \
+    return tfvm;                                                               \
+}                                                                              \
+                                                                               \
+                                                                               \
+template<>                                                                     \
+Foam::tmp<Foam::GeometricField<Foam::Type, Foam::fvPatchField, Foam::volMesh>> \
+Foam::fv::relaxedNonOrthoGaussLaplacianScheme<Foam::Type, Foam::scalar>::fvcLaplacian         \
+(                                                                              \
+    const GeometricField<scalar, fvsPatchField, surfaceMesh>& gamma,           \
+    const GeometricField<Type, fvPatchField, volMesh>& vf                      \
+)                                                                              \
+{                                                                              \
+    const fvMesh& mesh = this->mesh();                                         \
+                                                                               \
+    tmp<GeometricField<Type, fvPatchField, volMesh>> tLaplacian                \
+    (                                                                          \
+        fvc::div(gamma*this->tsnGradScheme_().snGrad(vf)*mesh.magSf())         \
+    );                                                                         \
+                                                                               \
+    tLaplacian.ref().rename                                                    \
+    (                                                                          \
+        "laplacian(" + gamma.name() + ',' + vf.name() + ')'                    \
+    );                                                                         \
+                                                                               \
+    return tLaplacian;                                                         \
+}
+
+
+declareFvmLaplacianScalarGamma(scalar);
+declareFvmLaplacianScalarGamma(vector);
+declareFvmLaplacianScalarGamma(sphericalTensor);
+declareFvmLaplacianScalarGamma(symmTensor);
+declareFvmLaplacianScalarGamma(tensor);
+
+
+// ************************************************************************* //

src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolationScheme/surfaceInterpolationScheme.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -68,7 +68,7 @@ Foam::surfaceInterpolationScheme<Type>::New
             "discretisation",
             schemeName,
             *MeshConstructorTablePtr_
-        ) << exit(FatalError);
+        ) << exit(FatalIOError);
     }
 
     return cstrIter()(mesh, schemeData);

src/fvOptions/sources/derived/actuationDiskSource/actuationDiskSourceTemplates.C

@@ -127,7 +127,7 @@ void Foam::fv::actuationDiskSource::calcFroudeMethod
         Ostream& os = file();
         writeCurrentTime(os);
 
-        os  << Uref << tab << Cp << tab << Ct << tab << a << tab << T << tab
+        os  << Uref << tab << Cp << tab << Ct << tab << a << tab << T
             << endl;
     }
 }
@@ -238,7 +238,7 @@ void Foam::fv::actuationDiskSource::calcVariableScalingMethod
         Ostream& os = file();
         writeCurrentTime(os);
 
-        os  << Uref << tab << Cp << tab << Ct
+        os  << Uref << tab << Cp << tab << Ct << tab
             << Udisk << tab << CpStar << tab << CtStar << tab << T << tab << P
             << endl;
     }

src/lumpedPointMotion/state/lumpedPointState.C

@@ -152,7 +152,7 @@ Foam::lumpedPointState::lumpedPointState
         FatalErrorInFunction
             << "Have " << points_.size() << " points but "
             << angles_.size() << " angles" << nl
-            exit(FatalError);
+            << exit(FatalError);
         #else
         WarningInFunction
             << "Have " << points_.size() << " points but "

src/mesh/blockMesh/gradingDescriptor/gradingDescriptor.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2015 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -51,10 +51,7 @@ Foam::gradingDescriptor::gradingDescriptor
     nDivFraction_(nDivFraction),
     expansionRatio_(expansionRatio)
 {
-    if (expansionRatio_ < 0)
-    {
-        expansionRatio_ = 1.0/(-expansionRatio_);
-    }
+    correct();
 }
 
 
@@ -67,10 +64,7 @@ Foam::gradingDescriptor::gradingDescriptor
     nDivFraction_(1),
     expansionRatio_(expansionRatio)
 {
-    if (expansionRatio_ < 0)
-    {
-        expansionRatio_ = 1.0/(-expansionRatio_);
-    }
+    correct();
 }
 
 
@@ -82,6 +76,15 @@ Foam::gradingDescriptor::gradingDescriptor(Istream& is)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+void Foam::gradingDescriptor::correct()
+{
+    if (expansionRatio_ < 0)
+    {
+        expansionRatio_ = 1.0/(-expansionRatio_);
+    }
+}
+
+
 Foam::gradingDescriptor Foam::gradingDescriptor::inv() const
 {
     return gradingDescriptor
@@ -129,6 +132,8 @@ Foam::Istream& Foam::operator>>(Istream& is, gradingDescriptor& gd)
         is.readEnd("gradingDescriptor");
     }
 
+    gd.correct();
+
     is.check(FUNCTION_NAME);
     return is;
 }

src/mesh/blockMesh/gradingDescriptor/gradingDescriptor.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2015 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -57,7 +57,7 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declarations
+// Forward Declarations
 class Istream;
 class Ostream;
 class gradingDescriptor;
@@ -81,12 +81,13 @@ class gradingDescriptor
 
 public:
 
+    // Friendship with gradingDescriptors
     friend class gradingDescriptors;
 
 
     // Constructors
 
-        //- Default constructor
+        //- Default construct (1, 1, 1)
         gradingDescriptor();
 
         //- Construct from components
@@ -110,8 +111,6 @@ public:
 
     // Member Functions
 
-    // Access
-
         scalar blockFraction() const
         {
             return blockFraction_;
@@ -127,7 +126,9 @@ public:
             return expansionRatio_;
         }
 
-    // Member Functions
+        //- Adjust expansion ratio.
+        //  Trap negative value and treat as its inverse.
+        void correct();
 
         //- Return the inverse gradingDescriptor with 1/expansionRatio
         gradingDescriptor inv() const;
@@ -143,8 +144,6 @@ public:
 
         friend Istream& operator>>(Istream&, gradingDescriptor&);
         friend Ostream& operator<<(Ostream&, const gradingDescriptor&);
-
-        friend Istream& operator>>(Istream&, gradingDescriptors&);
 };
 
 

src/mesh/blockMesh/gradingDescriptor/gradingDescriptors.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2015 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -32,7 +32,16 @@ License
 
 Foam::gradingDescriptors::gradingDescriptors()
 :
-    gradingDescriptors(gradingDescriptor())
+    List<gradingDescriptor>(1, gradingDescriptor())
+{}
+
+
+Foam::gradingDescriptors::gradingDescriptors
+(
+    const label len
+)
+:
+    List<gradingDescriptor>(len, gradingDescriptor())
 {}
 
 
@@ -44,6 +53,35 @@ Foam::gradingDescriptors::gradingDescriptors(const gradingDescriptor& gd)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+void Foam::gradingDescriptors::correct()
+{
+    for (gradingDescriptor& gd : *this)
+    {
+        gd.correct();
+    }
+}
+
+
+void Foam::gradingDescriptors::normalise()
+{
+    scalar sumBlockFraction = 0;
+    scalar sumNDivFraction = 0;
+
+    for (const gradingDescriptor& gd : *this)
+    {
+        sumBlockFraction += gd.blockFraction_;
+        sumNDivFraction += gd.nDivFraction_;
+    }
+
+    for (gradingDescriptor& gd : *this)
+    {
+        gd.blockFraction_ /= sumBlockFraction;
+        gd.nDivFraction_  /= sumNDivFraction;
+        gd.correct();
+    }
+}
+
+
 Foam::gradingDescriptors Foam::gradingDescriptors::inv() const
 {
     gradingDescriptors ret(*this);
@@ -67,6 +105,7 @@ Foam::Istream& Foam::operator>>(Istream& is, gradingDescriptors& gds)
     if (t.isNumber())
     {
         gds = gradingDescriptors(gradingDescriptor(t.number()));
+        gds.correct();
     }
     else
     {
@@ -75,28 +114,10 @@ Foam::Istream& Foam::operator>>(Istream& is, gradingDescriptors& gds)
         // Read the list for gradingDescriptors
         is >> static_cast<List<gradingDescriptor>&>(gds);
 
-        // Check state of Istream
-        is.check(FUNCTION_NAME);
-
-        // Normalize the blockFractions and nDivFractions
-        // of the list of gradingDescriptors
-
-        scalar sumBlockFraction = 0;
-        scalar sumNDivFraction = 0;
-
-        forAll(gds, i)
-        {
-            sumBlockFraction += gds[i].blockFraction_;
-            sumNDivFraction += gds[i].nDivFraction_;
-        }
-
-        forAll(gds, i)
-        {
-            gds[i].blockFraction_ /= sumBlockFraction;
-            gds[i].nDivFraction_ /= sumNDivFraction;
-        }
+        gds.normalise();
     }
 
+    is.check(FUNCTION_NAME);
     return is;
 }
 

src/mesh/blockMesh/gradingDescriptor/gradingDescriptors.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2015 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -46,10 +47,9 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declarations
-class gradingDescriptors;
-
+// Forward Declarations
 class Istream;
+class gradingDescriptors;
 Istream& operator>>(Istream& is, gradingDescriptors& gd);
 
 /*---------------------------------------------------------------------------*\
@@ -64,15 +64,26 @@ public:
 
     // Constructors
 
-        //- Default constructor
+        //- Default construct with a single default gradingDescriptor
         gradingDescriptors();
 
-        //- Construct from a gradingDescriptor
-        gradingDescriptors(const gradingDescriptor& gd);
+        //- Construct with specified number of default gradingDescriptor
+        explicit gradingDescriptors(const label len);
+
+        //- Construct from a single gradingDescriptor
+        explicit gradingDescriptors(const gradingDescriptor& gd);
 
 
     // Member Functions
 
+        //- Adjust expansion ratios.
+        //  Trap negative value and treat as its inverse.
+        void correct();
+
+        //- Normalize blockFractions and nDivFractions for the list
+        //- of gradingDescriptors, and call correct()
+        void normalise();
+
         //- Return the inverse gradingDescriptors with 1/expansionRatio
         gradingDescriptors inv() const;
 

src/meshTools/sets/faceSources/regionToFace/regionToFace.C

@@ -78,7 +78,7 @@ void Foam::regionToFace::markZone
     const indirectPrimitivePatch& patch,
     const label proci,
     const label facei,
-    const label zoneI,
+    const label zonei,
     labelList& faceZone
 ) const
 {
@@ -100,7 +100,7 @@ void Foam::regionToFace::markZone
                 edgeTopoDistanceData<label>
                 (
                     0,          // distance
-                    zoneI
+                    zonei
                 )
             );
         }
@@ -124,9 +124,13 @@ void Foam::regionToFace::markZone
 
     forAll(allFaceInfo, facei)
     {
-        if (allFaceInfo[facei].data() == zoneI)
+        if
+        (
+            allFaceInfo[facei].valid(calc.data())
+         && allFaceInfo[facei].data() == zonei
+        )
         {
-            faceZone[facei] = zoneI;
+            faceZone[facei] = zonei;
         }
     }
 }

src/parallel/decompose/Allwmake-mpi

@@ -1,6 +1,6 @@
 #!/bin/sh
 cd "${0%/*}" || exit                                # Run from this directory
-. ${WM_PROJECT_DIR:?}/wmake/scripts/AllwmakeParseArguments -no-recursion
+. ${WM_PROJECT_DIR:?}/wmake/scripts/AllwmakeParseArguments -no-recursion "$@"
 . ${WM_PROJECT_DIR:?}/wmake/scripts/wmakeFunctions  # Require wmake functions
 . ${WM_PROJECT_DIR:?}/wmake/scripts/have_scotch
 

tutorials/heatTransfer/chtMultiRegionSimpleFoam/jouleHeatingSolid/Allrun

@@ -7,6 +7,9 @@ cd "${0%/*}" || exit                                # Run from this directory
 
 runApplication $(getApplication)
 
-./createGraphs
+if notTest "$@"
+then
+    ./createGraphs
+fi
 
 #------------------------------------------------------------------------------

tutorials/incompressible/pimpleFoam/RAS/TJunctionFan/0.orig/p

@@ -27,7 +27,6 @@ boundaryField
         outOfBounds     clamp;
         direction       in;
         readerType      openFoam;
-        hasHeaderLine   true;
         file            "<constant>/FluxVsdP.dat";
         //nonDimensional  true;
         //rpm             300;

tutorials/incompressible/simpleFoam/pipeCyclic/constant/turbulenceProperties

@@ -24,6 +24,8 @@ RAS
     turbulence      on;
 
     printCoeffs     on;
+
+    implicitKEpsilonSource  true;
 }
 
 

tutorials/incompressible/simpleFoam/pipeCyclic/system/fvSolution

@@ -22,14 +22,14 @@ solvers
         solver          GAMG;
         smoother        GaussSeidel;
 
-        tolerance       1e-06;
+        tolerance       1e-09;
         relTol          0.05;
     }
 
     pFinal
     {
         $p;
-        tolerance       1e-06;
+        tolerance       1e-09;
         relTol          0;
     }
 
@@ -37,14 +37,14 @@ solvers
     {
         solver          smoothSolver;
         smoother        symGaussSeidel;
-        tolerance       1e-05;
+        tolerance       1e-015;
         relTol          0.1;
     }
 
     "(U|k|epsilon)Final"
     {
         $U;
-        tolerance       1e-05;
+        tolerance       1e-15;
         relTol          0;
     }
 }
@@ -61,12 +61,14 @@ PIMPLE
 SIMPLE
 {
     nNonOrthogonalCorrectors 0;
+    /*
     residualControl
     {
         p               1e-2;
         U               1e-3;
         "(k|epsilon)"   1e-3;
     }
+    */
 }
 
 relaxationFactors
@@ -79,7 +81,7 @@ relaxationFactors
     {
         U               0.7;
         k               0.7;
-        "epsilon.*"     0.7;
+        epsilon         0.7;
     }
 }
 

tutorials/incompressible/simpleFoam/pitzDaily/constant/turbulenceProperties

@@ -26,6 +26,8 @@ RAS
     turbulence      on;
 
     printCoeffs     on;
+
+    implicitKEpsilonSource  true;
 }
 
 

wmake/scripts/have_scotch

@@ -90,6 +90,7 @@ search_scotch()
     then
         header=$(findFirstFile \
             "$prefix/include/$localDir/$incName" \
+            "$prefix/include/scotch/$incName" \
             "$prefix/include/$incName" \
         )
         library=$(findExtLib "$libName")
@@ -201,12 +202,15 @@ search_ptscotch()
         header=$(findFirstFile  \
             "$prefix/include/$FOAM_MPI/$incName" \
             "$prefix/include/$localDir/$incName" \
+            "$prefix/include/ptscotch/$incName" \
+            "$prefix/include/scotch/$incName" \
             "$prefix/include/$incName"
         )
         library="$(findExtLib $FOAM_MPI/$libName $libName)"
     elif isSystem "$prefix"
     then
         header=$(findFirstFile \
+            "/usr/local/include/$localDir/$incName" \
             "/usr/local/include/ptscotch/$incName" \
             "/usr/local/include/scotch/$incName" \
             "/usr/local/include/$incName" \

wmake/scripts/wmakeFunctions

@@ -173,13 +173,18 @@ findObjectDir()
     # Default (local) build directory
     if [ -z "$objectsDir" ]
     then
-        relativeDir="$absdir"
-        appDir=.
-        [ -d Make ] || appDir=$(findTarget .) || exit 1 # Fatal
-        absdir=$(expandPath "$appDir"/.)
+        if [ -d "$absdir/Make" ]
+        then
+            objectsDir="${absdir}/Make/${WM_OPTIONS}"
+        else
+            relativeDir="$absdir"
+            appDir=.
+            [ -d Make ] || appDir=$(findTarget .) || exit 1 # Fatal
+            absdir=$(expandPath "$appDir"/.)
 
-        relativeDir="${relativeDir#${absdir}}"
-        objectsDir="${appDir}/Make/${WM_OPTIONS}${relativeDir}"
+            relativeDir="${relativeDir#${absdir}}"
+            objectsDir="${appDir}/Make/${WM_OPTIONS}${relativeDir}"
+        fi
     fi
 
     echo "$objectsDir"

wmake/wclean

@@ -313,16 +313,18 @@ fi
 # Clean the 'Make' directory if present
 #------------------------------------------------------------------------------
 
-if [ -d "$MakeDir" ]
+if [ -d "$MakeDir" ] && [ -n "$WM_OPTIONS" ]
 then
-    objectsDir="$MakeDir/$WM_OPTIONS"
-    case "$PWD" in
-    ("$WM_PROJECT_DIR"/*)
-        buildPath="$WM_PROJECT_DIR/build/${WM_OPTIONS}"
-        objectsDir=$buildPath$(echo $PWD | sed s%$WM_PROJECT_DIR%% )
-        ;;
-    esac
-    rm -rf "$objectsDir" 2>/dev/null
+    # Remove in-source directory (if any)
+    rm -rf "$MakeDir/$WM_OPTIONS" 2>/dev/null
+
+    # Remove out-of-source directory (if applicable)
+    relativeDir="${PWD#${WM_PROJECT_DIR}/}"
+    if [ "$relativeDir" != "$PWD" ]
+    then
+        objectsDir="${WM_PROJECT_DIR}/build/${WM_OPTIONS}/${relativeDir}"
+        rm -rf "$objectsDir" 2>/dev/null
+    fi
 fi
 
 #------------------------------------------------------------------------------

wmake/wmake

@@ -177,6 +177,8 @@ die()
 # Default make is the "make" in the path
 make="make"
 
+# Print compiler/system information (serial only)
+printInfo() { make --no-print-directory -f "$WM_DIR"/makefiles/info "$@"; }
 
 #------------------------------------------------------------------------------
 
@@ -268,11 +270,11 @@ do
         -show-api | -show-ext-so | \
         -show-compile-c   | -show-c   | -show-cflags   | -show-cflags-arch | \
         -show-compile-cxx | -show-cxx | -show-cxxflags | -show-cxxflags-arch )
-            $make -f $WM_DIR/makefiles/info "${1#-show-}"
+            printInfo "${1#-show-}"
             optShow=true
             ;;
         -show-path-c | -show-path-cxx )
-            command -v $($make -f $WM_DIR/makefiles/info "${1#-show-path-}")
+            command -v $(printInfo "${1#-show-path-}")
             optShow=true
             ;;
         -a | -all | all)
@@ -323,7 +325,7 @@ do
             : "${optAll:=all}" # implies 'all', unless previously set
             ;;
         -version | --version)
-            $make -f $WM_DIR/makefiles/info api
+            printInfo api
             optShow=true
             break;
             ;;