Merge branch 'master' of /home/dm4/OpenFOAM/OpenFOAM-dev

applications/solvers/incompressible/simpleFoam/porousSimpleFoam/Make/options

@@ -5,11 +5,16 @@ EXE_INC = \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/fieldSources/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
+
 
 EXE_LIBS = \
     -lincompressibleTurbulenceModel \
     -lincompressibleRASModels \
     -lincompressibleTransportModels \
     -lfiniteVolume \
-    -lmeshTools
+    -lmeshTools \
+    -lfieldSources \
+    -lsampling

applications/test/findSphereFeatureEdges-octree/Make/files

@@ -0,0 +1,3 @@
+Test-findSphereFeatureEdges-octree.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-findSphereFeatureEdges-octree

applications/test/findSphereFeatureEdges-octree/Make/options

@@ -0,0 +1,9 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/edgeMesh/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools \
+    -ledgeMesh

applications/test/findSphereFeatureEdges-octree/Test-findSphereFeatureEdges-octree.C

@@ -0,0 +1,134 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "Time.H"
+#include "fvMesh.H"
+#include "IStringStream.H"
+#include "indexedOctree.H"
+#include "treeDataEdge.H"
+#include "OFstream.H"
+#include "extendedFeatureEdgeMesh.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// Main program:
+
+int main(int argc, char *argv[])
+{
+
+#   include "setRootCase.H"
+#   include "createTime.H"
+#   include "createMesh.H"
+
+    fileName sFeatFileName("unit_cube_rotated.extendedFeatureEdgeMesh");
+
+    extendedFeatureEdgeMesh efem
+    (
+        IOobject
+        (
+            sFeatFileName,
+            runTime.time().constant(),
+            "extendedFeatureEdgeMesh",
+            runTime.time(),
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    );
+
+    // Slightly extended bb. Slightly off-centred just so on symmetric
+    // geometry there are less face/edge aligned items.
+    treeBoundBox bb
+    (
+        efem.points()
+    );
+
+    bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+    bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+
+    labelList allEdges(identity(efem.edges().size()));
+
+    indexedOctree<treeDataEdge> edgeTree
+    (
+        treeDataEdge
+        (
+            false,          // cachebb
+            efem.edges(),        // edges
+            efem.points(),       // points
+            allEdges        // selected edges
+        ),
+        bb,     // bb
+        8,      // maxLevel
+        10,     // leafsize
+        3.0     // duplicity
+    );
+
+    Info<< "Points: " << efem.points() << nl << endl;
+    Info<< "Edges: " << efem.edges() << nl << endl;
+
+    Info<< "Find edge labels within sphere from point (0, 0, 0):" << endl;
+
+    Info<< "    Radius = 0            : "
+        << edgeTree.findSphere(point(0, 0, 0), 0) << endl;
+
+    Info<< "    Radius = 1            : "
+        << edgeTree.findSphere(point(0, 0, 0), 1) << endl;
+
+    Info<< "    Radius = root(1.5)    : "
+        << edgeTree.findSphere(point(0, 0, 0), 1.5) << endl;
+
+    Info<< "    Radius = root(2)      : "
+        << edgeTree.findSphere(point(0, 0, 0), 2) << endl;
+
+    Info<< "    Radius = root(0.5)    : "
+        << edgeTree.findSphere(point(1, 0, 0), 0.5) << endl;
+
+    Info<< "    Radius = root(0.25)   : "
+        << edgeTree.findSphere(point(0, 0, 0.5), 0.25) << endl;
+
+    treeBoundBox tbb(point(0,0,0), point(0.1,0.1,0.1));
+    Info<< "    Box = " << tbb << "    : "
+        << edgeTree.findBox(tbb) << endl;
+
+    treeBoundBox tbb1(point(0,0,0), point(1,1,0.1));
+    Info<< "    Box = " << tbb1 << "        : "
+        << edgeTree.findBox(tbb1) << endl;
+
+    treeBoundBox tbb2(point(0.3,0,0), point(1,0.3,1));
+    Info<< "    Box = " << tbb2 << "      : "
+        << edgeTree.findBox(tbb2) << endl;
+
+    treeBoundBox tbb3(point(-0.2,0.5,0), point(0.3,0.9,1));
+    Info<< "    Box = " << tbb3 << " : "
+        << edgeTree.findBox(tbb3) << endl;
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/test/momentOfInertia/Test-momentOfInertia.C

@@ -2,7 +2,7 @@
  =========                   |
  \\      /   F ield          | OpenFOAM: The Open Source CFD Toolbox
   \\    /    O peration      |
-   \\  /     A nd            | Copyright (C) 2011 OpenFOAM Foundation
+   \\  /     A nd            | Copyright (C) 2012 OpenFOAM Foundation
     \\/      M anipulation   |
 -------------------------------------------------------------------------------
 License
@@ -35,7 +35,7 @@ Description
 #include "polyMesh.H"
 #include "ListOps.H"
 #include "face.H"
-#include "tetPointRef.H"
+#include "tetrahedron.H"
 #include "triFaceList.H"
 #include "OFstream.H"
 #include "meshTools.H"

applications/utilities/mesh/manipulation/createBaffles/createBaffles.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -290,18 +290,6 @@ int main(int argc, char *argv[])
 
     // Whether first use of face (modify) or consecutive (add)
     PackedBoolList modifiedFace(mesh.nFaces());
-    // Never modify coupled faces
-    forAll(patches, patchI)
-    {
-        const polyPatch& pp = patches[patchI];
-        if (pp.coupled())
-        {
-            forAll(pp, i)
-            {
-                modifiedFace[pp.start()+i] = 1;
-            }
-        }
-    }
     label nModified = 0;
 
     forAll(newMasterPatches, i)

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -254,5 +254,11 @@ Foam::UPstream::commsTypes Foam::UPstream::defaultCommsType
     commsTypeNames.read(debug::optimisationSwitches().lookup("commsType"))
 );
 
+// Number of polling cycles in processor updates
+int Foam::UPstream::nPollProcInterfaces
+(
+    debug::optimisationSwitch("nPollProcInterfaces", 0)
+);
+
 
 // ************************************************************************* //

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -242,6 +242,8 @@ public:
         //- Default commsType
         static commsTypes defaultCommsType;
 
+        //- Number of polling cycles in processor updates
+        static int nPollProcInterfaces;
 
     // Constructors
 
@@ -273,6 +275,9 @@ public:
             //- Wait until all requests (from start onwards) have finished.
             static void waitRequests(const label start = 0);
 
+            //- Wait until request i has finished.
+            static void waitRequest(const label i);
+
             //- Non-blocking comms: has request i finished?
             static bool finishedRequest(const label i);
 

src/OpenFOAM/global/argList/argList.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -772,10 +772,11 @@ Foam::argList::argList
                 Info<< "Roots  : " << roots << nl;
             }
             Info<< "Pstream initialized with:" << nl
-                << "    floatTransfer     : " << Pstream::floatTransfer << nl
+                << "    floatTransfer      : " << Pstream::floatTransfer << nl
-                << "    nProcsSimpleSum   : " << Pstream::nProcsSimpleSum << nl
+                << "    nProcsSimpleSum    : " << Pstream::nProcsSimpleSum << nl
-                << "    commsType         : "
+                << "    commsType          : "
-                << Pstream::commsTypeNames[Pstream::defaultCommsType]
+                << Pstream::commsTypeNames[Pstream::defaultCommsType] << nl
+                << "    polling iterations : " << Pstream::nPollProcInterfaces
                 << endl;
         }
     }

src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterfaceFields/lduInterfaceField/lduInterfaceField.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -60,6 +60,10 @@ class lduInterfaceField
         //- Reference to the coupled patch this field is defined for
         const lduInterface& interface_;
 
+        //- Update index used so that updateInterfaceMatrix is called only once
+        //  during the construction of the matrix
+        bool updatedMatrix_;
+
 
     // Private Member Functions
 
@@ -69,7 +73,6 @@ class lduInterfaceField
         //- Disallow default bitwise assignment
         void operator=(const lduInterfaceField&);
 
-
 public:
 
     //- Runtime type information
@@ -81,7 +84,8 @@ public:
         //- Construct given coupled patch
         lduInterfaceField(const lduInterface& patch)
         :
-            interface_(patch)
+            interface_(patch),
+            updatedMatrix_(false)
         {}
 
 
@@ -120,6 +124,24 @@ public:
             ) const
             {}
 
+            //- Whether matrix has been updated
+            bool updatedMatrix() const
+            {
+                return updatedMatrix_;
+            }
+
+            //- Whether matrix has been updated
+            bool& updatedMatrix()
+            {
+                return updatedMatrix_;
+            }
+
+            //- Is all data available
+            virtual bool ready() const
+            {
+                return true;
+            }
+
             //- Update result field based on interface functionality
             virtual void updateInterfaceMatrix
             (

src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixOperations.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -205,7 +205,15 @@ void Foam::lduMatrix::operator+=(const lduMatrix& A)
         if (debug > 1)
         {
             WarningIn("lduMatrix::operator+=(const lduMatrix& A)")
-                << "Unknown matrix type combination"
+                << "Unknown matrix type combination" << nl
+                << "    this :"
+                << " diagonal:" << diagonal()
+                << " symmetric:" << symmetric()
+                << " asymmetric:" << asymmetric() << nl
+                << "    A    :"
+                << " diagonal:" << A.diagonal()
+                << " symmetric:" << A.symmetric()
+                << " asymmetric:" << A.asymmetric()
                 << endl;
         }
     }
@@ -276,7 +284,15 @@ void Foam::lduMatrix::operator-=(const lduMatrix& A)
         if (debug > 1)
         {
             WarningIn("lduMatrix::operator-=(const lduMatrix& A)")
-                << "Unknown matrix type combination"
+                << "Unknown matrix type combination" << nl
+                << "    this :"
+                << " diagonal:" << diagonal()
+                << " symmetric:" << symmetric()
+                << " asymmetric:" << asymmetric() << nl
+                << "    A    :"
+                << " diagonal:" << A.diagonal()
+                << " symmetric:" << A.symmetric()
+                << " asymmetric:" << A.asymmetric()
                 << endl;
         }
     }

src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixUpdateMatrixInterfaces.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -87,7 +87,7 @@ void Foam::lduMatrix::initMatrixInterfaces
     }
     else
     {
-        FatalErrorIn("lduMatrix::initMatrixInterfaces")
+        FatalErrorIn("lduMatrix::initMatrixInterfaces(..)")
             << "Unsuported communications type "
             << Pstream::commsTypeNames[Pstream::defaultCommsType]
             << exit(FatalError);
@@ -104,22 +104,8 @@ void Foam::lduMatrix::updateMatrixInterfaces
     const direction cmpt
 ) const
 {
-    if
+    if (Pstream::defaultCommsType == Pstream::blocking)
-    (
-        Pstream::defaultCommsType == Pstream::blocking
-     || Pstream::defaultCommsType == Pstream::nonBlocking
-    )
     {
-        // Block until all sends/receives have been finished
-        if
-        (
-            Pstream::parRun()
-         && Pstream::defaultCommsType == Pstream::nonBlocking
-        )
-        {
-            UPstream::waitRequests();
-        }
-
         forAll(interfaces, interfaceI)
         {
             if (interfaces.set(interfaceI))
@@ -136,6 +122,87 @@ void Foam::lduMatrix::updateMatrixInterfaces
             }
         }
     }
+    else if (Pstream::defaultCommsType == Pstream::nonBlocking)
+    {
+        // Try and consume interfaces as they become available
+        bool allUpdated = false;
+
+        for (label i = 0; i < UPstream::nPollProcInterfaces; i++)
+        {
+            allUpdated = true;
+
+            forAll(interfaces, interfaceI)
+            {
+                if (interfaces.set(interfaceI))
+                {
+                    if (!interfaces[interfaceI].updatedMatrix())
+                    {
+                        if (interfaces[interfaceI].ready())
+                        {
+                            interfaces[interfaceI].updateInterfaceMatrix
+                            (
+                                psiif,
+                                result,
+                                *this,
+                                coupleCoeffs[interfaceI],
+                                cmpt,
+                                Pstream::defaultCommsType
+                            );
+                        }
+                        else
+                        {
+                            allUpdated = false;
+                        }
+                    }
+                }
+            }
+
+            if (allUpdated)
+            {
+                break;
+            }
+        }
+
+        // Block for everything
+        if (Pstream::parRun())
+        {
+            if (allUpdated)
+            {
+                // All received. Just remove all storage of requests
+                // Note that we don't know what starting number of requests
+                // was before start of sends and receives (since set from
+                // initMatrixInterfaces) so set to 0 and loose any in-flight
+                // requests.
+                UPstream::resetRequests(0);
+            }
+            else
+            {
+                // Block for all requests and remove storage
+                UPstream::waitRequests();
+            }
+        }
+
+        // Consume
+        forAll(interfaces, interfaceI)
+        {
+            if
+            (
+                interfaces.set(interfaceI)
+            && !interfaces[interfaceI].updatedMatrix()
+            )
+            {
+                interfaces[interfaceI].updateInterfaceMatrix
+                (
+                    psiif,
+                    result,
+                    *this,
+                    coupleCoeffs[interfaceI],
+                    cmpt,
+                    Pstream::defaultCommsType
+                );
+            }
+        }
+    }
     else if (Pstream::defaultCommsType == Pstream::scheduled)
     {
         const lduSchedule& patchSchedule = this->patchSchedule();
@@ -199,7 +266,7 @@ void Foam::lduMatrix::updateMatrixInterfaces
     }
     else
     {
-        FatalErrorIn("lduMatrix::updateMatrixInterfaces")
+        FatalErrorIn("lduMatrix::updateMatrixInterfaces(..)")
             << "Unsuported communications type "
             << Pstream::commsTypeNames[Pstream::defaultCommsType]
             << exit(FatalError);

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -80,11 +80,38 @@ void Foam::processorGAMGInterfaceField::initInterfaceMatrixUpdate
     const Pstream::commsTypes commsType
 ) const
 {
-    procInterface_.compressedSend
+    procInterface_.interfaceInternalField(psiInternal, scalarSendBuf_);
-    (
+
-        commsType,
+    if (commsType == Pstream::nonBlocking && !Pstream::floatTransfer)
-        procInterface_.interfaceInternalField(psiInternal)()
+    {
-    );
+        // Fast path.
+        scalarReceiveBuf_.setSize(scalarSendBuf_.size());
+        outstandingRecvRequest_ = UPstream::nRequests();
+        IPstream::read
+        (
+            Pstream::nonBlocking,
+            procInterface_.neighbProcNo(),
+            reinterpret_cast<char*>(scalarReceiveBuf_.begin()),
+            scalarReceiveBuf_.byteSize(),
+            procInterface_.tag()
+        );
+
+        outstandingSendRequest_ = UPstream::nRequests();
+        OPstream::write
+        (
+            Pstream::nonBlocking,
+            procInterface_.neighbProcNo(),
+            reinterpret_cast<const char*>(scalarSendBuf_.begin()),
+            scalarSendBuf_.byteSize(),
+            procInterface_.tag()
+        );
+    }
+    else
+    {
+        procInterface_.compressedSend(commsType, scalarSendBuf_);
+    }
+
+    const_cast<processorGAMGInterfaceField&>(*this).updatedMatrix() = false;
 }
 
 
@@ -98,18 +125,54 @@ void Foam::processorGAMGInterfaceField::updateInterfaceMatrix
     const Pstream::commsTypes commsType
 ) const
 {
-    scalarField pnf
+    if (updatedMatrix())
-    (
+    {
-        procInterface_.compressedReceive<scalar>(commsType, coeffs.size())
+        return;
-    );
+    }
-    transformCoupleField(pnf, cmpt);
 
     const labelUList& faceCells = procInterface_.faceCells();
 
-    forAll(faceCells, elemI)
+    if (commsType == Pstream::nonBlocking && !Pstream::floatTransfer)
     {
-        result[faceCells[elemI]] -= coeffs[elemI]*pnf[elemI];
+        // Fast path.
+        if
+        (
+            outstandingRecvRequest_ >= 0
+         && outstandingRecvRequest_ < Pstream::nRequests()
+        )
+        {
+            UPstream::waitRequest(outstandingRecvRequest_);
+        }
+        // Recv finished so assume sending finished as well.
+        outstandingSendRequest_ = -1;
+        outstandingRecvRequest_ = -1;
+
+        // Consume straight from scalarReceiveBuf_
+
+        // Transform according to the transformation tensor
+        transformCoupleField(scalarReceiveBuf_, cmpt);
+
+        // Multiply the field by coefficients and add into the result
+        forAll(faceCells, elemI)
+        {
+            result[faceCells[elemI]] -= coeffs[elemI]*scalarReceiveBuf_[elemI];
+        }
     }
+    else
+    {
+        scalarField pnf
+        (
+            procInterface_.compressedReceive<scalar>(commsType, coeffs.size())
+        );
+        transformCoupleField(pnf, cmpt);
+
+        forAll(faceCells, elemI)
+        {
+            result[faceCells[elemI]] -= coeffs[elemI]*pnf[elemI];
+        }
+    }
+
+    const_cast<processorGAMGInterfaceField&>(*this).updatedMatrix() = true;
 }
 
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/processorGAMGInterfaceField/processorGAMGInterfaceField.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -65,6 +65,22 @@ class processorGAMGInterfaceField
         int rank_;
 
 
+        // Sending and receiving
+
+            //- Outstanding request
+            mutable label outstandingSendRequest_;
+
+            //- Outstanding request
+            mutable label outstandingRecvRequest_;
+
+            //- Scalar send buffer
+            mutable Field<scalar> scalarSendBuf_;
+
+            //- Scalar receive buffer
+            mutable Field<scalar> scalarReceiveBuf_;
+
+
+
     // Private Member Functions
 
         //- Disallow default bitwise copy construct

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -188,6 +188,14 @@ public:
                 const UList<Type>& internalData
             ) const;
 
+            //- Get the interface internal field of the given field
+            template<class Type>
+            void interfaceInternalField
+            (
+                const UList<Type>& internalData,
+                List<Type>&
+            ) const;
+
             //- Return the values of the given internal data adjacent to
             //  the interface as a field
             virtual tmp<labelField> interfaceInternalField

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaces/GAMGInterface/GAMGInterfaceTemplates.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -34,14 +34,24 @@ Foam::tmp<Foam::Field<Type> > Foam::GAMGInterface::interfaceInternalField
 ) const
 {
     tmp<Field<Type> > tresult(new Field<Type>(size()));
-    Field<Type>& result = tresult();
+    interfaceInternalField(iF, tresult());
+    return tresult;
+}
+
+
+template<class Type>
+void Foam::GAMGInterface::interfaceInternalField
+(
+    const UList<Type>& iF,
+    List<Type>& result
+) const
+{
+    result.setSize(size());
 
     forAll(result, elemI)
     {
         result[elemI] = iF[faceCells_[elemI]];
     }
-
-    return tresult;
 }
 
 

src/OpenFOAM/meshes/meshShapes/tetCell/tetCell.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -43,7 +43,7 @@ SourceFiles
 #include "triFace.H"
 #include "edge.H"
 #include "pointField.H"
-#include "tetPointRef.H"
+#include "tetrahedron.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -42,7 +42,7 @@ SourceFiles
 #include "polyMesh.H"
 #include "coupledPolyPatch.H"
 #include "syncTools.H"
-#include "tetPointRef.H"
+#include "tetrahedron.H"
 #include "tetIndices.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -38,7 +38,7 @@ SourceFiles
 #define tetIndices_H
 
 #include "label.H"
-#include "tetPointRef.H"
+#include "tetrahedron.H"
 #include "triPointRef.H"
 #include "polyMesh.H"
 #include "triFace.H"

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -419,146 +419,148 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors
             << gAverage(half1Ctrs) << endl;
     }
 
-    switch (transform_)
+    if (half0Ctrs.size())
     {
-        case ROTATIONAL:
+        switch (transform_)
         {
-            vector n0 = findFaceMaxRadius(half0Ctrs);
+            case ROTATIONAL:
-            vector n1 = -findFaceMaxRadius(half1Ctrs);
-            n0 /= mag(n0) + VSMALL;
-            n1 /= mag(n1) + VSMALL;
-
-            if (debug)
             {
-                Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
+                vector n0 = findFaceMaxRadius(half0Ctrs);
-                    << " patch:" << name()
+                vector n1 = -findFaceMaxRadius(half1Ctrs);
-                    << " Specified rotation :"
+                n0 /= mag(n0) + VSMALL;
-                    << " n0:" << n0 << " n1:" << n1 << endl;
+                n1 /= mag(n1) + VSMALL;
-            }
 
-            // Extended tensor from two local coordinate systems calculated
-            // using normal and rotation axis
-            const tensor E0
-            (
-                rotationAxis_,
-                (n0 ^ rotationAxis_),
-                n0
-            );
-            const tensor E1
-            (
-                rotationAxis_,
-                (-n1 ^ rotationAxis_),
-                -n1
-            );
-            const tensor revT(E1.T() & E0);
-
-            // Rotation
-            forAll(half0Ctrs, faceI)
-            {
-                half0Ctrs[faceI] =
-                    Foam::transform
-                    (
-                        revT,
-                        half0Ctrs[faceI] - rotationCentre_
-                    )
-                  + rotationCentre_;
-                anchors0[faceI] =
-                    Foam::transform
-                    (
-                        revT,
-                        anchors0[faceI] - rotationCentre_
-                    )
-                  + rotationCentre_;
-            }
-
-            break;
-        }
-        case TRANSLATIONAL:
-        {
-            // Transform 0 points.
-
-            if (debug)
-            {
-                Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
-                    << " patch:" << name()
-                    << "Specified translation : " << separationVector_
-                    << endl;
-            }
-
-            // Note: getCentresAndAnchors gets called on the slave side
-            // so separationVector is owner-slave points.
-
-            half0Ctrs -= separationVector_;
-            anchors0 -= separationVector_;
-            break;
-        }
-        default:
-        {
-            // Assumes that cyclic is planar. This is also the initial
-            // condition for patches without faces.
-
-            // Determine the face with max area on both halves. These
-            // two faces are used to determine the transformation tensors
-            label max0I = findMaxArea(pp0.points(), pp0);
-            vector n0 = pp0[max0I].normal(pp0.points());
-            n0 /= mag(n0) + VSMALL;
-
-            label max1I = findMaxArea(pp1.points(), pp1);
-            vector n1 = pp1[max1I].normal(pp1.points());
-            n1 /= mag(n1) + VSMALL;
-
-            if (mag(n0 & n1) < 1-matchTolerance())
-            {
                 if (debug)
                 {
                     Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
                         << " patch:" << name()
-                        << " Detected rotation :"
+                        << " Specified rotation :"
                         << " n0:" << n0 << " n1:" << n1 << endl;
                 }
 
-                // Rotation (around origin)
+                // Extended tensor from two local coordinate systems calculated
-                const tensor revT(rotationTensor(n0, -n1));
+                // using normal and rotation axis
+                const tensor E0
+                (
+                    rotationAxis_,
+                    (n0 ^ rotationAxis_),
+                    n0
+                );
+                const tensor E1
+                (
+                    rotationAxis_,
+                    (-n1 ^ rotationAxis_),
+                    -n1
+                );
+                const tensor revT(E1.T() & E0);
 
                 // Rotation
                 forAll(half0Ctrs, faceI)
                 {
-                    half0Ctrs[faceI] = Foam::transform
+                    half0Ctrs[faceI] =
-                    (
+                        Foam::transform
-                        revT,
+                        (
-                        half0Ctrs[faceI]
+                            revT,
-                    );
+                            half0Ctrs[faceI] - rotationCentre_
-                    anchors0[faceI] = Foam::transform
+                        )
-                    (
+                      + rotationCentre_;
-                        revT,
+                    anchors0[faceI] =
-                        anchors0[faceI]
+                        Foam::transform
-                    );
+                        (
+                            revT,
+                            anchors0[faceI] - rotationCentre_
+                        )
+                      + rotationCentre_;
                 }
-            }
-            else
-            {
-                // Parallel translation. Get average of all used points.
 
-                const point ctr0(sum(pp0.localPoints())/pp0.nPoints());
+                break;
-                const point ctr1(sum(pp1.localPoints())/pp1.nPoints());
+            }
+            case TRANSLATIONAL:
+            {
+                // Transform 0 points.
 
                 if (debug)
                 {
                     Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
                         << " patch:" << name()
-                        << " Detected translation :"
+                        << "Specified translation : " << separationVector_
-                        << " n0:" << n0 << " n1:" << n1
+                        << endl;
-                        << " ctr0:" << ctr0 << " ctr1:" << ctr1 << endl;
                 }
 
-                half0Ctrs += ctr1 - ctr0;
+                // Note: getCentresAndAnchors gets called on the slave side
-                anchors0 += ctr1 - ctr0;
+                // so separationVector is owner-slave points.
+
+                half0Ctrs -= separationVector_;
+                anchors0 -= separationVector_;
+                break;
+            }
+            default:
+            {
+                // Assumes that cyclic is planar. This is also the initial
+                // condition for patches without faces.
+
+                // Determine the face with max area on both halves. These
+                // two faces are used to determine the transformation tensors
+                label max0I = findMaxArea(pp0.points(), pp0);
+                vector n0 = pp0[max0I].normal(pp0.points());
+                n0 /= mag(n0) + VSMALL;
+
+                label max1I = findMaxArea(pp1.points(), pp1);
+                vector n1 = pp1[max1I].normal(pp1.points());
+                n1 /= mag(n1) + VSMALL;
+
+                if (mag(n0 & n1) < 1-matchTolerance())
+                {
+                    if (debug)
+                    {
+                        Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
+                            << " patch:" << name()
+                            << " Detected rotation :"
+                            << " n0:" << n0 << " n1:" << n1 << endl;
+                    }
+
+                    // Rotation (around origin)
+                    const tensor revT(rotationTensor(n0, -n1));
+
+                    // Rotation
+                    forAll(half0Ctrs, faceI)
+                    {
+                        half0Ctrs[faceI] = Foam::transform
+                        (
+                            revT,
+                            half0Ctrs[faceI]
+                        );
+                        anchors0[faceI] = Foam::transform
+                        (
+                            revT,
+                            anchors0[faceI]
+                        );
+                    }
+                }
+                else
+                {
+                    // Parallel translation. Get average of all used points.
+
+                    const point ctr0(sum(pp0.localPoints())/pp0.nPoints());
+                    const point ctr1(sum(pp1.localPoints())/pp1.nPoints());
+
+                    if (debug)
+                    {
+                        Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
+                            << " patch:" << name()
+                            << " Detected translation :"
+                            << " n0:" << n0 << " n1:" << n1
+                            << " ctr0:" << ctr0 << " ctr1:" << ctr1 << endl;
+                    }
+
+                    half0Ctrs += ctr1 - ctr0;
+                    anchors0 += ctr1 - ctr0;
+                }
+                break;
             }
-            break;
         }
     }
 
-
     // Calculate typical distance per face
     tols = matchTolerance()*calcFaceTol(pp1, pp1.points(), half1Ctrs);
 }
@@ -1270,7 +1272,7 @@ bool Foam::cyclicPolyPatch::order
     rotation.setSize(pp.size());
     rotation = 0;
 
-    if (pp.empty() || transform_ == NOORDERING)
+    if (transform_ == NOORDERING)
     {
         // No faces, nothing to change.
         return false;

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,7 +25,7 @@ License
 
 #include "primitiveMesh.H"
 #include "pyramidPointFaceRef.H"
-#include "tetPointRef.H"
+#include "tetrahedron.H"
 #include "ListOps.H"
 #include "unitConversion.H"
 #include "SortableList.H"

src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetPointRef.H

@@ -1,52 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
-     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
-
-Typedef
-    Foam::tetPointRef
-
-Description
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef tetPointRef_H
-#define tetPointRef_H
-
-#include "point.H"
-#include "tetrahedron.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-typedef tetrahedron<point, const point&> tetPointRef;
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-#endif
-
-// ************************************************************************* //

src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetPoints.H

@@ -35,9 +35,9 @@ SourceFiles
 #ifndef tetPoints_H
 #define tetPoints_H
 
+#include "tetrahedron.H"
 #include "FixedList.H"
 #include "treeBoundBox.H"
-#include "tetPointRef.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -54,6 +54,8 @@ namespace Foam
 
 class Istream;
 class Ostream;
+class tetPoints;
+class plane;
 
 // Forward declaration of friend functions and operators
 
@@ -73,6 +75,7 @@ inline Ostream& operator<<
     const tetrahedron<Point, PointRef>&
 );
 
+typedef tetrahedron<point, const point&> tetPointRef;
 
 /*---------------------------------------------------------------------------*\
                            class tetrahedron Declaration
@@ -81,10 +84,71 @@ inline Ostream& operator<<
 template<class Point, class PointRef>
 class tetrahedron
 {
+public:
+
+        // Classes for use in sliceWithPlane. What to do with decomposition
+        // of tet.
+
+            //- Dummy
+            class dummyOp
+            {
+            public:
+                inline void operator()(const tetPoints&);
+            };
+
+            //- Sum resulting volumes
+            class sumVolOp
+            {
+            public:
+                scalar vol_;
+
+                inline sumVolOp();
+
+                inline void operator()(const tetPoints&);
+            };
+
+            //- Store resulting tets
+            class storeOp
+            {
+                FixedList<tetPoints, 200>& tets_;
+                label& nTets_;
+
+            public:
+                inline storeOp(FixedList<tetPoints, 200>&, label&);
+
+                inline void operator()(const tetPoints&);
+            };
+
+private:
+
     // Private data
 
         PointRef a_, b_, c_, d_;
 
+        inline static point planeIntersection
+        (
+            const FixedList<scalar, 4>&,
+            const tetPoints&,
+            const label,
+            const label
+        );
+
+        template<class TetOp>
+        inline static void decomposePrism
+        (
+            const FixedList<point, 6>& points,
+            TetOp& op
+        );
+
+        template<class AboveTetOp, class BelowTetOp>
+        inline static void tetSliceWithPlane
+        (
+            const plane& pl,
+            const tetPoints& tet,
+            AboveTetOp& aboveOp,
+            BelowTetOp& belowOp
+        );
+
 
 public:
 
@@ -184,6 +248,16 @@ public:
             //- Return true if point is inside tetrahedron
             inline bool inside(const point& pt) const;
 
+            //- Decompose tet into tets above and below plane
+            template<class AboveTetOp, class BelowTetOp>
+            inline void sliceWithPlane
+            (
+                const plane& pl,
+                AboveTetOp& aboveOp,
+                BelowTetOp& belowOp
+            ) const;
+
+
             //- Return (min)containment sphere, i.e. the smallest sphere with
             //  all points inside. Returns pointHit with:
             //  - hit         : if sphere is equal to circumsphere

src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,7 +25,8 @@ License
 
 #include "triangle.H"
 #include "IOstreams.H"
-#include "triPointRef.H"
+#include "tetPoints.H"
+#include "plane.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -492,6 +493,483 @@ bool Foam::tetrahedron<Point, PointRef>::inside(const point& pt) const
 }
 
 
+template<class Point, class PointRef>
+inline void Foam::tetrahedron<Point, PointRef>::dummyOp::operator()
+(
+    const tetPoints&
+)
+{}
+
+
+template<class Point, class PointRef>
+inline Foam::tetrahedron<Point, PointRef>::sumVolOp::sumVolOp()
+:
+    vol_(0.0)
+{}
+
+
+template<class Point, class PointRef>
+inline void Foam::tetrahedron<Point, PointRef>::sumVolOp::operator()
+(
+    const tetPoints& tet
+)
+{
+    vol_ += tet.tet().mag();
+}
+
+
+template<class Point, class PointRef>
+inline Foam::tetrahedron<Point, PointRef>::storeOp::storeOp
+(
+    FixedList<tetPoints, 200>& tets,
+    label& nTets
+)
+:
+    tets_(tets),
+    nTets_(nTets)
+{}
+
+
+template<class Point, class PointRef>
+inline void Foam::tetrahedron<Point, PointRef>::storeOp::operator()
+(
+    const tetPoints& tet
+)
+{
+    tets_[nTets_++] = tet;
+}
+
+
+template<class Point, class PointRef>
+inline Foam::point Foam::tetrahedron<Point, PointRef>::planeIntersection
+(
+    const FixedList<scalar, 4>& d,
+    const tetPoints& t,
+    const label negI,
+    const label posI
+)
+{
+    return
+        (d[posI]*t[negI] - d[negI]*t[posI])
+      / (-d[negI]+d[posI]);
+}
+
+
+template<class Point, class PointRef>
+template<class TetOp>
+inline void Foam::tetrahedron<Point, PointRef>::decomposePrism
+(
+    const FixedList<point, 6>& points,
+    TetOp& op
+)
+{
+    op(tetPoints(points[1], points[3], points[2], points[0]));
+    op(tetPoints(points[1], points[2], points[3], points[4]));
+    op(tetPoints(points[4], points[2], points[3], points[5]));
+}
+
+
+template<class Point, class PointRef>
+template<class AboveTetOp, class BelowTetOp>
+inline void Foam::tetrahedron<Point, PointRef>::
+tetSliceWithPlane
+(
+    const plane& pl,
+    const tetPoints& tet,
+    AboveTetOp& aboveOp,
+    BelowTetOp& belowOp
+)
+{
+    // Distance to plane
+    FixedList<scalar, 4> d;
+    label nPos = 0;
+    forAll(tet, i)
+    {
+        d[i] = ((tet[i]-pl.refPoint()) & pl.normal());
+        if (d[i] > 0)
+        {
+            nPos++;
+        }
+    }
+
+    if (nPos == 4)
+    {
+        aboveOp(tet);
+    }
+    else if (nPos == 3)
+    {
+        // Sliced into below tet and above prism. Prism gets split into
+        // two tets.
+
+        // Find the below tet
+        label i0 = -1;
+        forAll(d, i)
+        {
+            if (d[i] <= 0)
+            {
+                i0 = i;
+                break;
+            }
+        }
+
+        label i1 = d.fcIndex(i0);
+        label i2 = d.fcIndex(i1);
+        label i3 = d.fcIndex(i2);
+
+        point p01 = planeIntersection(d, tet, i0, i1);
+        point p02 = planeIntersection(d, tet, i0, i2);
+        point p03 = planeIntersection(d, tet, i0, i3);
+
+        // i0 = tetCell vertex 0: p01,p02,p03 outwards pointing triad
+        //          ,,         1 :     ,,     inwards pointing triad
+        //          ,,         2 :     ,,     outwards pointing triad
+        //          ,,         3 :     ,,     inwards pointing triad
+
+        //Pout<< "Split 3pos tet " << tet << " d:" << d << " into" << nl;
+
+        if (i0 == 0 || i0 == 2)
+        {
+            tetPoints t(tet[i0], p01, p02, p03);
+            //Pout<< "    belowtet:" << t << " around i0:" << i0 << endl;
+            //checkTet(t, "nPos 3, belowTet i0==0 or 2");
+            belowOp(t);
+
+            // Prism
+            FixedList<point, 6> p;
+            p[0] = tet[i1];
+            p[1] = tet[i3];
+            p[2] = tet[i2];
+            p[3] = p01;
+            p[4] = p03;
+            p[5] = p02;
+            //Pout<< "    aboveprism:" << p << endl;
+            decomposePrism(p, aboveOp);
+        }
+        else
+        {
+            tetPoints t(p01, p02, p03, tet[i0]);
+            //Pout<< "    belowtet:" << t << " around i0:" << i0 << endl;
+            //checkTet(t, "nPos 3, belowTet i0==1 or 3");
+            belowOp(t);
+
+            // Prism
+            FixedList<point, 6> p;
+            p[0] = tet[i3];
+            p[1] = tet[i1];
+            p[2] = tet[i2];
+            p[3] = p03;
+            p[4] = p01;
+            p[5] = p02;
+            //Pout<< "    aboveprism:" << p << endl;
+            decomposePrism(p, aboveOp);
+        }
+    }
+    else if (nPos == 2)
+    {
+        // Tet cut into two prisms. Determine the positive one.
+        label pos0 = -1;
+        label pos1 = -1;
+        label neg0 = -1;
+        label neg1 = -1;
+        forAll(d, i)
+        {
+            if (d[i] > 0)
+            {
+                if (pos0 == -1)
+                {
+                    pos0 = i;
+                }
+                else
+                {
+                    pos1 = i;
+                }
+            }
+            else
+            {
+                if (neg0 == -1)
+                {
+                    neg0 = i;
+                }
+                else
+                {
+                    neg1 = i;
+                }
+            }
+        }
+
+        //Pout<< "Split 2pos tet " << tet << " d:" << d
+        //    << " around pos0:" << pos0 << " pos1:" << pos1
+        //    << " neg0:" << neg0 << " neg1:" << neg1 << " into" << nl;
+
+        const edge posEdge(pos0, pos1);
+
+        if (posEdge == edge(0, 1))
+        {
+            point p02 = planeIntersection(d, tet, 0, 2);
+            point p03 = planeIntersection(d, tet, 0, 3);
+            point p12 = planeIntersection(d, tet, 1, 2);
+            point p13 = planeIntersection(d, tet, 1, 3);
+            // Split the resulting prism
+            {
+                FixedList<point, 6> p;
+                p[0] = tet[0];
+                p[1] = p02;
+                p[2] = p03;
+                p[3] = tet[1];
+                p[4] = p12;
+                p[5] = p13;
+                //Pout<< "    01 aboveprism:" << p << endl;
+                decomposePrism(p, aboveOp);
+            }
+            {
+                FixedList<point, 6> p;
+                p[0] = tet[2];
+                p[1] = p02;
+                p[2] = p12;
+                p[3] = tet[3];
+                p[4] = p03;
+                p[5] = p13;
+                //Pout<< "    01 belowprism:" << p << endl;
+                decomposePrism(p, belowOp);
+            }
+        }
+        else if (posEdge == edge(1, 2))
+        {
+            point p01 = planeIntersection(d, tet, 0, 1);
+            point p13 = planeIntersection(d, tet, 1, 3);
+            point p02 = planeIntersection(d, tet, 0, 2);
+            point p23 = planeIntersection(d, tet, 2, 3);
+            // Split the resulting prism
+            {
+                FixedList<point, 6> p;
+                p[0] = tet[1];
+                p[1] = p01;
+                p[2] = p13;
+                p[3] = tet[2];
+                p[4] = p02;
+                p[5] = p23;
+                //Pout<< "    12 aboveprism:" << p << endl;
+                decomposePrism(p, aboveOp);
+            }
+            {
+                FixedList<point, 6> p;
+                p[0] = tet[3];
+                p[1] = p23;
+                p[2] = p13;
+                p[3] = tet[0];
+                p[4] = p02;
+                p[5] = p01;
+                //Pout<< "    12 belowprism:" << p << endl;
+                decomposePrism(p, belowOp);
+            }
+        }
+        else if (posEdge == edge(2, 0))
+        {
+            point p01 = planeIntersection(d, tet, 0, 1);
+            point p03 = planeIntersection(d, tet, 0, 3);
+            point p12 = planeIntersection(d, tet, 1, 2);
+            point p23 = planeIntersection(d, tet, 2, 3);
+            // Split the resulting prism
+            {
+                FixedList<point, 6> p;
+                p[0] = tet[2];
+                p[1] = p12;
+                p[2] = p23;
+                p[3] = tet[0];
+                p[4] = p01;
+                p[5] = p03;
+                //Pout<< "    20 aboveprism:" << p << endl;
+                decomposePrism(p, aboveOp);
+            }
+            {
+                FixedList<point, 6> p;
+                p[0] = tet[1];
+                p[1] = p12;
+                p[2] = p01;
+                p[3] = tet[3];
+                p[4] = p23;
+                p[5] = p03;
+                //Pout<< "    20 belowprism:" << p << endl;
+                decomposePrism(p, belowOp);
+            }
+        }
+        else if (posEdge == edge(0, 3))
+        {
+            point p01 = planeIntersection(d, tet, 0, 1);
+            point p02 = planeIntersection(d, tet, 0, 2);
+            point p13 = planeIntersection(d, tet, 1, 3);
+            point p23 = planeIntersection(d, tet, 2, 3);
+            // Split the resulting prism
+            {
+                FixedList<point, 6> p;
+                p[0] = tet[3];
+                p[1] = p23;
+                p[2] = p13;
+                p[3] = tet[0];
+                p[4] = p02;
+                p[5] = p01;
+                //Pout<< "    03 aboveprism:" << p << endl;
+                decomposePrism(p, aboveOp);
+            }
+            {
+                FixedList<point, 6> p;
+                p[0] = tet[2];
+                p[1] = p23;
+                p[2] = p02;
+                p[3] = tet[1];
+                p[4] = p13;
+                p[5] = p01;
+                //Pout<< "    03 belowprism:" << p << endl;
+                decomposePrism(p, belowOp);
+            }
+        }
+        else if (posEdge == edge(1, 3))
+        {
+            point p01 = planeIntersection(d, tet, 0, 1);
+            point p12 = planeIntersection(d, tet, 1, 2);
+            point p03 = planeIntersection(d, tet, 0, 3);
+            point p23 = planeIntersection(d, tet, 2, 3);
+            // Split the resulting prism
+            {
+                FixedList<point, 6> p;
+                p[0] = tet[1];
+                p[1] = p12;
+                p[2] = p01;
+                p[3] = tet[3];
+                p[4] = p23;
+                p[5] = p03;
+                //Pout<< "    13 aboveprism:" << p << endl;
+                decomposePrism(p, aboveOp);
+            }
+            {
+                FixedList<point, 6> p;
+                p[0] = tet[2];
+                p[1] = p12;
+                p[2] = p23;
+                p[3] = tet[0];
+                p[4] = p01;
+                p[5] = p03;
+                //Pout<< "    13 belowprism:" << p << endl;
+                decomposePrism(p, belowOp);
+            }
+        }
+        else if (posEdge == edge(2, 3))
+        {
+            point p02 = planeIntersection(d, tet, 0, 2);
+            point p12 = planeIntersection(d, tet, 1, 2);
+            point p03 = planeIntersection(d, tet, 0, 3);
+            point p13 = planeIntersection(d, tet, 1, 3);
+            // Split the resulting prism
+            {
+                FixedList<point, 6> p;
+                p[0] = tet[2];
+                p[1] = p02;
+                p[2] = p12;
+                p[3] = tet[3];
+                p[4] = p03;
+                p[5] = p13;
+                //Pout<< "    23 aboveprism:" << p << endl;
+                decomposePrism(p, aboveOp);
+            }
+            {
+                FixedList<point, 6> p;
+                p[0] = tet[0];
+                p[1] = p02;
+                p[2] = p03;
+                p[3] = tet[1];
+                p[4] = p12;
+                p[5] = p13;
+                //Pout<< "    23 belowprism:" << p << endl;
+                decomposePrism(p, belowOp);
+            }
+        }
+        else
+        {
+            FatalErrorIn("tetSliceWithPlane(..)")
+                << "Missed edge:" << posEdge
+                << abort(FatalError);
+        }
+    }
+    else if (nPos == 1)
+    {
+        // Find the positive tet
+        label i0 = -1;
+        forAll(d, i)
+        {
+            if (d[i] > 0)
+            {
+                i0 = i;
+                break;
+            }
+        }
+
+        label i1 = d.fcIndex(i0);
+        label i2 = d.fcIndex(i1);
+        label i3 = d.fcIndex(i2);
+
+        point p01 = planeIntersection(d, tet, i0, i1);
+        point p02 = planeIntersection(d, tet, i0, i2);
+        point p03 = planeIntersection(d, tet, i0, i3);
+
+        //Pout<< "Split 1pos tet " << tet << " d:" << d << " into" << nl;
+
+        if (i0 == 0 || i0 == 2)
+        {
+            tetPoints t(tet[i0], p01, p02, p03);
+            //Pout<< "    abovetet:" << t << " around i0:" << i0 << endl;
+            //checkTet(t, "nPos 1, aboveTets i0==0 or 2");
+            aboveOp(t);
+
+            // Prism
+            FixedList<point, 6> p;
+            p[0] = tet[i1];
+            p[1] = tet[i3];
+            p[2] = tet[i2];
+            p[3] = p01;
+            p[4] = p03;
+            p[5] = p02;
+            //Pout<< "    belowprism:" << p << endl;
+            decomposePrism(p, belowOp);
+        }
+        else
+        {
+            tetPoints t(p01, p02, p03, tet[i0]);
+            //Pout<< "    abovetet:" << t << " around i0:" << i0 << endl;
+            //checkTet(t, "nPos 1, aboveTets i0==1 or 3");
+            aboveOp(t);
+
+            // Prism
+            FixedList<point, 6> p;
+            p[0] = tet[i3];
+            p[1] = tet[i1];
+            p[2] = tet[i2];
+            p[3] = p03;
+            p[4] = p01;
+            p[5] = p02;
+            //Pout<< "    belowprism:" << p << endl;
+            decomposePrism(p, belowOp);
+        }
+    }
+    else    // nPos == 0
+    {
+        belowOp(tet);
+    }
+}
+
+
+template<class Point, class PointRef>
+template<class AboveTetOp, class BelowTetOp>
+inline void Foam::tetrahedron<Point, PointRef>::sliceWithPlane
+(
+    const plane& pl,
+    AboveTetOp& aboveOp,
+    BelowTetOp& belowOp
+) const
+{
+    tetSliceWithPlane(pl, tetPoints(a_, b_, c_, d_), aboveOp, belowOp);
+}
+
+
 // * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
 
 template<class Point, class PointRef>

src/OpenFOAM/primitives/functions/DataEntry/TableFile/TableFile.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

src/Pstream/dummy/UPstream.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -74,6 +74,10 @@ void Foam::UPstream::waitRequests(const label start)
 {}
 
 
+void Foam::UPstream::waitRequest(const label i)
+{}
+
+
 bool Foam::UPstream::finishedRequest(const label i)
 {
     notImplemented("UPstream::finishedRequest()");

src/Pstream/mpi/UPstream.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -516,11 +516,55 @@ void Foam::UPstream::waitRequests(const label start)
 }
 
 
+void Foam::UPstream::waitRequest(const label i)
+{
+    if (debug)
+    {
+        Pout<< "UPstream::waitRequest : starting wait for request:" << i
+            << endl;
+    }
+
+    if (i >= PstreamGlobals::outstandingRequests_.size())
+    {
+        FatalErrorIn
+        (
+            "UPstream::waitRequest(const label)"
+        )   << "There are " << PstreamGlobals::outstandingRequests_.size()
+            << " outstanding send requests and you are asking for i=" << i
+            << nl
+            << "Maybe you are mixing blocking/non-blocking comms?"
+            << Foam::abort(FatalError);
+    }
+
+    int flag;
+    if
+    (
+        MPI_Wait
+        (
+           &PstreamGlobals::outstandingRequests_[i],
+            MPI_STATUS_IGNORE
+        )
+    )
+    {
+        FatalErrorIn
+        (
+            "UPstream::waitRequest()"
+        )   << "MPI_Wait returned with error" << Foam::endl;
+    }
+
+    if (debug)
+    {
+        Pout<< "UPstream::waitRequest : finished wait for request:" << i
+            << endl;
+    }
+}
+
+
 bool Foam::UPstream::finishedRequest(const label i)
 {
     if (debug)
     {
-        Pout<< "UPstream::waitRequests : starting wait for request:" << i
+        Pout<< "UPstream::waitRequests : checking finishedRequest:" << i
             << endl;
     }
 
@@ -546,7 +590,7 @@ bool Foam::UPstream::finishedRequest(const label i)
 
     if (debug)
     {
-        Pout<< "UPstream::waitRequests : finished wait for request:" << i
+        Pout<< "UPstream::waitRequests : finished finishedRequest:" << i
             << endl;
     }
 
@@ -554,6 +598,4 @@ bool Foam::UPstream::finishedRequest(const label i)
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 // ************************************************************************* //

src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,7 +26,7 @@ License
 #include "polyMeshGeometry.H"
 #include "polyMeshTetDecomposition.H"
 #include "pyramidPointFaceRef.H"
-#include "tetPointRef.H"
+#include "tetrahedron.H"
 #include "syncTools.H"
 #include "unitConversion.H"
 

src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C

@@ -756,6 +756,66 @@ void Foam::extendedFeatureEdgeMesh::nearestFeatureEdgeByType
 }
 
 
+void Foam::extendedFeatureEdgeMesh::allNearestFeatureEdges
+(
+    const point& sample,
+    const scalar searchRadiusSqr,
+    List<pointIndexHit>& info
+) const
+{
+    const PtrList<indexedOctree<treeDataEdge> >& edgeTrees = edgeTreesByType();
+
+    info.setSize(edgeTrees.size());
+
+    labelList sliceStarts(edgeTrees.size());
+
+    sliceStarts[0] = externalStart_;
+    sliceStarts[1] = internalStart_;
+    sliceStarts[2] = flatStart_;
+    sliceStarts[3] = openStart_;
+    sliceStarts[4] = multipleStart_;
+
+    DynamicList<pointIndexHit> dynEdgeHit;
+
+    // Loop over all the feature edge types
+    forAll(edgeTrees, i)
+    {
+        // Pick up all the edges that intersect the search sphere
+        labelList elems = edgeTrees[i].findSphere
+        (
+            sample,
+            searchRadiusSqr
+        );
+
+        forAll(elems, elemI)
+        {
+            label index = elems[elemI];
+            label edgeI = edgeTrees[i].shapes().edgeLabels()[index];
+            const edge& e = edges()[edgeI];
+
+            pointHit hitPoint = e.line(points()).nearestDist(sample);
+
+            label hitIndex = index + sliceStarts[i];
+
+            pointIndexHit nearHit;
+
+            if (!hitPoint.hit())
+            {
+                nearHit = pointIndexHit(true, hitPoint.missPoint(), hitIndex);
+            }
+            else
+            {
+                nearHit = pointIndexHit(true, hitPoint.hitPoint(), hitIndex);
+            }
+
+            dynEdgeHit.append(nearHit);
+        }
+    }
+
+    info.transfer(dynEdgeHit);
+}
+
+
 const Foam::indexedOctree<Foam::treeDataEdge>&
 Foam::extendedFeatureEdgeMesh::edgeTree() const
 {

src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.H

@@ -271,6 +271,14 @@ public:
                 List<pointIndexHit>& info
             ) const;
 
+            //- Find all the feature edges within searchDistSqr of sample
+            void allNearestFeatureEdges
+            (
+                const point& sample,
+                const scalar searchRadiusSqr,
+                List<pointIndexHit>& info
+            ) const;
+
         // Access
 
             //- Return the index of the start of the convex feature points

src/finiteVolume/fields/fvPatchFields/basic/sliced/slicedFvPatchField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -198,6 +198,17 @@ tmp<Field<Type> > slicedFvPatchField<Type>::patchInternalField() const
 }
 
 
+template<class Type>
+void slicedFvPatchField<Type>::patchInternalField(Field<Type>&) const
+{
+    notImplemented
+    (
+        "slicedFvPatchField<Type>::"
+        "patchInternalField(Field<Type>&)"
+    );
+}
+
+
 template<class Type>
 tmp<Field<Type> > slicedFvPatchField<Type>::patchNeighbourField
 (

src/finiteVolume/fields/fvPatchFields/basic/sliced/slicedFvPatchField.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -148,6 +148,9 @@ public:
             //- Return internal field next to patch as patch field
             virtual tmp<Field<Type> > patchInternalField() const;
 
+            //- Return internal field next to patch as patch field
+            virtual void patchInternalField(Field<Type>&) const;
+
             //- Return neighbour coupled given internal cell data
             virtual tmp<Field<Type> > patchNeighbourField
             (

src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -43,7 +43,12 @@ processorFvPatchField<Type>::processorFvPatchField
 )
 :
     coupledFvPatchField<Type>(p, iF),
-    procPatch_(refCast<const processorFvPatch>(p))
+    procPatch_(refCast<const processorFvPatch>(p)),
+    sendBuf_(0),
+    outstandingSendRequest_(-1),
+    outstandingRecvRequest_(-1),
+    scalarSendBuf_(0),
+    scalarReceiveBuf_(0)
 {}
 
 
@@ -56,7 +61,12 @@ processorFvPatchField<Type>::processorFvPatchField
 )
 :
     coupledFvPatchField<Type>(p, iF, f),
-    procPatch_(refCast<const processorFvPatch>(p))
+    procPatch_(refCast<const processorFvPatch>(p)),
+    sendBuf_(0),
+    outstandingSendRequest_(-1),
+    outstandingRecvRequest_(-1),
+    scalarSendBuf_(0),
+    scalarReceiveBuf_(0)
 {}
 
 
@@ -71,7 +81,12 @@ processorFvPatchField<Type>::processorFvPatchField
 )
 :
     coupledFvPatchField<Type>(ptf, p, iF, mapper),
-    procPatch_(refCast<const processorFvPatch>(p))
+    procPatch_(refCast<const processorFvPatch>(p)),
+    sendBuf_(0),
+    outstandingSendRequest_(-1),
+    outstandingRecvRequest_(-1),
+    scalarSendBuf_(0),
+    scalarReceiveBuf_(0)
 {
     if (!isA<processorFvPatch>(this->patch()))
     {
@@ -91,6 +106,12 @@ processorFvPatchField<Type>::processorFvPatchField
             << " in file " << this->dimensionedInternalField().objectPath()
             << exit(FatalIOError);
     }
+    if (debug && !ptf.ready())
+    {
+        FatalErrorIn("processorFvPatchField<Type>::processorFvPatchField(..)")
+            << "On patch " << procPatch_.name() << " outstanding request."
+            << abort(FatalError);
+    }
 }
 
 
@@ -103,7 +124,12 @@ processorFvPatchField<Type>::processorFvPatchField
 )
 :
     coupledFvPatchField<Type>(p, iF, dict),
-    procPatch_(refCast<const processorFvPatch>(p))
+    procPatch_(refCast<const processorFvPatch>(p)),
+    sendBuf_(0),
+    outstandingSendRequest_(-1),
+    outstandingRecvRequest_(-1),
+    scalarSendBuf_(0),
+    scalarReceiveBuf_(0)
 {
     if (!isA<processorFvPatch>(p))
     {
@@ -134,8 +160,20 @@ processorFvPatchField<Type>::processorFvPatchField
 :
     processorLduInterfaceField(),
     coupledFvPatchField<Type>(ptf),
-    procPatch_(refCast<const processorFvPatch>(ptf.patch()))
+    procPatch_(refCast<const processorFvPatch>(ptf.patch())),
-{}
+    sendBuf_(ptf.sendBuf_.xfer()),
+    outstandingSendRequest_(-1),
+    outstandingRecvRequest_(-1),
+    scalarSendBuf_(ptf.scalarSendBuf_.xfer()),
+    scalarReceiveBuf_(ptf.scalarReceiveBuf_.xfer())
+{
+    if (debug && !ptf.ready())
+    {
+        FatalErrorIn("processorFvPatchField<Type>::processorFvPatchField(..)")
+            << "On patch " << procPatch_.name() << " outstanding request."
+            << abort(FatalError);
+    }
+}
 
 
 template<class Type>
@@ -146,8 +184,20 @@ processorFvPatchField<Type>::processorFvPatchField
 )
 :
     coupledFvPatchField<Type>(ptf, iF),
-    procPatch_(refCast<const processorFvPatch>(ptf.patch()))
+    procPatch_(refCast<const processorFvPatch>(ptf.patch())),
-{}
+    sendBuf_(0),
+    outstandingSendRequest_(-1),
+    outstandingRecvRequest_(-1),
+    scalarSendBuf_(0),
+    scalarReceiveBuf_(0)
+{
+    if (debug && !ptf.ready())
+    {
+        FatalErrorIn("processorFvPatchField<Type>::processorFvPatchField(..)")
+            << "On patch " << procPatch_.name() << " outstanding request."
+            << abort(FatalError);
+    }
+}
 
 
 // * * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * //
@@ -162,6 +212,13 @@ processorFvPatchField<Type>::~processorFvPatchField()
 template<class Type>
 tmp<Field<Type> > processorFvPatchField<Type>::patchNeighbourField() const
 {
+    if (debug && !this->ready())
+    {
+        FatalErrorIn("processorFvPatchField<Type>::patchNeighbourField()")
+            << "On patch " << procPatch_.name()
+            << " outstanding request."
+            << abort(FatalError);
+    }
     return *this;
 }
 
@@ -174,7 +231,36 @@ void processorFvPatchField<Type>::initEvaluate
 {
     if (Pstream::parRun())
     {
-        procPatch_.compressedSend(commsType, this->patchInternalField()());
+        this->patchInternalField(sendBuf_);
+
+        if (commsType == Pstream::nonBlocking && !Pstream::floatTransfer)
+        {
+            // Fast path. Receive into *this
+            this->setSize(sendBuf_.size());
+            outstandingRecvRequest_ = UPstream::nRequests();
+            IPstream::read
+            (
+                Pstream::nonBlocking,
+                procPatch_.neighbProcNo(),
+                reinterpret_cast<char*>(this->begin()),
+                this->byteSize(),
+                procPatch_.tag()
+            );
+
+            outstandingSendRequest_ = UPstream::nRequests();
+            OPstream::write
+            (
+                Pstream::nonBlocking,
+                procPatch_.neighbProcNo(),
+                reinterpret_cast<const char*>(sendBuf_.begin()),
+                this->byteSize(),
+                procPatch_.tag()
+            );
+        }
+        else
+        {
+            procPatch_.compressedSend(commsType, sendBuf_);
+        }
     }
 }
 
@@ -187,7 +273,25 @@ void processorFvPatchField<Type>::evaluate
 {
     if (Pstream::parRun())
     {
-        procPatch_.compressedReceive<Type>(commsType, *this);
+        if (commsType == Pstream::nonBlocking && !Pstream::floatTransfer)
+        {
+            // Fast path. Received into *this
+
+            if
+            (
+                outstandingRecvRequest_ >= 0
+             && outstandingRecvRequest_ < Pstream::nRequests()
+            )
+            {
+                UPstream::waitRequest(outstandingRecvRequest_);
+            }
+            outstandingSendRequest_ = -1;
+            outstandingRecvRequest_ = -1;
+        }
+        else
+        {
+            procPatch_.compressedReceive<Type>(commsType, *this);
+        }
 
         if (doTransform())
         {
@@ -215,11 +319,49 @@ void processorFvPatchField<Type>::initInterfaceMatrixUpdate
     const Pstream::commsTypes commsType
 ) const
 {
-    procPatch_.compressedSend
+    this->patch().patchInternalField(psiInternal, scalarSendBuf_);
-    (
+
-        commsType,
+    if (commsType == Pstream::nonBlocking && !Pstream::floatTransfer)
-        this->patch().patchInternalField(psiInternal)()
+    {
-    );
+        // Fast path.
+        if (debug && !this->ready())
+        {
+            FatalErrorIn
+            (
+                "processorFvPatchField<Type>::initInterfaceMatrixUpdate(..)"
+            )   << "On patch " << procPatch_.name()
+                << " outstanding request."
+                << abort(FatalError);
+        }
+
+
+        scalarReceiveBuf_.setSize(scalarSendBuf_.size());
+        outstandingRecvRequest_ = UPstream::nRequests();
+        IPstream::read
+        (
+            Pstream::nonBlocking,
+            procPatch_.neighbProcNo(),
+            reinterpret_cast<char*>(scalarReceiveBuf_.begin()),
+            scalarReceiveBuf_.byteSize(),
+            procPatch_.tag()
+        );
+
+        outstandingSendRequest_ = UPstream::nRequests();
+        OPstream::write
+        (
+            Pstream::nonBlocking,
+            procPatch_.neighbProcNo(),
+            reinterpret_cast<const char*>(scalarSendBuf_.begin()),
+            scalarSendBuf_.byteSize(),
+            procPatch_.tag()
+        );
+    }
+    else
+    {
+        procPatch_.compressedSend(commsType, scalarSendBuf_);
+    }
+
+    const_cast<processorFvPatchField<Type>&>(*this).updatedMatrix() = false;
 }
 
 
@@ -234,22 +376,92 @@ void processorFvPatchField<Type>::updateInterfaceMatrix
     const Pstream::commsTypes commsType
 ) const
 {
-    scalarField pnf
+    if (this->updatedMatrix())
-    (
+    {
-        procPatch_.compressedReceive<scalar>(commsType, this->size())()
+        return;
-    );
+    }
-
-    // Transform according to the transformation tensor
-    transformCoupleField(pnf, cmpt);
-
-    // Multiply the field by coefficients and add into the result
 
     const labelUList& faceCells = this->patch().faceCells();
 
-    forAll(faceCells, elemI)
+    if (commsType == Pstream::nonBlocking && !Pstream::floatTransfer)
     {
-        result[faceCells[elemI]] -= coeffs[elemI]*pnf[elemI];
+        // Fast path.
+        if
+        (
+            outstandingRecvRequest_ >= 0
+         && outstandingRecvRequest_ < Pstream::nRequests()
+        )
+        {
+            UPstream::waitRequest(outstandingRecvRequest_);
+        }
+        // Recv finished so assume sending finished as well.
+        outstandingSendRequest_ = -1;
+        outstandingRecvRequest_ = -1;
+
+        // Consume straight from scalarReceiveBuf_
+
+        // Transform according to the transformation tensor
+        transformCoupleField(scalarReceiveBuf_, cmpt);
+
+        // Multiply the field by coefficients and add into the result
+        forAll(faceCells, elemI)
+        {
+            result[faceCells[elemI]] -= coeffs[elemI]*scalarReceiveBuf_[elemI];
+        }
     }
+    else
+    {
+        scalarField pnf
+        (
+            procPatch_.compressedReceive<scalar>(commsType, this->size())()
+        );
+
+        // Transform according to the transformation tensor
+        transformCoupleField(pnf, cmpt);
+
+        // Multiply the field by coefficients and add into the result
+        forAll(faceCells, elemI)
+        {
+            result[faceCells[elemI]] -= coeffs[elemI]*pnf[elemI];
+        }
+    }
+
+    const_cast<processorFvPatchField<Type>&>(*this).updatedMatrix() = true;
+}
+
+
+template<class Type>
+bool processorFvPatchField<Type>::ready() const
+{
+    if
+    (
+        outstandingSendRequest_ >= 0
+     && outstandingSendRequest_ < Pstream::nRequests()
+    )
+    {
+        bool finished = UPstream::finishedRequest(outstandingSendRequest_);
+        if (!finished)
+        {
+            return false;
+        }
+    }
+    outstandingSendRequest_ = -1;
+
+    if
+    (
+        outstandingRecvRequest_ >= 0
+     && outstandingRecvRequest_ < Pstream::nRequests()
+    )
+    {
+        bool finished = UPstream::finishedRequest(outstandingRecvRequest_);
+        if (!finished)
+        {
+            return false;
+        }
+    }
+    outstandingRecvRequest_ = -1;
+
+    return true;
 }
 
 

src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchField.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -59,6 +59,22 @@ class processorFvPatchField
         //- Local reference cast into the processor patch
         const processorFvPatch& procPatch_;
 
+        // Sending and receiving
+
+            //- Send buffer.
+            mutable Field<Type> sendBuf_;
+
+            //- Outstanding request
+            mutable label outstandingSendRequest_;
+
+            //- Outstanding request
+            mutable label outstandingRecvRequest_;
+
+            //- Scalar send buffer
+            mutable Field<scalar> scalarSendBuf_;
+
+            //- Scalar receive buffer
+            mutable Field<scalar> scalarReceiveBuf_;
 
 public:
 
@@ -179,6 +195,9 @@ public:
                 const Pstream::commsTypes commsType
             ) const;
 
+            //- Is all data available
+            virtual bool ready() const;
+
             //- Update result field based on interface functionality
             virtual void updateInterfaceMatrix
             (
@@ -221,6 +240,7 @@ public:
             {
                 return pTraits<Type>::rank;
             }
+
 };
 
 

src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchScalarField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -43,11 +43,53 @@ void processorFvPatchField<scalar>::initInterfaceMatrixUpdate
     const Pstream::commsTypes commsType
 ) const
 {
-    procPatch_.compressedSend
+    this->patch().patchInternalField(psiInternal, scalarSendBuf_);
+
+    if
     (
-        commsType,
+        Pstream::defaultCommsType == Pstream::nonBlocking
-        patch().patchInternalField(psiInternal)()
+    && !Pstream::floatTransfer
-    );
+    )
+    {
+        // Fast path.
+        if (debug && !this->ready())
+        {
+            FatalErrorIn
+            (
+                "processorFvPatchField<scalar>::initInterfaceMatrixUpdate(..)"
+            )   << "On patch " << procPatch_.name()
+                << " outstanding request."
+                << abort(FatalError);
+        }
+
+
+        scalarReceiveBuf_.setSize(scalarSendBuf_.size());
+        outstandingRecvRequest_ = UPstream::nRequests();
+        IPstream::read
+        (
+            Pstream::nonBlocking,
+            procPatch_.neighbProcNo(),
+            reinterpret_cast<char*>(scalarReceiveBuf_.begin()),
+            scalarReceiveBuf_.byteSize(),
+            procPatch_.tag()
+        );
+
+        outstandingSendRequest_ = UPstream::nRequests();
+        OPstream::write
+        (
+            Pstream::nonBlocking,
+            procPatch_.neighbProcNo(),
+            reinterpret_cast<const char*>(scalarSendBuf_.begin()),
+            scalarSendBuf_.byteSize(),
+            procPatch_.tag()
+        );
+    }
+    else
+    {
+        procPatch_.compressedSend(commsType, scalarSendBuf_);
+    }
+
+    const_cast<processorFvPatchField<scalar>&>(*this).updatedMatrix() = false;
 }
 
 
@@ -62,17 +104,48 @@ void processorFvPatchField<scalar>::updateInterfaceMatrix
     const Pstream::commsTypes commsType
 ) const
 {
-    scalarField pnf
+    if (this->updatedMatrix())
-    (
-        procPatch_.compressedReceive<scalar>(commsType, this->size())()
-    );
-
-    const labelUList& faceCells = patch().faceCells();
-
-    forAll(faceCells, facei)
     {
-        result[faceCells[facei]] -= coeffs[facei]*pnf[facei];
+        return;
     }
+
+    const labelUList& faceCells = this->patch().faceCells();
+
+    if (commsType == Pstream::nonBlocking && !Pstream::floatTransfer)
+    {
+        // Fast path.
+        if
+        (
+            outstandingRecvRequest_ >= 0
+         && outstandingRecvRequest_ < Pstream::nRequests()
+        )
+        {
+            UPstream::waitRequest(outstandingRecvRequest_);
+        }
+        outstandingSendRequest_ = -1;
+        outstandingRecvRequest_ = -1;
+
+
+        // Consume straight from scalarReceiveBuf_
+        forAll(faceCells, elemI)
+        {
+            result[faceCells[elemI]] -= coeffs[elemI]*scalarReceiveBuf_[elemI];
+        }
+    }
+    else
+    {
+        scalarField pnf
+        (
+            procPatch_.compressedReceive<scalar>(commsType, this->size())()
+        );
+
+        forAll(faceCells, elemI)
+        {
+            result[faceCells[elemI]] -= coeffs[elemI]*pnf[elemI];
+        }
+    }
+
+    const_cast<processorFvPatchField<scalar>&>(*this).updatedMatrix() = true;
 }
 
 

src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -190,6 +190,13 @@ Foam::fvPatchField<Type>::patchInternalField() const
 }
 
 
+template<class Type>
+void Foam::fvPatchField<Type>::patchInternalField(Field<Type>& pif) const
+{
+    patch_.patchInternalField(internalField_, pif);
+}
+
+
 template<class Type>
 void Foam::fvPatchField<Type>::autoMap
 (

src/finiteVolume/fields/fvPatchFields/fvPatchField/fvPatchField.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -352,6 +352,9 @@ public:
             //- Return internal field next to patch as patch field
             virtual tmp<Field<Type> > patchInternalField() const;
 
+            //- Return internal field next to patch as patch field
+            virtual void patchInternalField(Field<Type>&) const;
+
             //- Return patchField on the opposite patch of a coupled patch
             virtual tmp<Field<Type> > patchNeighbourField() const
             {

src/finiteVolume/fvMesh/fvPatches/fvPatch/fvPatch.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -239,6 +239,10 @@ public:
             template<class Type>
             tmp<Field<Type> > patchInternalField(const UList<Type>&) const;
 
+            //- Return given internal field next to patch as patch field
+            template<class Type>
+            void patchInternalField(const UList<Type>&, Field<Type>&) const;
+
             //- Return the corresponding patchField of the named field
             template<class GeometricField, class Type>
             const typename GeometricField::PatchFieldType& patchField

src/finiteVolume/fvMesh/fvPatches/fvPatch/fvPatchTemplates.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -47,6 +47,24 @@ Foam::tmp<Foam::Field<Type> > Foam::fvPatch::patchInternalField
 }
 
 
+template<class Type>
+void Foam::fvPatch::patchInternalField
+(
+    const UList<Type>& f,
+    Field<Type>& pif
+) const
+{
+    pif.setSize(size());
+
+    const labelUList& faceCells = this->faceCells();
+
+    forAll(pif, facei)
+    {
+        pif[facei] = f[faceCells[facei]];
+    }
+}
+
+
 template<class GeometricField, class Type>
 const typename GeometricField::PatchFieldType& Foam::fvPatch::patchField
 (

src/finiteVolume/interpolation/interpolation/interpolationCellPoint/cellPointWeight/cellPointWeight.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,7 +25,6 @@ License
 
 #include "cellPointWeight.H"
 #include "polyMesh.H"
-#include "tetPointRef.H"
 #include "polyMeshTetDecomposition.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

src/lagrangian/basic/Cloud/Cloud.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -41,9 +41,6 @@ SourceFiles
 #include "IOField.H"
 #include "CompactIOField.H"
 #include "polyMesh.H"
-#include "indexedOctree.H"
-#include "treeDataCell.H"
-#include "tetPointRef.H"
 #include "PackedBoolList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/lagrangian/basic/particle/particle.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,11 +35,10 @@ Description
 #include "vector.H"
 #include "Cloud.H"
 #include "IDLList.H"
-#include "labelList.H"
 #include "pointField.H"
 #include "faceList.H"
 #include "OFstream.H"
-#include "tetPointRef.H"
+#include "tetrahedron.H"
 #include "FixedList.H"
 #include "polyMeshTetDecomposition.H"
 

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C

@@ -1144,6 +1144,7 @@ void Foam::autoLayerDriver::syncPatchDisplacement
         );
 
         // Reset if differs
+        // 1. take max
         forAll(syncPatchNLayers, i)
         {
             if (syncPatchNLayers[i] != patchNLayers[i])
@@ -1174,6 +1175,7 @@ void Foam::autoLayerDriver::syncPatchDisplacement
         );
 
         // Reset if differs
+        // 2. take min
         forAll(syncPatchNLayers, i)
         {
             if (syncPatchNLayers[i] != patchNLayers[i])
@@ -2726,7 +2728,7 @@ void Foam::autoLayerDriver::addLayers
             Info<< "Writing shrunk mesh to " << meshRefiner_.timeName() << endl;
 
             // See comment in autoSnapDriver why we should not remove meshPhi
-            // using mesh.clearPout().
+            // using mesh.clearOut().
 
             meshRefiner_.write
             (

src/meshTools/indexedOctree/treeDataEdge.C

@@ -128,14 +128,37 @@ bool Foam::treeDataEdge::overlaps
     const treeBoundBox& cubeBb
 ) const
 {
-    if (cacheBb_)
+    const edge& e = edges_[edgeLabels_[index]];
+
+    const point& start = points_[e.start()];
+    const point& end = points_[e.end()];
+
+    point intersect;
+
+    return cubeBb.intersects(start, end, intersect);
+}
+
+
+// Check if any point on shape is inside sphere.
+bool Foam::treeDataEdge::overlaps
+(
+    const label index,
+    const point& centre,
+    const scalar radiusSqr
+) const
+{
+    const edge& e = edges_[edgeLabels_[index]];
+
+    const pointHit nearHit = e.line(points_).nearestDist(centre);
+
+    const scalar distSqr = sqr(nearHit.distance());
+
+    if (distSqr <= radiusSqr)
     {
-        return cubeBb.overlaps(bbs_[index]);
+        return true;
-    }
-    else
-    {
-        return cubeBb.overlaps(calcBb(edgeLabels_[index]));
     }
+
+    return false;
 }
 
 

src/meshTools/indexedOctree/treeDataEdge.H

@@ -148,6 +148,14 @@ public:
                 const treeBoundBox& sampleBb
             ) const;
 
+            //- Does (bb of) shape at index overlap bb
+            bool overlaps
+            (
+                const label index,
+                const point& centre,
+                const scalar radiusSqr
+            ) const;
+
             //- Calculates nearest (to sample) point in shape.
             //  Returns actual point and distance (squared)
             void findNearest

src/meshTools/momentOfInertia/momentOfInertia.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,7 +37,6 @@ SourceFiles
 #ifndef momentOfInertia_H
 #define momentOfInertia_H
 
-#include "tetPointRef.H"
 #include "triFaceList.H"
 #include "triSurface.H"
 #include "polyMesh.H"

src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,436 +24,27 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "tetOverlapVolume.H"
+#include "tetrahedron.H"
+#include "tetPoints.H"
+#include "polyMesh.H"
+#include "OFstream.H"
+#include "treeBoundBox.H"
+#include "indexedOctree.H"
+#include "treeDataCell.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 defineTypeNameAndDebug(Foam::tetOverlapVolume, 0);
 
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::tetOverlapVolume::tetOverlapVolume()
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::tetOverlapVolume::~tetOverlapVolume()
-{}
-
-
 // * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * * //
 
-Foam::point Foam::tetOverlapVolume::planeIntersection
-(
-    const FixedList<scalar, 4>& d,
-    const tetPoints& t,
-    const label negI,
-    const label posI
-)
-{
-    return (d[posI]*t[negI] - d[negI]*t[posI])/(-d[negI]+d[posI]);
-}
-
-
-template <class TetOp>
-inline void Foam::tetOverlapVolume::decomposePrism
-(
-    const FixedList<point, 6>& points,
-    TetOp& op
-)
-{
-    op(tetPoints(points[1], points[3], points[2], points[0]));
-    op(tetPoints(points[1], points[2], points[3], points[4]));
-    op(tetPoints(points[4], points[2], points[3], points[5]));
-}
-
-
-template <class AboveTetOp, class BelowTetOp>
-inline void Foam::tetOverlapVolume::tetSliceWithPlane
-(
-    const tetPoints& tet,
-    const plane& pl,
-
-    AboveTetOp& aboveOp,
-    BelowTetOp& belowOp
-)
-{
-    // Distance to plane
-    FixedList<scalar, 4> d;
-    label nPos = 0;
-    forAll(tet, i)
-    {
-        d[i] = ((tet[i]-pl.refPoint()) & pl.normal());
-        if (d[i] > 0)
-        {
-            nPos++;
-        }
-    }
-
-    if (nPos == 4)
-    {
-        aboveOp(tet);
-    }
-    else if (nPos == 3)
-    {
-        // Sliced into below tet and above prism. Prism gets split into
-        // two tets.
-
-        // Find the below tet
-        label i0 = -1;
-        forAll(d, i)
-        {
-            if (d[i] <= 0)
-            {
-                i0 = i;
-                break;
-            }
-        }
-
-        label i1 = d.fcIndex(i0);
-        label i2 = d.fcIndex(i1);
-        label i3 = d.fcIndex(i2);
-
-        point p01 = planeIntersection(d, tet, i0, i1);
-        point p02 = planeIntersection(d, tet, i0, i2);
-        point p03 = planeIntersection(d, tet, i0, i3);
-
-        // i0 = tetCell vertex 0: p01,p02,p03 outwards pointing triad
-        //          ,,         1 :     ,,     inwards pointing triad
-        //          ,,         2 :     ,,     outwards pointing triad
-        //          ,,         3 :     ,,     inwards pointing triad
-
-        //Pout<< "Split 3pos tet " << tet << " d:" << d << " into" << nl;
-
-        if (i0 == 0 || i0 == 2)
-        {
-            tetPoints t(tet[i0], p01, p02, p03);
-            //Pout<< "    belowtet:" << t << " around i0:" << i0 << endl;
-            //checkTet(t, "nPos 3, belowTet i0==0 or 2");
-            belowOp(t);
-
-            // Prism
-            FixedList<point, 6> p;
-            p[0] = tet[i1];
-            p[1] = tet[i3];
-            p[2] = tet[i2];
-            p[3] = p01;
-            p[4] = p03;
-            p[5] = p02;
-            //Pout<< "    aboveprism:" << p << endl;
-            decomposePrism(p, aboveOp);
-        }
-        else
-        {
-            tetPoints t(p01, p02, p03, tet[i0]);
-            //Pout<< "    belowtet:" << t << " around i0:" << i0 << endl;
-            //checkTet(t, "nPos 3, belowTet i0==1 or 3");
-            belowOp(t);
-
-            // Prism
-            FixedList<point, 6> p;
-            p[0] = tet[i3];
-            p[1] = tet[i1];
-            p[2] = tet[i2];
-            p[3] = p03;
-            p[4] = p01;
-            p[5] = p02;
-            //Pout<< "    aboveprism:" << p << endl;
-            decomposePrism(p, aboveOp);
-        }
-    }
-    else if (nPos == 2)
-    {
-        // Tet cut into two prisms. Determine the positive one.
-        label pos0 = -1;
-        label pos1 = -1;
-        label neg0 = -1;
-        label neg1 = -1;
-        forAll(d, i)
-        {
-            if (d[i] > 0)
-            {
-                if (pos0 == -1)
-                {
-                    pos0 = i;
-                }
-                else
-                {
-                    pos1 = i;
-                }
-            }
-            else
-            {
-                if (neg0 == -1)
-                {
-                    neg0 = i;
-                }
-                else
-                {
-                    neg1 = i;
-                }
-            }
-        }
-
-        //Pout<< "Split 2pos tet " << tet << " d:" << d
-        //    << " around pos0:" << pos0 << " pos1:" << pos1
-        //    << " neg0:" << neg0 << " neg1:" << neg1 << " into" << nl;
-
-        const edge posEdge(pos0, pos1);
-
-        if (posEdge == edge(0, 1))
-        {
-            point p02 = planeIntersection(d, tet, 0, 2);
-            point p03 = planeIntersection(d, tet, 0, 3);
-            point p12 = planeIntersection(d, tet, 1, 2);
-            point p13 = planeIntersection(d, tet, 1, 3);
-            // Split the resulting prism
-            {
-                FixedList<point, 6> p;
-                p[0] = tet[0];
-                p[1] = p02;
-                p[2] = p03;
-                p[3] = tet[1];
-                p[4] = p12;
-                p[5] = p13;
-                //Pout<< "    01 aboveprism:" << p << endl;
-                decomposePrism(p, aboveOp);
-            }
-            {
-                FixedList<point, 6> p;
-                p[0] = tet[2];
-                p[1] = p02;
-                p[2] = p12;
-                p[3] = tet[3];
-                p[4] = p03;
-                p[5] = p13;
-                //Pout<< "    01 belowprism:" << p << endl;
-                decomposePrism(p, belowOp);
-            }
-        }
-        else if (posEdge == edge(1, 2))
-        {
-            point p01 = planeIntersection(d, tet, 0, 1);
-            point p13 = planeIntersection(d, tet, 1, 3);
-            point p02 = planeIntersection(d, tet, 0, 2);
-            point p23 = planeIntersection(d, tet, 2, 3);
-            // Split the resulting prism
-            {
-                FixedList<point, 6> p;
-                p[0] = tet[1];
-                p[1] = p01;
-                p[2] = p13;
-                p[3] = tet[2];
-                p[4] = p02;
-                p[5] = p23;
-                //Pout<< "    12 aboveprism:" << p << endl;
-                decomposePrism(p, aboveOp);
-            }
-            {
-                FixedList<point, 6> p;
-                p[0] = tet[3];
-                p[1] = p23;
-                p[2] = p13;
-                p[3] = tet[0];
-                p[4] = p02;
-                p[5] = p01;
-                //Pout<< "    12 belowprism:" << p << endl;
-                decomposePrism(p, belowOp);
-            }
-        }
-        else if (posEdge == edge(2, 0))
-        {
-            point p01 = planeIntersection(d, tet, 0, 1);
-            point p03 = planeIntersection(d, tet, 0, 3);
-            point p12 = planeIntersection(d, tet, 1, 2);
-            point p23 = planeIntersection(d, tet, 2, 3);
-            // Split the resulting prism
-            {
-                FixedList<point, 6> p;
-                p[0] = tet[2];
-                p[1] = p12;
-                p[2] = p23;
-                p[3] = tet[0];
-                p[4] = p01;
-                p[5] = p03;
-                //Pout<< "    20 aboveprism:" << p << endl;
-                decomposePrism(p, aboveOp);
-            }
-            {
-                FixedList<point, 6> p;
-                p[0] = tet[1];
-                p[1] = p12;
-                p[2] = p01;
-                p[3] = tet[3];
-                p[4] = p23;
-                p[5] = p03;
-                //Pout<< "    20 belowprism:" << p << endl;
-                decomposePrism(p, belowOp);
-            }
-        }
-        else if (posEdge == edge(0, 3))
-        {
-            point p01 = planeIntersection(d, tet, 0, 1);
-            point p02 = planeIntersection(d, tet, 0, 2);
-            point p13 = planeIntersection(d, tet, 1, 3);
-            point p23 = planeIntersection(d, tet, 2, 3);
-            // Split the resulting prism
-            {
-                FixedList<point, 6> p;
-                p[0] = tet[3];
-                p[1] = p23;
-                p[2] = p13;
-                p[3] = tet[0];
-                p[4] = p02;
-                p[5] = p01;
-                //Pout<< "    03 aboveprism:" << p << endl;
-                decomposePrism(p, aboveOp);
-            }
-            {
-                FixedList<point, 6> p;
-                p[0] = tet[2];
-                p[1] = p23;
-                p[2] = p02;
-                p[3] = tet[1];
-                p[4] = p13;
-                p[5] = p01;
-                //Pout<< "    03 belowprism:" << p << endl;
-                decomposePrism(p, belowOp);
-            }
-        }
-        else if (posEdge == edge(1, 3))
-        {
-            point p01 = planeIntersection(d, tet, 0, 1);
-            point p12 = planeIntersection(d, tet, 1, 2);
-            point p03 = planeIntersection(d, tet, 0, 3);
-            point p23 = planeIntersection(d, tet, 2, 3);
-            // Split the resulting prism
-            {
-                FixedList<point, 6> p;
-                p[0] = tet[1];
-                p[1] = p12;
-                p[2] = p01;
-                p[3] = tet[3];
-                p[4] = p23;
-                p[5] = p03;
-                //Pout<< "    13 aboveprism:" << p << endl;
-                decomposePrism(p, aboveOp);
-            }
-            {
-                FixedList<point, 6> p;
-                p[0] = tet[2];
-                p[1] = p12;
-                p[2] = p23;
-                p[3] = tet[0];
-                p[4] = p01;
-                p[5] = p03;
-                //Pout<< "    13 belowprism:" << p << endl;
-                decomposePrism(p, belowOp);
-            }
-        }
-        else if (posEdge == edge(2, 3))
-        {
-            point p02 = planeIntersection(d, tet, 0, 2);
-            point p12 = planeIntersection(d, tet, 1, 2);
-            point p03 = planeIntersection(d, tet, 0, 3);
-            point p13 = planeIntersection(d, tet, 1, 3);
-            // Split the resulting prism
-            {
-                FixedList<point, 6> p;
-                p[0] = tet[2];
-                p[1] = p02;
-                p[2] = p12;
-                p[3] = tet[3];
-                p[4] = p03;
-                p[5] = p13;
-                //Pout<< "    23 aboveprism:" << p << endl;
-                decomposePrism(p, aboveOp);
-            }
-            {
-                FixedList<point, 6> p;
-                p[0] = tet[0];
-                p[1] = p02;
-                p[2] = p03;
-                p[3] = tet[1];
-                p[4] = p12;
-                p[5] = p13;
-                //Pout<< "    23 belowprism:" << p << endl;
-                decomposePrism(p, belowOp);
-            }
-        }
-        else
-        {
-            FatalErrorIn("tetSliceWithPlane(..)") << "Missed edge:" << posEdge
-                << abort(FatalError);
-        }
-    }
-    else if (nPos == 1)
-    {
-        // Find the positive tet
-        label i0 = -1;
-        forAll(d, i)
-        {
-            if (d[i] > 0)
-            {
-                i0 = i;
-                break;
-            }
-        }
-
-        label i1 = d.fcIndex(i0);
-        label i2 = d.fcIndex(i1);
-        label i3 = d.fcIndex(i2);
-
-        point p01 = planeIntersection(d, tet, i0, i1);
-        point p02 = planeIntersection(d, tet, i0, i2);
-        point p03 = planeIntersection(d, tet, i0, i3);
-
-        //Pout<< "Split 1pos tet " << tet << " d:" << d << " into" << nl;
-
-        if (i0 == 0 || i0 == 2)
-        {
-            tetPoints t(tet[i0], p01, p02, p03);
-            //Pout<< "    abovetet:" << t << " around i0:" << i0 << endl;
-            //checkTet(t, "nPos 1, aboveTets i0==0 or 2");
-            aboveOp(t);
-
-            // Prism
-            FixedList<point, 6> p;
-            p[0] = tet[i1];
-            p[1] = tet[i3];
-            p[2] = tet[i2];
-            p[3] = p01;
-            p[4] = p03;
-            p[5] = p02;
-            //Pout<< "    belowprism:" << p << endl;
-            decomposePrism(p, belowOp);
-        }
-        else
-        {
-            tetPoints t(p01, p02, p03, tet[i0]);
-            //Pout<< "    abovetet:" << t << " around i0:" << i0 << endl;
-            //checkTet(t, "nPos 1, aboveTets i0==1 or 3");
-            aboveOp(t);
-
-            // Prism
-            FixedList<point, 6> p;
-            p[0] = tet[i3];
-            p[1] = tet[i1];
-            p[2] = tet[i2];
-            p[3] = p03;
-            p[4] = p01;
-            p[5] = p02;
-            //Pout<< "    belowprism:" << p << endl;
-            decomposePrism(p, belowOp);
-        }
-    }
-    else    // nPos == 0
-    {
-        belowOp(tet);
-    }
-}
-
-
 void Foam::tetOverlapVolume::tetTetOverlap
 (
     const tetPoints& tetA,
@@ -462,16 +53,15 @@ void Foam::tetOverlapVolume::tetTetOverlap
     label& nInside,
     FixedList<tetPoints, 200>& outsideTets,
     label& nOutside
-)
+) const
 {
     // Work storage
     FixedList<tetPoints, 200> cutInsideTets;
     label nCutInside = 0;
 
-    storeTetOp inside(insideTets, nInside);
+    tetPointRef::storeOp inside(insideTets, nInside);
-    storeTetOp cutInside(cutInsideTets, nCutInside);
+    tetPointRef::storeOp cutInside(cutInsideTets, nCutInside);
-    dummyTetOp outside;
+    tetPointRef::dummyOp outside;
-
 
 
     // Cut tetA with all inwards pointing faces of tetB. Any tets remaining
@@ -484,7 +74,7 @@ void Foam::tetOverlapVolume::tetTetOverlap
         // Cut and insert subtets into cutInsideTets (either by getting
         // an index from freeSlots or by appending to insideTets) or
         // insert into outsideTets
-        tetSliceWithPlane(tetA, pl0, cutInside, outside);
+        tetA.tet().sliceWithPlane(pl0, cutInside, outside);
     }
 
     if (nCutInside == 0)
@@ -501,7 +91,7 @@ void Foam::tetOverlapVolume::tetTetOverlap
 
         for (label i = 0; i < nCutInside; i++)
         {
-            tetSliceWithPlane(cutInsideTets[i], pl1, inside, outside);
+            cutInsideTets[i].tet().sliceWithPlane(pl1, inside, outside);
         }
 
         if (nInside == 0)
@@ -518,7 +108,7 @@ void Foam::tetOverlapVolume::tetTetOverlap
 
         for (label i = 0; i < nInside; i++)
         {
-            tetSliceWithPlane(insideTets[i], pl2, cutInside, outside);
+            insideTets[i].tet().sliceWithPlane(pl2, cutInside, outside);
         }
 
         if (nCutInside == 0)
@@ -536,28 +126,27 @@ void Foam::tetOverlapVolume::tetTetOverlap
 
         for (label i = 0; i < nCutInside; i++)
         {
-            tetSliceWithPlane(cutInsideTets[i], pl3, inside, outside);
+            cutInsideTets[i].tet().sliceWithPlane(pl3, inside, outside);
         }
     }
 }
 
 
-inline Foam::scalar
+Foam::scalar Foam::tetOverlapVolume::tetTetOverlapVol
-Foam::tetOverlapVolume::tetTetOverlapVol
 (
     const tetPoints& tetA,
     const tetPoints& tetB
-)
+) const
 {
     FixedList<tetPoints, 200> insideTets;
     label nInside = 0;
     FixedList<tetPoints, 200> cutInsideTets;
     label nCutInside = 0;
 
-    storeTetOp inside(insideTets, nInside);
+    tetPointRef::storeOp inside(insideTets, nInside);
-    storeTetOp cutInside(cutInsideTets, nCutInside);
+    tetPointRef::storeOp cutInside(cutInsideTets, nCutInside);
-    sumTetVolOp volInside;
+    tetPointRef::sumVolOp volInside;
-    dummyTetOp outside;
+    tetPointRef::dummyOp outside;
 
     // face0
     plane pl0(tetB[1], tetB[3], tetB[2]);
@@ -605,12 +194,12 @@ Foam::tetOverlapVolume::tetTetOverlapVol
 }
 
 
-inline const Foam::treeBoundBox Foam::tetOverlapVolume::pyrBb
+Foam::treeBoundBox Foam::tetOverlapVolume::pyrBb
 (
     const pointField& points,
     const face& f,
     const point& fc
-)
+) const
 {
     treeBoundBox bb(fc, fc);
     forAll(f, fp)
@@ -750,8 +339,8 @@ Foam::scalar Foam::tetOverlapVolume::cellCellOverlapVolumeMinDecomp
 
 Foam::labelList Foam::tetOverlapVolume::overlappingCells
 (
-    const fvMesh& fromMesh,
+    const polyMesh& fromMesh,
-    const fvMesh& toMesh,
+    const polyMesh& toMesh,
     const label iTo
 ) const
 {
@@ -766,30 +355,4 @@ Foam::labelList Foam::tetOverlapVolume::overlappingCells
 }
 
 
-/*
-
-        forAll(cellsA, i)
-        {
-            label cellAI = cellsA[i];
-            treeBoundBox bbA
-            (
-                pointField(meshA.points(), meshA.cellPoints()[cellAI])
-            );
-
-            scalar v = cellCellOverlapVolumeMinDecomp
-            (
-                meshA,
-                cellAI,
-                bbA,
-                meshB,
-                cellBI,
-                bbB
-            );
-
-            overlapVol += v;
-            nOverlapTests++;
-        }
-
-*/
-
 // ************************************************************************* //

src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -36,21 +36,17 @@ SourceFiles
 #ifndef tetOverlapVolume_H
 #define tetOverlapVolume_H
 
-#include "tetrahedron.H"
+#include "FixedList.H"
-#include "fvMesh.H"
+#include "labelList.H"
-#include "plane.H"
+#include "treeBoundBox.H"
-#include "tetPointRef.H"
-#include "OFstream.H"
-#include "meshTools.H"
-#include "indexedOctree.H"
-#include "treeDataCell.H"
-#include "tetPoints.H"
-#include "tetCell.H"
-#include "EdgeMap.H"
 
 namespace Foam
 {
 
+class primitiveMesh;
+class polyMesh;
+class tetPoints;
+
 /*---------------------------------------------------------------------------*\
                          Class tetOverlapVolume Declaration
 \*---------------------------------------------------------------------------*/
@@ -59,82 +55,6 @@ class tetOverlapVolume
 {
     // Private member functions
 
-        //- Plane intersection
-        inline point planeIntersection
-        (
-            const FixedList<scalar, 4>& d,
-            const tetPoints& t,
-            const label negI,
-            const label posI
-        );
-
-
-        //- Decompose prism
-        template <class TetOp> inline void decomposePrism
-        (
-            const FixedList<point, 6>& points,
-            TetOp& op
-        );
-
-
-        //- Helping cľasses
-        class dummyTetOp
-        {
-         public:
-
-            inline void operator()(const tetPoints&){}
-        };
-
-
-        class sumTetVolOp
-        {
-         public:
-            scalar vol_;
-
-            inline sumTetVolOp()
-            :
-                vol_(0.0)
-            {}
-
-            inline void operator()(const tetPoints& tet)
-            {
-                vol_ += tet.tet().mag();
-            }
-        };
-
-
-        class storeTetOp
-        {
-            FixedList<tetPoints, 200>& tets_;
-            label& nTets_;
-
-         public:
-
-            inline storeTetOp(FixedList<tetPoints, 200>& tets, label& nTets)
-            :
-                tets_(tets),
-                nTets_(nTets)
-            {}
-
-            inline void operator()(const tetPoints& tet)
-            {
-                tets_[nTets_++] = tet;
-            }
-        };
-
-
-        //- Slice. Split tet into subtets above and below plane
-        template <class AboveTetOp, class BelowTetOp>
-        inline void tetSliceWithPlane
-        (
-            const tetPoints& tet,
-            const plane& pl,
-
-            AboveTetOp& aboveOp,
-            BelowTetOp& belowOp
-        );
-
-
         //- Tet overlap
         void tetTetOverlap
         (
@@ -144,31 +64,28 @@ class tetOverlapVolume
             label& nInside,
             FixedList<tetPoints, 200>& outsideTets,
             label& nOutside
-        );
+        ) const;
-
 
         //- Tet Overlap Vol
-        inline scalar tetTetOverlapVol
+        scalar tetTetOverlapVol
         (
             const tetPoints& tetA,
             const tetPoints& tetB
-        );
+        ) const;
-
 
         //- Return a const treeBoundBox
-        inline const treeBoundBox pyrBb
+        treeBoundBox pyrBb
         (
             const pointField& points,
             const face& f,
             const point& fc
-        );
+        ) const;
-
 
 public:
 
 
     //- Runtime type information
-    TypeName("tetOverlapVolume");
+    ClassName("tetOverlapVolume");
 
 
     // Constructors
@@ -177,18 +94,14 @@ public:
         tetOverlapVolume();
 
 
-    //- Destructor
-     virtual ~tetOverlapVolume();
-
-
     // Public members
 
         //- Return a list of cells in meshA which overlaps with cellBI in
         // meshB
         labelList overlappingCells
         (
-            const fvMesh& meshA,
+            const polyMesh& meshA,
-            const fvMesh& meshB,
+            const polyMesh& meshB,
             const label cellBI
         ) const;
 

src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetrahedron.C

@@ -1,341 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2012 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 3 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, see <http://www.gnu.org/licenses/>.
-
-Description
-    Calculation of shape function product for a tetrahedron
-
-\*---------------------------------------------------------------------------*/
-
-#include "tetrahedron.H"
-#include "triPointRef.H"
-#include "scalarField.H"
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-// (Probably very inefficient) minimum containment sphere calculation.
-// From http://www.imr.sandia.gov/papers/imr11/shewchuk2.pdf:
-// Sphere ctr is smallest one of
-// - tet circumcentre
-// - triangle circumcentre
-// - edge mids
-template<class Point, class PointRef>
-Foam::pointHit Foam::tetrahedron<Point, PointRef>::containmentSphere
-(
-    const scalar tol
-) const
-{
-    const scalar fac = 1 + tol;
-
-    // Halve order of tolerance for comparisons of sqr.
-    const scalar facSqr = Foam::sqrt(fac);
-
-
-    // 1. Circumcentre itself.
-
-    pointHit pHit(circumCentre());
-    pHit.setHit();
-    scalar minRadiusSqr = magSqr(pHit.rawPoint() - a_);
-
-
-    // 2. Try circumcentre of tet triangles. Create circumcircle for triFace and
-    // check if 4th point is inside.
-
-    {
-        point ctr = triPointRef(a_, b_, c_).circumCentre();
-        scalar radiusSqr = magSqr(ctr - a_);
-
-        if
-        (
-            radiusSqr < minRadiusSqr
-         && Foam::magSqr(d_-ctr) <= facSqr*radiusSqr
-        )
-        {
-            pHit.setMiss(false);
-            pHit.setPoint(ctr);
-            minRadiusSqr = radiusSqr;
-        }
-    }
-    {
-        point ctr = triPointRef(a_, b_, d_).circumCentre();
-        scalar radiusSqr = magSqr(ctr - a_);
-
-        if
-        (
-            radiusSqr < minRadiusSqr
-         && Foam::magSqr(c_-ctr) <= facSqr*radiusSqr
-        )
-        {
-            pHit.setMiss(false);
-            pHit.setPoint(ctr);
-            minRadiusSqr = radiusSqr;
-        }
-    }
-    {
-        point ctr = triPointRef(a_, c_, d_).circumCentre();
-        scalar radiusSqr = magSqr(ctr - a_);
-
-        if
-        (
-            radiusSqr < minRadiusSqr
-         && Foam::magSqr(b_-ctr) <= facSqr*radiusSqr
-        )
-        {
-            pHit.setMiss(false);
-            pHit.setPoint(ctr);
-            minRadiusSqr = radiusSqr;
-        }
-    }
-    {
-        point ctr = triPointRef(b_, c_, d_).circumCentre();
-        scalar radiusSqr = magSqr(ctr - b_);
-
-        if
-        (
-            radiusSqr < minRadiusSqr
-         && Foam::magSqr(a_-ctr) <= facSqr*radiusSqr
-        )
-        {
-            pHit.setMiss(false);
-            pHit.setPoint(ctr);
-            minRadiusSqr = radiusSqr;
-        }
-    }
-
-
-    // 3. Try midpoints of edges
-
-    // mid of edge A-B
-    {
-        point ctr = 0.5*(a_ + b_);
-        scalar radiusSqr = magSqr(a_ - ctr);
-        scalar testRadSrq = facSqr*radiusSqr;
-
-        if
-        (
-            radiusSqr < minRadiusSqr
-         && magSqr(c_-ctr) <= testRadSrq
-         && magSqr(d_-ctr) <= testRadSrq)
-        {
-            pHit.setMiss(false);
-            pHit.setPoint(ctr);
-            minRadiusSqr = radiusSqr;
-        }
-    }
-
-    // mid of edge A-C
-    {
-        point ctr = 0.5*(a_ + c_);
-        scalar radiusSqr = magSqr(a_ - ctr);
-        scalar testRadSrq = facSqr*radiusSqr;
-
-        if
-        (
-            radiusSqr < minRadiusSqr
-         && magSqr(b_-ctr) <= testRadSrq
-         && magSqr(d_-ctr) <= testRadSrq
-        )
-        {
-            pHit.setMiss(false);
-            pHit.setPoint(ctr);
-            minRadiusSqr = radiusSqr;
-        }
-    }
-
-    // mid of edge A-D
-    {
-        point ctr = 0.5*(a_ + d_);
-        scalar radiusSqr = magSqr(a_ - ctr);
-        scalar testRadSrq = facSqr*radiusSqr;
-
-        if
-        (
-            radiusSqr < minRadiusSqr
-         && magSqr(b_-ctr) <= testRadSrq
-         && magSqr(c_-ctr) <= testRadSrq
-        )
-        {
-            pHit.setMiss(false);
-            pHit.setPoint(ctr);
-            minRadiusSqr = radiusSqr;
-        }
-    }
-
-    // mid of edge B-C
-    {
-        point ctr = 0.5*(b_ + c_);
-        scalar radiusSqr = magSqr(b_ - ctr);
-        scalar testRadSrq = facSqr*radiusSqr;
-
-        if
-        (
-            radiusSqr < minRadiusSqr
-         && magSqr(a_-ctr) <= testRadSrq
-         && magSqr(d_-ctr) <= testRadSrq
-        )
-        {
-            pHit.setMiss(false);
-            pHit.setPoint(ctr);
-            minRadiusSqr = radiusSqr;
-        }
-    }
-
-    // mid of edge B-D
-    {
-        point ctr = 0.5*(b_ + d_);
-        scalar radiusSqr = magSqr(b_ - ctr);
-        scalar testRadSrq = facSqr*radiusSqr;
-
-        if
-        (
-            radiusSqr < minRadiusSqr
-         && magSqr(a_-ctr) <= testRadSrq
-         && magSqr(c_-ctr) <= testRadSrq)
-        {
-            pHit.setMiss(false);
-            pHit.setPoint(ctr);
-            minRadiusSqr = radiusSqr;
-        }
-    }
-
-    // mid of edge C-D
-    {
-        point ctr = 0.5*(c_ + d_);
-        scalar radiusSqr = magSqr(c_ - ctr);
-        scalar testRadSrq = facSqr*radiusSqr;
-
-        if
-        (
-            radiusSqr < minRadiusSqr
-         && magSqr(a_-ctr) <= testRadSrq
-         && magSqr(b_-ctr) <= testRadSrq
-        )
-        {
-            pHit.setMiss(false);
-            pHit.setPoint(ctr);
-            minRadiusSqr = radiusSqr;
-        }
-    }
-
-
-    pHit.setDistance(sqrt(minRadiusSqr));
-
-    return pHit;
-}
-
-
-template<class Point, class PointRef>
-void Foam::tetrahedron<Point, PointRef>::gradNiSquared
-(
-    scalarField& buffer
-) const
-{
-    // Change of sign between face area vector and gradient
-    // does not matter because of square
-
-    // Warning: Added a mag to produce positive coefficients even if
-    // the tetrahedron is twisted inside out.  This is pretty
-    // dangerous, but essential for mesh motion.
-    scalar magVol = Foam::mag(mag());
-
-    buffer[0] = (1.0/9.0)*magSqr(Sa())/magVol;
-    buffer[1] = (1.0/9.0)*magSqr(Sb())/magVol;
-    buffer[2] = (1.0/9.0)*magSqr(Sc())/magVol;
-    buffer[3] = (1.0/9.0)*magSqr(Sd())/magVol;
-}
-
-
-template<class Point, class PointRef>
-void Foam::tetrahedron<Point, PointRef>::gradNiDotGradNj
-(
-    scalarField& buffer
-) const
-{
-    // Warning. Ordering of edges needs to be the same for a tetrahedron
-    // class, a tetrahedron cell shape model and a tetCell
-
-    // Warning: Added a mag to produce positive coefficients even if
-    // the tetrahedron is twisted inside out.  This is pretty
-    // dangerous, but essential for mesh motion.
-
-    // Double change of sign between face area vector and gradient
-
-    scalar magVol = Foam::mag(mag());
-    vector sa = Sa();
-    vector sb = Sb();
-    vector sc = Sc();
-    vector sd = Sd();
-
-    buffer[0] = (1.0/9.0)*(sa & sb)/magVol;
-    buffer[1] = (1.0/9.0)*(sa & sc)/magVol;
-    buffer[2] = (1.0/9.0)*(sa & sd)/magVol;
-    buffer[3] = (1.0/9.0)*(sd & sb)/magVol;
-    buffer[4] = (1.0/9.0)*(sb & sc)/magVol;
-    buffer[5] = (1.0/9.0)*(sd & sc)/magVol;
-}
-
-
-template<class Point, class PointRef>
-void Foam::tetrahedron<Point, PointRef>::gradNiGradNi
-(
-    tensorField& buffer
-) const
-{
-    // Change of sign between face area vector and gradient
-    // does not matter because of square
-
-    scalar magVol = Foam::mag(mag());
-
-    buffer[0] = (1.0/9.0)*sqr(Sa())/magVol;
-    buffer[1] = (1.0/9.0)*sqr(Sb())/magVol;
-    buffer[2] = (1.0/9.0)*sqr(Sc())/magVol;
-    buffer[3] = (1.0/9.0)*sqr(Sd())/magVol;
-}
-
-
-template<class Point, class PointRef>
-void Foam::tetrahedron<Point, PointRef>::gradNiGradNj
-(
-    tensorField& buffer
-) const
-{
-    // Warning. Ordering of edges needs to be the same for a tetrahedron
-    // class, a tetrahedron cell shape model and a tetCell
-
-    // Double change of sign between face area vector and gradient
-
-    scalar magVol = Foam::mag(mag());
-    vector sa = Sa();
-    vector sb = Sb();
-    vector sc = Sc();
-    vector sd = Sd();
-
-    buffer[0] = (1.0/9.0)*(sa * sb)/magVol;
-    buffer[1] = (1.0/9.0)*(sa * sc)/magVol;
-    buffer[2] = (1.0/9.0)*(sa * sd)/magVol;
-    buffer[3] = (1.0/9.0)*(sd * sb)/magVol;
-    buffer[4] = (1.0/9.0)*(sb * sc)/magVol;
-    buffer[5] = (1.0/9.0)*(sd * sc)/magVol;
-}
-
-
-// ************************************************************************* //

src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetrahedron.H

@@ -1,315 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2012 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 3 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, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::tetrahedron
-
-Description
-    A tetrahedron primitive.
-
-    Ordering of edges needs to be the same for a tetrahedron
-    class, a tetrahedron cell shape model and a tetCell.
-
-SourceFiles
-    tetrahedronI.H
-    tetrahedron.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef tetrahedron_H
-#define tetrahedron_H
-
-#include "point.H"
-#include "primitiveFieldsFwd.H"
-#include "pointHit.H"
-#include "cachedRandom.H"
-#include "Random.H"
-#include "FixedList.H"
-#include "UList.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-class Istream;
-class Ostream;
-class tetPoints;
-class plane;
-
-// Forward declaration of friend functions and operators
-
-template<class Point, class PointRef> class tetrahedron;
-
-template<class Point, class PointRef>
-inline Istream& operator>>
-(
-    Istream&,
-    tetrahedron<Point, PointRef>&
-);
-
-template<class Point, class PointRef>
-inline Ostream& operator<<
-(
-    Ostream&,
-    const tetrahedron<Point, PointRef>&
-);
-
-
-/*---------------------------------------------------------------------------*\
-                           class tetrahedron Declaration
-\*---------------------------------------------------------------------------*/
-
-template<class Point, class PointRef>
-class tetrahedron
-{
-public:
-
-        // Classes for use in sliceWithPlane. What to do with decomposition
-        // of tet.
-
-            //- Dummy
-            class dummyOp
-            {
-            public:
-                inline void operator()(const tetPoints&);
-            };
-
-            //- Sum resulting volumes
-            class sumVolOp
-            {
-            public:
-                scalar vol_;
-
-                inline sumVolOp();
-
-                inline void operator()(const tetPoints&);
-            };
-
-            //- Store resulting tets
-            class storeOp
-            {
-                FixedList<tetPoints, 200>& tets_;
-                label& nTets_;
-
-            public:
-                inline storeOp(FixedList<tetPoints, 200>&, label&);
-
-                inline void operator()(const tetPoints&);
-            };
-
-private:
-
-    // Private data
-
-        PointRef a_, b_, c_, d_;
-
-        inline static point planeIntersection
-        (
-            const FixedList<scalar, 4>&,
-            const tetPoints&,
-            const label,
-            const label
-        );
-
-        template<class TetOp>
-        inline static void decomposePrism
-        (
-            const FixedList<point, 6>& points,
-            TetOp& op
-        );
-
-        template<class AboveTetOp, class BelowTetOp>
-        inline static void tetSliceWithPlane
-        (
-            const plane& pl,
-            const tetPoints& tet,
-            AboveTetOp& aboveOp,
-            BelowTetOp& belowOp
-        );
-
-
-public:
-
-    // Member constants
-
-        enum
-        {
-            nVertices = 4,  // Number of vertices in tetrahedron
-            nEdges = 6      // Number of edges in tetrahedron
-        };
-
-
-    // Constructors
-
-        //- Construct from points
-        inline tetrahedron
-        (
-            const Point& a,
-            const Point& b,
-            const Point& c,
-            const Point& d
-        );
-
-        //- Construct from four points in the list of points
-        inline tetrahedron
-        (
-            const UList<Point>&,
-            const FixedList<label, 4>& indices
-        );
-
-        //- Construct from Istream
-        inline tetrahedron(Istream&);
-
-
-    // Member Functions
-
-        // Access
-
-            //- Return vertices
-            inline const Point& a() const;
-
-            inline const Point& b() const;
-
-            inline const Point& c() const;
-
-            inline const Point& d() const;
-
-
-        // Properties
-
-            //- Return face normal
-            inline vector Sa() const;
-
-            inline vector Sb() const;
-
-            inline vector Sc() const;
-
-            inline vector Sd() const;
-
-            //- Return centre (centroid)
-            inline Point centre() const;
-
-            //- Return volume
-            inline scalar mag() const;
-
-            //- Return circum-centre
-            inline Point circumCentre() const;
-
-            //- Return circum-radius
-            inline scalar circumRadius() const;
-
-            //- Return quality: Ratio of tetrahedron and circum-sphere
-            //  volume, scaled so that a regular tetrahedron has a
-            //  quality of 1
-            inline scalar quality() const;
-
-            //- Return a random point in the tetrahedron from a
-            //  uniform distribution
-            inline Point randomPoint(Random& rndGen) const;
-
-            //- Return a random point in the tetrahedron from a
-            //  uniform distribution
-            inline Point randomPoint(cachedRandom& rndGen) const;
-
-            //- Calculate the barycentric coordinates of the given
-            //  point, in the same order as a, b, c, d.  Returns the
-            //  determinant of the solution.
-            inline scalar barycentric
-            (
-                const point& pt,
-                List<scalar>& bary
-            ) const;
-
-            //- Return nearest point to p on tetrahedron
-            inline pointHit nearestPoint(const point& p) const;
-
-            //- Return true if point is inside tetrahedron
-            inline bool inside(const point& pt) const;
-
-            //- Decompose tet into tets above and below plane
-            template<class AboveTetOp, class BelowTetOp>
-            inline void sliceWithPlane
-            (
-                const plane& pl,
-                AboveTetOp& aboveOp,
-                BelowTetOp& belowOp
-            ) const;
-
-
-            //- Return (min)containment sphere, i.e. the smallest sphere with
-            //  all points inside. Returns pointHit with:
-            //  - hit         : if sphere is equal to circumsphere
-            //                  (biggest sphere)
-            //  - point       : centre of sphere
-            //  - distance    : radius of sphere
-            //  - eligiblemiss: false
-            // Tol (small compared to 1, e.g. 1E-9) is used to determine
-            // whether point is inside: mag(pt - ctr) < (1+tol)*radius.
-            pointHit containmentSphere(const scalar tol) const;
-
-            //- Fill buffer with shape function products
-            void gradNiSquared(scalarField& buffer) const;
-
-            void gradNiDotGradNj(scalarField& buffer) const;
-
-            void gradNiGradNi(tensorField& buffer) const;
-
-            void gradNiGradNj(tensorField& buffer) const;
-
-
-    // IOstream operators
-
-        friend Istream& operator>> <Point, PointRef>
-        (
-            Istream&,
-            tetrahedron&
-        );
-
-        friend Ostream& operator<< <Point, PointRef>
-        (
-            Ostream&,
-            const tetrahedron&
-        );
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#include "tetrahedronI.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#ifdef NoRepository
-#   include "tetrahedron.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

tutorials/multiphase/interFoam/ras/damBreakPorousBaffle/system/fvSolution

@@ -68,11 +68,11 @@ PIMPLE
     cAlpha          2;
 }
 
-relaxation
+relaxationFactors
 {
-    U           1;
+    "U.*"           1;
-    k           1;
+    "k.*"           1;
-    epsilon     1;
+    "epsilon.*"     1;
 }
 
 // ************************************************************************* //