Initial autoMesh merge

applications/utilities/mesh/conversion/plot3dToFoam/hexBlock.C

@@ -87,7 +87,7 @@ void hexBlock::setHandedness()
 
     if (blockHandedness_ == noPoints)
     {
-        WarningIn("hexBlock::readPoints(const bool, Istream&)")
+        WarningIn("hexBlock::hexBlock::setHandedness()")
             << "Cannot determine orientation of block."
             << " Continuing as if right handed." << endl;
         blockHandedness_ = right;
@@ -105,35 +105,67 @@ hexBlock::hexBlock(const label nx, const label ny, const label nz)
     zDim_(nz - 1),
     blockHandedness_(noPoints),
     points_((xDim_ + 1)*(yDim_ + 1)*(zDim_ + 1))
-{
-    Pout<< "xDim:" << nx << " yDim:" << ny << " zDim:" << nz << endl;
-}
+{}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void hexBlock::readPoints(const bool readBlank, Istream& is)
+void hexBlock::readPoints
+(
+    const bool readBlank,
+    const scalar twoDThicknes,
+    Istream& is
+)
 {
     scalar iBlank;
 
-    forAll (points_, i)
+    label nPoints = points_.size();
+
+    if (twoDThicknes > 0)
+    {
+        nPoints /= 2;
+    }
+
+    Info<< "Reading " << nPoints << " x coordinates..." << endl;
+    for (label i=0; i < nPoints; i++)
     {
         is >> points_[i].x();
     }
 
-    forAll (points_, i)
+    Info<< "Reading " << nPoints << " y coordinates..." << endl;
+    for (label i=0; i < nPoints; i++)
     {
         is >> points_[i].y();
     }
 
-    forAll (points_, i)
+    if (twoDThicknes > 0)
     {
-        is >> points_[i].z();
+        Info<< "Extruding " << nPoints << " points in z direction..." << endl;
+        // Duplicate points
+        for (label i=0; i < nPoints; i++)
+        {
+            points_[i+nPoints] = points_[i];
+        }
+        for (label i=0; i < nPoints; i++)
+        {
+            points_[i].z() = 0;
+            points_[i+nPoints].z() = twoDThicknes;
+        }
     }
+    else
+    {
+        Info<< "Reading " << nPoints << " z coordinates..." << endl;
+        for (label i=0; i < nPoints; i++)
+        {
+            is >> points_[i].z();
+        }
+    }
+
 
     if (readBlank)
     {
-        forAll (points_, i)
+        Info<< "Reading " << nPoints << " blanks..." << endl;
+        for (label i=0; i < nPoints; i++)
         {
             is >> iBlank;
         }

applications/utilities/mesh/conversion/plot3dToFoam/hexBlock.H

@@ -145,7 +145,14 @@ public:
 
 
         //- Read block points either with or without blanking after every block.
-        void readPoints(const bool readBlank, Istream&);
+        //  If twoDThickness > 0 reads (half) the points and extrudes the
+        //  points in z direction.
+        void readPoints
+        (
+            const bool readBlank,
+            const scalar twoDThicknes,
+            Istream&
+        );
 };
 
 

applications/utilities/mesh/conversion/plot3dToFoam/plot3dToFoam.C

@@ -23,11 +23,12 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Description
-    Plot3d mesh (ascii format) converter.
+    Plot3d mesh (ascii/formatted format) converter.
 
     Work in progress! Handles ascii multiblock (and optionally singleBlock)
     format.
     By default expects blanking. Use -noBlank if none.
+    Use -2D <thickness> if 2D.
     Niklas Nordin has experienced a problem with lefthandedness of the blocks.
     The code should detect this automatically - see hexBlock::readPoints but
     if this goes wrong just set the blockHandedness_ variable to 'right'
@@ -59,6 +60,7 @@ int main(int argc, char *argv[])
     argList::validOptions.insert("scale", "scale factor");
     argList::validOptions.insert("noBlank", "");
     argList::validOptions.insert("singleBlock", "");
+    argList::validOptions.insert("2D", "thickness");
 
     argList args(argc, argv);
 
@@ -75,6 +77,13 @@ int main(int argc, char *argv[])
 
     bool readBlank = !args.options().found("noBlank");
     bool singleBlock = args.options().found("singleBlock");
+    scalar twoDThicknes = -1;
+    if (args.options().found("2D"))
+    {
+        twoDThicknes = readScalar(IStringStream(args.options()["2D"])());
+        Info<< "Reading 2D case by extruding points by " << twoDThicknes
+            << " in z direction." << nl << endl;
+    }
 
 
 #   include "createTime.H"
@@ -95,7 +104,7 @@ int main(int argc, char *argv[])
         plot3dFile >> nblock;
     }
 
-    Info << "Reading " << nblock << " blocks" << endl;
+    Info<< "Reading " << nblock << " blocks" << endl;
 
     PtrList<hexBlock> blocks(nblock);
 
@@ -104,20 +113,32 @@ int main(int argc, char *argv[])
 
         forAll (blocks, blockI)
         {
-            plot3dFile >> nx >> ny >> nz;
+            if (twoDThicknes > 0)
+            {
+                // Fake second set of points (done in readPoints below)
+                plot3dFile >> nx >> ny;
+                nz = 2;
+            }
+            else
+            {
+                plot3dFile >> nx >> ny >> nz;
+            }
+
+            Info<< "block " << blockI << " nx:" << nx
+                << " ny:" << ny << " nz:" << nz << endl;
 
             blocks.set(blockI, new hexBlock(nx, ny, nz));
         }
     }
 
-    Info << "Reading block points" << endl;
+    Info<< "Reading block points" << endl;
     label sumPoints(0);
     label nMeshCells(0);
 
     forAll (blocks, blockI)
     {
-        Info << "block " << blockI << ":" << nl;
-        blocks[blockI].readPoints(readBlank, plot3dFile);
+        Info<< "block " << blockI << ":" << nl;
+        blocks[blockI].readPoints(readBlank, twoDThicknes, plot3dFile);
         sumPoints += blocks[blockI].nBlockPoints();
         nMeshCells += blocks[blockI].nBlockCells();
         Info<< nl;
@@ -136,7 +157,6 @@ int main(int argc, char *argv[])
         }
     }
 
-
     // From old to new master point
     labelList oldToNew;
     pointField newPoints;
@@ -151,13 +171,17 @@ int main(int argc, char *argv[])
         newPoints
     );
 
+    Info<< "Merged points within " << SMALL << " distance. Merged from "
+        << oldToNew.size() << " down to " << newPoints.size()
+        << " points." << endl;
+
     // Scale the points
     if (scaleFactor > 1.0 + SMALL || scaleFactor < 1.0 - SMALL)
     {
         newPoints *= scaleFactor;
     }
 
-    Info << "Creating cells" << endl;
+    Info<< "Creating cells" << endl;
 
     cellShapeList cellShapes(nMeshCells);
 
@@ -190,7 +214,7 @@ int main(int argc, char *argv[])
         }
     }
 
-    Info << "Creating boundary patches" << endl;
+    Info<< "Creating boundary patches" << endl;
 
     faceListList boundary(0);
     wordList patchNames(0);
@@ -220,10 +244,10 @@ int main(int argc, char *argv[])
     // Set the precision of the points data to 10
     IOstream::defaultPrecision(10);
 
-    Info << "Writing polyMesh" << endl;
+    Info<< "Writing polyMesh" << endl;
     pShapeMesh.write();
 
-    Info << "End\n" << endl;
+    Info<< "End\n" << endl;
 
     return 0;
 }

applications/utilities/mesh/manipulation/Optional/setSet/setSet.C

@@ -42,6 +42,7 @@ Description
 #include "demandDrivenData.H"
 #include "writePatch.H"
 #include "writePointSet.H"
+#include "IOobjectList.H"
 
 #include <stdio.h>
 
@@ -233,8 +234,10 @@ void writeVTK
 
 void printHelp(Ostream& os)
 {
-    os  << "Please type 'help', 'quit', 'time ddd'"
+    os  << "Please type 'help', 'list', 'quit', 'time ddd'"
         << " or a set command after prompt." << endl
+        << "'list' will show all current cell/face/point sets." << endl
+        << "'time ddd' will change the current time." << endl
         << endl
         << "A set command should be of the following form" << endl
         << endl
@@ -272,6 +275,47 @@ void printHelp(Ostream& os)
 }
 
 
+void printAllSets(const polyMesh& mesh, Ostream& os)
+{
+    IOobjectList objects
+    (
+        mesh,
+        mesh.pointsInstance(),
+        polyMesh::meshSubDir/"sets"
+    );
+    IOobjectList cellSets(objects.lookupClass(cellSet::typeName));
+    if (cellSets.size() > 0)
+    {
+        os  << "cellSets:" << endl;
+        forAllConstIter(IOobjectList, cellSets, iter)
+        {
+            cellSet set(*iter());
+            os  << '\t' << set.name() << "\tsize:" << set.size() << endl;
+        }
+    }
+    IOobjectList faceSets(objects.lookupClass(faceSet::typeName));
+    if (faceSets.size() > 0)
+    {
+        os  << "faceSets:" << endl;
+        forAllConstIter(IOobjectList, faceSets, iter)
+        {
+            faceSet set(*iter());
+            os  << '\t' << set.name() << "\tsize:" << set.size() << endl;
+        }
+    }
+    IOobjectList pointSets(objects.lookupClass(pointSet::typeName));
+    if (pointSets.size() > 0)
+    {
+        os  << "pointSets:" << endl;
+        forAllConstIter(IOobjectList, pointSets, iter)
+        {
+            pointSet set(*iter());
+            os  << '\t' << set.name() << "\tsize:" << set.size() << endl;
+        }
+    }
+    os  << endl;
+}
+
 
 // Read command and execute. Return true if ok, false otherwise.
 bool doCommand
@@ -531,6 +575,12 @@ commandStatus parseType
 
         return INVALID;
     }
+    else if (setType == "list")
+    {
+        printAllSets(mesh, Pout);
+
+        return INVALID;
+    }
     else if (setType == "time")
     {
         scalar time = readScalar(is);
@@ -612,8 +662,9 @@ commandStatus parseType
         (
             "commandStatus parseType(Time&, polyMesh&, const word&"
             ", IStringStream&)"
-        )   << "Illegal set type " << setType << endl
-            << "Should be one of 'cellSet' 'faceSet' 'pointSet'"
+        )   << "Illegal command " << setType << endl
+            << "Should be one of 'help', 'list', 'time' or a set type :"
+            << " 'cellSet', 'faceSet', 'pointSet'"
             << endl;
 
         return INVALID;

applications/utilities/mesh/manipulation/couplePatches/FoamX/couplePatches.cfg

@@ -1,27 +0,0 @@
-/*---------------------------------------------------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.0                                   |
-|   \\  /    A nd           | Web:      http://www.openfoam.org               |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-// couplePatches tool definition
-
-description     "Utility to reorder face numbering for cyclic and processor patches";
-
-couplePatchesDict
-{
-    type            dictionary;
-    description     "couplePatches control dictionary";
-    dictionaryPath  "system";
-
-    entries
-    {
-        arguments
-        {
-            type        rootCaseArguments;
-        }
-    }
-}
-
-// ************************************************************************* //

applications/utilities/mesh/manipulation/couplePatches/Make/files

@@ -1,3 +0,0 @@
-couplePatches.C
-
-EXE = $(FOAM_APPBIN)/couplePatches

applications/utilities/mesh/manipulation/couplePatches/Make/options

@@ -1,7 +0,0 @@
-EXE_INC = \
-    -I$(LIB_SRC)/dynamicMesh/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
-
-EXE_LIBS = \
-    -ldynamicMesh \
-    -lmeshTools

applications/utilities/mesh/manipulation/couplePatches/couplePatches.C

@@ -1,160 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation; either version 2 of the License, or (at your
-    option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Description
-    Utility to reorder cyclic and processor patches.
-
-    Uses dummy morph to sort things out.
-
-    Is bit of hack since polyMesh constructor already checks for coupled
-    face areas so might bomb out. You might need to compile in a local
-    processorPolyPatch that gives a warning instead to make this utility
-    complete.
-
-\*---------------------------------------------------------------------------*/
-
-#include "argList.H"
-#include "polyMesh.H"
-#include "Time.H"
-#include "polyTopoChange.H"
-#include "mapPolyMesh.H"
-#include "OFstream.H"
-#include "coupledPolyPatch.H"
-#include "PstreamReduceOps.H"
-
-using namespace Foam;
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-// Main program:
-
-int main(int argc, char *argv[])
-{
-#   include "setRootCase.H"
-#   include "createTime.H"
-#   include "createPolyMesh.H"
-
-    Info<< "Using geometry to calculate face correspondence across"
-        << " coupled boundaries (processor, cyclic)" << nl
-        << "This will only work for cyclics if they are parallel or"
-        << " their rotation is defined across the origin" << nl
-        << endl;
-
-    const polyBoundaryMesh& patches = mesh.boundaryMesh();
-
-    bool hasCoupled = false;
-    forAll(patches, patchI)
-    {
-        const polyPatch& pp = patches[patchI];
-
-        if (pp.coupled())
-        {
-            hasCoupled = true;
-
-            break;
-        }
-    }
-
-    reduce(hasCoupled, orOp<bool>());
-
-    if (hasCoupled)
-    {
-        Info<< "Mesh has coupled patches ..." << nl << endl;
-
-        // Dummy topo changes container
-        polyTopoChange meshMod(mesh);
-
-        // Do all changes
-        Info<< "Doing dummy mesh morph to correct face ordering ..."
-            << endl;
-
-        runTime++;
-
-        faceList oldFaces(mesh.faces());
-
-        autoPtr<mapPolyMesh> morphMap = meshMod.changeMesh(mesh, false);
-
-        // Update fields
-        mesh.updateMesh(morphMap);
-
-        // Move mesh (since morphing does not do this)
-        if (morphMap().hasMotionPoints())
-        {
-            mesh.movePoints(morphMap().preMotionPoints());
-        }
-
-        // Find out if anything changed.
-        // Problem is that - if the topo change is only rotation -
-        // face::operator= does not detect a change so use labelList::operator=
-        bool meshChanged = false;
-
-        if (mesh.faces().size() != oldFaces.size())
-        {
-            meshChanged = true;
-        }
-        else
-        {
-            forAll(mesh.faces(), faceI)
-            {
-                const labelList& f = mesh.faces()[faceI];
-                const labelList& oldF = oldFaces[faceI];
-
-                if (f != oldF)
-                {
-                    meshChanged = true;
-                    break;
-                }
-            }
-        }
-
-        reduce(meshChanged, orOp<label>());
-
-        if (meshChanged)
-        {
-            // Set the precision of the points data to 10
-            IOstream::defaultPrecision(10);
-
-            // Write resulting mesh
-            Info << "Writing morphed mesh to time " << runTime.value() << endl;
-
-            mesh.write();
-        }
-        else
-        {
-            Info << "Mesh ordering ok. Nothing changed." << endl;
-        }
-    }
-    else
-    {
-        Info<< "Mesh has no coupled patches. Nothing changed ..." << nl << endl;
-    }
-        
-
-    Info << "End\n" << endl;
-
-    return 0;
-}
-
-
-// ************************************************************************* //

src/Allwmake

@@ -41,6 +41,8 @@ wmake libso randomProcesses
 
 ( cd postProcessing && ./Allwmake )
 
+wmake libso autoMesh
+
 wmake libso errorEstimation
 
 

src/OpenFOAM/containers/Lists/ListOps/ListOps.H

@@ -1,3 +1,4 @@
+
 /*---------------------------------------------------------------------------*\
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
@@ -26,7 +27,7 @@ InNamspace
     Foam
 
 Description
-    ListOps. Various functions to operate on Lists.
+    Various functions to operate on Lists.
 
 SourceFiles
     ListOps.C
@@ -66,6 +67,17 @@ template<class List>
 void inplaceReorder(const labelList& oldToNew, List&);
 
 
+// Variants to work with iterators and sparse tables. Need to have iterators
+// and insert()
+
+//- Map values. Do not map negative values.
+template<class Container>
+void inplaceMapValue(const labelList& oldToNew, Container&);
+//- Recreate with mapped keys. Remove elements with negative key.
+template<class Container>
+void inplaceMapKey(const labelList& oldToNew, Container&);
+
+
 //- Extract elements of List whose region is certain value. Use e.g.
 //  to extract all selected elements:
 //  subset<boolList, labelList>(selectedElems, true, lst);
@@ -103,11 +115,21 @@ labelList identity(const label len);
 //- Find first occurence of given element and return index,
 //  return -1 if not found. Linear search.
 template<class List>
-label findIndex(const List&, typename List::const_reference);
+label findIndex
+(
+    const List&,
+    typename List::const_reference,
+    const label start=0
+);
 
 //- Find all occurences of given element. Linear search.
 template<class List>
-labelList findIndices(const List&, typename List::const_reference);
+labelList findIndices
+(
+    const List&,
+    typename List::const_reference,
+    const label start=0
+);
 
 //- Opposite of findIndices: set values at indices to given value
 template<class List>
@@ -123,7 +145,7 @@ template<class List>
 List createWithValues
 (
     const label sz,
-    typename List::const_reference initValue,
+    const typename List::const_reference initValue,
     const labelList& indices,
     typename List::const_reference setValue
 );
@@ -131,25 +153,35 @@ List createWithValues
 //- Find index of max element (and larger than given element).
 //  return -1 if not found. Linear search.
 template<class List>
-label findMax(const List&);
+label findMax(const List&, const label start=0);
 
 
 //- Find index of min element (and less than given element).
 //  return -1 if not found. Linear search.
 template<class List>
-label findMin(const List&);
+label findMin(const List&, const label start=0);
 
 
 //- Find first occurence of given element in sorted list and return index,
 //  return -1 if not found. Binary search.
 template<class List>
-label findSortedIndex(const List&, typename List::const_reference);
+label findSortedIndex
+(
+    const List&,
+    typename List::const_reference,
+    const label start=0
+);
 
 
 //- Find last element < given value in sorted list and return index,
 //  return -1 if not found. Binary search.
 template<class List>
-label findLower(const List&, typename List::const_reference);
+label findLower
+(
+    const List&,
+    typename List::const_reference,
+    const label start=0
+);
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/OpenFOAM/containers/Lists/ListOps/ListOpsTemplates.C

@@ -118,6 +118,54 @@ void Foam::inplaceReorder
 }
 
 
+template<class Container>
+void Foam::inplaceMapValue
+(
+    const labelList& oldToNew,
+    Container& lst
+)
+{
+    for
+    (
+        typename Container::iterator iter = lst.begin();
+        iter != lst.end();
+        ++iter
+    )
+    {
+        if (iter() >= 0)
+        {
+            iter() = oldToNew[iter()];
+        }
+    }
+}
+
+
+template<class Container>
+void Foam::inplaceMapKey
+(
+    const labelList& oldToNew,
+    Container& lst
+)
+{
+    Container newLst(lst);
+
+    for
+    (
+        typename Container::iterator iter = lst.begin();
+        iter != lst.end();
+        ++iter
+    )
+    {
+        if (iter.key() >= 0)
+        {
+            newLst.insert(oldToNew[iter.key()], iter());
+        }
+    }
+    
+    lst.transfer(newLst);
+}
+
+
 template<class T, class List>
 List Foam::subset(const UList<T>& regions, const T& region, const List& lst)
 {
@@ -221,11 +269,16 @@ void Foam::invertManyToMany
 
 
 template<class List>
-Foam::label Foam::findIndex(const List& l, typename List::const_reference t)
+Foam::label Foam::findIndex
+(
+    const List& l,
+    typename List::const_reference t,
+    const label start
+)
 {
     label index = -1;
 
-    forAll(l, i)
+    for (label i = start; i < l.size(); i++)
     {
         if (l[i] == t)
         {
@@ -242,13 +295,14 @@ template<class List>
 Foam::labelList Foam::findIndices
 (
     const List& l,
-    typename List::const_reference t
+    typename List::const_reference t,
+    const label start
 )
 {
     // Count occurrences
     label n = 0;
 
-    forAll(l, i)
+    for (label i = start; i < l.size(); i++)
     {
         if (l[i] == t)
         {
@@ -260,7 +314,7 @@ Foam::labelList Foam::findIndices
     labelList indices(n);
     n = 0;
 
-    forAll(l, i)
+    for (label i = start; i < l.size(); i++)
     {
         if (l[i] == t)
         {
@@ -291,7 +345,7 @@ template<class List>
 List Foam::createWithValues
 (
     const label sz,
-    typename List::const_reference initValue,
+    const typename List::const_reference initValue,
     const labelList& indices,
     typename List::const_reference setValue
 )
@@ -303,16 +357,16 @@ List Foam::createWithValues
 
 
 template<class List>
-Foam::label Foam::findMax(const List& l)
+Foam::label Foam::findMax(const List& l, const label start)
 {
-    if (l.size() == 0)
+    if (start >= l.size())
     {
         return -1;
     }
 
-    label index = 0;
+    label index = start;
 
-    for (label i = 1; i < l.size(); i++)
+    for (label i = start+1; i < l.size(); i++)
     {
         if (l[i] > l[index])
         {
@@ -325,16 +379,16 @@ Foam::label Foam::findMax(const List& l)
 
 
 template<class List>
-Foam::label Foam::findMin(const List& l)
+Foam::label Foam::findMin(const List& l, const label start)
 {
-    if (l.size() == 0)
+    if (start >= l.size())
     {
         return -1;
     }
 
-    label index = 0;
+    label index = start;
 
-    for (label i = 1; i < l.size(); i++)
+    for (label i = start+1; i < l.size(); i++)
     {
         if (l[i] < l[index])
         {
@@ -350,15 +404,16 @@ template<class List>
 Foam::label Foam::findSortedIndex
 (
     const List& l,
-    typename List::const_reference t
+    typename List::const_reference t,
+    const label start
 )
 {
-    if (l.size() == 0)
+    if (start >= l.size())
     {
         return -1;
     }
 
-    label low = 0;
+    label low = start;
     label high = l.size() - 1;
 
     while (low <= high)
@@ -387,15 +442,16 @@ template<class List>
 Foam::label Foam::findLower
 (
     const List& l,
-    typename List::const_reference t
+    typename List::const_reference t,
+    const label start
 )
 {
-    if (l.size() == 0)
+    if (start >= l.size())
     {
         return -1;
     }
 
-    label low = 0;
+    label low = start;
     label high = l.size() - 1;
 
     while ((high - low) > 1)

src/autoMesh/Make/files

@@ -0,0 +1,18 @@
+autoHexMesh             = autoHexMesh
+
+$(autoHexMesh)/autoHexMeshDriver/autoHexMeshDriver.C
+$(autoHexMesh)/autoHexMeshDriver/autoHexMeshDriverLayers.C
+$(autoHexMesh)/autoHexMeshDriver/autoHexMeshDriverShrink.C
+$(autoHexMesh)/autoHexMeshDriver/autoHexMeshDriverSnap.C
+$(autoHexMesh)/autoHexMeshDriver/pointData/pointData.C
+$(autoHexMesh)/meshRefinement/meshRefinementBaffles.C
+$(autoHexMesh)/meshRefinement/meshRefinement.C
+$(autoHexMesh)/meshRefinement/meshRefinementMerge.C
+$(autoHexMesh)/meshRefinement/meshRefinementRefine.C
+$(autoHexMesh)/refinementSurfaces/refinementSurfaces.C
+$(autoHexMesh)/offsetTriSurfaceMesh/offsetTriSurfaceMesh.C
+$(autoHexMesh)/trackedParticle/trackedParticle.C
+$(autoHexMesh)/trackedParticle/trackedParticleCloud.C
+
+LIB = $(FOAM_LIBBIN)/libautoMesh
+

src/autoMesh/Make/options

@@ -0,0 +1,21 @@
+/* include /home/pinky2/mattijs/pub/Delaunay/CGAL-3.1/make/makefile_i686_Linux-2.6.4-52-default_g++-3.4.3 */
+
+EXE_INC = \
+    /* -g -DFULLDEBUG -O0 $(CGAL_CXXFLAGS) */ \
+    -I$(FOAM_UTILITIES)/parallelProcessing/decompositionMethods/decompositionMethods/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(FOAM_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/edgeMesh/lnInclude \
+    -I$(LIB_SRC)/triSurface/lnInclude \
+
+LIB_LIBS = \
+    -ldecompositionMethods \
+    -ldynamicMesh \
+    -lfiniteVolume \
+    -llagrangian \
+    -lmeshTools \
+    -ledgeMesh \
+    -ltriSurface \
+ /*   $(CGAL_LDFLAGS) */

src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.H

@@ -0,0 +1,810 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    autoHexMeshDriver
+
+Description
+    main meshing driver.
+
+SourceFiles
+    autoHexMeshDriver.C
+    autoHexMeshDriverSnap.C
+    autoHexMeshDriverLayers.C
+    autoHexMeshDriverShrink.C
+    autoHexMeshDriverTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef autoHexMeshDriver_H
+#define autoHexMeshDriver_H
+
+#include "autoPtr.H"
+#include "dictionary.H"
+#include "pointField.H"
+#include "boolList.H"
+#include "Switch.H"
+#include "PackedList.H"
+#include "wallPoint.H"
+#include "indirectPrimitivePatch.H"
+#include "featureEdgeMesh.H"
+#include "searchableSurface.H"
+#include "refinementSurfaces.H"
+#include "meshRefinement.H"
+#include "decompositionMethod.H"
+#include "fvMeshDistribute.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Class forward declarations
+class fvMesh;
+class pointMesh;
+class motionSmoother;
+class removePoints;
+class pointSet;
+class pointData;
+class faceSet;
+class addPatchCellLayer;
+class mapDistributePolyMesh;
+
+/*---------------------------------------------------------------------------*\
+                           Class autoHexMeshDriver Declaration
+\*---------------------------------------------------------------------------*/
+
+class autoHexMeshDriver
+{
+    // Static data members
+
+        //- Default angle for faces to be convcave
+        static const scalar defaultConcaveAngle;
+
+        //- Extrusion controls
+        enum extrudeMode
+        {
+            NOEXTRUDE,      // Do not extrude. No layers added.
+            EXTRUDE,        // Extrude
+            EXTRUDEREMOVE   // Extrude but afterwards remove added faces locally
+        };
+
+
+    // Private classes
+
+        //- Combine operator class for equalizing displacements.
+        class minMagEqOp
+        {
+        public:
+
+            void operator()(vector& x, const vector& y) const
+            {
+                if (magSqr(y) < magSqr(x))
+                {
+                    x = y;
+                }
+            }
+        };
+
+        // Combine operator class to combine normal with other normal.
+        class nomalsCombine
+        {
+        public:
+
+            void operator()(vector& x, const vector& y) const
+            {
+                if (y != wallPoint::greatPoint)
+                {
+                    if (x == wallPoint::greatPoint)
+                    {
+                        x = y;
+                    }
+                    else
+                    {
+                        x *= (x&y);
+                    }
+                }
+            }
+        };
+
+
+
+    // Private data
+
+        //- Reference to mesh
+        fvMesh& mesh_;
+
+        //- Input dictionarys
+        const dictionary dict_;
+
+        //- Debug level
+        const label debug_;
+
+        //- Total number of cells
+        const label maxGlobalCells_;
+
+        //- Per processor max number of cells
+        const label maxLocalCells_;
+
+        //- When to stop refining
+        const label minRefineCells_;
+
+        //- Curvature
+        const scalar curvature_;
+
+        //- Number of layers between different refinement levels
+        const label nBufferLayers_;
+
+        //- Areas to keep
+        const pointField keepPoints_;
+
+        //- Merge distance
+        const scalar mergeDist_;
+
+
+
+
+        //- Explicit features
+        PtrList<featureEdgeMesh> featureMeshes_;
+        //- Per feature the refinement level
+        labelList featureLevels_;
+
+        //- Shells
+        PtrList<searchableSurface> shells_;
+        //- Per shell the refinement level
+        labelList shellLevels_;
+        //- Per shell whether to refine inside or outside
+        boolList shellRefineInside_;
+
+        //- Surfaces with refinement levels built-in
+        autoPtr<refinementSurfaces> surfacesPtr_;
+
+        //- Per refinement surface region the patch
+        labelList globalToPatch_;
+
+        ////- Per refinement surface with a valid faceZone the cyclicpatch
+        ////  (or -1)
+        //labelList surfaceToCyclicPatch_;
+
+        //- Mesh refinement engine
+        autoPtr<meshRefinement> meshRefinerPtr_;
+
+        //- Decomposition engine
+        autoPtr<decompositionMethod> decomposerPtr_;
+
+        //- Mesh distribution engine
+        autoPtr<fvMeshDistribute> distributorPtr_;
+
+
+    // Private Member Functions
+
+        // Refinement
+
+            //- Get merge tolerance. Check against writing tolerance.
+            scalar getMergeDistance() const;
+
+            // Return per keeppoint -1 or the local cell label the point is in.
+            // Guaranteed to be only on one processor.
+            labelList findCells(const pointField&) const;
+
+            static void orientOutside(PtrList<searchableSurface>&);
+
+            //- Read feature edges
+            label readFeatureEdges();
+
+        // Snapping
+
+             //- Split surfaces into non-zoned and zones ones
+            void getZonedSurfaces(labelList&, labelList&) const;
+
+            //- Get faces to repatch. Returns map from face to patch.
+            Map<label> getZoneBafflePatches(const bool allowBoundary) const;
+
+            //- Calculates (geometric) shared points
+            static label getCollocatedPoints
+            (
+                const scalar tol,
+                const pointField&,
+                PackedList<1>&
+            );
+
+            //- Calculate displacement per patch point to smooth out patch.
+            //  Quite complicated in determining which points to move where.
+            pointField smoothPatchDisplacement(const motionSmoother&) const;
+
+            //- Check that face zones are synced
+            void checkCoupledFaceZones() const;
+
+            //- Per edge distance to patch
+            static tmp<scalarField> edgePatchDist
+            (
+                const pointMesh&,
+                const indirectPrimitivePatch&
+            );
+
+            //- Write displacement as .obj file.
+            static void dumpMove
+            (
+                const fileName&,
+                const pointField&,
+                const pointField&
+            );
+
+            //- Check displacement is outwards pointing
+            static bool outwardsDisplacement
+            (
+                const indirectPrimitivePatch&,
+                const vectorField&
+            );
+
+        // Face merging
+
+            //- Merge patch faces. Undo until no checkMesh errors.
+            label mergePatchFacesUndo
+            (
+                const scalar minCos,
+                const scalar concaveCos,
+                const dictionary&
+            );
+
+            //- Remove points.
+            autoPtr<mapPolyMesh> doRemovePoints
+            (
+                removePoints& pointRemover,
+                const boolList& pointCanBeDeleted
+            );
+
+            //- Restore faces (which contain removed points)
+            autoPtr<mapPolyMesh> doRestorePoints
+            (
+                removePoints& pointRemover,
+                const labelList& facesToRestore
+            );
+
+            //- Return candidateFaces that are also in set.
+            labelList collectFaces
+            (
+                const labelList& candidateFaces,
+                const labelHashSet& set
+            ) const;
+
+            //- Pick up faces of cells of faces in set.
+            labelList growFaceCellFace(const labelHashSet&) const;
+
+            //- Remove points not used by any face or points used by only
+            //  two faces where the edges are in line
+            label mergeEdgesUndo(const scalar minCos, const dictionary&);
+
+
+        // Layers
+
+            //- For debugging: Dump displacement to .obj files
+            static void dumpDisplacement
+            (
+                const fileName&,
+                const indirectPrimitivePatch&,
+                const vectorField&,
+                const List<extrudeMode>&
+            );
+
+            //- Check that primitivePatch is not multiply connected.
+            //  Collect non-manifold points in pointSet.
+            static void checkManifold
+            (
+                const indirectPrimitivePatch&,
+                pointSet& nonManifoldPoints
+            );
+
+
+            // Static extrusion setup
+
+                //- Unset extrusion on point. Returns true if anything unset.
+                static bool unmarkExtrusion
+                (
+                    const label patchPointI,
+                    pointField& patchDisp,
+                    labelList& patchNLayers,
+                    List<extrudeMode>& extrudeStatus
+                );
+
+                //- Unset extrusion on face. Returns true if anything unset.
+                static bool unmarkExtrusion
+                (
+                    const face& localFace,
+                    pointField& patchDisp,
+                    labelList& patchNLayers,
+                    List<extrudeMode>& extrudeStatus
+                );
+
+                //- No extrusion at non-manifold points.
+                void handleNonManifolds
+                (
+                    const indirectPrimitivePatch& pp,
+                    const labelList& meshEdges,
+                    pointField& patchDisp,
+                    labelList& patchNLayers,
+                    List<extrudeMode>& extrudeStatus
+                ) const;
+
+                //- No extrusion on feature edges. Assumes non-manifold
+                //  edges already handled.
+                void handleFeatureAngle
+                (
+                    const indirectPrimitivePatch& pp,
+                    const labelList& meshEdges,
+                    const scalar minCos,
+                    pointField& patchDisp,
+                    labelList& patchNLayers,
+                    List<extrudeMode>& extrudeStatus
+                ) const;
+
+                //- No extrusion on warped faces
+                void handleWarpedFaces
+                (
+                    const indirectPrimitivePatch& pp,
+                    const scalar faceRatio,
+                    const scalar edge0Len,
+                    const labelList& cellLevel,
+                    pointField& patchDisp,
+                    labelList& patchNLayers,
+                    List<extrudeMode>& extrudeStatus
+                ) const;
+
+                //- Determine the number of layers per point from the number of
+                //  layers per surface.
+                void setNumLayers
+                (
+                    const labelList& nLayers,
+                    const labelList& patchIDs,
+                    const indirectPrimitivePatch& pp,
+                    pointField& patchDisp,
+                    labelList& patchNLayers,
+                    List<extrudeMode>& extrudeStatus
+                ) const;
+
+                //- Grow no-extrusion layer.
+                static void growNoExtrusion
+                (
+                    const indirectPrimitivePatch& pp,
+                    pointField& patchDisp,
+                    labelList& patchNLayers,
+                    List<extrudeMode>& extrudeStatus
+                );
+
+                //- Calculate pointwise wanted and minimum thickness.
+                //  thickness: wanted thickness
+                //  minthickness: when to give up and not extrude
+                void calculateLayerThickness
+                (
+                    const scalar expansionRatio,
+                    const scalar finalLayerRatio,
+                    const scalar relMinThickness,
+                    const indirectPrimitivePatch& pp,
+                    const labelList& cellLevel,
+                    const labelList& patchNLayers,
+                    const scalar edge0Len,
+                    scalarField& thickness,
+                    scalarField& minThickness
+                ) const;
+
+
+            // Extrusion execution
+
+                //- Synchronize displacement among coupled patches.
+                void syncPatchDisplacement
+                (
+                    const motionSmoother& meshMover,
+                    const scalarField& minThickness,
+                    pointField& patchDisp,
+                    labelList& patchNLayers,
+                    List<extrudeMode>& extrudeStatus
+                ) const;
+
+                //- Get nearest point on surface to snap to
+                void getPatchDisplacement
+                (
+                    const motionSmoother& meshMover,
+                    const scalarField& thickness,
+                    const scalarField& minThickness,
+                    pointField& patchDisp,
+                    labelList& patchNLayers,
+                    List<extrudeMode>& extrudeStatus
+                ) const;
+
+                //- Truncates displacement
+                // - for all patchFaces in the faceset displacement gets set
+                //   to zero
+                // - all displacement < minThickness gets set to zero
+                label truncateDisplacement
+                (
+                    const motionSmoother& meshMover,
+                    const scalarField& minThickness,
+                    const faceSet& illegalPatchFaces,
+                    pointField& patchDisp,
+                    labelList& patchNLayers,
+                    List<extrudeMode>& extrudeStatus
+                ) const;
+
+                //- Setup layer information (at points and faces) to
+                //  modify mesh topology in
+                //  regions where layer mesh terminates. Guarantees an
+                //  optional slow decreasing of the number of layers.
+                //  Returns the number of layers per face and per point
+                //  to go into the actual layer addition engine.
+                void setupLayerInfoTruncation
+                (
+                    const motionSmoother& meshMover,
+                    const labelList& patchNLayers,
+                    const List<extrudeMode>& extrudeStatus,
+                    const label nBufferCellsNoExtrude,
+                    labelList& nPatchPointLayers,
+                    labelList& nPatchFaceLayers
+                ) const;
+
+                //- Does any of the cells use a face from faces?
+                static bool cellsUseFace
+                (
+                    const polyMesh& mesh,
+                    const labelList& cellLabels,
+                    const labelHashSet& faces
+                );
+
+                //- Checks the newly added cells and locally unmarks points
+                //  so they will not get extruded next time round. Returns
+                //  global number of unmarked points (0 if all was fine)
+                static label checkAndUnmark
+                (
+                    const addPatchCellLayer& addLayer,
+                    const dictionary& motionDict,
+                    const indirectPrimitivePatch& pp,
+                    const fvMesh&,
+
+                    pointField& patchDisp,
+                    labelList& patchNLayers,
+                    List<extrudeMode>& extrudeStatus
+                );
+
+                //- Count global number of extruded faces
+                static label countExtrusion
+                (
+                    const indirectPrimitivePatch& pp,
+                    const List<extrudeMode>& extrudeStatus
+                );
+
+                //- Collect layer faces and layer cells into bools
+                //  for ease of handling
+                static void getLayerCellsFaces
+                (
+                    const polyMesh&,
+                    const addPatchCellLayer&,
+                    boolList&,
+                    boolList&
+                );
+
+            // Mesh shrinking (to create space for layers)
+
+                //- Average field (over all subset of mesh points) by
+                //  summing contribution from edges. Global parallel since only
+                //  does master edges for coupled edges.
+                template<class Type>
+                void averageNeighbours
+                (
+                    const PackedList<1>& isMasterEdge,
+                    const labelList& meshEdges,
+                    const labelList& meshPoints,
+                    const edgeList& edges,
+                    const scalarField& invSumWeight,
+                    const Field<Type>& data,
+                    Field<Type>& average
+                ) const;
+
+                //- Calculate inverse sum of edge weights (currently always 1.0)
+                void sumWeights
+                (
+                    const PackedList<1>& isMasterEdge,
+                    const labelList& meshEdges,
+                    const labelList& meshPoints,
+                    const edgeList& edges,
+                    scalarField& invSumWeight
+                ) const;
+
+                //- Smooth scalar field on patch
+                void smoothField
+                (
+                    const motionSmoother& meshMover,
+                    const PackedList<1>& isMasterEdge,
+                    const labelList& meshEdges,
+                    const scalarField& fieldMin,
+                    const label& nSmoothDisp,
+                    scalarField& field
+                ) const;
+
+                //- Smooth normals on patch.
+                void smoothPatchNormals
+                (
+                    const motionSmoother& meshMover,
+                    const PackedList<1>& isMasterEdge,
+                    const labelList& meshEdges,
+                    const label nSmoothDisp,
+                    pointField& normals
+                ) const;
+
+                //- Smooth normals in interior.
+                void smoothNormals
+                (
+                    const label nSmoothDisp,
+                    const PackedList<1>& isMasterEdge,
+                    const labelList& fixedPoints,
+                    pointVectorField& normals
+                ) const;
+
+                bool isMaxEdge
+                (
+                    const List<pointData>&,
+                    const label edgeI,
+                    const scalar minCos
+                ) const;
+
+                //- Stop layer growth where mesh wraps around edge with a
+                //  large feature angle
+                void handleFeatureAngleLayerTerminations
+                (
+                    const indirectPrimitivePatch& pp,
+                    const scalar minCos,
+                    List<extrudeMode>& extrudeStatus,
+                    pointField& patchDisp,
+                    labelList& patchNLayers,
+                    label& nPointCounter
+                ) const;
+
+                //- Find isolated islands (points, edges and faces and
+                // layer terminations)
+                // in the layer mesh and stop any layer growth at these points.
+                void findIsolatedRegions
+                (
+                    const indirectPrimitivePatch& pp,
+                    const PackedList<1>& isMasterEdge,
+                    const labelList& meshEdges,
+                    const scalar minCosLayerTermination,
+                    scalarField& field,
+                    List<extrudeMode>& extrudeStatus,
+                    pointField& patchDisp,
+                    labelList& patchNLayers
+                ) const;
+
+                // Calculate medial axis fields
+                void medialAxisSmoothingInfo
+                (
+                    const motionSmoother& meshMover,
+                    const label nSmoothNormals,
+                    const label nSmoothSurfaceNormals,
+                    const scalar minMedianAxisAngleCos,
+
+                    pointVectorField& dispVec,
+                    pointScalarField& medialRatio,
+                    pointScalarField& medialDist
+                ) const;
+
+                //- Main routine to shrink mesh
+                void shrinkMeshMedialDistance
+                (
+                    motionSmoother& meshMover,
+                    const label nSmoothThickness,
+                    const scalar maxThicknessToMedialRatio,
+                    const label nAllowableErrors,
+                    const label nSnap,
+                    const scalar minCosLayerTermination,
+
+                    const scalarField& layerThickness,
+                    const scalarField& minThickness,
+
+                    const pointVectorField& dispVec,
+                    const pointScalarField& medialRatio,
+                    const pointScalarField& medialDist,
+
+                    List<extrudeMode>& extrudeStatus,
+                    pointField& patchDisp,
+                    labelList& patchNLayers
+                ) const;
+
+
+        //- Disallow default bitwise copy construct
+        autoHexMeshDriver(const autoHexMeshDriver&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const autoHexMeshDriver&);
+
+public:
+
+    //- Runtime type information
+    ClassName("autoHexMeshDriver");
+
+
+    // Constructors
+
+        //- Construct from dictionary and mesh to modify
+        autoHexMeshDriver
+        (
+            fvMesh& mesh,
+            const dictionary& meshDict,
+            const dictionary& decomposeDict
+        );
+
+
+    // Member Functions
+
+        // Access
+
+            //- reference to mesh
+            const fvMesh& mesh() const
+            {
+                return mesh_;
+            }
+            fvMesh& mesh()
+            {
+                return mesh_;
+            }
+
+            //- Surfaces to base refinement on
+            const refinementSurfaces& surfaces() const
+            {
+                return surfacesPtr_();
+            }
+
+            //- Per refinementsurface, per region the patch
+            const labelList& globalToPatch() const
+            {
+                return globalToPatch_;
+            }
+
+
+        // Refinement
+
+            //- Refine around explicit feature edges
+            label featureEdgeRefine
+            (
+                const label maxIter,
+                const label minRefine
+            );
+
+            //- Refine at surface intersections
+            label surfaceOnlyRefine(const label maxIter);
+
+            //- Remove cells not reachable from keep points
+            void removeInsideCells(const label nBufferLayers);
+
+            //- Refine volume regions (interior of shells)
+            label shellRefine(const label maxIter);
+
+            //- Introduce baffles; remove unreachable mesh parts
+            //  handleSnapProblems : whether to remove free floating cells
+            void baffleAndSplitMesh(const bool handleSnapProblems);
+
+            //- Move cells to zones
+            void zonify();
+
+            //- Split and recombine baffles to get rid of single face baffles.
+            void splitAndMergeBaffles(const bool handleSnapProblems);
+
+            //- Merge multiple boundary faces on single cell
+            void mergePatchFaces();
+
+            //- Redecompose according to cell count
+            //  keepZoneFaces : find all faceZones from zoned surfaces and keep
+            //                  owner and neighbour together
+            //  keepBaffles   : find all baffles and keep them together
+            autoPtr<mapDistributePolyMesh> balance
+            (
+                const bool keepZoneFaces,
+                const bool keepBaffles
+            );
+
+            //- Write mesh
+            void writeMesh(const string&) const;
+
+
+        // Snapping
+
+            //- Create baffles for faces straddling zoned surfaces. Return
+            //  baffles.
+            autoPtr<mapPolyMesh> createZoneBaffles(List<labelPair>&);
+
+            //- Merge baffles.
+            autoPtr<mapPolyMesh> mergeZoneBaffles(const List<labelPair>&);
+
+            //- Calculate edge length per patch point.
+            scalarField calcSnapDistance(const indirectPrimitivePatch&) const;
+
+            //- Get patches generated for surfaces.
+            labelList getSurfacePatches() const;
+
+            //- Smooth the mesh (patch and internal) to increase visibility
+            //  of surface points (on castellated mesh) w.r.t. surface.
+            void preSmoothPatch
+            (
+                const label nInitErrors,
+                const List<labelPair>& baffles,
+                motionSmoother&
+            ) const;
+
+            //- Per patch point calculate point on nearest surface. Set as
+            //  boundary conditions of motionSmoother displacement field. Return
+            //  displacement of patch points.
+            vectorField calcNearestSurface
+            (
+                const scalarField& snapDist,
+                motionSmoother& meshMover
+            ) const;
+
+            //- Smooth the displacement field to the internal.
+            void smoothDisplacement(motionSmoother&) const;
+
+            //- Do the hard work: move the mesh according to displacement,
+            //  locally relax the displacement.
+            void scaleMesh
+            (
+                const label nInitErrors,
+                const List<labelPair>& baffles,
+                motionSmoother&
+            );
+
+
+        // Layers
+
+            //- Merge patch faces on same cell.
+            void mergePatchFacesUndo();
+
+            //- Read layer per region
+            labelList readNumLayers() const;
+
+            //- Check that mesh outside is not multiply connected.
+            void checkMeshManifold() const;
+
+            //- Add cell layers
+            void addLayers(const scalar nAllowableErrors, motionSmoother&);
+
+
+        // Other
+
+            //- Do all : refine, snap, layers
+            void doMesh();
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "autoHexMeshDriverTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriverTemplates.C

@@ -0,0 +1,75 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "autoHexMeshDriver.H"
+#include "syncTools.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+void Foam::autoHexMeshDriver::averageNeighbours
+(
+    const PackedList<1>& isMasterEdge,
+    const labelList& meshEdges,
+    const labelList& meshPoints,
+    const edgeList& edges,
+    const scalarField& invSumWeight,
+    const Field<Type>& data,
+    Field<Type>& average
+) const
+{
+    average = pTraits<Type>::zero;
+
+    forAll(edges, edgeI)
+    {
+        if (isMasterEdge.get(meshEdges[edgeI]) == 1)
+        {
+            const edge& e = edges[edgeI];
+            //scalar eWeight = edgeWeights[edgeI];
+            scalar eWeight =  1.0;
+            label v0 = e[0];
+            label v1 = e[1];
+
+            average[v0] += eWeight*data[v1];
+            average[v1] += eWeight*data[v0];
+        }
+    }
+
+    syncTools::syncPointList
+    (
+        mesh_,
+        meshPoints,
+        average,
+        plusEqOp<Type>(),
+        pTraits<Type>::zero,    // null value
+        false                   // no separation
+    );
+
+    average *= invSumWeight;
+}
+
+
+// ************************************************************************* //

src/autoMesh/autoHexMesh/autoHexMeshDriver/pointData/pointData.C

@@ -0,0 +1,53 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+
+\*---------------------------------------------------------------------------*/
+
+#include "pointData.H"
+
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<(Ostream& os, const pointData& wDist)
+{
+    if (os.format() == IOstream::ASCII)
+    {
+        return os
+            << wDist.origin() << token::SPACE << wDist.distSqr()
+            << token::SPACE << wDist.s() << token::SPACE << wDist.v();
+    }
+    else
+    {
+        return os
+            << wDist.origin() << wDist.distSqr() << wDist.s() << wDist.v();
+    }
+}
+
+Foam::Istream& Foam::operator>>(Istream& is, pointData& wDist)
+{
+    return is >> wDist.origin_ >> wDist.distSqr_ >> wDist.s_ >> wDist.v_;
+}
+
+// ************************************************************************* //

src/autoMesh/autoHexMesh/autoHexMeshDriver/pointData/pointData.H

@@ -0,0 +1,228 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    pointData
+
+Description
+    Holds information regarding nearest wall point. Used in pointEdgeWave.
+    (so not standard meshWave)
+    To be used in wall distance calculation.
+
+SourceFiles
+    pointDataI.H
+    pointData.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef pointData_H
+#define pointData_H
+
+#include "point.H"
+#include "label.H"
+#include "scalar.H"
+#include "tensor.H"
+#include "pTraits.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Class forward declarations
+class polyPatch;
+class polyMesh;
+
+/*---------------------------------------------------------------------------*\
+                           Class pointData Declaration
+\*---------------------------------------------------------------------------*/
+
+class pointData
+{
+    // Private data
+
+        //- position of nearest wall center
+        point origin_;
+
+        //- normal distance (squared) from point to origin
+        scalar distSqr_;
+
+        //- additional information.
+        scalar s_;
+
+        //- additional information.
+        vector v_;
+
+
+    // Private Member Functions
+
+        //- Evaluate distance to point. Update distSqr, origin from whomever
+        //  is nearer pt. Return true if w2 is closer to point,
+        //  false otherwise.
+        inline bool update
+        (
+            const point&,
+            const pointData& w2,
+            const scalar tol
+        );
+
+        //- Combine current with w2. Update distSqr, origin if w2 has smaller
+        //  quantities and returns true.
+        inline bool update
+        (
+            const pointData& w2,
+            const scalar tol
+        );
+
+public:
+
+    // Constructors
+
+        //- Construct null
+        inline pointData();
+
+        //- Construct from origin, distance
+        inline pointData
+        (
+            const point& origin,
+            const scalar distSqr,
+            const scalar s,
+            const vector& v
+        );
+
+        //- Construct as copy
+        inline pointData(const pointData&);
+
+
+    // Member Functions
+
+        // Access
+
+            inline const point& origin() const;
+
+            inline scalar distSqr() const;
+
+            inline scalar s() const;
+
+            inline const vector& v() const;
+
+
+        // Needed by meshWave
+
+            //- Check whether origin has been changed at all or
+            //  still contains original (invalid) value.
+            inline bool valid() const;
+
+            //- Check for identical geometrical data. Used for cyclics checking.
+            inline bool sameGeometry(const pointData&, const scalar tol)
+             const;
+
+            //- Convert origin to relative vector to leaving point
+            //  (= point coordinate)
+            inline void leaveDomain
+            (
+                const polyPatch& patch,
+                const label patchPointI,
+                const point& pos
+            );
+
+            //- Convert relative origin to absolute by adding entering point
+            inline void enterDomain
+            (
+                const polyPatch& patch,
+                const label patchPointI,
+                const point& pos
+            );
+
+            //- Apply rotation matrix to origin
+            inline void transform(const tensor& rotTensor);
+
+            //- Influence of edge on point
+            inline bool updatePoint
+            (
+                const polyMesh& mesh,
+                const label pointI,
+                const label edgeI,
+                const pointData& edgeInfo,
+                const scalar tol
+            );
+
+            //- Influence of different value on same point.
+            //  Merge new and old info.
+            inline bool updatePoint
+            (
+                const polyMesh& mesh,
+                const label pointI,
+                const pointData& newPointInfo,
+                const scalar tol
+            );
+
+            //- Influence of different value on same point.
+            //  No information about current position whatsoever.
+            inline bool updatePoint
+            (
+                const pointData& newPointInfo,
+                const scalar tol
+            );
+
+            //- Influence of point on edge.
+            inline bool updateEdge
+            (
+                const polyMesh& mesh,
+                const label edgeI,
+                const label pointI,
+                const pointData& pointInfo,
+                const scalar tol
+            );
+
+
+    // Member Operators
+
+        //Note: Used to determine whether to call update.
+        inline bool operator==(const pointData&) const;
+
+        inline bool operator!=(const pointData&) const;
+
+
+    // IOstream Operators
+
+        friend Ostream& operator<<(Ostream&, const pointData&);
+        friend Istream& operator>>(Istream&, pointData&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "pointDataI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/autoMesh/autoHexMesh/autoHexMeshDriver/pointData/pointDataI.H

@@ -0,0 +1,356 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+
+\*---------------------------------------------------------------------------*/
+
+#include "polyMesh.H"
+#include "transform.H"
+#include "wallPoint.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+// Update this with w2 if w2 nearer to pt.
+inline bool Foam::pointData::update
+(
+    const point& pt,
+    const pointData& w2,
+    const scalar tol
+)
+{
+    scalar dist2 = magSqr(pt - w2.origin());
+
+    if (!valid())
+    {
+        distSqr_ = dist2;
+        origin_ = w2.origin();
+        s_ = w2.s();
+        v_ = w2.v();
+
+        return true;
+    }       
+
+
+//    if (v_ != w2.v())
+//    {
+//        return false;
+//    }
+
+
+    scalar diff = distSqr_ - dist2;
+
+    if (diff < 0)
+    {
+        // already nearer to pt
+        return false;
+    }
+
+    if ((diff < SMALL) || ((distSqr_ > SMALL) && (diff/distSqr_ < tol)))
+    {
+        // don't propagate small changes
+        return false;
+    }
+    else
+    {
+        // update with new values
+        distSqr_ = dist2;
+        origin_ = w2.origin();
+        s_ = w2.s();
+        v_ = w2.v();
+
+        return true;
+    }
+}
+
+    
+// Update this with w2 (information on same point)
+inline bool Foam::pointData::update
+(
+    const pointData& w2,
+    const scalar tol
+)
+{
+    if (!valid())
+    {
+        // current not yet set so use any value
+        distSqr_ = w2.distSqr();
+        origin_ = w2.origin();
+        s_ = w2.s();
+        v_ = w2.v();
+        return true;
+    }        
+
+
+//    if (v_ != w2.v())
+//    {
+//        return false;
+//    }
+
+
+    scalar diff = distSqr_ - w2.distSqr();
+
+    if (diff < 0)
+    {
+        // already nearer to pt
+        return false;
+    }
+
+    if ((diff < SMALL) || ((distSqr_ > SMALL) && (diff/distSqr_ < tol)))
+    {
+        // don't propagate small changes
+        return false;
+    }
+    else
+    {
+        // update with new values
+        distSqr_ =  w2.distSqr();
+        origin_ = w2.origin();
+        s_ = w2.s();
+        v_ = w2.v();
+
+        return true;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Null constructor
+inline Foam::pointData::pointData()
+:
+    origin_(wallPoint::greatPoint),
+    distSqr_(GREAT),
+    s_(GREAT),
+    v_(wallPoint::greatPoint)
+{}
+
+
+// Construct from origin, distance
+inline Foam::pointData::pointData
+(
+    const point& origin, 
+    const scalar distSqr,
+    const scalar s,
+    const vector& v
+)
+:
+    origin_(origin),
+    distSqr_(distSqr),
+    s_(s),
+    v_(v)
+{}
+
+
+// Construct as copy
+inline Foam::pointData::pointData(const pointData& wpt)
+:
+    origin_(wpt.origin()),
+    distSqr_(wpt.distSqr()),
+    s_(wpt.s()),
+    v_(wpt.v())
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline const Foam::point& Foam::pointData::origin() const
+{
+    return origin_;
+}
+
+
+inline Foam::scalar Foam::pointData::distSqr() const
+{
+    return distSqr_;
+}
+
+
+inline Foam::scalar Foam::pointData::s() const
+{
+    return s_;
+}
+
+
+inline const Foam::vector& Foam::pointData::v() const
+{
+    return v_;
+}
+
+
+inline bool Foam::pointData::valid() const
+{
+    return origin_ != wallPoint::greatPoint;
+}
+
+
+// Checks for cyclic points
+inline bool Foam::pointData::sameGeometry
+(
+    const pointData& w2,
+    const scalar tol
+) const
+{
+    scalar diff = Foam::mag(distSqr() - w2.distSqr());
+
+    if (diff < SMALL)
+    {
+        return true;
+    }
+    else
+    {
+        if ((distSqr() > SMALL) && ((diff/distSqr()) < tol))
+        {
+            return true;
+        }
+        else
+        {
+            return false;
+        }
+    }
+}
+
+
+inline void Foam::pointData::leaveDomain
+(
+    const polyPatch& patch,
+    const label patchPointI,
+    const point& coord
+)
+{
+    origin_ -= coord;
+}
+
+
+inline void Foam::pointData::transform(const tensor& rotTensor)
+{
+    origin_ = Foam::transform(rotTensor, origin_);
+}
+
+
+// Update absolute geometric quantities. Note that distance (distSqr_)
+// is not affected by leaving/entering domain.
+inline void Foam::pointData::enterDomain
+(
+    const polyPatch& patch,
+    const label patchPointI,
+    const point& coord
+)
+{
+    // back to absolute form
+    origin_ += coord;
+}
+
+
+// Update this with information from connected edge
+inline bool Foam::pointData::updatePoint
+(
+    const polyMesh& mesh,
+    const label pointI,
+    const label edgeI,
+    const pointData& edgeInfo,
+    const scalar tol
+)
+{
+    return 
+        update
+        (
+            mesh.points()[pointI],
+            edgeInfo,
+            tol
+        );
+    }    
+
+
+// Update this with new information on same point
+inline bool Foam::pointData::updatePoint
+(
+    const polyMesh& mesh,
+    const label pointI,
+    const pointData& newPointInfo,
+    const scalar tol
+)
+{
+    return
+        update
+        (
+            mesh.points()[pointI],
+            newPointInfo,
+            tol
+        );
+}    
+
+
+// Update this with new information on same point. No extra information.
+inline bool Foam::pointData::updatePoint
+(
+    const pointData& newPointInfo,
+    const scalar tol
+)
+{
+    return update(newPointInfo, tol);
+}    
+
+
+// Update this with information from connected point
+inline bool Foam::pointData::updateEdge
+(
+    const polyMesh& mesh,
+    const label edgeI,
+    const label pointI,
+    const pointData& pointInfo,
+    const scalar tol
+)
+{
+    const pointField& points = mesh.points();
+
+    const edge& e = mesh.edges()[edgeI];
+
+    const point edgeMid(0.5*(points[e[0]] + points[e[1]]));
+
+    return
+        update
+        (
+            edgeMid,
+            pointInfo,
+            tol
+        );
+}    
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+inline bool Foam::pointData::operator==(const pointData& rhs) const
+{
+    return origin() == rhs.origin();
+}
+
+
+inline bool Foam::pointData::operator!=(const pointData& rhs) const
+{
+    return !(*this == rhs);
+}
+
+
+// ************************************************************************* //

src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H

@@ -0,0 +1,667 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    meshRefinement
+
+Description
+    Helper class which maintains intersections of (changing) mesh with
+    (static) surfaces.
+    Maintains
+    - per face any intersections of this edge with any of the surfaces
+
+SourceFiles
+    meshRefinement.C
+    meshRefinementBaffles.C
+    meshRefinementMerge.C
+    meshRefinementRefine.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef meshRefinement_H
+#define meshRefinement_H
+
+#include "hexRef8.H"
+#include "mapPolyMesh.H"
+#include "autoPtr.H"
+#include "pointIOField.H"
+#include "labelPair.H"
+#include "indirectPrimitivePatch.H"
+#include "pointFieldsFwd.H"
+#include "Tuple2.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Class forward declarations
+class fvMesh;
+class mapDistributePolyMesh;
+class decompositionMethod;
+class refinementSurfaces;
+class removeCells;
+class featureEdgeMesh;
+class fvMeshDistribute;
+class searchableSurface;
+class regionSplit;
+
+/*---------------------------------------------------------------------------*\
+                           Class meshRefinement Declaration
+\*---------------------------------------------------------------------------*/
+
+class meshRefinement
+{
+
+public:
+
+    // Public data types
+
+        //- Enumeration for debug dumping
+        enum writeFlag
+        {
+            MESH = 1,
+            SCALARLEVELS = 2,
+            OBJINTERSECTIONS = 4
+        };
+
+
+        //- Enumeration for how the userdata is to be mapped upon refinement.
+        enum mapType
+        {
+            MASTERONLY = 1, // maintain master only
+            KEEPALL = 2,    // have slaves (upon refinement) from master
+            REMOVE = 4      // set value to -1 any face that has been refined
+        };
+
+
+private:
+
+    // Private data
+
+        //- Reference to mesh
+        fvMesh& mesh_;
+
+        //- tolerance used for sorting coordinates (used in 'less' routine)
+        const scalar tol_;
+
+        const refinementSurfaces& surfaces_;
+
+        //- refinement engine
+        hexRef8 meshCutter_;
+
+        //- per cc-cc vector the index of the surface hit
+        labelIOList surfaceIndex_;
+
+        //- user supplied face based data.
+        List<Tuple2<mapType, labelList> > userFaceData_;
+
+
+    // Private Member Functions
+
+        //- Reorder list according to map.
+        template<class T>
+        static void updateList
+        (
+            const labelList& newToOld,
+            const T& nullValue,
+            List<T>& elems
+        );
+
+        //- Add patchfield of given type to all fields on mesh
+        template<class GeoField>
+        static void addPatchFields(fvMesh&, const word& patchFieldType);
+
+        //- Reorder patchfields of all fields on mesh
+        template<class GeoField>
+        static void reorderPatchFields(fvMesh&, const labelList& oldToNew);
+
+        //- Find out which faces have changed given cells (old mesh labels)
+        //  that were marked for refinement.
+        static labelList getChangedFaces
+        (
+            const mapPolyMesh&,
+            const labelList& oldCellsToRefine
+        );
+
+        //- Calculate coupled boundary end vector and refinement level
+        void calcNeighbourData
+        (
+            labelList& neiLevel,
+            pointField& neiCc
+        ) const;
+        ////- Calculate coupled boundary end vector and refinement level. Sort
+        ////  so both sides of coupled face have data in same order.
+        //void calcCanonicalBoundaryData
+        //(
+        //    labelList& leftLevel,
+        //    pointField& leftCc,
+        //    labelList& rightLevel,
+        //    pointField& rightCc
+        //)  const;
+
+        //- Find any intersection of surface. Store in surfaceIndex_.
+        void updateIntersections(const labelList& changedFaces);
+
+        //- Set instance of all local IOobjects
+        void setInstance(const fileName&);
+
+        //- Remove cells. Put exposedFaces into exposedPatchIDs.
+        autoPtr<mapPolyMesh> doRemoveCells
+        (
+            const labelList& cellsToRemove,
+            const labelList& exposedFaces,
+            const labelList& exposedPatchIDs,
+            removeCells& cellRemover
+        );
+
+        // Get cells which are inside any closed surface. Note that
+        // all closed surfaces
+        // will have already been oriented to have keepPoint outside.
+        labelList getInsideCells(const word&) const;
+
+        // Do all to remove inside cells
+        autoPtr<mapPolyMesh> removeInsideCells
+        (
+            const string& msg,
+            const label exposedPatchI
+        );
+
+        //- Used in decomposeCombineRegions. Given global region per cell
+        //  determines master processor/cell for regions straddling
+        //  procboundaries.
+        void getRegionMaster
+        (
+            const boolList& blockedFace,
+            const regionSplit& globalRegion,
+            Map<label>& regionToMaster
+        ) const;
+
+
+        // Refinement candidate selection
+
+            //- Mark cell for refinement (if not already marked). Return false
+            // if refinelimit hit. Keeps running count (in nRefine) of cells
+            // marked for refinement
+            static bool markForRefine
+            (
+                const label markValue,
+                const label nAllowRefine,
+                label& cellValue,
+                label& nRefine
+            );
+
+            //- Calculate list of cells to refine based on intersection of
+            //  features.
+            label markFeatureRefinement
+            (
+                const point& keepPoint,
+                const PtrList<featureEdgeMesh>& featureMeshes,
+                const labelList& featureLevels,
+                const label nAllowRefine,
+
+                labelList& refineCell,
+                label& nRefine
+            ) const;
+
+            //- Mark cells for refinement-shells based refinement.
+            label markInternalRefinement
+            (
+                const PtrList<searchableSurface>&,
+                const labelList& shellLevels,
+                const boolList& shellRefineInside,
+                const label nAllowRefine,
+                labelList& refineCell,
+                label& nRefine
+            ) const;
+
+            //- Mark cells for surface intersection based refinement.
+            label markSurfaceRefinement
+            (
+                const label nAllowRefine,
+                const labelList& neiLevel,
+                const pointField& neiCc,
+                labelList& refineCell,
+                label& nRefine
+            ) const;
+
+            //- Mark cell if local curvature > curvature or
+            //  markDifferingRegions = true and intersections with different
+            //  regions.
+            bool checkCurvature
+            (
+                const labelList& globalToPatch,
+                const scalar curvature,
+                const bool markDifferingRegions,
+                const label nAllowRefine,
+                const label newSurfI,
+                const label newTriI,
+                const label cellI,
+
+                label& maxCellSurfI,
+                label& maxCellTriI,
+
+                labelList& refineCell,
+                label& nRefine
+            ) const;
+
+
+            //- Mark cells for surface curvature based refinement. Marks if
+            //  local curvature > curvature or if on different regions
+            //  (markDifferingRegions)
+            label markSurfaceCurvatureRefinement
+            (
+                const labelList& globalToPatch,
+                const scalar curvature,
+                const bool markDifferingRegions,
+                const label nAllowRefine,
+                const labelList& neiLevel,
+                const pointField& neiCc,
+                labelList& refineCell,
+                label& nRefine
+            ) const;
+
+        // Baffle handling
+
+            //- Determine patches for baffles
+            void getBafflePatches
+            (
+                const labelList& globalToPatch,
+                const labelList& neiLevel,
+                const pointField& neiCc,
+                labelList& ownPatch,
+                labelList& neiPatch
+            ) const;
+
+            //- Determine patch for baffle using some heuristic (and not
+            //  surface)
+            label getBafflePatch
+            (
+                const labelList& facePatch,
+                const label faceI
+            ) const;
+
+            //- Repatches external face or creates baffle for internal face
+            //  with user specified patches (might be different for both sides).
+            //  Returns label of added face.
+            label createBaffle
+            (
+                const label faceI,
+                const label ownPatch,
+                const label neiPatch,
+                polyTopoChange& meshMod
+            ) const;
+
+            //- Helper function to mark face as being on 'boundary'. Used by
+            //  markFacesOnProblemCells
+            void markBoundaryFace
+            (
+                const label faceI,
+                boolList& isBoundaryFace,
+                boolList& isBoundaryEdge,
+                boolList& isBoundaryPoint
+            ) const;
+
+            //- Returns list with for every internal face -1 or the patch
+            //  they should be baffled into.
+            labelList markFacesOnProblemCells
+            (
+                const labelList& globalToPatch
+            ) const;
+
+
+        // Baffle merging
+
+            //- Extract those baffles (duplicate) faces that are on the edge
+            //  of a baffle region. These are candidates for merging.
+            List<labelPair> filterDuplicateFaces(const List<labelPair>&) const;
+
+
+        // Zone handling
+
+            //- Finds zone per cell for cells inside closed named surfaces.
+            //  (uses geometric test for insideness)
+            void findCellZoneGeometric
+            (
+                const labelList& closedNamedSurfaces,
+                const labelList& namedSurfaceIndex,
+                const labelList& surfaceToCellZone,
+                labelList& cellToZone
+            ) const;
+
+            //- Finds zone per cell. Uses topological walk with all faces
+            //  marked in namedSurfaceIndex regarded as blocked.
+            void findCellZoneTopo
+            (
+                const point& keepPoint,
+                const labelList& namedSurfaceIndex,
+                const labelList& surfaceToCellZone,
+                labelList& cellToZone
+            ) const;
+
+
+        //- Disallow default bitwise copy construct
+        meshRefinement(const meshRefinement&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const meshRefinement&);
+
+public:
+
+    //- Runtime type information
+    ClassName("meshRefinement");
+
+
+    // Constructors
+
+        //- Construct from components
+        meshRefinement
+        (
+            fvMesh& mesh,
+            const scalar tol,
+            const refinementSurfaces&
+        );
+
+
+    // Member Functions
+
+        // Access
+
+            //- reference to mesh
+            const fvMesh& mesh() const
+            {
+                return mesh_;
+            }
+            fvMesh& mesh()
+            {
+                return mesh_;
+            }
+
+            //- reference to surface search engines
+            const refinementSurfaces& surfaces() const
+            {
+                return surfaces_;
+            }
+
+            //- reference to meshcutting engine
+            const hexRef8& meshCutter() const
+            {
+                return meshCutter_;
+            }
+
+            //- per start-end edge the index of the surface hit
+            const labelList& surfaceIndex() const
+            {
+                return surfaceIndex_;
+            }
+
+            labelList& surfaceIndex()
+            {
+                return surfaceIndex_;
+            }
+
+            //- Additional face data that is maintained across
+            //  topo changes. Every entry is a list over all faces.
+            //  Bit of a hack. Additional flag to say whether to maintain master
+            //  only (false) or increase set to account for face-from-face.
+            const List<Tuple2<mapType, labelList> >& userFaceData() const
+            {
+                return userFaceData_;
+            }
+
+            List<Tuple2<mapType, labelList> >& userFaceData()
+            {
+                return userFaceData_;
+            }
+
+
+        // Other
+
+            //- Count number of intersections (local)
+            label countHits() const;
+
+            //- Helper function to get decomposition such that all connected
+            //  regions get moved onto one processor. Used to prevent baffles
+            //  straddling processor boundaries. explicitConnections is to
+            //  keep pairs of non-coupled boundary faces together
+            //  (e.g. to keep baffles together)
+            labelList decomposeCombineRegions
+            (
+                const boolList& blockedFace,
+                const List<labelPair>& explicitConnections,
+                decompositionMethod&
+            ) const;
+
+            //- Get faces with intersection.
+            labelList intersectedFaces() const;
+
+            //- Get points on surfaces with intersection and boundary faces.
+            labelList intersectedPoints() const;
+
+            //- Get added patches (inverse of globalToPatch)
+            static labelList addedPatches(const labelList& globalToPatch);
+
+            //- Create patch from set of patches
+            static autoPtr<indirectPrimitivePatch> makePatch
+            (
+                const polyMesh&,
+                const labelList&
+            );
+
+            //- Helper function to make a pointVectorField with correct
+            //  bcs for mesh movement:
+            //  - adaptPatchIDs         : fixedValue
+            //  - processor             : calculated (so free to move)
+            //  - cyclic/wedge/symmetry : slip
+            //  - other                 : slip
+            static tmp<pointVectorField> makeDisplacementField
+            (
+                const pointMesh& pMesh,
+                const labelList& adaptPatchIDs
+            );
+
+        // Refinement
+
+            //- Calculate list of cells to refine.
+            labelList refineCandidates
+            (
+                const point& keepPoint,
+                const labelList& globalToPatch,
+                const scalar curvature,
+
+                const PtrList<featureEdgeMesh>& featureMeshes,
+                const labelList& featureLevels,
+
+                const PtrList<searchableSurface>&,
+                const labelList& shellLevels,
+                const boolList& shellRefineInside,
+
+                const bool featureRefinement,
+                const bool internalRefinement,
+                const bool surfaceRefinement,
+                const bool curvatureRefinement,
+                const label maxGlobalCells,
+                const label maxLocalCells
+            ) const;
+
+            //- Refine some cells
+            autoPtr<mapPolyMesh> refine(const labelList& cellsToRefine);
+
+            //- Refine some cells and rebalance
+            autoPtr<mapDistributePolyMesh> refineAndBalance
+            (
+                const string& msg,
+                decompositionMethod& decomposer,
+                fvMeshDistribute& distributor,
+                const labelList& cellsToRefine
+            );
+
+
+        // Baffle handling
+
+            //- Split off unreachable areas of mesh.
+            void baffleAndSplitMesh
+            (
+                const bool handleSnapProblems,
+                const bool mergeFreeStanding,
+                Time& runTime,
+                const labelList& globalToPatch,
+                const point& keepPoint
+            );
+
+            //- Split off (with optional buffer layers) unreachable areas
+            //  of mesh. Does not introduce baffles.
+            autoPtr<mapPolyMesh> splitMesh
+            (
+                const label nBufferLayers,
+                const labelList& globalToPatch,
+                const point& keepPoint
+            );
+
+            //- Find boundary points that connect to more than one cell
+            //  region and split them.
+            autoPtr<mapPolyMesh> dupNonManifoldPoints();
+
+            //- Create baffle for every internal face where ownPatch != -1
+            autoPtr<mapPolyMesh> createBaffles
+            (
+                const labelList& ownPatch,
+                const labelList& neiPatch
+            );
+
+            //- Return a list of coupled face pairs, i.e. faces that
+            //  use the same vertices.
+            List<labelPair> getDuplicateFaces(const labelList& testFaces) const;
+
+            //- Merge baffles. Gets pairs of faces.
+            autoPtr<mapPolyMesh> mergeBaffles(const List<labelPair>&);
+
+            //- Put faces/cells into zones according to surface specification.
+            //  Returns null if no zone surfaces present. Region containing
+            //  the keepPoint will not be put into a cellZone.
+            autoPtr<mapPolyMesh> zonify(const point& keepPoint);
+
+
+        // Other topo changes
+
+            //- Helper function to add patch to mesh
+            static label addPatch(fvMesh&, const word& name, const word& type);
+
+            //- Split mesh. Keep part containing point.
+            autoPtr<mapPolyMesh> splitMeshRegions(const point& keepPoint);
+
+            //- Update local numbering for mesh redistribution
+            void distribute(const mapDistributePolyMesh&);
+
+            //- Update for external change to mesh. changedFaces are in new mesh
+            //  face labels.
+            void updateMesh(const mapPolyMesh&, const labelList& changedFaces);
+
+
+            // Restoring : is where other processes delete and reinsert data.
+
+                //- Signal points/face/cells for which to store data
+                void storeData
+                (
+                    const labelList& pointsToStore,
+                    const labelList& facesToStore,
+                    const labelList& cellsToStore
+                );
+
+                //- Update local numbering + undo
+                //  Data to restore given as new pointlabel + stored pointlabel
+                //  (i.e. what was in pointsToStore)
+                void updateMesh
+                (
+                    const mapPolyMesh&,
+                    const labelList& changedFaces,
+                    const Map<label>& pointsToRestore,
+                    const Map<label>& facesToRestore,
+                    const Map<label>& cellsToRestore
+                );
+
+
+            //- Merge faces on the same patch (usually from exposing refinement)
+            //  Returns global number of faces merged.
+            label mergePatchFaces
+            (
+                const scalar minCos,
+                const scalar concaveCos
+                //const dictionary& motionDict,
+                //const labelList& globalToPatch
+            );
+
+            //- Remove points not used by any face or points used
+            // by only two faces where the edges are in line
+            autoPtr<mapPolyMesh> mergeEdges(const scalar minCos);
+
+
+        // Debug/IO
+
+            //- Debugging: check that all faces still obey start()>end()
+            void checkData();
+
+            //- Compare two lists over all boundary faces
+            template<class T>
+            void testSyncBoundaryFaceList
+            (
+                const scalar tol,
+                const string&,
+                const UList<T>&,
+                const UList<T>&
+            ) const;
+
+            //- Print some mesh stats.
+            void printMeshInfo(const bool, const string&) const;
+
+            //- Write mesh and all data
+            bool write() const;
+
+            //- Write refinement level as volScalarFields for postprocessing
+            void dumpRefinementLevel() const;
+
+            //- Debug: Write intersection information to OBJ format
+            void dumpIntersections(const fileName& prefix) const;
+
+            //- Do any one of above IO functions. flag is combination of
+            //  writeFlag values.
+            void write(const label flag, const fileName&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "meshRefinementTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/autoMesh/autoHexMesh/meshRefinement/meshRefinementMerge.C

@@ -0,0 +1,230 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    meshRefinement
+
+\*----------------------------------------------------------------------------*/
+
+#include "meshRefinement.H"
+#include "combineFaces.H"
+#include "polyTopoChange.H"
+#include "removePoints.H"
+
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+// Merge faces that are in-line.
+Foam::label Foam::meshRefinement::mergePatchFaces
+(
+    const scalar minCos,
+    const scalar concaveCos
+)
+{
+    // Patch face merging engine
+    combineFaces faceCombiner(mesh_);
+
+    // Get all sets of faces that can be merged
+    labelListList mergeSets(faceCombiner.getMergeSets(minCos, concaveCos));
+
+    label nFaceSets = returnReduce(mergeSets.size(), sumOp<label>());
+
+    Info<< "mergePatchFaces : Merging " << nFaceSets
+        << " sets of faces." << endl;
+
+    if (nFaceSets > 0)
+    {
+        // Topology changes container
+        polyTopoChange meshMod(mesh_);
+
+        // Merge all faces of a set into the first face of the set. Remove
+        // unused points.
+        faceCombiner.setRefinement(mergeSets, meshMod);
+
+        // Change the mesh (no inflation)
+        autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
+
+        // Update fields
+        mesh_.updateMesh(map);
+
+        // Move mesh (since morphing does not do this)
+        if (map().hasMotionPoints())
+        {
+            mesh_.movePoints(map().preMotionPoints());
+        }
+        else
+        {
+            // Delete mesh volumes. No other way to do this?
+            mesh_.clearOut();
+        }
+
+        faceCombiner.updateMesh(map);
+
+        // Get the kept faces that need to be recalculated.
+        // Merging two boundary faces might shift the cell centre
+        // (unless the faces are absolutely planar)
+        labelHashSet retestFaces(6*mergeSets.size());
+
+        forAll(mergeSets, setI)
+        {
+            label oldMasterI = mergeSets[setI][0];
+
+            label faceI = map().reverseFaceMap()[oldMasterI];
+
+            // faceI is always uncoupled boundary face
+            const cell& cFaces = mesh_.cells()[mesh_.faceOwner()[faceI]];
+
+            forAll(cFaces, i)
+            {
+                retestFaces.insert(cFaces[i]);
+            }
+        }
+        updateMesh(map, retestFaces.toc());
+    }
+
+
+    return nFaceSets;
+}
+
+
+// Remove points not used by any face or points used by only two faces where
+// the edges are in line
+Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::mergeEdges
+(
+    const scalar minCos
+)
+{
+    // Point removal analysis engine
+    removePoints pointRemover(mesh_);
+
+    // Count usage of points
+    boolList pointCanBeDeleted;
+    label nRemove = pointRemover.countPointUsage(minCos, pointCanBeDeleted);
+
+    Info<< "Removing " << nRemove
+        << " straight edge points." << endl;
+
+    autoPtr<mapPolyMesh> map;
+
+    if (nRemove > 0)
+    {
+        // Save my local faces that will change. These changed faces might
+        // cause a shift in the cell centre which needs to be retested.
+        // Have to do this before changing mesh since point will be removed.
+        labelHashSet retestOldFaces(nRemove / Pstream::nProcs());
+
+        {
+            const faceList& faces = mesh_.faces();
+
+            forAll(faces, faceI)
+            {
+                const face& f = faces[faceI];
+
+                forAll(f, fp)
+                {
+                    if (pointCanBeDeleted[f[fp]])
+                    {
+                        retestOldFaces.insert(faceI);
+                        break;
+                    }
+                }
+            }
+        }
+
+        // Topology changes container
+        polyTopoChange meshMod(mesh_);
+
+        pointRemover.setRefinement(pointCanBeDeleted, meshMod);
+
+        // Change the mesh (no inflation)
+        map = meshMod.changeMesh(mesh_, false, true);
+
+        // Update fields
+        mesh_.updateMesh(map);
+
+        // Move mesh (since morphing does not do this)
+        if (map().hasMotionPoints())
+        {
+            mesh_.movePoints(map().preMotionPoints());
+        }
+        else
+        {
+            // Delete mesh volumes. No other way to do this?
+            mesh_.clearOut();
+        }
+
+        pointRemover.updateMesh(map);
+
+        // Get the kept faces that need to be recalculated.
+        labelHashSet retestFaces(6*retestOldFaces.size());
+
+        const cellList& cells = mesh_.cells();
+
+        forAllConstIter(labelHashSet, retestOldFaces, iter)
+        {
+            label faceI = map().reverseFaceMap()[iter.key()];
+
+            const cell& ownFaces = cells[mesh_.faceOwner()[faceI]];
+
+            forAll(ownFaces, i)
+            {
+                retestFaces.insert(ownFaces[i]);
+            }
+
+            if (mesh_.isInternalFace(faceI))
+            {
+                const cell& neiFaces = cells[mesh_.faceNeighbour()[faceI]];
+
+                forAll(neiFaces, i)
+                {
+                    retestFaces.insert(neiFaces[i]);
+                }
+            }
+        }
+        updateMesh(map, retestFaces.toc());
+    }
+
+    return map;
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+
+// ************************************************************************* //

src/autoMesh/autoHexMesh/meshRefinement/meshRefinementTemplates.C

@@ -0,0 +1,195 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "meshRefinement.H"
+#include "fvMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+// Add a T entry
+template<class T> void meshRefinement::updateList
+(
+    const labelList& newToOld,
+    const T& nullValue,
+    List<T>& elems
+)
+{
+    List<T> newElems(newToOld.size(), nullValue);
+
+    forAll(newElems, i)
+    {
+        label oldI = newToOld[i];
+
+        if (oldI >= 0)
+        {
+            newElems[i] = elems[oldI];
+        }
+    }
+
+    elems.transfer(newElems);
+}
+
+
+// Compare two lists over all boundary faces
+template<class T>
+void meshRefinement::testSyncBoundaryFaceList
+(
+    const scalar tol,
+    const string& msg,
+    const UList<T>& faceData,
+    const UList<T>& syncedFaceData
+) const
+{
+    label nBFaces = mesh_.nFaces() - mesh_.nInternalFaces();
+
+    if (faceData.size() != nBFaces || syncedFaceData.size() != nBFaces)
+    {
+        FatalErrorIn
+        (
+            "meshRefinement::testSyncBoundaryFaceList"
+            "(const scalar, const string&, const List<T>&, const List<T>&)"
+        )   << "Boundary faces:" << nBFaces
+            << " faceData:" << faceData.size()
+            << " syncedFaceData:" << syncedFaceData.size()
+            << abort(FatalError);
+    }
+
+    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        label bFaceI = pp.start() - mesh_.nInternalFaces();
+
+        forAll(pp, i)
+        {
+            const T& data = faceData[bFaceI];
+            const T& syncData = syncedFaceData[bFaceI];
+
+            if (mag(data - syncData) > tol)
+            {
+                label faceI = pp.start()+i;
+
+                FatalErrorIn("testSyncFaces")
+                    << msg
+                    << "patchFace:" << i
+                    << " face:" << faceI
+                    << " fc:" << mesh_.faceCentres()[faceI]
+                    << " patch:" << pp.name()
+                    << " faceData:" << data
+                    << " syncedFaceData:" << syncData
+                    << " diff:" << mag(data - syncData)
+                    << abort(FatalError);
+            }
+
+            bFaceI++;
+        }
+    }
+}
+
+
+//template <class T, class Mesh>
+template<class GeoField>
+void meshRefinement::addPatchFields(fvMesh& mesh, const word& patchFieldType)
+{
+    HashTable<const GeoField*> flds
+    (
+        mesh.objectRegistry::lookupClass<GeoField>()
+    );
+
+    for
+    (
+        typename HashTable<const GeoField*>::const_iterator iter = flds.begin();
+        iter != flds.end();
+        ++iter
+    )
+    {
+        const GeoField& fld = *iter();
+
+        typename GeoField::GeometricBoundaryField& bfld =
+            const_cast<typename GeoField::GeometricBoundaryField&>
+            (
+                fld.boundaryField()
+            );
+
+        label sz = bfld.size();
+        bfld.setSize(sz+1);
+        bfld.set
+        (
+            sz,
+            GeoField::PatchFieldType::New
+            (
+                patchFieldType,
+                mesh.boundary()[sz],
+                fld.dimensionedInternalField()
+            )
+        );
+    }
+}
+
+
+// Reorder patch field
+template<class GeoField>
+void meshRefinement::reorderPatchFields(fvMesh& mesh, const labelList& oldToNew)
+{
+    HashTable<const GeoField*> flds
+    (
+        mesh.objectRegistry::lookupClass<GeoField>()
+    );
+
+    for
+    (
+        typename HashTable<const GeoField*>::const_iterator iter = flds.begin();
+        iter != flds.end();
+        ++iter
+    )
+    {
+        const GeoField& fld = *iter();
+
+        typename GeoField::GeometricBoundaryField& bfld =
+            const_cast<typename GeoField::GeometricBoundaryField&>
+            (
+                fld.boundaryField()
+            );
+
+        bfld.reorder(oldToNew);
+    }
+}
+
+
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/autoMesh/autoHexMesh/offsetTriSurfaceMesh/offsetTriSurfaceMesh.C

@@ -0,0 +1,199 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "offsetTriSurfaceMesh.H"
+#include "Random.H"
+#include "addToRunTimeSelectionTable.H"
+#include "triSurfaceTools.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(offsetTriSurfaceMesh, 0);
+addToRunTimeSelectionTable(searchableSurface, offsetTriSurfaceMesh, dict);
+
+}
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::offsetTriSurfaceMesh::offsetTriSurfaceMesh
+(
+    const IOobject& io,
+    const triSurface& s,
+    const scalar offset)
+:
+    triSurfaceMesh(io, s),
+    offset_(offset)
+{}
+
+
+Foam::offsetTriSurfaceMesh::offsetTriSurfaceMesh
+(
+    const IOobject& io,
+    const scalar offset
+)
+:
+    triSurfaceMesh(io),
+    offset_(offset)
+{}
+
+
+Foam::offsetTriSurfaceMesh::offsetTriSurfaceMesh
+(
+    const word& name,
+    const objectRegistry& obj,
+    const dictionary& dict
+)
+:
+    triSurfaceMesh(name, obj, dict),
+    offset_(readScalar(dict.lookup("offset")))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::offsetTriSurfaceMesh::~offsetTriSurfaceMesh()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::pointIndexHit Foam::offsetTriSurfaceMesh::findNearest
+(
+    const point& sample,
+    const scalar nearestDistSqr
+) const
+{
+    // Find nearest (add offset to search span)
+    pointIndexHit surfNearest = triSurfaceMesh::findNearest
+    (
+        sample,
+        nearestDistSqr + Foam::sqr(offset_)
+    );
+
+    // Shift back onto surface
+    if (surfNearest.hit())
+    {
+        vector n(sample-surfNearest.hitPoint());
+        n /= mag(n)+VSMALL;
+        surfNearest.setPoint(surfNearest.hitPoint() + offset_*n);
+    }
+    return surfNearest;
+}
+
+
+Foam::pointIndexHit Foam::offsetTriSurfaceMesh::findNearestOnEdge
+(
+    const point& sample,
+    const scalar nearestDistSqr
+) const
+{
+    // Find nearest (add offset to search span)
+    pointIndexHit surfNearest = triSurfaceMesh::findNearestOnEdge
+    (
+        sample,
+        nearestDistSqr + Foam::sqr(offset_)
+    );
+
+    // Shift back onto surface
+    if (surfNearest.hit())
+    {
+        vector n = sample-surfNearest.hitPoint();
+        n /= mag(n)+VSMALL;
+        surfNearest.setPoint(surfNearest.hitPoint() + offset_*n);
+    }
+    return surfNearest;
+}
+
+
+Foam::searchableSurface::volumeType Foam::offsetTriSurfaceMesh::getVolumeType
+(
+    const point& sample
+) const
+{
+    // Find the nearest point on background surface
+    pointIndexHit surfNearest = triSurfaceMesh::findNearest
+    (
+        sample,
+        Foam::sqr(GREAT)
+    );
+
+    if (!surfNearest.hit())
+    {
+        FatalErrorIn("offsetTriSurfaceMesh::getVolumeType(const point&)")
+            << "treeBb:" << tree().bb()
+            << " sample:" << sample
+            << " surfNearest:" << surfNearest
+            << abort(FatalError);
+    }
+
+    // Offset sample to the point.
+    vector n(surfNearest.hitPoint()-sample);
+    n /= mag(n)+VSMALL;
+
+    triSurfaceTools::sideType t = triSurfaceTools::surfaceSide
+    (
+        *this,
+        sample+offset_*n,
+        surfNearest.index(),
+        surfNearest.hitPoint()
+    );
+
+    if (t == triSurfaceTools::UNKNOWN)
+    {
+        return searchableSurface::UNKNOWN;
+    }
+    else if (t == triSurfaceTools::INSIDE)
+    {
+        return searchableSurface::INSIDE;
+    }
+    else if (t == triSurfaceTools::OUTSIDE)
+    {
+        return searchableSurface::OUTSIDE;
+    }
+    else
+    {
+        FatalErrorIn("offsetTriSurfaceMesh::getVolumeType(const point&)")
+            << "problem" << abort(FatalError);
+        return searchableSurface::UNKNOWN;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+
+// ************************************************************************* //

src/autoMesh/autoHexMesh/offsetTriSurfaceMesh/offsetTriSurfaceMesh.H

@@ -0,0 +1,167 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    offsetTriSurfaceMesh
+
+Description
+    triSurfaceMesh with offset. Used to define refinement boxes as a region
+    within certain distance to the refinement surface.
+    Note: reloads surface.
+
+SourceFiles
+    offsetTriSurfaceMesh.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef offsetTriSurfaceMesh_H
+#define offsetTriSurfaceMesh_H
+
+#include "triSurfaceMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class offsetTriSurfaceMesh Declaration
+\*---------------------------------------------------------------------------*/
+
+class offsetTriSurfaceMesh
+:
+    public triSurfaceMesh
+{
+    // Private data
+
+        const scalar offset_;
+
+public:
+
+    //- Runtime type information
+    TypeName("offsetTriSurfaceMesh");
+
+
+    // Constructors
+
+        //- Construct from triSurface and offset
+        offsetTriSurfaceMesh(const IOobject&, const triSurface&, const scalar);
+
+        //- Construct read and offset
+        offsetTriSurfaceMesh(const IOobject& io, const scalar);
+
+        //- Construct as searchableSurface
+        offsetTriSurfaceMesh
+        (
+            const word& name,
+            const objectRegistry& obj,
+            const dictionary& dict
+        );
+
+
+    // Destructor
+
+        virtual ~offsetTriSurfaceMesh();
+
+
+    // Member Functions
+
+        // searchableSurface implementation
+
+            //- Calculate nearest point on surface. Returns
+            //  - bool : any point found nearer than nearestDistSqr
+            //  - label: relevant index in surface
+            //  - point: actual nearest point found
+            virtual pointIndexHit findNearest
+            (
+                const point& sample,
+                const scalar nearestDistSqr
+            ) const;
+
+            //- Calculate nearest point on edge. Returns
+            //  - bool : any point found nearer than nearestDistSqr
+            //  - label: relevant index in surface
+            //  - point: actual nearest point found
+            virtual pointIndexHit findNearestOnEdge
+            (
+                const point& sample,
+                const scalar nearestDistSqr
+            ) const;
+
+            //- Find nearest to line. Returns
+            //  - bool : any point found?
+            //  - label: relevant index in shapes
+            //  - point: actual nearest point found
+            //  sets:
+            //  - tightest  : bounding box
+            //  - linePoint : corresponding nearest point on line
+            virtual pointIndexHit findNearest
+            (
+                const linePointRef& ln,
+                treeBoundBox& tightest,
+                point& linePoint
+            ) const
+            {
+                notImplemented("offsetTriSurfaceMesh::findNearest(..)");
+                return pointIndexHit();
+            }
+
+            //- Find nearest intersection of line between start and end.
+            virtual pointIndexHit findLine
+            (
+                const point& start,
+                const point& end
+            ) const
+            {
+                notImplemented("offsetTriSurfaceMesh::findLine(..)");
+                return pointIndexHit();
+            }
+
+            //- Find any intersection of line between start and end.
+            virtual pointIndexHit findLineAny
+            (
+                const point& start,
+                const point& end
+            ) const
+            {
+                notImplemented("offsetTriSurfaceMesh::findLine(..)");
+                return pointIndexHit();
+            }
+
+            //- Determine type (inside/outside/mixed) for point. unknown if
+            //  cannot be determined (e.g. non-manifold surface)
+            virtual volumeType getVolumeType(const point&) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C

@@ -0,0 +1,503 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    refinementSurfaces
+
+\*----------------------------------------------------------------------------*/
+
+#include "refinementSurfaces.H"
+#include "orientedSurface.H"
+#include "Time.H"
+#include "searchableSurface.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::fileNameList Foam::refinementSurfaces::extractFileNames
+(
+    const PtrList<dictionary>& surfaceDicts
+)
+{
+    fileNameList surfaceFileNames(surfaceDicts.size());
+
+    forAll(surfaceDicts, i)
+    {
+        surfaceFileNames[i] = fileName(surfaceDicts[i].lookup("file"));
+    }
+    return surfaceFileNames;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+Foam::refinementSurfaces::refinementSurfaces
+(
+    const IOobject& io,
+    const PtrList<dictionary>& surfaceDicts
+)
+:
+    triSurfaceMeshes(io, extractFileNames(surfaceDicts)),
+    names_(surfaceDicts.size()),
+    closed_(surfaceDicts.size(), true),
+    faceZoneNames_(surfaceDicts.size()),
+    cellZoneNames_(surfaceDicts.size()),
+    zoneInside_(surfaceDicts.size()),
+    regionOffset_(surfaceDicts.size())
+{
+    labelList globalMinLevel(surfaceDicts.size(), 0);
+    labelList globalMaxLevel(surfaceDicts.size(), 0);
+    List<HashTable<label> > regionMinLevel(surfaceDicts.size());
+    List<HashTable<label> > regionMaxLevel(surfaceDicts.size());
+
+    forAll(surfaceDicts, i)
+    {
+        const dictionary& dict = surfaceDicts[i];
+
+        names_[i] = word(dict.lookup("name"));
+
+        globalMinLevel[i] = readLabel(dict.lookup("minRefinementLevel"));
+        globalMaxLevel[i] = readLabel(dict.lookup("maxRefinementLevel"));
+
+        if
+        (
+            globalMinLevel[i] < 0
+         || globalMaxLevel[i] < globalMinLevel[i]
+        )
+        {
+            FatalErrorIn
+            (
+                "refinementSurfaces::refinementSurfaces"
+                "(const IOobject&, const PtrList<dictionary>&)"
+            )   << "Illegal level specification for surface " << names_[i]
+                << " : minLevel:" << globalMinLevel[i]
+                << " maxLevel:" << globalMaxLevel[i]
+                << exit(FatalError);
+        }
+
+        if (dict.found("faceZone"))
+        {
+            dict.lookup("faceZone") >> faceZoneNames_[i];
+            dict.lookup("cellZone") >> cellZoneNames_[i];
+            dict.lookup("zoneInside") >> zoneInside_[i];
+        }
+
+        if (dict.found("regions"))
+        {
+            PtrList<dictionary> regionDicts(dict.lookup("regions"));
+
+            forAll(regionDicts, dictI)
+            {
+                const dictionary& regionDict = regionDicts[dictI];
+
+                const word regionName(regionDict.lookup("name"));
+                label min = readLabel(regionDict.lookup("minRefinementLevel"));
+                label max = readLabel(regionDict.lookup("maxRefinementLevel"));
+
+                regionMinLevel[i].insert(regionName, min);
+                regionMaxLevel[i].insert(regionName, max);
+            }
+        }
+    }
+
+    // Calculate closedness
+    forAll(closed_, surfI)
+    {
+        const triSurface& s = operator[](surfI);
+
+        const labelListList& edgeFaces = s.edgeFaces();
+
+        forAll(edgeFaces, edgeI)
+        {
+            if (edgeFaces[edgeI].size() != 2)
+            {
+                closed_[surfI] = false;
+                break;
+            }
+        }
+    }
+
+
+    // Calculate local to global region offset
+    label nRegions = 0;
+
+    forAll(surfaceDicts, surfI)
+    {
+        regionOffset_[surfI] = nRegions;
+
+        nRegions += operator[](surfI).patches().size();
+    }
+
+    // From global region number to refinement level
+    minLevel_.setSize(nRegions);
+    minLevel_ = 0;
+    maxLevel_.setSize(nRegions);
+    maxLevel_ = 0;
+
+    forAll(surfaceDicts, surfI)
+    {
+        const geometricSurfacePatchList& regions = operator[](surfI).patches();
+
+        forAll(regions, i)
+        {
+            minLevel_[regionOffset_[surfI] + i] = globalMinLevel[surfI];
+            maxLevel_[regionOffset_[surfI] + i] = globalMaxLevel[surfI];
+        }
+
+        forAllConstIter(HashTable<label>, regionMinLevel[surfI], iter)
+        {
+            // Find the patch
+            forAll(regions, regionI)
+            {
+                if (regions[regionI].name() == iter.key())
+                {
+                    label globalRegionI = regionOffset_[surfI] + regionI;
+
+                    minLevel_[globalRegionI] = iter();
+                    maxLevel_[globalRegionI] =
+                        regionMaxLevel[surfI][iter.key()];
+                    break;
+                }
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+//- Get indices of named surfaces (surfaces with cellZoneNam)
+Foam::labelList Foam::refinementSurfaces::getNamedSurfaces() const
+{
+   labelList namedSurfaces(cellZoneNames_.size());
+
+    label namedI = 0;
+    forAll(cellZoneNames_, surfI)
+    {
+        if (cellZoneNames_[surfI].size() > 0)
+        {
+            namedSurfaces[namedI++] = surfI;
+        }
+    }
+    namedSurfaces.setSize(namedI);
+
+    return namedSurfaces;
+}
+
+
+bool Foam::refinementSurfaces::isSurfaceClosed(const triSurface& s)
+{
+    if (s.nEdges() == s.nInternalEdges())
+    {
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// Orient surface (if they're closed) before any searching is done.
+void Foam::refinementSurfaces::orientSurface
+(
+    const point& keepPoint,
+    triSurfaceMesh& s
+)
+{
+    // Flip surface so keepPoint is outside.
+    bool anyFlipped = orientedSurface::orient(s, keepPoint, true);
+
+    if (anyFlipped)
+    {
+        // orientedSurface will have done a clearOut of the surface.
+        // we could do a clearout of the triSurfaceMeshes::trees()
+        // but these aren't affected by orientation (except for cached
+        // sideness which should not be set at this point. !!Should
+        // check!)
+
+        Info<< "orientSurfaces : Flipped orientation of surface "
+            << s.searchableSurface::name()
+            << " so point " << keepPoint << " is outside." << endl;
+    }
+}
+
+
+// Count number of triangles per surface region
+Foam::labelList Foam::refinementSurfaces::countRegions(const triSurface& s)
+{
+    const geometricSurfacePatchList& regions = s.patches();
+
+    labelList nTris(regions.size(), 0);
+
+    forAll(s, triI)
+    {
+        nTris[s[triI].region()]++;
+    }
+    return nTris;
+}
+
+
+// Find intersections of edge. Return -1 or first surface with higher minLevel
+// number.
+Foam::label Foam::refinementSurfaces::findHigherIntersection
+(
+    const point& start,
+    const point& end,
+    const label currentLevel,   // current cell refinement level
+
+    pointIndexHit& hit
+) const
+{
+    forAll(*this, surfI)
+    {
+        hit = operator[](surfI).findLineAny(start, end);
+
+        if (hit.hit())
+        {
+            if (currentLevel == -1)
+            {
+                // Return any.
+                return surfI;
+            }
+            else
+            {
+                //const triSurface& s = trees()[surfI].shapes().surface();
+                //label region = globalRegion(surfI, s[hit.index()].region());
+                //if (minLevel_[region] > currentLevel)
+                if (minLevelFields_[surfI][hit.index()] > currentLevel)
+                {
+                    return surfI;
+                }
+            }
+        }
+    }
+
+    hit.setMiss();
+
+    return -1;
+}
+
+
+// Find intersection with max of edge. Return -1 or the surface
+// with the highest maxLevel above currentLevel
+Foam::label Foam::refinementSurfaces::findHighestIntersection
+(
+    const point& start,
+    const point& end,
+    const label currentLevel,   // current cell refinement level
+
+    pointIndexHit& maxHit
+) const
+{
+    // surface with highest maxlevel
+    label maxSurface = -1;
+    // maxLevel of maxSurface
+    label maxLevel = currentLevel;
+
+    forAll(*this, surfI)
+    {
+        pointIndexHit hit = operator[](surfI).findLineAny(start, end);
+
+        if (hit.hit())
+        {
+            const triSurface& s = operator[](surfI);
+
+            label region = globalRegion(surfI, s[hit.index()].region());
+
+            if (maxLevel_[region] > maxLevel)
+            {
+                maxSurface = surfI;
+                maxLevel = maxLevel_[region];
+                maxHit = hit;
+            }
+        }
+    }
+
+    if (maxSurface == -1)
+    {
+        // maxLevel unchanged. No interesting surface hit.
+        maxHit.setMiss();
+    }
+
+    return maxSurface;
+}
+
+
+Foam::label Foam::refinementSurfaces::insideZone
+(
+    const labelList& surfaces,
+    const point& pt
+) const
+{
+    forAll(surfaces, i)
+    {
+        label surfI = surfaces[i];
+
+        if (closed_[surfI])
+        {
+            const triSurfaceMesh& s = operator[](surfI);
+
+            triSurfaceMesh::volumeType t = s.getVolumeType(pt);
+
+            if
+            (
+                (t == triSurfaceMesh::INSIDE && zoneInside_[surfI])
+             || (t == triSurfaceMesh::OUTSIDE && !zoneInside_[surfI])
+            )
+            {
+                return surfI;
+            }
+        }
+    }
+    return -1;
+}
+
+
+Foam::label Foam::refinementSurfaces::markInsidePoints
+(
+    const pointField& points,
+    PackedList<1>& isInside
+) const
+{
+    isInside.setSize(points.size());
+    isInside = 0u;
+
+    label nPointInside = 0;
+
+    forAll(points, pointI)
+    {
+        forAll(*this, surfI)
+        {
+            if (closed()[surfI])
+            {
+                const triSurfaceMesh& s = operator[](surfI);
+
+                triSurfaceMesh::volumeType t = s.getVolumeType(points[pointI]);
+
+                if (t == triSurfaceMesh::INSIDE)
+                {
+                    isInside.set(pointI, 1u);
+                    nPointInside++;
+                    break;
+                }
+            }
+        }
+    }
+
+    return nPointInside;
+}
+
+
+void Foam::refinementSurfaces::setMinLevelFields
+(
+    const PtrList<searchableSurface>& shells,
+    const labelList& shellLevels,
+    const boolList& shellRefineInside
+)
+{
+    typedef Foam::DimensionedField<label, triSurfaceGeoMesh>
+        triSurfaceLabelField;
+
+    minLevelFields_.setSize(size());
+
+    forAll(*this, surfI)
+    {
+        const triSurfaceMesh& surfMesh = operator[](surfI);
+
+        minLevelFields_.set
+        (
+            surfI,
+            new triSurfaceLabelField
+            (
+                IOobject
+                (
+                    surfMesh.IOobject::name(),
+                    surfMesh.time().constant(), // directory
+                    "triSurface",               // instance
+                    surfMesh,
+                    IOobject::NO_READ,
+                    IOobject::AUTO_WRITE,
+                    false
+                ),
+                surfMesh,
+                dimless
+            )
+        );
+    }
+
+    // Initialise fields to region wise minLevel
+    forAll(*this, surfI)
+    {
+        const triSurface& s = operator[](surfI);
+
+        triSurfaceLabelField& minLevelField = minLevelFields_[surfI];
+
+        forAll(s, i)
+        {
+            minLevelField[i] = minLevel_[globalRegion(surfI, s[i].region())];
+        }
+    }
+
+    // Adapt for triangles inside shells
+    forAll(*this, surfI)
+    {
+        const triSurface& s = operator[](surfI);
+        triSurfaceLabelField& minLevelField = minLevelFields_[surfI];
+
+        forAll(s, i)
+        {
+            point fc = s[i].centre(s.points());
+
+            forAll(shells, shellI)
+            {
+                // Which side of shell is to be refined
+                searchableSurface::volumeType refSide =
+                (
+                    shellRefineInside[shellI]
+                  ? searchableSurface::INSIDE
+                  : searchableSurface::OUTSIDE
+                );
+
+                // Find whether point is inside or outside shell
+                searchableSurface::volumeType t =
+                    shells[shellI].getVolumeType(fc);
+
+                if (t == refSide)
+                {
+                    minLevelField[i] = max
+                    (
+                        minLevelField[i],
+                        shellLevels[shellI]
+                    );
+                }
+            }
+        }
+    }
+}
+
+
+// ************************************************************************* //

src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H

@@ -0,0 +1,279 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    refinementSurfaces
+
+Description
+    Container for triSurfaces used for surface-driven refinement.
+    These contain all the data about the level of refinement needed per
+    surface.
+
+SourceFiles
+    refinementSurfaces.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef refinementSurfaces_H
+#define refinementSurfaces_H
+
+#include "triSurfaceMeshes.H"
+#include "PackedList.H"
+#include "triSurfaceGeoMesh.H"
+#include "triSurfaceFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class searchableSurface;
+
+/*---------------------------------------------------------------------------*\
+                           Class refinementSurfaces Declaration
+\*---------------------------------------------------------------------------*/
+
+class refinementSurfaces
+:
+    public triSurfaceMeshes
+{
+    // Private data
+
+        //- Name (word)
+        wordList names_;
+
+        //- Per surface whether is closed
+        boolList closed_;
+
+        //- Per 'interface' surface : name of faceZone to put faces into
+        wordList faceZoneNames_;
+
+        //- Per 'interface' surface : name of cellZone to put cells into
+        wordList cellZoneNames_;
+
+        //- Per 'interface' surface : (only used if surface is closed)
+        //  whether to zone cells inside or outside surface.
+        boolList zoneInside_;
+
+
+        //- From local region number to global region number
+        labelList regionOffset_;
+
+        //- From global region number to refinement level
+        labelList minLevel_;
+
+        //- From global region number to refinement level
+        labelList maxLevel_;
+
+        //- Per surface refinement level adapted for shells.
+        PtrList<triSurfaceLabelField> minLevelFields_;
+
+
+    // Private Member Functions
+
+        static fileNameList extractFileNames(const PtrList<dictionary>&);
+
+        //- Disallow default bitwise copy construct
+        refinementSurfaces(const refinementSurfaces&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const refinementSurfaces&);
+
+
+public:
+
+    // Constructors
+
+        //- Construct from directory (io.instance()) and dictionaries
+        refinementSurfaces
+        (
+            const IOobject& io,
+            const PtrList<dictionary>&
+        );
+
+
+    // Member Functions
+
+        // Access
+
+            //- Names of surfaces
+            const wordList& names() const
+            {
+                return names_;
+            }
+
+            //- Per surface whether is closed
+            const boolList& closed() const
+            {
+                return closed_;
+            }
+
+            //- Per 'interface' surface : name of faceZone to put faces into
+            const wordList& faceZoneNames() const
+            {
+                return faceZoneNames_;
+            }
+
+            //- Per 'interface' surface : name of cellZone to put cells into
+            const wordList& cellZoneNames() const
+            {
+                return cellZoneNames_;
+            }
+
+            //- Per 'interface' surface : if closed: zonify cells inside surface
+            const boolList& zoneInside() const
+            {
+                return zoneInside_;
+            }
+
+            //- Get indices of named surfaces (surfaces with cellZoneName)
+            labelList getNamedSurfaces() const;
+
+
+            //- From local region number to global region number
+            const labelList& regionOffset() const
+            {
+                return regionOffset_;
+            }
+
+            //- From global region number to refinement level
+            const labelList& minLevel() const
+            {
+                return minLevel_;
+            }
+
+            //- From global region number to refinement level
+            const labelList& maxLevel() const
+            {
+                return maxLevel_;
+            }
+
+
+        // Helper
+
+            //- From surface and region on surface to global region
+            label globalRegion(const label surfI, const label regionI) const
+            {
+                return regionOffset_[surfI]+regionI;
+            }
+
+            //- From triangle on surface to global region
+            label triangleRegion(const label surfI, const label triI) const
+            {
+                const triSurface& s = operator[](surfI);
+
+                return globalRegion(surfI, s[triI].region());
+            }
+
+            //- Min level for surface and region on surface
+            label minLevel(const label surfI, const label regionI) const
+            {
+                return minLevel_[globalRegion(surfI, regionI)];
+            }
+
+            //- Max level for surface and region on surface
+            label maxLevel(const label surfI, const label regionI) const
+            {
+                return maxLevel_[globalRegion(surfI, regionI)];
+            }
+
+            label nRegions() const
+            {
+                return minLevel_.size();
+            }
+
+            //- Minlevel updated for refinement shells
+            const triSurfaceLabelField& minLevelField(const label surfI) const
+            {
+                return minLevelFields_[surfI];
+            }
+
+            //- Calculate minLevelFields
+            void setMinLevelFields
+            (
+                const PtrList<searchableSurface>& shells,
+                const labelList& shellLevels,
+                const boolList& shellRefineInside
+            );
+
+            //- Helper: is surface closed?
+            static bool isSurfaceClosed(const triSurface&);
+
+            //- Helper: orient (closed only) surfaces so keepPoint is outside.
+            static void orientSurface
+            (
+                const point& keepPoint,
+                triSurfaceMesh& surface
+            );
+
+            //- Helper: count number of triangles per region
+            static labelList countRegions(const triSurface&);
+
+
+        // Searching
+
+            //- Find intersection of edge. Return -1 or first surface
+            //  with higher (than currentLevel) minlevel.
+            //  Return surface number and hit info.
+            label findHigherIntersection
+            (
+                const point& start,
+                const point& end,
+                const label currentLevel,   // current cell refinement level
+                pointIndexHit&
+            ) const;
+
+            //- Find intersection with max level. Return -1 or the surface
+            //  with the highest maxLevel above currentLevel
+            label findHighestIntersection
+            (
+                const point& start,
+                const point& end,
+                const label currentLevel,   // current cell refinement level
+                pointIndexHit&
+            ) const;
+
+            //- Detect if a point is 'inside' (depending on zoneInside flag) a
+            //  zoneable surface. Returns -1 if not, returns first surface it
+            //  is.
+            label insideZone(const labelList& surfaces, const point& pt) const;
+
+
+            //- Mark for all points whether it is inside any closed surface
+            //  Return number of inside points.
+            label markInsidePoints(const pointField&, PackedList<1>& isInside)
+             const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/autoMesh/autoHexMesh/trackedParticle/ExactParticle.C

@@ -0,0 +1,263 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "ExactParticle.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class ParticleType>
+template<class TrackingData>
+Foam::label Foam::ExactParticle<ParticleType>::track
+(
+    const vector& endPosition,
+    TrackingData& td
+)
+{
+    this->facei_ = -1;
+
+    // Tracks to endPosition or stop on boundary
+    while(!this->onBoundary() && this->stepFraction_ < 1.0 - SMALL)
+    {
+        this->stepFraction_ +=
+            trackToFace(endPosition, td)
+           *(1.0 - this->stepFraction_);
+    }
+
+    return this->facei_;
+}
+
+
+template<class ParticleType>
+Foam::label Foam::ExactParticle<ParticleType>::track
+(
+    const vector& endPosition
+)
+{
+    int dummyTd;
+    return track(endPosition, dummyTd);
+}
+
+
+template<class ParticleType>
+template<class TrackingData>
+Foam::scalar Foam::ExactParticle<ParticleType>::trackToFace
+(
+    const vector& endPosition,
+    TrackingData& td
+)
+{
+    const polyMesh& mesh = this->cloud().pMesh();
+    const labelList& cFaces = mesh.cells()[this->celli_];
+
+    point intersection(vector::zero);
+    scalar trackFraction = VGREAT;
+    label hitFacei = -1;
+
+    const vector vec = endPosition-this->position_;
+
+    forAll(cFaces, i)
+    {
+        label facei = cFaces[i];
+
+        if (facei != this->face())
+        {
+            pointHit inter = mesh.faces()[facei].intersection
+            (
+                this->position_,
+                vec,
+                mesh.faceCentres()[facei],
+                mesh.points(),
+                intersection::HALF_RAY
+            );
+
+            if (inter.hit() && inter.distance() < trackFraction)
+            {
+                trackFraction = inter.distance();
+                hitFacei = facei;
+                intersection = inter.hitPoint();
+            }
+        }
+    }
+
+
+    if (hitFacei != -1)
+    {
+        if (trackFraction > 1.0)
+        {
+            // Nearest intersection beyond endPosition so we hit endPosition.
+            this->position_ = endPosition;
+            this->facei_ = -1;
+            return 1.0;
+        }
+        else
+        {
+            this->position_ = intersection;
+            this->facei_ = hitFacei;
+        }
+    }
+    else
+    {
+        // Did not find any intersection. Fall back to original approximate
+        // algorithm
+        return Particle<ParticleType>::trackToFace
+        (
+            endPosition,
+            td
+        );
+    }
+
+
+    // Normal situation (trackFraction 0..1)
+
+    bool internalFace = this->cloud().internalFace(this->facei_);
+
+    // change cell
+    if (internalFace) // Internal face
+    {
+        if (this->celli_ == mesh.faceOwner()[this->facei_])
+        {
+            this->celli_ = mesh.faceNeighbour()[this->facei_];
+        }
+        else if (this->celli_ == mesh.faceNeighbour()[this->facei_])
+        {
+            this->celli_ = mesh.faceOwner()[this->facei_];
+        }
+        else
+        {
+            FatalErrorIn
+            (
+                "ExactParticle::trackToFace"
+                "(const vector&, TrackingData&)"
+            )<< "addressing failure" << nl
+             << abort(FatalError);
+        }
+    }
+    else
+    {
+        ParticleType& p = static_cast<ParticleType&>(*this);
+
+        // Soft-sphere algorithm ignores the boundary
+        if (p.softImpact())
+        {
+            trackFraction = 1.0;
+            this->position_ = endPosition;
+        }
+
+        label patchi = patch(this->facei_);
+        const polyPatch& patch = mesh.boundaryMesh()[patchi];
+
+        if (isA<wedgePolyPatch>(patch))
+        {
+            p.hitWedgePatch
+            (
+                static_cast<const wedgePolyPatch&>(patch), td
+            );
+        }
+        else if (isA<symmetryPolyPatch>(patch))
+        {
+            p.hitSymmetryPatch
+            (
+                static_cast<const symmetryPolyPatch&>(patch), td
+            );
+        }
+        else if (isA<cyclicPolyPatch>(patch))
+        {
+            p.hitCyclicPatch
+            (
+                static_cast<const cyclicPolyPatch&>(patch), td
+            );
+        }
+        else if (isA<processorPolyPatch>(patch))
+        {
+            p.hitProcessorPatch
+            (
+                static_cast<const processorPolyPatch&>(patch), td
+            );
+        }
+        else if (isA<wallPolyPatch>(patch))
+        {
+            p.hitWallPatch
+            (
+                static_cast<const wallPolyPatch&>(patch), td
+            );
+        }
+        else if (isA<polyPatch>(patch))
+        {
+            p.hitPatch
+            (
+                static_cast<const polyPatch&>(patch), td
+            );
+        }
+        else
+        {
+            FatalErrorIn
+            (
+                "ExactParticle::trackToFace"
+                "(const vector& endPosition, scalar& trackFraction)"
+            )<< "patch type " << patch.type() << " not suported" << nl
+             << abort(FatalError);
+        }
+    }
+
+    // If the trackFraction = 0 something went wrong.
+    // Either the particle is flipping back and forth across a face perhaps
+    // due to velocity interpolation errors or it is in a "hole" in the mesh
+    // caused by face warpage.
+    // In both cases resolve the positional ambiguity by moving the particle
+    // slightly towards the cell-centre.
+    if (trackFraction < SMALL)
+    {
+        this->position_ +=
+            1.0e-6*(mesh.cellCentres()[this->celli_] - this->position_);
+    }
+
+    return trackFraction;
+}
+
+
+template<class ParticleType>
+Foam::scalar Foam::ExactParticle<ParticleType>::trackToFace
+(
+    const vector& endPosition
+)
+{
+    int dummyTd;
+    return trackToFace(endPosition, dummyTd);
+}
+
+
+template<class ParticleType>
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const ExactParticle<ParticleType>& p
+)
+{
+    return operator<<(os, static_cast<Particle<ParticleType> >(p));
+}
+
+
+// ************************************************************************* //

src/autoMesh/autoHexMesh/trackedParticle/ExactParticle.H

@@ -0,0 +1,193 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    ExactParticle
+
+Description
+    Special version of Particle to do tracking on non-convex cells.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef ExactParticle_H
+#define ExactParticle_H
+
+#include "face.H"
+#include "Particle.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+template<class ExactParticle>
+class Cloud;
+
+// Forward declaration of friend functions and operators
+
+template<class ParticleType>
+class ExactParticle;
+
+template<class ParticleType>
+Ostream& operator<<
+(
+    Ostream&,
+    const ExactParticle<ParticleType>&
+);
+
+
+/*---------------------------------------------------------------------------*\
+                           Class ExactParticle Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class ParticleType>
+class ExactParticle
+:
+    public Particle<ParticleType>
+{
+
+public:
+
+    friend class Cloud<ParticleType>;
+
+
+    // Constructors
+
+        //- Construct from components
+        ExactParticle
+        (
+            const Cloud<ParticleType>& cloud,
+            const vector& position,
+            const label celli
+        )
+        :
+            Particle<ParticleType>(cloud, position, celli)
+        {}
+            
+
+        //- Construct from Istream
+        ExactParticle
+        (
+            const Cloud<ParticleType>& cloud,
+            Istream& is,
+            bool readFields = true
+        )
+        :
+            Particle<ParticleType>(cloud, is, readFields)
+        {}
+
+
+        //- Factory class to read-construct particles used for parallel transfer
+        class iNew
+        {
+
+            // Private data
+
+            const Cloud<ParticleType>& cloud_;
+
+
+        public:
+
+            iNew(const Cloud<ParticleType>& cloud)
+            :
+                cloud_(cloud)
+            {}
+
+            autoPtr<ParticleType> operator()(Istream& is) const
+            {
+                return autoPtr<ParticleType>
+                (
+                    new ParticleType(cloud_, is)
+                );
+            }
+        };
+
+
+    // Destructor
+
+        virtual ~ExactParticle()
+        {}
+
+
+    // Member Functions
+
+            //- Track particle to end of trajectory
+            //  or until it hits the boundary.
+            //  On entry 'stepFraction()' should be set to the fraction of the
+            //  time-step at which the tracking starts and on exit it contains
+            //  the fraction of the time-step completed.
+            //  Returns the boundary face index if the track stops at the
+            //  boundary, -1 otherwise.
+            template<class TrackingData>
+            label track
+            (
+                const vector& endPosition,
+                TrackingData& td
+            );
+
+            //- Calls the templated track with dummy TrackingData
+            label track(const vector& endPosition);
+
+            //- Track particle to a given position and returns 1.0 if the
+            //  trajectory is completed without hitting a face otherwise
+            //  stops at the face and returns the fraction of the trajectory
+            //  completed.
+            //  on entry 'stepFraction()' should be set to the fraction of the
+            //  time-step at which the tracking starts.
+            template<class TrackingData>
+            scalar trackToFace
+            (
+                const vector& endPosition,
+                TrackingData& td
+            );
+
+            //- Calls the templated trackToFace with dummy TrackingData
+            scalar trackToFace(const vector& endPosition);
+
+
+    // Ostream Operator
+
+        friend Ostream& operator<< <ParticleType>
+        (
+            Ostream&,
+            const ExactParticle<ParticleType>&
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "ExactParticle.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/autoMesh/autoHexMesh/trackedParticle/Particle.H

@@ -0,0 +1,431 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Particle
+
+Description
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Particle_H
+#define Particle_H
+
+#include "vector.H"
+#include "IDLList.H"
+#include "labelList.H"
+#include "pointField.H"
+#include "typeInfo.H"
+#include "OFstream.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+template<class Particle>
+class Cloud;
+
+class wedgePolyPatch;
+class symmetryPolyPatch;
+class cyclicPolyPatch;
+class processorPolyPatch;
+class wallPolyPatch;
+class polyPatch;
+
+// Forward declaration of friend functions and operators
+
+template<class ParticleType>
+class Particle;
+
+template<class ParticleType>
+Ostream& operator<<
+(
+    Ostream&,
+    const Particle<ParticleType>&
+);
+
+
+/*---------------------------------------------------------------------------*\
+                           Class Particle Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class ParticleType>
+class Particle
+:
+    public IDLList<ParticleType>::link
+{
+
+protected:
+
+        //- Reference to the particle cloud
+        const Cloud<ParticleType>& cloud_;
+
+        //- Position of particle
+        vector position_;
+
+        //- Index of the cell it is in
+        label celli_;
+
+        //- Face index if the particle is on a face otherwise -1
+        label facei_;
+
+        //- Fraction of time-step completed
+        scalar stepFraction_;
+
+
+    // Private member functions
+
+        //- Return the 'lambda' value for the position, p, on the face,
+        // where, p = from + lamda*(to - from)
+        // for non-static meshes
+        inline scalar lambda
+        (
+            const vector& from,
+            const vector& to,
+            const label facei,
+            const scalar stepFraction
+        );
+
+        //- Return the 'lambda' value for the position, p, on the face,
+        // where, p = from + lamda*(to - from)
+        // for static meshes
+        inline scalar lambda
+        (
+            const vector& from,
+            const vector& to,
+            const label facei
+        );
+
+        //- Return the faces between position and cell centre
+        labelList findFaces
+        (
+            const vector& position
+        );
+
+        //- Return the faces between position and cell centre
+        labelList findFaces
+        (
+            const vector& position,
+            const label celli,
+            const scalar stepFraction
+        );
+
+
+        //- Overridable function to handle the particle hitting a wedgePatch
+        template<class TrackingData>
+        void hitWedgePatch
+        (
+            const wedgePolyPatch&,
+            TrackingData& td
+        );
+
+        //- Overridable function to handle the particle hitting a symmetryPatch
+        template<class TrackingData>
+        void hitSymmetryPatch
+        (
+            const symmetryPolyPatch&,
+            TrackingData& td
+        );
+
+        //- Overridable function to handle the particle hitting a cyclicPatch
+        template<class TrackingData>
+        void hitCyclicPatch
+        (
+            const cyclicPolyPatch&,
+            TrackingData& td
+        );
+
+        //- Overridable function to handle the particle hitting a processorPatch
+        template<class TrackingData>
+        void hitProcessorPatch
+        (
+            const processorPolyPatch&,
+            TrackingData& td
+        );
+
+        //- Overridable function to handle the particle hitting a wallPatch
+        template<class TrackingData>
+        void hitWallPatch
+        (
+            const wallPolyPatch&,
+            TrackingData& td
+        );
+
+        //- Overridable function to handle the particle hitting a general patch
+        template<class TrackingData>
+        void hitPatch
+        (
+            const polyPatch&,
+            TrackingData& td
+        );
+
+        //- Transform the position the particle
+        //  according to the given transformation tensor
+        void transformPosition(const tensor& T);
+
+        //- Transform the physical properties of the particle
+        //  according to the given transformation tensor
+        void transformProperties(const tensor& T);
+
+        //- Transform the physical properties of the particle
+        //  according to the given separation vector
+        void transformProperties(const vector& separation);
+
+        //- Convert global addressing to the processor patch
+        //  local equivalents
+        template<class TrackingData>
+        void prepareForParallelTransfer(const label patchi, TrackingData& td);
+
+        //- Convert processor patch addressing to the global equivalents
+        //  and set the celli to the face-neighbour
+        template<class TrackingData>
+        void correctAfterParallelTransfer(const label patchi, TrackingData& td);
+
+
+public:
+
+    friend class Cloud<ParticleType>;
+
+
+    //- Runtime type information
+    TypeName("Particle");
+
+
+    //- Class used to pass tracking data to the trackToFace function
+    class trackData
+    {
+
+        // Private data
+
+            //- Reference to the cloud containing this particle
+            Cloud<ParticleType>& cloud_;
+
+
+    public:
+
+        bool switchProcessor;
+        bool keepParticle;
+
+
+        // Constructors
+
+            inline trackData
+            (
+                Cloud<ParticleType>& cloud
+            );
+
+
+        // Member functions
+
+            inline Cloud<ParticleType>& cloud();
+    };
+
+
+    // Constructors
+
+        //- Construct from components
+        Particle
+        (
+            const Cloud<ParticleType>&,
+            const vector& position,
+            const label celli
+        );
+
+        //- Construct from Istream
+        Particle
+        (
+            const Cloud<ParticleType>&,
+            Istream&,
+            bool readFields = true
+        );
+
+        //- Factory class to read-construct particles used for parallel transfer
+        class iNew
+        {
+
+            // Private data
+
+            const Cloud<ParticleType>& cloud_;
+
+
+        public:
+
+            iNew(const Cloud<ParticleType>& cloud)
+            :
+                cloud_(cloud)
+            {}
+
+            autoPtr<ParticleType> operator()(Istream& is) const
+            {
+                return autoPtr<ParticleType>(new ParticleType(cloud_, is));
+            }
+        };
+
+
+    // Destructor
+
+        virtual ~Particle()
+        {}
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return true if particle is in cell
+            inline bool inCell();
+
+            //- Return true if position is in cell i
+            inline bool inCell
+            (
+                const vector& position,
+                const label celli,
+                const scalar stepFraction
+            );
+
+            //- Return current particle position
+            inline const vector& position() const;
+
+            //- Return current particle position
+            inline vector& position();
+
+            //- Return current cell particle is in
+            inline label cell() const;
+
+            //- Return current face particle is on otherwise -1
+            inline label face() const;
+
+            //- Return reference to the particle cloud
+            inline const Cloud<ParticleType>& cloud() const;
+
+            //- Return the impact model to be used, soft or hard (default).
+            inline bool softImpact() const;
+
+
+        // Check
+
+            //- Is the particle on the boundary/(or outside the domain)?
+            inline bool onBoundary() const;
+
+            //- Which patch is particle on
+            inline label patch(const label facei) const;
+
+            //- Which face of this patch is this particle on
+            inline label patchFace(const label patchi, const label facei) const;
+
+            //- The nearest distance to a wall that
+            //  the particle can be in the n direction
+            inline scalar wallImpactDistance(const vector& n) const;
+
+            //- Return the fraction of time-step completed
+            inline scalar& stepFraction();
+
+            //-  Return the fraction of time-step completed
+            inline scalar stepFraction() const;
+
+
+        // Track
+
+            //- Track particle to end of trajectory
+            //  or until it hits the boundary.
+            //  On entry 'stepFraction()' should be set to the fraction of the
+            //  time-step at which the tracking starts and on exit it contains
+            //  the fraction of the time-step completed.
+            //  Returns the boundary face index if the track stops at the
+            //  boundary, -1 otherwise.
+            template<class TrackingData>
+            label track
+            (
+                const vector& endPosition,
+                TrackingData& td
+            );
+
+            //- Calls the templated track with dummy TrackingData
+            label track(const vector& endPosition);
+
+            //- Track particle to a given position and returns 1.0 if the
+            //  trajectory is completed without hitting a face otherwise
+            //  stops at the face and returns the fraction of the trajectory
+            //  completed.
+            //  on entry 'stepFraction()' should be set to the fraction of the
+            //  time-step at which the tracking starts.
+            template<class TrackingData>
+            scalar trackToFace
+            (
+                const vector& endPosition,
+                TrackingData& td
+            );
+
+            //- Calls the templated trackToFace with dummy TrackingData
+            scalar trackToFace(const vector& endPosition);
+
+            //- Return the index of the face to be used in the interpolation routine
+            inline label faceInterpolation() const;
+
+
+    // I-O
+
+        //- Write the fields associated with the owner cloud
+        static void writeFields
+        (
+            const Cloud<ParticleType>& c
+        );
+
+
+    // Ostream Operator
+
+        friend Ostream& operator<< <ParticleType>
+        (
+            Ostream&,
+            const Particle<ParticleType>&
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "ParticleI.H"
+
+#define defineParticleTypeNameAndDebug(Type, DebugSwitch)                     \
+    template<>                                                                \
+    const Foam::word Particle<Type>::typeName(#Type);                         \
+    template<>                                                                \
+    int Particle<Type>::debug(Foam::debug::debugSwitch(#Type, DebugSwitch));
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "Particle.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/autoMesh/autoHexMesh/trackedParticle/trackedParticle.C

@@ -0,0 +1,261 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*----------------------------------------------------------------------------*/
+
+#include "trackedParticle.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+//- Construct from components
+Foam::trackedParticle::trackedParticle
+(
+    const Cloud<trackedParticle>& c,
+    const vector& position,
+    const label celli,
+    const point& end,
+    const label level,
+    const label i,
+    const label j
+)
+:
+    ExactParticle<trackedParticle>(c, position, celli),
+    end_(end),
+    level_(level),
+    i_(i),
+    j_(j)
+{}
+
+
+//- Construct from Istream
+Foam::trackedParticle::trackedParticle
+(
+    const Cloud<trackedParticle>& c,
+    Istream& is,
+    bool readFields
+)
+:
+    ExactParticle<trackedParticle>(c, is)
+{
+    if (readFields)
+    {
+        if (is.format() == IOstream::ASCII)
+        {
+            is >> end_;
+            level_ = readLabel(is);
+            i_ = readLabel(is);
+            j_ = readLabel(is);
+        }
+        else
+        {
+            is.read
+            (
+                reinterpret_cast<char*>(&end_),
+                sizeof(end_) + sizeof(level_) + sizeof(i_) + sizeof(j_)
+            );
+        }
+    }
+
+    // Check state of Istream
+    is.check
+    (
+        "trackedParticle::trackedParticle"
+        "(const Cloud<trackedParticle>&, Istream&, bool)"
+    );
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::trackedParticle::move(trackedParticle::trackData& td)
+{
+    td.switchProcessor = false;
+    td.keepParticle = true;
+
+    scalar deltaT = cloud().pMesh().time().deltaT().value();
+    scalar tEnd = (1.0 - stepFraction())*deltaT;
+    scalar dtMax = tEnd;
+
+    while (td.keepParticle && !td.switchProcessor && tEnd > SMALL)
+    {
+        // set the lagrangian time-step
+        scalar dt = min(dtMax, tEnd);
+
+        // mark visited cell with max level.
+        td.maxLevel()[cell()] = max(td.maxLevel()[cell()], level_);
+
+        dt *= trackToFace(end_, td);
+
+        tEnd -= dt;
+        stepFraction() = 1.0 - tEnd/deltaT;
+    }
+
+    return td.keepParticle;
+}
+
+
+void Foam::trackedParticle::hitWedgePatch
+(
+    const wedgePolyPatch&,
+    trackedParticle::trackData& td
+)
+{
+    // Remove particle
+    td.keepParticle = false;
+}
+
+
+void Foam::trackedParticle::hitWedgePatch
+(
+    const wedgePolyPatch&,
+    int&
+)
+{}
+
+
+void Foam::trackedParticle::hitSymmetryPatch
+(
+    const symmetryPolyPatch&,
+    trackedParticle::trackData& td
+)
+{
+    // Remove particle
+    td.keepParticle = false;
+}
+
+
+void Foam::trackedParticle::hitSymmetryPatch
+(
+    const symmetryPolyPatch&,
+    int&
+)
+{}
+
+
+void Foam::trackedParticle::hitCyclicPatch
+(
+    const cyclicPolyPatch&,
+    trackedParticle::trackData& td
+)
+{
+    // Remove particle
+    td.keepParticle = false;
+}
+
+
+void Foam::trackedParticle::hitCyclicPatch
+(
+    const cyclicPolyPatch&,
+    int&
+)
+{}
+
+
+void Foam::trackedParticle::hitProcessorPatch
+(
+    const processorPolyPatch&,
+    trackedParticle::trackData& td
+)
+{
+    // Remove particle
+    td.switchProcessor = true;
+}
+
+
+void Foam::trackedParticle::hitProcessorPatch
+(
+    const processorPolyPatch&,
+    int&
+)
+{}
+
+
+void Foam::trackedParticle::hitWallPatch
+(
+    const wallPolyPatch& wpp,
+    trackedParticle::trackData& td
+)
+{
+    // Remove particle
+    td.keepParticle = false;
+}
+
+
+void Foam::trackedParticle::hitWallPatch
+(
+    const wallPolyPatch& wpp,
+    int&
+)
+{}
+
+
+void Foam::trackedParticle::hitPatch
+(
+    const polyPatch& wpp,
+    trackedParticle::trackData& td
+)
+{
+    // Remove particle
+    td.keepParticle = false;
+}
+
+
+void Foam::trackedParticle::hitPatch
+(
+    const polyPatch& wpp,
+    int&
+)
+{}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<(Ostream& os, const trackedParticle& p)
+{
+    if (os.format() == IOstream::ASCII)
+    {
+        os  << static_cast<const ExactParticle<trackedParticle>&>(p)
+            << token::SPACE << p.end_
+            << token::SPACE << p.level_
+            << token::SPACE << p.i_
+            << token::SPACE << p.j_;
+    }
+    else
+    {
+        os  << static_cast<const ExactParticle<trackedParticle>&>(p);
+        os.write
+        (
+            reinterpret_cast<const char*>(&p.end_),
+            sizeof(p.end_) + sizeof(p.level_) + sizeof(p.i_) + sizeof(p.j_)
+        );
+    }
+
+    // Check state of Ostream
+    os.check("Ostream& operator<<(Ostream&, const trackedParticle&)");
+
+    return os;
+}
+
+
+// ************************************************************************* //

src/autoMesh/autoHexMesh/trackedParticle/trackedParticle.H

@@ -0,0 +1,276 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    trackedParticle
+
+Description
+    Particle class that marks cells it passes through. Used to mark cells
+    visited by feature edges. Uses ExactParticle tracking class so
+    will work on concave cells.
+
+SourceFiles
+    trackedParticle.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef trackedParticle_H
+#define trackedParticle_H
+
+#include "ExactParticle.H"
+#include "autoPtr.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class trackedParticleCloud;
+
+/*---------------------------------------------------------------------------*\
+                           Class trackedParticle Declaration
+\*---------------------------------------------------------------------------*/
+
+class trackedParticle
+:
+    public ExactParticle<trackedParticle>
+{
+    // Private data
+
+        //- end point to track to
+        point end_;
+
+        //- level of this particle
+        label level_;
+
+        //- passive label
+        label i_;
+
+        //- passive label
+        label j_;
+
+public:
+
+    friend class Cloud<trackedParticle>;
+
+    //- Class used to pass tracking data to the trackToFace function
+    class trackData
+    {
+        //- Reference to the cloud containing this particle
+        Cloud<trackedParticle>& cloud_;
+
+        labelList& maxLevel_;
+
+    public:
+
+        bool switchProcessor;
+        bool keepParticle;
+
+
+        // Constructors
+
+            trackData(Cloud<trackedParticle>& cloud, labelList& maxLevel)
+            :
+                cloud_(cloud),
+                maxLevel_(maxLevel)
+            {}
+
+
+        // Member functions
+
+            Cloud<trackedParticle>& cloud()
+            {
+                return cloud_;
+            }
+
+            labelList& maxLevel()
+            {
+                return maxLevel_;
+            }
+    };
+
+
+
+    // Constructors
+
+        //- Construct from components
+        trackedParticle
+        (
+            const Cloud<trackedParticle>& c,
+            const vector& position,
+            const label celli,
+            const point& end,
+            const label level,
+            const label i,
+            const label j
+        );
+
+        //- Construct from Istream
+        trackedParticle
+        (
+            const Cloud<trackedParticle>& c,
+            Istream& is,
+            bool readFields = true
+        );
+
+        //- Construct and return a clone
+        autoPtr<trackedParticle> clone() const
+        {
+            return autoPtr<trackedParticle>(new trackedParticle(*this));
+        }
+
+
+    // Member Functions
+
+        //- point to track to
+        point& end()
+        {
+            return end_;
+        }
+
+        //- transported label
+        label& i()
+        {
+            return i_;
+        }
+
+        //- transported label
+        label& j()
+        {
+            return j_;
+        }
+
+
+
+        // Tracking
+
+            //- Track all particles to their end point
+            bool move(trackData&);
+
+
+            //- Overridable function to handle the particle hitting a wedge
+            void hitWedgePatch
+            (
+                const wedgePolyPatch&,
+                trackedParticle::trackData& td
+            );
+            void hitWedgePatch
+            (
+                const wedgePolyPatch&,
+                int&
+            );
+
+            //- Overridable function to handle the particle hitting a
+            //  symmetryPlane
+            void hitSymmetryPatch
+            (
+                const symmetryPolyPatch&,
+                trackedParticle::trackData& td
+            );
+            void hitSymmetryPatch
+            (
+                const symmetryPolyPatch&,
+                int&
+            );
+
+            //- Overridable function to handle the particle hitting a cyclic
+            void hitCyclicPatch
+            (
+                const cyclicPolyPatch&,
+                trackedParticle::trackData& td
+            );
+            void hitCyclicPatch
+            (
+                const cyclicPolyPatch&,
+                int&
+            );
+
+            //- Overridable function to handle the particle hitting a
+            //- processorPatch
+            void hitProcessorPatch
+            (
+                const processorPolyPatch&,
+                trackedParticle::trackData& td
+            );
+            void hitProcessorPatch
+            (
+                const processorPolyPatch&,
+                int&
+            );
+
+            //- Overridable function to handle the particle hitting a wallPatch
+            void hitWallPatch
+            (
+                const wallPolyPatch&,
+                trackedParticle::trackData& td
+            );
+            void hitWallPatch
+            (
+                const wallPolyPatch&,
+                int&
+            );
+
+            //- Overridable function to handle the particle hitting a polyPatch
+            void hitPatch
+            (
+                const polyPatch&,
+                trackedParticle::trackData& td
+            );
+            void hitPatch
+            (
+                const polyPatch&,
+                int&
+            );
+
+
+    // Ostream Operator
+
+        friend Ostream& operator<<(Ostream&, const trackedParticle&);
+
+
+};
+
+
+template<>
+inline bool contiguous<trackedParticle>()
+{
+    return true;
+}
+
+//template<>
+//void Cloud<trackedParticle>::readFields();
+//
+//template<>
+//void Cloud<trackedParticle>::writeFields() const;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/autoMesh/autoHexMesh/trackedParticle/trackedParticleCloud.C

@@ -0,0 +1,44 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "trackedParticle.H"
+#include "Cloud.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineParticleTypeNameAndDebug(trackedParticle, 0);
+defineTemplateTypeNameAndDebug(Cloud<trackedParticle>, 0);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //