Merge branch 'master' into molecularDynamics

src/OpenFOAM/containers/HashTables/HashSet/HashSet.C

@@ -0,0 +1,90 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef HashSet_C
+#define HashSet_C
+
+#include "HashSet.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+template<class Key, class Hash>
+bool HashSet<Key, Hash>::operator==(const HashSet<Key, Hash>& ht) const
+{
+    const HashTable<empty, Key, Hash>& a = *this;
+
+    // Are all my elements in ht?
+    for
+    (
+        typename HashTable<empty, Key, Hash>::const_iterator iter = a.begin();
+        iter != a.end();
+        ++iter
+    )
+    {
+        if (!ht.found(iter.key()))
+        {
+            return false;
+        }
+    }
+
+    // Are all ht elements in me?
+    for
+    (
+        typename HashTable<empty, Key, Hash>::const_iterator iter = ht.begin();
+        iter != ht.end();
+        ++iter
+    )
+    {
+        if (!found(iter.key()))
+        {
+            return false;
+        }
+    }
+    return true;
+}
+
+
+template<class Key, class Hash>
+bool HashSet<Key, Hash>::operator!=(const HashSet<Key, Hash>& ht) const
+{
+    return !(operator==(ht));
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/containers/HashTables/HashSet/HashSet.H

@@ -100,6 +100,17 @@ public:
             {
                 return HashTable<empty, Key, Hash>::insert(key, empty());
             }
+
+
+    // Member Operators
+
+        //- Equality. Two hashtables are equal if all contents of first are
+        //  also in second and vice versa. So does not depend on table size or
+        //  order!
+        bool operator==(const HashSet<Key, Hash>&) const;
+
+        //- The opposite of the equality operation.
+        bool operator!=(const HashSet<Key, Hash>&) const;
 };
 
 
@@ -112,6 +123,12 @@ typedef HashSet<> wordHashSet;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#ifdef NoRepository
+#   include "HashSet.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //

src/OpenFOAM/containers/Lists/List/List.C

@@ -321,7 +321,7 @@ List<T>::List(const BiIndirectList<T>& idl)
 template<class T>
 List<T>::~List()
 {
-    if (this->size_) delete[] this->v_;
+    if (this->v_) delete[] this->v_;
 }
 
 
@@ -367,9 +367,8 @@ void List<T>::setSize(const label newSize)
                     register T* av = &nv[i];
                     while (i--) *--av = *--vv;
                 }
-
-                delete[] this->v_;
             }
+            if (this->v_) delete[] this->v_;
 
             this->size_ = newSize;
             this->v_ = nv;
@@ -400,7 +399,7 @@ void List<T>::setSize(const label newSize, const T& a)
 template<class T>
 void List<T>::clear()
 {
-    if (this->size_) delete[] this->v_;
+    if (this->v_) delete[] this->v_;
     this->size_ = 0;
     this->v_ = 0;
 }
@@ -411,7 +410,7 @@ void List<T>::clear()
 template<class T>
 void List<T>::transfer(List<T>& a)
 {
-    if (this->size_) delete[] this->v_;
+    if (this->v_) delete[] this->v_;
 
     this->size_ = a.size_;
     this->v_ = a.v_;
@@ -457,7 +456,8 @@ void List<T>::operator=(const UList<T>& a)
 {
     if (a.size_ != this->size_)
     {
-        if (this->size_) delete[] this->v_;
+        if (this->v_) delete[] this->v_;
+        this->v_ = 0;
         this->size_ = a.size_;
         if (this->size_) this->v_ = new T[this->size_];
     }
@@ -503,7 +503,8 @@ void List<T>::operator=(const SLList<T>& sll)
 {
     if (sll.size() != this->size_)
     {
-        if (this->size_) delete[] this->v_;
+        if (this->v_) delete[] this->v_;
+        this->v_ = 0;
         this->size_ = sll.size();
         if (this->size_) this->v_ = new T[this->size_];
     }
@@ -530,7 +531,8 @@ void List<T>::operator=(const IndirectList<T>& idl)
 {
     if (idl.size() != this->size_)
     {
-        if (this->size_) delete[] this->v_;
+        if (this->v_) delete[] this->v_;
+        this->v_ = 0;
         this->size_ = idl.size();
         if (this->size_) this->v_ = new T[this->size_];
     }
@@ -551,7 +553,8 @@ void List<T>::operator=(const BiIndirectList<T>& idl)
 {
     if (idl.size() != this->size_)
     {
-        if (this->size_) delete[] this->v_;
+        if (this->v_) delete[] this->v_;
+        this->v_ = 0;
         this->size_ = idl.size();
         if (this->size_) this->v_ = new T[this->size_];
     }

src/OpenFOAM/containers/Lists/List/List.H

@@ -137,6 +137,9 @@ public:
         //- Return a null List
         static const List<T>& null();
 
+        //- Return the number of elements in the UList.
+        inline label size() const;
+
 
         // Edit
 
@@ -156,6 +159,10 @@ public:
             //- Return subscript-checked element of UList.
             inline T& newElmt(const label);
 
+            //- Override size to be inconsistent with allocated storage.
+            //  Use with care.
+            inline label& size();
+
 
     // Member operators
 

src/OpenFOAM/containers/Lists/List/ListI.H

@@ -52,6 +52,20 @@ inline T& Foam::List<T>::newElmt(const label i)
 }
 
 
+template<class T>
+inline Foam::label Foam::List<T>::size() const
+{
+    return UList<T>::size_;
+}
+
+
+template<class T>
+inline Foam::label& Foam::List<T>::size()
+{
+    return UList<T>::size_;
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template<class T>

src/OpenFOAM/db/Time/TimeIO.C

@@ -86,13 +86,6 @@ void Foam::Time::readDict()
 
             purgeWrite_ = 0;
         }
-
-        if (writeControl_ != wcTimeStep && purgeWrite_ > 0)
-        {
-            FatalIOErrorIn("Time::readDict()", controlDict_)
-                << "writeControl must be set to timeStep for purgeWrite "
-                << exit(FatalIOError);
-        }
     }
 
     if (controlDict_.found("timeFormat"))

src/OpenFOAM/meshes/boundBox/boundBox.C

@@ -87,9 +87,47 @@ boundBox::boundBox(Istream& is)
 
 // * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
 
-Ostream& operator<<(Ostream& os, const boundBox& b)
+Ostream& operator<<(Ostream& os, const boundBox& bb)
 {
-    return os << b.min() << token::SPACE << b.max();
+    if (os.format() == IOstream::ASCII)
+    {
+        os << bb.min_ << token::SPACE << bb.max_;
+    }
+    else
+    {
+        os.write
+        (
+            reinterpret_cast<const char*>(&bb.min_),
+            sizeof(boundBox)
+        );
+    }
+
+    // Check state of Ostream
+    os.check("Ostream& operator<<(Ostream&, const boundBox&)");
+
+    return os;
+}
+
+
+Istream& operator>>(Istream& is, boundBox& bb)
+{
+    if (is.format() == IOstream::ASCII)
+    {
+        return is >> bb.min_ >>  bb.max_;
+    }
+    else
+    {
+        is.read
+        (
+            reinterpret_cast<char*>(&bb.min_),
+            sizeof(boundBox)
+        );
+    }
+
+    // Check state of Istream
+    is.check("Istream& operator>>(Istream&, boundBox&)");
+
+    return is;
 }
 
 

src/OpenFOAM/meshes/boundBox/boundBox.H

@@ -153,12 +153,31 @@ public:
             }
 
 
-    // Ostream operator
+    // Friend Operators
 
-        friend Ostream& operator<<(Ostream& os, const boundBox& b);
+        friend bool operator==(const boundBox& a, const boundBox& b)
+        {
+            return (a.min_ == b.min_) && (a.max_ == b.max_);
+        }
+
+        friend bool operator!=(const boundBox& a, const boundBox& b)
+        {
+            return !(a == b);
+        }
+
+
+    // IOstream operator
+
+        friend Istream& operator>>(Istream& is, boundBox&);
+        friend Ostream& operator<<(Ostream& os, const boundBox&);
 };
 
 
+//- Specify data associated with boundBox type is contiguous
+template<>
+inline bool contiguous<boundBox>() {return contiguous<point>();}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.C

@@ -282,8 +282,7 @@ const Foam::labelList& Foam::polyPatch::meshEdges() const
                 primitivePatch::meshEdges
                 (
                     boundaryMesh().mesh().edges(),
-                    boundaryMesh().mesh().cellEdges(),
-                    faceCells()
+                    boundaryMesh().mesh().pointEdges()
                 )
             );
     }

src/OpenFOAM/meshes/polyMesh/syncTools/syncTools.C

@@ -169,7 +169,7 @@ Foam::PackedList<1> Foam::syncTools::getMasterPoints(const polyMesh& mesh)
             else
             {
                 FatalErrorIn("syncTools::getMasterPoints(const polyMesh&)")
-                    << "Cannot handle patch " << patches[patchI].name()
+                    << "Cannot handle coupled patch " << patches[patchI].name()
                     << " of type " <<  patches[patchI].type()
                     << abort(FatalError);
             }
@@ -290,7 +290,7 @@ Foam::PackedList<1> Foam::syncTools::getMasterEdges(const polyMesh& mesh)
             else
             {
                 FatalErrorIn("syncTools::getMasterEdges(const polyMesh&)")
-                    << "Cannot handle patch " << patches[patchI].name()
+                    << "Cannot handle coupled patch " << patches[patchI].name()
                     << " of type " <<  patches[patchI].type()
                     << abort(FatalError);
             }
@@ -314,6 +314,54 @@ Foam::PackedList<1> Foam::syncTools::getMasterEdges(const polyMesh& mesh)
 }
 
 
+// Determines for every face whether it is coupled and if so sets only one.
+Foam::PackedList<1> Foam::syncTools::getMasterFaces(const polyMesh& mesh)
+{
+    PackedList<1> isMasterFace(mesh.nFaces(), 1);
+
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    forAll(patches, patchI)
+    {
+        if (patches[patchI].coupled())
+        {
+            if (Pstream::parRun() && isA<processorPolyPatch>(patches[patchI]))
+            {
+                const processorPolyPatch& pp =
+                    refCast<const processorPolyPatch>(patches[patchI]);
+
+                if (!pp.owner())
+                {
+                    forAll(pp, i)
+                    {
+                        isMasterFace.set(pp.start()+i, 0);
+                    }
+                }
+            }
+            else if (isA<cyclicPolyPatch>(patches[patchI]))
+            {
+                const cyclicPolyPatch& pp =
+                    refCast<const cyclicPolyPatch>(patches[patchI]);
+
+                for (label i = pp.size()/2; i < pp.size(); i++)
+                {
+                    isMasterFace.set(pp.start()+i, 0);
+                }
+            }
+            else
+            {
+                FatalErrorIn("syncTools::getMasterFaces(const polyMesh&)")
+                    << "Cannot handle coupled patch " << patches[patchI].name()
+                    << " of type " <<  patches[patchI].type()
+                    << abort(FatalError);
+            }
+        }
+    }
+
+    return isMasterFace;
+}
+
+
 template <>
 void Foam::syncTools::separateList
 (

src/OpenFOAM/meshes/polyMesh/syncTools/syncTools.H

@@ -226,9 +226,12 @@ public:
             //- Get per point whether is it master (of a coupled set of points)
             static PackedList<1> getMasterPoints(const polyMesh&);
 
-            //- Get per edge whether is it master (of a coupled set of edge)
+            //- Get per edge whether is it master (of a coupled set of edges)
             static PackedList<1> getMasterEdges(const polyMesh&);
 
+            //- Get per face whether is it master (of a coupled set of faces)
+            static PackedList<1> getMasterFaces(const polyMesh&);
+
 };
 
 

src/OpenFOAM/meshes/polyMesh/zones/faceZone/faceZone.C

@@ -362,16 +362,27 @@ const Foam::labelList& Foam::faceZone::meshEdges() const
 {
     if (!mePtr_)
     {
-        labelList faceCells(size());
-
-        const labelList& own = zoneMesh().mesh().faceOwner();
-
-        const labelList& faceLabels = *this;
-
-        forAll (faceCells, faceI)
-        {
-            faceCells[faceI] = own[faceLabels[faceI]];
-        }
+        //labelList faceCells(size());
+        //
+        //const labelList& own = zoneMesh().mesh().faceOwner();
+        //
+        //const labelList& faceLabels = *this;
+        //
+        //forAll (faceCells, faceI)
+        //{
+        //    faceCells[faceI] = own[faceLabels[faceI]];
+        //}
+        //
+        //mePtr_ =
+        //    new labelList
+        //    (
+        //        operator()().meshEdges
+        //        (
+        //            zoneMesh().mesh().edges(),
+        //            zoneMesh().mesh().cellEdges(),
+        //            faceCells
+        //        )
+        //    );
 
         mePtr_ =
             new labelList
@@ -379,8 +390,7 @@ const Foam::labelList& Foam::faceZone::meshEdges() const
                 operator()().meshEdges
                 (
                     zoneMesh().mesh().edges(),
-                    zoneMesh().mesh().cellEdges(),
-                    faceCells
+                    zoneMesh().mesh().pointEdges()
                 )
             );
     }

src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.H

@@ -348,7 +348,8 @@ public:
                 //  index in the edge list.  If the edge is not found, return -1
                 label whichEdge(const edge& e) const;
 
-                //- Return labels of patch edges in the global edge list
+                //- Return labels of patch edges in the global edge list using
+                //  cell addressing
                 labelList meshEdges
                 (
                     const edgeList& allEdges,
@@ -356,6 +357,14 @@ public:
                     const labelList& faceCells
                 ) const;
 
+                //- Return labels of patch edges in the global edge list using
+                //  basic edge addressing.
+                labelList meshEdges
+                (
+                    const edgeList& allEdges,
+                    const labelListList& pointEdges
+                ) const;
+
                 //- Return face normals for patch
                 const Field<PointType>& faceNormals() const;
 

src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchMeshEdges.C

@@ -111,6 +111,60 @@ labelList PrimitivePatch<Face, FaceList, PointField, PointType>::meshEdges
 }
 
 
+template
+<
+    class Face,
+    template<class> class FaceList,
+    class PointField,
+    class PointType
+>
+
+labelList PrimitivePatch<Face, FaceList, PointField, PointType>::meshEdges
+(
+    const edgeList& allEdges,
+    const labelListList& pointEdges
+) const
+{
+    if (debug)
+    {
+        Info<< "labelList PrimitivePatch<Face, FaceList, PointField, PointType>"
+            << "::meshEdges() : "
+            << "calculating labels of patch edges in mesh edge list"
+            << endl;
+    }
+
+    // get reference to the list of edges on the patch 
+    const edgeList& PatchEdges = edges();
+
+    // create the storage
+    labelList meshEdges(PatchEdges.size());
+
+    // get reference to the points on the patch
+    const labelList& pp = meshPoints();
+
+    // WARNING: Remember that local edges address into local point list;
+    // local-to-global point label translation is necessary
+    forAll (PatchEdges, edgeI)
+    {
+        const label globalPointI = pp[PatchEdges[edgeI].start()];
+        const edge curEdge(globalPointI, pp[PatchEdges[edgeI].end()]);
+
+        const labelList& pe = pointEdges[globalPointI];
+
+        forAll (pe, i)
+        {
+            if (allEdges[pe[i]] == curEdge)
+            {
+                meshEdges[edgeI] = pe[i];
+                break;
+            } 
+        }
+    }
+
+    return meshEdges;
+}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template

src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.C

@@ -66,6 +66,9 @@ primitiveMesh::primitiveMesh()
     ppPtr_(NULL),
     cpPtr_(NULL),
 
+    allocSize_(0),
+    labels_(0),
+
     cellCentresPtr_(NULL),
     faceCentresPtr_(NULL),
     cellVolumesPtr_(NULL),
@@ -106,6 +109,9 @@ primitiveMesh::primitiveMesh
     ppPtr_(NULL),
     cpPtr_(NULL),
 
+    allocSize_(0),
+    labels_(0),
+
     cellCentresPtr_(NULL),
     faceCentresPtr_(NULL),
     cellVolumesPtr_(NULL),

src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H

@@ -155,6 +155,17 @@ class primitiveMesh
             mutable labelListList* cpPtr_;
 
 
+        // On-the-fly edge addresing storage
+
+            //- Temporary storage for addressing. allocSize is the real size
+            //  of the labelList.
+            mutable label allocSize_;
+            mutable labelList labels_;
+
+            //- Temporary storage for addressing
+            mutable labelHashSet labelSet_;
+
+
         // Geometric data
 
             //- Cell centres
@@ -209,8 +220,14 @@ class primitiveMesh
             (
                 List<DynamicList<label> >&,
                 DynamicList<edge>&,
-                const label pA,
-                const label pB
+                const label,
+                const label
+            );
+            //- For on-the-fly addressing calculation
+            static label findFirstCommonElementFromSortedLists
+            ( 
+                const labelList&,
+                const labelList&
             );
 
 
@@ -667,6 +684,55 @@ public:
             //- Print a list of all the currently allocated mesh data
             void printAllocated() const;
 
+            // Per storage whether allocated
+            inline bool hasCellShapes() const;
+            inline bool hasEdges() const;
+            inline bool hasCellCells() const;
+            inline bool hasEdgeCells() const;
+            inline bool hasPointCells() const;
+            inline bool hasCells() const;
+            inline bool hasEdgeFaces() const;
+            inline bool hasPointFaces() const;
+            inline bool hasCellEdges() const;
+            inline bool hasFaceEdges() const;
+            inline bool hasPointEdges() const;
+            inline bool hasPointPoints() const;
+            inline bool hasCellPoints() const;
+            inline bool hasCellCentres() const;
+            inline bool hasFaceCentres() const;
+            inline bool hasCellVolumes() const;
+            inline bool hasFaceAreas() const;
+
+            // On-the-fly addressing calculation. These functions return either
+            // a reference to the full addressing (if already calculated) or
+            // a reference to member data labels_ so be careful when not storing
+            // result.
+
+            //- cellCells using cells
+            const labelList& cellCells(const label cellI) const;
+
+            //- cellPoints using cells
+            const labelList& cellPoints(const label cellI) const;
+
+            //- pointCells using pointFaces
+            const labelList& pointCells(const label pointI) const;
+
+            //- pointPoints using edges, pointEdges
+            const labelList& pointPoints(const label pointI) const;
+
+            //- faceEdges using pointFaces, edges, pointEdges
+            const labelList& faceEdges(const label faceI) const;
+
+            //- edgeFaces using pointFaces, edges, pointEdges
+            const labelList& edgeFaces(const label edgeI) const;
+
+            //- edgeCells using pointFaces, edges, pointEdges
+            const labelList& edgeCells(const label edgeI) const;
+
+            //- cellEdges using cells, pointFaces, edges, pointEdges
+            const labelList& cellEdges(const label cellI) const;
+
+
             //- Clear geometry
             void clearGeom();
 

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellCells.C

@@ -105,6 +105,53 @@ const labelListList& primitiveMesh::cellCells() const
 }
 
 
+const labelList& primitiveMesh::cellCells(const label cellI) const
+{
+    if (hasCellCells())
+    {
+        return cellCells()[cellI];
+    }
+    else
+    {
+        const labelList& own = faceOwner();
+        const labelList& nei = faceNeighbour();
+        const cell& cFaces = cells()[cellI];
+
+        labels_.size() = allocSize_;
+
+        if (cFaces.size() > allocSize_)
+        {
+            labels_.clear();
+            allocSize_ = cFaces.size();
+            labels_.setSize(allocSize_);
+        }
+
+        label n = 0;
+
+        forAll(cFaces, i)
+        {
+            label faceI = cFaces[i];
+
+            if (faceI < nInternalFaces())
+            {
+                if (own[faceI] == cellI)
+                {
+                    labels_[n++] = nei[faceI];
+                }
+                else
+                {
+                    labels_[n++] = own[faceI];
+                }
+            }
+        }
+
+        labels_.size() = n;
+
+        return labels_;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellPoints.C

@@ -53,6 +53,52 @@ const labelListList& primitiveMesh::cellPoints() const
 }
 
 
+const labelList& primitiveMesh::cellPoints(const label cellI) const
+{
+    if (hasCellPoints())
+    {
+        return cellPoints()[cellI];
+    }
+    else
+    {
+        const faceList& fcs = faces();
+        const labelList& cFaces = cells()[cellI];
+
+        labelSet_.clear();
+
+        forAll(cFaces, i)
+        {
+            const labelList& f = fcs[cFaces[i]];
+
+            forAll(f, fp)
+            {
+                labelSet_.insert(f[fp]);    
+            }
+        }
+
+        labels_.size() = allocSize_;
+
+        if (labelSet_.size() > allocSize_)
+        {
+            labels_.clear();
+            allocSize_ = labelSet_.size();
+            labels_.setSize(allocSize_);
+        }
+
+        label n = 0;
+
+        forAllConstIter(labelHashSet, labelSet_, iter)
+        {
+            labels_[n++] = iter.key();
+        }
+
+        labels_.size() = n;
+
+        return labels_;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C

@@ -771,11 +771,12 @@ bool primitiveMesh::checkPoints
         }
     }
 
-    const labelListList& pc = pointCells();
 
-    forAll (pc, pointI)
+    forAll (pf, pointI)
     {
-        if (pc[pointI].size() == 0)
+        const labelList& pc = pointCells(pointI);
+
+        if (pc.size() == 0)
         {
             if (setPtr)
             {

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdgeCells.C

@@ -51,6 +51,86 @@ const labelListList& primitiveMesh::edgeCells() const
 }
 
 
+const labelList& primitiveMesh::edgeCells(const label edgeI) const
+{
+    if (hasEdgeCells())
+    {
+        return edgeCells()[edgeI];
+    }
+    else
+    {
+        const labelList& own = faceOwner();
+        const labelList& nei = faceNeighbour();
+
+        // edge faces can either return labels_ or reference in edgeLabels.
+        labelList labelsCopy;
+        if (!hasEdgeFaces())
+        {
+            labelsCopy = edgeFaces(edgeI);
+        }
+
+        const labelList& eFaces =
+        (
+            hasEdgeFaces()
+          ? edgeFaces()[edgeI]
+          : labelsCopy
+        );
+
+        labels_.size() = allocSize_;
+
+        // labels_ should certainly be big enough for edge cells.
+        label n = 0;
+
+        // Do quadratic insertion.
+        forAll(eFaces, i)
+        {
+            label faceI = eFaces[i];
+
+            {
+                label ownCellI = own[faceI];
+
+                // Check if not already in labels_
+                for (label j = 0; j < n; j++)
+                {
+                    if (labels_[j] == ownCellI)
+                    {
+                        ownCellI = -1;
+                        break;
+                    }
+                }
+
+                if (ownCellI != -1)
+                {
+                    labels_[n++] = ownCellI;
+                }
+            }
+
+            if (isInternalFace(faceI))
+            {
+                label neiCellI = nei[faceI];
+
+                for (label j = 0; j < n; j++)
+                {
+                    if (labels_[j] == neiCellI)
+                    {
+                        neiCellI = -1;
+                        break;
+                    }
+                }
+
+                if (neiCellI != -1)
+                {
+                    labels_[n++] = neiCellI;
+                }
+            }
+        }
+
+        labels_.size() = n;
+
+        return labels_;
+    }
+}
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdgeFaces.C

@@ -51,6 +51,59 @@ const labelListList& primitiveMesh::edgeFaces() const
     return *efPtr_;
 }
 
+
+const labelList& primitiveMesh::edgeFaces(const label edgeI) const
+{
+    if (hasEdgeFaces())
+    {
+        return edgeFaces()[edgeI];
+    }
+    else
+    {
+        // Use the fact that pointEdges are sorted in incrementing edge order
+        const edge& e = edges()[edgeI];
+        const labelList& pFaces0 = pointFaces()[e[0]];
+        const labelList& pFaces1 = pointFaces()[e[1]];
+
+        label i0 = 0;
+        label i1 = 0;
+        label n = 0;
+
+        labels_.size() = allocSize_;
+
+        while (i0 < pFaces0.size() && i1 < pFaces1.size())
+        {
+            if (pFaces0[i0] < pFaces1[i1])
+            {
+                ++i0;
+            }
+            else if (pFaces0[i0] > pFaces1[i1])
+            {
+                ++i1;
+            }
+            else
+            {
+                // Equal. Append.
+                if (n == allocSize_)
+                {
+                    // Have setSize copy contents so far
+                    labels_.size() = n;
+                    allocSize_ = allocSize_*2 + 1;
+                    labels_.setSize(allocSize_);
+                }
+                labels_[n++] = pFaces0[i0];
+                ++i0;
+                ++i1;
+            }
+        }
+
+        labels_.size() = n;
+
+        return labels_;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdges.C

@@ -461,6 +461,46 @@ void primitiveMesh::calcEdges(const bool doFaceEdges) const
 }
 
 
+label primitiveMesh::findFirstCommonElementFromSortedLists
+( 
+    const labelList& list1,
+    const labelList& list2
+)
+{
+    label result = -1;
+
+    labelList::const_iterator iter1 = list1.begin();
+    labelList::const_iterator iter2 = list2.begin();
+
+    while (iter1 != list1.end() && iter2 != list2.end())
+    {
+        if( *iter1 < *iter2)
+        {
+            ++iter1;
+        }
+        else if (*iter1 > *iter2)
+        {
+            ++iter2;
+        }
+        else
+        {
+            result = *iter1;
+            break;
+        }
+    }
+    if (result == -1)
+    {
+        FatalErrorIn
+        (
+            "primitiveMesh::findFirstCommonElementFromSortedLists"
+            "(const labelList&, const labelList&)"
+        )   << "No common elements in lists " << list1 << " and " << list2
+            << abort(FatalError);        
+    }
+    return result;
+}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 const edgeList& primitiveMesh::edges() const
@@ -542,6 +582,91 @@ void primitiveMesh::clearOutEdges()
     deleteDemandDrivenData(edgesPtr_);
     deleteDemandDrivenData(pePtr_);
     deleteDemandDrivenData(fePtr_);
+    labels_.clear();
+    allocSize_ = 0;
+}
+
+
+const labelList& primitiveMesh::faceEdges(const label faceI) const
+{
+    if (hasFaceEdges())
+    {
+        return faceEdges()[faceI];
+    }
+    else
+    {
+        const labelListList& pointEs = pointEdges();
+        const face& f = faces()[faceI];
+
+        labels_.size() = allocSize_;
+
+        if (f.size() > allocSize_)
+        {
+            labels_.clear();
+            allocSize_ = f.size();
+            labels_.setSize(allocSize_);
+        }
+
+        label n = 0;
+
+        forAll(f, fp)
+        {
+            labels_[n++] = findFirstCommonElementFromSortedLists
+            (
+                pointEs[f[fp]],
+                pointEs[f.nextLabel(fp)]
+            );
+        }
+
+        labels_.size() = n;
+
+        return labels_;
+    }
+}
+
+
+const labelList& primitiveMesh::cellEdges(const label cellI) const
+{
+    if (hasCellEdges())
+    {
+        return cellEdges()[cellI];
+    }
+    else
+    {
+        const labelList& cFaces = cells()[cellI];
+
+        labelSet_.clear();
+
+        forAll(cFaces, i)
+        {
+            const labelList& fe = faceEdges(cFaces[i]);
+
+            forAll(fe, feI)
+            {
+                labelSet_.insert(fe[feI]);    
+            }
+        }
+
+        labels_.size() = allocSize_;
+
+        if (labelSet_.size() > allocSize_)
+        {
+            labels_.clear();
+            allocSize_ = labelSet_.size();
+            labels_.setSize(allocSize_);
+        }
+
+        label n =0;
+
+        forAllConstIter(labelHashSet, labelSet_, iter)
+        {
+            labels_[n++] = iter.key();
+        }
+
+        labels_.size() = n;
+
+        return labels_;
+    }
 }
 
 

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshI.H

@@ -104,6 +104,108 @@ inline bool primitiveMesh::isInternalFace(const label faceIndex) const
 }
 
 
+inline bool primitiveMesh::hasCellShapes() const
+{
+    return cellShapesPtr_;
+}
+
+
+inline bool primitiveMesh::hasEdges() const
+{
+    return edgesPtr_;
+}
+
+
+inline bool primitiveMesh::hasCellCells() const
+{
+    return ccPtr_;
+}
+
+
+inline bool primitiveMesh::hasEdgeCells() const
+{
+    return ecPtr_;
+}
+
+
+inline bool primitiveMesh::hasPointCells() const
+{
+    return pcPtr_;
+}
+
+
+inline bool primitiveMesh::hasCells() const
+{
+    return cfPtr_;
+}
+
+
+inline bool primitiveMesh::hasEdgeFaces() const
+{
+    return efPtr_;
+}
+
+
+inline bool primitiveMesh::hasPointFaces() const
+{
+    return pfPtr_;
+}
+
+
+inline bool primitiveMesh::hasCellEdges() const
+{
+    return cePtr_;
+}
+
+
+inline bool primitiveMesh::hasFaceEdges() const
+{
+    return fePtr_;
+}
+
+
+inline bool primitiveMesh::hasPointEdges() const
+{
+    return pePtr_;
+}
+
+
+inline bool primitiveMesh::hasPointPoints() const
+{
+    return ppPtr_;
+}
+
+
+inline bool primitiveMesh::hasCellPoints() const
+{
+    return cpPtr_;
+}
+
+
+inline bool primitiveMesh::hasCellCentres() const
+{
+    return cellCentresPtr_;
+}
+
+
+inline bool primitiveMesh::hasFaceCentres() const
+{
+    return faceCentresPtr_;
+}
+
+
+inline bool primitiveMesh::hasCellVolumes() const
+{
+    return cellVolumesPtr_;
+}
+
+
+inline bool primitiveMesh::hasFaceAreas() const
+{
+    return faceAreasPtr_;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshPointCells.C

@@ -114,6 +114,70 @@ const labelListList& primitiveMesh::pointCells() const
 }
 
 
+const labelList& primitiveMesh::pointCells(const label pointI) const
+{
+    if (hasPointCells())
+    {
+        return pointCells()[pointI];
+    }
+    else
+    {
+        const labelList& own = faceOwner();
+        const labelList& nei = faceNeighbour();
+        const labelList& pFaces = pointFaces()[pointI];
+
+        labels_.size() = allocSize_;
+
+        label n = 0;
+
+        forAll(pFaces, i)
+        {
+            const label faceI = pFaces[i];
+
+            // Append owner
+            if (n == allocSize_)
+            {
+                labels_.size() = n;
+                allocSize_ = allocSize_*2 + 1;
+                labels_.setSize(allocSize_);
+            }
+            labels_[n++] = own[faceI];
+
+            // Append neighbour
+            if (faceI < nInternalFaces())
+            {
+                if (n == allocSize_)
+                {
+                    labels_.size() = n;
+                    allocSize_ = allocSize_*2 + 1;
+                    labels_.setSize(allocSize_);
+                }
+                labels_[n++] = nei[faceI];
+            }
+        }
+        labels_.size() = n;
+
+
+        // Filter duplicates
+        sort(labels_);
+
+        n = 1;
+
+        for (label i = 1; i < labels_.size(); i++)
+        {
+            if (labels_[i] != labels_[i-1])
+            {
+                labels_[n++] = labels_[i];
+            }
+        }
+
+        labels_.size() = n;
+
+        return labels_;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshPointPoints.C

@@ -97,6 +97,40 @@ const labelListList& primitiveMesh::pointPoints() const
 }
 
 
+const labelList& primitiveMesh::pointPoints(const label pointI) const
+{
+    if (hasPointPoints())
+    {
+        return pointPoints()[pointI];
+    }
+    else
+    {
+        const edgeList& edges = this->edges();
+        const labelList& pEdges = pointEdges()[pointI];
+
+        labels_.size() = allocSize_;
+
+        if (pEdges.size() > allocSize_)
+        {
+            // Set size() so memory allocation behaves as normal.
+            labels_.clear();
+            allocSize_ = pEdges.size();
+            labels_.setSize(allocSize_);
+        }
+
+        label n = 0;
+
+        forAll(pEdges, i)
+        {
+            labels_[n++] = edges[pEdges[i]].otherVertex(pointI);
+        }
+
+        labels_.size() = n;
+
+        return labels_;
+    }
+}
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/OpenFOAM/primitives/contiguous/contiguous.H

@@ -47,34 +47,52 @@ namespace Foam
 
 // Assume the data associated with type T is not contiguous
 template<class T>
-inline bool contiguous()              {return false;}
+inline bool contiguous()                        {return false;}
 
 
 // Specify data associated with primitive types is contiguous
 
 template<>
-inline bool contiguous<bool>()        {return true;}
+inline bool contiguous<bool>()                  {return true;}
 
 template<>
-inline bool contiguous<char>()        {return true;}
+inline bool contiguous<char>()                  {return true;}
 
 template<>
-inline bool contiguous<short>()       {return true;}
+inline bool contiguous<unsigned char>()         {return true;}
 
 template<>
-inline bool contiguous<int>()         {return true;}
+inline bool contiguous<short>()                 {return true;}
 
 template<>
-inline bool contiguous<long>()        {return true;}
+inline bool contiguous<unsigned short>()        {return true;}
 
 template<>
-inline bool contiguous<float>()       {return true;}
+inline bool contiguous<int>()                   {return true;}
 
 template<>
-inline bool contiguous<double>()      {return true;}
+inline bool contiguous<unsigned int>()          {return true;}
 
 template<>
-inline bool contiguous<long double>() {return true;}
+inline bool contiguous<long>()                  {return true;}
+
+template<>
+inline bool contiguous<unsigned long>()         {return true;}
+
+template<>
+inline bool contiguous<long long>()             {return true;}
+
+template<>
+inline bool contiguous<unsigned long long>()    {return true;}
+
+template<>
+inline bool contiguous<float>()                 {return true;}
+
+template<>
+inline bool contiguous<double>()                {return true;}
+
+template<>
+inline bool contiguous<long double>()           {return true;}
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/autoMesh/autoHexMesh/shellSurfaces/shellSurfaces.C

@@ -30,6 +30,7 @@ License
 #include "triSurfaceMesh.H"
 #include "refinementSurfaces.H"
 #include "searchableSurfaces.H"
+#include "orientedSurface.H"
 #include "pointIndexHit.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -211,7 +212,27 @@ void Foam::shellSurfaces::orient()
                     refCast<const triSurfaceMesh>(s)
                 );
 
-                refinementSurfaces::orientSurface(outsidePt, shell);
+                // Flip surface so outsidePt is outside.
+                bool anyFlipped = orientedSurface::orient
+                (
+                    shell,
+                    outsidePt,
+                    true
+                );
+
+                if (anyFlipped)
+                {
+                    // orientedSurface will have done a clearOut of the surface.
+                    // we could do a clearout of the triSurfaceMeshes::trees()
+                    // but these aren't affected by orientation
+                    // (except for cached
+                    // sideness which should not be set at this point.
+                    // !!Should check!)
+
+                    Info<< "shellSurfaces : Flipped orientation of surface "
+                        << s.name()
+                        << " so point " << outsidePt << " is outside." << endl;
+                }
             }
         }
     }

src/conversion/polyDualMesh/polyDualMesh.C

@@ -746,8 +746,7 @@ void Foam::polyDualMesh::calcDual
         allBoundary.meshEdges
         (
             mesh.edges(),
-            mesh.cellEdges(),
-            SubList<label>(mesh.faceOwner(), allBoundary.size(), nIntFaces)
+            mesh.pointEdges()
         )
     );
 

src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C

@@ -29,8 +29,8 @@ License
 #include "matchPoints.H"
 #include "indirectPrimitivePatch.H"
 #include "meshTools.H"
-#include "octreeDataFace.H"
-#include "octree.H"
+#include "treeDataFace.H"
+#include "indexedOctree.H"
 #include "OFstream.H"
 #include "IndirectList.H"
 
@@ -1011,19 +1011,29 @@ void Foam::faceCoupleInfo::findSlavesCoveringMaster
 )
 {
     // Construct octree from all mesh0 boundary faces
-    octreeDataFace shapes(mesh0);
+    labelList bndFaces(mesh0.nFaces()-mesh0.nInternalFaces());
+    forAll(bndFaces, i)
+    {
+        bndFaces[i] = mesh0.nInternalFaces() + i;
+    }
 
     treeBoundBox overallBb(mesh0.points());
 
-    octree<octreeDataFace> tree
-    (
-        overallBb,  // overall search domain
-        shapes,     // all information needed to do checks on cells
-        1,          // min levels
-        20.0,       // maximum ratio of cubes v.s. cells
-        10.0
-    );
+    Random rndGen(123456);
 
+    indexedOctree<treeDataFace> tree
+    (
+        treeDataFace    // all information needed to search faces
+        (
+            false,                      // do not cache bb
+            mesh0,
+            bndFaces                    // boundary faces only
+        ),
+        overallBb.extend(rndGen, 1E-4), // overall search domain
+        8,                              // maxLevel
+        10,                             // leafsize
+        3.0                             // duplicity
+    );
 
     if (debug)
     {
@@ -1048,17 +1058,11 @@ void Foam::faceCoupleInfo::findSlavesCoveringMaster
         // Generate face centre (prevent cellCentres() reconstruction)
         point fc(f1.centre(mesh1.points()));
 
-        // Search in bounding box of face only.
-        treeBoundBox tightest(static_cast<const pointField&>(f1.points(mesh1.points())));
+        pointIndexHit nearInfo = tree.findNearest(fc, Foam::sqr(absTol));
 
-        scalar tightestDist = GREAT;
-
-        label index = tree.findNearest(fc, tightest, tightestDist);
-
-
-        if (index != -1)
+        if (nearInfo.hit())
         {
-            label mesh0FaceI = shapes.meshFaces()[index];
+            label mesh0FaceI = tree.shapes().faceLabels()[nearInfo.index()];
 
             // Check if points of f1 actually lie on top of mesh0 face
             // This is the bit that might fail since if f0 is severely warped

src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C

@@ -62,6 +62,7 @@ namespace Foam
     };
 }
 
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 void Foam::hexRef8::reorder
@@ -80,7 +81,7 @@ void Foam::hexRef8::reorder
 
         if (newI >= len)
         {
-            FatalErrorIn("hexRef8::reorder") << abort(FatalError);
+            FatalErrorIn("hexRef8::reorder(..)") << abort(FatalError);
         }
 
         if (newI >= 0)
@@ -557,22 +558,11 @@ Foam::label Foam::hexRef8::getAnchorCell
         }
 
         // Problem.
-
-        // Pick up points of the cell
-        const labelList cPoints(cellPoints(cellI));
-
-        Perr<< "cell:" << cellI << " points:" << endl;
-        forAll(cPoints, i)
-        {
-            label pointI = cPoints[i];
-
-            Perr<< "    " << pointI << " coord:" << mesh_.points()[pointI]
-                << nl;
-        }
+        dumpCell(cellI);
         Perr<< "cell:" << cellI << " anchorPoints:" << cellAnchorPoints[cellI]
             << endl;
 
-        FatalErrorIn("hexRef8::getAnchorCell")
+        FatalErrorIn("hexRef8::getAnchorCell(..)")
             << "Could not find point " << pointI
             << " in the anchorPoints for cell " << cellI << endl
             << "Does your original mesh obey the 2:1 constraint and"
@@ -690,9 +680,50 @@ Foam::label Foam::hexRef8::countAnchors
 }
 
 
+void Foam::hexRef8::dumpCell(const label cellI) const
+{
+    OFstream str(mesh_.time().path()/"cell_" + Foam::name(cellI) + ".obj");
+    Pout<< "hexRef8 : Dumping cell as obj to " << str.name() << endl;
+
+    const cell& cFaces = mesh_.cells()[cellI];
+
+    Map<label> pointToObjVert;
+    label objVertI = 0;
+
+    forAll(cFaces, i)
+    {
+        const face& f = mesh_.faces()[cFaces[i]];
+
+        forAll(f, fp)
+        {
+            if (pointToObjVert.insert(f[fp], objVertI))
+            {
+                meshTools::writeOBJ(str, mesh_.points()[f[fp]]);
+                objVertI++;
+            }
+        }
+    }
+
+    forAll(cFaces, i)
+    {
+        const face& f = mesh_.faces()[cFaces[i]];
+
+        forAll(f, fp)
+        {
+            label pointI = f[fp];
+            label nexPointI = f[f.fcIndex(fp)];
+
+            str << "l " << pointToObjVert[pointI]+1
+                << ' ' << pointToObjVert[nexPointI]+1 << nl;
+        }
+    }
+}
+
+
 // Find point with certain pointLevel. Skip any higher levels.
 Foam::label Foam::hexRef8::findLevel
 (
+    const label faceI,
     const face& f,
     const label startFp,
     const bool searchForward,
@@ -707,7 +738,13 @@ Foam::label Foam::hexRef8::findLevel
 
         if (pointLevel_[pointI] < wantedLevel)
         {
-            FatalErrorIn("hexRef8::findLevel")
+            dumpCell(mesh_.faceOwner()[faceI]);
+            if (mesh_.isInternalFace(faceI))
+            {
+                dumpCell(mesh_.faceNeighbour()[faceI]);
+            }
+
+            FatalErrorIn("hexRef8::findLevel(..)")
                 << "face:" << f
                 << " level:" << IndirectList<label>(pointLevel_, f)()
                 << " startFp:" << startFp
@@ -729,7 +766,13 @@ Foam::label Foam::hexRef8::findLevel
         }
     }
 
-    FatalErrorIn("hexRef8::findLevel")
+    dumpCell(mesh_.faceOwner()[faceI]);
+    if (mesh_.isInternalFace(faceI))
+    {
+        dumpCell(mesh_.faceNeighbour()[faceI]);
+    }
+
+    FatalErrorIn("hexRef8::findLevel(..)")
         << "face:" << f
         << " level:" << IndirectList<label>(pointLevel_, f)()
         << " startFp:" << startFp
@@ -763,15 +806,6 @@ Foam::label Foam::hexRef8::getAnchorLevel(const label faceI) const
         }
         else
         {
-            //FatalErrorIn("hexRef8::getAnchorLevel(const label) const")
-            //    << "face:" << faceI
-            //    << " centre:" << mesh_.faceCentres()[faceI]
-            //    << " verts:" << f
-            //    << " point levels:" << IndirectList<label>(pointLevel_, f)()
-            //    << " own:" << mesh_.faceOwner()[faceI]
-            //    << " ownLevel:" << cellLevel_[mesh_.faceOwner()[faceI]]
-            //    << abort(FatalError);
-
             return -1;
         }
     }
@@ -797,7 +831,7 @@ void Foam::hexRef8::checkInternalOrientation
 
     if ((dir & n) < 0)
     {
-        FatalErrorIn("checkInternalOrientation")
+        FatalErrorIn("checkInternalOrientation(..)")
             << "cell:" << cellI << " old face:" << faceI
             << " newFace:" << newFace << endl
             << " coords:" << compactPoints
@@ -812,7 +846,7 @@ void Foam::hexRef8::checkInternalOrientation
 
     if (s < 0.1 || s > 0.9)
     {
-        FatalErrorIn("checkInternalOrientation")
+        FatalErrorIn("checkInternalOrientation(..)")
             << "cell:" << cellI << " old face:" << faceI
             << " newFace:" << newFace << endl
             << " coords:" << compactPoints
@@ -843,7 +877,7 @@ void Foam::hexRef8::checkBoundaryOrientation
 
     if ((dir & n) < 0)
     {
-        FatalErrorIn("checkBoundaryOrientation")
+        FatalErrorIn("checkBoundaryOrientation(..)")
             << "cell:" << cellI << " old face:" << faceI
             << " newFace:" << newFace
             << " coords:" << compactPoints
@@ -858,7 +892,7 @@ void Foam::hexRef8::checkBoundaryOrientation
 
     if (s < 0.7 || s > 1.3)
     {
-        WarningIn("checkBoundaryOrientation")
+        WarningIn("checkBoundaryOrientation(..)")
             << "cell:" << cellI << " old face:" << faceI
             << " newFace:" << newFace
             << " coords:" << compactPoints
@@ -1191,8 +1225,22 @@ void Foam::hexRef8::createInternalFaces
             }
 
             // Now the face mid point is the second cLevel+1 point
-            label edgeMid = findLevel(f, f.fcIndex(anchorFp), true, cLevel+1);
-            label faceMid = findLevel(f, f.fcIndex(edgeMid), true, cLevel+1);
+            label edgeMid = findLevel
+            (
+                faceI,
+                f,
+                f.fcIndex(anchorFp),
+                true,
+                cLevel+1
+            );
+            label faceMid = findLevel
+            (
+                faceI,
+                f,
+                f.fcIndex(edgeMid),
+                true,
+                cLevel+1
+            );
 
             faceMidPointI = f[faceMid];
         }
@@ -1205,7 +1253,13 @@ void Foam::hexRef8::createInternalFaces
         }
         else
         {
-            FatalErrorIn("createInternalFaces")
+            dumpCell(mesh_.faceOwner()[faceI]);
+            if (mesh_.isInternalFace(faceI))
+            {
+                dumpCell(mesh_.faceNeighbour()[faceI]);
+            }
+
+            FatalErrorIn("createInternalFaces(..)")
                 << "nAnchors:" << nAnchors
                 << " faceI:" << faceI
                 << abort(FatalError);
@@ -1243,9 +1297,11 @@ void Foam::hexRef8::createInternalFaces
 
                     if (edgeMidPointI == -1)
                     {
+                        dumpCell(cellI);
+
                         const labelList cPoints(cellPoints(cellI));
 
-                        FatalErrorIn("createInternalFaces")
+                        FatalErrorIn("createInternalFaces(..)")
                             << "cell:" << cellI << " cLevel:" << cLevel
                             << " cell points:" << cPoints
                             << " pointLevel:"
@@ -1264,7 +1320,7 @@ void Foam::hexRef8::createInternalFaces
                 else
                 {
                     // Search forward in face to clevel+1
-                    label edgeMid = findLevel(f, fp1, true, cLevel+1);
+                    label edgeMid = findLevel(faceI, f, fp1, true, cLevel+1);
 
                     edgeMidPointI = f[edgeMid];
                 }
@@ -1314,9 +1370,11 @@ void Foam::hexRef8::createInternalFaces
 
                     if (edgeMidPointI == -1)
                     {
+                        dumpCell(cellI);
+
                         const labelList cPoints(cellPoints(cellI));
 
-                        FatalErrorIn("createInternalFaces")
+                        FatalErrorIn("createInternalFaces(..)")
                             << "cell:" << cellI << " cLevel:" << cLevel
                             << " cell points:" << cPoints
                             << " pointLevel:"
@@ -1335,7 +1393,14 @@ void Foam::hexRef8::createInternalFaces
                 else
                 {
                     // Search back to clevel+1
-                    label edgeMid = findLevel(f, fpMin1, false, cLevel+1);
+                    label edgeMid = findLevel
+                    (
+                        faceI,
+                        f,
+                        fpMin1,
+                        false,
+                        cLevel+1
+                    );
 
                     edgeMidPointI = f[edgeMid];
                 }
@@ -1591,9 +1656,11 @@ void Foam::hexRef8::checkWantedRefinementLevels
 
         if (mag(ownLevel-neiLevel) > 1)
         {
+            dumpCell(own);
+            dumpCell(nei);
             FatalErrorIn
             (
-                "hexRef8::checkRefinementLevels(const labelList&)"
+                "hexRef8::checkWantedRefinementLevels(const labelList&)"
             )   << "cell:" << own
                 << " current level:" << cellLevel_[own]
                 << " level after refinement:" << ownLevel
@@ -1633,9 +1700,10 @@ void Foam::hexRef8::checkWantedRefinementLevels
         {
             label patchI = mesh_.boundaryMesh().whichPatch(faceI);
 
+            dumpCell(own);
             FatalErrorIn
             (
-                "hexRef8::checkRefinementLevels(const labelList&)"
+                "hexRef8::checkWantedRefinementLevels(const labelList&)"
             )   << "Celllevel does not satisfy 2:1 constraint."
                 << " On coupled face "
                 << faceI
@@ -1730,6 +1798,25 @@ Foam::hexRef8::hexRef8(const polyMesh& mesh)
             << abort(FatalError);
     }
 
+    if
+    (
+        cellLevel_.size() != mesh_.nCells()
+     || pointLevel_.size() != mesh_.nPoints()
+    )
+    {
+        FatalErrorIn
+        (
+            "hexRef8::hexRef8(const polyMesh&)"
+        )   << "Restarted from inconsistent cellLevel or pointLevel files."
+            << endl
+            << "Number of cells in mesh:" << mesh_.nCells()
+            << " does not equal size of cellLevel:" << cellLevel_.size() << endl
+            << "Number of points in mesh:" << mesh_.nPoints()
+            << " does not equal size of pointLevel:" << pointLevel_.size()
+            << abort(FatalError);
+    }
+
+
     // Check refinement levels for consistency
     checkRefinementLevels(-1, labelList(0));
 
@@ -1806,7 +1893,24 @@ Foam::hexRef8::hexRef8
         )   << "History enabled but number of visible cells in it "
             << history_.visibleCells().size()
             << " is not equal to the number of cells in the mesh "
-            << mesh_.nCells()
+            << mesh_.nCells() << abort(FatalError);
+    }
+
+    if
+    (
+        cellLevel_.size() != mesh_.nCells()
+     || pointLevel_.size() != mesh_.nPoints()
+    )
+    {
+        FatalErrorIn
+        (
+            "hexRef8::hexRef8(const polyMesh&, const labelList&"
+            ", const labelList&, const refinementHistory&)"
+        )   << "Incorrect cellLevel or pointLevel size." << endl
+            << "Number of cells in mesh:" << mesh_.nCells()
+            << " does not equal size of cellLevel:" << cellLevel_.size() << endl
+            << "Number of points in mesh:" << mesh_.nPoints()
+            << " does not equal size of pointLevel:" << pointLevel_.size()
             << abort(FatalError);
     }
 
@@ -1877,6 +1981,24 @@ Foam::hexRef8::hexRef8
     savedPointLevel_(0),
     savedCellLevel_(0)
 {
+    if
+    (
+        cellLevel_.size() != mesh_.nCells()
+     || pointLevel_.size() != mesh_.nPoints()
+    )
+    {
+        FatalErrorIn
+        (
+            "hexRef8::hexRef8(const polyMesh&, const labelList&"
+            ", const labelList&)"
+        )   << "Incorrect cellLevel or pointLevel size." << endl
+            << "Number of cells in mesh:" << mesh_.nCells()
+            << " does not equal size of cellLevel:" << cellLevel_.size() << endl
+            << "Number of points in mesh:" << mesh_.nPoints()
+            << " does not equal size of pointLevel:" << pointLevel_.size()
+            << abort(FatalError);
+    }
+
     // Check refinement levels for consistency
     checkRefinementLevels(-1, labelList(0));
 
@@ -2356,9 +2478,13 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement
 
             if (mag(ownLevel-neiLevel) > 1)
             {
+                dumpCell(own);
+                dumpCell(nei);
                 FatalErrorIn
                 (
                     "hexRef8::consistentSlowRefinement"
+                    "(const label, const labelList&, const labelList&"
+                    ", const label, const labelList&)"
                 )   << "cell:" << own
                     << " current level:" << cellLevel_[own]
                     << " current refData:" << allCellInfo[own]
@@ -2412,11 +2538,14 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement
 
             if (mag(ownLevel - nbrLevel) > 1)
             {
+                dumpCell(own);
                 label patchI = mesh_.boundaryMesh().whichPatch(faceI);
 
                 FatalErrorIn
                 (
-                    "hexRef8::checkRefinementLevels(const labelList&)"
+                    "hexRef8::consistentSlowRefinement"
+                    "(const label, const labelList&, const labelList&"
+                    ", const label, const labelList&)"
                 )   << "Celllevel does not satisfy 2:1 constraint."
                     << " On coupled face "
                     << faceI
@@ -2805,7 +2934,6 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement2
             << cellsToRefine.size() << " to " << newCellsToRefine.size()
             << " cells to refine." << endl;
 
-//XXXX
         // Check that newCellsToRefine obeys at least 2:1.
 
         {
@@ -2861,6 +2989,7 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement2
                  && savedRefineCell.get(cellI) == 0
                 )
                 {
+                    dumpCell(cellI);
                     FatalErrorIn
                     (
                         "hexRef8::consistentSlowRefinement2"
@@ -2873,7 +3002,6 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement2
                 }
             }
         }
-//XXXX
     }
 
     return newCellsToRefine;
@@ -3342,6 +3470,7 @@ Foam::labelListList Foam::hexRef8::setRefinement
                 {
                     if (nAnchorPoints[cellI] == 8)
                     {
+                        dumpCell(cellI);
                         FatalErrorIn
                         (
                             "hexRef8::setRefinement(const labelList&"
@@ -3365,6 +3494,8 @@ Foam::labelListList Foam::hexRef8::setRefinement
             {
                 if (nAnchorPoints[cellI] != 8)
                 {
+                    dumpCell(cellI);
+
                     const labelList cPoints(cellPoints(cellI));
 
                     FatalErrorIn
@@ -3841,7 +3972,7 @@ Foam::labelListList Foam::hexRef8::setRefinement
 
         if (minPointI != labelMax && minPointI != mesh_.nPoints())
         {
-            FatalErrorIn("hexRef8::setRefinement")
+            FatalErrorIn("hexRef8::setRefinement(..)")
                 << "Added point labels not consecutive to existing mesh points."
                 << nl
                 << "mesh_.nPoints():" << mesh_.nPoints()
@@ -3995,13 +4126,6 @@ void Foam::hexRef8::updateMesh
 
                 if (oldCellI == -1)
                 {
-                    //FatalErrorIn("hexRef8::updateMesh(const mapPolyMesh&)")
-                    //    << "Problem : cell " << newCellI
-                    //    << " at " << mesh_.cellCentres()[newCellI]
-                    //    << " does not originate from another cell"
-                    //    << " (i.e. is inflated)." << nl
-                    //    << "Hence we cannot determine the new cellLevel"
-                    //    << " for it." << abort(FatalError);
                     newCellLevel[newCellI] = -1;
                 }
                 else
@@ -4036,8 +4160,8 @@ void Foam::hexRef8::updateMesh
         //{
         //    WarningIn("hexRef8::updateMesh(const mapPolyMesh&)")
         //        << "Problem : "
-        //        << "cellLevel_ contains illegal value -1 after mapping at cell "
-        //        << findIndex(cellLevel_, -1) << endl
+        //        << "cellLevel_ contains illegal value -1 after mapping
+        //        << " at cell " << findIndex(cellLevel_, -1) << endl
         //        << "This means that another program has inflated cells"
         //        << " (created cells out-of-nothing) and hence we don't know"
         //        << " their cell level. Continuing with illegal value."
@@ -4174,7 +4298,7 @@ void Foam::hexRef8::subset
 
         if (findIndex(cellLevel_, -1) != -1)
         {
-            FatalErrorIn("hexRef8::subset")
+            FatalErrorIn("hexRef8::subset(..)")
                 << "Problem : "
                 << "cellLevel_ contains illegal value -1 after mapping:"
                 << cellLevel_
@@ -4195,7 +4319,7 @@ void Foam::hexRef8::subset
 
         if (findIndex(pointLevel_, -1) != -1)
         {
-            FatalErrorIn("hexRef8::subset")
+            FatalErrorIn("hexRef8::subset(..)")
                 << "Problem : "
                 << "pointLevel_ contains illegal value -1 after mapping:"
                 << pointLevel_
@@ -4292,6 +4416,7 @@ void Foam::hexRef8::checkMesh() const
 
                     if (!cellToFace.insert(labelPair(own, nei[bFaceI]), faceI))
                     {
+                        dumpCell(own);
                         FatalErrorIn("hexRef8::checkMesh()")
                             << "Faces do not seem to be correct across coupled"
                             << " boundaries" << endl
@@ -4335,6 +4460,8 @@ void Foam::hexRef8::checkMesh() const
                 const face& f = mesh_.faces()[faceI];
                 label patchI = mesh_.boundaryMesh().whichPatch(faceI);
 
+                dumpCell(mesh_.faceOwner()[faceI]);
+
                 FatalErrorIn("hexRef8::checkMesh()")
                     << "Faces do not seem to be correct across coupled"
                     << " boundaries" << endl
@@ -4373,6 +4500,8 @@ void Foam::hexRef8::checkMesh() const
 
             if (f.size() != nVerts[i])
             {
+                dumpCell(mesh_.faceOwner()[faceI]);
+
                 label patchI = mesh_.boundaryMesh().whichPatch(faceI);
 
                 FatalErrorIn("hexRef8::checkMesh()")
@@ -4421,6 +4550,8 @@ void Foam::hexRef8::checkMesh() const
 
             if (mag(anchorVec - anchorPoints[i]) > smallDim)
             {
+                dumpCell(mesh_.faceOwner()[faceI]);
+
                 label patchI = mesh_.boundaryMesh().whichPatch(faceI);
 
                 FatalErrorIn("hexRef8::checkMesh()")
@@ -4487,6 +4618,9 @@ void Foam::hexRef8::checkRefinementLevels
 
             if (mag(cellLevel_[own] - cellLevel_[nei]) > 1)
             {
+                dumpCell(own);
+                dumpCell(nei);
+
                 FatalErrorIn
                 (
                     "hexRef8::checkRefinementLevels(const label)"
@@ -4520,6 +4654,8 @@ void Foam::hexRef8::checkRefinementLevels
 
             if (mag(cellLevel_[own] - neiLevel[i]) > 1)
             {
+                dumpCell(own);
+
                 label patchI = mesh_.boundaryMesh().whichPatch(faceI);
 
                 FatalErrorIn
@@ -4621,6 +4757,8 @@ void Foam::hexRef8::checkRefinementLevels
                   > maxPointDiff
                 )
                 {
+                    dumpCell(cellI);
+
                     FatalErrorIn
                     (
                         "hexRef8::checkRefinementLevels(const label)"
@@ -4637,6 +4775,72 @@ void Foam::hexRef8::checkRefinementLevels
             }
         }
     }
+
+
+    // Hanging points
+    {
+        // Any patches with points having only two edges.
+
+        boolList isHangingPoint(mesh_.nPoints(), false);
+
+        const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+
+        forAll(patches, patchI)
+        {
+            const polyPatch& pp = patches[patchI];
+
+            const labelList& meshPoints = pp.meshPoints();
+
+            forAll(meshPoints, i)
+            {
+                label pointI = meshPoints[i];
+
+                const labelList& pEdges = mesh_.pointEdges()[pointI];
+
+                if (pEdges.size() == 2)
+                {
+                    isHangingPoint[pointI] = true;
+                }
+            }
+        }
+
+        syncTools::syncPointList
+        (
+            mesh_,
+            isHangingPoint,
+            andEqOp<bool>(),        // only if all decide it is hanging point
+            true,                   // null
+            false                   // no separation
+        );
+
+        OFstream str(mesh_.time().path()/"hangingPoints.obj");
+
+        label nHanging = 0;
+
+        forAll(isHangingPoint, pointI)
+        {
+            if (isHangingPoint[pointI])
+            {
+                nHanging++;
+
+                Pout<< "Hanging boundary point " << pointI
+                    << " at " << mesh_.points()[pointI]
+                    << endl;
+
+                meshTools::writeOBJ(str, mesh_.points()[pointI]);
+            }
+        }
+
+        if (returnReduce(nHanging, sumOp<label>()) > 0)
+        {
+            FatalErrorIn
+            (
+                "hexRef8::checkRefinementLevels(const label)"
+            )   << "Detected a point used by two edges only (hanging point)"
+                << nl << "This is not allowed"
+                << abort(FatalError);
+        }
+    }
 }
 
 
@@ -4865,8 +5069,10 @@ Foam::labelList Foam::hexRef8::consistentUnrefinement
 
     if (maxSet)
     {
-        FatalErrorIn("hexRef8::consistentUnrefinement")
-            << "maxSet not implemented yet."
+        FatalErrorIn
+        (
+            "hexRef8::consistentUnrefinement(const labelList&, const bool"
+        )   << "maxSet not implemented yet."
             << abort(FatalError);
     }
 
@@ -4935,7 +5141,7 @@ Foam::labelList Foam::hexRef8::consistentUnrefinement
                 {
                     if (unrefineCell.get(own) == 0)
                     {
-                        FatalErrorIn("hexRef8::consistentUnrefinement")
+                        FatalErrorIn("hexRef8::consistentUnrefinement(..)")
                             << "problem" << abort(FatalError);
                     }
 
@@ -4953,7 +5159,7 @@ Foam::labelList Foam::hexRef8::consistentUnrefinement
                 {
                     if (unrefineCell.get(nei) == 0)
                     {
-                        FatalErrorIn("hexRef8::consistentUnrefinement")
+                        FatalErrorIn("hexRef8::consistentUnrefinement(..)")
                             << "problem" << abort(FatalError);
                     }
 
@@ -4989,7 +5195,7 @@ Foam::labelList Foam::hexRef8::consistentUnrefinement
                 {
                     if (unrefineCell.get(own) == 0)
                     {
-                        FatalErrorIn("hexRef8::consistentUnrefinement")
+                        FatalErrorIn("hexRef8::consistentUnrefinement(..)")
                             << "problem" << abort(FatalError);
                     }
 
@@ -5003,7 +5209,7 @@ Foam::labelList Foam::hexRef8::consistentUnrefinement
                 {
                     if (unrefineCell.get(own) == 1)
                     {
-                        FatalErrorIn("hexRef8::consistentUnrefinement")
+                        FatalErrorIn("hexRef8::consistentUnrefinement(..)")
                             << "problem" << abort(FatalError);
                     }
 
@@ -5086,8 +5292,10 @@ void Foam::hexRef8::setUnrefinement
 {
     if (!history_.active())
     {
-        FatalErrorIn("hexRef8::setUnrefinement()")
-            << "Only call if constructed with history capability"
+        FatalErrorIn
+        (
+            "hexRef8::setUnrefinement(const labelList&, polyTopoChange&)"
+        )   << "Only call if constructed with history capability"
             << abort(FatalError);
     }
 
@@ -5102,8 +5310,10 @@ void Foam::hexRef8::setUnrefinement
         {
             if (cellLevel_[cellI] < 0)
             {
-                FatalErrorIn("hexRef8::setUnrefinement()")
-                    << "Illegal cell level " << cellLevel_[cellI]
+                FatalErrorIn
+                (
+                    "hexRef8::setUnrefinement(const labelList&, polyTopoChange&)"
+                )   << "Illegal cell level " << cellLevel_[cellI]
                     << " for cell " << cellI
                     << abort(FatalError);
             }
@@ -5165,8 +5375,10 @@ void Foam::hexRef8::setUnrefinement
 
         if (facesToRemove.size() != splitFaces.size())
         {
-            FatalErrorIn("hexRef8::setUnrefinement")
-                << "Ininitial set of split points to unrefine does not"
+            FatalErrorIn
+            (
+                "hexRef8::setUnrefinement(const labelList&, polyTopoChange&)"
+            )   << "Ininitial set of split points to unrefine does not"
                 << " seem to be consistent or not mid points of refined cells"
                 << abort(FatalError);
         }
@@ -5187,8 +5399,10 @@ void Foam::hexRef8::setUnrefinement
         // Check
         if (pCells.size() != 8)
         {
-            FatalErrorIn("hexRef8::setUnrefinement")
-                << "splitPoint " << pointI
+            FatalErrorIn
+            (
+                "hexRef8::setUnrefinement(const labelList&, polyTopoChange&)"
+            )   << "splitPoint " << pointI
                 << " should have 8 cells using it. It has " << pCells
                 << abort(FatalError);
         }
@@ -5208,7 +5422,7 @@ void Foam::hexRef8::setUnrefinement
 
                 if (region == -1)
                 {
-                    FatalErrorIn("hexRef8::setUnrefinement")
+                    FatalErrorIn("hexRef8::setUnrefinement(..)")
                         << "Ininitial set of split points to unrefine does not"
                         << " seem to be consistent or not mid points"
                         << " of refined cells" << nl
@@ -5219,7 +5433,7 @@ void Foam::hexRef8::setUnrefinement
 
                 if (masterCellI != cellRegionMaster[region])
                 {
-                    FatalErrorIn("hexRef8::setUnrefinement")
+                    FatalErrorIn("hexRef8::setUnrefinement(..)")
                         << "cell:" << cellI << " on splitPoint:" << pointI
                         << " in region " << region
                         << " has master:" << cellRegionMaster[region]

src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.H

@@ -174,9 +174,12 @@ class hexRef8
         label findMaxLevel(const labelList& f) const;
         //- Count number of vertices <= anchorLevel
         label countAnchors(const labelList&, const label) const;
+        //- Debugging: dump cell as .obj file
+        void dumpCell(const label cellI) const;
         //- Find index of point with wantedLevel, starting from fp.
         label findLevel
         (
+            const label faceI,
             const face& f,
             const label startFp,
             const bool searchForward,

src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C

@@ -1295,7 +1295,7 @@ void Foam::polyTopoChange::calcFaceInflationMaps
                     selectFaces
                     (
                         mesh,
-                        mesh.edgeFaces()[iter()],
+                        mesh.edgeFaces(iter()),
                         true
                     )
                 );
@@ -1309,7 +1309,7 @@ void Foam::polyTopoChange::calcFaceInflationMaps
                     selectFaces
                     (
                         mesh,
-                        mesh.edgeFaces()[iter()],
+                        mesh.edgeFaces(iter()),
                         false
                     )
                 );

src/dynamicMesh/polyTopoChange/polyTopoChange/removeFaces.C

@@ -807,7 +807,7 @@ void Foam::removeFaces::setRefinement
     // Edges to remove
     labelHashSet edgesToRemove(faceLabels.size());
 
-    // Per face the region it is. -1 for removed faces, -2 for regions
+    // Per face the region it is in. -1 for removed faces, -2 for regions
     // consisting of single face only.
     labelList faceRegion(mesh_.nFaces(), -1);
 
@@ -1258,10 +1258,15 @@ void Foam::removeFaces::setRefinement
     // are only used by 2 unremoved edges.
     {
         // Usage of points by non-removed edges.
-        labelList nEdgesPerPoint(mesh_.nPoints(), labelMax);
+        labelList nEdgesPerPoint(mesh_.nPoints());
 
         const labelListList& pointEdges = mesh_.pointEdges();
 
+        forAll(pointEdges, pointI)
+        {
+            nEdgesPerPoint[pointI] = pointEdges[pointI].size();
+        }
+
         forAllConstIter(labelHashSet, edgesToRemove, iter)
         {
             // Edge will get removed.
@@ -1269,16 +1274,7 @@ void Foam::removeFaces::setRefinement
 
             forAll(e, i)
             {
-                label pointI = e[i];
-
-                if (nEdgesPerPoint[pointI] == labelMax)
-                {
-                    nEdgesPerPoint[pointI] = pointEdges[pointI].size()-1;
-                }
-                else
-                {
-                    nEdgesPerPoint[pointI]--;
-                }
+                nEdgesPerPoint[e[i]]--;
             }
         }
 

src/finiteVolume/Make/files

@@ -186,6 +186,8 @@ $(limitedSchemes)/MC/MC.C
 $(limitedSchemes)/Phi/Phi.C
 $(limitedSchemes)/filteredLinear/filteredLinear.C
 $(limitedSchemes)/filteredLinear2/filteredLinear2.C
+$(limitedSchemes)/filteredLinear3/filteredLinear3.C
+$(limitedSchemes)/limitWith/limitWith.C
 
 multivariateSchemes = $(surfaceInterpolation)/multivariateSchemes
 $(multivariateSchemes)/multivariateSurfaceInterpolationScheme/multivariateSurfaceInterpolationSchemes.C

src/finiteVolume/fields/fvPatchFields/derived/directMappedFixedValue/directMappedFixedValueFvPatchField.C

@@ -159,7 +159,7 @@ void directMappedFixedValueFvPatchField<Type>::updateCoeffs()
 
     // Get the scheduling information
     const List<labelPair>& schedule = mpp.schedule();
-    const labelListList& sendCellLabels = mpp.sendCellLabels();
+    const labelListList& sendLabels = mpp.sendLabels();
     const labelListList& receiveFaceLabels = mpp.receiveFaceLabels();
 
 
@@ -177,7 +177,7 @@ void directMappedFixedValueFvPatchField<Type>::updateCoeffs()
             toProc<< IndirectList<Type>
             (
                 this->internalField(),
-                sendCellLabels[recvProc]
+                sendLabels[recvProc]
             )();
         }
         else
@@ -204,7 +204,7 @@ void directMappedFixedValueFvPatchField<Type>::updateCoeffs()
         IndirectList<Type> fromFld
         (
             this->internalField(),
-            sendCellLabels[Pstream::myProcNo()]
+            sendLabels[Pstream::myProcNo()]
         );
     
         // Destination faces

src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/filteredLinear3/filteredLinear3.C

@@ -0,0 +1,48 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "LimitedScheme.H"
+#include "filteredLinear3.H"
+#include "filteredLinear3V.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    makeLimitedSurfaceInterpolationScheme
+    (
+        filteredLinear3,
+        filteredLinear3Limiter
+    )
+
+    makeLimitedVSurfaceInterpolationScheme
+    (
+        filteredLinear3V,
+        filteredLinear3VLimiter
+    )
+}
+
+// ************************************************************************* //

src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/filteredLinear3/filteredLinear3.H

@@ -0,0 +1,120 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::filteredLinear3Limiter
+
+Description
+    Class to generate weighting factors for the filtered-linear-3
+    differencing scheme.
+
+    The aim is to remove high-frequency modes with "staggering"
+    characteristics by comparing the face gradient with both neighbouring
+    cell gradients and introduce small amounts of upwind in order to damp
+    these modes.
+
+    Used in conjunction with the template class LimitedScheme.
+
+SourceFiles
+    filteredLinear3.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef filteredLinear3_H
+#define filteredLinear3_H
+
+#include "vector.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                     Class filteredLinear3Limiter Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class LimiterFunc>
+class filteredLinear3Limiter
+:
+    public LimiterFunc
+{
+    // Private data
+
+        // Scaling corefficient for the gradient ratio,
+        // 0 = linear
+        // 1 = fully limited
+        scalar k_;
+
+public:
+
+    filteredLinear3Limiter(Istream& is)
+    :
+        k_(readScalar(is))
+    {
+        if (k_ < 0 || k_ > 1)
+        {
+            FatalIOErrorIn("filteredLinear3Limiter(Istream& is)", is)
+                << "coefficient = " << k_
+                << " should be >= 0 and <= 1"
+                << exit(FatalIOError);
+        }    
+    }
+
+    scalar limiter
+    (
+        const scalar cdWeight,
+        const scalar faceFlux,
+        const typename LimiterFunc::phiType& phiP,
+        const typename LimiterFunc::phiType& phiN,
+        const typename LimiterFunc::gradPhiType& gradcP,
+        const typename LimiterFunc::gradPhiType& gradcN,
+        const vector& d
+    ) const
+    {
+        // Difference across face
+        scalar df = phiN - phiP;
+
+        // Twice the differences across face-neighbour cells
+        scalar dP = 2*(d & gradcP);
+        scalar dN = 2*(d & gradcN);
+
+        // Calculate the limiter
+        scalar limiter = 1 - k_*(dN - df)*(dP - df)/max(sqr(dN + dP), SMALL);
+
+        // Limit the limiter between linear and upwind
+        return max(min(limiter, 1), 0);
+    }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/filteredLinear3/filteredLinear3V.H

@@ -0,0 +1,123 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::filteredLinear3VLimiter
+
+Description
+    Class to generate weighting factors for the filteredLinear3V differencing
+    scheme.  The aim is to remove high-frequency modes with "staggering"
+    characteristics from vector fields by comparing the face gradient in the
+    direction of maximum gradient with both neighbouring cell gradients and
+    introduce small amounts of upwind in order to damp these modes.
+
+    Used in conjunction with the template class LimitedScheme.
+
+SourceFiles
+    filteredLinear3V.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef filteredLinear3V_H
+#define filteredLinear3V_H
+
+#include "vector.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                 Class filteredLinear3VLimiter Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class LimiterFunc>
+class filteredLinear3VLimiter
+:
+    public LimiterFunc
+{
+    // Private data
+
+        // Scaling corefficient for the gradient ratio,
+        // 0 = linear
+        // 1 = fully limited
+        scalar k_;
+
+public:
+
+    filteredLinear3VLimiter(Istream& is)
+    :
+        k_(readScalar(is))
+    {
+        if (k_ < 0 || k_ > 1)
+        {
+            FatalIOErrorIn("filteredLinear3VLimiter(Istream& is)", is)
+                << "coefficient = " << k_
+                << " should be >= 0 and <= 1"
+                << exit(FatalIOError);
+        }    
+    }
+
+    scalar limiter
+    (
+        const scalar cdWeight,
+        const scalar faceFlux,
+        const typename LimiterFunc::phiType& phiP,
+        const typename LimiterFunc::phiType& phiN,
+        const typename LimiterFunc::gradPhiType& gradcP,
+        const typename LimiterFunc::gradPhiType& gradcN,
+        const vector& d
+    ) const
+    {
+        // Difference across face
+        vector dfV = phiN - phiP;
+
+        // Scalar difference across the face
+        // in the direction in which the difference is largest
+        scalar df = dfV & dfV;
+
+        // Twice differences across face-neighbour cells
+        // in the direction in which the face-difference is largest
+        scalar dP = 2*(dfV & (d & gradcP));
+        scalar dN = 2*(dfV & (d & gradcN));
+
+        // Calculate the limiter
+        scalar limiter = 1 - k_*(dN - df)*(dP - df)/max(sqr(dN + dP), SMALL);
+
+        // Limit the limiter between linear and upwind
+        return max(min(limiter, 1), 0);
+    }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/limitWith/limitWith.C

@@ -0,0 +1,37 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvMesh.H"
+#include "limitWith.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    makeSurfaceInterpolationScheme(limitWith)
+}
+
+// ************************************************************************* //

src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/limitWith/limitWith.H

@@ -0,0 +1,161 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::limitWith
+
+Description
+    limitWith differencing scheme limits the specified scheme with the
+    specified limiter.
+
+SourceFiles
+    limitWith.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef limitWith_H
+#define limitWith_H
+
+#include "surfaceInterpolationScheme.H"
+#include "limitedSurfaceInterpolationScheme.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class limitWith Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class limitWith
+:
+    public surfaceInterpolationScheme<Type>
+{
+    // Private Member Functions
+
+        //- Interpolation scheme
+        tmp<surfaceInterpolationScheme<Type> > tInterp_;
+
+        //- Limiter
+        tmp<limitedSurfaceInterpolationScheme<Type> > tLimiter_;
+
+
+        //- Disallow default bitwise copy construct
+        limitWith(const limitWith&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const limitWith&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("limitWith");
+
+
+    // Constructors
+
+        //- Construct from mesh and Istream.
+        //  The name of the flux field is read from the Istream and looked-up
+        //  from the mesh objectRegistry
+        limitWith
+        (
+            const fvMesh& mesh,
+            Istream& is
+        )
+        :
+            surfaceInterpolationScheme<Type>(mesh),
+            tInterp_
+            (
+                surfaceInterpolationScheme<Type>::New(mesh, is)
+            ),
+            tLimiter_
+            (
+                limitedSurfaceInterpolationScheme<Type>::New(mesh, is)
+            )
+        {}
+
+        //- Construct from mesh, faceFlux and Istream
+        limitWith
+        (
+            const fvMesh& mesh,
+            const surfaceScalarField& faceFlux,
+            Istream& is
+        )
+        :
+            surfaceInterpolationScheme<Type>(mesh),
+            tInterp_
+            (
+                surfaceInterpolationScheme<Type>::New(mesh, faceFlux, is)
+            ),
+            tLimiter_
+            (
+                limitedSurfaceInterpolationScheme<Type>::New(mesh, faceFlux, is)
+            )
+        {}
+
+
+    // Member Functions
+
+        //- Return the interpolation weighting factors
+        virtual tmp<surfaceScalarField> weights
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& vf
+        ) const
+        {
+            return tLimiter_().weights
+            (
+                vf,
+                tInterp_().weights(vf),
+                tLimiter_().limiter(vf)
+            );
+        }
+
+        //- Return true if this scheme uses an explicit correction
+        virtual bool corrected() const
+        {
+            return tInterp_().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
+        {
+            return tLimiter_().limiter(vf)*tInterp_().correction(vf);
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/limitedSurfaceInterpolationScheme/limitedSurfaceInterpolationScheme.C

@@ -22,9 +22,6 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Description
-    Abstract base class for surface interpolation schemes.
-
 \*---------------------------------------------------------------------------*/
 
 #include "limitedSurfaceInterpolationScheme.H"
@@ -39,7 +36,6 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * //
 
-// Return weighting factors for scheme given by name in dictionary
 template<class Type>
 tmp<limitedSurfaceInterpolationScheme<Type> >
 limitedSurfaceInterpolationScheme<Type>::New
@@ -160,17 +156,14 @@ limitedSurfaceInterpolationScheme<Type>::~limitedSurfaceInterpolationScheme()
 template<class Type>
 tmp<surfaceScalarField> limitedSurfaceInterpolationScheme<Type>::weights
 (
-    const GeometricField<Type, fvPatchField, volMesh>& phi
+    const GeometricField<Type, fvPatchField, volMesh>& phi,
+    const surfaceScalarField& CDweights,
+    tmp<surfaceScalarField> tLimiter
 ) const
 {
-    const fvMesh& mesh = this->mesh();
-
     // Note that here the weights field is initialised as the limiter
     // from which the weight is calculated using the limiter value
-    tmp<surfaceScalarField> tWeights(this->limiter(phi));
-    surfaceScalarField& Weights = tWeights();
-
-    const surfaceScalarField& CDweights = mesh.surfaceInterpolation::weights();
+    surfaceScalarField& Weights = tLimiter();
 
     scalarField& pWeights = Weights.internalField();
 
@@ -199,9 +192,22 @@ tmp<surfaceScalarField> limitedSurfaceInterpolationScheme<Type>::weights
         }
     }
 
-    return tWeights;
+    return tLimiter;
 }
 
+template<class Type>
+tmp<surfaceScalarField> limitedSurfaceInterpolationScheme<Type>::weights
+(
+    const GeometricField<Type, fvPatchField, volMesh>& phi
+) const
+{
+    return this->weights
+    (
+        phi,
+        this->mesh().surfaceInterpolation::weights(),
+        this->limiter(phi)
+    );
+}
 
 template<class Type>
 tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >

src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/limitedSurfaceInterpolationScheme/limitedSurfaceInterpolationScheme.H

@@ -26,7 +26,7 @@ Class
     Foam::limitedSurfaceInterpolationScheme
 
 Description
-    Abstract base class for surface interpolation schemes.
+    Abstract base class for limited surface interpolation schemes.
 
 SourceFiles
     limitedSurfaceInterpolationScheme.C
@@ -119,7 +119,7 @@ public:
         {}
 
 
-        //- Construct from mesh and Istream. 
+        //- Construct from mesh and Istream.
         //  The name of the flux field is read from the Istream and looked-up
         //  from the mesh objectRegistry
         limitedSurfaceInterpolationScheme
@@ -170,6 +170,15 @@ public:
             const GeometricField<Type, fvPatchField, volMesh>&
         ) const = 0;
 
+        //- Return the interpolation weighting factors for the given field,
+        //  by limiting the given weights with the given limiter
+        tmp<surfaceScalarField> weights
+        (
+            const GeometricField<Type, fvPatchField, volMesh>&,
+            const surfaceScalarField& CDweights,
+            tmp<surfaceScalarField> tLimiter
+        ) const;
+
         //- Return the interpolation weighting factors for the given field
         virtual tmp<surfaceScalarField> weights
         (

src/lagrangian/basic/Cloud/CloudIO.C

@@ -118,7 +118,8 @@ void Foam::Cloud<ParticleType>::checkFieldIOobject
             "void Cloud<ParticleType>::checkFieldIOobject"
             "(Cloud<ParticleType>, IOField<DataType>)"
         )   << "Size of " << data.name()
-            << " field does not match the number of particles"
+            << " field " << data.size()
+            << " does not match the number of particles " << c.size()
             << abort(FatalError);
     }
 }

src/meshTools/Make/options

@@ -3,5 +3,5 @@ EXE_INC = \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude
 
 LIB_LIBS = \
-   -ltriSurface \
-   -llagrangian
+    -ltriSurface \
+    -llagrangian

src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.C

@@ -40,6 +40,18 @@ namespace Foam
 
     addToRunTimeSelectionTable(polyPatch, directMappedPolyPatch, word);
     addToRunTimeSelectionTable(polyPatch, directMappedPolyPatch, dictionary);
+
+
+    template<>
+    const char* NamedEnum<directMappedPolyPatch::sampleMode, 3>::names[] =
+    {
+        "nearestCell",
+        "nearestPatchFace",
+        "nearestFace"
+    };
+
+    const NamedEnum<directMappedPolyPatch::sampleMode, 3>
+        directMappedPolyPatch::sampleModeNames_;
 }
 
 
@@ -113,107 +125,236 @@ void Foam::directMappedPolyPatch::collectSamples
 void Foam::directMappedPolyPatch::findSamples
 (
     const pointField& samples,
-    labelList& sampleCellProcs,
-    labelList& sampleCells,
-    pointField& sampleCc
+    labelList& sampleProcs,
+    labelList& sampleIndices,
+    pointField& sampleLocations
 ) const
 {
-    sampleCellProcs.setSize(samples.size());
-    sampleCells.setSize(samples.size());
-    sampleCc.setSize(samples.size());
-    sampleCc = point(-GREAT, -GREAT, -GREAT);
+    const polyMesh& mesh = boundaryMesh().mesh();
 
+    // All the info for nearest. Construct to miss
+    List<nearInfo> nearest(samples.size());
+
+    switch (mode_)
     {
-        // Octree based search engine
-        meshSearch meshSearchEngine(boundaryMesh().mesh(), false);
-
-        forAll(samples, sampleI)
+        case NEARESTCELL:
         {
-            sampleCells[sampleI] = meshSearchEngine.findCell(samples[sampleI]);
+            if (samplePatch_.size() != 0 && samplePatch_ != "none")
+            {
+                FatalErrorIn
+                (
+                    "directMappedPolyPatch::findSamples(const pointField&,"
+                    " labelList&, labelList&, pointField&) const"
+                )   << "No need to supply a patch name when in "
+                    << sampleModeNames_[mode_] << " mode." << exit(FatalError);
+            }
 
-            if (sampleCells[sampleI] == -1)
+            // Octree based search engine
+            meshSearch meshSearchEngine(mesh, false);
+
+            forAll(samples, sampleI)
             {
-                sampleCellProcs[sampleI] = -1;
+                const point& sample = samples[sampleI];
+
+                label cellI = meshSearchEngine.findCell(sample);
+
+                if (cellI == -1)
+                {
+                    nearest[sampleI].second().first() = Foam::sqr(GREAT);
+                    nearest[sampleI].second().second() = Pstream::myProcNo();
+                }
+                else
+                {
+                    const point& cc = mesh.cellCentres()[cellI];
+
+                    nearest[sampleI].first() = pointIndexHit
+                    (
+                        true,
+                        cc,
+                        cellI
+                    );
+                    nearest[sampleI].second().first() = magSqr(cc-sample);
+                    nearest[sampleI].second().second() = Pstream::myProcNo();
+                }
             }
-            else
+            break;
+        }
+
+        case NEARESTPATCHFACE:
+        {
+            label patchI = boundaryMesh().findPatchID(samplePatch_);
+
+            if (patchI == -1)
             {
-                sampleCellProcs[sampleI] = Pstream::myProcNo();
-                sampleCc[sampleI] =
-                    boundaryMesh().mesh().cellCentres()[sampleCells[sampleI]];
+                FatalErrorIn("directMappedPolyPatch::findSamples(..)")
+                    << "Cannot find patch " << samplePatch_ << endl
+                    << "Valid patches are " << boundaryMesh().names()
+                    << exit(FatalError);
             }
+
+            Random rndGen(123456);
+
+            const polyPatch& pp = boundaryMesh()[patchI];
+
+            // patch faces
+            const labelList patchFaces(identity(pp.size()) + pp.start());
+
+            const treeBoundBox patchBb
+            (
+                treeBoundBox(pp.points(), pp.meshPoints()).extend
+                (
+                    rndGen,
+                    1E-4
+                )
+            );
+
+            autoPtr<indexedOctree<treeDataFace> > boundaryTree
+            (
+                new indexedOctree<treeDataFace>
+                (
+                    treeDataFace    // all information needed to search faces
+                    (
+                        false,                      // do not cache bb
+                        mesh,
+                        patchFaces                  // boundary faces only
+                    ),
+                    patchBb,   // overall search domain
+                    8,                              // maxLevel
+                    10,                             // leafsize
+                    3.0                             // duplicity
+                )
+            );
+
+            forAll(samples, sampleI)
+            {
+                const point& sample = samples[sampleI];
+
+                pointIndexHit& nearInfo = nearest[sampleI].first();
+                nearInfo = boundaryTree().findNearest
+                (
+                    sample,
+                    magSqr(patchBb.max()-patchBb.min())
+                );
+
+                if (!nearInfo.hit())
+                {
+                    nearest[sampleI].second().first() = Foam::sqr(GREAT);
+                    nearest[sampleI].second().second() = Pstream::myProcNo();
+                }
+                else
+                {
+                    point fc(pp[nearInfo.index()].centre(pp.points()));
+                    nearInfo.setPoint(fc);
+                    nearest[sampleI].second().first() = magSqr(fc-sample);
+                    nearest[sampleI].second().second() = Pstream::myProcNo();
+                }
+            }
+            break;
+        }
+
+        case NEARESTFACE:
+        {
+            if (samplePatch_.size() != 0 && samplePatch_ != "none")
+            {
+                FatalErrorIn
+                (
+                    "directMappedPolyPatch::findSamples(const pointField&,"
+                    " labelList&, labelList&, pointField&) const"
+                )   << "No need to supply a patch name when in "
+                    << sampleModeNames_[mode_] << " mode." << exit(FatalError);
+            }
+
+            // Octree based search engine
+            meshSearch meshSearchEngine(mesh, false);
+
+            forAll(samples, sampleI)
+            {
+                const point& sample = samples[sampleI];
+
+                label faceI = meshSearchEngine.findNearestFace(sample);
+
+                if (faceI == -1)
+                {
+                    nearest[sampleI].second().first() = Foam::sqr(GREAT);
+                    nearest[sampleI].second().second() = Pstream::myProcNo();
+                }
+                else
+                {
+                    const point& fc = mesh.faceCentres()[faceI];
+
+                    nearest[sampleI].first() = pointIndexHit
+                    (
+                        true,
+                        fc,
+                        faceI
+                    );
+                    nearest[sampleI].second().first() = magSqr(fc-sample);
+                    nearest[sampleI].second().second() = Pstream::myProcNo();
+                }
+            }
+            break;
+        }
+
+        default:
+        {
+            FatalErrorIn("directMappedPolyPatch::findSamples(..)")
+                << "problem." << abort(FatalError);
         }
     }
 
-    // Use only that processor that contains the sample
-    Pstream::listCombineGather(sampleCells, maxEqOp<label>());
-    Pstream::listCombineScatter(sampleCells);
 
-    labelList minSampleCellProcs(sampleCellProcs);
-    Pstream::listCombineGather(sampleCellProcs, maxEqOp<label>());
-    Pstream::listCombineScatter(sampleCellProcs);
-
-    Pstream::listCombineGather(sampleCc, maxEqOp<point>());
-    Pstream::listCombineScatter(sampleCc);
+    // Find nearest.
+    Pstream::listCombineGather(nearest, nearestEqOp());
+    Pstream::listCombineScatter(nearest);
 
     if (debug)
     {
         Info<< "directMappedPolyPatch::findSamples : " << endl;
-        forAll(sampleCells, sampleI)
+        forAll(nearest, sampleI)
         {
-            Info<< "    " << sampleI << " coord:" << samples[sampleI]
-                << " found on processor:" << sampleCellProcs[sampleI]
-                << " in cell:" << sampleCells[sampleI]
-                << " with cc:" << sampleCc[sampleI] << endl;
+            label procI = nearest[sampleI].second().second();
+            label localI = nearest[sampleI].first().index();
+
+            Info<< "    " << sampleI << " coord:"<< samples[sampleI]
+                << " found on processor:" << procI
+                << " in local cell/face:" << localI
+                << " with cc:" << nearest[sampleI].first().rawPoint() << endl;
         }
     }
 
     // Check for samples not being found
-    forAll(sampleCells, sampleI)
+    forAll(nearest, sampleI)
     {
-        if (sampleCells[sampleI] == -1)
+        if (!nearest[sampleI].first().hit())
         {
             FatalErrorIn
             (
-                "findSamples(const pointField&, labelList&, labelList&)"
+                "directMappedPolyPatch::findSamples"
+                "(const pointField&, labelList&"
+                ", labelList&, pointField&)"
             )   << "Did not find sample " << samples[sampleI]
                 << " on any processor." << exit(FatalError);
         }
     }
 
 
-    // Check for samples being found in multiple cells
-    {
-        forAll(minSampleCellProcs, sampleI)
-        {
-            if (minSampleCellProcs[sampleI] == -1)
-            {
-                minSampleCellProcs[sampleI] = labelMax;
-            }
-        }
-        Pstream::listCombineGather(minSampleCellProcs, minEqOp<label>());
-        Pstream::listCombineScatter(minSampleCellProcs);
+    // Convert back into proc+local index
+    sampleProcs.setSize(samples.size());
+    sampleIndices.setSize(samples.size());
+    sampleLocations.setSize(samples.size());
 
-        forAll(minSampleCellProcs, sampleI)
-        {
-            if (minSampleCellProcs[sampleI] != sampleCellProcs[sampleI])
-            {
-                FatalErrorIn
-                (
-                    "directMappedPolyPatch::findSamples"
-                    "(const pointField&, labelList&, labelList&)"
-                )   << "Found sample " << samples[sampleI]
-                    << " on processor " << minSampleCellProcs[sampleI]
-                    << " and on processor " << sampleCellProcs[sampleI] << endl
-                    << "This is illegal. {lease move your sampling plane."
-                    << exit(FatalError);
-            }
-        }
+    forAll(nearest, sampleI)
+    {
+        sampleProcs[sampleI] = nearest[sampleI].second().second();
+        sampleIndices[sampleI] = nearest[sampleI].first().index();
+        sampleLocations[sampleI] = nearest[sampleI].first().hitPoint();
     }
 }
 
+
 void Foam::directMappedPolyPatch::calcMapping() const
 {
-    if (sendCellLabelsPtr_.valid())
+    if (sendLabelsPtr_.valid())
     {
         FatalErrorIn("directMappedPolyPatch::calcMapping() const")
             << "Mapping already calculated" << exit(FatalError);
@@ -236,18 +377,18 @@ void Foam::directMappedPolyPatch::calcMapping() const
     pointField patchFc;
     collectSamples(samples, patchFaceProcs, patchFaces, patchFc);
 
-    // Find processor and cell samples are in
-    labelList sampleCellProcs;
-    labelList sampleCells;
-    pointField sampleCc;
-    findSamples(samples, sampleCellProcs, sampleCells, sampleCc);
+    // Find processor and cell/face samples are in and actual location.
+    labelList sampleProcs;
+    labelList sampleIndices;
+    pointField sampleLocations;
+    findSamples(samples, sampleProcs, sampleIndices, sampleLocations);
 
 
     // Now we have all the data we need:
     // - where sample originates from (so destination when mapping):
     //   patchFaces, patchFaceProcs.
-    // - cell sample is in (so source when mapping)
-    //   sampleCells, sampleCellProcs.
+    // - cell/face sample is in (so source when mapping)
+    //   sampleIndices, sampleProcs.
 
 
     if (debug && Pstream::master())
@@ -267,18 +408,18 @@ void Foam::directMappedPolyPatch::calcMapping() const
         {
             meshTools::writeOBJ(str, patchFc[i]);
             vertI++;
-            meshTools::writeOBJ(str, sampleCc[i]);
+            meshTools::writeOBJ(str, sampleLocations[i]);
             vertI++;
             str << "l " << vertI-1 << ' ' << vertI << nl;
         }
     }
 
 
-    // Check that actual offset vector (sampleCc - patchFc) is more or less
-    // constant.
+    // Check that actual offset vector (sampleLocations - patchFc) is more or
+    // less constant.
     if (Pstream::master())
     {
-        const scalarField magOffset(mag(sampleCc - patchFc));
+        const scalarField magOffset(mag(sampleLocations - patchFc));
         const scalar avgOffset(average(magOffset));
 
         forAll(magOffset, sampleI)
@@ -291,7 +432,7 @@ void Foam::directMappedPolyPatch::calcMapping() const
                     << "    on a single plane."
                     << " This might give numerical problems." << endl
                     << "    At patchface " << patchFc[sampleI]
-                    << " the sampled cell " << sampleCc[sampleI] << endl
+                    << " the sampled cell " << sampleLocations[sampleI] << endl
                     << "    is not on a plane " << avgOffset
                     << " offset from the patch." << endl
                     << "    You might want to shift your plane offset."
@@ -304,7 +445,7 @@ void Foam::directMappedPolyPatch::calcMapping() const
 
 
     // Determine schedule.
-    mapDistribute distMap(sampleCellProcs, patchFaceProcs);
+    mapDistribute distMap(sampleProcs, patchFaceProcs);
 
     // Rework the schedule to cell data to send, face data to receive.
     schedulePtr_.reset(new List<labelPair>(distMap.schedule()));
@@ -313,17 +454,21 @@ void Foam::directMappedPolyPatch::calcMapping() const
     const labelListList& constructMap = distMap.constructMap();
 
     // Extract the particular data I need to send and receive.
-    sendCellLabelsPtr_.reset(new labelListList(subMap.size()));
-    labelListList& sendCellLabels = sendCellLabelsPtr_();
+    sendLabelsPtr_.reset(new labelListList(subMap.size()));
+    labelListList& sendLabels = sendLabelsPtr_();
 
     forAll(subMap, procI)
     {
-        sendCellLabels[procI] = IndirectList<label>(sampleCells, subMap[procI]);
+        sendLabels[procI] = IndirectList<label>
+        (
+            sampleIndices,
+            subMap[procI]
+        );
 
         if (debug)
         {
-            Pout<< "To proc:" << procI << " sending values of cells:"
-                << sendCellLabels[procI] << endl;
+            Pout<< "To proc:" << procI << " sending values of cells/faces:"
+                << sendLabels[procI] << endl;
         }
     }
 
@@ -356,9 +501,11 @@ Foam::directMappedPolyPatch::directMappedPolyPatch
 )
 :
     polyPatch(name, size, start, index, bm),
+    mode_(NEARESTPATCHFACE),
+    samplePatch_("none"),
     offset_(vector::zero),
     schedulePtr_(NULL),
-    sendCellLabelsPtr_(NULL),
+    sendLabelsPtr_(NULL),
     receiveFaceLabelsPtr_(NULL)
 {}
 
@@ -372,9 +519,11 @@ Foam::directMappedPolyPatch::directMappedPolyPatch
 )
 :
     polyPatch(name, dict, index, bm),
+    mode_(sampleModeNames_.read(dict.lookup("sampleMode"))),
+    samplePatch_(dict.lookup("samplePatch")),
     offset_(dict.lookup("offset")),
     schedulePtr_(NULL),
-    sendCellLabelsPtr_(NULL),
+    sendLabelsPtr_(NULL),
     receiveFaceLabelsPtr_(NULL)
 {}
 
@@ -386,9 +535,11 @@ Foam::directMappedPolyPatch::directMappedPolyPatch
 )
 :
     polyPatch(pp, bm),
+    mode_(pp.mode_),
+    samplePatch_(pp.samplePatch_),
     offset_(pp.offset_),
     schedulePtr_(NULL),
-    sendCellLabelsPtr_(NULL),
+    sendLabelsPtr_(NULL),
     receiveFaceLabelsPtr_(NULL)
 {}
 
@@ -403,9 +554,11 @@ Foam::directMappedPolyPatch::directMappedPolyPatch
 )
 :
     polyPatch(pp, bm, index, newSize, newStart),
+    mode_(pp.mode_),
+    samplePatch_(pp.samplePatch_),
     offset_(pp.offset_),
     schedulePtr_(NULL),
-    sendCellLabelsPtr_(NULL),
+    sendLabelsPtr_(NULL),
     receiveFaceLabelsPtr_(NULL)
 {}
 
@@ -423,7 +576,7 @@ Foam::directMappedPolyPatch::~directMappedPolyPatch()
 void Foam::directMappedPolyPatch::clearOut()
 {
     schedulePtr_.clear();
-    sendCellLabelsPtr_.clear();
+    sendLabelsPtr_.clear();
     receiveFaceLabelsPtr_.clear();
 }
 
@@ -431,6 +584,10 @@ void Foam::directMappedPolyPatch::clearOut()
 void Foam::directMappedPolyPatch::write(Ostream& os) const
 {
     polyPatch::write(os);
+    os.writeKeyword("sampleMode") << sampleModeNames_[mode_]
+        << token::END_STATEMENT << nl;
+    os.writeKeyword("samplePatch") << samplePatch_
+        << token::END_STATEMENT << nl;
     os.writeKeyword("offset") << offset_ << token::END_STATEMENT << nl;
 }
 

src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.H

@@ -26,8 +26,8 @@ Class
     Foam::directMappedPolyPatch
 
 Description
-    Determines a mapping between patch face centres and mesh cell centres and
-    processors they're on.
+    Determines a mapping between patch face centres and mesh cell or face
+    centres and processors they're on.
 
 Note
     Storage is not optimal. It stores all face centres and cells on all
@@ -43,6 +43,9 @@ SourceFiles
 
 #include "polyPatch.H"
 #include "labelPair.H"
+#include "Tuple2.H"
+#include "pointIndexHit.H"
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -57,19 +60,41 @@ class directMappedPolyPatch
 :
     public polyPatch
 {
+
+public:
+
+        //- Mesh items to sample
+        enum sampleMode
+        {
+            NEARESTCELL,
+            NEARESTPATCHFACE,
+            NEARESTFACE
+        };
+
+private:
+
     // Private data
 
+        static const NamedEnum<sampleMode, 3> sampleModeNames_;
+
+        //- What to sample
+        sampleMode mode_;
+
+        //- Patch (only if NEARESTBOUNDARY)
+        const word samplePatch_;
+
         //- Offset vector
         const vector offset_;
 
 
+
         // Derived information
 
             //- Communication schedule.
             mutable autoPtr<List<labelPair> > schedulePtr_;
 
-            //- Cells to sample per processor
-            mutable autoPtr<labelListList> sendCellLabelsPtr_;
+            //- Cells/faces to sample per processor
+            mutable autoPtr<labelListList> sendLabelsPtr_;
 
             //- Patch faces to receive per processor
             mutable autoPtr<labelListList> receiveFaceLabelsPtr_;
@@ -86,19 +111,50 @@ class directMappedPolyPatch
             pointField& patchFc
         ) const;
 
-        //- Find cells containing samples
+        //- Find cells/faces containing samples
         void findSamples
         (
             const pointField&,
-            labelList& sampleCellProcs,
-            labelList& sampleCells,
-            pointField& sampleCc
+            labelList& sampleProcs,     // processor containing sample
+            labelList& sampleIndices,   // local index of cell/face
+            pointField& sampleLocations // actual representative location
         ) const;
 
         //- Calculate matching
         void calcMapping() const;
 
 
+    // Private class
+
+        //- 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;
+                    }
+                }
+            }
+        };
+
+
+
 protected:
 
         void clearOut();
@@ -231,21 +287,25 @@ public:
             return schedulePtr_();
         }
 
-        //- Cells to sample per processor
-        const labelListList& sendCellLabels() const
+        //- Cells/faces to sample per processor
+        const labelListList& sendLabels() const
         {
-            if (!sendCellLabelsPtr_.valid())
+            Pout<< "Asking for sendLabels." << endl;
+            if (!sendLabelsPtr_.valid())
             {
+                Pout<< "Calculating mapping." << endl;
                 calcMapping();
             }
-            return sendCellLabelsPtr_();
+            return sendLabelsPtr_();
         }
 
         //- Patch faces to receive per processor
         const labelListList& receiveFaceLabels() const
         {
+            Pout<< "Asking for receiveFaceLabels." << endl;
             if (!receiveFaceLabelsPtr_.valid())
             {
+                Pout<< "Calculating mapping." << endl;
                 calcMapping();
             }
             return receiveFaceLabelsPtr_();

src/meshTools/indexedOctree/indexedOctree.C

@@ -26,7 +26,7 @@ License
 
 #include "indexedOctree.H"
 #include "linePointRef.H"
-#include "triSurface.H"
+//#include "triSurface.H"
 #include "meshTools.H"
 #include "OFstream.H"
 

src/meshTools/indexedOctree/treeDataCell.C

@@ -26,7 +26,7 @@ License
 
 #include "treeDataCell.H"
 #include "indexedOctree.H"
-#include "polyMesh.H"
+#include "primitiveMesh.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -71,7 +71,7 @@ Foam::treeBoundBox Foam::treeDataCell::calcCellBb(const label cellI) const
 Foam::treeDataCell::treeDataCell
 (
     const bool cacheBb,
-    const polyMesh& mesh,
+    const primitiveMesh& mesh,
     const labelList& cellLabels
 )
 :
@@ -94,7 +94,7 @@ Foam::treeDataCell::treeDataCell
 Foam::treeDataCell::treeDataCell
 (
     const bool cacheBb,
-    const polyMesh& mesh
+    const primitiveMesh& mesh
 )
 :
     mesh_(mesh),

src/meshTools/indexedOctree/treeDataCell.H

@@ -47,7 +47,7 @@ namespace Foam
 {
 
 // Forward declaration of classes
-class polyMesh;
+class primitiveMesh;
 template<class Type> class indexedOctree;
 
 /*---------------------------------------------------------------------------*\
@@ -58,7 +58,7 @@ class treeDataCell
 {
     // Private data
 
-        const polyMesh& mesh_;
+        const primitiveMesh& mesh_;
 
         //- Subset of cells to work on
         const labelList cellLabels_;
@@ -87,12 +87,12 @@ public:
         treeDataCell
         (
             const bool cacheBb,
-            const polyMesh&,
+            const primitiveMesh&,
             const labelList&
         );
 
         //- Construct from mesh. Uses all cells in mesh.
-        treeDataCell(const bool cacheBb, const polyMesh&);
+        treeDataCell(const bool cacheBb, const primitiveMesh&);
 
 
     // Member Functions
@@ -104,7 +104,7 @@ public:
                 return cellLabels_;
             }
 
-            const polyMesh& mesh() const
+            const primitiveMesh& mesh() const
             {
                 return mesh_;
             }

src/meshTools/meshSearch/meshSearch.C

@@ -22,14 +22,11 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Description
-
 \*---------------------------------------------------------------------------*/
 
-#include "polyMesh.H"
 #include "meshSearch.H"
-#include "triPointRef.H"
-#include "octree.H"
+#include "polyMesh.H"
+#include "indexedOctree.H"
 #include "pointIndexHit.H"
 #include "DynamicList.H"
 #include "demandDrivenData.H"
@@ -43,14 +40,75 @@ Foam::scalar Foam::meshSearch::tol_ = 1E-3;
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+bool Foam::meshSearch::findNearer
+(
+    const point& sample,
+    const pointField& points,
+    label& nearestI,
+    scalar& nearestDistSqr
+)
+{
+    bool nearer = false;
+
+    forAll(points, pointI)
+    {
+        scalar distSqr = magSqr(points[pointI] - sample);
+
+        if (distSqr < nearestDistSqr)
+        {
+            nearestDistSqr = distSqr;
+            nearestI = pointI;
+            nearer = true;
+        }
+    }
+
+    return nearer;
+}
+
+
+bool Foam::meshSearch::findNearer
+(
+    const point& sample,
+    const pointField& points,
+    const labelList& indices,
+    label& nearestI,
+    scalar& nearestDistSqr
+)
+{
+    bool nearer = false;
+
+    forAll(indices, i)
+    {
+        label pointI = indices[i];
+
+        scalar distSqr = magSqr(points[pointI] - sample);
+
+        if (distSqr < nearestDistSqr)
+        {
+            nearestDistSqr = distSqr;
+            nearestI = pointI;
+            nearer = true;
+        }
+    }
+
+    return nearer;
+}
+
+
 // tree based searching
 Foam::label Foam::meshSearch::findNearestCellTree(const point& location) const
 {
-    treeBoundBox tightest(treeBoundBox::greatBox);
+    const indexedOctree<treeDataPoint>& tree = cellCentreTree();
 
-    scalar tightestDist = GREAT;
+    scalar span = mag(tree.bb().max() - tree.bb().min());
 
-    return cellCentreTree().findNearest(location, tightest, tightestDist);
+    pointIndexHit info = tree.findNearest(location, Foam::sqr(span));
+
+    if (!info.hit())
+    {
+        info = tree.findNearest(location, Foam::sqr(GREAT));
+    }
+    return info.index();
 }
 
 
@@ -59,21 +117,18 @@ Foam::label Foam::meshSearch::findNearestCellLinear(const point& location) const
 {
     const vectorField& centres = mesh_.cellCentres();
 
-    label nearestCelli = 0;
-    scalar minProximity = magSqr(centres[0] - location);
+    label nearestIndex = 0;
+    scalar minProximity = magSqr(centres[nearestIndex] - location);
 
-    forAll(centres, celli)
-    {
-        scalar proximity = magSqr(centres[celli] - location);
+    findNearer
+    (
+        location,
+        centres,
+        nearestIndex,
+        minProximity
+    );
 
-        if (proximity < minProximity)
-        {
-            nearestCelli = celli;
-            minProximity = proximity;
-        }
-    }
-
-    return nearestCelli;
+    return nearestIndex;
 }
 
 
@@ -92,40 +147,145 @@ Foam::label Foam::meshSearch::findNearestCellWalk
         )   << "illegal seedCell:" << seedCellI << exit(FatalError);
     }
 
-    const vectorField& centres = mesh_.cellCentres();
-    const labelListList& cc = mesh_.cellCells();
-
-
     // Walk in direction of face that decreases distance
 
-    label curCell = seedCellI;
-    scalar distanceSqr = magSqr(centres[curCell] - location);
+    label curCellI = seedCellI;
+    scalar distanceSqr = magSqr(mesh_.cellCentres()[curCellI] - location);
 
     bool closer;
 
     do
     {
-        closer = false;
-
-        // set the current list of neighbouring cells
-        const labelList& neighbours = cc[curCell];
-
-        forAll (neighbours, nI)
-        {
-            scalar curDistSqr = magSqr(centres[neighbours[nI]] - location);
-
-            // search through all the neighbours.
-            // If the cell is closer, reset current cell and distance
-            if (curDistSqr < distanceSqr)
-            {
-                distanceSqr = curDistSqr;
-                curCell = neighbours[nI];
-                closer = true;    // a closer neighbour has been found
-            }
-        }
+        // Try neighbours of curCellI
+        closer = findNearer
+        (
+            location,
+            mesh_.cellCentres(),
+            mesh_.cellCells()[curCellI],
+            curCellI,
+            distanceSqr
+        );
     } while (closer);
 
-    return curCell;
+    return curCellI;
+}
+
+
+// tree based searching
+Foam::label Foam::meshSearch::findNearestFaceTree(const point& location) const
+{
+    // Search nearest cell centre.
+    const indexedOctree<treeDataPoint>& tree = cellCentreTree();
+
+    scalar span = mag(tree.bb().max() - tree.bb().min());
+
+    // Search with decent span
+    pointIndexHit info = tree.findNearest(location, Foam::sqr(span));
+
+    if (!info.hit())
+    {
+        // Search with desparate span
+        info = tree.findNearest(location, Foam::sqr(GREAT));
+    }
+
+
+    // Now check any of the faces of the nearest cell
+    const vectorField& centres = mesh_.faceCentres();
+    const cell& ownFaces = mesh_.cells()[info.index()];
+
+    label nearestFaceI = ownFaces[0];
+    scalar minProximity = magSqr(centres[nearestFaceI] - location);
+
+    findNearer
+    (
+        location,
+        centres,
+        ownFaces,
+        nearestFaceI,
+        minProximity
+    );
+
+    return nearestFaceI;
+}
+
+
+// linear searching
+Foam::label Foam::meshSearch::findNearestFaceLinear(const point& location) const
+{
+    const vectorField& centres = mesh_.faceCentres();
+
+    label nearestFaceI = 0;
+    scalar minProximity = magSqr(centres[nearestFaceI] - location);
+
+    findNearer
+    (
+        location,
+        centres,
+        nearestFaceI,
+        minProximity
+    );
+
+    return nearestFaceI;
+}
+
+
+// walking from seed
+Foam::label Foam::meshSearch::findNearestFaceWalk
+(
+    const point& location,
+    const label seedFaceI
+) const
+{
+    if (seedFaceI < 0)
+    {
+        FatalErrorIn
+        (
+            "meshSearch::findNearestFaceWalk(const point&, const label)"
+        )   << "illegal seedFace:" << seedFaceI << exit(FatalError);
+    }
+
+    const vectorField& centres = mesh_.faceCentres();
+
+
+    // Walk in direction of face that decreases distance
+
+    label curFaceI = seedFaceI;
+    scalar distanceSqr = magSqr(centres[curFaceI] - location);
+
+    while (true)
+    {
+        label betterFaceI = curFaceI;
+
+        findNearer
+        (
+            location,
+            centres,
+            mesh_.cells()[mesh_.faceOwner()[curFaceI]],
+            betterFaceI,
+            distanceSqr
+        );
+
+        if (mesh_.isInternalFace(curFaceI))
+        {
+            findNearer
+            (
+                location,
+                centres,
+                mesh_.cells()[mesh_.faceNeighbour()[curFaceI]],
+                betterFaceI,
+                distanceSqr
+            );
+        }
+
+        if (betterFaceI == curFaceI)
+        {
+            break;
+        }
+
+        curFaceI = betterFaceI;
+    }
+
+    return curFaceI;
 }
 
 
@@ -180,12 +340,11 @@ Foam::label Foam::meshSearch::findNearestBoundaryFaceWalk
 
     const face& f =  mesh_.faces()[curFaceI];
 
-    scalar minDist =
-        f.nearestPoint
-        (
-            location,
-            mesh_.points()
-        ).distance();
+    scalar minDist = f.nearestPoint
+    (
+        location,
+        mesh_.points()
+    ).distance();
 
     bool closer;
 
@@ -202,8 +361,7 @@ Foam::label Foam::meshSearch::findNearestBoundaryFaceWalk
 
         forAll (myEdges, myEdgeI)
         {
-            const labelList& neighbours =
-                mesh_.edgeFaces()[myEdges[myEdgeI]];
+            const labelList& neighbours = mesh_.edgeFaces()[myEdges[myEdgeI]];
 
             // Check any face which uses edge, is boundary face and
             // is not curFaceI itself.
@@ -220,12 +378,11 @@ Foam::label Foam::meshSearch::findNearestBoundaryFaceWalk
                 {
                     const face& f =  mesh_.faces()[faceI];
 
-                    pointHit curHit =
-                        f.nearestPoint
-                        (
-                            location,
-                            mesh_.points()
-                        );
+                    pointHit curHit = f.nearestPoint
+                    (
+                        location,
+                        mesh_.points()
+                    );
 
                     // If the face is closer, reset current face and distance
                     if (curHit.distance() < minDist)
@@ -283,9 +440,11 @@ Foam::meshSearch::~meshSearch()
     clearOut();
 }
 
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-const Foam::octree<Foam::octreeDataFace>& Foam::meshSearch::boundaryTree() const
+const Foam::indexedOctree<Foam::treeDataFace>& Foam::meshSearch::boundaryTree()
+ const
 {
     if (!boundaryTreePtr_)
     {
@@ -294,17 +453,28 @@ const Foam::octree<Foam::octreeDataFace>& Foam::meshSearch::boundaryTree() const
         //
 
         // all boundary faces (not just walls)
-        octreeDataFace shapes(mesh_);
+        labelList bndFaces(mesh_.nFaces()-mesh_.nInternalFaces());
+        forAll(bndFaces, i)
+        {
+            bndFaces[i] = mesh_.nInternalFaces() + i;
+        }
 
         treeBoundBox overallBb(mesh_.points());
 
-        boundaryTreePtr_ = new octree<octreeDataFace>
+        Random rndGen(123456);
+
+        boundaryTreePtr_ = new indexedOctree<treeDataFace>
         (
-            overallBb,  // overall search domain
-            shapes,     // all information needed to do checks on cells
-            1,          // min levels
-            20.0,       // maximum ratio of cubes v.s. cells
-            10.0
+            treeDataFace    // all information needed to search faces
+            (
+                false,                      // do not cache bb
+                mesh_,
+                bndFaces                    // boundary faces only
+            ),
+            overallBb.extend(rndGen, 1E-3), // overall search domain
+            8,                              // maxLevel
+            10,                             // leafsize
+            3.0                             // duplicity
         );
     }
 
@@ -312,7 +482,8 @@ const Foam::octree<Foam::octreeDataFace>& Foam::meshSearch::boundaryTree() const
 }
 
 
-const Foam::octree<Foam::octreeDataCell>& Foam::meshSearch::cellTree() const
+const Foam::indexedOctree<Foam::treeDataCell>& Foam::meshSearch::cellTree()
+ const
 {
     if (!cellTreePtr_)
     {
@@ -320,17 +491,19 @@ const Foam::octree<Foam::octreeDataCell>& Foam::meshSearch::cellTree() const
         // Construct tree
         //
 
-        octreeDataCell shapes(mesh_);
-
         treeBoundBox overallBb(mesh_.points());
 
-        cellTreePtr_ = new octree<octreeDataCell>
+        cellTreePtr_ = new indexedOctree<treeDataCell>
         (
+            treeDataCell
+            (
+                false,  // not cache bb
+                mesh_
+            ),
             overallBb,  // overall search domain
-            shapes,     // all information needed to do checks on cells
-            1,          // min levels
-            20.0,       // maximum ratio of cubes v.s. cells
-            2.0
+            8,      // maxLevel
+            10,     // leafsize
+            3.0     // duplicity
         );
     }
 
@@ -339,8 +512,8 @@ const Foam::octree<Foam::octreeDataCell>& Foam::meshSearch::cellTree() const
 }
 
 
-const Foam::octree<Foam::octreeDataPoint>& Foam::meshSearch::cellCentreTree()
-    const
+const Foam::indexedOctree<Foam::treeDataPoint>&
+ Foam::meshSearch::cellCentreTree() const
 {
     if (!cellCentreTreePtr_)
     {
@@ -348,17 +521,14 @@ const Foam::octree<Foam::octreeDataPoint>& Foam::meshSearch::cellCentreTree()
         // Construct tree
         //
 
-        // shapes holds reference to cellCentres!
-        octreeDataPoint shapes(mesh_.cellCentres());
-
         treeBoundBox overallBb(mesh_.cellCentres());
 
-        cellCentreTreePtr_ = new octree<octreeDataPoint>
+        cellCentreTreePtr_ = new indexedOctree<treeDataPoint>
         (
+            treeDataPoint(mesh_.cellCentres()),
             overallBb,  // overall search domain
-            shapes,     // all information needed to do checks on cells
-            1,          // min levels
-            20.0,       // maximum ratio of cubes v.s. cells
+            10,         // max levels
+            10.0,       // maximum ratio of cubes v.s. cells
             100.0       // max. duplicity; n/a since no bounding boxes.
         );
     }
@@ -396,15 +566,14 @@ bool Foam::meshSearch::pointInCell(const point& p, label cellI) const
 
             forAll(f, fp)
             {
-                pointHit inter = 
-                    f.ray
-                    (
-                        ctr,
-                        dir,
-                        mesh_.points(),
-                        intersection::HALF_RAY,
-                        intersection::VECTOR
-                    );
+                pointHit inter = f.ray
+                (
+                    ctr,
+                    dir,
+                    mesh_.points(),
+                    intersection::HALF_RAY,
+                    intersection::VECTOR
+                );
 
                 if (inter.hit())
                 {
@@ -484,6 +653,31 @@ Foam::label Foam::meshSearch::findNearestCell
 }
 
 
+Foam::label Foam::meshSearch::findNearestFace
+(
+    const point& location,
+    const label seedFaceI,
+    const bool useTreeSearch
+) const
+{
+    if (seedFaceI == -1)
+    {
+        if (useTreeSearch)
+        {
+            return findNearestFaceTree(location);
+        }
+        else
+        {
+            return findNearestFaceLinear(location);
+        }
+    }
+    else
+    {
+        return findNearestFaceWalk(location, seedFaceI);
+    }
+}
+
+
 Foam::label Foam::meshSearch::findCell
 (
     const point& location,
@@ -607,19 +801,26 @@ Foam::label Foam::meshSearch::findNearestBoundaryFace
     {
         if (useTreeSearch)
         {
-            treeBoundBox tightest(treeBoundBox::greatBox);
+            const indexedOctree<treeDataFace>& tree =  boundaryTree();
 
-            scalar tightestDist = GREAT;
+            scalar span = mag(tree.bb().max() - tree.bb().min());
 
-            label index =
-                boundaryTree().findNearest
+            pointIndexHit info = boundaryTree().findNearest
+            (
+                location,
+                Foam::sqr(span)
+            );
+
+            if (!info.hit())
+            {
+                info = boundaryTree().findNearest
                 (
                     location,
-                    tightest,
-                    tightestDist
+                    Foam::sqr(GREAT)
                 );
+            }
 
-            return boundaryTree().shapes().meshFaces()[index];
+            return tree.shapes().faceLabels()[info.index()];
         }
         else
         {
@@ -670,7 +871,7 @@ Foam::pointIndexHit Foam::meshSearch::intersection
     if (curHit.hit())
     {
         // Change index into octreeData into face label
-        curHit.setIndex(boundaryTree().shapes().meshFaces()[curHit.index()]);
+        curHit.setIndex(boundaryTree().shapes().faceLabels()[curHit.index()]);
     }
     return curHit;
 }
@@ -733,7 +934,9 @@ Foam::List<Foam::pointIndexHit> Foam::meshSearch::intersections
 
 bool Foam::meshSearch::isInside(const point& p) const
 {
-    return boundaryTree().getSampleType(p) == octree<octreeDataFace>::INSIDE;
+    return
+        boundaryTree().getVolumeType(p)
+     == indexedOctree<treeDataFace>::INSIDE;
 }
 
 

src/meshTools/meshSearch/meshSearch.H

@@ -26,7 +26,8 @@ Class
     Foam::meshSearch
 
 Description
-    Various searches on polyMesh; uses (demand driven) octree to search.
+    Various (local, not parallel) searches on polyMesh;
+    uses (demand driven) octree to search.
 
 SourceFiles
     meshSearch.C
@@ -36,11 +37,10 @@ SourceFiles
 #ifndef meshSearch_H
 #define meshSearch_H
 
-#include "octreeDataCell.H"
-#include "octreeDataFace.H"
-#include "octreeDataPoint.H"
+#include "treeDataCell.H"
+#include "treeDataFace.H"
+#include "treeDataPoint.H"
 #include "pointIndexHit.H"
-#include "className.H"
 #include "Cloud.H"
 #include "passiveParticle.H"
 
@@ -71,46 +71,77 @@ class meshSearch
 
         //- demand driven octrees
 
-        mutable octree<octreeDataFace>* boundaryTreePtr_;
-        mutable octree<octreeDataCell>* cellTreePtr_;
-        mutable octree<octreeDataPoint>* cellCentreTreePtr_;
+        mutable indexedOctree<treeDataFace>* boundaryTreePtr_;
+        mutable indexedOctree<treeDataCell>* cellTreePtr_;
+        mutable indexedOctree<treeDataPoint>* cellCentreTreePtr_;
 
 
     // Private Member Functions
 
-        //- nearest cell centre using octree
-        label findNearestCellTree(const point& location) const;
-
-        //- nearest cell centre going through all cells
-        label findNearestCellLinear(const point& location) const;
-
-        //- walk from seed. Does not 'go around' boundary, just returns
-        //  last cell before boundary.
-        label findNearestCellWalk
+        //- Updates nearestI, nearestDistSqr from any closer ones.
+        static bool findNearer
         (
-            const point& location,
-            const label seedCellI
-        ) const;
+            const point& sample,
+            const pointField& points,
+            label& nearestI,
+            scalar& nearestDistSqr
+        );
 
-        //- cell containing location. Linear search.
-        label findCellLinear(const point& location) const;
-
-        //- walk from seed to find nearest boundary face. Gets stuck in
-        //  local minimum.
-        label findNearestBoundaryFaceWalk
+        //- Updates nearestI, nearestDistSqr from any selected closer ones.
+        static bool findNearer
         (
-            const point& location,
-            const label seedFaceI
-        ) const;
+            const point& sample,
+            const pointField& points,
+            const labelList& indices,
+            label& nearestI,
+            scalar& nearestDistSqr
+        );
 
-        //- Calculate offset vector in direction dir with as length a fraction
-        //  of the cell size (of the cell straddling boundary face)
-        vector offset
-        (
-            const point& bPoint,
-            const label bFaceI,
-            const vector& dir
-        ) const;
+
+        // Cells
+
+            //- nearest cell centre using octree
+            label findNearestCellTree(const point&) const;
+
+            //- nearest cell centre going through all cells
+            label findNearestCellLinear(const point&) const;
+
+            //- walk from seed. Does not 'go around' boundary, just returns
+            //  last cell before boundary.
+            label findNearestCellWalk(const point&, const label) const;
+
+            //- cell containing location. Linear search.
+            label findCellLinear(const point&) const;
+
+
+        // Cells
+
+            label findNearestFaceTree(const point&) const;
+
+            label findNearestFaceLinear(const point&) const;
+
+            label findNearestFaceWalk(const point&, const label) const;
+
+
+
+        // Boundary faces
+
+            //- walk from seed to find nearest boundary face. Gets stuck in
+            //  local minimum.
+            label findNearestBoundaryFaceWalk
+            (
+                const point& location,
+                const label seedFaceI
+            ) const;
+
+            //- Calculate offset vector in direction dir with as length a
+            //  fraction of the cell size (of the cell straddling boundary face)
+            vector offset
+            (
+                const point& bPoint,
+                const label bFaceI,
+                const vector& dir
+            ) const;
 
 
         //- Disallow default bitwise copy construct
@@ -154,13 +185,13 @@ public:
 
             //- Get (demand driven) reference to octree holding all
             //  boundary faces
-            const octree<octreeDataFace>& boundaryTree() const;
+            const indexedOctree<treeDataFace>& boundaryTree() const;
 
             //- Get (demand driven) reference to octree holding all cells
-            const octree<octreeDataCell>& cellTree() const;
+            const indexedOctree<treeDataCell>& cellTree() const;
 
             //- Get (demand driven) reference to octree holding all cell centres
-            const octree<octreeDataPoint>& cellCentreTree() const;
+            const indexedOctree<treeDataPoint>& cellCentreTree() const;
 
 
         // Queries
@@ -181,6 +212,13 @@ public:
                 const bool useTreeSearch = true
             ) const;
 
+            label findNearestFace
+            (
+                const point& location,
+                const label seedFaceI = -1,
+                const bool useTreeSearch = true
+            ) const;
+
             //- Find cell containing (using pointInCell) location.
             //  If seed provided walks and falls back to linear/tree search.
             //  (so handles holes correctly)s
@@ -206,7 +244,7 @@ public:
             //- Find first intersection of boundary in segment [pStart, pEnd]
             //  (so inclusive of endpoints). Always octree for now
             pointIndexHit intersection(const point& pStart, const point& pEnd)
-                const;
+            const;
 
             //- Find all intersections of boundary within segment pStart .. pEnd
             //  Always octree for now

src/meshTools/meshTools/meshTools.C

@@ -182,7 +182,7 @@ Foam::vectorField Foam::meshTools::calcBoxPointNormals(const primitivePatch& pp)
                     "Foam::meshTools::calcBoxPointNormals"
                     "(const primitivePatch& pp)"
                 )   << "No visible octant for point:" << pp.meshPoints()[pointI]
-                    << " cooord:" << pp.localPoints()[pointI] << nl
+                    << " cooord:" << pp.points()[pp.meshPoints()[pointI]] << nl
                     << "Normal set to " << pn[pointI] << endl;
             }
         }
@@ -299,7 +299,7 @@ bool Foam::meshTools::edgeOnCell
     const label edgeI
 )
 {
-    return findIndex(mesh.edgeCells()[edgeI], cellI) != -1;
+    return findIndex(mesh.edgeCells(edgeI), cellI) != -1;
 }        
 
 
@@ -310,7 +310,7 @@ bool Foam::meshTools::edgeOnFace
     const label edgeI
 )
 {
-    return findIndex(mesh.faceEdges()[faceI], edgeI) != -1;
+    return findIndex(mesh.faceEdges(faceI), edgeI) != -1;
 }        
 
 
@@ -403,7 +403,6 @@ Foam::label Foam::meshTools::getSharedEdge
     const labelList& f0Edges = mesh.faceEdges()[f0];
     const labelList& f1Edges = mesh.faceEdges()[f1];
 
-
     forAll(f0Edges, f0EdgeI)
     {
         label edge0 = f0Edges[f0EdgeI];
@@ -481,7 +480,7 @@ void Foam::meshTools::getEdgeFaces
     label& face1
 )
 {
-    const labelList& eFaces = mesh.edgeFaces()[edgeI];
+    const labelList& eFaces = mesh.edgeFaces(edgeI);
 
     face0 = -1;
     face1 = -1;
@@ -619,7 +618,7 @@ Foam::label Foam::meshTools::walkFace
     const label nEdges
 )
 {
-    const labelList& fEdges = mesh.faceEdges()[faceI];
+    const labelList& fEdges = mesh.faceEdges(faceI);
 
     label edgeI = startEdgeI;
 
@@ -790,13 +789,4 @@ Foam::label Foam::meshTools::cutDirToEdge
 }
 
 
-// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
-
-
 // ************************************************************************* //

src/meshTools/octree/treeBoundBox.H

@@ -314,6 +314,7 @@ public:
 
             //- Return slightly wider bounding box
             //  Extends all dimensions with s*span*Random::scalar01()
+            //  and guarantees in any direction s*mag(span) minimum width
             inline treeBoundBox extend(Random&, const scalar s) const;
 
     // Friend Operators

src/meshTools/octree/treeBoundBoxI.H

@@ -437,7 +437,15 @@ inline treeBoundBox treeBoundBox::extend(Random& rndGen, const scalar s) const
 {
     treeBoundBox bb(*this);
 
-    const vector span(bb.max() - bb.min());
+    vector span(bb.max() - bb.min());
+
+    // Make 3D
+    scalar magSpan = Foam::mag(span);
+
+    for (direction dir = 0; dir < vector::nComponents; dir++)
+    {
+        span[dir] = Foam::max(s*magSpan, span[dir]);
+    }
 
     bb.min() -= cmptMultiply(s*rndGen.vector01(), span);
     bb.max() += cmptMultiply(s*rndGen.vector01(), span);

src/meshTools/searchableSurface/triSurfaceMesh.C

@@ -230,7 +230,7 @@ const Foam::indexedOctree<Foam::treeDataTriSurface>&
             new indexedOctree<treeDataTriSurface>
             (
                 treeDataTriSurface(*this),
-                bb.extend(rndGen, 1E-3),    // slightly randomize bb
+                bb.extend(rndGen, 1E-4),    // slightly randomize bb
                 10,     // maxLevel
                 10,     // leafsize
                 3.0     // duplicity
@@ -274,7 +274,7 @@ const Foam::indexedOctree<Foam::treeDataEdge>&
                     localPoints(),  // points
                     bEdges          // selected edges
                 ),
-                bb.extend(rndGen, 1E-3),    // slightly randomize bb
+                bb.extend(rndGen, 1E-4),    // slightly randomize bb
                 8,      // maxLevel
                 10,     // leafsize
                 3.0     // duplicity

src/meshTools/sets/cellSources/surfaceToCell/surfaceToCell.C

@@ -31,8 +31,6 @@ Description
 #include "meshSearch.H"
 #include "triSurface.H"
 #include "triSurfaceSearch.H"
-#include "octree.H"
-#include "octreeDataTriSurface.H"
 #include "cellClassification.H"
 
 #include "addToRunTimeSelectionTable.H"

src/meshTools/sets/pointSources/surfaceToPoint/surfaceToPoint.C

@@ -30,8 +30,6 @@ Description
 #include "polyMesh.H"
 #include "meshSearch.H"
 #include "triSurfaceSearch.H"
-#include "octree.H"
-#include "octreeDataTriSurface.H"
 
 #include "addToRunTimeSelectionTable.H"
 

src/turbulenceModels/RAS/incompressible/kOmega/kOmega.C

@@ -44,7 +44,6 @@ addToRunTimeSelectionTable(RASModel, kOmega, dictionary);
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 kOmega::kOmega
 (
     const volVectorField& U,
@@ -72,6 +71,15 @@ kOmega::kOmega
             0.072
         )
     ),
+    alpha_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "alpha",
+            coeffDict_,
+            0.52
+        )
+    ),
     alphaK_
     (
         dimensioned<scalar>::lookupOrAddToDict
@@ -222,7 +230,7 @@ void kOmega::correct()
       - fvm::Sp(fvc::div(phi_), omega_)
       - fvm::laplacian(DomegaEff(), omega_)
      ==
-        G*omega_/k_
+        alpha_*G*omega_/k_
       - fvm::Sp(beta_*omega_, omega_)
     );
 

src/turbulenceModels/RAS/incompressible/kOmega/kOmega.H

@@ -35,6 +35,9 @@ Description
         D. C. Wilcox,
         DCW Industries, Inc., La Canada,
         California, 1998.
+
+        See also:
+        http://www.cfd-online.com/Wiki/Wilcox's_k-omega_model
     @endverbatim
 
     The default model coefficients correspond to the following:
@@ -42,6 +45,7 @@ Description
         kOmegaCoeffs
         {
             Cmu         0.09;  // Equivalent to betaStar
+            alpha       0.52;
             beta        0.072;
             alphak      0.5;
             alphaOmega  0.5;
@@ -81,6 +85,7 @@ class kOmega
 
             dimensionedScalar Cmu_;
             dimensionedScalar beta_;
+            dimensionedScalar alpha_;
             dimensionedScalar alphaK_;
             dimensionedScalar alphaOmega_;