Merge remote branch 'OpenCFD/master' into olesenm

Conflicts:
	applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C
	applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/createShellMesh.C
	applications/utilities/surface/surfaceCheck/surfaceCheck.C
	src/finiteVolume/fields/fvPatchFields/derived/advective/advectiveFvPatchField.C
	src/finiteVolume/fields/fvPatchFields/derived/waveTransmissive/waveTransmissiveFvPatchField.C
	src/meshTools/directMapped/directMappedPolyPatch/directMappedPatchBase.C

NOTE: also needed to strip trailing space/lines in various files

src/ODE/ODESolvers/KRR4/KRR4.C

@@ -58,17 +58,17 @@ const scalar
 Foam::KRR4::KRR4(const ODE& ode)
 :
     ODESolver(ode),
-    yTemp_(n_),
-    dydxTemp_(n_),
-    g1_(n_),
-    g2_(n_),
-    g3_(n_),
-    g4_(n_),
-    yErr_(n_),
-    dfdx_(n_),
-    dfdy_(n_, n_),
-    a_(n_, n_),
-    pivotIndices_(n_)
+    yTemp_(n_, 0.0),
+    dydxTemp_(n_, 0.0),
+    g1_(n_, 0.0),
+    g2_(n_, 0.0),
+    g3_(n_, 0.0),
+    g4_(n_, 0.0),
+    yErr_(n_, 0.0),
+    dfdx_(n_, 0.0),
+    dfdy_(n_, n_, 0.0),
+    a_(n_, n_, 0.0),
+    pivotIndices_(n_, 0.0)
 {}
 
 

src/ODE/ODESolvers/RK/RK.C

@@ -57,14 +57,14 @@ const scalar
 Foam::RK::RK(const ODE& ode)
 :
     ODESolver(ode),
-    yTemp_(n_),
-    ak2_(n_),
-    ak3_(n_),
-    ak4_(n_),
-    ak5_(n_),
-    ak6_(n_),
-    yErr_(n_),
-    yTemp2_(n_)
+    yTemp_(n_, 0.0),
+    ak2_(n_, 0.0),
+    ak3_(n_, 0.0),
+    ak4_(n_, 0.0),
+    ak5_(n_, 0.0),
+    ak6_(n_, 0.0),
+    yErr_(n_, 0.0),
+    yTemp2_(n_, 0.0)
 {}
 
 

src/ODE/ODESolvers/SIBS/SIBS.C

@@ -50,17 +50,17 @@ namespace Foam
 Foam::SIBS::SIBS(const ODE& ode)
 :
     ODESolver(ode),
-    a_(iMaxX_),
-    alpha_(kMaxX_, kMaxX_),
-    d_p_(n_, kMaxX_),
-    x_p_(kMaxX_),
-    err_(kMaxX_),
+    a_(iMaxX_, 0.0),
+    alpha_(kMaxX_, kMaxX_, 0.0),
+    d_p_(n_, kMaxX_, 0.0),
+    x_p_(kMaxX_, 0.0),
+    err_(kMaxX_, 0.0),
 
-    yTemp_(n_),
-    ySeq_(n_),
-    yErr_(n_),
-    dfdx_(n_),
-    dfdy_(n_, n_),
+    yTemp_(n_, 0.0),
+    ySeq_(n_, 0.0),
+    yErr_(n_, 0.0),
+    dfdx_(n_, 0.0),
+    dfdy_(n_, n_, 0.0),
     first_(1),
     epsOld_(-1.0)
 {}

src/OpenFOAM/algorithms/MeshWave/FaceCellWave.C

@@ -514,8 +514,8 @@ void Foam::FaceCellWave<Type, TrackingData>::handleProcPatches()
                 refCast<const processorPolyPatch>(patch);
 
             // Allocate buffers
-            labelList receiveFaces(patch.size());
-            List<Type> receiveFacesInfo(patch.size());
+            labelList receiveFaces;
+            List<Type> receiveFacesInfo;
 
             {
                 UIPstream fromNeighbour(procPatch.neighbProcNo(), pBufs);
@@ -673,8 +673,7 @@ Foam::FaceCellWave<Type, TrackingData>::FaceCellWave
     hasCyclicPatches_(hasCyclicPatch()),
     nEvals_(0),
     nUnvisitedCells_(mesh_.nCells()),
-    nUnvisitedFaces_(mesh_.nFaces()),
-    iter_(0)
+    nUnvisitedFaces_(mesh_.nFaces())
 {}
 
 
@@ -705,16 +704,15 @@ Foam::FaceCellWave<Type, TrackingData>::FaceCellWave
     hasCyclicPatches_(hasCyclicPatch()),
     nEvals_(0),
     nUnvisitedCells_(mesh_.nCells()),
-    nUnvisitedFaces_(mesh_.nFaces()),
-    iter_(0)
+    nUnvisitedFaces_(mesh_.nFaces())
 {
     // 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
         (
@@ -794,7 +792,7 @@ Foam::label Foam::FaceCellWave<Type, TrackingData>::faceToCell()
             );
         }
 
-        // Neighbour. Hack for check if face has neighbour.
+        // Neighbour.
         if (faceI < nInternalFaces)
         {
             cellI = neighbour[faceI];
@@ -925,11 +923,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;
@@ -961,10 +961,10 @@ 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
 
@@ -340,8 +337,8 @@ public:
             //  counted double)
             label cellToFace();
 
-            //- Iterate until no changes or maxIter reached. Returns number of
-            //  unset cells (see getUnsetCells)
+            //- Iterate until no changes or maxIter reached.  Returns actual
+            //  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
-        //  unset cells (see getUnsetCells)
+        //- Iterate until no changes or maxIter reached. Returns actual
+        //  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/fields/GeometricFields/GeometricTensorField/GeometricTensorField.C

@@ -36,6 +36,7 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+UNARY_FUNCTION(tensor, tensor, T, transform)
 UNARY_FUNCTION(scalar, tensor, tr, transform)
 UNARY_FUNCTION(sphericalTensor, tensor, sph, transform)
 UNARY_FUNCTION(symmTensor, tensor, symm, transform)

src/OpenFOAM/fields/GeometricFields/GeometricTensorField/GeometricTensorField.H

@@ -49,6 +49,7 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+UNARY_FUNCTION(tensor, tensor, T, transform)
 UNARY_FUNCTION(scalar, tensor, tr, transform)
 UNARY_FUNCTION(sphericalTensor, tensor, sph, transform)
 UNARY_FUNCTION(symmTensor, tensor, symm, transform)

src/OpenFOAM/global/argList/argList.C

@@ -736,7 +736,7 @@ Foam::argList::argList
         {
             Info<< "Slaves : " << slaveProcs << nl;
             if (roots.size())
-            {                
+            {
                 Info<< "Roots  : " << roots << nl;
             }
             Info<< "Pstream initialized with:" << nl
@@ -770,7 +770,8 @@ Foam::argList::argList
 
     if (bannerEnabled)
     {
-        Info<< "Monitoring run-time modified files using "
+        Info<< "fileModificationChecking : "
+            << "Monitoring run-time modified files using "
             << regIOobject::fileCheckTypesNames
                 [
                     regIOobject::fileModificationChecking

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

@@ -145,7 +145,7 @@ Foam::boundBox::boundBox
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-Foam::tmp<Foam::pointField> Foam::boundBox::corners() const
+Foam::tmp<Foam::pointField> Foam::boundBox::points() const
 {
     tmp<pointField> tPts = tmp<pointField>(new pointField(8));
     pointField& pt = tPts();
@@ -274,7 +274,7 @@ Foam::Istream& Foam::operator>>(Istream& is, boundBox& bb)
 {
     if (is.format() == IOstream::ASCII)
     {
-        return is >> bb.min_ >> bb.max_;
+        is >> bb.min_ >> bb.max_;
     }
     else
     {

src/OpenFOAM/meshes/boundBox/boundBox.H

@@ -158,7 +158,7 @@ public:
             inline scalar avgDim() const;
 
             //- Return corner points in an order corresponding to a 'hex' cell
-            tmp<pointField> corners() const;
+            tmp<pointField> points() const;
 
         // Query
 

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

@@ -165,6 +165,35 @@ const Foam::List<Foam::labelPair>& Foam::mapDistribute::schedule() const
 }
 
 
+void Foam::mapDistribute::checkReceivedSize
+(
+    const label procI,
+    const label expectedSize,
+    const label receivedSize
+)
+{
+    if (receivedSize != expectedSize)
+    {
+        FatalErrorIn
+        (
+            "template<class T>\n"
+            "void mapDistribute::distribute\n"
+            "(\n"
+            "    const Pstream::commsTypes commsType,\n"
+            "    const List<labelPair>& schedule,\n"
+            "    const label constructSize,\n"
+            "    const labelListList& subMap,\n"
+            "    const labelListList& constructMap,\n"
+            "    List<T>& field\n"
+            ")\n"
+        )   << "Expected from processor " << procI
+            << " " << expectedSize << " but received "
+            << receivedSize << " elements."
+            << abort(FatalError);
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 //- Construct from components

src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H

@@ -93,6 +93,15 @@ class mapDistribute
         mutable autoPtr<List<labelPair> > schedulePtr_;
 
 
+   // Private Member Functions
+
+        static void checkReceivedSize
+        (
+            const label procI,
+            const label expectedSize,
+            const label receivedSize
+        );
+
 public:
 
     // Constructors

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

@@ -111,25 +111,7 @@ void Foam::mapDistribute::distribute
                 IPstream fromNbr(Pstream::blocking, domain);
                 List<T> subField(fromNbr);
 
-                if (subField.size() != map.size())
-                {
-                    FatalErrorIn
-                    (
-                        "template<class T>\n"
-                        "void mapDistribute::distribute\n"
-                        "(\n"
-                        "    const Pstream::commsTypes commsType,\n"
-                        "    const List<labelPair>& schedule,\n"
-                        "    const label constructSize,\n"
-                        "    const labelListList& subMap,\n"
-                        "    const labelListList& constructMap,\n"
-                        "    List<T>& field\n"
-                        ")\n"
-                    )   << "Expected from processor " << domain
-                        << " " << map.size() << " but received "
-                        << subField.size() << " elements."
-                        << abort(FatalError);
-                }
+                checkReceivedSize(domain, map.size(), subField.size());
 
                 forAll(map, i)
                 {
@@ -160,46 +142,52 @@ void Foam::mapDistribute::distribute
         forAll(schedule, i)
         {
             const labelPair& twoProcs = schedule[i];
+            // twoProcs is a swap pair of processors. The first one is the
+            // one that needs to send first and then receive.
+
             label sendProc = twoProcs[0];
             label recvProc = twoProcs[1];
 
             if (Pstream::myProcNo() == sendProc)
             {
-                // I am sender. Send to recvProc.
-                OPstream toNbr(Pstream::scheduled, recvProc);
-                toNbr << UIndirectList<T>(field, subMap[recvProc]);
+                // I am send first, receive next
+                {
+                    OPstream toNbr(Pstream::scheduled, recvProc);
+                    toNbr << UIndirectList<T>(field, subMap[recvProc]);
+                }
+                {
+                    IPstream fromNbr(Pstream::scheduled, recvProc);
+                    List<T> subField(fromNbr);
+
+                    const labelList& map = constructMap[recvProc];
+
+                    checkReceivedSize(recvProc, map.size(), subField.size());
+
+                    forAll(map, i)
+                    {
+                        newField[map[i]] = subField[i];
+                    }
+                }
             }
             else
             {
-                // I am receiver. Receive from sendProc.
-                IPstream fromNbr(Pstream::scheduled, sendProc);
-                List<T> subField(fromNbr);
-
-                const labelList& map = constructMap[sendProc];
-
-                if (subField.size() != map.size())
+                // I am receive first, send next
                 {
-                    FatalErrorIn
-                    (
-                        "template<class T>\n"
-                        "void mapDistribute::distribute\n"
-                        "(\n"
-                        "    const Pstream::commsTypes commsType,\n"
-                        "    const List<labelPair>& schedule,\n"
-                        "    const label constructSize,\n"
-                        "    const labelListList& subMap,\n"
-                        "    const labelListList& constructMap,\n"
-                        "    List<T>& field\n"
-                        ")\n"
-                    )   << "Expected from processor " << sendProc
-                        << " " << map.size() << " but received "
-                        << subField.size() << " elements."
-                        << abort(FatalError);
+                    IPstream fromNbr(Pstream::scheduled, sendProc);
+                    List<T> subField(fromNbr);
+
+                    const labelList& map = constructMap[sendProc];
+
+                    checkReceivedSize(sendProc, map.size(), subField.size());
+
+                    forAll(map, i)
+                    {
+                        newField[map[i]] = subField[i];
+                    }
                 }
-
-                forAll(map, i)
                 {
-                    newField[map[i]] = subField[i];
+                    OPstream toNbr(Pstream::scheduled, sendProc);
+                    toNbr << UIndirectList<T>(field, subMap[sendProc]);
                 }
             }
         }
@@ -258,25 +246,7 @@ void Foam::mapDistribute::distribute
                     UIPstream str(domain, pBufs);
                     List<T> recvField(str);
 
-                    if (recvField.size() != map.size())
-                    {
-                        FatalErrorIn
-                        (
-                            "template<class T>\n"
-                            "void mapDistribute::distribute\n"
-                            "(\n"
-                            "    const Pstream::commsTypes commsType,\n"
-                            "    const List<labelPair>& schedule,\n"
-                            "    const label constructSize,\n"
-                            "    const labelListList& subMap,\n"
-                            "    const labelListList& constructMap,\n"
-                            "    List<T>& field\n"
-                            ")\n"
-                        )   << "Expected from processor " << domain
-                            << " " << map.size() << " but received "
-                            << recvField.size() << " elements."
-                            << abort(FatalError);
-                    }
+                    checkReceivedSize(domain, map.size(), recvField.size());
 
                     forAll(map, i)
                     {
@@ -379,29 +349,13 @@ void Foam::mapDistribute::distribute
 
                 if (domain != Pstream::myProcNo() && map.size())
                 {
-                    if (recvFields[domain].size() != map.size())
-                    {
-                        FatalErrorIn
-                        (
-                            "template<class T>\n"
-                            "void mapDistribute::distribute\n"
-                            "(\n"
-                            "    const Pstream::commsTypes commsType,\n"
-                            "    const List<labelPair>& schedule,\n"
-                            "    const label constructSize,\n"
-                            "    const labelListList& subMap,\n"
-                            "    const labelListList& constructMap,\n"
-                            "    List<T>& field\n"
-                            ")\n"
-                        )   << "Expected from processor " << domain
-                            << " " << map.size() << " but received "
-                            << recvFields[domain].size() << " elements."
-                            << abort(FatalError);
-                    }
+                    const List<T>& subField = recvFields[domain];
+
+                    checkReceivedSize(domain, map.size(), subField.size());
 
                     forAll(map, i)
                     {
-                        field[map[i]] = recvFields[domain][i];
+                        field[map[i]] = subField[i];
                     }
                 }
             }
@@ -502,25 +456,7 @@ void Foam::mapDistribute::distribute
                 IPstream fromNbr(Pstream::blocking, domain);
                 List<T> subField(fromNbr);
 
-                if (subField.size() != map.size())
-                {
-                    FatalErrorIn
-                    (
-                        "template<class T>\n"
-                        "void mapDistribute::distribute\n"
-                        "(\n"
-                        "    const Pstream::commsTypes commsType,\n"
-                        "    const List<labelPair>& schedule,\n"
-                        "    const label constructSize,\n"
-                        "    const labelListList& subMap,\n"
-                        "    const labelListList& constructMap,\n"
-                        "    List<T>& field\n"
-                        ")\n"
-                    )   << "Expected from processor " << domain
-                        << " " << map.size() << " but received "
-                        << subField.size() << " elements."
-                        << abort(FatalError);
-                }
+                checkReceivedSize(domain, map.size(), subField.size());
 
                 forAll(map, i)
                 {
@@ -551,46 +487,50 @@ void Foam::mapDistribute::distribute
         forAll(schedule, i)
         {
             const labelPair& twoProcs = schedule[i];
+            // twoProcs is a swap pair of processors. The first one is the
+            // one that needs to send first and then receive.
+
             label sendProc = twoProcs[0];
             label recvProc = twoProcs[1];
 
             if (Pstream::myProcNo() == sendProc)
             {
-                // I am sender. Send to recvProc.
-                OPstream toNbr(Pstream::scheduled, recvProc);
-                toNbr << UIndirectList<T>(field, subMap[recvProc]);
+                // I am send first, receive next
+                {
+                    OPstream toNbr(Pstream::scheduled, recvProc);
+                    toNbr << UIndirectList<T>(field, subMap[recvProc]);
+                }
+                {
+                    IPstream fromNbr(Pstream::scheduled, recvProc);
+                    List<T> subField(fromNbr);
+                    const labelList& map = constructMap[recvProc];
+
+                    checkReceivedSize(recvProc, map.size(), subField.size());
+
+                    forAll(map, i)
+                    {
+                        cop(newField[map[i]], subField[i]);
+                    }
+                }
             }
             else
             {
-                // I am receiver. Receive from sendProc.
-                IPstream fromNbr(Pstream::scheduled, sendProc);
-                List<T> subField(fromNbr);
-
-                const labelList& map = constructMap[sendProc];
-
-                if (subField.size() != map.size())
+                // I am receive first, send next
                 {
-                    FatalErrorIn
-                    (
-                        "template<class T>\n"
-                        "void mapDistribute::distribute\n"
-                        "(\n"
-                        "    const Pstream::commsTypes commsType,\n"
-                        "    const List<labelPair>& schedule,\n"
-                        "    const label constructSize,\n"
-                        "    const labelListList& subMap,\n"
-                        "    const labelListList& constructMap,\n"
-                        "    List<T>& field\n"
-                        ")\n"
-                    )   << "Expected from processor " << sendProc
-                        << " " << map.size() << " but received "
-                        << subField.size() << " elements."
-                        << abort(FatalError);
+                    IPstream fromNbr(Pstream::scheduled, sendProc);
+                    List<T> subField(fromNbr);
+                    const labelList& map = constructMap[sendProc];
+
+                    checkReceivedSize(sendProc, map.size(), subField.size());
+
+                    forAll(map, i)
+                    {
+                        cop(newField[map[i]], subField[i]);
+                    }
                 }
-
-                forAll(map, i)
                 {
-                    cop(newField[map[i]], subField[i]);
+                    OPstream toNbr(Pstream::scheduled, sendProc);
+                    toNbr << UIndirectList<T>(field, subMap[sendProc]);
                 }
             }
         }
@@ -649,25 +589,7 @@ void Foam::mapDistribute::distribute
                     UIPstream str(domain, pBufs);
                     List<T> recvField(str);
 
-                    if (recvField.size() != map.size())
-                    {
-                        FatalErrorIn
-                        (
-                            "template<class T>\n"
-                            "void mapDistribute::distribute\n"
-                            "(\n"
-                            "    const Pstream::commsTypes commsType,\n"
-                            "    const List<labelPair>& schedule,\n"
-                            "    const label constructSize,\n"
-                            "    const labelListList& subMap,\n"
-                            "    const labelListList& constructMap,\n"
-                            "    List<T>& field\n"
-                            ")\n"
-                        )   << "Expected from processor " << domain
-                            << " " << map.size() << " but received "
-                            << recvField.size() << " elements."
-                            << abort(FatalError);
-                    }
+                    checkReceivedSize(domain, map.size(), recvField.size());
 
                     forAll(map, i)
                     {
@@ -769,29 +691,13 @@ void Foam::mapDistribute::distribute
 
                 if (domain != Pstream::myProcNo() && map.size())
                 {
-                    if (recvFields[domain].size() != map.size())
-                    {
-                        FatalErrorIn
-                        (
-                            "template<class T>\n"
-                            "void mapDistribute::distribute\n"
-                            "(\n"
-                            "    const Pstream::commsTypes commsType,\n"
-                            "    const List<labelPair>& schedule,\n"
-                            "    const label constructSize,\n"
-                            "    const labelListList& subMap,\n"
-                            "    const labelListList& constructMap,\n"
-                            "    List<T>& field\n"
-                            ")\n"
-                        )   << "Expected from processor " << domain
-                            << " " << map.size() << " but received "
-                            << recvFields[domain].size() << " elements."
-                            << abort(FatalError);
-                    }
+                    const List<T>& subField = recvFields[domain];
+
+                    checkReceivedSize(domain, map.size(), subField.size());
 
                     forAll(map, i)
                     {
-                        cop(field[map[i]], recvFields[domain][i]);
+                        cop(field[map[i]], subField[i]);
                     }
                 }
             }

src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H

@@ -41,6 +41,7 @@ SourceFiles
 #include "tetPointRef.H"
 #include "triPointRef.H"
 #include "polyMesh.H"
+#include "triFace.H"
 #include "face.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -146,6 +147,10 @@ public:
             //  mesh face for this tet from the supplied mesh
             inline triPointRef faceTri(const polyMesh& mesh) const;
 
+            //- Return the point indices corresponding to the tri on the mesh
+            //  face for this tet from the supplied mesh
+            inline triFace faceTriIs(const polyMesh& mesh) const;
+
             //- Return the geometry corresponding to the tri on the
             //  mesh face for this tet from the supplied mesh using
             //  the old position

src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndicesI.H

@@ -122,6 +122,21 @@ Foam::triPointRef Foam::tetIndices::faceTri(const polyMesh& mesh) const
 }
 
 
+Foam::triFace Foam::tetIndices::faceTriIs(const polyMesh& mesh) const
+{
+    const faceList& pFaces = mesh.faces();
+
+    const Foam::face& f = pFaces[faceI_];
+
+    return triFace
+    (
+        f[faceBasePtI_],
+        f[facePtAI_],
+        f[facePtBI_]
+    );
+}
+
+
 Foam::triPointRef Foam::tetIndices::oldFaceTri(const polyMesh& mesh) const
 {
     const pointField& oldPPts = mesh.oldPoints();

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;
+        }
     }
 }
 
@@ -859,22 +883,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_;
 }
@@ -1014,22 +1022,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/meshes/primitiveShapes/tetrahedron/tetrahedronI.H

@@ -25,6 +25,7 @@ License
 
 #include "triangle.H"
 #include "IOstreams.H"
+#include "triPointRef.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -150,23 +151,63 @@ inline Point Foam::tetrahedron<Point, PointRef>::circumCentre() const
     vector ba = b ^ a;
     vector ca = c ^ a;
 
-    return a_ + 0.5*(a + (lambda*ba - mu*ca)/((c & ba) + ROOTVSMALL));
+    vector num = lambda*ba - mu*ca;
+    scalar denom = (c & ba);
+
+    if (Foam::mag(denom) < ROOTVSMALL)
+    {
+        WarningIn("Point tetrahedron<Point, PointRef>::circumCentre() const")
+            << "Degenerate tetrahedron:" << nl << *this << nl
+            <<"Returning centre instead of circumCentre."
+            << endl;
+
+        return centre();
+    }
+
+    return a_ + 0.5*(a + num/denom);
 }
 
 
 template<class Point, class PointRef>
 inline Foam::scalar Foam::tetrahedron<Point, PointRef>::circumRadius() const
 {
-    return Foam::mag(a_ - circumCentre());
+    vector a = b_ - a_;
+    vector b = c_ - a_;
+    vector c = d_ - a_;
+
+    scalar lambda = magSqr(c) - (a & c);
+    scalar mu = magSqr(b) - (a & b);
+
+    vector ba = b ^ a;
+    vector ca = c ^ a;
+
+    vector num = lambda*ba - mu*ca;
+    scalar denom = (c & ba);
+
+    if (Foam::mag(denom) < ROOTVSMALL)
+    {
+        WarningIn("Point tetrahedron<Point, PointRef>::circumCentre() const")
+            << "Degenerate tetrahedron:" << nl << *this << nl
+            << "Returning GREAT for circumRadius."
+            << endl;
+
+        return GREAT;
+    }
+
+    return Foam::mag(0.5*(a + num/denom));
 }
 
 
 template<class Point, class PointRef>
 inline Foam::scalar Foam::tetrahedron<Point, PointRef>::quality() const
 {
-    // Note: 8/(9*sqrt(3)) = 0.5132002393
-
-    return mag()/(0.5132002393*pow3(min(circumRadius(), GREAT)) + ROOTVSMALL);
+    return
+        mag()
+       /(
+            8.0/(9.0*sqrt(3.0))
+           *pow3(min(circumRadius(), GREAT))
+          + ROOTVSMALL
+        );
 }
 
 
@@ -238,7 +279,8 @@ Foam::scalar Foam::tetrahedron<Point, PointRef>::barycentric
                 "const point& pt"
             ") const"
         )
-            << "Degenerate tetrahedron - returning 1/4 barycentric coordinates."
+            << "Degenerate tetrahedron:" << nl << *this << nl
+            << "Returning 1/4 barycentric coordinates."
             << endl;
 
         bary = List<scalar>(4, 0.25);

src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H

@@ -109,9 +109,9 @@ inline Foam::vector Foam::triangle<Point, PointRef>::normal() const
 template<class Point, class PointRef>
 inline Point Foam::triangle<Point, PointRef>::circumCentre() const
 {
-    scalar d1 = (c_ - a_)&(b_ - a_);
-    scalar d2 = -(c_ - b_)&(b_ - a_);
-    scalar d3 = (c_ - a_)&(c_ - b_);
+    scalar d1 =  (c_ - a_) & (b_ - a_);
+    scalar d2 = -(c_ - b_) & (b_ - a_);
+    scalar d3 =  (c_ - a_) & (c_ - b_);
 
     scalar c1 = d2*d3;
     scalar c2 = d3*d1;
@@ -119,6 +119,16 @@ inline Point Foam::triangle<Point, PointRef>::circumCentre() const
 
     scalar c = c1 + c2 + c3;
 
+    if (Foam::mag(c) < ROOTVSMALL)
+    {
+        WarningIn("Point triangle<Point, PointRef>::circumCentre() const")
+            << "Degenerate triangle:" << nl << *this << nl
+            << "Returning centre instead of circumCentre."
+            << endl;
+
+        return centre();
+    }
+
     return
     (
         ((c2 + c3)*a_ + (c3 + c1)*b_ + (c1 + c2)*c_)/(2*c)
@@ -129,14 +139,19 @@ inline Point Foam::triangle<Point, PointRef>::circumCentre() const
 template<class Point, class PointRef>
 inline Foam::scalar Foam::triangle<Point, PointRef>::circumRadius() const
 {
-    scalar d1 = (c_ - a_) & (b_ - a_);
-    scalar d2 = - (c_ - b_) & (b_ - a_);
-    scalar d3 = (c_ - a_) & (c_ - b_);
+    scalar d1 =  (c_ - a_) & (b_ - a_);
+    scalar d2 = -(c_ - b_) & (b_ - a_);
+    scalar d3 =  (c_ - a_) & (c_ - b_);
 
     scalar denom = d2*d3 + d3*d1 + d1*d2;
 
     if (Foam::mag(denom) < VSMALL)
     {
+        WarningIn("scalar triangle<Point, PointRef>::circumRadius() const")
+            << "Degenerate triangle:" << nl << *this << nl
+            << "Returning GREAT for circumRadius."
+            << endl;
+
         return GREAT;
     }
     else
@@ -151,16 +166,7 @@ inline Foam::scalar Foam::triangle<Point, PointRef>::circumRadius() const
 template<class Point, class PointRef>
 inline Foam::scalar Foam::triangle<Point, PointRef>::quality() const
 {
-    // Note: 3*sqr(3)/(4*pi) = 0.4134966716
-
-    return
-        mag()
-      / (
-            constant::mathematical::pi
-           *Foam::sqr(circumRadius())
-           *0.4134966716
-          + VSMALL
-        );
+    return mag()/(Foam::sqr(circumRadius())*3.0*sqrt(3.0)/4.0 + VSMALL);
 }
 
 
@@ -267,7 +273,8 @@ Foam::scalar Foam::triangle<Point, PointRef>::barycentric
                 "const point& pt"
             ") const"
         )
-            << "Degenerate triangle - returning 1/3 barycentric coordinates."
+            << "Degenerate triangle:" << nl << *this << nl
+            << "Returning 1/3 barycentric coordinates."
             << endl;
 
         bary = List<scalar>(3, 1.0/3.0);
@@ -493,7 +500,7 @@ Foam::pointHit Foam::triangle<Point, PointRef>::nearestPointClassify
 ) const
 {
     // Adapted from:
-    // Real-time collision detection, Christer Ericson, 2005, 136-142
+    // Real-time collision detection, Christer Ericson, 2005, p136-142
 
     // Check if P in vertex region outside A
     vector ab = b_ - a_;
@@ -531,6 +538,27 @@ Foam::pointHit Foam::triangle<Point, PointRef>::nearestPointClassify
 
     if (vc <= 0.0 && d1 >= 0.0 && d3 <= 0.0)
     {
+        if ((d1 - d3) < ROOTVSMALL)
+        {
+            WarningIn
+            (
+                "pointHit triangle<Point, PointRef>::nearestPointClassify"
+                "("
+                    "const point& p,"
+                    "label& nearType,"
+                    "label& nearLabel"
+                ") const"
+            )
+                << "Degenerate triangle:" << nl << *this << nl
+                << "d1, d3: " << d1 << ", " << d3 << endl;
+
+            // For d1 = d3, a_ and b_ are likely coincident.
+
+            nearType = POINT;
+            nearLabel = 0;
+            return pointHit(false, a_, Foam::mag(a_ - p), true);
+        }
+
         // barycentric coordinates (1-v, v, 0)
         scalar v = d1/(d1 - d3);
 
@@ -559,6 +587,27 @@ Foam::pointHit Foam::triangle<Point, PointRef>::nearestPointClassify
 
     if (vb <= 0.0 && d2 >= 0.0 && d6 <= 0.0)
     {
+        if ((d2 - d6) < ROOTVSMALL)
+        {
+            WarningIn
+            (
+                "pointHit triangle<Point, PointRef>::nearestPointClassify"
+                "("
+                    "const point& p,"
+                    "label& nearType,"
+                    "label& nearLabel"
+                ") const"
+            )
+                << "Degenerate triangle:" << nl << *this << nl
+                << "d2, d6: " << d2 << ", " << d6 << endl;
+
+            // For d2 = d6, a_ and c_ are likely coincident.
+
+            nearType = POINT;
+            nearLabel = 0;
+            return pointHit(false, a_, Foam::mag(a_ - p), true);
+        }
+
         // barycentric coordinates (1-w, 0, w)
         scalar w = d2/(d2 - d6);
 
@@ -573,6 +622,28 @@ Foam::pointHit Foam::triangle<Point, PointRef>::nearestPointClassify
 
     if (va <= 0.0 && (d4 - d3) >= 0.0 && (d5 - d6) >= 0.0)
     {
+        if (((d4 - d3) + (d5 - d6)) < ROOTVSMALL)
+        {
+            WarningIn
+            (
+                "pointHit triangle<Point, PointRef>::nearestPointClassify"
+                "("
+                    "const point& p,"
+                    "label& nearType,"
+                    "label& nearLabel"
+                ") const"
+            )
+                << "Degenerate triangle:" << nl << *this << nl
+                << "(d4 - d3), (d6 - d5): " << (d4 - d3) << ", " << (d6 - d5)
+                << endl;
+
+            // For (d4 - d3) = (d6 - d5), b_ and c_ are likely coincident.
+
+            nearType = POINT;
+            nearLabel = 1;
+            return pointHit(false, b_, Foam::mag(b_ - p), true);
+        }
+
         // barycentric coordinates (0, 1-w, w)
         scalar w = (d4 - d3)/((d4 - d3) + (d5 - d6));
 
@@ -584,6 +655,28 @@ Foam::pointHit Foam::triangle<Point, PointRef>::nearestPointClassify
 
     // P inside face region. Compute Q through its barycentric
     // coordinates (u, v, w)
+
+    if ((va + vb + vc) < ROOTVSMALL)
+    {
+        WarningIn
+        (
+            "pointHit triangle<Point, PointRef>::nearestPointClassify"
+            "("
+                "const point& p,"
+                "label& nearType,"
+                "label& nearLabel"
+            ") const"
+        )
+            << "Degenerate triangle:" << nl << *this << nl
+            << "va, vb, vc: " << va << ", " << vb << ", " << vc
+            << endl;
+
+        point nearPt = centre();
+        nearType = NONE,
+        nearLabel = -1;
+        return pointHit(true, nearPt, Foam::mag(nearPt - p), false);
+    }
+
     scalar denom = 1.0/(va + vb + vc);
     scalar v = vb * denom;
     scalar w = vc * denom;

src/OpenFOAM/primitives/hashes/SHA1/SHA1Digest.C

@@ -32,7 +32,37 @@ License
 
 //! @cond fileScope
 const char hexChars[] = "0123456789abcdef";
-//! @endcond
+//! @endcond fileScope
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+unsigned char Foam::SHA1Digest::readHexDigit(Istream& is)
+{
+    // Takes into account that 'a' (or 'A') is 10
+    static const label alphaOffset = toupper('A') - 10;
+    // Takes into account that '0' is 0
+    static const label zeroOffset = int('0');
+
+    char c = 0;
+    is.read(c);
+
+    if (!isxdigit(c))
+    {
+        FatalIOErrorIn("SHA1Digest::readHexDigit(Istream&)", is)
+            << "Illegal hex digit: '" << c << "'"
+            << exit(FatalIOError);
+    }
+
+    if (isdigit(c))
+    {
+        return int(c) - zeroOffset;
+    }
+    else
+    {
+        return toupper(c) - alphaOffset;
+    }
+}
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -43,6 +73,12 @@ Foam::SHA1Digest::SHA1Digest()
 }
 
 
+Foam::SHA1Digest::SHA1Digest(Istream& is)
+{
+    is  >> *this;
+}
+
+
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 void Foam::SHA1Digest::clear()
@@ -75,6 +111,23 @@ bool Foam::SHA1Digest::operator!=(const SHA1Digest& rhs) const
 
 // * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //
 
+Foam::Istream& Foam::operator>>(Istream& is, SHA1Digest& dig)
+{
+    unsigned char *v = dig.v_;
+
+    for (unsigned i = 0; i < dig.length; ++i)
+    {
+        unsigned char c1 = SHA1Digest::readHexDigit(is);
+        unsigned char c2 = SHA1Digest::readHexDigit(is);
+
+        v[i] = (c1 << 4) + c2;
+    }
+
+    is.check("Istream& operator>>(Istream&, SHA1Digest&)");
+    return is;
+}
+
+
 Foam::Ostream& Foam::operator<<(Ostream& os, const SHA1Digest& dig)
 {
     const unsigned char *v = dig.v_;

src/OpenFOAM/primitives/hashes/SHA1/SHA1Digest.H

@@ -45,11 +45,13 @@ namespace Foam
 
 // Forward declaration of classes
 class Ostream;
+class Istream;
 
 // Forward declaration of friend functions and operators
 class SHA1;
 class SHA1Digest;
 Ostream& operator<<(Ostream&, const SHA1Digest&);
+Istream& operator>>(Istream&, SHA1Digest&);
 
 
 /*---------------------------------------------------------------------------*\
@@ -67,6 +69,9 @@ public:
     //- Construct a zero digest
     SHA1Digest();
 
+    //- Construct read a digest
+    SHA1Digest(Istream&);
+
     //- Reset the digest to zero
     void clear();
 
@@ -77,11 +82,14 @@ public:
     bool operator!=(const SHA1Digest&) const;
 
     friend Ostream& operator<<(Ostream&, const SHA1Digest&);
+    friend Istream& operator>>(Istream&, SHA1Digest&);
 
 private:
 
     //- The digest contents
     unsigned char v_[length];
+
+    static unsigned char readHexDigit(Istream&);
 };
 
 

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/strings/string/string.C

@@ -94,7 +94,7 @@ Foam::string& Foam::string::replaceAll
 
 
 // Expand all occurences of environment variables and initial tilde sequences
-Foam::string& Foam::string::expand()
+Foam::string& Foam::string::expand(const bool recurse)
 {
     size_type startEnvar = 0;
 
@@ -140,6 +140,10 @@ Foam::string& Foam::string::expand()
 
                 if (enVarString.size())
                 {
+                    if (recurse)
+                    {
+                        enVarString.expand();
+                    }
                     std::string::replace
                     (
                         startEnvar,

src/OpenFOAM/primitives/strings/string/string.H

@@ -182,7 +182,7 @@ public:
         //
         //  @sa
         //  Foam::findEtcFile
-        string& expand();
+        string& expand(const bool recurse=false);
 
         //- Remove repeated characters returning true if string changed
         bool removeRepeated(const char);

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
-      + (1 - (n1 & n2))*sqr(n1 ^ n2)/(magSqr(n1 ^ n2) + VSMALL)
+        s*I
+      + (1 - s)*sqr(n3)/(magSqr(n3) + VSMALL)
       + (n2*n1 - n1*n2);
 }
 

src/dynamicFvMesh/include/createDynamicFvMesh.H

@@ -9,7 +9,8 @@
             (
                 dynamicFvMesh::defaultRegion,
                 runTime.timeName(),
-                runTime
+                runTime,
+                IOobject::MUST_READ
             )
         )
     );

src/dynamicMesh/meshCut/directions/directionInfo/directionInfo.C

@@ -25,11 +25,8 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "directionInfo.H"
-//#include "hexMatcher.H"
-//#include "meshTools.H"
 #include "polyMesh.H"
 
-
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 // Find edge among edgeLabels that uses v0 and v1
@@ -209,13 +206,44 @@ Foam::Ostream& Foam::operator<<
     const Foam::directionInfo& wDist
 )
 {
-    return os << wDist.index_ << wDist.n_;
+    if (os.format() == IOstream::ASCII)
+    {
+        os << wDist.index_ << wDist.n_;
+    }
+    else
+    {
+        os.write
+        (
+            reinterpret_cast<const char*>(&wDist.index_),
+            sizeof(directionInfo)
+        );
+    }
+
+    // Check state of Ostream
+    os.check("Ostream& operator<<(Ostream&, const directionInfo&)");
+    return os;
+
 }
 
 
 Foam::Istream& Foam::operator>>(Foam::Istream& is, Foam::directionInfo& wDist)
 {
-    return is >> wDist.index_ >> wDist.n_;
+    if (is.format() == IOstream::ASCII)
+    {
+        is >> wDist.index_ >> wDist.n_;
+    }
+    else
+    {
+        is.read
+        (
+            reinterpret_cast<char*>(&wDist.index_),
+            sizeof(directionInfo)
+        );
+    }
+
+    // Check state of Istream
+    is.check("Istream& operator>>(Istream&, directionInfo&)");
+    return is;
 }
 
 // ************************************************************************* //

src/dynamicMesh/meshCut/directions/directionInfo/directionInfo.H

@@ -251,6 +251,14 @@ public:
 };
 
 
+//- Data associated with directionInfo type are contiguous
+template<>
+inline bool contiguous<directionInfo>()
+{
+    return true;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/dynamicMesh/meshCut/wallLayerCells/wallNormalInfo/wallNormalInfo.C

@@ -33,12 +33,42 @@ Foam::Ostream& Foam::operator<<
     const Foam::wallNormalInfo& wDist
 )
 {
-    return os << wDist.normal();
+    if (os.format() == IOstream::ASCII)
+    {
+        os << wDist.normal();
+    }
+    else
+    {
+        os.write
+        (
+            reinterpret_cast<const char*>(&wDist.normal_),
+            sizeof(vector)
+        );
+    }
+
+    // Check state of Ostream
+    os.check("Ostream& operator<<(Ostream&, const wallNormalInfo&)");
+    return os;
 }
 
 Foam::Istream& Foam::operator>>(Foam::Istream& is, Foam::wallNormalInfo& wDist)
 {
-    return is >> wDist.normal_;
+    if (is.format() == IOstream::ASCII)
+    {
+        is >> wDist.normal_;
+    }
+    else
+    {
+        is.read
+        (
+            reinterpret_cast<char*>(&wDist.normal_),
+            sizeof(vector)
+        );
+    }
+
+    // Check state of Istream
+    is.check("Istream& operator>>(Istream&, wallNormalInfo&)");
+    return is;
 }
 
 // ************************************************************************* //

src/dynamicMesh/meshCut/wallLayerCells/wallNormalInfo/wallNormalInfo.H

@@ -191,6 +191,14 @@ public:
 };
 
 
+//- Data associated with wallNormalInfo type are contiguous
+template<>
+inline bool contiguous<wallNormalInfo>()
+{
+    return true;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/dynamicMesh/polyTopoChange/polyTopoChange/refinementData.C

@@ -33,13 +33,43 @@ Foam::Ostream& Foam::operator<<
     const Foam::refinementData& wDist
 )
 {
-    return os << wDist.refinementCount_ << token::SPACE << wDist.count_;
+    if (os.format() == IOstream::ASCII)
+    {
+        os << wDist.refinementCount_ << token::SPACE << wDist.count_;
+    }
+    else
+    {
+        os.write
+        (
+            reinterpret_cast<const char*>(&wDist.refinementCount_),
+            sizeof(refinementData)
+        );
+    }
+
+    // Check state of Ostream
+    os.check("Ostream& operator<<(Ostream&, const refinementData&)");
+    return os;
 }
 
 
 Foam::Istream& Foam::operator>>(Foam::Istream& is, Foam::refinementData& wDist)
 {
-    return is >> wDist.refinementCount_ >> wDist.count_;
+    if (is.format() == IOstream::ASCII)
+    {
+        is >> wDist.refinementCount_ >> wDist.count_;
+    }
+    else
+    {
+        is.read
+        (
+            reinterpret_cast<char*>(&wDist.refinementCount_),
+            sizeof(refinementData)
+        );
+    }
+
+    // Check state of Istream
+    is.check("Istream& operator>>(Istream&, refinementData&)");
+    return is;
 }
 
 

src/dynamicMesh/polyTopoChange/polyTopoChange/refinementData.H

@@ -208,6 +208,14 @@ public:
 };
 
 
+//- Data associated with refinementData type are contiguous
+template<>
+inline bool contiguous<refinementData>()
+{
+    return true;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/dynamicMesh/polyTopoChange/polyTopoChange/refinementDistanceData.H

@@ -242,6 +242,14 @@ public:
 };
 
 
+//- Data associated with refinementDistanceData type are contiguous
+template<>
+inline bool contiguous<refinementDistanceData>()
+{
+    return true;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/finiteVolume/Make/files

@@ -227,6 +227,7 @@ $(schemes)/clippedLinear/clippedLinear.C
 $(schemes)/harmonic/harmonic.C
 $(schemes)/fixedBlended/fixedBlended.C
 $(schemes)/localBlended/localBlended.C
+$(schemes)/limiterBlended/limiterBlended.C
 $(schemes)/localMax/localMax.C
 $(schemes)/localMin/localMin.C
 

src/finiteVolume/cfdTools/general/fieldSources/basicSource/basicSource/IObasicSourceList.C

@@ -39,7 +39,7 @@ Foam::IObasicSourceList::IObasicSourceList
             "sourcesProperties",
             mesh.time().constant(),
             mesh,
-            IOobject::MUST_READ,
+            IOobject::MUST_READ_IF_MODIFIED,
             IOobject::NO_WRITE
         )
     ),

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

@@ -167,7 +167,6 @@ Foam::advectiveFvPatchField<Type>::advectionSpeed() const
             phiName_
         );
 
-
     if (phi.dimensions() == dimDensity*dimVelocity*dimArea)
     {
         const fvPatchScalarField& rhop =

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
-    (
-        phiName_,
-        reinterpret_cast<const surfaceScalarField*>(0),
-        reinterpret_cast<const scalar*>(0)
-    );
+    const Field<scalar>& phip =
+        this->patch().template lookupPatchField<surfaceScalarField, scalar>
+        (
+            phiName_
+        );
 
     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
-    (
-        "phi",
-        reinterpret_cast<const surfaceScalarField*>(0),
-        reinterpret_cast<const scalar*>(0)
-    );
+    const fvsPatchField<scalar>& phip =
+        this->patch().template lookupPatchField<surfaceScalarField, scalar>
+        (
+            phiName_
+        );
 
     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.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -129,7 +129,6 @@ selfContainedDirectMappedFixedValueFvPatchField
 {}
 
 
-
 template<class Type>
 selfContainedDirectMappedFixedValueFvPatchField<Type>::
 selfContainedDirectMappedFixedValueFvPatchField

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -48,7 +48,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-              Class selfContainedDirectMappedFixedValueFvPatch Declaration
+      Class selfContainedDirectMappedFixedValueFvPatchField 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
@@ -186,6 +188,7 @@ public:
             //- Update the coefficients associated with the patch field
             virtual void updateCoeffs();
 
+
         //- Write
         virtual void write(Ostream&) const;
 };

src/finiteVolume/fields/fvPatchFields/derived/selfContainedDirectMapped/selfContainedDirectMappedFixedValueFvPatchFields.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

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

@@ -107,32 +107,28 @@ Foam::tmp<Foam::scalarField>
 Foam::waveTransmissiveFvPatchField<Type>::advectionSpeed() const
 {
     // Lookup the velocity and compressibility of the patch
-    const fvPatchField<scalar>& psip = this->patch().lookupPatchField
-    (
-        psiName_,
-        reinterpret_cast<const volScalarField*>(0),
-        reinterpret_cast<const scalar*>(0)
-    );
+    const fvPatchField<scalar>& psip =
+        this->patch().template lookupPatchField<volScalarField, scalar>
+        (
+            psiName_
+        );
 
     const surfaceScalarField& phi =
-        this->db().objectRegistry::template lookupObject<surfaceScalarField>
-        (this->phiName_);
+        this->db().template lookupObject<surfaceScalarField>(this->phiName_);
 
-    fvsPatchField<scalar> phip = this->patch().lookupPatchField
-    (
-        this->phiName_,
-        reinterpret_cast<const surfaceScalarField*>(0),
-        reinterpret_cast<const scalar*>(0)
-    );
+    fvsPatchField<scalar> phip =
+        this->patch().template lookupPatchField<surfaceScalarField, scalar>
+        (
+            this->phiName_
+        );
 
     if (phi.dimensions() == dimDensity*dimVelocity*dimArea)
     {
-        const fvPatchScalarField& rhop = this->patch().lookupPatchField
-        (
-            this->rhoName_,
-            reinterpret_cast<const volScalarField*>(0),
-            reinterpret_cast<const scalar*>(0)
-        );
+        const fvPatchScalarField& rhop =
+            this->patch().template lookupPatchField<volScalarField, scalar>
+            (
+                this->rhoName_
+            );
 
         phip /= rhop;
     }
@@ -161,6 +157,7 @@ void Foam::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)

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/finiteVolume/finiteVolume/fvc/fvcSmooth/smoothData.H

@@ -219,6 +219,14 @@ public:
 };
 
 
+//- Data associated with smoothData type are contiguous
+template<>
+inline bool contiguous<smoothData>()
+{
+    return true;
+}
+
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/finiteVolume/finiteVolume/fvc/fvcSmooth/sweepData.H

@@ -210,6 +210,14 @@ public:
 };
 
 
+//- Data associated with sweepData type are contiguous
+template<>
+inline bool contiguous<sweepData>()
+{
+    return true;
+}
+
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/finiteVolume/fvMesh/wallDist/wallPointYPlus/wallPointYPlus.H

@@ -144,6 +144,14 @@ public:
 };
 
 
+//- Data associated with pointEdgePoint type as contiguous as underlying type
+template<>
+inline bool contiguous<wallPointYPlus>()
+{
+    return contiguous<wallPointData<scalar> >();
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 

src/finiteVolume/interpolation/interpolation/interpolationCellPoint/interpolationCellPoint.C

@@ -36,7 +36,10 @@ Foam::interpolationCellPoint<Type>::interpolationCellPoint
 :
     interpolation<Type>(psi),
     psip_(volPointInterpolation::New(psi.mesh()).interpolate(psi))
-{}
+{
+    // Uses cellPointWeight to do interpolation which needs tet decomposition
+    (void)psi.mesh().tetBasePtIs();
+}
 
 
 // ************************************************************************* //

src/finiteVolume/interpolation/surfaceInterpolation/schemes/limiterBlended/limiterBlended.C

@@ -0,0 +1,36 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvMesh.H"
+#include "limiterBlended.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    makeSurfaceInterpolationScheme(limiterBlended)
+}
+
+// ************************************************************************* //

src/finiteVolume/interpolation/surfaceInterpolation/schemes/limiterBlended/limiterBlended.H

@@ -0,0 +1,229 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::limiterBlended
+
+Description
+    Blends two specified schemes using the limiter function provided by a
+    limitedSurfaceInterpolationScheme.
+
+    The limited scheme is specified first followed by the scheme to be scaled
+    by the limiter and then the scheme scaled by 1 - limiter.
+
+SourceFiles
+    limiterBlended.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef limiterBlended_H
+#define limiterBlended_H
+
+#include "limitedSurfaceInterpolationScheme.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class limiterBlended Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class limiterBlended
+:
+    public surfaceInterpolationScheme<Type>
+{
+    // Private Member Functions
+
+        //- Limited scheme providing the limiter
+        tmp<limitedSurfaceInterpolationScheme<Type> > tLimitedScheme_;
+
+        //- Scheme 1
+        tmp<surfaceInterpolationScheme<Type> > tScheme1_;
+
+        //- Scheme 2
+        tmp<surfaceInterpolationScheme<Type> > tScheme2_;
+
+
+        //- Disallow default bitwise copy construct
+        limiterBlended(const limiterBlended&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const limiterBlended&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("limiterBlended");
+
+
+    // Constructors
+
+        //- Construct from mesh and Istream.
+        //  The name of the flux field is read from the Istream and looked-up
+        //  from the mesh objectRegistry
+        limiterBlended
+        (
+            const fvMesh& mesh,
+            Istream& is
+        )
+        :
+            surfaceInterpolationScheme<Type>(mesh),
+            tLimitedScheme_
+            (
+                limitedSurfaceInterpolationScheme<Type>::New(mesh, is)
+            ),
+            tScheme1_
+            (
+                surfaceInterpolationScheme<Type>::New(mesh, is)
+            ),
+            tScheme2_
+            (
+                surfaceInterpolationScheme<Type>::New(mesh, is)
+            )
+        {}
+
+        //- Construct from mesh, faceFlux and Istream
+        limiterBlended
+        (
+            const fvMesh& mesh,
+            const surfaceScalarField& faceFlux,
+            Istream& is
+        )
+        :
+            surfaceInterpolationScheme<Type>(mesh),
+            tLimitedScheme_
+            (
+                limitedSurfaceInterpolationScheme<Type>::New(mesh, faceFlux, is)
+            ),
+            tScheme1_
+            (
+                surfaceInterpolationScheme<Type>::New(mesh, faceFlux, is)
+            ),
+            tScheme2_
+            (
+                surfaceInterpolationScheme<Type>::New(mesh, faceFlux, is)
+            )
+        {}
+
+
+    // Member Functions
+
+        //- Return the interpolation weighting factors
+        tmp<surfaceScalarField> weights
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& vf
+        ) const
+        {
+            surfaceScalarField blendingFactor =
+                tLimitedScheme_().limiter(vf);
+
+            return
+                blendingFactor*tScheme1_().weights(vf)
+              + (scalar(1) - blendingFactor)*tScheme2_().weights(vf);
+        }
+
+        //- Return the face-interpolate of the given cell field
+        //  with explicit correction
+        tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
+        interpolate(const GeometricField<Type, fvPatchField, volMesh>& vf) const
+        {
+            surfaceScalarField blendingFactor =
+                tLimitedScheme_().limiter(vf);
+
+            return
+                blendingFactor*tScheme1_().interpolate(vf)
+              + (scalar(1) - blendingFactor)*tScheme2_().interpolate(vf);
+        }
+
+
+        //- Return true if this scheme uses an explicit correction
+        virtual bool corrected() const
+        {
+            return tScheme1_().corrected() || tScheme2_().corrected();
+        }
+
+
+        //- Return the explicit correction to the face-interpolate
+        //  for the given field
+        virtual tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
+        correction
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& vf
+        ) const
+        {
+            surfaceScalarField blendingFactor =
+                tLimitedScheme_().limiter(vf);
+
+            if (tScheme1_().corrected())
+            {
+                if (tScheme2_().corrected())
+                {
+                    return
+                    (
+                        blendingFactor
+                      * tScheme1_().correction(vf)
+                      + (scalar(1.0) - blendingFactor)
+                      * tScheme2_().correction(vf)
+                    );
+                }
+                else
+                {
+                    return
+                    (
+                        blendingFactor
+                      * tScheme1_().correction(vf)
+                    );
+                }
+            }
+            else if (tScheme2_().corrected())
+            {
+                return
+                (
+                    (scalar(1.0) - blendingFactor)
+                  * tScheme2_().correction(vf)
+                );
+            }
+            else
+            {
+                return tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
+                (
+                    NULL
+                );
+            }
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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.
-            inline bool valid() const;
+            //- Check whether origin has been changed at all or
+            //  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;
 
 
@@ -195,6 +218,14 @@ public:
 };
 
 
+//- Data associated with pointEdgeStructuredWalk type are contiguous
+template<>
+inline bool contiguous<pointEdgeStructuredWalk>()
+{
+    return true;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

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;
@@ -179,7 +180,7 @@ void Foam::InteractionLists<ParticleType>::buildInteractionLists()
 
         treeBoundBox tempTransformedBb
         (
-            transform.invTransform(cellBbsToExchange[bbI].corners())
+            transform.invTransform(cellBbsToExchange[bbI].points())
         );
 
         treeBoundBox extendedBb
@@ -396,7 +397,7 @@ void Foam::InteractionLists<ParticleType>::buildInteractionLists()
 
         treeBoundBox tempTransformedBb
         (
-            transform.invTransform(wallFaceBbsToExchange[bbI].corners())
+            transform.invTransform(wallFaceBbsToExchange[bbI].points())
         );
 
         treeBoundBox extendedBb
@@ -701,7 +702,7 @@ void Foam::InteractionLists<ParticleType>::findExtendedProcBbsInRange
                         (
                             transform.transform
                             (
-                                allExtendedProcBbs[procI].corners()
+                                allExtendedProcBbs[procI].points()
                             )
                         );
 
@@ -748,7 +749,7 @@ void Foam::InteractionLists<ParticleType>::findExtendedProcBbsInRange
                     (
                         transform.transform
                         (
-                            allExtendedProcBbs[procI].corners()
+                            allExtendedProcBbs[procI].points()
                         )
                     );
 
@@ -791,7 +792,7 @@ void Foam::InteractionLists<ParticleType>::findExtendedProcBbsInRange
                 (
                     transform.transform
                     (
-                        allExtendedProcBbs[procI].corners()
+                        allExtendedProcBbs[procI].points()
                     )
                 );
 

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/lagrangian/coalCombustion/Make/options

@@ -8,8 +8,8 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/liquids/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/liquidMixture/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/solids/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/solidMixture/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/pointSolids/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/pointSolidMixture/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
@@ -30,8 +30,8 @@ LIB_LIBS = \
     -lbasicThermophysicalModels \
     -lliquids \
     -lliquidMixture \
-    -lsolids \
-    -lsolidMixture \
+    -lpointSolids \
+    -lpointSolidMixture \
     -lreactionThermophysicalModels \
     -lSLGThermo \
     -lcompressibleRASModels \

src/lagrangian/intermediate/Make/options

@@ -7,8 +7,8 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/liquids/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/liquidMixture/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/solids/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/solidMixture/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/pointSolids/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/pointSolidMixture/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
@@ -27,8 +27,8 @@ LIB_LIBS = \
     -lpdf \
     -lliquids \
     -lliquidMixture \
-    -lsolids \
-    -lsolidMixture \
+    -lpointSolids \
+    -lpointSolidMixture \
     -lspecie \
     -lbasicThermophysicalModels \
     -lreactionThermophysicalModels \

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/PatchInjection/PatchInjection.C

@@ -38,7 +38,24 @@ Foam::label Foam::PatchInjection<CloudType>::parcelsToInject
 {
     if ((time0 >= 0.0) && (time0 < duration_))
     {
-        return round(fraction_*(time1 - time0)*parcelsPerSecond_);
+        scalar nParcels =fraction_*(time1 - time0)*parcelsPerSecond_;
+
+        cachedRandom& rnd = this->owner().rndGen();
+
+        label nParcelsToInject = floor(nParcels);
+
+        // Inject an additional parcel with a probability based on the
+        // remainder after the floor function
+        if
+        (
+            nParcelsToInject > 0
+         && (nParcels - scalar(nParcelsToInject) > rnd.position(0.0, 1.0))
+        )
+        {
+            ++nParcelsToInject;
+        }
+
+        return nParcelsToInject;
     }
     else
     {
@@ -173,6 +190,7 @@ void Foam::PatchInjection<CloudType>::setPositionAndCell
     if (cellOwners_.size() > 0)
     {
         cachedRandom& rnd = this->owner().rndGen();
+
         label cellI = rnd.position<label>(0, cellOwners_.size() - 1);
 
         cellOwner = cellOwners_[cellI];

src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.C

@@ -100,7 +100,7 @@ const Foam::liquidMixture& Foam::CompositionModel<CloudType>::liquids() const
 
 
 template<class CloudType>
-const Foam::solidMixture& Foam::CompositionModel<CloudType>::solids() const
+const Foam::pointSolidMixture& Foam::CompositionModel<CloudType>::solids() const
 {
     return thermo_.solids();
 }

src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.H

@@ -143,7 +143,7 @@ public:
                 const liquidMixture& liquids() const;
 
                 //- Return the global (additional) solids
-                const solidMixture& solids() const;
+                const pointSolidMixture& solids() const;
 
                 //- Return the list of phase properties
                 const phasePropertiesList& phaseProps() const;

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_.oldInstance());
-        }
+        newMesh.setInstance(meshRefiner_.timeName());
 
         // 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.
-    (so not standard meshWave)
-
-    To be used in wall distance calculation.
+    Variant of pointEdgePoint with some transported additional data.
+    WIP - should be templated on data like wallDistData.
+    Passive vector v_ is not a coordinate (so no enterDomain/leaveDomain
+    transformation needed)
 
 SourceFiles
     pointDataI.H
@@ -39,9 +39,7 @@ SourceFiles
 #ifndef pointData_H
 #define pointData_H
 
-#include "point.H"
-#include "label.H"
-#include "tensor.H"
+#include "pointEdgePoint.H"
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -49,51 +47,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 +87,6 @@ public:
 
         // Access
 
-            inline const point& origin() const;
-
-            inline scalar distSqr() const;
-
             inline scalar s() const;
 
             inline const vector& v() const;
@@ -129,81 +94,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&);
@@ -211,6 +155,14 @@ public:
 };
 
 
+//- Data associated with pointData as contiguous as pointEdgePoint
+template<>
+inline bool contiguous<pointData>()
+{
+    return contiguous<pointEdgePoint>();
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

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),
-    distSqr_(GREAT),
+    pointEdgePoint(),
     s_(GREAT),
     v_(point::max)
 {}
@@ -152,8 +46,7 @@ inline Foam::pointData::pointData
     const vector& v
 )
 :
-    origin_(origin),
-    distSqr_(distSqr),
+    pointEdgePoint(origin, 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()),
-    distSqr_(wpt.distSqr()),
+    pointEdgePoint(wpt),
     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
-{
-    return origin_ != point::max;
-}
-
-
-// Checks for cyclic points
-inline bool Foam::pointData::sameGeometry
+template <class TrackingData>
+inline void Foam::pointData::transform
 (
-    const pointData& w2,
-    const scalar tol
-) 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
+    const tensor& rotTensor,
+    TrackingData& td
 )
 {
-    origin_ -= coord;
-}
-
-
-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;
+    pointEdgePoint::transform(rotTensor, td);
+    v_ = Foam::transform(rotTensor, v_);
 }
 
 
 // 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
-        update
+    if
+    (
+        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
-        update
+    if
+    (
+        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();
-
-    const edge& e = mesh.edges()[edgeI];
-
-    const point edgeMid(0.5*(points[e[0]] + points[e[1]]));
-
-    return
-        update
+    if
+    (
+        pointEdgePoint::updateEdge
         (
-            edgeMid,
+            mesh,
+            edgeI,
+            pointI,
             pointInfo,
-            tol
-        );
-}
-
-
-// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
-
-inline bool Foam::pointData::operator==(const pointData& rhs) const
-{
-    return origin() == rhs.origin();
-}
-
-
-inline bool Foam::pointData::operator!=(const pointData& rhs) const
-{
-    return !(*this == rhs);
+            tol,
+            td
+        )
+    )
+    {
+        s_ = pointInfo.s_;
+        v_ = pointInfo.v_;
+        return true;
+    }
+    else
+    {
+        return false;
+    }
 }
 
 

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/Make/files

@@ -127,6 +127,7 @@ $(cellZoneSources)/setToCellZone/setToCellZone.C
 pointZoneSources = sets/pointZoneSources
 $(pointZoneSources)/setToPointZone/setToPointZone.C
 
+momentOfInertia/momentOfInertia.C
 
 surfaceSets/surfaceSets.C
 

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>
-Foam::scalar Foam::PointEdgeWave<Type>::propagationTol_ = 0.01;
+template <class Type, class TrackingData>
+Foam::scalar Foam::PointEdgeWave<Type, TrackingData>::propagationTol_ = 0.01;
 
-
-// Offset labelList. Used for transferring from one cyclic half to the other.
-template <class Type>
-void Foam::PointEdgeWave<Type>::offset(const label val, labelList& elems)
-{
-    forAll(elems, i)
-    {
-        elems[i] += val;
-    }
-}
+template <class Type, class TrackingData>
+Foam::label Foam::PointEdgeWave<Type, TrackingData>::dummyTrackData_ = 12345;
 
 
 // * * * * * * * * * * * * * 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>
-void Foam::PointEdgeWave<Type>::leaveDomain
+template <class Type, class TrackingData>
+void Foam::PointEdgeWave<Type, TrackingData>::leaveDomain
 (
-    const polyPatch& meshPatch,
-    const primitivePatch& patch,
+    const polyPatch& 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>
-void Foam::PointEdgeWave<Type>::enterDomain
+template <class Type, class TrackingData>
+void Foam::PointEdgeWave<Type, TrackingData>::enterDomain
 (
-    const polyPatch& meshPatch,
-    const primitivePatch& patch,
+    const polyPatch& 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>
-void Foam::PointEdgeWave<Type>::transform
+template <class Type, class TrackingData>
+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>&)"
-        )   << "Parallel cyclics not supported" << abort(FatalError);
+            "PointEdgeWave<Type, TrackingData>::transform"
+            "(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>
-bool Foam::PointEdgeWave<Type>::updatePoint
+template <class Type, class TrackingData>
+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>
-bool Foam::PointEdgeWave<Type>::updatePoint
+template <class Type, class TrackingData>
+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>
-bool Foam::PointEdgeWave<Type>::updateEdge
+template <class Type, class TrackingData>
+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>
-void Foam::PointEdgeWave<Type>::handleProcPatches()
+template <class Type, class TrackingData>
+void Foam::PointEdgeWave<Type, TrackingData>::handleProcPatches()
 {
-    // 1. Send all point info on processor patches. Send as
-    // face label + offset in face.
+    // 1. Send all point info on processor patches.
 
-    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
-
-            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)
+            label meshPointI = procPatch.meshPoints()[thisPointI];
+            if (changedPoint_[meshPointI])
             {
-                Pout<< "Processor patch " << patchI << ' ' << patch.name()
-                    << " communicating with " << procPatch.neighbProcNo()
-                    << "  Sending:" << patchInfo.size() << endl;
-            }
-
-            {
-                OPstream toNeighbour
-                (
-                    Pstream::blocking,
-                    procPatch.neighbProcNo()
-                );
-
-                toNeighbour << owner << ownerIndex << patchInfo;
+                patchInfo.append(allPointInfo_[meshPointI]);
+                thisPoints.append(thisPointI);
+                nbrPoints.append(neighbPoints[thisPointI]);
             }
         }
+
+        // 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 =
-                refCast<const processorPolyPatch>(patch);
+            UIPstream fromNeighbour(procPatch.neighbProcNo(), pBufs);
+            fromNeighbour >> patchPoints >> patchInfo;
+        }
 
-            List<Type> patchInfo;
-            labelList owner;
-            labelList ownerIndex;
+        if (debug)
+        {
+            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,
-                    procPatch.neighbProcNo()
+                    meshPointI,
+                    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
-        if (sharedData[pd.sharedPointAddr()[i]].valid())
+        // Retrieve my entries from the shared points.
+        const Type& nbrInfo = sharedData[pd.sharedPointAddr()[i]];
+
+        if (!allPointInfo_[meshPointI].equal(nbrInfo, td_))
         {
             updatePoint
             (
                 meshPointI,
-                sharedData[pd.sharedPointAddr()[i]],
-                propagationTol_,
+                nbrInfo,
                 allPointInfo_[meshPointI]
             );
         }
@@ -529,11 +433,14 @@ void Foam::PointEdgeWave<Type>::handleProcPatches()
 }
 
 
-template <class Type>
-void Foam::PointEdgeWave<Type>::handleCyclicPatches()
+template <class Type, class TrackingData>
+void Foam::PointEdgeWave<Type, TrackingData>::handleCyclicPatches()
 {
-    // 1. Send all point info on cyclic patches. Send as
-    // face label + offset in face.
+    // 1. Send all point info on cyclic patches.
+
+    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());
-            DynamicList<label> nbrPoints(nbrPatch.nPoints());
-            DynamicList<label> nbrOwner(nbrPatch.nPoints());
-            DynamicList<label> nbrIndex(nbrPatch.nPoints());
+            // Collect nbrPatch points that have changed
+            {
+                const cyclicPolyPatch& nbrPatch = cycPatch.neighbPatch();
+                const edgeList& pairs = cycPatch.coupledPoints();
+                const labelList& meshPoints = nbrPatch.meshPoints();
 
-            getChangedPatchPoints
-            (
-                nbrPatch,
-                nbrInfo,
-                nbrPoints,
-                nbrOwner,
-                nbrIndex
-            );
+                forAll(pairs, pairI)
+                {
+                    label thisPointI = pairs[pairI][0];
+                    label nbrPointI = pairs[pairI][1];
+                    label meshPointI = meshPoints[nbrPointI];
+
+                    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
-            updateFromPatchInfo
-            (
-                cycPatch,
-                cycPatch,
-                nbrOwner,
-                nbrIndex,
-                nbrInfo
-            );
+            // Adapt for entering domain
+            enterDomain(cycPatch, thisPoints, 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>
-Foam::PointEdgeWave<Type>::PointEdgeWave
+template <class Type, class TrackingData>
+Foam::PointEdgeWave<Type, TrackingData>::PointEdgeWave
 (
     const polyMesh& mesh,
     const labelList& changedPoints,
     const List<Type>& changedPointsInfo,
 
-    List<Type>& allPointInfo,
-    List<Type>& allEdgeInfo,
+    UList<Type>& allPointInfo,
+    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>
-Foam::PointEdgeWave<Type>::PointEdgeWave
+template <class Type, class TrackingData>
+Foam::PointEdgeWave<Type, TrackingData>::PointEdgeWave
 (
     const polyMesh& mesh,
-    List<Type>& allPointInfo,
-    List<Type>& allEdgeInfo
+    UList<Type>& allPointInfo,
+    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>
-Foam::PointEdgeWave<Type>::~PointEdgeWave()
+template <class Type, class TrackingData>
+Foam::PointEdgeWave<Type, TrackingData>::~PointEdgeWave()
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 
-template <class Type>
-Foam::label Foam::PointEdgeWave<Type>::getUnsetPoints() const
+template <class Type, class TrackingData>
+Foam::label Foam::PointEdgeWave<Type, TrackingData>::getUnsetPoints() const
 {
     return nUnvisitedPoints_;
 }
 
 
-template <class Type>
-Foam::label Foam::PointEdgeWave<Type>::getUnsetEdges() const
+template <class Type, class TrackingData>
+Foam::label Foam::PointEdgeWave<Type, TrackingData>::getUnsetEdges() const
 {
     return nUnvisitedEdges_;
 }
 
 
 // Copy point information into member data
-template <class Type>
-void Foam::PointEdgeWave<Type>::setPointInfo
+template <class Type, class TrackingData>
+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>
-Foam::label Foam::PointEdgeWave<Type>::edgeToPoint()
+template <class Type, class TrackingData>
+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>
-Foam::label Foam::PointEdgeWave<Type>::pointToEdge()
+template <class Type, class TrackingData>
+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>
-Foam::label Foam::PointEdgeWave<Type>::iterate(const label maxIter)
+template <class Type, class TrackingData>
+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 "tensor.H"
-#include "primitivePatch.H"
-#include "PtrList.H"
+#include "tensorField.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 primitivePatch& patch,
+            const polyPatch&,
             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 primitivePatch& patch,
+            const polyPatch&,
             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,
-            List<Type>& allEdgeInfo,
-            const label maxIter
+            UList<Type>& allPointInfo,
+            UList<Type>& allEdgeInfo,
+            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,
-            List<Type>& allEdgeInfo
+            UList<Type>& allPointInfo,
+            UList<Type>& allEdgeInfo,
+            TrackingData& td = dummyTrackData_
         );
 
 
@@ -291,18 +265,24 @@ public:
 
     // Member Functions
 
-        //- Get allPointInfo
-        const List<Type>& allPointInfo() const
+        //- Access allPointInfo
+        UList<Type>& allPointInfo() const
         {
             return allPointInfo_;
         }
 
-        //- Get allEdgeInfo
-        const List<Type>& allEdgeInfo() const
+        //- Access allEdgeInfo
+        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)
-             const;
+            template<class TrackingData>
+            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;
 
 
@@ -196,6 +226,14 @@ public:
 };
 
 
+//- Data associated with pointEdgePoint type are contiguous
+template<>
+inline bool contiguous<pointEdgePoint>()
+{
+    return true;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

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
-        update
-        (
-            mesh.points()[pointI],
-            edgeInfo,
-            tol
-        );
+    return update(mesh.points()[pointI], edgeInfo, tol, td);
 }
 
 
 // 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
-        update
-        (
-            mesh.points()[pointI],
-            newPointInfo,
-            tol
-        );
+    return update(mesh.points()[pointI], newPointInfo, tol, td);
 }
 
 
 // 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
-        update
-        (
-            edgeMid,
-            pointInfo,
-            tol
-        );
+template <class TrackingData>
+inline bool Foam::pointEdgePoint::equal
+(
+    const pointEdgePoint& rhs,
+    TrackingData& td
+) const
+{
+    return operator==(rhs);
 }
 
 

src/meshTools/cellClassification/cellInfo.H

@@ -201,6 +201,14 @@ public:
 };
 
 
+//- Data associated with cellInfo type are contiguous
+template<>
+inline bool contiguous<cellInfo>()
+{
+    return true;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/meshTools/cellDist/wallPoint/wallPoint.H

@@ -223,6 +223,14 @@ public:
 };
 
 
+//- Data associated with wallPoint type are contiguous
+template<>
+inline bool contiguous<wallPoint>()
+{
+    return true;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/meshTools/cellDist/wallPoint/wallPointData.H

@@ -160,6 +160,46 @@ public:
 };
 
 
+//- Data associated with wallPointData type are contiguous. List the usual
+//  ones.
+
+template <>
+inline bool contiguous<wallPointData<bool> >()
+{
+    return contiguous<wallPoint>();
+}
+template <>
+inline bool contiguous<wallPointData<label> >()
+{
+    return contiguous<wallPoint>();
+}
+template <>
+inline bool contiguous<wallPointData<scalar> >()
+{
+    return contiguous<wallPoint>();
+}
+template <>
+inline bool contiguous<wallPointData<vector> >()
+{
+    return contiguous<wallPoint>();
+}
+template <>
+inline bool contiguous<wallPointData<sphericalTensor> >()
+{
+    return contiguous<wallPoint>();
+}
+template <>
+inline bool contiguous<wallPointData<symmTensor> >()
+{
+    return contiguous<wallPoint>();
+}
+template <>
+inline bool contiguous<wallPointData<tensor> >()
+{
+    return contiguous<wallPoint>();
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 

src/meshTools/directMapped/directMappedPolyPatch/directMappedPatchBase.C

@@ -35,6 +35,7 @@ License
 #include "polyMesh.H"
 #include "polyPatch.H"
 #include "Time.H"
+#include "mapDistribute.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -65,32 +66,6 @@ namespace Foam
         "nonuniform",
         "normal"
     };
-
-    //- Private class for finding nearest
-    //  - point+local index
-    //  - sqr(distance)
-    //  - processor
-    typedef Tuple2<pointIndexHit, Tuple2<scalar, label> > nearInfo;
-
-    class nearestEqOp
-    {
-    public:
-
-        void operator()(nearInfo& x, const nearInfo& y) const
-        {
-            if (y.first().hit())
-            {
-                if (!x.first().hit())
-                {
-                    x = y;
-                }
-                else if (y.second().first() < x.second().first())
-                {
-                    x = y;
-                }
-            }
-        }
-    };
 }
 
 

src/meshTools/directMapped/directMappedPolyPatch/directMappedPatchBase.H

@@ -70,8 +70,6 @@ SourceFiles
 #include "pointField.H"
 #include "Tuple2.H"
 #include "pointIndexHit.H"
-#include "mapDistribute.H"
-
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -80,6 +78,7 @@ namespace Foam
 
 class polyPatch;
 class polyMesh;
+class mapDistribute;
 
 /*---------------------------------------------------------------------------*\
                       Class directMappedPatchBase Declaration
@@ -106,6 +105,34 @@ public:
             NORMAL              // use face normal + distance
         };
 
+
+    //- Helper class for finding nearest
+    //  - point+local index
+    //  - sqr(distance)
+    //  - processor
+    typedef Tuple2<pointIndexHit, Tuple2<scalar, label> > nearInfo;
+
+    class nearestEqOp
+    {
+
+    public:
+
+        void operator()(nearInfo& x, const nearInfo& y) const
+        {
+            if (y.first().hit())
+            {
+                if (!x.first().hit())
+                {
+                    x = y;
+                }
+                else if (y.second().first() < x.second().first())
+                {
+                    x = y;
+                }
+            }
+        }
+    };
+
 private:
 
     // Private data

src/meshTools/momentOfInertia/momentOfInertia.C

@@ -21,382 +21,328 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
-Class
-    momentOfInertia
-
-Description
-    Reimplementation of volInt.c by Brian Mirtich.
-        *  mirtich@cs.berkeley.edu                             *
-        *  http://www.cs.berkeley.edu/~mirtich                 *
-
--------------------------------------------------------------------------------
-*/
+\*---------------------------------------------------------------------------*/
 
 #include "momentOfInertia.H"
-//#include "pyramidPointFaceRef.H"
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 
-//Foam::tensor Foam::momentOfInertia
-//(
-//    const pointField& points,
-//    const faceList& faces,
-//    const cell& cFaces,
-//    const point& cc
-//)
-//{
-//    tensor t(tensor::zero);
-//
-//    forAll(cFaces, i)
-//    {
-//        const face& f = faces[cFaces[i]];
-//
-//        scalar pyrVol = pyramidPointFaceRef(f, cc).mag(points);
-//
-//        vector pyrCentre = pyramidPointFaceRef(f, cc).centre(points);
-//
-//        vector d = pyrCentre - cc;
-//
-//        t.xx() += pyrVol*(sqr(d.y()) + sqr(d.z()));
-//        t.yy() += pyrVol*(sqr(d.x()) + sqr(d.z()));
-//        t.zz() += pyrVol*(sqr(d.x()) + sqr(d.y()));
-//
-//        t.xy() -= pyrVol*d.x()*d.y();
-//        t.xz() -= pyrVol*d.x()*d.z();
-//        t.yz() -= pyrVol*d.y()*d.z();
-//    }
-//
-//    // Symmetric
-//    t.yx() = t.xy();
-//    t.zx() = t.xz();
-//    t.zy() = t.yz();
-//
-//    return t;
-//}
-
-
-#define sqr(x) ((x)*(x))
-#define pow3(x) ((x)*(x)*(x))
-
-// compute various integrations over projection of face
-void Foam::compProjectionIntegrals
+void Foam::momentOfInertia::massPropertiesSolid
 (
-    const pointField& points,
-    const face& f,
-    const direction A,
-    const direction B,
-
-    scalar& P1,
-    scalar& Pa,
-    scalar& Pb,
-    scalar& Paa,
-    scalar& Pab,
-    scalar& Pbb,
-    scalar& Paaa,
-    scalar& Paab,
-    scalar& Pabb,
-    scalar& Pbbb
-)
-{
-    P1 = Pa = Pb = Paa = Pab = Pbb = Paaa = Paab = Pabb = Pbbb = 0.0;
-
-    forAll(f, i)
-    {
-        scalar a0 = points[f[i]][A];
-        scalar b0 = points[f[i]][B];
-        scalar a1 = points[f[(i+1) % f.size()]][A];
-        scalar b1 = points[f[(i+1) % f.size()]][B];
-        scalar da = a1 - a0;
-        scalar db = b1 - b0;
-
-        scalar a0_2 = a0 * a0;
-        scalar a0_3 = a0_2 * a0;
-        scalar a0_4 = a0_3 * a0;
-
-        scalar b0_2 = b0 * b0;
-        scalar b0_3 = b0_2 * b0;
-        scalar b0_4 = b0_3 * b0;
-
-        scalar a1_2 = a1 * a1;
-        scalar a1_3 = a1_2 * a1;
-
-        scalar b1_2 = b1 * b1;
-        scalar b1_3 = b1_2 * b1;
-
-        scalar C1 = a1 + a0;
-
-        scalar Ca = a1*C1 + a0_2;
-        scalar Caa = a1*Ca + a0_3;
-        scalar Caaa = a1*Caa + a0_4;
-
-        scalar Cb = b1*(b1 + b0) + b0_2;
-        scalar Cbb = b1*Cb + b0_3;
-        scalar Cbbb = b1*Cbb + b0_4;
-
-        scalar Cab = 3*a1_2 + 2*a1*a0 + a0_2;
-        scalar Kab = a1_2 + 2*a1*a0 + 3*a0_2;
-
-        scalar Caab = a0*Cab + 4*a1_3;
-        scalar Kaab = a1*Kab + 4*a0_3;
-
-        scalar Cabb = 4*b1_3 + 3*b1_2*b0 + 2*b1*b0_2 + b0_3;
-        scalar Kabb = b1_3 + 2*b1_2*b0 + 3*b1*b0_2 + 4*b0_3;
-
-        P1 += db*C1;
-        Pa += db*Ca;
-        Paa += db*Caa;
-        Paaa += db*Caaa;
-        Pb += da*Cb;
-        Pbb += da*Cbb;
-        Pbbb += da*Cbbb;
-        Pab += db*(b1*Cab + b0*Kab);
-        Paab += db*(b1*Caab + b0*Kaab);
-        Pabb += da*(a1*Cabb + a0*Kabb);
-  }
-
-  P1 /= 2.0;
-  Pa /= 6.0;
-  Paa /= 12.0;
-  Paaa /= 20.0;
-  Pb /= -6.0;
-  Pbb /= -12.0;
-  Pbbb /= -20.0;
-  Pab /= 24.0;
-  Paab /= 60.0;
-  Pabb /= -60.0;
-}
-
-
-void Foam::compFaceIntegrals
-(
-    const pointField& points,
-    const face& f,
-    const vector& n,
-    const scalar w,
-    const direction A,
-    const direction B,
-    const direction C,
-
-    scalar& Fa,
-    scalar& Fb,
-    scalar& Fc,
-    scalar& Faa,
-    scalar& Fbb,
-    scalar& Fcc,
-    scalar& Faaa,
-    scalar& Fbbb,
-    scalar& Fccc,
-    scalar& Faab,
-    scalar& Fbbc,
-    scalar& Fcca
-)
-{
-    scalar P1, Pa, Pb, Paa, Pab, Pbb, Paaa, Paab, Pabb, Pbbb;
-
-    compProjectionIntegrals
-    (
-        points,
-        f,
-        A,
-        B,
-
-        P1,
-        Pa,
-        Pb,
-        Paa,
-        Pab,
-        Pbb,
-        Paaa,
-        Paab,
-        Pabb,
-        Pbbb
-    );
-
-    scalar k1 = 1 / n[C];
-    scalar k2 = k1 * k1;
-    scalar k3 = k2 * k1;
-    scalar k4 = k3 * k1;
-
-    Fa = k1 * Pa;
-    Fb = k1 * Pb;
-    Fc = -k2 * (n[A]*Pa + n[B]*Pb + w*P1);
-
-    Faa = k1 * Paa;
-    Fbb = k1 * Pbb;
-    Fcc = k3 * (sqr(n[A])*Paa + 2*n[A]*n[B]*Pab + sqr(n[B])*Pbb
-         + w*(2*(n[A]*Pa + n[B]*Pb) + w*P1));
-
-    Faaa = k1 * Paaa;
-    Fbbb = k1 * Pbbb;
-    Fccc = -k4 * (pow3(n[A])*Paaa + 3*sqr(n[A])*n[B]*Paab
-           + 3*n[A]*sqr(n[B])*Pabb + pow3(n[B])*Pbbb
-           + 3*w*(sqr(n[A])*Paa + 2*n[A]*n[B]*Pab + sqr(n[B])*Pbb)
-           + w*w*(3*(n[A]*Pa + n[B]*Pb) + w*P1));
-
-    Faab = k1 * Paab;
-    Fbbc = -k2 * (n[A]*Pabb + n[B]*Pbbb + w*Pbb);
-    Fcca = k3 * (sqr(n[A])*Paaa + 2*n[A]*n[B]*Paab + sqr(n[B])*Pabb
-         + w*(2*(n[A]*Paa + n[B]*Pab) + w*Pa));
-}
-
-
-void Foam::compVolumeIntegrals
-(
-    const pointField& points,
-    const faceList& faces,
-    const cell& cFaces,
-    const vectorField& fNorm,
-    const scalarField& fW,
-
-    scalar& T0,
-    vector& T1,
-    vector& T2,
-    vector& TP
-)
-{
-    T0 = 0;
-    T1 = vector::zero;
-    T2 = vector::zero;
-    TP = vector::zero;
-
-    forAll(cFaces, i)
-    {
-        const vector& n = fNorm[i];
-
-        scalar nx = mag(n[0]);
-        scalar ny = mag(n[1]);
-        scalar nz = mag(n[2]);
-
-        direction A, B, C;
-
-        if (nx > ny && nx > nz)
-        {
-            C = 0;
-        }
-        else
-        {
-            C = (ny > nz) ? 1 : 2;
-        }
-
-        A = (C + 1) % 3;
-        B = (A + 1) % 3;
-
-        scalar Fa, Fb, Fc, Faa, Fbb, Fcc, Faaa, Fbbb, Fccc, Faab, Fbbc, Fcca;
-        compFaceIntegrals
-        (
-            points,
-            faces[cFaces[i]],
-            n,
-            fW[i],
-            A,
-            B,
-            C,
-
-            Fa,
-            Fb,
-            Fc,
-            Faa,
-            Fbb,
-            Fcc,
-            Faaa,
-            Fbbb,
-            Fccc,
-            Faab,
-            Fbbc,
-            Fcca
-        );
-
-        T0 += n[0] * ((A == 0) ? Fa : ((B == 0) ? Fb : Fc));
-
-        T1[A] += n[A] * Faa;
-        T1[B] += n[B] * Fbb;
-        T1[C] += n[C] * Fcc;
-
-        T2[A] += n[A] * Faaa;
-        T2[B] += n[B] * Fbbb;
-        T2[C] += n[C] * Fccc;
-
-        TP[A] += n[A] * Faab;
-        TP[B] += n[B] * Fbbc;
-        TP[C] += n[C] * Fcca;
-    }
-
-    T1 /= 2;
-    T2 /= 3;
-    TP /= 2;
-}
-
-
-// Calculate
-// - r: centre of mass
-// - J: inertia around origin (point 0,0,0)
-void Foam::momentOfIntertia
-(
-    const pointField& points,
-    const faceList& faces,
-    const cell& cFaces,
-    point& r,
+    const pointField& pts,
+    const triFaceList& triFaces,
+    scalar density,
+    scalar& mass,
+    vector& cM,
     tensor& J
 )
 {
-    // Face normals
-    vectorField fNorm(cFaces.size());
-    scalarField fW(cFaces.size());
+    // Reimplemented from: Wm4PolyhedralMassProperties.cpp
+    // File Version: 4.10.0 (2009/11/18)
 
-    forAll(cFaces, i)
+    // Geometric Tools, LC
+    // Copyright (c) 1998-2010
+    // Distributed under the Boost Software License, Version 1.0.
+    // http://www.boost.org/LICENSE_1_0.txt
+    // http://www.geometrictools.com/License/Boost/LICENSE_1_0.txt
+
+    // Boost Software License - Version 1.0 - August 17th, 2003
+
+    // Permission is hereby granted, free of charge, to any person or
+    // organization obtaining a copy of the software and accompanying
+    // documentation covered by this license (the "Software") to use,
+    // reproduce, display, distribute, execute, and transmit the
+    // Software, and to prepare derivative works of the Software, and
+    // to permit third-parties to whom the Software is furnished to do
+    // so, all subject to the following:
+
+    // The copyright notices in the Software and this entire
+    // statement, including the above license grant, this restriction
+    // and the following disclaimer, must be included in all copies of
+    // the Software, in whole or in part, and all derivative works of
+    // the Software, unless such copies or derivative works are solely
+    // in the form of machine-executable object code generated by a
+    // source language processor.
+
+    // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+    // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+    // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE AND
+    // NON-INFRINGEMENT. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR
+    // ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE FOR ANY DAMAGES OR
+    // OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
+    // ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
+    // USE OR OTHER DEALINGS IN THE SOFTWARE.
+
+    const scalar r6 = 1.0/6.0;
+    const scalar r24 = 1.0/24.0;
+    const scalar r60 = 1.0/60.0;
+    const scalar r120 = 1.0/120.0;
+
+    // order:  1, x, y, z, x^2, y^2, z^2, xy, yz, zx
+    scalarField integrals(10, 0.0);
+
+    forAll(triFaces, i)
     {
-        label faceI = cFaces[i];
+        const triFace& tri(triFaces[i]);
 
-        const face& f = faces[faceI];
+        // vertices of triangle i
+        vector v0 = pts[tri[0]];
+        vector v1 = pts[tri[1]];
+        vector v2 = pts[tri[2]];
 
-        fNorm[i] = f.normal(points);
-        fNorm[i] /= mag(fNorm[i]) + VSMALL;
+        // cross product of edges
+        vector eA = v1 - v0;
+        vector eB = v2 - v0;
+        vector n = eA ^ eB;
 
-        fW[i] = - (fNorm[i] & points[f[0]]);
+        // compute integral terms
+        scalar tmp0, tmp1, tmp2;
+
+        scalar f1x, f2x, f3x, g0x, g1x, g2x;
+
+        tmp0 = v0.x() + v1.x();
+        f1x = tmp0 + v2.x();
+        tmp1 = v0.x()*v0.x();
+        tmp2 = tmp1 + v1.x()*tmp0;
+        f2x = tmp2 + v2.x()*f1x;
+        f3x = v0.x()*tmp1 + v1.x()*tmp2 + v2.x()*f2x;
+        g0x = f2x + v0.x()*(f1x + v0.x());
+        g1x = f2x + v1.x()*(f1x + v1.x());
+        g2x = f2x + v2.x()*(f1x + v2.x());
+
+        scalar f1y, f2y, f3y, g0y, g1y, g2y;
+
+        tmp0 = v0.y() + v1.y();
+        f1y = tmp0 + v2.y();
+        tmp1 = v0.y()*v0.y();
+        tmp2 = tmp1 + v1.y()*tmp0;
+        f2y = tmp2 + v2.y()*f1y;
+        f3y = v0.y()*tmp1 + v1.y()*tmp2 + v2.y()*f2y;
+        g0y = f2y + v0.y()*(f1y + v0.y());
+        g1y = f2y + v1.y()*(f1y + v1.y());
+        g2y = f2y + v2.y()*(f1y + v2.y());
+
+        scalar f1z, f2z, f3z, g0z, g1z, g2z;
+
+        tmp0 = v0.z() + v1.z();
+        f1z = tmp0 + v2.z();
+        tmp1 = v0.z()*v0.z();
+        tmp2 = tmp1 + v1.z()*tmp0;
+        f2z = tmp2 + v2.z()*f1z;
+        f3z = v0.z()*tmp1 + v1.z()*tmp2 + v2.z()*f2z;
+        g0z = f2z + v0.z()*(f1z + v0.z());
+        g1z = f2z + v1.z()*(f1z + v1.z());
+        g2z = f2z + v2.z()*(f1z + v2.z());
+
+        // update integrals
+        integrals[0] += n.x()*f1x;
+        integrals[1] += n.x()*f2x;
+        integrals[2] += n.y()*f2y;
+        integrals[3] += n.z()*f2z;
+        integrals[4] += n.x()*f3x;
+        integrals[5] += n.y()*f3y;
+        integrals[6] += n.z()*f3z;
+        integrals[7] += n.x()*(v0.y()*g0x + v1.y()*g1x + v2.y()*g2x);
+        integrals[8] += n.y()*(v0.z()*g0y + v1.z()*g1y + v2.z()*g2y);
+        integrals[9] += n.z()*(v0.x()*g0z + v1.x()*g1z + v2.x()*g2z);
     }
 
+    integrals[0] *= r6;
+    integrals[1] *= r24;
+    integrals[2] *= r24;
+    integrals[3] *= r24;
+    integrals[4] *= r60;
+    integrals[5] *= r60;
+    integrals[6] *= r60;
+    integrals[7] *= r120;
+    integrals[8] *= r120;
+    integrals[9] *= r120;
 
-    scalar T0;
-    vector T1, T2, TP;
+    // mass
+    mass = integrals[0];
 
-    compVolumeIntegrals
-    (
-        points,
-        faces,
-        cFaces,
-        fNorm,
-        fW,
+    // center of mass
+    cM = vector(integrals[1], integrals[2], integrals[3])/mass;
 
-        T0,
-        T1,
-        T2,
-        TP
-    );
+    // inertia relative to origin
+    J.xx() = integrals[5] + integrals[6];
+    J.xy() = -integrals[7];
+    J.xz() = -integrals[9];
+    J.yx() = J.xy();
+    J.yy() = integrals[4] + integrals[6];
+    J.yz() = -integrals[8];
+    J.zx() = J.xz();
+    J.zy() = J.yz();
+    J.zz() = integrals[4] + integrals[5];
 
-    const scalar density = 1.0;  /* assume unit density */
+    // inertia relative to center of mass
+    J -= mass*((cM & cM)*I - cM*cM);
 
-    scalar mass = density * T0;
-
-    /* compute center of mass */
-    r = T1 / T0;
-
-    /* compute inertia tensor */
-    J.xx() = density * (T2[1] + T2[2]);
-    J.yy() = density * (T2[2] + T2[0]);
-    J.zz() = density * (T2[0] + T2[1]);
-    J.xy() = J.yx() = - density * TP[0];
-    J.yz() = J.zy() = - density * TP[1];
-    J.zx() = J.xz() = - density * TP[2];
-
-    ///* translate inertia tensor to center of mass */
-    //J[XX] -= mass * (r[1]*r[1] + r[2]*r[2]);
-    //J[YY] -= mass * (r[2]*r[2] + r[0]*r[0]);
-    //J[ZZ] -= mass * (r[0]*r[0] + r[1]*r[1]);
-    //J[XY] = J[YX] += mass * r[0] * r[1];
-    //J[YZ] = J[ZY] += mass * r[1] * r[2];
-    //J[ZX] = J[XZ] += mass * r[2] * r[0];
+    // Apply density
+    mass *= density;
+    J *= density;
 }
 
 
+void Foam::momentOfInertia::massPropertiesShell
+(
+    const pointField& pts,
+    const triFaceList& triFaces,
+    scalar density,
+    scalar& mass,
+    vector& cM,
+    tensor& J
+)
+{
+    // Reset properties for accumulation
+
+    mass = 0.0;
+    cM = vector::zero;
+    J = tensor::zero;
+
+    // Find centre of mass
+
+    forAll(triFaces, i)
+    {
+        const triFace& tri(triFaces[i]);
+
+        triPointRef t
+        (
+            pts[tri[0]],
+            pts[tri[1]],
+            pts[tri[2]]
+        );
+
+        scalar triMag = t.mag();
+
+        cM +=  triMag*t.centre();
+
+        mass += triMag;
+    }
+
+    cM /= mass;
+
+    mass *= density;
+
+    // Find inertia around centre of mass
+
+    forAll(triFaces, i)
+    {
+        const triFace& tri(triFaces[i]);
+
+        J += triPointRef
+        (
+            pts[tri[0]],
+            pts[tri[1]],
+            pts[tri[2]]
+        ).inertia(cM, density);
+    }
+}
+
+
+void Foam::momentOfInertia::massPropertiesSolid
+(
+    const triSurface& surf,
+    scalar density,
+    scalar& mass,
+    vector& cM,
+    tensor& J
+)
+{
+    triFaceList faces(surf.size());
+
+    forAll(surf, i)
+    {
+        faces[i] = triFace(surf[i]);
+    }
+
+    massPropertiesSolid(surf.points(), faces, density, mass, cM, J);
+}
+
+
+void Foam::momentOfInertia::massPropertiesShell
+(
+    const triSurface& surf,
+    scalar density,
+    scalar& mass,
+    vector& cM,
+    tensor& J
+)
+{
+    triFaceList faces(surf.size());
+
+    forAll(surf, i)
+    {
+        faces[i] = triFace(surf[i]);
+    }
+
+    massPropertiesShell(surf.points(), faces, density, mass, cM, J);
+}
+
+
+Foam::tensor Foam::momentOfInertia::applyParallelAxisTheorem
+(
+    scalar mass,
+    const vector& cM,
+    const tensor& J,
+    const vector& refPt
+)
+{
+    // The displacement vector (refPt = cM) is the displacement of the
+    // new reference point from the centre of mass of the body that
+    // the inertia tensor applies to.
+
+    vector d = (refPt - cM);
+
+    return J + mass*((d & d)*I - d*d);
+}
+
+
+Foam::tmp<Foam::tensorField> Foam::momentOfInertia::meshInertia
+(
+    const polyMesh& mesh
+)
+{
+    tmp<tensorField> tTf = tmp<tensorField>(new tensorField(mesh.nCells()));
+
+    tensorField& tf = tTf();
+
+    forAll(tf, cI)
+    {
+        tf[cI] = meshInertia(mesh, cI);
+    }
+
+    return tTf;
+}
+
+
+Foam::tensor Foam::momentOfInertia::meshInertia
+(
+    const polyMesh& mesh,
+    label cellI
+)
+{
+    List<tetIndices> cellTets = polyMeshTetDecomposition::cellTetIndices
+    (
+        mesh,
+        cellI
+    );
+
+    triFaceList faces(cellTets.size());
+
+    forAll(cellTets, cTI)
+    {
+        faces[cTI] = cellTets[cTI].faceTriIs(mesh);
+    }
+
+    scalar m = 0.0;
+    vector cM = vector::zero;
+    tensor J = tensor::zero;
+
+    massPropertiesSolid(mesh.points(), faces, 1.0, m, cM, J);
+
+    return J;
+}
+
 
 // ************************************************************************* //

src/meshTools/momentOfInertia/momentOfInertia.H

@@ -25,6 +25,9 @@ Class
     momentOfInertia
 
 Description
+    Calculates the inertia tensor and principal axes and moments of a
+    polyhedra/cells/triSurfaces.  Inertia can either be of the solid body or
+    of a thin shell.
 
 SourceFiles
     momentOfInertia.H
@@ -34,34 +37,86 @@ SourceFiles
 #ifndef momentOfInertia_H
 #define momentOfInertia_H
 
-#include "tensor.H"
-#include "primitiveMesh.H"
+#include "tetPointRef.H"
+#include "triFaceList.H"
+#include "triSurface.H"
+#include "polyMesh.H"
+#include "polyMeshTetDecomposition.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
-////- Moment of inertia around cell centre for single cell.
-//tensor momentOfInertia
-//(
-//    const pointField&,
-//    const faceList&,
-//    const cell&,
-//    const point& cc
-//);
+/*---------------------------------------------------------------------------*\
+                     Class momentOfInertia Declaration
+\*---------------------------------------------------------------------------*/
+
+class momentOfInertia
+{
+
+public:
+
+    static void massPropertiesSolid
+    (
+        const pointField& pts,
+        const triFaceList& triFaces,
+        scalar density,
+        scalar& mass,
+        vector& cM,
+        tensor& J
+    );
+
+    static void massPropertiesShell
+    (
+        const pointField& pts,
+        const triFaceList& triFaces,
+        scalar density,
+        scalar& mass,
+        vector& cM,
+        tensor& J
+    );
+
+    static void massPropertiesSolid
+    (
+        const triSurface& surf,
+        scalar density,
+        scalar& mass,
+        vector& cM,
+        tensor& J
+    );
+
+    static void massPropertiesShell
+    (
+        const triSurface& surf,
+        scalar density,
+        scalar& mass,
+        vector& cM,
+        tensor& J
+    );
+
+    static tensor applyParallelAxisTheorem
+    (
+        scalar mass,
+        const vector& cM,
+        const tensor& J,
+        const vector& refPt
+    );
+
+    // Calculate the inertia tensor for all cells in the mesh
+    static tmp<tensorField> meshInertia
+    (
+        const polyMesh& mesh
+    );
+
+    // Calculate the inertia tensor the given cell
+    static tensor meshInertia
+    (
+        const polyMesh& mesh,
+        label cellI
+    );
+};
 
-// Calculate
-// - centre of mass
-// - inertia tensor around (0,0,0)
-void momentOfIntertia
-(
-    const pointField&,
-    const faceList&,
-    const cell&,
-    point& r,
-    tensor& Jorigin
-);
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/meshTools/octree/treeBoundBox.C

@@ -165,15 +165,19 @@ Foam::treeBoundBox::treeBoundBox
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::pointField Foam::treeBoundBox::points() const
+Foam::tmp<Foam::pointField> Foam::treeBoundBox::points() const
 {
-    pointField points(8);
+    tmp<pointField> tPts = tmp<pointField>(new pointField(8));
+
+    pointField& points = tPts();
+
     forAll(points, octant)
     {
         points[octant] = corner(octant);
 
     }
-    return points;
+
+    return tPts;
 }
 
 

src/meshTools/octree/treeBoundBox.H

@@ -203,7 +203,7 @@ public:
             inline scalar typDim() const;
 
             //- vertex coordinates. In octant coding.
-            pointField points() const;
+            tmp<pointField> points() const;
 
 
         // Check

src/parallel/decompose/decompositionMethods/structuredDecomp/topoDistanceData.H

@@ -195,6 +195,14 @@ public:
 };
 
 
+//- Data associated with topoDistanceData type are contiguous
+template<>
+inline bool contiguous<topoDistanceData>()
+{
+    return true;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

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

@@ -45,6 +45,7 @@ Foam::scalar Foam::streamLineParticle::calcSubCycleDeltaT
 ) const
 {
     streamLineParticle testParticle(*this);
+
     bool oldKeepParticle = td.keepParticle;
     bool oldSwitchProcessor = td.switchProcessor;
     scalar fraction = testParticle.trackToFace(position()+dt*U, td);
@@ -196,9 +197,7 @@ bool Foam::streamLineParticle::move
         // Cross cell in steps:
         // - at subiter 0 calculate dt to cross cell in nSubCycle steps
         // - at the last subiter do all of the remaining track
-        // - do a few more subiters than nSubCycle since velocity might
-        //   be decreasing
-        for (label subIter = 0; subIter < 2*td.nSubCycle_; subIter++)
+        for (label subIter = 0; subIter < td.nSubCycle_; subIter++)
         {
             --lifeTime_;
 
@@ -217,14 +216,14 @@ bool Foam::streamLineParticle::move
                 // Adapt dt to cross cell in a few steps
                 dt = calcSubCycleDeltaT(td, dt, U);
             }
-            else if (subIter == td.nSubCycle_-1)
+            else if (subIter == td.nSubCycle_ - 1)
             {
                 // Do full step on last subcycle
                 dt = min(dtMax, tEnd);
             }
 
 
-            scalar fraction = trackToFace(position()+dt*U, td);
+            scalar fraction = trackToFace(position() + dt*U, td);
             dt *= fraction;
 
             tEnd -= dt;
@@ -236,7 +235,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/sampling/Make/files

@@ -1,4 +1,5 @@
 probes/probes.C
+probes/patchProbes.C
 probes/probesGrouping.C
 probes/probesFunctionObject/probesFunctionObject.C
 

src/sampling/probes/IOprobes.H

@@ -33,6 +33,7 @@ Description
 #define IOprobes_H
 
 #include "probes.H"
+#include "patchProbes.H"
 #include "IOOutputFilter.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -40,6 +41,7 @@ Description
 namespace Foam
 {
     typedef IOOutputFilter<probes> IOprobes;
+    typedef IOOutputFilter<patchProbes> IOpatchProbes;
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/sampling/probes/patchProbes.C

@@ -0,0 +1,160 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "patchProbes.H"
+#include "volFields.H"
+#include "IOmanip.H"
+// For 'nearInfo' helper class only
+#include "directMappedPatchBase.H"
+#include "meshSearch.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(patchProbes, 0);
+}
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::patchProbes::findElements(const fvMesh& mesh)
+{
+    elementList_.clear();
+    elementList_.setSize(size());
+     // All the info for nearest. Construct to miss
+    List<directMappedPatchBase::nearInfo> nearest(this->size());
+
+    // Octree based search engine
+    meshSearch meshSearchEngine(mesh, false);
+
+    forAll(*this, probeI)
+    {
+        const vector& sample = operator[](probeI);
+        label faceI = meshSearchEngine.findNearestBoundaryFace(sample);
+        if (faceI == -1)
+        {
+            nearest[probeI].second().first() = Foam::sqr(GREAT);
+            nearest[probeI].second().second() = Pstream::myProcNo();
+        }
+        else
+        {
+            const point& fc = mesh.faceCentres()[faceI];
+            nearest[probeI].first() = pointIndexHit
+            (
+                true,
+                fc,
+                faceI
+            );
+            nearest[probeI].second().first() = magSqr(fc-sample);
+            nearest[probeI].second().second() = Pstream::myProcNo();
+        }
+    }
+
+
+    // Find nearest.
+    Pstream::listCombineGather(nearest, directMappedPatchBase::nearestEqOp());
+    Pstream::listCombineScatter(nearest);
+
+    if (debug)
+    {
+        Info<< "patchProbes::findElements" << " : " << endl;
+        forAll(nearest, sampleI)
+        {
+            label procI = nearest[sampleI].second().second();
+            label localI = nearest[sampleI].first().index();
+
+            Info<< "    " << sampleI << " coord:"<< operator[](sampleI)
+                << " found on processor:" << procI
+                << " in local cell/face:" << localI
+                << " with cc:" << nearest[sampleI].first().rawPoint() << endl;
+        }
+    }
+
+
+
+    // Check if all patchProbes have been found.
+    forAll(nearest, sampleI)
+    {
+        label localI = nearest[sampleI].first().index();
+
+        if (localI == -1)
+        {
+             if (Pstream::master())
+             {
+                WarningIn("patchProbes::findElements()")
+                    << "Did not find location "
+                    <<  nearest[sampleI].second().first()
+                    << " in any cell. Skipping location." << endl;
+             }
+        }
+        else
+        {
+            elementList_[sampleI] = localI;
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::patchProbes::patchProbes
+(
+    const word& name,
+    const objectRegistry& obr,
+    const dictionary& dict,
+    const bool loadFromFiles
+)
+:
+    probes(name, obr, dict, loadFromFiles)
+{
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::patchProbes::~patchProbes()
+{}
+
+
+void Foam::patchProbes::write()
+{
+    if (this->size() && prepare())
+    {
+        sampleAndWrite(scalarFields_);
+        sampleAndWrite(vectorFields_);
+        sampleAndWrite(sphericalTensorFields_);
+        sampleAndWrite(symmTensorFields_);
+        sampleAndWrite(tensorFields_);
+    }
+}
+
+void Foam::patchProbes::read(const dictionary& dict)
+{
+    probes::read(dict);
+}
+
+
+// ************************************************************************* //

src/sampling/probes/patchProbes.H

@@ -0,0 +1,147 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::patchProbes
+
+Description
+    Set of locations to sample.at patches
+
+    Call write() to sample and write files.
+
+SourceFiles
+    patchProbes.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef patchProbes_H
+#define patchProbes_H
+
+#include "probes.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class objectRegistry;
+class dictionary;
+class fvMesh;
+class mapPolyMesh;
+
+/*---------------------------------------------------------------------------*\
+                          Class patchProbes Declaration
+\*---------------------------------------------------------------------------*/
+
+class patchProbes
+:
+    public probes
+{
+
+    // Private Member Functions
+
+        //- Sample and write a particular volume field
+        template<class Type>
+        void sampleAndWrite
+        (
+            const GeometricField<Type, fvPatchField, volMesh>&
+        );
+
+
+        //- Sample and write all the fields of the given type
+        template <class Type>
+        void sampleAndWrite(const fieldGroup<Type>&);
+
+
+        //- Sample a volume field at all locations
+        template<class Type>
+        tmp<Field<Type> > sample
+        (
+            const GeometricField<Type, fvPatchField, volMesh>&
+        ) const;
+
+
+        //- Sample a single field on all sample locations
+        template <class Type>
+        tmp<Field<Type> > sample(const word& fieldName) const;
+
+
+        //- Disallow default bitwise copy construct
+        patchProbes(const patchProbes&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const patchProbes&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("patchProbes");
+
+
+    // Constructors
+
+        //- Construct for given objectRegistry and dictionary.
+        //  Allow the possibility to load fields from files
+        patchProbes
+        (
+            const word& name,
+            const objectRegistry&,
+            const dictionary&,
+            const bool loadFromFiles = false
+        );
+
+
+    //- Destructor
+    virtual ~patchProbes();
+
+    //- Public members
+
+        //- Sample and write
+        virtual void write();
+
+        //- Read
+        virtual void read(const dictionary&);
+
+        //- Find elements containing patchProbes
+        virtual void findElements(const fvMesh&);
+
+
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "patchProbesTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/sampling/probes/patchProbesTemplates.C

@@ -0,0 +1,162 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "patchProbes.H"
+#include "volFields.H"
+#include "IOmanip.H"
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+void Foam::patchProbes::sampleAndWrite
+(
+    const GeometricField<Type, fvPatchField, volMesh>& vField
+)
+{
+    Field<Type> values = sample(vField);
+
+    if (Pstream::master())
+    {
+        unsigned int w = IOstream::defaultPrecision() + 7;
+        OFstream& probeStream = *probeFilePtrs_[vField.name()];
+
+        probeStream << setw(w) << vField.time().value();
+
+        forAll(values, probeI)
+        {
+            probeStream << ' ' << setw(w) << values[probeI];
+        }
+        probeStream << endl;
+    }
+}
+
+
+template <class Type>
+void Foam::patchProbes::sampleAndWrite
+(
+    const fieldGroup<Type>& fields
+)
+{
+    forAll(fields, fieldI)
+    {
+        if (loadFromFiles_)
+        {
+            sampleAndWrite
+            (
+                GeometricField<Type, fvPatchField, volMesh>
+                (
+                    IOobject
+                    (
+                        fields[fieldI],
+                        mesh_.time().timeName(),
+                        mesh_,
+                        IOobject::MUST_READ,
+                        IOobject::NO_WRITE,
+                        false
+                    ),
+                    mesh_
+                )
+            );
+        }
+        else
+        {
+            objectRegistry::const_iterator iter = mesh_.find(fields[fieldI]);
+
+            if
+            (
+                iter != objectRegistry::end()
+             && iter()->type()
+             == GeometricField<Type, fvPatchField, volMesh>::typeName
+            )
+            {
+                sampleAndWrite
+                (
+                    mesh_.lookupObject
+                    <GeometricField<Type, fvPatchField, volMesh> >
+                    (
+                        fields[fieldI]
+                    )
+                );
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::tmp<Foam::Field<Type> >
+Foam::patchProbes::sample
+(
+    const GeometricField<Type, fvPatchField, volMesh>& vField
+) const
+{
+    const Type unsetVal(-VGREAT*pTraits<Type>::one);
+
+    tmp<Field<Type> > tValues
+    (
+        new Field<Type>(this->size(), unsetVal)
+    );
+
+    Field<Type>& values = tValues();
+
+    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+
+    forAll(*this, probeI)
+    {
+        label faceI = elementList_[probeI];
+
+        if (faceI >= 0)
+        {
+            label patchI = patches.whichPatch(faceI);
+            label localFaceI = patches[patchI].whichFace(faceI);
+            values[probeI] = vField.boundaryField()[patchI][localFaceI];
+        }
+    }
+
+    Pstream::listCombineGather(values, isNotEqOp<Type>());
+    Pstream::listCombineScatter(values);
+
+    return tValues;
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type> >
+Foam::patchProbes::sample(const word& fieldName) const
+{
+    return sample
+    (
+        mesh_.lookupObject<GeometricField<Type, fvPatchField, volMesh> >
+        (
+            fieldName
+        )
+    );
+}
+
+
+// ************************************************************************* //

src/sampling/probes/probes.C

@@ -36,30 +36,30 @@ defineTypeNameAndDebug(Foam::probes, 0);
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::probes::findCells(const fvMesh& mesh)
+void Foam::probes::findElements(const fvMesh& mesh)
 {
-    cellList_.clear();
-    cellList_.setSize(size());
+    elementList_.clear();
+    elementList_.setSize(size());
 
     forAll(*this, probeI)
     {
         const vector& location = operator[](probeI);
 
-        cellList_[probeI] = mesh.findCell(location);
+        elementList_[probeI] = mesh.findCell(location);
 
-        if (debug && cellList_[probeI] != -1)
+        if (debug && elementList_[probeI] != -1)
         {
             Pout<< "probes : found point " << location
-                << " in cell " << cellList_[probeI] << endl;
+                << " in cell " << elementList_[probeI] << endl;
         }
     }
 
 
     // Check if all probes have been found.
-    forAll(cellList_, probeI)
+    forAll(elementList_, probeI)
     {
         const vector& location = operator[](probeI);
-        label cellI = cellList_[probeI];
+        label cellI = elementList_[probeI];
 
         // Check at least one processor with cell.
         reduce(cellI, maxOp<label>());
@@ -76,12 +76,12 @@ void Foam::probes::findCells(const fvMesh& mesh)
         else
         {
             // Make sure location not on two domains.
-            if (cellList_[probeI] != -1 && cellList_[probeI] != cellI)
+            if (elementList_[probeI] != -1 && elementList_[probeI] != cellI)
             {
                 WarningIn("probes::read()")
                     << "Location " << location
                     << " seems to be on multiple domains:"
-                    << " cell " << cellList_[probeI]
+                    << " cell " << elementList_[probeI]
                     << " on my domain " << Pstream::myProcNo()
                         << " and cell " << cellI << " on some other domain."
                     << endl
@@ -249,7 +249,7 @@ void Foam::probes::read(const dictionary& dict)
     dict.lookup("fields") >> fieldSelection_;
 
     // redetermined all cell locations
-    findCells(mesh_);
+    findElements(mesh_);
     prepare();
 }