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

applications/utilities/parallelProcessing/decomposePar/decomposePar.C

@@ -99,6 +99,40 @@ Usage
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+const labelIOList& procAddressing
+(
+    const PtrList<fvMesh>& procMeshList,
+    const label procI,
+    const word& name,
+    PtrList<labelIOList>& procAddressingList
+)
+{
+    const fvMesh& procMesh = procMeshList[procI];
+
+    if (!procAddressingList.set(procI))
+    {
+        procAddressingList.set
+        (
+            procI,
+            new labelIOList
+            (
+                IOobject
+                (
+                    name,
+                    procMesh.facesInstance(),
+                    procMesh.meshSubDir,
+                    procMesh,
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE
+                )
+            )
+        );
+    }
+    return procAddressingList[procI];
+}
+
+
+
 int main(int argc, char *argv[])
 {
     argList::addNote
@@ -805,74 +839,29 @@ int main(int argc, char *argv[])
                 }
                 const fvMesh& procMesh = procMeshList[procI];
 
+                const labelIOList& faceProcAddressing = procAddressing
+                (
+                    procMeshList,
+                    procI,
+                    "faceProcAddressing",
+                    faceProcAddressingList
+                );
 
-                if (!faceProcAddressingList.set(procI))
+                const labelIOList& cellProcAddressing = procAddressing
-                {
+                (
-                    faceProcAddressingList.set
+                    procMeshList,
-                    (
+                    procI,
-                        procI,
+                    "cellProcAddressing",
-                        new labelIOList
+                    cellProcAddressingList
-                        (
+                );
-                            IOobject
-                            (
-                                "faceProcAddressing",
-                                procMesh.facesInstance(),
-                                procMesh.meshSubDir,
-                                procMesh,
-                                IOobject::MUST_READ,
-                                IOobject::NO_WRITE
-                            )
-                        )
-                    );
-                }
-                const labelIOList& faceProcAddressing =
-                    faceProcAddressingList[procI];
 
-
+                const labelIOList& boundaryProcAddressing = procAddressing
-                if (!cellProcAddressingList.set(procI))
+                (
-                {
+                    procMeshList,
-                    cellProcAddressingList.set
+                    procI,
-                    (
+                    "boundaryProcAddressing",
-                        procI,
+                    boundaryProcAddressingList
-                        new labelIOList
+                );
-                        (
-                            IOobject
-                            (
-                                "cellProcAddressing",
-                                procMesh.facesInstance(),
-                                procMesh.meshSubDir,
-                                procMesh,
-                                IOobject::MUST_READ,
-                                IOobject::NO_WRITE
-                            )
-                        )
-                    );
-                }
-                const labelIOList& cellProcAddressing =
-                    cellProcAddressingList[procI];
-
-
-                if (!boundaryProcAddressingList.set(procI))
-                {
-                    boundaryProcAddressingList.set
-                    (
-                        procI,
-                        new labelIOList
-                        (
-                            IOobject
-                            (
-                                "boundaryProcAddressing",
-                                procMesh.facesInstance(),
-                                procMesh.meshSubDir,
-                                procMesh,
-                                IOobject::MUST_READ,
-                                IOobject::NO_WRITE
-                            )
-                        )
-                    );
-                }
-                const labelIOList& boundaryProcAddressing =
-                    boundaryProcAddressingList[procI];
 
 
                 // FV fields
@@ -959,27 +948,13 @@ int main(int argc, char *argv[])
                  || pointTensorFields.size()
                 )
                 {
-                    if (!pointProcAddressingList.set(procI))
+                    const labelIOList& pointProcAddressing = procAddressing
-                    {
+                    (
-                        pointProcAddressingList.set
+                        procMeshList,
-                        (
+                        procI,
-                            procI,
+                        "pointProcAddressing",
-                            new labelIOList
+                        pointProcAddressingList
-                            (
+                    );
-                                IOobject
-                                (
-                                    "pointProcAddressing",
-                                    procMesh.facesInstance(),
-                                    procMesh.meshSubDir,
-                                    procMesh,
-                                    IOobject::MUST_READ,
-                                    IOobject::NO_WRITE
-                                )
-                            )
-                        );
-                    }
-                    const labelIOList& pointProcAddressing =
-                        pointProcAddressingList[procI];
 
                     const pointMesh& procPMesh = pointMesh::New(procMesh);
 

applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C

@@ -38,6 +38,7 @@ License
 #include "cellSet.H"
 #include "faceSet.H"
 #include "pointSet.H"
+#include "uniformDimensionedFields.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -195,6 +196,60 @@ bool Foam::domainDecomposition::writeDecomposition(const bool decomposeSets)
     }
 
 
+    autoPtr<labelIOList> cellLevelPtr;
+    {
+        IOobject io
+        (
+            "cellLevel",
+            facesInstance(),
+            polyMesh::meshSubDir,
+            *this,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        );
+        if (io.headerOk())
+        {
+            Info<< "Reading hexRef8 data : " << io.name() << endl;
+            cellLevelPtr.reset(new labelIOList(io));
+        }
+    }
+    autoPtr<labelIOList> pointLevelPtr;
+    {
+        IOobject io
+        (
+            "pointLevel",
+            facesInstance(),
+            polyMesh::meshSubDir,
+            *this,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        );
+        if (io.headerOk())
+        {
+            Info<< "Reading hexRef8 data : " << io.name() << endl;
+            pointLevelPtr.reset(new labelIOList(io));
+        }
+    }
+    autoPtr<uniformDimensionedScalarField> level0EdgePtr;
+    {
+        IOobject io
+        (
+            "level0Edge",
+            facesInstance(),
+            polyMesh::meshSubDir,
+            *this,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        );
+        if (io.headerOk())
+        {
+            Info<< "Reading hexRef8 data : " << io.name() << endl;
+            level0EdgePtr.reset(new uniformDimensionedScalarField(io));
+        }
+    }
+
+
+
     label maxProcCells = 0;
     label totProcFaces = 0;
     label maxProcPatches = 0;
@@ -767,8 +822,28 @@ bool Foam::domainDecomposition::writeDecomposition(const bool decomposeSets)
 
         procMesh.write();
 
+        // Write points if pointsInstance differing from facesInstance
+        if (facesInstancePointsPtr_.valid())
+        {
+            pointIOField pointsInstancePoints
+            (
+                IOobject
+                (
+                    "points",
+                    pointsInstance(),
+                    polyMesh::meshSubDir,
+                    procMesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                xferMove(procPoints)
+            );
+            pointsInstancePoints.write();
+        }
 
 
+        // Decompose any sets
         if (decomposeSets)
         {
             forAll(cellSets, i)
@@ -813,25 +888,67 @@ bool Foam::domainDecomposition::writeDecomposition(const bool decomposeSets)
         }
 
 
-        // Write points if pointsInstance differing from facesInstance
+        // hexRef8 data
-        if (facesInstancePointsPtr_.valid())
+        if (cellLevelPtr.valid())
         {
-            pointIOField pointsInstancePoints
+            labelIOList
             (
                 IOobject
                 (
-                    "points",
+                    cellLevelPtr().name(),
-                    pointsInstance(),
+                    facesInstance(),
                     polyMesh::meshSubDir,
                     procMesh,
                     IOobject::NO_READ,
-                    IOobject::NO_WRITE,
+                    IOobject::AUTO_WRITE
-                    false
                 ),
-                xferMove(procPoints)
+                UIndirectList<label>
-            );
+                (
-            pointsInstancePoints.write();
+                    cellLevelPtr(),
+                    procCellAddressing_[procI]
+                )()
+            ).write();
         }
+        if (pointLevelPtr.valid())
+        {
+            labelIOList
+            (
+                IOobject
+                (
+                    pointLevelPtr().name(),
+                    facesInstance(),
+                    polyMesh::meshSubDir,
+                    procMesh,
+                    IOobject::NO_READ,
+                    IOobject::AUTO_WRITE
+                ),
+                UIndirectList<label>
+                (
+                    pointLevelPtr(),
+                    procPointAddressing_[procI]
+                )()
+            ).write();
+        }
+        if (level0EdgePtr.valid())
+        {
+            uniformDimensionedScalarField
+            (
+                IOobject
+                (
+                    level0EdgePtr().name(),
+                    facesInstance(),
+                    polyMesh::meshSubDir,
+                    procMesh,
+                    IOobject::NO_READ,
+                    IOobject::AUTO_WRITE
+                ),
+                level0EdgePtr()
+            ).write();
+        }
+
+
+
+        // Statistics
 
         Info<< endl
             << "Processor " << procI << nl

src/OpenFOAM/fields/pointPatchFields/derived/timeVaryingUniformFixedValue/timeVaryingUniformFixedValuePointPatchField.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
@@ -55,7 +55,9 @@ timeVaryingUniformFixedValuePointPatchField
 :
     fixedValuePointPatchField<Type>(ptf, p, iF, mapper),
     timeSeries_(ptf.timeSeries_)
-{}
+{
+    updateCoeffs();
+}
 
 
 template<class Type>

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

@@ -251,27 +251,35 @@ public:
         {
         public:
 
-            template<class T>
+            template<class Type>
             void operator()
             (
                 const vectorTensorTransform& vt,
                 const bool forward,
-                List<T>& fld
+                List<Type>& fld
             ) const
             {
-                if (forward)
+                const tensor T(forward ? vt.R() : vt.R().T());
+                transformList(T, fld);
+            }
+
+            template<class Type>
+            void operator()
+            (
+                const vectorTensorTransform& vt,
+                const bool forward,
+                List<List<Type> >& flds
+            ) const
+            {
+                forAll(flds, i)
                 {
-                    transformList(vt.R(), fld);
+                    operator()(vt, forward, flds[i]);
-                }
-                else
-                {
-                    transformList(vt.R().T(), fld);
                 }
             }
 
             //- Transform patch-based field
-            template<class T>
+            template<class Type>
-            void operator()(const coupledPolyPatch& cpp, UList<T>& fld) const
+            void operator()(const coupledPolyPatch& cpp, UList<Type>& fld) const
             {
                 if (!cpp.parallel())
                 {
@@ -280,8 +288,8 @@ public:
             }
 
             //- Transform sparse field
-            template<class T, template<class> class Container>
+            template<class Type, template<class> class Container>
-            void operator()(const coupledPolyPatch& cpp, Container<T>& map)
+            void operator()(const coupledPolyPatch& cpp, Container<Type>& map)
             const
             {
                 if (!cpp.parallel())
@@ -313,6 +321,18 @@ public:
                     fld = vt.invTransformPosition(pfld);
                 }
             }
+            void operator()
+            (
+                const vectorTensorTransform& vt,
+                const bool forward,
+                List<List<point> >& flds
+            ) const
+            {
+                forAll(flds, i)
+                {
+                    operator()(vt, forward, flds[i]);
+                }
+            }
             //- Transform patch-based field
             void operator()(const coupledPolyPatch& cpp, pointField& fld) const
             {

src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsNormals.C

@@ -28,43 +28,6 @@ License
 #include "indirectPrimitivePatch.H"
 #include "globalMeshData.H"
 
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    //- Transformation
-    class listTransform
-    {
-    public:
-
-        void operator()
-        (
-            const vectorTensorTransform& vt,
-            const bool forward,
-            List<List<point> >& fld
-        ) const
-        {
-            const tensor T
-            (
-                forward
-              ? vt.R()
-              : vt.R().T()
-            );
-
-            forAll(fld, i)
-            {
-                List<point>& elems = fld[i];
-                forAll(elems, elemI)
-                {
-                    elems[elemI] = transform(T, elems[elemI]);
-                }
-            }
-        }
-    };
-}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template
@@ -141,7 +104,7 @@ Foam::PatchTools::pointNormals
     (
         transforms,
         pointFaceNormals,
-        listTransform()
+        mapDistribute::transform()
     );
 
 

src/dynamicMesh/polyTopoChange/polyTopoChange/removePoints.C

@@ -219,7 +219,7 @@ Foam::label Foam::removePoints::countPointUsage
 
     pointCanBeDeleted.setSize(mesh_.nPoints());
     pointCanBeDeleted = false;
-    label nDeleted = 0;
+    //label nDeleted = 0;
 
     forAll(edge0, pointI)
     {
@@ -243,14 +243,14 @@ Foam::label Foam::removePoints::countPointUsage
             if ((e0Vec & e1Vec) > minCos)
             {
                 pointCanBeDeleted[pointI] = true;
-                nDeleted++;
+                //nDeleted++;
             }
         }
         else if (edge0[pointI] == -1)
         {
             // point not used at all
             pointCanBeDeleted[pointI] = true;
-            nDeleted++;
+            //nDeleted++;
         }
     }
     edge0.clear();
@@ -300,6 +300,15 @@ Foam::label Foam::removePoints::countPointUsage
         true                // null value
     );
 
+    label nDeleted = 0;
+    forAll(pointCanBeDeleted, pointI)
+    {
+        if (pointCanBeDeleted[pointI])
+        {
+            nDeleted++;
+        }
+    }
+
     return returnReduce(nDeleted, sumOp<label>());
 }
 

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

@@ -33,6 +33,7 @@ Description
 #include "fvMesh.H"
 #include "Time.H"
 #include "OFstream.H"
+#include "OBJstream.H"
 #include "mapPolyMesh.H"
 #include "pointEdgePoint.H"
 #include "PointEdgeWave.H"
@@ -1193,12 +1194,29 @@ void Foam::autoSnapDriver::detectNearSurfaces
 
     const pointField avgCc(avgCellCentres(mesh, pp));
 
-    // Construct rays from localPoints to beyond cell centre
+    // Construct rays through localPoints to beyond cell centre
+    pointField start(pp.nPoints());
     pointField end(pp.nPoints());
     forAll(localPoints, pointI)
     {
         const point& pt = localPoints[pointI];
-        end[pointI] = pt + 2*(avgCc[pointI]-pt);
+        const vector d = 2*(avgCc[pointI]-pt);
+        start[pointI] = pt - d;
+        end[pointI] = pt + d;
+    }
+
+
+    autoPtr<OBJstream> gapStr;
+    if (debug&meshRefinement::OBJINTERSECTIONS)
+    {
+        gapStr.reset
+        (
+            new OBJstream
+            (
+                mesh.time().path()
+              / "detectNearSurfaces_" + meshRefiner_.timeName() + ".obj"
+            )
+        );
     }
 
 
@@ -1226,7 +1244,7 @@ void Foam::autoSnapDriver::detectNearSurfaces
         surfaces.findNearestIntersection
         (
             unzonedSurfaces,
-            localPoints,
+            start,
             end,
 
             surface1,
@@ -1248,38 +1266,67 @@ void Foam::autoSnapDriver::detectNearSurfaces
 
             bool override = false;
 
-            if (hit1[pointI].hit())
+            //if (hit1[pointI].hit())
+            //{
+            //    if
+            //    (
+            //        meshRefiner_.isGap
+            //        (
+            //            planarCos,
+            //            nearestPoint[pointI],
+            //            nearestNormal[pointI],
+            //            hit1[pointI].hitPoint(),
+            //            normal1[pointI]
+            //        )
+            //    )
+            //    {
+            //        disp[pointI] = hit1[pointI].hitPoint()-pt;
+            //        override = true;
+            //    }
+            //}
+            //if (hit2[pointI].hit())
+            //{
+            //    if
+            //    (
+            //        meshRefiner_.isGap
+            //        (
+            //            planarCos,
+            //            nearestPoint[pointI],
+            //            nearestNormal[pointI],
+            //            hit2[pointI].hitPoint(),
+            //            normal2[pointI]
+            //        )
+            //    )
+            //    {
+            //        disp[pointI] = hit2[pointI].hitPoint()-pt;
+            //        override = true;
+            //    }
+            //}
+
+            if (hit1[pointI].hit() && hit2[pointI].hit())
             {
                 if
                 (
                     meshRefiner_.isGap
                     (
                         planarCos,
-                        nearestPoint[pointI],
-                        nearestNormal[pointI],
                         hit1[pointI].hitPoint(),
-                        normal1[pointI]
+                        normal1[pointI],
-                    )
-                )
-                {
-                    disp[pointI] = hit1[pointI].hitPoint()-pt;
-                    override = true;
-                }
-            }
-            if (hit2[pointI].hit())
-            {
-                if
-                (
-                    meshRefiner_.isGap
-                    (
-                        planarCos,
-                        nearestPoint[pointI],
-                        nearestNormal[pointI],
                         hit2[pointI].hitPoint(),
                         normal2[pointI]
                     )
                 )
                 {
+                    // TBD: check if the attraction (to nearest) would attract
+                    // good enough and not override attraction
+
+                    if (gapStr.valid())
+                    {
+                        const point& intPt = hit2[pointI].hitPoint();
+                        gapStr().write(linePointRef(pt, intPt));
+                    }
+
+                    // Choose hit2 : nearest to end point (so inside the domain)
                     disp[pointI] = hit2[pointI].hitPoint()-pt;
                     override = true;
                 }
@@ -1337,7 +1384,7 @@ void Foam::autoSnapDriver::detectNearSurfaces
             surfaces.findNearestIntersection
             (
                 surfacesToTest,
-                pointField(localPoints, zonePointIndices),
+                pointField(start, zonePointIndices),
                 pointField(end, zonePointIndices),
 
                 surface1,
@@ -1363,38 +1410,63 @@ void Foam::autoSnapDriver::detectNearSurfaces
 
                 bool override = false;
 
-                if (hit1[i].hit())
+                //if (hit1[i].hit())
+                //{
+                //    if
+                //    (
+                //        meshRefiner_.isGap
+                //        (
+                //            planarCos,
+                //            nearestPoint[pointI],
+                //            nearestNormal[pointI],
+                //            hit1[i].hitPoint(),
+                //            normal1[i]
+                //        )
+                //    )
+                //    {
+                //        disp[pointI] = hit1[i].hitPoint()-pt;
+                //        override = true;
+                //    }
+                //}
+                //if (hit2[i].hit())
+                //{
+                //    if
+                //    (
+                //        meshRefiner_.isGap
+                //        (
+                //            planarCos,
+                //            nearestPoint[pointI],
+                //            nearestNormal[pointI],
+                //            hit2[i].hitPoint(),
+                //            normal2[i]
+                //        )
+                //    )
+                //    {
+                //        disp[pointI] = hit2[i].hitPoint()-pt;
+                //        override = true;
+                //    }
+                //}
+
+                if (hit1[i].hit() && hit2[i].hit())
                 {
                     if
                     (
                         meshRefiner_.isGap
                         (
                             planarCos,
-                            nearestPoint[pointI],
-                            nearestNormal[pointI],
                             hit1[i].hitPoint(),
-                            normal1[i]
+                            normal1[i],
-                        )
-                    )
-                    {
-                        disp[pointI] = hit1[i].hitPoint()-pt;
-                        override = true;
-                    }
-                }
-                if (hit2[i].hit())
-                {
-                    if
-                    (
-                        meshRefiner_.isGap
-                        (
-                            planarCos,
-                            nearestPoint[pointI],
-                            nearestNormal[pointI],
                             hit2[i].hitPoint(),
                             normal2[i]
                         )
                     )
                     {
+                        if (gapStr.valid())
+                        {
+                            const point& intPt = hit2[i].hitPoint();
+                            gapStr().write(linePointRef(pt, intPt));
+                        }
+
                         disp[pointI] = hit2[i].hitPoint()-pt;
                         override = true;
                     }

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C

@@ -623,6 +623,149 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::doRemoveCells
 }
 
 
+// Split faces
+Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitFaces
+(
+    const labelList& splitFaces,
+    const labelPairList& splits
+)
+{
+    polyTopoChange meshMod(mesh_);
+
+    forAll(splitFaces, i)
+    {
+        label faceI = splitFaces[i];
+        const face& f = mesh_.faces()[faceI];
+
+
+        // Split as start and end index in face
+        const labelPair& split = splits[i];
+
+        label nVerts = split[1]-split[0];
+        if (nVerts < 0)
+        {
+            nVerts += f.size();
+        }
+        nVerts += 1;
+
+
+        // Split into f0, f1
+        face f0(nVerts);
+
+        label fp = split[0];
+        forAll(f0, i)
+        {
+            f0[i] = f[fp];
+            fp = f.fcIndex(fp);
+        }
+
+        face f1(f.size()-f0.size()+2);
+        fp = split[1];
+        forAll(f1, i)
+        {
+            f1[i] = f[fp];
+            fp = f.fcIndex(fp);
+        }
+
+
+        // Determine face properties
+        label own = mesh_.faceOwner()[faceI];
+        label nei = -1;
+        label patchI = -1;
+        if (faceI >= mesh_.nInternalFaces())
+        {
+            patchI = mesh_.boundaryMesh().whichPatch(faceI);
+        }
+        else
+        {
+            nei = mesh_.faceNeighbour()[faceI];
+        }
+
+        label zoneI = mesh_.faceZones().whichZone(faceI);
+        bool zoneFlip = false;
+        if (zoneI != -1)
+        {
+            const faceZone& fz = mesh_.faceZones()[zoneI];
+            zoneFlip = fz.flipMap()[fz.whichFace(faceI)];
+        }
+
+
+Pout<< "face:" << faceI << " verts:" << f
+    << " split into f0:" << f0
+    << " f1:" << f1 << endl;
+
+        // Change/add faces
+        meshMod.modifyFace
+        (
+            f0,                         // modified face
+            faceI,                      // label of face
+            own,                        // owner
+            nei,                        // neighbour
+            false,                      // face flip
+            patchI,                     // patch for face
+            zoneI,                      // zone for face
+            zoneFlip                    // face flip in zone
+        );
+        meshMod.addFace
+        (
+            f1,                         // modified face
+            own,                        // owner
+            nei,                        // neighbour
+            -1,                         // master point
+            -1,                         // master edge
+            faceI,                      // master face
+            false,                      // face flip
+            patchI,                     // patch for face
+            zoneI,                      // zone for face
+            zoneFlip                    // face flip in zone
+        );
+    }
+
+
+
+    // Change the mesh (no inflation)
+    autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
+
+    // Update fields
+    mesh_.updateMesh(map);
+
+    // Move mesh (since morphing might not do this)
+    if (map().hasMotionPoints())
+    {
+        mesh_.movePoints(map().preMotionPoints());
+    }
+    else
+    {
+        // Delete mesh volumes. No other way to do this?
+        mesh_.clearOut();
+    }
+
+    // Reset the instance for if in overwrite mode
+    mesh_.setInstance(timeName());
+    setInstance(mesh_.facesInstance());
+
+    // Update local mesh data
+    const labelList& oldToNew = map().reverseFaceMap();
+    labelList newSplitFaces(renumber(oldToNew, splitFaces));
+    // Add added faces (every splitFaces becomes two faces)
+    label sz = newSplitFaces.size();
+    newSplitFaces.setSize(2*sz);
+    forAll(map().faceMap(), faceI)
+    {
+        label oldFaceI = map().faceMap()[faceI];
+        if (oldToNew[oldFaceI] != faceI)
+        {
+            // Added face
+            newSplitFaces[sz++] = faceI;
+        }
+    }
+
+    updateMesh(map, newSplitFaces);
+
+    return map;
+}
+
+
 //// Determine for multi-processor regions the lowest numbered cell on the
 //// lowest numbered processor.
 //void Foam::meshRefinement::getCoupledRegionMaster

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H

@@ -866,6 +866,13 @@ public:
                 const point& keepPoint
             );
 
+            //- Split faces into two
+            autoPtr<mapPolyMesh> splitFaces
+            (
+                const labelList& splitFaces,
+                const labelPairList& splits
+            );
+
             //- Update local numbering for mesh redistribution
             void distribute(const mapDistributePolyMesh&);
 
@@ -1004,6 +1011,26 @@ public:
             //- Do any one of above IO functions. flag is combination of
             //  writeFlag values.
             void write(const label flag, const fileName&) const;
+
+            //- Helper: calculate average
+            template<class T>
+            static T gAverage
+            (
+                const polyMesh& mesh,
+                const PackedBoolList& isMasterElem,
+                const UList<T>& values
+            );
+
+            //- Helper: calculate average over selected elements
+            template<class T>
+            static T gAverage
+            (
+                const polyMesh& mesh,
+                const PackedBoolList& isMasterElem,
+                const labelList& meshPoints,
+                const UList<T>& values
+            );
+
 };
 
 

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C

@@ -267,7 +267,12 @@ void Foam::meshRefinement::markFeatureCellLevel
 
 
     // Find all start cells of features. Is done by tracking from keepPoint.
-    Cloud<trackedParticle> cloud(mesh_, IDLList<trackedParticle>());
+    Cloud<trackedParticle> startPointCloud
+    (
+        mesh_,
+        "startPointCloud",
+        IDLList<trackedParticle>()
+    );
 
 
     // Features are identical on all processors. Number them so we know
@@ -315,7 +320,7 @@ void Foam::meshRefinement::markFeatureCellLevel
                     }
 
                     // Non-manifold point. Create particle.
-                    cloud.addParticle
+                    startPointCloud.addParticle
                     (
                         new trackedParticle
                         (
@@ -359,7 +364,7 @@ void Foam::meshRefinement::markFeatureCellLevel
                     }
 
                     // Non-manifold point. Create particle.
-                    cloud.addParticle
+                    startPointCloud.addParticle
                     (
                         new trackedParticle
                         (
@@ -384,12 +389,13 @@ void Foam::meshRefinement::markFeatureCellLevel
     maxFeatureLevel = labelList(mesh_.nCells(), -1);
 
     // Database to pass into trackedParticle::move
-    trackedParticle::trackingData td(cloud, maxFeatureLevel);
+    trackedParticle::trackingData td(startPointCloud, maxFeatureLevel);
 
     // Track all particles to their end position (= starting feature point)
     // Note that the particle might have started on a different processor
     // so this will transfer across nicely until we can start tracking proper.
-    cloud.move(td, GREAT);
+    startPointCloud.move(td, GREAT);
+
 
     // Reset level
     maxFeatureLevel = -1;
@@ -403,8 +409,64 @@ void Foam::meshRefinement::markFeatureCellLevel
         featureEdgeVisited[featI] = 0u;
     }
 
+
+    Cloud<trackedParticle> cloud
+    (
+        mesh_,
+        "featureCloud",
+        IDLList<trackedParticle>()
+    );
+
+    forAllIter(Cloud<trackedParticle>, startPointCloud, iter)
+    {
+        const trackedParticle& startTp = iter();
+
+        label featI = startTp.i();
+        label pointI = startTp.j();
+
+        const featureEdgeMesh& featureMesh = features_[featI];
+        const labelList& pEdges = featureMesh.pointEdges()[pointI];
+
+        // Now shoot particles down all pEdges.
+        forAll(pEdges, pEdgeI)
+        {
+            label edgeI = pEdges[pEdgeI];
+
+            if (featureEdgeVisited[featI].set(edgeI, 1u))
+            {
+                // Unvisited edge. Make the particle go to the other point
+                // on the edge.
+
+                const edge& e = featureMesh.edges()[edgeI];
+                label otherPointI = e.otherVertex(pointI);
+
+                trackedParticle* tp(new trackedParticle(startTp));
+                tp->end() = featureMesh.points()[otherPointI];
+                tp->j() = otherPointI;
+
+                if (debug&meshRefinement::FEATURESEEDS)
+                {
+                    Pout<< "Adding particle for point:" << pointI
+                        << " coord:" << tp->position()
+                        << " feature:" << featI
+                        << " to track to:" << tp->end()
+                        << endl;
+                }
+
+                cloud.addParticle(tp);
+            }
+        }
+    }
+
+    startPointCloud.clear();
+
+
     while (true)
     {
+        // Track all particles to their end position.
+        cloud.move(td, GREAT);
+
+
         label nParticles = 0;
 
         // Make particle follow edge.
@@ -460,10 +522,23 @@ void Foam::meshRefinement::markFeatureCellLevel
         {
             break;
         }
-
-        // Track all particles to their end position.
-        cloud.move(td, GREAT);
     }
+
+
+
+    //if (debug&meshRefinement::FEATURESEEDS)
+    //{
+    //    forAll(maxFeatureLevel, cellI)
+    //    {
+    //        if (maxFeatureLevel[cellI] != -1)
+    //        {
+    //            Pout<< "Feature went through cell:" << cellI
+    //                << " coord:" << mesh_.cellCentres()[cellI]
+    //                << " leve:" << maxFeatureLevel[cellI]
+    //                << endl;
+    //        }
+    //    }
+    //}
 }
 
 

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementTemplates.C

@@ -57,6 +57,108 @@ template<class T> void meshRefinement::updateList
 }
 
 
+template<class T>
+T meshRefinement::gAverage
+(
+    const polyMesh& mesh,
+    const PackedBoolList& isMasterElem,
+    const UList<T>& values
+)
+{
+    if (values.size() != isMasterElem.size())
+    {
+        FatalErrorIn
+        (
+            "meshRefinement::gAverage\n"
+            "(\n"
+            "    const polyMesh&,\n"
+            "    const PackedBoolList& isMasterElem,\n"
+            "    const UList<T>& values\n"
+            ")\n"
+        )   << "Number of elements in list " << values.size()
+            << " does not correspond to number of elements in isMasterElem "
+            << isMasterElem.size()
+            << exit(FatalError);
+    }
+
+    T sum = pTraits<T>::zero;
+    label n = 0;
+
+    forAll(values, i)
+    {
+        if (isMasterElem[i])
+        {
+            sum += values[i];
+            n++;
+        }
+    }
+
+    reduce(sum, sumOp<T>());
+    reduce(n, sumOp<label>());
+
+    if (n > 0)
+    {
+        return sum/n;
+    }
+    else
+    {
+        return pTraits<T>::max;
+    }
+}
+
+
+template<class T>
+T meshRefinement::gAverage
+(
+    const polyMesh& mesh,
+    const PackedBoolList& isMasterElem,
+    const labelList& meshPoints,
+    const UList<T>& values
+)
+{
+    if (values.size() != meshPoints.size())
+    {
+        FatalErrorIn
+        (
+            "meshRefinement::gAverage\n"
+            "(\n"
+            "    const polyMesh&,\n"
+            "    const labelList&,\n"
+            "    const PackedBoolList& isMasterElem,\n"
+            "    const UList<T>& values\n"
+            ")\n"
+        )   << "Number of elements in list " << values.size()
+            << " does not correspond to number of elements in meshPoints "
+            << meshPoints.size()
+            << exit(FatalError);
+    }
+
+    T sum = pTraits<T>::zero;
+    label n = 0;
+
+    forAll(values, i)
+    {
+        if (isMasterElem[meshPoints[i]])
+        {
+            sum += values[i];
+            n++;
+        }
+    }
+
+    reduce(sum, sumOp<T>());
+    reduce(n, sumOp<label>());
+
+    if (n > 0)
+    {
+        return sum/n;
+    }
+    else
+    {
+        return pTraits<T>::max;
+    }
+}
+
+
 // Compare two lists over all boundary faces
 template<class T>
 void meshRefinement::testSyncBoundaryFaceList

src/mesh/autoMesh/autoHexMesh/trackedParticle/trackedParticle.C

@@ -169,7 +169,7 @@ void Foam::trackedParticle::hitProcessorPatch
     trackingData& td
 )
 {
-    // Remove particle
+    // Move to different processor
     td.switchProcessor = true;
 }
 

src/mesh/autoMesh/autoHexMesh/trackedParticle/trackedParticle.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
@@ -163,12 +163,24 @@ public:
             return end_;
         }
 
+        //- transported label
+        label i() const
+        {
+            return i_;
+        }
+
         //- transported label
         label& i()
         {
             return i_;
         }
 
+        //- transported label
+        label j() const
+        {
+            return j_;
+        }
+
         //- transported label
         label& j()
         {

src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C

@@ -135,145 +135,12 @@ Foam::sampledTriSurfaceMesh::nonCoupledboundaryTree() const
 }
 
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+bool Foam::sampledTriSurfaceMesh::update(const meshSearch& meshSearcher)
-
-Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
-(
-    const word& name,
-    const polyMesh& mesh,
-    const word& surfaceName,
-    const samplingSource sampleSource
-)
-:
-    sampledSurface(name, mesh),
-    surface_
-    (
-        IOobject
-        (
-            surfaceName,
-            mesh.time().constant(), // instance
-            "triSurface",           // local
-            mesh,                   // registry
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE,
-            false
-        )
-    ),
-    sampleSource_(sampleSource),
-    needsUpdate_(true),
-    sampleElements_(0),
-    samplePoints_(0)
-{}
-
-
-Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
-(
-    const word& name,
-    const polyMesh& mesh,
-    const dictionary& dict
-)
-:
-    sampledSurface(name, mesh, dict),
-    surface_
-    (
-        IOobject
-        (
-            dict.lookup("surface"),
-            mesh.time().constant(), // instance
-            "triSurface",           // local
-            mesh,                   // registry
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE,
-            false
-        )
-    ),
-    sampleSource_(samplingSourceNames_[dict.lookup("source")]),
-    needsUpdate_(true),
-    sampleElements_(0),
-    samplePoints_(0)
-{}
-
-
-Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
-(
-    const word& name,
-    const polyMesh& mesh,
-    const triSurface& surface,
-    const word& sampleSourceName
-)
-:
-    sampledSurface(name, mesh),
-    surface_
-    (
-        IOobject
-        (
-            name,
-            mesh.time().constant(), // instance
-            "triSurface",           // local
-            mesh,                  // registry
-            IOobject::NO_READ,
-            IOobject::NO_WRITE,
-            false
-        ),
-        surface
-    ),
-    sampleSource_(samplingSourceNames_[sampleSourceName]),
-    needsUpdate_(true),
-    sampleElements_(0),
-    samplePoints_(0)
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::sampledTriSurfaceMesh::~sampledTriSurfaceMesh()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-bool Foam::sampledTriSurfaceMesh::needsUpdate() const
 {
-    return needsUpdate_;
-}
-
-
-bool Foam::sampledTriSurfaceMesh::expire()
-{
-    // already marked as expired
-    if (needsUpdate_)
-    {
-        return false;
-    }
-
-    sampledSurface::clearGeom();
-    MeshStorage::clear();
-
-    boundaryTreePtr_.clear();
-    sampleElements_.clear();
-    samplePoints_.clear();
-
-    needsUpdate_ = true;
-    return true;
-}
-
-
-bool Foam::sampledTriSurfaceMesh::update()
-{
-    if (!needsUpdate_)
-    {
-        return false;
-    }
-
-
     // Find the cells the triangles of the surface are in.
     // Does approximation by looking at the face centres only
     const pointField& fc = surface_.faceCentres();
 
-    // Mesh search engine, no triangulation of faces.
-    meshSearch meshSearcher(mesh(), polyMesh::FACEPLANES);
-
-
     List<nearInfo> nearest(fc.size());
 
     // Global numbering for cells/faces - only used to uniquely identify local
@@ -611,6 +478,157 @@ bool Foam::sampledTriSurfaceMesh::update()
 }
 
 
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
+(
+    const word& name,
+    const polyMesh& mesh,
+    const word& surfaceName,
+    const samplingSource sampleSource
+)
+:
+    sampledSurface(name, mesh),
+    surface_
+    (
+        IOobject
+        (
+            surfaceName,
+            mesh.time().constant(), // instance
+            "triSurface",           // local
+            mesh,                   // registry
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE,
+            false
+        )
+    ),
+    sampleSource_(sampleSource),
+    needsUpdate_(true),
+    sampleElements_(0),
+    samplePoints_(0)
+{}
+
+
+Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
+(
+    const word& name,
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+:
+    sampledSurface(name, mesh, dict),
+    surface_
+    (
+        IOobject
+        (
+            dict.lookup("surface"),
+            mesh.time().constant(), // instance
+            "triSurface",           // local
+            mesh,                   // registry
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE,
+            false
+        )
+    ),
+    sampleSource_(samplingSourceNames_[dict.lookup("source")]),
+    needsUpdate_(true),
+    sampleElements_(0),
+    samplePoints_(0)
+{}
+
+
+Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
+(
+    const word& name,
+    const polyMesh& mesh,
+    const triSurface& surface,
+    const word& sampleSourceName
+)
+:
+    sampledSurface(name, mesh),
+    surface_
+    (
+        IOobject
+        (
+            name,
+            mesh.time().constant(), // instance
+            "triSurface",           // local
+            mesh,                  // registry
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        surface
+    ),
+    sampleSource_(samplingSourceNames_[sampleSourceName]),
+    needsUpdate_(true),
+    sampleElements_(0),
+    samplePoints_(0)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::sampledTriSurfaceMesh::~sampledTriSurfaceMesh()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::sampledTriSurfaceMesh::needsUpdate() const
+{
+    return needsUpdate_;
+}
+
+
+bool Foam::sampledTriSurfaceMesh::expire()
+{
+    // already marked as expired
+    if (needsUpdate_)
+    {
+        return false;
+    }
+
+    sampledSurface::clearGeom();
+    MeshStorage::clear();
+
+    boundaryTreePtr_.clear();
+    sampleElements_.clear();
+    samplePoints_.clear();
+
+    needsUpdate_ = true;
+    return true;
+}
+
+
+bool Foam::sampledTriSurfaceMesh::update()
+{
+    if (!needsUpdate_)
+    {
+        return false;
+    }
+
+    // Mesh search engine, no triangulation of faces.
+    meshSearch meshSearcher(mesh(), polyMesh::FACEPLANES);
+
+    return update(meshSearcher);
+}
+
+
+bool Foam::sampledTriSurfaceMesh::update(const treeBoundBox& bb)
+{
+    if (!needsUpdate_)
+    {
+        return false;
+    }
+
+    // Mesh search engine on subset, no triangulation of faces.
+    meshSearch meshSearcher(mesh(), bb, polyMesh::FACEPLANES);
+
+    return update(meshSearcher);
+}
+
+
 Foam::tmp<Foam::scalarField> Foam::sampledTriSurfaceMesh::sample
 (
     const volScalarField& vField

src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.H

@@ -72,6 +72,7 @@ namespace Foam
 {
 
 class treeDataFace;
+class meshSearch;
 
 /*---------------------------------------------------------------------------*\
                        Class sampledTriSurfaceMesh Declaration
@@ -136,6 +137,8 @@ private:
         tmp<Field<Type> >
         interpolateField(const interpolation<Type>&) const;
 
+        bool update(const meshSearch& meshSearcher);
+
 public:
 
     //- Runtime type information
@@ -189,6 +192,10 @@ public:
         //  Do nothing (and return false) if no update was needed
         virtual bool update();
 
+        //- Update the surface using a bound box to limit the searching.
+        //  For direct use, i.e. not through sample.
+        //  Do nothing (and return false) if no update was needed
+        bool update(const treeBoundBox&);
 
         //- Points of surface
         virtual const pointField& points() const