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

src/OpenFOAM/Make/files

@@ -45,7 +45,9 @@ primitives/Tensor/lists/symmTensorList.C
 primitives/Tensor/lists/tensorList.C
 
 primitives/Vector/complexVector/complexVector.C
+#if !defined(WM_SP)
 primitives/Vector/floatVector/floatVector.C
+#endif
 primitives/Vector/labelVector/labelVector.C
 primitives/Vector/vector/vector.C
 primitives/Vector/lists/vectorList.C

src/OpenFOAM/containers/LinkedLists/linkTypes/DLListBase/DLListBase.C

@@ -94,23 +94,23 @@ bool Foam::DLListBase::swapUp(DLListBase::link* a)
         if (ap == first_)
         {
             first_ = a;
+            ap->prev_ = a;
+        }
+        else
+        {
+            ap->prev_->next_ = a;
         }
 
         if (a == last_)
         {
             last_ = ap;
+            a->next_ = ap;
         }
-
-        if (a->next_)
+        else
         {
             a->next_->prev_ = ap;
         }
 
-        if (ap->prev_)
-        {
-            ap->prev_->next_ = a;
-        }
-
         a->prev_ = ap->prev_;
         ap->prev_ = a;
 
@@ -135,19 +135,19 @@ bool Foam::DLListBase::swapDown(DLListBase::link* a)
         if (a == first_)
         {
             first_ = an;
+            a->prev_ = an;
+        }
+        else
+        {
+            a->prev_->next_ = an;
         }
 
         if (an == last_)
         {
             last_ = a;
+            an->next_ = a;
         }
-
-        if (a->prev_)
-        {
-            a->prev_->next_ = an;
-        }
-
-        if (an->next_)
+        else
         {
             an->next_->prev_ = a;
         }

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/db/Time/Time.C

@@ -203,10 +203,15 @@ void Foam::Time::setControls()
         )
     );
 
-    if (timeDict.readIfPresent("deltaT", deltaT_))
+    // Read and set the deltaT only if time-step adjustment is active
+    // otherwise use the deltaT from the controlDict
+    if (controlDict_.lookupOrDefault<Switch>("adjustTimeStep", false))
     {
-        deltaTSave_ = deltaT_;
-        deltaT0_ = deltaT_;
+        if (timeDict.readIfPresent("deltaT", deltaT_))
+        {
+            deltaTSave_ = deltaT_;
+            deltaT0_ = deltaT_;
+        }
     }
 
     timeDict.readIfPresent("deltaT0", deltaT0_);
@@ -984,6 +989,17 @@ Foam::Time& Foam::Time::operator++()
             << " to " << precision_
             << " to distinguish between timeNames at time " << value()
             << endl;
+
+        if (precision_ == 100 && precision_ != oldPrecision)
+        {
+            // Reached limit.
+            WarningIn("Time::operator++()")
+                << "Current time name " << dimensionedScalar::name()
+                << " is the old as the previous one " << oldTimeName
+                << endl
+                << "    This might result in overwriting old results."
+                << endl;
+        }
     }
 
 

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
-                << "    nProcsSimpleSum   : " << Pstream::nProcsSimpleSum << nl
-                << "    commsType         : "
-                << Pstream::commsTypeNames[Pstream::defaultCommsType]
+                << "    floatTransfer      : " << Pstream::floatTransfer << nl
+                << "    nProcsSimpleSum    : " << Pstream::nProcsSimpleSum << nl
+                << "    commsType          : "
+                << 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
-    (
-        Pstream::defaultCommsType == Pstream::blocking
-     || Pstream::defaultCommsType == Pstream::nonBlocking
-    )
+    if (Pstream::defaultCommsType == Pstream::blocking)
     {
-        // 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
-    (
-        commsType,
-        procInterface_.interfaceInternalField(psiInternal)()
-    );
+    procInterface_.interfaceInternalField(psiInternal, scalarSendBuf_);
+
+    if (commsType == Pstream::nonBlocking && !Pstream::floatTransfer)
+    {
+        // 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
-    (
-        procInterface_.compressedReceive<scalar>(commsType, coeffs.size())
-    );
-    transformCoupleField(pnf, cmpt);
+    if (updatedMatrix())
+    {
+        return;
+    }
 
     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);
-            vector n1 = -findFaceMaxRadius(half1Ctrs);
-            n0 /= mag(n0) + VSMALL;
-            n1 /= mag(n1) + VSMALL;
-
-            if (debug)
+            case ROTATIONAL:
             {
-                Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
-                    << " patch:" << name()
-                    << " Specified rotation :"
-                    << " n0:" << n0 << " n1:" << n1 << endl;
-            }
+                vector n0 = findFaceMaxRadius(half0Ctrs);
+                vector n1 = -findFaceMaxRadius(half1Ctrs);
+                n0 /= mag(n0) + VSMALL;
+                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)
-                const tensor revT(rotationTensor(n0, -n1));
+                // 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]
-                    );
-                    anchors0[faceI] = Foam::transform
-                    (
-                        revT,
-                        anchors0[faceI]
-                    );
+                    half0Ctrs[faceI] =
+                        Foam::transform
+                        (
+                            revT,
+                            half0Ctrs[faceI] - rotationCentre_
+                        )
+                      + rotationCentre_;
+                    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());
-                const point ctr1(sum(pp1.localPoints())/pp1.nPoints());
+                break;
+            }
+            case TRANSLATIONAL:
+            {
+                // Transform 0 points.
 
                 if (debug)
                 {
                     Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
                         << " patch:" << name()
-                        << " Detected translation :"
-                        << " n0:" << n0 << " n1:" << n1
-                        << " ctr0:" << ctr0 << " ctr1:" << ctr1 << endl;
+                        << "Specified translation : " << separationVector_
+                        << endl;
                 }
 
-                half0Ctrs += ctr1 - ctr0;
-                anchors0 += ctr1 - ctr0;
+                // 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 :"
+                            << " 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/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/Vector/floatVector/floatVector.H

@@ -47,11 +47,11 @@ namespace Foam
 
 typedef Vector<float> floatVector;
 
-
 //- Data associated with floatVector type are contiguous
+#if !defined(WM_SP)
 template<>
 inline bool contiguous<floatVector>() {return true;}
-
+#endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

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/conversion/ensight/file/ensightFile.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) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "ensightFile.H"
+#include <sstream>
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -31,6 +32,53 @@ bool Foam::ensightFile::allowUndef_ = false;
 
 Foam::scalar Foam::ensightFile::undefValue_ = Foam::floatScalarVGREAT;
 
+// default is width 8
+Foam::string Foam::ensightFile::mask_ = "********";
+
+Foam::string Foam::ensightFile::dirFmt_ = "%08d";
+
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+Foam::string Foam::ensightFile::mask()
+{
+    return mask_;
+}
+
+
+Foam::string Foam::ensightFile::subDir(const label n)
+{
+    char buf[32];
+
+    sprintf(buf, dirFmt_.c_str(), n);
+    return buf;
+}
+
+
+void Foam::ensightFile::subDirWidth(const label n)
+{
+    // enforce max limit to avoid buffer overflow in subDir()
+    if (n < 1 || n > 31)
+    {
+        return;
+    }
+
+    // appropriate printf format
+    std::ostringstream oss;
+    oss << "%0" << n << "d";
+    dirFmt_ = oss.str();
+
+    // set mask accordingly
+    mask_.resize(n, '*');
+}
+
+
+Foam::label Foam::ensightFile::subDirWidth()
+{
+    return mask_.size();
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::ensightFile::ensightFile
@@ -247,22 +295,4 @@ Foam::Ostream& Foam::ensightFile::writeBinaryHeader()
 }
 
 
-// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
-
-Foam::string Foam::ensightFile::mask()
-{
-    char buf[16] = "********";
-    return buf;
-}
-
-
-Foam::string Foam::ensightFile::subDir(const label n)
-{
-    char buf[16];
-
-    sprintf(buf, "%08d", n);
-    return buf;
-}
-
-
 // ************************************************************************* //

src/conversion/ensight/file/ensightFile.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) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -57,6 +57,12 @@ class ensightFile
         //- value to represent undef in results
         static scalar undefValue_;
 
+        //- The '*' mask appropriate for subDir
+        static string mask_;
+
+        //- The printf format for zero-padded subdirectory numbers
+        static string dirFmt_;
+
 
     // Private Member Functions
 
@@ -88,12 +94,19 @@ public:
         //- Return setting for whether 'undef' values are allowed in results
         static bool allowUndef();
 
-        //- '*' mask appropriate for subDir
+        //- The '*' mask appropriate for subDir
         static string mask();
 
-        //- consistent zero-padded numbers for subdirectories
+        //- Consistent zero-padded numbers for subdirectories
         static string subDir(const label);
 
+        //- Set width of subDir and mask. Default width is 8 digits.
+        //  Max width is 31 digits.
+        static void subDirWidth(const label);
+
+        //- Return current width of subDir and mask.
+        static label subDirWidth();
+
 
     // Edit
 

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/edgeFormats/emesh/EMESHedgeFormat.C

@@ -176,7 +176,7 @@ void Foam::fileFormats::EMESHedgeFormat::write
         << "    version     " << os.version() << ";\n"
         << "    format      " << os.format() << ";\n"
         << "    class       " << "featureEdgeMesh" << ";\n"
-        << "    note        " << "written " + clock::dateTime() << nl
+        << "    note        " << "written " + clock::dateTime() << ";\n"
         << "    object      " << filename.name() << ";\n"
         << "}" << nl;
 

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/fieldSources/basicSource/basicSource/IObasicSourceList.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) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -45,7 +45,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                        Class IObasicSourceList Declaration
+                      Class IObasicSourceList Declaration
 \*---------------------------------------------------------------------------*/
 
 class IObasicSourceList
@@ -80,6 +80,8 @@ public:
         {}
 
 
+    // Member Functions
+
         //- Read dictionary
         virtual bool read();
 };

src/fieldSources/basicSource/basicSource/basicSource.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) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,6 +26,7 @@ License
 #include "basicSource.H"
 #include "fvMesh.H"
 #include "fvMatrices.H"
+#include "volFields.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -161,7 +162,7 @@ void Foam::basicSource::setCellSet()
         }
         case smMapRegion:
         {
-            if(active_)
+            if (active_)
             {
                 Info<< indent << "- selecting inter region mapping" << endl;
                 const fvMesh& secondaryMesh =
@@ -170,7 +171,6 @@ void Foam::basicSource::setCellSet()
                 const boundBox secondaryBB = secondaryMesh.bounds();
                 if (secondaryBB.overlaps(primaryBB))
                 {
-
                     // Dummy patches
                     wordList cuttingPatches;
                     HashTable<word> patchMap;
@@ -218,7 +218,7 @@ void Foam::basicSource::setCellSet()
     }
 
     // Set volume information
-    if(selectionMode_ != smMapRegion)
+    if (selectionMode_ != smMapRegion)
     {
         V_ = 0.0;
         forAll(cells_, i)
@@ -303,6 +303,7 @@ Foam::autoPtr<Foam::basicSource> Foam::basicSource::New
     return autoPtr<basicSource>(cstrIter()(name, modelType, coeffs, mesh));
 }
 
+
 Foam::basicSource::~basicSource()
 {
     if (!secondaryToPrimaryInterpPtr_.empty())
@@ -364,6 +365,36 @@ void Foam::basicSource::checkApplied() const
 }
 
 
+void Foam::basicSource::correct(volScalarField& fld)
+{
+    // do nothing
+}
+
+
+void Foam::basicSource::correct(volVectorField& fld)
+{
+    // do nothing
+}
+
+
+void Foam::basicSource::correct(volSphericalTensorField& fld)
+{
+    // do nothing
+}
+
+
+void Foam::basicSource::correct(volSymmTensorField& fld)
+{
+    // do nothing
+}
+
+
+void Foam::basicSource::correct(volTensorField& fld)
+{
+    // do nothing
+}
+
+
 void Foam::basicSource::addSup(fvMatrix<scalar>& eqn, const label fieldI)
 {
     // do nothing

src/fieldSources/basicSource/basicSource/basicSource.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) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -47,6 +47,7 @@ SourceFiles
 #define basicSource_H
 
 #include "fvMatricesFwd.H"
+#include "volFieldsFwd.H"
 #include "cellSet.H"
 #include "autoPtr.H"
 #include "meshToMesh.H"
@@ -60,7 +61,6 @@ namespace Foam
 
 class fvMesh;
 
-
 /*---------------------------------------------------------------------------*\
                          Class basicSource Declaration
 \*---------------------------------------------------------------------------*/
@@ -330,6 +330,24 @@ public:
 
         // Evaluation
 
+            // Correct
+
+                //- Scalar
+                virtual void correct(volScalarField& fld);
+
+                //- Vector
+                virtual void correct(volVectorField& fld);
+
+                //- Spherical tensor
+                virtual void correct(volSphericalTensorField& fld);
+
+                //- Symmetric tensor
+                virtual void correct(volSymmTensorField& fld);
+
+                //- Tensor
+                virtual void correct(volTensorField& fld);
+
+
             // Add explicit and implicit contributions
 
                 //- Scalar

src/fieldSources/basicSource/basicSource/basicSourceI.H

@@ -135,4 +135,5 @@ inline const Foam::word Foam::basicSource::mapRegionName() const
     return mapRegionName_;
 }
 
+
 // ************************************************************************* //

src/fieldSources/basicSource/basicSource/basicSourceList.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) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -95,6 +95,11 @@ public:
 
     // Member Functions
 
+        //- Correct
+        template<class Type>
+        void correct(GeometricField<Type, fvPatchField, volMesh>& fld);
+
+
         // Sources
 
             //- Return source for equation

src/fieldSources/basicSource/basicSource/basicSourceListTemplates.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) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,6 +25,39 @@ License
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class Type>
+void Foam::basicSourceList::correct
+(
+    GeometricField<Type, fvPatchField, volMesh>& fld
+)
+{
+    const word& fieldName = fld.name();
+
+    forAll(*this, i)
+    {
+        basicSource& source = this->operator[](i);
+
+        label fieldI = source.applyToField(fieldName);
+
+        if (fieldI != -1)
+        {
+            source.setApplied(fieldI);
+
+            if (source.isActive())
+            {
+                if (debug)
+                {
+                    Info<< "Correcting source " << source.name()
+                        << " for field " << fieldName << endl;
+                }
+
+                source.correct(fld);
+            }
+        }
+    }
+}
+
+
 template<class Type>
 Foam::tmp<Foam::fvMatrix<Type> > Foam::basicSourceList::operator()
 (

src/fieldSources/basicSource/pressureGradientExplicitSource/pressureGradientExplicitSource.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) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -69,53 +69,6 @@ void Foam::pressureGradientExplicitSource::writeGradP() const
 }
 
 
-void Foam::pressureGradientExplicitSource::update(fvMatrix<vector>& eqn)
-{
-    volVectorField& U = const_cast<volVectorField&>(eqn.psi());
-
-    const volScalarField& rAU =
-        mesh_.lookupObject<volScalarField>("(1|A(" + U.name() + "))");
-
-    // Integrate flow variables over cell set
-    scalar magUbarAve = 0.0;
-    scalar rAUave = 0.0;
-    const scalarField& cv = mesh_.V();
-    forAll(cells_, i)
-    {
-        label cellI = cells_[i];
-        scalar volCell = cv[cellI];
-        magUbarAve += (flowDir_ & U[cellI])*volCell;
-        rAUave += rAU[cellI]*volCell;
-    }
-
-    // Collect across all processors
-    reduce(magUbarAve, sumOp<scalar>());
-
-    // Volume averages
-    magUbarAve /= V_;
-    rAUave /= V_;
-
-    // Calculate the pressure gradient increment needed to adjust the average
-    // flow-rate to the desired value
-    scalar gradPplus = (mag(Ubar_) - magUbarAve)/rAUave;
-
-    // Apply correction to velocity field
-    forAll(cells_, i)
-    {
-        label cellI = cells_[i];
-        U[cellI] += flowDir_*rAU[cellI]*gradPplus;
-    }
-
-    // Update pressure gradient
-    gradP_.value() += gradPplus;
-
-    Info<< "Uncorrected Ubar = " << magUbarAve << tab
-        << "Pressure gradient = " << gradP_.value() << endl;
-
-    writeGradP();
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::pressureGradientExplicitSource::pressureGradientExplicitSource
@@ -123,7 +76,7 @@ Foam::pressureGradientExplicitSource::pressureGradientExplicitSource
     const word& sourceName,
     const word& modelType,
     const dictionary& dict,
-    const fvMesh& mesh    
+    const fvMesh& mesh
 )
 :
     basicSource(sourceName, modelType, dict, mesh),
@@ -133,6 +86,23 @@ Foam::pressureGradientExplicitSource::pressureGradientExplicitSource
     flowDir_(Ubar_/mag(Ubar_))
 {
     coeffs_.lookup("fieldNames") >> fieldNames_;
+
+    if (fieldNames_.size() != 1)
+    {
+        FatalErrorIn
+        (
+            "Foam::pressureGradientExplicitSource::"
+            "pressureGradientExplicitSource"
+            "("
+                "onst word&, "
+                "const word&, "
+                "const dictionary&, "
+                "const fvMesh&"
+            ")"
+        )   << "Source can only be applied to a single field.  Current "
+            << "settings are:" << fieldNames_ << exit(FatalError);
+    }
+
     applied_.setSize(fieldNames_.size(), false);
 
     // Read the initial pressure gradient from file if it exists
@@ -160,8 +130,6 @@ void Foam::pressureGradientExplicitSource::addSup
     const label fieldI
 )
 {
-    update(eqn);
-
     DimensionedField<vector, volMesh> Su
     (
         IOobject
@@ -178,7 +146,53 @@ void Foam::pressureGradientExplicitSource::addSup
 
     UIndirectList<vector>(Su, cells_) = flowDir_*gradP_.value();
 
-    eqn -= Su;
+    eqn += Su;
+}
+
+
+void Foam::pressureGradientExplicitSource::correct(volVectorField& U)
+{
+    const volScalarField& rAU =
+        mesh_.lookupObject<volScalarField>("(1|A(" + U.name() + "))");
+
+    // Integrate flow variables over cell set
+    scalar magUbarAve = 0.0;
+    scalar rAUave = 0.0;
+    const scalarField& cv = mesh_.V();
+    forAll(cells_, i)
+    {
+        label cellI = cells_[i];
+        scalar volCell = cv[cellI];
+        magUbarAve += (flowDir_ & U[cellI])*volCell;
+        rAUave += rAU[cellI]*volCell;
+    }
+
+    // Collect across all processors
+    reduce(magUbarAve, sumOp<scalar>());
+    reduce(rAUave, sumOp<scalar>());
+
+    // Volume averages
+    magUbarAve /= V_;
+    rAUave /= V_;
+
+    // Calculate the pressure gradient increment needed to adjust the average
+    // flow-rate to the desired value
+    scalar gradPplus = (mag(Ubar_) - magUbarAve)/rAUave;
+
+    // Apply correction to velocity field
+    forAll(cells_, i)
+    {
+        label cellI = cells_[i];
+        U[cellI] += flowDir_*rAU[cellI]*gradPplus;
+    }
+
+    // Update pressure gradient
+    gradP_.value() += gradPplus;
+
+    Info<< "Pressure gradient source: uncorrected Ubar = " << magUbarAve
+        << ", pressure gradient = " << gradP_.value() << endl;
+
+    writeGradP();
 }
 
 

src/fieldSources/basicSource/pressureGradientExplicitSource/pressureGradientExplicitSource.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) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -118,7 +118,10 @@ public:
 
     // Member Functions
 
-        // Access
+        // Evaluate
+
+            //- Correct the pressure gradient
+            virtual void correct(volVectorField& U);
 
             //- Add explicit contribution to equation
             virtual void addSup(fvMatrix<vector>& eqn, const label fieldI);

src/finiteVolume/Make/files

@@ -170,7 +170,7 @@ $(derivedFvPatchFields)/swirlFlowRateInletVelocity/swirlFlowRateInletVelocityFvP
 $(derivedFvPatchFields)/cylindricalInletVelocity/cylindricalInletVelocityFvPatchVectorField.C
 $(derivedFvPatchFields)/outletMappedUniformInlet/outletMappedUniformInletFvPatchFields.C
 $(derivedFvPatchFields)/waveSurfacePressure/waveSurfacePressureFvPatchScalarField.C
-
+$(derivedFvPatchFields)/phaseHydrostaticPressure/phaseHydrostaticPressureFvPatchScalarField.C
 
 fvsPatchFields = fields/fvsPatchFields
 $(fvsPatchFields)/fvsPatchField/fvsPatchFields.C

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
-    (
-        commsType,
-        this->patch().patchInternalField(psiInternal)()
-    );
+    this->patch().patchInternalField(psiInternal, scalarSendBuf_);
+
+    if (commsType == Pstream::nonBlocking && !Pstream::floatTransfer)
+    {
+        // 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
-    (
-        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
+    if (this->updatedMatrix())
+    {
+        return;
+    }
 
     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,
-        patch().patchInternalField(psiInternal)()
-    );
+        Pstream::defaultCommsType == Pstream::nonBlocking
+    && !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
-    (
-        procPatch_.compressedReceive<scalar>(commsType, this->size())()
-    );
-
-    const labelUList& faceCells = patch().faceCells();
-
-    forAll(faceCells, facei)
+    if (this->updatedMatrix())
     {
-        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/derived/phaseHydrostaticPressure/phaseHydrostaticPressureFvPatchScalarField.C

@@ -0,0 +1,202 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "phaseHydrostaticPressureFvPatchScalarField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+#include "uniformDimensionedFields.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::phaseHydrostaticPressureFvPatchScalarField::
+phaseHydrostaticPressureFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    mixedFvPatchScalarField(p, iF),
+    phaseName_("alpha"),
+    rho_(0.0),
+    pRefValue_(0.0),
+    pRefPoint_(vector::zero)
+{
+    this->refValue() = 0.0;
+    this->refGrad() = 0.0;
+    this->valueFraction() = 0.0;
+}
+
+
+Foam::phaseHydrostaticPressureFvPatchScalarField::
+phaseHydrostaticPressureFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    mixedFvPatchScalarField(p, iF),
+    phaseName_(dict.lookupOrDefault<word>("phaseName", "alpha")),
+    rho_(readScalar(dict.lookup("rho"))),
+    pRefValue_(readScalar(dict.lookup("pRefValue"))),
+    pRefPoint_(dict.lookup("pRefPoint"))
+{
+    this->refValue() = pRefValue_;
+
+    if (dict.found("value"))
+    {
+        fvPatchScalarField::operator=
+        (
+            scalarField("value", dict, p.size())
+        );
+    }
+    else
+    {
+        fvPatchScalarField::operator=(this->refValue());
+    }
+
+    this->refGrad() = 0.0;
+    this->valueFraction() = 0.0;
+}
+
+
+Foam::phaseHydrostaticPressureFvPatchScalarField::
+phaseHydrostaticPressureFvPatchScalarField
+(
+    const phaseHydrostaticPressureFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    mixedFvPatchScalarField(ptf, p, iF, mapper),
+    phaseName_(ptf.phaseName_),
+    rho_(ptf.rho_),
+    pRefValue_(ptf.pRefValue_),
+    pRefPoint_(ptf.pRefPoint_)
+{}
+
+
+Foam::phaseHydrostaticPressureFvPatchScalarField::
+phaseHydrostaticPressureFvPatchScalarField
+(
+    const phaseHydrostaticPressureFvPatchScalarField& ptf
+)
+:
+    mixedFvPatchScalarField(ptf),
+    phaseName_(ptf.phaseName_)
+{}
+
+
+Foam::phaseHydrostaticPressureFvPatchScalarField::
+phaseHydrostaticPressureFvPatchScalarField
+(
+    const phaseHydrostaticPressureFvPatchScalarField& ptf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    mixedFvPatchScalarField(ptf, iF),
+    phaseName_(ptf.phaseName_),
+    rho_(ptf.rho_),
+    pRefValue_(ptf.pRefValue_),
+    pRefPoint_(ptf.pRefPoint_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::phaseHydrostaticPressureFvPatchScalarField::updateCoeffs()
+{
+    if (this->updated())
+    {
+        return;
+    }
+
+    const scalarField& alphap =
+        patch().lookupPatchField<volScalarField, scalar>
+        (
+            phaseName_
+        );
+
+    const uniformDimensionedVectorField& g =
+        db().lookupObject<uniformDimensionedVectorField>("g");
+
+    // scalar rhor = 1000;
+    // scalarField alphap1 = max(min(alphap, 1.0), 0.0);
+    // valueFraction() = alphap1/(alphap1 + rhor*(1.0 - alphap1));
+    valueFraction() = max(min(alphap, 1.0), 0.0);
+
+    refValue() =
+        pRefValue_
+      + rho_*((g.value() & patch().Cf()) - (g.value() & pRefPoint_));
+
+    mixedFvPatchScalarField::updateCoeffs();
+}
+
+
+void Foam::phaseHydrostaticPressureFvPatchScalarField::write(Ostream& os) const
+{
+    fvPatchScalarField::write(os);
+    if (phaseName_ != "alpha")
+    {
+        os.writeKeyword("phaseName")
+            << phaseName_ << token::END_STATEMENT << nl;
+    }
+    os.writeKeyword("rho") << rho_ << token::END_STATEMENT << nl;
+    os.writeKeyword("pRefValue") << pRefValue_ << token::END_STATEMENT << nl;
+    os.writeKeyword("pRefPoint") << pRefPoint_ << token::END_STATEMENT << nl;
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+void Foam::phaseHydrostaticPressureFvPatchScalarField::operator=
+(
+    const fvPatchScalarField& ptf
+)
+{
+    fvPatchScalarField::operator=
+    (
+        valueFraction()*refValue()
+        + (1 - valueFraction())*ptf
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    makePatchTypeField
+    (
+        fvPatchScalarField,
+        phaseHydrostaticPressureFvPatchScalarField
+    );
+}
+
+// ************************************************************************* //

src/finiteVolume/fields/fvPatchFields/derived/phaseHydrostaticPressure/phaseHydrostaticPressureFvPatchScalarField.H

@@ -0,0 +1,223 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Class
+    Foam::phaseHydrostaticPressureFvPatchScalarField
+
+Description
+    Phase hydrostatic pressure boundary condition calculated as
+
+    pRefValue + rho*g.(x - pRefPoint)
+
+    where rho is provided and assumed uniform
+
+    applied according to the phase-fraction field provided:
+        1 -> fix value to that provided
+        0 -> zero-gradient
+
+SourceFiles
+    phaseHydrostaticPressureFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef phaseHydrostaticPressureFvPatchScalarField_H
+#define phaseHydrostaticPressureFvPatchScalarField_H
+
+#include "mixedFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+              Class phaseHydrostaticPressureFvPatch Declaration
+\*---------------------------------------------------------------------------*/
+
+class phaseHydrostaticPressureFvPatchScalarField
+:
+    public mixedFvPatchScalarField
+{
+
+protected:
+
+    // Protected data
+
+        //- Name of phase-fraction field
+        word phaseName_;
+
+        //- Constant density in the far-field
+        scalar rho_;
+
+        //- Reference pressure
+        scalar pRefValue_;
+
+        //- Reference pressure location
+        vector pRefPoint_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("phaseHydrostaticPressure");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        phaseHydrostaticPressureFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        phaseHydrostaticPressureFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        //  phaseHydrostaticPressureFvPatchScalarField onto a new patch
+        phaseHydrostaticPressureFvPatchScalarField
+        (
+            const phaseHydrostaticPressureFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        phaseHydrostaticPressureFvPatchScalarField
+        (
+            const phaseHydrostaticPressureFvPatchScalarField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField >
+            (
+                new phaseHydrostaticPressureFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        phaseHydrostaticPressureFvPatchScalarField
+        (
+            const phaseHydrostaticPressureFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new phaseHydrostaticPressureFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        // Access
+
+            //- Return the phaseName
+            const word& phaseName() const
+            {
+                return phaseName_;
+            }
+
+            //- Return reference to the phaseName to allow adjustment
+            word& phaseName()
+            {
+                return phaseName_;
+            }
+
+            //- Return the constant density in the far-field
+            scalar rho() const
+            {
+                return rho_;
+            }
+
+            //- Return reference to the constant density in the far-field
+            //  to allow adjustment
+            scalar& rho()
+            {
+                return rho_;
+            }
+
+            //- Return the reference pressure
+            scalar pRefValue() const
+            {
+                return pRefValue_;
+            }
+
+            //- Return reference to the reference pressure to allow adjustment
+            scalar& pRefValue()
+            {
+                return pRefValue_;
+            }
+
+            //- Return the pressure reference location
+            const vector& pRefPoint() const
+            {
+                return pRefPoint_;
+            }
+
+            //- Return reference to the pressure reference location
+            //  to allow adjustment
+            vector& pRefPoint()
+            {
+                return pRefPoint_;
+            }
+
+
+        //- Update the coefficients associated with the patch field
+        virtual void updateCoeffs();
+
+        //- Write
+        virtual void write(Ostream&) const;
+
+
+    // Member operators
+
+        virtual void operator=(const fvPatchScalarField& pvf);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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/constraint/cyclicAMI/cyclicAMIFvPatch.C

@@ -75,8 +75,8 @@ Foam::tmp<Foam::vectorField> Foam::cyclicAMIFvPatch::delta() const
     {
         forAll(patchD, faceI)
         {
-            vector ddi = patchD[faceI];
-            vector dni = nbrPatchD[faceI];
+            const vector& ddi = patchD[faceI];
+            const vector& dni = nbrPatchD[faceI];
 
             pdv[faceI] = ddi - dni;
         }
@@ -85,8 +85,8 @@ Foam::tmp<Foam::vectorField> Foam::cyclicAMIFvPatch::delta() const
     {
         forAll(patchD, faceI)
         {
-            vector ddi = patchD[faceI];
-            vector dni = nbrPatchD[faceI];
+            const vector& ddi = patchD[faceI];
+            const vector& dni = nbrPatchD[faceI];
 
             pdv[faceI] = ddi - transform(forwardT()[0], dni);
         }

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/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelIO.C

@@ -128,7 +128,7 @@ void Foam::ReactingParcel<ParcelType>::readFields
     wordList stateLabels(nPhases, "");
     if (compModel.nPhase() == 1)
     {
-        stateLabels = compModel.stateLabels();
+        stateLabels = compModel.stateLabels()[0];
     }
 
 
@@ -198,7 +198,7 @@ void Foam::ReactingParcel<ParcelType>::writeFields
         wordList stateLabels(phaseTypes.size(), "");
         if (compModel.nPhase() == 1)
         {
-            stateLabels = compModel.stateLabels();
+            stateLabels = compModel.stateLabels()[0];
         }
 
         forAll(phaseTypes, j)

src/lagrangian/intermediate/submodels/CloudFunctionObjects/PatchPostProcessing/PatchPostProcessing.C

@@ -26,6 +26,7 @@ License
 #include "PatchPostProcessing.H"
 #include "Pstream.H"
 #include "stringListOps.H"
+#include "ListOps.H"
 #include "ListListOps.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
@@ -55,9 +56,12 @@ void Foam::PatchPostProcessing<CloudType>::write()
 {
     forAll(patchData_, i)
     {
+        List<List<scalar> > procTimes(Pstream::nProcs());
+        procTimes[Pstream::myProcNo()] = times_[i];
+        Pstream::gatherList(procTimes);
+
         List<List<string> > procData(Pstream::nProcs());
         procData[Pstream::myProcNo()] = patchData_[i];
-
         Pstream::gatherList(procData);
 
         if (Pstream::master())
@@ -100,18 +104,33 @@ void Foam::PatchPostProcessing<CloudType>::write()
                 procData,
                 accessOp<List<string> >()
             );
-            sort(globalData);
+
+            List<scalar> globalTimes;
+            globalTimes = ListListOps::combine<List<scalar> >
+            (
+                procTimes,
+                accessOp<List<scalar> >()
+            );
+
+            labelList indices;
+            sortedOrder(globalTimes, indices);
 
             string header("# Time currentProc " + parcelType::propHeader);
             patchOutFile<< header.c_str() << nl;
 
-            forAll(globalData, dataI)
+            forAll(globalTimes, i)
             {
-                patchOutFile<< globalData[dataI].c_str() << nl;
+                label dataI = indices[i];
+
+                patchOutFile
+                    << globalTimes[dataI] << ' '
+                    << globalData[dataI].c_str()
+                    << nl;
             }
         }
 
         patchData_[i].clearStorage();
+        times_[i].clearStorage();
     }
 }
 
@@ -128,6 +147,7 @@ Foam::PatchPostProcessing<CloudType>::PatchPostProcessing
     CloudFunctionObject<CloudType>(dict, owner, typeName),
     maxStoredParcels_(readScalar(this->coeffDict().lookup("maxStoredParcels"))),
     patchIDs_(),
+    times_(),
     patchData_()
 {
     const wordList allPatchNames = owner.mesh().boundaryMesh().names();
@@ -167,6 +187,7 @@ Foam::PatchPostProcessing<CloudType>::PatchPostProcessing
     }
 
     patchData_.setSize(patchIDs_.size());
+    times_.setSize(patchIDs_.size());
 }
 
 
@@ -179,6 +200,7 @@ Foam::PatchPostProcessing<CloudType>::PatchPostProcessing
     CloudFunctionObject<CloudType>(ppm),
     maxStoredParcels_(ppm.maxStoredParcels_),
     patchIDs_(ppm.patchIDs_),
+    times_(ppm.times_),
     patchData_(ppm.patchData_)
 {}
 
@@ -203,9 +225,11 @@ void Foam::PatchPostProcessing<CloudType>::postPatch
     const label localPatchI = applyToPatch(patchI);
     if (localPatchI != -1 && patchData_[localPatchI].size() < maxStoredParcels_)
     {
+        times_[localPatchI].append(this->owner().time().value());
+
         OStringStream data;
-        data<< this->owner().time().timeName() << ' ' << Pstream::myProcNo()
-            << ' ' << p;
+        data<< Pstream::myProcNo() << ' ' << p;
+
         patchData_[localPatchI].append(data.str());
     }
 }

src/lagrangian/intermediate/submodels/CloudFunctionObjects/PatchPostProcessing/PatchPostProcessing.H

@@ -61,6 +61,9 @@ class PatchPostProcessing
         //- List of patch indices to post-process
         labelList patchIDs_;
 
+        //- List of time for each data record
+        List<DynamicList<scalar> > times_;
+
         //- List of output data per patch
         List<DynamicList<string> > patchData_;
 

src/lagrangian/molecularDynamics/molecule/reducedUnits/reducedUnits.H

@@ -76,7 +76,7 @@ class reducedUnits
 
                 scalar refPressure_;
 
-               scalar refMassDensity_;
+                scalar refMassDensity_;
 
                 scalar refNumberDensity_;
 
@@ -103,9 +103,9 @@ public:
     // Constructors
 
         //- Construct with no argument, uses default values:
-        // length  = 1nm
-        // mass = 1.660538782e−27kg (unified atomic mass unit)
-        // temperature = 1K (therefore, energy = 1*kb)
+        //  length  = 1nm
+        //  mass = 1.660538782e-27kg (unified atomic mass unit)
+        //  temperature = 1K (therefore, energy = 1*kb)
         reducedUnits();
 
         //- Construct from components

src/lagrangian/molecularDynamics/old/reducedUnits/reducedUnits.H

@@ -76,7 +76,7 @@ class reducedUnits
 
                 scalar refPressure_;
 
-               scalar refMassDensity_;
+                scalar refMassDensity_;
 
                 scalar refNumberDensity_;
 
@@ -103,9 +103,9 @@ public:
     // Constructors
 
         //- Construct with no argument, uses default values:
-        // length  = 1nm
-        // mass = 1.660538782e−27kg (unified atomic mass unit)
-        // temperature = 1K (therefore, energy = 1*kb)
+        //  length  = 1nm
+        //  mass = 1.660538782e-27kg (unified atomic mass unit)
+        //  temperature = 1K (therefore, energy = 1*kb)
         reducedUnits();
 
         //- Construct from components

src/lagrangian/spray/submodels/BreakupModel/PilchErdman/PilchErdman.C

@@ -101,25 +101,28 @@ bool Foam::PilchErdman<CloudType>::update
         // We > 1335, wave crest stripping
         scalar taubBar = 5.5;
 
-        if (We > 175.0)
+        if (We < 1335)
         {
-            // sheet stripping
-            taubBar = 0.766*pow(2.0*We - 12.0, 0.25);
-        }
-        else if (We > 22.0)
-        {
-            // Bag-and-stamen breakup
-            taubBar = 14.1*pow(2.0*We - 12.0, -0.25);
-        }
-        else if (We > 9.0)
-        {
-            // Bag breakup
-            taubBar = 2.45*pow(2.0*We - 12.0, 0.25);
-        }
-        else if (We > 6.0)
-        {
-            // Vibrational breakup
-            taubBar = 6.0*pow(2.0*We - 12.0, -0.25);
+            if (We > 175.0)
+            {
+                // sheet stripping
+                taubBar = 0.766*pow(2.0*We - 12.0, 0.25);
+            }
+            else if (We > 22.0)
+            {
+                // Bag-and-stamen breakup
+                taubBar = 14.1*pow(2.0*We - 12.0, -0.25);
+            }
+            else if (We > 9.0)
+            {
+                // Bag breakup
+                taubBar = 2.45*pow(2.0*We - 12.0, 0.25);
+            }
+            else if (We > 6.0)
+            {
+                // Vibrational breakup
+                taubBar = 6.0*pow(2.0*We - 12.0, -0.25);
+            }
         }
 
         scalar rho12 = sqrt(rhoc/rho);

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/AMIInterpolation/patches/cyclic/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C

@@ -366,7 +366,7 @@ Foam::cyclicAMIPolyPatch::cyclicAMIPolyPatch
     AMIPtr_(NULL),
     AMIReverse_(false),
     surfPtr_(NULL),
-    surfDict_(dictionary::null)
+    surfDict_(fileName("surface"))
 {
     // Neighbour patch might not be valid yet so no transformation
     // calculation possible
@@ -474,7 +474,7 @@ Foam::cyclicAMIPolyPatch::cyclicAMIPolyPatch
     AMIPtr_(NULL),
     AMIReverse_(pp.AMIReverse_),
     surfPtr_(NULL),
-    surfDict_(dictionary::null)
+    surfDict_(pp.surfDict_)
 {
     // Neighbour patch might not be valid yet so no transformation
     // calculation possible
@@ -501,7 +501,7 @@ Foam::cyclicAMIPolyPatch::cyclicAMIPolyPatch
     AMIPtr_(NULL),
     AMIReverse_(pp.AMIReverse_),
     surfPtr_(NULL),
-    surfDict_(dictionary::null)
+    surfDict_(pp.surfDict_)
 {
     if (nbrPatchName_ == name())
     {
@@ -825,7 +825,7 @@ void Foam::cyclicAMIPolyPatch::write(Ostream& os) const
             << token::END_STATEMENT << nl;
     }
 
-    if (surfDict_ != dictionary::null)
+    if (!surfDict_.empty())
     {
         os.writeKeyword(surfDict_.dictName());
         os  << surfDict_;

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]);
-    }
-    else
-    {
-        return cubeBb.overlaps(calcBb(edgeLabels_[index]));
+        return true;
     }
+
+    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/mappedPatches/mappedPolyPatch/mappedPatchBase.C

@@ -798,7 +798,7 @@ Foam::mappedPatchBase::mappedPatchBase
     AMIPtr_(NULL),
     AMIReverse_(false),
     surfPtr_(NULL),
-    surfDict_(dictionary::null)
+    surfDict_(fileName("surface"))
 {}
 
 
@@ -824,7 +824,7 @@ Foam::mappedPatchBase::mappedPatchBase
     AMIPtr_(NULL),
     AMIReverse_(false),
     surfPtr_(NULL),
-    surfDict_(dictionary::null)
+    surfDict_(fileName("surface"))
 {}
 
 
@@ -850,7 +850,7 @@ Foam::mappedPatchBase::mappedPatchBase
     AMIPtr_(NULL),
     AMIReverse_(false),
     surfPtr_(NULL),
-    surfDict_(dictionary::null)
+    surfDict_(fileName("surface"))
 {}
 
 
@@ -876,7 +876,7 @@ Foam::mappedPatchBase::mappedPatchBase
     AMIPtr_(NULL),
     AMIReverse_(false),
     surfPtr_(NULL),
-    surfDict_(dictionary::null)
+    surfDict_(fileName("surface"))
 {}
 
 
@@ -1216,8 +1216,11 @@ void Foam::mappedPatchBase::write(Ostream& os) const
                 << token::END_STATEMENT << nl;
         }
 
-        os.writeKeyword(surfDict_.dictName());
-        os  << surfDict_;
+        if (!surfDict_.empty())
+        {
+            os.writeKeyword(surfDict_.dictName());
+            os  << surfDict_;
+        }
     }
 }
 

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/meshTools/regionSplit/regionSplit.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
@@ -128,7 +128,7 @@ void Foam::regionSplit::fillSeedMask
     changedFaces.setSize(nFaces);
 
 
-    // Loop over changed faces. MeshWave in small.
+    // Loop over changed faces. FaceCellWave in small.
 
     while (changedFaces.size())
     {
@@ -260,7 +260,7 @@ void Foam::regionSplit::fillSeedMask
 }
 
 
-Foam::label Foam::regionSplit::calcRegionSplit
+Foam::label Foam::regionSplit::calcLocalRegionSplit
 (
     const boolList& blockedFace,
     const List<labelPair>& explicitConnections,
@@ -282,7 +282,7 @@ Foam::label Foam::regionSplit::calcRegionSplit
                 {
                     FatalErrorIn
                     (
-                        "regionSplit::calcRegionSplit(..)"
+                        "regionSplit::calcLocalRegionSplit(..)"
                     )   << "Face " << faceI << " not synchronised. My value:"
                         << blockedFace[faceI] << "  coupled value:"
                         << syncBlockedFace[faceI]
@@ -313,7 +313,7 @@ Foam::label Foam::regionSplit::calcRegionSplit
     // ~~~~~~~~~~~~~~~~~~~~
 
     // Start with region 0
-    label nRegions = 0;
+    label nLocalRegions = 0;
 
     label unsetCellI = 0;
 
@@ -340,11 +340,11 @@ Foam::label Foam::regionSplit::calcRegionSplit
             cellRegion,
             faceRegion,
             unsetCellI,
-            nRegions
+            nLocalRegions
         );
 
         // Current unsetCell has now been handled. Go to next region.
-        nRegions++;
+        nLocalRegions++;
         unsetCellI++;
     }
     while (true);
@@ -356,7 +356,7 @@ Foam::label Foam::regionSplit::calcRegionSplit
         {
             if (cellRegion[cellI] < 0)
             {
-                FatalErrorIn("regionSplit::calcRegionSplit(..)")
+                FatalErrorIn("regionSplit::calcLocalRegionSplit(..)")
                     << "cell:" << cellI << " region:" << cellRegion[cellI]
                     << abort(FatalError);
             }
@@ -366,34 +366,62 @@ Foam::label Foam::regionSplit::calcRegionSplit
         {
             if (faceRegion[faceI] == -1)
             {
-                FatalErrorIn("regionSplit::calcRegionSplit(..)")
+                FatalErrorIn("regionSplit::calcLocalRegionSplit(..)")
                     << "face:" << faceI << " region:" << faceRegion[faceI]
                     << abort(FatalError);
             }
         }
     }
 
+    return nLocalRegions;
+}
 
 
-    // Assign global regions
-    // ~~~~~~~~~~~~~~~~~~~~~
+Foam::autoPtr<Foam::globalIndex> Foam::regionSplit::calcRegionSplit
+(
+    const boolList& blockedFace,
+    const List<labelPair>& explicitConnections,
+
+    labelList& cellRegion
+) const
+{
+    // See header in regionSplit.H
+
+    // 1. Do local analysis
+    label nLocalRegions = calcLocalRegionSplit
+    (
+        blockedFace,
+        explicitConnections,
+        cellRegion
+    );
+
+    if (!Pstream::parRun())
+    {
+        return autoPtr<globalIndex>(new globalIndex(nLocalRegions));
+    }
+
+
+    // 2. Assign global regions
+    // ~~~~~~~~~~~~~~~~~~~~~~~~
     // Offset local regions to create unique global regions.
 
-    globalIndex globalRegions(nRegions);
+    globalIndex globalRegions(nLocalRegions);
 
 
-    // Merge global regions
-    // ~~~~~~~~~~~~~~~~~~~~
+    // Convert regions to global ones
+    forAll(cellRegion, cellI)
+    {
+        cellRegion[cellI] = globalRegions.toGlobal(cellRegion[cellI]);
+    }
+
+
+    // 3. Merge global regions
+    // ~~~~~~~~~~~~~~~~~~~~~~~
     // Regions across non-blocked proc patches get merged.
     // This will set merged global regions to be the min of both.
     // (this will create gaps in the global region list so they will get
     // merged later on)
 
-    // Map from global to merged global
-    labelList mergedGlobal(identity(globalRegions.size()));
-
-
-    // See if any regions get merged. Only nessecary for parallel
     while (Pstream::parRun())
     {
         if (debug)
@@ -401,123 +429,81 @@ Foam::label Foam::regionSplit::calcRegionSplit
             Pout<< nl << "-- Starting Iteration --" << endl;
         }
 
+
         const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
-        // Send global regions across (or -2 if blocked face)
+        labelList nbrRegion(mesh_.nFaces()-mesh_.nInternalFaces(), -1);
         forAll(patches, patchI)
         {
             const polyPatch& pp = patches[patchI];
 
-            if (isA<processorPolyPatch>(pp))
+            if (pp.coupled())
             {
-                labelList myGlobalRegions(pp.size());
-
-                label faceI = pp.start();
-
-                forAll(pp, i)
-                {
-                    if (faceRegion[faceI] < 0)
-                    {
-                        myGlobalRegions[i] = faceRegion[faceI];
-                    }
-                    else
-                    {
-                        myGlobalRegions[i] = mergedGlobal
-                        [globalRegions.toGlobal(faceRegion[faceI])];
-                    }
-
-                    faceI++;
-                }
-
-                OPstream toProcNbr
+                const labelList& patchCells = pp.faceCells();
+                SubList<label> patchNbrRegion
                 (
-                    Pstream::blocking,
-                    refCast<const processorPolyPatch>(pp).neighbProcNo()
+                    nbrRegion,
+                    pp.size(),
+                    pp.start()-mesh_.nInternalFaces()
                 );
 
-                toProcNbr << myGlobalRegions;
+                forAll(patchCells, i)
+                {
+                    label faceI = pp.start()+i;
+                    if (!blockedFace.size() || !blockedFace[faceI])
+                    {
+                        patchNbrRegion[i] = cellRegion[patchCells[i]];
+                    }
+                }
             }
         }
+        syncTools::swapBoundaryFaceList(mesh_, nbrRegion);
 
-
-        // Receive global regions
-
-        label nMerged = 0;
+        Map<label> globalToMerged(mesh_.nFaces()-mesh_.nInternalFaces());
 
         forAll(patches, patchI)
         {
             const polyPatch& pp = patches[patchI];
 
-            if (isA<processorPolyPatch>(pp))
+            if (pp.coupled())
             {
-                const processorPolyPatch& procPp =
-                    refCast<const processorPolyPatch>(pp);
+                const labelList& patchCells = pp.faceCells();
+                SubList<label> patchNbrRegion
+                (
+                    nbrRegion,
+                    pp.size(),
+                    pp.start()-mesh_.nInternalFaces()
+                );
 
-                IPstream fromProcNbr(Pstream::blocking, procPp.neighbProcNo());
-
-                labelList nbrRegions(fromProcNbr);
-
-
-                // Compare with my regions to see which get merged.
-
-                label faceI = pp.start();
-
-                forAll(pp, i)
+                forAll(patchCells, i)
                 {
-                    if
-                    (
-                        faceRegion[faceI] < 0
-                     || nbrRegions[i] < 0
-                    )
+                    label faceI = pp.start()+i;
+
+                    if (!blockedFace.size() || !blockedFace[faceI])
                     {
-                        if (faceRegion[faceI] != nbrRegions[i])
+                        if (patchNbrRegion[i] < cellRegion[patchCells[i]])
                         {
-                            FatalErrorIn("regionSplit::calcRegionSplit(..)")
-                                << "On patch:" << pp.name()
-                                << " face:" << faceI
-                                << " my local region:" << faceRegion[faceI]
-                                << " neighbouring region:"
-                                << nbrRegions[i] << nl
-                                << "Maybe your blockedFaces are not"
-                                << " synchronized across coupled faces?"
-                                << abort(FatalError);
+                            //Pout<< "on patch:" << pp.name()
+                            //    << " cell:" << patchCells[i]
+                            //    << " at:"
+                            // << mesh_.cellCentres()[patchCells[i]]
+                            //    << " was:" << cellRegion[patchCells[i]]
+                            //    << " nbr:" << patchNbrRegion[i]
+                            //    << endl;
+
+                            globalToMerged.insert
+                            (
+                                cellRegion[patchCells[i]],
+                                patchNbrRegion[i]
+                            );
                         }
                     }
-                    else
-                    {
-                        label uncompactGlobal =
-                            globalRegions.toGlobal(faceRegion[faceI]);
-
-                        label myGlobal = mergedGlobal[uncompactGlobal];
-
-                        if (myGlobal != nbrRegions[i])
-                        {
-                            label minRegion = min(myGlobal, nbrRegions[i]);
-
-                            if (debug)
-                            {
-                                Pout<< "Merging region " << myGlobal
-                                    << " (on proc " << Pstream::myProcNo()
-                                    << ") and region " << nbrRegions[i]
-                                    << " (on proc " << procPp.neighbProcNo()
-                                    << ") into region " << minRegion << endl;
-                            }
-
-                            mergedGlobal[uncompactGlobal] = minRegion;
-                            mergedGlobal[myGlobal] = minRegion;
-                            mergedGlobal[nbrRegions[i]] = minRegion;
-
-                            nMerged++;
-                        }
-                    }
-
-                    faceI++;
                 }
             }
         }
 
 
-        reduce(nMerged, sumOp<label>());
+        label nMerged = returnReduce(globalToMerged.size(), sumOp<label>());
 
         if (debug)
         {
@@ -529,51 +515,176 @@ Foam::label Foam::regionSplit::calcRegionSplit
             break;
         }
 
-        // Merge the compacted regions.
-        Pstream::listCombineGather(mergedGlobal, minEqOp<label>());
-        Pstream::listCombineScatter(mergedGlobal);
+        // Renumber the regions according to the globalToMerged
+        forAll(cellRegion, cellI)
+        {
+            label regionI = cellRegion[cellI];
+            Map<label>::const_iterator iter = globalToMerged.find(regionI);
+            if (iter != globalToMerged.end())
+            {
+                 cellRegion[cellI] = iter();
+            }
+        }
     }
 
 
-    // Compact global regions
-    // ~~~~~~~~~~~~~~~~~~~~~~
+    // Now our cellRegion will have non-local elements in it. So compact
+    // it.
 
-    // All procs will have the same global mergedGlobal region.
-    // There might be gaps in it however so compact.
-
-    labelList mergedToCompacted(globalRegions.size(), -1);
-
-    label compactI = 0;
-
-    forAll(mergedGlobal, i)
+    // 4a: count. Use a labelHashSet to count regions only once.
+    label nCompact = 0;
     {
-        label merged = mergedGlobal[i];
-
-        if (mergedToCompacted[merged] == -1)
+        labelHashSet localRegion(mesh_.nFaces()-mesh_.nInternalFaces());
+        forAll(cellRegion, cellI)
         {
-            mergedToCompacted[merged] = compactI++;
+            if
+            (
+                globalRegions.isLocal(cellRegion[cellI])
+             && localRegion.insert(cellRegion[cellI])
+            )
+            {
+                nCompact++;
+            }
         }
     }
 
     if (debug)
     {
-        Pout<< "Compacted down to " << compactI << " regions." << endl;
+        Pout<< "Compacted from " << nLocalRegions
+            << " down to " << nCompact << " local regions." << endl;
+    }
+
+
+    // Global numbering for compacted local regions
+    autoPtr<globalIndex> globalCompactPtr(new globalIndex(nCompact));
+    const globalIndex& globalCompact = globalCompactPtr();
+
+
+    // 4b: renumber
+    // Renumber into compact indices. Note that since we've already made
+    // all regions global we now need a Map to store the compacting information
+    // instead of a labelList - otherwise we could have used a straight
+    // labelList.
+
+    // Local compaction map
+    Map<label> globalToCompact(2*nCompact);
+    // Remote regions we want the compact number for
+    List<labelHashSet> nonLocal(Pstream::nProcs());
+    forAll(nonLocal, procI)
+    {
+        nonLocal[procI].resize((nLocalRegions-nCompact)/Pstream::nProcs());
     }
 
-    // Renumber cellRegion to be global regions
     forAll(cellRegion, cellI)
     {
         label region = cellRegion[cellI];
-
-        if (region >= 0)
+        if (globalRegions.isLocal(region))
         {
-            label merged = mergedGlobal[globalRegions.toGlobal(region)];
-
-            cellRegion[cellI] = mergedToCompacted[merged];
+            Map<label>::const_iterator iter = globalToCompact.find(region);
+            if (iter == globalToCompact.end())
+            {
+                label compactRegion = globalCompact.toGlobal
+                (
+                    globalToCompact.size()
+                );
+                globalToCompact.insert(region, compactRegion);
+            }
+        }
+        else
+        {
+            nonLocal[globalRegions.whichProcID(region)].insert(region);
         }
     }
 
-    return compactI;
+
+    // Now we have all the local regions compacted. Now we need to get the
+    // non-local ones from the processors to whom they are local.
+    // Convert the nonLocal (labelHashSets) to labelLists.
+
+    labelListList sendNonLocal(Pstream::nProcs());
+    labelList nNonLocal(Pstream::nProcs(), 0);
+    forAll(sendNonLocal, procI)
+    {
+        sendNonLocal[procI].setSize(nonLocal[procI].size());
+        forAllConstIter(labelHashSet, nonLocal[procI], iter)
+        {
+            sendNonLocal[procI][nNonLocal[procI]++] = iter.key();
+        }
+    }
+
+    if (debug)
+    {
+        forAll(nNonLocal, procI)
+        {
+            Pout<< "    from processor " << procI
+                << " want " << nNonLocal[procI]
+                << " region numbers."
+                << endl;
+        }
+        Pout<< endl;
+    }
+
+
+    // Get the wanted region labels into recvNonLocal
+    labelListList recvNonLocal;
+    labelListList sizes;
+    Pstream::exchange<labelList, label>
+    (
+        sendNonLocal,
+        recvNonLocal,
+        sizes
+    );
+
+    // Now we have the wanted compact region labels that procI wants in
+    // recvNonLocal[procI]. Construct corresponding list of compact
+    // region labels to send back.
+
+    labelListList sendWantedLocal(Pstream::nProcs());
+    forAll(recvNonLocal, procI)
+    {
+        const labelList& nonLocal = recvNonLocal[procI];
+        sendWantedLocal[procI].setSize(nonLocal.size());
+
+        forAll(nonLocal, i)
+        {
+            sendWantedLocal[procI][i] = globalToCompact[nonLocal[i]];
+        }
+    }
+
+
+    // Send back (into recvNonLocal)
+    recvNonLocal.clear();
+    sizes.clear();
+    Pstream::exchange<labelList, label>
+    (
+        sendWantedLocal,
+        recvNonLocal,
+        sizes
+    );
+    sendWantedLocal.clear();
+
+    // Now recvNonLocal contains for every element in setNonLocal the
+    // corresponding compact number. Insert these into the local compaction
+    // map.
+
+    forAll(recvNonLocal, procI)
+    {
+        const labelList& wantedRegions = sendNonLocal[procI];
+        const labelList& compactRegions = recvNonLocal[procI];
+
+        forAll(wantedRegions, i)
+        {
+            globalToCompact.insert(wantedRegions[i], compactRegions[i]);
+        }
+    }
+
+    // Finally renumber the regions
+    forAll(cellRegion, cellI)
+    {
+        cellRegion[cellI] = globalToCompact[cellRegion[cellI]];
+    }
+
+    return globalCompactPtr;
 }
 
 
@@ -582,9 +693,15 @@ Foam::label Foam::regionSplit::calcRegionSplit
 Foam::regionSplit::regionSplit(const polyMesh& mesh)
 :
     labelList(mesh.nCells(), -1),
-    mesh_(mesh),
-    nRegions_(calcRegionSplit(boolList(0, false), List<labelPair>(0), *this))
-{}
+    mesh_(mesh)
+{
+    globalNumberingPtr_ = calcRegionSplit
+    (
+        boolList(0, false), //blockedFaces
+        List<labelPair>(0), //explicitConnections,
+        *this
+    );
+}
 
 
 Foam::regionSplit::regionSplit
@@ -594,9 +711,15 @@ Foam::regionSplit::regionSplit
 )
 :
     labelList(mesh.nCells(), -1),
-    mesh_(mesh),
-    nRegions_(calcRegionSplit(blockedFace, List<labelPair>(0), *this))
-{}
+    mesh_(mesh)
+{
+    globalNumberingPtr_ = calcRegionSplit
+    (
+        blockedFace,        //blockedFaces
+        List<labelPair>(0), //explicitConnections,
+        *this
+    );
+}
 
 
 Foam::regionSplit::regionSplit
@@ -607,12 +730,20 @@ Foam::regionSplit::regionSplit
 )
 :
     labelList(mesh.nCells(), -1),
-    mesh_(mesh),
-    nRegions_(calcRegionSplit(blockedFace, explicitConnections, *this))
-{}
+    mesh_(mesh)
+{
+    globalNumberingPtr_ = calcRegionSplit
+    (
+        blockedFace,            //blockedFaces
+        explicitConnections,    //explicitConnections,
+        *this
+    );
+}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 
+
+
 // ************************************************************************* //

src/meshTools/regionSplit/regionSplit.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
@@ -26,7 +26,66 @@ Class
 
 Description
     This class separates the mesh into distinct unconnected regions,
-    each of which is then given a label.
+    each of which is then given a label according to globalNumbering().
+
+
+    Say 6 cells, 3 processors, with single baffle on proc1.
+
+              baffle
+                |
+    +---+---+---+---+---+---+
+    |   |   |   |   |   |   |
+    +---+---+---+---+---+---+
+      proc0 | proc1 | proc2
+
+
+
+    1: determine local regions (uncoupled)
+
+    +---+---+---+---+---+---+
+    | 0 | 0 | 0 | 1 | 0 | 0 |
+    +---+---+---+---+---+---+
+      proc0 | proc1 | proc2
+
+
+
+    2: make global
+
+    +---+---+---+---+---+---+
+    | 0 | 0 | 1 | 2 | 3 | 3 |
+    +---+---+---+---+---+---+
+      proc0 | proc1 | proc2
+
+
+
+    3: merge connected across procs
+
+    +---+---+---+---+---+---+
+    | 0 | 0 | 0 | 2 | 2 | 2 |
+    +---+---+---+---+---+---+
+      proc0 | proc1 | proc2
+
+
+
+    4. determine locally owner regions. determine compact numbering for the
+    local regions and send these to all processors that need them:
+
+    proc0 uses regions:
+        - 0 which is local to it.
+    proc1 uses regions
+        - 0 which originates from proc0
+        - 2 which is local to it
+    proc2 uses regions
+        - 2 which originates from proc1
+
+    So proc1 needs to get the compact number for region 0 from proc0 and proc2
+    needs to get the compact number for region 2 from proc1:
+
+    +---+---+---+---+---+---+
+    | 0 | 0 | 0 | 1 | 1 | 1 |
+    +---+---+---+---+---+---+
+      proc0 | proc1 | proc2
+
 
 SourceFiles
     regionSplit.C
@@ -36,14 +95,17 @@ SourceFiles
 #ifndef regionSplit_H
 #define regionSplit_H
 
-#include "polyMesh.H"
-#include "demandDrivenData.H"
+#include "globalIndex.H"
+#include "labelPair.H"
+#include "boolList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
+class polyMesh;
+
 /*---------------------------------------------------------------------------*\
                            Class regionSplit Declaration
 \*---------------------------------------------------------------------------*/
@@ -57,8 +119,7 @@ class regionSplit
         //- Reference to mesh
         const polyMesh& mesh_;
 
-        //- Number of regions
-        label nRegions_;
+        mutable autoPtr<globalIndex> globalNumberingPtr_;
 
 
     // Private Member Functions
@@ -84,8 +145,16 @@ class regionSplit
             const label markValue
         ) const;
 
-        //- Calculate region split. Return number of regions.
-        label calcRegionSplit
+        //- Calculate local region split. Return number of regions.
+        label calcLocalRegionSplit
+        (
+            const boolList& blockedFace,
+            const List<labelPair>& explicitConnections,
+            labelList& cellRegion
+        ) const;
+
+        //- Calculate global region split. Return globalIndex.
+        autoPtr<globalIndex> calcRegionSplit
         (
             const boolList& blockedFace,
             const List<labelPair>& explicitConnections,
@@ -118,11 +187,18 @@ public:
 
     // Member Functions
 
-        //- Return number of regions
+        //- Return global region numbering
+        const globalIndex& globalNumbering() const
+        {
+            return globalNumberingPtr_();
+        }
+
+        //- Return total number of regions
         label nRegions() const
         {
-            return nRegions_;
+            return globalNumbering().size();
         }
+
 };
 
 

src/regionModels/regionModel/regionProperties/regionProperties.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
@@ -52,4 +52,5 @@ Foam::regionProperties::regionProperties(const Time& runTime)
 Foam::regionProperties::~regionProperties()
 {}
 
+
 // ************************************************************************* //

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);
-    storeTetOp cutInside(cutInsideTets, nCutInside);
-    dummyTetOp outside;
-
+    tetPointRef::storeOp inside(insideTets, nInside);
+    tetPointRef::storeOp cutInside(cutInsideTets, nCutInside);
+    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::tetOverlapVolume::tetTetOverlapVol
+Foam::scalar 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);
-    storeTetOp cutInside(cutInsideTets, nCutInside);
-    sumTetVolOp volInside;
-    dummyTetOp outside;
+    tetPointRef::storeOp inside(insideTets, nInside);
+    tetPointRef::storeOp cutInside(cutInsideTets, nCutInside);
+    tetPointRef::sumVolOp volInside;
+    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 fvMesh& toMesh,
+    const polyMesh& fromMesh,
+    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 "fvMesh.H"
-#include "plane.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"
+#include "FixedList.H"
+#include "labelList.H"
+#include "treeBoundBox.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 fvMesh& meshB,
+            const polyMesh& meshA,
+            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
-
-// ************************************************************************* //

src/thermophysicalModels/basic/Make/files

@@ -2,6 +2,7 @@ mixtures/basicMixture/basicMixture.C
 mixtures/basicMixture/basicMixtures.C
 
 basicThermo/basicThermo.C
+basicThermo/basicThermoNew.C
 
 psiThermo/basicPsiThermo/basicPsiThermo.C
 psiThermo/basicPsiThermo/basicPsiThermoNew.C

src/thermophysicalModels/basic/basicThermo/basicThermo.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
@@ -39,6 +39,7 @@ License
 namespace Foam
 {
     defineTypeNameAndDebug(basicThermo, 0);
+    defineRunTimeSelectionTable(basicThermo, fvMesh);
 }
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //

src/thermophysicalModels/basic/basicThermo/basicThermo.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
@@ -29,7 +29,7 @@ Description
 
 SourceFiles
     basicThermo.C
-    newBasicThermo.C
+    basicThermoNew.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -109,12 +109,25 @@ public:
     TypeName("basicThermo");
 
 
+    //- Declare run-time constructor selection table
+    declareRunTimeSelectionTable
+    (
+        autoPtr,
+        basicThermo,
+        fvMesh,
+        (const fvMesh& mesh),
+        (mesh)
+    );
+
     // Constructors
 
         //- Construct from mesh
         basicThermo(const fvMesh&);
 
 
+    //- Selector
+    static autoPtr<basicThermo> New(const fvMesh&);
+
     //- Destructor
     virtual ~basicThermo();
 

src/thermophysicalModels/basic/basicThermo/basicThermoNew.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
@@ -21,32 +21,51 @@ License
     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"
+#include "basicThermo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-namespace Foam
+Foam::autoPtr<Foam::basicThermo> Foam::basicThermo::New
+(
+    const fvMesh& mesh
+)
 {
+    // get model name, but do not register the dictionary
+    // otherwise it is registered in the database twice
+    const word modelType
+    (
+        IOdictionary
+        (
+            IOobject
+            (
+                "thermophysicalProperties",
+                mesh.time().constant(),
+                mesh,
+                IOobject::MUST_READ_IF_MODIFIED,
+                IOobject::NO_WRITE,
+                false
+            )
+        ).lookup("thermoType")
+    );
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    Info<< "Selecting thermodynamics package " << modelType << endl;
 
-typedef tetrahedron<point, const point&> tetPointRef;
+    fvMeshConstructorTable::iterator cstrIter =
+        fvMeshConstructorTablePtr_->find(modelType);
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    if (cstrIter == fvMeshConstructorTablePtr_->end())
+    {
+        FatalErrorIn("basicThermo::New(const fvMesh&)")
+            << "Unknown basicThermo type " << modelType << nl << nl
+            << "Valid basicThermo types are:" << nl
+            << fvMeshConstructorTablePtr_->sortedToc() << nl
+            << exit(FatalError);
+    }
 
-} // End namespace Foam
+    return autoPtr<basicThermo>(cstrIter()(mesh));
+}
 
-#endif
 
 // ************************************************************************* //

src/thermophysicalModels/basic/psiThermo/basicPsiThermo/makeBasicPsiThermo.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
@@ -55,7 +55,14 @@ addToRunTimeSelectionTable                                                    \
     basicPsiThermo,                                                           \
     Cthermo##Mixture##Transport##Thermo##EqnOfState,                          \
     fvMesh                                                                    \
-)
+);                                                                            \
+                                                                              \
+addToRunTimeSelectionTable                                                    \
+(                                                                             \
+    basicThermo,                                                              \
+    Cthermo##Mixture##Transport##Thermo##EqnOfState,                          \
+    fvMesh                                                                    \
+);
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/thermophysicalModels/basic/rhoThermo/basicRhoThermo/makeBasicRhoThermo.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
@@ -55,7 +55,14 @@ addToRunTimeSelectionTable                                                    \
     basicRhoThermo,                                                           \
     Cthermo##Mixture##Transport##Thermo##EqnOfState,                          \
     fvMesh                                                                    \
-)
+);                                                                            \
+                                                                              \
+addToRunTimeSelectionTable                                                    \
+(                                                                             \
+    basicThermo,                                                              \
+    Cthermo##Mixture##Transport##Thermo##EqnOfState,                          \
+    fvMesh                                                                    \
+);
 
 
 #define makeBasicRhoPolyThermo(Cthermo,Mixture,Order)                         \
@@ -88,7 +95,14 @@ addToRunTimeSelectionTable                                                    \
     basicRhoThermo,                                                           \
     Cthermo##Mixture##icoPoly##Order##ThermoPhysics,                          \
     fvMesh                                                                    \
-)
+);                                                                            \
+                                                                              \
+addToRunTimeSelectionTable                                                    \
+(                                                                             \
+    basicThermo,                                                              \
+    Cthermo##Mixture##icoPoly##Order##ThermoPhysics,                          \
+    fvMesh                                                                    \
+);
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/turbulenceModels/compressible/LES/Smagorinsky/Smagorinsky.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) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -113,7 +113,13 @@ public:
         //  calculated from the given velocity gradient
         tmp<volScalarField> k(const tmp<volTensorField>& gradU) const
         {
-            return (2.0*ck_/ce_)*sqr(delta())*magSqr(dev(symm(gradU)));
+            volSymmTensorField D(symm(gradU));
+
+            volScalarField a(ce_/delta());
+            volScalarField b((2.0/3.0)*tr(D));
+            volScalarField c(2*ck_*delta()*(dev(D) && D));
+
+            return sqr((-b + sqrt(sqr(b) + 4*a*c))/(2*a));
         }
 
         //- Return SGS kinetic energy

src/turbulenceModels/compressible/RAS/SpalartAllmaras/SpalartAllmaras.H

@@ -33,7 +33,7 @@ Description
         "A One-Equation Turbulence Model for Aerodynamic Flows"
         P.R. Spalart,
         S.R. Allmaras,
-        La Recherche A´rospatiale, No. 1, 1994, pp. 5–21.
+        La Recherche Aerospatiale, No. 1, 1994, pp. 5-21.
 
         Extended according to:
 

src/turbulenceModels/incompressible/RAS/SpalartAllmaras/SpalartAllmaras.H

@@ -33,7 +33,7 @@ Description
         "A One-Equation Turbulence Model for Aerodynamic Flows"
         P.R. Spalart,
         S.R. Allmaras,
-        La Recherche A´rospatiale, No. 1, 1994, pp. 5–21.
+        La Recherche Aerospatiale, No. 1, 1994, pp. 5-21.
 
         Extended according to: