Merge commit 'OpenCFD/master' into olesenm

applications/solvers/discreteMethods/dsmc/dsmcFoam/dsmcFoam.C

@@ -91,7 +91,7 @@ int main(int argc, char *argv[])
 
         runTime.write();
 
-        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+        Info<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
             << nl << endl;
     }

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 =
+    volScalarField DpDt
-        fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+    (
+        "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/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
+        //- Construct by mapping given
-        //  onto a new patch
+        //  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/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/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 (newPatches[0] == -1)
-
-    if (wantedPatchI == -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:
+    // Do the actual changes
-    // - 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);
-
 
     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];
+            createBaffles
-            label zoneID = faceZones.whichZone(faceI);
-            bool zoneFlip = false;
-            if (zoneID >= 0)
-            {
-                const faceZone& fZone = faceZones[zoneID];
-                zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
-            }
-
-            meshMod.setAction
             (
-                polyModifyFace
+                mesh,
-                (
+                faceI,
-                    f,                          // modified face
+                -1,                             // oldPatchI,
-                    faceI,                      // label of face
+                newPatches,
-                    mesh.faceOwner()[faceI],    // owner
+                doneFace,
-                    -1,                         // neighbour
+                meshMod
-                    false,                      // face flip
-                    newPatch[faceI],            // patch for face
-                    false,                      // remove from zone
-                    zoneID,                     // zone for face
-                    zoneFlip                    // face flip in zone
-                )
             );
-
-            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/postProcessing/sampling/sample/sample.C

@@ -97,10 +97,11 @@ using namespace Foam;
 int main(int argc, char *argv[])
 {
     timeSelector::addOptions();
+#   include "addRegionOption.H"
 #   include "setRootCase.H"
 #   include "createTime.H"
     instantList timeDirs = timeSelector::select0(runTime, args);
-#   include "createMesh.H"
+#   include "createNamedMesh.H"
 
     IOsampledSets sSets
     (

applications/utilities/postProcessing/sampling/sample/sampleDict

@@ -105,6 +105,14 @@ sets
         end         (2 0.51  0.005);
         nPoints     10;
     }
+
+    somePoints
+    {
+        type    cloud;
+        axis    xyz;
+        points  ((0.049 0.049 0.005)(0.051 0.049 0.005));
+    }
+
 );
 
 

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;
 
@@ -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) : "

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/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

@@ -58,7 +58,13 @@ void Foam::mapDistribute::distribute
         }
 
         // Subset myself
-        UIndirectList<T> subField(field, subMap[Pstream::myProcNo()]);
+        const labelList& mySubMap = subMap[Pstream::myProcNo()];
+
+        List<T> subField(mySubMap.size());
+        forAll(mySubMap, i)
+        {
+            subField[i] = field[mySubMap[i]];
+        }
 
         // Receive sub field from myself (subField)
         const labelList& map = constructMap[Pstream::myProcNo()];
@@ -355,7 +361,13 @@ void Foam::mapDistribute::distribute
         }
 
         // Subset myself
-        UIndirectList<T> subField(field, subMap[Pstream::myProcNo()]);
+        const labelList& mySubMap = subMap[Pstream::myProcNo()];
+
+        List<T> subField(mySubMap.size());
+        forAll(mySubMap, i)
+        {
+            subField[i] = field[mySubMap[i]];
+        }
 
         // Receive sub field from myself (subField)
         const labelList& map = constructMap[Pstream::myProcNo()];

src/OpenFOAM/meshes/polyMesh/polyMeshInitMesh.C

@@ -82,10 +82,10 @@ void Foam::polyMesh::initMesh()
     );
 
     string meshInfo =
-        "nPoints: " + Foam::name(nPoints())
+        "nPoints:" + Foam::name(nPoints())
-      + " nCells: " + Foam::name(this->nCells())
+      + "  nCells:" + Foam::name(this->nCells())
-      + " nFaces: " + Foam::name(nFaces())
+      + "  nFaces:" + Foam::name(nFaces())
-      + " nInternalFaces: " + Foam::name(nInternalFaces());
+      + "  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& p0 = p.points()[f[0]];
-        const point p1(p.points()[f[1]]);
+        const point& p1 = p.points()[f[1]];
-        const point p2(p.points()[f[f.size()-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"
+                    << "Normal calculated from points inconsistent"
-                    << nl
+                    << " 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/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/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
+    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
-    const directMappedPolyPatch& mpp = refCast<const directMappedPolyPatch>
+
+    // 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 =
+            const label nbrPatchID = nbrMesh.boundaryMesh().findPatchID
-                 this->patch().patch().boundaryMesh().findPatchID
             (
                 mpp.samplePatch()
             );
-            if (patchID < 0)
+            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;
+            Field<Type> allValues(nbrMesh.nFaces(), pTraits<Type>::zero);
-            const word& fieldName = this->dimensionedInternalField().name();
-            const fieldType& sendField =
-                this->db().objectRegistry::lookupObject<fieldType>(fieldName);
 
-            Field<Type> allValues
+            const fieldType& nbrField = nbrMesh.lookupObject<fieldType>
             (
-                this->patch().patch().boundaryMesh().mesh().nFaces(),
+                fldName
-                pTraits<Type>::zero
             );
-
+            forAll(nbrField.boundaryField(), patchI)
-            forAll(sendField.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
+    // Get the directMappedPatchBase
-    const directMappedPolyPatch& mpp = refCast<const directMappedPolyPatch>
+    const directMappedPatchBase& mpp = refCast<const directMappedPatchBase>
     (
         directMappedVelocityFluxFixedValueFvPatchField::patch().patch()
     );
-
+    const mapDistribute& distMap = mpp.map();
-    vectorField newUValues(size());
+    const fvMesh& nbrMesh = refCast<const fvMesh>(mpp.sampleMesh());
-    scalarField newPhiValues(size());
-
     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
+            vectorField allUValues(nbrMesh.nFaces(), vector::zero);
-            (
+            scalarField allPhiValues(nbrMesh.nFaces(), 0.0);
-                patch().patch().boundaryMesh().mesh().nFaces(),
-                vector::zero
-            );
-            scalarField allPhiValues
-            (
-                patch().patch().boundaryMesh().mesh().nFaces(),
-                0.0
-            );
 
             forAll(UField.boundaryField(), patchI)
             {
@@ -273,34 +197,58 @@ void directMappedVelocityFluxFixedValueFvPatchField::updateCoeffs()
                 }
             }
 
-            getNewValues
+            mapDistribute::distribute
             (
-                mpp,
+                Pstream::defaultCommsType,
-                allUValues,
+                distMap.schedule(),
-                allPhiValues,
+                distMap.constructSize(),
-                newUValues,
+                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 =
+            const label nbrPatchID = nbrMesh.boundaryMesh().findPatchID
-                patch().patch().boundaryMesh().findPatchID
             (
                 mpp.samplePatch()
             );
 
-            getNewValues
+            newUValues = UField.boundaryField()[nbrPatchID];
+
+            mapDistribute::distribute
             (
-                mpp,
+                Pstream::defaultCommsType,
-                UField.boundaryField()[patchID],
+                distMap.schedule(),
-                phiField.boundaryField()[patchID],
+                distMap.constructSize(),
-                newUValues,
+                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/lagrangian/basic/Particle/Particle.C

@@ -390,7 +390,7 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
     // slightly towards the cell-centre.
     if (trackFraction < SMALL)
     {
-        position_ += 1.0e-6*(mesh.cellCentres()[celli_] - position_);
+        position_ += 1.0e-3*(mesh.cellCentres()[celli_] - position_);
     }
 
     return trackFraction;

src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C

@@ -37,6 +37,7 @@ Foam::scalar Foam::DsmcCloud<ParcelType>::kb = 1.380650277e-23;
 template<class ParcelType>
 Foam::scalar Foam::DsmcCloud<ParcelType>::Tref = 273;
 
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 template<class ParcelType>
@@ -113,14 +114,57 @@ void Foam::DsmcCloud<ParcelType>::initialise
 
     numberDensities /= nParticle_;
 
-    scalar x0 = mesh_.bounds().min().x();
+    forAll(mesh_.cells(), cell)
-    scalar xR = mesh_.bounds().max().x() - x0;
+    {
+        const vector& cC = mesh_.cellCentres()[cell];
+        const labelList& cellFaces = mesh_.cells()[cell];
+        const scalar cV = mesh_.cellVolumes()[cell];
 
-    scalar y0 = mesh_.bounds().min().y();
+        label nTets = 0;
-    scalar yR = mesh_.bounds().max().y() - y0;
 
-    scalar z0 = mesh_.bounds().min().z();
+        // Each face is split into nEdges (or nVertices) - 2 tets.
-    scalar zR = mesh_.bounds().max().z() - z0;
+        forAll(cellFaces, face)
+        {
+            nTets += mesh_.faces()[cellFaces[face]].size() - 2;
+        }
+
+        // Calculate the cumulative tet volumes circulating around the cell and
+        // record the vertex labels of each.
+        scalarList cTetVFracs(nTets, 0.0);
+
+        List<labelList> tetPtIs(nTets, labelList(3,-1));
+
+        // Keep track of which tet this is.
+        label tet = 0;
+
+        forAll(cellFaces, face)
+        {
+            const labelList& facePoints = mesh_.faces()[cellFaces[face]];
+
+            label pointI = 1;
+            while ((pointI + 1) < facePoints.size())
+            {
+
+                const vector& pA = mesh_.points()[facePoints[0]];
+                const vector& pB = mesh_.points()[facePoints[pointI]];
+                const vector& pC = mesh_.points()[facePoints[pointI + 1]];
+
+                cTetVFracs[tet] =
+                    mag(((pA - cC) ^ (pB - cC)) & (pC - cC))/(cV*6.0)
+                  + cTetVFracs[max((tet - 1),0)];
+
+                tetPtIs[tet][0] = facePoints[0];
+                tetPtIs[tet][1] = facePoints[pointI];
+                tetPtIs[tet][2] = facePoints[pointI + 1];
+
+                pointI++;
+                tet++;
+            }
+        }
+
+        // Force the last volume fraction value to 1.0 to avoid any
+        // rounding/non-flat face errors giving a value < 1.0
+        cTetVFracs[nTets - 1] = 1.0;
 
         forAll(molecules, i)
         {
@@ -135,38 +179,73 @@ void Foam::DsmcCloud<ParcelType>::initialise
                     << abort(FatalError);
             }
 
-        const typename ParcelType::constantProperties& cP = constProps(typeId);
+            const typename ParcelType::constantProperties& cP =
+                constProps(typeId);
 
             scalar numberDensity = numberDensities[i];
 
-        scalar spacing = pow(numberDensity,-(1.0/3.0));
+            // Calculate the number of particles required
+            scalar particlesRequired = numberDensity*mesh_.cellVolumes()[cell];
 
-        int ni = label(xR/spacing) + 1;
+            // Only integer numbers of particles can be inserted
-        int nj = label(yR/spacing) + 1;
+            label nParticlesToInsert = label(particlesRequired);
-        int nk = label(zR/spacing) + 1;
 
-        vector delta(xR/ni, yR/nj, zR/nk);
+            // Add another particle with a probability proportional to the
-
+            // remainder of taking the integer part of particlesRequired
-        scalar pert = spacing;
+            if ((particlesRequired - nParticlesToInsert) > rndGen_.scalar01())
-
-        for (int i = 0; i < ni; i++)
             {
-            for (int j = 0; j < nj; j++)
+                nParticlesToInsert++;
-            {
+            }
-                for (int k = 0; k < nk; k++)
-                {
-                    point p
-                    (
-                        x0 + (i + 0.5)*delta.x(),
-                        y0 + (j + 0.5)*delta.y(),
-                        z0 + (k + 0.5)*delta.z()
-                    );
 
-                    p.x() += pert*(rndGen_.scalar01() - 0.5);
+            for (label pI = 0; pI < nParticlesToInsert; pI++)
-                    p.y() += pert*(rndGen_.scalar01() - 0.5);
+            {
-                    p.z() += pert*(rndGen_.scalar01() - 0.5);
+                // Choose a random point in a generic tetrahedron
 
-                    label cell = mesh_.findCell(p);
+                scalar s = rndGen_.scalar01();
+                scalar t = rndGen_.scalar01();
+                scalar u = rndGen_.scalar01();
+
+                if (s + t > 1.0)
+                {
+                    s = 1.0 - s;
+                    t = 1.0 - t;
+                }
+
+                if (t + u > 1.0)
+                {
+                    scalar tmp = u;
+                    u = 1.0 - s - t;
+                    t = 1.0 - tmp;
+                }
+                else if (s + t + u > 1.0)
+                {
+                    scalar tmp = u;
+                    u = s + t + u - 1.0;
+                    s = 1.0 - t - tmp;
+                }
+
+                // Choose a tetrahedron to insert in, based on their relative
+                // volumes
+                scalar tetSelection = rndGen_.scalar01();
+
+                // Selected tetrahedron
+                label sTet = -1;
+
+                forAll(cTetVFracs, tet)
+                {
+                    sTet = tet;
+
+                    if (cTetVFracs[tet] >= tetSelection)
+                    {
+                        break;
+                    }
+                }
+
+                vector p =
+                    (1 - s - t - u)*cC
+                  + s*mesh_.points()[tetPtIs[sTet][0]]
+                  + t*mesh_.points()[tetPtIs[sTet][1]]
+                  + u*mesh_.points()[tetPtIs[sTet][2]];
 
                 vector U = equipartitionLinearVelocity
                 (
@@ -182,8 +261,6 @@ void Foam::DsmcCloud<ParcelType>::initialise
 
                 U += velocity;
 
-                    if (cell >= 0)
-                    {
                 addNewParcel
                 (
                     p,
@@ -195,8 +272,6 @@ void Foam::DsmcCloud<ParcelType>::initialise
             }
         }
     }
-        }
-    }
 
     // Initialise the sigmaTcRMax_ field to the product of the cross section of
     // the most abundant species and the most probable thermal speed (Bird,
@@ -276,7 +351,7 @@ void Foam::DsmcCloud<ParcelType>::collisions()
 
             scalar selectedPairs =
                 collisionSelectionRemainder_[celli]
-              + 0.5*nC*(nC-1)*nParticle_*sigmaTcRMax*deltaT
+              + 0.5*nC*(nC - 1)*nParticle_*sigmaTcRMax*deltaT
                /mesh_.cellVolumes()[celli];
 
             label nCandidates(selectedPairs);
@@ -551,6 +626,13 @@ Foam::DsmcCloud<ParcelType>::DsmcCloud
     buildConstProps();
 
     buildCellOccupancy();
+
+    // Initialise the collision selection remainder to a random value between 0
+    // and 1.
+    forAll(collisionSelectionRemainder_, i)
+    {
+        collisionSelectionRemainder_[i] = rndGen_.scalar01();
+    }
 }
 
 
@@ -739,8 +821,12 @@ void Foam::DsmcCloud<ParcelType>::info() const
 
     Info<< "Cloud name: " << this->name() << nl
         << "    Number of dsmc particles        = "
-        << nDsmcParticles << nl
+        << nDsmcParticles
-        << "    Number of molecules             = "
+        << endl;
+
+    if (nDsmcParticles)
+    {
+        Info<< "    Number of molecules             = "
             << nMol << nl
             << "    Mass in system                  = "
             << returnReduce(massInSystem(), sumOp<scalar>()) << nl
@@ -753,8 +839,9 @@ void Foam::DsmcCloud<ParcelType>::info() const
             << "    Average internal energy         = "
             << internalEnergy/nMol << nl
             << "    Average total energy            = "
-        << (internalEnergy + linearKineticEnergy)/nMol << nl
+            << (internalEnergy + linearKineticEnergy)/nMol
             << endl;
+    }
 }
 
 
@@ -785,7 +872,11 @@ Foam::scalar Foam::DsmcCloud<ParcelType>::equipartitionInternalEnergy
 {
     scalar Ei = 0.0;
 
-    if (iDof < 2.0 + SMALL && iDof > 2.0 - SMALL)
+    if (iDof < SMALL)
+    {
+        return Ei;
+    }
+    else if (iDof < 2.0 + SMALL && iDof > 2.0 - SMALL)
     {
         // Special case for iDof = 2, i.e. diatomics;
         Ei = -log(rndGen_.scalar01())*kb*temperature;
@@ -796,7 +887,7 @@ Foam::scalar Foam::DsmcCloud<ParcelType>::equipartitionInternalEnergy
 
         scalar energyRatio;
 
-        scalar P;
+        scalar P = -1;
 
         do
         {

src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H

@@ -299,6 +299,12 @@ public:
                     scalar mass
                 ) const;
 
+                inline scalarField maxwellianAverageSpeed
+                (
+                    scalarField temperature,
+                    scalar mass
+                ) const;
+
                 //- RMS particle speed
                 inline scalar maxwellianRMSSpeed
                 (
@@ -306,6 +312,12 @@ public:
                     scalar mass
                 ) const;
 
+                inline scalarField maxwellianRMSSpeed
+                (
+                    scalarField temperature,
+                    scalar mass
+                ) const;
+
                 //- Most probable speed
                 inline scalar maxwellianMostProbableSpeed
                 (
@@ -313,6 +325,11 @@ public:
                     scalar mass
                 ) const;
 
+                inline scalarField maxwellianMostProbableSpeed
+                (
+                    scalarField temperature,
+                    scalar mass
+                ) const;
 
             // Sub-models
 

src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H

@@ -302,6 +302,17 @@ inline Foam::scalar Foam::DsmcCloud<ParcelType>::maxwellianAverageSpeed
 }
 
 
+template<class ParcelType>
+inline Foam::scalarField Foam::DsmcCloud<ParcelType>::maxwellianAverageSpeed
+(
+    scalarField temperature,
+    scalar mass
+) const
+{
+    return 2.0*sqrt(2.0*kb*temperature/(mathematicalConstant::pi*mass));
+}
+
+
 template<class ParcelType>
 inline Foam::scalar Foam::DsmcCloud<ParcelType>::maxwellianRMSSpeed
 (
@@ -313,6 +324,17 @@ inline Foam::scalar Foam::DsmcCloud<ParcelType>::maxwellianRMSSpeed
 }
 
 
+template<class ParcelType>
+inline Foam::scalarField Foam::DsmcCloud<ParcelType>::maxwellianRMSSpeed
+(
+    scalarField temperature,
+    scalar mass
+) const
+{
+    return sqrt(3.0*kb*temperature/mass);
+}
+
+
 template<class ParcelType>
 inline Foam::scalar
 Foam::DsmcCloud<ParcelType>::maxwellianMostProbableSpeed
@@ -325,6 +347,18 @@ Foam::DsmcCloud<ParcelType>::maxwellianMostProbableSpeed
 }
 
 
+template<class ParcelType>
+inline Foam::scalarField
+Foam::DsmcCloud<ParcelType>::maxwellianMostProbableSpeed
+(
+    scalarField temperature,
+    scalar mass
+) const
+{
+    return sqrt(2.0*kb*temperature/mass);
+}
+
+
 template<class ParcelType>
 inline const Foam::tmp<Foam::volScalarField>
 Foam::DsmcCloud<ParcelType>::rhoN() const
@@ -343,7 +377,7 @@ Foam::DsmcCloud<ParcelType>::rhoN() const
                 false
             ),
             mesh_,
-            dimensionedScalar("zero", dimensionSet(0, -3, 0, 0, 0), 0.0)
+            dimensionedScalar("zero", dimensionSet(0, -3, 0, 0, 0), VSMALL)
         )
     );
 
@@ -380,7 +414,7 @@ Foam::DsmcCloud<ParcelType>::rhoM() const
                 false
             ),
             mesh_,
-            dimensionedScalar("zero", dimensionSet(1, -3, 0, 0, 0), 0.0)
+            dimensionedScalar("zero", dimensionSet(1, -3, 0, 0, 0), VSMALL)
         )
     );
 
@@ -568,7 +602,7 @@ Foam::DsmcCloud<ParcelType>::iDof() const
                 false
             ),
             mesh_,
-            dimensionedScalar("zero", dimensionSet(0, -3, 0, 0, 0), 0.0)
+            dimensionedScalar("zero", dimensionSet(0, -3, 0, 0, 0), VSMALL)
         )
     );
 

src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C

@@ -36,31 +36,26 @@ Foam::FreeStream<CloudType>::FreeStream
 )
 :
     InflowBoundaryModel<CloudType>(dict, cloud, typeName),
-    patchIndex_(),
+    patches_(),
-    temperature_(readScalar(this->coeffDict().lookup("temperature"))),
-    velocity_(this->coeffDict().lookup("velocity")),
     moleculeTypeIds_(),
     numberDensities_(),
     particleFluxAccumulators_()
 {
-    word patchName = this->coeffDict().lookup("patch");
+    // Identify which patches to use
 
-    patchIndex_ = cloud.mesh().boundaryMesh().findPatchID(patchName);
+    DynamicList<label> patches;
 
-    const polyPatch& patch = cloud.mesh().boundaryMesh()[patchIndex_];
+    forAll(cloud.mesh().boundaryMesh(), p)
-
-    if (patchIndex_ == -1)
     {
-        FatalErrorIn
+        const polyPatch& patch = cloud.mesh().boundaryMesh()[p];
-        (
+
-            "Foam::FreeStream<CloudType>::FreeStream"
+        if (patch.type() == polyPatch::typeName)
-            "("
+        {
-                "const dictionary&, "
+            patches.append(p);
-                "CloudType&"
-            ")"
-        )   << "patch " << patchName << " not found." << nl
-            << abort(FatalError);
         }
+    }
+
+    patches_.transfer(patches);
 
     const dictionary& numberDensitiesDict
     (
@@ -69,10 +64,24 @@ Foam::FreeStream<CloudType>::FreeStream
 
     List<word> molecules(numberDensitiesDict.toc());
 
-    numberDensities_.setSize(molecules.size());
+    // Initialise the particleFluxAccumulators_
+    particleFluxAccumulators_.setSize(patches_.size());
+
+    forAll(patches_, p)
+    {
+        const polyPatch& patch = cloud.mesh().boundaryMesh()[patches_[p]];
+
+        particleFluxAccumulators_[p] = List<Field<scalar> >
+        (
+            molecules.size(),
+            Field<scalar>(patch.size(), 0.0)
+        );
+    }
 
     moleculeTypeIds_.setSize(molecules.size());
 
+    numberDensities_.setSize(molecules.size());
+
     forAll(molecules, i)
     {
         numberDensities_[i] = readScalar
@@ -97,12 +106,6 @@ Foam::FreeStream<CloudType>::FreeStream
     }
 
     numberDensities_ /= cloud.nParticle();
-
-    particleFluxAccumulators_.setSize
-    (
-        molecules.size(),
-        Field<scalar>(patch.size(), 0)
-    );
 }
 
 
@@ -127,17 +130,65 @@ void Foam::FreeStream<CloudType>::inflow()
 
     Random& rndGen(cloud.rndGen());
 
-    const polyPatch& patch = mesh.boundaryMesh()[patchIndex_];
+    scalar sqrtPi = sqrt(mathematicalConstant::pi);
 
     label particlesInserted = 0;
 
+    const volScalarField::GeometricBoundaryField& boundaryT
+    (
+        cloud.T().boundaryField()
+    );
+
+    const volVectorField::GeometricBoundaryField& boundaryU
+    (
+        cloud.U().boundaryField()
+    );
+
+    forAll(patches_, p)
+    {
+        label patchI = patches_[p];
+
+        const polyPatch& patch = mesh.boundaryMesh()[patchI];
+
         // Add mass to the accumulators.  negative face area dotted with the
         // velocity to point flux into the domain.
 
-    forAll(particleFluxAccumulators_, i)
+        // Take a reference to the particleFluxAccumulator for this patch
+        List<Field<scalar> >& pFA = particleFluxAccumulators_[p];
+
+        forAll(pFA, i)
         {
-        particleFluxAccumulators_[i] +=
+            label typeId = moleculeTypeIds_[i];
-           -patch.faceAreas() & (velocity_*numberDensities_[i]*deltaT);
+
+            scalar mass = cloud.constProps(typeId).mass();
+
+            scalarField mostProbableSpeed
+            (
+                cloud.maxwellianMostProbableSpeed
+                (
+                    boundaryT[patchI],
+                    mass
+                )
+            );
+
+            // Dotting boundary velocity with the face unit normal (which points
+            // out of the domain, so it must be negated), dividing by the most
+            // probable speed to form molecularSpeedRatio * cosTheta
+
+            scalarField sCosTheta =
+                boundaryU[patchI]
+              & -patch.faceAreas()/mag(patch.faceAreas())
+               /mostProbableSpeed;
+
+            // From Bird eqn 4.22
+
+            pFA[i] +=
+                mag(patch.faceAreas())*numberDensities_[i]*deltaT
+               *mostProbableSpeed
+               *(
+                   exp(-sqr(sCosTheta)) + sqrtPi*sCosTheta*(1 + erf(sCosTheta))
+                )
+               /(2.0*sqrtPi);
         }
 
         forAll(patch, f)
@@ -157,7 +208,7 @@ void Foam::FreeStream<CloudType>::inflow()
 
             scalar fA = mag(patch.faceAreas()[f]);
 
-        // Cummulative triangle area fractions
+            // Cumulative triangle area fractions
             List<scalar> cTriAFracs(nVertices);
 
             for (label v = 0; v < nVertices - 1; v++)
@@ -177,35 +228,47 @@ void Foam::FreeStream<CloudType>::inflow()
 
             // Normal unit vector *negative* so normal is pointing into the
             // domain
-        vector nw = patch.faceAreas()[f];
+            vector n = patch.faceAreas()[f];
-        nw /= -mag(nw);
+            n /= -mag(n);
 
             // Wall tangential unit vector. Use the direction between the
             // face centre and the first vertex in the list
-        vector tw1 = fC - (mesh.points()[faceVertices[0]]);
+            vector t1 = fC - (mesh.points()[faceVertices[0]]);
-        tw1 /= mag(tw1);
+            t1 /= mag(t1);
 
             // Other tangential unit vector.  Rescaling in case face is not
-        // flat and nw and tw1 aren't perfectly orthogonal
+            // flat and n and t1 aren't perfectly orthogonal
-        vector tw2 = nw^tw1;
+            vector t2 = n^t1;
-        tw2 /= mag(tw2);
+            t2 /= mag(t2);
 
-        forAll(particleFluxAccumulators_, i)
+            scalar faceTemperature = boundaryT[patchI][f];
+
+            const vector& faceVelocity = boundaryU[patchI][f];
+
+            forAll(pFA, i)
             {
-            scalar& faceAccumulator = particleFluxAccumulators_[i][f];
+                scalar& faceAccumulator = pFA[i][f];
 
-            // Number of particles to insert
+                // Number of whole particles to insert
                 label nI = max(label(faceAccumulator), 0);
 
+                // Add another particle with a probability proportional to the
+                // remainder of taking the integer part of faceAccumulator
+                if ((faceAccumulator - nI) > rndGen.scalar01())
+                {
+                    nI++;
+                }
+
                 faceAccumulator -= nI;
 
                 label typeId = moleculeTypeIds_[i];
 
                 scalar mass = cloud.constProps(typeId).mass();
 
-            for (label n = 0; n < nI; n++)
+                for (label i = 0; i < nI; i++)
                 {
-                // Choose a triangle to insert on, based on their relative area
+                    // Choose a triangle to insert on, based on their relative
+                    // area
 
                     scalar triSelection = rndGen.scalar01();
 
@@ -239,19 +302,76 @@ void Foam::FreeStream<CloudType>::inflow()
 
                     point p = (1 - t)*A + (1 - s)*t*B + s*t*C;
 
-                vector U =
+                    // Velocity generation
-                    sqrt(CloudType::kb*temperature_/mass)
+
-                   *(
+                    scalar mostProbableSpeed
-                        rndGen.GaussNormal()*tw1
+                    (
-                      + rndGen.GaussNormal()*tw2
+                        cloud.maxwellianMostProbableSpeed
-                      - sqrt(-2.0*log(max(1 - rndGen.scalar01(), VSMALL)))*nw
+                        (
+                            faceTemperature,
+                            mass
+                        )
                     );
 
-                U += velocity_;
+                    scalar sCosTheta = (faceVelocity & n)/mostProbableSpeed;
+
+                    // Coefficients required for Bird eqn 12.5
+                    scalar uNormProbCoeffA =
+                        sCosTheta + sqrt(sqr(sCosTheta) + 2.0);
+
+                    scalar uNormProbCoeffB =
+                        0.5*
+                        (
+                            1.0
+                          + sCosTheta*(sCosTheta - sqrt(sqr(sCosTheta) + 2.0))
+                        );
+
+                    // Equivalent to the QA value in Bird's DSMC3.FOR
+                    scalar randomScaling = 3.0;
+
+                    if (sCosTheta < -3)
+                    {
+                        randomScaling = mag(sCosTheta) + 1;
+                    }
+
+                    scalar P = -1;
+
+                    // Normalised candidates for the normal direction velocity
+                    // component
+                    scalar uNormal;
+                    scalar uNormalThermal;
+
+                    // Select a velocity using Bird eqn 12.5
+                    do
+                    {
+                        uNormalThermal =
+                            randomScaling*(2.0*rndGen.scalar01() - 1);
+
+                        uNormal = uNormalThermal + sCosTheta;
+
+                        if (uNormal < 0.0)
+                        {
+                            P = -1;
+                        }
+                        else
+                        {
+                            P = 2.0*uNormal/uNormProbCoeffA
+                               *exp(uNormProbCoeffB - sqr(uNormalThermal));
+                        }
+
+                    } while (P < rndGen.scalar01());
+
+                    vector U =
+                        sqrt(CloudType::kb*faceTemperature/mass)
+                       *(
+                            rndGen.GaussNormal()*t1
+                          + rndGen.GaussNormal()*t2
+                        )
+                      + mostProbableSpeed*uNormal*n;
 
                     scalar Ei = cloud.equipartitionInternalEnergy
                     (
-                    temperature_,
+                        faceTemperature,
                         cloud.constProps(typeId).internalDegreesOfFreedom()
                     );
 
@@ -268,6 +388,7 @@ void Foam::FreeStream<CloudType>::inflow()
                 }
             }
         }
+    }
 
     reduce(particlesInserted, sumOp<label>());
 

src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.H

@@ -26,8 +26,10 @@ Class
     Foam::FreeStream
 
 Description
-    Inserting new particles across the faces of a specified patch for a free
+    Inserting new particles across the faces of a all patched of type
-    stream.  Uniform values of temperature, velocity and number densities
+    "patch" for a free stream.  Uniform values number density, temperature
+    and velocity sourced face-by-face from the boundaryT and boundaryU fields
+    of the cloud.
 
 \*---------------------------------------------------------------------------*/
 
@@ -52,14 +54,8 @@ class FreeStream
 {
     // Private data
 
-        //- Index of patch to introduce particles across
+        //- The indices of patches to introduce molecules across
-        label patchIndex_;
+        labelList patches_;
-
-        //- Temperature of the free stream
-        scalar temperature_;
-
-        //- Velocity of the free stream
-        vector velocity_;
 
         //- The molecule types to be introduced
         List<label> moleculeTypeIds_;
@@ -67,10 +63,13 @@ class FreeStream
         //- The number density of the species in the inflow
         Field<scalar> numberDensities_;
 
-        //- A List of Fields, one Field for every species to be introduced, each
+        //- A List of Lists of Fields specifying carry-over of mass flux from
-        // field entry corresponding to a face on the patch to be injected
+        // one timestep to the next
-        // across.
+        // + Outer List - one inner List for each patch
-        List<Field<scalar> > particleFluxAccumulators_;
+        // + Inner List - one Field for every species to be introduced
+        // + Each field entry corresponding to a face to be injected across
+        //   with a particular species
+        List<List<Field<scalar> > > particleFluxAccumulators_;
 
 
 public:

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,14 +26,6 @@ License
 
 #include "directMappedPolyPatch.H"
 #include "addToRunTimeSelectionTable.H"
-#include "ListListOps.H"
-#include "meshSearch.H"
-#include "mapDistribute.H"
-#include "meshTools.H"
-#include "OFstream.H"
-#include "Random.H"
-#include "treeDataFace.H"
-#include "indexedOctree.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -43,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  * * * * * * * * * * * * * * * * * * //
 
@@ -519,12 +53,7 @@ Foam::directMappedPolyPatch::directMappedPolyPatch
 )
 :
     polyPatch(name, size, start, index, bm),
-    mode_(NEARESTPATCHFACE),
+    directMappedPatchBase(static_cast<const polyPatch&>(*this))
-    samplePatch_("none"),
-    offset_(vector::zero),
-    schedulePtr_(NULL),
-    sendLabelsPtr_(NULL),
-    receiveFaceLabelsPtr_(NULL)
 {}
 
 
@@ -537,12 +66,7 @@ Foam::directMappedPolyPatch::directMappedPolyPatch
 )
 :
     polyPatch(name, dict, index, bm),
-    mode_(sampleModeNames_.read(dict.lookup("sampleMode"))),
+    directMappedPatchBase(*this, dict)
-    samplePatch_(dict.lookup("samplePatch")),
-    offset_(dict.lookup("offset")),
-    schedulePtr_(NULL),
-    sendLabelsPtr_(NULL),
-    receiveFaceLabelsPtr_(NULL)
 {}
 
 
@@ -553,12 +77,7 @@ Foam::directMappedPolyPatch::directMappedPolyPatch
 )
 :
     polyPatch(pp, bm),
-    mode_(pp.mode_),
+    directMappedPatchBase(*this, pp)
-    samplePatch_(pp.samplePatch_),
-    offset_(pp.offset_),
-    schedulePtr_(NULL),
-    sendLabelsPtr_(NULL),
-    receiveFaceLabelsPtr_(NULL)
 {}
 
 
@@ -572,12 +91,7 @@ Foam::directMappedPolyPatch::directMappedPolyPatch
 )
 :
     polyPatch(pp, bm, index, newSize, newStart),
-    mode_(pp.mode_),
+    directMappedPatchBase(*this, pp)
-    samplePatch_(pp.samplePatch_),
-    offset_(pp.offset_),
-    schedulePtr_(NULL),
-    sendLabelsPtr_(NULL),
-    receiveFaceLabelsPtr_(NULL)
 {}
 
 
@@ -585,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();
+    polyPatch::initGeometry();
-    sendLabelsPtr_.clear();
+    directMappedPatchBase::clearOut();
-    receiveFaceLabelsPtr_.clear();
+}
+
+//- 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_]
+    directMappedPatchBase::write(os);
-        << token::END_STATEMENT << nl;
-    os.writeKeyword("samplePatch") << samplePatch_
-        << token::END_STATEMENT << nl;
-    os.writeKeyword("offset") << offset_ << token::END_STATEMENT << nl;
 }
 
 

src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.H

@@ -42,9 +42,7 @@ SourceFiles
 #define directMappedPolyPatch_H
 
 #include "polyPatch.H"
-#include "labelPair.H"
+#include "directMappedPatchBase.H"
-#include "Tuple2.H"
-#include "pointIndexHit.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()
+        virtual void initGeometry();
-        {
-            clearOut();
-        }
 
         //- Calculate the patch geometry
-        virtual void calcGeometry()
+        virtual void calcGeometry();
-        {
-            clearOut();
-        }
 
         //- Initialise the patches for moving points
-        virtual void initMovePoints(const pointField&)
+        virtual void initMovePoints(const pointField&);
-        {
-            clearOut();
-        }
 
         //- Correct patches after moving points
-        virtual void movePoints(const pointField& p)
+        virtual void movePoints(const pointField&);
-        {
-            clearOut();
-        }
 
         //- Initialise the update of the patch topology
-        virtual void initUpdateMesh()
+        virtual void initUpdateMesh();
-        {
-            clearOut();
-        }
 
         //- Update of the patch topology
-        virtual void updateMesh()
+        virtual void updateMesh();
-        {
-            clearOut();
-        }
 
 
 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/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& f0 = s.localFaces()[face0];
-            const labelledTri& f1 = s[face1];
+            const labelledTri& f1 = s.localFaces()[face1];
 
             if (flip[face0] == UNVISITED)
             {
@@ -349,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/postProcessing/functionObjects/utilities/dsmcFields/dsmcFields.C

@@ -137,7 +137,7 @@ void Foam::dsmcFields::write()
             iDofMeanName
         );
 
-        if (min(rhoNMean).value() > VSMALL)
+        if (min(mag(rhoNMean)).value() > VSMALL)
         {
             Info<< "Calculating dsmcFields." << endl;
 
@@ -223,10 +223,9 @@ void Foam::dsmcFields::write()
         }
         else
         {
-            Info<< "Small or negative value (" << min(rhoNMean)
+            Info<< "Small value (" << min(mag(rhoNMean))
                 << ") found in rhoNMean field. "
-                << "Not calculating dsmcFields to avoid division by zero "
+                << "Not calculating dsmcFields to avoid division by zero."
-                << "or invalid results."
                 << endl;
         }
     }

src/sampling/sampledSet/sampledSets/sampledSets.C

@@ -236,8 +236,6 @@ Foam::sampledSets::sampledSets
     loadFromFiles_(loadFromFiles),
     outputPath_(fileName::null),
     searchEngine_(mesh_, true),
-//    pMeshPtr_(NULL),
-//    pInterpPtr_(NULL),
     fieldNames_(),
     interpolationScheme_(word::null),
     writeFormat_(word::null)
@@ -250,6 +248,10 @@ Foam::sampledSets::sampledSets
     {
         outputPath_ = mesh_.time().path()/name_;
     }
+    if (mesh_.name() != fvMesh::defaultRegion)
+    {
+        outputPath_ = outputPath_/mesh_.name();
+    }
 
     read(dict);
 }

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
+            //- Check triply (or more) connected edges.
-            //  sharing these 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);
 

tutorials/compressible/rhoPisoFoam/ras/cavity/constant/polyMesh/boundary

@@ -1,8 +1,8 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.5                                   |
+|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
 FoamFile
@@ -10,32 +10,31 @@ FoamFile
     version     2.0;
     format      ascii;
     class       polyBoundaryMesh;
+    location    "constant/polyMesh";
     object      boundary;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 3
 (
-movingWall
+    movingWall
-{
+    {
         type            wall;
         nFaces          20;
         startFace       760;
-}
+    }
-
+    fixedWalls
-fixedWalls
+    {
-{
         type            wall;
         nFaces          60;
         startFace       780;
-}
+    }
-
+    frontAndBack
-frontAndBack
+    {
-{
         type            empty;
         nFaces          800;
         startFace       840;
-}
+    }
 )
 
 // ************************************************************************* //

tutorials/heatTransfer/buoyantPisoFoam/hotRoom/0/epsilon

@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-dimensions      [ 0 2 -3 0 0 0 0 ];
+dimensions      [0 2 -3 0 0 0 0];
 
 internalField   uniform 0.01;
 
@@ -23,20 +23,18 @@ boundaryField
 {
     floor
     {
-        type            epsilonWallFunction;
+        type            compressible::epsilonWallFunction;
-        value           uniform 0;
+        value           uniform 0.01;
     }
-
     ceiling
     {
-        type            epsilonWallFunction;
+        type            compressible::epsilonWallFunction;
-        value           uniform 0;
+        value           uniform 0.01;
     }
-
     fixedWalls
     {
-        type            epsilonWallFunction;
+        type            compressible::epsilonWallFunction;
-        value           uniform 0;
+        value           uniform 0.01;
     }
 }
 

tutorials/heatTransfer/buoyantPisoFoam/hotRoom/0/k

@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-dimensions      [ 0 2 -2 0 0 0 0 ];
+dimensions      [0 2 -2 0 0 0 0];
 
 internalField   uniform 0.1;
 
@@ -23,20 +23,18 @@ boundaryField
 {
     floor
     {
-        type            kQRWallFunction;
+        type            compressible::kQRWallFunction;
-        value           uniform 0;
+        value           uniform 0.1;
     }
-
     ceiling
     {
-        type            kQRWallFunction;
+        type            compressible::kQRWallFunction;
-        value           uniform 0;
+        value           uniform 0.1;
     }
-
     fixedWalls
     {
-        type            kQRWallFunction;
+        type            compressible::kQRWallFunction;
-        value           uniform 0;
+        value           uniform 0.1;
     }
 }
 

tutorials/heatTransfer/buoyantPisoFoam/hotRoom/0/mut

@@ -1,44 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       volScalarField;
-    location    "0";
-    object      mut;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [ 1 -1 -1 0 0 0 0 ];
-
-internalField   uniform 0;
-
-boundaryField
-{
-    floor
-    {
-        type            mutWallFunction;
-        value           uniform 0;
-    }
-
-    ceiling
-    {
-        type            mutWallFunction;
-        value           uniform 0;
-    }
-
-    fixedWalls
-    {
-        type            mutWallFunction;
-        value           uniform 0;
-    }
-}
-
-
-// ************************************************************************* //

tutorials/heatTransfer/buoyantPisoFoam/hotRoom/constant/polyMesh/boundary

@@ -1,8 +1,8 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.5                                   |
+|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
 FoamFile
@@ -10,6 +10,7 @@ FoamFile
     version     2.0;
     format      ascii;
     class       polyBoundaryMesh;
+    location    "constant/polyMesh";
     object      boundary;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

tutorials/heatTransfer/buoyantPisoFoam/hotRoom/system/controlDict

@@ -43,5 +43,8 @@ timePrecision   6;
 
 runTimeModifiable yes;
 
+adjustTimeStep  no;
+
+maxCo           0.5;
 
 // ************************************************************************* //

tutorials/heatTransfer/buoyantPisoFoam/hotRoom/system/fvSchemes

@@ -33,6 +33,7 @@ divSchemes
     div(phi,k)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
     div(phi,R)      Gauss upwind;
+    div(phiU,p)     Gauss linear;
     div(R)          Gauss linear;
     div((muEff*dev2(grad(U).T()))) Gauss linear;
 }

tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/0/epsilon

@@ -1,8 +1,8 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.5                                   |
+|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
 FoamFile
@@ -23,17 +23,17 @@ boundaryField
 {
     floor
     {
-        type            epsilonWallFunction;
+        type            compressible::epsilonWallFunction;
         value           uniform 0.01;
     }
     ceiling
     {
-        type            epsilonWallFunction;
+        type            compressible::epsilonWallFunction;
         value           uniform 0.01;
     }
     fixedWalls
     {
-        type            epsilonWallFunction;
+        type            compressible::epsilonWallFunction;
         value           uniform 0.01;
     }
 }

tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/0/k

@@ -1,8 +1,8 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.5                                   |
+|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
 FoamFile
@@ -23,17 +23,17 @@ boundaryField
 {
     floor
     {
-        type            kQRWallFunction;
+        type            compressible::kQRWallFunction;
         value           uniform 0.1;
     }
     ceiling
     {
-        type            kQRWallFunction;
+        type            compressible::kQRWallFunction;
         value           uniform 0.1;
     }
     fixedWalls
     {
-        type            kQRWallFunction;
+        type            compressible::kQRWallFunction;
         value           uniform 0.1;
     }
 }

tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/0/mut

@@ -1,42 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.5                                   |
-|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       volScalarField;
-    location    "0";
-    object      mut;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [1 -1 -1 0 0 0 0];
-
-internalField   uniform 0;
-
-boundaryField
-{
-    floor
-    {
-        type            mutWallFunction;
-        value           uniform 0;
-    }
-    ceiling
-    {
-        type            mutWallFunction;
-        value           uniform 0;
-    }
-    fixedWalls
-    {
-        type            mutWallFunction;
-        value           uniform 0;
-    }
-}
-
-
-// ************************************************************************* //

tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/Allrun

@@ -4,7 +4,7 @@
 
 application="buoyantSimpleFoam"
 
-compileApplication ../../buoyantFoam/hotRoom/setHotRoom
+compileApplication ../../buoyantPisoFoam/hotRoom/setHotRoom
 runApplication blockMesh
 runApplication setHotRoom
 runApplication $application

tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/constant/polyMesh/boundary

@@ -1,8 +1,8 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.5                                   |
+|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
 FoamFile

tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/0/mut

@@ -1,50 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       volScalarField;
-    location    "0";
-    object      mut;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [ 1 -1 -1 0 0 0 0 ];
-
-internalField   uniform 0;
-
-boundaryField
-{
-    box
-    {
-        type            mutWallFunction;
-        value           uniform 0;
-    }
-
-    floor
-    {
-        type            mutWallFunction;
-        value           uniform 0;
-    }
-
-    ceiling
-    {
-        type            mutWallFunction;
-        value           uniform 0;
-    }
-
-    fixedWalls
-    {
-        type            mutWallFunction;
-        value           uniform 0;
-    }
-}
-
-
-// ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/0/alphat

@@ -21,6 +21,11 @@ internalField   uniform 0;
 
 boundaryField
 {
+    minY
+    {
+        type            alphatWallFunction;
+        value           uniform 0;
+    }
     maxY
     {
         type            alphatWallFunction;

tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/0/epsilon

@@ -21,6 +21,11 @@ internalField   uniform 0.01;
 
 boundaryField
 {
+    minY
+    {
+        type            epsilonWallFunction;
+        value           uniform 0.01;
+    }
     maxY
     {
         type            epsilonWallFunction;

tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/0/k

@@ -21,6 +21,11 @@ internalField   uniform 0.1;
 
 boundaryField
 {
+    minY
+    {
+        type            kQRWallFunction;
+        value           uniform 0.1;
+    }
     maxY
     {
         type            kQRWallFunction;

tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/0/mut

@@ -1,67 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       volScalarField;
-    location    "0.001";
-    object      mut;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [1 -1 -1 0 0 0 0];
-
-internalField   uniform 0;
-
-boundaryField
-{
-    maxY
-    {
-        type            mutWallFunction;
-        value           uniform 0;
-    }
-    minX
-    {
-        type            calculated;
-        value           uniform 0;
-    }
-    maxX
-    {
-        type            calculated;
-        value           uniform 0;
-    }
-    minZ
-    {
-        type            mutWallFunction;
-        value           uniform 0;
-    }
-    maxZ
-    {
-        type            mutWallFunction;
-        value           uniform 0;
-    }
-    topAir_to_leftSolid
-    {
-        type            calculated;
-        value           uniform 0;
-    }
-    topAir_to_heater
-    {
-        type            calculated;
-        value           uniform 0;
-    }
-    topAir_to_rightSolid
-    {
-        type            calculated;
-        value           uniform 0;
-    }
-}
-
-
-// ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/bottomAir/fvSolution

@@ -16,26 +16,30 @@ FoamFile
 
 solvers
 {
-    rho PCG
+    rho
     {
+        solver          PCG
         preconditioner  DIC;
         tolerance       1e-6;
         relTol          0;
     };
-//    pd PCG
+//    pd
 //    {
+//        solver           PCG
 //        preconditioner   DIC;
 //        tolerance        1e-6;
 //        relTol           0.1;
 //    };
-//    pdFinal PCG
+//    pdFinal
 //    {
+//        solver           PCG;
 //        preconditioner   DIC;
 //        tolerance        1e-08;
 //        relTol           0;
 //    };
-    pd GAMG
+    pd
     {
+        solver           GAMG;
         tolerance        1e-6;
         relTol           0.1;
 
@@ -46,8 +50,9 @@ solvers
         agglomerator     faceAreaPair;
         mergeLevels      1;
     };
-    pdFinal GAMG
+    pdFinal
     {
+        solver           GAMG;
         tolerance        1e-6;
         relTol           0;
 
@@ -58,32 +63,37 @@ solvers
         agglomerator     faceAreaPair;
         mergeLevels      1;
     };
-    U PBiCG
+    U
     {
+        solver           PBiCG;
         preconditioner   DILU;
         tolerance        1e-08;
         relTol           0;
     };
-    h PBiCG
+    h
     {
+        solver           PBiCG;
         preconditioner   DILU;
         tolerance        1e-06;
         relTol           0;
     };
-    k PBiCG
+    k
     {
+        solver           PBiCG;
         preconditioner   DILU;
         tolerance        1e-06;
         relTol           0;
     };
-    epsilon PBiCG
+    epsilon
     {
+        solver           PBiCG;
         preconditioner   DILU;
         tolerance        1e-06;
         relTol           0;
     };
-    R PBiCG
+    R
     {
+        solver           PBiCG;
         preconditioner   DILU;
         tolerance        1e-06;
         relTol           0;

tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/fvSchemes

@@ -14,5 +14,33 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+ddtSchemes
+{
+}
+
+gradSchemes
+{
+}
+
+divSchemes
+{
+}
+
+laplacianSchemes
+{
+}
+
+interpolationSchemes
+{
+}
+
+snGradSchemes
+{
+}
+
+fluxRequired
+{
+}
+
 
 // ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/heater/fvSolution

@@ -16,8 +16,9 @@ FoamFile
 
 solvers
 {
-    T PCG
+    T
     {
+        solver           PCG;
         preconditioner   DIC;
         tolerance        1E-06;
         relTol           0;

tutorials/incompressible/pisoFoam/les/pitzDailyDirectMapped/constant/polyMesh/boundary

@@ -20,7 +20,7 @@ FoamFile
 (
 inlet
 {
-    type            directMappedPatch;
+    type            directMapped;
     nFaces          30;
     startFace       27238;
     sampleMode      nearestCell;

tutorials/incompressible/pisoFoam/les/pitzDailyDirectMapped/system/changeDictionaryDict

@@ -23,6 +23,7 @@ dictionaryReplacement
         {
             type            directMappedPatch;
             offset          ( 0.0495 0 0 );
+            sampleRegion    region0;
             sampleMode      nearestCell;
             samplePatch     none;
         }

wmake/rules/General/bison

@@ -1,8 +1,8 @@
 .SUFFIXES: .y .Y
 
-ytoo = bison -v $(YYPREFIX) -d $$SOURCE ; mv $$SOURCE_DIR/*.tab.c $*.c ; mv $$SOURCE_DIR/*.tab.h $*.h ; $(cc) $(cFLAGS) -c $*.c -o $@ 
+ytoo = bison -v -d -y $$SOURCE ; mv y.tab.c $*.c ; mv y.tab.h $*.h ; $(cc) $(cFLAGS) -c $*.c -o $@ 
 
-Ytoo = bison -v $(YYPREFIX) -d $$SOURCE ; mv $$SOURCE_DIR/*.tab.c $*.C ; mv $$SOURCE_DIR/*.tab.h $*.H ; $(CC) $(c++FLAGS) -c $*.C -o $@ 
+Ytoo = bison -v -d -y $$SOURCE ; mv y.tab.c $*.C ; mv y.tab.h $*.H ; $(CC) $(c++FLAGS) -c $*.C -o $@ 
 
 .y.dep:
 	$(MAKE_DEP)

wmake/rules/General/flex

@@ -1,6 +1,6 @@
 .SUFFIXES: .l
 
-ltoo = flex $$SOURCE ; mv lex.yy.c $*.C ; $(CC) $(c++FLAGS) -c $*.C -o $@ 
+ltoo = flex $$SOURCE ; mv lex.yy.c $*.c ; $(cc) $(cFLAGS) -c $*.c -o $@ 
 
 .l.dep:
 	$(MAKE_DEP)

wmake/rules/General/standard

@@ -5,6 +5,7 @@ include $(GENERAL_RULES)/sourceToDep
 include $(GENERAL_RULES)/java
 include $(GENERAL_RULES)/flex
 include $(GENERAL_RULES)/flex++
-include $(GENERAL_RULES)/byacc
+#include $(GENERAL_RULES)/byacc
-include $(GENERAL_RULES)/btyacc++
+#include $(GENERAL_RULES)/btyacc++
+include $(GENERAL_RULES)/bison
 include $(GENERAL_RULES)/moc