Parallel GGI performance improvements

applications/utilities/surface/surfaceMeshConvertTesting/Make/options

@@ -1,5 +1,9 @@
 EXE_INC = \
+    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/triSurface/lnInclude \
     -I$(LIB_SRC)/surfMesh/lnInclude
 
-EXE_LIBS = -ltriSurface -lsurfMesh
+EXE_LIBS = \
+    -lmeshTools \
+    -ltriSurface \
+    -lsurfMesh

etc/controlDict

@@ -471,7 +471,7 @@ DebugSwitches
     genericPatch        0;
     geomCellLooper      0;
     geometricSurfacePatch   0;
-    ggi                 0;
+    GGIInterpolation    0;
     global              0;
     globalMeshData      0;
     globalPoints        0;

etc/cshrc

@@ -108,8 +108,8 @@ if ( ! $?WM_OSTYPE ) setenv WM_OSTYPE POSIX
 # Compiler: set to Gcc, Gcc43 or Icc (for Intel's icc)
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 #if ( ! $?WM_COMPILER ) setenv WM_COMPILER Gcc
-setenv WM_COMPILER Gcc
-#setenv WM_COMPILER Icc
+#setenv WM_COMPILER Gcc
+setenv WM_COMPILER Icc
 
 setenv WM_COMPILER_ARCH
 setenv WM_COMPILER_LIB_ARCH
@@ -127,8 +127,8 @@ if ( ! $?WM_PRECISION_OPTION ) setenv WM_PRECISION_OPTION DP
 
 # WM_COMPILE_OPTION = Opt | Debug | Prof
 #if ( ! $?WM_COMPILE_OPTION ) setenv WM_COMPILE_OPTION Opt
-setenv WM_COMPILE_OPTION Opt
-#setenv WM_COMPILE_OPTION Debug
+#setenv WM_COMPILE_OPTION Opt
+setenv WM_COMPILE_OPTION Debug
 
 # WM_MPLIB = | OPENMPI | MPICH | MPICH-GM | HPMPI | GAMMA | MPI | QSMPI
 if ( ! $?WM_MPLIB ) setenv WM_MPLIB OPENMPI

src/Allwmake

@@ -15,6 +15,7 @@ set -x
 wmakePrintBuild -check || /bin/rm -f OpenFOAM/Make/$WM_OPTIONS/global.? 2>/dev/null
 
 wmakeLnInclude OpenFOAM
+wmakeLnInclude meshTools
 wmakeLnInclude OSspecific/$WM_OSTYPE
 Pstream/Allwmake
 

src/OpenFOAM/Make/files

@@ -308,6 +308,7 @@ primitiveShapes = meshes/primitiveShapes
 
 $(primitiveShapes)/line/line.C
 $(primitiveShapes)/plane/plane.C
+$(primitiveShapes)/triangle/triangleFuncs.C
 $(primitiveShapes)/triangle/intersection.C
 
 meshShapes = meshes/meshShapes
@@ -495,6 +496,7 @@ $(pointBoundaryMesh)/pointBoundaryMesh.C
 meshes/boundBox/boundBox.C
 
 meshTools = meshes/meshTools
+$(meshTools)/meshTools.C
 $(meshTools)/matchPoints.C
 $(meshTools)/mergePoints.C
 
@@ -563,12 +565,33 @@ $(interpolations)/RBFInterpolation/RBFFunctions/Gauss/Gauss.C
 $(interpolations)/RBFInterpolation/RBFFunctions/TPS/TPS.C
 $(interpolations)/RBFInterpolation/RBFFunctions/IMQB/IMQB.C
 
+
 algorithms/MeshWave/MeshWaveName.C
 algorithms/MeshWave/FaceCellWaveName.C
 algorithms/polygon/clipping/SutherlandHodgman.C
 algorithms/polygon/pointInPolygon/HormannAgathos.C
 algorithms/rotation/RodriguesRotation.C
 
+octree = algorithms/octree/octree
+$(octree)/octreeName.C
+$(octree)/octreeDataPoint.C
+$(octree)/octreeDataPointTreeLeaf.C
+$(octree)/octreeDataEdges.C
+$(octree)/octreeDataCell.C
+$(octree)/octreeDataFace.C
+$(octree)/octreeDataBoundBox.C
+$(octree)/treeBoundBox.C
+$(octree)/treeNodeName.C
+$(octree)/treeLeafName.C
+$(octree)/pointIndexHitIOList.C
+
+indexedOctree = algorithms/octree/indexedOctree
+$(indexedOctree)/indexedOctreeName.C
+$(indexedOctree)/treeDataCell.C
+$(indexedOctree)/treeDataEdge.C
+$(indexedOctree)/treeDataFace.C
+$(indexedOctree)/treeDataPoint.C
+
 graph/curve/curve.C
 graph/graph.C
 

src/OpenFOAM/algorithms/octree/indexedOctree/treeDataCell.H

@@ -120,8 +120,8 @@ public:
 
         // Search
 
-            //- Get type (inside,outside,mixed,unknown) of point w.r.t. surface.
-            //  Only makes sense for closed surfaces.
+            //- Get type (inside,outside,mixed,unknown) of point w.r.t. surface
+            //  Only makes sense for closed surfaces
             label getVolumeType
             (
                 const indexedOctree<treeDataCell>&,

src/OpenFOAM/algorithms/octree/indexedOctree/treeDataPoint.C

@@ -77,8 +77,8 @@ bool Foam::treeDataPoint::overlaps
 }
 
 
-// Calculate nearest point to sample. Updates (if any) nearestDistSqr, minIndex,
-// nearestPoint.
+// Calculate nearest point to sample. Updates (if any) nearestDistSqr,
+//  minIndex, nearestPoint.
 void Foam::treeDataPoint::findNearest
 (
     const labelList& indices,

src/OpenFOAM/algorithms/octree/octree/octree.C

@@ -338,7 +338,7 @@ Foam::label Foam::octree<Type>::findNearest
 
 
 template <class Type>
-Foam::labelList Foam::octree<Type>::findBox(const boundBox& bb) const
+Foam::labelList Foam::octree<Type>::findBox(const treeBoundBox& bb) const
 {
     // Storage for labels of shapes inside bb. Size estimate.
     labelHashSet elements(100);

src/OpenFOAM/algorithms/octree/octree/octree.H

@@ -339,7 +339,7 @@ public:
 
             //- Find (in no particular order) indices of all shapes inside or
             //  overlapping bounding box (i.e. all shapes not outside box)
-            labelList findBox(const boundBox& bb) const;
+            labelList findBox(const treeBoundBox& bb) const;
 
             //- Find intersected shape along line. pointIndexHit contains index
             //  of nearest shape intersected and the intersection point. Uses

src/OpenFOAM/algorithms/octree/octree/octreeDataBoundBox.C

@@ -0,0 +1,210 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     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 "octreeDataBoundBox.H"
+#include "labelList.H"
+#include "polyMesh.H"
+#include "octree.H"
+#include "polyPatch.H"
+#include "triangleFuncs.H"
+#include "linePointRef.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(Foam::octreeDataBoundBox, 0);
+
+Foam::scalar Foam::octreeDataBoundBox::tol
+(
+    debug::tolerances("octreeDataBoundBoxTol", 1e-6)
+);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::octreeDataBoundBox::octreeDataBoundBox
+(
+    const treeBoundBoxList& bbL
+)
+:
+    allBb_(bbL)
+{}
+
+
+Foam::octreeDataBoundBox::octreeDataBoundBox(const octreeDataBoundBox& shapes)
+:
+    allBb_(shapes.allBb())
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::octreeDataBoundBox::~octreeDataBoundBox()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::label Foam::octreeDataBoundBox::getSampleType
+(
+    const octree<octreeDataBoundBox>& oc,
+    const point& sample
+) const
+{
+    return octree<octreeDataBoundBox>::UNKNOWN;
+}
+
+
+bool Foam::octreeDataBoundBox::overlaps
+(
+    const label index,
+    const treeBoundBox& sampleBb
+) const
+{
+    return sampleBb.overlaps(allBb_[index]);
+}
+
+
+bool Foam::octreeDataBoundBox::contains(const label, const point&) const
+{
+    notImplemented
+    (
+        "octreeDataBoundBox::contains(const label, const point&)"
+    );
+
+    return false;
+}
+
+
+bool Foam::octreeDataBoundBox::intersects
+(
+    const label index,
+    const point& start,
+    const point& end,
+    point& intersectionPoint
+) const
+{
+    notImplemented
+    (
+        "octreeDataBoundBox::intersects(const label, const point&)"
+    );
+    return false;
+}
+
+
+bool Foam::octreeDataBoundBox::findTightest
+(
+    const label index,
+    const point& sample,
+    treeBoundBox& tightest
+) const
+{
+    // Get furthest away vertex
+    point myNear, myFar;
+    allBb_[index].calcExtremities(sample, myNear, myFar);
+
+    const point dist = myFar - sample;
+    scalar myFarDist = mag(dist);
+
+    point tightestNear, tightestFar;
+    tightest.calcExtremities(sample, tightestNear, tightestFar);
+
+    scalar tightestFarDist = mag(tightestFar - sample);
+
+    if (tightestFarDist < myFarDist)
+    {
+        // Keep current tightest.
+        return false;
+    }
+    else
+    {
+        // Construct bb around sample and myFar
+        const point dist2(fabs(dist.x()), fabs(dist.y()), fabs(dist.z()));
+
+        tightest.min() = sample - dist2;
+        tightest.max() = sample + dist2;
+
+        return true;
+    }
+}
+
+
+// Determine numerical value of sign of sample compared to shape at index
+Foam::scalar Foam::octreeDataBoundBox::calcSign
+(
+    const label index,
+    const point& sample,
+    point& n
+) const
+{
+    n = vector::zero;
+
+    return 1;
+}
+
+
+// Calculate nearest point on/in shapei
+Foam::scalar Foam::octreeDataBoundBox::calcNearest
+(
+    const label index,
+    const point& sample,
+    point& nearest
+) const
+{
+    notImplemented
+    (
+        "octreeDataBoundBox::calcNearest"
+        "(const label index, const point& sample, point& nearest)"
+    );
+
+    return GREAT;
+}
+
+
+// Calculate nearest point on/in shapei
+Foam::scalar Foam::octreeDataBoundBox::calcNearest
+(
+    const label index,
+    const linePointRef& ln,
+    point& linePt,
+    point& shapePt
+) const
+{
+    notImplemented
+    (
+        "octreeDataBoundBox::calcNearest"
+        "(const label, const linePointRef&, point&, point&)"
+    );
+    return GREAT;
+}
+
+
+void Foam::octreeDataBoundBox::write(Ostream& os, const label index) const
+{
+    os << allBb_[index];
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/algorithms/octree/octree/octreeDataBoundBox.H

@@ -0,0 +1,213 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     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::octreeDataBoundBox
+
+Description
+    Holds data for octree to work on simple bounding boxes
+
+    For example, calculate (in calcNearest) the correct intersection point
+    with a bounding box.
+
+    Using bounding boxes simplify the usage of octree when only
+    primitivePatches entities are available (eg: from inside ggi interpolators)
+
+SourceFiles
+    octreeDataBoundBox.C
+
+Author
+    Martin Beaudoin, Hydro-Quebec, (2010)
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef octreeDataBoundBox_H
+#define octreeDataBoundBox_H
+
+#include "treeBoundBoxList.H"
+#include "faceList.H"
+#include "point.H"
+#include "className.H"
+#include "linePointRef.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class primitiveMesh;
+template<class Type> class octree;
+class polyPatch;
+
+/*---------------------------------------------------------------------------*\
+                           Class octreeDataBoundBox Declaration
+\*---------------------------------------------------------------------------*/
+
+class octreeDataBoundBox
+{
+    // Static data
+
+        //- Tolerance on linear dimensions
+        static scalar tol;
+
+
+    // Private data
+
+        //- Bounding boxes
+        treeBoundBoxList allBb_;
+
+
+    // Private Member Functions
+
+
+public:
+
+    // Declare name of the class and its debug switch
+    ClassName("octreeDataBoundBox");
+
+
+    // Constructors
+
+        //- Construct from list of precomputed of bound boxes
+        explicit octreeDataBoundBox
+        (
+            const treeBoundBoxList& bbL
+        );
+
+        //- Construct as copy
+        octreeDataBoundBox(const octreeDataBoundBox&);
+
+
+    // Destructor
+
+        ~octreeDataBoundBox();
+
+
+    // Member Functions
+
+        // Access
+
+            const treeBoundBoxList& allBb() const
+            {
+                return allBb_;
+            }
+
+            label size() const
+            {
+                return allBb_.size();
+            }
+
+
+        // Search
+
+            //- Get type of sample
+            label getSampleType
+            (
+                const octree<octreeDataBoundBox>&,
+                const point&
+            ) const;
+
+            //- Does (bb of) shape at index overlap bb
+            bool overlaps
+            (
+                const label index,
+                const treeBoundBox& sampleBb
+            ) const;
+
+            //- Does shape at index contain sample
+            bool contains
+            (
+                const label index,
+                const point& sample
+            ) const;
+
+            //- Segment (from start to end) intersection with shape
+            //  at index. If intersects returns true and sets intersectionPoint
+            bool intersects
+            (
+                const label index,
+                const point& start,
+                const point& end,
+                point& intersectionPoint
+            ) const;
+
+            //- Sets newTightest to bounding box (and returns true) if
+            //  nearer to sample than tightest bounding box. Otherwise
+            //  returns false.
+            bool findTightest
+            (
+                const label index,
+                const point& sample,
+                treeBoundBox& tightest
+            ) const;
+
+            //- Given index get unit normal and calculate (numerical) sign
+            //  of sample.
+            //  Used to determine accuracy of calcNearest or inside/outside.
+            scalar calcSign
+            (
+                const label index,
+                const point& sample,
+                vector& n
+            ) const;
+
+            //- Calculates nearest (to sample) point in shape.
+            //  Returns point and mag(nearest - sample). Returns GREAT if
+            //  sample does not project onto (triangle decomposition) of face.
+            scalar calcNearest
+            (
+                const label index,
+                const point& sample,
+                point& nearest
+            ) const;
+
+            //- Calculates nearest (to line segment) point in shape.
+            //  Returns distance and both point.
+            scalar calcNearest
+            (
+                const label index,
+                const linePointRef& ln,
+                point& linePt,          // nearest point on line
+                point& shapePt          // nearest point on shape
+            ) const;
+
+
+        // Write
+
+            //- Write shape at index
+            void write(Ostream& os, const label index) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/algorithms/octree/octree/octreeDataFace.C

@@ -545,6 +545,7 @@ bool Foam::octreeDataFace::overlaps
             return true;
         }
     }
+
     return false;
 }
 
@@ -555,6 +556,7 @@ bool Foam::octreeDataFace::contains(const label, const point&) const
     (
         "octreeDataFace::contains(const label, const point&)"
     );
+
     return false;
 }
 

src/OpenFOAM/algorithms/octree/octree/treeBoundBox.H

@@ -260,13 +260,13 @@ public:
             ) const;
 
             //- Overlaps other boundingbox?
-            inline bool overlaps(const treeBoundBox&) const;
+            inline bool overlaps(const boundBox&) const;
 
             //- Overlaps boundingSphere (centre + sqr(radius))?
             bool overlaps(const point&, const scalar radiusSqr) const;
 
-            //- Intersects segment; set point to intersection position and face,
-            //  return true if intersection found.
+            //- Intersects segment; set point to intersection position and face
+            //  Will return true if intersection found.
             //  (pt argument used during calculation even if not intersecting).
             //  Calculates intersections from outside supplied vector
             //  (overallStart, overallVec). This is so when
@@ -333,6 +333,12 @@ public:
             //  and guarantees in any direction s*mag(span) minimum width
             inline treeBoundBox extend(Random&, const scalar s) const;
 
+            //- Return slightly wider bounding box
+            //  Extends all dimensions with s*span
+            //  MB, 28/May/2011
+            inline treeBoundBox extend(const scalar s) const;
+
+
     // Friend Operators
 
         friend bool operator==(const treeBoundBox&, const treeBoundBox&);

src/OpenFOAM/algorithms/octree/octree/treeBoundBoxI.H

@@ -309,7 +309,7 @@ inline void Foam::treeBoundBox::searchOrder
 
 // true if bb's intersect or overlap.
 // Note: <= to make sure we catch all.
-inline bool Foam::treeBoundBox::overlaps(const treeBoundBox& bb) const
+inline bool Foam::treeBoundBox::overlaps(const boundBox& bb) const
 {
     return boundBox::overlaps(bb);
 }
@@ -347,6 +347,28 @@ inline Foam::treeBoundBox Foam::treeBoundBox::extend
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+//- Return slightly wider bounding box, without a random gen param
+//  So we are sure to stay away from very small offsets like rndGen == 0
+//  MB, 28/May/2011
+inline Foam::treeBoundBox Foam::treeBoundBox::extend(const scalar s) const
+{
+    treeBoundBox bb(*this);
+
+    vector span(bb.max() - bb.min());
+
+    // Make 3D
+    scalar magSpan = Foam::mag(span);
+
+    for (direction dir = 0; dir < vector::nComponents; dir++)
+    {
+        span[dir] = Foam::max(s*magSpan, span[dir]);
+    }
+
+    bb.min() -= cmptMultiply(s*vector::one, span);
+    bb.max() += cmptMultiply(s*vector::one, span);
+
+    return bb;
+}
+
 
 // ************************************************************************* //

src/OpenFOAM/algorithms/octree/octree/treeLeaf.C

@@ -334,7 +334,7 @@ template <class Type>
 bool Foam::treeLeaf<Type>::findBox
 (
     const Type& shapes,
-    const boundBox& box,
+    const treeBoundBox& box,
     labelHashSet& elements
 ) const
 {

src/OpenFOAM/algorithms/octree/octree/treeLeaf.H

@@ -235,7 +235,7 @@ public:
             bool findBox
             (
                 const Type& shapes,
-                const boundBox& bb,
+                const treeBoundBox& bb,
                 labelHashSet& elements
             ) const;
 
@@ -276,7 +276,7 @@ public:
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #include "octreeDataPointTreeLeaf.H"
-#include "octreeDataTriSurfaceTreeLeaf.H"
+//#include "octreeDataTriSurfaceTreeLeaf.H"  HJ, header moved to triSurfaceMesh.H
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/OpenFOAM/algorithms/octree/octree/treeNode.C

@@ -764,7 +764,7 @@ bool Foam::treeNode<Type>::findTightest
     // Estimate for best place to start searching
     label sampleOctant = this->bb().subOctant(midpoint(), sample, onEdge);
 
-    // Go into all suboctants (one containing sample first) and update tightest.
+    // Go into all suboctants (one containing sample first) and update tightest
     // Order of visiting is if e.g. sampleOctant = 5:
     //  5 1 2 3 4 0 6 7
     for (label octantI=0; octantI<8; octantI++)
@@ -846,7 +846,7 @@ bool Foam::treeNode<Type>::findNearest
     // Estimate for best place to start searching
     label sampleOctant = this->bb().subOctant(midpoint(), sample, onEdge);
 
-    // Go into all suboctants (one containing sample first) and update tightest.
+    // Go into all suboctants (one containing sample first) and update tightest
     // Order of visiting is if e.g. sampleOctant = 5:
     //  5 1 2 3 4 0 6 7
     for (label octantI=0; octantI<8; octantI++)
@@ -1004,7 +1004,7 @@ template <class Type>
 bool Foam::treeNode<Type>::findBox
 (
     const Type& shapes,
-    const boundBox& box,
+    const treeBoundBox& box,
     labelHashSet& elements
 ) const
 {
@@ -1018,7 +1018,7 @@ bool Foam::treeNode<Type>::findBox
         onEdge
     );
 
-    // Go into all suboctants (one containing sample first) and update tightest.
+    // Go into all suboctants (one containing sample first) and update tightest
     // Order of visiting is if e.g. sampleOctant = 5:
     //  5 1 2 3 4 0 6 7
     for (label octantI=0; octantI<8; octantI++)

src/OpenFOAM/algorithms/octree/octree/treeNode.H

@@ -271,7 +271,7 @@ public:
             bool findBox
             (
                 const Type& shapes,
-                const boundBox& bb,
+                const treeBoundBox& bb,
                 labelHashSet& elements
             ) const;
 

src/OpenFOAM/db/IOstreams/Pstreams/Pstream.H

@@ -279,7 +279,7 @@ public:
         //- Am I the master process
         static bool master()
         {
-            return myProcNo_ == 0;
+            return myProcNo_ == masterNo();
         }
 
         //- Process index of the master

src/OpenFOAM/db/IOstreams/Pstreams/gatherScatterList.C

@@ -25,9 +25,10 @@ License
 Description
     Gather data from all processors onto single processor according to some
     communication schedule (usually linear-to-master or tree-to-master).
-    The gathered data will be a list with element procID the data from processor
-    procID. Before calling every processor should insert its value into
-    Values[Pstream::myProcNo()].
+    The gathered data will be a list with element procID the data from
+     processor procID. Before calling every processor should insert
+    its value into Values[Pstream::myProcNo()].
+
     Note: after gather every processor only knows its own data and that of the
     processors below it. Only the 'master' of the communication schedule holds
     a fully filled List. Use scatter to distribute the data.

src/OpenFOAM/interpolations/GGIInterpolation/GGIInterpolate.C

@@ -60,6 +60,38 @@ void GGIInterpolation<MasterPatch, SlavePatch>::interpolate
 }
 
 
+template<class MasterPatch, class SlavePatch>
+template<class Type>
+void GGIInterpolation<MasterPatch, SlavePatch>::maskedInterpolate
+(
+    const Field<Type>& ff,
+    Field<Type>& result,
+    const labelList& mask,
+    const labelListList& addr,
+    const scalarListList& weights
+)
+{
+    forAll (mask, maskI)
+    {
+        // Pick the masked face
+        const label faceI = mask[maskI];
+
+        const labelList& curAddr = addr[faceI];
+
+        const scalarList& curWeights = weights[faceI];
+
+        // Clear condensed list entry: masked faces only
+        result[maskI] = pTraits<Type>::zero;
+
+        forAll (curAddr, i)
+        {
+            // Put the result into condensed list: masked faces only
+            result[maskI] += ff[curAddr[i]]*curWeights[i];
+        }
+    }
+}
+
+
 template<class MasterPatch, class SlavePatch>
 template<class Type>
 void GGIInterpolation<MasterPatch, SlavePatch>::bridge
@@ -76,6 +108,44 @@ void GGIInterpolation<MasterPatch, SlavePatch>::bridge
 }
 
 
+template<class MasterPatch, class SlavePatch>
+template<class Type>
+void GGIInterpolation<MasterPatch, SlavePatch>::maskedBridge
+(
+    const Field<Type>& bridgeField,
+    Field<Type>& ff,
+    const labelList& mask,
+    const labelList& addr
+)
+{
+    // Note: tricky algorithm
+    // In order for a face to be bridged it needs to be both in the
+    // mask and in selection of faces that are bridged (addr).
+    // This implies an n-squared search, but we can use the fact that
+    // both lists are ordered.
+
+    label curAddrI = 0;
+
+    forAll (mask, maskI)
+    {
+        // Pick the masked face
+        const label faceI = mask[maskI];
+
+        for (; curAddrI < addr.size(); curAddrI++)
+        {
+            if (faceI == addr[curAddrI])
+            {
+                // Found masked bridged face
+                // Put the result into condensed list: masked faces only
+                ff[maskI] = bridgeField[faceI];
+
+                break;
+            }
+        }
+    }
+}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class MasterPatch, class SlavePatch>
@@ -161,6 +231,84 @@ GGIInterpolation<MasterPatch, SlavePatch>::masterToSlave
 }
 
 
+template<class MasterPatch, class SlavePatch>
+template<class Type>
+void GGIInterpolation<MasterPatch, SlavePatch>::maskedMasterToSlave
+(
+    const Field<Type>& ff,
+    Field<Type>& result,
+    const labelList& mask
+) const
+{
+    if (ff.size() != masterPatch_.size())
+    {
+        FatalErrorIn
+        (
+            "void GGIInterpolation::maskedMasterToSlave\n"
+            "(\n"
+            "    const Field<Type>& ff,\n"
+            "    Field<Type>& result,\n"
+            "    const labelList& mask\n"
+            ") const"
+        )   << "given field does not correspond to patch. Patch size: "
+            << masterPatch_.size() << " field size: " << ff.size()
+            << abort(FatalError);
+    }
+
+    if (result.size() != mask.size())
+    {
+        FatalErrorIn
+        (
+            "void GGIInterpolation::maskedMasterToSlave\n"
+            "(\n"
+            "    const Field<Type>& ff,\n"
+            "    Field<Type>& result,\n"
+            "    const labelList& mask\n"
+            ") const"
+        )   << "result field does not correspond to mask. Field size: "
+            << result.size() << " mask size: " << mask.size()
+            << abort(FatalError);
+    }
+
+    if (this->doTransform() && pTraits<Type>::rank > 0)
+    {
+        // Transform master data to slave
+        Field<Type> transformFF;
+
+        if (reverseT_.size() == 1)
+        {
+            // Constant transform
+            transformFF = transform(reverseT_[0], ff);
+        }
+        else
+        {
+            // Full patch transform
+            transformFF = transform(reverseT_, ff);
+        }
+
+        GGIInterpolation<MasterPatch, SlavePatch>::maskedInterpolate
+        (
+            transformFF,
+            result,
+            mask,
+            this->slaveAddr(),
+            this->slaveWeights()
+        );
+    }
+    else
+    {
+        GGIInterpolation<MasterPatch, SlavePatch>::maskedInterpolate
+        (
+            ff,
+            result,
+            mask,
+            this->slaveAddr(),
+            this->slaveWeights()
+        );
+    }
+}
+
+
 template<class MasterPatch, class SlavePatch>
 template<class Type>
 tmp<Field<Type> >
@@ -244,6 +392,83 @@ GGIInterpolation<MasterPatch, SlavePatch>::slaveToMaster
 }
 
 
+template<class MasterPatch, class SlavePatch>
+template<class Type>
+void GGIInterpolation<MasterPatch, SlavePatch>::maskedSlaveToMaster
+(
+    const Field<Type>& ff,
+    Field<Type>& result,
+    const labelList& mask
+) const
+{
+    if (ff.size() != slavePatch_.size())
+    {
+        FatalErrorIn
+        (
+            "void GGIInterpolation::maskedSlaveToMaster"
+            "(\n"
+            "    const Field<Type>& ff,\n"
+            "    Field<Type>& result,\n"
+            "    const labelList& mask\n"
+            ") const"
+        )   << "given field does not correspond to patch. Patch size: "
+            << slavePatch_.size() << " field size: " << ff.size()
+            << abort(FatalError);
+    }
+
+    if (result.size() != mask.size())
+    {
+        FatalErrorIn
+        (
+            "void GGIInterpolation::maskedSlaveToMaster\n"
+            "(\n"
+            "    const Field<Type>& ff,\n"
+            "    Field<Type>& result,\n"
+            "    const labelList& mask\n"
+            ") const"
+        )   << "result field does not correspond to mask. Field size: "
+            << result.size() << " mask size: " << mask.size()
+            << abort(FatalError);
+    }
+
+    if (this->doTransform() && pTraits<Type>::rank > 0)
+    {
+        // Transform slave data to master
+        Field<Type> transformFF;
+        if (forwardT_.size() == 1)
+        {
+            // Constant transform
+            transformFF = transform(forwardT_[0], ff);
+        }
+        else
+        {
+            // Full patch transform
+            transformFF = transform(forwardT_, ff);
+        }
+
+        GGIInterpolation<MasterPatch, SlavePatch>::maskedInterpolate
+        (
+            transformFF,
+            result,
+            mask,
+            this->masterAddr(),
+            this->masterWeights()
+        );
+    }
+    else
+    {
+        GGIInterpolation<MasterPatch, SlavePatch>::maskedInterpolate
+        (
+            ff,
+            result,
+            mask,
+            this->masterAddr(),
+            this->masterWeights()
+        );
+    }
+}
+
+
 template<class MasterPatch, class SlavePatch>
 template<class Type>
 void GGIInterpolation<MasterPatch, SlavePatch>::bridgeMaster
@@ -275,6 +500,42 @@ void GGIInterpolation<MasterPatch, SlavePatch>::bridgeMaster
 }
 
 
+template<class MasterPatch, class SlavePatch>
+template<class Type>
+void GGIInterpolation<MasterPatch, SlavePatch>::maskedBridgeMaster
+(
+    const Field<Type>& bridgeField,
+    Field<Type>& ff,
+    const labelList& mask
+) const
+{
+    if
+    (
+        bridgeField.size() != masterPatch_.size()
+     || ff.size() != mask.size()
+    )
+    {
+        FatalErrorIn
+        (
+            "void GGIInterpolation<MasterPatch, SlavePatch>::"
+            "maskedBridgeMaster\n"
+            "(\n"
+            "    const Field<Type>& bridgeField,\n"
+            "    Field<Type>& ff,\n"
+            "    const labelList& mask\n"
+            ") const"
+        )   << "given field does not correspond to patch. Patch size: "
+            << masterPatch_.size()
+            << " bridge field size: " << bridgeField.size()
+            << " field size: " << ff.size()
+            << " mask size: " << mask.size()
+            << abort(FatalError);
+    }
+
+    maskedBridge(bridgeField, ff, mask, uncoveredMasterFaces());
+}
+
+
 template<class MasterPatch, class SlavePatch>
 template<class Type>
 void GGIInterpolation<MasterPatch, SlavePatch>::bridgeSlave
@@ -286,14 +547,16 @@ void GGIInterpolation<MasterPatch, SlavePatch>::bridgeSlave
     if
     (
         bridgeField.size() != slavePatch_.size()
-     || ff.size() != slavePatch_.size())
+     || ff.size() != slavePatch_.size()
+    )
     {
         FatalErrorIn
         (
-            "void GGIInterpolation<MasterPatch, SlavePatch>::bridgeSlave\n"
+            "void GGIInterpolation<MasterPatch, SlavePatch>::"
+            "bridgeSlave\n"
             "(\n"
             "    const Field<Type>& bridgeField,\n"
-            "    Field<Type>& ff\n"
+            "    Field<Type>& ff"
             ") const"
         )   << "given field does not correspond to patch. Patch size: "
             << slavePatch_.size()
@@ -306,6 +569,42 @@ void GGIInterpolation<MasterPatch, SlavePatch>::bridgeSlave
 }
 
 
+template<class MasterPatch, class SlavePatch>
+template<class Type>
+void GGIInterpolation<MasterPatch, SlavePatch>::maskedBridgeSlave
+(
+    const Field<Type>& bridgeField,
+    Field<Type>& ff,
+    const labelList& mask
+) const
+{
+    if
+    (
+        bridgeField.size() != slavePatch_.size()
+     || ff.size() != mask.size()
+    )
+    {
+        FatalErrorIn
+        (
+            "void GGIInterpolation<MasterPatch, SlavePatch>::"
+            "maskedBridgeSlave\n"
+            "(\n"
+            "    const Field<Type>& bridgeField,\n"
+            "    Field<Type>& ff\n,"
+            "    const labelList& mask\n"
+            ") const"
+        )   << "given field does not correspond to patch. Patch size: "
+            << slavePatch_.size()
+            << " bridge field size: " << bridgeField.size()
+            << " field size: " << ff.size()
+            << " mask size: " << mask.size()
+            << abort(FatalError);
+    }
+
+    maskedBridge(bridgeField, ff, mask, uncoveredSlaveFaces());
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/interpolations/GGIInterpolation/GGIInterpolation.H

@@ -85,6 +85,7 @@ public:
         {
             THREE_D_DISTANCE,
             AABB,
+            BB_OCTREE,
             N_SQUARED
         };
 
@@ -192,6 +193,39 @@ class GGIInterpolation
         //- Facet normal featureCos criteria
         static const scalar featureCosTol_;
 
+       //- Facet bound box extension factor
+       static const scalar faceBoundBoxExtendSpanFraction_;
+
+       //- The next 3 attributes are parameters controlling the creation
+       //  of an octree search engine for the GGI facets neighbourhood
+       //  determination.
+       //
+       //  Tests using GGI patches of up to ~100K facets have shown that
+       //  the following values gave the best compromise between the
+       //  "octree creation time" vs "octree search time":
+       //
+       //     octreeSearchMinNLevel_     : 3
+       //     octreeSearchMaxLeafRatio_  : 3
+       //     octreeSearchMaxShapeRatio_ : 1
+       //
+       //  For GGI patches larger than ~100K facets, your mileage may vary.
+       //  So these 3 control parameters are adjustable using the following
+       //  global optimization switches:
+       //
+       //     GGIOctreeSearchMinNLevel
+       //     GGIOctreeSearchMaxLeafRatio
+       //     GGIOctreeSearchMaxShapeRatio
+       //
+
+       //- Octree search: minNlevel parameter for octree constructor
+       static const label octreeSearchMinNLevel_;
+
+       //- Octree search: maxLeafRatio parameter for octree constructor
+       static const scalar octreeSearchMaxLeafRatio_;
+
+       //- Octree search: maxShapeRatio parameter for octree constructor
+       static const scalar octreeSearchMaxShapeRatio_;
+
 
     // Private Member Functions
 
@@ -210,6 +244,10 @@ class GGIInterpolation
         //- Evaluate faces neighborhood based of faces Axis Aligned BB
         void findNeighboursAABB(labelListList& result) const;
 
+        //- Evaluate faces neighborhood based of faces BB and an octree
+        //  search engine
+        void findNeighboursBBOctree(labelListList& result) const;
+
         //- Projects a list of points onto a plane located at
         //  planeOrig, oriented along planeNormal
         tmp<pointField> projectPointsOnPlane
@@ -326,6 +364,18 @@ class GGIInterpolation
             const scalarListList& weights
         );
 
+        //- Interpolate given source and target, addressing and weights
+        //  for masked faces only
+        template<class Type>
+        static void maskedInterpolate
+        (
+            const Field<Type>& ff,
+            Field<Type>& result,
+            const labelList& mask,
+            const labelListList& addr,
+            const scalarListList& weights
+        );
+
         //- Bridge uncovered faces given addressing
         template<class Type>
         static void bridge
@@ -335,6 +385,17 @@ class GGIInterpolation
             const labelList& addr
         );
 
+        //- Bridge uncovered faces given addressing
+        //  for masked faces only
+        template<class Type>
+        static void maskedBridge
+        (
+            const Field<Type>& bridgeField,
+            Field<Type>& ff,
+            const labelList& mask,
+            const labelList& addr
+        );
+
         //- Is a transform required?
         inline bool doTransform() const
         {
@@ -360,8 +421,8 @@ public:
             const tensorField& forwardT,
             const tensorField& reverseT,
             const vectorField& forwardSep,
-            const scalar masterFaceNonOverlapFaceTol = 0.0,
-            const scalar slaveFaceNonOverlapFaceTol = 0.0,
+            const scalar masterFaceNonOverlapFaceTol = 0,
+            const scalar slaveFaceNonOverlapFaceTol = 0,
             const bool rescaleGGIWeightingFactors = true,
             const quickReject reject = AABB
         );
@@ -399,22 +460,45 @@ public:
 
         //- Interpolate from master to slave
         template<class Type>
-        tmp<Field<Type> > masterToSlave(const Field<Type>& pf) const;
+        tmp<Field<Type> > masterToSlave(const Field<Type>& ff) const;
 
         template<class Type>
         tmp<Field<Type> > masterToSlave
         (
-            const tmp<Field<Type> >& tpf
+            const tmp<Field<Type> >& tff
         ) const;
 
+        //- Interpolate from master to slave, only for marked master faces
+        //  Addressing picks the faces from patch that are selected
+        //  for collection into the result field
+        template<class Type>
+        void maskedMasterToSlave
+        (
+            const Field<Type>& ff,
+            Field<Type>& result,
+            const labelList& mask
+        ) const;
+
+
         //- Interpolate from slave to master
         template<class Type>
-        tmp<Field<Type> > slaveToMaster(const Field<Type>& pf) const;
+        tmp<Field<Type> > slaveToMaster(const Field<Type>& ff) const;
 
         template<class Type>
         tmp<Field<Type> > slaveToMaster
         (
-            const tmp<Field<Type> >& tpf
+            const tmp<Field<Type> >& tff
+        ) const;
+
+        //- Interpolate from slave to master, only for marked slave faces
+        //  Addressing picks the faces from patch that are selected
+        //  for collection into the result field
+        template<class Type>
+        void maskedSlaveToMaster
+        (
+            const Field<Type>& ff,
+            Field<Type>& result,
+            const labelList& mask
         ) const;
 
 
@@ -426,6 +510,15 @@ public:
             Field<Type>& ff
         ) const;
 
+        //- Bridge uncovered master patch field, only for marked master faces
+        template<class Type>
+        void maskedBridgeMaster
+        (
+            const Field<Type>& bridgeField,
+            Field<Type>& ff,
+            const labelList& mask
+        ) const;
+
         //- Bridge uncovered slave patch field
         template<class Type>
         void bridgeSlave
@@ -434,6 +527,15 @@ public:
             Field<Type>& ff
         ) const;
 
+        //- Bridge uncovered slave patch field, only for marked slave faces
+        template<class Type>
+        void maskedBridgeSlave
+        (
+            const Field<Type>& bridgeField,
+            Field<Type>& ff,
+            const labelList& mask
+        ) const;
+
 
     // Edit
 

src/OpenFOAM/interpolations/GGIInterpolation/GGIInterpolationQuickRejectTests.C

@@ -34,12 +34,45 @@ Author
 #include "boundBox.H"
 #include "plane.H"
 #include "transformField.H"
+#include "octree.H"
+#include "octreeDataBoundBox.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+template<class MasterPatch, class SlavePatch>
+const scalar
+GGIInterpolation<MasterPatch, SlavePatch>::faceBoundBoxExtendSpanFraction_
+(
+    debug::tolerances("GGIFaceBoundBoxExtendSpanFraction", 1.0e-2)
+);
+
+template<class MasterPatch, class SlavePatch>
+const label
+ GGIInterpolation<MasterPatch, SlavePatch>::octreeSearchMinNLevel_
+(
+    debug::optimisationSwitch("GGIOctreeSearchMinNLevel", 3)
+);
+
+template<class MasterPatch, class SlavePatch>
+const scalar
+GGIInterpolation<MasterPatch, SlavePatch>::octreeSearchMaxLeafRatio_
+(
+    debug::optimisationSwitch("GGIOctreeSearchMaxLeafRatio", 3)
+);
+
+template<class MasterPatch, class SlavePatch>
+const scalar
+GGIInterpolation<MasterPatch, SlavePatch>::octreeSearchMaxShapeRatio_
+(
+    debug::optimisationSwitch("GGIOctreeSearchMaxShapeRatio", 1)
+);
+
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 // From: http://www.gamasutra.com/features/20000330/bobic_02.htm
@@ -382,6 +415,148 @@ void GGIInterpolation<MasterPatch, SlavePatch>::findNeighboursAABB
     }
 }
 
+// This algorithm find the faces in proximity of another face based
+// on the face BB (Bounding Box) and an octree of bounding boxes.
+template<class MasterPatch, class SlavePatch>
+void GGIInterpolation<MasterPatch, SlavePatch>::findNeighboursBBOctree
+(
+    labelListList& result
+) const
+{
+    List<DynamicList<label> > candidateMasterNeighbors(masterPatch_.size());
+
+    // Initialize the list of master patch faces bounding box
+    treeBoundBoxList lmasterFaceBB(masterPatch_.size());
+
+    forAll (masterPatch_, faceMi)
+    {
+        pointField facePoints
+        (
+            masterPatch_[faceMi].points(masterPatch_.points())
+        );
+
+        // Construct face BB with an extension of face span defined by the
+        //  global tolerance factor faceBoundBoxExtendSpanFraction_
+        // (1% by default)
+        treeBoundBox bbFaceMaster(facePoints);
+
+        lmasterFaceBB[faceMi] =
+            bbFaceMaster.extend(faceBoundBoxExtendSpanFraction_);
+    }
+
+    // Initialize the list of slave patch faces bounding box
+    treeBoundBoxList lslaveFaceBB(slavePatch_.size());
+
+    forAll (slavePatch_, faceSi)
+    {
+        pointField facePoints
+        (
+            slavePatch_[faceSi].points(slavePatch_.points())
+        );
+
+        // possible transformation and separation for cyclic patches
+        if (doTransform())
+        {
+            if (forwardT_.size() == 1)
+            {
+                transform(facePoints, forwardT_[0], facePoints);
+            }
+            else
+            {
+                transform(facePoints, forwardT_[faceSi], facePoints);
+            }
+        }
+
+        if (doSeparation())
+        {
+            if (forwardSep_.size() == 1)
+            {
+                facePoints += forwardSep_[0];
+            }
+            else
+            {
+                facePoints += forwardSep_[faceSi];
+            }
+        }
+
+        // Construct face BB with an extension of face span defined by the
+        //  global tolerance factor faceBoundBoxExtendSpanFraction_
+        // (1% by default)
+        treeBoundBox bbFaceSlave(facePoints);
+
+        lslaveFaceBB[faceSi] =
+            bbFaceSlave.extend(faceBoundBoxExtendSpanFraction_);
+    }
+
+    // Create the slave octreeData, using the boundBox flavor
+    octreeDataBoundBox slaveDataBB(lslaveFaceBB);
+
+    // Overall slave patch BB
+    treeBoundBox slaveOverallBB(slavePatch_.points());
+
+    // Create the slave patch octree
+
+    octree<octreeDataBoundBox> slavePatchOctree
+    (
+        slaveOverallBB,              // overall search domain
+        slaveDataBB,
+        octreeSearchMinNLevel_,      // min number of levels
+        octreeSearchMaxLeafRatio_,   // max avg. size of leaves
+        octreeSearchMaxShapeRatio_   // max avg. duplicity.
+    );
+
+    const vectorField& masterFaceNormals = masterPatch_.faceNormals();
+    vectorField slaveNormals = slavePatch_.faceNormals();
+
+    // Transform slave normals to master plane if needed
+    if (doTransform())
+    {
+        if (forwardT_.size() == 1)
+        {
+            transform(slaveNormals, forwardT_[0], slaveNormals);
+        }
+        else
+        {
+            transform(slaveNormals, forwardT_, slaveNormals);
+        }
+    }
+
+    // Visit each master patch face BB and find the potential neighbours
+    // using the slave patch octree
+
+    forAll (lmasterFaceBB, faceMi)
+    {
+        // List of candidate neighbours
+        labelList overlappedFaces  =
+            slavePatchOctree.findBox(lmasterFaceBB[faceMi]);
+
+        forAll (overlappedFaces, ovFi)
+        {
+            label faceSi = overlappedFaces[ovFi];
+
+            // Compute and verify featureCos between the two face normals
+            // before adding to the list of candidates
+            scalar featureCos =
+                masterFaceNormals[faceMi] & slaveNormals[faceSi];
+
+            if (mag(featureCos) > featureCosTol_)
+            {
+                candidateMasterNeighbors[faceMi].append(faceSi);
+            }
+        }
+    }
+
+    // Repack the list
+    result.setSize(masterPatch_.size());
+
+    forAll (result, i)
+    {
+        result[i].transfer(candidateMasterNeighbors[i].shrink());
+    }
+
+    return;
+}
+
 
 
 // Projects a list of points onto a plane located at planeOrig,
@@ -664,8 +839,8 @@ bool GGIInterpolation<MasterPatch, SlavePatch>::detect2dPolygonsOverlap
         // accept to discard an intersecting area roughly 10e-6 times the
         // surface of the smallest polygon, so 1 PPM.
         //
-        // Which means also that our GGI weighting factors will never be smaller
-        // than roughly 10e-6, because this is the fraction of
+        // Which means also that our GGI weighting factors will never be
+        //  smaller than roughly 10e-6, because this is the fraction of
         // the intersection surface area we choose to discard.
 
         scalar _epsilon = tolFactor*
@@ -691,12 +866,14 @@ bool GGIInterpolation<MasterPatch, SlavePatch>::detect2dPolygonsOverlap
     {
         // We have not found any separating axes by exploring from
         // poly1, let's switch by exploring from poly2 instead
-        isOverlapping = detect2dPolygonsOverlap(poly2, poly1, tolFactor, false);
+        isOverlapping =
+            detect2dPolygonsOverlap(poly2, poly1, tolFactor, false);
     }
 
     return isOverlapping;
 }
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/interpolations/GGIInterpolation/GGIInterpolationWeights.C

@@ -87,7 +87,11 @@ void GGIInterpolation<MasterPatch, SlavePatch>::calcAddressing() const
 
     if (debug)
     {
-        Info << "Evaluation of GGI weighting factors:" << endl;
+        InfoIn
+        (
+            "void GGIInterpolation<MasterPatch, SlavePatch>::"
+            "calcAddressing() const"
+        )   << "Evaluation of GGI weighting factors:" << endl;
     }
 
     // Create the dynamic lists to hold the addressing
@@ -112,14 +116,19 @@ void GGIInterpolation<MasterPatch, SlavePatch>::calcAddressing() const
     // The candidates master neighbours
     // Choice of algorithm:
     // 1) Axis-aligned bounding box
-    // 2) 3-D vector distance
-    // 3) n-Squared search
+    // 2) Octree search with bounding box
+    // 3) 3-D vector distance
+    // 4) n-Squared search
     labelListList candidateMasterNeighbors;
 
     if (reject_ == AABB)
     {
          findNeighboursAABB(candidateMasterNeighbors);
     }
+    else if (reject_ == BB_OCTREE)
+    {
+         findNeighboursBBOctree(candidateMasterNeighbors);
+    }
     else if (reject_ == THREE_D_DISTANCE)
     {
          findNeighbours3D(candidateMasterNeighbors);
@@ -500,7 +509,7 @@ void GGIInterpolation<MasterPatch, SlavePatch>::calcAddressing() const
 // other.  A good case: the ercoftac conical diffuser, Case0...  GGI
 // located right between the cylindrical and conical parts, rotate the
 // cylindrical by 15 degrees.  For this case, we will need to devise a
-// decent strategy in order to intelligently take care of this
+// decent strategy in order to intelligently take care of these
 // "missing weights"
 //
 // The purpose of the ::rescaleWeightingFactors() method is mainly for
@@ -584,6 +593,7 @@ void GGIInterpolation<MasterPatch, SlavePatch>::rescaleWeightingFactors() const
     }
 }
 
+
 // Find non-overlapping faces from both master and slave patches
 // The default non-overlapping criteria is total absence of neighbours.
 // Later on, ths criteria will be based on minimum surface intersection, or

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/ggiGAMGInterfaceField/ggiGAMGInterfaceField.C

@@ -107,12 +107,29 @@ void Foam::ggiGAMGInterfaceField::updateInterfaceMatrix
 
     const unallocLabelList& faceCells = ggiInterface_.faceCells();
 
+    // New treatment.  HJ, 26/Jun/2011
+    if (zonePnf.size() != faceCells.size())
+    {
+        FatalErrorIn("ggiGAMGInterfaceField::updateInterfaceMatrix")
+            << "Bananas!!!"
+            << abort(FatalError);
+    }
+
+    forAll(faceCells, elemI)
+    {
+        result[faceCells[elemI]] -= coeffs[elemI]*zonePnf[elemI];
+    }
+
+
+    // Old treatment
+#if(0)
     const labelList& za = ggiInterface_.zoneAddressing();
 
     forAll(faceCells, elemI)
     {
         result[faceCells[elemI]] -= coeffs[elemI]*zonePnf[za[elemI]];
     }
+#endif
 }
 
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaces/ggiGAMGInterface/ggiGAMGInterface.C

@@ -41,6 +41,49 @@ namespace Foam
 }
 
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::ggiGAMGInterface::initFastReduce() const
+{
+    // Init should be executed only once
+    initReduce_ = true;
+
+    // Establish parallel comms pattern
+    if (!localParallel() && Pstream::parRun())
+    {
+        // Only master handles communication
+        if (Pstream::master())
+        {
+            receiveAddr_.setSize(Pstream::nProcs());
+            sendAddr_.setSize(Pstream::nProcs());
+
+            receiveAddr_[0] = zoneAddressing();
+
+            for (label procI = 1; procI < Pstream::nProcs(); procI++)
+            {
+                // Note: must use normal comms because the size of the
+                // communicated lists is unknown on the receiving side
+                // HJ, 4/Jun/2011
+
+                // Opt: reconsider mode of communication
+                IPstream ip(Pstream::scheduled, procI);
+
+                receiveAddr_[procI] = labelList(ip);
+
+                sendAddr_[procI] = labelList(ip);
+            }
+        }
+        else
+        {
+            // Opt: reconsider mode of communication
+            OPstream op(Pstream::scheduled, Pstream::masterNo());
+
+            op << zoneAddressing() << shadowInterface().zoneAddressing();
+        }
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::ggiGAMGInterface::ggiGAMGInterface
@@ -54,7 +97,10 @@ Foam::ggiGAMGInterface::ggiGAMGInterface
     GAMGInterface(lduMesh),
     fineGgiInterface_(refCast<const ggiLduInterface>(fineInterface)),
     zoneSize_(0),
-    zoneAddressing_()
+    zoneAddressing_(),
+    initReduce_(false),
+    receiveAddr_(),
+    sendAddr_()
 {
     // Note.
     // All processors will do the same coarsening and then filter
@@ -644,6 +690,7 @@ const Foam::labelListList& Foam::ggiGAMGInterface::addressing() const
     FatalErrorIn("const labelListList& ggiGAMGInterface::addressing() const")
         << "Requested fine addressing at coarse level"
         << abort(FatalError);
+
     return labelListList::null();
 }
 
@@ -659,6 +706,7 @@ const Foam::scalarListList& Foam::ggiGAMGInterface::weights() const
     FatalErrorIn("const labelListList& ggiGAMGInterface::weights() const")
         << "Requested fine addressing at coarse level"
         << abort(FatalError);
+
     return scalarListList::null();
 }
 
@@ -681,9 +729,8 @@ void Foam::ggiGAMGInterface::initTransfer
     const unallocLabelList& interfaceData
 ) const
 {
-    // Label transfer is straight local
+    // Label transfer is local
     labelTransferBuffer_ = interfaceData;
-//     labelTransferBuffer_ = expand(interfaceData);
 }
 
 
@@ -706,9 +753,6 @@ void Foam::ggiGAMGInterface::initInternalFieldTransfer
 {
     // Label transfer is local without global reduction
     labelTransferBuffer_ = interfaceInternalField(iF);
-
-//     labelList pif = interfaceInternalField(iF);
-//     labelTransferBuffer_ = expand(pif);
 }
 
 
@@ -730,8 +774,14 @@ void Foam::ggiGAMGInterface::initInternalFieldTransfer
 {
     scalarField pif = interfaceInternalField(iF);
 
+    // New treatment.  HJ, 26/Jun/2011
+    fieldTransferBuffer_ = fastReduce(pif);
+
+    // Old treatment
+#if(0)
     // Expand the field, executing a reduce operation in init
     fieldTransferBuffer_ = expand(pif);
+#endif
 }
 
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaces/ggiGAMGInterface/ggiGAMGInterface.H

@@ -69,6 +69,18 @@ class ggiGAMGInterface
         labelList zoneAddressing_;
 
 
+        // Parallel communication optimisation, stored on master processor
+
+            //- Has reduce been initialised?
+            mutable bool initReduce_;
+
+            //- List of zone faces indices received from each processor
+            mutable labelListList receiveAddr_;
+
+            //- List of zone faces indices to send to each processor
+            mutable labelListList sendAddr_;
+
+
     // Private Member Functions
 
         //- Disallow default bitwise copy construct
@@ -77,8 +89,19 @@ class ggiGAMGInterface
         //- Disallow default bitwise assignment
         void operator=(const ggiGAMGInterface&);
 
+
+        //- Init fast reduce
+        void initFastReduce() const;
+
+        //- Fast reduce
+        //  Note: contains global communications
+        //  New, HJ, 24/Jun/2011
+        template<class Type>
+        tmp<Field<Type> > fastReduce(const UList<Type>&) const;
+
         //- Expand data to zone size
         //  Note: contains global reduce.  HJ, 15/May/2009
+        //  Obsolete.  HJ, 24/Jun/2011
         template<class Type>
         tmp<Field<Type> > expand(const UList<Type>&) const;
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaces/ggiGAMGInterface/ggiGAMGInterfaceTemplates.C

@@ -33,6 +33,174 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class Type>
+tmp<Field<Type> > ggiGAMGInterface::fastReduce(const UList<Type>& ff) const
+{
+    // Algorithm
+    // Local processor contains faceCells part of the zone and requires
+    // zoneAddressing part.
+    // For fast communications, each processor will send the faceCells and
+    // zoneAddressing to the master.  Master will assemble global zone
+    // and send off messages to all processors containing only
+    // the required data
+    // HJ, 24/Jun/2011
+
+    if (ff.size() != this->size())
+    {
+        FatalErrorIn
+        (
+            "tmp<Field<Type> > ggiGAMGInterface::fastReduce"
+            "("
+            "    const UList<Type>& ff"
+            ") const"
+        )   << "Wrong field size.  ff: " << ff.size()
+            << " interface: " << this->size()
+            << abort(FatalError);
+    }
+
+    if (localParallel() || !Pstream::parRun())
+    {
+        // Field remains identical: no parallel communications required
+        tmp<Field<Type> > tresult(new Field<Type>(ff));
+
+        return tresult;
+    }
+
+    // Execute reduce if not already done
+    if (!initReduce_)
+    {
+        initFastReduce();
+    }
+
+    if (Pstream::master())
+    {
+        // Master collects information and distributes data.
+        Field<Type> expandField(zoneSize(), pTraits<Type>::zero);
+
+        // Insert master processor
+        const labelList& za = zoneAddressing();
+
+        forAll (za, i)
+        {
+            expandField[za[i]] = ff[i];
+        }
+
+        // Master receives and inserts data from all processors for which
+        // receiveAddr contains entries
+        for (label procI = 1; procI < Pstream::nProcs(); procI++)
+        {
+            const labelList& curRAddr = receiveAddr_[procI];
+
+            if (!curRAddr.empty())
+            {
+                Field<Type> receiveBuf(curRAddr.size());
+
+                // Opt: reconsider mode of communication
+                IPstream::read
+                (
+                    Pstream::blocking,
+                    procI,
+                    reinterpret_cast<char*>(receiveBuf.begin()),
+                    receiveBuf.byteSize()
+                );
+
+                // Insert received information
+                forAll (curRAddr, i)
+                {
+                    expandField[curRAddr[i]] = receiveBuf[i];
+                }
+            }
+        }
+
+        // Expanded field complete, send required data to other processors
+        for (label procI = 1; procI < Pstream::nProcs(); procI++)
+        {
+            const labelList& curSAddr = sendAddr_[procI];
+
+            if (!curSAddr.empty())
+            {
+                Field<Type> sendBuf(curSAddr.size());
+
+                forAll (curSAddr, i)
+                {
+                    sendBuf[i] = expandField[curSAddr[i]];
+                }
+
+                // Opt: reconsider mode of communication
+                OPstream::write
+                (
+                    Pstream::blocking,
+                    procI,
+                    reinterpret_cast<const char*>(sendBuf.begin()),
+                    sendBuf.byteSize()
+                );
+            }
+        }
+
+        // Note: different from ggi patch: field reduction happens within
+        // fastReduce.  HJ, 26/Jun/2011
+        const labelList& sza = shadowInterface().zoneAddressing();
+
+        tmp<Field<Type> > tredField
+        (
+            new Field<Type>(sza.size(), pTraits<Type>::zero)
+        );
+        Field<Type>& redField = tredField();
+
+        // Select elements from shadow zone addressing
+        forAll (sza, i)
+        {
+            redField[i] = expandField[sza[i]];
+        }
+
+        return tredField;
+    }
+    else
+    {
+        // Send local data to master and receive remote data
+        // If patch is empty, communication is avoided
+        // HJ, 4/Jun/2011
+        if (size())
+        {
+            // Opt: reconsider mode of communication
+            OPstream::write
+            (
+                Pstream::blocking,
+                Pstream::masterNo(),
+                reinterpret_cast<const char*>(ff.begin()),
+                ff.byteSize()
+            );
+        }
+
+        // Prepare to receive remote data
+        const labelList& za = zoneAddressing();
+
+        tmp<Field<Type> > treceiveBuf
+        (
+            new Field<Type>(za.size(), pTraits<Type>::zero)
+        );
+        Field<Type>& receiveBuf = treceiveBuf();
+
+        if (!za.empty())
+        {
+            // Opt: reconsider mode of communication
+            IPstream::read
+            (
+                Pstream::blocking,
+                Pstream::masterNo(),
+                reinterpret_cast<char*>(receiveBuf.begin()),
+                receiveBuf.byteSize()
+            );
+
+            // Note: different from ggi patch: field reduction happens within
+            // fastReduce.  HJ, 26/Jun/2011
+        }
+
+        return treceiveBuf;
+    }
+}
+
+
 template<class Type>
 tmp<Field<Type> > ggiGAMGInterface::expand(const UList<Type>& pd) const
 {

src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C

@@ -144,6 +144,19 @@ void Foam::polyBoundaryMesh::clearAddressing()
 
 void Foam::polyBoundaryMesh::calcGeometry()
 {
+    // Calculation of addressing, with communication
+    // HJ, 12/Jun/2011
+    forAll(*this, patchi)
+    {
+        operator[](patchi).initAddressing();
+    }
+
+    forAll(*this, patchi)
+    {
+        operator[](patchi).calcAddressing();
+    }
+
+    // Calculation of geometry with communications
     forAll(*this, patchi)
     {
         operator[](patchi).initGeometry();

src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H

@@ -78,7 +78,10 @@ class polyBoundaryMesh
         //- Create identity map
         static labelList ident(const label len);
 
-        //- Calculate the geometry for the patches (transformation tensors etc.)
+        //- Calculate the geometry for the patches,
+        //  eg transformation tensors etc.
+        //  Enhanced for parallel addressing calculations
+        //  HJ, 11/Jun/2011
         void calcGeometry();
 
         //- Disallow construct as copy

src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.H

@@ -86,6 +86,12 @@ protected:
             const scalar absTol = matchTol_
         ) const;
 
+        //- Initialise the calculation of the patch addressing
+        virtual void initAddressing() = 0;
+
+        //- Calculate the patch addressing
+        virtual void calcAddressing() = 0;
+
         //- Initialise the calculation of the patch geometry
         virtual void initGeometry() = 0;
 

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

@@ -950,6 +950,18 @@ Foam::cyclicPolyPatch::~cyclicPolyPatch()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+void Foam::cyclicPolyPatch::initAddressing()
+{
+    polyPatch::initAddressing();
+}
+
+
+void Foam::cyclicPolyPatch::calcAddressing()
+{
+    polyPatch::calcAddressing();
+}
+
+
 void Foam::cyclicPolyPatch::initGeometry()
 {
     polyPatch::initGeometry();
@@ -959,6 +971,7 @@ void Foam::cyclicPolyPatch::initGeometry()
 void Foam::cyclicPolyPatch::calcGeometry()
 {
     polyPatch::calcGeometry();
+
     calcTransforms();
 }
 

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

@@ -178,6 +178,12 @@ protected:
 
     // Protected Member functions
 
+        //- Initialise the calculation of the patch addressing
+        virtual void initAddressing();
+
+        //- Calculate the patch addressing
+        virtual void calcAddressing();
+
         //- Initialise the calculation of the patch geometry
         virtual void initGeometry();
 

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.C

@@ -51,6 +51,15 @@ namespace Foam
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+bool Foam::ggiPolyPatch::active() const
+{
+    polyPatchID shadow(shadowName_, boundaryMesh());
+    faceZoneID zone(zoneName_, boundaryMesh().mesh().faceZones());
+
+    return shadow.active() && zone.active();
+}
+
+
 void Foam::ggiPolyPatch::calcZoneAddressing() const
 {
     // Calculate patch-to-zone addressing
@@ -61,18 +70,24 @@ void Foam::ggiPolyPatch::calcZoneAddressing() const
             << abort(FatalError);
     }
 
+    if (debug)
+    {
+        Pout<< "ggiPolyPatch::calcZoneAddressing() const for patch "
+            << index() << endl;
+    }
+
     // Calculate patch-to-zone addressing
     zoneAddressingPtr_ = new labelList(size());
-    labelList& addr = *zoneAddressingPtr_;
+    labelList& zAddr = *zoneAddressingPtr_;
     const faceZone& myZone = zone();
 
     for (label i = 0; i < size(); i++)
     {
-        addr[i] = myZone.whichFace(start() + i);
+        zAddr[i] = myZone.whichFace(start() + i);
     }
 
     // Check zone addressing
-    if (addr.size() > 0 && min(addr) < 0)
+    if (zAddr.size() > 0 && min(zAddr) < 0)
     {
         FatalErrorIn("void ggiPolyPatch::calcZoneAddressing() const")
             << "Problem with patch-to zone addressing: some patch faces "
@@ -82,6 +97,85 @@ void Foam::ggiPolyPatch::calcZoneAddressing() const
 }
 
 
+void Foam::ggiPolyPatch::calcRemoteZoneAddressing() const
+{
+    // Calculate patch-to-zone addressing
+    if (remoteZoneAddressingPtr_)
+    {
+        FatalErrorIn("void ggiPolyPatch::calcRemoteZoneAddressing() const")
+            << "Patch to remote zone addressing already calculated"
+            << abort(FatalError);
+    }
+
+    if (debug)
+    {
+        Pout<< "ggiPolyPatch::calcRemoteZoneAddressing() const for patch "
+            << index() << endl;
+    }
+
+    // Once zone addressing is established, visit the opposite side and find
+    // out which face data is needed for interpolation
+    boolList usedShadows(shadow().zone().size(), false);
+
+    const labelList& zAddr = zoneAddressing();
+
+    if (master())
+    {
+        const labelListList& addr = patchToPatch().masterAddr();
+
+        forAll (zAddr, mfI)
+        {
+            const labelList& nbrs = addr[zAddr[mfI]];
+
+            forAll (nbrs, nbrI)
+            {
+                usedShadows[nbrs[nbrI]] = true;
+            }
+        }
+    }
+    else
+    {
+        const labelListList& addr = patchToPatch().slaveAddr();
+
+        forAll (zAddr, mfI)
+        {
+            const labelList& nbrs = addr[zAddr[mfI]];
+
+            forAll (nbrs, nbrI)
+            {
+                usedShadows[nbrs[nbrI]] = true;
+            }
+        }
+    }
+
+    // Count and pick up shadow indices
+    label nShadows = 0;
+
+    forAll (usedShadows, sI)
+    {
+        if (usedShadows[sI])
+        {
+            nShadows++;
+        }
+    }
+
+    remoteZoneAddressingPtr_ = new labelList(nShadows);
+    labelList& rza = *remoteZoneAddressingPtr_;
+
+    // Reset counter for re-use
+    nShadows = 0;
+
+    forAll (usedShadows, sI)
+    {
+        if (usedShadows[sI])
+        {
+            rza[nShadows] = sI;
+            nShadows++;
+        }
+    }
+}
+
+
 void Foam::ggiPolyPatch::calcPatchToPatch() const
 {
     // Create patch-to-patch interpolation
@@ -98,31 +192,29 @@ void Foam::ggiPolyPatch::calcPatchToPatch() const
         patchToPatchPtr_ =
             new ggiZoneInterpolation
             (
-                zone()(),
-                shadow().zone()(),
+                zone()(),           // This zone reference
+                shadow().zone()(),  // This shadow zone reference
                 forwardT(),
                 reverseT(),
                 shadow().separation(), // Slave-to-master separation. Bug fix
                 0,             // Non-overlapping face tolerances
                 0,             // HJ, 24/Oct/2008
                 true,          // Rescale weighting factors.  Bug fix, MB.
-                ggiInterpolation::AABB
+//                 ggiInterpolation::AABB
+                ggiInterpolation::BB_OCTREE  // Octree search, MB.
             );
 
-        // Abort immediatly if uncovered faces are present and the option
+        // Abort immediately if uncovered faces are present and the option
         // bridgeOverlap is not set.
         if
         (
             (
                 patchToPatch().uncoveredMasterFaces().size() > 0
-                &&
-                !bridgeOverlap()
+            && !bridgeOverlap()
             )
-            ||
-            (
+         || (
                 patchToPatch().uncoveredSlaveFaces().size() > 0
-                &&
-                !shadow().bridgeOverlap()
+            && !shadow().bridgeOverlap()
             )
         )
         {
@@ -144,43 +236,6 @@ void Foam::ggiPolyPatch::calcPatchToPatch() const
 }
 
 
-void Foam::ggiPolyPatch::calcLocalParallel() const
-{
-    // Calculate patch-to-zone addressing
-    if (localParallelPtr_)
-    {
-        FatalErrorIn("void ggiPolyPatch::calcLocalParallel() const")
-            << "Local parallel switch already calculated"
-            << abort(FatalError);
-    }
-
-    localParallelPtr_ = new bool(false);
-    bool& emptyOrComplete = *localParallelPtr_;
-
-    // Calculate localisation on master and shadow
-    emptyOrComplete =
-        (zone().size() == size() && shadow().zone().size() == shadow().size())
-     || (size() == 0 && shadow().size() == 0);
-
-    reduce(emptyOrComplete, andOp<bool>());
-
-    if (debug && Pstream::parRun())
-    {
-        Info<< "GGI patch Master: " << name()
-            << " Slave: " << shadowName() << " is ";
-
-        if (emptyOrComplete)
-        {
-           Info<< "local parallel" << endl;
-        }
-        else
-        {
-            Info<< "split between multiple processors" << endl;
-        }
-    }
-}
-
-
 void Foam::ggiPolyPatch::calcReconFaceCellCentres() const
 {
     // Create neighbouring face centres using interpolation
@@ -209,12 +264,162 @@ void Foam::ggiPolyPatch::calcReconFaceCellCentres() const
 }
 
 
+void Foam::ggiPolyPatch::calcLocalParallel() const
+{
+    // Calculate patch-to-zone addressing
+    if (localParallelPtr_)
+    {
+        FatalErrorIn("void ggiPolyPatch::calcLocalParallel() const")
+            << "Local parallel switch already calculated"
+            << abort(FatalError);
+    }
+
+    localParallelPtr_ = new bool(false);
+    bool& emptyOrComplete = *localParallelPtr_;
+
+    // If running in serial, all GGIs are local parallel
+    // HJ, 1/Jun/2011
+    if (!Pstream::parRun())
+    {
+        return;
+    }
+
+    // Calculate localisation on master and shadow
+    emptyOrComplete =
+        (zone().size() == size() && shadow().zone().size() == shadow().size())
+     || (size() == 0 && shadow().size() == 0);
+
+    reduce(emptyOrComplete, andOp<bool>());
+
+    if (debug && Pstream::parRun())
+    {
+        Info<< "GGI patch Master: " << name()
+            << " Slave: " << shadowName() << " is ";
+
+        if (emptyOrComplete)
+        {
+           Info<< "local parallel" << endl;
+        }
+        else
+        {
+            Info<< "split between multiple processors" << endl;
+        }
+    }
+}
+
+
+void Foam::ggiPolyPatch::calcSendReceive() const
+{
+    // Note: all processors will execute calcSendReceive but only master will
+    // hold the information.  Therefore, pointers on slave processors
+    // will remain meaningless, but for purposes of consistency
+    // (of the calc-call) they will be set to zero-sized array
+    // HJ, 4/Jun/2011
+
+    if (receiveAddrPtr_ || sendAddrPtr_)
+    {
+        FatalErrorIn("void ggiPolyPatch::calcSendReceive() const")
+            << "Send-receive addressing already calculated"
+            << abort(FatalError);
+    }
+
+    if (debug)
+    {
+        Pout<< "ggiPolyPatch::calcSendReceive() const for patch "
+            << index() << endl;
+    }
+
+    if (!Pstream::parRun())
+    {
+        FatalErrorIn("void ggiPolyPatch::calcSendReceive() const")
+            << "Requested calculation of send-receive addressing for a "
+            << "serial run.  This is not allowed"
+            << abort(FatalError);
+    }
+
+    // Master will receive and store the maps
+    if (Pstream::master())
+    {
+        receiveAddrPtr_ = new labelListList(Pstream::nProcs());
+        labelListList& rAddr = *receiveAddrPtr_;
+
+        sendAddrPtr_ = new labelListList(Pstream::nProcs());
+        labelListList& sAddr = *sendAddrPtr_;
+
+        // Insert master
+        rAddr[0] = zoneAddressing();
+
+        for (label procI = 1; procI < Pstream::nProcs(); procI++)
+        {
+            // Note: must use normal comms because the size of the
+            // communicated lists is unknown on the receiving side
+            // HJ, 4/Jun/2011
+
+            // Opt: reconsider mode of communication
+            IPstream ip(Pstream::scheduled, procI);
+
+            rAddr[procI] = labelList(ip);
+
+            sAddr[procI] = labelList(ip);
+        }
+    }
+    else
+    {
+        // Create dummy pointers: only master processor stores maps
+        receiveAddrPtr_ = new labelListList();
+        sendAddrPtr_ = new labelListList();
+
+        // Send information to master
+        const labelList& za = zoneAddressing();
+        const labelList& ra = remoteZoneAddressing();
+
+        // Note: must use normal comms because the size of the
+        // communicated lists is unknown on the receiving side
+        // HJ, 4/Jun/2011
+
+        // Opt: reconsider mode of communication
+        OPstream op(Pstream::scheduled, Pstream::masterNo());
+
+        // Send local and remote addressing to master
+        op << za << ra;
+    }
+}
+
+
+const Foam::labelListList& Foam::ggiPolyPatch::receiveAddr() const
+{
+    if (!receiveAddrPtr_)
+    {
+        calcSendReceive();
+    }
+
+    return *receiveAddrPtr_;
+}
+
+
+const Foam::labelListList& Foam::ggiPolyPatch::sendAddr() const
+{
+    if (!sendAddrPtr_)
+    {
+        calcSendReceive();
+    }
+
+    return *sendAddrPtr_;
+}
+
+
 void Foam::ggiPolyPatch::clearGeom()
 {
-    deleteDemandDrivenData(patchToPatchPtr_);
-    deleteDemandDrivenData(zoneAddressingPtr_);
-
     deleteDemandDrivenData(reconFaceCellCentresPtr_);
+
+    // Remote addressing and send-receive maps depend on the local
+    // position.  Therefore, it needs to be recalculated at mesh motion.
+    // Local zone addressing does not change with mesh motion
+    // HJ, 23/Jun/2011
+    deleteDemandDrivenData(remoteZoneAddressingPtr_);
+
+    deleteDemandDrivenData(receiveAddrPtr_);
+    deleteDemandDrivenData(sendAddrPtr_);
 }
 
 
@@ -222,6 +427,8 @@ void Foam::ggiPolyPatch::clearOut()
 {
     clearGeom();
 
+    deleteDemandDrivenData(zoneAddressingPtr_);
+    deleteDemandDrivenData(patchToPatchPtr_);
     deleteDemandDrivenData(localParallelPtr_);
 }
 
@@ -245,8 +452,11 @@ Foam::ggiPolyPatch::ggiPolyPatch
     zoneIndex_(-1),
     patchToPatchPtr_(NULL),
     zoneAddressingPtr_(NULL),
+    remoteZoneAddressingPtr_(NULL),
+    reconFaceCellCentresPtr_(NULL),
     localParallelPtr_(NULL),
-    reconFaceCellCentresPtr_(NULL)
+    receiveAddrPtr_(NULL),
+    sendAddrPtr_(NULL)
 {}
 
 
@@ -270,8 +480,11 @@ Foam::ggiPolyPatch::ggiPolyPatch
     zoneIndex_(-1),
     patchToPatchPtr_(NULL),
     zoneAddressingPtr_(NULL),
+    remoteZoneAddressingPtr_(NULL),
+    reconFaceCellCentresPtr_(NULL),
     localParallelPtr_(NULL),
-    reconFaceCellCentresPtr_(NULL)
+    receiveAddrPtr_(NULL),
+    sendAddrPtr_(NULL)
 {}
 
 
@@ -291,8 +504,11 @@ Foam::ggiPolyPatch::ggiPolyPatch
     zoneIndex_(-1),
     patchToPatchPtr_(NULL),
     zoneAddressingPtr_(NULL),
+    remoteZoneAddressingPtr_(NULL),
+    reconFaceCellCentresPtr_(NULL),
     localParallelPtr_(NULL),
-    reconFaceCellCentresPtr_(NULL)
+    receiveAddrPtr_(NULL),
+    sendAddrPtr_(NULL)
 {}
 
 
@@ -310,8 +526,11 @@ Foam::ggiPolyPatch::ggiPolyPatch
     zoneIndex_(-1),
     patchToPatchPtr_(NULL),
     zoneAddressingPtr_(NULL),
+    remoteZoneAddressingPtr_(NULL),
+    reconFaceCellCentresPtr_(NULL),
     localParallelPtr_(NULL),
-    reconFaceCellCentresPtr_(NULL)
+    receiveAddrPtr_(NULL),
+    sendAddrPtr_(NULL)
 {}
 
 
@@ -333,8 +552,11 @@ Foam::ggiPolyPatch::ggiPolyPatch
     zoneIndex_(-1),
     patchToPatchPtr_(NULL),
     zoneAddressingPtr_(NULL),
+    remoteZoneAddressingPtr_(NULL),
+    reconFaceCellCentresPtr_(NULL),
     localParallelPtr_(NULL),
-    reconFaceCellCentresPtr_(NULL)
+    receiveAddrPtr_(NULL),
+    sendAddrPtr_(NULL)
 {}
 
 
@@ -437,6 +659,17 @@ const Foam::labelList& Foam::ggiPolyPatch::zoneAddressing() const
 }
 
 
+const Foam::labelList& Foam::ggiPolyPatch::remoteZoneAddressing() const
+{
+    if (!remoteZoneAddressingPtr_)
+    {
+        calcRemoteZoneAddressing();
+    }
+
+    return *remoteZoneAddressingPtr_;
+}
+
+
 bool Foam::ggiPolyPatch::localParallel() const
 {
     // Calculate patch-to-zone addressing
@@ -483,8 +716,45 @@ const Foam::vectorField& Foam::ggiPolyPatch::reconFaceCellCentres() const
 }
 
 
+void Foam::ggiPolyPatch::initAddressing()
+{
+    if (active())
+    {
+        // Force zone addressing first
+        zoneAddressing();
+        remoteZoneAddressing();
+
+        if (Pstream::parRun() && !localParallel())
+        {
+            sendAddr();
+        }
+    }
+
+    polyPatch::initAddressing();
+}
+
+
+void Foam::ggiPolyPatch::calcAddressing()
+{
+    polyPatch::calcAddressing();
+}
+
+
 void Foam::ggiPolyPatch::initGeometry()
 {
+    // Communication is allowed either before or after processor
+    // patch comms.  HJ, 11/Jul/2011
+    if (active())
+    {
+        calcTransforms();
+
+        // Note: Only master calculates recon; slave uses master interpolation
+        if (master())
+        {
+            reconFaceCellCentres();
+        }
+    }
+
     polyPatch::initGeometry();
 }
 
@@ -497,18 +767,40 @@ void Foam::ggiPolyPatch::calcGeometry()
     // reconFaceCellCentres in order to correctly set the transformation
     // in the interpolation routines.
     // HJ, 3/Jul/2009
-    calcTransforms();
-
-    // Reconstruct the cell face centres
-    if (patchToPatchPtr_ && master())
-    {
-        reconFaceCellCentres();
-    }
 }
 
 
 void Foam::ggiPolyPatch::initMovePoints(const pointField& p)
 {
+    clearGeom();
+
+    // Calculate transforms on mesh motion?
+    calcTransforms();
+
+    // Update interpolation for new relative position of GGI interfaces
+    if (patchToPatchPtr_)
+    {
+        patchToPatchPtr_->movePoints();
+    }
+
+    // Recalculate send and receive maps
+    if (active())
+    {
+        // Force zone addressing first
+        zoneAddressing();
+        remoteZoneAddressing();
+
+        if (Pstream::parRun() && !localParallel())
+        {
+            sendAddr();
+        }
+    }
+
+    if (active() && master())
+    {
+        reconFaceCellCentres();
+    }
+
     polyPatch::initMovePoints(p);
 }
 
@@ -516,7 +808,6 @@ void Foam::ggiPolyPatch::initMovePoints(const pointField& p)
 void Foam::ggiPolyPatch::movePoints(const pointField& p)
 {
     polyPatch::movePoints(p);
-    clearGeom();
 }
 
 

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatch.H

@@ -85,25 +85,43 @@ class ggiPolyPatch
         //- Patch-to-zone addressing
         mutable labelList* zoneAddressingPtr_;
 
-        //- Is the patch localised on a single processor
-        // (single processor in a parallel run)?
-        //  Used for parallel optimisation
-        mutable bool* localParallelPtr_;
+        //- Remote zone addressing: data needed for interpolation
+        mutable labelList* remoteZoneAddressingPtr_;
 
         //- Reconstructed patch neighbour cell centres
         mutable vectorField* reconFaceCellCentresPtr_;
 
 
+        // Parallel communication optimisation, stored on master processor
+
+            //- Is the patch localised on a single processor
+            // (single processor in a parallel run)?
+            //  Used for parallel optimisation
+            mutable bool* localParallelPtr_;
+
+            //- List of zone faces indices received from each processor
+            mutable labelListList* receiveAddrPtr_;
+
+            //- List of zone faces indices to send to each processor
+            mutable labelListList* sendAddrPtr_;
+
+
     // Private member functions
 
+        //- Is the GGI active& (zone and shadow present)
+        bool active() const;
+
         //- Calculate patch-to-zone addressing
         virtual void calcZoneAddressing() const;
 
+        //- Calculate remote patch-to-zone addressing
+        virtual void calcRemoteZoneAddressing() const;
+
         //- Calculate interpolation
         virtual void calcPatchToPatch() const;
 
-        //- Calculate local parallel switch
-        void calcLocalParallel() const;
+
+        // Geometry
 
             //- Calculate reconstructed cell centres
             void calcReconFaceCellCentres() const;
@@ -111,6 +129,25 @@ class ggiPolyPatch
             //- Force calculation of transformation tensors
             virtual void calcTransforms();
 
+
+        // Parallel communication optimisation, stored on master processor
+
+            //- Calculate local parallel switch
+            void calcLocalParallel() const;
+
+            //- Calculate send and receive addressing
+            void calcSendReceive() const;
+
+
+            //- Return receive addressing
+            const labelListList& receiveAddr() const;
+
+            //- Return send addressing
+            const labelListList& sendAddr() const;
+
+
+        // Memory management
+
             //- Clear geometry
              void clearGeom();
 
@@ -122,6 +159,12 @@ protected:
 
     // Protected Member functions
 
+        //- Initialise the calculation of the patch addressing
+        virtual void initAddressing();
+
+        //- Calculate the patch addressing
+        virtual void calcAddressing();
+
         //- Initialise the calculation of the patch geometry
         virtual void initGeometry();
 
@@ -281,6 +324,9 @@ public:
         //- Return patch-to-zone addressing
         const labelList& zoneAddressing() const;
 
+        //- Return remote patch-to-zone addressing
+        const labelList& remoteZoneAddressing() const;
+
         //- Is the patch localised on a single processor
         bool localParallel() const;
 
@@ -289,13 +335,21 @@ public:
         const ggiZoneInterpolation& patchToPatch() const;
 
         //- Expand face field to full zone size, including reduction
+        //  New.  HJ, 12/Jun/2011
+        template<class Type>
+        tmp<Field<Type> > fastExpand(const Field<Type>& pf) const;
+
+        //- Expand face field to full zone size, including reduction
+        //  Obsolete.  HJ, 12/Jun/2011
         template<class Type>
         tmp<Field<Type> > expand(const Field<Type>& pf) const;
 
         //- Filter zone field to patch size
+        //  Obsolete.  HJ, 12/Jun/2011
         template<class Type>
         tmp<Field<Type> > filter(const Field<Type>& zf) const;
 
+
         //- Interpolate face field: given field on a the shadow side,
         //  create an interpolated field on this side
         template<class Type>

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/ggi/ggiPolyPatchTemplates.C

@@ -28,26 +28,193 @@ Author
 \*---------------------------------------------------------------------------*/
 
 #include "ggiPolyPatch.H"
+#include "OSspecific.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+// New.  HJ, 12/Jun/2011
 template<class Type>
-Foam::tmp<Foam::Field<Type> > Foam::ggiPolyPatch::expand
+Foam::tmp<Foam::Field<Type> > Foam::ggiPolyPatch::fastExpand
 (
-    const Field<Type>& pf
+    const Field<Type>& ff
 ) const
 {
     // Check and expand the field from patch size to zone size
-    if (pf.size() != size())
+    // with communication
+
+    // Algorithm:
+    // 1) Master processor holds maps of all zone addressing (data provided)
+    // and all remote zone addressing (data required)
+    // 2) Each processor will send the locally active data to the master
+    // 3) Master assembles all the data
+    // 4) Master sends to all processors the data they need to receive
+    //
+    // Notes:
+    // A) If the size of zone addressing is zero, data is not sent
+    // B) Communicated data on each processor has the size of live faces
+    // C) Expanded data will be equal to actual data fronm other processors
+    //    only for the faces marked in remote; for other faces, it will be
+    //    equal to zero
+    // D) On processor zero, complete data is available
+    // HJ, 4/Jun/2011
+    if (ff.size() != size())
+    {
+        FatalErrorIn
+        (
+            "tmp<Field<Type> > ggiPolyPatch::fastExpand"
+            "("
+            "    const Field<Type>& ff"
+            ") const"
+        )   << "Incorrect patch field size.  Field size: "
+            << ff.size() << " patch size: " << size()
+            << abort(FatalError);
+    }
+
+    if (localParallel())
+    {
+        FatalErrorIn
+        (
+            "tmp<Field<Type> > ggiPolyPatch::fastExpand"
+            "("
+            "    const Field<Type>& ff"
+            ") const"
+        )   << "Requested expand on local parallel.  This is not allowed"
+            << abort(FatalError);
+    }
+
+    // Expand the field to zone size
+    tmp<Field<Type> > texpandField
+    (
+        new Field<Type>(zone().size(), pTraits<Type>::zero)
+    );
+
+    Field<Type>& expandField = texpandField();
+
+    if (Pstream::master())
+    {
+        // Insert master processor
+        const labelList& za = zoneAddressing();
+
+        forAll (za, i)
+        {
+            expandField[za[i]] = ff[i];
+        }
+
+        // Master receives and inserts data from all processors for which
+        // receiveAddr contains entries
+        for (label procI = 1; procI < Pstream::nProcs(); procI++)
+        {
+            const labelList& curRAddr = receiveAddr()[procI];
+
+            if (!curRAddr.empty())
+            {
+                Field<Type> receiveBuf(curRAddr.size());
+
+                // Opt: reconsider mode of communication
+                IPstream::read
+                (
+                    Pstream::blocking,
+                    procI,
+                    reinterpret_cast<char*>(receiveBuf.begin()),
+                    receiveBuf.byteSize()
+                );
+
+                // Insert received information
+                forAll (curRAddr, i)
+                {
+                    expandField[curRAddr[i]] = receiveBuf[i];
+                }
+            }
+        }
+
+        // Expanded field complete, send required data to other processors
+        for (label procI = 1; procI < Pstream::nProcs(); procI++)
+        {
+            const labelList& curSAddr = sendAddr()[procI];
+
+            if (!curSAddr.empty())
+            {
+                Field<Type> sendBuf(curSAddr.size());
+
+                forAll (curSAddr, i)
+                {
+                    sendBuf[i] = expandField[curSAddr[i]];
+                }
+
+                // Opt: reconsider mode of communication
+                OPstream::write
+                (
+                    Pstream::blocking,
+                    procI,
+                    reinterpret_cast<const char*>(sendBuf.begin()),
+                    sendBuf.byteSize()
+                );
+            }
+        }
+    }
+    else
+    {
+        // Send local data to master and receive remote data
+        // If patch is empty, communication is avoided
+        // HJ, 4/Jun/2011
+        if (size())
+        {
+            // Opt: reconsider mode of communication
+            OPstream::write
+            (
+                Pstream::blocking,
+                Pstream::masterNo(),
+                reinterpret_cast<const char*>(ff.begin()),
+                ff.byteSize()
+            );
+        }
+
+        // Prepare to receive remote data
+        const labelList& rza = remoteZoneAddressing();
+
+        if (!rza.empty())
+        {
+            Field<Type> receiveBuf(rza.size());
+
+            // Opt: reconsider mode of communication
+            IPstream::read
+            (
+                Pstream::blocking,
+                Pstream::masterNo(),
+                reinterpret_cast<char*>(receiveBuf.begin()),
+                receiveBuf.byteSize()
+            );
+
+            // Insert the data into expanded field
+            forAll (rza, i)
+            {
+                expandField[rza[i]] = receiveBuf[i];
+            }
+        }
+    }
+
+    return texpandField;
+}
+
+
+// Obsolete.  HJ, 12/Jun/2011
+template<class Type>
+Foam::tmp<Foam::Field<Type> > Foam::ggiPolyPatch::expand
+(
+    const Field<Type>& ff
+) const
+{
+    // Check and expand the field from patch size to zone size
+    if (ff.size() != size())
     {
         FatalErrorIn
         (
             "tmp<Field<Type> > ggiPolyPatch::expand"
             "("
-            "    const Field<Type>& pf"
+            "    const Field<Type>& ff"
             ") const"
         )   << "Incorrect patch field size.  Field size: "
-            << pf.size() << " patch size: " << size()
+            << ff.size() << " patch size: " << size()
             << abort(FatalError);
     }
 
@@ -63,7 +230,7 @@ Foam::tmp<Foam::Field<Type> > Foam::ggiPolyPatch::expand
 
     forAll (addr, i)
     {
-        expandField[addr[i]] = pf[i];
+        expandField[addr[i]] = ff[i];
     }
 
     // Parallel data exchange: update surface field on all processors
@@ -79,6 +246,7 @@ Foam::tmp<Foam::Field<Type> > Foam::ggiPolyPatch::expand
 }
 
 
+// Obsolete.  HJ, 12/Jun/2011
 template<class Type>
 Foam::tmp<Foam::Field<Type> > Foam::ggiPolyPatch::filter
 (
@@ -121,25 +289,82 @@ Foam::tmp<Foam::Field<Type> > Foam::ggiPolyPatch::filter
 template<class Type>
 Foam::tmp<Foam::Field<Type> > Foam::ggiPolyPatch::interpolate
 (
-    const Field<Type>& pf
+    const Field<Type>& ff
 ) const
 {
     // Check and expand the field from patch size to zone size
-    if (pf.size() != shadow().size())
+    if (ff.size() != shadow().size())
     {
         FatalErrorIn
         (
             "tmp<Field<Type> > ggiPolyPatch::interpolate"
             "("
-            "    const Field<Type>& pf"
+            "    const Field<Type>& ff"
             ") const"
         )   << "Incorrect slave patch field size.  Field size: "
-            << pf.size() << " patch size: " << shadow().size()
+            << ff.size() << " patch size: " << shadow().size()
             << abort(FatalError);
     }
 
+#   if 1
+
+    // New.  HJ, 12/Jun/2011
+    if (localParallel())
+    {
+        // No expansion or filtering needed.  HJ, 4/Jun/2011
+
+        // Interpolate field
+        if (master())
+        {
+            return patchToPatch().slaveToMaster(ff);
+        }
+        else
+        {
+            return patchToPatch().masterToSlave(ff);
+        }
+    }
+    else
+    {
+        // Note: fast expand is always done on the local side
+        // HJ, 24/Jun/2011
+        Field<Type> expandField = fastExpand(ff);
+
+        tmp<Field<Type> > tresult(new Field<Type>(size()));
+        Field<Type>& result = tresult();
+
+        if (master())
+        {
+            patchToPatch().maskedSlaveToMaster
+            (
+                expandField,
+                result,
+                zoneAddressing()
+            );
+        }
+        else
+        {
+            patchToPatch().maskedMasterToSlave
+            (
+                expandField,
+                result,
+                zoneAddressing()
+            );
+        }
+
+        return tresult;
+    }
+
+#   else
+
+    // Old treatment
+    // Obsolete.  HJ, 12/Jun/2011
+
     // Expand the field to zone size
-    Field<Type> expandField = shadow().expand(pf);
+    // Note: with full fields it is the shadow side that does
+    // the expand.  This is different than fastExpand because
+    // the addressing is stored remotely.
+    // HJ, 24/Jun/2011
+    Field<Type> expandField = shadow().expand(ff);
 
     Field<Type> zoneField;
 
@@ -154,17 +379,19 @@ Foam::tmp<Foam::Field<Type> > Foam::ggiPolyPatch::interpolate
     }
 
     return this->filter(zoneField);
+
+#   endif
 }
 
 
 template<class Type>
 Foam::tmp<Foam::Field<Type> > Foam::ggiPolyPatch::interpolate
 (
-    const tmp<Field<Type> >& tpf
+    const tmp<Field<Type> >& tff
 ) const
 {
-    tmp<Field<Type> > tint = interpolate(tpf());
-    tpf.clear();
+    tmp<Field<Type> > tint = interpolate(tff());
+    tff.clear();
     return tint;
 }
 
@@ -178,7 +405,43 @@ void Foam::ggiPolyPatch::bridge
 {
     if (bridgeOverlap())
     {
-        // Parallelisation
+#       ifdef FAST_PARALLEL_GGI
+
+        if (localParallel())
+        {
+            if (master())
+            {
+                patchToPatch().bridgeMaster(bridgeField, ff);
+            }
+            else
+            {
+                patchToPatch().bridgeSlave(bridgeField, ff);
+            }
+        }
+        else
+        {
+            // Note: since bridging is only a local operation
+            if (master())
+            {
+                patchToPatch().maskedBridgeMaster
+                (
+                    bridgeField,
+                    ff,
+                    zoneAddressing()
+                );
+            }
+            else
+            {
+                patchToPatch().maskedBridgeSlave
+                (
+                    bridgeField,
+                    ff,
+                    zoneAddressing()
+                );
+            }
+        }
+
+#       else
 
         // Expand the field to zone size
         Field<Type> expandBridge = this->expand(bridgeField);
@@ -200,15 +463,8 @@ void Foam::ggiPolyPatch::bridge
         {
             ff[i] = expandField[addr[i]];
         }
-    }
-    else
-    {
-        // HJ, temporary
-        InfoIn
-        (
-            "void bridge(const Field<Type>& bridgeField, "
-            "Field<Type>& ff) const"
-        )   << "Bridging is switched off " << endl;
+
+#       endif
     }
 }
 

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/overlapGgi/overlapGgiPolyPatch.H

@@ -147,6 +147,12 @@ protected:
 
     // Protected Member functions
 
+        //- Initialise the calculation of the patch addressing
+        virtual void initAddressing();
+
+        //- Calculate the patch addressing
+        virtual void calcAddressing();
+
         //- Initialise the calculation of the patch geometry
         virtual void initGeometry();
 

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/overlapGgi/overlapGgiPolyPatchGeometry.C

@@ -175,8 +175,7 @@ void Foam::overlapGgiPolyPatch::calcPatchToPatch() const
         if
         (
             patchToPatch().uncoveredMasterFaces().size() > 0
-            ||
-            patchToPatch().uncoveredSlaveFaces().size() > 0
+         || patchToPatch().uncoveredSlaveFaces().size() > 0
         )
         {
             FatalErrorIn("void overlapGgiPolyPatch::calcPatchToPatch() const")
@@ -309,6 +308,18 @@ Foam::overlapGgiPolyPatch::reconFaceCellCentres() const
 }
 
 
+void Foam::overlapGgiPolyPatch::initAddressing()
+{
+    polyPatch::initAddressing();
+}
+
+
+void Foam::overlapGgiPolyPatch::calcAddressing()
+{
+    polyPatch::calcAddressing();
+}
+
+
 void Foam::overlapGgiPolyPatch::initGeometry()
 {
     polyPatch::initGeometry();

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

@@ -135,6 +135,18 @@ Foam::processorPolyPatch::~processorPolyPatch()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+void Foam::processorPolyPatch::initAddressing()
+{
+    polyPatch::initAddressing();
+}
+
+
+void Foam::processorPolyPatch::calcAddressing()
+{
+    polyPatch::calcAddressing();
+}
+
+
 void Foam::processorPolyPatch::initGeometry()
 {
     if (Pstream::parRun())

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

@@ -92,6 +92,12 @@ protected:
 
     // Protected Member functions
 
+        //- Initialise the calculation of the patch addressing
+        virtual void initAddressing();
+
+        //- Calculate the patch addressing
+        virtual void calcAddressing();
+
         //- Initialise the calculation of the patch geometry
         void initGeometry();
 

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/regionCouple/regionCouplePolyPatch.C

@@ -79,7 +79,8 @@ void Foam::regionCouplePolyPatch::calcInterpolation() const
     {
         FatalErrorIn("void regionCouplePolyPatch::calcInterpolation() const")
             << "Shadow of regionCouple patch " << name()
-            << " named " << shadowPatchName() << " is not a regionCouple." << nl
+            << " named " << shadowPatchName()
+            << " is not a regionCouple." << nl
             << "This is not allowed.  Please check your mesh definition."
             << abort(FatalError);
     }
@@ -153,10 +154,17 @@ void Foam::regionCouplePolyPatch::calcReconFaceCellCentres() const
 }
 
 
+void Foam::regionCouplePolyPatch::clearGeom() const
+{
+    deleteDemandDrivenData(reconFaceCellCentresPtr_);
+}
+
+
 void Foam::regionCouplePolyPatch::clearOut() const
 {
+    clearGeom();
+
     deleteDemandDrivenData(patchToPatchPtr_);
-    deleteDemandDrivenData(reconFaceCellCentresPtr_);
 }
 
 
@@ -173,14 +181,14 @@ Foam::regionCouplePolyPatch::regionCouplePolyPatch
     const word& shadowRegionName,
     const word& shadowPatchName,
     const bool attached,
-    const bool attachedWalls
+    const bool isWall
 )
 :
     coupledPolyPatch(name, size, start, index, bm),
     shadowRegionName_(shadowRegionName),
     shadowPatchName_(shadowPatchName),
     attached_(attached),
-    attachedWalls_(attachedWalls),
+    isWall_(isWall),
     shadowIndex_(-1),
     patchToPatchPtr_(NULL),
     reconFaceCellCentresPtr_(NULL)
@@ -200,7 +208,7 @@ Foam::regionCouplePolyPatch::regionCouplePolyPatch
     shadowRegionName_(dict.lookup("shadowRegion")),
     shadowPatchName_(dict.lookup("shadowPatch")),
     attached_(dict.lookup("attached")),
-    attachedWalls_(dict.lookup("attachedWalls")),
+    isWall_(dict.lookup("isWall")),
     shadowIndex_(-1),
     patchToPatchPtr_(NULL),
     reconFaceCellCentresPtr_(NULL)
@@ -218,7 +226,7 @@ Foam::regionCouplePolyPatch::regionCouplePolyPatch
     shadowRegionName_(pp.shadowRegionName_),
     shadowPatchName_(pp.shadowPatchName_),
     attached_(pp.attached_),
-    attachedWalls_(pp.attachedWalls_),
+    isWall_(pp.isWall_),
     shadowIndex_(-1),
     patchToPatchPtr_(NULL),
     reconFaceCellCentresPtr_(NULL)
@@ -239,7 +247,7 @@ Foam::regionCouplePolyPatch::regionCouplePolyPatch
     shadowRegionName_(pp.shadowRegionName_),
     shadowPatchName_(pp.shadowPatchName_),
     attached_(pp.attached_),
-    attachedWalls_(pp.attachedWalls_),
+    isWall_(pp.isWall_),
     shadowIndex_(-1),
     patchToPatchPtr_(NULL),
     reconFaceCellCentresPtr_(NULL)
@@ -278,7 +286,11 @@ void Foam::regionCouplePolyPatch::attach() const
     {
         attached_ = true;
         shadow().attach();
-        clearOut();
+
+        // Patch-to-patch interpolation does not need to be cleared,
+        // only face/cell centres and interpolation factors
+        // HJ, 6/Jun/2011
+        clearGeom();
     }
 }
 
@@ -289,7 +301,11 @@ void Foam::regionCouplePolyPatch::detach() const
     {
         attached_ = false;
         shadow().detach();
-        clearOut();
+
+        // Patch-to-patch interpolation does not need to be cleared,
+        // only face/cell centres and interpolation factors
+        // HJ, 6/Jun/2011
+        clearGeom();
     }
 }
 
@@ -336,6 +352,18 @@ Foam::regionCouplePolyPatch::reconFaceCellCentres() const
 }
 
 
+void Foam::regionCouplePolyPatch::initAddressing()
+{
+    polyPatch::initAddressing();
+}
+
+
+void Foam::regionCouplePolyPatch::calcAddressing()
+{
+    polyPatch::calcAddressing();
+}
+
+
 void Foam::regionCouplePolyPatch::initGeometry()
 {
     polyPatch::initGeometry();
@@ -344,6 +372,7 @@ void Foam::regionCouplePolyPatch::initGeometry()
 
 void Foam::regionCouplePolyPatch::calcGeometry()
 {
+    polyPatch::calcGeometry();
     // Reconstruct the cell face centres
 //     reconFaceCellCentres();
 }
@@ -432,6 +461,8 @@ void Foam::regionCouplePolyPatch::write(Ostream& os) const
         << shadowPatchName_ << token::END_STATEMENT << nl;
     os.writeKeyword("attached")
         << attached_ << token::END_STATEMENT << nl;
+    os.writeKeyword("isWall")
+        << isWall_ << token::END_STATEMENT << nl;
 }
 
 

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/regionCouple/regionCouplePolyPatch.H

@@ -70,7 +70,7 @@ class regionCouplePolyPatch
         mutable Switch attached_;
 
         //- Are the region attached walls
-        mutable Switch attachedWalls_;
+        Switch isWall_;
 
         //- Shadow patch index.  Delayed evaluation for construction
         mutable label shadowIndex_;
@@ -93,6 +93,12 @@ protected:
 
     // Protected Member functions
 
+        //- Initialise the calculation of the patch addressing
+        virtual void initAddressing();
+
+        //- Calculate the patch addressing
+        virtual void calcAddressing();
+
         //- Initialise the calculation of the patch geometry
         virtual void initGeometry();
 
@@ -121,9 +127,13 @@ protected:
         //- Calculate reconstructed cell centres
         void calcReconFaceCellCentres() const;
 
+        //- Clear geometry
+        void clearGeom() const;
+
         //- Clear out
         void clearOut() const;
 
+
 public:
 
     //- Runtime type information
@@ -237,10 +247,10 @@ public:
             return attached_;
         }
 
-        //AJ: read from dictionary if coupled patch is also a wall
+        // Return true if patch is considered a wall
         bool isWall() const
         {
-            return attachedWalls_;
+            return isWall_;
         }
 
         //- Attach regions

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

@@ -115,6 +115,14 @@ protected:
         // The polyPatch geometry initialisation is called by polyBoundaryMesh
         friend class polyBoundaryMesh;
 
+        //- Initialise the calculation of the patch addressing
+        virtual void initAddressing()
+        {}
+
+        //- Calculate the patch addressing
+        virtual void calcAddressing()
+        {}
+
         //- Initialise the calculation of the patch geometry
         virtual void initGeometry()
         {}

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

@@ -119,7 +119,7 @@ public:
         const treeBoundBox& cubeBb
     );
 
-    //- Does triangle intersect plane. Return bool and set intersection segment.
+    //- Does triangle intersect plane. Return bool and set intersection segment
     static bool intersect
     (
         const point& va0,
@@ -133,7 +133,7 @@ public:
         point& pInter1
     );
 
-    //- Do triangles intersect. Return bool and set intersection segment.
+    //- Do triangles intersect. Return bool and set intersection segment
     static bool intersect
     (
         const point& va0,

src/dynamicMesh/dynamicFvMesh/mixerGgiFvMesh/mixerGgiFvMesh.C

@@ -203,7 +203,7 @@ void Foam::mixerGgiFvMesh::calcMovingMasks() const
     forAll (staticPatches, patchI)
     {
         const label staticSliderID =
-            boundaryMesh().findPatchID(movingPatches[patchI]);
+            boundaryMesh().findPatchID(staticPatches[patchI]);
 
         if (staticSliderID < 0)
         {

src/dynamicMesh/dynamicMesh/polyMeshModifiers/slidingInterface/coupleSlidingInterface.C

@@ -1764,7 +1764,9 @@ void Foam::slidingInterface::coupleInterface(polyTopoChange& ref) const
                     "polyTopoChange& ref) const"
                 )   << "Face " << curFaceID << " reduced to less than "
                     << "3 points.  Topological/cutting error B." << nl
-                    << "Old face: " << oldFace << " new face: " << newFaceLabels
+                    << "Old face: " << oldFace << " new face: "
+                    << newFaceLabels << nl
+                    << "old points: " << oldFace.points(points)
                     << abort(FatalError);
             }
 

src/finiteVolume/cfdTools/general/SRF/derivedFvPatchFields/SRFVelocityFvPatchVectorField/SRFVelocityFvPatchVectorField.C

@@ -44,7 +44,7 @@ SRFVelocityFvPatchVectorField::SRFVelocityFvPatchVectorField
 )
 :
     fixedValueFvPatchVectorField(p, iF),
-    relative_(0),
+    relative_(false),
     inletValue_(p.size(), vector::zero)
 {}
 

src/finiteVolume/cfdTools/incompressible/movingMeshContinuityErrs.H

@@ -30,6 +30,7 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
+if (mesh.moving())
 {
     volScalarField contErr = fvc::div(phi + fvc::meshPhi(U));
 
@@ -47,5 +48,10 @@ Description
         << ", cumulative = " << cumulativeContErr
         << endl;
 }
+else
+{
+#   include "continuityErrs.H"
+}
+
 
 // ************************************************************************* //

src/finiteVolume/cfdTools/incompressible/movingMeshRhoUContinuityErrs.H

@@ -53,4 +53,5 @@ else
 #   include "continuityErrs.H"
 }
 
+
 // ************************************************************************* //

src/finiteVolume/fields/fvPatchFields/derived/totalPressure/totalPressureFvPatchScalarField.C

@@ -30,7 +30,6 @@ License
 #include "volFields.H"
 #include "surfaceFields.H"
 
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::totalPressureFvPatchScalarField::totalPressureFvPatchScalarField

src/meshTools/Make/files

@@ -16,31 +16,11 @@ primitiveMeshGeometry/primitiveMeshGeometry.C
 
 meshSearch/meshSearch.C
 
-meshTools/meshTools.C
-
 PointEdgeWave/PointEdgeWaveName.C
 PointEdgeWave/pointEdgePoint.C
 
 regionSplit/regionSplit.C
 
-octree/octreeName.C
-octree/octreeDataPoint.C
-octree/octreeDataPointTreeLeaf.C
-octree/octreeDataEdges.C
-octree/octreeDataCell.C
-octree/octreeDataFace.C
-octree/treeBoundBox.C
-octree/treeNodeName.C
-octree/treeLeafName.C
-octree/pointIndexHitIOList.C
-
-indexedOctree/indexedOctreeName.C
-indexedOctree/treeDataCell.C
-indexedOctree/treeDataEdge.C
-indexedOctree/treeDataFace.C
-indexedOctree/treeDataPoint.C
-indexedOctree/treeDataTriSurface.C
-
 searchableSurface = searchableSurface
 $(searchableSurface)/distributedTriSurfaceMesh.C
 $(searchableSurface)/searchableBox.C
@@ -127,7 +107,6 @@ $(intersectedSurface)/edgeSurface.C
 triSurface/triSurfaceSearch/triSurfaceSearch.C
 triSurface/octreeData/octreeDataTriSurface.C
 triSurface/octreeData/octreeDataTriSurfaceTreeLeaf.C
-triSurface/triangleFuncs/triangleFuncs.C
 triSurface/surfaceFeatures/surfaceFeatures.C
 triSurface/triSurfaceTools/triSurfaceTools.C
 triSurface/triSurfaceTools/geompack/geompack.C

src/meshTools/searchableSurface/triSurfaceMesh.H

@@ -46,6 +46,7 @@ SourceFiles
 #include "objectRegistry.H"
 #include "indexedOctree.H"
 #include "treeDataTriSurface.H"
+#include "octreeDataTriSurfaceTreeLeaf.H"
 #include "treeDataEdge.H"
 #include "EdgeMap.H"
 

src/meshTools/sets/cellSources/rotatedBoxToCell/rotatedBoxToCell.C

@@ -129,7 +129,7 @@ Foam::rotatedBoxToCell::rotatedBoxToCell
 )
 :
     topoSetSource(mesh),
-    origin_(origin_),
+    origin_(origin),
     i_(i),
     j_(j),
     k_(k)

src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.H

@@ -73,6 +73,7 @@ class triSurfaceTools
                 RED,
                 GREEN
             };
+
             static void calcRefineStatus
             (
                 const triSurface& surf,
@@ -118,8 +119,8 @@ class triSurfaceTools
                 label edgeI
             );
 
-            // Return value of faceUsed for faces using vertI (local numbering).
-            // Used internally.
+            // Return value of faceUsed for faces using vertI (local numbering)
+            // Used internally
             static label vertexUsesFace
             (
                 const triSurface& surf,

src/tetDecompositionFiniteElement/Make.cellDecomp/options

@@ -1,2 +1,6 @@
 EXE_INC = \
+    -I$(LIB_SRC)/meshTools/lnInclude \
     -DCELL_DECOMP
+
+EXE_LIBS = \
+    -lmeshTools

src/tetDecompositionFiniteElement/Make.faceDecomp/options

@@ -1 +1,6 @@
-EXE_INC = -DFACE_DECOMP
+EXE_INC = \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -DFACE_DECOMP
+
+EXE_LIBS = \
+    -lmeshTools

src/triSurface/Make/files

@@ -7,6 +7,7 @@ meshTriangulation/meshTriangulation.C
 triSurface/triSurface.C
 triSurface/triSurfaceAddressing.C
 triSurface/stitchTriangles.C
+triSurface/treeDataTriSurface/treeDataTriSurface.C
 
 interfaces = triSurface/interfaces
 $(interfaces)/STL/writeSTL.C

src/triSurface/Make/options

@@ -1,3 +1,4 @@
 EXE_INC = \
+    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(WM_THIRD_PARTY_DIR)/zlib-1.2.3
 

tutorials/basic/scalarTransportFoam/swirlTest/constant/polyMesh/boundary

@@ -1,6 +1,6 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM Extend Project: Open Source CFD        |
+| \\      /  F ield         | OpenFOAM Extend Project: Open source CFD        |
 |  \\    /   O peration     | Version:  1.6-ext                               |
 |   \\  /    A nd           | Web:      www.extend-project.de                 |
 |    \\/     M anipulation  |                                                 |
@@ -10,6 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       polyBoundaryMesh;
+    location    "constant/polyMesh";
     object      boundary;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -22,21 +23,18 @@ inlet
         nFaces          20;
         startFace       1540;
     }
-
     outlet
     {
         type            patch;
         nFaces          20;
         startFace       1560;
     }
-
     fixedWalls
     {
         type            wall;
         nFaces          80;
         startFace       1580;
     }
-
     defaultFaces
     {
         type            empty;

tutorials/coupled/conjugateHeatFoam/conjugateCavity/boundary

@@ -39,6 +39,7 @@ right
     shadowRegion    solid;
     shadowPatch     left;
     attached        on;
+    isWall          yes;
 }
 
 bottom

tutorials/coupled/conjugateHeatFoam/heatedBlock/boundary

@@ -32,6 +32,7 @@ left
     shadowRegion    region0;
     shadowPatch     right;
     attached        on;
+    isWall          yes;  // Not relevant
 }
 
 right