Merge branch 'master' into cvm

ReleaseNotes-dev

@@ -184,6 +184,8 @@
     + =setSet=: allows time range (e.g. 0:100) in combination with -batch argument
     to execute the commands for multiple times.
 * Post-processing
+  + =paraFoam=, =foamToVTK=: full support for polyhedral cell type in recent
+     Paraview versions.
   + =foamToEnsight=: parallel continuous data. new =-nodeValues= option to generate and output nodal
     field data.
   + =singleCellMesh=: new utility to convert mesh and fields to a single cell

applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C

@@ -218,7 +218,7 @@ void Foam::kineticTheoryModel::solve(const volTensorField& gradUat)
     //  The solution is higly unstable close to the packing limit.
     gs0_ = radialModel_->g0
     (
-        min(max(alpha_, 1e-6), alphaMax_ - 0.01),
+        min(max(alpha_, scalar(1e-6)), alphaMax_ - 0.01),
         alphaMax_
     );
 
@@ -255,7 +255,7 @@ void Foam::kineticTheoryModel::solve(const volTensorField& gradUat)
     volScalarField J1 = 3.0*betaPrim;
     volScalarField J2 =
         0.25*sqr(betaPrim)*da_*sqr(Ur)
-       /(max(alpha_, 1e-6)*rhoa_*sqrtPi*(ThetaSqrt + TsmallSqrt));
+       /(max(alpha_, scalar(1e-6))*rhoa_*sqrtPi*(ThetaSqrt + TsmallSqrt));
 
     // bulk viscosity  p. 45 (Lun et al. 1984).
     lambda_ = (4.0/3.0)*sqr(alpha_)*rhoa_*da_*gs0_*(1.0+e_)*ThetaSqrt/sqrtPi;
@@ -309,7 +309,11 @@ void Foam::kineticTheoryModel::solve(const volTensorField& gradUat)
         volScalarField t1 = K1*alpha_ + rhoa_;
         volScalarField l1 = -t1*trD;
         volScalarField l2 = sqr(t1)*tr2D;
-        volScalarField l3 = 4.0*K4*max(alpha_, 1e-6)*(2.0*K3*trD2 + K2*tr2D);
+        volScalarField l3 =
+            4.0
+           *K4
+           *max(alpha_, scalar(1e-6))
+           *(2.0*K3*trD2 + K2*tr2D);
 
         Theta_ = sqr((l1 + sqrt(l2 + l3))/(2.0*(alpha_ + 1.0e-4)*K4));
     }

applications/utilities/parallelProcessing/redistributeMeshPar/redistributeMeshPar.C

@@ -561,10 +561,6 @@ int main(int argc, char *argv[])
         mkDir(args.path());
     }
 
-    // Switch timeStamp checking to one which does not do any
-    // parallel sync for same reason
-    regIOobject::fileModificationChecking = regIOobject::timeStamp;
-
 #   include "createTime.H"
 
 

etc/controlDict

@@ -595,6 +595,7 @@ DebugSwitches
     muSgsSpalartAllmarasWallFunction 0;
     multiDirRefinement  0;
     multiHoleInjector   0;
+    multiLevel          1;
     multivariateSelection 0;
     mutRoughWallFunction 0;
     mutSpalartAllmarasStandardRoughWallFunction  0;

src/OpenFOAM/algorithms/MeshWave/FaceCellWave.C

@@ -40,8 +40,7 @@ const Foam::scalar
 Foam::FaceCellWave<Type, TrackingData>::geomTol_ = 1e-6;
 
 template <class Type, class TrackingData>
-const Foam::scalar
+Foam::scalar Foam::FaceCellWave<Type, TrackingData>::propagationTol_ = 0.01;
-Foam::FaceCellWave<Type, TrackingData>::propagationTol_ = 0.01;
 
 template <class Type, class TrackingData>
 Foam::label Foam::FaceCellWave<Type, TrackingData>::dummyTrackData_ = 12345;
@@ -516,8 +515,8 @@ void Foam::FaceCellWave<Type, TrackingData>::handleProcPatches()
                 refCast<const processorPolyPatch>(patch);
 
             // Allocate buffers
-            labelList receiveFaces(patch.size());
+            labelList receiveFaces;
-            List<Type> receiveFacesInfo(patch.size());
+            List<Type> receiveFacesInfo;
 
             {
                 UIPstream fromNeighbour(procPatch.neighbProcNo(), pBufs);
@@ -675,8 +674,7 @@ Foam::FaceCellWave<Type, TrackingData>::FaceCellWave
     hasCyclicPatches_(hasCyclicPatch()),
     nEvals_(0),
     nUnvisitedCells_(mesh_.nCells()),
-    nUnvisitedFaces_(mesh_.nFaces()),
+    nUnvisitedFaces_(mesh_.nFaces())
-    iter_(0)
 {}
 
 
@@ -707,16 +705,15 @@ Foam::FaceCellWave<Type, TrackingData>::FaceCellWave
     hasCyclicPatches_(hasCyclicPatch()),
     nEvals_(0),
     nUnvisitedCells_(mesh_.nCells()),
-    nUnvisitedFaces_(mesh_.nFaces()),
+    nUnvisitedFaces_(mesh_.nFaces())
-    iter_(0)
 {
     // Copy initial changed faces data
     setFaceInfo(changedFaces, changedFacesInfo);
 
     // Iterate until nothing changes
-    iterate(maxIter);
+    label iter = iterate(maxIter);
 
-    if ((maxIter > 0) && (iter_ >= maxIter))
+    if ((maxIter > 0) && (iter >= maxIter))
     {
         FatalErrorIn
         (
@@ -796,7 +793,7 @@ Foam::label Foam::FaceCellWave<Type, TrackingData>::faceToCell()
             );
         }
 
-        // Neighbour. Hack for check if face has neighbour.
+        // Neighbour.
         if (faceI < nInternalFaces)
         {
             cellI = neighbour[faceI];
@@ -927,11 +924,13 @@ Foam::label Foam::FaceCellWave<Type, TrackingData>::iterate(const label maxIter)
         handleProcPatches();
     }
 
-    while (iter_ < maxIter)
+    label iter = 0;
+
+    while (iter < maxIter)
     {
         if (debug)
         {
-            Pout<< " Iteration " << iter_ << endl;
+            Pout<< " Iteration " << iter << endl;
         }
 
         nEvals_ = 0;
@@ -963,10 +962,11 @@ Foam::label Foam::FaceCellWave<Type, TrackingData>::iterate(const label maxIter)
             break;
         }
 
-        ++iter_;
+        ++iter;
     }
 
-    return nUnvisitedCells_;
+    return iter;
 }
 
+
 // ************************************************************************* //

src/OpenFOAM/algorithms/MeshWave/FaceCellWave.H

@@ -114,9 +114,6 @@ class FaceCellWave
         label nUnvisitedCells_;
         label nUnvisitedFaces_;
 
-        //- Iteration counter
-        label iter_;
-
 
     // Private Member Functions
 
@@ -232,7 +229,7 @@ class FaceCellWave
       // Private static data
 
             static const scalar geomTol_;
-            static const scalar propagationTol_;
+            static scalar propagationTol_;
 
             //- Used as default trackdata value to satisfy default template
             //  argument.
@@ -340,8 +337,8 @@ public:
             //  counted double)
             label cellToFace();
 
-            //- Iterate until no changes or maxIter reached. Returns number of
+            //- Iterate until no changes or maxIter reached.  Returns actual
-            //  unset cells (see getUnsetCells)
+            //  number of iterations.
             label iterate(const label maxIter);
 
 };

src/OpenFOAM/algorithms/MeshWave/MeshWave.H

@@ -134,8 +134,8 @@ public:
             return calc_.data();
         }
 
-        //- Iterate until no changes or maxIter reached. Returns number of
+        //- Iterate until no changes or maxIter reached. Returns actual
-        //  unset cells (see getUnsetCells)
+        //  number of iterations.
         label iterate(const label maxIter)
         {
             return calc_.iterate(maxIter);

src/OpenFOAM/db/IOobjects/IOdictionary/IOdictionary.C

@@ -66,6 +66,10 @@ void Foam::IOdictionary::readFile(const bool masterOnly)
           : Pstream::treeCommunication()
         );
 
+        // Master reads headerclassname from file. Make sure this gets
+        // transfered as well as contents.
+        Pstream::scatter(comms, const_cast<word&>(headerClassName()));
+        Pstream::scatter(comms, note());
 
         // Get my communication order
         const Pstream::commsStruct& myComm = comms[Pstream::myProcNo()];

src/OpenFOAM/db/regIOobject/regIOobjectRead.C

@@ -201,6 +201,11 @@ bool Foam::regIOobject::read()
           : Pstream::treeCommunication()
         );
 
+        // Master reads headerclassname from file. Make sure this gets
+        // transfered as well as contents.
+        Pstream::scatter(comms, const_cast<word&>(headerClassName()));
+        Pstream::scatter(comms, note());
+
 
         // Get my communication order
         const Pstream::commsStruct& myComm = comms[Pstream::myProcNo()];

src/OpenFOAM/interpolations/interpolationTable/interpolationTable.C

@@ -284,10 +284,8 @@ Type Foam::interpolationTable<Type>::rateOfChange(const scalar value) const
             case interpolationTable::REPEAT:
             {
                 // adjust lookupValue to >= minLimit
-                while (lookupValue < minLimit)
+                scalar span = maxLimit-minLimit;
-                {
+                lookupValue = fmod(lookupValue-minLimit, span) + minLimit;
-                    lookupValue += maxLimit;
-                }
                 break;
             }
         }
@@ -325,10 +323,8 @@ Type Foam::interpolationTable<Type>::rateOfChange(const scalar value) const
             case interpolationTable::REPEAT:
             {
                 // adjust lookupValue <= maxLimit
-                while (lookupValue > maxLimit)
+                scalar span = maxLimit-minLimit;
-                {
+                lookupValue = fmod(lookupValue-minLimit, span) + minLimit;
-                    lookupValue -= maxLimit;
-                }
                 break;
             }
         }
@@ -540,11 +536,9 @@ Type Foam::interpolationTable<Type>::operator()(const scalar value) const
             }
             case interpolationTable::REPEAT:
             {
-                // adjust lookupValue to >= minLimin
+                // adjust lookupValue to >= minLimit
-                while (lookupValue < minLimit)
+                scalar span = maxLimit-minLimit;
-                {
+                lookupValue = fmod(lookupValue-minLimit, span) + minLimit;
-                    lookupValue += maxLimit;
-                }
                 break;
             }
         }
@@ -582,10 +576,8 @@ Type Foam::interpolationTable<Type>::operator()(const scalar value) const
             case interpolationTable::REPEAT:
             {
                 // adjust lookupValue <= maxLimit
-                while (lookupValue > maxLimit)
+                scalar span = maxLimit-minLimit;
-                {
+                lookupValue = fmod(lookupValue-minLimit, span) + minLimit;
-                    lookupValue -= maxLimit;
-                }
                 break;
             }
         }

src/OpenFOAM/matrices/solution/solution.C

@@ -84,7 +84,11 @@ Foam::solution::solution
             dictName,
             obr.time().system(),
             obr,
-            IOobject::MUST_READ_IF_MODIFIED,
+            (
+                obr.readOpt() == IOobject::MUST_READ
+              ? IOobject::MUST_READ_IF_MODIFIED
+              : obr.readOpt()
+            ),
             IOobject::NO_WRITE
         )
     ),
@@ -94,7 +98,14 @@ Foam::solution::solution
     defaultRelaxationFactor_(0),
     solvers_(ITstream("solvers", tokenList())())
 {
-    read(solutionDict());
+    if
+    (
+        readOpt() == IOobject::MUST_READ
+     || readOpt() == IOobject::MUST_READ_IF_MODIFIED
+    )
+    {
+        read(solutionDict());
+    }
 }
 
 

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C

@@ -273,6 +273,30 @@ void Foam::cyclicPolyPatch::calcTransforms
             half1Normals,
             half0Tols
         );
+
+        if (transform_ == ROTATIONAL && !parallel() && forwardT().size() > 1)
+        {
+            const_cast<tensorField&>(forwardT()).setSize(1);
+            const_cast<tensorField&>(reverseT()).setSize(1);
+            const_cast<boolList&>(collocated()).setSize(1);
+
+            WarningIn
+            (
+                "cyclicPolyPatch::calcTransforms\n"
+                "    (\n"
+                "        const primitivePatch&,\n"
+                "        const UList<point>&,\n"
+                "        const UList<point>&,\n"
+                "        const UList<point>&,\n"
+                "        const UList<point>&\n"
+                "    )"
+            )   << "For patch " << name()
+                << " calculated non-uniform transform tensor even though"
+                << " the transform type is " << transformTypeNames[transform_]
+                << ". Setting the transformation tensor to be a uniform"
+                << " rotation."
+                << endl;
+        }
     }
 }
 
@@ -851,22 +875,6 @@ const Foam::edgeList& Foam::cyclicPolyPatch::coupledPoints() const
                 str<< "l " << vertI-1 << ' ' << vertI << nl;
             }
         }
-
-        // Remove any addressing calculated for the coupled edges calculation
-        const_cast<primitivePatch&>
-        (
-            static_cast<const primitivePatch&>
-            (
-                *this
-            )
-        ).clearOut();
-        const_cast<primitivePatch&>
-        (
-            static_cast<const primitivePatch&>
-            (
-                neighbPatch()
-            )
-        ).clearOut();
     }
     return *coupledPointsPtr_;
 }
@@ -1006,22 +1014,6 @@ const Foam::edgeList& Foam::cyclicPolyPatch::coupledEdges() const
                 str<< "l " << vertI-1 << ' ' << vertI << nl;
             }
         }
-
-        // Remove any addressing calculated for the coupled edges calculation
-        const_cast<primitivePatch&>
-        (
-            static_cast<const primitivePatch&>
-            (
-                *this
-            )
-        ).clearOut();
-        const_cast<primitivePatch&>
-        (
-            static_cast<const primitivePatch&>
-            (
-                neighbPatch
-            )
-        ).clearOut();
     }
     return *coupledEdgesPtr_;
 }

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.C

@@ -409,6 +409,7 @@ void Foam::processorPolyPatch::updateMesh(PstreamBuffers& pBufs)
         }
 
         // Remove any addressing used for shared points/edges calculation
+        // since mostly not needed.
         primitivePatch::clearOut();
     }
 }

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processorCyclic/processorCyclicPolyPatch.C

@@ -174,8 +174,6 @@ void Foam::processorCyclicPolyPatch::calcGeometry(PstreamBuffers& pBufs)
         // Update underlying cyclic
         coupledPolyPatch& pp = const_cast<coupledPolyPatch&>(referPatch());
 
-        Pout<< "updating geometry on refered patch:" << pp.name() << endl;
-
         pp.calcGeometry
         (
             *this,

src/OpenFOAM/primitives/ops/ops.H

@@ -76,6 +76,8 @@ EqOp(eqMag, x = mag(y))
 EqOp(plusEqMagSqr, x += magSqr(y))
 EqOp(maxEq, x = max(x, y))
 EqOp(minEq, x = min(x, y))
+EqOp(minMagSqrEq, x = (magSqr(x)<=magSqr(y) ? x : y))
+EqOp(maxMagSqrEq, x = (magSqr(x)>=magSqr(y) ? x : y))
 EqOp(andEq, x = (x && y))
 EqOp(orEq, x  = (x || y))
 

src/OpenFOAM/primitives/transform/transform.H

@@ -48,9 +48,11 @@ inline tensor rotationTensor
     const vector& n2
 )
 {
+    const scalar s = n1 & n2;
+    const vector n3 = n1 ^ n2;
     return
-        (n1 & n2)*I
+        s*I
-      + (1 - (n1 & n2))*sqr(n1 ^ n2)/(magSqr(n1 ^ n2) + VSMALL)
+      + (1 - s)*sqr(n3)/(magSqr(n3) + VSMALL)
       + (n2*n1 - n1*n2);
 }
 

src/finiteVolume/fields/fvPatchFields/derived/advective/advectiveFvPatchField.C

@@ -163,21 +163,19 @@ tmp<scalarField> advectiveFvPatchField<Type>::advectionSpeed() const
     const surfaceScalarField& phi =
         this->db().objectRegistry::lookupObject<surfaceScalarField>(phiName_);
 
-    fvsPatchField<scalar> phip = this->patch().lookupPatchField
+    fvsPatchField<scalar> phip =
-    (
+        this->patch().template lookupPatchField<surfaceScalarField, scalar>
-        phiName_,
+        (
-        reinterpret_cast<const surfaceScalarField*>(0),
+            phiName_
-        reinterpret_cast<const scalar*>(0)
+        );
-    );
 
     if (phi.dimensions() == dimDensity*dimVelocity*dimArea)
     {
-        const fvPatchScalarField& rhop = this->patch().lookupPatchField
+        const fvPatchScalarField& rhop =
-        (
+            this->patch().template lookupPatchField<volScalarField, scalar>
-            rhoName_,
+            (
-            reinterpret_cast<const volScalarField*>(0),
+                rhoName_
-            reinterpret_cast<const scalar*>(0)
+            );
-        );
 
         return phip/(rhop*this->patch().magSf());
     }

src/finiteVolume/fields/fvPatchFields/derived/inletOutlet/inletOutletFvPatchField.C

@@ -25,15 +25,10 @@ License
 
 #include "inletOutletFvPatchField.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class Type>
-inletOutletFvPatchField<Type>::inletOutletFvPatchField
+Foam::inletOutletFvPatchField<Type>::inletOutletFvPatchField
 (
     const fvPatch& p,
     const DimensionedField<Type, volMesh>& iF
@@ -49,7 +44,7 @@ inletOutletFvPatchField<Type>::inletOutletFvPatchField
 
 
 template<class Type>
-inletOutletFvPatchField<Type>::inletOutletFvPatchField
+Foam::inletOutletFvPatchField<Type>::inletOutletFvPatchField
 (
     const inletOutletFvPatchField<Type>& ptf,
     const fvPatch& p,
@@ -63,7 +58,7 @@ inletOutletFvPatchField<Type>::inletOutletFvPatchField
 
 
 template<class Type>
-inletOutletFvPatchField<Type>::inletOutletFvPatchField
+Foam::inletOutletFvPatchField<Type>::inletOutletFvPatchField
 (
     const fvPatch& p,
     const DimensionedField<Type, volMesh>& iF,
@@ -93,7 +88,7 @@ inletOutletFvPatchField<Type>::inletOutletFvPatchField
 
 
 template<class Type>
-inletOutletFvPatchField<Type>::inletOutletFvPatchField
+Foam::inletOutletFvPatchField<Type>::inletOutletFvPatchField
 (
     const inletOutletFvPatchField<Type>& ptf
 )
@@ -104,7 +99,7 @@ inletOutletFvPatchField<Type>::inletOutletFvPatchField
 
 
 template<class Type>
-inletOutletFvPatchField<Type>::inletOutletFvPatchField
+Foam::inletOutletFvPatchField<Type>::inletOutletFvPatchField
 (
     const inletOutletFvPatchField<Type>& ptf,
     const DimensionedField<Type, volMesh>& iF
@@ -118,19 +113,18 @@ inletOutletFvPatchField<Type>::inletOutletFvPatchField
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type>
-void inletOutletFvPatchField<Type>::updateCoeffs()
+void Foam::inletOutletFvPatchField<Type>::updateCoeffs()
 {
     if (this->updated())
     {
         return;
     }
 
-    const Field<scalar>& phip = this->patch().lookupPatchField
+    const Field<scalar>& phip =
-    (
+        this->patch().template lookupPatchField<surfaceScalarField, scalar>
-        phiName_,
+        (
-        reinterpret_cast<const surfaceScalarField*>(0),
+            phiName_
-        reinterpret_cast<const scalar*>(0)
+        );
-    );
 
     this->valueFraction() = 1.0 - pos(phip);
 
@@ -139,7 +133,7 @@ void inletOutletFvPatchField<Type>::updateCoeffs()
 
 
 template<class Type>
-void inletOutletFvPatchField<Type>::write(Ostream& os) const
+void Foam::inletOutletFvPatchField<Type>::write(Ostream& os) const
 {
     fvPatchField<Type>::write(os);
     if (phiName_ != "phi")
@@ -154,7 +148,7 @@ void inletOutletFvPatchField<Type>::write(Ostream& os) const
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template<class Type>
-void inletOutletFvPatchField<Type>::operator=
+void Foam::inletOutletFvPatchField<Type>::operator=
 (
     const fvPatchField<Type>& ptf
 )
@@ -167,8 +161,4 @@ void inletOutletFvPatchField<Type>::operator=
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //

src/finiteVolume/fields/fvPatchFields/derived/inletOutlet/inletOutletFvPatchField.H

@@ -56,6 +56,7 @@ protected:
 
     // Protected data
 
+        //- Name of flux field
         word phiName_;
 
 

src/finiteVolume/fields/fvPatchFields/derived/outletInlet/outletInletFvPatchField.C

@@ -25,21 +25,17 @@ License
 
 #include "outletInletFvPatchField.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class Type>
-outletInletFvPatchField<Type>::outletInletFvPatchField
+Foam::outletInletFvPatchField<Type>::outletInletFvPatchField
 (
     const fvPatch& p,
     const DimensionedField<Type, volMesh>& iF
 )
 :
-    mixedFvPatchField<Type>(p, iF)
+    mixedFvPatchField<Type>(p, iF),
+    phiName_("phi")
 {
     this->refValue() = *this;
     this->refGrad() = pTraits<Type>::zero;
@@ -48,7 +44,7 @@ outletInletFvPatchField<Type>::outletInletFvPatchField
 
 
 template<class Type>
-outletInletFvPatchField<Type>::outletInletFvPatchField
+Foam::outletInletFvPatchField<Type>::outletInletFvPatchField
 (
     const outletInletFvPatchField<Type>& ptf,
     const fvPatch& p,
@@ -56,19 +52,21 @@ outletInletFvPatchField<Type>::outletInletFvPatchField
     const fvPatchFieldMapper& mapper
 )
 :
-    mixedFvPatchField<Type>(ptf, p, iF, mapper)
+    mixedFvPatchField<Type>(ptf, p, iF, mapper),
+    phiName_(ptf.phiName_)
 {}
 
 
 template<class Type>
-outletInletFvPatchField<Type>::outletInletFvPatchField
+Foam::outletInletFvPatchField<Type>::outletInletFvPatchField
 (
     const fvPatch& p,
     const DimensionedField<Type, volMesh>& iF,
     const dictionary& dict
 )
 :
-    mixedFvPatchField<Type>(p, iF)
+    mixedFvPatchField<Type>(p, iF),
+    phiName_(dict.lookupOrDefault<word>("phi", "phi"))
 {
     this->refValue() = Field<Type>("outletValue", dict, p.size());
 
@@ -90,42 +88,43 @@ outletInletFvPatchField<Type>::outletInletFvPatchField
 
 
 template<class Type>
-outletInletFvPatchField<Type>::outletInletFvPatchField
+Foam::outletInletFvPatchField<Type>::outletInletFvPatchField
 (
     const outletInletFvPatchField<Type>& ptf
 )
 :
-    mixedFvPatchField<Type>(ptf)
+    mixedFvPatchField<Type>(ptf),
+    phiName_(ptf.phiName_)
 {}
 
 
 template<class Type>
-outletInletFvPatchField<Type>::outletInletFvPatchField
+Foam::outletInletFvPatchField<Type>::outletInletFvPatchField
 (
     const outletInletFvPatchField<Type>& ptf,
     const DimensionedField<Type, volMesh>& iF
 )
 :
-    mixedFvPatchField<Type>(ptf, iF)
+    mixedFvPatchField<Type>(ptf, iF),
+    phiName_(ptf.phiName_)
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type>
-void outletInletFvPatchField<Type>::updateCoeffs()
+void Foam::outletInletFvPatchField<Type>::updateCoeffs()
 {
     if (this->updated())
     {
         return;
     }
 
-    const fvsPatchField<scalar>& phip = this->patch().lookupPatchField
+    const fvsPatchField<scalar>& phip =
-    (
+        this->patch().template lookupPatchField<surfaceScalarField, scalar>
-        "phi",
+        (
-        reinterpret_cast<const surfaceScalarField*>(0),
+            phiName_
-        reinterpret_cast<const scalar*>(0)
+        );
-    );
 
     this->valueFraction() = pos(phip);
 
@@ -134,16 +133,16 @@ void outletInletFvPatchField<Type>::updateCoeffs()
 
 
 template<class Type>
-void outletInletFvPatchField<Type>::write(Ostream& os) const
+void Foam::outletInletFvPatchField<Type>::write(Ostream& os) const
 {
     fvPatchField<Type>::write(os);
+    if (phiName_ != "phi")
+    {
+        os.writeKeyword("phi") << phiName_ << token::END_STATEMENT << nl;
+    }
     this->refValue().writeEntry("outletValue", os);
     this->writeEntry("value", os);
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //

src/finiteVolume/fields/fvPatchFields/derived/outletInlet/outletInletFvPatchField.H

@@ -52,6 +52,14 @@ class outletInletFvPatchField
     public mixedFvPatchField<Type>
 {
 
+protected:
+
+    // Protected data
+
+        //- Name of flux field
+        word phiName_;
+
+
 public:
 
     //- Runtime type information

src/finiteVolume/fields/fvPatchFields/derived/selfContainedDirectMapped/selfContainedDirectMappedFixedValueFvPatchField.H

@@ -48,7 +48,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-              Class selfContainedDirectMappedFixedValueFvPatch Declaration
+        Class selfContainedDirectMappedFixedValueFvPatch Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class Type>
@@ -75,6 +75,7 @@ class selfContainedDirectMappedFixedValueFvPatchField
 
         mutable autoPtr<interpolation<Type> > interpolator_;
 
+
     // Private Member Functions
 
         //- Field to sample. Either on my or nbr mesh
@@ -83,6 +84,7 @@ class selfContainedDirectMappedFixedValueFvPatchField
         //- Access the interpolation method
         const interpolation<Type>& interpolator() const;
 
+
 public:
 
     //- Runtime type information

src/finiteVolume/fields/fvPatchFields/derived/waveTransmissive/waveTransmissiveFvPatchField.C

@@ -31,15 +31,10 @@ License
 #include "CrankNicholsonDdtScheme.H"
 #include "backwardDdtScheme.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class Type>
-waveTransmissiveFvPatchField<Type>::waveTransmissiveFvPatchField
+Foam::waveTransmissiveFvPatchField<Type>::waveTransmissiveFvPatchField
 (
     const fvPatch& p,
     const DimensionedField<Type, volMesh>& iF
@@ -52,7 +47,7 @@ waveTransmissiveFvPatchField<Type>::waveTransmissiveFvPatchField
 
 
 template<class Type>
-waveTransmissiveFvPatchField<Type>::waveTransmissiveFvPatchField
+Foam::waveTransmissiveFvPatchField<Type>::waveTransmissiveFvPatchField
 (
     const waveTransmissiveFvPatchField& ptf,
     const fvPatch& p,
@@ -67,7 +62,7 @@ waveTransmissiveFvPatchField<Type>::waveTransmissiveFvPatchField
 
 
 template<class Type>
-waveTransmissiveFvPatchField<Type>::waveTransmissiveFvPatchField
+Foam::waveTransmissiveFvPatchField<Type>::waveTransmissiveFvPatchField
 (
     const fvPatch& p,
     const DimensionedField<Type, volMesh>& iF,
@@ -81,7 +76,7 @@ waveTransmissiveFvPatchField<Type>::waveTransmissiveFvPatchField
 
 
 template<class Type>
-waveTransmissiveFvPatchField<Type>::waveTransmissiveFvPatchField
+Foam::waveTransmissiveFvPatchField<Type>::waveTransmissiveFvPatchField
 (
     const waveTransmissiveFvPatchField& ptpsf
 )
@@ -93,7 +88,7 @@ waveTransmissiveFvPatchField<Type>::waveTransmissiveFvPatchField
 
 
 template<class Type>
-waveTransmissiveFvPatchField<Type>::waveTransmissiveFvPatchField
+Foam::waveTransmissiveFvPatchField<Type>::waveTransmissiveFvPatchField
 (
     const waveTransmissiveFvPatchField& ptpsf,
     const DimensionedField<Type, volMesh>& iF
@@ -108,35 +103,32 @@ waveTransmissiveFvPatchField<Type>::waveTransmissiveFvPatchField
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type>
-tmp<scalarField> waveTransmissiveFvPatchField<Type>::advectionSpeed() const
+Foam::tmp<Foam::scalarField>
+Foam::waveTransmissiveFvPatchField<Type>::advectionSpeed() const
 {
     // Lookup the velocity and compressibility of the patch
-    const fvPatchField<scalar>& psip = this->patch().lookupPatchField
+    const fvPatchField<scalar>& psip =
-    (
+        this->patch().template lookupPatchField<volScalarField, scalar>
-        psiName_,
+        (
-        reinterpret_cast<const volScalarField*>(0),
+            psiName_
-        reinterpret_cast<const scalar*>(0)
+        );
-    );
 
     const surfaceScalarField& phi =
-        this->db().objectRegistry::lookupObject<surfaceScalarField>
+        this->db().template lookupObject<surfaceScalarField>(this->phiName_);
-        (this->phiName_);
 
-    fvsPatchField<scalar> phip = this->patch().lookupPatchField
+    fvsPatchField<scalar> phip =
-    (
+        this->patch().template lookupPatchField<surfaceScalarField, scalar>
-        this->phiName_,
+        (
-        reinterpret_cast<const surfaceScalarField*>(0),
+            this->phiName_
-        reinterpret_cast<const scalar*>(0)
+        );
-    );
 
     if (phi.dimensions() == dimDensity*dimVelocity*dimArea)
     {
-        const fvPatchScalarField& rhop = this->patch().lookupPatchField
+        const fvPatchScalarField& rhop =
-        (
+            this->patch().template lookupPatchField<volScalarField, scalar>
-            this->rhoName_,
+            (
-            reinterpret_cast<const volScalarField*>(0),
+                this->rhoName_
-            reinterpret_cast<const scalar*>(0)
+            );
-        );
 
         phip /= rhop;
     }
@@ -149,7 +141,7 @@ tmp<scalarField> waveTransmissiveFvPatchField<Type>::advectionSpeed() const
 
 
 template<class Type>
-void waveTransmissiveFvPatchField<Type>::write(Ostream& os) const
+void Foam::waveTransmissiveFvPatchField<Type>::write(Ostream& os) const
 {
     fvPatchField<Type>::write(os);
 
@@ -165,6 +157,7 @@ void waveTransmissiveFvPatchField<Type>::write(Ostream& os) const
     {
         os.writeKeyword("psi") << psiName_ << token::END_STATEMENT << nl;
     }
+
     os.writeKeyword("gamma") << gamma_ << token::END_STATEMENT << nl;
 
     if (this->lInf_ > SMALL)
@@ -179,8 +172,4 @@ void waveTransmissiveFvPatchField<Type>::write(Ostream& os) const
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //

src/finiteVolume/finiteVolume/fvSchemes/fvSchemes.C

@@ -256,7 +256,11 @@ Foam::fvSchemes::fvSchemes(const objectRegistry& obr)
             "fvSchemes",
             obr.time().system(),
             obr,
-            IOobject::MUST_READ_IF_MODIFIED,
+            (
+                obr.readOpt() == IOobject::MUST_READ
+              ? IOobject::MUST_READ_IF_MODIFIED
+              : obr.readOpt()
+            ),
             IOobject::NO_WRITE
         )
     ),
@@ -364,7 +368,14 @@ Foam::fvSchemes::fvSchemes(const objectRegistry& obr)
     // persistent settings across reads is incorrect
     clear();
 
-    read(schemesDict());
+    if
+    (
+        readOpt() == IOobject::MUST_READ
+     || readOpt() == IOobject::MUST_READ_IF_MODIFIED
+    )
+    {
+        read(schemesDict());
+    }
 }
 
 

src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/pointEdgeStructuredWalk.H

@@ -37,11 +37,7 @@ SourceFiles
 #define pointEdgeStructuredWalk_H
 
 #include "point.H"
-#include "label.H"
-#include "scalar.H"
 #include "tensor.H"
-#include "pTraits.H"
-
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -75,10 +71,12 @@ class pointEdgeStructuredWalk
     // Private Member Functions
 
         //- Evaluate distance to point.
+        template<class TrackingData>
         inline bool update
         (
             const pointEdgeStructuredWalk& w2,
-            const scalar tol
+            const scalar tol,
+            TrackingData& td
         );
 
 public:
@@ -111,80 +109,105 @@ public:
 
         // Needed by meshWave
 
-            //- Check whether still contains original (invalid) value.
+            //- Check whether origin has been changed at all or
-            inline bool valid() const;
+            //  still contains original (invalid) value.
+            template<class TrackingData>
+            inline bool valid(TrackingData& td) const;
 
             //- Check for identical geometrical data. Used for cyclics checking.
+            template<class TrackingData>
             inline bool sameGeometry
             (
                 const pointEdgeStructuredWalk&,
-                const scalar tol
+                const scalar tol,
+                TrackingData& td
             ) const;
 
             //- Convert origin to relative vector to leaving point
             //  (= point coordinate)
+            template<class TrackingData>
             inline void leaveDomain
             (
                 const polyPatch& patch,
                 const label patchPointI,
-                const point& pos
+                const point& pos,
+                TrackingData& td
             );
 
             //- Convert relative origin to absolute by adding entering point
+            template<class TrackingData>
             inline void enterDomain
             (
                 const polyPatch& patch,
                 const label patchPointI,
-                const point& pos
+                const point& pos,
+                TrackingData& td
             );
 
             //- Apply rotation matrix to origin
-            inline void transform(const tensor& rotTensor);
+            template<class TrackingData>
+            inline void transform
+            (
+                const tensor& rotTensor,
+                TrackingData& td
+            );
 
             //- Influence of edge on point
+            template<class TrackingData>
             inline bool updatePoint
             (
                 const polyMesh& mesh,
                 const label pointI,
                 const label edgeI,
                 const pointEdgeStructuredWalk& edgeInfo,
-                const scalar tol
+                const scalar tol,
+                TrackingData& td
             );
 
             //- Influence of different value on same point.
             //  Merge new and old info.
+            template<class TrackingData>
             inline bool updatePoint
             (
                 const polyMesh& mesh,
                 const label pointI,
                 const pointEdgeStructuredWalk& newPointInfo,
-                const scalar tol
+                const scalar tol,
+                TrackingData& td
             );
 
             //- Influence of different value on same point.
             //  No information about current position whatsoever.
+            template<class TrackingData>
             inline bool updatePoint
             (
                 const pointEdgeStructuredWalk& newPointInfo,
-                const scalar tol
+                const scalar tol,
+                TrackingData& td
             );
 
             //- Influence of point on edge.
+            template<class TrackingData>
             inline bool updateEdge
             (
                 const polyMesh& mesh,
                 const label edgeI,
                 const label pointI,
                 const pointEdgeStructuredWalk& pointInfo,
-                const scalar tol
+                const scalar tol,
+                TrackingData& td
             );
 
+            //- Same (like operator==)
+            template<class TrackingData>
+            inline bool equal(const pointEdgeStructuredWalk&, TrackingData&)
+            const;
+
 
     // Member Operators
 
         //Note: Used to determine whether to call update.
         inline bool operator==(const pointEdgeStructuredWalk&) const;
-
         inline bool operator!=(const pointEdgeStructuredWalk&) const;
 
 

src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/pointEdgeStructuredWalkI.H

@@ -29,13 +29,15 @@ License
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 // Update this with w2.
+template<class TrackingData>
 inline bool Foam::pointEdgeStructuredWalk::update
 (
     const pointEdgeStructuredWalk& w2,
-    const scalar tol
+    const scalar tol,
+    TrackingData& td
 )
 {
-    if (!valid())
+    if (!valid(td))
     {
         // current not yet set. Walked from w2 to here (=point0)
         dist_ = w2.dist_ + mag(point0_-w2.previousPoint_);
@@ -105,17 +107,20 @@ inline const Foam::vector& Foam::pointEdgeStructuredWalk::data() const
 }
 
 
-inline bool Foam::pointEdgeStructuredWalk::valid() const
+template<class TrackingData>
+inline bool Foam::pointEdgeStructuredWalk::valid(TrackingData& td) const
 {
     return previousPoint_ != vector::max;
 }
 
 
 // Checks for cyclic points
+template<class TrackingData>
 inline bool Foam::pointEdgeStructuredWalk::sameGeometry
 (
     const pointEdgeStructuredWalk& w2,
-    const scalar tol
+    const scalar tol,
+    TrackingData& td
 ) const
 {
     scalar diff = Foam::mag(dist() - w2.dist());
@@ -138,18 +143,25 @@ inline bool Foam::pointEdgeStructuredWalk::sameGeometry
 }
 
 
+template<class TrackingData>
 inline void Foam::pointEdgeStructuredWalk::leaveDomain
 (
     const polyPatch& patch,
     const label patchPointI,
-    const point& coord
+    const point& coord,
+    TrackingData& td
 )
 {
     previousPoint_ -= coord;
 }
 
 
-inline void Foam::pointEdgeStructuredWalk::transform(const tensor& rotTensor)
+template<class TrackingData>
+inline void Foam::pointEdgeStructuredWalk::transform
+(
+    const tensor& rotTensor,
+    TrackingData& td
+)
 {
     previousPoint_ = Foam::transform(rotTensor, previousPoint_);
 }
@@ -157,11 +169,13 @@ inline void Foam::pointEdgeStructuredWalk::transform(const tensor& rotTensor)
 
 // Update absolute geometric quantities. Note that distance (dist_)
 // is not affected by leaving/entering domain.
+template<class TrackingData>
 inline void Foam::pointEdgeStructuredWalk::enterDomain
 (
     const polyPatch& patch,
     const label patchPointI,
-    const point& coord
+    const point& coord,
+    TrackingData& td
 )
 {
     // back to absolute form
@@ -170,18 +184,20 @@ inline void Foam::pointEdgeStructuredWalk::enterDomain
 
 
 // Update this with information from connected edge
+template<class TrackingData>
 inline bool Foam::pointEdgeStructuredWalk::updatePoint
 (
     const polyMesh& mesh,
     const label pointI,
     const label edgeI,
     const pointEdgeStructuredWalk& edgeInfo,
-    const scalar tol
+    const scalar tol,
+    TrackingData& td
 )
 {
     if (inZone())
     {
-        return update(edgeInfo, tol);
+        return update(edgeInfo, tol, td);
     }
     else
     {
@@ -191,17 +207,19 @@ inline bool Foam::pointEdgeStructuredWalk::updatePoint
 
 
 // Update this with new information on same point
+template<class TrackingData>
 inline bool Foam::pointEdgeStructuredWalk::updatePoint
 (
     const polyMesh& mesh,
     const label pointI,
     const pointEdgeStructuredWalk& newPointInfo,
-    const scalar tol
+    const scalar tol,
+    TrackingData& td
 )
 {
     if (inZone())
     {
-        return update(newPointInfo, tol);
+        return update(newPointInfo, tol, td);
     }
     else
     {
@@ -211,29 +229,33 @@ inline bool Foam::pointEdgeStructuredWalk::updatePoint
 
 
 // Update this with new information on same point. No extra information.
+template<class TrackingData>
 inline bool Foam::pointEdgeStructuredWalk::updatePoint
 (
     const pointEdgeStructuredWalk& newPointInfo,
-    const scalar tol
+    const scalar tol,
+    TrackingData& td
 )
 {
-    return update(newPointInfo, tol);
+    return update(newPointInfo, tol, td);
 }
 
 
 // Update this with information from connected point
+template<class TrackingData>
 inline bool Foam::pointEdgeStructuredWalk::updateEdge
 (
     const polyMesh& mesh,
     const label edgeI,
     const label pointI,
     const pointEdgeStructuredWalk& pointInfo,
-    const scalar tol
+    const scalar tol,
+    TrackingData& td
 )
 {
     if (inZone())
     {
-        return update(pointInfo, tol);
+        return update(pointInfo, tol, td);
     }
     else
     {
@@ -242,6 +264,17 @@ inline bool Foam::pointEdgeStructuredWalk::updateEdge
 }
 
 
+template <class TrackingData>
+inline bool Foam::pointEdgeStructuredWalk::equal
+(
+    const pointEdgeStructuredWalk& rhs,
+    TrackingData& td
+) const
+{
+    return operator==(rhs);
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 inline bool Foam::pointEdgeStructuredWalk::operator==

src/fvMotionSolver/motionDiffusivity/inversePointDistance/inversePointDistanceDiffusivity.C

@@ -91,6 +91,9 @@ void Foam::inversePointDistanceDiffusivity::correct()
     List<pointEdgePoint> pointWallDist(mesh.nPoints());
     List<pointEdgePoint> edgeWallDist(mesh.nEdges());
 
+    int dummyTrackData = 0;
+
+
     {
         // Seeds
         List<pointEdgePoint> seedInfo(nPatchEdges);
@@ -108,7 +111,7 @@ void Foam::inversePointDistanceDiffusivity::correct()
             {
                 label pointI = meshPoints[i];
 
-                if (!pointWallDist[pointI].valid())
+                if (!pointWallDist[pointI].valid(dummyTrackData))
                 {
                     // Not yet seeded
                     seedInfo[nPatchEdges] = pointEdgePoint
@@ -135,7 +138,8 @@ void Foam::inversePointDistanceDiffusivity::correct()
 
             pointWallDist,
             edgeWallDist,
-            mesh.globalData().nTotalPoints() // max iterations
+            mesh.globalData().nTotalPoints(),// max iterations
+            dummyTrackData
         );
     }
 

src/lagrangian/basic/InteractionLists/InteractionLists.C

@@ -86,7 +86,8 @@ void Foam::InteractionLists<ParticleType>::buildInteractionLists()
     // referred cells that they interact with.
     PackedBoolList cellInRangeOfCoupledPatch(mesh_.nCells(), false);
 
-    // IAndT: index and transform
+    // IAndT: index (=local cell index) and transform (from
+    // globalIndexAndTransform)
     DynamicList<labelPair> cellIAndTToExchange;
 
     DynamicList<treeBoundBox> cellBbsToExchange;

src/lagrangian/basic/Particle/Particle.C

@@ -371,10 +371,21 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
         {
             // Apply tracking correction towards tet centre
 
+            if (debug)
+            {
+                Pout<< "tracking rescue using tetCentre from " << position();
+            }
+
             position_ +=
                 Cloud<ParticleType>::trackingCorrectionTol
                *(tet.centre() - position_);
 
+            if (debug)
+            {
+                Pout<< " to " << position() << " due to "
+                    << (tet.centre() - position_) << endl;
+            }
+
             cloud_.trackingRescue();
 
             return trackFraction;
@@ -639,6 +650,12 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
 
         tetPointRef tet = currentTet();
 
+        if (debug)
+        {
+            Pout<< "tracking rescue for lambdaMin:" << lambdaMin
+                << "from " << position();
+        }
+
         position_ +=
             Cloud<ParticleType>::trackingCorrectionTol
            *(tet.centre() - position_);
@@ -692,6 +709,11 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
             }
         }
 
+        if (debug)
+        {
+            Pout<< " to " << position() << endl;
+        }
+
         cloud_.trackingRescue();
     }
 

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C

@@ -159,11 +159,6 @@ Foam::label Foam::autoLayerDriver::mergePatchFacesUndo
             mesh.clearOut();
         }
 
-        if (meshRefiner_.overwrite())
-        {
-            mesh.setInstance(meshRefiner_.oldInstance());
-        }
-
         faceCombiner.updateMesh(map);
 
         meshRefiner_.updateMesh(map, labelList(0));
@@ -324,11 +319,6 @@ Foam::label Foam::autoLayerDriver::mergePatchFacesUndo
                 mesh.clearOut();
             }
 
-            if (meshRefiner_.overwrite())
-            {
-                mesh.setInstance(meshRefiner_.oldInstance());
-            }
-
             faceCombiner.updateMesh(map);
 
             // Renumber restore maps
@@ -396,11 +386,6 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::autoLayerDriver::doRemovePoints
         mesh.clearOut();
     }
 
-    if (meshRefiner_.overwrite())
-    {
-        mesh.setInstance(meshRefiner_.oldInstance());
-    }
-
     pointRemover.updateMesh(map);
     meshRefiner_.updateMesh(map, labelList(0));
 
@@ -454,11 +439,6 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::autoLayerDriver::doRestorePoints
         mesh.clearOut();
     }
 
-    if (meshRefiner_.overwrite())
-    {
-        mesh.setInstance(meshRefiner_.oldInstance());
-    }
-
     pointRemover.updateMesh(map);
     meshRefiner_.updateMesh(map, labelList(0));
 
@@ -3100,10 +3080,7 @@ void Foam::autoLayerDriver::addLayers
         //?neccesary? Update fields
         newMesh.updateMesh(map);
 
-        if (meshRefiner_.overwrite())
+        newMesh.setInstance(meshRefiner_.timeName());
-        {
-            newMesh.setInstance(meshRefiner_.oldInstance());
-        }
 
         // Update numbering on addLayer:
         // - cell/point labels to be newMesh.
@@ -3229,11 +3206,6 @@ void Foam::autoLayerDriver::addLayers
         mesh.clearOut();
     }
 
-    if (meshRefiner_.overwrite())
-    {
-        mesh.setInstance(meshRefiner_.oldInstance());
-    }
-
     meshRefiner_.updateMesh(map, labelList(0));
 
 

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C

@@ -756,6 +756,9 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo
     // Distance to wall
     List<pointData> pointWallDist(mesh.nPoints());
 
+    // Dummy additional info for PointEdgeWave
+    int dummyTrackData = 0;
+
 
     // 1. Calculate distance to points where displacement is specified.
     {
@@ -783,7 +786,8 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo
             wallInfo,
             pointWallDist,
             edgeWallDist,
-            mesh.globalData().nTotalPoints()    // max iterations
+            mesh.globalData().nTotalPoints(),   // max iterations
+            dummyTrackData
         );
     }
 
@@ -813,7 +817,7 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo
                 {
                     label pointI = e[ep];
 
-                    if (!pointMedialDist[pointI].valid())
+                    if (!pointMedialDist[pointI].valid(dummyTrackData))
                     {
                         maxPoints.append(pointI);
                         maxInfo.append
@@ -857,7 +861,7 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo
                 {
                     label pointI = meshPoints[i];
 
-                    if (!pointMedialDist[pointI].valid())
+                    if (!pointMedialDist[pointI].valid(dummyTrackData))
                     {
                         maxPoints.append(pointI);
                         maxInfo.append
@@ -888,7 +892,8 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo
 
             pointMedialDist,
             edgeMedialDist,
-            mesh.globalData().nTotalPoints()    // max iterations
+            mesh.globalData().nTotalPoints(),   // max iterations
+            dummyTrackData
         );
 
         // Extract medial axis distance as pointScalarField

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/pointData/pointData.C

@@ -32,19 +32,20 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const pointData& wDist)
     if (os.format() == IOstream::ASCII)
     {
         return os
-            << wDist.origin() << token::SPACE << wDist.distSqr()
+            << static_cast<const pointEdgePoint&>(wDist)
             << token::SPACE << wDist.s() << token::SPACE << wDist.v();
     }
     else
     {
         return os
-            << wDist.origin() << wDist.distSqr() << wDist.s() << wDist.v();
+            << static_cast<const pointEdgePoint&>(wDist)
+            << wDist.s() << wDist.v();
     }
 }
 
 Foam::Istream& Foam::operator>>(Istream& is, pointData& wDist)
 {
-    return is >> wDist.origin_ >> wDist.distSqr_ >> wDist.s_ >> wDist.v_;
+    return is >> static_cast<pointEdgePoint&>(wDist) >> wDist.s_ >> wDist.v_;
 }
 
 // ************************************************************************* //

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/pointData/pointData.H

@@ -25,10 +25,10 @@ Class
     Foam::pointData
 
 Description
-    Holds information regarding nearest wall point. Used in pointEdgeWave.
+    Variant of pointEdgePoint with some transported additional data. 
-    (so not standard meshWave)
+    WIP - should be templated on data like wallDistData.
-
+    Passive vector v_ is not a coordinate (so no enterDomain/leaveDomain
-    To be used in wall distance calculation.
+    transformation needed)
 
 SourceFiles
     pointDataI.H
@@ -39,9 +39,10 @@ SourceFiles
 #ifndef pointData_H
 #define pointData_H
 
-#include "point.H"
+#include "pointEdgePoint.H"
-#include "label.H"
+//#include "point.H"
-#include "tensor.H"
+//#include "label.H"
+//#include "tensor.H"
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -49,51 +50,22 @@ SourceFiles
 namespace Foam
 {
 
-// Class forward declarations
-class polyPatch;
-class polyMesh;
-
 /*---------------------------------------------------------------------------*\
                            Class pointData Declaration
 \*---------------------------------------------------------------------------*/
 
 class pointData
+:
+    public pointEdgePoint
 {
     // Private data
 
-        //- position of nearest wall center
-        point origin_;
-
-        //- normal distance (squared) from point to origin
-        scalar distSqr_;
-
         //- additional information.
         scalar s_;
 
         //- additional information.
         vector v_;
 
-
-    // Private Member Functions
-
-        //- Evaluate distance to point. Update distSqr, origin from whomever
-        //  is nearer pt. Return true if w2 is closer to point,
-        //  false otherwise.
-        inline bool update
-        (
-            const point&,
-            const pointData& w2,
-            const scalar tol
-        );
-
-        //- Combine current with w2. Update distSqr, origin if w2 has smaller
-        //  quantities and returns true.
-        inline bool update
-        (
-            const pointData& w2,
-            const scalar tol
-        );
-
 public:
 
     // Constructors
@@ -118,10 +90,6 @@ public:
 
         // Access
 
-            inline const point& origin() const;
-
-            inline scalar distSqr() const;
-
             inline scalar s() const;
 
             inline const vector& v() const;
@@ -129,81 +97,60 @@ public:
 
         // Needed by meshWave
 
-            //- Check whether origin has been changed at all or
-            //  still contains original (invalid) value.
-            inline bool valid() const;
-
-            //- Check for identical geometrical data. Used for cyclics checking.
-            inline bool sameGeometry(const pointData&, const scalar tol)
-             const;
-
-            //- Convert origin to relative vector to leaving point
-            //  (= point coordinate)
-            inline void leaveDomain
-            (
-                const polyPatch& patch,
-                const label patchPointI,
-                const point& pos
-            );
-
-            //- Convert relative origin to absolute by adding entering point
-            inline void enterDomain
-            (
-                const polyPatch& patch,
-                const label patchPointI,
-                const point& pos
-            );
-
             //- Apply rotation matrix to origin
-            inline void transform(const tensor& rotTensor);
+            template<class TrackingData>
+            inline void transform
+            (
+                const tensor& rotTensor,
+                TrackingData& td
+            );
 
             //- Influence of edge on point
+            template<class TrackingData>
             inline bool updatePoint
             (
                 const polyMesh& mesh,
                 const label pointI,
                 const label edgeI,
                 const pointData& edgeInfo,
-                const scalar tol
+                const scalar tol,
+                TrackingData& td
             );
 
             //- Influence of different value on same point.
             //  Merge new and old info.
+            template<class TrackingData>
             inline bool updatePoint
             (
                 const polyMesh& mesh,
                 const label pointI,
                 const pointData& newPointInfo,
-                const scalar tol
+                const scalar tol,
+                TrackingData& td
             );
 
             //- Influence of different value on same point.
             //  No information about current position whatsoever.
+            template<class TrackingData>
             inline bool updatePoint
             (
                 const pointData& newPointInfo,
-                const scalar tol
+                const scalar tol,
+                TrackingData& td
             );
 
             //- Influence of point on edge.
+            template<class TrackingData>
             inline bool updateEdge
             (
                 const polyMesh& mesh,
                 const label edgeI,
                 const label pointI,
                 const pointData& pointInfo,
-                const scalar tol
+                const scalar tol,
+                TrackingData& td
             );
 
-
-    // Member Operators
-
-        //Note: Used to determine whether to call update.
-        inline bool operator==(const pointData&) const;
-
-        inline bool operator!=(const pointData&) const;
-
-
     // IOstream Operators
 
         friend Ostream& operator<<(Ostream&, const pointData&);

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/pointData/pointDataI.H

@@ -26,118 +26,12 @@ License
 #include "polyMesh.H"
 #include "transform.H"
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-// Update this with w2 if w2 nearer to pt.
-inline bool Foam::pointData::update
-(
-    const point& pt,
-    const pointData& w2,
-    const scalar tol
-)
-{
-    scalar dist2 = magSqr(pt - w2.origin());
-
-    if (!valid())
-    {
-        distSqr_ = dist2;
-        origin_ = w2.origin();
-        s_ = w2.s();
-        v_ = w2.v();
-
-        return true;
-    }
-
-
-//    if (v_ != w2.v())
-//    {
-//        return false;
-//    }
-
-
-    scalar diff = distSqr_ - dist2;
-
-    if (diff < 0)
-    {
-        // already nearer to pt
-        return false;
-    }
-
-    if ((diff < SMALL) || ((distSqr_ > SMALL) && (diff/distSqr_ < tol)))
-    {
-        // don't propagate small changes
-        return false;
-    }
-    else
-    {
-        // update with new values
-        distSqr_ = dist2;
-        origin_ = w2.origin();
-        s_ = w2.s();
-        v_ = w2.v();
-
-        return true;
-    }
-}
-
-
-// Update this with w2 (information on same point)
-inline bool Foam::pointData::update
-(
-    const pointData& w2,
-    const scalar tol
-)
-{
-    if (!valid())
-    {
-        // current not yet set so use any value
-        distSqr_ = w2.distSqr();
-        origin_ = w2.origin();
-        s_ = w2.s();
-        v_ = w2.v();
-        return true;
-    }
-
-
-//    if (v_ != w2.v())
-//    {
-//        return false;
-//    }
-
-
-    scalar diff = distSqr_ - w2.distSqr();
-
-    if (diff < 0)
-    {
-        // already nearer to pt
-        return false;
-    }
-
-    if ((diff < SMALL) || ((distSqr_ > SMALL) && (diff/distSqr_ < tol)))
-    {
-        // don't propagate small changes
-        return false;
-    }
-    else
-    {
-        // update with new values
-        distSqr_ =  w2.distSqr();
-        origin_ = w2.origin();
-        s_ = w2.s();
-        v_ = w2.v();
-
-        return true;
-    }
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Null constructor
 inline Foam::pointData::pointData()
 :
-    origin_(point::max),
+    pointEdgePoint(),
-    distSqr_(GREAT),
     s_(GREAT),
     v_(point::max)
 {}
@@ -152,8 +46,7 @@ inline Foam::pointData::pointData
     const vector& v
 )
 :
-    origin_(origin),
+    pointEdgePoint(origin, distSqr),
-    distSqr_(distSqr),
     s_(s),
     v_(v)
 {}
@@ -162,8 +55,7 @@ inline Foam::pointData::pointData
 // Construct as copy
 inline Foam::pointData::pointData(const pointData& wpt)
 :
-    origin_(wpt.origin()),
+    pointEdgePoint(wpt),
-    distSqr_(wpt.distSqr()),
     s_(wpt.s()),
     v_(wpt.v())
 {}
@@ -171,18 +63,6 @@ inline Foam::pointData::pointData(const pointData& wpt)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline const Foam::point& Foam::pointData::origin() const
-{
-    return origin_;
-}
-
-
-inline Foam::scalar Foam::pointData::distSqr() const
-{
-    return distSqr_;
-}
-
-
 inline Foam::scalar Foam::pointData::s() const
 {
     return s_;
@@ -195,157 +75,143 @@ inline const Foam::vector& Foam::pointData::v() const
 }
 
 
-inline bool Foam::pointData::valid() const
+template <class TrackingData>
-{
+inline void Foam::pointData::transform
-    return origin_ != point::max;
-}
-
-
-// Checks for cyclic points
-inline bool Foam::pointData::sameGeometry
 (
-    const pointData& w2,
+    const tensor& rotTensor,
-    const scalar tol
+    TrackingData& td
-) const
-{
-    scalar diff = Foam::mag(distSqr() - w2.distSqr());
-
-    if (diff < SMALL)
-    {
-        return true;
-    }
-    else
-    {
-        if ((distSqr() > SMALL) && ((diff/distSqr()) < tol))
-        {
-            return true;
-        }
-        else
-        {
-            return false;
-        }
-    }
-}
-
-
-inline void Foam::pointData::leaveDomain
-(
-    const polyPatch& patch,
-    const label patchPointI,
-    const point& coord
 )
 {
-    origin_ -= coord;
+    pointEdgePoint::transform(rotTensor, td);
-}
+    v_ = Foam::transform(rotTensor, v_);
-
-
-inline void Foam::pointData::transform(const tensor& rotTensor)
-{
-    origin_ = Foam::transform(rotTensor, origin_);
-}
-
-
-// Update absolute geometric quantities. Note that distance (distSqr_)
-// is not affected by leaving/entering domain.
-inline void Foam::pointData::enterDomain
-(
-    const polyPatch& patch,
-    const label patchPointI,
-    const point& coord
-)
-{
-    // back to absolute form
-    origin_ += coord;
 }
 
 
 // Update this with information from connected edge
+template <class TrackingData>
 inline bool Foam::pointData::updatePoint
 (
     const polyMesh& mesh,
     const label pointI,
     const label edgeI,
     const pointData& edgeInfo,
-    const scalar tol
+    const scalar tol,
+    TrackingData& td
 )
 {
-    return
+    if
-        update
+    (
+        pointEdgePoint::updatePoint
         (
-            mesh.points()[pointI],
+            mesh,
+            pointI,
+            edgeI,
             edgeInfo,
-            tol
+            tol,
-        );
+            td
+        )
+    )
+    {
+        s_ = edgeInfo.s_;
+        v_ = edgeInfo.v_;
+        return true;
     }
-
+    else
+    {
+        return false;
+    }
+}
 
 // Update this with new information on same point
+template <class TrackingData>
 inline bool Foam::pointData::updatePoint
 (
     const polyMesh& mesh,
     const label pointI,
     const pointData& newPointInfo,
-    const scalar tol
+    const scalar tol,
+    TrackingData& td
 )
 {
-    return
+    if
-        update
+    (
+        pointEdgePoint::updatePoint
         (
-            mesh.points()[pointI],
+            mesh,
+            pointI,
             newPointInfo,
-            tol
+            tol,
-        );
+            td
+        )
+    )
+    {
+        s_ = newPointInfo.s_;
+        v_ = newPointInfo.v_;
+        return true;
+    }
+    else
+    {
+        return false;
+    }
 }
 
 
 // Update this with new information on same point. No extra information.
+template <class TrackingData>
 inline bool Foam::pointData::updatePoint
 (
     const pointData& newPointInfo,
-    const scalar tol
+    const scalar tol,
+    TrackingData& td
 )
 {
-    return update(newPointInfo, tol);
+    if (pointEdgePoint::updatePoint(newPointInfo, tol, td))
+    {
+        s_ = newPointInfo.s_;
+        v_ = newPointInfo.v_;
+        return true;
+    }
+    else
+    {
+        return false;
+    }
 }
 
 
 // Update this with information from connected point
+template <class TrackingData>
 inline bool Foam::pointData::updateEdge
 (
     const polyMesh& mesh,
     const label edgeI,
     const label pointI,
     const pointData& pointInfo,
-    const scalar tol
+    const scalar tol,
+    TrackingData& td
+
 )
 {
-    const pointField& points = mesh.points();
+    if
-
+    (
-    const edge& e = mesh.edges()[edgeI];
+        pointEdgePoint::updateEdge
-
-    const point edgeMid(0.5*(points[e[0]] + points[e[1]]));
-
-    return
-        update
         (
-            edgeMid,
+            mesh,
+            edgeI,
+            pointI,
             pointInfo,
-            tol
+            tol,
-        );
+            td
-}
+        )
-
+    )
-
+    {
-// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+        s_ = pointInfo.s_;
-
+        v_ = pointInfo.v_;
-inline bool Foam::pointData::operator==(const pointData& rhs) const
+        return true;
-{
+    }
-    return origin() == rhs.origin();
+    else
-}
+    {
-
+        return false;
-
+    }
-inline bool Foam::pointData::operator!=(const pointData& rhs) const
-{
-    return !(*this == rhs);
 }
 
 

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

@@ -481,11 +481,6 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::doRemoveCells
         mesh_.clearOut();
     }
 
-    if (overwrite_)
-    {
-        mesh_.setInstance(oldInstance_);
-    }
-
     // Update local mesh data
     cellRemover.updateMesh(map);
 
@@ -1903,9 +1898,14 @@ void Foam::meshRefinement::distribute(const mapDistributePolyMesh& map)
             pointMap.clear();
         }
     }
+
+    // If necessary reset the instance
+    mesh_.setInstance(timeName());
+    setInstance(mesh_.facesInstance());
 }
 
 
+// Update local data for a mesh change
 void Foam::meshRefinement::updateMesh
 (
     const mapPolyMesh& map,
@@ -2027,6 +2027,10 @@ void Foam::meshRefinement::updateMesh
             data.transfer(newFaceData);
         }
     }
+
+    // If necessary reset the instance
+    mesh_.setInstance(timeName());
+    setInstance(mesh_.facesInstance());
 }
 
 
@@ -2137,7 +2141,7 @@ void Foam::meshRefinement::dumpRefinementLevel() const
         IOobject
         (
             "cellLevel",
-            timeName(),
+            mesh_.time().timeName(),// Dump to current time, not to mesh inst
             mesh_,
             IOobject::NO_READ,
             IOobject::AUTO_WRITE,
@@ -2165,7 +2169,7 @@ void Foam::meshRefinement::dumpRefinementLevel() const
         IOobject
         (
             "pointLevel",
-            timeName(),
+            mesh_.time().timeName(),// Dump to current time, not to mesh inst
             mesh_,
             IOobject::NO_READ,
             IOobject::NO_WRITE,

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

@@ -517,11 +517,6 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createBaffles
         mesh_.clearOut();
     }
 
-    if (overwrite())
-    {
-        mesh_.setInstance(oldInstance());
-    }
-
     //- Redo the intersections on the newly create baffle faces. Note that
     //  this changes also the cell centre positions.
     faceSet baffledFacesSet(mesh_, "baffledFacesSet", 2*nBaffles);
@@ -965,11 +960,6 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::mergeBaffles
         mesh_.clearOut();
     }
 
-    if (overwrite())
-    {
-        mesh_.setInstance(oldInstance());
-    }
-
     // Update intersections. Recalculate intersections on merged faces since
     // this seems to give problems? Note: should not be nessecary since
     // baffles preserve intersections from when they were created.
@@ -2249,11 +2239,6 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::dupNonManifoldPoints()
         mesh_.clearOut();
     }
 
-    if (overwrite())
-    {
-        mesh_.setInstance(oldInstance());
-    }
-
     // Update intersections. Is mapping only (no faces created, positions stay
     // same) so no need to recalculate intersections.
     updateMesh(map, labelList(0));
@@ -2783,11 +2768,6 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
         mesh_.clearOut();
     }
 
-    if (overwrite())
-    {
-        mesh_.setInstance(oldInstance());
-    }
-
     // Print some stats (note: zones are synchronised)
     if (mesh_.cellZones().size() > 0)
     {

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

@@ -103,11 +103,6 @@ Foam::label Foam::meshRefinement::mergePatchFaces
             mesh_.clearOut();
         }
 
-        if (overwrite())
-        {
-            mesh_.setInstance(oldInstance());
-        }
-
         faceCombiner.updateMesh(map);
 
         // Get the kept faces that need to be recalculated.
@@ -203,11 +198,6 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::mergeEdges
             mesh_.clearOut();
         }
 
-        if (overwrite())
-        {
-            mesh_.setInstance(oldInstance());
-        }
-
         pointRemover.updateMesh(map);
 
         // Get the kept faces that need to be recalculated.

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

@@ -1236,11 +1236,6 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::refine
         mesh_.clearOut();
     }
 
-    if (overwrite())
-    {
-        mesh_.setInstance(oldInstance());
-    }
-
     // Update intersection info
     updateMesh(map, getChangedFaces(map, cellsToRefine));
 

src/meshTools/PointEdgeWave/PointEdgeWave.C

@@ -32,36 +32,25 @@ License
 #include "PstreamCombineReduceOps.H"
 #include "debug.H"
 #include "typeInfo.H"
+#include "globalMeshData.H"
+#include "pointFields.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-template <class Type>
+template <class Type, class TrackingData>
-Foam::scalar Foam::PointEdgeWave<Type>::propagationTol_ = 0.01;
+Foam::scalar Foam::PointEdgeWave<Type, TrackingData>::propagationTol_ = 0.01;
 
-
+template <class Type, class TrackingData>
-// Offset labelList. Used for transferring from one cyclic half to the other.
+Foam::label Foam::PointEdgeWave<Type, TrackingData>::dummyTrackData_ = 12345;
-template <class Type>
-void Foam::PointEdgeWave<Type>::offset(const label val, labelList& elems)
-{
-    forAll(elems, i)
-    {
-        elems[i] += val;
-    }
-}
 
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-// Gets point-point correspondence. Is
-// - list of halfA points (in cyclic patch points)
-// - list of halfB points (can overlap with A!)
-// - for every patchPoint its corresponding point
 // Handle leaving domain. Implementation referred to Type
-template <class Type>
+template <class Type, class TrackingData>
-void Foam::PointEdgeWave<Type>::leaveDomain
+void Foam::PointEdgeWave<Type, TrackingData>::leaveDomain
 (
-    const polyPatch& meshPatch,
+    const polyPatch& patch,
-    const primitivePatch& patch,
     const labelList& patchPointLabels,
     List<Type>& pointInfo
 ) const
@@ -74,17 +63,16 @@ void Foam::PointEdgeWave<Type>::leaveDomain
 
         const point& pt = patch.points()[meshPoints[patchPointI]];
 
-        pointInfo[i].leaveDomain(meshPatch, patchPointI, pt);
+        pointInfo[i].leaveDomain(patch, patchPointI, pt, td_);
     }
 }
 
 
 // Handle entering domain. Implementation referred to Type
-template <class Type>
+template <class Type, class TrackingData>
-void Foam::PointEdgeWave<Type>::enterDomain
+void Foam::PointEdgeWave<Type, TrackingData>::enterDomain
 (
-    const polyPatch& meshPatch,
+    const polyPatch& patch,
-    const primitivePatch& patch,
     const labelList& patchPointLabels,
     List<Type>& pointInfo
 ) const
@@ -97,15 +85,16 @@ void Foam::PointEdgeWave<Type>::enterDomain
 
         const point& pt = patch.points()[meshPoints[patchPointI]];
 
-        pointInfo[i].enterDomain(meshPatch, patchPointI, pt);
+        pointInfo[i].enterDomain(patch, patchPointI, pt, td_);
     }
 }
 
 
 // Transform. Implementation referred to Type
-template <class Type>
+template <class Type, class TrackingData>
-void Foam::PointEdgeWave<Type>::transform
+void Foam::PointEdgeWave<Type, TrackingData>::transform
 (
+    const polyPatch& patch,
     const tensorField& rotTensor,
     List<Type>& pointInfo
 ) const
@@ -116,19 +105,22 @@ void Foam::PointEdgeWave<Type>::transform
 
         forAll(pointInfo, i)
         {
-            pointInfo[i].transform(T);
+            pointInfo[i].transform(T, td_);
         }
     }
     else
     {
         FatalErrorIn
         (
-            "PointEdgeWave<Type>::transform(const tensorField&, List<Type>&)"
+            "PointEdgeWave<Type, TrackingData>::transform"
-        )   << "Parallel cyclics not supported" << abort(FatalError);
+            "(const tensorField&, List<Type>&)"
+        )   << "Non-uniform transformation on patch " << patch.name()
+            << " not supported for point fields"
+            << abort(FatalError);
 
         forAll(pointInfo, i)
         {
-            pointInfo[i].transform(rotTensor[i]);
+            pointInfo[i].transform(rotTensor[i], td_);
         }
     }
 }
@@ -139,19 +131,18 @@ void Foam::PointEdgeWave<Type>::transform
 // Updates:
 //      - changedPoint_, changedPoints_, nChangedPoints_,
 //      - statistics: nEvals_, nUnvisitedPoints_
-template <class Type>
+template <class Type, class TrackingData>
-bool Foam::PointEdgeWave<Type>::updatePoint
+bool Foam::PointEdgeWave<Type, TrackingData>::updatePoint
 (
     const label pointI,
     const label neighbourEdgeI,
     const Type& neighbourInfo,
-    const scalar tol,
     Type& pointInfo
 )
 {
     nEvals_++;
 
-    bool wasValid = pointInfo.valid();
+    bool wasValid = pointInfo.valid(td_);
 
     bool propagate =
         pointInfo.updatePoint
@@ -160,7 +151,8 @@ bool Foam::PointEdgeWave<Type>::updatePoint
             pointI,
             neighbourEdgeI,
             neighbourInfo,
-            tol
+            propagationTol_,
+            td_
         );
 
     if (propagate)
@@ -172,7 +164,7 @@ bool Foam::PointEdgeWave<Type>::updatePoint
         }
     }
 
-    if (!wasValid && pointInfo.valid())
+    if (!wasValid && pointInfo.valid(td_))
     {
         --nUnvisitedPoints_;
     }
@@ -186,18 +178,17 @@ bool Foam::PointEdgeWave<Type>::updatePoint
 // Updates:
 //      - changedPoint_, changedPoints_, nChangedPoints_,
 //      - statistics: nEvals_, nUnvisitedPoints_
-template <class Type>
+template <class Type, class TrackingData>
-bool Foam::PointEdgeWave<Type>::updatePoint
+bool Foam::PointEdgeWave<Type, TrackingData>::updatePoint
 (
     const label pointI,
     const Type& neighbourInfo,
-    const scalar tol,
     Type& pointInfo
 )
 {
     nEvals_++;
 
-    bool wasValid = pointInfo.valid();
+    bool wasValid = pointInfo.valid(td_);
 
     bool propagate =
         pointInfo.updatePoint
@@ -205,7 +196,8 @@ bool Foam::PointEdgeWave<Type>::updatePoint
             mesh_,
             pointI,
             neighbourInfo,
-            tol
+            propagationTol_,
+            td_
         );
 
     if (propagate)
@@ -217,7 +209,7 @@ bool Foam::PointEdgeWave<Type>::updatePoint
         }
     }
 
-    if (!wasValid && pointInfo.valid())
+    if (!wasValid && pointInfo.valid(td_))
     {
         --nUnvisitedPoints_;
     }
@@ -231,19 +223,18 @@ bool Foam::PointEdgeWave<Type>::updatePoint
 // Updates:
 //      - changedEdge_, changedEdges_, nChangedEdges_,
 //      - statistics: nEvals_, nUnvisitedEdge_
-template <class Type>
+template <class Type, class TrackingData>
-bool Foam::PointEdgeWave<Type>::updateEdge
+bool Foam::PointEdgeWave<Type, TrackingData>::updateEdge
 (
     const label edgeI,
     const label neighbourPointI,
     const Type& neighbourInfo,
-    const scalar tol,
     Type& edgeInfo
 )
 {
     nEvals_++;
 
-    bool wasValid = edgeInfo.valid();
+    bool wasValid = edgeInfo.valid(td_);
 
     bool propagate =
         edgeInfo.updateEdge
@@ -252,7 +243,8 @@ bool Foam::PointEdgeWave<Type>::updateEdge
             edgeI,
             neighbourPointI,
             neighbourInfo,
-            tol
+            propagationTol_,
+            td_
         );
 
     if (propagate)
@@ -264,7 +256,7 @@ bool Foam::PointEdgeWave<Type>::updateEdge
         }
     }
 
-    if (!wasValid && edgeInfo.valid())
+    if (!wasValid && edgeInfo.valid(td_))
     {
         --nUnvisitedEdges_;
     }
@@ -274,9 +266,9 @@ bool Foam::PointEdgeWave<Type>::updateEdge
 
 
 // Check if patches of given type name are present
-template <class Type>
+template <class Type, class TrackingData>
 template <class PatchType>
-Foam::label Foam::PointEdgeWave<Type>::countPatchType() const
+Foam::label Foam::PointEdgeWave<Type, TrackingData>::countPatchType() const
 {
     label nPatches = 0;
 
@@ -291,199 +283,110 @@ Foam::label Foam::PointEdgeWave<Type>::countPatchType() const
 }
 
 
-// Collect changed patch points
-template <class Type>
-void Foam::PointEdgeWave<Type>::getChangedPatchPoints
-(
-    const primitivePatch& patch,
-
-    DynamicList<Type>& patchInfo,
-    DynamicList<label>& patchPoints,
-    DynamicList<label>& owner,
-    DynamicList<label>& ownerIndex
-) const
-{
-    const labelList& meshPoints = patch.meshPoints();
-    const faceList& localFaces = patch.localFaces();
-    const labelListList& pointFaces = patch.pointFaces();
-
-    forAll(meshPoints, patchPointI)
-    {
-        label meshPointI = meshPoints[patchPointI];
-
-        if (changedPoint_[meshPointI])
-        {
-            patchInfo.append(allPointInfo_[meshPointI]);
-            patchPoints.append(patchPointI);
-
-            label patchFaceI = pointFaces[patchPointI][0];
-
-            const face& f = localFaces[patchFaceI];
-
-            label index = findIndex(f, patchPointI);
-
-            owner.append(patchFaceI);
-            ownerIndex.append(index);
-        }
-    }
-
-    patchInfo.shrink();
-    patchPoints.shrink();
-    owner.shrink();
-    ownerIndex.shrink();
-}
-
-
-// Update overall for changed patch points
-template <class Type>
-void Foam::PointEdgeWave<Type>::updateFromPatchInfo
-(
-    const polyPatch& meshPatch,
-    const primitivePatch& patch,
-    const labelList& owner,
-    const labelList& ownerIndex,
-    List<Type>& patchInfo
-)
-{
-    const faceList& localFaces = patch.localFaces();
-    const labelList& meshPoints = patch.meshPoints();
-
-    // Get patch and mesh points.
-    labelList changedPatchPoints(patchInfo.size());
-    labelList changedMeshPoints(patchInfo.size());
-
-    forAll(owner, i)
-    {
-        label faceI = owner[i];
-
-        const face& f = localFaces[faceI];
-
-        label index = (f.size() - ownerIndex[i]) % f.size();
-
-        changedPatchPoints[i] = f[index];
-        changedMeshPoints[i] = meshPoints[f[index]];
-    }
-
-    // Adapt for entering domain
-    enterDomain(meshPatch, patch, changedPatchPoints, patchInfo);
-
-    // Merge received info
-    forAll(patchInfo, i)
-    {
-        updatePoint
-        (
-            changedMeshPoints[i],
-            patchInfo[i],
-            propagationTol_,
-            allPointInfo_[changedMeshPoints[i]]
-        );
-    }
-}
-
-
-
 // Transfer all the information to/from neighbouring processors
-template <class Type>
+template <class Type, class TrackingData>
-void Foam::PointEdgeWave<Type>::handleProcPatches()
+void Foam::PointEdgeWave<Type, TrackingData>::handleProcPatches()
 {
-    // 1. Send all point info on processor patches. Send as
+    // 1. Send all point info on processor patches.
-    // face label + offset in face.
 
-    forAll(mesh_.boundaryMesh(), patchI)
+    PstreamBuffers pBufs(Pstream::nonBlocking);
+
+    DynamicList<Type> patchInfo;
+    DynamicList<label> thisPoints;
+    DynamicList<label> nbrPoints;
+
+    forAll(mesh_.globalData().processorPatches(), i)
     {
-        const polyPatch& patch = mesh_.boundaryMesh()[patchI];
+        label patchI = mesh_.globalData().processorPatches()[i];
+        const processorPolyPatch& procPatch =
+            refCast<const processorPolyPatch>(mesh_.boundaryMesh()[patchI]);
 
-        if (isA<processorPolyPatch>(patch))
+        patchInfo.clear();
+        patchInfo.reserve(procPatch.nPoints());
+        thisPoints.clear();
+        thisPoints.reserve(procPatch.nPoints());
+        nbrPoints.clear();
+        nbrPoints.reserve(procPatch.nPoints());
+
+        // Get all changed points in reverse order
+        const labelList& neighbPoints = procPatch.neighbPoints();
+        forAll(neighbPoints, thisPointI)
         {
-            // Get all changed points in relative addressing
+            label meshPointI = procPatch.meshPoints()[thisPointI];
-
+            if (changedPoint_[meshPointI])
-            DynamicList<Type> patchInfo(patch.nPoints());
-            DynamicList<label> patchPoints(patch.nPoints());
-            DynamicList<label> owner(patch.nPoints());
-            DynamicList<label> ownerIndex(patch.nPoints());
-
-            getChangedPatchPoints
-            (
-                patch,
-                patchInfo,
-                patchPoints,
-                owner,
-                ownerIndex
-            );
-
-            // Adapt for leaving domain
-            leaveDomain(patch, patch, patchPoints, patchInfo);
-
-            const processorPolyPatch& procPatch =
-                refCast<const processorPolyPatch>(patch);
-
-            if (debug)
             {
-                Pout<< "Processor patch " << patchI << ' ' << patch.name()
+                patchInfo.append(allPointInfo_[meshPointI]);
-                    << " communicating with " << procPatch.neighbProcNo()
+                thisPoints.append(thisPointI);
-                    << "  Sending:" << patchInfo.size() << endl;
+                nbrPoints.append(neighbPoints[thisPointI]);
-            }
-
-            {
-                OPstream toNeighbour
-                (
-                    Pstream::blocking,
-                    procPatch.neighbProcNo()
-                );
-
-                toNeighbour << owner << ownerIndex << patchInfo;
             }
         }
+
+        // Adapt for leaving domain
+        leaveDomain(procPatch, thisPoints, patchInfo);
+
+        if (debug)
+        {
+            Pout<< "Processor patch " << patchI << ' ' << procPatch.name()
+                << " communicating with " << procPatch.neighbProcNo()
+                << "  Sending:" << patchInfo.size() << endl;
+        }
+
+        UOPstream toNeighbour(procPatch.neighbProcNo(), pBufs);
+        toNeighbour << nbrPoints << patchInfo;
     }
 
 
+    pBufs.finishedSends();
+
     //
     // 2. Receive all point info on processor patches.
     //
 
-    forAll(mesh_.boundaryMesh(), patchI)
+    forAll(mesh_.globalData().processorPatches(), i)
     {
-        const polyPatch& patch = mesh_.boundaryMesh()[patchI];
+        label patchI = mesh_.globalData().processorPatches()[i];
+        const processorPolyPatch& procPatch =
+            refCast<const processorPolyPatch>(mesh_.boundaryMesh()[patchI]);
+
+        List<Type> patchInfo;
+        labelList patchPoints;
 
-        if (isA<processorPolyPatch>(patch))
         {
-            const processorPolyPatch& procPatch =
+            UIPstream fromNeighbour(procPatch.neighbProcNo(), pBufs);
-                refCast<const processorPolyPatch>(patch);
+            fromNeighbour >> patchPoints >> patchInfo;
+        }
 
-            List<Type> patchInfo;
+        if (debug)
-            labelList owner;
+        {
-            labelList ownerIndex;
+            Pout<< "Processor patch " << patchI << ' ' << procPatch.name()
+                << " communicating with " << procPatch.neighbProcNo()
+                << "  Received:" << patchInfo.size() << endl;
+        }
+
+        // Apply transform to received data for non-parallel planes
+        if (!procPatch.parallel())
+        {
+            transform(procPatch, procPatch.forwardT(), patchInfo);
+        }
+
+        // Adapt for entering domain
+        enterDomain(procPatch, patchPoints, patchInfo);
+
+        // Merge received info
+        const labelList& meshPoints = procPatch.meshPoints();
+        forAll(patchInfo, i)
+        {
+            label meshPointI = meshPoints[patchPoints[i]];
+
+            if (!allPointInfo_[meshPointI].equal(patchInfo[i], td_))
             {
-                IPstream fromNeighbour
+                updatePoint
                 (
-                    Pstream::blocking,
+                    meshPointI,
-                    procPatch.neighbProcNo()
+                    patchInfo[i],
+                    allPointInfo_[meshPointI]
                 );
-
-                fromNeighbour >> owner >> ownerIndex >> patchInfo;
             }
-
-            if (debug)
-            {
-                Pout<< "Processor patch " << patchI << ' ' << patch.name()
-                    << " communicating with " << procPatch.neighbProcNo()
-                    << "  Received:" << patchInfo.size() << endl;
-            }
-
-            // Apply transform to received data for non-parallel planes
-            if (!procPatch.parallel())
-            {
-                transform(procPatch.forwardT(), patchInfo);
-            }
-
-            updateFromPatchInfo
-            (
-                patch,
-                patch,
-                owner,
-                ownerIndex,
-                patchInfo
-            );
         }
     }
 
@@ -508,20 +411,21 @@ void Foam::PointEdgeWave<Type>::handleProcPatches()
 
     // Combine on master. Reduce operator has to handle a list and call
     // Type.updatePoint for all elements
-    combineReduce(sharedData, listUpdateOp<Type>());
+    combineReduce(sharedData, listUpdateOp<Type>(propagationTol_, td_));
 
     forAll(pd.sharedPointLabels(), i)
     {
         label meshPointI = pd.sharedPointLabels()[i];
 
-        // Retrieve my entries from the shared points
+        // Retrieve my entries from the shared points.
-        if (sharedData[pd.sharedPointAddr()[i]].valid())
+        const Type& nbrInfo = sharedData[pd.sharedPointAddr()[i]];
+
+        if (!allPointInfo_[meshPointI].equal(nbrInfo, td_))
         {
             updatePoint
             (
                 meshPointI,
-                sharedData[pd.sharedPointAddr()[i]],
+                nbrInfo,
-                propagationTol_,
                 allPointInfo_[meshPointI]
             );
         }
@@ -529,11 +433,14 @@ void Foam::PointEdgeWave<Type>::handleProcPatches()
 }
 
 
-template <class Type>
+template <class Type, class TrackingData>
-void Foam::PointEdgeWave<Type>::handleCyclicPatches()
+void Foam::PointEdgeWave<Type, TrackingData>::handleCyclicPatches()
 {
-    // 1. Send all point info on cyclic patches. Send as
+    // 1. Send all point info on cyclic patches.
-    // face label + offset in face.
+
+    DynamicList<Type> nbrInfo;
+    DynamicList<label> nbrPoints;
+    DynamicList<label> thisPoints;
 
     forAll(mesh_.boundaryMesh(), patchI)
     {
@@ -544,31 +451,44 @@ void Foam::PointEdgeWave<Type>::handleCyclicPatches()
             const cyclicPolyPatch& cycPatch =
                 refCast<const cyclicPolyPatch>(patch);
 
-            const cyclicPolyPatch& nbrPatch = cycPatch.neighbPatch();
+            nbrInfo.clear();
+            nbrInfo.reserve(cycPatch.nPoints());
+            nbrPoints.clear();
+            nbrPoints.reserve(cycPatch.nPoints());
+            thisPoints.clear();
+            thisPoints.reserve(cycPatch.nPoints());
 
-            DynamicList<Type> nbrInfo(nbrPatch.nPoints());
+            // Collect nbrPatch points that have changed
-            DynamicList<label> nbrPoints(nbrPatch.nPoints());
+            {
-            DynamicList<label> nbrOwner(nbrPatch.nPoints());
+                const cyclicPolyPatch& nbrPatch = cycPatch.neighbPatch();
-            DynamicList<label> nbrIndex(nbrPatch.nPoints());
+                const edgeList& pairs = cycPatch.coupledPoints();
+                const labelList& meshPoints = nbrPatch.meshPoints();
 
-            getChangedPatchPoints
+                forAll(pairs, pairI)
-            (
+                {
-                nbrPatch,
+                    label thisPointI = pairs[pairI][0];
-                nbrInfo,
+                    label nbrPointI = pairs[pairI][1];
-                nbrPoints,
+                    label meshPointI = meshPoints[nbrPointI];
-                nbrOwner,
+
-                nbrIndex
+                    if (changedPoint_[meshPointI])
-            );
+                    {
+                        nbrInfo.append(allPointInfo_[meshPointI]);
+                        nbrPoints.append(nbrPointI);
+                        thisPoints.append(thisPointI);
+                    }
+                }
+
+                // nbr : adapt for leaving domain
+                leaveDomain(nbrPatch, nbrPoints, nbrInfo);
+            }
 
-            // nbr : adapt for leaving domain
-            leaveDomain(nbrPatch, nbrPatch, nbrPoints, nbrInfo);
 
             // Apply rotation for non-parallel planes
 
             if (!cycPatch.parallel())
             {
                 // received data from half1
-                transform(cycPatch.forwardT(), nbrInfo);
+                transform(cycPatch, cycPatch.forwardT(), nbrInfo);
             }
 
             if (debug)
@@ -578,19 +498,24 @@ void Foam::PointEdgeWave<Type>::handleCyclicPatches()
                     << endl;
             }
 
-            // Half1: update with data from halfB
+            // Adapt for entering domain
-            updateFromPatchInfo
+            enterDomain(cycPatch, thisPoints, nbrInfo);
-            (
-                cycPatch,
-                cycPatch,
-                nbrOwner,
-                nbrIndex,
-                nbrInfo
-            );
 
-            if (debug)
+            // Merge received info
+            const labelList& meshPoints = cycPatch.meshPoints();
+            forAll(nbrInfo, i)
             {
-                //checkCyclic(patch);
+                label meshPointI = meshPoints[thisPoints[i]];
+
+                if (!allPointInfo_[meshPointI].equal(nbrInfo[i], td_))
+                {
+                    updatePoint
+                    (
+                        meshPointI,
+                        nbrInfo[i],
+                        allPointInfo_[meshPointI]
+                    );
+                }
             }
         }
     }
@@ -601,22 +526,24 @@ void Foam::PointEdgeWave<Type>::handleCyclicPatches()
 
 // Iterate, propagating changedPointsInfo across mesh, until no change (or
 // maxIter reached). Initial point values specified.
-template <class Type>
+template <class Type, class TrackingData>
-Foam::PointEdgeWave<Type>::PointEdgeWave
+Foam::PointEdgeWave<Type, TrackingData>::PointEdgeWave
 (
     const polyMesh& mesh,
     const labelList& changedPoints,
     const List<Type>& changedPointsInfo,
 
-    List<Type>& allPointInfo,
+    UList<Type>& allPointInfo,
-    List<Type>& allEdgeInfo,
+    UList<Type>& allEdgeInfo,
 
-    const label maxIter
+    const label maxIter,
+    TrackingData& td
 )
 :
     mesh_(mesh),
     allPointInfo_(allPointInfo),
     allEdgeInfo_(allEdgeInfo),
+    td_(td),
     changedPoint_(mesh_.nPoints(), false),
     changedPoints_(mesh_.nPoints()),
     nChangedPoints_(0),
@@ -632,7 +559,7 @@ Foam::PointEdgeWave<Type>::PointEdgeWave
     {
         FatalErrorIn
         (
-            "PointEdgeWave<Type>::PointEdgeWave"
+            "PointEdgeWave<Type, TrackingData>::PointEdgeWave"
             "(const polyMesh&, const labelList&, const List<Type>,"
             " List<Type>&, List<Type>&, const label maxIter)"
         )   << "size of pointInfo work array is not equal to the number"
@@ -645,7 +572,7 @@ Foam::PointEdgeWave<Type>::PointEdgeWave
     {
         FatalErrorIn
         (
-            "PointEdgeWave<Type>::PointEdgeWave"
+            "PointEdgeWave<Type, TrackingData>::PointEdgeWave"
             "(const polyMesh&, const labelList&, const List<Type>,"
             " List<Type>&, List<Type>&, const label maxIter)"
         )   << "size of edgeInfo work array is not equal to the number"
@@ -671,7 +598,7 @@ Foam::PointEdgeWave<Type>::PointEdgeWave
     {
         FatalErrorIn
         (
-            "PointEdgeWave<Type>::PointEdgeWave"
+            "PointEdgeWave<Type, TrackingData>::PointEdgeWave"
             "(const polyMesh&, const labelList&, const List<Type>,"
             " List<Type>&, List<Type>&, const label maxIter)"
         )   << "Maximum number of iterations reached. Increase maxIter." << endl
@@ -683,17 +610,19 @@ Foam::PointEdgeWave<Type>::PointEdgeWave
 }
 
 
-template <class Type>
+template <class Type, class TrackingData>
-Foam::PointEdgeWave<Type>::PointEdgeWave
+Foam::PointEdgeWave<Type, TrackingData>::PointEdgeWave
 (
     const polyMesh& mesh,
-    List<Type>& allPointInfo,
+    UList<Type>& allPointInfo,
-    List<Type>& allEdgeInfo
+    UList<Type>& allEdgeInfo,
+    TrackingData& td
 )
 :
     mesh_(mesh),
     allPointInfo_(allPointInfo),
     allEdgeInfo_(allEdgeInfo),
+    td_(td),
     changedPoint_(mesh_.nPoints(), false),
     changedPoints_(mesh_.nPoints()),
     nChangedPoints_(0),
@@ -709,31 +638,31 @@ Foam::PointEdgeWave<Type>::PointEdgeWave
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-template <class Type>
+template <class Type, class TrackingData>
-Foam::PointEdgeWave<Type>::~PointEdgeWave()
+Foam::PointEdgeWave<Type, TrackingData>::~PointEdgeWave()
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 
-template <class Type>
+template <class Type, class TrackingData>
-Foam::label Foam::PointEdgeWave<Type>::getUnsetPoints() const
+Foam::label Foam::PointEdgeWave<Type, TrackingData>::getUnsetPoints() const
 {
     return nUnvisitedPoints_;
 }
 
 
-template <class Type>
+template <class Type, class TrackingData>
-Foam::label Foam::PointEdgeWave<Type>::getUnsetEdges() const
+Foam::label Foam::PointEdgeWave<Type, TrackingData>::getUnsetEdges() const
 {
     return nUnvisitedEdges_;
 }
 
 
 // Copy point information into member data
-template <class Type>
+template <class Type, class TrackingData>
-void Foam::PointEdgeWave<Type>::setPointInfo
+void Foam::PointEdgeWave<Type, TrackingData>::setPointInfo
 (
     const labelList& changedPoints,
     const List<Type>& changedPointsInfo
@@ -743,13 +672,13 @@ void Foam::PointEdgeWave<Type>::setPointInfo
     {
         label pointI = changedPoints[changedPointI];
 
-        bool wasValid = allPointInfo_[pointI].valid();
+        bool wasValid = allPointInfo_[pointI].valid(td_);
 
         // Copy info for pointI
         allPointInfo_[pointI] = changedPointsInfo[changedPointI];
 
         // Maintain count of unset points
-        if (!wasValid && allPointInfo_[pointI].valid())
+        if (!wasValid && allPointInfo_[pointI].valid(td_))
         {
             --nUnvisitedPoints_;
         }
@@ -766,8 +695,8 @@ void Foam::PointEdgeWave<Type>::setPointInfo
 
 
 // Propagate information from edge to point. Return number of points changed.
-template <class Type>
+template <class Type, class TrackingData>
-Foam::label Foam::PointEdgeWave<Type>::edgeToPoint()
+Foam::label Foam::PointEdgeWave<Type, TrackingData>::edgeToPoint()
 {
     for
     (
@@ -780,7 +709,7 @@ Foam::label Foam::PointEdgeWave<Type>::edgeToPoint()
 
         if (!changedEdge_[edgeI])
         {
-            FatalErrorIn("PointEdgeWave<Type>::edgeToPoint()")
+            FatalErrorIn("PointEdgeWave<Type, TrackingData>::edgeToPoint()")
                 << "edge " << edgeI
                 << " not marked as having been changed" << nl
                 << "This might be caused by multiple occurences of the same"
@@ -797,14 +726,13 @@ Foam::label Foam::PointEdgeWave<Type>::edgeToPoint()
         {
             Type& currentWallInfo = allPointInfo_[e[eI]];
 
-            if (currentWallInfo != neighbourWallInfo)
+            if (!currentWallInfo.equal(neighbourWallInfo, td_))
             {
                 updatePoint
                 (
                     e[eI],
                     edgeI,
                     neighbourWallInfo,
-                    propagationTol_,
                     currentWallInfo
                 );
             }
@@ -843,8 +771,8 @@ Foam::label Foam::PointEdgeWave<Type>::edgeToPoint()
 
 
 // Propagate information from point to edge. Return number of edges changed.
-template <class Type>
+template <class Type, class TrackingData>
-Foam::label Foam::PointEdgeWave<Type>::pointToEdge()
+Foam::label Foam::PointEdgeWave<Type, TrackingData>::pointToEdge()
 {
     const labelListList& pointEdges = mesh_.pointEdges();
 
@@ -859,7 +787,7 @@ Foam::label Foam::PointEdgeWave<Type>::pointToEdge()
 
         if (!changedPoint_[pointI])
         {
-            FatalErrorIn("PointEdgeWave<Type>::pointToEdge()")
+            FatalErrorIn("PointEdgeWave<Type, TrackingData>::pointToEdge()")
                 << "Point " << pointI
                 << " not marked as having been changed" << nl
                 << "This might be caused by multiple occurences of the same"
@@ -877,14 +805,13 @@ Foam::label Foam::PointEdgeWave<Type>::pointToEdge()
 
             Type& currentWallInfo = allEdgeInfo_[edgeI];
 
-            if (currentWallInfo != neighbourWallInfo)
+            if (!currentWallInfo.equal(neighbourWallInfo, td_))
             {
                 updateEdge
                 (
                     edgeI,
                     pointI,
                     neighbourWallInfo,
-                    propagationTol_,
                     currentWallInfo
                 );
             }
@@ -912,8 +839,11 @@ Foam::label Foam::PointEdgeWave<Type>::pointToEdge()
 
 
 // Iterate
-template <class Type>
+template <class Type, class TrackingData>
-Foam::label Foam::PointEdgeWave<Type>::iterate(const label maxIter)
+Foam::label Foam::PointEdgeWave<Type, TrackingData>::iterate
+(
+    const label maxIter
+)
 {
     if (nCyclicPatches_ > 0)
     {
@@ -975,4 +905,5 @@ Foam::label Foam::PointEdgeWave<Type>::iterate(const label maxIter)
     return iter;
 }
 
+
 // ************************************************************************* //

src/meshTools/PointEdgeWave/PointEdgeWave.H

@@ -58,13 +58,9 @@ SourceFiles
 #ifndef PointEdgeWave_H
 #define PointEdgeWave_H
 
-#include "label.H"
 #include "boolList.H"
 #include "scalarField.H"
-#include "pointFields.H"
+#include "tensorField.H"
-#include "tensor.H"
-#include "primitivePatch.H"
-#include "PtrList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -73,7 +69,7 @@ namespace Foam
 
 // Forward declaration of classes
 class polyMesh;
-
+class polyPatch;
 
 /*---------------------------------------------------------------------------*\
                         Class PointEdgeWaveName Declaration
@@ -86,7 +82,7 @@ TemplateName(PointEdgeWave);
                            Class PointEdgeWave Declaration
 \*---------------------------------------------------------------------------*/
 
-template <class Type>
+template <class Type, class TrackingData = int>
 class PointEdgeWave
 :
     public PointEdgeWaveName
@@ -98,6 +94,10 @@ class PointEdgeWave
         //  up to Type implementation)
         static scalar propagationTol_;
 
+        //- Used as default trackdata value to satisfy default template
+        //  argument.
+        static label dummyTrackData_;
+
 
     // Private data
 
@@ -105,10 +105,13 @@ class PointEdgeWave
         const polyMesh& mesh_;
 
         //- Wall information for all points
-        List<Type>& allPointInfo_;
+        UList<Type>& allPointInfo_;
 
         //- Information on all mesh edges
-        List<Type>& allEdgeInfo_;
+        UList<Type>& allEdgeInfo_;
+
+        //- Additional data to be passed into container
+        TrackingData& td_;
 
         //- Has point changed
         boolList changedPoint_;
@@ -127,9 +130,6 @@ class PointEdgeWave
         //- Number of cyclic patches
         label nCyclicPatches_;
 
-        //- For every cyclic patch two primitivePatches
-        PtrList<primitivePatch> cycHalves_;
-
         //- Number of evaluations
         label nEvals_;
 
@@ -140,15 +140,10 @@ class PointEdgeWave
 
     // Private Member Functions
 
-        //- Add value to all elements of labelList
-        static void offset(const label val, labelList& elems);
-
-
         //- Adapt pointInfo for leaving domain
         void leaveDomain
         (
-            const polyPatch& meshPatch,
+            const polyPatch&,
-            const primitivePatch& patch,
             const List<label>& patchPointLabels,
             List<Type>& pointInfo
         ) const;
@@ -156,8 +151,7 @@ class PointEdgeWave
         //- Adapt pointInfo for entering domain
         void enterDomain
         (
-            const polyPatch& meshPatch,
+            const polyPatch&,
-            const primitivePatch& patch,
             const List<label>& patchPointLabels,
             List<Type>& pointInfo
         ) const;
@@ -165,6 +159,7 @@ class PointEdgeWave
         //- Transform. Implementation referred to Type
         void transform
         (
+            const polyPatch& patch,
             const tensorField& rotTensor,
             List<Type>& pointInfo
         ) const;
@@ -176,7 +171,6 @@ class PointEdgeWave
             const label pointI,
             const label neighbourEdgeI,
             const Type& neighbourInfo,
-            const scalar tol,
             Type& pointInfo
         );
 
@@ -186,7 +180,6 @@ class PointEdgeWave
         (
             const label pointI,
             const Type& neighbourInfo,
-            const scalar tol,
             Type& pointInfo
         );
 
@@ -197,7 +190,6 @@ class PointEdgeWave
             const label edgeI,
             const label neighbourPointI,
             const Type& neighbourInfo,
-            const scalar tol,
             Type& edgeInfo
         );
 
@@ -207,26 +199,6 @@ class PointEdgeWave
             template <class PatchType>
             label countPatchType() const;
 
-            //- Get info on patch points
-            void getChangedPatchPoints
-            (
-                const primitivePatch& patch,
-                DynamicList<Type>& patchInfo,
-                DynamicList<label>& patchPoints,
-                DynamicList<label>& owner,
-                DynamicList<label>& ownerIndex
-            ) const;
-
-            //- Merge data from patch into overall data
-            void updateFromPatchInfo
-            (
-                const polyPatch& meshPatch,
-                const primitivePatch& patch,
-                const labelList& owner,
-                const labelList& ownerIndex,
-                List<Type>& patchInfo
-            );
-
             //- Merge data from across processor boundaries
             void handleProcPatches();
 
@@ -270,9 +242,10 @@ public:
             const polyMesh& mesh,
             const labelList& initialPoints,
             const List<Type>& initialPointsInfo,
-            List<Type>& allPointInfo,
+            UList<Type>& allPointInfo,
-            List<Type>& allEdgeInfo,
+            UList<Type>& allEdgeInfo,
-            const label maxIter
+            const label maxIter,
+            TrackingData& td = dummyTrackData_
         );
 
         //- Construct from mesh. Use setPointInfo and iterate() to do
@@ -280,8 +253,9 @@ public:
         PointEdgeWave
         (
             const polyMesh& mesh,
-            List<Type>& allPointInfo,
+            UList<Type>& allPointInfo,
-            List<Type>& allEdgeInfo
+            UList<Type>& allEdgeInfo,
+            TrackingData& td = dummyTrackData_
         );
 
 
@@ -291,18 +265,24 @@ public:
 
     // Member Functions
 
-        //- Get allPointInfo
+        //- Access allPointInfo
-        const List<Type>& allPointInfo() const
+        UList<Type>& allPointInfo() const
         {
             return allPointInfo_;
         }
 
-        //- Get allEdgeInfo
+        //- Access allEdgeInfo
-        const List<Type>& allEdgeInfo() const
+        UList<Type>& allEdgeInfo() const
         {
             return allEdgeInfo_;
         }
 
+        //- Additional data to be passed into container
+        const TrackingData& data() const
+        {
+            return td_;
+        }
+
         //- Get number of unvisited edges, i.e. edges that were not (yet)
         //  reached from walking across mesh. This can happen from
         //  - not enough iterations done
@@ -340,19 +320,29 @@ public:
 \*---------------------------------------------------------------------------*/
 
 //- List update operation
-template <class Type>
+template <class Type, class TrackingData = int>
 class listUpdateOp
 {
+    //- Additional data to be passed into container
+
+    const scalar tol_;
+
+    TrackingData& td_;
 
 public:
+    listUpdateOp(const scalar tol, TrackingData& td)
+    :
+        tol_(tol),
+        td_(td)
+    {}
 
     void operator()(List<Type>& x, const List<Type>& y) const
     {
         forAll(x, i)
         {
-            if (y[i].valid())
+            if (y[i].valid(td_))
             {
-                x[i].updatePoint(y[i], PointEdgeWave<Type>::propagationTol());
+                x[i].updatePoint(y[i], tol_, td_);
             }
         }
     }

src/meshTools/PointEdgeWave/pointEdgePoint.H

@@ -73,19 +73,23 @@ class pointEdgePoint
         //- Evaluate distance to point. Update distSqr, origin from whomever
         //  is nearer pt. Return true if w2 is closer to point,
         //  false otherwise.
+        template<class TrackingData>
         inline bool update
         (
             const point&,
             const pointEdgePoint& w2,
-            const scalar tol
+            const scalar tol,
+            TrackingData& td
         );
 
         //- Combine current with w2. Update distSqr, origin if w2 has smaller
         //  quantities and returns true.
+        template<class TrackingData>
         inline bool update
         (
             const pointEdgePoint& w2,
-            const scalar tol
+            const scalar tol,
+            TrackingData& td
         );
 
 
@@ -97,7 +101,7 @@ public:
         inline pointEdgePoint();
 
         //- Construct from origin, distance
-        inline pointEdgePoint(const point& origin, const scalar distSqr);
+        inline pointEdgePoint(const point&, const scalar);
 
         //- Construct as copy
         inline pointEdgePoint(const pointEdgePoint&);
@@ -116,76 +120,102 @@ public:
 
             //- Check whether origin has been changed at all or
             //  still contains original (invalid) value.
-            inline bool valid() const;
+            template<class TrackingData>
+            inline bool valid(TrackingData& td) const;
 
             //- Check for identical geometrical data. Used for cyclics checking.
-            inline bool sameGeometry(const pointEdgePoint&, const scalar tol)
+            template<class TrackingData>
-             const;
+            inline bool sameGeometry
+            (
+                const pointEdgePoint&,
+                const scalar tol,
+                TrackingData& td
+            ) const;
 
             //- Convert origin to relative vector to leaving point
             //  (= point coordinate)
+            template<class TrackingData>
             inline void leaveDomain
             (
                 const polyPatch& patch,
                 const label patchPointI,
-                const point& pos
+                const point& pos,
+                TrackingData& td
             );
 
             //- Convert relative origin to absolute by adding entering point
+            template<class TrackingData>
             inline void enterDomain
             (
                 const polyPatch& patch,
                 const label patchPointI,
-                const point& pos
+                const point& pos,
+                TrackingData& td
             );
 
             //- Apply rotation matrix to origin
-            inline void transform(const tensor& rotTensor);
+            template<class TrackingData>
+            inline void transform
+            (
+                const tensor& rotTensor,
+                TrackingData& td
+            );
 
             //- Influence of edge on point
+            template<class TrackingData>
             inline bool updatePoint
             (
                 const polyMesh& mesh,
                 const label pointI,
                 const label edgeI,
                 const pointEdgePoint& edgeInfo,
-                const scalar tol
+                const scalar tol,
+                TrackingData& td
             );
 
             //- Influence of different value on same point.
             //  Merge new and old info.
+            template<class TrackingData>
             inline bool updatePoint
             (
                 const polyMesh& mesh,
                 const label pointI,
                 const pointEdgePoint& newPointInfo,
-                const scalar tol
+                const scalar tol,
+                TrackingData& td
             );
 
             //- Influence of different value on same point.
             //  No information about current position whatsoever.
+            template<class TrackingData>
             inline bool updatePoint
             (
                 const pointEdgePoint& newPointInfo,
-                const scalar tol
+                const scalar tol,
+                TrackingData& td
             );
 
             //- Influence of point on edge.
+            template<class TrackingData>
             inline bool updateEdge
             (
                 const polyMesh& mesh,
                 const label edgeI,
                 const label pointI,
                 const pointEdgePoint& pointInfo,
-                const scalar tol
+                const scalar tol,
+                TrackingData& td
             );
 
+            //- Same (like operator==)
+            template<class TrackingData>
+            inline bool equal(const pointEdgePoint&, TrackingData& td) const;
+
 
     // Member Operators
 
-        //Note: Used to determine whether to call update.
+        // Needed for List IO
         inline bool operator==(const pointEdgePoint&) const;
-
         inline bool operator!=(const pointEdgePoint&) const;
 
 

src/meshTools/PointEdgeWave/pointEdgePointI.H

@@ -29,16 +29,18 @@ License
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 // Update this with w2 if w2 nearer to pt.
+template<class TrackingData>
 inline bool Foam::pointEdgePoint::update
 (
     const point& pt,
     const pointEdgePoint& w2,
-    const scalar tol
+    const scalar tol,
+    TrackingData& td
 )
 {
     scalar dist2 = magSqr(pt - w2.origin());
 
-    if (!valid())
+    if (!valid(td))
     {
         // current not yet set so use any value
         distSqr_ = dist2;
@@ -72,13 +74,15 @@ inline bool Foam::pointEdgePoint::update
 
 
 // Update this with w2 (information on same point)
+template<class TrackingData>
 inline bool Foam::pointEdgePoint::update
 (
     const pointEdgePoint& w2,
-    const scalar tol
+    const scalar tol,
+    TrackingData& td
 )
 {
-    if (!valid())
+    if (!valid(td))
     {
         // current not yet set so use any value
         distSqr_ = w2.distSqr();
@@ -155,17 +159,20 @@ inline Foam::scalar Foam::pointEdgePoint::distSqr() const
 }
 
 
-inline bool Foam::pointEdgePoint::valid() const
+template<class TrackingData>
+inline bool Foam::pointEdgePoint::valid(TrackingData& td) const
 {
     return origin_ != point::max;
 }
 
 
 // Checks for cyclic points
+template<class TrackingData>
 inline bool Foam::pointEdgePoint::sameGeometry
 (
     const pointEdgePoint& w2,
-    const scalar tol
+    const scalar tol,
+    TrackingData& td
 ) const
 {
     scalar diff = Foam::mag(distSqr() - w2.distSqr());
@@ -188,18 +195,25 @@ inline bool Foam::pointEdgePoint::sameGeometry
 }
 
 
+template<class TrackingData>
 inline void Foam::pointEdgePoint::leaveDomain
 (
     const polyPatch& patch,
     const label patchPointI,
-    const point& coord
+    const point& coord,
+    TrackingData& td
 )
 {
     origin_ -= coord;
 }
 
 
-inline void Foam::pointEdgePoint::transform(const tensor& rotTensor)
+template<class TrackingData>
+inline void Foam::pointEdgePoint::transform
+(
+    const tensor& rotTensor,
+    TrackingData& td
+)
 {
     origin_ = Foam::transform(rotTensor, origin_);
 }
@@ -207,11 +221,13 @@ inline void Foam::pointEdgePoint::transform(const tensor& rotTensor)
 
 // Update absolute geometric quantities. Note that distance (distSqr_)
 // is not affected by leaving/entering domain.
+template<class TrackingData>
 inline void Foam::pointEdgePoint::enterDomain
 (
     const polyPatch& patch,
     const label patchPointI,
-    const point& coord
+    const point& coord,
+    TrackingData& td
 )
 {
     // back to absolute form
@@ -220,78 +236,74 @@ inline void Foam::pointEdgePoint::enterDomain
 
 
 // Update this with information from connected edge
+template<class TrackingData>
 inline bool Foam::pointEdgePoint::updatePoint
 (
     const polyMesh& mesh,
     const label pointI,
     const label edgeI,
     const pointEdgePoint& edgeInfo,
-    const scalar tol
+    const scalar tol,
+    TrackingData& td
 )
 {
-    return
+    return update(mesh.points()[pointI], edgeInfo, tol, td);
-        update
-        (
-            mesh.points()[pointI],
-            edgeInfo,
-            tol
-        );
 }
 
 
 // Update this with new information on same point
+template<class TrackingData>
 inline bool Foam::pointEdgePoint::updatePoint
 (
     const polyMesh& mesh,
     const label pointI,
     const pointEdgePoint& newPointInfo,
-    const scalar tol
+    const scalar tol,
+    TrackingData& td
 )
 {
-    return
+    return update(mesh.points()[pointI], newPointInfo, tol, td);
-        update
-        (
-            mesh.points()[pointI],
-            newPointInfo,
-            tol
-        );
 }
 
 
 // Update this with new information on same point. No extra information.
+template<class TrackingData>
 inline bool Foam::pointEdgePoint::updatePoint
 (
     const pointEdgePoint& newPointInfo,
-    const scalar tol
+    const scalar tol,
+    TrackingData& td
 )
 {
-    return update(newPointInfo, tol);
+    return update(newPointInfo, tol, td);
 }
 
 
 // Update this with information from connected point
+template<class TrackingData>
 inline bool Foam::pointEdgePoint::updateEdge
 (
     const polyMesh& mesh,
     const label edgeI,
     const label pointI,
     const pointEdgePoint& pointInfo,
-    const scalar tol
+    const scalar tol,
+    TrackingData& td
 )
 {
-    const pointField& points = mesh.points();
-
     const edge& e = mesh.edges()[edgeI];
+    return update(e.centre(mesh.points()), pointInfo, tol, td);
+}
 
-    const point edgeMid(0.5*(points[e[0]] + points[e[1]]));
 
-    return
+template <class TrackingData>
-        update
+inline bool Foam::pointEdgePoint::equal
-        (
+(
-            edgeMid,
+    const pointEdgePoint& rhs,
-            pointInfo,
+    TrackingData& td
-            tol
+) const
-        );
+{
+    return operator==(rhs);
 }
 
 

src/postProcessing/functionObjects/field/streamLine/streamLineParticle.C

@@ -240,7 +240,13 @@ bool Foam::streamLineParticle::move
                 lifeTime_ = 0;
             }
 
-            if (!td.keepParticle || td.switchProcessor || lifeTime_ == 0)
+            if
+            (
+                face() != -1
+            || !td.keepParticle
+            ||  td.switchProcessor
+            ||  lifeTime_ == 0
+            )
             {
                 break;
             }

src/turbulenceModels/LES/LESdeltas/Make/files

@@ -2,6 +2,6 @@ LESdelta/LESdelta.C
 cubeRootVolDelta/cubeRootVolDelta.C
 PrandtlDelta/PrandtlDelta.C
 smoothDelta/smoothDelta.C
-maxhxhyhzDelta/maxhxhyhzDelta.C
+maxDeltaxyz/maxDeltaxyz.C
 
 LIB = $(FOAM_LIBBIN)/libLESdeltas

src/turbulenceModels/LES/LESdeltas/maxDeltaxyz/maxDeltaxyz.C

@@ -23,7 +23,7 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "maxhxhyhzDelta.H"
+#include "maxDeltaxyz.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -33,12 +33,12 @@ namespace Foam
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-defineTypeNameAndDebug(maxhxhyhzDelta, 0);
+defineTypeNameAndDebug(maxDeltaxyz, 0);
-addToRunTimeSelectionTable(LESdelta, maxhxhyhzDelta, dictionary);
+addToRunTimeSelectionTable(LESdelta, maxDeltaxyz, dictionary);
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void maxhxhyhzDelta::calcDelta()
+void maxDeltaxyz::calcDelta()
 {
     label nD = mesh().nGeometricD();
 
@@ -86,7 +86,7 @@ void maxhxhyhzDelta::calcDelta()
     }
     else if (nD == 2)
     {
-        WarningIn("maxhxhyhzDelta::calcDelta()")
+        WarningIn("maxDeltaxyz::calcDelta()")
             << "Case is 2D, LES is not strictly applicable\n"
             << endl;
 
@@ -94,7 +94,7 @@ void maxhxhyhzDelta::calcDelta()
     }
     else
     {
-        FatalErrorIn("maxhxhyhzDelta::calcDelta()")
+        FatalErrorIn("maxDeltaxyz::calcDelta()")
             << "Case is not 3D or 2D, LES is not applicable"
             << exit(FatalError);
     }
@@ -103,7 +103,7 @@ void maxhxhyhzDelta::calcDelta()
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-maxhxhyhzDelta::maxhxhyhzDelta
+maxDeltaxyz::maxDeltaxyz
 (
     const word& name,
     const fvMesh& mesh,
@@ -119,14 +119,14 @@ maxhxhyhzDelta::maxhxhyhzDelta
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void maxhxhyhzDelta::read(const dictionary& dd)
+void maxDeltaxyz::read(const dictionary& dd)
 {
     dd.subDict(type() + "Coeffs").lookup("deltaCoeff") >> deltaCoeff_;
     calcDelta();
 }
 
 
-void maxhxhyhzDelta::correct()
+void maxDeltaxyz::correct()
 {
     if (mesh_.changing())
     {

src/turbulenceModels/LES/LESdeltas/maxDeltaxyz/maxDeltaxyz.H

@@ -22,20 +22,20 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::maxhxhyhzDelta
+    Foam::maxDeltaxyz
 
 Description
-    maxhxhyhzDelta takes the maximum of the three dimensions per cell:
+    maxDeltaxyz takes the maximum of the three dimensions per cell:
     max(hx, hy, hz). Valid for structures hexahedral cells only.
 
 
 SourceFiles
-    maxhxhyhzDelta.C
+    maxDeltaxyz.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef maxhxhyhzDeltaDelta_H
+#ifndef maxDeltaxyzDelta_H
-#define maxhxhyhzDeltaDelta_H
+#define maxDeltaxyzDelta_H
 
 #include "LESdelta.H"
 
@@ -45,10 +45,10 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class maxhxhyhzDelta Declaration
+                           Class maxDeltaxyz Declaration
 \*---------------------------------------------------------------------------*/
 
-class maxhxhyhzDelta
+class maxDeltaxyz
 :
     public LESdelta
 {
@@ -60,8 +60,8 @@ class maxhxhyhzDelta
     // Private Member Functions
 
         //- Disallow default bitwise copy construct and assignment
-        maxhxhyhzDelta(const maxhxhyhzDelta&);
+        maxDeltaxyz(const maxDeltaxyz&);
-        void operator=(const maxhxhyhzDelta&);
+        void operator=(const maxDeltaxyz&);
 
         // Calculate the delta values
         void calcDelta();
@@ -70,13 +70,13 @@ class maxhxhyhzDelta
 public:
 
     //- Runtime type information
-    TypeName("maxhxhyhzDelta");
+    TypeName("maxDeltaxyz");
 
 
     // Constructors
 
         //- Construct from name, mesh and IOdictionary
-        maxhxhyhzDelta
+        maxDeltaxyz
         (
             const word& name,
             const fvMesh& mesh,
@@ -85,7 +85,7 @@ public:
 
 
     //- Destructor
-    virtual ~maxhxhyhzDelta()
+    virtual ~maxDeltaxyz()
     {}
 
 

tutorials/heatTransfer/buoyantBoussinesqPimpleFoam/hotRoom/system/fvSchemes

@@ -34,7 +34,7 @@ divSchemes
     div(phi,epsilon) Gauss upwind;
     div(phi,R)      Gauss upwind;
     div(R)          Gauss linear;
-    div((nuEff*dev(grad(U).T()))) Gauss linear;
+    div((nuEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes

tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/hotRoom/system/fvSchemes

@@ -32,7 +32,7 @@ divSchemes
     div(phi,T)      Gauss upwind;
     div(phi,k)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
-    div((nuEff*dev(grad(U).T()))) Gauss linear;
+    div((nuEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes

tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/fvSchemes

@@ -34,7 +34,7 @@ divSchemes
     div(phi,epsilon) Gauss upwind;
     div(phi,R)      Gauss upwind;
     div(R)          Gauss linear;
-    div((nuEff*dev(grad(U).T()))) Gauss linear;
+    div((nuEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes

tutorials/incompressible/MRFSimpleFoam/mixerVessel2D/system/fvSchemes

@@ -33,7 +33,7 @@ divSchemes
     div(phi,U)      Gauss limitedLinearV 1;
     div(phi,k)      Gauss limitedLinear 1;
     div(phi,epsilon) Gauss limitedLinear 1;
-    div((nuEff*dev(grad(U).T()))) Gauss linear;
+    div((nuEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes

tutorials/incompressible/adjointShapeOptimizationFoam/pitzDaily/system/fvSchemes

@@ -35,7 +35,7 @@ divSchemes
     div(phi,U)      Gauss upwind;
     div(phi,k)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
-    div((nuEff*dev(grad(U).T()))) Gauss linear;
+    div((nuEff*dev(T(grad(U))))) Gauss linear;
 
     div(-phi,Ua)    Gauss upwind;
     div((nuEff*dev(grad(Ua).T()))) Gauss linear;

tutorials/incompressible/boundaryFoam/boundaryLaunderSharma/system/fvSchemes

@@ -32,7 +32,7 @@ divSchemes
     div(phi,epsilon) Gauss linear;
     div(phi,R)      Gauss linear;
     div(phi,nuTilda) Gauss linear;
-    div((nuEff*dev(grad(U).T()))) Gauss linear;
+    div((nuEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes

tutorials/incompressible/boundaryFoam/boundaryWallFunctions/system/fvSchemes

@@ -32,7 +32,7 @@ divSchemes
     div(phi,epsilon) Gauss linear;
     div(phi,R)      Gauss linear;
     div(phi,nuTilda) Gauss linear;
-    div((nuEff*dev(grad(U).T()))) Gauss linear;
+    div((nuEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes

tutorials/incompressible/boundaryFoam/boundaryWallFunctionsProfile/system/fvSchemes

@@ -33,7 +33,7 @@ divSchemes
     div(phi,omega)  Gauss linear;
     div(phi,R)      Gauss linear;
     div(phi,nuTilda) Gauss linear;
-    div((nuEff*dev(grad(U).T()))) Gauss linear;
+    div((nuEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes

tutorials/incompressible/channelFoam/channel395/system/fvSchemes

@@ -35,7 +35,7 @@ divSchemes
     div(phi,B)      Gauss limitedLinear 1;
     div(B)          Gauss linear;
     div(phi,nuTilda) Gauss limitedLinear 1;
-    div((nuEff*dev(grad(U).T()))) Gauss linear;
+    div((nuEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes

tutorials/incompressible/pimpleDyMFoam/movingCone/system/fvSchemes

@@ -30,7 +30,7 @@ divSchemes
 {
     default         none;
     div(phi,U)      Gauss linear;
-    div((nuEff*dev(grad(U).T()))) Gauss linear;
+    div((nuEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes

tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion2D_pimpleDyMFoam/system/fvSchemes

@@ -32,7 +32,7 @@ divSchemes
     div(phi,U)      Gauss linearUpwind grad(U);
     div(phi,k)      Gauss limitedLinear 1;
     div(phi,omega)  Gauss limitedLinear 1;
-    div((nuEff*dev(grad(U).T()))) Gauss linear;
+    div((nuEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes

tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion2D_simpleFoam/system/fvSchemes

@@ -32,7 +32,7 @@ divSchemes
     div(phi,U)      Gauss linearUpwind grad(U);
     div(phi,k)      Gauss upwind;
     div(phi,omega)  Gauss upwind;
-    div((nuEff*dev(grad(U).T()))) Gauss linear;
+    div((nuEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes

tutorials/incompressible/pimpleFoam/t-junction-with-fan/system/fvSchemes

@@ -36,7 +36,7 @@ divSchemes
     div(phi,R)      Gauss limitedLinear 1;
     div(R)          Gauss linear;
     div(phi,nuTilda) Gauss limitedLinear 1;
-    div((nuEff*dev(grad(U).T()))) Gauss linear;
+    div((nuEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes

tutorials/incompressible/pimpleFoam/t-junction/system/fvSchemes

@@ -36,7 +36,7 @@ divSchemes
     div(phi,R)      Gauss limitedLinear 1;
     div(R)          Gauss linear;
     div(phi,nuTilda) Gauss limitedLinear 1;
-    div((nuEff*dev(grad(U).T()))) Gauss linear;
+    div((nuEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes

tutorials/incompressible/pisoFoam/les/pitzDaily/system/fvSchemes

@@ -35,7 +35,7 @@ divSchemes
     div(phi,B)      Gauss limitedLinear 1;
     div(phi,nuTilda) Gauss limitedLinear 1;
     div(B)          Gauss linear;
-    div((nuEff*dev(grad(U).T()))) Gauss linear;
+    div((nuEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes

tutorials/incompressible/pisoFoam/les/pitzDailyDirectMapped/system/fvSchemes

@@ -35,7 +35,7 @@ divSchemes
     div(phi,B)      Gauss limitedLinear 1;
     div(phi,nuTilda) Gauss limitedLinear 1;
     div(B)          Gauss linear;
-    div((nuEff*dev(grad(U).T()))) Gauss linear;
+    div((nuEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes

tutorials/incompressible/pisoFoam/ras/cavity/system/fvSchemes

@@ -36,7 +36,7 @@ divSchemes
     div(phi,R)      Gauss limitedLinear 1;
     div(R)          Gauss linear;
     div(phi,nuTilda) Gauss limitedLinear 1;
-    div((nuEff*dev(grad(U).T()))) Gauss linear;
+    div((nuEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes

tutorials/incompressible/porousSimpleFoam/angledDuctExplicit/system/fvSchemes

@@ -30,7 +30,7 @@ gradSchemes
 divSchemes
 {
     div(phi,U)      Gauss upwind;
-    div((nuEff*dev(grad(U).T()))) Gauss linear;
+    div((nuEff*dev(T(grad(U))))) Gauss linear;
     div(phi,epsilon) Gauss upwind;
     div(phi,k)      Gauss upwind;
 }

tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/system/fvSchemes

@@ -30,7 +30,7 @@ gradSchemes
 divSchemes
 {
     div(phi,U)      Gauss upwind;
-    div((nuEff*dev(grad(U).T()))) Gauss linear;
+    div((nuEff*dev(T(grad(U))))) Gauss linear;
     div(phi,epsilon) Gauss upwind;
     div(phi,k)      Gauss upwind;
 }

tutorials/incompressible/simpleFoam/airFoil2D/system/fvSchemes

@@ -32,7 +32,7 @@ divSchemes
     default         none;
     div(phi,U)      Gauss linearUpwind Gauss linear;
     div(phi,nuTilda) Gauss linearUpwind Gauss linear;
-    div((nuEff*dev(grad(U).T()))) Gauss linear;
+    div((nuEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes

tutorials/incompressible/simpleFoam/motorBike/system/fvSchemes

@@ -32,7 +32,7 @@ divSchemes
     div(phi,U)      Gauss linearUpwind grad(U);
     div(phi,k)      Gauss upwind;
     div(phi,omega)  Gauss upwind;
-    div((nuEff*dev(grad(U).T()))) Gauss linear;
+    div((nuEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes

tutorials/incompressible/simpleFoam/pitzDailyExptInlet/system/fvSchemes

@@ -36,7 +36,7 @@ divSchemes
     div(phi,R)      Gauss upwind;
     div(R)          Gauss linear;
     div(phi,nuTilda) Gauss upwind;
-    div((nuEff*dev(grad(U).T()))) Gauss linear;
+    div((nuEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes

tutorials/incompressible/simpleFoam/windTurbineTerrain/0/U

@@ -46,13 +46,15 @@ boundaryField
 
     "terrain_.*"
     {
-        type            fixedValue;
+        type            uniformFixedValue;
+        uniformValue    (0 0 0);
         value           uniform (0 0 0);
     }
 
     ground
     {
-        type            fixedValue;
+        type            uniformFixedValue;
+        uniformValue    (0 0 0);
         value           uniform (0 0 0);
     }
 

tutorials/incompressible/simpleFoam/windTurbineTerrain/0/epsilon

@@ -23,7 +23,7 @@ internalField   uniform $turbulentEpsilon;
 
 boundaryField
 {
-     #include "include/ABLConditions"
+    #include "include/ABLConditions"
 
     "terrain_.*"
     {

tutorials/incompressible/simpleFoam/windTurbineTerrain/0/k

@@ -32,8 +32,9 @@ boundaryField
     }
     inlet
     {
-        type            fixedValue;
+        type            uniformFixedValue;
-        value           uniform $turbulentKE;
+        uniformValue    $turbulentKE;
+        value           $turbulentKE;
     }
     "terrain_.*"
     {

tutorials/incompressible/simpleFoam/windTurbineTerrain/0/p

@@ -29,8 +29,9 @@ boundaryField
 
     outlet
     {
-        type            fixedValue;
+        type            uniformFixedValue;
-        value           $internalField;
+        value           uniform $pressure;
+        uniformValue    $pressure;
     }
 
     "terrain_.*"

tutorials/incompressible/simpleFoam/windTurbineTerrain/Allclean

@@ -19,4 +19,7 @@ rm -rf constant/polyMesh/sets
 #rm -rf constant/refinementHistory
 #rm -rf constant/surfaceIndex
 
+# Reset decomposeParDict
+cp system/decomposeParDict-nonPar system/decomposeParDict
+
 # ----------------------------------------------------------------- end-of-file

tutorials/incompressible/simpleFoam/windTurbineTerrain/Allrun

@@ -4,10 +4,29 @@ cd ${0%/*} || exit 1    # run from this directory
 # Source tutorial run functions
 . $WM_PROJECT_DIR/bin/tools/RunFunctions
 
+compileApplication windSimpleFoam
+
 runApplication blockMesh
-runApplication snappyHexMesh -overwrite
+cp system/decomposeParDict-nonPar system/decomposeParDict
-runApplication setSet -batch makeZones
+runApplication decomposePar
-runApplication setsToZones -noFlipMap
+
-runApplication windSimpleFoam
+#runApplication snappyHexMesh -overwrite
+#runApplication setSet -batch makeZones
+#runApplication setsToZones -noFlipMap
+#runApplication windSimpleFoam
+
+cp system/decomposeParDict-par system/decomposeParDict
+runParallel snappyHexMesh 2 -overwrite
+
+# Add wildcard entries for meshes patches since not preserved
+# by decomposePar. Notice -literalRE option to add wildcard itself
+# without evaluation.
+runParallel changeDictionary 2 -literalRE
+
+runParallel setSet 2 -batch makeZones
+runParallel windSimpleFoam 2
+
+runApplication reconstructParMesh -constant
+runApplication reconstructPar
 
 # ----------------------------------------------------------------- end-of-file

tutorials/incompressible/simpleFoam/windTurbineTerrain/makeZones

@@ -1,2 +1,4 @@
 cellSet actuationDisk1 new boxToCell (581850.5 4785805 1061) (581850.8 4785815 1071)
-cellSet actuationDisk2 new boxToCell (581754 4785658 1065) (581754.4 4785668 1075)
+cellZoneSet actuationDisk1 new setToCellZone actuationDisk1
+cellSet actuationDisk2 new boxToCell (581754 4785658 1065) (581754.4 4785668 1075)
+cellZoneSet actuationDisk2 new setToCellZone actuationDisk2

tutorials/incompressible/simpleFoam/windTurbineTerrain/system/changeDictionaryDict

@@ -0,0 +1,201 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                     |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      changeDictionaryDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Use absolute path to make sure it also works in parallel
+#include        "$FOAM_CASE/0/include/initialConditions"
+
+dictionaryReplacement
+{
+
+    // Specify 
+    // - all fvPatchFields with potential non-uniform values
+    // - all fvPatchFields originating from meshing
+    // - all fvPatchFields originating from mesh-redistribution
+
+    p
+    {
+        boundaryField
+        {
+            outlet
+            {
+                type            uniformFixedValue;
+                value           $pressure;
+                uniformValue    $pressure;
+            }
+            inlet
+            {
+                type            zeroGradient;
+            }
+            "terrain_.*"
+            {
+                type            zeroGradient;
+            }
+            ground
+            {
+                type            zeroGradient;
+            }
+            #include "$FOAM_CASE/0/include/sideAndTopPatches"
+            "procBoundary.*"
+            {
+                type            processor;
+            }
+        }
+    }
+
+    k
+    {
+        boundaryField
+        {
+            outlet
+            {
+                type            inletOutlet;
+                inletValue      uniform 0.0;
+                value           uniform $turbulentKE;
+            }
+            inlet
+            {
+                type            uniformFixedValue;
+                uniformValue    $turbulentKE;
+            }
+            "terrain_.*"
+            {
+                type            kqRWallFunction;
+                value           uniform 0.0;
+            }
+            ground
+            {
+                type            zeroGradient;
+            }
+            #include "$FOAM_CASE/0/include/sideAndTopPatches"
+            "procBoundary.*"
+            {
+                type            processor;
+            }
+        }
+    }
+
+    U
+    {
+        boundaryField
+        {
+            outlet
+            {
+                type            inletOutlet;
+                inletValue      uniform (0 0 0);
+                value           uniform $flowVelocity;
+            }
+            inlet
+            {
+                type            atmBoundaryLayerInletVelocity;
+                Uref            10.0;
+                Href            20;
+                n               (1 0 0);
+                z               (0 0 1);
+                z0              0.1;
+                zGround         935.0;
+                value           uniform $flowVelocity;
+            }
+            "terrain_.*"
+            {
+                type            uniformFixedValue;
+                uniformValue    $flowVelocity;
+            }
+            ground
+            {
+                type            uniformFixedValue;
+                uniformValue    $flowVelocity;
+            }
+            #include "$FOAM_CASE/0/include/sideAndTopPatches"
+            "procBoundary.*"
+            {
+                type            processor;
+            }
+        }
+    }
+
+    nut
+    {
+        boundaryField
+        {
+            outlet
+            {
+                type            calculated;
+                value           uniform 0;
+            }
+            inlet
+            {
+                type            calculated;
+                value           uniform 0;
+            }
+            "terrain_.*"
+            {
+                type            nutkRoughWallFunction;
+                Ks              uniform 0.2; //Ks = 20 Z0
+                Cs              uniform 0.5;
+                value           uniform 0.0;
+            }
+            ground
+            {
+                type            calculated;
+                value           uniform 0;
+            }
+            #include "$FOAM_CASE/0/include/sideAndTopPatches"
+            "procBoundary.*"
+            {
+                type            processor;
+            }
+        }
+    }
+
+    epsilon
+    {
+        boundaryField
+        {
+            outlet
+            {
+                type            zeroGradient;
+            }
+            inlet
+            {
+                type            atmBoundaryLayerInletEpsilon;
+                Ustar           0.82;
+                z               (0 0 1);
+                z0              0.1;
+                value           uniform $turbulentEpsilon;
+                zGround         935.0;
+            }
+            "terrain_.*"
+            {
+                type            epsilonWallFunction;
+                Cmu             0.09;
+                kappa           0.4;
+                E               9.8;
+                value           uniform $turbulentEpsilon;
+            }
+            ground
+            {
+                type            zeroGradient;
+            }
+            #include "$FOAM_CASE/0/include/sideAndTopPatches"
+            "procBoundary.*"
+            {
+                type            processor;
+            }
+        }
+    }
+}
+
+// ************************************************************************* //

tutorials/incompressible/simpleFoam/windTurbineTerrain/system/decomposeParDict

@@ -15,35 +15,15 @@ FoamFile
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-numberOfSubdomains 1;
+numberOfSubdomains 2;
 
 method          hierarchical;
-// method          metis;
-// method          ptscotch;
-
-simpleCoeffs
-{
-    n               (4 1 1);
-    delta           0.001;
-}
 
 hierarchicalCoeffs
 {
-    n               (1 1 1);
+    n           (2 1 1);
-    delta           0.001;
+    delta       0.001;
-    order           xyz;
+    order       xyz;
 }
 
-manualCoeffs
-{
-    dataFile        "cellDecomposition";
-}
-
-metisCoeffs
-{
-    //n                   (5 1 1);
-    //cellWeightsFile     "constant/cellWeightsFile";
-}
-
-
 // ************************************************************************* //

tutorials/incompressible/simpleFoam/windTurbineTerrain/system/decomposeParDict-nonPar

@@ -0,0 +1,29 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      decomposeParDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains 2;
+
+method          hierarchical;
+
+hierarchicalCoeffs
+{
+    n           (2 1 1);
+    delta       0.001;
+    order       xyz;
+}
+
+// ************************************************************************* //

tutorials/incompressible/simpleFoam/windTurbineTerrain/system/decomposeParDict-par

@@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      decomposeParDict;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains 2;
+
+method          ptscotch;
+
+// ************************************************************************* //

tutorials/incompressible/simpleFoam/windTurbineTerrain/system/fvSchemes

@@ -30,7 +30,7 @@ divSchemes
 {
     default             none;
     div(phi,U)          Gauss upwind grad(U);
-    div((nuEff*dev(grad(U).T())))   Gauss linear;
+    div((nuEff*dev(T(grad(U)))))    Gauss linear;
     div(phi,epsilon)    Gauss upwind;
     div(phi,k)          Gauss upwind;
 }

tutorials/incompressible/simpleFoam/windTurbineTerrain/system/fvSolution

@@ -18,49 +18,49 @@ solvers
 {
     p
     {
-        solver      smoothSolver;
+        solver           GAMG;
-        smoother    GaussSeidel;
+        tolerance        1e-7;
-        tolerance   1e-6;
+        relTol           0.1;
-        relTol      0.01;
+        smoother         GaussSeidel;
-        maxIter     300;
+        nPreSweeps       0;
-    };
+        nPostSweeps      2;
+        cacheAgglomeration on;
+        agglomerator     faceAreaPair;
+        nCellsInCoarsestLevel 10;
+        mergeLevels      1;
+    }
 
     U
     {
-        solver      smoothSolver;
+        solver           smoothSolver;
-        smoother    GaussSeidel;
+        smoother         GaussSeidel;
-        tolerance   1e-6;
+        tolerance        1e-8;
-        relTol      0.01;
+        relTol           0.1;
-    };
+        nSweeps          1;
+    }
 
     k
     {
-        solver      smoothSolver;
+        solver           smoothSolver;
-        smoother    GaussSeidel;
+        smoother         GaussSeidel;
-        tolerance   1e-6;
+        tolerance        1e-8;
-        relTol      0.01;
+        relTol           0.1;
-    };
+        nSweeps          1;
+    }
 
     epsilon
     {
-        solver      smoothSolver;
+        solver           smoothSolver;
-        smoother    GaussSeidel;
+        smoother         GaussSeidel;
-        tolerance   1e-6;
+        tolerance        1e-8;
-        relTol      0.01;
+        relTol           0.1;
-    };
+        nSweeps          1;
-
+    }
-    omega
-    {
-        solver      smoothSolver;
-        smoother    GaussSeidel;
-        tolerance   1e-6;
-        relTol      0.1;
-    };
 }
 
 SIMPLE
 {
-    nNonOrthogonalCorrectors     1;
+    nNonOrthogonalCorrectors     0;
     convergence                  1e-3;
 }
 

tutorials/multiphase/compressibleInterFoam/les/depthCharge2D/system/fvSchemes

@@ -32,7 +32,7 @@ divSchemes
     div(phirb,alpha) Gauss interfaceCompression 1;
     div(phi,p_rgh)  Gauss upwind;
     div(phi,k)      Gauss vanLeer;
-    div((nuEff*dev(grad(U).T()))) Gauss linear;
+    div((nuEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes

tutorials/multiphase/compressibleInterFoam/les/depthCharge3D/system/fvSchemes

@@ -32,7 +32,7 @@ divSchemes
     div(phirb,alpha) Gauss interfaceCompression 1;
     div(phi,p_rgh)  Gauss upwind;
     div(phi,k)      Gauss vanLeer;
-    div((nuEff*dev(grad(U).T()))) Gauss linear;
+    div((nuEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes

tutorials/multiphase/interFoam/les/nozzleFlow2D/system/fvSchemes

@@ -34,7 +34,7 @@ divSchemes
     div(phi,B)      Gauss limitedLinear 1;
     div(B)          Gauss linear;
     div(phi,nuTilda) Gauss limitedLinear 1;
-    div((nuEff*dev(grad(U).T()))) Gauss linear;
+    div((nuEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes

tutorials/multiphase/interFoam/ras/damBreak/system/fvSchemes

@@ -35,7 +35,7 @@ divSchemes
     div(phi,R)      Gauss upwind;
     div(R)          Gauss linear;
     div(phi,nuTilda) Gauss upwind;
-    div((nuEff*dev(grad(U).T()))) Gauss linear;
+    div((nuEff*dev(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes