Merge branch 'master' into dsmc

applications/solvers/heatTransfer/buoyantPisoFoam/Make/options

@@ -1,7 +1,6 @@
 EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
-    -I../XiFoam \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
 EXE_LIBS = \

applications/solvers/heatTransfer/buoyantPisoFoam/buoyantPisoFoam.C

@@ -70,11 +70,12 @@ int main(int argc, char *argv[])
 
         #include "UEqn.H"
 
+        #include "hEqn.H"
+
         // --- PISO loop
 
         for (int corr=0; corr<nCorr; corr++)
         {
-            #include "hEqn.H"
             #include "pEqn.H"
         }
 

applications/solvers/heatTransfer/buoyantPisoFoam/createFields.H

@@ -53,12 +53,15 @@
     );
 
     Info<< "Creating field DpDt\n" << endl;
-    volScalarField DpDt =
-        fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+    volScalarField DpDt
+    (
+        "DpDt",
+        fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p)
+    );
 
     Info<< "Calculating field g.h\n" << endl;
     volScalarField gh("gh", g & mesh.C());
-    surfaceScalarField ghf("gh", g & mesh.Cf());
+    surfaceScalarField ghf("ghf", g & mesh.Cf());
 
     dimensionedScalar pRef("pRef", p.dimensions(), thermo->lookup("pRef"));
 

applications/solvers/heatTransfer/buoyantPisoFoam/pEqn.H

@@ -12,8 +12,8 @@
     (
         fvc::interpolate(rho)
        *(
-           (fvc::interpolate(U) & mesh.Sf())
-         + fvc::ddtPhiCorr(rUA, rho, U, phi)
+            (fvc::interpolate(U) & mesh.Sf())
+          + fvc::ddtPhiCorr(rUA, rho, U, phi)
         )
     );
 

applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H

@@ -53,7 +53,7 @@
 
     Info<< "Calculating field g.h\n" << endl;
     volScalarField gh("gh", g & mesh.C());
-    surfaceScalarField ghf("gh", g & mesh.Cf());
+    surfaceScalarField ghf("ghf", g & mesh.Cf());
 
     dimensionedScalar pRef("pRef", p.dimensions(), thermo->lookup("pRef"));
 

applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.H

@@ -103,8 +103,8 @@ public:
             const dictionary&
         );
 
-        //- Construct by mapping given solidWallMixedTemperatureCoupledFvPatchScalarField
-        //  onto a new patch
+        //- Construct by mapping given
+        //  solidWallMixedTemperatureCoupledFvPatchScalarField onto a new patch
         solidWallMixedTemperatureCoupledFvPatchScalarField
         (
             const solidWallMixedTemperatureCoupledFvPatchScalarField&,

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H

@@ -4,7 +4,7 @@
     volScalarField& rho = rhoFluid[i];
     volScalarField& K = KFluid[i];
     volVectorField& U = UFluid[i];
-    surfaceScalarField phi = phiFluid[i];
+    surfaceScalarField& phi = phiFluid[i];
     compressible::turbulenceModel& turb = turbulence[i];
     volScalarField& DpDt = DpDtFluid[i];
     const volScalarField& gh = ghFluid[i];

applications/solvers/multiphase/interFoam/createFields.H

@@ -85,7 +85,7 @@
 
     Info<< "Calculating field g.h\n" << endl;
     volScalarField gh("gh", g & mesh.C());
-    surfaceScalarField ghf("gh", g & mesh.Cf());
+    surfaceScalarField ghf("ghf", g & mesh.Cf());
 
 
     volScalarField p

applications/solvers/multiphase/multiphaseInterFoam/createFields.H

@@ -47,7 +47,7 @@
 
     Info<< "Calculating field g.h\n" << endl;
     volScalarField gh("gh", g & mesh.C());
-    surfaceScalarField ghf("gh", g & mesh.Cf());
+    surfaceScalarField ghf("ghf", g & mesh.Cf());
 
 
     volScalarField p

applications/solvers/multiphase/settlingFoam/createFields.H

@@ -339,4 +339,4 @@
     );
 
     Info<< "Calculating field (g.h)f\n" << endl;
-    surfaceScalarField ghf = surfaceScalarField("gh", g & mesh.Cf());
+    surfaceScalarField ghf = surfaceScalarField("ghf", g & mesh.Cf());

applications/test/List/ListTest.C

@@ -43,31 +43,33 @@ using namespace Foam;
 
 int main(int argc, char *argv[])
 {
-    List<vector> list(IStringStream("1 ((0 1 2))")());
-    Info<< list << endl;
+    List<vector> list1(IStringStream("1 ((0 1 2))")());
+    Info<< "list1: " << list1 << endl;
 
     List<vector> list2(IStringStream("((0 1 2) (3 4 5) (6 7 8))")());
-    Info<< list2 << endl;
+    Info<< "list2: " << list2 << endl;
+
+    list1.append(list2);
+    Info<< "list1.append(list2): " << list1 << endl;
 
     Info<< findIndex(list2, vector(3, 4, 5)) << endl;
 
     list2.setSize(10, vector(1, 2, 3));
-    Info<< list2 << endl;
+    Info<< "list2: " << list2 << endl;
 
     List<vector> list3(list2.xfer());
     Info<< "Transferred via the xfer() method" << endl;
-    Info<< list2 << nl
-        << list3 << endl;
+    Info<< "list2: " << list2 << nl
+        << "list3: " << list3 << endl;
 
 
     // Subset
     const labelList map(IStringStream("2 (0 2)")());
     List<vector> subList3(list3, map);
     Info<< "Elements " << map << " out of " << list3
-        << " : " << subList3 << endl;
+        << " => " << subList3 << endl;
 
     return 0;
 }
 
-
 // ************************************************************************* //

applications/utilities/mesh/conversion/writeMeshObj/writeMeshObj.C

@@ -344,6 +344,7 @@ int main(int argc, char *argv[])
     argList::validOptions.insert("point", "pointI");
     argList::validOptions.insert("cellSet", "setName");
     argList::validOptions.insert("faceSet", "setName");
+#   include "addRegionOption.H"
 
 #   include "setRootCase.H"
 #   include "createTime.H"
@@ -364,7 +365,7 @@ int main(int argc, char *argv[])
 
     instantList timeDirs = timeSelector::select0(runTime, args);
 
-#   include "createPolyMesh.H"
+#   include "createNamedPolyMesh.H"
 
     forAll(timeDirs, timeI)
     {

applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict

@@ -274,6 +274,14 @@ addLayersControls
 
     // Create buffer region for new layer terminations
     nBufferCellsNoExtrude 0;
+
+
+    // Overall max number of layer addition iterations
+    nLayerIter 50;
+
+    // Max number of iterations after which relaxed meshQuality controls
+    // get used.
+    nRelaxedIter 20;
 }
 
 
@@ -335,6 +343,16 @@ meshQualityControls
     nSmoothScale 4;
     //- amount to scale back displacement at error points
     errorReduction 0.75;
+
+
+
+    // Optional : some meshing phases allow usage of relaxed rules.
+    // See e.g. addLayersControls::nRelaxedIter.
+    relaxed
+    {
+        //- Maximum non-orthogonality allowed. Set to 180 to disable.
+        maxNonOrtho 75;
+    }
 }
 
 

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

@@ -48,12 +48,164 @@ using namespace Foam;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+void createBaffles
+(
+    const polyMesh& mesh,
+    const label faceI,
+    const label oldPatchI,
+    const labelList& newPatches,
+    PackedBoolList& doneFace,
+    polyTopoChange& meshMod
+)
+{
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+    const faceZoneMesh& faceZones = mesh.faceZones();
+
+    const face& f = mesh.faces()[faceI];
+    label zoneID = faceZones.whichZone(faceI);
+    bool zoneFlip = false;
+    if (zoneID >= 0)
+    {
+        const faceZone& fZone = faceZones[zoneID];
+        zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
+    }
+    label nei = -1;
+    if (oldPatchI == -1)
+    {
+        nei = mesh.faceNeighbour()[faceI];
+    }
+
+
+    face revFace(f.reverseFace());
+
+    forAll(newPatches, i)
+    {
+        if (oldPatchI == -1)
+        {
+            // Internal face
+            if (doneFace.set(faceI))
+            {
+                // First usage of face. Modify.
+                meshMod.setAction
+                (
+                    polyModifyFace
+                    (
+                        f,                          // modified face
+                        faceI,                      // label of face
+                        mesh.faceOwner()[faceI],    // owner
+                        -1,                         // neighbour
+                        false,                      // face flip
+                        newPatches[i],              // patch for face
+                        false,                      // remove from zone
+                        zoneID,                     // zone for face
+                        zoneFlip                    // face flip in zone
+                    )
+                );
+            }
+            else
+            {
+                // Second or more usage of face. Add.
+                meshMod.setAction
+                (
+                    polyAddFace
+                    (
+                        f,                          // modified face
+                        mesh.faceOwner()[faceI],    // owner
+                        -1,                         // neighbour
+                        -1,                         // master point
+                        -1,                         // master edge
+                        faceI,                      // master face
+                        false,                      // face flip
+                        newPatches[i],              // patch for face
+                        zoneID,                     // zone for face
+                        zoneFlip                    // face flip in zone
+                    )
+                );
+            }
+
+            meshMod.setAction
+            (
+                polyAddFace
+                (
+                    revFace,                    // modified face
+                    nei,                        // owner
+                    -1,                         // neighbour
+                    -1,                         // masterPointID
+                    -1,                         // masterEdgeID
+                    faceI,                      // masterFaceID,
+                    true,                       // face flip
+                    newPatches[i],              // patch for face
+                    zoneID,                     // zone for face
+                    zoneFlip                    // face flip in zone
+                )
+            );
+        }
+        else if
+        (
+            patches[oldPatchI].coupled()
+         && patches[newPatches[i]].coupled()
+        )
+        {
+            // Do not allow coupled patches to be moved to different coupled
+            // patches.
+            //WarningIn("createBaffles()")
+            //    << "Not moving face from coupled patch "
+            //    << patches[oldPatchI].name()
+            //    << " to another coupled patch " << patches[newPatches[i]]
+            //    << endl;
+        }
+        else
+        {
+            if (doneFace.set(faceI))
+            {
+                meshMod.setAction
+                (
+                    polyModifyFace
+                    (
+                        f,                          // modified face
+                        faceI,                      // label of face
+                        mesh.faceOwner()[faceI],    // owner
+                        -1,                         // neighbour
+                        false,                      // face flip
+                        newPatches[i],              // patch for face
+                        false,                      // remove from zone
+                        zoneID,                     // zone for face
+                        zoneFlip                    // face flip in zone
+                    )
+                );
+            }
+            else
+            {
+                meshMod.setAction
+                (
+                    polyAddFace
+                    (
+                        f,                          // modified face
+                        mesh.faceOwner()[faceI],    // owner
+                        -1,                         // neighbour
+                        -1,                         // master point
+                        -1,                         // master edge
+                        faceI,                      // master face
+                        false,                      // face flip
+                        newPatches[i],              // patch for face
+                        zoneID,                     // zone for face
+                        zoneFlip                    // face flip in zone
+                    )
+                );
+            }
+        }
+    }
+}
+
+
+
 // Main program:
 
 int main(int argc, char *argv[])
 {
     argList::validArgs.append("set");
     argList::validArgs.append("patch");
+    argList::validOptions.insert("additionalPatches", "(patch2 .. patchN)");
     argList::validOptions.insert("overwrite", "");
 
 #   include "setRootCase.H"
@@ -63,29 +215,64 @@ int main(int argc, char *argv[])
     const word oldInstance = mesh.pointsInstance();
 
     const polyBoundaryMesh& patches = mesh.boundaryMesh();
-    const faceZoneMesh& faceZones = mesh.faceZones();
 
     // Faces to baffle
     word setName(args.additionalArgs()[0]);
     Pout<< "Reading faceSet from " << setName << nl << endl;
     faceSet facesToSplit(mesh, setName);
+    // Make sure set is synchronised across couples
+    facesToSplit.sync(mesh);
     Pout<< "Read " << facesToSplit.size() << " faces from " << setName
         << nl << endl;
 
-    // Patch to put them into
+
+    // Patches to put baffles into
+    labelList newPatches(1);
+
     word patchName(args.additionalArgs()[1]);
-    label wantedPatchI = patches.findPatchID(patchName);
+    newPatches[0] = patches.findPatchID(patchName);
+    Pout<< "Using patch " << patchName
+        << " at index " << newPatches[0] << endl;
 
-    Pout<< "Using patch " << patchName << " at index " << wantedPatchI << endl;
-
-    if (wantedPatchI == -1)
+    if (newPatches[0] == -1)
     {
         FatalErrorIn(args.executable())
-            << "Cannot find patch " << patchName << exit(FatalError);
+            << "Cannot find patch " << patchName << endl
+            << "Valid patches are " << patches.names() << exit(FatalError);
     }
 
+
+    // Additional patches
+    if (args.options().found("additionalPatches"))
+    {
+        const wordList patchNames
+        (
+            IStringStream(args.options()["additionalPatches"])()
+        );
+
+        forAll(patchNames, i)
+        {
+            label patchI = patches.findPatchID(patchNames[i]);
+            Info<< "Using additional patch " << patchNames[i]
+                << " at index " << patchI << endl;
+
+            if (patchI == -1)
+            {
+                FatalErrorIn(args.executable())
+                    << "Cannot find patch " << patchNames[i] << endl
+                    << "Valid patches are " << patches.names()
+                    << exit(FatalError);
+            }
+
+            newPatches.append(patchI);
+        }
+    }
+
+
     bool overwrite = args.options().found("overwrite");
 
+
+
     // Read objects in time directory
     IOobjectList objects(mesh, runTime.timeName());
 
@@ -128,107 +315,53 @@ int main(int argc, char *argv[])
     polyTopoChange meshMod(mesh);
 
 
-    // Creating baffles:
-    // - coupled boundary faces : become the patch specified
-    // - non-coupled ,,         : illegal
-    // - internal faces         : converted into boundary faces.
-
-    labelList newPatch(mesh.nFaces(), -1);
-
-    forAllConstIter(faceSet, facesToSplit, iter)
-    {
-        label faceI = iter.key();
-
-        label patchI = patches.whichPatch(faceI);
-
-        if (patchI == -1)
-        {
-            newPatch[faceI] = wantedPatchI;
-        }
-        else
-        {
-            if (patches[patchI].coupled())
-            {
-                if (patchI != wantedPatchI)
-                {
-                    newPatch[faceI] = wantedPatchI;
-                }
-            }
-            else
-            {
-                FatalErrorIn(args.executable())
-                    << "Can only create baffles from internal faces"
-                    << " or coupled boundary faces." << endl
-                    << "Face " << faceI << " is a boundary face on patch "
-                    << patches[patchI].name() << exit(FatalError);
-            }
-        }
-    }
-
-
-    // If one side of a coupled boundary is marked for baffling, make sure to
-    // also do the other side.
-
-    syncTools::syncFaceList(mesh, newPatch, maxEqOp<label>(), false);
-
+    // Do the actual changes
 
     label nBaffled = 0;
+    PackedBoolList doneFace(mesh.nFaces());
 
-    forAll(newPatch, faceI)
+    for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
     {
-        if (newPatch[faceI] != -1)
+        if (facesToSplit.found(faceI))
         {
-            const face& f = mesh.faces()[faceI];
-            label zoneID = faceZones.whichZone(faceI);
-            bool zoneFlip = false;
-            if (zoneID >= 0)
-            {
-                const faceZone& fZone = faceZones[zoneID];
-                zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
-            }
-
-            meshMod.setAction
+            createBaffles
             (
-                polyModifyFace
-                (
-                    f,                          // modified face
-                    faceI,                      // label of face
-                    mesh.faceOwner()[faceI],    // owner
-                    -1,                         // neighbour
-                    false,                      // face flip
-                    newPatch[faceI],            // patch for face
-                    false,                      // remove from zone
-                    zoneID,                     // zone for face
-                    zoneFlip                    // face flip in zone
-                )
+                mesh,
+                faceI,
+                -1,                             // oldPatchI,
+                newPatches,
+                doneFace,
+                meshMod
             );
-
-            if (mesh.isInternalFace(faceI))
-            {
-                meshMod.setAction
-                (
-                    polyAddFace
-                    (
-                        f.reverseFace(),            // modified face
-                        mesh.faceNeighbour()[faceI],// owner
-                        -1,                         // neighbour
-                        -1,                         // masterPointID
-                        -1,                         // masterEdgeID
-                        faceI,                      // masterFaceID,
-                        false,                      // face flip
-                        newPatch[faceI],            // patch for face
-                        zoneID,                     // zone for face
-                        zoneFlip                    // face flip in zone
-                    )
-                );
-            }
-
             nBaffled++;
         }
     }
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        forAll(pp, i)
+        {
+            label faceI = pp.start()+i;
+
+            if (facesToSplit.found(faceI))
+            {
+                createBaffles
+                (
+                    mesh,
+                    faceI,
+                    patchI,                     // oldPatchI,
+                    newPatches,
+                    doneFace,
+                    meshMod
+                );
+                nBaffled++;
+            }
+        }
+    }
 
 
-    Pout<< "Converted locally " << nBaffled
+    Info<< "Converted " << returnReduce(nBaffled, sumOp<label>())
         << " faces into boundary faces on patch " << patchName << nl << endl;
 
     if (!overwrite)
@@ -252,7 +385,7 @@ int main(int argc, char *argv[])
     {
         mesh.setInstance(oldInstance);
     }
-    Pout<< "Writing mesh to " << runTime.timeName() << endl;
+    Info<< "Writing mesh to " << runTime.timeName() << endl;
 
     mesh.write();
 

applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C

@@ -57,6 +57,7 @@ Description
 #include "EdgeMap.H"
 #include "syncTools.H"
 #include "ReadFields.H"
+#include "directMappedWallPolyPatch.H"
 
 using namespace Foam;
 
@@ -1022,13 +1023,13 @@ EdgeMap<label> addRegionPatches
             (
                 mesh,
                 regionNames[e[0]] + "_to_" + regionNames[e[1]],
-                polyPatch::typeName
+                directMappedWallPolyPatch::typeName
             );
             addPatch
             (
                 mesh,
                 regionNames[e[1]] + "_to_" + regionNames[e[0]],
-                polyPatch::typeName
+                directMappedWallPolyPatch::typeName
             );
 
             Info<< "For interface between region " << e[0]
@@ -1100,7 +1101,6 @@ EdgeMap<label> addRegionPatches
 //}
 
 
-//XXXXXXXXX
 // Find region that covers most of cell zone
 label findCorrespondingRegion
 (
@@ -1152,7 +1152,6 @@ label findCorrespondingRegion
 
     return regionI;
 }
-//XXXXXXXXX
 
 
 //// Checks if cellZone has corresponding cellRegion.

applications/utilities/parallelProcessing/decomposePar/decomposePar.C

@@ -83,6 +83,7 @@ Usage
 int main(int argc, char *argv[])
 {
     argList::noParallel();
+#   include "addRegionOption.H"
     argList::validOptions.insert("cellDist", "");
     argList::validOptions.insert("copyUniform", "");
     argList::validOptions.insert("fields", "");
@@ -92,6 +93,17 @@ int main(int argc, char *argv[])
 
 #   include "setRootCase.H"
 
+    word regionName = fvMesh::defaultRegion;
+    word regionDir = word::null;
+
+    if (args.options().found("region"))
+    {
+        regionName = args.options()["region"];
+        regionDir = regionName;
+        Info<< "Decomposing mesh " << regionName << nl << endl;
+    }
+
+
     bool writeCellDist(args.options().found("cellDist"));
     bool copyUniform(args.options().found("copyUniform"));
     bool decomposeFieldsOnly(args.options().found("fields"));
@@ -105,7 +117,17 @@ int main(int argc, char *argv[])
 
     // determine the existing processor count directly
     label nProcs = 0;
-    while (isDir(runTime.path()/(word("processor") + name(nProcs))))
+    while
+    (
+        isDir
+        (
+            runTime.path()
+           /(word("processor") + name(nProcs))
+           /runTime.constant()
+           /regionDir
+           /polyMesh::meshSubDir
+        )
+    )
     {
         ++nProcs;
     }
@@ -119,6 +141,7 @@ int main(int argc, char *argv[])
             (
                 "decomposeParDict",
                 runTime.time().system(),
+                regionDir,          // use region if non-standard
                 runTime,
                 IOobject::MUST_READ,
                 IOobject::NO_WRITE,
@@ -196,7 +219,7 @@ int main(int argc, char *argv[])
     (
         IOobject
         (
-            domainDecomposition::defaultRegion,
+            regionName,
             runTime.timeName(),
             runTime
         )
@@ -219,7 +242,7 @@ int main(int argc, char *argv[])
             (
                 runTime.path()
               / mesh.facesInstance()
-              / polyMesh::defaultRegion
+              / regionName
               / "cellDecomposition"
             );
 
@@ -383,7 +406,12 @@ int main(int argc, char *argv[])
 
             label i = 0;
 
-            forAllIter(Cloud<indexedParticle>, lagrangianPositions[cloudI], iter)
+            forAllIter
+            (
+                Cloud<indexedParticle>,
+                lagrangianPositions[cloudI],
+                iter
+            )
             {
                 iter().index() = i++;
 
@@ -405,7 +433,8 @@ int main(int argc, char *argv[])
 
                 if (!cellParticles[cloudI][celli])
                 {
-                    cellParticles[cloudI][celli] = new SLList<indexedParticle*>();
+                    cellParticles[cloudI][celli] = new SLList<indexedParticle*>
+                    ();
                 }
 
                 cellParticles[cloudI][celli]->append(&iter());
@@ -513,7 +542,7 @@ int main(int argc, char *argv[])
         (
             IOobject
             (
-                fvMesh::defaultRegion,
+                regionName,
                 processorDb.timeName(),
                 processorDb
             )

applications/utilities/parallelProcessing/decomposePar/decomposeParDict

@@ -19,11 +19,12 @@ FoamFile
 
 numberOfSubdomains  4;
 
-// preservePatches (inlet);
+//- Keep owner and neighbour on same processor for faces in zones:
 // preserveFaceZones (heater solid1 solid3);
 
-method          simple;
+method          scotch;
 // method          hierarchical;
+// method          simple;
 // method          metis;
 // method          manual;
 
@@ -53,6 +54,9 @@ metisCoeffs
   */
 }
 
+scotchCoeffs
+{}
+
 manualCoeffs
 {
     dataFile    "decompositionData";

applications/utilities/parallelProcessing/decomposePar/distributeCells.C

@@ -45,35 +45,6 @@ void domainDecomposition::distributeCells()
 
     labelHashSet sameProcFaces;
 
-    if (decompositionDict_.found("preservePatches"))
-    {
-        wordList pNames(decompositionDict_.lookup("preservePatches"));
-
-        Info<< "Keeping owner and neighbour of faces in patches " << pNames
-            << " on same processor" << endl;
-
-        const polyBoundaryMesh& patches = boundaryMesh();
-
-        forAll(pNames, i)
-        {
-            label patchI = patches.findPatchID(pNames[i]);
-
-            if (patchI == -1)
-            {
-                FatalErrorIn("domainDecomposition::distributeCells()")
-                    << "Unknown preservePatch " << pNames[i]
-                    << endl << "Valid patches are " << patches.names()
-                    << exit(FatalError);
-            }
-
-            const polyPatch& pp = patches[patchI];
-
-            forAll(pp, i)
-            {
-                sameProcFaces.insert(pp.start() + i);
-            }
-        }
-    }
     if (decompositionDict_.found("preserveFaceZones"))
     {
         wordList zNames(decompositionDict_.lookup("preserveFaceZones"));
@@ -104,14 +75,6 @@ void domainDecomposition::distributeCells()
         }
     }
 
-    if (sameProcFaces.size())
-    {
-        Info<< "Selected " << sameProcFaces.size()
-            << " faces whose owner and neighbour cell should be kept on the"
-            << " same processor" << endl;
-    }
-
-
 
     // Construct decomposition method and either do decomposition on
     // cell centres or on agglomeration
@@ -129,6 +92,10 @@ void domainDecomposition::distributeCells()
     }
     else
     {
+        Info<< "Selected " << sameProcFaces.size()
+            << " faces whose owner and neighbour cell should be kept on the"
+            << " same processor" << endl;
+
         // Faces where owner and neighbour are not 'connected' (= all except
         // sameProcFaces)
         boolList blockedFace(nFaces(), true);

applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C

@@ -268,7 +268,7 @@ bool domainDecomposition::writeDecomposition()
         (
             IOobject
             (
-                polyMesh::defaultRegion,
+                this->polyMesh::name(),  // region name of undecomposed mesh
                 "constant",
                 processorDb
             ),

applications/utilities/parallelProcessing/reconstructPar/processorMeshes.C

@@ -200,12 +200,6 @@ void Foam::processorMeshes::reconstructPoints(fvMesh& mesh)
     // Read the field for all the processors
     PtrList<pointIOField> procsPoints(meshes_.size());
 
-    fileName regionPrefix = "";
-    if (meshName_ != fvMesh::defaultRegion)
-    {
-        regionPrefix = meshName_;
-    }
-
     forAll (meshes_, procI)
     {
         procsPoints.set
@@ -217,7 +211,7 @@ void Foam::processorMeshes::reconstructPoints(fvMesh& mesh)
                 (
                     "points",
                     meshes_[procI].time().timeName(),
-                    regionPrefix/polyMesh::meshSubDir,
+                    polyMesh::meshSubDir,
                     meshes_[procI],
                     IOobject::MUST_READ,
                     IOobject::NO_WRITE

applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C

@@ -116,6 +116,13 @@ int main(int argc, char *argv[])
     {
         regionPrefix = regionName;
     }
+
+    // Set all times on processor meshes equal to reconstructed mesh
+    forAll (databases, procI)
+    {
+        databases[procI].setTime(runTime.timeName(), runTime.timeIndex());
+    }
+
     // Read all meshes and addressing to reconstructed mesh
     processorMeshes procMeshes(databases, regionName);
 

applications/utilities/surface/surfaceCheck/surfaceCheck.C

@@ -32,6 +32,7 @@ License
 #include "OFstream.H"
 #include "surfaceIntersection.H"
 #include "SortableList.H"
+#include "PatchTools.H"
 
 using namespace Foam;
 
@@ -171,11 +172,11 @@ int main(int argc, char *argv[])
 
     argList::validArgs.clear();
     argList::validArgs.append("surface file");
-    argList::validOptions.insert("noSelfIntersection", "");
+    argList::validOptions.insert("checkSelfIntersection", "");
     argList::validOptions.insert("verbose", "");
     argList args(argc, argv);
 
-    bool checkSelfIntersection = !args.options().found("noSelfIntersection");
+    bool checkSelfIntersection = args.options().found("checkSelfIntersection");
     bool verbose = args.options().found("verbose");
 
     fileName surfFileName(args.additionalArgs()[0]);
@@ -596,13 +597,14 @@ int main(int argc, char *argv[])
     // Check orientation
     // ~~~~~~~~~~~~~~~~~
 
-    boolList borderEdge(surf.checkOrientation(false));
+    labelHashSet borderEdge(surf.size()/1000);
+    PatchTools::checkOrientation(surf, false, &borderEdge);
 
     //
     // Colour all faces into zones using borderEdge
     //
     labelList normalZone;
-    label numNormalZones = surf.markZones(borderEdge, normalZone);
+    label numNormalZones = PatchTools::markZones(surf, borderEdge, normalZone);
 
     Pout<< endl
         << "Number of zones (connected area with consistent normal) : "

applications/utilities/surface/surfaceConvert/surfaceConvert.C

@@ -115,15 +115,13 @@ int main(int argc, char *argv[])
     }
 
     Info<< "writing " << exportName;
-    if (scaleFactor <= 0)
+    if (scaleFactor > 0)
     {
-        Info<< " without scaling" << endl;
-    }
-    else
-    {
-        Info<< " with scaling " << scaleFactor << endl;
+        Info<< " with scaling " << scaleFactor;
         surf.scalePoints(scaleFactor);
     }
+    Info<< endl;
+
     surf.write(exportName, sortByRegion);
 
     Info<< "\nEnd\n" << endl;

applications/utilities/surface/surfaceMeshConvert/surfaceMeshConvert.C

@@ -73,7 +73,7 @@ int main(int argc, char *argv[])
     argList::noParallel();
     argList::validArgs.append("inputFile");
     argList::validArgs.append("outputFile");
-    argList::validOptions.insert("clean",  "scale");
+    argList::validOptions.insert("clean", "");
     argList::validOptions.insert("scaleIn",  "scale");
     argList::validOptions.insert("scaleOut", "scale");
     argList::validOptions.insert("dict", "coordinateSystemsDict");

applications/utilities/surface/surfaceMeshExport/surfaceMeshExport.C

@@ -76,7 +76,7 @@ int main(int argc, char *argv[])
     argList::noParallel();
     argList::validArgs.append("outputFile");
     argList::validOptions.insert("name",  "name");
-    argList::validOptions.insert("clean",  "scale");
+    argList::validOptions.insert("clean", "");
     argList::validOptions.insert("scaleIn",  "scale");
     argList::validOptions.insert("scaleOut", "scale");
     argList::validOptions.insert("dict", "coordinateSystemsDict");

applications/utilities/surface/surfaceMeshImport/surfaceMeshImport.C

@@ -76,7 +76,7 @@ int main(int argc, char *argv[])
     argList::noParallel();
     argList::validArgs.append("inputFile");
     argList::validOptions.insert("name",  "name");
-    argList::validOptions.insert("clean",  "scale");
+    argList::validOptions.insert("clean", "");
     argList::validOptions.insert("scaleIn",  "scale");
     argList::validOptions.insert("scaleOut", "scale");
     argList::validOptions.insert("dict", "coordinateSystemsDict");

etc/controlDict

@@ -360,7 +360,7 @@ DebugSwitches
     diagonal            0;
     dictionary          0;
     dimensionSet        1;
-    directMapped        0;
+    directMappedBase    0;
     directMappedPatch   0;
     directMappedVelocityFlux 0;
     directionMixed      0;

src/OpenFOAM/containers/Lists/List/List.C

@@ -404,6 +404,61 @@ void Foam::List<T>::clear()
 }
 
 
+template<class T>
+void Foam::List<T>::append(const UList<T>& lst)
+{
+    if (this == &lst)
+    {
+        FatalErrorIn
+        (
+            "List<T>::append(const UList<T>&)"
+        )   << "attempted appending to self" << abort(FatalError);
+    }
+
+    label nextFree = this->size_;
+    setSize(nextFree + lst.size());
+
+    forAll(lst, elemI)
+    {
+        this->operator[](nextFree++) = lst[elemI];
+    }
+}
+
+
+template<class T>
+void Foam::List<T>::append(const UIndirectList<T>& lst)
+{
+    label nextFree = this->size_;
+    setSize(nextFree + lst.size());
+
+    forAll(lst, elemI)
+    {
+        this->operator[](nextFree++) = lst[elemI];
+    }
+}
+
+
+template<class T>
+void Foam::List<T>::append(const SLList<T>& lst)
+{
+    if (lst.size())
+    {
+        label nextFree = this->size_;
+        setSize(nextFree + lst.size());
+
+        for
+        (
+            typename SLList<T>::const_iterator iter = lst.begin();
+            iter != lst.end();
+            ++iter
+        )
+        {
+            this->operator[](nextFree++) = iter();
+        }
+    }
+}
+
+
 // Transfer the contents of the argument List into this List
 // and anull the argument list
 template<class T>
@@ -559,12 +614,9 @@ void Foam::List<T>::operator=(const IndirectList<T>& lst)
         if (this->size_) this->v_ = new T[this->size_];
     }
 
-    if (this->size_)
+    forAll(*this, i)
     {
-        forAll(*this, i)
-        {
-            this->operator[](i) = lst[i];
-        }
+        this->operator[](i) = lst[i];
     }
 }
 
@@ -581,12 +633,9 @@ void Foam::List<T>::operator=(const UIndirectList<T>& lst)
         if (this->size_) this->v_ = new T[this->size_];
     }
 
-    if (this->size_)
+    forAll(*this, i)
     {
-        forAll(*this, i)
-        {
-            this->operator[](i) = lst[i];
-        }
+        this->operator[](i) = lst[i];
     }
 }
 
@@ -603,12 +652,9 @@ void Foam::List<T>::operator=(const BiIndirectList<T>& lst)
         if (this->size_) this->v_ = new T[this->size_];
     }
 
-    if (this->size_)
+    forAll(*this, i)
     {
-        forAll(*this, i)
-        {
-            this->operator[](i) = lst[i];
-        }
+        this->operator[](i) = lst[i];
     }
 }
 

src/OpenFOAM/containers/Lists/List/List.H

@@ -109,7 +109,7 @@ public:
         List(const List<T>&);
 
         //- Construct by transferring the parameter contents
-        List(const Xfer<List<T> >&);
+        List(const Xfer< List<T> >&);
 
         //- Construct as copy or re-use as specified.
         List(List<T>&, bool reUse);
@@ -181,6 +181,15 @@ public:
             //- Clear the list, i.e. set size to zero.
             void clear();
 
+            //- Append a List at the end of this list
+            void append(const UList<T>&);
+
+            //- Append a UIndirectList at the end of this list
+            void append(const UIndirectList<T>&);
+
+            //- Append a SLList at the end of this list
+            void append(const SLList<T>&);
+
             //- Transfer the contents of the argument List into this List
             //  and annull the argument list.
             void transfer(List<T>&);
@@ -195,7 +204,7 @@ public:
             void transfer(SortableList<T>&);
 
             //- Transfer contents to the Xfer container
-            inline Xfer<List<T> > xfer();
+            inline Xfer< List<T> > xfer();
 
             //- Return subscript-checked element of UList.
             inline T& newElmt(const label);

src/OpenFOAM/containers/Lists/PackedList/PackedList.H

@@ -26,7 +26,7 @@ Class
     Foam::PackedList
 
 Description
-    A Dynamically allocatable list of packed unsigned ints.
+    A dynamically allocatable list of packed unsigned integers.
 
     The list resizing is similar to DynamicList, thus the methods clear()
     and setSize() behave like their DynamicList counterparts and the methods
@@ -38,7 +38,7 @@ Note
     In a const context, the '[]' operator simply returns the stored value,
     with out-of-range elements returned as zero.
     In a non-const context, the '[]' operator returns an iteratorBase, which
-    may not have a valid reference for out-of-range elements.
+    might not have a valid reference for out-of-range elements.
     The iteratorBase class handles the assignment of new values.
 
     Using the iteratorBase as a proxy allows assignment of values
@@ -50,11 +50,11 @@ Note
     @endcode
 
     Using get() or the '[]' operator are similarly fast. Looping and reading
-    with an iterator is approx. 15% slower, but can be more flexible.
+    via an iterator is approx. 15% slower, but can be more flexible.
 
     Using the set() operator (and the '[]' operator) are marginally slower
-    (approx. 5%) than using an iterator, but the set() method has an
-    advantage that it also returns a bool if the value changed.  This can be
+    (approx. 5%) than using an iterator, but the set() method has the
+    advantage of also returning a bool if the value changed.  This can be
     useful for branching on changed values.
 
     @code
@@ -65,7 +65,7 @@ Note
 
     The lazy evaluation used means that reading an out-of-range element
     returns zero, but does not affect the list size.  Even in a non-const
-    context, only the assigment causes the element to be created.
+    context, only the assigment itself causes the element to be created.
     For example,
     @code
         list.resize(4);
@@ -171,7 +171,7 @@ public:
         inline PackedList(const PackedList<nBits>&);
 
         //- Construct by transferring the parameter contents
-        inline PackedList(const Xfer<PackedList<nBits> >&);
+        inline PackedList(const Xfer< PackedList<nBits> >&);
 
         //- Construct from a list of labels
         PackedList(const UList<label>&);
@@ -240,7 +240,6 @@ public:
 
         //- Reserve allocation space for at least this size.
         //  Never shrinks the allocated size.
-        //  Optionally provide an initialization value for new elements.
         inline void reserve(const label);
 
         //- Clear the list, i.e. set addressable size to zero.
@@ -258,7 +257,7 @@ public:
         inline void transfer(PackedList<nBits>&);
 
         //- Transfer contents to the Xfer container
-        inline Xfer<PackedList<nBits> > xfer();
+        inline Xfer< PackedList<nBits> > xfer();
 
 
     // Member operators
@@ -413,7 +412,7 @@ public:
         //- iterator set to the beginning of the PackedList
         inline iterator begin();
 
-        //- iterator set to beyond the end of the HashTable
+        //- iterator set to beyond the end of the PackedList
         inline iterator end();
 
 

src/OpenFOAM/matrices/solution/solution.C

@@ -134,7 +134,7 @@ Foam::label Foam::solution::upgradeSolverDict
 
 
             // write out information to help people adjust to the new syntax
-            if (verbose)
+            if (verbose && Pstream::master())
             {
                 Info<< "// using new solver syntax:\n"
                     << iter().keyword() << subdict << endl;

src/OpenFOAM/meshes/boundBox/boundBox.H

@@ -174,7 +174,7 @@ public:
 
         // Query
 
-            //- Completely contains other boundingBox? (inside or on edge)
+            //- Overlaps/touches boundingBox?
             bool overlaps(const boundBox& bb) const
             {
                 return

src/OpenFOAM/meshes/meshShapes/face/face.H

@@ -222,16 +222,16 @@ public:
         ) const;
 
         //- Fast intersection with a ray.
-        //  For a hit, the pointHit.distance() is the line parameter t :
-        //  intersection=p+t*q. Only defined for FULL_RAY or
-        //  HALF_RAY.
+        //  Does face-center decomposition and returns triangle intersection
+        //  point closest to p. See triangle::intersection for details.
         pointHit intersection
         (
             const point& p,
             const vector& q,
             const point& ctr,
             const pointField& meshPoints,
-            const intersection::algorithm alg
+            const intersection::algorithm alg,
+            const scalar tol = 0.0
         ) const;
 
         //- Return nearest point to face

src/OpenFOAM/meshes/meshShapes/face/faceIntersection.C

@@ -136,7 +136,8 @@ Foam::pointHit Foam::face::intersection
     const vector& q,
     const point& ctr,
     const pointField& meshPoints,
-    const intersection::algorithm alg
+    const intersection::algorithm alg,
+    const scalar tol
 ) const
 {
     scalar nearestHitDist = VGREAT;
@@ -154,7 +155,7 @@ Foam::pointHit Foam::face::intersection
             meshPoints[f[pI]],
             meshPoints[f[fcIndex(pI)]],
             ctr
-        ).intersection(p, q, alg);
+        ).intersection(p, q, alg, tol);
 
         if (curHit.hit())
         {

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

@@ -205,7 +205,7 @@ Foam::mapDistribute::mapDistribute
     const labelList& recvProcs
 )
 :
-    constructSize_(sendProcs.size()),
+    constructSize_(0),
     schedulePtr_()
 {
     if (sendProcs.size() != recvProcs.size())
@@ -266,6 +266,8 @@ Foam::mapDistribute::mapDistribute
         {
             // I am the receiver.
             constructMap_[sendProc][nRecv[sendProc]++] = sampleI;
+            // Largest entry inside constructMap
+            constructSize_ = sampleI+1;
         }
     }
 }

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

@@ -131,18 +131,36 @@ public:
                 return constructSize_;
             }
 
+            //- Constructed data size
+            label& constructSize()
+            {
+                return constructSize_;
+            }
+
             //- From subsetted data back to original data
             const labelListList& subMap() const
             {
                 return subMap_;
             }
 
+            //- From subsetted data back to original data
+            labelListList& subMap()
+            {
+                return subMap_;
+            }
+
             //- From subsetted data to new reconstructed data
             const labelListList& constructMap() const
             {
                 return constructMap_;
             }
 
+            //- From subsetted data to new reconstructed data
+            labelListList& constructMap()
+            {
+                return constructMap_;
+            }
+
             //- Calculate a schedule. See above.
             static List<labelPair> schedule
             (

src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C

@@ -52,13 +52,8 @@ void Foam::mapDistribute::distribute
 
             if (domain != Pstream::myProcNo() && map.size())
             {
-                List<T> subField(map.size());
-                forAll(map, i)
-                {
-                    subField[i] = field[map[i]];
-                }
                 OPstream toNbr(Pstream::blocking, domain);
-                toNbr << subField;
+                toNbr << UIndirectList<T>(field, map);
             }
         }
 
@@ -126,13 +121,7 @@ void Foam::mapDistribute::distribute
         List<T> newField(constructSize);
 
         // Subset myself
-        const labelList& mySubMap = subMap[Pstream::myProcNo()];
-
-        List<T> subField(mySubMap.size());
-        forAll(mySubMap, i)
-        {
-            subField[i] = field[mySubMap[i]];
-        }
+        UIndirectList<T> subField(field, subMap[Pstream::myProcNo()]);
 
         // Receive sub field from myself (subField)
         const labelList& map = constructMap[Pstream::myProcNo()];
@@ -152,16 +141,8 @@ void Foam::mapDistribute::distribute
             if (Pstream::myProcNo() == sendProc)
             {
                 // I am sender. Send to recvProc.
-                const labelList& map = subMap[recvProc];
-
-                List<T> subField(map.size());
-                forAll(map, i)
-                {
-                    subField[i] = field[map[i]];
-                }
-
                 OPstream toNbr(Pstream::scheduled, recvProc);
-                toNbr << subField;
+                toNbr << UIndirectList<T>(field, subMap[recvProc]);
             }
             else
             {
@@ -374,13 +355,8 @@ void Foam::mapDistribute::distribute
 
             if (domain != Pstream::myProcNo() && map.size())
             {
-                List<T> subField(map.size());
-                forAll(map, i)
-                {
-                    subField[i] = field[map[i]];
-                }
                 OPstream toNbr(Pstream::blocking, domain);
-                toNbr << subField;
+                toNbr << UIndirectList<T>(field, map);
             }
         }
 
@@ -449,13 +425,7 @@ void Foam::mapDistribute::distribute
         List<T> newField(constructSize, nullValue);
 
         // Subset myself
-        const labelList& mySubMap = subMap[Pstream::myProcNo()];
-
-        List<T> subField(mySubMap.size());
-        forAll(mySubMap, i)
-        {
-            subField[i] = field[mySubMap[i]];
-        }
+        UIndirectList<T> subField(field, subMap[Pstream::myProcNo()]);
 
         // Receive sub field from myself (subField)
         const labelList& map = constructMap[Pstream::myProcNo()];
@@ -475,16 +445,8 @@ void Foam::mapDistribute::distribute
             if (Pstream::myProcNo() == sendProc)
             {
                 // I am sender. Send to recvProc.
-                const labelList& map = subMap[recvProc];
-
-                List<T> subField(map.size());
-                forAll(map, i)
-                {
-                    subField[i] = field[map[i]];
-                }
-
                 OPstream toNbr(Pstream::scheduled, recvProc);
-                toNbr << subField;
+                toNbr << UIndirectList<T>(field, subMap[recvProc]);
             }
             else
             {

src/OpenFOAM/meshes/polyMesh/polyMeshInitMesh.C

@@ -82,10 +82,10 @@ void Foam::polyMesh::initMesh()
     );
 
     string meshInfo =
-        "nPoints: " + Foam::name(nPoints())
-      + " nCells: " + Foam::name(this->nCells())
-      + " nFaces: " + Foam::name(nFaces())
-      + " nInternalFaces: " + Foam::name(nInternalFaces());
+        "nPoints:" + Foam::name(nPoints())
+      + "  nCells:" + Foam::name(this->nCells())
+      + "  nFaces:" + Foam::name(nFaces())
+      + "  nInternalFaces:" + Foam::name(nInternalFaces());
 
     owner_.note() = meshInfo;
     neighbour_.note() = meshInfo;

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

@@ -58,7 +58,7 @@ Foam::PatchTools::checkOrientation
             {
                 Info<< "Face[" << faceI << "] " << p[faceI]
                     << " has fewer than 3 edges. Edges: " << edgeLabels
-                    << nl;
+                    << endl;
             }
             valid = false;
         }
@@ -70,10 +70,12 @@ Foam::PatchTools::checkOrientation
                 {
                     if (report)
                     {
-                        Info<< "edge number " << edgeLabels[i] << " on face " << faceI
+                        Info<< "edge number " << edgeLabels[i]
+                            << " on face " << faceI
                             << " out-of-range\n"
                             << "This usually means the input surface has "
-                            << "edges with more than 2 faces connected." << nl;
+                            << "edges with more than 2 faces connected."
+                            << endl;
                     }
                     valid = false;
                 }
@@ -91,9 +93,9 @@ Foam::PatchTools::checkOrientation
         //- Compute normal from 3 points, use the first as the origin
         // minor warpage should not be a problem
         const Face& f = p[faceI];
-        const point p0(p.points()[f[0]]);
-        const point p1(p.points()[f[1]]);
-        const point p2(p.points()[f[f.size()-1]]);
+        const point& p0 = p.points()[f[0]];
+        const point& p1 = p.points()[f[1]];
+        const point& p2 = p.points()[f[f.size()-1]];
 
         const vector pointNormal((p1 - p0) ^ (p2 - p0));
         if ((pointNormal & p.faceNormals()[faceI]) < 0)
@@ -103,12 +105,12 @@ Foam::PatchTools::checkOrientation
             if (report)
             {
                 Info
-                    << "Normal calculated from points inconsistent with faceNormal"
-                    << nl
+                    << "Normal calculated from points inconsistent"
+                    << " with faceNormal" << nl
                     << "face: " << f << nl
                     << "points: " << p0 << ' ' << p1 << ' ' << p2 << nl
                     << "pointNormal:" << pointNormal << nl
-                    << "faceNormal:" << p.faceNormals()[faceI] << nl;
+                    << "faceNormal:" << p.faceNormals()[faceI] << endl;
             }
         }
     }

src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H

@@ -153,7 +153,7 @@ public:
             inline scalar sweptVol(const triangle& t) const;
 
             //- Return point intersection with a ray.
-            //  For a hit, the distance is signed.  Positive number
+            //  For a hit, the distance is signed. Positive number
             //  represents the point in front of triangle.
             //  In case of miss pointHit is set to nearest point
             //  on triangle and its distance to the distance between
@@ -169,12 +169,14 @@ public:
             //- Fast intersection with a ray.
             //  For a hit, the pointHit.distance() is the line parameter t :
             //  intersection=p+t*q. Only defined for VISIBLE, FULL_RAY or
-            //  HALF_RAY.
+            //  HALF_RAY. tol increases the virtual size of the triangle
+            // by a relative factor.
             inline pointHit intersection
             (
                 const point& p,
                 const vector& q,
-                const intersection::algorithm alg
+                const intersection::algorithm alg,
+                const scalar tol = 0.0
             ) const;
 
             //- Return nearest point to p on triangle

src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H

@@ -439,7 +439,8 @@ inline pointHit triangle<Point, PointRef>::intersection
 (
     const point& orig,
     const vector& dir,
-    const intersection::algorithm alg
+    const intersection::algorithm alg,
+    const scalar tol
 ) const
 {
     const vector edge1 = b_ - a_;
@@ -457,99 +458,61 @@ inline pointHit triangle<Point, PointRef>::intersection
     if (alg == intersection::VISIBLE)
     {
         // Culling branch
-        if (det < SMALL)
+        if (det < ROOTVSMALL)
         {
-            // return miss
+            // Ray on wrong side of triangle. Return miss
             return intersection;
         }
-        /* calculate distance from a_ to ray origin */
-        const vector tVec = orig-a_;
-
-        /* calculate U parameter and test bounds */
-        scalar u = tVec & pVec;
-
-        if (u < 0.0 || u > det)
-        {
-            // return miss
-            return intersection;
-        }
-
-        /* prepare to test V parameter */
-        const vector qVec = tVec ^ edge1;
-
-        /* calculate V parameter and test bounds */
-        scalar v = dir & qVec;
-
-        if (v < 0.0 || u + v > det)
-        {
-            // return miss
-            return intersection;
-        }
-
-        /* calculate t, scale parameters, ray intersects triangle */
-        scalar t = edge2 & qVec;
-        scalar inv_det = 1.0 / det;
-        t *= inv_det;
-        u *= inv_det;
-        v *= inv_det;
-
-        intersection.setHit();
-        intersection.setPoint(a_ + u*edge1 + v*edge2);
-        intersection.setDistance(t);
     }
     else if (alg == intersection::HALF_RAY || alg == intersection::FULL_RAY)
     {
         // Non-culling branch
-        if (det > -SMALL && det < SMALL)
+        if (det > -ROOTVSMALL && det < ROOTVSMALL)
         {
-            // return miss
+            // Ray parallel to triangle. Return miss
             return intersection;
         }
-        const scalar inv_det = 1.0 / det;
-
-        /* calculate distance from a_ to ray origin */
-        const vector tVec = orig - a_;
-        /* calculate U parameter and test bounds */
-        const scalar u = (tVec & pVec)*inv_det;
-
-        if (u < 0.0 || u > 1.0)
-        {
-            // return miss
-            return intersection;
-        }
-        /* prepare to test V parameter */
-        const vector qVec = tVec ^ edge1;
-        /* calculate V parameter and test bounds */
-        const scalar v = (dir & qVec) * inv_det;
-
-        if (v < 0.0 || u + v > 1.0)
-        {
-            // return miss
-            return intersection;
-        }
-        /* calculate t, ray intersects triangle */
-        const scalar t = (edge2 & qVec) * inv_det;
-
-        if (alg == intersection::HALF_RAY && t < 0)
-        {
-            // return miss
-            return intersection;
-        }
-
-        intersection.setHit();
-        intersection.setPoint(a_ + u*edge1 + v*edge2);
-        intersection.setDistance(t);
     }
-    else
+
+    const scalar inv_det = 1.0 / det;
+
+    /* calculate distance from a_ to ray origin */
+    const vector tVec = orig-a_;
+
+    /* calculate U parameter and test bounds */
+    const scalar u = (tVec & pVec)*inv_det;
+
+    if (u < -tol || u > 1.0+tol)
     {
-        FatalErrorIn
-        (
-            "triangle<Point, PointRef>::intersection(const point&"
-            ", const vector&, const intersection::algorithm)"
-        )   << "intersection only defined for VISIBLE, FULL_RAY or HALF_RAY"
-            << abort(FatalError);
+        // return miss
+        return intersection;
     }
 
+    /* prepare to test V parameter */
+    const vector qVec = tVec ^ edge1;
+
+    /* calculate V parameter and test bounds */
+    const scalar v = (dir & qVec) * inv_det;
+
+    if (v < -tol || u + v > 1.0+tol)
+    {
+        // return miss
+        return intersection;
+    }
+
+    /* calculate t, scale parameters, ray intersects triangle */
+    const scalar t = (edge2 & qVec) * inv_det;
+
+    if (alg == intersection::HALF_RAY && t < -tol)
+    {
+        // Wrong side of orig. Return miss
+        return intersection;
+    }
+
+    intersection.setHit();
+    intersection.setPoint(a_ + u*edge1 + v*edge2);
+    intersection.setDistance(t);
+
     return intersection;
 }
 
@@ -622,7 +585,7 @@ inline bool triangle<Point, PointRef>::classify
 
     bool hit = false;
 
-    if (Foam::mag(u1) < SMALL)
+    if (Foam::mag(u1) < ROOTVSMALL)
     {
         beta = u0/u2;
 

src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.H

@@ -38,7 +38,6 @@ SourceFiles
 
 #include "autoPtr.H"
 #include "dictionary.H"
-#include "boolList.H"
 #include "wallPoint.H"
 #include "searchableSurfaces.H"
 #include "refinementSurfaces.H"

src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C

@@ -45,6 +45,7 @@ Description
 #include "OFstream.H"
 #include "layerParameters.H"
 #include "combineFaces.H"
+#include "IOmanip.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -1414,8 +1415,11 @@ void Foam::autoLayerDriver::calculateLayerThickness
     const indirectPrimitivePatch& pp,
     const labelList& patchIDs,
     const scalarField& patchExpansionRatio,
-    const scalarField& patchFinalLayerRatio,
-    const scalarField& patchRelMinThickness,
+
+    const bool relativeSizes,
+    const scalarField& patchFinalLayerThickness,
+    const scalarField& patchMinThickness,
+
     const labelList& cellLevel,
     const labelList& patchNLayers,
     const scalar edge0Len,
@@ -1428,22 +1432,100 @@ void Foam::autoLayerDriver::calculateLayerThickness
     const fvMesh& mesh = meshRefiner_.mesh();
     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 
-    if (min(patchRelMinThickness) < 0 || max(patchRelMinThickness) > 2)
+
+    // Rework patch-wise layer parameters into minimum per point
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Reuse input fields
+    expansionRatio.setSize(pp.nPoints());
+    expansionRatio = GREAT;
+    thickness.setSize(pp.nPoints());
+    thickness = GREAT;
+    minThickness.setSize(pp.nPoints());
+    minThickness = GREAT;
+
+    forAll(patchIDs, i)
     {
-        FatalErrorIn("calculateLayerThickness(..)")
-            << "Thickness should be factor of local undistorted cell size."
-            << " Valid values are [0..2]." << nl
-            << " minThickness:" << patchRelMinThickness
-            << exit(FatalError);
+        label patchI = patchIDs[i];
+
+        const labelList& meshPoints = patches[patchI].meshPoints();
+
+        forAll(meshPoints, patchPointI)
+        {
+            label ppPointI = pp.meshPointMap()[meshPoints[patchPointI]];
+
+            expansionRatio[ppPointI] = min
+            (
+                expansionRatio[ppPointI],
+                patchExpansionRatio[patchI]
+            );
+            thickness[ppPointI] = min
+            (
+                thickness[ppPointI],
+                patchFinalLayerThickness[patchI]
+            );
+            minThickness[ppPointI] = min
+            (
+                minThickness[ppPointI],
+                patchMinThickness[patchI]
+            );
+        }
     }
 
+    syncTools::syncPointList
+    (
+        mesh,
+        pp.meshPoints(),
+        expansionRatio,
+        minEqOp<scalar>(),
+        GREAT,              // null value
+        false               // no separation
+    );
+    syncTools::syncPointList
+    (
+        mesh,
+        pp.meshPoints(),
+        thickness,
+        minEqOp<scalar>(),
+        GREAT,              // null value
+        false               // no separation
+    );
+    syncTools::syncPointList
+    (
+        mesh,
+        pp.meshPoints(),
+        minThickness,
+        minEqOp<scalar>(),
+        GREAT,              // null value
+        false               // no separation
+    );
 
-    // Per point the max cell level of connected cells
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    labelList maxPointLevel(pp.nPoints(), labelMin);
+    // Now the thicknesses are set according to the minimum of connected
+    // patches.
 
+
+    // Rework relative thickness into absolute
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // by multiplying with the internal cell size.
+
+    if (relativeSizes)
     {
+        if (min(patchMinThickness) < 0 || max(patchMinThickness) > 2)
+        {
+            FatalErrorIn("calculateLayerThickness(..)")
+                << "Thickness should be factor of local undistorted cell size."
+                << " Valid values are [0..2]." << nl
+                << " minThickness:" << patchMinThickness
+                << exit(FatalError);
+        }
+
+
+        // Determine per point the max cell level of connected cells
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        labelList maxPointLevel(pp.nPoints(), labelMin);
+
         forAll(pp, i)
         {
             label ownLevel = cellLevel[mesh.faceOwner()[pp.addressing()[i]]];
@@ -1465,113 +1547,44 @@ void Foam::autoLayerDriver::calculateLayerThickness
             labelMin,           // null value
             false               // no separation
         );
-    }
 
 
-    // Rework patch-wise layer parameters into minimum per point
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    expansionRatio.setSize(pp.nPoints());
-    expansionRatio = GREAT;
-    scalarField finalLayerRatio(pp.nPoints(), GREAT);
-    scalarField relMinThickness(pp.nPoints(), GREAT);
-
-    {
-        forAll(patchIDs, i)
+        forAll(maxPointLevel, pointI)
         {
-            label patchI = patchIDs[i];
-
-            const labelList& meshPoints = patches[patchI].meshPoints();
-
-            forAll(meshPoints, patchPointI)
-            {
-                label ppPointI = pp.meshPointMap()[meshPoints[patchPointI]];
-
-                expansionRatio[ppPointI] = min
-                (
-                    expansionRatio[ppPointI],
-                    patchExpansionRatio[patchI]
-                );
-                finalLayerRatio[ppPointI] = min
-                (
-                    finalLayerRatio[ppPointI],
-                    patchFinalLayerRatio[patchI]
-                );
-                relMinThickness[ppPointI] = min
-                (
-                    relMinThickness[ppPointI],
-                    patchRelMinThickness[patchI]
-                );
-            }
+            // Find undistorted edge size for this level.
+            scalar edgeLen = edge0Len/(1<<maxPointLevel[pointI]);
+            thickness[pointI] *= edgeLen;
+            minThickness[pointI] *= edgeLen;
         }
-
-        syncTools::syncPointList
-        (
-            mesh,
-            pp.meshPoints(),
-            expansionRatio,
-            minEqOp<scalar>(),
-            GREAT,              // null value
-            false               // no separation
-        );
-        syncTools::syncPointList
-        (
-            mesh,
-            pp.meshPoints(),
-            finalLayerRatio,
-            minEqOp<scalar>(),
-            GREAT,              // null value
-            false               // no separation
-        );
-        syncTools::syncPointList
-        (
-            mesh,
-            pp.meshPoints(),
-            relMinThickness,
-            minEqOp<scalar>(),
-            GREAT,              // null value
-            false               // no separation
-        );
     }
 
 
 
-    // Per mesh point the expansion parameters
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Rework thickness (of final layer) into overall thickness of all layers
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    thickness.setSize(pp.nPoints());
-    minThickness.setSize(pp.nPoints());
-
-    forAll(maxPointLevel, pointI)
+    forAll(thickness, pointI)
     {
-        // Find undistorted edge size for this level.
-        scalar edgeLen = edge0Len/(1<<maxPointLevel[pointI]);
-
         // Calculate layer thickness based on expansion ratio
         // and final layer height
         if (expansionRatio[pointI] == 1)
         {
-            thickness[pointI] =
-                finalLayerRatio[pointI]
-              * patchNLayers[pointI]
-              * edgeLen;
-            minThickness[pointI] = relMinThickness[pointI]*edgeLen;
+            thickness[pointI] *= patchNLayers[pointI];
         }
         else
         {
+
             scalar invExpansion = 1.0 / expansionRatio[pointI];
             label nLay = patchNLayers[pointI];
-            thickness[pointI] =
-                finalLayerRatio[pointI]
-              * edgeLen
-              * (1.0 - pow(invExpansion, nLay))
+            thickness[pointI] *=
+                (1.0 - pow(invExpansion, nLay))
               / (1.0 - invExpansion);
-            minThickness[pointI] = relMinThickness[pointI]*edgeLen;
         }
     }
 
-    Info<< "calculateLayerThickness : min:" << gMin(thickness)
-        << " max:" << gMax(thickness) << endl;
+
+    //Info<< "calculateLayerThickness : min:" << gMin(thickness)
+    //    << " max:" << gMax(thickness) << endl;
 }
 
 
@@ -2292,7 +2305,7 @@ bool Foam::autoLayerDriver::cellsUseFace
 Foam::label Foam::autoLayerDriver::checkAndUnmark
 (
     const addPatchCellLayer& addLayer,
-    const dictionary& motionDict,
+    const dictionary& meshQualityDict,
     const indirectPrimitivePatch& pp,
     const fvMesh& newMesh,
 
@@ -2304,7 +2317,7 @@ Foam::label Foam::autoLayerDriver::checkAndUnmark
     // Check the resulting mesh for errors
     Info<< nl << "Checking mesh with layer ..." << endl;
     faceSet wrongFaces(newMesh, "wrongFaces", newMesh.nFaces()/1000);
-    motionSmoother::checkMesh(false, newMesh, motionDict, wrongFaces);
+    motionSmoother::checkMesh(false, newMesh, meshQualityDict, wrongFaces);
     Info<< "Detected " << returnReduce(wrongFaces.size(), sumOp<label>())
         << " illegal faces"
         << " (concave, zero area or negative cell pyramid volume)"
@@ -2474,8 +2487,8 @@ void Foam::autoLayerDriver::mergePatchFacesUndo
         << "      (cos:" << minCos << ')' << nl
         << "    - as long as the resulting face doesn't become concave"
         << " by more than "
-        << layerParams.concaveAngle()
-        << " degrees (0=straight, 180=fully concave)" << nl
+        << layerParams.concaveAngle() << " degrees" << nl
+        << "      (0=straight, 180=fully concave)" << nl
         << endl;
 
     label nChanged = mergePatchFacesUndo(minCos, concaveCos, motionDict);
@@ -2621,8 +2634,11 @@ void Foam::autoLayerDriver::addLayers
         pp,
         meshMover.adaptPatchIDs(),
         layerParams.expansionRatio(),
-        layerParams.finalLayerRatio(),
-        layerParams.minThickness(),
+
+        layerParams.relativeSizes(),        // thickness relative to cellsize?
+        layerParams.finalLayerThickness(),  // wanted thicknes
+        layerParams.minThickness(),         // minimum thickness
+
         cellLevel,
         patchNLayers,
         edge0Len,
@@ -2632,6 +2648,79 @@ void Foam::autoLayerDriver::addLayers
         expansionRatio
     );
 
+
+    // Print a bit
+    {
+        const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+        Info<< nl
+            << "patch               faces    layers avg thickness[m]" << nl
+            << "                                    near-wall overall" << nl
+            << "-----               -----    ------ --------- -------" << endl;
+
+        forAll(meshMover.adaptPatchIDs(), i)
+        {
+            label patchI = meshMover.adaptPatchIDs()[i];
+
+            const labelList& meshPoints = patches[patchI].meshPoints();
+
+            //scalar maxThickness = -VGREAT;
+            //scalar minThickness = VGREAT;
+            scalar sumThickness = 0;
+            scalar sumNearWallThickness = 0;
+
+            forAll(meshPoints, patchPointI)
+            {
+                label ppPointI = pp.meshPointMap()[meshPoints[patchPointI]];
+
+                //maxThickness = max(maxThickness, thickness[ppPointI]);
+                //minThickness = min(minThickness, thickness[ppPointI]);
+                sumThickness += thickness[ppPointI];
+
+                label nLay = patchNLayers[ppPointI];
+                if (nLay > 0)
+                {
+                    if (expansionRatio[ppPointI] == 1)
+                    {
+                        sumNearWallThickness += thickness[ppPointI]/nLay;
+                    }
+                    else
+                    {
+                        scalar s =
+                            (1.0-pow(expansionRatio[ppPointI], nLay))
+                          / (1.0-expansionRatio[ppPointI]);
+                        sumNearWallThickness += thickness[ppPointI]/s;
+                    }
+                }
+            }
+
+            label totNPoints = returnReduce(meshPoints.size(), sumOp<label>());
+
+            //reduce(maxThickness, maxOp<scalar>());
+            //reduce(minThickness, minOp<scalar>());
+            scalar avgThickness =
+                returnReduce(sumThickness, sumOp<scalar>())
+              / totNPoints;
+            scalar avgNearWallThickness =
+                returnReduce(sumNearWallThickness, sumOp<scalar>())
+              / totNPoints;
+
+            Info<< setf(ios_base::left) << setw(19) << patches[patchI].name();
+            //Sout.unsetf(ios_base::left);
+            Info<< setprecision(3)
+                << " " << setw(8)
+                << returnReduce(patches[patchI].size(), sumOp<scalar>())
+                << " " << setw(6) << layerParams.numLayers()[patchI]
+                << " " << setw(8) << avgNearWallThickness
+                << "  " << setw(8) << avgThickness
+                //<< " " << setw(8) << minThickness
+                //<< " " << setw(8) << maxThickness
+                << endl;
+        }
+        Info<< endl;
+    }
+
+
     // Calculate wall to medial axis distance for smoothing displacement
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -2709,8 +2798,12 @@ void Foam::autoLayerDriver::addLayers
     boolList flaggedCells;
     boolList flaggedFaces;
 
-    while (true)
+    for (label iteration = 0; iteration < layerParams.nLayerIter(); iteration++)
     {
+        Info<< nl
+            << "Layer addition iteration " << iteration << nl
+            << "--------------------------" << endl;
+
         // Make sure displacement is equal on both sides of coupled patches.
         syncPatchDisplacement
         (
@@ -2831,7 +2924,8 @@ void Foam::autoLayerDriver::addLayers
             nPatchFaceLayers
         );
 
-        // Calculate displacement for first layer for addPatchLayer
+        // Calculate displacement for first layer for addPatchLayer.
+        // (first layer = layer of cells next to the original mesh)
         vectorField firstDisp(patchNLayers.size(), vector::zero);
 
         forAll(patchNLayers, i)
@@ -2847,9 +2941,9 @@ void Foam::autoLayerDriver::addLayers
                     label nLay = nPatchPointLayers[i];
                     scalar h =
                         pow(expansionRatio[i], nLay - 1)
-                      * (mag(patchDisp[i])*(1.0 - expansionRatio[i]))
+                      * (1.0 - expansionRatio[i])
                       / (1.0 - pow(expansionRatio[i], nLay));
-                    firstDisp[i] = h/mag(patchDisp[i])*patchDisp[i];
+                    firstDisp[i] = h*patchDisp[i];
                 }
             }
         }
@@ -2864,7 +2958,7 @@ void Foam::autoLayerDriver::addLayers
             pp,
             nPatchFaceLayers,   // layers per face
             nPatchPointLayers,  // layers per point
-            firstDisp,          // thickness of first layer
+            firstDisp,          // thickness of layer nearest internal mesh
             meshMod
         );
 
@@ -2951,10 +3045,23 @@ void Foam::autoLayerDriver::addLayers
         }
 
         // Unset the extrusion at the pp.
+        const dictionary& meshQualityDict =
+        (
+            iteration < layerParams.nRelaxedIter()
+          ? motionDict
+          : motionDict.subDict("relaxed")
+        );
+
+        if (iteration >= layerParams.nRelaxedIter())
+        {
+            Info<< "Switched to relaxed meshQuality constraints." << endl;
+        }
+
+
         label nTotChanged = checkAndUnmark
         (
             addLayer,
-            motionDict,
+            meshQualityDict,
             pp,
             newMesh,
 
@@ -2976,6 +3083,8 @@ void Foam::autoLayerDriver::addLayers
         // Reset mesh points and start again
         meshMover.movePoints(oldPoints);
         meshMover.correct();
+
+        Info<< endl;
     }
 
 

src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H

@@ -255,9 +255,12 @@ class autoLayerDriver
                 (
                     const indirectPrimitivePatch& pp,
                     const labelList& patchIDs,
+
                     const scalarField& patchExpansionRatio,
-                    const scalarField& patchFinalLayerRatio,
-                    const scalarField& patchRelMinThickness,
+                    const bool relativeSizes,
+                    const scalarField& patchFinalLayerThickness,
+                    const scalarField& patchMinThickness,
+
                     const labelList& cellLevel,
                     const labelList& patchNLayers,
                     const scalar edge0Len,

src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C

@@ -95,7 +95,7 @@ Foam::label Foam::autoRefineDriver::readFeatureEdges
         Info<< "Refinement level " << featureLevels[i]
             << " for all cells crossed by feature " << featFileName
             << " (" << featureMeshes[i].points().size() << " points, "
-            << featureMeshes[i].edges().size() << ")." << endl;
+            << featureMeshes[i].edges().size() << " edges)." << endl;
     }
 
     Info<< "Read feature lines in = "

src/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C

@@ -163,7 +163,8 @@ Foam::layerParameters::layerParameters
         numLayers_.size(),
         readScalar(dict.lookup("expansionRatio"))
     ),
-    finalLayerRatio_
+    relativeSizes_(false),
+    finalLayerThickness_
     (
         numLayers_.size(),
         readScalar(dict.lookup("finalLayerRatio"))
@@ -211,8 +212,24 @@ Foam::layerParameters::layerParameters
     (
         readLabel(dict.lookup("nBufferCellsNoExtrude"))
     ),
-    nSnap_(readLabel(dict.lookup("nSnap")))
-{}
+    nSnap_(readLabel(dict.lookup("nSnap"))),
+    nLayerIter_(readLabel(dict.lookup("nLayerIter"))),
+    nRelaxedIter_(labelMax)
+{
+    if (dict.found("nRelaxedIter"))
+    {
+        dict.lookup("nRelaxedIter") >> nRelaxedIter_;
+    }
+
+    if (nLayerIter_ < 0 || nRelaxedIter_ < 0)
+    {
+        FatalErrorIn("layerParameters::layerParameters(..)")
+            << "Layer iterations should be >= 0." << endl
+            << "nLayerIter:" << nLayerIter_
+            << " nRelaxedIter:" << nRelaxedIter_
+            << exit(FatalError);
+    }
+}
 
 
 // Construct from dictionary
@@ -228,10 +245,11 @@ Foam::layerParameters::layerParameters
         boundaryMesh.size(),
         readScalar(dict.lookup("expansionRatio"))
     ),
-    finalLayerRatio_
+    relativeSizes_(dict.lookup("relativeSizes")),
+    finalLayerThickness_
     (
         boundaryMesh.size(),
-        readScalar(dict.lookup("finalLayerRatio"))
+        readScalar(dict.lookup("finalLayerThickness"))
     ),
     minThickness_
     (
@@ -276,8 +294,24 @@ Foam::layerParameters::layerParameters
     (
         readLabel(dict.lookup("nBufferCellsNoExtrude"))
     ),
-    nSnap_(readLabel(dict.lookup("nRelaxIter")))
+    nSnap_(readLabel(dict.lookup("nRelaxIter"))),
+    nLayerIter_(readLabel(dict.lookup("nLayerIter"))),
+    nRelaxedIter_(labelMax)
 {
+    if (dict.found("nRelaxedIter"))
+    {
+        dict.lookup("nRelaxedIter") >> nRelaxedIter_;
+    }
+    if (nLayerIter_ < 0 || nRelaxedIter_ < 0)
+    {
+        FatalErrorIn("layerParameters::layerParameters(..)")
+            << "Layer iterations should be >= 0." << endl
+            << "nLayerIter:" << nLayerIter_
+            << " nRelaxedIter:" << nRelaxedIter_
+            << exit(FatalError);
+    }
+
+
     const dictionary& layersDict = dict.subDict("layers");
 
     forAll(boundaryMesh, patchI)
@@ -291,14 +325,21 @@ Foam::layerParameters::layerParameters
             numLayers_[patchI] =
                 readLabel(layerDict.lookup("nSurfaceLayers"));
 
-            //- Patch-wise layer parameters disabled for now. Just remove
-            //  settings in initialiser list and uncomment below.
-            //expansionRatio_[patchI] =
-            //    readScalar(layerDict.lookup("expansionRatio"));
-            //finalLayerRatio_[patchI] =
-            //    readScalar(layerDict.lookup("finalLayerRatio"));
-            //minThickness_[patchI] =
-            //    readScalar(layerDict.lookup("minThickness"));
+            layerDict.readIfPresent
+            (
+                "expansionRatio",
+                expansionRatio_[patchI]
+            );
+            layerDict.readIfPresent
+            (
+                "finalLayerThickness",
+                finalLayerThickness_[patchI]
+            );
+            layerDict.readIfPresent
+            (
+                "minThickness",
+                minThickness_[patchI]
+            );
         }
     }
 

src/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.H

@@ -39,6 +39,7 @@ SourceFiles
 #include "dictionary.H"
 #include "scalarField.H"
 #include "labelList.H"
+#include "Switch.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -70,14 +71,10 @@ class layerParameters
 
             scalarField expansionRatio_;
 
-            //- Wanted thickness of final added cell layer. If multiple layers
-            //  is the thickness of the layer furthest away from the wall.
-            //  Relative to undistorted size of cell outside layer.
-            scalarField finalLayerRatio_;
+            Switch relativeSizes_;
+
+            scalarField finalLayerThickness_;
 
-            //- Minimum thickness of cell layer. If for any reason layer
-            //  cannot be above minThickness do not add layer.
-            //  Relative to undistorted size of cell outside layer.
             scalarField minThickness_;
 
 
@@ -105,6 +102,10 @@ class layerParameters
 
         label nSnap_;
 
+        label nLayerIter_;
+
+        label nRelaxedIter_;
+
 
     // Private Member Functions
 
@@ -160,18 +161,27 @@ public:
                     return expansionRatio_;
                 }
 
-                //- Wanted thickness of final added cell layer. If multiple
-                //  layers
-                //  is the thickness of the layer furthest away from the wall.
-                //  Relative to undistorted size of cell outside layer.
-                const scalarField& finalLayerRatio() const
+                //- Are size parameters relative to inner cell size or
+                //  absolute distances.
+                bool relativeSizes() const
                 {
-                    return finalLayerRatio_;
+                    return relativeSizes_;
+                }
+
+                //- Wanted thickness of final added cell layer. If multiple
+                //  layers is the thickness of the layer furthest away
+                //  from the wall (i.e. nearest the original mesh)
+                //  If relativeSize() this number is relative to undistorted
+                //  size of the cell outside layer.
+                const scalarField& finalLayerThickness() const
+                {
+                    return finalLayerThickness_;
                 }
 
                 //- Minimum thickness of cell layer. If for any reason layer
                 //  cannot be above minThickness do not add layer.
-                //  Relative to undistorted size of cell outside layer.
+                //  If relativeSize() this number is relative to undistorted
+                //  size of the cell outside layer.
                 const scalarField& minThickness() const
                 {
                     return minThickness_;
@@ -196,26 +206,20 @@ public:
                 return nGrow_;
             }
 
-            // Number of smoothing iterations of surface normals 
+            //- Number of smoothing iterations of surface normals 
             label nSmoothSurfaceNormals() const
             {
                 return nSmoothSurfaceNormals_;
             }
 
-            // Number of smoothing iterations of interior mesh movement
-            // direction  
+            //- Number of smoothing iterations of interior mesh movement
+            //  direction  
             label nSmoothNormals() const
             {
                 return nSmoothNormals_;
             }
 
-            // Smooth layer thickness over surface patches
-            label nSmoothThickness() const
-            {
-                return nSmoothThickness_;
-            }
-
-            // Stop layer growth on highly warped cells 
+            //- Stop layer growth on highly warped cells 
             scalar maxFaceThicknessRatio() const
             {
                 return maxFaceThicknessRatio_;
@@ -226,20 +230,26 @@ public:
                 return layerTerminationCos_;
             }
 
-            // Reduce layer growth where ratio thickness to medial 
-            // distance is large 
+            //- Smooth layer thickness over surface patches
+            label nSmoothThickness() const
+            {
+                return nSmoothThickness_;
+            }
+
+            //- Reduce layer growth where ratio thickness to medial 
+            //  distance is large 
             scalar maxThicknessToMedialRatio() const
             {
                 return maxThicknessToMedialRatio_;
             }
 
-            // Angle used to pick up medial axis points
+            //- Angle used to pick up medial axis points
             scalar minMedianAxisAngleCos() const
             {
                 return minMedianAxisAngleCos_;
             }
 
-            // Create buffer region for new layer terminations
+            //- Create buffer region for new layer terminations
             label nBufferCellsNoExtrude() const
             {
                 return nBufferCellsNoExtrude_;
@@ -249,6 +259,24 @@ public:
             {
                 return nSnap_;
             }
+
+            // Overall
+
+                //- Number of overall layer addition iterations
+                label nLayerIter() const
+                {
+                    return nLayerIter_;
+                }
+
+                //- Number of iterations after which relaxed motion rules
+                //  are to be used.
+                label nRelaxedIter() const
+                {
+                    return nRelaxedIter_;
+                }
+
+
+
 };
 
 

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

@@ -42,9 +42,7 @@ SourceFiles
 
 #include "point.H"
 #include "label.H"
-#include "scalar.H"
 #include "tensor.H"
-#include "pTraits.H"
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

@@ -22,8 +22,6 @@ 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"

src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H

@@ -47,7 +47,6 @@ SourceFiles
 #include "hexRef8.H"
 #include "mapPolyMesh.H"
 #include "autoPtr.H"
-#include "pointIOField.H"
 #include "labelPair.H"
 #include "indirectPrimitivePatch.H"
 #include "pointFieldsFwd.H"

src/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C

@@ -328,7 +328,7 @@ Foam::label Foam::meshRefinement::markFeatureRefinement
     // Database to pass into trackedParticle::move
     trackedParticle::trackData td(cloud, maxFeatureLevel);
 
-    // Track all particles to their end position.
+    // Track all particles to their end position (= starting feature point)
     cloud.move(td);
 
     // Reset level

src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H

@@ -38,7 +38,6 @@ SourceFiles
 #ifndef refinementSurfaces_H
 #define refinementSurfaces_H
 
-#include "PackedBoolList.H"
 #include "triSurfaceGeoMesh.H"
 #include "triSurfaceFields.H"
 #include "vectorList.H"

src/autoMesh/autoHexMesh/trackedParticle/ExactParticle.C

@@ -102,23 +102,7 @@ Foam::scalar Foam::ExactParticle<ParticleType>::trackToFace
         }
     }
 
-
-    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
+    if (hitFacei == -1)
     {
         // Did not find any intersection. Fall back to original approximate
         // algorithm
@@ -129,8 +113,23 @@ Foam::scalar Foam::ExactParticle<ParticleType>::trackToFace
         );
     }
 
+    if (trackFraction >= (1.0-SMALL))
+    {
+        // Nearest intersection beyond endPosition so we hit endPosition.
+        trackFraction = 1.0;
+        this->position_ = endPosition;
+        this->facei_ = -1;
+        return 1.0;
+    }
+    else
+    {
+        this->position_ = intersection;
+        this->facei_ = hitFacei;
+    }
 
-    // Normal situation (trackFraction 0..1)
+
+    // Normal situation (trackFraction 0..1). Straight copy
+    // of Particle::trackToFace.
 
     bool internalFace = this->cloud().internalFace(this->facei_);
 

src/decompositionMethods/decompositionMethods/Make/files

@@ -2,7 +2,10 @@ decompositionMethod/decompositionMethod.C
 geomDecomp/geomDecomp.C
 simpleGeomDecomp/simpleGeomDecomp.C
 hierarchGeomDecomp/hierarchGeomDecomp.C
-metisDecomp/metisDecomp.C
 manualDecomp/manualDecomp.C
 
+scotchDecomp/scotchDecomp.C
+
+metisDecomp/metisDecomp.C
+
 LIB = $(FOAM_LIBBIN)/libdecompositionMethods

src/decompositionMethods/decompositionMethods/Make/options

@@ -1,6 +1,8 @@
 EXE_INC = \
+    -I$(WM_THIRD_PARTY_DIR)/scotch_5.1/src/libscotch/lnInclude \
     -I$(WM_THIRD_PARTY_DIR)/metis-5.0pre2/include
 
 LIB_LIBS = \
+    -lscotch \
     -lmetis \
     -lGKlib

src/decompositionMethods/decompositionMethods/scotchDecomp/scotchDecomp.C

@@ -0,0 +1,425 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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
+
+S       Strat=b{job=t,map=t,poli=S,sep=(m{asc=b{bnd=d{pass=40,dif=1,rem=1}f{move=80,pass=-1,bal=0.005},org=f{move=80,pass=-1,bal=0.005},width=3},low=h{pass=10}f{move=80,pass=-1,bal=0.0005},type=h,vert=80,rat=0.8}|m{asc=b{bnd=d{pass=40,dif=1,rem=1}f{move=80,pass=-1,bal=0.005},org=f{move=80,pass=-1,bal=0.005},width=3},low=h{pass=10}f{move=80,pass=-1,bal=0.0005},type=h,vert=80,rat=0.8})}
+\*---------------------------------------------------------------------------*/
+
+#include "scotchDecomp.H"
+#include "addToRunTimeSelectionTable.H"
+#include "floatScalar.H"
+#include "IFstream.H"
+#include "Time.H"
+#include "cyclicPolyPatch.H"
+
+extern "C"
+{
+#define OMPI_SKIP_MPICXX
+#include "module.h"
+#include "common.h"
+#include "scotch.h"
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(scotchDecomp, 0);
+
+    addToRunTimeSelectionTable
+    (
+        decompositionMethod,
+        scotchDecomp,
+        dictionaryMesh
+    );
+}
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::scotchDecomp::check(const int retVal, const char* str)
+{
+    if (retVal)
+    {
+        FatalErrorIn("scotchDecomp::decompose(..)")
+            << "Called to scotch routine " << str << " failed."
+            << exit(FatalError);
+    }
+}
+
+
+// Call Metis with options from dictionary.
+Foam::label Foam::scotchDecomp::decompose
+(
+    const List<int>& adjncy,
+    const List<int>& xadj,
+
+    List<int>& finalDecomp
+)
+{
+    // Strategy
+    // ~~~~~~~~
+    // Default.
+    SCOTCH_Strat stradat;
+    check(SCOTCH_stratInit(&stradat), "SCOTCH_stratInit");
+    //SCOTCH_stratGraphMap(&stradat, &argv[i][2]);
+    //fprintf(stdout, "S\tStrat=");
+    //SCOTCH_stratSave(&stradat, stdout);
+    //fprintf(stdout, "\n");
+
+
+    // Graph
+    // ~~~~~
+
+    SCOTCH_Graph grafdat;
+    check(SCOTCH_graphInit(&grafdat), "SCOTCH_graphInit");
+    check
+    (
+        SCOTCH_graphBuild
+        (
+            &grafdat,
+            0,                      // baseval, c-style numbering
+            xadj.size()-1,          // vertnbr, nCells
+            xadj.begin(),           // verttab, start index per cell into adjncy
+            &xadj[1],               // vendtab, end index  ,,
+            NULL,                   // velotab, vertex weights
+            NULL,                   // vlbltab
+            adjncy.size(),          // edgenbr, number of arcs
+            adjncy.begin(),         // edgetab
+            NULL                    // edlotab, edge weights
+        ),
+        "SCOTCH_graphBuild"
+    );
+    check(SCOTCH_graphCheck(&grafdat), "SCOTCH_graphCheck");
+
+
+
+    //// Architecture
+    //// ~~~~~~~~~~~~
+    //// (fully connected network topology since using switch)
+    //
+    //SCOTCH_Arch archdat;
+    //check(SCOTCH_archInit(&archdat), "SCOTCH_archInit");
+    //check
+    //(
+    //    // SCOTCH_archCmpltw for weighted.
+    //    SCOTCH_archCmplt(&archdat, nProcessors_),
+    //    "SCOTCH_archCmplt"
+    //);
+
+
+
+
+    //SCOTCH_Mapping mapdat;
+    //SCOTCH_graphMapInit(&grafdat, &mapdat, &archdat, NULL);
+    //SCOTCH_graphMapCompute(&grafdat, &mapdat, &stradat); /* Perform mapping */
+
+
+    finalDecomp.setSize(xadj.size()-1);
+    check
+    (
+        SCOTCH_graphPart
+        (
+            &grafdat,
+            nProcessors_,       // partnbr
+            &stradat,           // const SCOTCH_Strat * 
+            finalDecomp.begin() // parttab
+        ),
+        "SCOTCH_graphPart"
+    );
+
+    // Release storage for graph
+    SCOTCH_graphExit(&grafdat);
+    SCOTCH_stratExit(&stradat);
+    //// Release storage for network topology
+    //SCOTCH_archExit(&archdat);
+
+    return 0;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::scotchDecomp::scotchDecomp
+(
+    const dictionary& decompositionDict,
+    const polyMesh& mesh
+)
+:
+    decompositionMethod(decompositionDict),
+    mesh_(mesh)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::labelList Foam::scotchDecomp::decompose(const pointField& points)
+{
+    if (points.size() != mesh_.nCells())
+    {
+        FatalErrorIn("scotchDecomp::decompose(const pointField&)")
+            << "Can use this decomposition method only for the whole mesh"
+            << endl
+            << "and supply one coordinate (cellCentre) for every cell." << endl
+            << "The number of coordinates " << points.size() << endl
+            << "The number of cells in the mesh " << mesh_.nCells()
+            << exit(FatalError);
+    }
+
+    // Make Metis CSR (Compressed Storage Format) storage
+    //   adjncy      : contains neighbours (= edges in graph)
+    //   xadj(celli) : start of information in adjncy for celli
+
+    List<int> xadj(mesh_.nCells()+1);
+
+    // Initialise the number of internal faces of the cells to twice the
+    // number of internal faces
+    label nInternalFaces = 2*mesh_.nInternalFaces();
+
+    // Check the boundary for coupled patches and add to the number of
+    // internal faces
+    const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
+
+    forAll(pbm, patchi)
+    {
+        if (isA<cyclicPolyPatch>(pbm[patchi]))
+        {
+            nInternalFaces += pbm[patchi].size();
+        }
+    }
+
+    // Create the adjncy array the size of the total number of internal and
+    // coupled faces
+    List<int> adjncy(nInternalFaces);
+
+    // Fill in xadj
+    // ~~~~~~~~~~~~
+    label freeAdj = 0;
+
+    for (label cellI = 0; cellI < mesh_.nCells(); cellI++)
+    {
+        xadj[cellI] = freeAdj;
+
+        const labelList& cFaces = mesh_.cells()[cellI];
+
+        forAll(cFaces, i)
+        {
+            label faceI = cFaces[i];
+
+            if
+            (
+                mesh_.isInternalFace(faceI)
+             || isA<cyclicPolyPatch>(pbm[pbm.whichPatch(faceI)])
+            )
+            {
+                freeAdj++;
+            }
+        }
+    }
+    xadj[mesh_.nCells()] = freeAdj;
+
+
+    // Fill in adjncy
+    // ~~~~~~~~~~~~~~
+
+    labelList nFacesPerCell(mesh_.nCells(), 0);
+
+    // Internal faces
+    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
+    {
+        label own = mesh_.faceOwner()[faceI];
+        label nei = mesh_.faceNeighbour()[faceI];
+
+        adjncy[xadj[own] + nFacesPerCell[own]++] = nei;
+        adjncy[xadj[nei] + nFacesPerCell[nei]++] = own;
+    }
+
+    // Coupled faces. Only cyclics done.
+    forAll(pbm, patchi)
+    {
+        if (isA<cyclicPolyPatch>(pbm[patchi]))
+        {
+            const unallocLabelList& faceCells = pbm[patchi].faceCells();
+
+            label sizeby2 = faceCells.size()/2;
+
+            for (label facei=0; facei<sizeby2; facei++)
+            {
+                label own = faceCells[facei];
+                label nei = faceCells[facei + sizeby2];
+
+                adjncy[xadj[own] + nFacesPerCell[own]++] = nei;
+                adjncy[xadj[nei] + nFacesPerCell[nei]++] = own;
+            }
+        }
+    }
+
+    // Decompose using default weights
+    List<int> finalDecomp;
+    decompose(adjncy, xadj, finalDecomp);
+
+    // Copy back to labelList
+    labelList decomp(finalDecomp.size());
+    forAll(decomp, i)
+    {
+        decomp[i] = finalDecomp[i];
+    }
+    return decomp;
+}
+
+
+// From cell-cell connections to Metis format (like CompactListList)
+void Foam::scotchDecomp::calcMetisCSR
+(
+    const labelListList& cellCells,
+    List<int>& adjncy,
+    List<int>& xadj
+)
+{
+    // Count number of internal faces
+    label nConnections = 0;
+
+    forAll(cellCells, coarseI)
+    {
+        nConnections += cellCells[coarseI].size();
+    }
+
+    // Create the adjncy array as twice the size of the total number of
+    // internal faces
+    adjncy.setSize(nConnections);
+
+    xadj.setSize(cellCells.size()+1);
+
+
+    // Fill in xadj
+    // ~~~~~~~~~~~~
+    label freeAdj = 0;
+
+    forAll(cellCells, coarseI)
+    {
+        xadj[coarseI] = freeAdj;
+
+        const labelList& cCells = cellCells[coarseI];
+
+        forAll(cCells, i)
+        {
+            adjncy[freeAdj++] = cCells[i];
+        }
+    }
+    xadj[cellCells.size()] = freeAdj;
+}
+
+
+Foam::labelList Foam::scotchDecomp::decompose
+(
+    const labelList& agglom,
+    const pointField& agglomPoints
+)
+{
+    if (agglom.size() != mesh_.nCells())
+    {
+        FatalErrorIn
+        (
+            "parMetisDecomp::decompose(const labelList&, const pointField&)"
+        )   << "Size of cell-to-coarse map " << agglom.size()
+            << " differs from number of cells in mesh " << mesh_.nCells()
+            << exit(FatalError);
+    }
+
+    // Make Metis CSR (Compressed Storage Format) storage
+    //   adjncy      : contains neighbours (= edges in graph)
+    //   xadj(celli) : start of information in adjncy for celli
+    List<int> adjncy;
+    List<int> xadj;
+    {
+        // Get cellCells on coarse mesh.
+        labelListList cellCells;
+        calcCellCells
+        (
+            mesh_,
+            agglom,
+            agglomPoints.size(),
+            cellCells
+        );
+
+        calcMetisCSR(cellCells, adjncy, xadj);
+    }
+
+    // Decompose using default weights
+    List<int> finalDecomp;
+    decompose(adjncy, xadj, finalDecomp);
+
+
+    // Rework back into decomposition for original mesh_
+    labelList fineDistribution(agglom.size());
+
+    forAll(fineDistribution, i)
+    {
+        fineDistribution[i] = finalDecomp[agglom[i]];
+    }
+
+    return fineDistribution;
+}
+
+
+Foam::labelList Foam::scotchDecomp::decompose
+(
+    const labelListList& globalCellCells,
+    const pointField& cellCentres
+)
+{
+    if (cellCentres.size() != globalCellCells.size())
+    {
+        FatalErrorIn
+        (
+            "scotchDecomp::decompose(const pointField&, const labelListList&)"
+        )   << "Inconsistent number of cells (" << globalCellCells.size()
+            << ") and number of cell centres (" << cellCentres.size()
+            << ")." << exit(FatalError);
+    }
+
+
+    // Make Metis CSR (Compressed Storage Format) storage
+    //   adjncy      : contains neighbours (= edges in graph)
+    //   xadj(celli) : start of information in adjncy for celli
+
+    List<int> adjncy;
+    List<int> xadj;
+    calcMetisCSR(globalCellCells, adjncy, xadj);
+
+
+    // Decompose using default weights
+    List<int> finalDecomp;
+    decompose(adjncy, xadj, finalDecomp);
+
+    // Copy back to labelList
+    labelList decomp(finalDecomp.size());
+    forAll(decomp, i)
+    {
+        decomp[i] = finalDecomp[i];
+    }
+    return decomp;
+}
+
+
+// ************************************************************************* //

src/decompositionMethods/decompositionMethods/scotchDecomp/scotchDecomp.H

@@ -0,0 +1,150 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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
+    Foam::scotchDecomp
+
+Description
+    Scotch domain decomposition
+
+SourceFiles
+    scotchDecomp.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef scotchDecomp_H
+#define scotchDecomp_H
+
+#include "decompositionMethod.H"
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class scotchDecomp Declaration
+\*---------------------------------------------------------------------------*/
+
+class scotchDecomp
+:
+    public decompositionMethod
+{
+    // Private data
+
+        const polyMesh& mesh_;
+
+
+    // Private Member Functions
+
+        //- Check and print error message
+        static void check(const int, const char*);
+
+        label decompose
+        (
+            const List<int>& adjncy,
+            const List<int>& xadj,
+            List<int>& finalDecomp
+        );
+
+        //- Disallow default bitwise copy construct and assignment
+        void operator=(const scotchDecomp&);
+        scotchDecomp(const scotchDecomp&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("scotch");
+
+
+    // Constructors
+
+        //- Construct given the decomposition dictionary and mesh
+        scotchDecomp
+        (
+            const dictionary& decompositionDict,
+            const polyMesh& mesh
+        );
+
+
+    // Destructor
+
+        virtual ~scotchDecomp()
+        {}
+
+
+    // Member Functions
+
+        virtual bool parallelAware() const
+        {
+            // Metis does not know about proc boundaries
+            return false;
+        }
+
+        //- Return for every coordinate the wanted processor number. Use the
+        //  mesh connectivity (if needed)
+        virtual labelList decompose(const pointField&);
+
+        //- Return for every coordinate the wanted processor number. Gets
+        //  passed agglomeration map (from fine to coarse cells) and coarse cell
+        //  location. Can be overridden by decomposers that provide this
+        //  functionality natively.
+        virtual labelList decompose
+        (
+            const labelList& agglom,
+            const pointField&
+        );
+
+        //- Return for every coordinate the wanted processor number. Explicitly
+        //  provided mesh connectivity.
+        //  The connectivity is equal to mesh.cellCells() except for
+        //  - in parallel the cell numbers are global cell numbers (starting
+        //    from 0 at processor0 and then incrementing all through the
+        //    processors)
+        //  - the connections are across coupled patches
+        virtual labelList decompose
+        (
+            const labelListList& globalCellCells,
+            const pointField& cc
+        );
+
+        //- Helper to convert cellcells into compact row storage
+        static void calcMetisCSR
+        (
+            const labelListList& globalCellCells,
+            List<int>& adjncy,
+            List<int>& xadj
+        );
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/dynamicMesh/motionSmoother/motionSmoother.C

@@ -825,6 +825,28 @@ bool Foam::motionSmoother::scaleMesh
     const bool smoothMesh,
     const label nAllowableErrors
 )
+{
+    return scaleMesh
+    (
+        checkFaces,
+        baffles,
+        paramDict_,
+        paramDict_,
+        smoothMesh,
+        nAllowableErrors
+    );
+}
+
+
+bool Foam::motionSmoother::scaleMesh
+(
+    labelList& checkFaces,
+    const List<labelPair>& baffles,
+    const dictionary& paramDict,
+    const dictionary& meshQualityDict,
+    const bool smoothMesh,
+    const label nAllowableErrors
+)
 {
     if (!smoothMesh && adaptPatchIDs_.empty())
     {
@@ -859,9 +881,9 @@ bool Foam::motionSmoother::scaleMesh
     }
 
     const scalar errorReduction =
-        readScalar(paramDict_.lookup("errorReduction"));
+        readScalar(paramDict.lookup("errorReduction"));
     const label nSmoothScale =
-        readLabel(paramDict_.lookup("nSmoothScale"));
+        readLabel(paramDict.lookup("nSmoothScale"));
 
 
     // Note: displacement_ should already be synced already from setDisplacement
@@ -921,7 +943,7 @@ bool Foam::motionSmoother::scaleMesh
 
     // Check. Returns parallel number of incorrect faces.
     faceSet wrongFaces(mesh_, "wrongFaces", mesh_.nFaces()/100+100);
-    checkMesh(false, mesh_, paramDict_, checkFaces, baffles, wrongFaces);
+    checkMesh(false, mesh_, meshQualityDict, checkFaces, baffles, wrongFaces);
 
     if (returnReduce(wrongFaces.size(), sumOp<label>()) <= nAllowableErrors)
     {

src/dynamicMesh/motionSmoother/motionSmoother.H

@@ -52,7 +52,7 @@ Description
     @endverbatim
 
 Note
-    Shared points (parallel): a processor can have points which are part of
+    - Shared points (parallel): a processor can have points which are part of
     pp on another processor but have no pp itself (i.e. it has points
     and/or edges but no faces of pp). Hence we have to be careful when e.g.
     synchronising displacements that the value from the processor which has
@@ -61,10 +61,13 @@ Note
     else. The combine operator used will give preference to non-zero
     values.
 
-    Various routines take baffles. These are sets of boundary faces that
+    - Various routines take baffles. These are sets of boundary faces that
     are treated as a single internal face. This is a hack used to apply
     movement to internal faces.
 
+    - Mesh constraints are looked up from the supplied dictionary. (uses
+    recursive lookup)
+
 SourceFiles
     motionSmoother.C
     motionSmootherTemplates.C
@@ -395,6 +398,17 @@ public:
                 const label nAllow = 0
             );
 
+            //- Move mesh with externally provided mesh constraints
+            bool scaleMesh
+            (
+                labelList& checkFaces,
+                const List<labelPair>& baffles,
+                const dictionary& paramDict,
+                const dictionary& meshQualityDict,
+                const bool smoothMesh = true,
+                const label nAllow = 0
+            );
+
             //- Update topology
             void updateMesh();
 

src/dynamicMesh/motionSmoother/motionSmootherCheck.C

@@ -28,19 +28,6 @@ License
 #include "polyMeshGeometry.H"
 #include "IOmanip.H"
 
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 bool Foam::motionSmoother::checkMesh
@@ -74,17 +61,50 @@ bool Foam::motionSmoother::checkMesh
     labelHashSet& wrongFaces
 )
 {
-    const scalar maxNonOrtho(readScalar(dict.lookup("maxNonOrtho")));
-    const scalar minVol(readScalar(dict.lookup("minVol")));
-    const scalar maxConcave(readScalar(dict.lookup("maxConcave")));
-    const scalar minArea(readScalar(dict.lookup("minArea")));
-    const scalar maxIntSkew(readScalar(dict.lookup("maxInternalSkewness")));
-    const scalar maxBounSkew(readScalar(dict.lookup("maxBoundarySkewness")));
-    const scalar minWeight(readScalar(dict.lookup("minFaceWeight")));
-    const scalar minVolRatio(readScalar(dict.lookup("minVolRatio")));
-    const scalar minTwist(readScalar(dict.lookup("minTwist")));
-    const scalar minTriangleTwist(readScalar(dict.lookup("minTriangleTwist")));
-    const scalar minDet(readScalar(dict.lookup("minDeterminant")));
+    const scalar maxNonOrtho
+    (
+        readScalar(dict.lookup("maxNonOrtho", true))
+    );
+    const scalar minVol
+    (
+        readScalar(dict.lookup("minVol", true))
+    );
+    const scalar maxConcave
+    (
+        readScalar(dict.lookup("maxConcave", true))
+    );
+    const scalar minArea
+    (
+        readScalar(dict.lookup("minArea", true))
+    );
+    const scalar maxIntSkew
+    (
+        readScalar(dict.lookup("maxInternalSkewness", true))
+    );
+    const scalar maxBounSkew
+    (
+        readScalar(dict.lookup("maxBoundarySkewness", true))
+    );
+    const scalar minWeight
+    (
+        readScalar(dict.lookup("minFaceWeight", true))
+    );
+    const scalar minVolRatio
+    (
+        readScalar(dict.lookup("minVolRatio", true))
+    );
+    const scalar minTwist
+    (
+        readScalar(dict.lookup("minTwist", true))
+    );
+    const scalar minTriangleTwist
+    (
+        readScalar(dict.lookup("minTriangleTwist", true))
+    );
+    const scalar minDet
+    (
+        readScalar(dict.lookup("minDeterminant", true))
+    );
 
     label nWrongFaces = 0;
 
@@ -389,17 +409,50 @@ bool Foam::motionSmoother::checkMesh
     labelHashSet& wrongFaces
 )
 {
-    const scalar maxNonOrtho(readScalar(dict.lookup("maxNonOrtho")));
-    const scalar minVol(readScalar(dict.lookup("minVol")));
-    const scalar maxConcave(readScalar(dict.lookup("maxConcave")));
-    const scalar minArea(readScalar(dict.lookup("minArea")));
-    const scalar maxIntSkew(readScalar(dict.lookup("maxInternalSkewness")));
-    const scalar maxBounSkew(readScalar(dict.lookup("maxBoundarySkewness")));
-    const scalar minWeight(readScalar(dict.lookup("minFaceWeight")));
-    const scalar minVolRatio(readScalar(dict.lookup("minVolRatio")));
-    const scalar minTwist(readScalar(dict.lookup("minTwist")));
-    const scalar minTriangleTwist(readScalar(dict.lookup("minTriangleTwist")));
-    const scalar minDet(readScalar(dict.lookup("minDeterminant")));
+    const scalar maxNonOrtho
+    (
+        readScalar(dict.lookup("maxNonOrtho", true))
+    );
+    const scalar minVol
+    (
+        readScalar(dict.lookup("minVol", true))
+    );
+    const scalar maxConcave
+    (
+        readScalar(dict.lookup("maxConcave", true))
+    );
+    const scalar minArea
+    (
+        readScalar(dict.lookup("minArea", true))
+    );
+    const scalar maxIntSkew
+    (
+        readScalar(dict.lookup("maxInternalSkewness", true))
+    );
+    const scalar maxBounSkew
+    (
+        readScalar(dict.lookup("maxBoundarySkewness", true))
+    );
+    const scalar minWeight
+    (
+        readScalar(dict.lookup("minFaceWeight", true))
+    );
+    const scalar minVolRatio
+    (
+        readScalar(dict.lookup("minVolRatio", true))
+    );
+    const scalar minTwist
+    (
+        readScalar(dict.lookup("minTwist", true))
+    );
+    const scalar minTriangleTwist
+    (
+        readScalar(dict.lookup("minTriangleTwist", true))
+    );
+    const scalar minDet
+    (
+        readScalar(dict.lookup("minDeterminant", true))
+    );
 
     label nWrongFaces = 0;
 

src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.H

@@ -294,9 +294,10 @@ public:
             //  - nPointLayers : number of layers per (patch)point
             //  - nFaceLayers : number of layers per (patch) face
             //  - firstDisplacement : displacement per point for first
-            //    layer of points. If zero do not add point.
+            //    layer of points (i.e. nearest to original mesh). If zero
+            //    do not add point.
             //  Layer thicknesses are calculated to constant geometric
-            //  expansion. Use expansionRatio for constant size.
+            //  expansion. Use expansionRatio 1 for constant size.
             //  Sets addedPoints_ which is per pp point a list of points
             //  added.
             //  Note: firstDisplacement has to be parallel synchronised before

src/finiteVolume/Make/files

@@ -24,6 +24,7 @@ $(constraintFvPatches)/processor/processorFvPatch.C
 derivedFvPatches = $(fvPatches)/derived
 $(derivedFvPatches)/wall/wallFvPatch.C
 $(derivedFvPatches)/directMapped/directMappedFvPatch.C
+$(derivedFvPatches)/directMapped/directMappedWallFvPatch.C
 
 wallDist = fvMesh/wallDist
 $(wallDist)/wallPointYPlus/wallPointYPlus.C

src/finiteVolume/fields/fvPatchFields/derived/directMappedFixedValue/directMappedFixedValueFvPatchField.C

@@ -25,7 +25,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "directMappedFixedValueFvPatchField.H"
-#include "directMappedFvPatch.H"
+#include "directMappedPatchBase.H"
 #include "volFields.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -61,7 +61,7 @@ directMappedFixedValueFvPatchField<Type>::directMappedFixedValueFvPatchField
     setAverage_(ptf.setAverage_),
     average_(ptf.average_)
 {
-    if (!isType<directMappedFvPatch>(this->patch()))
+    if (!isA<directMappedPatchBase>(this->patch().patch()))
     {
         FatalErrorIn
         (
@@ -74,7 +74,7 @@ directMappedFixedValueFvPatchField<Type>::directMappedFixedValueFvPatchField
             "    const fvPatchFieldMapper&\n"
             ")\n"
         )   << "\n    patch type '" << p.type()
-            << "' not type '" << typeName << "'"
+            << "' not type '" << directMappedPatchBase::typeName << "'"
             << "\n    for patch " << p.name()
             << " of field " << this->dimensionedInternalField().name()
             << " in file " << this->dimensionedInternalField().objectPath()
@@ -95,7 +95,7 @@ directMappedFixedValueFvPatchField<Type>::directMappedFixedValueFvPatchField
     setAverage_(readBool(dict.lookup("setAverage"))),
     average_(pTraits<Type>(dict.lookup("average")))
 {
-    if (!isType<directMappedFvPatch>(this->patch()))
+    if (!isA<directMappedPatchBase>(this->patch().patch()))
     {
         FatalErrorIn
         (
@@ -107,12 +107,19 @@ directMappedFixedValueFvPatchField<Type>::directMappedFixedValueFvPatchField
             "    const dictionary& dict\n"
             ")\n"
         )   << "\n    patch type '" << p.type()
-            << "' not type '" << typeName << "'"
+            << "' not type '" << directMappedPatchBase::typeName << "'"
             << "\n    for patch " << p.name()
             << " of field " << this->dimensionedInternalField().name()
             << " in file " << this->dimensionedInternalField().objectPath()
             << exit(FatalError);
     }
+
+    // Force calculation of schedule (uses parallel comms)
+    const directMappedPatchBase& mpp = refCast<const directMappedPatchBase>
+    (
+        this->patch().patch()
+    );
+    (void)mpp.map().schedule();
 }
 
 
@@ -143,70 +150,6 @@ directMappedFixedValueFvPatchField<Type>::directMappedFixedValueFvPatchField
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class Type>
-void directMappedFixedValueFvPatchField<Type>::getNewValues
-(
-    const directMappedPolyPatch& mpp,
-    const Field<Type>& sendValues,
-    Field<Type>& newValues
-) const
-{
-    // Get the scheduling information
-    const List<labelPair>& schedule = mpp.schedule();
-    const labelListList& sendLabels = mpp.sendLabels();
-    const labelListList& receiveFaceLabels = mpp.receiveFaceLabels();
-
-    forAll(schedule, i)
-    {
-        const labelPair& twoProcs = schedule[i];
-        label sendProc = twoProcs[0];
-        label recvProc = twoProcs[1];
-
-        if (Pstream::myProcNo() == sendProc)
-        {
-            OPstream toProc(Pstream::scheduled, recvProc);
-            toProc<< UIndirectList<Type>(sendValues, sendLabels[recvProc])();
-        }
-        else
-        {
-            // I am receiver. Receive from sendProc.
-            IPstream fromProc(Pstream::scheduled, sendProc);
-
-            Field<Type> fromFld(fromProc);
-
-            // Destination faces
-            const labelList& faceLabels = receiveFaceLabels[sendProc];
-
-            forAll(fromFld, i)
-            {
-                label patchFaceI = faceLabels[i];
-
-                newValues[patchFaceI] = fromFld[i];
-            }
-        }
-    }
-
-    // Do data from myself
-    {
-        UIndirectList<Type> fromFld
-        (
-            sendValues,
-            sendLabels[Pstream::myProcNo()]
-        );
-
-        // Destination faces
-        const labelList& faceLabels = receiveFaceLabels[Pstream::myProcNo()];
-
-        forAll(fromFld, i)
-        {
-            label patchFaceI = faceLabels[i];
-
-            newValues[patchFaceI] = fromFld[i];
-        }
-    }
-}
-
-
 template<class Type>
 void directMappedFixedValueFvPatchField<Type>::updateCoeffs()
 {
@@ -215,65 +158,94 @@ void directMappedFixedValueFvPatchField<Type>::updateCoeffs()
         return;
     }
 
-    // Get the directMappedPolyPatch
-    const directMappedPolyPatch& mpp = refCast<const directMappedPolyPatch>
+    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
+
+    // Get the scheduling information from the directMappedPatchBase
+    const directMappedPatchBase& mpp = refCast<const directMappedPatchBase>
     (
         directMappedFixedValueFvPatchField<Type>::patch().patch()
     );
+    const mapDistribute& distMap = mpp.map();
+    const fvMesh& nbrMesh = refCast<const fvMesh>(mpp.sampleMesh());
+    const word& fldName = this->dimensionedInternalField().name();
 
-    Field<Type> newValues(this->size());
+    // Result of obtaining remote values
+    Field<Type> newValues;
 
     switch (mpp.mode())
     {
-        case directMappedPolyPatch::NEARESTCELL:
+        case directMappedPatchBase::NEARESTCELL:
         {
-            getNewValues(mpp, this->internalField(), newValues);
+            if (mpp.sameRegion())
+            {
+                newValues = this->internalField();
+            }
+            else
+            {
+                newValues = nbrMesh.lookupObject<fieldType>
+                (
+                    fldName
+                ).internalField();
+            }
+            mapDistribute::distribute
+            (
+                Pstream::defaultCommsType,
+                distMap.schedule(),
+                distMap.constructSize(),
+                distMap.subMap(),
+                distMap.constructMap(),
+                newValues
+            );
 
             break;
         }
-        case directMappedPolyPatch::NEARESTPATCHFACE:
+        case directMappedPatchBase::NEARESTPATCHFACE:
         {
-            const label patchID =
-                 this->patch().patch().boundaryMesh().findPatchID
-                 (
-                    mpp.samplePatch()
-                 );
-            if (patchID < 0)
+            const label nbrPatchID = nbrMesh.boundaryMesh().findPatchID
+            (
+                mpp.samplePatch()
+            );
+            if (nbrPatchID < 0)
             {
                 FatalErrorIn
                 (
                     "void directMappedFixedValueFvPatchField<Type>::"
                     "updateCoeffs()"
                 )<< "Unable to find sample patch " << mpp.samplePatch()
+                 << " in region " << mpp.sampleRegion()
                  << " for patch " << this->patch().name() << nl
                  << abort(FatalError);
             }
-            typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
-            const word& fieldName = this->dimensionedInternalField().name();
-            const fieldType& sendField =
-                this->db().objectRegistry::lookupObject<fieldType>(fieldName);
 
-            getNewValues(mpp, sendField.boundaryField()[patchID], newValues);
+            const fieldType& nbrField = nbrMesh.lookupObject<fieldType>
+            (
+                fldName
+            );
+            newValues = nbrField.boundaryField()[nbrPatchID];
+            mapDistribute::distribute
+            (
+                Pstream::defaultCommsType,
+                distMap.schedule(),
+                distMap.constructSize(),
+                distMap.subMap(),
+                distMap.constructMap(),
+                newValues
+            );
 
             break;
         }
-        case directMappedPolyPatch::NEARESTFACE:
+        case directMappedPatchBase::NEARESTFACE:
         {
-            typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
-            const word& fieldName = this->dimensionedInternalField().name();
-            const fieldType& sendField =
-                this->db().objectRegistry::lookupObject<fieldType>(fieldName);
+            Field<Type> allValues(nbrMesh.nFaces(), pTraits<Type>::zero);
 
-            Field<Type> allValues
+            const fieldType& nbrField = nbrMesh.lookupObject<fieldType>
             (
-                this->patch().patch().boundaryMesh().mesh().nFaces(),
-                pTraits<Type>::zero
+                fldName
             );
-
-            forAll(sendField.boundaryField(), patchI)
+            forAll(nbrField.boundaryField(), patchI)
             {
                 const fvPatchField<Type>& pf =
-                    sendField.boundaryField()[patchI];
+                    nbrField.boundaryField()[patchI];
                 label faceStart = pf.patch().patch().start();
 
                 forAll(pf, faceI)
@@ -282,9 +254,17 @@ void directMappedFixedValueFvPatchField<Type>::updateCoeffs()
                 }
             }
 
-            getNewValues(mpp, allValues, newValues);
+            mapDistribute::distribute
+            (
+                Pstream::defaultCommsType,
+                distMap.schedule(),
+                distMap.constructSize(),
+                distMap.subMap(),
+                distMap.constructMap(),
+                allValues
+            );
 
-            newValues = this->patch().patchSlice(newValues);
+            newValues = this->patch().patchSlice(allValues);
 
             break;
         }
@@ -316,6 +296,16 @@ void directMappedFixedValueFvPatchField<Type>::updateCoeffs()
 
     this->operator==(newValues);
 
+    if (debug)
+    {
+        Info<< "directMapped on field:" << fldName
+            << " patch:" << this->patch().name()
+            << "  avg:" << gAverage(*this)
+            << "  min:" << gMin(*this)
+            << "  max:" << gMax(*this)
+            << endl;
+    }
+
     fixedValueFvPatchField<Type>::updateCoeffs();
 }
 

src/finiteVolume/fields/fvPatchFields/derived/directMappedFixedValue/directMappedFixedValueFvPatchField.H

@@ -62,18 +62,6 @@ class directMappedFixedValueFvPatchField
         //  setAverage_ is set true
         Type average_;
 
-
-   // Private member functions
-
-        // Helper function to return the new field values
-        void getNewValues
-        (
-            const directMappedPolyPatch& mpp,
-            const Field<Type>& sendValues,
-            Field<Type>& newValues
-        ) const;
-
-
 public:
 
     //- Runtime type information

src/finiteVolume/fields/fvPatchFields/derived/directMappedVelocityFluxFixedValue/directMappedVelocityFluxFixedValueFvPatchField.C

@@ -26,7 +26,7 @@ License
 
 #include "directMappedVelocityFluxFixedValueFvPatchField.H"
 #include "fvPatchFieldMapper.H"
-#include "directMappedFvPatch.H"
+#include "directMappedPatchBase.H"
 #include "volFields.H"
 #include "surfaceFields.H"
 #include "addToRunTimeSelectionTable.H"
@@ -62,12 +62,12 @@ directMappedVelocityFluxFixedValueFvPatchField
     fixedValueFvPatchVectorField(ptf, p, iF, mapper),
     phiName_(ptf.phiName_)
 {
-    if (!isType<directMappedFvPatch>(patch()))
+    if (!isType<directMappedPatchBase>(this->patch().patch()))
     {
         FatalErrorIn
         (
             "directMappedVelocityFluxFixedValueFvPatchField::"
-            "directMappedFixedValueFvPatchField\n"
+            "directMappedVelocityFluxFixedValueFvPatchField\n"
             "(\n"
             "    const directMappedVelocityFluxFixedValueFvPatchField&,\n"
             "    const fvPatch&,\n"
@@ -75,7 +75,7 @@ directMappedVelocityFluxFixedValueFvPatchField
             "    const fvPatchFieldMapper&\n"
             ")\n"
         )   << "\n    patch type '" << p.type()
-            << "' not type '" << typeName << "'"
+            << "' not type '" << directMappedPatchBase::typeName << "'"
             << "\n    for patch " << p.name()
             << " of field " << dimensionedInternalField().name()
             << " in file " << dimensionedInternalField().objectPath()
@@ -95,24 +95,31 @@ directMappedVelocityFluxFixedValueFvPatchField
     fixedValueFvPatchVectorField(p, iF, dict),
     phiName_(dict.lookup("phi"))
 {
-    if (!isType<directMappedFvPatch>(patch()))
+    if (!isType<directMappedPatchBase>(this->patch().patch()))
     {
         FatalErrorIn
         (
             "directMappedVelocityFluxFixedValueFvPatchField::"
-            "directMappedFixedValueFvPatchField\n"
+            "directMappedVelocityFluxFixedValueFvPatchField\n"
             "(\n"
             "    const fvPatch& p,\n"
             "    const DimensionedField<vector, volMesh>& iF,\n"
             "    const dictionary& dict\n"
             ")\n"
         )   << "\n    patch type '" << p.type()
-            << "' not type '" << typeName << "'"
+            << "' not type '" << directMappedPatchBase::typeName << "'"
             << "\n    for patch " << p.name()
             << " of field " << dimensionedInternalField().name()
             << " in file " << dimensionedInternalField().objectPath()
             << exit(FatalError);
     }
+
+    // Force calculation of schedule (uses parallel comms)
+    const directMappedPolyPatch& mpp = refCast<const directMappedPolyPatch>
+    (
+        patch().patch()
+    );
+    (void)mpp.map().schedule();
 }
 
 
@@ -139,85 +146,6 @@ directMappedVelocityFluxFixedValueFvPatchField
 {}
 
 
-void directMappedVelocityFluxFixedValueFvPatchField::getNewValues
-(
-    const directMappedPolyPatch& mpp,
-    const vectorField& sendUValues,
-    const scalarField& sendPhiValues,
-    vectorField& newUValues,
-    scalarField& newPhiValues
-) const
-{
-    // Get the scheduling information
-    const List<labelPair>& schedule = mpp.schedule();
-    const labelListList& sendLabels = mpp.sendLabels();
-    const labelListList& receiveFaceLabels = mpp.receiveFaceLabels();
-
-    forAll(schedule, i)
-    {
-        const labelPair& twoProcs = schedule[i];
-        label sendProc = twoProcs[0];
-        label recvProc = twoProcs[1];
-
-        if (Pstream::myProcNo() == sendProc)
-        {
-            OPstream toProc(Pstream::scheduled, recvProc);
-            toProc<< UIndirectList<vector>(sendUValues, sendLabels[recvProc])();
-            toProc<< UIndirectList<scalar>
-                (
-                    sendPhiValues,
-                    sendLabels[recvProc]
-                )();
-        }
-        else
-        {
-            // I am receiver. Receive from sendProc.
-            IPstream fromProc(Pstream::scheduled, sendProc);
-
-            vectorField fromUFld(fromProc);
-            scalarField fromPhiFld(fromProc);
-
-            // Destination faces
-            const labelList& faceLabels = receiveFaceLabels[sendProc];
-
-            forAll(fromUFld, i)
-            {
-                label patchFaceI = faceLabels[i];
-
-                newUValues[patchFaceI] = fromUFld[i];
-                newPhiValues[patchFaceI] = fromPhiFld[i];
-            }
-        }
-    }
-
-    // Do data from myself
-    {
-        UIndirectList<vector> fromUFld
-            (
-                sendUValues,
-                sendLabels[Pstream::myProcNo()]
-            );
-
-        UIndirectList<scalar> fromPhiFld
-            (
-                sendPhiValues,
-                sendLabels[Pstream::myProcNo()]
-            );
-
-        // Destination faces
-        const labelList& faceLabels = receiveFaceLabels[Pstream::myProcNo()];
-
-        forAll(fromUFld, i)
-        {
-            label patchFaceI = faceLabels[i];
-
-            newUValues[patchFaceI] = fromUFld[i];
-            newPhiValues[patchFaceI] = fromPhiFld[i];
-        }
-    }
-}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 void directMappedVelocityFluxFixedValueFvPatchField::updateCoeffs()
@@ -227,37 +155,33 @@ void directMappedVelocityFluxFixedValueFvPatchField::updateCoeffs()
         return;
     }
 
-    // Get the directMappedPolyPatch
-    const directMappedPolyPatch& mpp = refCast<const directMappedPolyPatch>
+    // Get the directMappedPatchBase
+    const directMappedPatchBase& mpp = refCast<const directMappedPatchBase>
     (
         directMappedVelocityFluxFixedValueFvPatchField::patch().patch()
     );
-
-    vectorField newUValues(size());
-    scalarField newPhiValues(size());
-
+    const mapDistribute& distMap = mpp.map();
+    const fvMesh& nbrMesh = refCast<const fvMesh>(mpp.sampleMesh());
     const word& fieldName = dimensionedInternalField().name();
-    const volVectorField& UField = db().lookupObject<volVectorField>(fieldName);
+    const volVectorField& UField = nbrMesh.lookupObject<volVectorField>
+    (
+        fieldName
+    );
 
     surfaceScalarField& phiField = const_cast<surfaceScalarField&>
     (
-        db().lookupObject<surfaceScalarField>(phiName_)
+        nbrMesh.lookupObject<surfaceScalarField>(phiName_)
     );
 
+    vectorField newUValues;
+    scalarField newPhiValues;
+
     switch (mpp.mode())
     {
         case directMappedPolyPatch::NEARESTFACE:
         {
-            vectorField allUValues
-            (
-                patch().patch().boundaryMesh().mesh().nFaces(),
-                vector::zero
-            );
-            scalarField allPhiValues
-            (
-                patch().patch().boundaryMesh().mesh().nFaces(),
-                0.0
-            );
+            vectorField allUValues(nbrMesh.nFaces(), vector::zero);
+            scalarField allPhiValues(nbrMesh.nFaces(), 0.0);
 
             forAll(UField.boundaryField(), patchI)
             {
@@ -273,34 +197,58 @@ void directMappedVelocityFluxFixedValueFvPatchField::updateCoeffs()
                 }
             }
 
-            getNewValues
+            mapDistribute::distribute
             (
-                mpp,
-                allUValues,
-                allPhiValues,
-                newUValues,
+                Pstream::defaultCommsType,
+                distMap.schedule(),
+                distMap.constructSize(),
+                distMap.subMap(),
+                distMap.constructMap(),
+                allUValues
+            );
+            newUValues = patch().patchSlice(newUValues);
+            
+            mapDistribute::distribute
+            (
+                Pstream::defaultCommsType,
+                distMap.schedule(),
+                distMap.constructSize(),
+                distMap.subMap(),
+                distMap.constructMap(),
                 newPhiValues
             );
-
-            newUValues = patch().patchSlice(newUValues);
             newPhiValues = patch().patchSlice(newPhiValues);
 
             break;
         }
         case directMappedPolyPatch::NEARESTPATCHFACE:
         {
-            const label patchID =
-                patch().patch().boundaryMesh().findPatchID
-                (
-                    mpp.samplePatch()
-                );
-
-            getNewValues
+            const label nbrPatchID = nbrMesh.boundaryMesh().findPatchID
             (
-                mpp,
-                UField.boundaryField()[patchID],
-                phiField.boundaryField()[patchID],
-                newUValues,
+                mpp.samplePatch()
+            );
+
+            newUValues = UField.boundaryField()[nbrPatchID];
+
+            mapDistribute::distribute
+            (
+                Pstream::defaultCommsType,
+                distMap.schedule(),
+                distMap.constructSize(),
+                distMap.subMap(),
+                distMap.constructMap(),
+                newUValues
+            );
+
+            newPhiValues = phiField.boundaryField()[nbrPatchID];
+
+            mapDistribute::distribute
+            (
+                Pstream::defaultCommsType,
+                distMap.schedule(),
+                distMap.constructSize(),
+                distMap.subMap(),
+                distMap.constructMap(),
                 newPhiValues
             );
 

src/finiteVolume/fields/fvPatchFields/derived/directMappedVelocityFluxFixedValue/directMappedVelocityFluxFixedValueFvPatchField.H

@@ -57,20 +57,6 @@ class directMappedVelocityFluxFixedValueFvPatchField
         //- Name of flux field
         word phiName_;
 
-
-   // Private member functions
-
-        // Helper function to return the new field values
-        void getNewValues
-        (
-            const directMappedPolyPatch& mpp,
-            const vectorField& sendUValues,
-            const scalarField& sendPhiValues,
-            vectorField& newUValues,
-            scalarField& newPhiValues
-        ) const;
-
-
 public:
 
     //- Runtime type information

src/finiteVolume/fields/fvPatchFields/derived/inletOutletTotalTemperature/inletOutletTotalTemperatureFvPatchScalarField.C

@@ -189,7 +189,8 @@ void Foam::inletOutletTotalTemperatureFvPatchScalarField::updateCoeffs()
 }
 
 
-void Foam::inletOutletTotalTemperatureFvPatchScalarField::write(Ostream& os) const
+void Foam::inletOutletTotalTemperatureFvPatchScalarField::write(Ostream& os)
+const
 {
     fvPatchScalarField::write(os);
     if (UName_ != "U")

src/finiteVolume/fvMesh/fvPatches/derived/directMapped/directMappedWallFvPatch.C

@@ -0,0 +1,38 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 "directMappedWallFvPatch.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(directMappedWallFvPatch, 0);
+    addToRunTimeSelectionTable(fvPatch, directMappedWallFvPatch, polyPatch);
+}
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/fvPatches/derived/directMapped/directMappedWallFvPatch.H

@@ -0,0 +1,84 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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
+    Foam::directMappedWallFvPatch
+
+Description
+    Foam::directMappedWallFvPatch
+
+SourceFiles
+    directMappedWallFvPatch.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef directMappedWallFvPatch_H
+#define directMappedWallFvPatch_H
+
+#include "fvPatch.H"
+#include "directMappedWallPolyPatch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class directMappedWallFvPatch Declaration
+\*---------------------------------------------------------------------------*/
+
+class directMappedWallFvPatch
+:
+    public fvPatch
+{
+
+public:
+
+    //- Runtime type information
+    TypeName(directMappedWallPolyPatch::typeName_());
+
+
+    // Constructors
+
+        //- Construct from components
+        directMappedWallFvPatch
+        (
+            const polyPatch& patch,
+            const fvBoundaryMesh& bm
+        )
+        :
+            fvPatch(patch, bm)
+        {}
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/meshTools/Make/files

@@ -141,7 +141,10 @@ triSurface/triSurfaceTools/geompack/geompack.C
 
 twoDPointCorrector/twoDPointCorrector.C
 
+directMapped/directMappedPolyPatch/directMappedPatchBase.C
 directMapped/directMappedPolyPatch/directMappedPolyPatch.C
+directMapped/directMappedPolyPatch/directMappedWallPolyPatch.C
 directMapped/directMappedPointPatch/directMappedPointPatch.C
+directMapped/directMappedPointPatch/directMappedWallPointPatch.C
 
 LIB = $(FOAM_LIBBIN)/libmeshTools

src/meshTools/directMapped/directMappedPointPatch/directMappedWallPointPatch.C

@@ -0,0 +1,52 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 "directMappedWallPointPatch.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(directMappedWallPointPatch, 0);
+
+// Add the patch constructor functions to the hash tables
+addToRunTimeSelectionTable
+(
+    facePointPatch,
+    directMappedWallPointPatch,
+    polyPatch
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/meshTools/directMapped/directMappedPointPatch/directMappedWallPointPatch.H

@@ -0,0 +1,84 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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
+    Foam::directMappedWallPointPatch
+
+Description
+    DirectMapped patch.
+
+SourceFiles
+    directMappedWallPointPatch.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef directMappedWallPointPatch_H
+#define directMappedWallPointPatch_H
+
+#include "facePointPatch.H"
+#include "directMappedWallPolyPatch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class directMappedWallPointPatch Declaration
+\*---------------------------------------------------------------------------*/
+
+class directMappedWallPointPatch
+:
+    public facePointPatch
+{
+
+public:
+
+    //- Runtime type information
+    TypeName(directMappedWallPolyPatch::typeName_());
+
+
+    // Constructors
+
+        //- Construct from polyPatch
+        directMappedWallPointPatch
+        (
+            const polyPatch& patch,
+            const pointBoundaryMesh& bm
+        )
+        :
+            facePointPatch(patch, bm)
+        {}
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/meshTools/directMapped/directMappedPolyPatch/directMappedPatchBase.C

@@ -0,0 +1,657 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 "directMappedPatchBase.H"
+#include "addToRunTimeSelectionTable.H"
+#include "ListListOps.H"
+#include "meshSearch.H"
+#include "meshTools.H"
+#include "OFstream.H"
+#include "Random.H"
+#include "treeDataFace.H"
+#include "indexedOctree.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(directMappedPatchBase, 0);
+
+    template<>
+    const char* NamedEnum<directMappedPatchBase::sampleMode, 3>::names[] =
+    {
+        "nearestCell",
+        "nearestPatchFace",
+        "nearestFace"
+    };
+
+    const NamedEnum<directMappedPatchBase::sampleMode, 3>
+        directMappedPatchBase::sampleModeNames_;
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::directMappedPatchBase::collectSamples
+(
+    pointField& samples,
+    labelList& patchFaceProcs,
+    labelList& patchFaces,
+    pointField& patchFc
+) const
+{
+
+    // Collect all sample points and the faces they come from.
+    List<pointField> globalFc(Pstream::nProcs());
+    List<pointField> globalSamples(Pstream::nProcs());
+    labelListList globalFaces(Pstream::nProcs());
+
+    globalFc[Pstream::myProcNo()] = patch_.faceCentres();
+    globalSamples[Pstream::myProcNo()] = globalFc[Pstream::myProcNo()]+offset_;
+    globalFaces[Pstream::myProcNo()] = identity(patch_.size());
+
+    // Distribute to all processors
+    Pstream::gatherList(globalSamples);
+    Pstream::scatterList(globalSamples);
+    Pstream::gatherList(globalFaces);
+    Pstream::scatterList(globalFaces);
+    Pstream::gatherList(globalFc);
+    Pstream::scatterList(globalFc);
+
+    // Rework into straight list
+    samples = ListListOps::combine<pointField>
+    (
+        globalSamples,
+        accessOp<pointField>()
+    );
+    patchFaces = ListListOps::combine<labelList>
+    (
+        globalFaces,
+        accessOp<labelList>()
+    );
+    patchFc = ListListOps::combine<pointField>
+    (
+        globalFc,
+        accessOp<pointField>()
+    );
+
+    patchFaceProcs.setSize(patchFaces.size());
+    labelList nPerProc
+    (
+        ListListOps::subSizes
+        (
+            globalFaces,
+            accessOp<labelList>()
+        )
+    );
+    label sampleI = 0;
+    forAll(nPerProc, procI)
+    {
+        for (label i = 0; i < nPerProc[procI]; i++)
+        {
+            patchFaceProcs[sampleI++] = procI;
+        }
+    }
+}
+
+
+// Find the processor/cell containing the samples. Does not account
+// for samples being found in two processors.
+void Foam::directMappedPatchBase::findSamples
+(
+    const pointField& samples,
+    labelList& sampleProcs,
+    labelList& sampleIndices,
+    pointField& sampleLocations
+) const
+{
+    // Lookup the correct region
+    const polyMesh& mesh = sampleMesh();
+
+    // All the info for nearest. Construct to miss
+    List<nearInfo> nearest(samples.size());
+
+    switch (mode_)
+    {
+        case NEARESTCELL:
+        {
+            if (samplePatch_.size() && samplePatch_ != "none")
+            {
+                FatalErrorIn
+                (
+                    "directMappedPatchBase::findSamples(const pointField&,"
+                    " labelList&, labelList&, pointField&) const"
+                )   << "No need to supply a patch name when in "
+                    << sampleModeNames_[mode_] << " mode." << exit(FatalError);
+            }
+
+            // Octree based search engine
+            meshSearch meshSearchEngine(mesh, false);
+
+            forAll(samples, sampleI)
+            {
+                const point& sample = samples[sampleI];
+
+                label cellI = meshSearchEngine.findCell(sample);
+
+                if (cellI == -1)
+                {
+                    nearest[sampleI].second().first() = Foam::sqr(GREAT);
+                    nearest[sampleI].second().second() = Pstream::myProcNo();
+                }
+                else
+                {
+                    const point& cc = mesh.cellCentres()[cellI];
+
+                    nearest[sampleI].first() = pointIndexHit
+                    (
+                        true,
+                        cc,
+                        cellI
+                    );
+                    nearest[sampleI].second().first() = magSqr(cc-sample);
+                    nearest[sampleI].second().second() = Pstream::myProcNo();
+                }
+            }
+            break;
+        }
+
+        case NEARESTPATCHFACE:
+        {
+            Random rndGen(123456);
+
+            const polyPatch& pp = samplePolyPatch();
+
+            if (pp.empty())
+            {
+                forAll(samples, sampleI)
+                {
+                    nearest[sampleI].second().first() = Foam::sqr(GREAT);
+                    nearest[sampleI].second().second() = Pstream::myProcNo();
+                }
+            }
+            else
+            {
+                // patch faces
+                const labelList patchFaces(identity(pp.size()) + pp.start());
+
+                treeBoundBox patchBb
+                (
+                    treeBoundBox(pp.points(), pp.meshPoints()).extend
+                    (
+                        rndGen,
+                        1E-4
+                    )
+                );
+                patchBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+                patchBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+
+                indexedOctree<treeDataFace> boundaryTree
+                (
+                    treeDataFace    // all information needed to search faces
+                    (
+                        false,      // do not cache bb
+                        mesh,
+                        patchFaces  // boundary faces only
+                    ),
+                    patchBb,        // overall search domain
+                    8,              // maxLevel
+                    10,             // leafsize
+                    3.0             // duplicity
+                );
+
+                forAll(samples, sampleI)
+                {
+                    const point& sample = samples[sampleI];
+
+                    pointIndexHit& nearInfo = nearest[sampleI].first();
+                    nearInfo = boundaryTree.findNearest
+                    (
+                        sample,
+                        magSqr(patchBb.span())
+                    );
+
+                    if (!nearInfo.hit())
+                    {
+                        nearest[sampleI].second().first() = Foam::sqr(GREAT);
+                        nearest[sampleI].second().second() =
+                            Pstream::myProcNo();
+                    }
+                    else
+                    {
+                        point fc(pp[nearInfo.index()].centre(pp.points()));
+                        nearInfo.setPoint(fc);
+                        nearest[sampleI].second().first() = magSqr(fc-sample);
+                        nearest[sampleI].second().second() =
+                            Pstream::myProcNo();
+                    }
+                }
+            }
+            break;
+        }
+
+        case NEARESTFACE:
+        {
+            if (samplePatch_.size() && samplePatch_ != "none")
+            {
+                FatalErrorIn
+                (
+                    "directMappedPatchBase::findSamples(const pointField&,"
+                    " labelList&, labelList&, pointField&) const"
+                )   << "No need to supply a patch name when in "
+                    << sampleModeNames_[mode_] << " mode." << exit(FatalError);
+            }
+
+            // Octree based search engine
+            meshSearch meshSearchEngine(mesh, false);
+
+            forAll(samples, sampleI)
+            {
+                const point& sample = samples[sampleI];
+
+                label faceI = meshSearchEngine.findNearestFace(sample);
+
+                if (faceI == -1)
+                {
+                    nearest[sampleI].second().first() = Foam::sqr(GREAT);
+                    nearest[sampleI].second().second() = Pstream::myProcNo();
+                }
+                else
+                {
+                    const point& fc = mesh.faceCentres()[faceI];
+
+                    nearest[sampleI].first() = pointIndexHit
+                    (
+                        true,
+                        fc,
+                        faceI
+                    );
+                    nearest[sampleI].second().first() = magSqr(fc-sample);
+                    nearest[sampleI].second().second() = Pstream::myProcNo();
+                }
+            }
+            break;
+        }
+
+        default:
+        {
+            FatalErrorIn("directMappedPatchBase::findSamples(..)")
+                << "problem." << abort(FatalError);
+        }
+    }
+
+
+    // Find nearest.
+    Pstream::listCombineGather(nearest, nearestEqOp());
+    Pstream::listCombineScatter(nearest);
+
+    if (debug)
+    {
+        Info<< "directMappedPatchBase::findSamples on mesh " << sampleRegion_
+            << " : " << endl;
+        forAll(nearest, sampleI)
+        {
+            label procI = nearest[sampleI].second().second();
+            label localI = nearest[sampleI].first().index();
+
+            Info<< "    " << sampleI << " coord:"<< samples[sampleI]
+                << " found on processor:" << procI
+                << " in local cell/face:" << localI
+                << " with cc:" << nearest[sampleI].first().rawPoint() << endl;
+        }
+    }
+
+    // Check for samples not being found
+    forAll(nearest, sampleI)
+    {
+        if (!nearest[sampleI].first().hit())
+        {
+            FatalErrorIn
+            (
+                "directMappedPatchBase::findSamples"
+                "(const pointField&, labelList&"
+                ", labelList&, pointField&)"
+            )   << "Did not find sample " << samples[sampleI]
+                << " on any processor of region" << sampleRegion_
+                << exit(FatalError);
+        }
+    }
+
+
+    // Convert back into proc+local index
+    sampleProcs.setSize(samples.size());
+    sampleIndices.setSize(samples.size());
+    sampleLocations.setSize(samples.size());
+
+    forAll(nearest, sampleI)
+    {
+        sampleProcs[sampleI] = nearest[sampleI].second().second();
+        sampleIndices[sampleI] = nearest[sampleI].first().index();
+        sampleLocations[sampleI] = nearest[sampleI].first().hitPoint();
+    }
+}
+
+
+void Foam::directMappedPatchBase::calcMapping() const
+{
+    if (mapPtr_.valid())
+    {
+        FatalErrorIn("directMappedPatchBase::calcMapping() const")
+            << "Mapping already calculated" << exit(FatalError);
+    }
+
+    if
+    (
+        offset_ == vector::zero
+     && sampleRegion_ == patch_.boundaryMesh().mesh().name()
+    )
+    {
+        FatalErrorIn("directMappedPatchBase::calcMapping() const")
+            << "Invalid offset " << offset_ << endl
+            << "Offset is the vector added to the patch face centres to"
+            << " find the cell supplying the data."
+            << exit(FatalError);
+    }
+
+
+    // Get global list of all samples and the processor and face they come from.
+    pointField samples;
+    labelList patchFaceProcs;
+    labelList patchFaces;
+    pointField patchFc;
+    collectSamples(samples, patchFaceProcs, patchFaces, patchFc);
+
+    // Find processor and cell/face samples are in and actual location.
+    labelList sampleProcs;
+    labelList sampleIndices;
+    pointField sampleLocations;
+    findSamples(samples, sampleProcs, sampleIndices, sampleLocations);
+
+
+    // Now we have all the data we need:
+    // - where sample originates from (so destination when mapping):
+    //   patchFaces, patchFaceProcs.
+    // - cell/face sample is in (so source when mapping)
+    //   sampleIndices, sampleProcs.
+
+    //forAll(samples, i)
+    //{
+    //    Info<< i << " need data in region "
+    //        << patch_.boundaryMesh().mesh().name()
+    //        << " for proc:" << patchFaceProcs[i]
+    //        << " face:" << patchFaces[i]
+    //        << " at:" << patchFc[i] << endl
+    //        << "Found data in region " << sampleRegion_
+    //        << " at proc:" << sampleProcs[i]
+    //        << " face:" << sampleIndices[i]
+    //        << " at:" << sampleLocations[i]
+    //        << nl << endl;
+    //}
+
+
+
+    if (debug && Pstream::master())
+    {
+        OFstream str
+        (
+            patch_.boundaryMesh().mesh().time().path()
+          / patch_.name()
+          + "_directMapped.obj"
+        );
+        Pout<< "Dumping mapping as lines from patch faceCentres to"
+            << " sampled cellCentres to file " << str.name() << endl;
+
+        label vertI = 0;
+
+        forAll(patchFc, i)
+        {
+            meshTools::writeOBJ(str, patchFc[i]);
+            vertI++;
+            meshTools::writeOBJ(str, sampleLocations[i]);
+            vertI++;
+            str << "l " << vertI-1 << ' ' << vertI << nl;
+        }
+    }
+
+
+    // Check that actual offset vector (sampleLocations - patchFc) is more or
+    // less constant.
+    if (Pstream::master())
+    {
+        const scalarField magOffset(mag(sampleLocations - patchFc));
+        const scalar avgOffset(average(magOffset));
+
+        forAll(magOffset, sampleI)
+        {
+            if (mag(magOffset[sampleI]-avgOffset) > max(SMALL, 0.001*avgOffset))
+            {
+                WarningIn("directMappedPatchBase::calcMapping() const")
+                    << "The actual cell/face centres picked up using offset "
+                    << offset_ << " are not" << endl
+                    << "    on a single plane."
+                    << " This might give numerical problems." << endl
+                    << "    At patchface " << patchFc[sampleI]
+                    << " the sampled cell/face " << sampleLocations[sampleI]
+                    << endl
+                    << "    is not on a plane " << avgOffset
+                    << " offset from the patch." << endl
+                    << "    You might want to shift your plane offset."
+                    << " Set the debug flag to get a dump of sampled cells."
+                    << endl;
+                break;
+            }
+        }
+    }
+
+
+    // Determine schedule.
+    mapPtr_.reset(new mapDistribute(sampleProcs, patchFaceProcs));
+
+    // Rework the schedule from indices into samples to cell data to send,
+    // face data to receive.
+
+    labelListList& subMap = mapPtr_().subMap();
+    labelListList& constructMap = mapPtr_().constructMap();
+
+    forAll(subMap, procI)
+    {
+        subMap[procI] = UIndirectList<label>
+        (
+            sampleIndices,
+            subMap[procI]
+        );
+        constructMap[procI] = UIndirectList<label>
+        (
+            patchFaces,
+            constructMap[procI]
+        );
+
+        if (debug)
+        {
+            Pout<< "To proc:" << procI << " sending values of cells/faces:"
+                << subMap[procI] << endl;
+            Pout<< "From proc:" << procI << " receiving values of patch faces:"
+                << constructMap[procI] << endl;
+        }
+    }
+
+    // Redo constructSize
+    mapPtr_().constructSize() = patch_.size();
+
+    if (debug)
+    {
+        // Check that all elements get a value.
+        PackedBoolList used(patch_.size());
+        forAll(constructMap, procI)
+        {
+            const labelList& map = constructMap[procI];
+
+            forAll(map, i)
+            {
+                label faceI = map[i];
+
+                if (used[faceI] == 0)
+                {
+                    used[faceI] = 1;
+                }
+                else
+                {
+                    FatalErrorIn("directMappedPatchBase::calcMapping() const")
+                        << "On patch " << patch_.name()
+                        << " patchface " << faceI
+                        << " is assigned to more than once."
+                        << abort(FatalError);
+                }
+            }
+        }
+        forAll(used, faceI)
+        {
+            if (used[faceI] == 0)
+            {
+                FatalErrorIn("directMappedPatchBase::calcMapping() const")
+                    << "On patch " << patch_.name()
+                    << " patchface " << faceI
+                    << " is never assigned to."
+                    << abort(FatalError);
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * * * * * //
+
+Foam::directMappedPatchBase::directMappedPatchBase
+(
+     const polyPatch& pp
+)
+:
+    patch_(pp),
+    sampleRegion_(patch_.boundaryMesh().mesh().name()),
+    mode_(NEARESTPATCHFACE),
+    samplePatch_("none"),
+    offset_(vector::zero),
+    sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
+    mapPtr_(NULL)
+{}
+
+
+Foam::directMappedPatchBase::directMappedPatchBase
+(
+    const polyPatch& pp,
+    const dictionary& dict
+)
+:
+    patch_(pp),
+    sampleRegion_
+    (
+        dict.lookupOrDefault
+        (
+            "sampleRegion",
+            patch_.boundaryMesh().mesh().name()
+        )
+    ),
+    mode_(sampleModeNames_.read(dict.lookup("sampleMode"))),
+    samplePatch_(dict.lookup("samplePatch")),
+    offset_(dict.lookup("offset")),
+    sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
+    mapPtr_(NULL)
+{}
+
+
+Foam::directMappedPatchBase::directMappedPatchBase
+(
+    const polyPatch& pp,
+    const directMappedPatchBase& dmp
+)
+:
+    patch_(pp),
+    sampleRegion_(dmp.sampleRegion_),
+    mode_(dmp.mode_),
+    samplePatch_(dmp.samplePatch_),
+    offset_(dmp.offset_),
+    sameRegion_(dmp.sameRegion_),
+    mapPtr_(NULL)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::directMappedPatchBase::~directMappedPatchBase()
+{
+    clearOut();
+}
+
+
+void Foam::directMappedPatchBase::clearOut()
+{
+    mapPtr_.clear();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::polyMesh& Foam::directMappedPatchBase::sampleMesh() const
+{
+    return patch_.boundaryMesh().mesh().time().lookupObject<polyMesh>
+    (
+        sampleRegion_
+    );
+}
+
+
+const Foam::polyPatch& Foam::directMappedPatchBase::samplePolyPatch() const
+{
+    const polyMesh& nbrMesh = sampleMesh();
+
+    label patchI = nbrMesh.boundaryMesh().findPatchID(samplePatch_);
+
+    if (patchI == -1)
+    {
+        FatalErrorIn("directMappedPatchBase::samplePolyPatch() ")
+            << "Cannot find patch " << samplePatch_
+            << " in region " << sampleRegion_ << endl
+            << "Valid patches are " << nbrMesh.boundaryMesh().names()
+            << exit(FatalError);
+    }
+
+    return nbrMesh.boundaryMesh()[patchI];
+}
+
+
+void Foam::directMappedPatchBase::write(Ostream& os) const
+{
+    os.writeKeyword("sampleMode") << sampleModeNames_[mode_]
+        << token::END_STATEMENT << nl;
+    os.writeKeyword("sampleRegion") << sampleRegion_
+        << token::END_STATEMENT << nl;
+    os.writeKeyword("samplePatch") << samplePatch_
+        << token::END_STATEMENT << nl;
+    os.writeKeyword("offset") << offset_ << token::END_STATEMENT << nl;
+}
+
+
+// ************************************************************************* //

src/meshTools/directMapped/directMappedPolyPatch/directMappedPatchBase.H

@@ -0,0 +1,248 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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
+    Foam::directMappedPatchBase
+
+Description
+    Determines a mapping between patch face centres and mesh cell or face
+    centres and processors they're on.
+
+Note
+    Storage is not optimal. It temporary collects all (patch)face centres
+    on all processors to keep the addressing calculation simple.
+
+SourceFiles
+    directMappedPatchBase.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef directMappedPatchBase_H
+#define directMappedPatchBase_H
+
+#include "pointField.H"
+#include "Tuple2.H"
+#include "pointIndexHit.H"
+#include "mapDistribute.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class polyPatch;
+class polyMesh;
+
+/*---------------------------------------------------------------------------*\
+                      Class directMappedPatchBase Declaration
+\*---------------------------------------------------------------------------*/
+
+class directMappedPatchBase
+{
+
+public:
+
+        //- Mesh items to sample
+        enum sampleMode
+        {
+            NEARESTCELL,
+            NEARESTPATCHFACE,
+            NEARESTFACE
+        };
+
+private:
+
+    // Private data
+
+        static const NamedEnum<sampleMode, 3> sampleModeNames_;
+
+        //- Patch to sample
+        const polyPatch& patch_;
+
+        //- Region to sample
+        const word sampleRegion_;
+
+        //- What to sample
+        const sampleMode mode_;
+
+        //- Patch (only if NEARESTBOUNDARY)
+        const word samplePatch_;
+
+        //- Offset vector
+        const vector offset_;
+
+        //- Same region
+        const bool sameRegion_;
+
+
+        // Derived information
+
+            //- Communication schedule:
+            //  - Cells/faces to sample per processor
+            //  - Patch faces to receive per processor
+            //  - schedule
+            mutable autoPtr<mapDistribute> mapPtr_;
+
+    // Private Member Functions
+
+        //- Collect single list of samples and originating processor+face.
+        void collectSamples
+        (
+            pointField&,
+            labelList& patchFaceProcs,
+            labelList& patchFaces,
+            pointField& patchFc
+        ) const;
+
+        //- Find cells/faces containing samples
+        void findSamples
+        (
+            const pointField&,
+            labelList& sampleProcs,     // processor containing sample
+            labelList& sampleIndices,   // local index of cell/face
+            pointField& sampleLocations // actual representative location
+        ) const;
+
+        //- Calculate matching
+        void calcMapping() const;
+
+
+    // Private class
+
+        //- Private class for finding nearest
+        //  - point+local index
+        //  - sqr(distance)
+        //  - processor
+        typedef Tuple2<pointIndexHit, Tuple2<scalar, label> > nearInfo;
+
+        class nearestEqOp
+        {
+
+        public:
+
+            void operator()(nearInfo& x, const nearInfo& y) const
+            {
+                if (y.first().hit())
+                {
+                    if (!x.first().hit())
+                    {
+                        x = y;
+                    }
+                    else if (y.second().first() < x.second().first())
+                    {
+                        x = y;
+                    }
+                }
+            }
+        };
+
+
+public:
+
+    //- Runtime type information
+    ClassName("directMappedPatchBase");
+
+
+    // Constructors
+
+        //- Construct from components
+        directMappedPatchBase(const polyPatch&);
+
+        //- Construct from dictionary
+        directMappedPatchBase(const polyPatch&, const dictionary&);
+
+        //- Construct as copy, resetting patch
+        directMappedPatchBase(const polyPatch&, const directMappedPatchBase&);
+
+
+    // Destructor
+
+        ~directMappedPatchBase();
+
+        void clearOut();
+
+
+    // Member functions
+
+        //- What to sample
+        const sampleMode& mode() const
+        {
+            return mode_;
+        }
+
+        //- Region to sample
+        const word& sampleRegion() const
+        {
+            return sampleRegion_;
+        }
+
+        //- Patch (only if NEARESTBOUNDARY)
+        const word& samplePatch() const
+        {
+            return samplePatch_;
+        }
+
+        //- Offset vector (from patch faces to destination mesh objects)
+        const vector& offset() const
+        {
+            return offset_;
+        }
+
+        //- Return reference to the parallel distribution map
+        const mapDistribute& map() const
+        {
+            if (mapPtr_.empty())
+            {
+                calcMapping();
+            }
+            return mapPtr_();
+        }
+
+        //- Cached sampleRegion != mesh.name()
+        bool sameRegion() const
+        {
+            return sameRegion_;
+        }
+
+        //- Get the region mesh
+        const polyMesh& sampleMesh() const;
+
+        //- Get the patch on the region
+        const polyPatch& samplePolyPatch() const;
+
+        //- Write as a dictionary
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.C

@@ -26,11 +26,6 @@ License
 
 #include "directMappedPolyPatch.H"
 #include "addToRunTimeSelectionTable.H"
-#include "ListListOps.H"
-#include "meshSearch.H"
-#include "mapDistribute.H"
-#include "meshTools.H"
-#include "OFstream.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -40,469 +35,11 @@ namespace Foam
 
     addToRunTimeSelectionTable(polyPatch, directMappedPolyPatch, word);
     addToRunTimeSelectionTable(polyPatch, directMappedPolyPatch, dictionary);
-
-
-    template<>
-    const char* NamedEnum<directMappedPolyPatch::sampleMode, 3>::names[] =
-    {
-        "nearestCell",
-        "nearestPatchFace",
-        "nearestFace"
-    };
-
-    const NamedEnum<directMappedPolyPatch::sampleMode, 3>
-        directMappedPolyPatch::sampleModeNames_;
 }
 
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::directMappedPolyPatch::collectSamples
-(
-    pointField& samples,
-    labelList& patchFaceProcs,
-    labelList& patchFaces,
-    pointField& patchFc
-) const
-{
-
-    // Collect all sample points and the faces they come from.
-    List<pointField> globalFc(Pstream::nProcs());
-    List<pointField> globalSamples(Pstream::nProcs());
-    labelListList globalFaces(Pstream::nProcs());
-
-    globalFc[Pstream::myProcNo()] = this->faceCentres();
-    globalSamples[Pstream::myProcNo()] = globalFc[Pstream::myProcNo()]+offset_;
-    globalFaces[Pstream::myProcNo()] = identity(size());
-
-    // Distribute to all processors
-    Pstream::gatherList(globalSamples);
-    Pstream::scatterList(globalSamples);
-    Pstream::gatherList(globalFaces);
-    Pstream::scatterList(globalFaces);
-    Pstream::gatherList(globalFc);
-    Pstream::scatterList(globalFc);
-
-    // Rework into straight list
-    samples = ListListOps::combine<pointField>
-    (
-        globalSamples,
-        accessOp<pointField>()
-    );
-    patchFaces = ListListOps::combine<labelList>
-    (
-        globalFaces,
-        accessOp<labelList>()
-    );
-    patchFc = ListListOps::combine<pointField>
-    (
-        globalFc,
-        accessOp<pointField>()
-    );
-
-    patchFaceProcs.setSize(patchFaces.size());
-    labelList nPerProc
-    (
-        ListListOps::subSizes
-        (
-            globalFaces,
-            accessOp<labelList>()
-        )
-    );
-    label sampleI = 0;
-    forAll(nPerProc, procI)
-    {
-        for (label i = 0; i < nPerProc[procI]; i++)
-        {
-            patchFaceProcs[sampleI++] = procI;
-        }
-    }
-}
-
-
-// Find the processor/cell containing the samples. Does not account
-// for samples being found in two processors.
-void Foam::directMappedPolyPatch::findSamples
-(
-    const pointField& samples,
-    labelList& sampleProcs,
-    labelList& sampleIndices,
-    pointField& sampleLocations
-) const
-{
-    const polyMesh& mesh = boundaryMesh().mesh();
-
-    // All the info for nearest. Construct to miss
-    List<nearInfo> nearest(samples.size());
-
-    switch (mode_)
-    {
-        case NEARESTCELL:
-        {
-            if (samplePatch_.size() && samplePatch_ != "none")
-            {
-                FatalErrorIn
-                (
-                    "directMappedPolyPatch::findSamples(const pointField&,"
-                    " labelList&, labelList&, pointField&) const"
-                )   << "No need to supply a patch name when in "
-                    << sampleModeNames_[mode_] << " mode." << exit(FatalError);
-            }
-
-            // Octree based search engine
-            meshSearch meshSearchEngine(mesh, false);
-
-            forAll(samples, sampleI)
-            {
-                const point& sample = samples[sampleI];
-
-                label cellI = meshSearchEngine.findCell(sample);
-
-                if (cellI == -1)
-                {
-                    nearest[sampleI].second().first() = Foam::sqr(GREAT);
-                    nearest[sampleI].second().second() = Pstream::myProcNo();
-                }
-                else
-                {
-                    const point& cc = mesh.cellCentres()[cellI];
-
-                    nearest[sampleI].first() = pointIndexHit
-                    (
-                        true,
-                        cc,
-                        cellI
-                    );
-                    nearest[sampleI].second().first() = magSqr(cc-sample);
-                    nearest[sampleI].second().second() = Pstream::myProcNo();
-                }
-            }
-            break;
-        }
-
-        case NEARESTPATCHFACE:
-        {
-            label patchI = boundaryMesh().findPatchID(samplePatch_);
-
-            if (patchI == -1)
-            {
-                FatalErrorIn("directMappedPolyPatch::findSamples(..)")
-                    << "Cannot find patch " << samplePatch_ << endl
-                    << "Valid patches are " << boundaryMesh().names()
-                    << exit(FatalError);
-            }
-
-            Random rndGen(123456);
-
-            const polyPatch& pp = boundaryMesh()[patchI];
-
-            if (pp.empty())
-            {
-                forAll(samples, sampleI)
-                {
-                    nearest[sampleI].second().first() = Foam::sqr(GREAT);
-                    nearest[sampleI].second().second() = Pstream::myProcNo();
-                }
-            }
-            else
-            {
-                // patch faces
-                const labelList patchFaces(identity(pp.size()) + pp.start());
-
-                treeBoundBox patchBb
-                (
-                    treeBoundBox(pp.points(), pp.meshPoints()).extend
-                    (
-                        rndGen,
-                        1E-4
-                    )
-                );
-                patchBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
-                patchBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
-
-                autoPtr<indexedOctree<treeDataFace> > boundaryTree
-                (
-                    new indexedOctree<treeDataFace>
-                    (
-                        treeDataFace  // all information needed to search faces
-                        (
-                            false,                      // do not cache bb
-                            mesh,
-                            patchFaces                  // boundary faces only
-                        ),
-                        patchBb,   // overall search domain
-                        8,                              // maxLevel
-                        10,                             // leafsize
-                        3.0                             // duplicity
-                    )
-                );
-
-                forAll(samples, sampleI)
-                {
-                    const point& sample = samples[sampleI];
-
-                    pointIndexHit& nearInfo = nearest[sampleI].first();
-                    nearInfo = boundaryTree().findNearest
-                    (
-                        sample,
-                        magSqr(patchBb.span())
-                    );
-
-                    if (!nearInfo.hit())
-                    {
-                        nearest[sampleI].second().first() = Foam::sqr(GREAT);
-                        nearest[sampleI].second().second() =
-                            Pstream::myProcNo();
-                    }
-                    else
-                    {
-                        point fc(pp[nearInfo.index()].centre(pp.points()));
-                        nearInfo.setPoint(fc);
-                        nearest[sampleI].second().first() = magSqr(fc-sample);
-                        nearest[sampleI].second().second() =
-                            Pstream::myProcNo();
-                    }
-                }
-            }
-            break;
-        }
-
-        case NEARESTFACE:
-        {
-            if (samplePatch_.size() && samplePatch_ != "none")
-            {
-                FatalErrorIn
-                (
-                    "directMappedPolyPatch::findSamples(const pointField&,"
-                    " labelList&, labelList&, pointField&) const"
-                )   << "No need to supply a patch name when in "
-                    << sampleModeNames_[mode_] << " mode." << exit(FatalError);
-            }
-
-            // Octree based search engine
-            meshSearch meshSearchEngine(mesh, false);
-
-            forAll(samples, sampleI)
-            {
-                const point& sample = samples[sampleI];
-
-                label faceI = meshSearchEngine.findNearestFace(sample);
-
-                if (faceI == -1)
-                {
-                    nearest[sampleI].second().first() = Foam::sqr(GREAT);
-                    nearest[sampleI].second().second() = Pstream::myProcNo();
-                }
-                else
-                {
-                    const point& fc = mesh.faceCentres()[faceI];
-
-                    nearest[sampleI].first() = pointIndexHit
-                    (
-                        true,
-                        fc,
-                        faceI
-                    );
-                    nearest[sampleI].second().first() = magSqr(fc-sample);
-                    nearest[sampleI].second().second() = Pstream::myProcNo();
-                }
-            }
-            break;
-        }
-
-        default:
-        {
-            FatalErrorIn("directMappedPolyPatch::findSamples(..)")
-                << "problem." << abort(FatalError);
-        }
-    }
-
-
-    // Find nearest.
-    Pstream::listCombineGather(nearest, nearestEqOp());
-    Pstream::listCombineScatter(nearest);
-
-    if (debug)
-    {
-        Info<< "directMappedPolyPatch::findSamples : " << endl;
-        forAll(nearest, sampleI)
-        {
-            label procI = nearest[sampleI].second().second();
-            label localI = nearest[sampleI].first().index();
-
-            Info<< "    " << sampleI << " coord:"<< samples[sampleI]
-                << " found on processor:" << procI
-                << " in local cell/face:" << localI
-                << " with cc:" << nearest[sampleI].first().rawPoint() << endl;
-        }
-    }
-
-    // Check for samples not being found
-    forAll(nearest, sampleI)
-    {
-        if (!nearest[sampleI].first().hit())
-        {
-            FatalErrorIn
-            (
-                "directMappedPolyPatch::findSamples"
-                "(const pointField&, labelList&"
-                ", labelList&, pointField&)"
-            )   << "Did not find sample " << samples[sampleI]
-                << " on any processor." << exit(FatalError);
-        }
-    }
-
-
-    // Convert back into proc+local index
-    sampleProcs.setSize(samples.size());
-    sampleIndices.setSize(samples.size());
-    sampleLocations.setSize(samples.size());
-
-    forAll(nearest, sampleI)
-    {
-        sampleProcs[sampleI] = nearest[sampleI].second().second();
-        sampleIndices[sampleI] = nearest[sampleI].first().index();
-        sampleLocations[sampleI] = nearest[sampleI].first().hitPoint();
-    }
-}
-
-
-void Foam::directMappedPolyPatch::calcMapping() const
-{
-    if (sendLabelsPtr_.valid())
-    {
-        FatalErrorIn("directMappedPolyPatch::calcMapping() const")
-            << "Mapping already calculated" << exit(FatalError);
-    }
-
-    if (offset_ == vector::zero)
-    {
-        FatalErrorIn("directMappedPolyPatch::calcMapping() const")
-            << "Invalid offset " << offset_ << endl
-            << "Offset is the vector added to the patch face centres to"
-            << " find the cell supplying the data."
-            << exit(FatalError);
-    }
-
-
-    // Get global list of all samples and the processor and face they come from.
-    pointField samples;
-    labelList patchFaceProcs;
-    labelList patchFaces;
-    pointField patchFc;
-    collectSamples(samples, patchFaceProcs, patchFaces, patchFc);
-
-    // Find processor and cell/face samples are in and actual location.
-    labelList sampleProcs;
-    labelList sampleIndices;
-    pointField sampleLocations;
-    findSamples(samples, sampleProcs, sampleIndices, sampleLocations);
-
-
-    // Now we have all the data we need:
-    // - where sample originates from (so destination when mapping):
-    //   patchFaces, patchFaceProcs.
-    // - cell/face sample is in (so source when mapping)
-    //   sampleIndices, sampleProcs.
-
-
-    if (debug && Pstream::master())
-    {
-        OFstream str
-        (
-            boundaryMesh().mesh().time().path()
-          / name()
-          + "_directMapped.obj"
-        );
-        Pout<< "Dumping mapping as lines from patch faceCentres to"
-            << " sampled cellCentres to file " << str.name() << endl;
-
-        label vertI = 0;
-
-        forAll(patchFc, i)
-        {
-            meshTools::writeOBJ(str, patchFc[i]);
-            vertI++;
-            meshTools::writeOBJ(str, sampleLocations[i]);
-            vertI++;
-            str << "l " << vertI-1 << ' ' << vertI << nl;
-        }
-    }
-
-
-    // Check that actual offset vector (sampleLocations - patchFc) is more or
-    // less constant.
-    if (Pstream::master())
-    {
-        const scalarField magOffset(mag(sampleLocations - patchFc));
-        const scalar avgOffset(average(magOffset));
-
-        forAll(magOffset, sampleI)
-        {
-            if (mag(magOffset[sampleI]-avgOffset) > 0.001*avgOffset)
-            {
-                WarningIn("directMappedPolyPatch::calcMapping() const")
-                    << "The actual cell centres picked up using offset "
-                    << offset_ << " are not" << endl
-                    << "    on a single plane."
-                    << " This might give numerical problems." << endl
-                    << "    At patchface " << patchFc[sampleI]
-                    << " the sampled cell " << sampleLocations[sampleI] << endl
-                    << "    is not on a plane " << avgOffset
-                    << " offset from the patch." << endl
-                    << "    You might want to shift your plane offset."
-                    << " Set the debug flag to get a dump of sampled cells."
-                    << endl;
-                break;
-            }
-        }
-    }
-
-
-    // Determine schedule.
-    mapDistribute distMap(sampleProcs, patchFaceProcs);
-
-    // Rework the schedule to cell data to send, face data to receive.
-    schedulePtr_.reset(new List<labelPair>(distMap.schedule()));
-
-    const labelListList& subMap = distMap.subMap();
-    const labelListList& constructMap = distMap.constructMap();
-
-    // Extract the particular data I need to send and receive.
-    sendLabelsPtr_.reset(new labelListList(subMap.size()));
-    labelListList& sendLabels = sendLabelsPtr_();
-
-    forAll(subMap, procI)
-    {
-        sendLabels[procI] = UIndirectList<label>
-        (
-            sampleIndices,
-            subMap[procI]
-        );
-
-        if (debug)
-        {
-            Pout<< "To proc:" << procI << " sending values of cells/faces:"
-                << sendLabels[procI] << endl;
-        }
-    }
-
-    receiveFaceLabelsPtr_.reset(new labelListList(constructMap.size()));
-    labelListList& receiveFaceLabels = receiveFaceLabelsPtr_();
-
-    forAll(constructMap, procI)
-    {
-        receiveFaceLabels[procI] =
-            UIndirectList<label>(patchFaces, constructMap[procI]);
-
-        if (debug)
-        {
-            Pout<< "From proc:" << procI << " receiving values of patch faces:"
-                << receiveFaceLabels[procI] << endl;
-        }
-    }
-}
-
 
 // * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * * * * * //
 
@@ -516,12 +53,7 @@ Foam::directMappedPolyPatch::directMappedPolyPatch
 )
 :
     polyPatch(name, size, start, index, bm),
-    mode_(NEARESTPATCHFACE),
-    samplePatch_("none"),
-    offset_(vector::zero),
-    schedulePtr_(NULL),
-    sendLabelsPtr_(NULL),
-    receiveFaceLabelsPtr_(NULL)
+    directMappedPatchBase(static_cast<const polyPatch&>(*this))
 {}
 
 
@@ -534,12 +66,7 @@ Foam::directMappedPolyPatch::directMappedPolyPatch
 )
 :
     polyPatch(name, dict, index, bm),
-    mode_(sampleModeNames_.read(dict.lookup("sampleMode"))),
-    samplePatch_(dict.lookup("samplePatch")),
-    offset_(dict.lookup("offset")),
-    schedulePtr_(NULL),
-    sendLabelsPtr_(NULL),
-    receiveFaceLabelsPtr_(NULL)
+    directMappedPatchBase(*this, dict)
 {}
 
 
@@ -550,12 +77,7 @@ Foam::directMappedPolyPatch::directMappedPolyPatch
 )
 :
     polyPatch(pp, bm),
-    mode_(pp.mode_),
-    samplePatch_(pp.samplePatch_),
-    offset_(pp.offset_),
-    schedulePtr_(NULL),
-    sendLabelsPtr_(NULL),
-    receiveFaceLabelsPtr_(NULL)
+    directMappedPatchBase(*this, pp)
 {}
 
 
@@ -569,12 +91,7 @@ Foam::directMappedPolyPatch::directMappedPolyPatch
 )
 :
     polyPatch(pp, bm, index, newSize, newStart),
-    mode_(pp.mode_),
-    samplePatch_(pp.samplePatch_),
-    offset_(pp.offset_),
-    schedulePtr_(NULL),
-    sendLabelsPtr_(NULL),
-    receiveFaceLabelsPtr_(NULL)
+    directMappedPatchBase(*this, pp)
 {}
 
 
@@ -582,28 +99,59 @@ Foam::directMappedPolyPatch::directMappedPolyPatch
 
 Foam::directMappedPolyPatch::~directMappedPolyPatch()
 {
-    clearOut();
+    directMappedPatchBase::clearOut();
 }
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::directMappedPolyPatch::clearOut()
+//- Initialise the calculation of the patch geometry
+void Foam::directMappedPolyPatch::initGeometry()
 {
-    schedulePtr_.clear();
-    sendLabelsPtr_.clear();
-    receiveFaceLabelsPtr_.clear();
+    polyPatch::initGeometry();
+    directMappedPatchBase::clearOut();
+}
+
+//- Calculate the patch geometry
+void Foam::directMappedPolyPatch::calcGeometry()
+{
+    polyPatch::calcGeometry();
+    directMappedPatchBase::clearOut();
+}
+
+//- Initialise the patches for moving points
+void Foam::directMappedPolyPatch::initMovePoints(const pointField& p)
+{
+    polyPatch::initMovePoints(p);
+    directMappedPatchBase::clearOut();
+}
+
+//- Correct patches after moving points
+void Foam::directMappedPolyPatch::movePoints(const pointField& p)
+{
+    polyPatch::movePoints(p);
+    directMappedPatchBase::clearOut();
+}
+
+//- Initialise the update of the patch topology
+void Foam::directMappedPolyPatch::initUpdateMesh()
+{
+    polyPatch::initUpdateMesh();
+    directMappedPatchBase::clearOut();
+}
+
+//- Update of the patch topology
+void Foam::directMappedPolyPatch::updateMesh()
+{
+    polyPatch::updateMesh();
+    directMappedPatchBase::clearOut();
 }
 
 
 void Foam::directMappedPolyPatch::write(Ostream& os) const
 {
     polyPatch::write(os);
-    os.writeKeyword("sampleMode") << sampleModeNames_[mode_]
-        << token::END_STATEMENT << nl;
-    os.writeKeyword("samplePatch") << samplePatch_
-        << token::END_STATEMENT << nl;
-    os.writeKeyword("offset") << offset_ << token::END_STATEMENT << nl;
+    directMappedPatchBase::write(os);
 }
 
 

src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.H

@@ -42,9 +42,7 @@ SourceFiles
 #define directMappedPolyPatch_H
 
 #include "polyPatch.H"
-#include "labelPair.H"
-#include "Tuple2.H"
-#include "pointIndexHit.H"
+#include "directMappedPatchBase.H"
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -52,148 +50,37 @@ SourceFiles
 namespace Foam
 {
 
+class polyMesh;
+
 /*---------------------------------------------------------------------------*\
                       Class directMappedPolyPatch Declaration
 \*---------------------------------------------------------------------------*/
 
 class directMappedPolyPatch
 :
-    public polyPatch
+    public polyPatch,
+    public directMappedPatchBase
 {
 
-public:
-
-        //- Mesh items to sample
-        enum sampleMode
-        {
-            NEARESTCELL,
-            NEARESTPATCHFACE,
-            NEARESTFACE
-        };
-
-private:
-
-    // Private data
-
-        static const NamedEnum<sampleMode, 3> sampleModeNames_;
-
-        //- What to sample
-        sampleMode mode_;
-
-        //- Patch (only if NEARESTBOUNDARY)
-        const word samplePatch_;
-
-        //- Offset vector
-        const vector offset_;
-
-
-
-        // Derived information
-
-            //- Communication schedule.
-            mutable autoPtr<List<labelPair> > schedulePtr_;
-
-            //- Cells/faces to sample per processor
-            mutable autoPtr<labelListList> sendLabelsPtr_;
-
-            //- Patch faces to receive per processor
-            mutable autoPtr<labelListList> receiveFaceLabelsPtr_;
-
-
-    // Private Member Functions
-
-        //- Collect single list of samples and originating processor+face.
-        void collectSamples
-        (
-            pointField&,
-            labelList& patchFaceProcs,
-            labelList& patchFaces,
-            pointField& patchFc
-        ) const;
-
-        //- Find cells/faces containing samples
-        void findSamples
-        (
-            const pointField&,
-            labelList& sampleProcs,     // processor containing sample
-            labelList& sampleIndices,   // local index of cell/face
-            pointField& sampleLocations // actual representative location
-        ) const;
-
-        //- Calculate matching
-        void calcMapping() const;
-
-
-    // Private class
-
-        //- Private class for finding nearest
-        //  - point+local index
-        //  - sqr(distance)
-        //  - processor
-        typedef Tuple2<pointIndexHit, Tuple2<scalar, label> > nearInfo;
-
-        class nearestEqOp
-        {
-
-        public:
-
-            void operator()(nearInfo& x, const nearInfo& y) const
-            {
-                if (y.first().hit())
-                {
-                    if (!x.first().hit())
-                    {
-                        x = y;
-                    }
-                    else if (y.second().first() < x.second().first())
-                    {
-                        x = y;
-                    }
-                }
-            }
-        };
-
-
-
 protected:
 
-        void clearOut();
-
         //- Initialise the calculation of the patch geometry
-        virtual void initGeometry()
-        {
-            clearOut();
-        }
+        virtual void initGeometry();
 
         //- Calculate the patch geometry
-        virtual void calcGeometry()
-        {
-            clearOut();
-        }
+        virtual void calcGeometry();
 
         //- Initialise the patches for moving points
-        virtual void initMovePoints(const pointField&)
-        {
-            clearOut();
-        }
+        virtual void initMovePoints(const pointField&);
 
         //- Correct patches after moving points
-        virtual void movePoints(const pointField& p)
-        {
-            clearOut();
-        }
+        virtual void movePoints(const pointField&);
 
         //- Initialise the update of the patch topology
-        virtual void initUpdateMesh()
-        {
-            clearOut();
-        }
+        virtual void initUpdateMesh();
 
         //- Update of the patch topology
-        virtual void updateMesh()
-        {
-            clearOut();
-        }
+        virtual void updateMesh();
 
 
 public:
@@ -271,72 +158,6 @@ public:
 
     // Member functions
 
-        //- What to sample
-        const sampleMode& mode() const
-        {
-            return mode_;
-        }
-
-        //- Patch (only if NEARESTBOUNDARY)
-        const word& samplePatch() const
-        {
-            return samplePatch_;
-        }
-
-        //- Offset vector (from patch faces to destination mesh objects)
-        const vector& offset() const
-        {
-            return offset_;
-        }
-
-        //- Access to communication.
-        const List<labelPair>& schedule() const
-        {
-            if (schedulePtr_.empty())
-            {
-                calcMapping();
-            }
-            return schedulePtr_();
-        }
-
-        //- Cells/faces to sample per processor
-        const labelListList& sendLabels() const
-        {
-            if (debug)
-            {
-                Pout<< "Asking for sendLabels." << endl;
-            }
-
-            if (sendLabelsPtr_.empty())
-            {
-                if (debug)
-                {
-                    Pout<< "Calculating mapping." << endl;
-                    calcMapping();
-                }
-            }
-            return sendLabelsPtr_();
-        }
-
-        //- Patch faces to receive per processor
-        const labelListList& receiveFaceLabels() const
-        {
-            if (debug)
-            {
-                Pout<< "Asking for receiveFaceLabels." << endl;
-            }
-
-            if (receiveFaceLabelsPtr_.empty())
-            {
-                if (debug)
-                {
-                    Pout<< "Calculating mapping." << endl;
-                    calcMapping();
-                }
-            }
-            return receiveFaceLabelsPtr_();
-        }
-
         //- Write the polyPatch data as a dictionary
         virtual void write(Ostream&) const;
 };

src/meshTools/directMapped/directMappedPolyPatch/directMappedWallPolyPatch.C

@@ -0,0 +1,163 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 "directMappedWallPolyPatch.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(directMappedWallPolyPatch, 0);
+
+    addToRunTimeSelectionTable(polyPatch, directMappedWallPolyPatch, word);
+    addToRunTimeSelectionTable
+    (
+        polyPatch,
+        directMappedWallPolyPatch,
+        dictionary
+    );
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * * * * * //
+
+Foam::directMappedWallPolyPatch::directMappedWallPolyPatch
+(
+    const word& name,
+    const label size,
+    const label start,
+    const label index,
+    const polyBoundaryMesh& bm
+)
+:
+    wallPolyPatch(name, size, start, index, bm),
+    directMappedPatchBase(static_cast<const polyPatch&>(*this))
+{}
+
+
+Foam::directMappedWallPolyPatch::directMappedWallPolyPatch
+(
+    const word& name,
+    const dictionary& dict,
+    const label index,
+    const polyBoundaryMesh& bm
+)
+:
+    wallPolyPatch(name, dict, index, bm),
+    directMappedPatchBase(*this, dict)
+{}
+
+
+Foam::directMappedWallPolyPatch::directMappedWallPolyPatch
+(
+    const directMappedWallPolyPatch& pp,
+    const polyBoundaryMesh& bm
+)
+:
+    wallPolyPatch(pp, bm),
+    directMappedPatchBase(*this, pp)
+{}
+
+
+Foam::directMappedWallPolyPatch::directMappedWallPolyPatch
+(
+    const directMappedWallPolyPatch& pp,
+    const polyBoundaryMesh& bm,
+    const label index,
+    const label newSize,
+    const label newStart
+)
+:
+    wallPolyPatch(pp, bm, index, newSize, newStart),
+    directMappedPatchBase(*this, pp)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::directMappedWallPolyPatch::~directMappedWallPolyPatch()
+{
+    directMappedPatchBase::clearOut();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+//- Initialise the calculation of the patch geometry
+void Foam::directMappedWallPolyPatch::initGeometry()
+{
+    wallPolyPatch::initGeometry();
+    directMappedPatchBase::clearOut();
+}
+
+//- Calculate the patch geometry
+void Foam::directMappedWallPolyPatch::calcGeometry()
+{
+    wallPolyPatch::calcGeometry();
+    directMappedPatchBase::clearOut();
+}
+
+//- Initialise the patches for moving points
+void Foam::directMappedWallPolyPatch::initMovePoints(const pointField& p)
+{
+    wallPolyPatch::initMovePoints(p);
+    directMappedPatchBase::clearOut();
+}
+
+//- Correct patches after moving points
+void Foam::directMappedWallPolyPatch::movePoints(const pointField& p)
+{
+    wallPolyPatch::movePoints(p);
+    directMappedPatchBase::clearOut();
+}
+
+//- Initialise the update of the patch topology
+void Foam::directMappedWallPolyPatch::initUpdateMesh()
+{
+    wallPolyPatch::initUpdateMesh();
+    directMappedPatchBase::clearOut();
+}
+
+//- Update of the patch topology
+void Foam::directMappedWallPolyPatch::updateMesh()
+{
+    wallPolyPatch::updateMesh();
+    directMappedPatchBase::clearOut();
+}
+
+
+void Foam::directMappedWallPolyPatch::write(Ostream& os) const
+{
+    wallPolyPatch::write(os);
+    directMappedPatchBase::write(os);
+}
+
+
+// ************************************************************************* //

src/meshTools/directMapped/directMappedPolyPatch/directMappedWallPolyPatch.H

@@ -0,0 +1,181 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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
+    Foam::directMappedWallPolyPatch
+
+Description
+    Determines a mapping between patch face centres and mesh cell or face
+    centres and processors they're on.
+
+Note
+    Storage is not optimal. It stores all face centres and cells on all
+    processors to keep the addressing calculation simple.
+
+SourceFiles
+    directMappedWallPolyPatch.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef directMappedWallPolyPatch_H
+#define directMappedWallPolyPatch_H
+
+#include "wallPolyPatch.H"
+#include "directMappedPatchBase.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class polyMesh;
+
+/*---------------------------------------------------------------------------*\
+                      Class directMappedWallPolyPatch Declaration
+\*---------------------------------------------------------------------------*/
+
+class directMappedWallPolyPatch
+:
+    public wallPolyPatch,
+    public directMappedPatchBase
+{
+
+protected:
+
+        //- Initialise the calculation of the patch geometry
+        virtual void initGeometry();
+
+        //- Calculate the patch geometry
+        virtual void calcGeometry();
+
+        //- Initialise the patches for moving points
+        virtual void initMovePoints(const pointField&);
+
+        //- Correct patches after moving points
+        virtual void movePoints(const pointField&);
+
+        //- Initialise the update of the patch topology
+        virtual void initUpdateMesh();
+
+        //- Update of the patch topology
+        virtual void updateMesh();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("directMappedWall");
+
+
+    // Constructors
+
+        //- Construct from components
+        directMappedWallPolyPatch
+        (
+            const word& name,
+            const label size,
+            const label start,
+            const label index,
+            const polyBoundaryMesh& bm
+        );
+
+        //- Construct from dictionary
+        directMappedWallPolyPatch
+        (
+            const word& name,
+            const dictionary& dict,
+            const label index,
+            const polyBoundaryMesh& bm
+        );
+
+        //- Construct as copy, resetting the boundary mesh
+        directMappedWallPolyPatch
+        (
+            const directMappedWallPolyPatch&,
+            const polyBoundaryMesh&
+        );
+
+        //- Construct given the original patch and resetting the
+        //  face list and boundary mesh information
+        directMappedWallPolyPatch
+        (
+            const directMappedWallPolyPatch& pp,
+            const polyBoundaryMesh& bm,
+            const label index,
+            const label newSize,
+            const label newStart
+        );
+
+        //- Construct and return a clone, resetting the boundary mesh
+        virtual autoPtr<polyPatch> clone(const polyBoundaryMesh& bm) const
+        {
+            return autoPtr<polyPatch>(new directMappedWallPolyPatch(*this, bm));
+        }
+
+        //- Construct and return a clone, resetting the face list
+        //  and boundary mesh
+        virtual autoPtr<polyPatch> clone
+        (
+            const polyBoundaryMesh& bm,
+            const label index,
+            const label newSize,
+            const label newStart
+        ) const
+        {
+            return autoPtr<polyPatch>
+            (
+                new directMappedWallPolyPatch
+                (
+                    *this,
+                    bm,
+                    index,
+                    newSize,
+                    newStart
+                )
+            );
+        }
+
+
+    // Destructor
+
+        ~directMappedWallPolyPatch();
+
+
+    // Member functions
+
+        //- Write the polyPatch data as a dictionary
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/meshTools/indexedOctree/indexedOctree.H

@@ -128,6 +128,15 @@ public:
 
 private:
 
+    // Static data
+
+        //- Relative peturbation tolerance. Determines when point is
+        //  considered to be close to face/edge of bb of node.
+        //  The tolerance is relative to the bounding box of the smallest
+        //  node.
+        static scalar perturbTol_;
+
+
     // Private data
 
         //- Underlying shapes for geometric queries.
@@ -144,8 +153,9 @@ private:
 
     // Private Member Functions
 
-        //- Like above but now bb is implicitly provided as parent bb + mid
-        //  + octant
+        //- Helper: does bb intersect a sphere around sample? Or is any
+        //  corner point of bb closer than nearestDistSqr to sample.
+        //  (bb is implicitly provided as parent bb + octant)
         static bool overlaps
         (
             const treeBoundBox& parentBb,
@@ -215,6 +225,60 @@ private:
                 point& nearestPoint
             ) const;
 
+            //- Return bbox of octant
+            treeBoundBox subBbox
+            (
+                const label parentNodeI,
+                const direction octant
+            ) const;
+
+            //- Helper: take a point on/close to face of bb and push it
+            //  inside or outside of bb.
+            static point pushPoint
+            (
+                const treeBoundBox&,
+                const point&,
+                const bool pushInside
+            );
+
+            //- Helper: take a point on face of bb and push it
+            //  inside or outside of bb.
+            static point pushPoint
+            (
+                const treeBoundBox&,
+                const direction,
+                const point&,
+                const bool pushInside
+            );
+
+            //- Helper: take point on face(s) of bb and push it away from
+            //  edges of face.
+            static point pushPointIntoFace
+            (
+                const treeBoundBox& bb,
+                const vector& dir,          // end-start
+                const point& pt
+            );
+
+            ////- Push point on multiple faces away from any corner/edge.
+            //static void checkMultipleFaces
+            //(
+            //    const treeBoundBox& bb,
+            //    const vector& dir,          // end-start
+            //    pointIndexHit& faceHitInfo,
+            //    direction& faceID
+            //);
+
+            //- Walk to parent of node+octant.
+            bool walkToParent
+            (
+                const label nodeI,
+                const direction octant,
+
+                label& parentNodeI,
+                label& parentOctant
+            ) const;
+
             //- Walk tree to neighbouring node. Return false if edge of tree
             //  hit.
             bool walkToNeighbour
@@ -225,8 +289,8 @@ private:
                 direction& octant
             ) const;
 
-            //- Categorize point on bounding box.
-            static direction getFace(const treeBoundBox&, const point&);
+            //- Debug: return verbose the bounding box faces
+            static word faceString(const direction faceID);
 
             //- Traverse a node. If intersects a triangle return first
             // intersection point.
@@ -235,9 +299,11 @@ private:
             void traverseNode
             (
                 const bool findAny,
+                const point& treeStart,
+                const vector& treeVec,
+
                 const point& start,
                 const point& end,
-                const vector& dir,
                 const label nodeI,
                 const direction octantI,
 
@@ -252,7 +318,8 @@ private:
                 const point& treeStart,
                 const point& treeEnd,
                 const label startNodeI,
-                const direction startOctantI
+                const direction startOctantI,
+                const bool verbose = false
             ) const;
 
             //- Find any or nearest intersection of line between start and end.
@@ -310,6 +377,10 @@ private:
 
 public:
 
+        //- Get the perturbation tolerance
+        static scalar& perturbTol();
+
+
     // Constructors
 
         //- Construct null
@@ -505,6 +576,7 @@ public:
                 const point& sample
             );
 
+
         // Write
 
             //- Print tree. Either print all indices (printContent = true) or

src/meshTools/indexedOctree/treeDataCell.H

@@ -37,9 +37,7 @@ SourceFiles
 #ifndef treeDataCell_H
 #define treeDataCell_H
 
-#include "treeBoundBox.H"
 #include "treeBoundBoxList.H"
-#include "labelList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/meshTools/indexedOctree/treeDataEdge.C

@@ -22,16 +22,15 @@ 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 "treeDataEdge.H"
 #include "indexedOctree.H"
-#include "polyMesh.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
+defineTypeNameAndDebug(Foam::treeDataEdge, 0);
+
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 

src/meshTools/indexedOctree/treeDataEdge.H

@@ -36,8 +36,6 @@ SourceFiles
 #ifndef treeDataEdge_H
 #define treeDataEdge_H
 
-#include "treeBoundBox.H"
-#include "pointField.H"
 #include "treeBoundBoxList.H"
 #include "linePointRef.H"
 

src/meshTools/indexedOctree/treeDataPoint.C

@@ -32,6 +32,11 @@ Description
 #include "polyMesh.H"
 #include "triangleFuncs.H"
 
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(Foam::treeDataPoint, 0);
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct from components

src/meshTools/indexedOctree/treeDataTriSurface.C

@@ -22,8 +22,6 @@ 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 "treeDataTriSurface.H"
@@ -226,7 +224,7 @@ Foam::label Foam::treeDataTriSurface::getVolumeType
 
     if (!pHit.hit())
     {
-        FatalErrorIn("treeDataTriSurface::getVolumeType")
+        FatalErrorIn("treeDataTriSurface::getVolumeType(..)")
             << "treeBb:" << treeBb
             << " sample:" << sample
             << " pHit:" << pHit
@@ -238,7 +236,8 @@ Foam::label Foam::treeDataTriSurface::getVolumeType
         surface_,
         sample,
         pHit.index(),
-        pHit.hitPoint()
+        pHit.hitPoint(),
+        indexedOctree<treeDataTriSurface>::perturbTol()
     );
 
     if (t == triSurfaceTools::UNKNOWN)
@@ -255,12 +254,13 @@ Foam::label Foam::treeDataTriSurface::getVolumeType
     }
     else
     {
-        FatalErrorIn("treeDataTriSurface::getVolumeType")
+        FatalErrorIn("treeDataTriSurface::getVolumeType(..)")
             << "problem" << abort(FatalError);
         return indexedOctree<treeDataTriSurface>::UNKNOWN;
     }
 }
 
+
 // Check if any point on triangle is inside cubeBb.
 bool Foam::treeDataTriSurface::overlaps
 (
@@ -305,6 +305,7 @@ bool Foam::treeDataTriSurface::overlaps
     // Now we have the difficult case: all points are outside but connecting
     // edges might go through cube. Use fast intersection of bounding box.
 
+    //return triangleFuncs::intersectBbExact(p0, p1, p2, cubeBb);
     return triangleFuncs::intersectBb(p0, p1, p2, cubeBb);
 }
 
@@ -445,7 +446,15 @@ bool Foam::treeDataTriSurface::intersects
 
     const vector dir(end - start);
 
-    pointHit inter = tri.intersection(start, dir, intersection::HALF_RAY);
+    // Use relative tolerance (from octree) to determine intersection.
+
+    pointHit inter = tri.intersection
+    (
+        start,
+        dir,
+        intersection::HALF_RAY,
+        indexedOctree<treeDataTriSurface>::perturbTol()
+    );
 
     if (inter.hit() && inter.distance() <= 1)
     {
@@ -459,6 +468,16 @@ bool Foam::treeDataTriSurface::intersects
     {
         return false;
     }
+
+
+    //- Using exact intersections
+    //pointHit info = f.tri(points).intersectionExact(start, end);
+    //
+    //if (info.hit())
+    //{
+    //    intersectionPoint = info.hitPoint();
+    //}
+    //return info.hit();
 }
 
 

src/meshTools/meshSearch/meshSearch.C

@@ -27,9 +27,11 @@ License
 #include "meshSearch.H"
 #include "polyMesh.H"
 #include "indexedOctree.H"
-#include "pointIndexHit.H"
 #include "DynamicList.H"
 #include "demandDrivenData.H"
+#include "treeDataCell.H"
+#include "treeDataFace.H"
+#include "treeDataPoint.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 

src/meshTools/meshSearch/meshSearch.H

@@ -37,9 +37,6 @@ SourceFiles
 #ifndef meshSearch_H
 #define meshSearch_H
 
-#include "treeDataCell.H"
-#include "treeDataFace.H"
-#include "treeDataPoint.H"
 #include "pointIndexHit.H"
 #include "Cloud.H"
 #include "passiveParticle.H"
@@ -51,6 +48,10 @@ namespace Foam
 
 // Forward declaration of classes
 class polyMesh;
+class treeDataCell;
+class treeDataFace;
+class treeDataPoint;
+template<class Type> class indexedOctree;
 
 /*---------------------------------------------------------------------------*\
                            Class meshSearch Declaration

src/meshTools/octree/treeBoundBox.C

@@ -309,12 +309,14 @@ bool Foam::treeBoundBox::overlaps
 //     This makes sure that posBits tests 'inside'
 bool Foam::treeBoundBox::intersects
 (
+    const point& overallStart,
+    const vector& overallVec,
     const point& start,
     const point& end,
-    point& pt
+    point& pt,
+    direction& ptOnFaces
 ) const
 {
-    const vector vec(end - start);
     const direction endBits = posBits(end);
     pt = start;
 
@@ -325,80 +327,106 @@ bool Foam::treeBoundBox::intersects
         if (ptBits == 0)
         {
             // pt inside bb
+            ptOnFaces = faceBits(pt);
             return true;
         }
 
         if ((ptBits & endBits) != 0)
         {
             // pt and end in same block outside of bb
+            ptOnFaces = faceBits(pt);
             return false;
         }
 
         if (ptBits & LEFTBIT)
         {
             // Intersect with plane V=min, n=-1,0,0
-            if (Foam::mag(vec.x()) > VSMALL)
+            if (Foam::mag(overallVec.x()) > VSMALL)
             {
-                scalar s = (min().x() - pt.x())/vec.x();
+                scalar s = (min().x() - overallStart.x())/overallVec.x();
+                pt.x() = min().x();
+                pt.y() = overallStart.y() + overallVec.y()*s;
+                pt.z() = overallStart.z() + overallVec.z()*s;
+            }
+            else
+            {
+                // Vector not in x-direction. But still intersecting bb planes.
+                // So must be close - just snap to bb.
                 pt.x() = min().x();
-                pt.y() = pt.y() + vec.y()*s;
-                pt.z() = pt.z() + vec.z()*s;
             }
         }
-        if (ptBits & RIGHTBIT)
+        else if (ptBits & RIGHTBIT)
         {
             // Intersect with plane V=max, n=1,0,0
-            if (Foam::mag(vec.x()) > VSMALL)
+            if (Foam::mag(overallVec.x()) > VSMALL)
+            {
+                scalar s = (max().x() - overallStart.x())/overallVec.x();
+                pt.x() = max().x();
+                pt.y() = overallStart.y() + overallVec.y()*s;
+                pt.z() = overallStart.z() + overallVec.z()*s;
+            }
+            else
             {
-                scalar s = (max().x() - pt.x())/vec.x();
                 pt.x() = max().x();
-                pt.y() = pt.y() + vec.y()*s;
-                pt.z() = pt.z() + vec.z()*s;
             }
         }
-
-        if (ptBits & BOTTOMBIT)
+        else if (ptBits & BOTTOMBIT)
         {
             // Intersect with plane V=min, n=0,-1,0
-            if (Foam::mag(vec.y()) > VSMALL)
+            if (Foam::mag(overallVec.y()) > VSMALL)
             {
-                scalar s = (min().y() - pt.y())/vec.y();
-                pt.x() = pt.x() + vec.x()*s;
+                scalar s = (min().y() - overallStart.y())/overallVec.y();
+                pt.x() = overallStart.x() + overallVec.x()*s;
                 pt.y() = min().y();
-                pt.z() = pt.z() + vec.z()*s;
+                pt.z() = overallStart.z() + overallVec.z()*s;
+            }
+            else
+            {
+                pt.x() = min().y();
             }
         }
-        if (ptBits & TOPBIT)
+        else if (ptBits & TOPBIT)
         {
             // Intersect with plane V=max, n=0,1,0
-            if (Foam::mag(vec.y()) > VSMALL)
+            if (Foam::mag(overallVec.y()) > VSMALL)
+            {
+                scalar s = (max().y() - overallStart.y())/overallVec.y();
+                pt.x() = overallStart.x() + overallVec.x()*s;
+                pt.y() = max().y();
+                pt.z() = overallStart.z() + overallVec.z()*s;
+            }
+            else
             {
-                scalar s = (max().y() - pt.y())/vec.y();
-                pt.x() = pt.x() + vec.x()*s;
                 pt.y() = max().y();
-                pt.z() = pt.z() + vec.z()*s;
             }
         }
-
-        if (ptBits & BACKBIT)
+        else if (ptBits & BACKBIT)
         {
             // Intersect with plane V=min, n=0,0,-1
-            if (Foam::mag(vec.z()) > VSMALL)
+            if (Foam::mag(overallVec.z()) > VSMALL)
+            {
+                scalar s = (min().z() - overallStart.z())/overallVec.z();
+                pt.x() = overallStart.x() + overallVec.x()*s;
+                pt.y() = overallStart.y() + overallVec.y()*s;
+                pt.z() = min().z();
+            }
+            else
             {
-                scalar s = (min().z() - pt.z())/vec.z();
-                pt.x() = pt.x() + vec.x()*s;
-                pt.y() = pt.y() + vec.y()*s;
                 pt.z() = min().z();
             }
         }
-        if (ptBits & FRONTBIT)
+        else if (ptBits & FRONTBIT)
         {
             // Intersect with plane V=max, n=0,0,1
-            if (Foam::mag(vec.z()) > VSMALL)
+            if (Foam::mag(overallVec.z()) > VSMALL)
+            {
+                scalar s = (max().z() - overallStart.z())/overallVec.z();
+                pt.x() = overallStart.x() + overallVec.x()*s;
+                pt.y() = overallStart.y() + overallVec.y()*s;
+                pt.z() = max().z();
+            }
+            else
             {
-                scalar s = (max().z() - pt.z())/vec.z();
-                pt.x() = pt.x() + vec.x()*s;
-                pt.y() = pt.y() + vec.y()*s;
                 pt.z() = max().z();
             }
         }
@@ -406,6 +434,18 @@ bool Foam::treeBoundBox::intersects
 }
 
 
+bool Foam::treeBoundBox::intersects
+(
+    const point& start,
+    const point& end,
+    point& pt
+) const
+{
+    direction ptBits;
+    return intersects(start, end-start, start, end, pt, ptBits);
+}
+
+
 // this.bb fully contains bb
 bool Foam::treeBoundBox::contains(const treeBoundBox& bb) const
 {
@@ -452,6 +492,40 @@ bool Foam::treeBoundBox::contains(const vector& dir, const point& pt) const
 }
 
 
+// Code position of pt on bounding box faces
+Foam::direction Foam::treeBoundBox::faceBits(const point& pt) const
+{
+    direction faceBits = 0;
+    if (pt.x() == min().x())
+    {
+        faceBits |= LEFTBIT;
+    }
+    else if (pt.x() == max().x())
+    {
+        faceBits |= RIGHTBIT;
+    }
+
+    if (pt.y() == min().y())
+    {
+        faceBits |= BOTTOMBIT;
+    }
+    else if (pt.y() == max().y())
+    {
+        faceBits |= TOPBIT;
+    }
+
+    if (pt.z() == min().z())
+    {
+        faceBits |= BACKBIT;
+    }
+    else if (pt.z() == max().z())
+    {
+        faceBits |= FRONTBIT;
+    }
+    return faceBits;
+}
+
+
 // Code position of point relative to box
 Foam::direction Foam::treeBoundBox::posBits(const point& pt) const
 {
@@ -461,7 +535,7 @@ Foam::direction Foam::treeBoundBox::posBits(const point& pt) const
     {
         posBits |= LEFTBIT;
     }
-    if (pt.x() > max().x())
+    else if (pt.x() > max().x())
     {
         posBits |= RIGHTBIT;
     }
@@ -470,7 +544,7 @@ Foam::direction Foam::treeBoundBox::posBits(const point& pt) const
     {
         posBits |= BOTTOMBIT;
     }
-    if (pt.y() > max().y())
+    else if (pt.y() > max().y())
     {
         posBits |= TOPBIT;
     }
@@ -479,7 +553,7 @@ Foam::direction Foam::treeBoundBox::posBits(const point& pt) const
     {
         posBits |= BACKBIT;
     }
-    if (pt.z() > max().z())
+    else if (pt.z() > max().z())
     {
         posBits |= FRONTBIT;
     }

src/meshTools/octree/treeBoundBox.H

@@ -120,12 +120,12 @@ public:
         enum faceBit
         {
             NOFACE    = 0,
-            LEFTBIT   = 0x1 << LEFT,
-            RIGHTBIT  = 0x1 << RIGHT,
-            BOTTOMBIT = 0x1 << BOTTOM,
-            TOPBIT    = 0x1 << TOP,
-            BACKBIT   = 0x1 << BACK,
-            FRONTBIT  = 0x1 << FRONT,
+            LEFTBIT   = 0x1 << LEFT,    //1
+            RIGHTBIT  = 0x1 << RIGHT,   //2
+            BOTTOMBIT = 0x1 << BOTTOM,  //4
+            TOPBIT    = 0x1 << TOP,     //8
+            BACKBIT   = 0x1 << BACK,    //16
+            FRONTBIT  = 0x1 << FRONT,   //32
         };
 
         //- Edges codes.
@@ -265,9 +265,25 @@ public:
             //- Overlaps boundingSphere (centre + sqr(radius))?
             bool overlaps(const point&, const scalar radiusSqr) const;
 
-            //- Intersects segment; set point to intersection position,
+            //- Intersects segment; set point to intersection position and face,
             //  return true if intersection found.
-            // (intPt argument used during calculation even if not intersecting)
+            //  (pt argument used during calculation even if not intersecting).
+            //  Calculates intersections from outside supplied vector
+            //  (overallStart, overallVec). This is so when
+            //  e.g. tracking through lots of consecutive boxes
+            // (typical octree) we're not accumulating truncation errors. Set
+            // to start, (end-start) if not used.
+            bool intersects
+            (
+                const point& overallStart,
+                const vector& overallVec,
+                const point& start,
+                const point& end,
+                point& pt,
+                direction& ptBits
+            ) const;
+
+            //- Like above but does not return faces point is on
             bool intersects
             (
                 const point& start,
@@ -285,6 +301,9 @@ public:
             //  dir would cause it to go inside.
             bool contains(const vector& dir, const point&) const;
 
+            //- Code position of pt on bounding box faces
+            direction faceBits(const point& pt) const;
+
             //- Position of point relative to bounding box
             direction posBits(const point&) const;
 

src/meshTools/searchableSurface/triSurfaceMesh.C

@@ -237,6 +237,7 @@ Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io, const triSurface& s)
         )
     ),
     triSurface(s),
+    tolerance_(indexedOctree<treeDataTriSurface>::perturbTol()),
     surfaceClosed_(-1)
 {}
 
@@ -279,6 +280,7 @@ Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io)
             searchableSurface::objectPath()
         )
     ),
+    tolerance_(indexedOctree<treeDataTriSurface>::perturbTol()),
     surfaceClosed_(-1)
 {}
 
@@ -324,6 +326,7 @@ Foam::triSurfaceMesh::triSurfaceMesh
             searchableSurface::objectPath()
         )
     ),
+    tolerance_(indexedOctree<treeDataTriSurface>::perturbTol()),
     surfaceClosed_(-1)
 {
     scalar scaleFactor = 0;
@@ -332,8 +335,18 @@ Foam::triSurfaceMesh::triSurfaceMesh
     // eg, CAD geometries are often done in millimeters
     if (dict.readIfPresent("scale", scaleFactor) && scaleFactor > 0)
     {
+        Info<< searchableSurface::name() << " : using scale " << scaleFactor
+            << endl;
         triSurface::scalePoints(scaleFactor);
     }
+
+
+    // Have optional non-standard search tolerance for gappy surfaces.
+    if (dict.readIfPresent("tolerance", tolerance_) && tolerance_ > 0)
+    {
+        Info<< searchableSurface::name() << " : using intersection tolerance "
+            << tolerance_ << endl;
+    }
 }
 
 
@@ -380,6 +393,9 @@ const Foam::indexedOctree<Foam::treeDataTriSurface>&
         bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
         bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
 
+        scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
+        indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;
+
         tree_.reset
         (
             new indexedOctree<treeDataTriSurface>
@@ -391,6 +407,8 @@ const Foam::indexedOctree<Foam::treeDataTriSurface>&
                 3.0     // duplicity
             )
         );
+
+        indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
     }
 
     return tree_();
@@ -491,6 +509,9 @@ void Foam::triSurfaceMesh::findNearest
 
     info.setSize(samples.size());
 
+    scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
+    indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;
+
     forAll(samples, i)
     {
         static_cast<pointIndexHit&>(info[i]) = octree.findNearest
@@ -499,6 +520,8 @@ void Foam::triSurfaceMesh::findNearest
             nearestDistSqr[i]
         );
     }
+
+    indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
 }
 
 
@@ -513,6 +536,9 @@ void Foam::triSurfaceMesh::findLine
 
     info.setSize(start.size());
 
+    scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
+    indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;
+
     forAll(start, i)
     {
         static_cast<pointIndexHit&>(info[i]) = octree.findLine
@@ -521,6 +547,8 @@ void Foam::triSurfaceMesh::findLine
             end[i]
         );
     }
+
+    indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
 }
 
 
@@ -535,6 +563,9 @@ void Foam::triSurfaceMesh::findLineAny
 
     info.setSize(start.size());
 
+    scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
+    indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;
+
     forAll(start, i)
     {
         static_cast<pointIndexHit&>(info[i]) = octree.findLineAny
@@ -543,6 +574,8 @@ void Foam::triSurfaceMesh::findLineAny
             end[i]
         );
     }
+
+    indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
 }
 
 
@@ -557,6 +590,9 @@ void Foam::triSurfaceMesh::findLineAll
 
     info.setSize(start.size());
 
+    scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
+    indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;
+
     // Work array
     DynamicList<pointIndexHit, 1, 1> hits;
 
@@ -570,7 +606,7 @@ void Foam::triSurfaceMesh::findLineAll
     const scalarField magSqrDirVec(magSqr(dirVec));
     const vectorField smallVec
     (
-        Foam::sqrt(SMALL)*dirVec
+        indexedOctree<treeDataTriSurface>::perturbTol()*dirVec
       + vector(ROOTVSMALL,ROOTVSMALL,ROOTVSMALL)
     );
 
@@ -613,6 +649,8 @@ void Foam::triSurfaceMesh::findLineAll
             info[pointI].clear();
         }
     }
+
+    indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
 }
 
 
@@ -649,7 +687,13 @@ void Foam::triSurfaceMesh::getNormal
     {
         if (info[i].hit())
         {
-            normal[i] = faceNormals()[info[i].index()];
+            label triI = info[i].index();
+            //- Cached:
+            //normal[i] = faceNormals()[triI];
+
+            //- Uncached
+            normal[i] = triSurface::operator[](triI).normal(points());
+            normal[i] /= mag(normal[i]) + VSMALL;
         }
         else
         {
@@ -691,6 +735,9 @@ void Foam::triSurfaceMesh::getVolumeType
 {
     volType.setSize(points.size());
 
+    scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
+    indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;
+
     forAll(points, pointI)
     {
         const point& pt = points[pointI];
@@ -699,6 +746,8 @@ void Foam::triSurfaceMesh::getVolumeType
         // - cheat conversion since same values
         volType[pointI] = static_cast<volumeType>(tree().getVolumeType(pt));
     }
+
+    indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
 }
 
 

src/meshTools/searchableSurface/triSurfaceMesh.H

@@ -28,6 +28,11 @@ Class
 Description
     IOoject and searching on triSurface
 
+    Note: when constructing from dictionary has optional parameters:
+        - scale     : scaling factor.
+        - tolerance : relative tolerance for doing intersections
+                      (see triangle::intersection)
+
 SourceFiles
     triSurfaceMesh.C
 
@@ -63,6 +68,9 @@ private:
 
     // Private member data
 
+        //- Optional tolerance to use in searches
+        scalar tolerance_;
+
         //- Search tree (triangles)
         mutable autoPtr<indexedOctree<treeDataTriSurface> > tree_;
 

src/meshTools/triSurface/orientedSurface/orientedSurface.C

@@ -121,8 +121,8 @@ Foam::labelList Foam::orientedSurface::edgeToFace
             label face0 = eFaces[0];
             label face1 = eFaces[1];
 
-            const labelledTri& f0 = s[face0];
-            const labelledTri& f1 = s[face1];
+            const labelledTri& f0 = s.localFaces()[face0];
+            const labelledTri& f1 = s.localFaces()[face1];
 
             if (flip[face0] == UNVISITED)
             {
@@ -238,7 +238,8 @@ void Foam::orientedSurface::propagateOrientation
         s,
         samplePoint,
         nearestFaceI,
-        nearestPt
+        nearestPt,
+        10*SMALL
     );
 
     if (side == triSurfaceTools::UNKNOWN)
@@ -348,7 +349,10 @@ Foam::orientedSurface::orientedSurface
 :
     triSurface(surf)
 {
-    point outsidePoint = 2 * treeBoundBox(localPoints()).span();
+    // BoundBox calculation without localPoints
+    treeBoundBox bb(surf.points(), surf.meshPoints());
+
+    point outsidePoint = bb.max() + bb.span();
 
     orient(*this, outsidePoint, orientOutside);
 }

src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C

@@ -2203,7 +2203,8 @@ Foam::triSurfaceTools::sideType Foam::triSurfaceTools::surfaceSide
     const triSurface& surf,
     const point& sample,
     const label nearestFaceI,   // nearest face
-    const point& nearestPoint   // nearest point on nearest face
+    const point& nearestPoint,  // nearest point on nearest face
+    const scalar tol
 )
 {
     const labelledTri& f = surf[nearestFaceI];
@@ -2217,7 +2218,7 @@ Foam::triSurfaceTools::sideType Foam::triSurfaceTools::surfaceSide
         points[f[0]],
         points[f[1]],
         points[f[2]]
-    ).classify(nearestPoint, 1E-6, nearType, nearLabel);
+    ).classify(nearestPoint, tol, nearType, nearLabel);
 
     if (nearType == triPointRef::NONE)
     {

src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.H

@@ -454,13 +454,14 @@ public:
 
         //- Given nearest point (to sample) on surface determines which side
         //  sample is. Uses either face normal, edge normal or point normal
-        //  (non-trivial). Feed in output of e.g. triangle::classify.
+        //  (non-trivial). Uses triangle::classify.
         static sideType surfaceSide
         (
             const triSurface& surf,
             const point& sample,
             const label nearestFaceI,   // nearest face
-            const point& nearestPt      // nearest point on nearest face
+            const point& nearestPt,     // nearest point on nearest face
+            const scalar tol            // tolerance for nearness test.
         );
 
     // Triangulation of faces

src/meshTools/triSurface/triangleFuncs/triangleFuncs.C

@@ -22,8 +22,6 @@ 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 "triangleFuncs.H"
@@ -31,8 +29,6 @@ Description
 #include "treeBoundBox.H"
 #include "SortableList.H"
 #include "boolList.H"
-#include "scalar.H"
-
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -170,7 +166,6 @@ bool Foam::triangleFuncs::intersectAxesBundle
 }
 
 
-
 // Intersect triangle with bounding box. Return true if
 // any of the faces of bb intersect triangle.
 // Note: so returns false if triangle inside bb.
@@ -262,6 +257,272 @@ bool Foam::triangleFuncs::intersectBb
 }
 
 
+//// Intersect triangle with bounding box. Return true if
+//// any of the faces of bb intersect triangle.
+//// Note: so returns false if triangle inside bb.
+//bool Foam::triangleFuncs::intersectBbExact
+//(
+//    const point& p0,
+//    const point& p1,
+//    const point& p2,
+//    const treeBoundBox& cubeBb
+//)
+//{
+//    const point& min = cubeBb.min();
+//    const point& max = cubeBb.max();
+//
+//    const point& cube0 = min;
+//    const point  cube1(min.x(), min.y(), max.z());
+//    const point  cube2(max.x(), min.y(), max.z());
+//    const point  cube3(max.x(), min.y(), min.z());
+//
+//    const point  cube4(min.x(), max.y(), min.z());
+//    const point  cube5(min.x(), max.y(), max.z());
+//    const point& cube6 = max;
+//    const point  cube7(max.x(), max.y(), min.z());
+//
+//    // Test intersection of triangle with twelve edges of box.
+//    {
+//        triPointRef tri(p0, p1, p2);
+//        if (tri.intersectionExact(cube0, cube1).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(cube1, cube2).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(cube2, cube3).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(cube3, cube0).hit())
+//        {   
+//            return true;
+//        }
+//
+//        if (tri.intersectionExact(cube4, cube5).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(cube5, cube6).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(cube6, cube7).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(cube7, cube4).hit())
+//        {
+//            return true;
+//        }
+//
+//        if (tri.intersectionExact(cube0, cube4).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(cube1, cube5).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(cube2, cube6).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(cube3, cube7).hit())
+//        {
+//            return true;
+//        }
+//    }
+//    // Test intersection of triangle edges with bounding box
+//    {
+//        triPointRef tri(cube0, cube1, cube2);
+//        if (tri.intersectionExact(p0, p1).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p1, p2).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p2, p0).hit())
+//        {
+//            return true;
+//        }
+//    }
+//    {
+//        triPointRef tri(cube2, cube3, cube0);
+//        if (tri.intersectionExact(p0, p1).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p1, p2).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p2, p0).hit())
+//        {
+//            return true;
+//        }
+//    }
+//    {
+//        triPointRef tri(cube4, cube5, cube6);
+//        if (tri.intersectionExact(p0, p1).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p1, p2).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p2, p0).hit())
+//        {
+//            return true;
+//        }
+//    }
+//    {
+//        triPointRef tri(cube6, cube7, cube4);
+//        if (tri.intersectionExact(p0, p1).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p1, p2).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p2, p0).hit())
+//        {
+//            return true;
+//        }
+//    }
+//
+//
+//    {
+//        triPointRef tri(cube4, cube5, cube1);
+//        if (tri.intersectionExact(p0, p1).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p1, p2).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p2, p0).hit())
+//        {
+//            return true;
+//        }
+//    }
+//    {
+//        triPointRef tri(cube1, cube0, cube4);
+//        if (tri.intersectionExact(p0, p1).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p1, p2).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p2, p0).hit())
+//        {
+//            return true;
+//        }
+//    }
+//    {
+//        triPointRef tri(cube7, cube6, cube2);
+//        if (tri.intersectionExact(p0, p1).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p1, p2).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p2, p0).hit())
+//        {
+//            return true;
+//        }
+//    }
+//    {
+//        triPointRef tri(cube2, cube3, cube7);
+//        if (tri.intersectionExact(p0, p1).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p1, p2).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p2, p0).hit())
+//        {
+//            return true;
+//        }
+//    }
+//
+//    {
+//        triPointRef tri(cube0, cube4, cube7);
+//        if (tri.intersectionExact(p0, p1).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p1, p2).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p2, p0).hit())
+//        {
+//            return true;
+//        }
+//    }
+//    {
+//        triPointRef tri(cube7, cube3, cube0);
+//        if (tri.intersectionExact(p0, p1).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p1, p2).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p2, p0).hit())
+//        {
+//            return true;
+//        }
+//    }
+//    {
+//        triPointRef tri(cube1, cube5, cube6);
+//        if (tri.intersectionExact(p0, p1).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p1, p2).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p2, p0).hit())
+//        {
+//            return true;
+//        }
+//    }
+//    {
+//        triPointRef tri(cube6, cube2, cube1);
+//        if (tri.intersectionExact(p0, p1).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p1, p2).hit())
+//        {
+//            return true;
+//        }
+//        if (tri.intersectionExact(p2, p0).hit())
+//        {
+//            return true;
+//        }
+//    }
+//    return false;
+//}
+
+
 bool Foam::triangleFuncs::intersect
 (
     const point& va0,
@@ -470,145 +731,4 @@ bool Foam::triangleFuncs::intersect
 }
 
 
-bool Foam::triangleFuncs::classify
-(
-    const point& baseVertex,
-    const vector& E0,
-    const vector& E1,
-    const vector& n,
-    const point& pInter,
-    const scalar tol,
-    label& nearType,
-    label& nearLabel
-)
-{
-    // Get largest component of normal
-    scalar magX = Foam::mag(n.x());
-    scalar magY = Foam::mag(n.y());
-    scalar magZ = Foam::mag(n.z());
-
-    label i0 = -1;
-    if ((magX >= magY) && (magX >= magZ))
-    {
-        i0 = 0;
-    }
-    else if ((magY >= magX) && (magY >= magZ))
-    {
-        i0 = 1;
-    }
-    else
-    {
-        i0 = 2;
-    }
-
-    // Get other components
-    label i1 = (i0 + 1) % 3;
-    label i2 = (i1 + 1) % 3;
-
-    scalar u1 = E0[i1];
-    scalar v1 = E0[i2];
-
-    scalar u2 = E1[i1];
-    scalar v2 = E1[i2];
-
-    scalar det = v2*u1 - u2*v1;
-
-    scalar u0 = pInter[i1] - baseVertex[i1];
-    scalar v0 = pInter[i2] - baseVertex[i2];
-
-    scalar alpha = 0;
-    scalar beta = 0;
-
-    bool hit = false;
-    
-    if (Foam::mag(u1) < ROOTVSMALL)
-    {
-        beta = u0/u2;
-
-        alpha = (v0 - beta*v2)/v1;
-
-        hit = ((alpha >= 0) && ((alpha + beta) <= 1));
-    }
-    else
-    {
-        beta = (v0*u1 - u0*v1)/det;
-
-        alpha = (u0 - beta*u2)/u1;
-
-        hit = ((alpha >= 0) && ((alpha + beta) <= 1));
-    }
-
-    //
-    // Now alpha, beta are the coordinates in the triangle local coordinate
-    // system E0, E1
-    //
-
-    nearType = NONE;
-    nearLabel = -1;
-
-    if (mag(alpha+beta-1) < tol)
-    {
-        // On edge between vert 1-2 (=edge 1)
-
-        if (mag(alpha) < tol)
-        {
-            nearType = POINT;
-            nearLabel = 2;
-        }
-        else if (mag(beta) < tol)
-        {
-            nearType = POINT;
-            nearLabel = 1;
-        }
-        else if ((alpha >= 0) && (alpha <= 1) && (beta >= 0) && (beta <= 1))
-        {
-            nearType = EDGE;
-            nearLabel = 1;
-        }
-    }
-    else if (mag(alpha) < tol)
-    {
-        // On edge between vert 2-0 (=edge 2)
-
-        if (mag(beta) < tol)
-        {
-            nearType = POINT;
-            nearLabel = 0;
-        }
-        else if (mag(beta-1) < tol)
-        {
-            nearType = POINT;
-            nearLabel = 2;
-        }
-        else if ((beta >= 0) && (beta <= 1))
-        {
-            nearType = EDGE;
-            nearLabel = 2;
-        }
-    }
-    else if (mag(beta) < tol)
-    {
-        // On edge between vert 0-1 (= edge 0)
-
-        if (mag(alpha) < tol)
-        {
-            nearType = POINT;
-            nearLabel = 0;
-        }
-        else if (mag(alpha-1) < tol)
-        {
-            nearType = POINT;
-            nearLabel = 1;
-        }
-        else if ((alpha >= 0) && (alpha <= 1))
-        {
-            nearType = EDGE;
-            nearLabel = 0;
-        }
-    }
-
-    return hit;
-}
-
-
 // ************************************************************************* //

src/meshTools/triSurface/triangleFuncs/triangleFuncs.H

@@ -147,25 +147,6 @@ public:
         point& pInter0,
         point& pInter1
     );
-
-
-    //- Classify point on triangle plane w.r.t. triangle edges.
-    //  - inside/outside (returned)
-    //  - near point (0, 1, 2)
-    //  - near edge (0, 1, 2). Note: edges are counted from starting vertex
-    //    so e.g. edge 2 is from f[2] to f[0]
-    static bool classify
-    (
-        const point& baseVertex,
-        const vector& E0,
-        const vector& E1,
-        const vector& n,
-        const point& pInter,
-        const scalar tol,
-        label& nearType,
-        label& nearLabel
-    );
-
 };
 
 

src/triSurface/triSurface/triSurface.C

@@ -347,127 +347,6 @@ void Foam::triSurface::checkEdges(const bool verbose)
 }
 
 
-// Check normals and orientation
-Foam::boolList Foam::triSurface::checkOrientation(const bool verbose)
-{
-    const edgeList& es = edges();
-    const labelListList& faceEs = faceEdges();
-
-    // Check edge normals, face normals, point normals.
-    forAll(faceEs, facei)
-    {
-        const labelList& edgeLabels = faceEs[facei];
-
-        if (edgeLabels.size() != 3)
-        {
-            FatalErrorIn("triSurface::checkOrientation(bool)")
-                << "triangle " << (*this)[facei]
-                << " does not have 3 edges. Edges:" << edgeLabels
-                << exit(FatalError);
-        }
-
-        bool valid = true;
-        forAll(edgeLabels, i)
-        {
-            if (edgeLabels[i] < 0 || edgeLabels[i] >= nEdges())
-            {
-                WarningIn
-                (
-                    "triSurface::checkOrientation(bool)"
-                )   << "edge number " << edgeLabels[i] << " on face " << facei
-                    << " out-of-range\n"
-                    << "This usually means that the input surface has "
-                    << "edges with more than 2 triangles connected.\n"
-                    << endl;
-                valid = false;
-            }
-        }
-        if (! valid)
-        {
-            continue;
-        }
-
-
-        //
-        //- Compute normal from triangle points.
-        //
-
-        const labelledTri& tri = (*this)[facei];
-        const point pa(points()[tri[0]]);
-        const point pb(points()[tri[1]]);
-        const point pc(points()[tri[2]]);
-
-        const vector pointNormal((pc - pb) ^ (pa - pb));
-        if ((pointNormal & faceNormals()[facei]) < 0)
-        {
-            FatalErrorIn("triSurface::checkOrientation(bool)")
-                << "Normal calculated from points not consistent with"
-                " faceNormal" << endl
-                << "triangle:" << tri << endl
-                << "points:" << pa << ' ' << pb << ' ' << pc << endl
-                << "pointNormal:" << pointNormal << endl
-                << "faceNormal:" << faceNormals()[facei]
-                << exit(FatalError);
-        }
-    }
-
-
-    const labelListList& eFaces = edgeFaces();
-
-    // Storage for holding status of edge. True if normal flips across this
-    // edge
-    boolList borderEdge(nEdges(), false);
-
-    forAll(es, edgeI)
-    {
-        const edge& e = es[edgeI];
-        const labelList& neighbours = eFaces[edgeI];
-
-        if (neighbours.size() == 2)
-        {
-            const labelledTri& faceA = (*this)[neighbours[0]];
-            const labelledTri& faceB = (*this)[neighbours[1]];
-
-            // The edge cannot be going in the same direction if both faces
-            // are oriented counterclockwise.
-            // Thus the next face point *must* different between the faces.
-            if
-            (
-                faceA[faceA.fcIndex(findIndex(faceA, e.start()))]
-             == faceB[faceB.fcIndex(findIndex(faceB, e.start()))]
-            )
-            {
-                borderEdge[edgeI] = true;
-                if (verbose)
-                {
-                    WarningIn("PrimitivePatchExtra::checkOrientation(bool)")
-                        << "face orientation incorrect." << nl
-                        << "edge[" << edgeI << "] " << e
-                        << " between faces " << neighbours << ":" << nl
-                        << "face[" << neighbours[0] << "] " << faceA << nl
-                        << "face[" << neighbours[1] << "] " << faceB << endl;
-                }
-            }
-        }
-        else if (neighbours.size() != 1)
-        {
-            if (verbose)
-            {
-                WarningIn("triSurface::checkOrientation(bool)")
-                    << "Wrong number of edge neighbours." << endl
-                    << "edge[" << edgeI << "] " << e
-                    << "with points:" << localPoints()[e.start()]
-                    << ' ' << localPoints()[e.end()]
-                    << " has neighbours:" << neighbours << endl;
-            }
-            borderEdge[edgeI] = true;
-        }
-    }
-
-    return borderEdge;
-}
-
-
 // Read triangles, points from Istream
 bool Foam::triSurface::read(Istream& is)
 {

src/triSurface/triSurface/triSurface.H

@@ -326,17 +326,12 @@ public:
             //- Scale points. A non-positive factor is ignored
             virtual void scalePoints(const scalar&);
 
-            //- Check/fix duplicate/degenerate triangles
+            //- Check/remove duplicate/degenerate triangles
             void checkTriangles(const bool verbose);
 
-            //- Check triply (or more) connected edges. Return list of faces
-            //  sharing these edges.
+            //- Check triply (or more) connected edges.
             void checkEdges(const bool verbose);
 
-            //- Check orientation (normals) and normals of neighbouring
-            //  triangles
-            boolList checkOrientation(const bool verbose);
-
             //- Remove non-valid triangles
             void cleanup(const bool verbose);
 

src/turbulenceModels/compressible/RAS/backwardsCompatibilityWallFunctions/backwardsCompatibilityWallFunctions.C

@@ -64,8 +64,8 @@ tmp<volScalarField> autoCreateAlphat
     }
     else
     {
-        Info<< "--> Upgrading " << fieldName << " to employ run-time "
-            << "selectable wall functions" << endl;
+        Info<< "--> Creating " << fieldName
+            << " to employ run-time selectable wall functions" << endl;
 
         const fvBoundaryMesh& bm = mesh.boundary();
 
@@ -104,7 +104,7 @@ tmp<volScalarField> autoCreateAlphat
             )
         );
 
-        Info<< "    Writing updated " << fieldName << endl;
+        Info<< "    Writing new " << fieldName << endl;
         alphat().write();
 
         return alphat;
@@ -134,8 +134,8 @@ tmp<volScalarField> autoCreateMut
     }
     else
     {
-        Info<< "--> Upgrading " << fieldName << " to employ run-time "
-            << "selectable wall functions" << endl;
+        Info<< "--> Creating " << fieldName
+            << " to employ run-time selectable wall functions" << endl;
 
         const fvBoundaryMesh& bm = mesh.boundary();
 
@@ -174,7 +174,7 @@ tmp<volScalarField> autoCreateMut
             )
         );
 
-        Info<< "    Writing updated " << fieldName << endl;
+        Info<< "    Writing new " << fieldName << endl;
         mut().write();
 
         return mut;

src/turbulenceModels/compressible/RAS/backwardsCompatibilityWallFunctions/backwardsCompatibilityWallFunctionsTemplates.C

@@ -26,6 +26,7 @@ License
 
 #include "backwardsCompatibilityWallFunctions.H"
 #include "Time.H"
+#include "OSspecific.H"
 
 #include "wallPolyPatch.H"
 
@@ -77,27 +78,35 @@ autoCreateWallFunctionField
     }
     else
     {
-        Info<< "--> Upgrading " << fieldName << " to employ run-time "
-            << "selectable wall functions" << endl;
+        Info<< "--> Upgrading " << fieldName
+            << " to employ run-time selectable wall functions" << endl;
+
+        // Read existing field
+        IOobject ioObj
+        (
+            fieldName,
+            mesh.time().timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE,
+            false
+        );
 
-        // Read existing epsilon field
         tmp<fieldType> fieldOrig
         (
             new fieldType
             (
-                IOobject
-                (
-                    fieldName,
-                    mesh.time().timeName(),
-                    mesh,
-                    IOobject::MUST_READ,
-                    IOobject::NO_WRITE,
-                    false
-                ),
+                ioObj,
                 mesh
             )
         );
 
+        // rename file
+        Info<< "    Backup original " << fieldName << " to "
+            << fieldName << ".old" << endl;
+        mv(ioObj.objectPath(), ioObj.objectPath() + ".old");
+
+
         PtrList<fvPatchField<Type> > newPatchFields(mesh.boundary().size());
 
         forAll(newPatchFields, patchI)
@@ -145,11 +154,6 @@ autoCreateWallFunctionField
             )
         );
 
-        Info<< "    Writing backup of original " << fieldName << " to "
-            << fieldName << ".old" << endl;
-        fieldOrig().rename(fieldName + ".old");
-        fieldOrig().write();
-
         Info<< "    Writing updated " << fieldName << endl;
         fieldNew().write();