Initial mattijsLib merge

src/Allwmake

 #!/bin/sh
 set -x
 
-( cd $FOAM_SRC/other && ./Allwmake )
+#( cd $FOAM_SRC/other && ./Allwmake )
 
 ( cd $FOAM_SRC/OpenFOAM && wmakeLnInclude . )
 

src/OpenFOAM/Make/files

@@ -400,7 +400,8 @@ mapPolyMesh = $(polyMesh)/mapPolyMesh
 $(mapPolyMesh)/mapPolyMesh.C
 $(mapPolyMesh)/faceMapper/faceMapper.C
 $(mapPolyMesh)/cellMapper/cellMapper.C
-$(mapPolyMesh)/mapDistributePolyMesh/mapDistributePolyMesh.C
+$(mapPolyMesh)/mapDistribute/mapDistribute.C
+$(mapPolyMesh)/mapDistribute/mapDistributePolyMesh.C
 $(mapPolyMesh)/mapAddedPolyMesh.C
 
 PrimitivePatch = $(primitiveMesh)/PrimitivePatch

src/OpenFOAM/algorithms/MeshWave/FaceCellWave.C

@@ -527,7 +527,7 @@ void Foam::FaceCellWave<Type>::handleProcPatches()
     {
         const polyPatch& patch = mesh_.boundaryMesh()[patchI];
 
-        if (Pstream::parRun() && isA<processorPolyPatch>(patch))
+        if (isA<processorPolyPatch>(patch))
         {
             // Allocate buffers
             label nSendFaces;
@@ -580,7 +580,7 @@ void Foam::FaceCellWave<Type>::handleProcPatches()
     {
         const polyPatch& patch = mesh_.boundaryMesh()[patchI];
 
-        if (Pstream::parRun() && isA<processorPolyPatch>(patch))
+        if (isA<processorPolyPatch>(patch))
         {
             const processorPolyPatch& procPatch =
                 refCast<const processorPolyPatch>(patch);
@@ -812,7 +812,6 @@ Foam::FaceCellWave<Type>::FaceCellWave
     changedCells_(mesh_.nCells()),
     nChangedCells_(0),
     hasCyclicPatches_(hasPatchType(cyclicPolyPatch::typeName)),
-    hasProcPatches_(hasPatchType(processorPolyPatch::typeName)),
     nEvals_(0),
     nUnvisitedCells_(mesh_.nCells()),
     nUnvisitedFaces_(mesh_.nFaces()),
@@ -843,7 +842,6 @@ Foam::FaceCellWave<Type>::FaceCellWave
     changedCells_(mesh_.nCells()),
     nChangedCells_(0),
     hasCyclicPatches_(hasPatchType(cyclicPolyPatch::typeName)),
-    hasProcPatches_(hasPatchType(processorPolyPatch::typeName)),
     nEvals_(0),
     nUnvisitedCells_(mesh_.nCells()),
     nUnvisitedFaces_(mesh_.nFaces()),
@@ -1031,7 +1029,7 @@ Foam::label Foam::FaceCellWave<Type>::cellToFace()
         // Transfer changed faces across cyclic halves
         handleCyclicPatches();
     }
-    if (hasProcPatches_)
+    if (Pstream::parRun())
     {
         // Transfer changed faces from neighbouring processors.
         handleProcPatches();
@@ -1060,7 +1058,7 @@ Foam::label Foam::FaceCellWave<Type>::iterate(const label maxIter)
         // Transfer changed faces across cyclic halves
         handleCyclicPatches();
     }
-    if (hasProcPatches_)
+    if (Pstream::parRun())
     {
         // Transfer changed faces from neighbouring processors.
         handleProcPatches();

src/OpenFOAM/algorithms/MeshWave/FaceCellWave.H

@@ -32,8 +32,7 @@ Description
 
     Handles parallel and cyclics and non-parallel cyclics.
 
-Note
-    Whether to propagate depends on the return value of Type::update
+    Note: whether to propagate depends on the return value of Type::update
     which returns true (i.e. propagate) if the value changes by more than a
     certain tolerance.
     This tolerance can be very strict for normal face-cell and parallel
@@ -106,9 +105,6 @@ class FaceCellWave
         //- Contains cyclics
         bool hasCyclicPatches_;
 
-        //- Contains processor patches
-        bool hasProcPatches_;
-
         //- Number of evaluations
         label nEvals_;
 

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

@@ -164,9 +164,8 @@ public:
             //- Get value at index I
             inline unsigned int get(const label i) const;
 
-
-            //- Set value at index I
-            inline void set(const label i, const unsigned int val);
+            //- Set value at index I. Return true if value changed.
+            inline bool set(const label i, const unsigned int val);
 
 
     // Member operators

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

@@ -151,7 +151,7 @@ inline unsigned int PackedList<nBits>::operator[](const label i) const
 
 // Set value at i
 template<int nBits>
-inline void PackedList<nBits>::set(const label i, const unsigned int val)
+inline bool PackedList<nBits>::set(const label i, const unsigned int val)
 {
 #   ifdef DEBUGList
     checkIndex(i);
@@ -175,7 +175,11 @@ inline void PackedList<nBits>::set(const label i, const unsigned int val)
 
     unsigned int& elem = List<unsigned int>::operator[](intIndex(i));
 
+    unsigned int oldElem = elem;
+
     elem = (elem & ~shiftedMask) | shiftedVal;
+
+    return elem != oldElem;
 }
 
 

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

@@ -438,6 +438,38 @@ public:
             template <class T>
             static void listCombineScatter(List<T>& Value);
 
+        // Combine variants working on whole map at a time. Container needs to
+        // have iterators and find() defined.
+
+            template <class Container, class CombineOp>
+            static void mapCombineGather
+            (
+                const List<commsStruct>& comms,
+                Container& Values,
+                const CombineOp& cop
+            );
+
+            //- Like above but switches between linear/tree communication
+            template <class Container, class CombineOp>
+            static void mapCombineGather
+            (
+                Container& Values,
+                const CombineOp& cop
+            );
+
+            //- Scatter data. Reverse of combineGather
+            template <class Container>
+            static void mapCombineScatter
+            (
+                const List<commsStruct>& comms,
+                Container& Values
+            );
+
+            //- Like above but switches between linear/tree communication
+            template <class Container>
+            static void mapCombineScatter(Container& Values);
+            
+
 
         // Gather/scatter keeping the individual processor data separate.
         // Values is a List of size Pstream::nProcs() where

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

@@ -409,6 +409,147 @@ void Pstream::listCombineScatter(List<T>& Values)
 
 
 
+// Same thing but for sparse list (map)
+// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+
+template <class Container, class CombineOp>
+void Pstream::mapCombineGather
+(
+    const List<Pstream::commsStruct>& comms,
+    Container& Values,
+    const CombineOp& cop
+)
+{
+    if (Pstream::parRun())
+    {
+        // Get my communication order
+        const commsStruct& myComm = comms[Pstream::myProcNo()];
+
+        // Receive from my downstairs neighbours
+        forAll(myComm.below(), belowI)
+        {
+            label belowID = myComm.below()[belowI];
+
+            IPstream fromBelow(Pstream::scheduled, belowID);
+            Container receivedValues(fromBelow);
+
+            if (debug & 2)
+            {
+                Pout<< " received from "
+                    << belowID << " data:" << receivedValues << endl;
+            }
+
+            for
+            (
+                typename Container::const_iterator slaveIter =
+                    receivedValues.begin();
+                slaveIter != receivedValues.end();
+                ++slaveIter
+            )
+            {
+                typename Container::iterator
+                    masterIter = Values.find(slaveIter.key());
+
+                if (masterIter != Values.end())
+                {
+                    cop(masterIter(), slaveIter());
+                }
+                else
+                {
+                    Values.insert(slaveIter.key(), slaveIter());
+                }
+            }
+        }
+
+        // Send up Value
+        if (myComm.above() != -1)
+        {
+            if (debug & 2)
+            {
+                Pout<< " sending to " << myComm.above()
+                    << " data:" << Values << endl;
+            }
+
+            OPstream toAbove(Pstream::scheduled, myComm.above());
+            toAbove << Values;
+        }
+    }
+}
+
+
+template <class Container, class CombineOp>
+void Pstream::mapCombineGather(Container& Values, const CombineOp& cop)
+{
+    if (Pstream::nProcs() < Pstream::nProcsSimpleSum)
+    {
+        mapCombineGather(Pstream::linearCommunication(), Values, cop);
+    }
+    else
+    {
+        mapCombineGather(Pstream::treeCommunication(), Values, cop);
+    }
+}
+
+
+template <class Container>
+void Pstream::mapCombineScatter
+(
+    const List<Pstream::commsStruct>& comms,
+    Container& Values
+)
+{
+    if (Pstream::parRun())
+    {
+        // Get my communication order
+        const Pstream::commsStruct& myComm = comms[Pstream::myProcNo()];
+
+        // Reveive from up
+        if (myComm.above() != -1)
+        {
+            IPstream fromAbove(Pstream::scheduled, myComm.above());
+            fromAbove >> Values;
+
+            if (debug & 2)
+            {
+                Pout<< " received from "
+                    << myComm.above() << " data:" << Values << endl;
+            }
+        }
+
+        // Send to my downstairs neighbours
+        forAll(myComm.below(), belowI)
+        {
+            label belowID = myComm.below()[belowI];
+
+            if (debug & 2)
+            {
+                Pout<< " sending to " << belowID << " data:" << Values << endl;
+            }
+
+            OPstream toBelow(Pstream::scheduled, belowID);
+            toBelow << Values;
+        }
+    }
+}
+
+
+template <class Container>
+void Pstream::mapCombineScatter(Container& Values)
+{
+    if (Pstream::nProcs() < Pstream::nProcsSimpleSum)
+    {
+        mapCombineScatter(Pstream::linearCommunication(), Values);
+    }
+    else
+    {
+        mapCombineScatter(Pstream::treeCommunication(), Values);
+    }
+}
+
+
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/meshes/meshShapes/face/face.H

@@ -188,7 +188,7 @@ public:
         //- Return potential intersection with face with a ray starting
         //  at p, direction n (does not need to be normalized)
         //  Does face-center decomposition and returns triangle intersection
-        //  point closest to p.
+        //  point closest to p. Face-center is calculated from point average.
         //  For a hit, the distance is signed.  Positive number
         //  represents the point in front of triangle
         //  In case of miss the point is the nearest point on the face
@@ -206,6 +206,19 @@ public:
             const intersection::direction dir = intersection::VECTOR
         ) const;
 
+        //- Fast intersection with a ray.
+        //  For a hit, the pointHit.distance() is the line parameter t :
+        //  intersection=p+t*q. Only defined for FULL_RAY or
+        //  HALF_RAY.
+        pointHit intersection
+        (
+            const point& p,
+            const vector& q,
+            const point& ctr,
+            const pointField& meshPoints,
+            const intersection::algorithm alg
+        ) const;
+
         //- Return nearest point to face
         pointHit nearestPoint
         (

src/OpenFOAM/meshes/meshShapes/face/faceIntersection.C

@@ -133,6 +133,52 @@ pointHit face::ray
 }
 
 
+pointHit face::intersection
+(
+    const point& p,
+    const vector& q,
+    const point& ctr,
+    const pointField& meshPoints,
+    const intersection::algorithm alg
+) const
+{
+    scalar nearestHitDist = VGREAT;
+
+    // Initialize to miss, distance = GREAT
+    pointHit nearest(p);
+
+    const labelList& f = *this;
+
+    forAll(f, pI)
+    {
+        // Note: for best accuracy, centre point always comes last
+        pointHit curHit = triPointRef
+        (
+            meshPoints[f[pI]],
+            meshPoints[f[fcIndex(pI)]],
+            ctr
+        ).intersection(p, q, alg);
+
+        if (curHit.hit())
+        {
+            if (Foam::mag(curHit.distance()) < nearestHitDist)
+            {
+                nearestHitDist = Foam::mag(curHit.distance());
+                nearest.setHit();
+                nearest.setPoint(curHit.hitPoint());
+            }
+        }
+    }
+
+    if (nearest.hit())
+    {
+        nearest.setDistance(nearestHitDist);
+    }
+
+    return nearest;
+}
+
+
 pointHit face::nearestPoint
 (
     const point& p,

src/OpenFOAM/meshes/meshTools/matchPoints.C

@@ -64,9 +64,11 @@ bool Foam::matchPoints
             startI = 0;
         }
 
-        label face1I = -1;
 
-        // Go through range of equal mag and find equal vector.
+        // Go through range of equal mag and find nearest vector.
+        scalar minDistSqr = VGREAT;
+        label minFaceI = -1;
+    
         for
         (
             label j = startI;
@@ -78,17 +80,17 @@ bool Foam::matchPoints
         )
         {
             label faceI = pts1MagSqr.indices()[j];
-
             // Compare actual vectors
-            if (magSqr(pts0[face0I] - pts1[faceI]) <= sqr(matchDist))
-            {
-                face1I = faceI;
+            scalar distSqr = magSqr(pts0[face0I] - pts1[faceI]);
 
-                break;
+            if (distSqr <= sqr(matchDist) && distSqr < minDistSqr)
+            {
+                minDistSqr = distSqr;
+                minFaceI = faceI;
             }
         }
 
-        if (face1I == -1)
+        if (minFaceI == -1)
         {
             fullMatch = false;
 
@@ -120,7 +122,7 @@ bool Foam::matchPoints
             }
         }
 
-        from0To1[face0I] = face1I;
+        from0To1[face0I] = minFaceI;
     }
 
     return fullMatch;

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

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "mapDistribute.H"
+#include "commSchedule.H"
+#include "HashSet.H"
+#include "ListOps.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::mapDistribute::calcSchedule() const
+{
+    // Communications: send and receive processor
+    List<labelPair> allComms;
+
+    {
+        HashSet<labelPair, labelPair::Hash<> > commsSet(Pstream::nProcs());
+
+        // Find what communication is required
+        forAll(subMap_, procI)
+        {
+            if (procI != Pstream::myProcNo())
+            {
+                if (subMap_[procI].size() > 0)
+                {
+                    // I need to send to procI
+                    commsSet.insert(labelPair(Pstream::myProcNo(), procI));
+                }
+                if (constructMap_[procI].size() > 0)
+                {
+                    // I need to receive from procI
+                    commsSet.insert(labelPair(procI, Pstream::myProcNo()));
+                }
+            }
+        }
+        allComms = commsSet.toc();
+    }
+
+
+    // Reduce
+    if (Pstream::master())
+    {
+        // Receive and merge
+        for
+        (
+            int slave=Pstream::firstSlave();
+            slave<=Pstream::lastSlave();
+            slave++
+        )
+        {
+            IPstream fromSlave(Pstream::blocking, slave);
+            List<labelPair> nbrData(fromSlave);
+
+            forAll(nbrData, i)
+            {
+                if (findIndex(allComms, nbrData[i]) == -1)
+                {
+                    label sz = allComms.size();
+                    allComms.setSize(sz+1);
+                    allComms[sz] = nbrData[i];
+                }
+            }
+        }
+        // Send back
+        for
+        (
+            int slave=Pstream::firstSlave();
+            slave<=Pstream::lastSlave();
+            slave++
+        )
+        {
+            OPstream toSlave(Pstream::blocking, slave);
+            toSlave << allComms;
+        }
+    }
+    else
+    {
+        {
+            OPstream toMaster(Pstream::blocking, Pstream::masterNo());
+            toMaster << allComms;
+        }
+        {
+            IPstream fromMaster(Pstream::blocking, Pstream::masterNo());
+            fromMaster >> allComms;
+        }
+    }
+
+
+    // Determine my schedule.
+    labelList mySchedule
+    (
+        commSchedule
+        (
+            Pstream::nProcs(),
+            allComms
+        ).procSchedule()[Pstream::myProcNo()]
+    );
+
+    // Processors involved in my schedule
+    schedulePtr_.reset
+    (
+        new List<labelPair>
+        (
+            IndirectList<labelPair>(allComms, mySchedule)
+        )
+    );
+
+
+    //if (debug)
+    //{
+    //    Pout<< "I need to:" << endl;
+    //    const List<labelPair>& comms = schedule();
+    //    forAll(comms, i)
+    //    {
+    //        const labelPair& twoProcs = comms[i];
+    //        label sendProc = twoProcs[0];
+    //        label recvProc = twoProcs[1];
+    //
+    //        if (recvProc == Pstream::myProcNo())
+    //        {
+    //            Pout<< "    receive from " << sendProc << endl;
+    //        }
+    //        else
+    //        {
+    //            Pout<< "    send to " << recvProc << endl;
+    //        }
+    //    }
+    //}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+//- Construct from components
+Foam::mapDistribute::mapDistribute
+(
+    const label constructSize,