Merge branch 'master' into AMI

bin/tools/pre-commit-hook

@@ -273,7 +273,7 @@ checkCopyright()
 
     for f in $fileList
     do
-        sYear=`grep "Copyright.*OpenFOAM" $f | sed 's/[^0-9]//g' | cut -c 5-9`
+        sYear=`grep "Copyright.*OpenFOAM" $f | sed 's/[^0-9]//g'`
         if [ "$year" != "" ] && [ "$year" != "$sYear" ]; then
             echo "Updated copyright for: $f"
             sed -i "s/$sYear OpenFOAM/$year OpenFOAM/g" $f

src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchFieldNew.C

@@ -1,25 +1,25 @@
 /*---------------------------------------------------------------------------*\
   =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+  \\      /  F ield         |2011 OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
-    This file is part of OpenFOAM.
+    This file is part of2011 OpenFOAM.
 
-    OpenFOAM is free software: you can redistribute it and/or modify it
+   2011 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
+   2011 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/>.
+    along with2011 OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 \*---------------------------------------------------------------------------*/
 
@@ -169,7 +169,7 @@ Foam::autoPtr<Foam::pointPatchField<Type> > Foam::pointPatchField<Type>::New
             typename dictionaryConstructorTable::iterator patchTypeCstrIter
                 = dictionaryConstructorTablePtr_->find(p.type());
 
-            if (patchTypeCstrIter == pointPatchConstructorTablePtr_->end())
+            if (patchTypeCstrIter == dictionaryConstructorTablePtr_->end())
             {
                 FatalIOErrorIn
                 (

src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.C

@@ -375,7 +375,7 @@ Foam::primitiveMesh::cellTree() const
                 overallBb,
                 8,              // maxLevel
                 10,             // leafsize
-                3.0             // duplicity
+                5.0             // duplicity
             );
     }
 

src/lagrangian/basic/InteractionLists/InteractionLists.C

@@ -167,7 +167,7 @@ void Foam::InteractionLists<ParticleType>::buildInteractionLists()
         procBbRndExt,
         8,              // maxLevel,
         10,             // leafSize,
-        100.0
+        100.0           // duplicity
     );
 
     ril_.setSize(cellBbsToExchange.size());
@@ -386,7 +386,7 @@ void Foam::InteractionLists<ParticleType>::buildInteractionLists()
         procBbRndExt,
         8,              // maxLevel,
         10,             // leafSize,
-        100.0
+        100.0           // duplicity
     );
 
     rwfil_.setSize(wallFaceBbsToExchange.size());

src/meshTools/meshSearch/meshSearch.C

@@ -30,7 +30,6 @@ License
 #include "demandDrivenData.H"
 #include "treeDataCell.H"
 #include "treeDataFace.H"
-#include "treeDataPoint.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -508,6 +507,21 @@ Foam::meshSearch::meshSearch(const polyMesh& mesh, const bool faceDecomp)
 {}
 
 
+// Construct with a custom bounding box
+Foam::meshSearch::meshSearch
+(
+    const polyMesh& mesh,
+    const treeBoundBox& bb,
+    const bool faceDecomp
+)
+:
+    mesh_(mesh),
+    faceDecomp_(faceDecomp)
+{
+    overallBbPtr_.reset(new treeBoundBox(bb));
+}
+
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::meshSearch::~meshSearch()
@@ -527,6 +541,21 @@ const Foam::indexedOctree<Foam::treeDataFace>& Foam::meshSearch::boundaryTree()
         // Construct tree
         //
 
+        if (!overallBbPtr_.valid())
+        {
+            Random rndGen(261782);
+            overallBbPtr_.reset
+            (
+                new treeBoundBox(mesh_.points())
+            );
+
+            treeBoundBox& overallBb = overallBbPtr_();
+            // Extend slightly and make 3D
+            overallBb = overallBb.extend(rndGen, 1E-4);
+            overallBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+            overallBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+        }
+
         // all boundary faces (not just walls)
         labelList bndFaces(mesh_.nFaces()-mesh_.nInternalFaces());
         forAll(bndFaces, i)
@@ -534,12 +563,6 @@ const Foam::indexedOctree<Foam::treeDataFace>& Foam::meshSearch::boundaryTree()
             bndFaces[i] = mesh_.nInternalFaces() + i;
         }
 
-        treeBoundBox overallBb(mesh_.points());
-        Random rndGen(123456);
-        overallBb = overallBb.extend(rndGen, 1E-4);
-        overallBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
-        overallBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
-
         boundaryTreePtr_.reset
         (
             new indexedOctree<treeDataFace>
@@ -550,7 +573,7 @@ const Foam::indexedOctree<Foam::treeDataFace>& Foam::meshSearch::boundaryTree()
                     mesh_,
                     bndFaces                    // boundary faces only
                 ),
-                overallBb,                      // overall search domain
+                overallBbPtr_(),                // overall search domain
                 8,                              // maxLevel
                 10,                             // leafsize
                 3.0                             // duplicity
@@ -567,13 +590,24 @@ const
 {
     if (!cellTreePtr_.valid())
     {
-        treeBoundBox overallBb(mesh_.points());
+        //
+        // Construct tree
+        //
 
+        if (!overallBbPtr_.valid())
+        {
             Random rndGen(261782);
+            overallBbPtr_.reset
+            (
+                new treeBoundBox(mesh_.points())
+            );
 
+            treeBoundBox& overallBb = overallBbPtr_();
+            // Extend slightly and make 3D
             overallBb = overallBb.extend(rndGen, 1E-4);
             overallBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
             overallBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+        }
 
         cellTreePtr_.reset
         (
@@ -584,10 +618,10 @@ const
                     false,  // not cache bb
                     mesh_
                 ),
-                overallBb,
+                overallBbPtr_(),
                 8,              // maxLevel
                 10,             // leafsize
-                3.0             // duplicity
+                6.0             // duplicity
             )
         );
     }
@@ -904,6 +938,7 @@ void Foam::meshSearch::clearOut()
 {
     boundaryTreePtr_.clear();
     cellTreePtr_.clear();
+    overallBbPtr_.clear();
 }
 
 

src/meshTools/meshSearch/meshSearch.H

@@ -49,6 +49,7 @@ class polyMesh;
 class treeDataCell;
 class treeDataFace;
 template<class Type> class indexedOctree;
+class treeBoundBox;
 
 /*---------------------------------------------------------------------------*\
                            Class meshSearch Declaration
@@ -64,8 +65,10 @@ class meshSearch
         //- Whether to use face decomposition for all geometric tests
         const bool faceDecomp_;
 
-        //- demand driven octrees
+        //- data bounding box
+        mutable autoPtr<treeBoundBox> overallBbPtr_;
 
+        //- demand driven octrees
         mutable autoPtr<indexedOctree<treeDataFace> > boundaryTreePtr_;
         mutable autoPtr<indexedOctree<treeDataCell> > cellTreePtr_;
 
@@ -163,9 +166,19 @@ public:
 
     // Constructors
 
-        //- Construct from components
+        //- Construct from components. Constructs bb slightly bigger than
+        //  mesh points bb.
         meshSearch(const polyMesh& mesh, const bool faceDecomp = true);
 
+        //- Construct with a custom bounding box. Any mesh element outside
+        //  bb will not be found. Up to user to make sure bb
+        //  extends slightly beyond wanted elements.
+        meshSearch
+        (
+            const polyMesh& mesh,
+            const treeBoundBox& bb,
+            const bool faceDecomp = true
+        );
 
     //- Destructor
     ~meshSearch();

src/postProcessing/functionObjects/field/fieldAverage/controlDict

@@ -68,6 +68,7 @@ functions
                 mean        on;
                 prime2Mean  on;
                 base        time;
+                window      0.01; // optional averaging window
             }
 
             p
@@ -75,6 +76,7 @@ functions
                 mean        on;
                 prime2Mean  on;
                 base        time;
+                window      0.01; // optional averaging window
             }
         );
     }

src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverage.C

@@ -377,8 +377,7 @@ void Foam::fieldAverage::execute()
 
 
 void Foam::fieldAverage::end()
-{
+{}
-}
 
 
 void Foam::fieldAverage::write()

src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverageTemplates.C

@@ -149,7 +149,7 @@ const
 {
     typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
 
-    const scalar dt = obr_.time().deltaTValue();
+    scalar dt = obr_.time().deltaTValue();
 
     forAll(faItems_, i)
     {
@@ -163,17 +163,24 @@ const
                 obr_.lookupObject<fieldType>(meanFieldList[i])
             );
 
-            scalar alpha = 0.0;
+            scalar Dt = totalTime_[i];
-            scalar beta = 0.0;
+            if (faItems_[i].iterBase())
-            if (faItems_[i].timeBase())
             {
-                 alpha = (totalTime_[i] - dt)/totalTime_[i];
+                dt = 1.0;
-                 beta = dt/totalTime_[i];
+                Dt = scalar(totalIter_[i]);
             }
-            else
+
+            scalar alpha = (Dt - dt)/Dt;
+            scalar beta = dt/Dt;
+            if (faItems_[i].window() > 0)
             {
-                alpha = scalar(totalIter_[i] - 1)/scalar(totalIter_[i]);
+                const scalar w = faItems_[i].window();
-                beta = 1.0/scalar(totalIter_[i]);
+
+                if (Dt - dt >= w)
+                {
+                    alpha = (w - dt)/w;
+                    beta = dt/w;
+                }
             }
 
             meanField = alpha*meanField + beta*baseField;
@@ -192,7 +199,7 @@ void Foam::fieldAverage::calculatePrime2MeanFields
     typedef GeometricField<Type1, fvPatchField, volMesh> fieldType1;
     typedef GeometricField<Type2, fvPatchField, volMesh> fieldType2;
 
-    const scalar dt = obr_.time().deltaTValue();
+    scalar dt = obr_.time().deltaTValue();
 
     forAll(faItems_, i)
     {
@@ -213,17 +220,24 @@ void Foam::fieldAverage::calculatePrime2MeanFields
                 obr_.lookupObject<fieldType2>(prime2MeanFieldList[i])
             );
 
-            scalar alpha = 0.0;
+            scalar Dt = totalTime_[i];
-            scalar beta = 0.0;
+            if (faItems_[i].iterBase())
-            if (faItems_[i].timeBase())
             {
-                alpha = (totalTime_[i] - dt)/totalTime_[i];
+                dt = 1.0;
-                beta = dt/totalTime_[i];
+                Dt = scalar(totalIter_[i]);
             }
-            else
+
+            scalar alpha = (Dt - dt)/Dt;
+            scalar beta = dt/Dt;
+            if (faItems_[i].window() > 0)
             {
-                alpha = scalar(totalIter_[i] - 1)/scalar(totalIter_[i]);
+                const scalar w = faItems_[i].window();
-                beta = 1.0/scalar(totalIter_[i]);
+
+                if (Dt - dt >= w)
+                {
+                    alpha = (w - dt)/w;
+                    beta = dt/w;
+                }
             }
 
             prime2MeanField =

src/postProcessing/functionObjects/field/fieldAverage/fieldAverageItem/fieldAverageItem.C

@@ -53,7 +53,8 @@ Foam::fieldAverageItem::fieldAverageItem()
     fieldName_("unknown"),
     mean_(0),
     prime2Mean_(0),
-    base_(ITER)
+    base_(ITER),
+    window_(-1.0)
 {}
 
 
@@ -62,7 +63,8 @@ Foam::fieldAverageItem::fieldAverageItem(const fieldAverageItem& faItem)
     fieldName_(faItem.fieldName_),
     mean_(faItem.mean_),
     prime2Mean_(faItem.prime2Mean_),
-    base_(faItem.base_)
+    base_(faItem.base_),
+    window_(faItem.window_)
 {}
 
 
@@ -91,6 +93,7 @@ void Foam::fieldAverageItem::operator=(const fieldAverageItem& rhs)
     mean_ = rhs.mean_;
     prime2Mean_ = rhs.prime2Mean_;
     base_ = rhs.base_;
+    window_ = rhs.window_;
 }
 
 

src/postProcessing/functionObjects/field/fieldAverage/fieldAverageItem/fieldAverageItem.H

@@ -33,9 +33,13 @@ Description
             mean            on;
             prime2Mean      on;
             base            time; // iteration
+            window          200;  // optional averaging window
         }
     \endverbatim
 
+    The averaging window corresponds to the averaging interval (iters or time)
+    If not specified, the averaging is over 'all iters/time'
+
 SourceFiles
     fieldAverageItem.C
     fieldAverageItemIO.C
@@ -100,6 +104,9 @@ private:
         //- Averaging base type
         baseType base_;
 
+        //- Averaging window - defaults to -1 for 'all iters/time'
+        scalar window_;
+
 
 public:
 
@@ -148,7 +155,7 @@ public:
             }
 
             //- Return true if base is ITER
-            Switch ITERBase() const
+            Switch iterBase() const
             {
                 return base_ == ITER;
             }
@@ -159,6 +166,11 @@ public:
                 return base_ == TIME;
             }
 
+            scalar window() const
+            {
+                return window_;
+            }
+
 
     // Member Operators
 
@@ -177,7 +189,8 @@ public:
                 a.fieldName_ == b.fieldName_
              && a.mean_ == b.mean_
              && a.prime2Mean_ == b.prime2Mean_
-             && a.base_ == b.base_;
+             && a.base_ == b.base_
+             && a.window_ == b.window_;
         }
 
         friend bool operator!=

src/postProcessing/functionObjects/field/fieldAverage/fieldAverageItem/fieldAverageItemIO.C

@@ -33,7 +33,9 @@ Foam::fieldAverageItem::fieldAverageItem(Istream& is)
 :
     fieldName_("unknown"),
     mean_(0),
-    prime2Mean_(0)
+    prime2Mean_(0),
+    base_(ITER),
+    window_(-1.0)
 {
     is.check("Foam::fieldAverageItem::fieldAverageItem(Foam::Istream&)");
 
@@ -43,6 +45,7 @@ Foam::fieldAverageItem::fieldAverageItem(Istream& is)
     entry.lookup("mean") >> mean_;
     entry.lookup("prime2Mean") >> prime2Mean_;
     base_ = baseTypeNames_[entry.lookup("base")];
+    window_ = entry.lookupOrDefault<scalar>("window", -1.0);
 }
 
 
@@ -62,6 +65,7 @@ Foam::Istream& Foam::operator>>(Istream& is, fieldAverageItem& faItem)
     entry.lookup("mean") >> faItem.mean_;
     entry.lookup("prime2Mean") >> faItem.prime2Mean_;
     faItem.base_ = faItem.baseTypeNames_[entry.lookup("base")];
+    faItem.window_ = entry.lookupOrDefault<scalar>("window", -1.0);
 
     return is;
 }
@@ -80,7 +84,15 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const fieldAverageItem& faItem)
     os.writeKeyword("prime2Mean") << faItem.mean_
         << token::END_STATEMENT << nl;
     os.writeKeyword("base") << faItem.baseTypeNames_[faItem.base_]
-        << token::END_STATEMENT << nl << token::END_BLOCK << nl;
+        << token::END_STATEMENT << nl;
+
+    if (faItem.window_ > 0)
+    {
+        os.writeKeyword("window") << faItem.window_
+            << token::END_STATEMENT << nl;
+    }
+
+    os  << token::END_BLOCK << nl;
 
     os.check
     (

src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSource.C

@@ -145,23 +145,8 @@ void Foam::fieldValues::cellSource::initialise(const dictionary& dict)
     if (operation_ == opWeightedAverage)
     {
         dict.lookup("weightField") >> weightFieldName_;
-        if
-        (
-            obr().foundObject<volScalarField>(weightFieldName_)
-        )
-        {
         Info<< "    weight field = " << weightFieldName_;
     }
-        else
-        {
-            FatalErrorIn("cellSource::initialise()")
-                << type() << " " << name_ << ": "
-                << sourceTypeNames_[source_] << "(" << sourceName_ << "):"
-                << nl << "    Weight field " << weightFieldName_
-                << " must be a " << volScalarField::typeName
-                << nl << exit(FatalError);
-        }
-    }
 
     Info<< nl << endl;
 }

src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSource.H

@@ -157,7 +157,8 @@ protected:
         template<class Type>
         tmp<Field<Type> > setFieldValues
         (
-            const word& fieldName
+            const word& fieldName,
+            const bool mustGet = false
         ) const;
 
         //- Apply the 'operation' to the values

src/postProcessing/functionObjects/field/fieldValues/cellSource/cellSourceTemplates.C

@@ -45,7 +45,8 @@ bool Foam::fieldValues::cellSource::validField(const word& fieldName) const
 template<class Type>
 Foam::tmp<Foam::Field<Type> > Foam::fieldValues::cellSource::setFieldValues
 (
-    const word& fieldName
+    const word& fieldName,
+    const bool mustGet
 ) const
 {
     typedef GeometricField<Type, fvPatchField, volMesh> vf;
@@ -55,6 +56,20 @@ Foam::tmp<Foam::Field<Type> > Foam::fieldValues::cellSource::setFieldValues
         return filterField(obr_.lookupObject<vf>(fieldName));
     }
 
+    if (mustGet)
+    {
+        FatalErrorIn
+        (
+            "Foam::tmp<Foam::Field<Type> > "
+            "Foam::fieldValues::cellSource::setFieldValues"
+            "("
+                "const word&, "
+                "const bool"
+            ") const"
+        )   << "Field " << fieldName << " not found in database"
+            << abort(FatalError);
+    }
+
     return tmp<Field<Type> >(new Field<Type>(0.0));
 }
 
@@ -125,7 +140,13 @@ bool Foam::fieldValues::cellSource::writeValues(const word& fieldName)
         scalarField V(filterField(mesh().V()));
         combineFields(V);
 
-        scalarField weightField(setFieldValues<scalar>(weightFieldName_));
+        scalarField weightField;
+
+        if (operation_ == opWeightedAverage)
+        {
+            weightField = setFieldValues<scalar>(weightFieldName_, true);
+        }
+
         combineFields(weightField);
 
         if (Pstream::master())

src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.C

@@ -284,24 +284,8 @@ void Foam::fieldValues::faceSource::initialise(const dictionary& dict)
     if (operation_ == opWeightedAverage)
     {
         dict.lookup("weightField") >> weightFieldName_;
-        if
-        (
-            obr().foundObject<volScalarField>(weightFieldName_)
-         || obr().foundObject<surfaceScalarField>(weightFieldName_)
-        )
-        {
         Info<< "    weight field = " << weightFieldName_;
     }
-        else
-        {
-            FatalErrorIn("faceSource::initialise()")
-                << type() << " " << name_ << ": "
-                << sourceTypeNames_[source_] << "(" << sourceName_ << "):"
-                << nl << "    Weight field " << weightFieldName_
-                << " must be either a " << volScalarField::typeName << " or "
-                << surfaceScalarField::typeName << nl << exit(FatalError);
-        }
-    }
 
     Info<< nl << endl;
 }

src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSource.H

@@ -200,7 +200,11 @@ protected:
 
         //- Return field values by looking up field name
         template<class Type>
-        tmp<Field<Type> > getFieldValues(const word& fieldName) const;
+        tmp<Field<Type> > getFieldValues
+        (
+            const word& fieldName,
+            const bool mustGet = false
+        ) const;
 
         //- Apply the 'operation' to the values
         template<class Type>

src/postProcessing/functionObjects/field/fieldValues/faceSource/faceSourceTemplates.C

@@ -52,7 +52,8 @@ bool Foam::fieldValues::faceSource::validField(const word& fieldName) const
 template<class Type>
 Foam::tmp<Foam::Field<Type> > Foam::fieldValues::faceSource::getFieldValues
 (
-    const word& fieldName
+    const word& fieldName,
+    const bool mustGet
 ) const
 {
     typedef GeometricField<Type, fvsPatchField, surfaceMesh> sf;
@@ -74,6 +75,20 @@ Foam::tmp<Foam::Field<Type> > Foam::fieldValues::faceSource::getFieldValues
         }
     }
 
+    if (mustGet)
+    {
+        FatalErrorIn
+        (
+            "Foam::tmp<Foam::Field<Type> > "
+            "Foam::fieldValues::faceSource::getFieldValues"
+            "("
+                "const word&, "
+                "const bool"
+            ") const"
+        )   << "Field " << fieldName << " not found in database"
+            << abort(FatalError);
+    }
+
     return tmp<Field<Type> >(new Field<Type>(0));
 }
 
@@ -139,7 +154,13 @@ bool Foam::fieldValues::faceSource::writeValues(const word& fieldName)
     if (ok)
     {
         Field<Type> values(getFieldValues<Type>(fieldName));
-        scalarField weightField(getFieldValues<scalar>(weightFieldName_));
+        scalarField weightField;
+
+        if (operation_ == opWeightedAverage)
+        {
+            weightField = getFieldValues<scalar>(weightFieldName_, true);
+        }
+
         scalarField magSf;
 
         if (surfacePtr_.valid())

src/sampling/meshToMeshInterpolation/meshToMesh/calculateMeshToMeshAddressing.C

@@ -106,7 +106,7 @@ void Foam::meshToMesh::calcAddressing()
         shiftedBb,      // overall bounding box
         8,              // maxLevel
         10,             // leafsize
-        3.0             // duplicity
+        6.0             // duplicity
     );
 
     if (debug)

src/sampling/sampledSet/sampledSet/sampledSet.C

@@ -28,6 +28,7 @@ License
 #include "primitiveMesh.H"
 #include "meshSearch.H"
 #include "writer.H"
+#include "particle.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -183,26 +184,58 @@ Foam::point Foam::sampledSet::pushIn
 ) const
 {
     label cellI = mesh().faceOwner()[faceI];
+    const point& cC = mesh().cellCentres()[cellI];
 
-    const point& cellCtr = mesh().cellCentres()[cellI];
+    point newPosition = facePt;
 
-    point newSample =
+    // Taken from particle::initCellFacePt()
-        facePt + tol*(cellCtr - facePt);
+    label tetFaceI;
+    label tetPtI;
+    mesh().findTetFacePt(cellI, facePt, tetFaceI, tetPtI);
 
-    if (!searchEngine().pointInCell(newSample, cellI))
+    if (tetFaceI == -1 || tetPtI == -1)
+    {
+        newPosition = facePt;
+
+        label trap(1.0/particle::trackingCorrectionTol + 1);
+
+        label iterNo = 0;
+
+        do
+        {
+            newPosition += particle::trackingCorrectionTol*(cC - facePt);
+
+            mesh().findTetFacePt
+            (
+                cellI,
+                newPosition,
+                tetFaceI,
+                tetPtI
+            );
+
+            iterNo++;
+
+        } while (tetFaceI < 0  && iterNo <= trap);
+    }
+
+    if (tetFaceI == -1)
     {
         FatalErrorIn
         (
             "sampledSet::pushIn(const point&, const label)"
-        )   << "After pushing " << facePt << " to " << newSample
+        )   << "After pushing " << facePt << " to " << newPosition
-            << " it is still outside faceI " << faceI << endl
+            << " it is still outside face " << faceI
+            << " at " << mesh().faceCentres()[faceI]
+            << " of cell " << cellI
+            << " at " << cC << endl
             << "Please change your starting point"
             << abort(FatalError);
     }
-    //Info<< "pushIn : moved " << facePt << " to " << newSample
+
+    //Info<< "pushIn : moved " << facePt << " to " << newPosition
     //    << endl;
 
-    return newSample;
+    return newPosition;
 }