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

applications/solvers/combustion/chemFoam/readInitialConditions.H

@@ -87,7 +87,7 @@
     forAll(Y, i)
     {
         Y[i] = Y0[i];
-        h0 += Y0[i]*specieData[i].Hs(p[i], T0);
+        h0 += Y0[i]*specieData[i].Hs(p[0], T0);
     }
 
     thermo.he() = dimensionedScalar("h", dimEnergy/dimMass, h0);

applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/Make/options

@@ -1,5 +1,8 @@
 EXE_INC = \
     -I../buoyantBoussinesqSimpleFoam \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/turbulenceModels \
     -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/lnInclude \
@@ -9,6 +12,8 @@ EXE_INC = \
 
 EXE_LIBS = \
     -lfiniteVolume \
+    -lfvOptions \
+    -lsampling \
     -lmeshTools \
     -lincompressibleTurbulenceModel \
     -lincompressibleRASModels \

applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/TEqn.H

@@ -11,12 +11,18 @@
       - fvm::laplacian(alphaEff, T)
      ==
         radiation->ST(rhoCpRef, T)
+      + fvOptions(T)
     );
 
     TEqn.relax();
+
+    fvOptions.constrain(TEqn);
+
     TEqn.solve();
 
     radiation->correct();
 
+    fvOptions.correct(T);
+
     rhok = 1.0 - beta*(T - TRef);
 }

applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/UEqn.H

@@ -5,10 +5,14 @@
         fvm::ddt(U)
       + fvm::div(phi, U)
       + turbulence->divDevReff(U)
+     ==
+        fvOptions(U)
     );
 
     UEqn.relax();
 
+    fvOptions.constrain(UEqn);
+
     if (pimple.momentumPredictor())
     {
         solve
@@ -23,4 +27,6 @@
                 )*mesh.magSf()
             )
         );
+
+        fvOptions.correct(U);
     }

applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C

@@ -48,8 +48,9 @@ Description
 #include "fvCFD.H"
 #include "singlePhaseTransportModel.H"
 #include "RASModel.H"
-#include "pimpleControl.H"
 #include "radiationModel.H"
+#include "fvIOoptionList.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -61,6 +62,7 @@ int main(int argc, char *argv[])
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "createIncompressibleRadiationModel.H"
+    #include "createFvOptions.H"
     #include "initContinuityErrs.H"
     #include "readTimeControls.H"
     #include "CourantNo.H"

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/Make/options

@@ -1,5 +1,8 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude \
     -I$(LIB_SRC)/turbulenceModels \
     -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/lnInclude \
     -I$(LIB_SRC)/transportModels \
@@ -7,7 +10,9 @@ EXE_INC = \
 
 EXE_LIBS = \
     -lfiniteVolume \
+    -lsampling \
     -lmeshTools \
+    -lfvOptions \
     -lincompressibleTurbulenceModel \
     -lincompressibleRASModels \
     -lincompressibleTransportModels

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/TEqn.H

@@ -8,10 +8,17 @@
     (
         fvm::div(phi, T)
       - fvm::laplacian(alphaEff, T)
+     ==
+        fvOptions(T)
     );
 
     TEqn.relax();
+
+    fvOptions.constrain(TEqn);
+
     TEqn.solve();
 
+    fvOptions.correct(T);
+
     rhok = 1.0 - beta*(T - TRef);
 }

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H

@@ -4,16 +4,20 @@
     (
         fvm::div(phi, U)
       + turbulence->divDevReff(U)
+     ==
+        fvOptions(U)
     );
 
     UEqn().relax();
 
+    fvOptions.constrain(UEqn());
+
     if (simple.momentumPredictor())
     {
         solve
         (
             UEqn()
-            ==
+          ==
             fvc::reconstruct
             (
                 (
@@ -22,4 +26,6 @@
                 )*mesh.magSf()
             )
         );
+
+        fvOptions.correct(U);
     }

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C

@@ -48,6 +48,7 @@ Description
 #include "fvCFD.H"
 #include "singlePhaseTransportModel.H"
 #include "RASModel.H"
+#include "fvIOoptionList.H"
 #include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -59,6 +60,7 @@ int main(int argc, char *argv[])
     #include "createMesh.H"
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
+    #include "createFvOptions.H"
     #include "initContinuityErrs.H"
 
     simpleControl simple(mesh);

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C

@@ -26,21 +26,18 @@ Application
 
 Description
     Combination of heatConductionFoam and buoyantFoam for conjugate heat
-    transfer between a solid region and fluid region. It includes
+    transfer between solid regions and fluid regions. Both regions include
-    porous media in the primary fluid region treated explicitly.
+    the fvOptions framework.
 
     It handles secondary fluid or solid circuits which can be coupled
     thermally with the main fluid region. i.e radiators, etc.
 
-    The secondary fluid region is
-
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
 #include "rhoThermo.H"
 #include "turbulenceModel.H"
 #include "fixedGradientFvPatchFields.H"
-#include "zeroGradientFvPatchFields.H"
 #include "regionProperties.H"
 #include "compressibleCourantNo.H"
 #include "solidRegionDiffNo.H"

applications/solvers/incompressible/pimpleFoam/pimpleFoam.C

@@ -30,7 +30,7 @@ Description
 
     Sub-models include:
     - turbulence modelling, i.e. laminar, RAS or LES
-    - run-time selectable finitie volume options, e.g. MRF, explicit porosity
+    - run-time selectable finite volume options, e.g. MRF, explicit porosity
 
 \*---------------------------------------------------------------------------*/
 

applications/test/GAMGAgglomeration/Test-GAMGAgglomeration.C

@@ -35,116 +35,6 @@ Description
 #include "meshTools.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-void checkConnectedAgglomeration
-(
-    const lduMesh& mesh,
-    const labelUList& restrict,
-    const label nCoarse
-)
-{
-    if (mesh.lduAddr().size() != restrict.size())
-    {
-        FatalErrorIn
-        (
-            "checkConnectedAgglomeration(const lduMesh&, const labelList&)"
-        )   << "nCells:" << mesh.lduAddr().size()
-            << " agglom:" << restrict.size()
-            << abort(FatalError);
-    }
-
-    // Seed (master) for every region
-    labelList regionToMaster(nCoarse, -1);
-    labelList master(mesh.lduAddr().size(), -1);
-    forAll(restrict, cellI)
-    {
-        label region = restrict[cellI];
-        if (regionToMaster[region] == -1)
-        {
-            // Set cell to be master for region
-            //Pout<< "For region " << region
-            //    << " allocating local master " << cellI
-            //    << endl;
-            regionToMaster[region] = cellI;
-            master[cellI] = cellI;
-        }
-    }
-
-    // Now loop and transport master through region
-    const labelUList& lower = mesh.lduAddr().lowerAddr();
-    const labelUList& upper = mesh.lduAddr().upperAddr();
-
-    while (true)
-    {
-        label nChanged = 0;
-
-        forAll(lower, faceI)
-        {
-            label own = lower[faceI];
-            label nei = upper[faceI];
-
-            if (restrict[own] == restrict[nei])
-            {
-                // Region-internal face
-
-                if (master[own] != -1)
-                {
-                    if (master[nei] == -1)
-                    {
-                        master[nei] = master[own];
-                        nChanged++;
-                    }
-                    else if (master[nei] != master[own])
-                    {
-                        FatalErrorIn("checkConnectedAgglomeration(..)")
-                            << "problem" << abort(FatalError);
-                    }
-                }
-                else if (master[nei] != -1)
-                {
-                    master[own] = master[nei];
-                    nChanged++;
-                }
-            }
-        }
-
-        reduce(nChanged, sumOp<label>());
-
-        if (nChanged == 0)
-        {
-            break;
-        }
-    }
-
-    // Check that master is set for all cells
-    boolList singleRegion(nCoarse, true);
-    label nSet = nCoarse;
-    forAll(master, cellI)
-    {
-        if (master[cellI] == -1)
-        {
-            label region = restrict[cellI];
-            if (singleRegion[region] == true)
-            {
-                singleRegion[region] = false;
-                nSet--;
-            }
-        }
-    }
-
-    label totalNCoarse = returnReduce(nCoarse, sumOp<label>());
-    label totalNVisited = returnReduce(nSet, sumOp<label>());
-
-    if (totalNVisited < totalNCoarse)
-    {
-        WarningIn("checkConnectedAgglomeration(..)")
-            << "out of " << totalNCoarse
-            << " agglomerated cells have " << totalNCoarse-totalNVisited
-            << " cells that are not a single connected region" << endl;
-    }
-}
-
-
 // Main program:
 
 int main(int argc, char *argv[])
@@ -227,12 +117,28 @@ int main(int argc, char *argv[])
             << "    agglomerated size : "
             << returnReduce(coarseSize, sumOp<label>()) << endl;
 
-        checkConnectedAgglomeration
+        labelList newAddr;
+        label newCoarseSize = 0;
+        bool ok = GAMGAgglomeration::checkRestriction
         (
-            agglom.meshLevel(level),
+            newAddr,
+            newCoarseSize,
+
+            agglom.meshLevel(level).lduAddr(),
             addr,
             coarseSize
         );
+        if (!ok)
+        {
+            WarningIn(args.executable())
+                << "At level " << level
+                << " there are " << coarseSize
+                << " agglomerated cells but " << newCoarseSize
+                << " disconnected regions" << endl
+                << "    This means that some agglomerations (coarse cells)"
+                << "    consist of multiple disconnected regions."
+                << endl;
+        }
 
 
         forAll(addr, fineI)

applications/test/Matrix/Test-Matrix.C

@@ -116,17 +116,19 @@ int main(int argc, char *argv[])
         squareMatrix[2][1] = -43;
         squareMatrix[2][2] = 98;
 
+        const scalarSquareMatrix squareMatrixCopy = squareMatrix;
         Info<< nl << "Square Matrix = " << squareMatrix << endl;
 
-        scalarDiagonalMatrix rhs(3, 0);
+        Info<< "det = " << det(squareMatrixCopy) << endl;
-        rhs[0] = 1;
-        rhs[1] = 2;
-        rhs[2] = 3;
 
-        LUsolve(squareMatrix, rhs);
+        labelList rhs(3, 0);
+        label sign;
+        LUDecompose(squareMatrix, rhs, sign);
 
         Info<< "Decomposition = " << squareMatrix << endl;
-        Info<< "Solution = " << rhs << endl;
+        Info<< "Pivots = " << rhs << endl;
+        Info<< "Sign = " << sign << endl;
+        Info<< "det = " << detDecomposed(squareMatrix, sign) << endl;
     }
 
     Info<< "\nEnd\n" << endl;

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DistributedDelaunayMesh.C

@@ -254,6 +254,60 @@ void Foam::DistributedDelaunayMesh<Triangulation>::findProcessorBoundaryCells
        /Pstream::nProcs()
     );
 
+//    std::list<Cell_handle> infinite_cells;
+//    Triangulation::incident_cells
+//    (
+//        Triangulation::infinite_vertex(),
+//        std::back_inserter(infinite_cells)
+//    );
+//
+//    for
+//    (
+//        typename std::list<Cell_handle>::iterator vcit
+//            = infinite_cells.begin();
+//        vcit != infinite_cells.end();
+//        ++vcit
+//    )
+//    {
+//        Cell_handle cit = *vcit;
+//
+//        // Index of infinite vertex in this cell.
+//        int i = cit->index(Triangulation::infinite_vertex());
+//
+//        Cell_handle c = cit->neighbor(i);
+//
+//        if (c->unassigned())
+//        {
+//            c->cellIndex() = this->getNewCellIndex();
+//
+//            if (checkProcBoundaryCell(c, circumsphereOverlaps))
+//            {
+//                cellToCheck.insert(c->cellIndex());
+//            }
+//        }
+//    }
+//
+//
+//    for
+//    (
+//        Finite_cells_iterator cit = Triangulation::finite_cells_begin();
+//        cit != Triangulation::finite_cells_end();
+//        ++cit
+//    )
+//    {
+//        if (cit->parallelDualVertex())
+//        {
+//            if (cit->unassigned())
+//            {
+//                if (checkProcBoundaryCell(cit, circumsphereOverlaps))
+//                {
+//                    cellToCheck.insert(cit->cellIndex());
+//                }
+//            }
+//        }
+//    }
+
+
     for
     (
         All_cells_iterator cit = Triangulation::all_cells_begin();
@@ -315,12 +369,20 @@ void Foam::DistributedDelaunayMesh<Triangulation>::findProcessorBoundaryCells
                     continue;
                 }
 
-                checkProcBoundaryCell
+                if
                 (
-                    citNeighbor,
+                    checkProcBoundaryCell
-                    circumsphereOverlaps
+                    (
-                );
+                        citNeighbor,
+                        circumsphereOverlaps
+                    )
+                )
+                {
+                    cellToCheck.insert(citNeighbor->cellIndex());
+                }
             }
+
+            cellToCheck.unset(cit->cellIndex());
         }
     }
 }
@@ -521,7 +583,6 @@ Foam::label Foam::DistributedDelaunayMesh<Triangulation>::referVertices
                         << originalParallelVertices[vI].procIndex()
                         << " " << originalParallelVertices[vI].index() << endl;
                 }
-
             }
         }
     }

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C

@@ -194,12 +194,24 @@ void Foam::conformalVoronoiMesh::insertPoints
         }
     }
 
+    label preReinsertionSize(number_of_vertices());
+
     rangeInsertWithInfo
     (
         vertices.begin(),
         vertices.end(),
-        true
+        false
     );
+
+    const label nReinserted = returnReduce
+    (
+        label(number_of_vertices()) - preReinsertionSize,
+        sumOp<label>()
+    );
+
+    Info<< "    Reinserted " << nReinserted << " vertices out of "
+        << returnReduce(vertices.size(), sumOp<label>())
+        << endl;
 }
 
 
@@ -787,8 +799,8 @@ Foam::face Foam::conformalVoronoiMesh::buildDualFace
                 << "Dual face uses circumcenter defined by a "
                 << "Delaunay tetrahedron with no internal "
                 << "or boundary points.  Defining Delaunay edge ends: "
-                << topoint(vA->point()) << " "
+                << vA->info() << " "
-                << topoint(vB->point()) << nl
+                << vB->info() << nl
                 << exit(FatalError);
         }
 
@@ -1034,6 +1046,18 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
         )
     ),
     decomposition_()
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::conformalVoronoiMesh::~conformalVoronoiMesh()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::conformalVoronoiMesh::initialiseForMotion()
 {
     if (foamyHexMeshControls().objOutput())
     {
@@ -1049,7 +1073,10 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
                 runTime_,
                 rndGen_,
                 geometryToConformTo_,
-                foamyHexMeshDict.subDict("backgroundMeshDecomposition")
+                foamyHexMeshControls().foamyHexMeshDict().subDict
+                (
+                    "backgroundMeshDecomposition"
+                )
             )
         );
     }
@@ -1113,13 +1140,53 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+void Foam::conformalVoronoiMesh::initialiseForConformation()
+{
+    if (Pstream::parRun())
+    {
+        decomposition_.reset
+        (
+            new backgroundMeshDecomposition
+            (
+                runTime_,
+                rndGen_,
+                geometryToConformTo_,
+                foamyHexMeshControls().foamyHexMeshDict().subDict
+                (
+                    "backgroundMeshDecomposition"
+                )
+            )
+        );
+    }
 
-Foam::conformalVoronoiMesh::~conformalVoronoiMesh()
+    insertInitialPoints();
-{}
 
+    insertFeaturePoints();
+
+    // Improve the guess that the backgroundMeshDecomposition makes with the
+    // initial positions.  Use before building the surface conformation to
+    // better balance the surface conformation load.
+    distributeBackground(*this);
+
+    buildSurfaceConformation();
+
+    // The introduction of the surface conformation may have distorted the
+    // balance of vertices, distribute if necessary.
+    distributeBackground(*this);
+
+    if (Pstream::parRun())
+    {
+        sync(decomposition_().procBounds());
+    }
+
+    cellSizeMeshOverlapsBackground();
+
+    if (foamyHexMeshControls().printVertexInfo())
+    {
+        printVertexInfo(Info);
+    }
+}
 
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 void Foam::conformalVoronoiMesh::move()
 {
@@ -1630,6 +1697,8 @@ void Foam::conformalVoronoiMesh::move()
             );
         }
 
+        DynamicList<Vertex_handle> pointsToRemove;
+
         for
         (
             Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
@@ -1640,15 +1709,18 @@ void Foam::conformalVoronoiMesh::move()
             if
             (
                 (vit->internalPoint() || vit->internalBoundaryPoint())
-             && !vit->referred()
+             //&& !vit->referred()
             )
             {
-                bool inside = geometryToConformTo_.inside
+                const Foam::point& pt = topoint(vit->point());
-                (
-                    topoint(vit->point())
-                );
 
-                if (!inside)
+                bool inside = geometryToConformTo_.inside(pt);
+
+                if
+                (
+                    !inside
+                 || !geometryToConformTo_.globalBounds().contains(pt)
+                )
                 {
                     if
                     (
@@ -1658,13 +1730,16 @@ void Foam::conformalVoronoiMesh::move()
                     {
                         str().write(topoint(vit->point()));
                     }
-                    remove(vit);
+
+                    pointsToRemove.append(vit);
                     internalPtIsOutside++;
                 }
             }
         }
 
-        Info<< "    " << internalPtIsOutside
+        remove(pointsToRemove.begin(), pointsToRemove.end());
+
+        Info<< "    " << returnReduce(internalPtIsOutside, sumOp<label>())
             << " internal points were inserted outside the domain. "
             << "They have been removed." << endl;
     }

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H

@@ -317,7 +317,7 @@ private:
         void insertPoints
         (
             List<Vb>& vertices,
-            bool distribute = true
+            bool distribute
         );
 
         //- Create a point-pair at a ppDist distance either side of
@@ -764,6 +764,7 @@ private:
             pointIndexHitAndFeatureDynList& featureEdgeHits,
             DynamicList<label>& surfaceToTreeShape,
             DynamicList<label>& edgeToTreeShape,
+            Map<scalar>& surfacePtToEdgePtDist,
             bool firstPass
         ) const;
 
@@ -1037,6 +1038,10 @@ public:
 
     // Member Functions
 
+        void initialiseForMotion();
+
+        void initialiseForConformation();
+
         //- Move the vertices according to the controller, re-conforming to the
         //  surface as required
         void move();

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C

@@ -1787,6 +1787,11 @@ void Foam::conformalVoronoiMesh::indexDualVertices
         }
     }
 
+    OBJstream snapping1("snapToSurface1.obj");
+    OBJstream snapping2("snapToSurface2.obj");
+    OFstream tetToSnapTo("tetsToSnapTo.obj");
+    label offset = 0;
+
     for
     (
         Delaunay::Finite_cells_iterator cit = finite_cells_begin();
@@ -1892,6 +1897,88 @@ void Foam::conformalVoronoiMesh::indexDualVertices
                 }
             }
 
+            {
+                // Snapping points far outside
+                if (cit->boundaryDualVertex() && !cit->parallelDualVertex())
+                {
+                    pointFromPoint dual = cit->dual();
+
+                    pointIndexHit hitInfo;
+                    label surfHit;
+
+                    // Find nearest surface point
+                    geometryToConformTo_.findSurfaceNearest
+                    (
+                        dual,
+                        sqr(targetCellSize(dual)),
+                        hitInfo,
+                        surfHit
+                    );
+
+                    if (!hitInfo.hit())
+                    {
+                        // Project dual to nearest point on tet
+
+                        tetPointRef tet
+                        (
+                            topoint(cit->vertex(0)->point()),
+                            topoint(cit->vertex(1)->point()),
+                            topoint(cit->vertex(2)->point()),
+                            topoint(cit->vertex(3)->point())
+                        );
+
+                        pointFromPoint nearestPointOnTet =
+                            tet.nearestPoint(dual).rawPoint();
+
+                        // Get nearest point on surface from tet.
+                        geometryToConformTo_.findSurfaceNearest
+                        (
+                            nearestPointOnTet,
+                            sqr(targetCellSize(nearestPointOnTet)),
+                            hitInfo,
+                            surfHit
+                        );
+
+                        vector snapDir = nearestPointOnTet - dual;
+                        snapDir /= mag(snapDir) + SMALL;
+
+                        drawDelaunayCell(tetToSnapTo, cit, offset);
+                        offset += 1;
+
+                        vectorField norm(1);
+                        allGeometry_[surfHit].getNormal
+                        (
+                            List<pointIndexHit>(1, hitInfo),
+                            norm
+                        );
+                        norm[0] /= mag(norm[0]) + SMALL;
+
+                        if
+                        (
+                            hitInfo.hit()
+                         && (mag(snapDir & norm[0]) > 0.5)
+                        )
+                        {
+                            snapping1.write
+                            (
+                                linePointRef(dual, nearestPointOnTet)
+                            );
+
+                            snapping2.write
+                            (
+                                linePointRef
+                                (
+                                    nearestPointOnTet,
+                                    hitInfo.hitPoint()
+                                )
+                            );
+
+                            pts[cit->cellIndex()] = hitInfo.hitPoint();
+                        }
+                    }
+                }
+            }
+
             if (cit->boundaryDualVertex())
             {
                 if (cit->featureEdgeDualVertex())
@@ -2458,6 +2545,54 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
          || (vB->internalOrBoundaryPoint() && !vB->referred())
         )
         {
+            if
+            (
+                (vA->internalPoint() && vB->externalBoundaryPoint())
+             || (vB->internalPoint() && vA->externalBoundaryPoint())
+            )
+            {
+                Cell_circulator ccStart = incident_cells(*eit);
+                Cell_circulator cc1 = ccStart;
+                Cell_circulator cc2 = cc1;
+
+                cc2++;
+
+                bool skipEdge = false;
+
+                do
+                {
+                    if
+                    (
+                        cc1->hasFarPoint() || cc2->hasFarPoint()
+                     || is_infinite(cc1) || is_infinite(cc2)
+                    )
+                    {
+                        Pout<< "Ignoring edge between internal and external: "
+                            << vA->info()
+                            << vB->info();
+
+                        skipEdge = true;
+                        break;
+                    }
+
+                    cc1++;
+                    cc2++;
+
+                } while (cc1 != ccStart);
+
+
+                // Do not create faces if the internal point is outside!
+                // This occurs because the internal point is not determined to
+                // be outside in the inside/outside test. This is most likely
+                // due to the triangle.nearestPointClassify test not returning
+                // edge/point as the nearest type.
+
+                if (skipEdge)
+                {
+                    continue;
+                }
+            }
+
             face newDualFace = buildDualFace(eit);
 
             if (newDualFace.size() >= 3)

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C

@@ -298,6 +298,8 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation()
         DynamicList<label> edgeToTreeShape(AtoV/4);
         DynamicList<label> surfaceToTreeShape(AtoV);
 
+        Map<scalar> surfacePtToEdgePtDist(AtoV/4);
+
         for
         (
             Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
@@ -332,6 +334,7 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation()
                         featureEdgeHits,
                         surfaceToTreeShape,
                         edgeToTreeShape,
+                        surfacePtToEdgePtDist,
                         true
                     );
                 }
@@ -463,6 +466,8 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation()
         DynamicList<label> surfaceToTreeShape(AtoV/2);
         DynamicList<label> edgeToTreeShape(AtoV/4);
 
+        Map<scalar> surfacePtToEdgePtDist;
+
         for
         (
             Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
@@ -507,6 +512,7 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation()
                         featureEdgeHits,
                         surfaceToTreeShape,
                         edgeToTreeShape,
+                        surfacePtToEdgePtDist,
                         false
                     );
                 }
@@ -552,6 +558,7 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation()
                         featureEdgeHits,
                         surfaceToTreeShape,
                         edgeToTreeShape,
+                        surfacePtToEdgePtDist,
                         false
                     );
                 }
@@ -569,7 +576,20 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation()
             Vertex_handle vA = c->vertex(eit->second);
             Vertex_handle vB = c->vertex(eit->third);
 
-            if (vA->referred() || vB->referred())
+            if
+            (
+                vA->referred()
+             || vB->referred()
+            )
+            {
+                continue;
+            }
+
+            if
+            (
+                (vA->internalPoint() && vA->referred())
+             || (vB->internalPoint() && vB->referred())
+            )
             {
                 continue;
             }
@@ -617,6 +637,7 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation()
                         featureEdgeHits,
                         surfaceToTreeShape,
                         edgeToTreeShape,
+                        surfacePtToEdgePtDist,
                         false
                     );
                 }
@@ -2076,6 +2097,7 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits
     pointIndexHitAndFeatureDynList& featureEdgeHits,
     DynamicList<label>& surfaceToTreeShape,
     DynamicList<label>& edgeToTreeShape,
+    Map<scalar>& surfacePtToEdgePtDist,
     bool firstPass
 ) const
 {
@@ -2182,48 +2204,35 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits
                             (
                                 pointIndexHitAndFeature(edHit, featureHit)
                             );
+
+//                            Info<< "Add " << existingEdgeLocations_.size() - 1
+//                                << " " << magSqr(edPt - surfPt) << endl;
+
+                            surfacePtToEdgePtDist.insert
+                            (
+                                existingEdgeLocations_.size() - 1,
+                                magSqr(edPt - surfPt)
+                            );
                         }
                         else if (firstPass)
                         {
                             label hitIndex = nearestEdgeHit.index();
 
-                            // Calc new edge location
+//                            Info<< "Close to " << nearestEdgeHit << endl;
-                            Foam::point newPt =
-                                0.5
-                               *(
-                                    nearestEdgeHit.hitPoint()
-                                  + edHit.hitPoint()
-                                );
-
-                            pointIndexHit pHitOld =
-                                edgeLocationTreePtr_().findNearest
-                                (
-                                    nearestEdgeHit.hitPoint(), GREAT
-                                );
-
-                            pointIndexHit pHitNew =
-                                edgeLocationTreePtr_().findNearest
-                                (
-                                    edHit.hitPoint(), GREAT
-                                );
 
                             if
                             (
-                                magSqr(pHitNew.hitPoint() - edHit.hitPoint())
+                                magSqr(edPt - surfPt)
-                              < magSqr
+                              < surfacePtToEdgePtDist[hitIndex]
-                                (
-                                    pHitOld.hitPoint()
-                                  - nearestEdgeHit.hitPoint()
-                                )
                             )
                             {
-                                edHit.setPoint(edHit.hitPoint());
-
                                 featureEdgeHits[hitIndex] =
                                    pointIndexHitAndFeature(edHit, featureHit);
 
                                 existingEdgeLocations_[hitIndex] =
                                     edHit.hitPoint();
+                                surfacePtToEdgePtDist[hitIndex] =
+                                    magSqr(edPt - surfPt);
 
                                 // Change edge location in featureEdgeHits
                                 // remove index from edge tree

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C

@@ -873,17 +873,7 @@ void Foam::conformalVoronoiMesh::reinsertFeaturePoints(bool distribute)
 {
     Info<< nl << "Reinserting stored feature points" << endl;
 
-    label preReinsertionSize(number_of_vertices());
-
     insertPoints(featureVertices_, distribute);
-
-    const label nReinserted = returnReduce
-    (
-        label(number_of_vertices()) - preReinsertionSize,
-        sumOp<label>()
-    );
-
-    Info<< "    Reinserted " << nReinserted << " vertices" << endl;
 }
 
 
@@ -1458,8 +1448,6 @@ void Foam::conformalVoronoiMesh::createFeaturePoints(DynamicList<Vb>& pts)
 
     forAll(feMeshes, i)
     {
-        Info<< indent << "Edge mesh = " << feMeshes[i].name() << nl << endl;
-
         const extendedFeatureEdgeMesh& feMesh(feMeshes[i]);
 
         for

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C

@@ -101,7 +101,9 @@ void Foam::conformalVoronoiMesh::drawDelaunayCell
     // Supply offset as tet number
     offset *= 4;
 
-    os  << "# cell index: " << label(c->cellIndex()) << endl;
+    os  << "# cell index: " << label(c->cellIndex())
+        << " INT_MIN = " << INT_MIN
+        << endl;
 
     os  << "# circumradius "
         << mag(c->dual() - topoint(c->vertex(0)->point()))
@@ -112,7 +114,15 @@ void Foam::conformalVoronoiMesh::drawDelaunayCell
         os  << "# index / type / procIndex: "
             << label(c->vertex(i)->index()) << " "
             << label(c->vertex(i)->type()) << " "
-            << label(c->vertex(i)->procIndex()) << endl;
+            << label(c->vertex(i)->procIndex())
+            << (is_infinite(c->vertex(i)) ? " # This vertex is infinite!" : "")
+            <<
+                (
+                    c->vertex(i)->uninitialised()
+                  ? " # This vertex is uninitialised!"
+                  : ""
+                )
+            << endl;
 
         meshTools::writeOBJ(os, topoint(c->vertex(i)->point()));
     }

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C

@@ -706,8 +706,8 @@ Foam::Field<bool> Foam::conformationSurfaces::wellInside
 //            Info<< surface.name() << " = "
 //                << volumeType::names[surfaceVolumeTests[s][i]] << endl;
 
-            //if (surfaceVolumeTests[s][i] == volumeType::OUTSIDE)
+            if (surfaceVolumeTests[s][i] == volumeType::OUTSIDE)
-            if (surfaceVolumeTests[s][i] != volumeType::INSIDE)
+//            if (surfaceVolumeTests[s][i] != volumeType::INSIDE)
             {
                 insidePoint[i] = false;
 

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/pointFile/pointFile.C

@@ -58,7 +58,7 @@ List<Vb::Point> pointFile::initialPoints() const
         IOobject
         (
             pointFileName_.name(),
-            foamyHexMesh_.time().constant(),
+            foamyHexMesh_.time().timeName(),
             foamyHexMesh_.time(),
             IOobject::MUST_READ,
             IOobject::NO_WRITE

applications/utilities/mesh/generation/foamyHexMesh/foamyHexMesh.C

@@ -45,12 +45,22 @@ int main(int argc, char *argv[])
         "check all surface geometry for quality"
     );
 
+    Foam::argList::addBoolOption
+    (
+        "conformationOnly",
+        "conform to the initial points without any point motion"
+    );
+
+    #include "addOverwriteOption.H"
+
     #include "setRootCase.H"
     #include "createTime.H"
 
     runTime.functionObjects().off();
 
     const bool checkGeometry = args.optionFound("checkGeometry");
+    const bool conformationOnly = args.optionFound("conformationOnly");
+    const bool overwrite = args.optionFound("overwrite");
 
     IOdictionary foamyHexMeshDict
     (
@@ -104,16 +114,33 @@ int main(int argc, char *argv[])
     conformalVoronoiMesh mesh(runTime, foamyHexMeshDict);
 
 
-    while (runTime.loop())
+    if (conformationOnly)
     {
-        Info<< nl << "Time = " << runTime.timeName() << endl;
+        mesh.initialiseForConformation();
 
-        mesh.move();
+        if (!overwrite)
+        {
+            runTime++;
+        }
 
-        Info<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+        mesh.writeMesh(runTime.timeName());
-            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
-            << endl;
     }
+    else
+    {
+        mesh.initialiseForMotion();
+
+        while (runTime.loop())
+        {
+            Info<< nl << "Time = " << runTime.timeName() << endl;
+
+            mesh.move();
+
+            Info<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+                << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+                << endl;
+        }
+    }
+
 
     Info<< nl << "End" << nl << endl;
 

applications/utilities/postProcessing/velocityField/Pe/Pe.C

@@ -130,7 +130,7 @@ void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
                             mesh.magSf()
                           * mesh.surfaceInterpolation::deltaCoeffs()
                           * fvc::interpolate(RASModel->nuEff())
-                         )
+                        )
                     )
                 );
             }
@@ -192,8 +192,12 @@ void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
                             mesh,
                             IOobject::NO_READ
                         ),
-                        mesh.surfaceInterpolation::deltaCoeffs()
+                        mag(phi)
-                      * (mag(phi)/mesh.magSf())*(runTime.deltaT()/nu)
+                       /(
+                            mesh.magSf()
+                          * mesh.surfaceInterpolation::deltaCoeffs()
+                          * nu
+                        )
                     )
                 );
             }
@@ -317,8 +321,12 @@ void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
                             mesh,
                             IOobject::NO_READ
                         ),
-                        mesh.surfaceInterpolation::deltaCoeffs()
+                        mag(phi)
-                      * (mag(phi)/(mesh.magSf()))*(runTime.deltaT()/mu)
+                       /(
+                            mesh.magSf()
+                          * mesh.surfaceInterpolation::deltaCoeffs()
+                          * mu
+                        )
                     )
                 );
             }
@@ -330,26 +338,6 @@ void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
                     << abort(FatalError);
         }
 
-
-        // can also check how many cells exceed a particular Pe limit
-        /*
-        {
-            label count = 0;
-            label PeLimit = 200;
-            forAll(PePtr(), i)
-            {
-                if (PePtr()[i] > PeLimit)
-                {
-                    count++;
-                }
-
-            }
-
-            Info<< "Fraction > " << PeLimit << " = "
-                << scalar(count)/Pe.size() << endl;
-        }
-        */
-
         Info<< "Pe max : " << max(PePtr()).value() << endl;
 
         if (writeResults)

applications/utilities/surface/surfaceCheck/surfaceCheck.C

@@ -653,6 +653,36 @@ int main(int argc, char *argv[])
         Info<< "Checking self-intersection." << endl;
 
         triSurfaceSearch querySurf(surf);
+
+        //{
+        //    OBJstream intStream("selfInter2.obj");
+        //    const indexedOctree<treeDataTriSurface>& tree = querySurf.tree();
+        //    forAll(surf.edges(), edgeI)
+        //    {
+        //        const edge& e = surf.edges()[edgeI];
+        //
+        //        pointIndexHit hitInfo
+        //        (
+        //            tree.findLine
+        //            (
+        //                surf.points()[surf.meshPoints()[e[0]]],
+        //                surf.points()[surf.meshPoints()[e[1]]],
+        //                treeDataTriSurface::findSelfIntersectOp
+        //                (
+        //                    tree,
+        //                    edgeI
+        //                )
+        //            )
+        //        );
+        //
+        //        if (hitInfo.hit())
+        //        {
+        //            Pout<< "Found hit:" << hitInfo.hitPoint() << endl;
+        //            intStream.write(hitInfo.hitPoint());
+        //        }
+        //    }
+        //}
+
         surfaceIntersection inter(querySurf);
 
         if (inter.cutEdges().empty() && inter.cutPoints().empty())

etc/caseDicts/foamyHexMeshDict

@@ -78,6 +78,7 @@ motionControl
 
     objOutput                   no;
     timeChecks                  no;
+    printVertexInfo             off;
 
     maxLoadUnbalance            0.2;
 

etc/codeTemplates/dynamicCode/functionObjectTemplate.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -123,6 +123,19 @@ void ${typeName}FunctionObject::end()
 }
 
 
+void ${typeName}FunctionObject::timeSet()
+{
+    if (${verbose:-false})
+    {
+        Info<<"timeSet ${typeName} sha1: ${SHA1sum}\n";
+    }
+
+//{{{ begin codeTime
+    ${codeTimeSet}
+//}}} end code
+}
+
+
 void ${typeName}FunctionObject::write()
 {
     if (${verbose:-false})

etc/codeTemplates/dynamicCode/functionObjectTemplate.H

@@ -127,6 +127,9 @@ public:
         //- Write, execute the "writeCalls"
         virtual void write();
 
+        //- Called when time was set at the end of the Time::operator++
+        virtual void timeSet();
+
         //- Update for changes of mesh
         virtual void updateMesh(const mapPolyMesh&)
         {}

src/OpenFOAM/db/Time/Time.C

@@ -1220,6 +1220,7 @@ Foam::Time& Foam::Time::operator++()
             writeOnce_ = false;
         }
 
+        functionObjects_.timeSet();
     }
 
     return *this;

src/OpenFOAM/db/functionObjects/OutputFilterFunctionObject/OutputFilterFunctionObject.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -202,6 +202,18 @@ bool Foam::OutputFilterFunctionObject<OutputFilter>::end()
 }
 
 
+template<class OutputFilter>
+bool Foam::OutputFilterFunctionObject<OutputFilter>::timeSet()
+{
+    if (active())
+    {
+        ptr_->timeSet();
+    }
+
+    return true;
+}
+
+
 template<class OutputFilter>
 bool Foam::OutputFilterFunctionObject<OutputFilter>::read
 (

src/OpenFOAM/db/functionObjects/OutputFilterFunctionObject/OutputFilterFunctionObject.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -201,6 +201,8 @@ public:
             //- Called when Time::run() determines that the time-loop exits
             virtual bool end();
 
+            //- Called when time was set at the end of the Time::operator++
+            virtual bool timeSet();
 
             //- Read and set the function object if its data have changed
             virtual bool read(const dictionary&);

src/OpenFOAM/db/functionObjects/functionObject/functionObject.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -120,6 +120,12 @@ bool Foam::functionObject::end()
 }
 
 
+bool Foam::functionObject::timeSet()
+{
+    return false;
+}
+
+
 Foam::autoPtr<Foam::functionObject> Foam::functionObject::iNew::operator()
 (
     const word& name,

src/OpenFOAM/db/functionObjects/functionObject/functionObject.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -157,6 +157,9 @@ public:
         //  By default it simply calls execute().
         virtual bool end();
 
+        //- Called when time was set at the end of the Time::operator++
+        virtual bool timeSet();
+
         //- Read and set the function object if its data have changed
         virtual bool read(const dictionary&) = 0;
 

src/OpenFOAM/db/functionObjects/functionObjectList/functionObjectList.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -190,6 +190,27 @@ bool Foam::functionObjectList::end()
 }
 
 
+bool Foam::functionObjectList::timeSet()
+{
+    bool ok = true;
+
+    if (execution_)
+    {
+        if (!updated_)
+        {
+            read();
+        }
+
+        forAll(*this, objectI)
+        {
+            ok = operator[](objectI).timeSet() && ok;
+        }
+    }
+
+    return ok;
+}
+
+
 bool Foam::functionObjectList::read()
 {
     bool ok = true;

src/OpenFOAM/db/functionObjects/functionObjectList/functionObjectList.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -163,6 +163,9 @@ public:
         //- Called when Time::run() determines that the time-loop exits
         virtual bool end();
 
+        //- Called when time was set at the end of the Time::operator++
+        virtual bool timeSet();
+
         //- Read and set the function objects if their data have changed
         virtual bool read();
 

src/OpenFOAM/fields/Fields/scalarField/scalarField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -27,6 +27,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "scalarField.H"
+#include "unitConversion.H"
 
 #define TEMPLATE
 #include "FieldFunctionsM.C"
@@ -155,6 +156,11 @@ UNARY_FUNCTION(scalar, scalar, j1)
 UNARY_FUNCTION(scalar, scalar, y0)
 UNARY_FUNCTION(scalar, scalar, y1)
 
+UNARY_FUNCTION(scalar, scalar, degToRad)
+UNARY_FUNCTION(scalar, scalar, radToDeg)
+UNARY_FUNCTION(scalar, scalar, atmToPa)
+UNARY_FUNCTION(scalar, scalar, paToAtm)
+
 
 #define BesselFunc(func)                                                      \
 void func(scalarField& res, const int n, const UList<scalar>& sf)             \

src/OpenFOAM/fields/Fields/scalarField/scalarField.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -130,6 +130,10 @@ UNARY_FUNCTION(scalar, scalar, j1)
 UNARY_FUNCTION(scalar, scalar, y0)
 UNARY_FUNCTION(scalar, scalar, y1)
 
+UNARY_FUNCTION(scalar, scalar, degToRad)
+UNARY_FUNCTION(scalar, scalar, radToDeg)
+UNARY_FUNCTION(scalar, scalar, atmToPa)
+UNARY_FUNCTION(scalar, scalar, paToAtm)
 
 #define BesselFunc(func)                                            \
 void func(scalarField& Res, const int n, const UList<scalar>& sf);  \

src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.C

@@ -0,0 +1,73 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 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 "SquareMatrix.H"
+#include "labelList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::scalar Foam::detDecomposed
+(
+    const SquareMatrix<Type>& matrix,
+    const label sign
+)
+{
+    scalar diagProduct = 1.0;
+
+    for (label i = 0; i < matrix.n(); ++i)
+    {
+        diagProduct *= matrix[i][i];
+    }
+
+    return sign*diagProduct;
+}
+
+
+template<class Type>
+Foam::scalar Foam::det(const SquareMatrix<Type>& matrix)
+{
+    SquareMatrix<Type> matrixTmp = matrix;
+
+    labelList pivotIndices(matrix.n());
+    label sign;
+    LUDecompose(matrixTmp, pivotIndices, sign);
+
+    return detDecomposed(matrixTmp, sign);
+}
+
+
+template<class Type>
+Foam::scalar Foam::det(SquareMatrix<Type>& matrix)
+{
+    labelList pivotIndices(matrix.n());
+    label sign;
+    LUDecompose(matrix, pivotIndices, sign);
+
+    return detDecomposed(matrix, sign);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -81,6 +81,21 @@ public:
 };
 
 
+// Global functions
+
+//- Return the LU decomposed SquareMatrix det
+template<class Type>
+scalar detDecomposed(const SquareMatrix<Type>&, const label sign);
+
+//- Return the SquareMatrix det
+template<class Type>
+scalar det(const SquareMatrix<Type>&);
+
+//- Return the SquareMatrix det and the LU decomposition in the original matrix
+template<class Type>
+scalar det(SquareMatrix<Type>&);
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
@@ -91,6 +106,12 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#ifdef NoRepository
+#   include "SquareMatrix.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.C

@@ -535,4 +535,123 @@ const Foam::labelListListList& Foam::GAMGAgglomeration::boundaryFaceMap
 }
 
 
+bool Foam::GAMGAgglomeration::checkRestriction
+(
+    labelList& newRestrict,
+    label& nNewCoarse,
+    const lduAddressing& fineAddressing,
+    const labelUList& restrict,
+    const label nCoarse
+)
+{
+    if (fineAddressing.size() != restrict.size())
+    {
+        FatalErrorIn
+        (
+            "checkRestriction(..)"
+        )   << "nCells:" << fineAddressing.size()
+            << " agglom:" << restrict.size()
+            << abort(FatalError);
+    }
+
+    // Seed (master) for every region
+    labelList master(identity(fineAddressing.size()));
+
+    // Now loop and transport master through region
+    const labelUList& lower = fineAddressing.lowerAddr();
+    const labelUList& upper = fineAddressing.upperAddr();
+
+    while (true)
+    {
+        label nChanged = 0;
+
+        forAll(lower, faceI)
+        {
+            label own = lower[faceI];
+            label nei = upper[faceI];
+
+            if (restrict[own] == restrict[nei])
+            {
+                // coarse-mesh-internal face
+
+                if (master[own] < master[nei])
+                {
+                    master[nei] = master[own];
+                    nChanged++;
+                }
+                else if (master[own] > master[nei])
+                {
+                    master[own] = master[nei];
+                    nChanged++;
+                }
+            }
+        }
+
+        reduce(nChanged, sumOp<label>());
+
+        if (nChanged == 0)
+        {
+            break;
+        }
+    }
+
+
+    // Count number of regions/masters per coarse cell
+    labelListList coarseToMasters(nCoarse);
+    nNewCoarse = 0;
+    forAll(restrict, cellI)
+    {
+        labelList& masters = coarseToMasters[restrict[cellI]];
+
+        if (findIndex(masters, master[cellI]) == -1)
+        {
+            masters.append(master[cellI]);
+            nNewCoarse++;
+        }
+    }
+
+    if (nNewCoarse > nCoarse)
+    {
+        //WarningIn("GAMGAgglomeration::checkRestriction(..)")
+        //    << "Have " << nCoarse
+        //    << " agglomerated cells but " << nNewCoarse
+        //    << " disconnected regions" << endl;
+
+        // Keep coarseToMasters[0] the original coarse, allocate new ones
+        // for the others
+        labelListList coarseToNewCoarse(coarseToMasters.size());
+
+        nNewCoarse = nCoarse;
+
+        forAll(coarseToMasters, coarseI)
+        {
+            const labelList& masters = coarseToMasters[coarseI];
+
+            labelList& newCoarse = coarseToNewCoarse[coarseI];
+            newCoarse.setSize(masters.size());
+            newCoarse[0] = coarseI;
+            for (label i = 1; i < newCoarse.size(); i++)
+            {
+                newCoarse[i] = nNewCoarse++;
+            }
+        }
+
+        newRestrict.setSize(fineAddressing.size());
+        forAll(restrict, cellI)
+        {
+            label coarseI = restrict[cellI];
+
+            label index = findIndex(coarseToMasters[coarseI], master[cellI]);
+            newRestrict[cellI] = coarseToNewCoarse[coarseI][index];
+        }
+
+        return false;
+    }
+    else
+    {
+        return true;
+    }
+}
+
+
 // ************************************************************************* //

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.H

@@ -474,6 +474,16 @@ public:
             const labelListListList& boundaryFaceMap(const label fineLeveli)
             const;
 
+        //- Given restriction determines if coarse cells are connected.
+        //  Return ok is so, otherwise creates new restriction that is
+        static bool checkRestriction
+        (
+            labelList& newRestrict,
+            label& nNewCoarse,
+            const lduAddressing& fineAddressing,
+            const labelUList& restrict,
+            const label nCoarse
+        );
 };
 
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/pairGAMGAgglomeration/pairGAMGAgglomerate.C

@@ -26,177 +26,7 @@ License
 #include "pairGAMGAgglomeration.H"
 #include "lduAddressing.H"
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+// * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
-
-Foam::tmp<Foam::labelField> Foam::pairGAMGAgglomeration::agglomerate
-(
-    label& nCoarseCells,
-    const lduAddressing& fineMatrixAddressing,
-    const scalarField& faceWeights
-)
-{
-    const label nFineCells = fineMatrixAddressing.size();
-
-    const labelUList& upperAddr = fineMatrixAddressing.upperAddr();
-    const labelUList& lowerAddr = fineMatrixAddressing.lowerAddr();
-
-    // For each cell calculate faces
-    labelList cellFaces(upperAddr.size() + lowerAddr.size());
-    labelList cellFaceOffsets(nFineCells + 1);
-
-    // memory management
-    {
-        labelList nNbrs(nFineCells, 0);
-
-        forAll(upperAddr, facei)
-        {
-            nNbrs[upperAddr[facei]]++;
-        }
-
-        forAll(lowerAddr, facei)
-        {
-            nNbrs[lowerAddr[facei]]++;
-        }
-
-        cellFaceOffsets[0] = 0;
-        forAll(nNbrs, celli)
-        {
-            cellFaceOffsets[celli+1] = cellFaceOffsets[celli] + nNbrs[celli];
-        }
-
-        // reset the whole list to use as counter
-        nNbrs = 0;
-
-        forAll(upperAddr, facei)
-        {
-            cellFaces
-            [
-                cellFaceOffsets[upperAddr[facei]] + nNbrs[upperAddr[facei]]
-            ] = facei;
-
-            nNbrs[upperAddr[facei]]++;
-        }
-
-        forAll(lowerAddr, facei)
-        {
-            cellFaces
-            [
-                cellFaceOffsets[lowerAddr[facei]] + nNbrs[lowerAddr[facei]]
-            ] = facei;
-
-            nNbrs[lowerAddr[facei]]++;
-        }
-    }
-
-
-    // go through the faces and create clusters
-
-    tmp<labelField> tcoarseCellMap(new labelField(nFineCells, -1));
-    labelField& coarseCellMap = tcoarseCellMap();
-
-    nCoarseCells = 0;
-
-    for (label celli=0; celli<nFineCells; celli++)
-    {
-        if (coarseCellMap[celli] < 0)
-        {
-            label matchFaceNo = -1;
-            scalar maxFaceWeight = -GREAT;
-
-            // check faces to find ungrouped neighbour with largest face weight
-            for
-            (
-                label faceOs=cellFaceOffsets[celli];
-                faceOs<cellFaceOffsets[celli+1];
-                faceOs++
-            )
-            {
-                label facei = cellFaces[faceOs];
-
-                // I don't know whether the current cell is owner or neighbour.
-                // Therefore I'll check both sides
-                if
-                (
-                    coarseCellMap[upperAddr[facei]] < 0
-                 && coarseCellMap[lowerAddr[facei]] < 0
-                 && faceWeights[facei] > maxFaceWeight
-                )
-                {
-                    // Match found. Pick up all the necessary data
-                    matchFaceNo = facei;
-                    maxFaceWeight = faceWeights[facei];
-                }
-            }
-
-            if (matchFaceNo >= 0)
-            {
-                // Make a new group
-                coarseCellMap[upperAddr[matchFaceNo]] = nCoarseCells;
-                coarseCellMap[lowerAddr[matchFaceNo]] = nCoarseCells;
-                nCoarseCells++;
-            }
-            else
-            {
-                // No match. Find the best neighbouring cluster and
-                // put the cell there
-                label clusterMatchFaceNo = -1;
-                scalar clusterMaxFaceCoeff = -GREAT;
-
-                for
-                (
-                    label faceOs=cellFaceOffsets[celli];
-                    faceOs<cellFaceOffsets[celli+1];
-                    faceOs++
-                )
-                {
-                    label facei = cellFaces[faceOs];
-
-                    if (faceWeights[facei] > clusterMaxFaceCoeff)
-                    {
-                        clusterMatchFaceNo = facei;
-                        clusterMaxFaceCoeff = faceWeights[facei];
-                    }
-                }
-
-                if (clusterMatchFaceNo >= 0)
-                {
-                    // Add the cell to the best cluster
-                    coarseCellMap[celli] = max
-                    (
-                        coarseCellMap[upperAddr[clusterMatchFaceNo]],
-                        coarseCellMap[lowerAddr[clusterMatchFaceNo]]
-                    );
-                }
-            }
-        }
-    }
-
-
-    // Check that all cells are part of clusters,
-    // if not create single-cell "clusters" for each
-    for (label celli=0; celli<nFineCells; celli++)
-    {
-        if (coarseCellMap[celli] < 0)
-        {
-            coarseCellMap[celli] = nCoarseCells;
-            nCoarseCells++;
-        }
-    }
-
-    // Reverse the map ordering to improve the next level of agglomeration
-    // (doesn't always help and is sometimes detrimental)
-    nCoarseCells--;
-
-    forAll(coarseCellMap, celli)
-    {
-        coarseCellMap[celli] = nCoarseCells - coarseCellMap[celli];
-    }
-
-    nCoarseCells++;
-
-    return tcoarseCellMap;
-}
-
 
 void Foam::pairGAMGAgglomeration::agglomerate
 (
@@ -285,4 +115,187 @@ void Foam::pairGAMGAgglomeration::agglomerate
 }
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::labelField> Foam::pairGAMGAgglomeration::agglomerate
+(
+    label& nCoarseCells,
+    const lduAddressing& fineMatrixAddressing,
+    const scalarField& faceWeights
+)
+{
+    const label nFineCells = fineMatrixAddressing.size();
+
+    const labelUList& upperAddr = fineMatrixAddressing.upperAddr();
+    const labelUList& lowerAddr = fineMatrixAddressing.lowerAddr();
+
+    // For each cell calculate faces
+    labelList cellFaces(upperAddr.size() + lowerAddr.size());
+    labelList cellFaceOffsets(nFineCells + 1);
+
+    // memory management
+    {
+        labelList nNbrs(nFineCells, 0);
+
+        forAll(upperAddr, facei)
+        {
+            nNbrs[upperAddr[facei]]++;
+        }
+
+        forAll(lowerAddr, facei)
+        {
+            nNbrs[lowerAddr[facei]]++;
+        }
+
+        cellFaceOffsets[0] = 0;
+        forAll(nNbrs, celli)
+        {
+            cellFaceOffsets[celli+1] = cellFaceOffsets[celli] + nNbrs[celli];
+        }
+
+        // reset the whole list to use as counter
+        nNbrs = 0;
+
+        forAll(upperAddr, facei)
+        {
+            cellFaces
+            [
+                cellFaceOffsets[upperAddr[facei]] + nNbrs[upperAddr[facei]]
+            ] = facei;
+
+            nNbrs[upperAddr[facei]]++;
+        }
+
+        forAll(lowerAddr, facei)
+        {
+            cellFaces
+            [
+                cellFaceOffsets[lowerAddr[facei]] + nNbrs[lowerAddr[facei]]
+            ] = facei;
+
+            nNbrs[lowerAddr[facei]]++;
+        }
+    }
+
+
+    // go through the faces and create clusters
+
+    tmp<labelField> tcoarseCellMap(new labelField(nFineCells, -1));
+    labelField& coarseCellMap = tcoarseCellMap();
+
+    nCoarseCells = 0;
+    label celli;
+
+    for (label cellfi=0; cellfi<nFineCells; cellfi++)
+    {
+        // Change cell ordering depending on direction for this level
+        celli = forward_ ? cellfi : nFineCells - cellfi - 1;
+
+        if (coarseCellMap[celli] < 0)
+        {
+            label matchFaceNo = -1;
+            scalar maxFaceWeight = -GREAT;
+
+            // check faces to find ungrouped neighbour with largest face weight
+            for
+            (
+                label faceOs=cellFaceOffsets[celli];
+                faceOs<cellFaceOffsets[celli+1];
+                faceOs++
+            )
+            {
+                label facei = cellFaces[faceOs];
+
+                // I don't know whether the current cell is owner or neighbour.
+                // Therefore I'll check both sides
+                if
+                (
+                    coarseCellMap[upperAddr[facei]] < 0
+                 && coarseCellMap[lowerAddr[facei]] < 0
+                 && faceWeights[facei] > maxFaceWeight
+                )
+                {
+                    // Match found. Pick up all the necessary data
+                    matchFaceNo = facei;
+                    maxFaceWeight = faceWeights[facei];
+                }
+            }
+
+            if (matchFaceNo >= 0)
+            {
+                // Make a new group
+                coarseCellMap[upperAddr[matchFaceNo]] = nCoarseCells;
+                coarseCellMap[lowerAddr[matchFaceNo]] = nCoarseCells;
+                nCoarseCells++;
+            }
+            else
+            {
+                // No match. Find the best neighbouring cluster and
+                // put the cell there
+                label clusterMatchFaceNo = -1;
+                scalar clusterMaxFaceCoeff = -GREAT;
+
+                for
+                (
+                    label faceOs=cellFaceOffsets[celli];
+                    faceOs<cellFaceOffsets[celli+1];
+                    faceOs++
+                )
+                {
+                    label facei = cellFaces[faceOs];
+
+                    if (faceWeights[facei] > clusterMaxFaceCoeff)
+                    {
+                        clusterMatchFaceNo = facei;
+                        clusterMaxFaceCoeff = faceWeights[facei];
+                    }
+                }
+
+                if (clusterMatchFaceNo >= 0)
+                {
+                    // Add the cell to the best cluster
+                    coarseCellMap[celli] = max
+                    (
+                        coarseCellMap[upperAddr[clusterMatchFaceNo]],
+                        coarseCellMap[lowerAddr[clusterMatchFaceNo]]
+                    );
+                }
+            }
+        }
+    }
+
+    // Check that all cells are part of clusters,
+    // if not create single-cell "clusters" for each
+    for (label cellfi=0; cellfi<nFineCells; cellfi++)
+    {
+        // Change cell ordering depending on direction for this level
+        celli = forward_ ? cellfi : nFineCells - cellfi - 1;
+
+        if (coarseCellMap[celli] < 0)
+        {
+            coarseCellMap[celli] = nCoarseCells;
+            nCoarseCells++;
+        }
+    }
+
+    if (!forward_)
+    {
+        nCoarseCells--;
+
+        forAll(coarseCellMap, celli)
+        {
+            coarseCellMap[celli] = nCoarseCells - coarseCellMap[celli];
+        }
+
+        nCoarseCells++;
+    }
+
+    // Reverse the map ordering for the next level
+    // to improve the next level of agglomeration
+    forward_ = !forward_;
+
+    return tcoarseCellMap;
+}
+
+
 // ************************************************************************* //

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/pairGAMGAgglomeration/pairGAMGAgglomeration.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -29,7 +29,8 @@ License
 
 namespace Foam
 {
-defineTypeNameAndDebug(pairGAMGAgglomeration, 0);
+    defineTypeNameAndDebug(pairGAMGAgglomeration, 0);
+    bool pairGAMGAgglomeration::forward_(true);
 }
 
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/pairGAMGAgglomeration/pairGAMGAgglomeration.H

@@ -56,6 +56,9 @@ class pairGAMGAgglomeration
         //- Number of levels to merge, 1 = don't merge, 2 = merge pairs etc.
         label mergeLevels_;
 
+        //- Direction of cell loop for the current level
+        static bool forward_;
+
 
 protected:
 
@@ -97,7 +100,6 @@ public:
             const lduAddressing& fineMatrixAddressing,
             const scalarField& faceWeights
         );
-
 };
 
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.H

@@ -230,7 +230,8 @@ class GAMGSolver
             const lduMatrix& m,
             const FieldField<Field, scalar>& interfaceBouCoeffs,
             const lduInterfaceFieldPtrsList& interfaces,
-            const scalarField& source,
+            const labelList& restrictAddressing,
+            const scalarField& psiC,
             const direction cmpt
         ) const;
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverInterpolate.C

@@ -34,7 +34,8 @@ void Foam::GAMGSolver::interpolate
     const lduMatrix& m,
     const FieldField<Field, scalar>& interfaceBouCoeffs,
     const lduInterfaceFieldPtrsList& interfaces,
-    const scalarField& source,
+    const labelList& restrictAddressing,
+    const scalarField& psiC,
     const direction cmpt
 ) const
 {
@@ -80,6 +81,30 @@ void Foam::GAMGSolver::interpolate
     {
         psiPtr[celli] = -ApsiPtr[celli]/(diagPtr[celli]);
     }
+
+    register const label nCCells = psiC.size();
+    scalarField corrC(nCCells, 0);
+    scalarField diagC(nCCells, 0);
+
+    for (register label celli=0; celli<nCells; celli++)
+    {
+        corrC[restrictAddressing[celli]] += diagPtr[celli]*psiPtr[celli];
+        diagC[restrictAddressing[celli]] += diagPtr[celli];
+        //corrC[restrictAddressing[celli]] += psiPtr[celli]/diagPtr[celli];
+        //diagC[restrictAddressing[celli]] += 1.0/diagPtr[celli];
+        //corrC[restrictAddressing[celli]] += psiPtr[celli];
+        //diagC[restrictAddressing[celli]] += 1.0;
+    }
+
+    for (register label ccelli=0; ccelli<nCCells; ccelli++)
+    {
+        corrC[ccelli] = psiC[ccelli] - corrC[ccelli]/diagC[ccelli];
+    }
+
+    for (register label celli=0; celli<nCells; celli++)
+    {
+        psiPtr[celli] += corrC[restrictAddressing[celli]];
+    }
 }
 
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C

@@ -313,7 +313,7 @@ void Foam::GAMGSolver::Vcycle
             scalarField& ACfRef =
                 const_cast<scalarField&>(ACf.operator const scalarField&());
 
-            if (interpolateCorrection_)
+            if (interpolateCorrection_ && leveli < coarsestLevel - 2)
             {
                 interpolate
                 (
@@ -322,7 +322,8 @@ void Foam::GAMGSolver::Vcycle
                     matrixLevels_[leveli],
                     interfaceLevelsBouCoeffs_[leveli],
                     interfaceLevels_[leveli],
-                    coarseSources[leveli],
+                    agglomeration_.restrictAddressing(leveli + 1),
+                    coarseCorrFields[leveli + 1],
                     cmpt
                 );
             }
@@ -382,7 +383,8 @@ void Foam::GAMGSolver::Vcycle
             matrix_,
             interfaceBouCoeffs_,
             interfaces_,
-            finestResidual,
+            agglomeration_.restrictAddressing(0),
+            coarseCorrFields[0],
             cmpt
         );
     }

src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C

@@ -33,9 +33,22 @@ void Foam::LUDecompose
     scalarSquareMatrix& matrix,
     labelList& pivotIndices
 )
+{
+    label sign;
+    LUDecompose(matrix, pivotIndices, sign);
+}
+
+
+void Foam::LUDecompose
+(
+    scalarSquareMatrix& matrix,
+    labelList& pivotIndices,
+    label& sign
+)
 {
     label n = matrix.n();
     scalar vv[n];
+    sign = 1;
 
     for (register label i=0; i<n; i++)
     {
@@ -113,6 +126,7 @@ void Foam::LUDecompose
                 Swap(matrixj[k], matrixiMax[k]);
             }
 
+            sign *= -1;
             vv[iMax] = vv[j];
         }
 

src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.H

@@ -79,6 +79,15 @@ void LUDecompose
     labelList& pivotIndices
 );
 
+//- LU decompose the matrix with pivoting.
+//  sign is -1 for odd number of row interchanges and 1 for even number.
+void LUDecompose
+(
+    scalarSquareMatrix& matrix,
+    labelList& pivotIndices,
+    label& sign
+);
+
 //- LU decompose the matrix into a lower (L) and upper (U) part. U = L.T()
 void LUDecompose(scalarSymmetricSquareMatrix& matrix);
 

src/OpenFOAM/memory/autoPtr/autoPtr.H

@@ -63,11 +63,11 @@ public:
         //- Store object pointer
         inline explicit autoPtr(T* = 0);
 
-        //- Construct as copy by transfering pointer to this autoPtr and
+        //- Construct as copy by transferring pointer to this autoPtr and
         //  setting the arguments pointer to NULL
         inline autoPtr(const autoPtr<T>&);
 
-        //- Construct either by transfering pointer or cloning. Should
+        //- Construct either by transferring pointer or cloning. Should
         //  only be called with type that supports cloning.
         inline autoPtr(const autoPtr<T>&, const bool reUse);
 

src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.H

@@ -63,16 +63,6 @@ public:
         const vectorField& cellCtrs
     );
 
-//    static tmp<scalarField> faceOrthogonality(const polyMesh& mesh)
-//    {
-//        return faceOrthogonality
-//        (
-//            mesh,
-//            mesh.faceAreas(),
-//            mesh.cellCentres()
-//        );
-//    }
-
     //- Generate skewness field
     static tmp<scalarField> faceSkewness
     (
@@ -83,18 +73,6 @@ public:
         const vectorField& cellCtrs
     );
 
-//    static tmp<scalarField> faceSkewness(const polyMesh& mesh)
-//    {
-//        return faceSkewness
-//        (
-//            mesh,
-//            mesh.points(),
-//            mesh.faceCentres(),
-//            mesh.faceAreas(),
-//            mesh.cellCentres()
-//        );
-//    }
-
     //- Generate interpolation factors field
     static tmp<scalarField> faceWeights
     (

src/dynamicMesh/polyMeshFilter/polyMeshFilter.C

@@ -550,6 +550,16 @@ Foam::label Foam::polyMeshFilter::filter(const label nOriginalBadFaces)
     // Maintain the number of times a point has been part of a bad face
     labelList pointErrorCount(mesh_.nPoints(), 0);
 
+    PackedBoolList newErrorPoint(mesh_.nPoints());
+    edgeCollapser::checkMeshQuality
+    (
+        mesh_,
+        meshQualityCoeffDict_,
+        newErrorPoint
+    );
+
+    bool newBadFaces = true;
+
     // Main loop
     // ~~~~~~~~~
     // It tries and do some collapses, checks the resulting mesh and
@@ -559,7 +569,8 @@ Foam::label Foam::polyMeshFilter::filter(const label nOriginalBadFaces)
     while
     (
         nOuterIterations < maxIterations_
-     && nBadFaces > nOriginalBadFaces
+     //&& nBadFaces > nOriginalBadFaces
+     && newBadFaces
     )
     {
         Info<< nl << "Outer Iteration = " << nOuterIterations++ << nl
@@ -902,6 +913,21 @@ Foam::label Foam::polyMeshFilter::filter(const label nOriginalBadFaces)
                 isErrorPoint,
                 pointErrorCount
             );
+
+            newBadFaces = false;
+            forAll(mesh_.points(), pI)
+            {
+                if (isErrorPoint[origToCurrentPointMap[pI]])
+                {
+                    if (!newErrorPoint[pI])
+                    {
+                        newBadFaces = true;
+                        break;
+                    }
+                }
+            }
+
+            reduce(newBadFaces, orOp<bool>());
         }
         else
         {

src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMeshTemplates.C

@@ -45,7 +45,7 @@ void Foam::extendedFeatureEdgeMesh::sortPointsAndEdges
     const edgeList& sFeatEds(surf.edges());
     const labelListList edgeFaces = PatchTools::sortedEdgeFaces(surf);
     const vectorField& faceNormals = surf.faceNormals();
-    const labelListList& pointEdges = PatchTools::sortedPointEdges(surf);
+    const labelListList pointEdges = PatchTools::sortedPointEdges(surf);
 
     // Extract and reorder the data from surfaceFeatures
 

src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.C

@@ -88,7 +88,17 @@ Foam::porosityModel::porosityModel
     const word& cellZoneName
 )
 :
-    MeshObject<fvMesh, Foam::UpdateableMeshObject, porosityModel>(mesh),
+    regIOobject
+    (
+        IOobject
+        (
+            name,
+            mesh.time().timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        )
+    ),
     name_(name),
     mesh_(mesh),
     dict_(dict),

src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.H

@@ -36,7 +36,6 @@ SourceFiles
 #ifndef porosityModel_H
 #define porosityModel_H
 
-#include "MeshObject.H"
 #include "fvMesh.H"
 #include "dictionary.H"
 #include "fvMatricesFwd.H"
@@ -56,7 +55,7 @@ namespace Foam
 
 class porosityModel
 :
-    public MeshObject<fvMesh, UpdateableMeshObject, porosityModel>
+    public regIOobject
 {
 private:
 

src/fvAgglomerationMethods/MGridGenGamgAgglomeration/MGridGenGAMGAgglomerate.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-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -30,11 +30,6 @@ Description
 #include "fvMesh.H"
 #include "syncTools.H"
 
-//extern "C"
-//{
-//#   include "mgridgen.h"
-//}
-
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 void Foam::MGridGenGAMGAgglomeration::
@@ -156,6 +151,26 @@ Foam::tmp<Foam::labelField> Foam::MGridGenGAMGAgglomeration::agglomerate
         finalAgglom.begin()
     );
 
+    {
+        label nNewCoarseCells = 0;
+        labelList newRestrictAddr;
+        bool ok = checkRestriction
+        (
+            newRestrictAddr,
+            nNewCoarseCells
+,
+            fineAddressing,
+            finalAgglom,
+            nCoarseCells
+        );
+
+        if (!ok)
+        {
+            nCoarseCells = nNewCoarseCells;
+            finalAgglom.transfer(newRestrictAddr);
+        }
+    }
+
     return tmp<labelField>(new labelField(finalAgglom));
 }
 

src/fvAgglomerationMethods/MGridGenGamgAgglomeration/MGridGenGAMGAgglomeration.C

@@ -59,7 +59,11 @@ Foam::MGridGenGAMGAgglomeration::MGridGenGAMGAgglomeration
 
     // Start geometric agglomeration from the cell volumes and areas of the mesh
     scalarField* VPtr = const_cast<scalarField*>(&fvMesh_.cellVolumes());
-    SubField<vector> Sf(fvMesh_.faceAreas(), fvMesh_.nInternalFaces());
+
+    vectorField magFaceAreas(vector::one*mag(fvMesh_.faceAreas()));
+    SubField<vector> Sf(magFaceAreas, fvMesh_.nInternalFaces());
+    //SubField<vector> Sf(fvMesh_.faceAreas(), fvMesh_.nInternalFaces());
+
     vectorField* SfPtr = const_cast<vectorField*>
     (
         &Sf.operator const vectorField&()

src/fvMotionSolver/fvMotionSolvers/displacement/interpolation/displacementInterpolationMotionSolver.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -36,9 +36,6 @@ Description
 
     Tables are in the \a constant/tables directory.
 
-Note
-    could be a motionSolver - does not use any fvMesh structure.
-
 SourceFiles
     displacementInterpolationMotionSolver.C
 

src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/displacementLayeredMotionMotionSolver.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -253,6 +253,18 @@ Foam::displacementLayeredMotionMotionSolver::faceZoneEvaluate
         // Only on boundary faces - follow boundary conditions
         fld = vectorField(pointDisplacement_, meshPoints);
     }
+    else if (type == "uniformFollow")
+    {
+        // Reads name of name of patch. Then get average point dislacement on
+        // patch. That becomes the value of fld.
+        const word patchName(dict.lookup("patch"));
+        label patchID = mesh().boundaryMesh().findPatchID(patchName);
+        pointField pdf
+        (
+            pointDisplacement_.boundaryField()[patchID].patchInternalField()
+        );
+        fld = gAverage(pdf);
+    }
     else
     {
         FatalIOErrorIn
@@ -399,22 +411,6 @@ Info<< "For cellZone:" << cellZoneI
         // Implement real bc.
         patchDisp[patchI].correctBoundaryConditions();
 
-
-//Info<< "Writing displacement for faceZone " << fz.name()
-//    << " to " << patchDisp[patchI].name() << endl;
-//patchDisp[patchI].write();
-
-//        // Copy into pointDisplacement for other fields to use
-//        forAll(isZonePoint, pointI)
-//        {
-//            if (isZonePoint[pointI])
-//            {
-//                pointDisplacement_[pointI] = patchDisp[patchI][pointI];
-//            }
-//        }
-//        pointDisplacement_.correctBoundaryConditions();
-
-
         patchI++;
     }
 
@@ -423,37 +419,40 @@ Info<< "For cellZone:" << cellZoneI
     // ~~~~~
     // solving the interior is just interpolating
 
-//    // Get normalised distance
+    if (debug)
-//    pointScalarField distance
+    {
-//    (
+        // Get normalised distance
-//        IOobject
+        pointScalarField distance
-//        (
+        (
-//            "distance",
+            IOobject
-//            mesh().time().timeName(),
+            (
-//            mesh(),
+                "distance",
-//            IOobject::NO_READ,
+                mesh().time().timeName(),
-//            IOobject::NO_WRITE,
+                mesh(),
-//            false
+                IOobject::NO_READ,
-//        ),
+                IOobject::NO_WRITE,
-//        pointMesh::New(mesh()),
+                false
-//        dimensionedScalar("distance", dimLength, 0.0)
+            ),
-//    );
+            pointMesh::New(mesh()),
-//    forAll(distance, pointI)
+            dimensionedScalar("distance", dimLength, 0.0)
-//    {
+        );
-//        if (isZonePoint[pointI])
+        forAll(distance, pointI)
-//        {
+        {
-//            scalar d1 = patchDist[0][pointI];
+            if (isZonePoint[pointI])
-//            scalar d2 = patchDist[1][pointI];
+            {
-//            if (d1+d2 > SMALL)
+                scalar d1 = patchDist[0][pointI];
-//            {
+                scalar d2 = patchDist[1][pointI];
-//                scalar s = d1/(d1+d2);
+                if (d1+d2 > SMALL)
-//                distance[pointI] = s;
+                {
-//            }
+                    scalar s = d1/(d1+d2);
-//        }
+                    distance[pointI] = s;
-//    }
+                }
-//    Info<< "Writing distance pointScalarField to " << mesh().time().timeName()
+            }
-//        << endl;
+        }
-//    distance.write();
+        Info<< "Writing distance pointScalarField to "
+            << mesh().time().timeName() << endl;
+        distance.write();
+    }
 
     // Average
     forAll(pointDisplacement_, pointI)
@@ -470,7 +469,6 @@ Info<< "For cellZone:" << cellZoneI
               + s*patchDisp[1][pointI];
         }
     }
-    pointDisplacement_.correctBoundaryConditions();
 }
 
 
@@ -484,9 +482,7 @@ displacementLayeredMotionMotionSolver
 )
 :
     displacementMotionSolver(mesh, dict, typeName)
-{
+{}
-    pointDisplacement_.correctBoundaryConditions();
-}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
@@ -518,6 +514,9 @@ void Foam::displacementLayeredMotionMotionSolver::solve()
     // the motionSolver accordingly
     movePoints(mesh().points());
 
+    // Apply boundary conditions
+    pointDisplacement_.boundaryField().updateCoeffs();
+
     // Apply all regions (=cellZones)
 
     const dictionary& regionDicts = coeffDict().subDict("regions");
@@ -544,6 +543,9 @@ void Foam::displacementLayeredMotionMotionSolver::solve()
 
         cellZoneSolve(zoneI, regionDict);
     }
+
+    // Update pointDisplacement for solved values
+    pointDisplacement_.correctBoundaryConditions();
 }
 
 

src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/displacementLayeredMotionMotionSolver.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -28,11 +28,6 @@ Description
     Mesh motion solver for an (multi-block) extruded fvMesh. Gets given the
     structure of the mesh blocks and boundary conditions on these blocks.
 
-    Note: should not be an fvMotionSolver but just a motionSolver. Only here
-    so we can reuse displacementFvMotionSolver functionality (e.g. surface
-    following boundary conditions)
-
-
     The displacementLayeredMotionCoeffs subdict of dynamicMeshDict specifies
     per region (=cellZone) the boundary conditions on two opposing patches
     (=faceZones). It then interpolates the boundary values using topological
@@ -44,6 +39,9 @@ Description
             Use this for faceZones on boundary faces (so it uses the
             proper boundary conditions on the pointDisplacement).
 
+    uniformFollow: like 'follow' but takes the average value of
+            a specified 'patch' (which is not necessarily colocated)
+
     fixedValue: fixed value.
 
     timeVaryingUniformFixedValue: table-driven fixed value.

src/meshTools/algorithms/PointEdgeWave/PointEdgeWave.C

@@ -415,46 +415,10 @@ void Foam::PointEdgeWave<Type, TrackingData>::handleProcPatches()
         }
     }
 
-
+    // Collocated points should be handled by face based transfer
-
+    // (since that is how connectivity is worked out)
-    //
+    // They are also explicitly equalised in handleCollocatedPoints to
-    // 3. Handle all shared points
+    // guarantee identical values.
-    //    (Note:irrespective if changed or not for now)
-    //
-
-    const globalMeshData& pd = mesh_.globalData();
-
-    List<Type> sharedData(pd.nGlobalPoints());
-
-    forAll(pd.sharedPointLabels(), i)
-    {
-        label meshPointI = pd.sharedPointLabels()[i];
-
-        // Fill my entries in the shared points
-        sharedData[pd.sharedPointAddr()[i]] = allPointInfo_[meshPointI];
-    }
-
-    // Combine on master. Reduce operator has to handle a list and call
-    // Type.updatePoint for all elements
-    combineReduce(sharedData, listUpdateOp<Type>(propagationTol_, td_));
-
-    forAll(pd.sharedPointLabels(), i)
-    {
-        label meshPointI = pd.sharedPointLabels()[i];
-
-        // Retrieve my entries from the shared points.
-        const Type& nbrInfo = sharedData[pd.sharedPointAddr()[i]];
-
-        if (!allPointInfo_[meshPointI].equal(nbrInfo, td_))
-        {
-            updatePoint
-            (
-                meshPointI,
-                nbrInfo,
-                allPointInfo_[meshPointI]
-            );
-        }
-    }
 }
 
 

src/meshTools/indexedOctree/treeDataEdge.C

@@ -109,8 +109,6 @@ Foam::treeDataEdge::findIntersectOp::findIntersectOp
 (
     const indexedOctree<treeDataEdge>& tree
 )
-:
-    tree_(tree)
 {}
 
 

src/meshTools/indexedOctree/treeDataEdge.H

@@ -121,8 +121,6 @@ public:
 
     class findIntersectOp
     {
-        const indexedOctree<treeDataEdge>& tree_;
-
     public:
 
         findIntersectOp(const indexedOctree<treeDataEdge>& tree);

src/meshTools/indexedOctree/treeDataPoint.C

@@ -70,8 +70,6 @@ Foam::treeDataPoint::findIntersectOp::findIntersectOp
 (
     const indexedOctree<treeDataPoint>& tree
 )
-:
-    tree_(tree)
 {}
 
 

src/meshTools/indexedOctree/treeDataPoint.H

@@ -104,8 +104,6 @@ public:
 
     class findIntersectOp
     {
-        const indexedOctree<treeDataPoint>& tree_;
-
     public:
 
         findIntersectOp(const indexedOctree<treeDataPoint>& tree);

src/meshTools/indexedOctree/treeDataPrimitivePatch.C

@@ -66,76 +66,6 @@ void Foam::treeDataPrimitivePatch<PatchType>::update()
 }
 
 
-template<class PatchType>
-bool Foam::treeDataPrimitivePatch<PatchType>::findIntersection
-(
-    const indexedOctree<treeDataPrimitivePatch<PatchType> >& tree,
-    const label index,
-    const point& start,
-    const point& end,
-    point& intersectionPoint
-)
-{
-    const treeDataPrimitivePatch<PatchType>& shape = tree.shapes();
-    const PatchType& patch = shape.patch();
-
-    const pointField& points = patch.points();
-    const typename PatchType::FaceType& f = patch[index];
-
-    // Do quick rejection test
-    if (shape.cacheBb_)
-    {
-        const treeBoundBox& faceBb = shape.bbs_[index];
-
-        if ((faceBb.posBits(start) & faceBb.posBits(end)) != 0)
-        {
-            // start and end in same block outside of faceBb.
-            return false;
-        }
-    }
-
-    const vector dir(end - start);
-    pointHit inter;
-
-    if (f.size() == 3)
-    {
-        inter = triPointRef
-        (
-            points[f[0]],
-            points[f[1]],
-            points[f[2]]
-        ).intersection(start, dir, intersection::HALF_RAY, shape.planarTol_);
-    }
-    else
-    {
-        const pointField& faceCentres = patch.faceCentres();
-
-        inter = f.intersection
-        (
-            start,
-            dir,
-            faceCentres[index],
-            points,
-            intersection::HALF_RAY,
-            shape.planarTol_
-        );
-    }
-
-    if (inter.hit() && inter.distance() <= 1)
-    {
-        // Note: no extra test on whether intersection is in front of us
-        // since using half_ray
-        intersectionPoint = inter.hitPoint();
-
-        return true;
-    }
-    else
-    {
-        return false;
-    }
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct from components
@@ -187,6 +117,19 @@ Foam::treeDataPrimitivePatch<PatchType>::findAllIntersectOp::findAllIntersectOp
 {}
 
 
+template<class PatchType>
+Foam::treeDataPrimitivePatch<PatchType>::
+findSelfIntersectOp::findSelfIntersectOp
+(
+    const indexedOctree<treeDataPrimitivePatch<PatchType> >& tree,
+    const label edgeID
+)
+:
+    tree_(tree),
+    edgeID_(edgeID)
+{}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class PatchType>
@@ -645,4 +588,116 @@ bool Foam::treeDataPrimitivePatch<PatchType>::findAllIntersectOp::operator()
 }
 
 
+template<class PatchType>
+bool Foam::treeDataPrimitivePatch<PatchType>::findSelfIntersectOp::operator()
+(
+    const label index,
+    const point& start,
+    const point& end,
+    point& intersectionPoint
+) const
+{
+    if (edgeID_ == -1)
+    {
+        FatalErrorIn
+        (
+            "findSelfIntersectOp::operator()\n"
+            "(\n"
+            "    const label index,\n"
+            "    const point& start,\n"
+            "    const point& end,\n"
+            "    point& intersectionPoint\n"
+            ") const"
+        )   << "EdgeID not set. Please set edgeID to the index of"
+            << " the edge you are testing"
+            << exit(FatalError);
+    }
+
+    const treeDataPrimitivePatch<PatchType>& shape = tree_.shapes();
+    const PatchType& patch = shape.patch();
+
+    const typename PatchType::FaceType& f = patch.localFaces()[index];
+    const edge& e = patch.edges()[edgeID_];
+
+    if (findIndex(f, e[0]) == -1 && findIndex(f, e[1]) == -1)
+    {
+        return findIntersection(tree_, index, start, end, intersectionPoint);
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+template<class PatchType>
+bool Foam::treeDataPrimitivePatch<PatchType>::findIntersection
+(
+    const indexedOctree<treeDataPrimitivePatch<PatchType> >& tree,
+    const label index,
+    const point& start,
+    const point& end,
+    point& intersectionPoint
+)
+{
+    const treeDataPrimitivePatch<PatchType>& shape = tree.shapes();
+    const PatchType& patch = shape.patch();
+
+    const pointField& points = patch.points();
+    const typename PatchType::FaceType& f = patch[index];
+
+    // Do quick rejection test
+    if (shape.cacheBb_)
+    {
+        const treeBoundBox& faceBb = shape.bbs_[index];
+
+        if ((faceBb.posBits(start) & faceBb.posBits(end)) != 0)
+        {
+            // start and end in same block outside of faceBb.
+            return false;
+        }
+    }
+
+    const vector dir(end - start);
+    pointHit inter;
+
+    if (f.size() == 3)
+    {
+        inter = triPointRef
+        (
+            points[f[0]],
+            points[f[1]],
+            points[f[2]]
+        ).intersection(start, dir, intersection::HALF_RAY, shape.planarTol_);
+    }
+    else
+    {
+        const pointField& faceCentres = patch.faceCentres();
+
+        inter = f.intersection
+        (
+            start,
+            dir,
+            faceCentres[index],
+            points,
+            intersection::HALF_RAY,
+            shape.planarTol_
+        );
+    }
+
+    if (inter.hit() && inter.distance() <= 1)
+    {
+        // Note: no extra test on whether intersection is in front of us
+        // since using half_ray
+        intersectionPoint = inter.hitPoint();
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
 // ************************************************************************* //

src/meshTools/indexedOctree/treeDataPrimitivePatch.H

@@ -86,16 +86,6 @@ class treeDataPrimitivePatch
         //- Initialise all member data
         void update();
 
-        //- Find intersection of line with shapes
-        static bool findIntersection
-        (
-            const indexedOctree<treeDataPrimitivePatch<PatchType> >& tree,
-            const label index,
-            const point& start,
-            const point& end,
-            point& intersectionPoint
-        );
-
 
 public:
 
@@ -141,8 +131,8 @@ public:
 
         findIntersectOp(const indexedOctree<treeDataPrimitivePatch>& tree);
 
-        //- Calculate intersection of triangle with ray. Sets result
+        //- Calculate intersection of any face with ray. Sets result
-        //  accordingly
+        //  accordingly. Used to find first intersection.
         bool operator()
         (
             const label index,
@@ -167,8 +157,34 @@ public:
             DynamicList<label>& shapeMask
         );
 
-        //- Calculate intersection of triangle with ray. Sets result
+        //- Calculate intersection of unique face with ray. Sets result
-        //  accordingly
+        //  accordingly. Used to find all faces.
+        bool operator()
+        (
+            const label index,
+            const point& start,
+            const point& end,
+            point& intersectionPoint
+        ) const;
+    };
+
+
+    class findSelfIntersectOp
+    {
+        const indexedOctree<treeDataPrimitivePatch>& tree_;
+
+        const label edgeID_;
+
+    public:
+
+        findSelfIntersectOp
+        (
+            const indexedOctree<treeDataPrimitivePatch>& tree,
+            const label edgeID
+        );
+
+        //- Calculate intersection of face with edge of patch. Excludes
+        //  faces that use edgeID. Used to find self intersection.
         bool operator()
         (
             const label index,
@@ -234,6 +250,16 @@ public:
                 const point& centre,
                 const scalar radiusSqr
             ) const;
+
+            //- Helper: find intersection of line with shapes
+            static bool findIntersection
+            (
+                const indexedOctree<treeDataPrimitivePatch<PatchType> >& tree,
+                const label index,
+                const point& start,
+                const point& end,
+                point& intersectionPoint
+            );
 };
 
 

src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C

@@ -74,6 +74,26 @@ const Foam::NamedEnum<Foam::distributedTriSurfaceMesh::distributionType, 3>
 // Read my additional data from the dictionary
 bool Foam::distributedTriSurfaceMesh::read()
 {
+
+    if
+    (
+        Pstream::parRun()
+     &&
+     (
+        regIOobject::fileModificationChecking == timeStampMaster
+     || regIOobject::fileModificationChecking == inotifyMaster
+     )
+    )
+    {
+        FatalErrorIn("Foam::distributedTriSurfaceMesh::read()")
+            << "    distributedTriSurfaceMesh is being constructed\n"
+            << "    using 'timeStampMaster' or 'inotifyMaster.'\n"
+            << "    Modify the entry fileModificationChecking\n"
+            << "    in the etc/controlDict.\n"
+            << "    Use 'timeStamp' instead."
+            << exit(FatalError);
+    }
+
     // Get bb of all domains.
     procBb_.setSize(Pstream::nProcs());
 
@@ -183,10 +203,6 @@ void Foam::distributedTriSurfaceMesh::distributeSegment
     List<DynamicList<label> >& sendMap
 ) const
 {
-    // Work points
-    point clipPt;
-
-
     // 1. Fully local already handled outside. Note: retest is cheap.
     if (isLocal(procBb_[Pstream::myProcNo()], start, end))
     {
@@ -227,6 +243,9 @@ void Foam::distributedTriSurfaceMesh::distributeSegment
             // Scheme a: any processor that intersects the segment gets
             // the segment.
 
+            // Intersection point
+            point clipPt;
+
             if (bb.intersects(start, end, clipPt))
             {
                 sendMap[procI].append(allSegments.size());

src/postProcessing/functionObjects/IO/partialWrite/partialWrite.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) 2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -87,13 +87,33 @@ void Foam::partialWrite::read(const dictionary& dict)
             << ". It should be >= 1."
             << exit(FatalIOError);
     }
+
+    // Clear out any previously loaded fields
+    vsf_.clear();
+    vvf_.clear();
+    vSpheretf_.clear();
+    vSymmtf_.clear();
+    vtf_.clear();
+
+    ssf_.clear();
+    svf_.clear();
+    sSpheretf_.clear();
+    sSymmtf_.clear();
+    stf_.clear();
+
+    forAllConstIter(HashSet<word>, objectNames_, iter)
+    {
+        loadField<scalar>(iter.key(), vsf_, ssf_);
+        loadField<vector>(iter.key(), vvf_, svf_);
+        loadField<sphericalTensor>(iter.key(), vSpheretf_, sSpheretf_);
+        loadField<symmTensor>(iter.key(), vSymmtf_, sSymmtf_);
+        loadField<tensor>(iter.key(), vtf_, stf_);
+    }
 }
 
 
 void Foam::partialWrite::execute()
 {
-    //Pout<< "execute at time " << obr_.time().timeName()
-    //    << " index:" << obr_.time().timeIndex() << endl;
 }
 
 
@@ -104,99 +124,61 @@ void Foam::partialWrite::end()
 }
 
 
-void Foam::partialWrite::write()
+void Foam::partialWrite::timeSet()
 {
-    //Pout<< "write at time " << obr_.time().timeName() << endl;
     if (obr_.time().outputTime())
     {
-        // Above check so it can be used both with
-        //  outputControl   timeStep;
-        //  outputInterval  1;
-        // or with
-        //  outputControl   outputTime;
-
         writeInstance_++;
 
         if (writeInstance_ == writeInterval_)
         {
-            // Normal dump
+            // Next overall dump corresponf to partial write. Change
+            // write options to AUTO_WRITE
             writeInstance_ = 0;
+
+            changeWriteOptions<scalar>(vsf_, ssf_, IOobject::AUTO_WRITE);
+            changeWriteOptions<vector>(vvf_, svf_, IOobject::AUTO_WRITE);
+            changeWriteOptions<sphericalTensor>
+            (
+                vSpheretf_,
+                sSpheretf_,
+                IOobject::AUTO_WRITE
+            );
+            changeWriteOptions<symmTensor>
+            (
+                vSymmtf_,
+                sSymmtf_,
+                IOobject::AUTO_WRITE
+            );
+            changeWriteOptions<tensor>(vtf_, stf_, IOobject::AUTO_WRITE);
         }
         else
         {
-            // Delete all but marked objects
+            changeWriteOptions<scalar>(vsf_, ssf_, IOobject::NO_WRITE);
-            fileName dbDir;
+            changeWriteOptions<vector>(vvf_, svf_, IOobject::NO_WRITE);
-            if (isA<polyMesh>(obr_))
+            changeWriteOptions<sphericalTensor>
-            {
-                dbDir = dynamic_cast<const polyMesh&>(obr_).dbDir();
-            }
-
-            IOobjectList objects(obr_, obr_.time().timeName());
-
-            if (debug)
-            {
-                Pout<< "For region:" << obr_.name() << endl;
-            }
-
-            forAllConstIter(IOobjectList, objects, iter)
-            {
-                if (!objectNames_.found(iter()->name()))
-                {
-                    const fileName f =
-                        obr_.time().timePath()
-                       /dbDir
-                       /iter()->name();
-                    if (debug)
-                    {
-                        Pout<< "   rm " << f << endl;
-                    }
-                    rm(f);
-                }
-            }
-
-            // Do the lagrangian files as well.
-            fileNameList cloudDirs
             (
-                readDir
+                vSpheretf_,
-                (
+                sSpheretf_,
-                    obr_.time().timePath()/dbDir/cloud::prefix,
+                IOobject::NO_WRITE
-                    fileName::DIRECTORY
-                )
             );
-            forAll(cloudDirs, i)
+            changeWriteOptions<symmTensor>
-            {
+            (
-                if (debug)
+                vSymmtf_,
-                {
+                sSymmtf_,
-                    Pout<< "For cloud:" << cloudDirs[i] << endl;
+                IOobject::NO_WRITE
-                }
+            );
-
+            changeWriteOptions<tensor>(vtf_, stf_, IOobject::NO_WRITE);
-                IOobjectList sprayObjs
-                (
-                    obr_,
-                    obr_.time().timeName(),
-                    cloud::prefix/cloudDirs[i]
-                );
-                forAllConstIter(IOobjectList, sprayObjs, iter)
-                {
-                    if (!objectNames_.found(iter()->name()))
-                    {
-                        const fileName f =
-                            obr_.time().timePath()
-                           /dbDir
-                           /cloud::prefix
-                           /cloudDirs[i]
-                           /iter()->name();
-                        if (debug)
-                        {
-                            Pout<< "   rm " << f << endl;
-                        }
-                        rm(f);
-                    }
-                }
-            }
         }
     }
 }
 
 
+void Foam::partialWrite::write()
+{
+    // Do nothing. The fields get written through the
+    // standard dump
+}
+
+
 // ************************************************************************* //

src/postProcessing/functionObjects/IO/partialWrite/partialWrite.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -29,8 +29,8 @@ Group
 
 Description
     This function object allows user-selected fields/registered objects to be
-    written at a custom write interval.  It operates by deleting all entries
+    written at a custom write interval. The interval is given in terms of
-    except those selected after writing.
+    number of overall dumps
 
     Example of function object specification:
     \verbatim
@@ -66,8 +66,9 @@ SourceFiles
 #define partialWrite_H
 
 #include "HashSet.H"
-#include "DynamicList.H"
 #include "runTimeSelectionTables.H"
+#include "volFields.H"
+#include "surfaceFields.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -95,6 +96,18 @@ protected:
 
         const objectRegistry& obr_;
 
+        //- Loaded fields
+        UPtrList<volScalarField> vsf_;
+        UPtrList<volVectorField> vvf_;
+        UPtrList<volSphericalTensorField> vSpheretf_;
+        UPtrList<volSymmTensorField> vSymmtf_;
+        UPtrList<volTensorField> vtf_;
+
+        UPtrList<surfaceScalarField> ssf_;
+        UPtrList<surfaceVectorField> svf_;
+        UPtrList<surfaceSphericalTensorField> sSpheretf_;
+        UPtrList<surfaceSymmTensorField> sSymmtf_;
+        UPtrList<surfaceTensorField> stf_;
 
         // Read from dictionary
 
@@ -105,6 +118,7 @@ protected:
             label writeInterval_;
 
 
+
         //- Current dump instance. If reaches writeInterval do a full write.
         label writeInstance_;
 
@@ -118,6 +132,24 @@ protected:
         void operator=(const partialWrite&);
 
 
+        //- Load objects in the objectNames
+        template<class Type>
+        void loadField
+        (
+            const word&,
+            UPtrList<GeometricField<Type, fvPatchField, volMesh> >&,
+            UPtrList<GeometricField<Type, fvsPatchField, surfaceMesh> >&
+        ) const;
+
+        template<class Type>
+        void changeWriteOptions
+        (
+            UPtrList<GeometricField<Type, fvPatchField, volMesh> >&,
+            UPtrList<GeometricField<Type, fvsPatchField, surfaceMesh> >&,
+            const IOobject::writeOption
+        ) const;
+
+
 public:
 
     //- Runtime type information
@@ -152,12 +184,15 @@ public:
         //- Read the partialWrite data
         virtual void read(const dictionary&);
 
-        //- Execute, currently does nothing
+        //- Execute
         virtual void execute();
 
         //- Execute at the final time-loop, currently does nothing
         virtual void end();
 
+        //- Called when time was set at the end of the Time::operator++
+        virtual void timeSet();
+
         //- Write the partialWrite
         virtual void write();
 
@@ -177,6 +212,12 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#ifdef NoRepository
+#   include "partialWriteTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //

src/postProcessing/functionObjects/IO/partialWrite/partialWriteTemplates.C

@@ -0,0 +1,105 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 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 "partialWrite.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+#include "Time.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+void Foam::partialWrite::loadField
+(
+    const word& fieldName,
+    UPtrList<GeometricField<Type, fvPatchField, volMesh> >& vflds,
+    UPtrList<GeometricField<Type, fvsPatchField, surfaceMesh> >& sflds
+) const
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> vfType;
+    typedef GeometricField<Type, fvsPatchField, surfaceMesh> sfType;
+
+    if (obr_.foundObject<vfType>(fieldName))
+    {
+        if (debug)
+        {
+            Info<< "partialWrite : Field "
+                << fieldName << " found in database" << endl;
+        }
+
+        vfType& vField =
+            const_cast<vfType&>
+            (
+                obr_.lookupObject<vfType>(fieldName)
+            );
+
+
+        const unsigned int sz = vflds.size();
+        vflds.setSize(sz + 1);
+        vflds.set(sz, &vField);
+    }
+    else if (obr_.foundObject<sfType>(fieldName))
+    {
+        if (debug)
+        {
+            Info<< "partialWrite : Field " << fieldName
+                << " found in database" << endl;
+        }
+
+         sfType& sField =
+            const_cast<sfType&>
+            (
+                obr_.lookupObject<sfType>(fieldName)
+            );
+
+
+        const unsigned int sz = sflds.size();
+        sflds.setSize(sz + 1);
+        sflds.set(sz, &sField);
+    }
+}
+
+
+template<class Type>
+void Foam::partialWrite::changeWriteOptions
+(
+    UPtrList<GeometricField<Type, fvPatchField, volMesh> >& vflds,
+    UPtrList<GeometricField<Type, fvsPatchField, surfaceMesh> >& sflds,
+    const IOobject::writeOption wOption
+) const
+{
+    forAll(vflds , i)
+    {
+        vflds[i].writeOpt() = wOption;
+    }
+
+    forAll(sflds , i)
+    {
+        sflds[i].writeOpt() = wOption;
+    }
+}
+
+
+// ************************************************************************* //

src/postProcessing/functionObjects/IO/removeRegisteredObject/removeRegisteredObject.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -92,6 +92,12 @@ void Foam::removeRegisteredObject::end()
 }
 
 
+void Foam::removeRegisteredObject::timeSet()
+{
+    // Do nothing - only valid on execute
+}
+
+
 void Foam::removeRegisteredObject::write()
 {
     // Do nothing - only valid on execute

src/postProcessing/functionObjects/IO/removeRegisteredObject/removeRegisteredObject.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -146,6 +146,9 @@ public:
         //- Execute at the final time-loop, currently does nothing
         virtual void end();
 
+        //- Called when time was set at the end of the Time::operator++
+        virtual void timeSet();
+
         //- Write the removeRegisteredObject
         virtual void write();
 

src/postProcessing/functionObjects/IO/writeDictionary/writeDictionary.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -199,6 +199,12 @@ void Foam::writeDictionary::end()
 }
 
 
+void Foam::writeDictionary::timeSet()
+{
+    // do nothing
+}
+
+
 void Foam::writeDictionary::write()
 {
     // do nothing

src/postProcessing/functionObjects/IO/writeDictionary/writeDictionary.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -134,6 +134,9 @@ public:
         //- Execute at the final time-loop, currently does nothing
         virtual void end();
 
+        //- Called when time was set at the end of the Time::operator++
+        virtual void timeSet();
+
         //- Write the writeDictionary
         virtual void write();
 

src/postProcessing/functionObjects/IO/writeRegisteredObject/writeRegisteredObject.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -79,6 +79,12 @@ void Foam::writeRegisteredObject::end()
 }
 
 
+void Foam::writeRegisteredObject::timeSet()
+{
+    // Do nothing - only valid on write
+}
+
+
 void Foam::writeRegisteredObject::write()
 {
     forAll(objectNames_, i)

src/postProcessing/functionObjects/IO/writeRegisteredObject/writeRegisteredObject.H

@@ -147,6 +147,9 @@ public:
         //- Execute at the final time-loop, currently does nothing
         virtual void end();
 
+        //- Called when time was set at the end of the Time::operator++
+        virtual void timeSet();
+
         //- Write the writeRegisteredObject
         virtual void write();
 

src/postProcessing/functionObjects/cloud/cloudInfo/cloudInfo.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -106,6 +106,12 @@ void Foam::cloudInfo::end()
 }
 
 
+void Foam::cloudInfo::timeSet()
+{
+    // Do nothing
+}
+
+
 void Foam::cloudInfo::write()
 {
     if (active_)

src/postProcessing/functionObjects/cloud/cloudInfo/cloudInfo.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -162,6 +162,9 @@ public:
         //- Execute at the final time-loop, currently does nothing
         virtual void end();
 
+        //- Called when time was set at the end of the Time::operator++
+        virtual void timeSet();
+
         //- Write
         virtual void write();
 

src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverage.C

@@ -393,6 +393,10 @@ void Foam::fieldAverage::end()
 {}
 
 
+void Foam::fieldAverage::timeSet()
+{}
+
+
 void Foam::fieldAverage::write()
 {
     if (active_)

src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverage.H

@@ -318,6 +318,9 @@ public:
         //- Execute the averaging at the final time-loop, currently does nothing
         virtual void end();
 
+        //- Called when time was set at the end of the Time::operator++
+        virtual void timeSet();
+
         //- Calculate the field average data and write
         virtual void write();
 

src/postProcessing/functionObjects/field/fieldCoordinateSystemTransform/fieldCoordinateSystemTransform.C

@@ -104,6 +104,12 @@ void Foam::fieldCoordinateSystemTransform::end()
 }
 
 
+void Foam::fieldCoordinateSystemTransform::timeSet()
+{
+    // Do nothing
+}
+
+
 void Foam::fieldCoordinateSystemTransform::write()
 {
     forAll(fieldSet_, fieldI)

src/postProcessing/functionObjects/field/fieldCoordinateSystemTransform/fieldCoordinateSystemTransform.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -173,6 +173,9 @@ public:
         //- Execute at the final time-loop, currently does nothing
         virtual void end();
 
+        //- Called when time was set at the end of the Time::operator++
+        virtual void timeSet();
+
         //- Write
         virtual void write();
 

src/postProcessing/functionObjects/field/fieldMinMax/fieldMinMax.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -139,6 +139,12 @@ void Foam::fieldMinMax::end()
 }
 
 
+void Foam::fieldMinMax::timeSet()
+{
+    // Do nothing - only valid on write
+}
+
+
 void Foam::fieldMinMax::write()
 {
     if (active_)

src/postProcessing/functionObjects/field/fieldMinMax/fieldMinMax.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -189,6 +189,9 @@ public:
         //- Execute at the final time-loop, currently does nothing
         virtual void end();
 
+        //- Called when time was set at the end of the Time::operator++
+        virtual void timeSet();
+
         //- Calculate the field min/max
         template<class Type>
         void calcMinMaxFields

src/postProcessing/functionObjects/field/fieldValues/fieldValue/fieldValue.C

@@ -131,6 +131,12 @@ void Foam::fieldValue::end()
 }
 
 
+void Foam::fieldValue::timeSet()
+{
+    // Do nothing
+}
+
+
 void Foam::fieldValue::updateMesh(const mapPolyMesh&)
 {
     // Do nothing

src/postProcessing/functionObjects/field/fieldValues/fieldValue/fieldValue.H

@@ -191,6 +191,9 @@ public:
             //- Execute the at the final time-loop, currently does nothing
             virtual void end();
 
+            //- Called when time was set at the end of the Time::operator++
+            virtual void timeSet();
+
             //- Update mesh
             virtual void updateMesh(const mapPolyMesh&);
 

src/postProcessing/functionObjects/field/fieldValues/fieldValueDelta/fieldValueDelta.C

@@ -202,6 +202,12 @@ void Foam::fieldValues::fieldValueDelta::end()
 }
 
 
+void Foam::fieldValues::fieldValueDelta::timeSet()
+{
+    // Do nothing
+}
+
+
 void Foam::fieldValues::fieldValueDelta::updateMesh(const mapPolyMesh&)
 {
     // Do nothing

src/postProcessing/functionObjects/field/fieldValues/fieldValueDelta/fieldValueDelta.H

@@ -192,6 +192,9 @@ public:
             //- Execute the at the final time-loop, currently does nothing
             virtual void end();
 
+             //- Called when time was set at the end of the Time::operator++
+            virtual void timeSet();
+
             //- Update mesh
             virtual void updateMesh(const mapPolyMesh&);
 

src/postProcessing/functionObjects/field/nearWallFields/nearWallFields.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -163,6 +163,11 @@ void Foam::nearWallFields::end()
 }
 
 
+void Foam::nearWallFields::timeSet()
+{
+}
+
+
 void Foam::nearWallFields::write()
 {
     if (debug)

src/postProcessing/functionObjects/field/nearWallFields/nearWallFields.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -191,6 +191,9 @@ public:
         //- Execute at the final time-loop, currently does nothing
         virtual void end();
 
+        //- Called when time was set at the end of the Time::operator++
+        virtual void timeSet();
+
         //- Write
         virtual void write();
 

src/postProcessing/functionObjects/field/processorField/processorField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -96,6 +96,12 @@ void Foam::processorField::end()
 }
 
 
+void Foam::processorField::timeSet()
+{
+    // Do nothing
+}
+
+
 void Foam::processorField::write()
 {
     if (active_)

src/postProcessing/functionObjects/field/processorField/processorField.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -145,6 +145,9 @@ public:
         //- Execute at the final time-loop, currently does nothing
         virtual void end();
 
+        //- Called when time was set at the end of the Time::operator++
+        virtual void timeSet();
+
         //- Write
         virtual void write();
 

src/postProcessing/functionObjects/field/readFields/readFields.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -122,6 +122,12 @@ void Foam::readFields::end()
 }
 
 
+void Foam::readFields::timeSet()
+{
+    // Do nothing
+}
+
+
 void Foam::readFields::write()
 {
     // Do nothing

src/postProcessing/functionObjects/field/readFields/readFields.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -173,6 +173,9 @@ public:
         //- Execute at the final time-loop, currently does nothing
         virtual void end();
 
+        //- Called when time was set at the end of the Time::operator++
+        virtual void timeSet();
+
         //- Write
         virtual void write();
 

src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -384,6 +384,12 @@ void Foam::regionSizeDistribution::end()
 }
 
 
+void Foam::regionSizeDistribution::timeSet()
+{
+    // Do nothing - only valid on write
+}
+
+
 void Foam::regionSizeDistribution::write()
 {
     if (active_)

src/postProcessing/functionObjects/field/regionSizeDistribution/regionSizeDistribution.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -266,6 +266,9 @@ public:
         //- Execute at the final time-loop, currently does nothing
         virtual void end();
 
+        //- Called when time was set at the end of the Time::operator++
+        virtual void timeSet();
+
         //- Calculate the regionSizeDistribution and write
         virtual void write();
 

src/postProcessing/functionObjects/field/streamLine/streamLine.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -517,6 +517,10 @@ void Foam::streamLine::end()
 {}
 
 
+void Foam::streamLine::timeSet()
+{}
+
+
 void Foam::streamLine::write()
 {
     if (active_)

src/postProcessing/functionObjects/field/streamLine/streamLine.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -266,6 +266,9 @@ public:
         //- Execute the averaging at the final time-loop, currently does nothing
         virtual void end();
 
+        //- Called when time was set at the end of the Time::operator++
+        virtual void timeSet();
+
         //- Calculate the field average data and write
         virtual void write();
 

src/postProcessing/functionObjects/field/surfaceInterpolateFields/surfaceInterpolateFields.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) 2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -112,6 +112,12 @@ void Foam::surfaceInterpolateFields::end()
 }
 
 
+void Foam::surfaceInterpolateFields::timeSet()
+{
+    // Do nothing
+}
+
+
 void Foam::surfaceInterpolateFields::write()
 {
     if (active_)