Merge branch 'master' of ssh://hunt/home/hunt2/OpenFOAM/OpenFOAM-dev

applications/solvers/multiphase/lesInterFoam/lesInterFoam.C

@@ -85,6 +85,8 @@ int main(int argc, char *argv[])
 
         #include "continuityErrs.H"
 
+        p = pd + rho*gh;
+
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

applications/utilities/mesh/advanced/refineWallLayer/refineWallLayer.C

@@ -184,10 +184,6 @@ int main(int argc, char *argv[])
         << endl;
 
     // Transfer DynamicLists to straight ones.
-    labelList cutEdges;
-    cutEdges.transfer(allCutEdges);
-    allCutEdges.clear();
-
     scalarField cutEdgeWeights;
     cutEdgeWeights.transfer(allCutEdgeWeights);
     allCutEdgeWeights.clear();
@@ -199,7 +195,7 @@ int main(int argc, char *argv[])
         mesh,
         cutCells.toc(),     // cells candidate for cutting
         labelList(0),       // cut vertices
-        cutEdges,           // cut edges
+        allCutEdges,        // cut edges
         cutEdgeWeights      // weight on cut edges
     );
 

applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C

@@ -433,7 +433,7 @@ int main(int argc, char *argv[])
     Info<< "Finished meshing in = "
         << runTime.elapsedCpuTime() << " s." << endl;
 
-    Pout<< "End\n" << endl;
+    Info<< "End\n" << endl;
 
     return(0);
 }

applications/utilities/parallelProcessing/decomposePar/decomposeMesh.C

@@ -35,11 +35,8 @@ Description
 #include "IOstreams.H"
 #include "SLPtrList.H"
 #include "boolList.H"
-#include "cellList.H"
 #include "primitiveMesh.H"
-#include "cyclicFvPatch.H"
+#include "cyclicPolyPatch.H"
-#include "fvPatchList.H"
-#include "DynamicList.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 

applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C

@@ -33,7 +33,7 @@ License
 #include "OSspecific.H"
 #include "Map.H"
 #include "globalMeshData.H"
-
+#include "DynamicList.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 

applications/utilities/postProcessing/dataConversion/foamToFieldview9/foamToFieldview9.C

@@ -45,7 +45,6 @@ Description
 #include "IOobjectList.H"
 #include "boolList.H"
 #include "stringList.H"
-#include "DynamicList.H"
 #include "cellModeller.H"
 
 #include "floatScalar.H"

applications/utilities/postProcessing/dataConversion/foamToVTK/internalWriter.C

@@ -27,9 +27,6 @@ License
 #include "internalWriter.H"
 #include "writeFuns.H"
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct from components

applications/utilities/postProcessing/dataConversion/foamToVTK/lagrangianWriter.C

@@ -29,9 +29,6 @@ License
 #include "Cloud.H"
 #include "passiveParticle.H"
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct from components

applications/utilities/postProcessing/dataConversion/foamToVTK/writePatchGeom.C

@@ -78,6 +78,7 @@ void writePatchGeom
     writeFuns::write(pStream, binary, vertLabels);
 }
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

applications/utilities/postProcessing/dataConversion/foamToVTK/writePointSet.C

@@ -103,6 +103,7 @@ void writePointSet
     writeFuns::write(pStream, binary, pointIDs);
 }
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/containers/Lists/DynamicList/DynamicList.C

@@ -33,7 +33,7 @@ template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
 Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::DynamicList(Istream& is)
 :
     List<T>(is),
-    nextFree_(List<T>::size())
+    allocSize_(List<T>::size())
 {}
 
 
@@ -44,9 +44,6 @@ Foam::Ostream& Foam::operator<<
     const Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>& DL
 )
 {
-    const_cast<DynamicList<T, SizeInc, SizeMult, SizeDiv>&>(DL)
-        .setSize(DL.nextFree_);
-
     os << static_cast<const List<T>&>(DL);
     return os;
 }
@@ -60,7 +57,7 @@ Foam::Istream& Foam::operator>>
 )
 {
     is >> static_cast<List<T>&>(DL);
-    DL.nextFree_ = DL.List<T>::size();
+    DL.allocSize_ = DL.List<T>::size();
 
     return is;
 }

src/OpenFOAM/containers/Lists/DynamicList/DynamicList.H

@@ -81,24 +81,29 @@ class DynamicList
 {
     // Private data
 
-        //- Number of next free element
+        //- Allocated size for underlying List.
-        label nextFree_;
+        label allocSize_;
 
 
 public:
 
+    // Related types
+
+        //- Declare friendship with the List class
+        friend class List<T>;
+
     // Constructors
 
         //- Construct null
         inline DynamicList();
 
-        //- Construct given size
+        //- Construct given size.
         explicit inline DynamicList(const label);
 
-        //- Construct from UList. nextFree_ set to size().
+        //- Construct from UList. Size set to UList size.
         explicit inline DynamicList(const UList<T>&);
 
-        //- Construct from Istream. nextFree_ set to size().
+        //- Construct from Istream. Size set to size of read list.
         explicit DynamicList(Istream&);
 
 
@@ -106,22 +111,24 @@ public:
 
         // Access
 
-            //- Size of the active part of the list.
+            //- Size of the underlying storage.
-            //  Direct over-ride of list size member function
+            inline label allocSize() const;
-            inline label size() const;
 
 
         // Edit
 
             //- Reset size of List.
-            void setSize(const label);
+            inline void setSize(const label);
 
             //- Reset size of List and value for new elements.
-            void setSize(const label, const T&);
+            inline void setSize(const label, const T&);
 
-            //- Clear the list, i.e. set next free to zero.
+            //- Clear the list, i.e. set the size to zero.
             //  Allocated size does not change
-            void clear();
+            inline void clear();
+
+            //- Clear the list and delete storage.
+            inline void clearStorage();
 
             //- Shrink the List<T> to the number of elements used
             inline DynamicList<T, SizeInc, SizeMult, SizeDiv>& shrink();
@@ -130,11 +137,11 @@ public:
             //  and annull the argument list. Is same as List::transfer except
             //  checks that you're not changing the underlying list to something
             //  smaller than nextFree_.
-            void transfer(List<T>&);
+            inline void transfer(List<T>&);
 
             //- Transfer the contents of the argument DynamicList into this
             //  DynamicList and annull the argument list.
-            void transfer(DynamicList<T, SizeInc, SizeMult, SizeDiv>&);
+            inline void transfer(DynamicList<T, SizeInc, SizeMult, SizeDiv>&);
 
 
     // Member Operators
@@ -152,9 +159,15 @@ public:
         //- Assignment of all entries to the given value
         inline void operator=(const T&);
 
-        //- Assignment to List<T>
+        //- Assignment from List<T>
         inline void operator=(const List<T>&);
 
+        //- Assignment from DynamicList<T>
+        inline void operator=
+        (
+            const DynamicList<T, SizeInc, SizeMult, SizeDiv>&
+        );
+
 
     // IOstream operators
 

src/OpenFOAM/containers/Lists/DynamicList/DynamicListI.H

@@ -31,8 +31,10 @@ template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
 inline Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::DynamicList()
 :
     List<T>(SizeInc),
-    nextFree_(0)
+    allocSize_(SizeInc)
-{}
+{
+    List<T>::size() = 0;
+}
 
 
 //- Construct given size
@@ -43,8 +45,10 @@ inline Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::DynamicList
 )
 :
     List<T>(s),
-    nextFree_(0)
+    allocSize_(s)
-{}
+{
+    List<T>::size() = 0;
+}
 
 
 //- Construct given size
@@ -55,17 +59,17 @@ inline Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::DynamicList
 )
 :
     List<T>(s),
-    nextFree_(s.size())
+    allocSize_(s.size())
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
-inline Foam::label Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::size()
+inline Foam::label Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::allocSize()
 const
 {
-    return nextFree_;
+    return allocSize_;
 }
 
 
@@ -75,13 +79,14 @@ inline void Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::setSize
     const label s
 )
 {
-    if (s < nextFree_)
+    if (s < List<T>::size())
     {
-        nextFree_ = s;
+        List<T>::size() = s;
     }
     else
     {
-        List<T>::setSize(s);
+        allocSize_ = s;
+        List<T>::setSize(allocSize_);
     }
 }
 
@@ -93,14 +98,14 @@ inline void Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::setSize
     const T& t
 )
 {
-    if (s < nextFree_)
+    if (s < List<T>::size())
     {
-        nextFree_ = s;
+        List<T>::size() = s;
     }
     else
     {
-        List<T>::setSize(s, t);
+        allocSize_ = s;
-        nextFree_ = s;
+        List<T>::setSize(allocSize_, t);
     }
 }
 
@@ -108,7 +113,16 @@ inline void Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::setSize
 template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
 inline void Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::clear()
 {
-    nextFree_ = 0;
+    List<T>::size() = 0;
+}
+
+
+template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+inline void Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::clearStorage()
+{
+    List<T>::size() = allocSize_;    // make List<T> consistent
+    List<T>::clear();
+    allocSize_ = 0;
 }
 
 
@@ -116,7 +130,8 @@ template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
 inline Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>&
 Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::shrink()
 {
-    List<T>::setSize(nextFree_);
+    allocSize_ = List<T>::size();
+    List<T>::setSize(allocSize_);
     return *this;
 }
 
@@ -125,20 +140,20 @@ template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
 inline void
 Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::transfer(List<T>& l)
 {
-    if (l.size() < nextFree_)
+    if (l.size() < List<T>::size())
     {
         FatalErrorIn
         (
             "void DynamicList<T, SizeInc, SizeMult"
             ", SizeDiv>::transfer(List<T>&)"
         )   << "Cannot replace the underlying storage of this DynamicList"
-            << " of which " << nextFree_ << " elements are used" << nl
+            << " of which " << List<T>::size() << " elements are used" << nl
             << "with a List of size " << l.size() << abort(FatalError);
     }
     else
     {
+        allocSize_ = l.size();
         List<T>::transfer(l);   // take over storage
-        l.clear();              // set nextFree of l to 0
     }
 }
 
@@ -150,40 +165,40 @@ Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::transfer
     DynamicList<T, SizeInc, SizeMult, SizeDiv>& l
 )
 {
+    allocSize_ = l.allocSize();
     List<T>::transfer(l);       // take over storage
-    nextFree_ = l.size();       // take over used size
+    l.allocSize_ = 0;
-    l.clear();                  // set nextFree of l to 0
 }
 
 
 template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
-inline void Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::append
+inline void Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::append(const T& e)
-(
-    const T& e
-)
 {
-    nextFree_++;
+    // Work on copy free index since gets overwritten by setSize
+    label nextFree = List<T>::size();
 
-    if (nextFree_ > List<T>::size())
+    nextFree++;
+
+    if (nextFree > allocSize_)
     {
-        List<T>::setSize
+        allocSize_ = max
         (
-            max
+            nextFree,
-            (
+            label(SizeMult*allocSize_/SizeDiv + SizeInc)
-                nextFree_,
-                label(SizeMult*List<T>::size()/SizeDiv + SizeInc)
-            )
         );
+        List<T>::setSize(allocSize_);
     }
 
-    this->operator[](nextFree_ - 1) = e;
+    this->operator[](nextFree - 1) = e;
+
+    List<T>::size() = nextFree;
 }
 
 
 template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
 inline T Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::remove()
 {
-    if (nextFree_ == 0)
+    if (List<T>::size() == 0)
     {
         FatalErrorIn
         (
@@ -191,7 +206,7 @@ inline T Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::remove()
         )   << "List is empty" << abort(FatalError);
     }
 
-    return List<T>::operator[](--nextFree_);
+    return List<T>::operator[](--List<T>::size());
 }
 
 
@@ -203,20 +218,22 @@ inline T& Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::operator()
     const label i
 )
 {
-    nextFree_ = max(nextFree_, i + 1);
+    label nextFree = List<T>::size();
 
-    if (nextFree_ > List<T>::size())
+    nextFree = max(nextFree, i + 1);
+
+    if (nextFree > allocSize_)
     {
-        List<T>::setSize
+       allocSize_ = max
         (
-            max
+            nextFree,
-            (
+            label(SizeMult*allocSize_/SizeDiv + SizeInc)
-                nextFree_,
-                label(SizeMult*List<T>::size()/SizeDiv + SizeInc)
-            )
         );
+        List<T>::setSize(allocSize_);
     }
 
+    List<T>::size() = nextFree;
+
     return this->operator[](i);
 }
 
@@ -228,7 +245,7 @@ inline void Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::operator=
 )
 {
     List<T>::operator=(t);
-    nextFree_ = List<T>::size();
+    allocSize_ = List<T>::size();
 }
 
 
@@ -239,7 +256,18 @@ inline void Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::operator=
 )
 {
     List<T>::operator=(l);
-    nextFree_ = l.size();
+    allocSize_ = List<T>::size();
+}
+
+
+template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+inline void Foam::DynamicList<T, SizeInc, SizeMult, SizeDiv>::operator=
+(
+    const DynamicList<T, SizeInc, SizeMult, SizeDiv>& l
+)
+{
+    List<T>::operator=(l);
+    allocSize_ = l.allocSize();
 }
 
 

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

@@ -420,6 +420,23 @@ void List<T>::transfer(List<T>& a)
 }
 
 
+// Transfer the contents of the argument DynamicList into this List
+// and anull the argument list
+template<class T>
+template<unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+void List<T>::transfer(DynamicList<T, SizeInc, SizeMult, SizeDiv>& a)
+{
+    if (this->v_) delete[] this->v_;
+
+    this->size_ = a.size_;
+    this->v_ = a.v_;
+
+    a.size_ = 0;
+    a.allocSize_ = 0;
+    a.v_ = 0;
+}
+
+
 template<class T>
 void sort(List<T>& a)
 {

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

@@ -61,6 +61,8 @@ template<class T> Istream& operator>>(Istream&, List<T>&);
 template<class T, label Size> class FixedList;
 template<class T> class PtrList;
 template<class T> class SLList;
+template<class T, unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+ class DynamicList;
 template<class T> class IndirectList;
 template<class T> class BiIndirectList;
 
@@ -156,6 +158,11 @@ public:
             //  and annull the argument list.
             void transfer(List<T>&);
 
+            //- Transfer the contents of the argument List into this List
+            //  and annull the argument list.
+            template<unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
+            void transfer(DynamicList<T, SizeInc, SizeMult, SizeDiv>&);
+
             //- Return subscript-checked element of UList.
             inline T& newElmt(const label);
 

src/OpenFOAM/containers/Lists/PackedList/PackedList.C

@@ -40,10 +40,7 @@ PackedList<nBits>::PackedList(const label size, const unsigned int val)
     List<unsigned int>(intSize(size)),
     size_(size)
 {
-    for (label i = 0; i < size; i++)
+    operator=(val);
-    {
-        set(i, val);
-    }
 }
 
 

src/OpenFOAM/containers/Lists/PackedList/PackedList.H

@@ -160,13 +160,15 @@ public:
             //- Number of packed elements
             inline label size() const;
 
-
             //- Get value at index I
             inline unsigned int get(const label i) const;
 
             //- Set value at index I. Return true if value changed.
             inline bool set(const label i, const unsigned int val);
 
+            //- Underlying storage
+            inline List<unsigned int>& storage();
+
 
     // Member operators
 
@@ -179,7 +181,8 @@ public:
             //- Assignment operator. Takes linear time.
             void operator=(const PackedList<nBits>&);
 
-            //- Assignment of all entries to the given value
+            //- Assignment of all entries to the given value. Does set on all
+            //  elements.
             inline void operator=(const unsigned int val);
 
             //- Return as labelList

src/OpenFOAM/containers/Lists/PackedList/PackedListI.H

@@ -83,7 +83,7 @@ inline void PackedList<nBits>::checkValue(const unsigned int val) const
 {
     if (val>=(1u << nBits))
     {
-        FatalErrorIn("PackedList<T>::set(const unsigned int)")
+        FatalErrorIn("PackedList<T>::checkValue(const unsigned int)")
             << "value " << label(val) << " out of range 0 ... "
             << label((1u << nBits)-1)
             << " representable by " << nBits << " bits"
@@ -183,6 +183,13 @@ inline bool PackedList<nBits>::set(const label i, const unsigned int val)
 }
 
 
+template<int nBits>
+inline List<unsigned int>& PackedList<nBits>::storage()
+{
+    return static_cast<List<unsigned int>&>(*this);
+}
+
+
 template<int nBits>
 inline ::Foam::reference PackedList<nBits>::operator[](const label i)
 {
@@ -213,9 +220,21 @@ inline void PackedList<nBits>::operator=(const unsigned int val)
 #   ifdef DEBUGList
     checkValue(val);
 #   endif
-    List<unsigned int>::operator=(val);
+
+    if (val == 0)
+    {
+        List<unsigned int>::operator=(val);
+    }
+    else
+    {
+        for (label i = 0; i < size_; i++)
+        {
+            set(i, val);
+        }
+    }
 }
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/containers/Lists/PackedList/PackedListIO.C

@@ -1,152 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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 "PackedList.H"
-#include "Ostream.H"
-#include "token.H"
-#include "contiguous.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * Ostream Operator *  * * * * * * * * * * * * //
-
-template<class T>
-void PackedList<T>::writeEntry(Ostream& os) const
-{
-    if
-    (
-        size()
-     && token::compound::isCompound
-        (
-            "List<" + word(pTraits<T>::typeName) + '>'
-        )
-    )
-    {
-        os  << word("List<" + word(pTraits<T>::typeName) + '>') << " ";
-    }
-    
-    os << *this;
-}
-
-
-template<class T>
-void PackedList<T>::writeEntry(const word& keyword, Ostream& os) const
-{
-    os.writeKeyword(keyword);
-    writeEntry(os);
-    os << token::END_STATEMENT << endl;
-}
-
-
-template<class T>
-Ostream& operator<<(Ostream& os, const PackedList<T>& L)
-{
-    // Write list contents depending on data format
-    if (os.format() == IOstream::ASCII || !contiguous<T>())
-    {
-        bool uniform = false;
-
-        if (L.size() > 1 && contiguous<T>())
-        {
-            uniform = true;
-
-            unsigned int L0 = L.get(0);
-
-            forAll(L, i)
-            {
-                if (L.get(i) != L0)
-                {
-                    uniform = false;
-                    break;
-                }
-            }
-        }
-
-        if (uniform)
-        {
-            // Write size of list and start contents delimiter
-            os << L.size() << token::BEGIN_BLOCK;
-
-            // Write list contents
-            os << L.get(0);
-
-            // Write end of contents delimiter
-            os << token::END_BLOCK;
-        }
-        else if (L.size() < 11 && contiguous<T>())
-        {
-            // Write size of list and start contents delimiter
-            os << L.size() << token::BEGIN_LIST;
-
-            // Write list contents
-            forAll(L, i)
-            {
-                if (i > 0) os << token::SPACE;
-                os << L.get(i);
-            }
-
-            // Write end of contents delimiter
-            os << token::END_LIST;
-        }
-        else
-        {
-            // Write size of list and start contents delimiter
-            os << nl << L.size() << nl << token::BEGIN_LIST;
-
-            // Write list contents
-            forAll(L, i)
-            {
-                os << nl << L.get(i);
-            }
-
-            // Write end of contents delimiter
-            os << nl << token::END_LIST << nl;
-        }
-    }
-    else
-    {
-        os << nl << L.size() << nl;
-        if (L.size())
-        {
-            os.write((const char*)L.v_, L.size()*sizeof(T));
-        }
-    }
-
-    // Check state of IOstream
-    os.check("Ostream& operator<<(Ostream&, const PackedList&)");
-
-    return os;
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// ************************************************************************* //

src/OpenFOAM/meshes/boundBox/boundBox.C

@@ -79,10 +79,9 @@ boundBox::boundBox(const pointField& points, const bool doReduce)
 
 
 boundBox::boundBox(Istream& is)
-:
+{
-    min_(is),
+    operator>>(is, *this);
-    max_(is)
+}
-{}
 
 
 // * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
@@ -113,7 +112,7 @@ Istream& operator>>(Istream& is, boundBox& bb)
 {
     if (is.format() == IOstream::ASCII)
     {
-        return is >> bb.min_ >>  bb.max_;
+        return is >> bb.min_ >> bb.max_;
     }
     else
     {

src/OpenFOAM/meshes/meshShapes/edge/edge.H

@@ -141,6 +141,9 @@ inline label Hash<edge>::operator()(const edge& e) const
     return e[0]*e[1] + e[0]+e[1];
 }
 
+template<>
+inline bool contiguous<edge>()  {return true;}
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/OpenFOAM/meshes/meshShapes/triFace/triFace.H

@@ -166,6 +166,9 @@ inline label Hash<triFace>::operator()(const triFace& t) const
     return (t[0]*t[1]*t[2] + t[0]+t[1]+t[2]);
 }
 
+template<>
+inline bool contiguous<triFace>()  {return true;}
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndex.C

@@ -26,12 +26,6 @@ License
 
 #include "globalIndex.H"
 
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::globalIndex::globalIndex(const label localSize)
@@ -68,15 +62,6 @@ Foam::globalIndex::globalIndex(Istream& is)
 }
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
-
-
 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 
 Foam::Istream& Foam::operator>>(Istream& is, globalIndex& gi)

src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C

@@ -397,13 +397,11 @@ void Foam::globalMeshData::calcSharedEdges() const
     sharedEdgeLabelsPtr_ = new labelList();
     labelList& sharedEdgeLabels = *sharedEdgeLabelsPtr_;
     sharedEdgeLabels.transfer(dynSharedEdgeLabels);
-    dynSharedEdgeLabels.clear();
 
     dynSharedEdgeAddr.shrink();
     sharedEdgeAddrPtr_ = new labelList();
     labelList& sharedEdgeAddr = *sharedEdgeAddrPtr_;
     sharedEdgeAddr.transfer(dynSharedEdgeAddr);
-    dynSharedEdgeAddr.clear();
 
     if (debug)
     {

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

@@ -84,16 +84,32 @@ void Foam::mapDistribute::distribute
         // Receive sub field from neighbour
         for (label domain = 0; domain < Pstream::nProcs(); domain++)
         {
-            if
+            const labelList& map = constructMap[domain];
-            (
+
-                domain != Pstream::myProcNo()
+            if (domain != Pstream::myProcNo() && map.size() > 0)
-             && constructMap[domain].size() > 0
-            )
             {
                 IPstream fromNbr(Pstream::blocking, domain);
                 List<T> subField(fromNbr);
 
-                const labelList& map = constructMap[domain];
+                if (subField.size() != map.size())
+                {
+                    FatalErrorIn
+                    (
+                        "template<class T>\n"
+                        "void mapDistribute::distribute\n"
+                        "(\n"
+                        "    const Pstream::commsTypes commsType,\n"
+                        "    const List<labelPair>& schedule,\n"
+                        "    const label constructSize,\n"
+                        "    const labelListList& subMap,\n"
+                        "    const labelListList& constructMap,\n"
+                        "    List<T>& field\n"
+                        ")\n"
+                    )   << "Expected from processor " << domain
+                        << " " << map.size() << " but received "
+                        << subField.size() << " elements."
+                        << abort(FatalError);
+                }
 
                 forAll(map, i)
                 {
@@ -126,6 +142,7 @@ void Foam::mapDistribute::distribute
             newField[map[i]] = subField[i];
         }
 
+        // Schedule will already have pruned 0-sized comms
         forAll(schedule, i)
         {
             const labelPair& twoProcs = schedule[i];
@@ -154,6 +171,26 @@ void Foam::mapDistribute::distribute
 
                 const labelList& map = constructMap[sendProc];
 
+                if (subField.size() != map.size())
+                {
+                    FatalErrorIn
+                    (
+                        "template<class T>\n"
+                        "void mapDistribute::distribute\n"
+                        "(\n"
+                        "    const Pstream::commsTypes commsType,\n"
+                        "    const List<labelPair>& schedule,\n"
+                        "    const label constructSize,\n"
+                        "    const labelListList& subMap,\n"
+                        "    const labelListList& constructMap,\n"
+                        "    List<T>& field\n"
+                        ")\n"
+                    )   << "Expected from processor " << sendProc
+                        << " " << map.size() << " but received "
+                        << subField.size() << " elements."
+                        << abort(FatalError);
+                }
+
                 forAll(map, i)
                 {
                     newField[map[i]] = subField[i];
@@ -164,7 +201,74 @@ void Foam::mapDistribute::distribute
     }
     else if (commsType == Pstream::nonBlocking)
     {
-        List<T> newField(constructSize);
+        if (!contiguous<T>())
+        {
+            FatalErrorIn
+            (
+                "template<class T>\n"
+                "void mapDistribute::distribute\n"
+                "(\n"
+                "    const Pstream::commsTypes commsType,\n"
+                "    const List<labelPair>& schedule,\n"
+                "    const label constructSize,\n"
+                "    const labelListList& subMap,\n"
+                "    const labelListList& constructMap,\n"
+                "    List<T>& field\n"
+                ")\n"
+            )   << "Non-blocking only supported for contiguous data."
+                << exit(FatalError);
+        }
+
+        // Set up sends to neighbours
+
+        List<List<T > > sendFields(Pstream::nProcs());
+
+        for (label domain = 0; domain < Pstream::nProcs(); domain++)
+        {
+            const labelList& map = subMap[domain];
+
+            if (domain != Pstream::myProcNo() && map.size() > 0)
+            {
+                List<T>& subField = sendFields[domain];
+                subField.setSize(map.size());
+                forAll(map, i)
+                {
+                    subField[i] = field[map[i]];
+                }
+
+                OPstream::write
+                (
+                    Pstream::nonBlocking,
+                    domain,
+                    reinterpret_cast<const char*>(subField.begin()),
+                    subField.size()
+                );
+            }
+        }
+
+        // Set up receives from neighbours
+
+        List<List<T > > recvFields(Pstream::nProcs());
+
+        for (label domain = 0; domain < Pstream::nProcs(); domain++)
+        {
+            const labelList& map = constructMap[domain];
+
+            if (domain != Pstream::myProcNo() && map.size() > 0)
+            {
+                recvFields[domain].setSize(map.size());
+                IPstream::read
+                (
+                    Pstream::nonBlocking,
+                    domain,
+                    reinterpret_cast<char*>(recvFields[domain].begin()),
+                    recvFields[domain].size()
+                );
+            }
+        }
+
+
+        // Combine bits. Note that can reuse field storage
 
         // Subset myself
         const labelList& mySubMap = subMap[Pstream::myProcNo()];
@@ -175,54 +279,56 @@ void Foam::mapDistribute::distribute
             subField[i] = field[mySubMap[i]];
         }
 
+        field.setSize(constructSize);
+
         // Receive sub field from myself (subField)
         const labelList& map = constructMap[Pstream::myProcNo()];
 
         forAll(map, i)
         {
-            newField[map[i]] = subField[i];
+            field[map[i]] = subField[i];
-        }
-
-        // Send sub field to neighbour
-        for (label domain = 0; domain < Pstream::nProcs(); domain++)
-        {
-            const labelList& map = subMap[domain];
-
-            if (domain != Pstream::myProcNo() && map.size() > 0)
-            {
-
-                List<T> subField(map.size());
-                forAll(map, i)
-                {
-                    subField[i] = field[map[i]];
-                }
-
-                OPstream toNbr(Pstream::nonBlocking, domain);
-                toNbr << subField;
-            }
         }
 
 
-        // Receive sub field from neighbour
+        // Wait for all to finish
+
+        OPstream::waitRequests();
+        IPstream::waitRequests();
+        
+        // Collect neighbour fields
+
         for (label domain = 0; domain < Pstream::nProcs(); domain++)
         {
             const labelList& map = constructMap[domain];
 
             if (domain != Pstream::myProcNo() && map.size() > 0)
             {
-                IPstream fromNbr(Pstream::nonBlocking, domain);
+                if (recvFields[domain].size() != map.size())
-                List<T> subField(fromNbr);
+                {
+                    FatalErrorIn
+                    (
+                        "template<class T>\n"
+                        "void mapDistribute::distribute\n"
+                        "(\n"
+                        "    const Pstream::commsTypes commsType,\n"
+                        "    const List<labelPair>& schedule,\n"
+                        "    const label constructSize,\n"
+                        "    const labelListList& subMap,\n"
+                        "    const labelListList& constructMap,\n"
+                        "    List<T>& field\n"
+                        ")\n"
+                    )   << "Expected from processor " << domain
+                        << " " << map.size() << " but received "
+                        << recvFields[domain].size() << " elements."
+                        << abort(FatalError);
+                }
 
                 forAll(map, i)
                 {
-                    newField[map[i]] = subField[i];
+                    field[map[i]] = recvFields[domain][i];
                 }
             }
         }
-        OPstream::waitRequests();
-        IPstream::waitRequests();
-
-        field.transfer(newField);
     }
     else
     {

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

@@ -221,6 +221,41 @@ public:
                 const bool applySeparation
             );
 
+        // PackedList versions
+
+            template <int nBits, class CombineOp>
+            static void syncFaceList
+            (
+                const polyMesh& mesh,
+                PackedList<nBits>& faceValues,
+                const CombineOp& cop
+            );
+
+            template <int nBits>
+            static void swapFaceList
+            (
+                const polyMesh& mesh,
+                PackedList<nBits>& faceValues
+            );
+
+            template <int nBits, class CombineOp>
+            static void syncPointList
+            (
+                const polyMesh& mesh,
+                PackedList<nBits>& pointValues,
+                const CombineOp& cop,
+                const unsigned int nullValue
+            );
+
+            template <int nBits, class CombineOp>
+            static void syncEdgeList
+            (
+                const polyMesh& mesh,
+                PackedList<nBits>& edgeValues,
+                const CombineOp& cop,
+                const unsigned int nullValue
+            );
+
         // Other
 
             //- Get per point whether is it master (of a coupled set of points)

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

@@ -940,17 +940,7 @@ void Foam::syncTools::syncPointList
                     fromNbr >> nbrPatchInfo;
                 }
                 // Null any value which is not on neighbouring processor
-                label nbrSize = nbrPatchInfo.size();
+                nbrPatchInfo.setSize(procPatch.nPoints(), nullValue);
-                nbrPatchInfo.setSize(procPatch.nPoints());
-                for
-                (
-                    label i = nbrSize;
-                    i < procPatch.nPoints();
-                    i++
-                )
-                {
-                    nbrPatchInfo[i] = nullValue;
-                }
 
                 if (!procPatch.parallel())
                 {
@@ -1225,18 +1215,7 @@ void Foam::syncTools::syncEdgeList
                     fromNeighb >> nbrPatchInfo;
                 }
                 // Null any value which is not on neighbouring processor
-                label nbrSize = nbrPatchInfo.size();
+                nbrPatchInfo.setSize(procPatch.nEdges(), nullValue);
-                nbrPatchInfo.setSize(procPatch.nEdges());
-                for
-                (
-                    label i = nbrSize;
-                    i < procPatch.nEdges();
-                    i++
-                )
-                {
-                    nbrPatchInfo[i] = nullValue;
-                }
-
 
                 if (!procPatch.parallel())
                 {
@@ -1591,4 +1570,469 @@ void Foam::syncTools::swapFaceList
 }
 
 
+template <int nBits, class CombineOp>
+void Foam::syncTools::syncFaceList
+(
+    const polyMesh& mesh,
+    PackedList<nBits>& faceValues,
+    const CombineOp& cop
+)
+{
+    if (faceValues.size() != mesh.nFaces())
+    {
+        FatalErrorIn
+        (
+            "syncTools<int nBits, class CombineOp>::syncFaceList"
+            "(const polyMesh&, PackedList<nBits>&, const CombineOp&)"
+        )   << "Number of values " << faceValues.size()
+            << " is not equal to the number of faces in the mesh "
+            << mesh.nFaces() << abort(FatalError);
+    }
+
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    if (!hasCouples(patches))
+    {
+        return;
+    }
+
+    // Patch data (proc patches only).
+    List<List<unsigned int> > patchValues(patches.size());
+
+    if (Pstream::parRun())
+    {
+        // Send
+
+        forAll(patches, patchI)
+        {
+            if
+            (
+                isA<processorPolyPatch>(patches[patchI])
+             && patches[patchI].size() > 0
+            )
+            {
+                const processorPolyPatch& procPatch =
+                    refCast<const processorPolyPatch>(patches[patchI]);
+
+                patchValues[patchI].setSize(procPatch.size());
+                forAll(procPatch, i)
+                {
+                    patchValues[patchI][i] =
+                        faceValues.get(procPatch.start()+i);
+                }
+
+                OPstream toNbr(Pstream::blocking, procPatch.neighbProcNo());
+                toNbr << patchValues[patchI];
+            }
+        }
+
+
+        // Receive and combine.
+
+        forAll(patches, patchI)
+        {
+            if
+            (
+                isA<processorPolyPatch>(patches[patchI])
+             && patches[patchI].size() > 0
+            )
+            {
+                const processorPolyPatch& procPatch =
+                    refCast<const processorPolyPatch>(patches[patchI]);
+
+                {
+                    IPstream fromNbr
+                    (
+                        Pstream::blocking,
+                        procPatch.neighbProcNo()
+                    );
+                    fromNbr >> patchValues[patchI];
+                }
+
+                // Combine (bitwise)
+                forAll(procPatch, i)
+                {
+                    unsigned int patchVal = patchValues[patchI][i];
+                    label meshFaceI = procPatch.start()+i;
+                    unsigned int faceVal = faceValues.get(meshFaceI);
+                    cop(faceVal, patchVal);
+                    faceValues.set(meshFaceI, faceVal);
+                }
+            }
+        }
+    }
+
+    // Do the cyclics.
+    forAll(patches, patchI)
+    {
+        if (isA<cyclicPolyPatch>(patches[patchI]))
+        {
+            const cyclicPolyPatch& cycPatch =
+                refCast<const cyclicPolyPatch>(patches[patchI]);
+
+            label half = cycPatch.size()/2;
+
+            for (label i = 0; i < half; i++)
+            {
+                label meshFace0 = cycPatch.start()+i;
+                unsigned int val0 = faceValues.get(meshFace0);
+                label meshFace1 = meshFace0 + half;
+                unsigned int val1 = faceValues.get(meshFace1);
+
+                unsigned int t = val0;
+                cop(t, val1);
+                faceValues.set(meshFace0, t);
+
+                cop(val1, val0);                
+                faceValues.set(meshFace1, val1);
+            }
+        }
+    }
+}
+
+
+template <int nBits>
+void Foam::syncTools::swapFaceList
+(
+    const polyMesh& mesh,
+    PackedList<nBits>& faceValues
+)
+{
+    syncFaceList(mesh, faceValues, eqOp<unsigned int>());
+}
+
+
+template <int nBits, class CombineOp>
+void Foam::syncTools::syncPointList
+(
+    const polyMesh& mesh,
+    PackedList<nBits>& pointValues,
+    const CombineOp& cop,
+    const unsigned int nullValue
+)
+{
+    if (pointValues.size() != mesh.nPoints())
+    {
+        FatalErrorIn
+        (
+            "syncTools<int nBits, class CombineOp>::syncPointList"
+            "(const polyMesh&, PackedList<nBits>&, const CombineOp&"
+            ", const unsigned int&)"
+        )   << "Number of values " << pointValues.size()
+            << " is not equal to the number of points in the mesh "
+            << mesh.nPoints() << abort(FatalError);
+    }
+
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    if (!hasCouples(patches))
+    {
+        return;
+    }
+
+    // Patch data (proc patches only).
+    List<List<unsigned int> > patchValues(patches.size());
+
+    if (Pstream::parRun())
+    {
+        // Send
+
+        forAll(patches, patchI)
+        {
+            if
+            (
+                isA<processorPolyPatch>(patches[patchI])
+             && patches[patchI].nPoints() > 0
+            )
+            {
+                const processorPolyPatch& procPatch =
+                    refCast<const processorPolyPatch>(patches[patchI]);
+
+                patchValues[patchI].setSize(procPatch.nPoints());
+                patchValues[patchI] = nullValue;
+
+                const labelList& meshPts = procPatch.meshPoints();
+                const labelList& nbrPts = procPatch.neighbPoints();
+
+                forAll(nbrPts, pointI)
+                {
+                    label nbrPointI = nbrPts[pointI];
+                    if (nbrPointI >= 0 && nbrPointI < procPatch.nPoints())
+                    {
+                        patchValues[patchI][nbrPointI] =
+                            pointValues.get(meshPts[pointI]);
+                    }
+                }
+
+                OPstream toNbr(Pstream::blocking, procPatch.neighbProcNo());
+                toNbr << patchValues[patchI];
+            }
+        }
+
+
+        // Receive and combine.
+
+        forAll(patches, patchI)
+        {
+            if
+            (
+                isA<processorPolyPatch>(patches[patchI])
+             && patches[patchI].nPoints() > 0
+            )
+            {
+                const processorPolyPatch& procPatch =
+                    refCast<const processorPolyPatch>(patches[patchI]);
+
+                {
+                    // We do not know the number of points on the other side
+                    // so cannot use Pstream::read.
+                    IPstream fromNbr
+                    (
+                        Pstream::blocking,
+                        procPatch.neighbProcNo()
+                    );
+                    fromNbr >> patchValues[patchI];
+                }
+
+                // Null any value which is not on neighbouring processor
+                patchValues[patchI].setSize(procPatch.nPoints(), nullValue);
+
+                const labelList& meshPts = procPatch.meshPoints();
+
+                forAll(meshPts, pointI)
+                {
+                    label meshPointI = meshPts[pointI];
+                    unsigned int pointVal = pointValues.get(meshPointI);
+                    cop(pointVal, patchValues[patchI][pointI]);
+                    pointValues.set(meshPointI, pointVal);
+                }
+            }
+        }
+    }
+
+    // Do the cyclics.
+    forAll(patches, patchI)
+    {
+        if (isA<cyclicPolyPatch>(patches[patchI]))
+        {
+            const cyclicPolyPatch& cycPatch =
+                refCast<const cyclicPolyPatch>(patches[patchI]);
+
+            const edgeList& coupledPoints = cycPatch.coupledPoints();
+            const labelList& meshPts = cycPatch.meshPoints();
+
+            forAll(coupledPoints, i)
+            {
+                const edge& e = coupledPoints[i];
+
+                label point0 = meshPts[e[0]];
+                label point1 = meshPts[e[1]];
+
+                unsigned int val0 = pointValues.get(point0);
+                unsigned int t = val0;
+                unsigned int val1 = pointValues.get(point1);
+
+                cop(t, val1);
+                pointValues.set(point0, t);
+                cop(val1, val0);
+                pointValues.set(point1, val1);
+            }
+        }
+    }
+
+    // Synchronize multiple shared points.
+    const globalMeshData& pd = mesh.globalData();
+
+    if (pd.nGlobalPoints() > 0)
+    {
+        // Values on shared points. Use unpacked storage for ease!
+        List<unsigned int> sharedPts(pd.nGlobalPoints(), nullValue);
+
+        forAll(pd.sharedPointLabels(), i)
+        {
+            label meshPointI = pd.sharedPointLabels()[i];
+            // Fill my entries in the shared points
+            sharedPts[pd.sharedPointAddr()[i]] = pointValues.get(meshPointI);
+        }
+
+        // Combine on master.
+        Pstream::listCombineGather(sharedPts, cop);
+        Pstream::listCombineScatter(sharedPts);
+
+        // Now we will all have the same information. Merge it back with
+        // my local information.
+        forAll(pd.sharedPointLabels(), i)
+        {
+            label meshPointI = pd.sharedPointLabels()[i];
+            pointValues.set(meshPointI, sharedPts[pd.sharedPointAddr()[i]]);
+        }
+    }
+}
+
+
+template <int nBits, class CombineOp>
+void Foam::syncTools::syncEdgeList
+(
+    const polyMesh& mesh,
+    PackedList<nBits>& edgeValues,
+    const CombineOp& cop,
+    const unsigned int nullValue
+)
+{
+    if (edgeValues.size() != mesh.nEdges())
+    {
+        FatalErrorIn
+        (
+            "syncTools<int nBits, class CombineOp>::syncEdgeList"
+            "(const polyMesh&, PackedList<nBits>&, const CombineOp&"
+            ", const unsigned int&)"
+        )   << "Number of values " << edgeValues.size()
+            << " is not equal to the number of edges in the mesh "
+            << mesh.nEdges() << abort(FatalError);
+    }
+
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    if (!hasCouples(patches))
+    {
+        return;
+    }
+
+    // Patch data (proc patches only).
+    List<List<unsigned int> > patchValues(patches.size());
+
+    if (Pstream::parRun())
+    {
+        // Send
+
+        forAll(patches, patchI)
+        {
+            if
+            (
+                isA<processorPolyPatch>(patches[patchI])
+             && patches[patchI].nEdges() > 0
+            )
+            {
+                const processorPolyPatch& procPatch =
+                    refCast<const processorPolyPatch>(patches[patchI]);
+
+                patchValues[patchI].setSize(procPatch.nEdges(), nullValue);
+
+                const labelList& meshEdges = procPatch.meshEdges();
+                const labelList& neighbEdges = procPatch.neighbEdges();
+
+                forAll(neighbEdges, edgeI)
+                {
+                    label nbrEdgeI = neighbEdges[edgeI];
+                    if (nbrEdgeI >= 0 && nbrEdgeI < procPatch.nEdges())
+                    {
+                        patchValues[patchI][nbrEdgeI] =
+                            edgeValues.get(meshEdges[edgeI]);
+                    }
+                }
+
+                OPstream toNbr(Pstream::blocking, procPatch.neighbProcNo());
+                toNbr << patchValues[patchI];
+            }
+        }
+
+
+        // Receive and combine.
+
+        forAll(patches, patchI)
+        {
+            if
+            (
+                isA<processorPolyPatch>(patches[patchI])
+             && patches[patchI].nEdges() > 0
+            )
+            {
+                const processorPolyPatch& procPatch =
+                    refCast<const processorPolyPatch>(patches[patchI]);
+
+                {
+                    IPstream fromNeighb
+                    (
+                        Pstream::blocking,
+                        procPatch.neighbProcNo()
+                    );
+                    fromNeighb >> patchValues[patchI];
+                }
+
+                patchValues[patchI].setSize(procPatch.nEdges(), nullValue);
+
+                const labelList& meshEdges = procPatch.meshEdges();
+
+                forAll(meshEdges, edgeI)
+                {
+                    unsigned int patchVal = patchValues[patchI][edgeI];
+                    label meshEdgeI = meshEdges[edgeI];
+                    unsigned int edgeVal = edgeValues.get(meshEdgeI);
+                    cop(edgeVal, patchVal);
+                    edgeValues.set(meshEdgeI, edgeVal);
+                }
+            }
+        }
+    }
+
+    // Do the cyclics.
+    forAll(patches, patchI)
+    {
+        if (isA<cyclicPolyPatch>(patches[patchI]))
+        {
+            const cyclicPolyPatch& cycPatch =
+                refCast<const cyclicPolyPatch>(patches[patchI]);
+
+            const edgeList& coupledEdges = cycPatch.coupledEdges();
+            const labelList& meshEdges = cycPatch.meshEdges();
+
+            forAll(coupledEdges, i)
+            {
+                const edge& e = coupledEdges[i];
+
+                label edge0 = meshEdges[e[0]];
+                label edge1 = meshEdges[e[1]];
+
+                unsigned int val0 = edgeValues.get(edge0);
+                unsigned int t = val0;
+                unsigned int val1 = edgeValues.get(edge1);
+
+                cop(t, val1);
+                edgeValues.set(edge0, t);
+                cop(val1, val0);
+                edgeValues.set(edge1, val1);
+            }
+        }
+    }
+
+    // Synchronize multiple shared edges.
+    const globalMeshData& pd = mesh.globalData();
+
+    if (pd.nGlobalEdges() > 0)
+    {
+        // Values on shared edges. Use unpacked storage for ease!
+        List<unsigned int> sharedPts(pd.nGlobalEdges(), nullValue);
+
+        forAll(pd.sharedEdgeLabels(), i)
+        {
+            label meshEdgeI = pd.sharedEdgeLabels()[i];
+            // Fill my entries in the shared edges
+            sharedPts[pd.sharedEdgeAddr()[i]] = edgeValues.get(meshEdgeI);
+        }
+
+        // Combine on master.
+        Pstream::listCombineGather(sharedPts, cop);
+        Pstream::listCombineScatter(sharedPts);
+
+        // Now we will all have the same information. Merge it back with
+        // my local information.
+        forAll(pd.sharedEdgeLabels(), i)
+        {
+            label meshEdgeI = pd.sharedEdgeLabels()[i];
+            edgeValues.set(meshEdgeI, sharedPts[pd.sharedEdgeAddr()[i]]);
+        }
+    }
+}
+
+
 // ************************************************************************* //

src/OpenFOAM/meshes/primitiveMesh/PrimitivePatchExtra/PrimitivePatchExtraSearch.C

@@ -102,8 +102,9 @@ void Foam::PrimitivePatchExtra<Face, ListType, PointField, PointType>::markZone
             break;
         }
 
-        changedFaces.transfer(newChangedFaces.shrink());
+        // New dynamiclist: can leave dynamicList unshrunk
-        newChangedFaces.clear();
+        //changedFaces.transfer(newChangedFaces.shrink());
+        changedFaces.transfer(newChangedFaces);
     }
 }
 

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdges.C

@@ -441,11 +441,10 @@ void primitiveMesh::calcEdges(const bool doFaceEdges) const
         forAll(pe, pointI)
         {
             DynamicList<label>& pEdges = pe[pointI];
-            inplaceRenumber(oldToNew, pEdges);
             pEdges.shrink();
+            inplaceRenumber(oldToNew, pEdges);
             pointEdges[pointI].transfer(pEdges);
             Foam::sort(pointEdges[pointI]);
-            pEdges.clear();
         }
 
         // faceEdges

src/OpenFOAM/meshes/primitiveMesh/primitivePatch/primitivePatch.H

@@ -26,7 +26,7 @@ Typedef
     Foam::primitivePatch
 
 Description
-    Foam::primitivePatch
+    Addressing for a faceList slice.
 
 \*---------------------------------------------------------------------------*/
 

src/OpenFOAM/meshes/primitiveShapes/line/linePointRef.H

@@ -26,6 +26,7 @@ Typedef
     Foam::linePointRef
 
 Description
+    Line using referred points
 
 \*---------------------------------------------------------------------------*/
 

src/OpenFOAM/primitives/contiguous/contiguous.H

@@ -38,6 +38,8 @@ Description
 #ifndef contiguous_H
 #define contiguous_H
 
+#include "label.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -45,54 +47,117 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+// Forward declaration of friend functions and operators
+template<class T, label Size> class FixedList;
+template<class T> class Pair;
+
+
 // Assume the data associated with type T is not contiguous
 template<class T>
 inline bool contiguous()                        {return false;}
 
 
-// Specify data associated with primitive types is contiguous
+
+// Specify data associated with primitive types (and simple fixed size
+// containers - only size 2 defined here) is contiguous
 
 template<>
-inline bool contiguous<bool>()                  {return true;}
+inline bool contiguous<bool>()                              {return true;}
+template<>
+inline bool contiguous<FixedList<bool, 2> >()               {return true;}
+template<>
+inline bool contiguous<Pair<bool> >()                       {return true;}
 
 template<>
-inline bool contiguous<char>()                  {return true;}
+inline bool contiguous<char>()                              {return true;}
+template<>
+inline bool contiguous<FixedList<char, 2> >()               {return true;}
+template<>
+inline bool contiguous<Pair<char> >()                       {return true;}
 
 template<>
-inline bool contiguous<unsigned char>()         {return true;}
+inline bool contiguous<unsigned char>()                     {return true;}
+template<>
+inline bool contiguous<FixedList<unsigned char, 2> >()      {return true;}
+template<>
+inline bool contiguous<Pair<unsigned char> >()              {return true;}
 
 template<>
-inline bool contiguous<short>()                 {return true;}
+inline bool contiguous<short>()                             {return true;}
+template<>
+inline bool contiguous<FixedList<short, 2> >()              {return true;}
+template<>
+inline bool contiguous<Pair<short> >()                      {return true;}
 
 template<>
-inline bool contiguous<unsigned short>()        {return true;}
+inline bool contiguous<unsigned short>()                    {return true;}
+template<>
+inline bool contiguous<FixedList<unsigned short, 2> >()     {return true;}
+template<>
+inline bool contiguous<Pair<unsigned short> >()             {return true;}
 
 template<>
-inline bool contiguous<int>()                   {return true;}
+inline bool contiguous<int>()                               {return true;}
+template<>
+inline bool contiguous<FixedList<int, 2> >()                {return true;}
+template<>
+inline bool contiguous<Pair<int> >()                        {return true;}
 
 template<>
-inline bool contiguous<unsigned int>()          {return true;}
+inline bool contiguous<unsigned int>()                      {return true;}
+template<>
+inline bool contiguous<FixedList<unsigned int, 2> >()       {return true;}
+template<>
+inline bool contiguous<Pair<unsigned int> >()               {return true;}
 
 template<>
-inline bool contiguous<long>()                  {return true;}
+inline bool contiguous<long>()                              {return true;}
+template<>
+inline bool contiguous<FixedList<long, 2> >()               {return true;}
+template<>
+inline bool contiguous<Pair<long> >()                       {return true;}
 
 template<>
-inline bool contiguous<unsigned long>()         {return true;}
+inline bool contiguous<unsigned long>()                     {return true;}
+template<>
+inline bool contiguous<FixedList<unsigned long, 2> >()      {return true;}
+template<>
+inline bool contiguous<Pair<unsigned long> >()              {return true;}
 
 template<>
-inline bool contiguous<long long>()             {return true;}
+inline bool contiguous<long long>()                         {return true;}
+template<>
+inline bool contiguous<FixedList<long long, 2> >()          {return true;}
+template<>
+inline bool contiguous<Pair<long long> >()                  {return true;}
 
 template<>
-inline bool contiguous<unsigned long long>()    {return true;}
+inline bool contiguous<unsigned long long>()                {return true;}
+template<>
+inline bool contiguous<FixedList<unsigned long long, 2> >() {return true;}
+template<>
+inline bool contiguous<Pair<unsigned long long> >()         {return true;}
 
 template<>
-inline bool contiguous<float>()                 {return true;}
+inline bool contiguous<float>()                             {return true;}
+template<>
+inline bool contiguous<FixedList<float, 2> >()              {return true;}
+template<>
+inline bool contiguous<Pair<float> >()                       {return true;}
 
 template<>
-inline bool contiguous<double>()                {return true;}
+inline bool contiguous<double>()                            {return true;}
+template<>
+inline bool contiguous<FixedList<double, 2> >()             {return true;}
+template<>
+inline bool contiguous<Pair<double> >()                     {return true;}
 
 template<>
-inline bool contiguous<long double>()           {return true;}
+inline bool contiguous<long double>()                       {return true;}
+template<>
+inline bool contiguous<FixedList<long double, 2> >()        {return true;}
+template<>
+inline bool contiguous<Pair<long double> >()                {return true;}
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C

@@ -25,6 +25,7 @@ License
 \*----------------------------------------------------------------------------*/
 
 #include "autoRefineDriver.H"
+#include "meshRefinement.H"
 #include "fvMesh.H"
 #include "Time.H"
 #include "boundBox.H"

src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H

@@ -35,7 +35,7 @@ SourceFiles
 #ifndef autoRefineDriver_H
 #define autoRefineDriver_H
 
-#include "meshRefinement.H"
+#include "treeBoundBox.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -45,6 +45,9 @@ namespace Foam
 // Forward declaration of classes
 class featureEdgeMesh;
 class refinementParameters;
+class meshRefinement;
+class decompositionMethod;
+class fvMeshDistribute;
 
 /*---------------------------------------------------------------------------*\
                            Class autoRefineDriver Declaration

src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C

@@ -51,6 +51,9 @@ License
 #include "meshTools.H"
 #include "OFstream.H"
 #include "geomDecomp.H"
+#include "Random.H"
+#include "searchableSurfaces.H"
+#include "treeBoundBox.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -122,13 +125,36 @@ void Foam::meshRefinement::updateIntersections(const labelList& changedFaces)
 {
     const pointField& cellCentres = mesh_.cellCentres();
 
-    label nTotEdges = returnReduce(surfaceIndex_.size(), sumOp<label>());
+    // Stats on edges to test. Count proc faces only once.
-    label nChangedFaces = returnReduce(changedFaces.size(), sumOp<label>());
+    PackedList<1> isMasterFace(syncTools::getMasterFaces(mesh_));
+
+    {
+        label nMasterFaces = 0;
+        forAll(isMasterFace, faceI)
+        {
+            if (isMasterFace.get(faceI) == 1)
+            {
+                nMasterFaces++;
+            }
+        }
+        reduce(nMasterFaces, sumOp<label>());
+
+        label nChangedFaces = 0;
+        forAll(changedFaces, i)
+        {
+            if (isMasterFace.get(changedFaces[i]) == 1)
+            {
+                nChangedFaces++;
+            }
+        }
+        reduce(nChangedFaces, sumOp<label>());
+
+        Info<< "Edge intersection testing:" << nl
+            << "    Number of edges             : " << nMasterFaces << nl
+            << "    Number of edges to retest   : " << nChangedFaces
+            << endl;
+    }
 
-    Info<< "Edge intersection testing:" << nl
-        << "    Number of edges             : " << nTotEdges << nl
-        << "    Number of edges to retest   : " << nChangedFaces
-        << endl;
 
     // Get boundary face centre and level. Coupled aware.
     labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
@@ -838,11 +864,14 @@ Foam::meshRefinement::meshRefinement
 
 Foam::label Foam::meshRefinement::countHits() const
 {
+    // Stats on edges to test. Count proc faces only once.
+    PackedList<1> isMasterFace(syncTools::getMasterFaces(mesh_));
+
     label nHits = 0;
 
     forAll(surfaceIndex_, faceI)
     {
-        if (surfaceIndex_[faceI] >= 0)
+        if (surfaceIndex_[faceI] >= 0 && isMasterFace.get(faceI) == 1)
         {
             nHits++;
         }
@@ -996,11 +1025,6 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
 
     if (Pstream::parRun())
     {
-        //Info<< nl
-        //    << "Doing final balancing" << nl
-        //    << "---------------------" << nl
-        //    << endl;
-        //
         //if (debug_)
         //{
         //    const_cast<Time&>(mesh_.time())++;
@@ -1014,17 +1038,12 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
             // Faces where owner and neighbour are not 'connected' so can
             // go to different processors.
             boolList blockedFace(mesh_.nFaces(), true);
+            label nUnblocked = 0;
             // Pairs of baffles
             List<labelPair> couples;
 
             if (keepZoneFaces)
             {
-                label nNamed = surfaces().getNamedSurfaces().size();
-
-                Info<< "Found " << nNamed << " surfaces with faceZones."
-                    << " Applying special decomposition to keep those together."
-                    << endl;
-
                 // Determine decomposition to keep/move surface zones
                 // on one processor. The reason is that snapping will make these
                 // into baffles, move and convert them back so if they were
@@ -1039,8 +1058,6 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
                 // Get faces whose owner and neighbour should stay together,
                 // i.e. they are not 'blocked'.
 
-                label nZoned = 0;
-
                 forAll(fzNames, surfI)
                 {
                     if (fzNames[surfI].size() > 0)
@@ -1055,14 +1072,18 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
                             if (blockedFace[fZone[i]])
                             {
                                 blockedFace[fZone[i]] = false;
-                                nZoned++;
+                                nUnblocked++;
                             }
                         }
                     }
                 }
-                Info<< "Found " << returnReduce(nZoned, sumOp<label>())
+            }
-                    << " zoned faces to keep together."
+            reduce(nUnblocked, sumOp<label>());
-                    << endl;
+
+            if (keepZoneFaces)
+            {
+                Info<< "Found " << nUnblocked
+                    << " zoned faces to keep together." << endl;
             }
 
             if (keepBaffles)
@@ -1073,29 +1094,43 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
                     identity(mesh_.nFaces()-mesh_.nInternalFaces())
                    +mesh_.nInternalFaces()
                 );
+            }
+            label nCouples = returnReduce(couples.size(), sumOp<label>());
 
-                Info<< "Found " << returnReduce(couples.size(), sumOp<label>())
+            if (keepBaffles)
-                    << " baffles to keep together."
+            {
+                Info<< "Found " << nCouples << " baffles to keep together."
                     << endl;
             }
 
-            distribution = decomposeCombineRegions
+            if (nUnblocked > 0 || nCouples > 0)
-            (
-                blockedFace,
-                couples,
-                decomposer
-            );
-
-            labelList nProcCells(distributor.countCells(distribution));
-            Pstream::listCombineGather(nProcCells, plusEqOp<label>());
-            Pstream::listCombineScatter(nProcCells);
-
-            Info<< "Calculated decomposition:" << endl;
-            forAll(nProcCells, procI)
             {
-                Info<< "    " << procI << '\t' << nProcCells[procI] << endl;
+                Info<< "Applying special decomposition to keep baffles"
+                    << " and zoned faces together." << endl;
+
+                distribution = decomposeCombineRegions
+                (
+                    blockedFace,
+                    couples,
+                    decomposer
+                );
+
+                labelList nProcCells(distributor.countCells(distribution));
+                Pstream::listCombineGather(nProcCells, plusEqOp<label>());
+                Pstream::listCombineScatter(nProcCells);
+
+                Info<< "Calculated decomposition:" << endl;
+                forAll(nProcCells, procI)
+                {
+                    Info<< "    " << procI << '\t' << nProcCells[procI] << endl;
+                }
+                Info<< endl;
+            }
+            else
+            {
+                // Normal decomposition
+                distribution = decomposer.decompose(mesh_.cellCentres());
             }
-            Info<< endl;
         }
         else
         {
@@ -1715,6 +1750,36 @@ void Foam::meshRefinement::distribute(const mapDistributePolyMesh& map)
     {
         map.distributeFaceData(userFaceData_[i].second());
     }
+
+    // Redistribute surface and any fields on it.
+    {
+        Random rndGen(653213);
+
+        // Get local mesh bounding box. Single box for now.
+        List<treeBoundBox> meshBb(1);
+        treeBoundBox& bb = meshBb[0];
+        bb = boundBox(mesh_.points(), false);
+        bb = bb.extend(rndGen, 1E-4);
+
+        // Distribute all geometry (so refinementSurfaces and shellSurfaces)
+        searchableSurfaces& geometry =
+            const_cast<searchableSurfaces&>(surfaces_.geometry());
+
+        forAll(geometry, i)
+        {
+            autoPtr<mapDistribute> faceMap;
+            autoPtr<mapDistribute> pointMap;
+            geometry[i].distribute
+            (
+                meshBb,
+                false,          // do not keep outside triangles
+                faceMap,
+                pointMap
+            );
+            faceMap.clear();
+            pointMap.clear();
+        }
+    }
 }
 
 

src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C

@@ -386,7 +386,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createBaffles
         FatalErrorIn
         (
             "meshRefinement::createBaffles"
-            "(const label, const labelList&, const labelList&)"
+            "(const labelList&, const labelList&)"
         )   << "Illegal size :"
             << " ownPatch:" << ownPatch.size()
             << " neiPatch:" << neiPatch.size()
@@ -412,7 +412,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createBaffles
                 FatalErrorIn
                 (
                     "meshRefinement::createBaffles"
-                    "(const label, const labelList&, const labelList&)"
+                    "(const labelList&, const labelList&)"
                 )   << "Non synchronised at face:" << faceI
                     << " on patch:" << mesh_.boundaryMesh().whichPatch(faceI)
                     << " fc:" << mesh_.faceCentres()[faceI] << endl
@@ -461,7 +461,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createBaffles
 
     //- Redo the intersections on the newly create baffle faces. Note that
     //  this changes also the cell centre positions.
-    labelHashSet baffledFacesSet(2*nBaffles);
+    faceSet baffledFacesSet(mesh_, "baffledFacesSet", 2*nBaffles);
 
     const labelList& reverseFaceMap = map().reverseFaceMap();
     const labelList& faceMap = map().faceMap();
@@ -496,6 +496,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createBaffles
             }
         }
     }
+    baffledFacesSet.sync(mesh_);
 
     updateMesh(map, baffledFacesSet.toc());
 
@@ -1886,15 +1887,15 @@ void Foam::meshRefinement::baffleAndSplitMesh
 
         labelList facePatch
         (
-            //markFacesOnProblemCells
+            markFacesOnProblemCells
-            //(
-            //    globalToPatch
-            //)
-            markFacesOnProblemCellsGeometric
             (
-                motionDict,
                 globalToPatch
             )
+            //markFacesOnProblemCellsGeometric
+            //(
+            //    motionDict,
+            //    globalToPatch
+            //)
         );
         Info<< "Analyzed problem cells in = "
             << runTime.cpuTimeIncrement() << " s\n" << nl << endl;

src/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C

@@ -1091,8 +1091,8 @@ Foam::labelList Foam::meshRefinement::refineCandidates
         //    << "    local allowable refinement:" << nAllowRefine << nl
         //    << endl;
 
-        //- Disable refinement shortcut
+        //- Disable refinement shortcut. nAllowRefine is per processor limit.
-        label nAllowRefine = labelMax;
+        label nAllowRefine = labelMax / Pstream::nProcs();
 
         // Marked for refinement (>= 0) or not (-1). Actual value is the
         // index of the surface it intersects.

src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C

@@ -25,19 +25,13 @@ License
 \*----------------------------------------------------------------------------*/
 
 #include "refinementSurfaces.H"
-#include "orientedSurface.H"
 #include "Time.H"
 #include "searchableSurfaces.H"
 #include "shellSurfaces.H"
 #include "triSurfaceMesh.H"
 #include "labelPair.H"
 #include "searchableSurfacesQueries.H"
-
+#include "UPtrList.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -438,31 +432,6 @@ Foam::labelList Foam::refinementSurfaces::getClosedNamedSurfaces() const
 }
 
 
-// Orient surface (if they're closed) before any searching is done.
-void Foam::refinementSurfaces::orientSurface
-(
-    const point& keepPoint,
-    triSurfaceMesh& s
-)
-{
-    // Flip surface so keepPoint is outside.
-    bool anyFlipped = orientedSurface::orient(s, keepPoint, 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<< "orientSurfaces : Flipped orientation of surface "
-            << s.searchableSurface::name()
-            << " so point " << keepPoint << " is outside." << endl;
-    }
-}
-
-
 // Count number of triangles per surface region
 Foam::labelList Foam::refinementSurfaces::countRegions(const triSurface& s)
 {
@@ -485,7 +454,7 @@ void Foam::refinementSurfaces::setMinLevelFields
     const shellSurfaces& shells
 )
 {
-    minLevelFields_.setSize(surfaces_.size());
+    //minLevelFields_.setSize(surfaces_.size());
 
     forAll(surfaces_, surfI)
     {
@@ -495,26 +464,24 @@ void Foam::refinementSurfaces::setMinLevelFields
         {
             const triSurfaceMesh& triMesh = refCast<const triSurfaceMesh>(geom);
 
-            minLevelFields_.set
+            autoPtr<triSurfaceLabelField> minLevelFieldPtr
             (
-                surfI,
                 new triSurfaceLabelField
                 (
                     IOobject
                     (
-                        triMesh.searchableSurface::name(),
+                        "minLevel",
                         triMesh.objectRegistry::time().constant(),// directory
                         "triSurface",               // instance
                         triMesh,
                         IOobject::NO_READ,
-                        IOobject::AUTO_WRITE,
+                        IOobject::AUTO_WRITE
-                        false
                     ),
                     triMesh,
                     dimless
                 )
             );
-            triSurfaceLabelField& minLevelField = minLevelFields_[surfI];
+            triSurfaceLabelField& minLevelField = minLevelFieldPtr();
 
             const triSurface& s = static_cast<const triSurface&>(triMesh);
 
@@ -542,6 +509,9 @@ void Foam::refinementSurfaces::setMinLevelFields
                     shellLevel[triI]
                 );
             }
+
+            // Store field on triMesh
+            minLevelFieldPtr.ptr()->store();
         }
     }
 }
@@ -569,17 +539,38 @@ void Foam::refinementSurfaces::findHigherIntersection
         return;
     }
 
+
+    // Precalculate per surface whether it has a minlevelfield
+    UPtrList<triSurfaceLabelField> minLevelFields(surfaces_.size());
+    forAll(surfaces_, surfI)
+    {
+        const searchableSurface& geom = allGeometry_[surfaces_[surfI]];
+
+        if (isA<triSurfaceMesh>(geom))
+        {
+            const triSurfaceMesh& triMesh = refCast<const triSurfaceMesh>(geom);
+            minLevelFields.set
+            (
+                surfI,
+               &const_cast<triSurfaceLabelField&>
+                (
+                    triMesh.lookupObject<triSurfaceLabelField>("minLevel")
+                )
+            );
+        }
+    }
+
     // Work arrays
     labelList hitMap(identity(start.size()));
     pointField p0(start);
     pointField p1(end);
     List<pointIndexHit> hitInfo(start.size());
 
-
     forAll(surfaces_, surfI)
     {
         allGeometry_[surfaces_[surfI]].findLineAny(p0, p1, hitInfo);
 
+
         // Copy all hits into arguments, continue with misses
         label newI = 0;
         forAll(hitInfo, hitI)
@@ -592,12 +583,12 @@ void Foam::refinementSurfaces::findHigherIntersection
                 // Check if minLevelField for this surface.
                 if
                 (
-                    minLevelFields_.set(surfI)
+                    minLevelFields.set(surfI)
-                 && minLevelFields_[surfI].size() > 0
+                 && minLevelFields[surfI].size() > 0
                 )
                 {
                     minLocalLevel =
-                        minLevelFields_[surfI][hitInfo[hitI].index()];
+                        minLevelFields[surfI][hitInfo[hitI].index()];
                 }
                 else
                 {
@@ -668,30 +659,57 @@ void Foam::refinementSurfaces::findAllHigherIntersections
     {
         allGeometry_[surfaces_[surfI]].findLineAll(start, end, hitInfo);
 
+        // Repack hits for surface into flat list
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        // To avoid overhead of calling getRegion for every point
+
+        label n = 0;
+        forAll(hitInfo, pointI)
+        {
+            n += hitInfo[pointI].size();
+        }
+
+        List<pointIndexHit> surfInfo(n);
+        labelList pointMap(n);
+        n = 0;
+
         forAll(hitInfo, pointI)
         {
             const List<pointIndexHit>& pHits = hitInfo[pointI];
-            allGeometry_[surfaces_[surfI]].getRegion(pHits, pRegions);
-            allGeometry_[surfaces_[surfI]].getNormal(pHits, pNormals);
 
-            // Extract those hits that are on higher levelled surfaces.
+            forAll(pHits, i)
-            // Note: should move extraction of region, normal outside of loop
-            // below for if getRegion/getNormal have high overhead.
-
-            forAll(pHits, pHitI)
             {
-                label region = globalRegion(surfI, pRegions[pHitI]);
+                surfInfo[n] = pHits[i];
+                pointMap[n] = pointI;
+                n++;
+            }
+        }
 
-                if (maxLevel_[region] > currentLevel[pointI])
+        labelList surfRegion(n);
-                {
+        vectorField surfNormal(n);
-                    // Append to pointI info
+        allGeometry_[surfaces_[surfI]].getRegion(surfInfo, surfRegion);
-                    label sz = surfaceNormal[pointI].size();
+        allGeometry_[surfaces_[surfI]].getNormal(surfInfo, surfNormal);
-                    surfaceNormal[pointI].setSize(sz+1);
-                    surfaceNormal[pointI][sz] = pNormals[pHitI];
 
-                    surfaceLevel[pointI].setSize(sz+1);
+        surfInfo.clear();
-                    surfaceLevel[pointI][sz] = maxLevel_[region];
+
-                }
+
+        // Extract back into pointwise
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        forAll(surfRegion, i)
+        {
+            label region = globalRegion(surfI, surfRegion[i]);
+            label pointI = pointMap[i];
+
+            if (maxLevel_[region] > currentLevel[pointI])
+            {
+                // Append to pointI info
+                label sz = surfaceNormal[pointI].size();
+                surfaceNormal[pointI].setSize(sz+1);
+                surfaceNormal[pointI][sz] = surfNormal[i];
+
+                surfaceLevel[pointI].setSize(sz+1);
+                surfaceLevel[pointI][sz] = maxLevel_[region];
             }
         }
     }

src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H

@@ -90,9 +90,6 @@ class refinementSurfaces
         //- From global region number to refinement level
         labelList maxLevel_;
 
-        //- Per surface refinement level adapted for shells.
-        PtrList<triSurfaceLabelField> minLevelFields_;
-
 
     // Private Member Functions
 
@@ -213,13 +210,6 @@ public:
                 const shellSurfaces& shells
             );
 
-            //- Helper: orient (closed only) surfaces so keepPoint is outside.
-            static void orientSurface
-            (
-                const point& keepPoint,
-                triSurfaceMesh& surface
-            );
-
             //- Helper: count number of triangles per region
             static labelList countRegions(const triSurface&);
 

src/dynamicMesh/meshCut/cellLooper/hexCellLooper.C

@@ -29,7 +29,6 @@ License
 #include "polyMesh.H"
 #include "cellModeller.H"
 #include "mathematicalConstants.H"
-#include "DynamicList.H"
 #include "plane.H"
 #include "ListOps.H"
 #include "meshTools.H"

src/dynamicMesh/meshCut/meshModifiers/meshCutAndRemove/meshCutAndRemove.H

@@ -40,7 +40,6 @@ SourceFiles
 #include "boolList.H"
 #include "labelList.H"
 #include "typeInfo.H"
-#include "DynamicList.H"
 #include "Map.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/dynamicMesh/meshCut/meshModifiers/meshCutter/meshCutter.H

@@ -113,10 +113,8 @@ SourceFiles
 #define meshCutter_H
 
 #include "edgeVertex.H"
-#include "boolList.H"
 #include "labelList.H"
 #include "typeInfo.H"
-#include "DynamicList.H"
 #include "Map.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/dynamicMesh/meshCut/meshModifiers/refinementIterator/refinementIterator.C

@@ -26,12 +26,10 @@ License
 
 #include "refinementIterator.H"
 #include "polyMesh.H"
-#include "polyTopoChanger.H"
 #include "Time.H"
 #include "refineCell.H"
 #include "undoableMeshCutter.H"
 #include "polyTopoChange.H"
-#include "DynamicList.H"
 #include "mapPolyMesh.H"
 #include "cellCuts.H"
 #include "OFstream.H"

src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C

@@ -595,6 +595,7 @@ void Foam::edgeCollapser::updateMesh(const mapPolyMesh& map)
 {
     pointRegion_.setSize(mesh_.nPoints());
     pointRegion_ = -1;
+    // Reset count, do not remove underlying storage
     pointRegionMaster_.clear();
     freeRegions_.clear();
 }

src/dynamicMesh/polyTopoChange/polyTopoChange/faceCollapser.H

@@ -54,7 +54,6 @@ SourceFiles
 #define faceCollapser_H
 
 #include "labelList.H"
-#include "DynamicList.H"
 #include "point.H"
 #include "Map.H"
 #include "labelHashSet.H"

src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C

@@ -3021,7 +3021,9 @@ Foam::labelListList Foam::hexRef8::setRefinement
             << " Checking initial mesh just to make sure" << endl;
 
         checkMesh();
-        checkRefinementLevels(-1, labelList(0));
+        // Cannot call checkRefinementlevels since hanging points might
+        // get triggered by the mesher after subsetting.
+        //checkRefinementLevels(-1, labelList(0));
     }
 
     // Clear any saved point/cell data.
@@ -4813,7 +4815,7 @@ void Foam::hexRef8::checkRefinementLevels
             false                   // no separation
         );
 
-        OFstream str(mesh_.time().path()/"hangingPoints.obj");
+        //OFstream str(mesh_.time().path()/"hangingPoints.obj");
 
         label nHanging = 0;
 
@@ -4826,8 +4828,7 @@ void Foam::hexRef8::checkRefinementLevels
                 Pout<< "Hanging boundary point " << pointI
                     << " at " << mesh_.points()[pointI]
                     << endl;
-
+                //meshTools::writeOBJ(str, mesh_.points()[pointI]);
-                meshTools::writeOBJ(str, mesh_.points()[pointI]);
             }
         }
 

src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChangeTemplates.C

@@ -96,13 +96,5 @@ void Foam::polyTopoChange::renumberKey
     elems.transfer(newElems);
 }
 
-// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
-
 
 // ************************************************************************* //

src/edgeMesh/edgeMesh.C

@@ -156,7 +156,7 @@ label edgeMesh::regions(labelList& edgeRegion) const
                 }
             }
 
-            edgesToVisit.transfer(newEdgesToVisit.shrink());
+            edgesToVisit.transfer(newEdgesToVisit);
         }
 
         currentRegion++;

src/finiteVolume/finiteVolume/ddtSchemes/CrankNicholsonDdtScheme/CrankNicholsonDdtScheme.C

@@ -1027,8 +1027,7 @@ CrankNicholsonDdtScheme<Type>::fvcDdtPhiCorr
             {
                 dUdt0 = rDtCoef0_(dUdt0)*
                 (
-                    rho.oldTime()*U.oldTime()
+                    U.oldTime() - U.oldTime().oldTime()
-                  - rho.oldTime().oldTime()*U.oldTime().oldTime()
                 ) - offCentre_(dUdt0());
             }
 
@@ -1054,11 +1053,8 @@ CrankNicholsonDdtScheme<Type>::fvcDdtPhiCorr
                       - (
                             fvc::interpolate
                             (
-                                rA*
+                                rA*rho.oldTime()
-                                (
+                               *(rDtCoef*U.oldTime() + offCentre_(dUdt0()))
-                                    rDtCoef*rho.oldTime()*U.oldTime()
-                                    + offCentre_(dUdt0())
-                                )
                             ) & mesh().Sf()
                         )
                     )

src/meshTools/Make/files

@@ -54,12 +54,15 @@ indexedOctree/treeDataPoint.C
 indexedOctree/treeDataTriSurface.C
 
 searchableSurface = searchableSurface
+$(searchableSurface)/distributedTriSurfaceMesh.C
 $(searchableSurface)/searchableBox.C
+$(searchableSurface)/searchablePlate.C
 $(searchableSurface)/searchableSphere.C
 $(searchableSurface)/searchableSurface.C
 $(searchableSurface)/searchableSurfaces.C
 $(searchableSurface)/searchableSurfacesQueries.C
 $(searchableSurface)/triSurfaceMesh.C
+$(searchableSurface)/searchableSurfaceWithGaps.C
 
 topoSets = sets/topoSets
 $(topoSets)/cellSet.C

src/meshTools/cellFeatures/cellFeatures.C

@@ -371,7 +371,7 @@ void Foam::cellFeatures::calcSuperFaces() const
                 {
                     superFace.shrink();
 
-                    faces[superFaceI] = face(superFace);
+                    faces[superFaceI].transfer(superFace);
                 }
             }
         }

src/meshTools/meshSearch/meshSearch.C

@@ -921,14 +921,7 @@ Foam::List<Foam::pointIndexHit> Foam::meshSearch::intersections
 
     hits.shrink();
 
-    // Copy into straight list
+    return hits;
-    List<pointIndexHit> allHits(hits.size());
-
-    forAll(hits, hitI)
-    {
-        allHits[hitI] = hits[hitI];
-    }
-    return allHits;
 }
 
 

src/meshTools/octree/pointIndexHit.H

@@ -45,7 +45,7 @@ typedef PointIndexHit<point> pointIndexHit;
 
 //- Specify data associated with pointIndexHit type is contiguous
 template<>
-inline bool contiguous<pointIndexHit>() {return true;}
+inline bool contiguous<pointIndexHit>() {return contiguous<point>();}
 
 }
 

src/meshTools/octree/treeBoundBox.C

@@ -679,7 +679,11 @@ Foam::label Foam::treeBoundBox::distanceCmp
 
 bool Foam::operator==(const treeBoundBox& a, const treeBoundBox& b)
 {
-    return (a.min() == b.min()) && (a.max() == b.max());
+    return operator==
+    (
+        static_cast<const boundBox&>(a),
+        static_cast<const boundBox&>(b)
+    );
 }
 
 
@@ -689,12 +693,17 @@ bool Foam::operator!=(const treeBoundBox& a, const treeBoundBox& b)
 }
 
 
-// * * * * * * * * * * * * * * * IOstream Operator  * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<(Ostream& os, const treeBoundBox& bb)
+{
+    return os << static_cast<const boundBox&>(bb);
+}
+
 
 Foam::Istream& Foam::operator>>(Istream& is, treeBoundBox& bb)
 {
-    is >> bb.min() >> bb.max();
+    return is >> static_cast<boundBox&>(bb);
-    return is;
 }
 
 

src/meshTools/octree/treeBoundBox.H

@@ -322,13 +322,19 @@ public:
         friend bool operator==(const treeBoundBox&, const treeBoundBox&);
         friend bool operator!=(const treeBoundBox&, const treeBoundBox&);
 
-
     // IOstream operator
 
-        friend Istream& operator>>(Istream&, treeBoundBox&);
+        friend Istream& operator>>(Istream& is, treeBoundBox&);
+        friend Ostream& operator<<(Ostream& os, const treeBoundBox&);
+
 };
 
 
+//- Specify data associated with boundBox type is contiguous
+template<>
+inline bool contiguous<treeBoundBox>() {return contiguous<boundBox>();}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/meshTools/searchableSurface/distributedTriSurfaceMesh.H

@@ -0,0 +1,459 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::distributedTriSurfaceMesh
+
+Description
+    IOoject and searching on triSurface
+
+SourceFiles
+    distributedTriSurfaceMesh.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef distributedTriSurfaceMesh_H
+#define distributedTriSurfaceMesh_H
+
+#include "triSurfaceMesh.H"
+#include "IOdictionary.H"
+#include "Pair.H"
+#include "globalIndex.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class mapDistribute;
+
+// Typedefs
+typedef Pair<point> segment;
+template<>
+inline bool contiguous<segment>() {return contiguous<point>();}
+
+
+/*---------------------------------------------------------------------------*\
+                           Class distributedTriSurfaceMesh Declaration
+\*---------------------------------------------------------------------------*/
+
+class distributedTriSurfaceMesh
+:
+    public triSurfaceMesh
+{
+private:
+
+    // Static data
+
+        //- Merging distance
+        static scalar mergeDist_;
+
+
+    // Private member data
+
+        //- Bounding box settings
+        IOdictionary dict_;
+
+        //- bounding boxes of all processors
+        List<List<treeBoundBox> > procBb_;
+
+        //- Global triangle numbering
+        mutable autoPtr<globalIndex> globalTris_;
+
+
+    // Private Member Functions
+
+        // Read
+
+            //- Read my additional data
+            bool read();
+
+
+        // Line intersection
+
+            static bool isLocal
+            (
+                const List<treeBoundBox>& myBbs,
+                const point& start,
+                const point& end
+            );
+
+            //- Split segment into subsegments overlapping the processor
+            //  bounding boxes. Sort per processor.
+            void splitSegment
+            (
+                const label,
+                const point& start,
+                const point& end,
+
+                DynamicList<segment>&,
+                DynamicList<label>&,
+                List<DynamicList<label> >&
+            ) const;
+
+            //- Divide edges into local and remote segments. Construct map to
+            //  distribute and collect data.
+            autoPtr<mapDistribute> constructSegments
+            (
+                const pointField& start,
+                const pointField& end,
+
+                List<segment>& allSegments,
+                List<label>& allSegmentMap
+            ) const;
+
+            //- Split edges, distribute, test and collect.
+            void findLine
+            (
+                const bool nearestIntersection,
+                const pointField& start,
+                const pointField& end,
+                List<pointIndexHit>& info
+            ) const;
+
+
+        // Triangle index
+
+            //- Obtains global indices from pointIndexHit and swaps them back
+            //  to their original processor. Used to calculate local region
+            //  and normal.
+            autoPtr<mapDistribute> calcLocalQueries
+            (
+                const List<pointIndexHit>&,
+                labelList& triangleIndex
+            ) const;
+
+
+        // Nearest
+
+            label calcOverlappingProcs
+            (
+                const point& centre,
+                const scalar radiusSqr,
+                boolList& overlaps
+            ) const;
+
+            autoPtr<mapDistribute> calcLocalQueries
+            (
+                const pointField& centres,
+                const scalarField& radiusSqr,
+
+                pointField& allCentres,
+                scalarField& allRadiusSqr,
+                labelList& allSegmentMap
+            ) const;
+
+
+        // Surface redistribution
+
+            //- Calculate surface bounding box
+            void calcBounds(boundBox& bb, label& nPoints) const;
+
+            //- Does any part of triangle overlap bb.
+            static bool overlaps
+            (
+                const List<treeBoundBox>& bb,
+                const point& p0,
+                const point& p1,
+                const point& p2
+            );
+
+            //- Find points used in subset
+            static void subsetMeshMap
+            (
+                const triSurface& s,
+                const boolList& include,
+                const label nIncluded,
+                labelList& newToOldPoints,
+                labelList& oldToNewPoints,
+                labelList& newToOldFaces
+            );
+
+            //- Construct subsetted surface
+            static triSurface subsetMesh
+            (
+                const triSurface& s,
+                const labelList& newToOldPoints,
+                const labelList& oldToNewPoints,
+                const labelList& newToOldFaces
+            );
+
+            //- Subset given marked faces
+            static triSurface subsetMesh
+            (
+                const triSurface& s,
+                const boolList& include,
+                labelList& newToOldPoints,
+                labelList& newToOldFaces
+            );
+
+            //- Subset given marked faces
+            static triSurface subsetMesh
+            (
+                const triSurface& s,
+                const labelList& newToOldFaces,
+                labelList& newToOldPoints
+            );
+
+            //- Find triangle otherF in allFaces.
+            static label findTriangle
+            (
+                const List<labelledTri>& allFaces,
+                const labelListList& allPointFaces,
+                const labelledTri& otherF
+            );
+
+            //- Merge triSurface (subTris, subPoints) into allTris, allPoints.
+            static void merge
+            (
+                const scalar mergeDist,
+                const List<labelledTri>& subTris,
+                const pointField& subPoints,
+
+                List<labelledTri>& allTris,
+                pointField& allPoints,
+
+                labelList& faceConstructMap,
+                labelList& pointConstructMap
+            );
+
+            //- Distribute stored fields
+            template<class Type>
+            void distributeFields(const mapDistribute& map);
+
+
+        //- Disallow default bitwise copy construct
+        distributedTriSurfaceMesh(const distributedTriSurfaceMesh&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const distributedTriSurfaceMesh&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("distributedTriSurfaceMesh");
+
+
+    // Constructors
+
+        //- Construct from triSurface
+        distributedTriSurfaceMesh
+        (
+            const IOobject&,
+            const triSurface&,
+            const dictionary& dict
+        );
+
+        //- Construct read
+        distributedTriSurfaceMesh(const IOobject& io);
+
+        //- Construct from dictionary (used by searchableSurface)
+        distributedTriSurfaceMesh
+        (
+            const IOobject& io,
+            const dictionary& dict
+        );
+
+
+    // Destructor
+
+        virtual ~distributedTriSurfaceMesh();
+
+        //- Clear storage
+        void clearOut();
+
+
+    // Member Functions
+
+        //- Return merge distance
+        static scalar mergeDist()
+        {
+            return mergeDist_;
+        }
+
+        //- Set the merge distance, returning the previous value
+        static scalar setMergeDist(const scalar t)
+        {
+            if (t < -VSMALL)
+            {
+                FatalErrorIn
+                (
+                    "scalar distributedTriSurfaceMesh::setMergeDist"
+                    "(const scalar t)"
+                )   << "Negative merge distance tolerance. This is not allowed."
+                    << abort(FatalError);
+            }
+
+            scalar oldTol = mergeDist_;
+            mergeDist_ = t;
+
+            return oldTol;
+        }
+
+
+        //- Triangle indexing (demand driven)
+        const globalIndex& globalTris() const;
+
+
+        // searchableSurface implementation
+
+            //- Whether supports volume type below. I.e. whether is closed.
+            //  Not supported.
+            virtual bool hasVolumeType() const
+            {
+                return false;
+            }
+
+            //- Range of global indices that can be returned.
+            virtual label globalSize() const
+            {
+                const labelList& offsets = globalTris().offsets();
+                return offsets[offsets.size()-1];
+            }
+
+            virtual void findNearest
+            (
+                const pointField& sample,
+                const scalarField& nearestDistSqr,
+                List<pointIndexHit>&
+            ) const;
+
+            virtual void findLine
+            (
+                const pointField& start,
+                const pointField& end,
+                List<pointIndexHit>&
+            ) const;
+
+            virtual void findLineAny
+            (
+                const pointField& start,
+                const pointField& end,
+                List<pointIndexHit>&
+            ) const;
+
+            //- Get all intersections in order from start to end.
+            virtual void findLineAll
+            (
+                const pointField& start,
+                const pointField& end,
+                List<List<pointIndexHit> >&
+            ) const;
+
+            //- From a set of points and indices get the region
+            virtual void getRegion
+            (
+                const List<pointIndexHit>&,
+                labelList& region
+            ) const;
+
+            //- From a set of points and indices get the normal
+            virtual void getNormal
+            (
+                const List<pointIndexHit>&,
+                vectorField& normal
+            ) const;
+
+            //- Determine type (inside/outside/mixed) for point. unknown if
+            //  cannot be determined (e.g. non-manifold surface)
+            virtual void getVolumeType
+            (
+                const pointField&,
+                List<volumeType>&
+            ) const;
+
+            //- Set bounds of surface. Bounds currently set as list of
+            //  bounding boxes. Will do redistribution of surface to locally
+            //  have all triangles overlapping bounds.
+            //  Larger bounds: more triangles (memory), more fully local tests
+            //  (quick).
+            //  keepNonLocal = true : keep triangles that do not overlap
+            //  any processor bounds.
+            //  Should really be split into a routine to determine decomposition
+            //  and one that does actual distribution but determining
+            //  decomposition with duplicate triangle merging requires
+            //  same amoun as work as actual distribution.
+            virtual void distribute
+            (
+                const List<treeBoundBox>&,
+                const bool keepNonLocal,
+                autoPtr<mapDistribute>& faceMap,
+                autoPtr<mapDistribute>& pointMap
+            );
+
+
+        // Other
+
+            //- Specific to triSurfaceMesh: from a set of hits (points and
+            //  indices) get the specified field. Misses do not get set.
+            virtual void getField
+            (
+                const word& fieldName,
+                const List<pointIndexHit>&,
+                labelList& values
+            ) const;
+
+            //- Subset the part of surface that is overlapping bounds.
+            static triSurface overlappingSurface
+            (
+                const triSurface&,
+                const List<treeBoundBox>&,
+                labelList& subPointMap,
+                labelList& subFaceMap
+            );
+
+            //- Print some stats. Parallel aware version of
+            //  triSurface::writeStats.
+            void writeStats(Ostream& os) const;
+
+
+        // regIOobject implementation
+
+            //- Write using given format, version and compression
+            virtual bool writeObject
+            (
+                IOstream::streamFormat fmt,
+                IOstream::versionNumber ver,
+                IOstream::compressionType cmp
+            ) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "distributedTriSurfaceMeshTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/meshTools/searchableSurface/distributedTriSurfaceMeshTemplates.C

@@ -0,0 +1,139 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "distributedTriSurfaceMesh.H"
+#include "triSurfaceFields.H"
+#include "mapDistribute.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+//template<class Type>
+//void Foam::distributedTriSurfaceMesh::getField
+//(
+//    const word& fieldName,
+//    const List<pointIndexHit>& info,
+//    List<Type>& values
+//) const
+//{
+//    typedef DimensionedField<Type, triSurfaceGeoMesh> DimensionedSurfField;
+//
+//
+//    // Get query data (= local index of triangle)
+//    // ~~~~~~~~~~~~~~
+//
+//    labelList triangleIndex(info.size());
+//    autoPtr<mapDistribute> mapPtr
+//    (
+//        calcLocalQueries
+//        (
+//            info,
+//            triangleIndex
+//        )
+//    );
+//    const mapDistribute& map = mapPtr();
+//
+//
+//    // Do my tests
+//    // ~~~~~~~~~~~
+//
+//    const DimensionedSurfField& fld = lookupObject<DimensionedSurfField>
+//    (
+//        fieldName
+//    );
+//    const triSurface& s = static_cast<const triSurface&>(*this);
+//
+//    values.setSize(triangleIndex.size());
+//
+//    forAll(triangleIndex, i)
+//    {
+//        label triI = triangleIndex[i];
+//        values[i] = fld[triI];
+//    }
+//
+//
+//    // Send back results
+//    // ~~~~~~~~~~~~~~~~~
+//
+//    map.distribute
+//    (
+//        Pstream::nonBlocking,
+//        List<labelPair>(0),
+//        info.size(),
+//        map.constructMap(),     // what to send
+//        map.subMap(),           // what to receive
+//        values
+//    );
+//}
+
+
+template<class Type>
+void Foam::distributedTriSurfaceMesh::distributeFields
+(
+    const mapDistribute& map
+)
+{
+    typedef DimensionedField<Type, triSurfaceGeoMesh> DimensionedSurfField;
+
+    HashTable<const DimensionedSurfField*> fields
+    (
+        objectRegistry::lookupClass
+        <DimensionedSurfField >()
+    );
+
+    for
+    (
+        typename HashTable<const DimensionedSurfField*>::iterator fieldIter =
+            fields.begin();
+        fieldIter != fields.end();
+        ++fieldIter
+    )
+    {
+        DimensionedSurfField& field =
+            const_cast<DimensionedSurfField&>(*fieldIter());
+
+        label oldSize = field.size();
+
+        map.distribute
+        (
+            Pstream::nonBlocking,
+            List<labelPair>(0),
+            map.constructSize(),
+            map.subMap(),
+            map.constructMap(),
+            field
+        );
+
+        if (debug)
+        {
+            Info<< "Mapped " << field.typeName << ' ' << field.name()
+                << " from size " << oldSize << " to size " << field.size()
+                << endl;
+        }
+    }
+}
+
+
+// ************************************************************************* //

src/meshTools/searchableSurface/searchableBox.C

@@ -484,7 +484,6 @@ void Foam::searchableBox::findLineAll
 
             hits.shrink();
             info[pointI].transfer(hits);
-            hits.clear();
         }
         else
         {

src/meshTools/searchableSurface/searchableBox.H

@@ -121,6 +121,12 @@ public:
             return true;
         }
 
+        //- Range of local indices that can be returned.
+        virtual label size() const
+        {
+            return 6;
+        }
+
 
         // Single point queries.
 

src/meshTools/searchableSurface/searchablePlate.C

@@ -0,0 +1,366 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "searchablePlate.H"
+#include "addToRunTimeSelectionTable.H"
+#include "SortableList.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(searchablePlate, 0);
+addToRunTimeSelectionTable(searchableSurface, searchablePlate, dict);
+
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::direction Foam::searchablePlate::calcNormal(const point& span)
+{
+    direction normalDir = 3;
+
+    for (direction dir = 0; dir < vector::nComponents; dir++)
+    {
+        if (span[dir] < 0)
+        {
+            FatalErrorIn("searchablePlate::calcNormal()")
+                << "Span should have two positive and one zero entry. Now:"
+                << span << exit(FatalError);
+        }
+        else if (span[dir] < VSMALL)
+        {
+            if (normalDir == 3)
+            {
+                normalDir = dir;
+            }
+            else
+            {
+                // Multiple zero entries. Flag and exit.
+                normalDir = 3;
+                break;
+            }
+        }
+    }
+
+    if (normalDir == 3)
+    {
+        FatalErrorIn("searchablePlate::calcNormal()")
+            << "Span should have one and only zero entry. Now:" << span
+            << exit(FatalError);
+    }
+
+    return normalDir;
+}
+
+
+// Returns miss or hit with face (always 0)
+Foam::pointIndexHit Foam::searchablePlate::findNearest
+(
+    const point& sample,
+    const scalar nearestDistSqr
+) const
+{
+    // For every component direction can be
+    // left of min, right of max or inbetween.
+    // - outside points: project first one x plane (either min().x()
+    // or max().x()), then onto y plane and finally z. You should be left
+    // with intersection point
+    // - inside point: find nearest side (compare to mid point). Project onto
+    //   that.
+
+    // Project point on plane.
+    pointIndexHit info(true, sample, 0);
+    info.rawPoint()[normalDir_] = origin_[normalDir_];
+
+    // Clip to edges if outside
+    for (direction dir = 0; dir < vector::nComponents; dir++)
+    {
+        if (dir != normalDir_)
+        {
+            if (info.rawPoint()[dir] < origin_[dir])
+            {
+                info.rawPoint()[dir] = origin_[dir];
+            }
+            else if (info.rawPoint()[dir] > origin_[dir]+span_[dir])
+            {
+                info.rawPoint()[dir] = origin_[dir]+span_[dir];
+            }
+        }
+    }
+
+    // Check if outside. Optimisation: could do some checks on distance already
+    // on components above
+    if (magSqr(info.rawPoint() - sample) > nearestDistSqr)
+    {
+        info.setMiss();
+        info.setIndex(-1);
+    }
+
+    return info;
+}
+
+
+Foam::pointIndexHit Foam::searchablePlate::findLine
+(
+    const point& start,
+    const point& end
+) const
+{
+    pointIndexHit info
+    (
+        true,
+        vector::zero,
+        0
+    );
+
+    const vector dir(end-start);
+
+    if (mag(dir[normalDir_]) < VSMALL)
+    {
+        info.setMiss();
+        info.setIndex(-1);
+    }
+    else
+    {
+        scalar t = (origin_[normalDir_]-start[normalDir_]) / dir[normalDir_];
+
+        if (t < 0 || t > 1)
+        {
+            info.setMiss();
+            info.setIndex(-1);
+        }
+        else
+        {
+            info.rawPoint() = start+t*dir;
+            info.rawPoint()[normalDir_] = origin_[normalDir_];
+
+            // Clip to edges
+            for (direction dir = 0; dir < vector::nComponents; dir++)
+            {
+                if (dir != normalDir_)
+                {
+                    if (info.rawPoint()[dir] < origin_[dir])
+                    {
+                        info.setMiss();
+                        info.setIndex(-1);
+                        break;
+                    }
+                    else if (info.rawPoint()[dir] > origin_[dir]+span_[dir])
+                    {
+                        info.setMiss();
+                        info.setIndex(-1);
+                        break;
+                    }
+                }
+            }
+        }
+    }
+
+    // Debug
+    if (info.hit())
+    {
+        treeBoundBox bb(origin_, origin_+span_);
+        bb.min()[normalDir_] -= 1E-6;
+        bb.max()[normalDir_] += 1E-6;
+
+        if (!bb.contains(info.hitPoint()))
+        {
+            FatalErrorIn("searchablePlate::findLine(..)")
+                << "bb:" << bb << endl
+                << "origin_:" << origin_ << endl
+                << "span_:" << span_ << endl
+                << "normalDir_:" << normalDir_ << endl
+                << "hitPoint:" << info.hitPoint()
+                << abort(FatalError);
+        }
+    }
+
+    return info;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::searchablePlate::searchablePlate
+(
+    const IOobject& io,
+    const point& origin,
+    const vector& span
+)
+:
+    searchableSurface(io),
+    origin_(origin),
+    span_(span),
+    normalDir_(calcNormal(span_))
+{}
+
+
+Foam::searchablePlate::searchablePlate
+(
+    const IOobject& io,
+    const dictionary& dict
+)
+:
+    searchableSurface(io),
+    origin_(dict.lookup("origin")),
+    span_(dict.lookup("span")),
+    normalDir_(calcNormal(span_))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::searchablePlate::~searchablePlate()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::wordList& Foam::searchablePlate::regions() const
+{
+    if (regions_.size() == 0)
+    {
+        regions_.setSize(1);
+        regions_[0] = "region0";
+    }
+    return regions_;
+}
+
+
+void Foam::searchablePlate::findNearest
+(
+    const pointField& samples,
+    const scalarField& nearestDistSqr,
+    List<pointIndexHit>& info
+) const
+{
+    info.setSize(samples.size());
+
+    forAll(samples, i)
+    {
+        info[i] = findNearest(samples[i], nearestDistSqr[i]);
+    }
+}
+
+
+void Foam::searchablePlate::findLine
+(
+    const pointField& start,
+    const pointField& end,
+    List<pointIndexHit>& info
+) const
+{
+    info.setSize(start.size());
+
+    forAll(start, i)
+    {
+        info[i] = findLine(start[i], end[i]);
+    }
+}
+
+
+void Foam::searchablePlate::findLineAny
+(
+    const pointField& start,
+    const pointField& end,
+    List<pointIndexHit>& info
+) const
+{
+    findLine(start, end, info);
+}
+
+
+void Foam::searchablePlate::findLineAll
+(
+    const pointField& start,
+    const pointField& end,
+    List<List<pointIndexHit> >& info
+) const
+{
+    List<pointIndexHit> nearestInfo;
+    findLine(start, end, nearestInfo);
+
+    info.setSize(start.size());
+    forAll(info, pointI)
+    {
+        if (nearestInfo[pointI].hit())
+        {
+            info[pointI].setSize(1);
+            info[pointI][0] = nearestInfo[pointI];
+        }
+        else
+        {
+            info[pointI].clear();
+        }
+    }
+}
+
+
+void Foam::searchablePlate::getRegion
+(
+    const List<pointIndexHit>& info,
+    labelList& region
+) const
+{
+    region.setSize(info.size());
+    region = 0;
+}
+
+
+void Foam::searchablePlate::getNormal
+(
+    const List<pointIndexHit>& info,
+    vectorField& normal
+) const
+{
+    normal.setSize(info.size());
+    normal = vector::zero;
+    forAll(normal, i)
+    {
+        normal[i][normalDir_] = 1.0;
+    }
+}
+
+
+void Foam::searchablePlate::getVolumeType
+(
+    const pointField& points,
+    List<volumeType>& volType
+) const
+{
+    FatalErrorIn
+    (
+        "searchableCollection::getVolumeType(const pointField&"
+        ", List<volumeType>&) const"
+    )   << "Volume type not supported for plate."
+        << exit(FatalError);
+}
+
+
+// ************************************************************************* //

src/meshTools/searchableSurface/searchablePlate.H

@@ -0,0 +1,215 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::searchablePlate
+
+Description
+    Searching on finite plate. Plate has to be aligned with coordinate
+    axes!
+
+SourceFiles
+    searchablePlate.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef searchablePlate_H
+#define searchablePlate_H
+
+#include "searchableSurface.H"
+#include "treeBoundBox.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+
+/*---------------------------------------------------------------------------*\
+                           Class searchablePlate Declaration
+\*---------------------------------------------------------------------------*/
+
+class searchablePlate
+:
+    public searchableSurface
+{
+private:
+
+    // Private Member Data
+
+        const point origin_;
+
+        const vector span_;
+
+        //- Coordinate direction which is normal
+        const direction normalDir_;
+
+        mutable wordList regions_;
+
+
+    // Private Member Functions
+
+        //- Calculate normal direction from span
+        static direction calcNormal(const point&);
+
+        pointIndexHit findNearest
+        (
+            const point& sample,
+            const scalar nearestDistSqr
+        ) const;
+
+        pointIndexHit findLine
+        (
+            const point& start,
+            const point& end
+        ) const;
+
+
+        //- Disallow default bitwise copy construct
+        searchablePlate(const searchablePlate&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const searchablePlate&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("searchablePlate");
+
+
+    // Constructors
+
+        //- Construct from components
+        searchablePlate
+        (
+            const IOobject& io,
+            const point& origin,
+            const point& span
+        );
+
+        //- Construct from dictionary (used by searchableSurface)
+        searchablePlate
+        (
+            const IOobject& io,
+            const dictionary& dict
+        );
+
+    // Destructor
+
+        virtual ~searchablePlate();
+
+
+    // Member Functions
+
+        virtual const wordList& regions() const;
+
+        //- Whether supports volume type below
+        virtual bool hasVolumeType() const
+        {
+            return false;
+        }
+
+        //- Range of local indices that can be returned.
+        virtual label size() const
+        {
+            return 1;
+        }
+
+
+        // Multiple point queries.
+
+            virtual void findNearest
+            (
+                const pointField& sample,
+                const scalarField& nearestDistSqr,
+                List<pointIndexHit>&
+            ) const;
+
+            virtual void findLine
+            (
+                const pointField& start,
+                const pointField& end,
+                List<pointIndexHit>&
+            ) const;
+
+            virtual void findLineAny
+            (
+                const pointField& start,
+                const pointField& end,
+                List<pointIndexHit>&
+            ) const;
+
+            //- Get all intersections in order from start to end.
+            virtual void findLineAll
+            (
+                const pointField& start,
+                const pointField& end,
+                List<List<pointIndexHit> >&
+            ) const;
+
+            //- From a set of points and indices get the region
+            virtual void getRegion
+            (
+                const List<pointIndexHit>&,
+                labelList& region
+            ) const;
+
+            //- From a set of points and indices get the normal
+            virtual void getNormal
+            (
+                const List<pointIndexHit>&,
+                vectorField& normal
+            ) const;
+
+            //- Determine type (inside/outside/mixed) for point. unknown if
+            //  cannot be determined (e.g. non-manifold surface)
+            virtual void getVolumeType
+            (
+                const pointField&,
+                List<volumeType>&
+            ) const;
+
+
+        // regIOobject implementation
+
+            bool writeData(Ostream&) const
+            {
+                notImplemented("searchablePlate::writeData(Ostream&) const");
+                return false;
+            }
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/meshTools/searchableSurface/searchableSphere.H

@@ -36,8 +36,8 @@ SourceFiles
 #ifndef searchableSphere_H
 #define searchableSphere_H
 
-#include "searchableSurface.H"
 #include "treeBoundBox.H"
+#include "searchableSurface.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -127,6 +127,12 @@ public:
             return true;
         }
 
+        //- Range of local indices that can be returned.
+        virtual label size() const
+        {
+            return 1;
+        }
+
 
         // Multiple point queries.
 

src/meshTools/searchableSurface/searchableSurface.H

@@ -56,6 +56,7 @@ namespace Foam
 
 // Forward declaration of classes
 class objectRegistry;
+class mapDistribute;
 class treeBoundBox;
 
 /*---------------------------------------------------------------------------*\
@@ -180,6 +181,15 @@ public:
         //- Whether supports volume type below.
         virtual bool hasVolumeType() const = 0;
 
+        //- Range of local indices that can be returned.
+        virtual label size() const = 0;
+
+        //- Range of global indices that can be returned.
+        virtual label globalSize() const
+        {
+            return size();
+        }
+
 
         // Single point queries.
 
@@ -296,16 +306,18 @@ public:
 
         // Other
 
-            ////- Get bounding box.
+            //- Set bounds of surface. Bounds currently set as list of
-            //const boundBox& bounds() const = 0;
+            //  bounding boxes. The bounds are hints to the surface as for
-
+            //  the range of queries it can expect. faceMap/pointMap can be
-            ////- Set bounding box.
+            //  set if the surface has done any redistribution.
-            //void setBounds
+            virtual void distribute
-            //(
+            (
-            //    const boundBox&,
+                const List<treeBoundBox>&,
-            //    autoPtr<mapDistribute>& faceMap,
+                const bool keepNonLocal,
-            //    autoPtr<mapDistribute>& pointMap
+                autoPtr<mapDistribute>& faceMap,
-            //) = 0;
+                autoPtr<mapDistribute>& pointMap
+            )
+            {}
 
 
 };

src/meshTools/searchableSurface/searchableSurfaceWithGaps.C

@@ -0,0 +1,321 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "searchableSurfaceWithGaps.H"
+#include "addToRunTimeSelectionTable.H"
+#include "SortableList.H"
+#include "Time.H"
+#include "ListOps.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(searchableSurfaceWithGaps, 0);
+addToRunTimeSelectionTable(searchableSurface, searchableSurfaceWithGaps, dict);
+
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::Pair<Foam::vector> Foam::searchableSurfaceWithGaps::offsetVecs
+(
+    const point& start,
+    const point& end
+) const
+{
+    Pair<vector> offsets(vector::zero, vector::zero);
+
+    vector n(end-start);
+
+    scalar magN = mag(n);
+
+    if (magN > SMALL)
+    {
+        n /= magN;
+
+        // Do first offset vector. Is the coordinate axes with the smallest
+        // component along the vector n.
+        scalar minMag = GREAT;
+        direction minCmpt = 0;
+
+        for (direction cmpt = 0; cmpt < vector::nComponents; cmpt++)
+        {
+            if (mag(n[cmpt]) < minMag)
+            {
+                minMag = mag(n[cmpt]);
+                minCmpt = cmpt;
+            }
+        }
+
+        offsets[0][minCmpt] = 1.0;
+        // Orthogonalise
+        offsets[0] -= n[minCmpt]*n;
+        // Scale
+        offsets[0] *= gap_/mag(offsets[0]);
+
+
+        // Do second offset vector perp to original edge and first offset vector
+        offsets[1] = n ^ offsets[0];
+        offsets[1] *= gap_/mag(offsets[1]);
+    }
+
+    return offsets;    
+}
+
+
+void Foam::searchableSurfaceWithGaps::offsetVecs
+(
+    const pointField& start,
+    const pointField& end,
+    pointField& offset0,
+    pointField& offset1
+) const
+{
+    offset0.setSize(start.size());
+    offset1.setSize(start.size());
+
+    forAll(start, i)
+    {
+        const Pair<vector> offsets(offsetVecs(start[i], end[i]));
+        offset0[i] = offsets[0];
+        offset1[i] = offsets[1];
+    }
+}
+
+
+Foam::label Foam::searchableSurfaceWithGaps::countMisses
+(
+    const List<pointIndexHit>& info,
+    labelList& missMap
+)
+{
+    label nMiss = 0;
+    forAll(info, i)
+    {
+        if (!info[i].hit())
+        {
+            nMiss++;
+        }
+    }
+
+    missMap.setSize(nMiss);
+    nMiss = 0;
+
+    forAll(info, i)
+    {
+        if (!info[i].hit())
+        {
+            missMap[nMiss++] = i;
+        }
+    }
+
+    return nMiss;
+}
+
+
+// Anything not a hit in both counts as a hit
+Foam::label Foam::searchableSurfaceWithGaps::countMisses
+(
+    const List<pointIndexHit>& plusInfo,
+    const List<pointIndexHit>& minInfo,
+    labelList& missMap
+)
+{
+    label nMiss = 0;
+    forAll(plusInfo, i)
+    {
+        if (!plusInfo[i].hit() || !minInfo[i].hit())
+        {
+            nMiss++;
+        }
+    }
+
+    missMap.setSize(nMiss);
+    nMiss = 0;
+
+    forAll(plusInfo, i)
+    {
+        if (!plusInfo[i].hit() || !minInfo[i].hit())
+        {
+            missMap[nMiss++] = i;
+        }
+    }
+
+    return nMiss;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::searchableSurfaceWithGaps::searchableSurfaceWithGaps
+(
+    const IOobject& io,
+    const dictionary& dict
+)
+:
+    searchableSurface(io),
+    gap_(readScalar(dict.lookup("gap"))),
+    subGeom_(1)
+{
+    const word subGeomName(dict.lookup("surface"));
+
+    const searchableSurface& s =
+        io.db().lookupObject<searchableSurface>(subGeomName);
+
+    subGeom_.set(0, &const_cast<searchableSurface&>(s));
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::searchableSurfaceWithGaps::~searchableSurfaceWithGaps()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::searchableSurfaceWithGaps::findLine
+(
+    const pointField& start,
+    const pointField& end,
+    List<pointIndexHit>& info
+) const
+{
+    surface().findLine(start, end, info);
+
+    // Count number of misses. Determine map
+    labelList compactMap;
+    label nMiss = countMisses(info, compactMap);
+
+    if (returnReduce(nMiss, sumOp<label>()) > 0)
+    {
+        pointField compactStart(start, compactMap);
+        pointField compactEnd(end, compactMap);
+
+        // Calculate offset vector
+        pointField offset0, offset1;
+        offsetVecs
+        (
+            compactStart,
+            compactEnd,
+            offset0,
+            offset1
+        );
+
+        // test in pairs: only if both perturbations hit something
+        // do we accept the hit.
+
+        List<pointIndexHit> plusInfo;
+        surface().findLine(start+offset0, end+offset0, plusInfo);
+        List<pointIndexHit> minInfo;
+        surface().findLine(start-offset0, end-offset0, minInfo);
+
+        // Extract any hits
+        forAll(plusInfo, i)
+        {
+            if (plusInfo[i].hit() && minInfo[i].hit())
+            {
+                info[compactMap[i]] = plusInfo[i].hitPoint()-offset0[i];
+            }
+        }
+
+        labelList plusMissMap;
+        nMiss = countMisses(plusInfo, minInfo, plusMissMap);
+
+        if (returnReduce(nMiss, sumOp<label>()) > 0)
+        {
+            // Extract (inplace possible because of order)
+            forAll(plusMissMap, i)
+            {
+                label mapI = plusMissMap[i];
+                compactStart[i] = compactStart[mapI];
+                compactEnd[i] = compactEnd[mapI];
+                compactMap[i] = compactMap[mapI];
+            }
+            compactStart.setSize(plusMissMap.size());
+            compactEnd.setSize(plusMissMap.size());
+            compactMap.setSize(plusMissMap.size());
+
+
+            surface().findLine(start+offset1, end+offset1, plusInfo);
+            surface().findLine(start-offset1, end-offset1, minInfo);
+
+            // Extract any hits
+            forAll(plusInfo, i)
+            {
+                if (plusInfo[i].hit() && minInfo[i].hit())
+                {
+                    info[compactMap[i]] = plusInfo[i].hitPoint()-offset1[i];
+                }
+            }
+        }
+    }
+}
+
+
+void Foam::searchableSurfaceWithGaps::findLineAny
+(
+    const pointField& start,
+    const pointField& end,
+    List<pointIndexHit>& info
+) const
+{
+    // To be done ...
+    findLine(start, end, info);
+}
+
+
+void Foam::searchableSurfaceWithGaps::findLineAll
+(
+    const pointField& start,
+    const pointField& end,
+    List<List<pointIndexHit> >& info
+) const
+{
+    // To be done. Assume for now only one intersection.
+    List<pointIndexHit> nearestInfo;
+    findLine(start, end, nearestInfo);
+
+    info.setSize(start.size());
+    forAll(info, pointI)
+    {
+        if (nearestInfo[pointI].hit())
+        {
+            info[pointI].setSize(1);
+            info[pointI][0] = nearestInfo[pointI];
+        }
+        else
+        {
+            info[pointI].clear();
+        }
+    }
+}
+
+
+// ************************************************************************* //

src/meshTools/searchableSurface/searchableSurfaceWithGaps.H

@@ -0,0 +1,238 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::searchableSurfaceWithGaps
+
+Description
+    searchableSurface using multiple slightly shifted underlying surfaces
+    to make sure pierces don't go through gaps. 
+
+SourceFiles
+    searchableSurfaceWithGaps.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef searchableSurfaceWithGaps_H
+#define searchableSurfaceWithGaps_H
+
+#include "searchableSurface.H"
+#include "UPtrList.H"
+#include "Pair.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+
+/*---------------------------------------------------------------------------*\
+                           Class searchableSurfaceWithGaps Declaration
+\*---------------------------------------------------------------------------*/
+
+class searchableSurfaceWithGaps
+:
+    public searchableSurface
+{
+private:
+
+    // Private Member Data
+
+        //- gap size in meter
+        const scalar gap_;
+
+        //- Underlying geometry (size 1)
+        UPtrList<searchableSurface> subGeom_;
+
+
+    // Private Member Functions
+
+        Pair<vector> offsetVecs(const point&, const point&) const;
+
+        void offsetVecs
+        (
+            const pointField& start,
+            const pointField& end,
+            pointField& offset0,
+            pointField& offset1
+        ) const;
+
+        static label countMisses
+        (
+            const List<pointIndexHit>& info,
+            labelList& missMap
+        );
+
+        static label countMisses
+        (
+            const List<pointIndexHit>& plusInfo,
+            const List<pointIndexHit>& minInfo,
+            labelList& missMap
+        );
+
+
+        //- Disallow default bitwise copy construct
+        searchableSurfaceWithGaps(const searchableSurfaceWithGaps&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const searchableSurfaceWithGaps&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("searchableSurfaceWithGaps");
+
+
+    // Constructors
+
+        //- Construct from dictionary (used by searchableSurface)
+        searchableSurfaceWithGaps
+        (
+            const IOobject& io,
+            const dictionary& dict
+        );
+
+    // Destructor
+
+        virtual ~searchableSurfaceWithGaps();
+
+
+    // Member Functions
+
+        const searchableSurface& surface() const
+        {
+            return subGeom_[0];
+        }
+
+
+        virtual const wordList& regions() const
+        {
+            return surface().regions();
+        }
+
+        //- Whether supports volume type below
+        virtual bool hasVolumeType() const
+        {
+            return surface().hasVolumeType();
+        }
+
+        //- Range of local indices that can be returned.
+        virtual label size() const
+        {
+            return surface().size();
+        }
+
+
+        // Multiple point queries.
+
+            virtual void findNearest
+            (
+                const pointField& sample,
+                const scalarField& nearestDistSqr,
+                List<pointIndexHit>& info
+            ) const
+            {
+                surface().findNearest
+                (
+                    sample,
+                    nearestDistSqr,
+                    info
+                );
+            }
+
+            virtual void findLine
+            (
+                const pointField& start,
+                const pointField& end,
+                List<pointIndexHit>&
+            ) const;
+
+            virtual void findLineAny
+            (
+                const pointField& start,
+                const pointField& end,
+                List<pointIndexHit>&
+            ) const;
+
+            //- Get all intersections in order from start to end.
+            virtual void findLineAll
+            (
+                const pointField& start,
+                const pointField& end,
+                List<List<pointIndexHit> >&
+            ) const;
+
+            //- From a set of points and indices get the region
+            virtual void getRegion
+            (
+                const List<pointIndexHit>& info,
+                labelList& region
+            ) const
+            {
+                surface().getRegion(info, region);
+            }
+
+            //- From a set of points and indices get the normal
+            virtual void getNormal
+            (
+                const List<pointIndexHit>& info,
+                vectorField& normal
+            ) const
+            {
+                surface().getNormal(info, normal);
+            }
+
+            //- Determine type (inside/outside/mixed) for point. unknown if
+            //  cannot be determined (e.g. non-manifold surface)
+            virtual void getVolumeType
+            (
+                const pointField& samples,
+                List<volumeType>& info
+            ) const
+            {
+                surface().getVolumeType(samples, info);
+            }
+
+
+        // regIOobject implementation
+
+            bool writeData(Ostream& os) const
+            {
+                return surface().writeData(os);
+            }
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/meshTools/searchableSurface/searchableSurfaces.H

@@ -149,7 +149,7 @@ public:
                 const pointField& start,
                 const pointField& end,
                 labelListList& surfaces,
-                List<List<pointIndexHit> >& surfaceHits
+                List<List<pointIndexHit> >&
             ) const;
 
             //- Find nearest. Return -1 (and a miss()) or surface and nearest

src/meshTools/searchableSurface/triSurfaceMesh.C

@@ -28,6 +28,8 @@ License
 #include "Random.H"
 #include "addToRunTimeSelectionTable.H"
 #include "EdgeMap.H"
+#include "triSurfaceFields.H"
+#include "Time.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -299,57 +301,6 @@ const Foam::wordList& Foam::triSurfaceMesh::regions() const
 }
 
 
-//Foam::pointIndexHit Foam::triSurfaceMesh::findNearest
-//(
-//    const point& sample,
-//    const scalar nearestDistSqr
-//) const
-//{
-//    return tree().findNearest(sample, nearestDistSqr);
-//}
-//
-//
-//Foam::pointIndexHit Foam::triSurfaceMesh::findNearestOnEdge
-//(
-//    const point& sample,
-//    const scalar nearestDistSqr
-//) const
-//{
-//    return = edgeTree().findNearest(sample, nearestDistSqr);
-//}
-//
-//
-//Foam::pointIndexHit Foam::triSurfaceMesh::findNearest
-//(
-//    const linePointRef& ln,
-//    treeBoundBox& tightest,
-//    point& linePoint
-//) const
-//{
-//    return tree().findNearest(ln, tightest, linePoint);
-//}
-//
-//
-//Foam::pointIndexHit Foam::triSurfaceMesh::findLine
-//(
-//    const point& start,
-//    const point& end
-//) const
-//{
-//    return tree().findLine(start, end);
-//}
-//
-//
-//Foam::pointIndexHit Foam::triSurfaceMesh::findLineAny
-//(
-//    const point& start,
-//    const point& end
-//) const
-//{
-//    return tree().findLineAny(start, end);
-//}
-
-
 // Find out if surface is closed.
 bool Foam::triSurfaceMesh::hasVolumeType() const
 {
@@ -491,7 +442,6 @@ void Foam::triSurfaceMesh::findLineAll
 
             hits.shrink();
             info[pointI].transfer(hits);
-            hits.clear();
         }
         else
         {
@@ -545,6 +495,29 @@ void Foam::triSurfaceMesh::getNormal
 }
 
 
+void Foam::triSurfaceMesh::getField
+(
+    const word& fieldName,
+    const List<pointIndexHit>& info,
+    labelList& values
+) const
+{
+    const triSurfaceLabelField& fld = lookupObject<triSurfaceLabelField>
+    (
+        fieldName
+    );
+
+    values.setSize(info.size());
+    forAll(info, i)
+    {
+        if (info[i].hit())
+        {
+            values[i] = fld[info[i].index()];
+        }
+    }
+}
+
+
 void Foam::triSurfaceMesh::getVolumeType
 (
     const pointField& points,
@@ -564,4 +537,31 @@ void Foam::triSurfaceMesh::getVolumeType
 }
 
 
+//- Write using given format, version and compression
+bool Foam::triSurfaceMesh::writeObject
+(
+    IOstream::streamFormat fmt,
+    IOstream::versionNumber ver,
+    IOstream::compressionType cmp
+) const
+{
+    fileName fullPath(searchableSurface::objectPath());
+
+    if (!mkDir(fullPath.path()))
+    {
+        return false;
+    }
+
+    triSurface::write(fullPath);
+
+    if (!file(fullPath))
+    {
+        return false;
+    }
+
+    //return objectRegistry::writeObject(fmt, ver, cmp);
+    return true;
+}
+
+
 // ************************************************************************* //

src/meshTools/searchableSurface/triSurfaceMesh.H

@@ -36,6 +36,7 @@ SourceFiles
 #ifndef triSurfaceMesh_H
 #define triSurfaceMesh_H
 
+#include "treeBoundBox.H"
 #include "searchableSurface.H"
 #include "objectRegistry.H"
 #include "indexedOctree.H"
@@ -142,8 +143,11 @@ public:
             //- Whether supports volume type below. I.e. whether is closed.
             virtual bool hasVolumeType() const;
 
-
+            //- Range of local indices that can be returned.
-        // Multiple point queries.
+            virtual label size() const
+            {
+                return triSurface::size();
+            }
 
             virtual void findNearest
             (
@@ -196,6 +200,30 @@ public:
                 List<volumeType>&
             ) const;
 
+            //- Set bounds of surface. Bounds currently set as list of
+            //  bounding boxes. The bounds are hints to the surface as for
+            //  the range of queries it can expect. faceMap/pointMap can be
+            //  set if the surface has done any redistribution.
+            virtual void distribute
+            (
+                const List<treeBoundBox>&,
+                const bool keepNonLocal,
+                autoPtr<mapDistribute>& faceMap,
+                autoPtr<mapDistribute>& pointMap
+            )
+            {}
+
+        // Other
+
+            //- Specific to triSurfaceMesh: from a set of hits (points and
+            //  indices) get the specified field. Misses do not get set.
+            virtual void getField
+            (
+                const word& fieldName,
+                const List<pointIndexHit>&,
+                labelList& values
+            ) const;
+
 
         // regIOobject implementation
 
@@ -205,6 +233,14 @@ public:
                 return false;
             }
 
+            //- Write using given format, version and compression
+            virtual bool writeObject
+            (
+                IOstream::streamFormat fmt,
+                IOstream::versionNumber ver,
+                IOstream::compressionType cmp
+            ) const;
+
 };
 
 

src/meshTools/triSurface/booleanOps/intersectedSurface/edgeSurface.C

@@ -288,11 +288,9 @@ Foam::edgeSurface::edgeSurface
 
     // Transfer.
     allEdges.shrink();
-    allEdges.clear();
     edges_.transfer(allEdges);
 
     allParentEdges.shrink();
-    allParentEdges.clear();
     parentEdges_.transfer(allParentEdges);
 
     forAll(allFaceEdges, faceI)
@@ -300,7 +298,6 @@ Foam::edgeSurface::edgeSurface
         DynamicList<label>& allFEdges = allFaceEdges[faceI];
 
         allFEdges.shrink();
-        allFEdges.clear();
         faceEdges_[faceI].transfer(allFEdges);
     }
 

src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersectionTemplates.C

@@ -40,7 +40,6 @@ void Foam::surfaceIntersection::transfer
 {
     dList.shrink();
     lList.transfer(dList);
-    dList.clear();
 }
 
 

src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C

@@ -199,7 +199,6 @@ void Foam::surfaceFeatures::calcFeatPoints(const List<edgeStatus>& edgeStat)
     }
     featurePoints.shrink();
     featurePoints_.transfer(featurePoints);
-    featurePoints.clear();
 }
 
 

src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.H

@@ -36,7 +36,6 @@ SourceFiles
 #ifndef triSurfaceSearch_H
 #define triSurfaceSearch_H
 
-#include "DynamicList.H"
 #include "pointField.H"
 #include "boolList.H"
 #include "pointIndexHit.H"

src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C

@@ -300,11 +300,7 @@ Foam::triSurface Foam::triSurfaceTools::doRefine
     allPoints.transfer(newPoints);
     newPoints.clear();
 
-    List<labelledTri> allFaces;
+    return triSurface(newFaces, surf.patches(), allPoints, true);
-    allFaces.transfer(newFaces);
-    newFaces.clear();
-
-    return triSurface(allFaces, surf.patches(), allPoints);
 }
 
 
@@ -1988,7 +1984,7 @@ Foam::triSurface Foam::triSurfaceTools::greenRefine
     newFaces.shrink();
     newPoints.setSize(newPointI);
 
-    return triSurface(newFaces, surf.patches(), newPoints);
+    return triSurface(newFaces, surf.patches(), newPoints, true);
 }
 
 

src/meshTools/triSurface/triangleFuncs/triangleFuncs.C

@@ -101,8 +101,6 @@ bool Foam::triangleFuncs::intersectAxesBundle
     // Since direction is coordinate axis there is no need to do projection,
     // we can directly check u,v components for inclusion in triangle.
 
-    scalar localScale = max(max(magSqr(V10), magSqr(V20)), 1.0);
-
     // Get other components
     label i1 = (i0 + 1) % 3;
     label i2 = (i1 + 1) % 3;
@@ -113,6 +111,8 @@ bool Foam::triangleFuncs::intersectAxesBundle
     scalar u2 = V20[i1];
     scalar v2 = V20[i2];
 
+    scalar localScale = mag(u1)+mag(v1)+mag(u2)+mag(v2);
+
     scalar det = v2*u1 - u2*v1;
 
     // Fix for  V0:(-31.71428 0 -15.10714)
@@ -136,7 +136,7 @@ bool Foam::triangleFuncs::intersectAxesBundle
         scalar beta = 0;
         bool inter = false;
 
-        if (Foam::mag(u1)/localScale < SMALL)
+        if (Foam::mag(u1) < SMALL)
         {
             beta = u0/u2;
             if ((beta >= 0) && (beta <= 1))

src/sampling/cuttingPlane/cuttingPlane.C

@@ -154,7 +154,6 @@ Foam::labelList Foam::cuttingPlane::intersectEdges
 
     dynCuttingPoints.shrink();
     cutPoints_.transfer(dynCuttingPoints);
-    dynCuttingPoints.clear();
 
     return edgePoint;
 }

src/sampling/sampledSet/midPoint/midPointSet.C

@@ -25,7 +25,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "midPointSet.H"
-#include "DynamicList.H"
 #include "polyMesh.H"
 #include "addToRunTimeSelectionTable.H"
 

src/sampling/sampledSet/midPointAndFace/midPointAndFaceSet.C

@@ -25,7 +25,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "midPointAndFaceSet.H"
-#include "DynamicList.H"
 #include "polyMesh.H"
 #include "addToRunTimeSelectionTable.H"
 

src/thermophysicalModels/combustion/chemistryReaders/foamChemistryReader/foamChemistryReader.H

@@ -41,7 +41,6 @@ SourceFiles
 #include "typeInfo.H"
 #include "HashPtrTable.H"
 #include "SLPtrList.H"
-#include "DynamicList.H"
 #include "labelList.H"
 #include "speciesTable.H"
 #include "atomicWeights.H"

src/triSurface/tools/labelledTri/labelledTri.H

@@ -105,6 +105,10 @@ public:
 };
 
 
+template<>
+inline bool contiguous<labelledTri>()  {return true;}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/triSurface/tools/labelledTri/labelledTriI.H

@@ -58,16 +58,7 @@ inline labelledTri::labelledTri
 
 inline labelledTri::labelledTri(Istream& is)
 {
-    // Read beginning of labelledTri point pair
+    operator>>(is, *this);
-    is.readBegin("labelledTri");
-
-    is >> static_cast<triFace&>(*this) >> region_;
-
-    // Read end of labelledTri point pair
-    is.readEnd("labelledTri");
-
-    // Check state of Istream
-    is.check("labelledTri::labelledTri(Istream& is)");
 }
 
 
@@ -88,13 +79,20 @@ inline label& labelledTri::region()
 
 inline Istream& operator>>(Istream& is, labelledTri& t)
 {
-    // Read beginning of labelledTri point pair
+    if (is.format() == IOstream::ASCII)
-    is.readBegin("labelledTri");
+    {
+        // Read beginning of labelledTri point pair
+        is.readBegin("labelledTri");
 
-    is  >> static_cast<triFace&>(t) >> t.region_;
+        is  >> static_cast<triFace&>(t) >> t.region_;
 
-    // Read end of labelledTri point pair
+        // Read end of labelledTri point pair
-    is.readEnd("labelledTri");
+        is.readEnd("labelledTri");
+    }
+    else
+    {
+        is.read(reinterpret_cast<char*>(&t), sizeof(labelledTri));
+    }
 
     // Check state of Ostream
     is.check("Istream& operator>>(Istream&, labelledTri&)");
@@ -105,9 +103,24 @@ inline Istream& operator>>(Istream& is, labelledTri& t)
 
 inline Ostream& operator<<(Ostream& os, const labelledTri& t)
 {
-    os  << token::BEGIN_LIST
+    if (os.format() == IOstream::ASCII)
-        << static_cast<const triFace&>(t) << token::SPACE << t.region_
+    {
-        << token::END_LIST;
+        os  << token::BEGIN_LIST
+            << static_cast<const triFace&>(t) << token::SPACE << t.region_
+            << token::END_LIST;
+    }
+    else
+    {
+        os.write
+        (
+            reinterpret_cast<const char*>(&t),
+            sizeof(labelledTri)
+        );
+    }
+
+    // Check state of Ostream
+    os.check("Ostream& operator<<(Ostream&, const labelledTri&)");
+
 
     return os;
 }

src/triSurface/triSurface/interfaces/AC3D/readAC.C

@@ -339,11 +339,7 @@ bool triSurface::readAC(const fileName& ACfileName)
     allPoints.transfer(points);
     points.clear();
 
-    List<labelledTri> allFaces;
+    *this = triSurface(faces, patches, allPoints, true);
-    allFaces.transfer(faces);
-    faces.clear();
-
-    *this = triSurface(allFaces, patches, allPoints);
 
     stitchTriangles(allPoints);
 

src/triSurface/triSurface/interfaces/NAS/readNAS.C

@@ -350,7 +350,7 @@ bool triSurface::readNAS(const fileName& OBJfileName)
     points.clear();
 
     // Create triSurface
-    *this = triSurface(faces, patches, allPoints);
+    *this = triSurface(faces, patches, allPoints, true);
 
     return true;
 }

src/triSurface/triSurface/interfaces/OBJ/readOBJ.C

@@ -195,14 +195,9 @@ bool triSurface::readOBJ(const fileName& OBJfileName)
     // Transfer DynamicLists to straight ones.
     pointField allPoints;
     allPoints.transfer(points);
-    points.clear();
-
-    List<labelledTri> allFaces;
-    allFaces.transfer(faces);
-    faces.clear();
 
     // Create triSurface
-    *this = triSurface(allFaces, patches, allPoints);
+    *this = triSurface(faces, patches, allPoints, true);
 
     return true;
 }

src/triSurface/triSurface/triSurface.C

@@ -634,7 +634,8 @@ surfacePatchList triSurface::calcPatches(labelList& faceMap) const
     // Compact regions
 
     // Get last region
-    label maxRegion = 0;
+    //label maxRegion = 0;                  // for compacted regions
+    label maxRegion = patches_.size()-1;    // for non-compacted regions
 
     if (faceMap.size() > 0)
     {