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

.gitignore

@@ -42,6 +42,7 @@ linuxming*/
 SiCortex*Gcc*/
 solaris*Gcc*/
 SunOS*Gcc*/
+platforms/
 
 # reinstate wmake/rules that might look like build folders
 !wmake/rules/*/

ReleaseNotes-dev

@@ -104,6 +104,11 @@
       taking film into account
     + Parallel aware
 *** *New* ptscotch decomposition method
+*** *Updated* scotch decomposition method to run in parallel by doing
+    decomposition on the master. Unfortunately scotch and ptscotch cannot
+    be linked in to the same executable.
+*** *Updated* simple decomposition method to run in parallel by doing
+    decomposition on the master.
 *** *Updated* decomposePar maps polyPatches instead of recreating them so
     polyPatches holding data can map the data.
 *** *Updated* particle tracking algorithm
@@ -183,6 +188,9 @@
       e.g. pitzDailyDirectMapped tutorial.
     + =setSet=: allows time range (e.g. 0:100) in combination with -batch argument
     to execute the commands for multiple times.
+    + =extrudeMesh=:
+            - option to add extrusion to existing mesh.
+            - works in parallel
 * Post-processing
   + =paraFoam=, =foamToVTK=: full support for polyhedral cell type in recent
      Paraview versions.

applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H

@@ -2,8 +2,8 @@
     word scheme("div(phi,alpha)");
     word schemer("div(phir,alpha)");
 
-    surfaceScalarField phic(phi);
+    surfaceScalarField phic("phic", phi);
-    surfaceScalarField phir(phia - phib);
+    surfaceScalarField phir("phir", phia - phib);
 
     if (g0.value() > 0.0)
     {

applications/test/globalMeshData/Test-globalMeshData.C

@@ -70,7 +70,12 @@ int main(int argc, char *argv[])
         );
 
         // Exchange data. Apply positional transforms.
-        globalPointSlavesMap.distribute(transforms, coords, true);
+        globalPointSlavesMap.distribute
+        (
+            transforms,
+            coords,
+            mapDistribute::transformPosition()
+        );
 
         // Print
         forAll(slaves, pointI)
@@ -127,7 +132,12 @@ int main(int argc, char *argv[])
         }
 
         // Exchange data Apply positional transforms.
-        globalEdgeSlavesMap.distribute(transforms, ec, true);
+        globalEdgeSlavesMap.distribute
+        (
+            transforms,
+            ec,
+            mapDistribute::transformPosition()
+        );
 
         // Print
         forAll(slaves, edgeI)
@@ -163,82 +173,131 @@ int main(int argc, char *argv[])
     }
 
 
-    //// Test: (collocated) point to faces addressing
+    // Test: point to faces addressing
-    //{
+    {
-    //    const labelListList& globalPointBoundaryFaces =
+        const mapDistribute& globalPointBoundaryFacesMap =
-    //        globalData.globalPointBoundaryFaces();
+            globalData.globalPointBoundaryFacesMap();
-    //    const mapDistribute& globalPointBoundaryFacesMap =
+        const labelListList& slaves =
-    //        globalData.globalPointBoundaryFacesMap();
+            globalData.globalPointBoundaryFaces();
-    //
+        const labelListList& transformedSlaves =
-    //    label nBnd = mesh.nFaces()-mesh.nInternalFaces();
+            globalData.globalPointTransformedBoundaryFaces();
-    //
+
-    //    pointField fc(globalPointBoundaryFacesMap.constructSize());
+        label nBnd = mesh.nFaces()-mesh.nInternalFaces();
-    //    SubList<point>(fc, nBnd).assign
+
-    //    (
+        pointField fc(globalPointBoundaryFacesMap.constructSize());
-    //        primitivePatch
+        SubList<point>(fc, nBnd).assign
-    //        (
+        (
-    //            SubList<face>
+            primitivePatch
-    //            (
+            (
-    //                mesh.faces(),
+                SubList<face>
-    //                nBnd,
+                (
-    //                mesh.nInternalFaces()
+                    mesh.faces(),
-    //            ),
+                    nBnd,
-    //            mesh.points()
+                    mesh.nInternalFaces()
-    //        ).faceCentres()
+                ),
-    //    );
+                mesh.points()
-    //
+            ).faceCentres()
-    //    // Exchange data
+        );
-    //    globalPointBoundaryFacesMap.distribute(fc);
+
-    //
+        // Exchange data
-    //    // Print
+        globalPointBoundaryFacesMap.distribute
-    //    forAll(globalPointBoundaryFaces, pointI)
+        (
-    //    {
+            transforms,
-    //        const labelList& bFaces = globalPointBoundaryFaces[pointI];
+            fc,
-    //
+            mapDistribute::transformPosition()
-    //        Pout<< "Point:" << pointI
+        );
-    //            << " at:" << coupledPatch.localPoints()[pointI]
+
-    //            << " connected to faces:" << endl;
+        // Print
-    //
+        forAll(slaves, pointI)
-    //        forAll(bFaces, i)
+        {
-    //        {
+            const labelList& slaveFaces = slaves[pointI];
-    //            Pout<< "    " << fc[bFaces[i]] << endl;
+
-    //        }
+            if (slaveFaces.size() > 0)
-    //    }
+            {
-    //}
+                Pout<< "Master point:" << pointI
-    //
+                    << " at:" << coupledPatch.localPoints()[pointI]
-    //
+                    << " connected to " << slaveFaces.size()
-    //// Test:(collocated) point to cells addressing
+                    << " untransformed faces:" << endl;
-    //{
+
-    //    const labelList& boundaryCells = globalData.boundaryCells();
+                forAll(slaveFaces, i)
-    //    const labelListList& globalPointBoundaryCells =
+                {
-    //        globalData.globalPointBoundaryCells();
+                    Pout<< "    " << fc[slaveFaces[i]] << endl;
-    //    const mapDistribute& globalPointBoundaryCellsMap =
+                }
-    //        globalData.globalPointBoundaryCellsMap();
+            }
-    //
+
-    //    pointField cc(globalPointBoundaryCellsMap.constructSize());
+            const labelList& transformedSlaveFaces = transformedSlaves[pointI];
-    //    forAll(boundaryCells, i)
+
-    //    {
+            if (transformedSlaveFaces.size() > 0)
-    //        cc[i] = mesh.cellCentres()[boundaryCells[i]];
+            {
-    //    }
+                Pout<< "Master point:" << pointI
-    //
+                    << " connected to " << transformedSlaveFaces.size()
-    //    // Exchange data
+                    << " transformed faces:" << endl;
-    //    globalPointBoundaryCellsMap.distribute(cc);
+
-    //
+                forAll(transformedSlaveFaces, i)
-    //    // Print
+                {
-    //    forAll(globalPointBoundaryCells, pointI)
+                    Pout<< "    " << fc[transformedSlaveFaces[i]] << endl;
-    //    {
+                }
-    //        const labelList& bCells = globalPointBoundaryCells[pointI];
+            }
-    //
+        }
-    //        Pout<< "Point:" << pointI
+    }
-    //            << " at:" << coupledPatch.localPoints()[pointI]
+
-    //            << " connected to cells:" << endl;
+
-    //
+    // Test: point to cells addressing
-    //        forAll(bCells, i)
+    {
-    //        {
+        const labelList& boundaryCells = globalData.boundaryCells();
-    //            Pout<< "    " << cc[bCells[i]] << endl;
+        const mapDistribute& globalPointBoundaryCellsMap =
-    //        }
+            globalData.globalPointBoundaryCellsMap();
-    //    }
+        const labelListList& slaves = globalData.globalPointBoundaryCells();
-    //}
+        const labelListList& transformedSlaves =
+            globalData.globalPointTransformedBoundaryCells();
+
+        pointField cc(globalPointBoundaryCellsMap.constructSize());
+        forAll(boundaryCells, i)
+        {
+            cc[i] = mesh.cellCentres()[boundaryCells[i]];
+        }
+
+        // Exchange data
+        globalPointBoundaryCellsMap.distribute
+        (
+            transforms,
+            cc,
+            mapDistribute::transformPosition()
+        );
+
+        // Print
+        forAll(slaves, pointI)
+        {
+            const labelList& pointCells = slaves[pointI];
+
+            if (pointCells.size() > 0)
+            {
+                Pout<< "Master point:" << pointI
+                    << " at:" << coupledPatch.localPoints()[pointI]
+                    << " connected to " << pointCells.size()
+                    << " untransformed boundaryCells:" << endl;
+
+                forAll(pointCells, i)
+                {
+                    Pout<< "    " << cc[pointCells[i]] << endl;
+                }
+            }
+
+            const labelList& transformPointCells = transformedSlaves[pointI];
+
+            if (transformPointCells.size() > 0)
+            {
+                Pout<< "Master point:" << pointI
+                    << " connected to " << transformPointCells.size()
+                    << " transformed boundaryCells:" << endl;
+
+                forAll(transformPointCells, i)
+                {
+                    Pout<< "    " << cc[transformPointCells[i]] << endl;
+                }
+            }
+        }
+    }
 
 
     Info<< "End\n" << endl;

applications/test/syncTools/Test-syncTools.C

@@ -30,12 +30,12 @@ Description
 \*---------------------------------------------------------------------------*/
 
 
+#include "syncTools.H"
 #include "argList.H"
 #include "polyMesh.H"
 #include "Time.H"
 #include "Random.H"
 #include "PackedList.H"
-#include "syncTools.H"
 
 using namespace Foam;
 
@@ -210,7 +210,7 @@ void testSparseData(const polyMesh& mesh, Random& rndGen)
     {
         // Create some data. Use slightly perturbed positions.
         Map<point> sparseData;
-        pointField fullData(mesh.nPoints(), point::max);
+        pointField fullData(mesh.nPoints(), point(GREAT, GREAT, GREAT));
 
         forAll(localPoints, i)
         {
@@ -236,7 +236,7 @@ void testSparseData(const polyMesh& mesh, Random& rndGen)
             mesh,
             fullData,
             minMagSqrEqOp<point>(),
-            point::max
+            point(GREAT, GREAT, GREAT)
             // true                    // apply separation
         );
 
@@ -246,7 +246,7 @@ void testSparseData(const polyMesh& mesh, Random& rndGen)
         {
             const point& fullPt = fullData[meshPointI];
 
-            if (fullPt != point::max)
+            if (fullPt != point(GREAT, GREAT, GREAT))
             {
                 const point& sparsePt = sparseData[meshPointI];
 
@@ -285,8 +285,8 @@ void testSparseData(const polyMesh& mesh, Random& rndGen)
 
     {
         // Create some data. Use slightly perturbed positions.
-        EdgeMap<vector> sparseData;
+        EdgeMap<point> sparseData;
-        pointField fullData(mesh.nEdges(), point::max);
+        pointField fullData(mesh.nEdges(), point(GREAT, GREAT, GREAT));
 
         const edgeList& edges = allBoundary.edges();
         const labelList meshEdges = allBoundary.meshEdges
@@ -313,14 +313,14 @@ void testSparseData(const polyMesh& mesh, Random& rndGen)
         (
             mesh,
             sparseData,
-            minEqOp<vector>()
+            minMagSqrEqOp<point>()
         );
         syncTools::syncEdgeList
         (
             mesh,
             fullData,
-            minEqOp<vector>(),
+            minMagSqrEqOp<point>(),
-            point::max
+            point(GREAT, GREAT, GREAT)
         );
 
         // Compare.
@@ -329,7 +329,7 @@ void testSparseData(const polyMesh& mesh, Random& rndGen)
         {
             const point& fullPt = fullData[meshEdgeI];
 
-            if (fullPt != point::max)
+            if (fullPt != point(GREAT, GREAT, GREAT))
             {
                 const point& sparsePt = sparseData[mesh.edges()[meshEdgeI]];
 
@@ -350,7 +350,7 @@ void testSparseData(const polyMesh& mesh, Random& rndGen)
         {
             const edge& e = mesh.edges()[meshEdgeI];
 
-            EdgeMap<vector>::const_iterator iter = sparseData.find(e);
+            EdgeMap<point>::const_iterator iter = sparseData.find(e);
 
             if (iter != sparseData.end())
             {
@@ -385,7 +385,7 @@ void testPointSync(const polyMesh& mesh, Random& rndGen)
             mesh,
             syncedPoints,
             minMagSqrEqOp<point>(),
-            point::max
+            point(GREAT, GREAT, GREAT)
         );
 
         forAll(syncedPoints, pointI)
@@ -461,7 +461,7 @@ void testEdgeSync(const polyMesh& mesh, Random& rndGen)
             mesh,
             syncedMids,
             minMagSqrEqOp<point>(),
-            point::max
+            point(GREAT, GREAT, GREAT)
         );
 
         forAll(syncedMids, edgeI)

applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C

@@ -43,6 +43,7 @@ Description
 #include "addPatchCellLayer.H"
 #include "fvMesh.H"
 #include "MeshedSurfaces.H"
+#include "globalIndex.H"
 
 #include "extrudedMesh.H"
 #include "extrudeModel.H"
@@ -261,6 +262,12 @@ int main(int argc, char *argv[])
         sourceCasePath.expand();
         fileName sourceRootDir = sourceCasePath.path();
         fileName sourceCaseDir = sourceCasePath.name();
+        if (Pstream::parRun())
+        {
+            sourceCaseDir =
+                sourceCaseDir
+               /"processor" + Foam::name(Pstream::myProcNo());
+        }
         wordList sourcePatches;
         dict.lookup("sourcePatches") >> sourcePatches;
 
@@ -279,54 +286,173 @@ int main(int argc, char *argv[])
             sourceRootDir,
             sourceCaseDir
         );
+
         #include "createMesh.H"
 
         const polyBoundaryMesh& patches = mesh.boundaryMesh();
 
 
-        // Topo change container. Either copy an existing mesh or start
-        // with empty storage (number of patches only needed for checking)
-        autoPtr<polyTopoChange> meshMod
-        (
-            (
-                mode == MESH
-              ? new polyTopoChange(mesh)
-              : new polyTopoChange(patches.size())
-            )
-        );
-
         // Extrusion engine. Either adding to existing mesh or
         // creating separate mesh.
         addPatchCellLayer layerExtrude(mesh, (mode == MESH));
 
+        const labelList meshFaces(patchFaces(patches, sourcePatches));
         indirectPrimitivePatch extrudePatch
         (
             IndirectList<face>
             (
                 mesh.faces(),
-                patchFaces(patches, sourcePatches)
+                meshFaces
             ),
             mesh.points()
         );
 
+        // Determine extrudePatch normal
+        pointField extrudePatchPointNormals
+        (
+            PatchTools::pointNormals    //calcNormals
+            (
+                mesh,
+                extrudePatch,
+                meshFaces
+            )
+        );
+
+
+        // Precalculate mesh edges for pp.edges.
+        const labelList meshEdges
+        (
+            extrudePatch.meshEdges
+            (
+                mesh.edges(),
+                mesh.pointEdges()
+            )
+        );
+
+        // Global face indices engine
+        const globalIndex globalFaces(mesh.nFaces());
+
+        // Determine extrudePatch.edgeFaces in global numbering (so across
+        // coupled patches)
+        labelListList edgeGlobalFaces
+        (
+            addPatchCellLayer::globalEdgeFaces
+            (
+                mesh,
+                globalFaces,
+                extrudePatch
+            )
+        );
+
+
+        // Determine what patches boundary edges need to get extruded into.
+        // This might actually cause edge-connected processors to become
+        // face-connected so might need to introduce new processor boundaries.
+        // Calculates:
+        //  - per pp.edge the patch to extrude into
+        //  - any additional processor boundaries (patchToNbrProc = map from
+        //    new patchID to neighbour processor)
+        //  - number of new patches (nPatches)
+
+        labelList sidePatchID;
+        label nPatches;
+        Map<label> nbrProcToPatch;
+        Map<label> patchToNbrProc;
+        addPatchCellLayer::calcSidePatch
+        (
+            mesh,
+            globalFaces,
+            edgeGlobalFaces,
+            extrudePatch,
+
+            sidePatchID,
+            nPatches,
+            nbrProcToPatch,
+            patchToNbrProc
+        );
+
+
+        // Add any patches.
+
+        label nAdded = nPatches - mesh.boundaryMesh().size();
+        reduce(nAdded, sumOp<label>());
+
+        Info<< "Adding overall " << nAdded << " processor patches." << endl;
+
+        if (nAdded > 0)
+        {
+            DynamicList<polyPatch*> newPatches(nPatches);
+            forAll(mesh.boundaryMesh(), patchI)
+            {
+                newPatches.append
+                (
+                    mesh.boundaryMesh()[patchI].clone
+                    (
+                        mesh.boundaryMesh()
+                    ).ptr()
+                );
+            }
+            for
+            (
+                label patchI = mesh.boundaryMesh().size();
+                patchI < nPatches;
+                patchI++
+            )
+            {
+                label nbrProcI = patchToNbrProc[patchI];
+                word name =
+                        "procBoundary"
+                      + Foam::name(Pstream::myProcNo())
+                      + "to"
+                      + Foam::name(nbrProcI);
+
+                Pout<< "Adding patch " << patchI
+                    << " name:" << name
+                    << " between " << Pstream::myProcNo()
+                    << " and " << nbrProcI
+                    << endl;
+
+
+                newPatches.append
+                (
+                    new processorPolyPatch
+                    (
+                        name,
+                        0,                  // size
+                        mesh.nFaces(),      // start
+                        patchI,             // index
+                        mesh.boundaryMesh(),// polyBoundaryMesh
+                        Pstream::myProcNo(),// myProcNo
+                        nbrProcI            // neighbProcNo
+                    )
+                );
+            }
+
+            // Add patches. Do no parallel updates.
+            mesh.removeFvBoundary();
+            mesh.addFvPatches(newPatches, true);
+        }
+
+
+
         // Only used for addPatchCellLayer into new mesh
-        labelList exposedPatchIDs;
+        labelList exposedPatchID;
         if (mode == PATCH)
         {
             dict.lookup("exposedPatchName") >> backPatchName;
-            exposedPatchIDs.setSize
+            exposedPatchID.setSize
             (
                 extrudePatch.size(),
                 findPatchID(patches, backPatchName)
             );
         }
 
-
+        // Determine points and extrusion
         pointField layer0Points(extrudePatch.nPoints());
         pointField displacement(extrudePatch.nPoints());
         forAll(displacement, pointI)
         {
-            const vector& patchNormal = extrudePatch.pointNormals()[pointI];
+            const vector& patchNormal = extrudePatchPointNormals[pointI];
 
             // layer0 point
             layer0Points[pointI] = model()
@@ -373,15 +499,31 @@ int main(int argc, char *argv[])
         // Layers per point
         labelList nPointLayers(extrudePatch.nPoints(), model().nLayers());
         // Displacement for first layer
-        vectorField firstLayerDisp(displacement*model().sumThickness(1));
+        vectorField firstLayerDisp = displacement*model().sumThickness(1);
+
         // Expansion ratio not used.
         scalarField ratio(extrudePatch.nPoints(), 1.0);
 
+        // Topo change container. Either copy an existing mesh or start
+        // with empty storage (number of patches only needed for checking)
+        autoPtr<polyTopoChange> meshMod
+        (
+            (
+                mode == MESH
+              ? new polyTopoChange(mesh)
+              : new polyTopoChange(patches.size())
+            )
+        );
+
         layerExtrude.setRefinement
         (
+            globalFaces,
+            edgeGlobalFaces,
+
             ratio,              // expansion ratio
             extrudePatch,       // patch faces to extrude
-            exposedPatchIDs,    // if new mesh: patches for exposed faces
+            sidePatchID,        // if boundary edge: patch to use
+            exposedPatchID,     // if new mesh: patches for exposed faces
             nFaceLayers,
             nPointLayers,
             firstLayerDisp,
@@ -401,7 +543,7 @@ int main(int argc, char *argv[])
                     model()
                     (
                         extrudePatch.localPoints()[pointI],
-                        extrudePatch.pointNormals()[pointI],
+                        extrudePatchPointNormals[pointI],
                         pPointI+1       // layer
                     )
                 );

applications/utilities/miscellaneous/foamDebugSwitches/Make/options

@@ -1,25 +1,32 @@
 EXE_LIBS = \
     -lautoMesh \
     -lbarotropicCompressibilityModel \
+    -lbasicSolidThermo \
     -lbasicThermophysicalModels \
     -lblockMesh \
     -lchemistryModel \
     -lcoalCombustion \
+    -lcombustionModels \
     -lcompressibleLESModels \
     -lcompressibleRASModels \
     -lcompressibleTurbulenceModel \
     -lconversion \
     -ldecompositionMethods \
     -ldieselSpray \
+    -ldistributed \
+    -ldistributionModels \
     -ldsmc \
     -ldynamicFvMesh \
     -ldynamicMesh \
     -ledgeMesh \
     -lengine \
     -lerrorEstimation \
+    -lEulerianInterfacialModels \
+    -lextrudeModel \
     -lfieldFunctionObjects \
     -lfileFormats \
     -lfiniteVolume \
+    -lfoamCalcFunctions \
     -lforces \
     -lfvMotionSolvers \
     -lgenericPatchFields \
@@ -29,6 +36,8 @@ EXE_LIBS = \
     -lincompressibleTurbulenceModel \
     -linterfaceProperties \
     -lIOFunctionObjects \
+    -ljobControl \
+    -lkineticTheoryModel \
     -llagrangian \
     -llagrangianIntermediate \
     -llaminarFlameSpeedModels \
@@ -37,22 +46,29 @@ EXE_LIBS = \
     -lliquidMixtureProperties \
     -lliquidProperties \
     -lmeshTools \
+    -lMGridGenGAMGAgglomeration \
     -lmolecularMeasurements \
     -lmolecule \
+    -lmultiphaseInterFoam \
     -lODE \
     -lOpenFOAM \
-    -ldistributionModels \
+    -lphaseModel \
     -lpotential \
     -lradiationModels \
     -lrandomProcesses \
     -lreactionThermophysicalModels \
+    -lreconstruct \
     -lsampling \
+    -lSLGThermo \
     -lsolidMixtureProperties \
     -lsolidParticle \
     -lsolidProperties \
+    -lsolid \
     -lspecie \
+    -lsurfaceFilmModels \
     -lsurfMesh \
     -lsystemCall \
+    -ltabulatedWallFunctions \
     -lthermalPorousZone \
     -lthermophysicalFunctions \
     -ltopoChangerFvMesh \

applications/utilities/parallelProcessing/redistributeMeshPar/redistributeMeshPar.C

@@ -73,8 +73,8 @@ autoPtr<fvMesh> createMesh
     const bool haveMesh
 )
 {
-    Pout<< "Create mesh for time = "
+    //Pout<< "Create mesh for time = "
-        << runTime.timeName() << nl << endl;
+    //    << runTime.timeName() << nl << endl;
 
     IOobject io
     (
@@ -98,12 +98,12 @@ autoPtr<fvMesh> createMesh
             xferCopy(labelList()),
             false
         );
-        Pout<< "Writing dummy mesh to " << dummyMesh.polyMesh::objectPath()
+        //Pout<< "Writing dummy mesh to " << dummyMesh.polyMesh::objectPath()
-            << endl;
+        //    << endl;
         dummyMesh.write();
     }
 
-    Pout<< "Reading mesh from " << io.objectPath() << endl;
+    //Pout<< "Reading mesh from " << io.objectPath() << endl;
     autoPtr<fvMesh> meshPtr(new fvMesh(io));
     fvMesh& mesh = meshPtr();
 
@@ -201,9 +201,8 @@ autoPtr<fvMesh> createMesh
                     break;
                 }
 
-                Pout<< "Adding patch:" << nPatches
+                //Pout<< "Adding patch:" << nPatches
-                    << " name:" << name
+                //    << " name:" << name << " type:" << type << endl;
-                    << " type:" << type << endl;
 
                 dictionary patchDict(e.dict());
                 patchDict.remove("nFaces");
@@ -282,8 +281,7 @@ autoPtr<fvMesh> createMesh
     if (!haveMesh)
     {
         // We created a dummy mesh file above. Delete it.
-        Pout<< "Removing dummy mesh " << io.objectPath()
+        //Pout<< "Removing dummy mesh " << io.objectPath() << endl;
-            << endl;
         rmDir(io.objectPath());
     }
 
@@ -351,17 +349,19 @@ scalar getMergeDistance
 }
 
 
-void printMeshData(Ostream& os, const polyMesh& mesh)
+//void printMeshData(Ostream& os, const polyMesh& mesh)
-{
+//{
-    os  << "Number of points:           " << mesh.points().size() << nl
+//    os  << "Number of points:           " << mesh.points().size() << nl
-        << "          faces:            " << mesh.faces().size() << nl
+//        << "          faces:            " << mesh.faces().size() << nl
-        << "          internal faces:   " << mesh.faceNeighbour().size() << nl
+//        << "          internal faces:   " << mesh.faceNeighbour().size() << nl
-        << "          cells:            " << mesh.cells().size() << nl
+//        << "          cells:            " << mesh.cells().size() << nl
-        << "          boundary patches: " << mesh.boundaryMesh().size() << nl
+//        << "          boundary patches: " << mesh.boundaryMesh().size() << nl
-        << "          point zones:      " << mesh.pointZones().size() << nl
+//        << "          point zones:      " << mesh.pointZones().size() << nl
-        << "          face zones:       " << mesh.faceZones().size() << nl
+//        << "          face zones:       " << mesh.faceZones().size() << nl
-        << "          cell zones:       " << mesh.cellZones().size() << nl;
+//        << "          cell zones:       " << mesh.cellZones().size() << nl;
-}
+//}
+
+
 void printMeshData(const polyMesh& mesh)
 {
     // Collect all data on master
@@ -525,6 +525,7 @@ void readFields
 
             // Receive field
             IPstream fromMaster(Pstream::blocking, Pstream::masterNo());
+            dictionary fieldDict(fromMaster);
 
             fields.set
             (
@@ -540,7 +541,7 @@ void readFields
                         IOobject::AUTO_WRITE
                     ),
                     mesh,
-                    fromMaster
+                    fieldDict
                 )
             );
 
@@ -916,8 +917,8 @@ int main(int argc, char *argv[])
     // Mesh distribution engine
     fvMeshDistribute distributor(mesh, tolDim);
 
-    Pout<< "Wanted distribution:"
+    //Pout<< "Wanted distribution:"
-        << distributor.countCells(finalDecomp) << nl << endl;
+    //    << distributor.countCells(finalDecomp) << nl << endl;
 
     // Do actual sending/receiving of mesh
     autoPtr<mapDistributePolyMesh> map = distributor.distribute(finalDecomp);
@@ -987,7 +988,7 @@ int main(int argc, char *argv[])
     Info<< endl;
 
 
-    Pout<< "End\n" << endl;
+    Info<< "End\n" << endl;
 
     return 0;
 }

applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -36,6 +36,7 @@ Description
 #include "Time.H"
 #include "surfaceFeatures.H"
 #include "featureEdgeMesh.H"
+#include "extendedFeatureEdgeMesh.H"
 #include "treeBoundBox.H"
 #include "meshTools.H"
 #include "OFstream.H"
@@ -298,20 +299,46 @@ int main(int argc, char *argv[])
 
     // Extracting and writing a featureEdgeMesh
 
-    Pout<< nl << "Writing featureEdgeMesh to constant/featureEdgeMesh."
+    extendedFeatureEdgeMesh feMesh
-        << endl;
-
-    featureEdgeMesh feMesh
     (
         newSet,
         runTime,
         surfFileName.lessExt().name() + ".featureEdgeMesh"
     );
 
+    Info<< nl << "Writing extendedFeatureEdgeMesh to " << feMesh.objectPath()
+        << endl;
+
+
     feMesh.writeObj(surfFileName.lessExt().name());
 
     feMesh.write();
 
+
+    // Write a featureEdgeMesh for backwards compatibility
+    {
+        featureEdgeMesh bfeMesh
+        (
+            IOobject
+            (
+                surfFileName.lessExt().name() + ".eMesh",   // name
+                runTime.constant(), // instance
+                "triSurface",
+                runTime,            // registry
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE,
+                false
+            ),
+            feMesh.points(),
+            feMesh.edges()
+        );
+
+        Info<< nl << "Writing featureEdgeMesh to "
+            << bfeMesh.objectPath() << endl;
+
+        bfeMesh.regIOobject::write();
+    }
+
     Info<< "End\n" << endl;
 
     return 0;

doc/codingStyleGuide.org

@@ -24,23 +24,23 @@
     + stream output
       + =<<= is always four characters after the start of the stream,
         so that the =<<= symbols align, i.e.
-        #+BEGIN_EXAMPLE
+        #+BEGIN_SRC c++
         Info<< ...
         os  << ...
-        #+END_EXAMPLE
+        #+END_SRC
         so
-        #+BEGIN_EXAMPLE
+        #+BEGIN_SRC C++
         WarningIn("className::functionName()")
             << "Warning message"
-        #+END_EXAMPLE
+        #+END_SRC
-        NOT
+        *not*
-        #+BEGIN_EXAMPLE
+        #+BEGIN_SRC C++
         WarningIn("className::functionName()")
         << "Warning message"
-        #+END_EXAMPLE
+        #+END_SRC
 
     + no unnecessary class section headers, i.e. remove
-#+BEGIN_EXAMPLE
+#+BEGIN_SRC C++
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
     // Check
@@ -48,37 +48,96 @@
     // Edit
 
     // Write
-#+END_EXAMPLE
+#+END_SRC
       if they contain nothing, even if planned for 'future use'
 
     + class titles are centred
-#+BEGIN_EXAMPLE
+#+BEGIN_SRC C++
 /*---------------------------------------------------------------------------*\
                         Class exampleClass Declaration
 \*---------------------------------------------------------------------------*/
-#+END_EXAMPLE
+#+END_SRC
 
-      NOT
+      *not*
 
-#+BEGIN_EXAMPLE
+#+BEGIN_SRC C++
 /*---------------------------------------------------------------------------*\
                 Class exampleClass Declaration
 \*---------------------------------------------------------------------------*/
-#+END_EXAMPLE
+#+END_SRC
 
 *** The /.H/ Files
     + header file spacing
       + Leave two empty lines between sections
         (as per functions in the /.C/ file etc)
-    + use =//- Comment= comments in header file
+    + use =//- Comment= comments in header file to add descriptions to class
-      + add descriptions to class data and functions
+      data and functions do be included in the Doxygen documentation:
+      + text on the line starting with =//-= becomes the Doxygen brief
+        description;
+      + text on subsequent lines becomes the Doxygen detailed description /e.g./
+        #+BEGIN_SRC C++
+        //- A function which returns a thing
+        //  This is a detailed description of the function
+        //  which processes stuff and returns other stuff
+        //  depending on things.
+        thing function(stuff1, stuff2);
+        #+END_SRC
+      + list entries start with =-= or =-#= for numbered lists but cannot start
+        on the line immediately below the brief description so
+        #+BEGIN_SRC C++
+        //- Compare triFaces
+        //  Returns:
+        //  -  0: different
+        //  - +1: identical
+        //  - -1: same face, but different orientation
+        static inline int compare(const triFace&, const triFace&);
+        #+END_SRC
+        or
+        #+BEGIN_SRC C++
+        //- Compare triFaces returning 0, +1 or -1
+        //
+        //  -  0: different
+        //  - +1: identical
+        //  - -1: same face, but different orientation
+        static inline int compare(const triFace&, const triFace&);
+        #+END_SRC
+        *not*
+        #+BEGIN_SRC C++
+        //- Compare triFaces returning 0, +1 or -1
+        //  -  0: different
+        //  - +1: identical
+        //  - -1: same face, but different orientation
+        static inline int compare(const triFace&, const triFace&);
+        #+END_SRC
+      + list can be nested for example
+        #+BEGIN_SRC C++
+        //- Search for \em name
+        //  in the following hierarchy:
+        //  -# personal settings:
+        //    - ~/.OpenFOAM/\<VERSION\>/
+        //      <em>for version-specific files</em>
+        //    - ~/.OpenFOAM/
+        //      <em>for version-independent files</em>
+        //  -# site-wide settings:
+        //    - $WM_PROJECT_INST_DIR/site/\<VERSION\>
+        //      <em>for version-specific files</em>
+        //    - $WM_PROJECT_INST_DIR/site/
+        //      <em>for version-independent files</em>
+        //  -# shipped settings:
+        //    - $WM_PROJECT_DIR/etc/
+        //
+        //  \return the full path name or fileName() if the name cannot be found
+        //  Optionally abort if the file cannot be found
+        fileName findEtcFile(const fileName&, bool mandatory=false);
+        #+END_SRC
+      + for more details see the Doxygen documentation.
     + destructor
       + If adding a comment to the destructor -
         use =//-= and code as a normal function:
-        #+BEGIN_EXAMPLE
+        #+BEGIN_SRC C++
         //- Destructor
         ~className();
-        #+END_EXAMPLE
+        #+END_SRC
     + inline functions
       + Use inline functions where appropriate in a separate /classNameI.H/
         file.  Avoid cluttering the header file with function bodies.
@@ -86,23 +145,23 @@
 *** The /.C/ Files
     + Do not open/close namespaces in a /.C/ file
       + Fully scope the function name, i.e.
-        #+BEGIN_EXAMPLE
+        #+BEGIN_SRC C++
         Foam::returnType Foam::className::functionName()
-        #+END_EXAMPLE
+        #+END_SRC
-        NOT
+        *not*
-        #+BEGIN_EXAMPLE
+        #+BEGIN_SRC C++
         namespace Foam
         {
             ...
             returnType className::functionName()
             ...
         }
-        #+END_EXAMPLE
+        #+END_SRC
         EXCEPTION
 
         When there are multiple levels of namespace, they may be used in the
         /.C/ file, i.e.
-        #+BEGIN_EXAMPLE
+        #+BEGIN_SRC C++
         namespace Foam
         {
         namespace compressible
@@ -113,7 +172,7 @@
         } // End namespace RASModels
         } // End namespace compressible
         } // End namespace Foam
-        #+END_EXAMPLE
+        #+END_SRC
 
     + Use two empty lines between functions
 
@@ -123,25 +182,25 @@
     + const
       + Use everywhere it is applicable.
     + variable initialisation using
-      #+BEGIN_EXAMPLE
+      #+BEGIN_SRC C++
       const className& variableName = otherClass.data();
-      #+END_EXAMPLE
+      #+END_SRC
-      NOT
+      *not*
-      #+BEGIN_EXAMPLE
+      #+BEGIN_SRC C++
       const className& variableName(otherClass.data());
-      #+END_EXAMPLE
+      #+END_SRC
     + virtual functions
       + If a class is virtual, make all derived classes virtual.
 
 *** Conditional Statements
-    #+BEGIN_EXAMPLE
+    #+BEGIN_SRC C++
     if (condition)
     {
         code;
     }
-    #+END_EXAMPLE
+    #+END_SRC
     OR
-    #+BEGIN_EXAMPLE
+    #+BEGIN_SRC C++
     if
     (
        long condition
@@ -149,24 +208,24 @@
     {
         code;
     }
-    #+END_EXAMPLE
+    #+END_SRC
-    NOT (no space between =if= and =(= used)
+    *not* (no space between =if= and =(= used)
-    #+BEGIN_EXAMPLE
+    #+BEGIN_SRC C++
     if(condition)
     {
         code;
     }
-    #+END_EXAMPLE
+    #+END_SRC
 
 *** =for= and =while= Loops
-    #+BEGIN_EXAMPLE
+    #+BEGIN_SRC C++
     for (i = 0; i < maxI; i++)
     {
         code;
     }
-    #+END_EXAMPLE
+    #+END_SRC
     OR
-    #+BEGIN_EXAMPLE
+    #+BEGIN_SRC C++
     for
     (
         i = 0;
@@ -176,26 +235,26 @@
     {
         code;
     }
-    #+END_EXAMPLE
+    #+END_SRC
-    NOT this (no space between =for= and =(= used)
+    *not* this (no space between =for= and =(= used)
-    #+BEGIN_EXAMPLE
+    #+BEGIN_SRC C++
     for(i = 0; i < maxI; i++)
     {
         code;
     }
-    #+END_EXAMPLE
+    #+END_SRC
     Note that when indexing through iterators, it is often slightly more
     efficient to use the pre-increment form. Eg, =++iter= instead of =iter++=
 
 *** =forAll=, =forAllIter=, =forAllConstIter=, etc. loops
     like =for= loops, but
-    #+BEGIN_EXAMPLE
+    #+BEGIN_SRC C++
     forAll(
-    #+END_EXAMPLE
+    #+END_SRC
-    NOT
+    *not*
-    #+BEGIN_EXAMPLE
+    #+BEGIN_SRC C++
     forAll (
-    #+END_EXAMPLE
+    #+END_SRC
     Using the =forAllIter= and =forAllConstIter= macros is generally
     advantageous - less typing, easier to find later.  However, since
     they are macros, they will fail if the iterated object contains
@@ -203,9 +262,9 @@
 
     The following will FAIL!:
 
-    #+BEGIN_EXAMPLE
+    #+BEGIN_SRC C++
     forAllIter(HashTable<labelPair, edge, Hash<edge> >, foo, iter)
-    #+END_EXAMPLE
+    #+END_SRC
     These convenience macros are also generally avoided in other
     container classes and OpenFOAM primitive classes.
 
@@ -214,106 +273,106 @@
       + split initially after the function return type and left align
       + do not put =const= onto its own line - use a split to keep it with
         the function name and arguments.
-        #+BEGIN_EXAMPLE
+        #+BEGIN_SRC C++
         const Foam::longReturnTypeName&
         Foam::longClassName::longFunctionName const
-        #+END_EXAMPLE
+        #+END_SRC
-        NOT
+        *not*
-        #+BEGIN_EXAMPLE
+        #+BEGIN_SRC C++
         const Foam::longReturnTypeName&
             Foam::longClassName::longFunctionName const
-        #+END_EXAMPLE
+        #+END_SRC
-        NOR
+        *nor*
-        #+BEGIN_EXAMPLE
+        #+BEGIN_SRC C++
         const Foam::longReturnTypeName& Foam::longClassName::longFunctionName
         const
-        #+END_EXAMPLE
+        #+END_SRC
-        NOR
+        *nor*
-        #+BEGIN_EXAMPLE
+        #+BEGIN_SRC C++
         const Foam::longReturnTypeName& Foam::longClassName::
         longFunctionName const
-        #+END_EXAMPLE
+        #+END_SRC
       + if it needs to be split again, split at the function name (leaving
         behind the preceding scoping =::=s), and again, left align, i.e.
-        #+BEGIN_EXAMPLE
+        #+BEGIN_SRC C++
         const Foam::longReturnTypeName&
         Foam::veryveryveryverylongClassName::
         veryveryveryverylongFunctionName const
-        #+END_EXAMPLE
+        #+END_SRC
 
 ***** Splitting long lines at an "="
      Indent after split
-     #+BEGIN_EXAMPLE
+     #+BEGIN_SRC C++
      variableName =
          longClassName.longFunctionName(longArgument);
-     #+END_EXAMPLE
+     #+END_SRC
      OR (where necessary)
-     #+BEGIN_EXAMPLE
+     #+BEGIN_SRC C++
      variableName =
          longClassName.longFunctionName
          (
              longArgument1,
              longArgument2
          );
-     #+END_EXAMPLE
+     #+END_SRC
-     NOT
+     *not*
-     #+BEGIN_EXAMPLE
+     #+BEGIN_SRC C++
      variableName =
      longClassName.longFunctionName(longArgument);
-     #+END_EXAMPLE
+     #+END_SRC
-     NOR
+     *nor*
-     #+BEGIN_EXAMPLE
+     #+BEGIN_SRC C++
      variableName = longClassName.longFunctionName
      (
          longArgument1,
          longArgument2
      );
-     #+END_EXAMPLE
+     #+END_SRC
 
 *** Maths and Logic
     + operator spacing
-      #+BEGIN_EXAMPLE
+      #+BEGIN_SRC C++
       a + b, a - b
       a*b, a/b
       a & b, a ^ b
       a = b, a != b
       a < b, a > b, a >= b, a <= b
       a || b, a && b
-      #+END_EXAMPLE
+      #+END_SRC
 
     + splitting formulae over several lines
 
       Split and indent as per "splitting long lines at an ="
       with the operator on the lower line.  Align operator so that first
       variable, function or bracket on the next line is 4 spaces indented i.e.
-      #+BEGIN_EXAMPLE
+      #+BEGIN_SRC C++
       variableName =
           a*(a + b)
          *exp(c/d)
          *(k + t);
-      #+END_EXAMPLE
+      #+END_SRC
       This is sometimes more legible when surrounded by extra parentheses:
 
-      #+BEGIN_EXAMPLE
+      #+BEGIN_SRC C++
       variableName =
       (
           a*(a + b)
          *exp(c/d)
          *(k + t)
       );
-      #+END_EXAMPLE
+      #+END_SRC
 
     + splitting logical tests over several lines
 
       outdent the operator so that the next variable to test is aligned with
       the four space indentation, i.e.
-      #+BEGIN_EXAMPLE
+      #+BEGIN_SRC C++
       if
       (
           a == true
        && b == c
       )
-      #+END_EXAMPLE
+      #+END_SRC
 
 ** Documentation
 *** General
@@ -331,7 +390,7 @@
       description.
 
       For example,
-      #+BEGIN_EXAMPLE
+      #+BEGIN_SRC C++
       Class
           Foam::myClass
 
@@ -340,7 +399,7 @@
 
           The class is implemented as a set of recommendations that may
           sometimes be useful.
-      #+END_EXAMPLE
+      #+END_SRC
 
     + The class name must be qualified by its namespace, otherwise Doxygen
       will think you are documenting some other class.
@@ -348,33 +407,33 @@
     + If you don't have anything to say about the class (at the moment), use
       the namespace-qualified class name for the description. This aids with
       finding these under-documented classes later.
-      #+BEGIN_EXAMPLE
+      #+BEGIN_SRC C++
       Class
           Foam::myUnderDocumentedClass
 
       Description
           Foam::myUnderDocumentedClass
-      #+END_EXAMPLE
+      #+END_SRC
 
     + Use 'Class' and 'Namespace' tags in the header files.
       The Description block then applies to documenting the class.
 
     + Use 'InClass' and 'InNamespace' in the source files.
       The Description block then applies to documenting the file itself.
-      #+BEGIN_EXAMPLE
+      #+BEGIN_SRC C++
       InClass
           Foam::myClass
 
       Description
           Implements the read and writing of files.
-      #+END_EXAMPLE
+      #+END_SRC
 
 *** Doxygen Special Commands
     Doxygen has a large number of special commands with a =\= prefix.
 
     Since the filtering removes the leading spaces within the blocks, the
     Doxygen commmands can be inserted within the block without problems.
-    #+BEGIN_EXAMPLE
+    #+BEGIN_SRC C++
     InClass
         Foam::myClass
 
@@ -398,7 +457,7 @@
                 ...  // some operation
             }
         \endcode
-    #+END_EXAMPLE
+    #+END_SRC
 
 *** HTML Special Commands
     Since Doxygen also handles HTML tags to a certain extent, the angle
@@ -418,16 +477,16 @@
     + If the namespaces is used to hold sub-models, the namespace can be
       documented in the same file as the class with the model selector.
       eg,
-      #+BEGIN_EXAMPLE
+      #+BEGIN_SRC C++
       documented namespace 'Foam::functionEntries' within the
       class 'Foam::functionEntry'
-      #+END_EXAMPLE
+      #+END_SRC
 
     + If nothing else helps, find some sensible header.
       eg,
-      #+BEGIN_EXAMPLE
+      #+BEGIN_SRC C++
       namespace 'Foam' is documented in the foamVersion.H file
-      #+END_EXAMPLE
+      #+END_SRC
 
 *** Documenting typedefs and classes defined via macros
     ... not yet properly resolved
@@ -456,9 +515,9 @@
     used as a variable or class method name, it is probably better to use
     '-ize', which is considered the main form by the Oxford University
     Press. Eg,
-    #+BEGIN_EXAMPLE
+    #+BEGIN_SRC C++
     myClass.initialize()
-    #+END_EXAMPLE
+    #+END_SRC
 
     The word "its" (possesive) vs. "it's" (colloquial for "it is" or "it has")
     seems to confuse non-native (and some native) English speakers.

src/OpenFOAM/fields/Fields/transformList/transformList.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -46,11 +46,17 @@ Foam::List<T> Foam::transform
 
 
 template <class T>
-void Foam::transformList
+void Foam::transformList(const tensor& rotTensor, UList<T>& field)
-(
+{
-    const tensorField& rotTensor,
+    forAll(field, i)
-    UList<T>& field
+    {
-)
+        field[i] = transform(rotTensor, field[i]);
+    }
+}
+
+
+template <class T>
+void Foam::transformList(const tensorField& rotTensor, UList<T>& field)
 {
     if (rotTensor.size() == 1)
     {
@@ -79,11 +85,17 @@ void Foam::transformList
 
 
 template <class T>
-void Foam::transformList
+void Foam::transformList(const tensor& rotTensor, Map<T>& field)
-(
+{
-    const tensorField& rotTensor,
+    forAllIter(typename Map<T>, field, iter)
-    Map<T>& field
+    {
-)
+        iter() = transform(rotTensor[0], iter());
+    }
+}
+
+
+template <class T>
+void Foam::transformList(const tensorField& rotTensor, Map<T>& field)
 {
     if (rotTensor.size() == 1)
     {
@@ -105,11 +117,17 @@ void Foam::transformList
 
 
 template <class T>
-void Foam::transformList
+void Foam::transformList(const tensor& rotTensor, EdgeMap<T>& field)
-(
+{
-    const tensorField& rotTensor,
+    forAllIter(typename EdgeMap<T>, field, iter)
-    EdgeMap<T>& field
+    {
-)
+        iter() = transform(rotTensor[0], iter());
+    }
+}
+
+
+template <class T>
+void Foam::transformList(const tensorField& rotTensor, EdgeMap<T>& field)
 {
     if (rotTensor.size() == 1)
     {

src/OpenFOAM/fields/Fields/transformList/transformList.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -60,33 +60,57 @@ List<T> transform
 //- Apply transformation to list. Either single transformation tensor
 //  or one tensor per element.
 template<class T>
+void transformList(const tensor&, UList<T>&);
+template<class T>
 void transformList(const tensorField&, UList<T>&);
 
+template<class T>
+void transformList(const tensor&, Map<T>&);
 template<class T>
 void transformList(const tensorField&, Map<T>&);
 
+template<class T>
+void transformList(const tensor&, EdgeMap<T>&);
 template<class T>
 void transformList(const tensorField&, EdgeMap<T>&);
 
 
+template<>
+inline void transformList(const tensor&, labelUList&)
+{}
 template<>
 inline void transformList(const tensorField&, labelUList&)
 {}
 template<>
+inline void transformList(const tensor&, Map<label>&)
+{}
+template<>
 inline void transformList(const tensorField&, Map<label>&)
 {}
 template<>
+inline void transformList(const tensor&, EdgeMap<label>&)
+{}
+template<>
 inline void transformList(const tensorField&, EdgeMap<label>&)
 {}
 
 
+template<>
+inline void transformList(const tensor&, UList<scalar>&)
+{}
 template<>
 inline void transformList(const tensorField&, UList<scalar>&)
 {}
 template<>
+inline void transformList(const tensor&, Map<scalar>&)
+{}
+template<>
 inline void transformList(const tensorField&, Map<scalar>&)
 {}
 template<>
+inline void transformList(const tensor&, EdgeMap<scalar>&)
+{}
+template<>
 inline void transformList(const tensorField&, EdgeMap<scalar>&)
 {}
 

src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H

@@ -194,18 +194,22 @@ class globalMeshData
             mutable autoPtr<mapDistribute> globalEdgeSlavesMapPtr_;
 
 
-            //// Coupled point to boundary faces
+            // Coupled point to boundary faces
-            //
+
-            //mutable autoPtr<globalIndex> globalBoundaryFaceNumberingPtr_;
+            mutable autoPtr<globalIndex> globalBoundaryFaceNumberingPtr_;
-            //mutable autoPtr<labelListList> globalPointBoundaryFacesPtr_;
+            mutable autoPtr<labelListList> globalPointBoundaryFacesPtr_;
-            //mutable autoPtr<mapDistribute> globalPointBoundaryFacesMapPtr_;
+            mutable autoPtr<labelListList>
-            //
+                globalPointTransformedBoundaryFacesPtr_;
-            //// Coupled point to collocated boundary cells
+            mutable autoPtr<mapDistribute> globalPointBoundaryFacesMapPtr_;
-            //
+
-            //mutable autoPtr<labelList> boundaryCellsPtr_;
+            // Coupled point to collocated boundary cells
-            //mutable autoPtr<globalIndex> globalBoundaryCellNumberingPtr_;
+
-            //mutable autoPtr<labelListList> globalPointBoundaryCellsPtr_;
+            mutable autoPtr<labelList> boundaryCellsPtr_;
-            //mutable autoPtr<mapDistribute> globalPointBoundaryCellsMapPtr_;
+            mutable autoPtr<globalIndex> globalBoundaryCellNumberingPtr_;
+            mutable autoPtr<labelListList> globalPointBoundaryCellsPtr_;
+            mutable autoPtr<labelListList>
+                globalPointTransformedBoundaryCellsPtr_;
+            mutable autoPtr<mapDistribute> globalPointBoundaryCellsMapPtr_;
 
 
         // Globally shared point addressing
@@ -288,6 +292,18 @@ class globalMeshData
             void calcGlobalEdgeSlaves() const;
 
 
+        // Global boundary face/cell addressing
+
+            //- Calculate coupled point to uncoupled boundary faces. Local only.
+            void calcPointBoundaryFaces(labelListList&) const;
+
+            //- Calculate global point to global boundary face addressing.
+            void calcGlobalPointBoundaryFaces() const;
+
+            //- Calculate global point to global boundary cell addressing.
+            void calcGlobalPointBoundaryCells() const;
+
+
         //- Disallow default bitwise copy construct
         globalMeshData(const globalMeshData&);
 
@@ -458,7 +474,7 @@ public:
             const globalIndexAndTransform& globalTransforms() const;
 
             //- Helper: synchronise data with transforms
-            template<class Type, class CombineOp>
+            template<class Type, class CombineOp, class TransformOp>
             static void syncData
             (
                 List<Type>& pointData,
@@ -467,7 +483,7 @@ public:
                 const mapDistribute& slavesMap,
                 const globalIndexAndTransform&,
                 const CombineOp& cop,
-                const bool isPosition
+                const TransformOp& top
             );
 
             //- Helper: synchronise data without transforms
@@ -490,13 +506,13 @@ public:
                 const labelListList& globalPointSlaves() const;
                 const labelListList& globalPointTransformedSlaves() const;
                 const mapDistribute& globalPointSlavesMap() const;
-                //- Helper to synchronise mesh point data
+                //- Helper to synchronise coupled patch point data
-                template<class Type, class CombineOp>
+                template<class Type, class CombineOp, class TransformOp>
                 void syncPointData
                 (
                     List<Type>& pointData,
                     const CombineOp& cop,
-                    const bool isPosition
+                    const TransformOp& top
                 ) const;
 
             // Coupled edge to coupled edges.
@@ -506,6 +522,28 @@ public:
                 const labelListList& globalEdgeTransformedSlaves() const;
                 const mapDistribute& globalEdgeSlavesMap() const;
 
+            // Coupled point to boundary faces. These are uncoupled boundary
+            // faces only but include empty patches.
+
+                //- Numbering of boundary faces is face-mesh.nInternalFaces()
+                const globalIndex& globalBoundaryFaceNumbering() const;
+                const labelListList& globalPointBoundaryFaces() const;
+                const labelListList& globalPointTransformedBoundaryFaces()
+                const;
+                const mapDistribute& globalPointBoundaryFacesMap() const;
+
+            // Coupled point to boundary cell
+
+                //- From boundary cell to mesh cell
+                const labelList& boundaryCells() const;
+
+                //- Numbering of boundary cells is according to boundaryCells()
+                const globalIndex& globalBoundaryCellNumbering() const;
+                const labelListList& globalPointBoundaryCells() const;
+                const labelListList& globalPointTransformedBoundaryCells()
+                const;
+                const mapDistribute& globalPointBoundaryCellsMap() const;
+
 
             // Other
 

src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshDataTemplates.C

@@ -29,7 +29,7 @@ License
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class Type, class CombineOp>
+template<class Type, class CombineOp, class TransformOp>
 void Foam::globalMeshData::syncData
 (
     List<Type>& elems,
@@ -38,11 +38,11 @@ void Foam::globalMeshData::syncData
     const mapDistribute& slavesMap,
     const globalIndexAndTransform& transforms,
     const CombineOp& cop,
-    const bool isPosition
+    const TransformOp& top
 )
 {
     // Pull slave data onto master
-    slavesMap.distribute(transforms, elems, isPosition);
+    slavesMap.distribute(transforms, elems, top);
 
     // Combine master data with slave data
     forAll(slaves, i)
@@ -85,7 +85,7 @@ void Foam::globalMeshData::syncData
         transforms,
         elems.size(),
         elems,
-        isPosition
+        top
     );
 }
 
@@ -143,12 +143,12 @@ void Foam::globalMeshData::syncData
 }
 
 
-template<class Type, class CombineOp>
+template<class Type, class CombineOp, class TransformOp>
 void Foam::globalMeshData::syncPointData
 (
     List<Type>& pointData,
     const CombineOp& cop,
-    const bool isPosition
+    const TransformOp& top
 ) const
 {
     if (pointData.size() != mesh_.nPoints())
@@ -171,7 +171,7 @@ void Foam::globalMeshData::syncPointData
         globalPointSlavesMap(),
         globalTransforms(),
         cop,
-        isPosition
+        top
     );
 
     // Extract back onto mesh

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

@@ -37,6 +37,99 @@ defineTypeNameAndDebug(Foam::mapDistribute, 0);
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+template<>
+void Foam::mapDistribute::transform::operator()
+(
+    const vectorTensorTransform&,
+    const bool,
+    List<label>&
+) const
+{}
+template<>
+void Foam::mapDistribute::transform::operator()
+(
+    const coupledPolyPatch&,
+    UList<label>&
+) const
+{}
+template<>
+void Foam::mapDistribute::transform::operator()
+(
+    const coupledPolyPatch&,
+    Map<label>&
+) const
+{}
+template<>
+void Foam::mapDistribute::transform::operator()
+(
+    const coupledPolyPatch&,
+    EdgeMap<label>&
+) const
+{}
+
+
+template<>
+void Foam::mapDistribute::transform::operator()
+(
+    const vectorTensorTransform&,
+    const bool,
+    List<scalar>&
+) const
+{}
+template<>
+void Foam::mapDistribute::transform::operator()
+(
+    const coupledPolyPatch&,
+    UList<scalar>&
+) const
+{}
+template<>
+void Foam::mapDistribute::transform::operator()
+(
+    const coupledPolyPatch&,
+    Map<scalar>&
+) const
+{}
+template<>
+void Foam::mapDistribute::transform::operator()
+(
+    const coupledPolyPatch&,
+    EdgeMap<scalar>&
+) const
+{}
+
+
+template<>
+void Foam::mapDistribute::transform::operator()
+(
+    const vectorTensorTransform&,
+    const bool,
+    List<bool>&
+) const
+{}
+template<>
+void Foam::mapDistribute::transform::operator()
+(
+    const coupledPolyPatch&,
+    UList<bool>&
+) const
+{}
+template<>
+void Foam::mapDistribute::transform::operator()
+(
+    const coupledPolyPatch&,
+    Map<bool>&
+) const
+{}
+template<>
+void Foam::mapDistribute::transform::operator()
+(
+    const coupledPolyPatch&,
+    EdgeMap<bool>&
+) const
+{}
+
+
 Foam::List<Foam::labelPair> Foam::mapDistribute::schedule
 (
     const labelListList& subMap,
@@ -535,101 +628,6 @@ void Foam::mapDistribute::exchangeAddressing
 }
 
 
-template<>
-void Foam::mapDistribute::applyTransforms
-(
-    const globalIndexAndTransform& globalTransforms,
-    List<point>& field,
-    const bool isPosition
-) const
-{
-    const List<vectorTensorTransform>& totalTransform =
-        globalTransforms.transformPermutations();
-
-    forAll(totalTransform, trafoI)
-    {
-        const vectorTensorTransform& vt = totalTransform[trafoI];
-        const labelList& elems = transformElements_[trafoI];
-        label n = transformStart_[trafoI];
-
-        // Could be optimised to avoid memory allocations
-
-        if (isPosition)
-        {
-            Field<point> transformFld
-            (
-                vt.transformPosition(Field<point>(field, elems))
-            );
-            forAll(transformFld, i)
-            {
-                //cop(field[n++], transformFld[i]);
-                field[n++] = transformFld[i];
-            }
-        }
-        else
-        {
-            Field<point> transformFld
-            (
-                transform(vt.R(), Field<point>(field, elems))
-            );
-
-            forAll(transformFld, i)
-            {
-                //cop(field[n++], transformFld[i]);
-                field[n++] = transformFld[i];
-            }
-        }
-    }
-}
-
-
-template<>
-void Foam::mapDistribute::applyInverseTransforms
-(
-    const globalIndexAndTransform& globalTransforms,
-    List<point>& field,
-    const bool isPosition
-) const
-{
-    const List<vectorTensorTransform>& totalTransform =
-        globalTransforms.transformPermutations();
-
-    forAll(totalTransform, trafoI)
-    {
-        const vectorTensorTransform& vt = totalTransform[trafoI];
-        const labelList& elems = transformElements_[trafoI];
-        label n = transformStart_[trafoI];
-
-        // Could be optimised to avoid memory allocations
-
-        if (isPosition)
-        {
-            Field<point> transformFld
-            (
-                vt.invTransformPosition
-                (
-                    SubField<point>(field, elems.size(), n)
-                )
-            );
-            forAll(transformFld, i)
-            {
-                field[elems[i]] = transformFld[i];
-            }
-        }
-        else
-        {
-            Field<point> transformFld(SubField<point>(field, elems.size(), n));
-            transform(transformFld, vt.R().T(), transformFld);
-
-            forAll(transformFld, i)
-            {
-                field[elems[i]] = transformFld[i];
-            }
-        }
-    }
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 //- Construct null

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

@@ -127,12 +127,14 @@ SourceFiles
 #ifndef mapDistribute_H
 #define mapDistribute_H
 
+#include "transformList.H"
 #include "labelList.H"
 #include "labelPair.H"
 #include "Pstream.H"
 #include "boolList.H"
 #include "Map.H"
-#include "point.H"
+#include "vectorTensorTransform.H"
+#include "coupledPolyPatch.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -217,31 +219,124 @@ class mapDistribute
         template<class T>
         void applyDummyTransforms(List<T>& field) const;
 
-        template<class T>   //, class CombineOp>
+        template<class T, class TransformOp>
         void applyTransforms
         (
             const globalIndexAndTransform& globalTransforms,
             List<T>& field,
-            const bool isPosition
+            const TransformOp& top
-            //const CombineOp& cop
         ) const;
 
         //- Helper function: copy transformElements without transformation
         template<class T>
         void applyDummyInverseTransforms(List<T>& field) const;
 
-        template<class T>   //, class CombineOp>
+        template<class T, class TransformOp>
         void applyInverseTransforms
         (
             const globalIndexAndTransform& globalTransforms,
             List<T>& field,
-            const bool isPosition
+            const TransformOp& top
-            //const CombineOp& cop
         ) const;
 
 
 public:
 
+    // Public classes
+
+        //- Default transformation behaviour
+        class transform
+        {
+        public:
+
+            template<class T>
+            void operator()
+            (
+                const vectorTensorTransform& vt,
+                const bool forward,
+                List<T>& fld
+            ) const
+            {
+                if (forward)
+                {
+                    transformList(vt.R(), fld);
+                }
+                else
+                {
+                    transformList(vt.R().T(), fld);
+                }
+            }
+
+            //- Transform patch-based field
+            template<class T>
+            void operator()(const coupledPolyPatch& cpp, UList<T>& fld) const
+            {
+                if (!cpp.parallel())
+                {
+                    transformList(cpp.forwardT(), fld);
+                }
+            }
+
+            //- Transform sparse field
+            template<class T, template<class> class Container>
+            void operator()(const coupledPolyPatch& cpp, Container<T>& map)
+            const
+            {
+                if (!cpp.parallel())
+                {
+                    transformList(cpp.forwardT(), map);
+                }
+            }
+        };
+
+        //- Default transformation behaviour for position
+        class transformPosition
+        {
+        public:
+
+            void operator()
+            (
+                const vectorTensorTransform& vt,
+                const bool forward,
+                List<point>& fld
+            ) const
+            {
+                pointField pfld(fld.xfer());
+                if (forward)
+                {
+                    fld = vt.transformPosition(pfld);
+                }
+                else
+                {
+                    fld = vt.invTransformPosition(pfld);
+                }
+            }
+            //- Transform patch-based field
+            void operator()(const coupledPolyPatch& cpp, pointField& fld) const
+            {
+                cpp.transformPosition(fld);
+            }
+            template<template<class> class Container>
+            void operator()(const coupledPolyPatch& cpp, Container<point>& map)
+            const
+            {
+                Field<point> fld(map.size());
+                label i = 0;
+                forAllConstIter(typename Container<point>, map, iter)
+                {
+                    fld[i++] = iter();
+                }
+                cpp.transformPosition(fld);
+                i = 0;
+                forAllIter(typename Container<point>, map, iter)
+                {
+                    iter() = fld[i++];
+                }
+            }
+        };
+
+
+
     // Declare name of the class and its debug switch
     ClassName("mapDistribute");
 
@@ -456,12 +551,12 @@ public:
             const;
 
             //- Same but with transforms
-            template<class T>
+            template<class T, class TransformOp>
             void distribute
             (
                 const globalIndexAndTransform&,
                 List<T>& fld,
-                const bool isPosition
+                const TransformOp& top
             ) const;
 
             //- Reverse distribute data using default commsType.
@@ -474,13 +569,13 @@ public:
             ) const;
 
             //- Same but with transforms
-            template<class T>
+            template<class T, class TransformOp>
             void reverseDistribute
             (
                 const globalIndexAndTransform&,
                 const label constructSize,
                 List<T>& fld,
-                const bool isPosition
+                const TransformOp& top
             ) const;
 
             //- Reverse distribute data using default commsType.
@@ -496,14 +591,14 @@ public:
             ) const;
 
             //- Same but with transforms
-            template<class T>
+            template<class T, class TransformOp>
             void reverseDistribute
             (
                 const globalIndexAndTransform&,
                 const label constructSize,
                 const T& nullValue,
                 List<T>& fld,
-                const bool isPosition
+                const TransformOp& top
             ) const;
 
             //- Do all sends using PstreamBuffers
@@ -531,24 +626,85 @@ public:
 
 };
 
-//- Specialisation for transforms that can apply positional transform
+
 template<>
-void mapDistribute::applyTransforms
+void mapDistribute::transform::operator()
 (
-    const globalIndexAndTransform& globalTransforms,
+    const vectorTensorTransform&,
-    List<point>& field,
+    const bool,
-    const bool isPosition
+    List<label>&
-    //const CombineOp& cop
+) const;
+template<>
+void mapDistribute::transform::operator()
+(
+    const coupledPolyPatch&,
+    UList<label>&
+) const;
+template<>
+void mapDistribute::transform::operator()
+(
+    const coupledPolyPatch&,
+    Map<label>&
+) const;
+template<>
+void mapDistribute::transform::operator()
+(
+    const coupledPolyPatch&,
+    EdgeMap<label>&
 ) const;
 
-template<>   //, class CombineOp>
+template<>
-void mapDistribute::applyInverseTransforms
+void mapDistribute::transform::operator()
 (
-    const globalIndexAndTransform& globalTransforms,
+    const coupledPolyPatch&,
-    List<point>& field,
+    UList<scalar>&
-    const bool isPosition
-    //const CombineOp& cop
 ) const;
+template<>
+void mapDistribute::transform::operator()
+(
+    const vectorTensorTransform&,
+    const bool,
+    List<scalar>&
+) const;
+template<>
+void mapDistribute::transform::operator()
+(
+    const coupledPolyPatch&,
+    Map<scalar>&
+) const;
+template<>
+void mapDistribute::transform::operator()
+(
+    const coupledPolyPatch&,
+    EdgeMap<scalar>&
+) const;
+
+template<>
+void mapDistribute::transform::operator()
+(
+    const coupledPolyPatch& cpp,
+    UList<bool>& fld
+) const;
+template<>
+void mapDistribute::transform::operator()
+(
+    const vectorTensorTransform&,
+    const bool,
+    List<bool>&
+) const;
+template<>
+void mapDistribute::transform::operator()
+(
+    const coupledPolyPatch&,
+    Map<bool>&
+) const;
+template<>
+void mapDistribute::transform::operator()
+(
+    const coupledPolyPatch&,
+    EdgeMap<bool>&
+) const;
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

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

@@ -812,30 +812,14 @@ void Foam::mapDistribute::applyDummyInverseTransforms(List<T>& field) const
 
 
 // Calculate transformed elements.
-template<class T>   //, class CombineOp>
+template<class T, class TransformOp>   //, class CombineOp>
 void Foam::mapDistribute::applyTransforms
 (
     const globalIndexAndTransform& globalTransforms,
     List<T>& field,
-    const bool isPosition
+    const TransformOp& top
-    //const CombineOp& cop
 ) const
 {
-    if (isPosition)
-    {
-        FatalErrorIn
-        (
-            "mapDistribute::applyTransforms\n"
-            "(\n"
-            "    const globalIndexAndTransform&,\n"
-            "    List<T>&,\n"
-            "    const bool\n"
-            ") const\n"
-        )   << "It does not make sense to apply position transformation"
-            << " for anything else than pointFields."
-            << abort(FatalError);
-    }
-
     const List<vectorTensorTransform>& totalTransform =
         globalTransforms.transformPermutations();
 
@@ -846,7 +830,8 @@ void Foam::mapDistribute::applyTransforms
         label n = transformStart_[trafoI];
 
         // Could be optimised to avoid memory allocations
-        Field<T> transformFld(transform(vt.R(), Field<T>(field, elems)));
+        List<T> transformFld(UIndirectList<T>(field, elems));
+        top(vt, true, transformFld);
 
         forAll(transformFld, i)
         {
@@ -858,30 +843,14 @@ void Foam::mapDistribute::applyTransforms
 
 
 // Calculate transformed elements.
-template<class T>   //, class CombineOp>
+template<class T, class TransformOp>   //, class CombineOp>
 void Foam::mapDistribute::applyInverseTransforms
 (
     const globalIndexAndTransform& globalTransforms,
     List<T>& field,
-    const bool isPosition
+    const TransformOp& top
-    //const CombineOp& cop
 ) const
 {
-    if (isPosition)
-    {
-        FatalErrorIn
-        (
-            "mapDistribute::applyInverseTransforms\n"
-            "(\n"
-            "    const globalIndexAndTransform&,\n"
-            "    List<T>&,\n"
-            "    const bool\n"
-            ") const\n"
-        )   << "It does not make sense to apply position transformation"
-            << " for anything else than pointFields."
-            << abort(FatalError);
-    }
-
     const List<vectorTensorTransform>& totalTransform =
         globalTransforms.transformPermutations();
 
@@ -892,8 +861,8 @@ void Foam::mapDistribute::applyInverseTransforms
         label n = transformStart_[trafoI];
 
         // Could be optimised to avoid memory allocations
-        Field<T> transformFld(SubField<T>(field, elems.size(), n));
+        List<T> transformFld(SubList<T>(field, elems.size(), n));
-        transform(transformFld, vt.R().T(), transformFld);
+        top(vt, false, transformFld);
 
         forAll(transformFld, i)
         {
@@ -1073,54 +1042,54 @@ void Foam::mapDistribute::reverseDistribute
 
 
 //- Distribute data using default commsType.
-template<class T>
+template<class T, class TransformOp>
 void Foam::mapDistribute::distribute
 (
     const globalIndexAndTransform& git,
     List<T>& fld,
-    const bool isPosition
+    const TransformOp& top
 ) const
 {
     // Distribute. Leave out dummy transforms since we're doing them ourselves
     distribute(fld, false);
     // Do transforms
-    applyTransforms(git, fld, isPosition);   //, eqOp<T>());
+    applyTransforms(git, fld, top);
 }
 
 
-template<class T>
+template<class T, class TransformOp>
 void Foam::mapDistribute::reverseDistribute
 (
     const globalIndexAndTransform& git,
     const label constructSize,
     List<T>& fld,
-    const bool isPosition
+    const TransformOp& top
 ) const
 {
     // Fill slots with reverse-transformed data. Note that it also copies
     // back into the non-remote part of fld even though these values are not
     // used.
-    applyInverseTransforms(git, fld, isPosition);   //, eqOp<T>());
+    applyInverseTransforms(git, fld, top);
 
     // And send back (the remote slots). Disable dummy transformations.
     reverseDistribute(constructSize, fld, false);
 }
 
 
-template<class T>
+template<class T, class TransformOp>
 void Foam::mapDistribute::reverseDistribute
 (
     const globalIndexAndTransform& git,
     const label constructSize,
     const T& nullValue,
     List<T>& fld,
-    const bool isPosition
+    const TransformOp& top
 ) const
 {
     // Fill slots with reverse-transformed data Note that it also copies
     // back into the non-remote part of fld even though these values are not
     // used.
-    applyInverseTransforms(git, fld, isPosition);   //, eqOp<T>());
+    applyInverseTransforms(git, fld, top);   //, eqOp<T>());
 
     // And send back (the remote slots) Disable dummy transformations.
     reverseDistribute(constructSize, nullValue, fld, false);

src/OpenFOAM/meshes/polyMesh/syncTools/syncTools.C

@@ -25,77 +25,6 @@ License
 
 #include "syncTools.H"
 
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-template<>
-void Foam::syncTools::transform::operator()
-(
-    const coupledPolyPatch&,
-    Field<label>&
-) const
-{}
-template<>
-void Foam::syncTools::transform::operator()
-(
-    const coupledPolyPatch&,
-    Map<label>&
-) const
-{}
-template<>
-void Foam::syncTools::transform::operator()
-(
-    const coupledPolyPatch&,
-    EdgeMap<label>&
-) const
-{}
-
-
-template<>
-void Foam::syncTools::transform::operator()
-(
-    const coupledPolyPatch&,
-    Field<scalar>&
-) const
-{}
-template<>
-void Foam::syncTools::transform::operator()
-(
-    const coupledPolyPatch&,
-    Map<scalar>&
-) const
-{}
-template<>
-void Foam::syncTools::transform::operator()
-(
-    const coupledPolyPatch&,
-    EdgeMap<scalar>&
-) const
-{}
-
-
-template<>
-void Foam::syncTools::transform::operator()
-(
-    const coupledPolyPatch&,
-    Field<bool>&
-) const
-{}
-template<>
-void Foam::syncTools::transform::operator()
-(
-    const coupledPolyPatch&,
-    Map<bool>&
-) const
-{}
-template<>
-void Foam::syncTools::transform::operator()
-(
-    const coupledPolyPatch&,
-    EdgeMap<bool>&
-) const
-{}
-
-
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 // Determines for every point whether it is coupled and if so sets only one.

src/OpenFOAM/meshes/polyMesh/syncTools/syncTools.H

@@ -47,7 +47,7 @@ SourceFiles
 #include "PackedBoolList.H"
 #include "polyMesh.H"
 #include "coupledPolyPatch.H"
-#include "transformList.H"
+#include "mapDistribute.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -87,59 +87,6 @@ class syncTools
 
 public:
 
-    // Public classes
-
-        class transform
-        {
-        public:
-            //- Transform patch-based field
-            template<class T>
-            void operator()(const coupledPolyPatch& cpp, Field<T>& fld) const
-            {
-                if (!cpp.parallel())
-                {
-                    transformList(cpp.forwardT(), fld);
-                }
-            }
-
-            //- Transform sparse field
-            template<class T, template<class> class Container>
-            void operator()(const coupledPolyPatch& cpp, Container<T>& map)
-            const
-            {
-                if (!cpp.parallel())
-                {
-                    transformList(cpp.forwardT(), map);
-                }
-            }
-        };
-
-        class transformPosition
-        {
-        public:
-            void operator()(const coupledPolyPatch& cpp, pointField& fld) const
-            {
-                cpp.transformPosition(fld);
-            }
-            template<template<class> class Container>
-            void operator()(const coupledPolyPatch& cpp, Container<point>& map)
-            const
-            {
-                Field<point> fld(map.size());
-                label i = 0;
-                forAllConstIter(typename Container<point>, map, iter)
-                {
-                    fld[i++] = iter();
-                }
-                cpp.transformPosition(fld);
-                i = 0;
-                forAllIter(typename Container<point>, map, iter)
-                {
-                    iter() = fld[i++];
-                }
-            }
-        };
-
 
         // Basic routines with user-supplied transformation. Preferably
         // use specialisations below.
@@ -163,41 +110,41 @@ public:
                 const CombineOp& cop,
                 const TransformOp& top
             );
-//
+
-//            //- Synchronize values on all mesh points.
+            //- Synchronize values on all mesh points.
-//            template <class T, class CombineOp, class TransformOp>
+            template <class T, class CombineOp, class TransformOp>
-//            static void syncPointList
+            static void syncPointList
-//            (
+            (
-//                const polyMesh&,
+                const polyMesh&,
-//                UList<T>&,
+                List<T>&,
-//                const CombineOp& cop,
+                const CombineOp& cop,
-//                const T& nullValue,
+                const T& nullValue,
-//                const TransformOp& top
+                const TransformOp& top
-//            );
+            );
-//
+
-//            //- Synchronize values on selected mesh points.
+            //- Synchronize values on selected mesh points.
-//            template <class T, class CombineOp, class TransformOp>
+            template <class T, class CombineOp, class TransformOp>
-//            static void syncPointList
+            static void syncPointList
-//            (
+            (
-//                const polyMesh&,
+                const polyMesh&,
-//                const labelList& meshPoints,
+                const labelList& meshPoints,
-//                UList<T>&,
+                List<T>&,
-//                const CombineOp& cop,
+                const CombineOp& cop,
-//                const T& nullValue,
+                const T& nullValue,
-//                const TransformOp& top
+                const TransformOp& top
-//            );
+            );
-//
+
-//            //- Synchronize values on all mesh edges.
+            //- Synchronize values on all mesh edges.
-//            template <class T, class CombineOp, class TransformOp>
+            template <class T, class CombineOp, class TransformOp>
-//            static void syncEdgeList
+            static void syncEdgeList
-//            (
+            (
-//                const polyMesh&,
+                const polyMesh&,
-//                UList<T>&,
+                List<T>&,
-//                const CombineOp& cop,
+                const CombineOp& cop,
-//                const T& nullValue,
+                const T& nullValue,
-//                const TransformOp& top
+                const TransformOp& top
-//            );
+            );
-//
+
             //- Synchronize values on boundary faces only.
             template <class T, class CombineOp, class TransformOp>
             static void syncBoundaryFaceList
@@ -219,7 +166,17 @@ public:
                 List<T>& l,
                 const CombineOp& cop,
                 const T& nullValue
+            )
+            {
+                syncPointList
+                (
+                    mesh,
+                    l,
+                    cop,
+                    nullValue,
+                    mapDistribute::transform()
                 );
+            }
 
             //- Synchronize locations on all mesh points.
             template <class CombineOp>
@@ -229,7 +186,17 @@ public:
                 List<point>& l,
                 const CombineOp& cop,
                 const point& nullValue
+            )
+            {
+                syncPointList
+                (
+                    mesh,
+                    l,
+                    cop,
+                    nullValue,
+                    mapDistribute::transformPosition()
                 );
+            }
 
             //- Synchronize values on selected mesh points.
             template <class T, class CombineOp>
@@ -237,10 +204,21 @@ public:
             (
                 const polyMesh& mesh,
                 const labelList& meshPoints,
-                UList<T>& l,
+                List<T>& l,
                 const CombineOp& cop,
                 const T& nullValue
+            )
+            {
+                syncPointList
+                (
+                    mesh,
+                    meshPoints,
+                    l,
+                    cop,
+                    nullValue,
+                    mapDistribute::transform()
                 );
+            }
 
             //- Synchronize locations on selected mesh points.
             template <class CombineOp>
@@ -248,10 +226,21 @@ public:
             (
                 const polyMesh& mesh,
                 const labelList& meshPoints,
-                UList<point>& l,
+                List<point>& l,
                 const CombineOp& cop,
                 const point& nullValue
+            )
+            {
+                syncPointList
+                (
+                    mesh,
+                    meshPoints,
+                    l,
+                    cop,
+                    nullValue,
+                    mapDistribute::transformPosition()
                 );
+            }
 
 
         // Synchronise edge-wise data
@@ -261,20 +250,40 @@ public:
             static void syncEdgeList
             (
                 const polyMesh& mesh,
-                UList<T>& l,
+                List<T>& l,
                 const CombineOp& cop,
                 const T& nullValue
+            )
+            {
+                syncEdgeList
+                (
+                    mesh,
+                    l,
+                    cop,
+                    nullValue,
+                    mapDistribute::transform()
                 );
+            }
 
             //- Synchronize values on all mesh edges.
             template <class CombineOp>
             static void syncEdgePositions
             (
                 const polyMesh& mesh,
-                UList<point>& l,
+                List<point>& l,
                 const CombineOp& cop,
                 const point& nullValue
+            )
+            {
+                syncEdgeList
+                (
+                    mesh,
+                    l,
+                    cop,
+                    nullValue,
+                    mapDistribute::transformPosition()
                 );
+            }
 
 
         // Synchronise face-wise data
@@ -288,7 +297,7 @@ public:
                 const CombineOp& cop
             )
             {
-                syncBoundaryFaceList(mesh, l, cop, transform());
+                syncBoundaryFaceList(mesh, l, cop, mapDistribute::transform());
             }
 
             //- Synchronize locations on boundary faces only.
@@ -300,7 +309,13 @@ public:
                 const CombineOp& cop
             )
             {
-                syncBoundaryFaceList(mesh, l, cop, transformPosition());
+                syncBoundaryFaceList
+                (
+                    mesh,
+                    l,
+                    cop,
+                    mapDistribute::transformPosition()
+                );
             }
 
             //- Synchronize values on all mesh faces.
@@ -319,7 +334,13 @@ public:
                     mesh.nInternalFaces()
                 );
 
-                syncBoundaryFaceList(mesh, bndValues, cop, transform());
+                syncBoundaryFaceList
+                (
+                    mesh,
+                    bndValues,
+                    cop,
+                    mapDistribute::transform()
+                );
             }
 
             //- Synchronize locations on all mesh faces.
@@ -337,7 +358,13 @@ public:
                     mesh.nFaces()-mesh.nInternalFaces(),
                     mesh.nInternalFaces()
                 );
-                syncBoundaryFaceList(mesh, bndValues, cop, transformPosition());
+                syncBoundaryFaceList
+                (
+                    mesh,
+                    bndValues,
+                    cop,
+                    mapDistribute::transformPosition()
+                );
             }
 
             //- Swap coupled boundary face values.
@@ -348,7 +375,13 @@ public:
                 UList<T>& l
             )
             {
-                syncBoundaryFaceList(mesh, l, eqOp<T>(), transform());
+                syncBoundaryFaceList
+                (
+                    mesh,
+                    l,
+                    eqOp<T>(),
+                    mapDistribute::transform()
+                );
             }
 
              //- Swap coupled positions.
@@ -359,7 +392,13 @@ public:
                 UList<T>& l
             )
             {
-                syncBoundaryFaceList(mesh, l, eqOp<T>(), transformPosition());
+                syncBoundaryFaceList
+                (
+                    mesh,
+                    l,
+                    eqOp<T>(),
+                    mapDistribute::transformPosition()
+                );
             }
 
             //- Swap coupled face values.
@@ -376,7 +415,13 @@ public:
                     mesh.nFaces()-mesh.nInternalFaces(),
                     mesh.nInternalFaces()
                 );
-                syncBoundaryFaceList(mesh, bndValues, eqOp<T>(), transform());
+                syncBoundaryFaceList
+                (
+                    mesh,
+                    bndValues,
+                    eqOp<T>(),
+                    mapDistribute::transform()
+                );
             }
 
         // Sparse versions
@@ -390,7 +435,7 @@ public:
                 const CombineOp& cop
             )
             {
-                syncPointMap(mesh, l, cop, transform());
+                syncPointMap(mesh, l, cop, mapDistribute::transform());
             }
 
             //- Synchronize locations on selected points.
@@ -402,7 +447,7 @@ public:
                 const CombineOp& cop
             )
             {
-                syncPointMap(mesh, l, cop, transformPosition());
+                syncPointMap(mesh, l, cop, mapDistribute::transformPosition());
             }
 
             //- Synchronize values on selected edges. Edges are represented
@@ -416,7 +461,7 @@ public:
                 const CombineOp& cop
             )
             {
-                syncEdgeMap(mesh, l, cop, transform());
+                syncEdgeMap(mesh, l, cop, mapDistribute::transform());
             }
 
             //- Synchronize locations on selected edges.
@@ -428,7 +473,7 @@ public:
                 const CombineOp& cop
             )
             {
-                syncEdgeMap(mesh, l, cop, transformPosition());
+                syncEdgeMap(mesh, l, cop, mapDistribute::transformPosition());
             }
 
         // PackedList versions
@@ -480,65 +525,6 @@ public:
 };
 
 
-template<>
-void syncTools::transform::operator()
-(
-    const coupledPolyPatch&,
-    Field<label>&
-) const;
-template<>
-void syncTools::transform::operator()
-(
-    const coupledPolyPatch&,
-    Map<label>&
-) const;
-template<>
-void syncTools::transform::operator()
-(
-    const coupledPolyPatch&,
-    EdgeMap<label>&
-) const;
-
-
-template<>
-void syncTools::transform::operator()
-(
-    const coupledPolyPatch&,
-    Field<scalar>&
-) const;
-template<>
-void syncTools::transform::operator()
-(
-    const coupledPolyPatch&,
-    Map<scalar>&
-) const;
-template<>
-void syncTools::transform::operator()
-(
-    const coupledPolyPatch&,
-    EdgeMap<scalar>&
-) const;
-
-template<>
-void syncTools::transform::operator()
-(
-    const coupledPolyPatch& cpp,
-    Field<bool>& fld
-) const;
-template<>
-void syncTools::transform::operator()
-(
-    const coupledPolyPatch&,
-    Map<bool>&
-) const;
-template<>
-void syncTools::transform::operator()
-(
-    const coupledPolyPatch&,
-    EdgeMap<bool>&
-) const;
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/meshes/polyMesh/syncTools/syncToolsTemplates.C

@@ -764,7 +764,7 @@ void Foam::syncTools::syncEdgeMap
 //void Foam::syncTools::syncPointList
 //(
 //    const polyMesh& mesh,
-//    UList<T>& pointValues,
+//    List<T>& pointValues,
 //    const CombineOp& cop,
 //    const T& nullValue,
 //    const TransformOp& top
@@ -775,7 +775,7 @@ void Foam::syncTools::syncEdgeMap
 //        FatalErrorIn
 //        (
 //            "syncTools<class T, class CombineOp>::syncPointList"
-//            "(const polyMesh&, UList<T>&, const CombineOp&, const T&"
+//            "(const polyMesh&, List<T>&, const CombineOp&, const T&"
 //            ", const bool)"
 //        )   << "Number of values " << pointValues.size()
 //            << " is not equal to the number of points in the mesh "
@@ -941,7 +941,7 @@ void Foam::syncTools::syncEdgeMap
 //(
 //    const polyMesh& mesh,
 //    const labelList& meshPoints,
-//    UList<T>& pointValues,
+//    List<T>& pointValues,
 //    const CombineOp& cop,
 //    const T& nullValue,
 //    const TransformOp& top
@@ -952,7 +952,7 @@ void Foam::syncTools::syncEdgeMap
 //        FatalErrorIn
 //        (
 //            "syncTools<class T, class CombineOp>::syncPointList"
-//            "(const polyMesh&, const labelList&, UList<T>&, const CombineOp&"
+//            "(const polyMesh&, const labelList&, List<T>&, const CombineOp&"
 //            ", const T&, const bool)"
 //        )   << "Number of values " << pointValues.size()
 //            << " is not equal to the number of points "
@@ -986,13 +986,14 @@ void Foam::syncTools::syncEdgeMap
 //    }
 //}
 
-template <class T, class CombineOp>
+template <class T, class CombineOp, class TransformOp>
 void Foam::syncTools::syncPointList
 (
     const polyMesh& mesh,
     List<T>& pointValues,
     const CombineOp& cop,
-    const T& nullValue
+    const T& nullValue,
+    const TransformOp& top
 )
 {
     if (pointValues.size() != mesh.nPoints())
@@ -1000,48 +1001,50 @@ void Foam::syncTools::syncPointList
         FatalErrorIn
         (
             "syncTools<class T, class CombineOp>::syncPointList"
-            "(const polyMesh&, UList<T>&, const CombineOp&, const T&)"
+            "(const polyMesh&, List<T>&, const CombineOp&, const T&"
+            ", const bool)"
         )   << "Number of values " << pointValues.size()
             << " is not equal to the number of points in the mesh "
             << mesh.nPoints() << abort(FatalError);
     }
 
-    mesh.globalData().syncPointData(pointValues, cop, false);
+    mesh.globalData().syncPointData(pointValues, cop, top);
 }
 
 
-template <class CombineOp>
+//template <class CombineOp>
-void Foam::syncTools::syncPointPositions
+//void Foam::syncTools::syncPointPositions
-(
+//(
-    const polyMesh& mesh,
+//    const polyMesh& mesh,
-    List<point>& pointValues,
+//    List<point>& pointValues,
-    const CombineOp& cop,
+//    const CombineOp& cop,
-    const point& nullValue
+//    const point& nullValue
-)
+//)
-{
+//{
-    if (pointValues.size() != mesh.nPoints())
+//    if (pointValues.size() != mesh.nPoints())
-    {
+//    {
-        FatalErrorIn
+//        FatalErrorIn
-        (
+//        (
-            "syncTools<class CombineOp>::syncPointPositions"
+//            "syncTools<class CombineOp>::syncPointPositions"
-            "(const polyMesh&, List<point>&, const CombineOp&, const point&)"
+//            "(const polyMesh&, List<point>&, const CombineOp&, const point&)"
-        )   << "Number of values " << pointValues.size()
+//        )   << "Number of values " << pointValues.size()
-            << " is not equal to the number of points in the mesh "
+//            << " is not equal to the number of points in the mesh "
-            << mesh.nPoints() << abort(FatalError);
+//            << mesh.nPoints() << abort(FatalError);
-    }
+//    }
-
+//
-    mesh.globalData().syncPointData(pointValues, cop, true);
+//    mesh.globalData().syncPointData(pointValues, cop, true);
-}
+//}
 
 
-template <class T, class CombineOp>
+template <class T, class CombineOp, class TransformOp>
 void Foam::syncTools::syncPointList
 (
     const polyMesh& mesh,
     const labelList& meshPoints,
-    UList<T>& pointValues,
+    List<T>& pointValues,
     const CombineOp& cop,
-    const T& nullValue
+    const T& nullValue,
+    const TransformOp& top
 )
 {
     if (pointValues.size() != meshPoints.size())
@@ -1049,7 +1052,7 @@ void Foam::syncTools::syncPointList
         FatalErrorIn
         (
             "syncTools<class T, class CombineOp>::syncPointList"
-            "(const polyMesh&, UList<T>&, const CombineOp&, const T&)"
+            "(const polyMesh&, List<T>&, const CombineOp&, const T&)"
         )   << "Number of values " << pointValues.size()
             << " is not equal to the number of meshPoints "
             << meshPoints.size() << abort(FatalError);
@@ -1078,7 +1081,7 @@ void Foam::syncTools::syncPointList
         gd.globalPointSlavesMap(),
         gd.globalTransforms(),
         cop,
-        false       //position?
+        top
     );
 
     forAll(meshPoints, i)
@@ -1093,72 +1096,74 @@ void Foam::syncTools::syncPointList
 }
 
 
-template <class CombineOp>
+//template <class CombineOp>
-void Foam::syncTools::syncPointPositions
+//void Foam::syncTools::syncPointPositions
-(
+//(
-    const polyMesh& mesh,
+//    const polyMesh& mesh,
-    const labelList& meshPoints,
+//    const labelList& meshPoints,
-    UList<point>& pointValues,
+//    List<point>& pointValues,
-    const CombineOp& cop,
+//    const CombineOp& cop,
-    const point& nullValue
+//    const point& nullValue
-)
+//)
-{
+//{
-    if (pointValues.size() != meshPoints.size())
+//    if (pointValues.size() != meshPoints.size())
-    {
+//    {
-        FatalErrorIn
+//        FatalErrorIn
-        (
+//        (
-            "syncTools<class CombineOp>::syncPointList"
+//            "syncTools<class CombineOp>::syncPointList"
-            "(const polyMesh&, UList<point>&, const CombineOp&, const point&)"
+//            "(const polyMesh&, List<point>&, const CombineOp&, const point&)"
-        )   << "Number of values " << pointValues.size()
+//        )   << "Number of values " << pointValues.size()
-            << " is not equal to the number of meshPoints "
+//            << " is not equal to the number of meshPoints "
-            << meshPoints.size() << abort(FatalError);
+//            << meshPoints.size() << abort(FatalError);
-    }
+//    }
-    const globalMeshData& gd = mesh.globalData();
+//    const globalMeshData& gd = mesh.globalData();
-    const indirectPrimitivePatch& cpp = gd.coupledPatch();
+//    const indirectPrimitivePatch& cpp = gd.coupledPatch();
-    const Map<label>& mpm = cpp.meshPointMap();
+//    const Map<label>& mpm = cpp.meshPointMap();
-
+//
-    List<point> cppFld(cpp.nPoints(), nullValue);
+//    List<point> cppFld(cpp.nPoints(), nullValue);
-
+//
-    forAll(meshPoints, i)
+//    forAll(meshPoints, i)
-    {
+//    {
-        label pointI = meshPoints[i];
+//        label pointI = meshPoints[i];
-        Map<label>::const_iterator iter = mpm.find(pointI);
+//        Map<label>::const_iterator iter = mpm.find(pointI);
-        if (iter != mpm.end())
+//        if (iter != mpm.end())
-        {
+//        {
-            cppFld[iter()] = pointValues[i];
+//            cppFld[iter()] = pointValues[i];
-        }
+//        }
-    }
+//    }
-
+//
-    globalMeshData::syncData
+//    globalMeshData::syncData
-    (
+//    (
-        cppFld,
+//        cppFld,
-        gd.globalPointSlaves(),
+//        gd.globalPointSlaves(),
-        gd.globalPointTransformedSlaves(),
+//        gd.globalPointTransformedSlaves(),
-        gd.globalPointSlavesMap(),
+//        gd.globalPointSlavesMap(),
-        gd.globalTransforms(),
+//        gd.globalTransforms(),
-        cop,
+//        cop,
-        true    //position?
+//        true,   //position?
-    );
+//        mapDistribute::transform()  // not used
-
+//    );
-    forAll(meshPoints, i)
+//
-    {
+//    forAll(meshPoints, i)
-        label pointI = meshPoints[i];
+//    {
-        Map<label>::const_iterator iter = mpm.find(pointI);
+//        label pointI = meshPoints[i];
-        if (iter != mpm.end())
+//        Map<label>::const_iterator iter = mpm.find(pointI);
-        {
+//        if (iter != mpm.end())
-            pointValues[i] = cppFld[iter()];
+//        {
-        }
+//            pointValues[i] = cppFld[iter()];
-    }
+//        }
-}
+//    }
+//}
 
 
-template <class T, class CombineOp>
+template <class T, class CombineOp, class TransformOp>
 void Foam::syncTools::syncEdgeList
 (
     const polyMesh& mesh,
-    UList<T>& edgeValues,
+    List<T>& edgeValues,
     const CombineOp& cop,
-    const T& nullValue
+    const T& nullValue,
+    const TransformOp& top
 )
 {
     if (edgeValues.size() != mesh.nEdges())
@@ -1166,7 +1171,7 @@ void Foam::syncTools::syncEdgeList
         FatalErrorIn
         (
             "syncTools<class T, class CombineOp>::syncEdgeList"
-            "(const polyMesh&, UList<T>&, const CombineOp&, const T&)"
+            "(const polyMesh&, List<T>&, const CombineOp&, const T&)"
         )   << "Number of values " << edgeValues.size()
             << " is not equal to the number of edges in the mesh "
             << mesh.nEdges() << abort(FatalError);
@@ -1187,7 +1192,7 @@ void Foam::syncTools::syncEdgeList
         edgeMap,
         git,
         cop,
-        false       //position?
+        top
     );
 
     // Extract back onto mesh
@@ -1198,50 +1203,51 @@ void Foam::syncTools::syncEdgeList
 }
 
 
-template <class CombineOp>
+//template <class CombineOp>
-void Foam::syncTools::syncEdgePositions
+//void Foam::syncTools::syncEdgePositions
-(
+//(
-    const polyMesh& mesh,
+//    const polyMesh& mesh,
-    UList<point>& edgeValues,
+//    List<point>& edgeValues,
-    const CombineOp& cop,
+//    const CombineOp& cop,
-    const point& nullValue
+//    const point& nullValue
-)
+//)
-{
+//{
-    if (edgeValues.size() != mesh.nEdges())
+//    if (edgeValues.size() != mesh.nEdges())
-    {
+//    {
-        FatalErrorIn
+//        FatalErrorIn
-        (
+//        (
-            "syncTools<class CombineOp>::syncEdgePositions"
+//            "syncTools<class CombineOp>::syncEdgePositions"
-            "(const polyMesh&, UList<point>&, const CombineOp&, const point&)"
+//            "(const polyMesh&, List<point>&, const CombineOp&, const point&)"
-        )   << "Number of values " << edgeValues.size()
+//        )   << "Number of values " << edgeValues.size()
-            << " is not equal to the number of edges in the mesh "
+//            << " is not equal to the number of edges in the mesh "
-            << mesh.nEdges() << abort(FatalError);
+//            << mesh.nEdges() << abort(FatalError);
-    }
+//    }
-
+//
-    const globalMeshData& gd = mesh.globalData();
+//    const globalMeshData& gd = mesh.globalData();
-    const labelList& meshEdges = gd.coupledPatchMeshEdges();
+//    const labelList& meshEdges = gd.coupledPatchMeshEdges();
-    const globalIndexAndTransform& git = gd.globalTransforms();
+//    const globalIndexAndTransform& git = gd.globalTransforms();
-    const mapDistribute& map = gd.globalEdgeSlavesMap();
+//    const mapDistribute& map = gd.globalEdgeSlavesMap();
-
+//
-    List<point> cppFld(UIndirectList<point>(edgeValues, meshEdges));
+//    List<point> cppFld(UIndirectList<point>(edgeValues, meshEdges));
-
+//
-    globalMeshData::syncData
+//    globalMeshData::syncData
-    (
+//    (
-        cppFld,
+//        cppFld,
-        gd.globalEdgeSlaves(),
+//        gd.globalEdgeSlaves(),
-        gd.globalEdgeTransformedSlaves(),
+//        gd.globalEdgeTransformedSlaves(),
-        map,
+//        map,
-        git,
+//        git,
-        cop,
+//        cop,
-        true        //position?
+//        true,       //position?
-    );
+//        mapDistribute::transform()  // not used
-
+//    );
-    // Extract back onto mesh
+//
-    forAll(meshEdges, i)
+//    // Extract back onto mesh
-    {
+//    forAll(meshEdges, i)
-        edgeValues[meshEdges[i]] = cppFld[i];
+//    {
-    }
+//        edgeValues[meshEdges[i]] = cppFld[i];
-}
+//    }
+//}
 
 
 template <class T, class CombineOp, class TransformOp>

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -29,6 +29,7 @@ License
 #include "PatchToolsEdgeOwner.C"
 #include "PatchToolsSearch.C"
 #include "PatchToolsSortEdges.C"
+#include "PatchToolsNormals.C"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchTools.H

@@ -35,6 +35,7 @@ SourceFiles
     PatchToolsEdgeOwner.C
     PatchToolsSearch.C
     PatchToolsSortEdges.C
+    PatchToolsNormals.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -49,6 +50,8 @@ SourceFiles
 namespace Foam
 {
 
+class polyMesh;
+
 /*---------------------------------------------------------------------------*\
                          Class PatchTools Declaration
 \*---------------------------------------------------------------------------*/
@@ -166,6 +169,22 @@ public:
     );
 
 
+    //- Return parallel consistent point normals for patches (on boundary faces)
+    //  using mesh points.
+    template
+    <
+        class Face,
+        template<class> class FaceList,
+        class PointField,
+        class PointType
+    >
+    static tmp<pointField> pointNormals
+    (
+        const polyMesh&,
+        const PrimitivePatch<Face, FaceList, PointField, PointType>&,
+        const labelList& meshFaces
+    );
+
 };
 
 

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

@@ -0,0 +1,227 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 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 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "PatchTools.H"
+#include "polyMesh.H"
+#include "indirectPrimitivePatch.H"
+#include "globalMeshData.H"
+
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    //- Transformation
+    class listTransform
+    {
+    public:
+
+        void operator()
+        (
+            const vectorTensorTransform& vt,
+            const bool forward,
+            List<List<point> >& fld
+        ) const
+        {
+            const tensor T
+            (
+                forward
+              ? vt.R()
+              : vt.R().T()
+            );
+
+            forAll(fld, i)
+            {
+                List<point>& elems = fld[i];
+                forAll(elems, elemI)
+                {
+                    elems[elemI] = transform(T, elems[elemI]);
+                }
+            }
+        }
+    };
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template
+<
+    class Face,
+    template<class> class FaceList,
+    class PointField,
+    class PointType
+>
+
+Foam::tmp<Foam::pointField>
+Foam::PatchTools::pointNormals
+(
+    const polyMesh& mesh,
+    const PrimitivePatch<Face, FaceList, PointField, PointType>& p,
+    const labelList& meshFaces
+)
+{
+    // Assume patch is smaller than the globalData().coupledPatch() (?) so
+    // loop over patch meshPoints.
+
+    const globalMeshData& globalData = mesh.globalData();
+    const indirectPrimitivePatch& coupledPatch = globalData.coupledPatch();
+    const Map<label>& coupledPatchMP = coupledPatch.meshPointMap();
+    const mapDistribute& map = globalData.globalPointSlavesMap();
+    const globalIndexAndTransform& transforms =
+        globalData.globalTransforms();
+
+
+    // 1. Start off with local normals (note:without calculating pointNormals
+    //    to avoid them being stored)
+
+    tmp<pointField> textrudeN(new pointField(p.nPoints(), vector::zero));
+    pointField& extrudeN = textrudeN();
+    {
+        const faceList& localFaces = p.localFaces();
+        const vectorField& faceNormals = p.faceNormals();
+
+        forAll(localFaces, faceI)
+        {
+            const face& f = localFaces[faceI];
+            const vector& n = faceNormals[faceI];
+            forAll(f, fp)
+            {
+                extrudeN[f[fp]] += n;
+            }
+        }
+        extrudeN /= mag(extrudeN)+VSMALL;
+    }
+
+
+    // Collect local pointFaces
+    List<List<point> > pointFaceNormals(map.constructSize());
+    forAll(p.meshPoints(), patchPointI)
+    {
+        label meshPointI = p.meshPoints()[patchPointI];
+        Map<label>::const_iterator fnd = coupledPatchMP.find(meshPointI);
+        if (fnd != coupledPatchMP.end())
+        {
+            label coupledPointI = fnd();
+
+            List<point>& pNormals = pointFaceNormals[coupledPointI];
+            const labelList& pFaces = p.pointFaces()[patchPointI];
+            pNormals.setSize(pFaces.size());
+            forAll(pFaces, i)
+            {
+                pNormals[i] = p.faceNormals()[pFaces[i]];
+            }
+        }
+    }
+
+
+    // Pull remote data into local slots
+    map.distribute
+    (
+        transforms,
+        pointFaceNormals,
+        listTransform()
+    );
+
+
+    // Combine normals
+    const labelListList& slaves = globalData.globalPointSlaves();
+    const labelListList& transformedSlaves =
+        globalData.globalPointTransformedSlaves();
+
+
+    pointField coupledPointNormals(map.constructSize(), vector::zero);
+
+    forAll(p.meshPoints(), patchPointI)
+    {
+        label meshPointI = p.meshPoints()[patchPointI];
+        Map<label>::const_iterator fnd = coupledPatchMP.find(meshPointI);
+        if (fnd != coupledPatchMP.end())
+        {
+            label coupledPointI = fnd();
+            const labelList& slaveSlots =
+                slaves[coupledPointI];
+            const labelList& transformedSlaveSlots =
+                transformedSlaves[coupledPointI];
+
+            label nFaces = slaveSlots.size()+transformedSlaveSlots.size();
+            if (nFaces > 0)
+            {
+                // Combine
+                point& n = coupledPointNormals[coupledPointI];
+
+                n += sum(pointFaceNormals[coupledPointI]);
+
+                forAll(slaveSlots, i)
+                {
+                    n += sum(pointFaceNormals[slaveSlots[i]]);
+                }
+                forAll(transformedSlaveSlots, i)
+                {
+                    n += sum(pointFaceNormals[transformedSlaveSlots[i]]);
+                }
+                n /= mag(n)+VSMALL;
+
+                // Put back into slave slots
+                forAll(slaveSlots, i)
+                {
+                    coupledPointNormals[slaveSlots[i]] = n;
+                }
+                forAll(transformedSlaveSlots, i)
+                {
+                    coupledPointNormals[transformedSlaveSlots[i]] = n;
+                }
+            }
+        }
+    }
+
+
+    // Send back
+    map.reverseDistribute
+    (
+        transforms,
+        coupledPointNormals.size(),
+        coupledPointNormals,
+        mapDistribute::transform()
+    );
+
+
+    // Override patch normals
+    forAll(p.meshPoints(), patchPointI)
+    {
+        label meshPointI = p.meshPoints()[patchPointI];
+        Map<label>::const_iterator fnd = coupledPatchMP.find(meshPointI);
+        if (fnd != coupledPatchMP.end())
+        {
+            label coupledPointI = fnd();
+            extrudeN[patchPointI] = coupledPointNormals[coupledPointI];
+        }
+    }
+
+    return textrudeN;
+}
+
+
+// ************************************************************************* //

src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C

@@ -43,75 +43,6 @@ defineTypeNameAndDebug(Foam::addPatchCellLayer, 0);
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-// Calculate global faces per pp edge.
-Foam::labelListList Foam::addPatchCellLayer::calcGlobalEdgeFaces
-(
-    const polyMesh& mesh,
-    const globalIndex& globalFaces,
-    const indirectPrimitivePatch& pp,
-    const labelList& meshEdges
-)
-{
-    //// Determine coupled edges just so we don't have to have storage
-    //// for all non-coupled edges.
-    //
-    //PackedBoolList isCoupledEdge(mesh.nEdges());
-    //
-    //const polyBoundaryMesh& patches = mesh.boundaryMesh();
-    //
-    //forAll(patches, patchI)
-    //{
-    //    const polyPatch& pp = patches[patchI];
-    //
-    //    if (pp.coupled())
-    //    {
-    //        const labelList& meshEdges = pp.meshEdges();
-    //
-    //        forAll(meshEdges, i)
-    //        {
-    //            isCoupledEdge.set(meshEdges[i], 1);
-    //        }
-    //    }
-    //}
-
-    // From mesh edge to global face labels. Sized only for pp edges.
-    labelListList globalEdgeFaces(mesh.nEdges());
-
-    const labelListList& edgeFaces = pp.edgeFaces();
-
-    forAll(edgeFaces, edgeI)
-    {
-        label meshEdgeI = meshEdges[edgeI];
-
-        //if (isCoupledEdge.get(meshEdgeI) == 1)
-        {
-            const labelList& eFaces = edgeFaces[edgeI];
-
-            // Store face and processor as unique tag.
-            labelList& globalEFaces = globalEdgeFaces[meshEdgeI];
-            globalEFaces.setSize(eFaces.size());
-            forAll(eFaces, i)
-            {
-                globalEFaces[i] =
-                    globalFaces.toGlobal(pp.addressing()[eFaces[i]]);
-            }
-        }
-    }
-
-    // Synchronise across coupled edges.
-    syncTools::syncEdgeList
-    (
-        mesh,
-        globalEdgeFaces,
-        uniqueEqOp(),
-        labelList()     // null value
-    );
-
-    // Extract pp part
-    return labelListList(UIndirectList<labelList>(globalEdgeFaces, meshEdges));
-}
-
-
 Foam::label Foam::addPatchCellLayer::nbrFace
 (
     const labelListList& edgeFaces,
@@ -316,12 +247,12 @@ Foam::labelPair Foam::addPatchCellLayer::getEdgeString
 Foam::label Foam::addPatchCellLayer::addSideFace
 (
     const indirectPrimitivePatch& pp,
-    const labelList& patchID,           // prestored patch per pp face
     const labelListList& addedCells,    // per pp face the new extruded cell
     const face& newFace,
+    const label newPatchID,
+
     const label ownFaceI,               // pp face that provides owner
     const label nbrFaceI,
-    const label patchEdgeI,             // edge to add to
     const label meshEdgeI,              // corresponding mesh edge
     const label layerI,                 // layer
     const label numEdgeFaces,           // number of layers for edge
@@ -329,8 +260,9 @@ Foam::label Foam::addPatchCellLayer::addSideFace
     polyTopoChange& meshMod
 ) const
 {
-    // Edge to 'inflate' from
+    // Face or edge to 'inflate' from
     label inflateEdgeI = -1;
+    label inflateFaceI = -1;
 
     // Check mesh faces using edge
     if (addToMesh_)
@@ -346,8 +278,6 @@ Foam::label Foam::addPatchCellLayer::addSideFace
         }
     }
 
-    // Get my mesh face and its zone.
-    label meshFaceI = pp.addressing()[ownFaceI];
     // Zone info comes from any side patch face. Otherwise -1 since we
     // don't know what to put it in - inherit from the extruded faces?
     label zoneI = -1;   //mesh_.faceZones().whichZone(meshFaceI);
@@ -358,14 +288,15 @@ Foam::label Foam::addPatchCellLayer::addSideFace
     // Is patch edge external edge of indirectPrimitivePatch?
     if (nbrFaceI == -1)
     {
-        // External edge so external face. Patch id is obtained from
+        // External edge so external face.
-        // any other patch connected to edge.
 
         const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
         // Loop over all faces connected to edge to inflate and
-        // see if any boundary face (but not meshFaceI)
+        // see if we can find a face that is otherPatchID
-        label otherPatchID = patchID[ownFaceI];
+
+        // Get my mesh face and its zone.
+        label meshFaceI = pp.addressing()[ownFaceI];
 
         forAll(meshFaces, k)
         {
@@ -373,11 +304,14 @@ Foam::label Foam::addPatchCellLayer::addSideFace
 
             if
             (
-                faceI != meshFaceI
+                (faceI != meshFaceI)
-             && !mesh_.isInternalFace(faceI)
+             && (patches.whichPatch(faceI) == newPatchID)
             )
             {
-                otherPatchID = patches.whichPatch(faceI);
+                // Found the patch face. Use it to inflate from
+                inflateEdgeI = -1;
+                inflateFaceI = faceI;
+
                 zoneI = mesh_.faceZones().whichZone(faceI);
                 if (zoneI != -1)
                 {
@@ -414,7 +348,7 @@ Foam::label Foam::addPatchCellLayer::addSideFace
 
         //Pout<< "Added boundary face:" << newFace
         //    << " own:" << addedCells[ownFaceI][layerOwn]
-        //    << " patch:" << otherPatchID
+        //    << " patch:" << newPatchID
         //    << endl;
 
         addedFaceI = meshMod.setAction
@@ -426,9 +360,9 @@ Foam::label Foam::addPatchCellLayer::addSideFace
                 -1,                         // neighbour
                 -1,                         // master point
                 inflateEdgeI,               // master edge
-                -1,                         // master face
+                inflateFaceI,               // master face
                 false,                      // flux flip
-                otherPatchID,               // patch for face
+                newPatchID,                 // patch for face
                 zoneI,                      // zone for face
                 flip                        // face zone flip
             )
@@ -510,6 +444,51 @@ Foam::label Foam::addPatchCellLayer::addSideFace
 }
 
 
+Foam::label Foam::addPatchCellLayer::findProcPatch
+(
+    const polyMesh& mesh,
+    const label nbrProcID
+)
+{
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    forAll(mesh.globalData().processorPatches(), i)
+    {
+        label patchI = mesh.globalData().processorPatches()[i];
+
+        if
+        (
+            refCast<const processorPolyPatch>(patches[patchI]).neighbProcNo()
+         == nbrProcID
+        )
+        {
+            return patchI;
+        }
+    }
+    return -1;
+}
+
+
+void Foam::addPatchCellLayer::setFaceProps
+(
+    const polyMesh& mesh,
+    const label faceI,
+
+    label& patchI,
+    label& zoneI,
+    bool& zoneFlip
+)
+{
+    patchI = mesh.boundaryMesh().whichPatch(faceI);
+    zoneI = mesh.faceZones().whichZone(faceI);
+    if (zoneI != -1)
+    {
+        label index = mesh.faceZones()[zoneI].whichFace(faceI);
+        zoneFlip = mesh.faceZones()[zoneI].flipMap()[index];
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct from mesh
@@ -561,10 +540,251 @@ Foam::labelListList Foam::addPatchCellLayer::addedCells() const
 }
 
 
+// Calculate global faces per pp edge.
+Foam::labelListList Foam::addPatchCellLayer::globalEdgeFaces
+(
+    const polyMesh& mesh,
+    const globalIndex& globalFaces,
+    const indirectPrimitivePatch& pp
+)
+{
+    // Precalculate mesh edges for pp.edges.
+    const labelList meshEdges(pp.meshEdges(mesh.edges(), mesh.pointEdges()));
+
+    // From mesh edge to global face labels. Non-empty sublists only for
+    // pp edges.
+    labelListList globalEdgeFaces(mesh.nEdges());
+
+    const labelListList& edgeFaces = pp.edgeFaces();
+
+    forAll(edgeFaces, edgeI)
+    {
+        label meshEdgeI = meshEdges[edgeI];
+
+        const labelList& eFaces = edgeFaces[edgeI];
+
+        // Store face and processor as unique tag.
+        labelList& globalEFaces = globalEdgeFaces[meshEdgeI];
+        globalEFaces.setSize(eFaces.size());
+        forAll(eFaces, i)
+        {
+            globalEFaces[i] = globalFaces.toGlobal(pp.addressing()[eFaces[i]]);
+        }
+    }
+
+    // Synchronise across coupled edges.
+    syncTools::syncEdgeList
+    (
+        mesh,
+        globalEdgeFaces,
+        uniqueEqOp(),
+        labelList()             // null value
+    );
+
+    // Extract pp part
+    return labelListList(UIndirectList<labelList>(globalEdgeFaces, meshEdges));
+}
+
+
+void Foam::addPatchCellLayer::calcSidePatch
+(
+    const polyMesh& mesh,
+    const globalIndex& globalFaces,
+    const labelListList& globalEdgeFaces,
+    const indirectPrimitivePatch& pp,
+
+    labelList& sidePatchID,
+    label& nPatches,
+    Map<label>& nbrProcToPatch,
+    Map<label>& patchToNbrProc
+)
+{
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    // Precalculate mesh edges for pp.edges.
+    const labelList meshEdges(pp.meshEdges(mesh.edges(), mesh.pointEdges()));
+
+    sidePatchID.setSize(pp.nEdges());
+    sidePatchID = -1;
+
+    // These also get determined but not (yet) exported:
+    // - whether face is created from other face or edge
+    // - what zone&orientation face should have
+
+    labelList inflateEdgeI(pp.nEdges(), -1);
+    labelList inflateFaceI(pp.nEdges(), -1);
+    labelList sideZoneID(pp.nEdges(), -1);
+    boolList sideFlip(pp.nEdges(), false);
+
+    nPatches = patches.size();
+
+    forAll(globalEdgeFaces, edgeI)
+    {
+        const labelList& eGlobalFaces = globalEdgeFaces[edgeI];
+        if
+        (
+            eGlobalFaces.size() == 2
+         && pp.edgeFaces()[edgeI].size() == 1
+        )
+        {
+            // Locally but not globally a boundary edge. Hence a coupled
+            // edge. Find the patch to use if on different
+            // processors.
+
+            label f0 = eGlobalFaces[0];
+            label f1 = eGlobalFaces[1];
+
+            label otherProcI = -1;
+            if (globalFaces.isLocal(f0) && !globalFaces.isLocal(f1))
+            {
+                otherProcI = globalFaces.whichProcID(f1);
+            }
+            else if (!globalFaces.isLocal(f0) && globalFaces.isLocal(f1))
+            {
+                otherProcI = globalFaces.whichProcID(f0);
+            }
+
+
+            if (otherProcI != -1)
+            {
+                sidePatchID[edgeI] = findProcPatch(mesh, otherProcI);
+                if (sidePatchID[edgeI] == -1)
+                {
+                    // Cannot find a patch to processor. See if already
+                    // marked for addition
+                    if (nbrProcToPatch.found(otherProcI))
+                    {
+                        sidePatchID[edgeI] = nbrProcToPatch[otherProcI];
+                    }
+                    else
+                    {
+                        sidePatchID[edgeI] = nPatches;
+                        nbrProcToPatch.insert(otherProcI, nPatches);
+                        patchToNbrProc.insert(nPatches, otherProcI);
+                        nPatches++;
+                    }
+                }
+            }
+        }
+    }
+
+
+
+    // Determine face properties for all other boundary edges
+    // ------------------------------------------------------
+
+    const labelListList& edgeFaces = pp.edgeFaces();
+    forAll(edgeFaces, edgeI)
+    {
+        if (edgeFaces[edgeI].size() == 1 && sidePatchID[edgeI] == -1)
+        {
+            // Proper, uncoupled patch edge.
+
+            label myFaceI = pp.addressing()[edgeFaces[edgeI][0]];
+
+            // Pick up any boundary face on this edge and use its properties
+            label meshEdgeI = meshEdges[edgeI];
+            const labelList& meshFaces = mesh.edgeFaces()[meshEdgeI];
+
+            forAll(meshFaces, k)
+            {
+                label faceI = meshFaces[k];
+
+                if (faceI != myFaceI && !mesh.isInternalFace(faceI))
+                {
+                    setFaceProps
+                    (
+                        mesh,
+                        faceI,
+
+                        sidePatchID[edgeI],
+                        sideZoneID[edgeI],
+                        sideFlip[edgeI]
+                    );
+                    inflateFaceI[edgeI] = faceI;
+                    inflateEdgeI[edgeI] = -1;
+
+                    break;
+                }
+            }
+        }
+    }
+
+
+
+    // Now hopefully every boundary edge has a side patch. Check
+    forAll(edgeFaces, edgeI)
+    {
+        if (edgeFaces[edgeI].size() == 1 && sidePatchID[edgeI] == -1)
+        {
+            const edge& e = pp.edges()[edgeI];
+            FatalErrorIn("addPatchCellLayer::calcSidePatch(..)")
+                << "Have no sidePatchID for edge " << edgeI << " points "
+                << pp.points()[pp.meshPoints()[e[0]]]
+                << pp.points()[pp.meshPoints()[e[1]]]
+                << abort(FatalError);
+        }
+    }
+
+
+
+    // Now we have sidepatch see if we have patchface or edge to inflate
+    // from.
+    forAll(edgeFaces, edgeI)
+    {
+        if (edgeFaces[edgeI].size() == 1 && inflateFaceI[edgeI] == -1)
+        {
+            // 1. Do we have a boundary face to inflate from
+
+            label myFaceI = pp.addressing()[edgeFaces[edgeI][0]];
+
+            // Pick up any boundary face on this edge and use its properties
+            label meshEdgeI = meshEdges[edgeI];
+            const labelList& meshFaces = mesh.edgeFaces()[meshEdgeI];
+
+            forAll(meshFaces, k)
+            {
+                label faceI = meshFaces[k];
+
+                if (faceI != myFaceI)
+                {
+                    if (mesh.isInternalFace(faceI))
+                    {
+                        inflateEdgeI[edgeI] = meshEdgeI;
+                    }
+                    else
+                    {
+                        if (patches.whichPatch(faceI) == sidePatchID[edgeI])
+                        {
+                            setFaceProps
+                            (
+                                mesh,
+                                faceI,
+
+                                sidePatchID[edgeI],
+                                sideZoneID[edgeI],
+                                sideFlip[edgeI]
+                            );
+                            inflateFaceI[edgeI] = faceI;
+                            inflateEdgeI[edgeI] = -1;
+
+                            break;
+                        }
+                    }
+                }
+            }
+        }
+    }
+}
+
+
 void Foam::addPatchCellLayer::setRefinement
 (
+    const globalIndex& globalFaces,
+    const labelListList& globalEdgeFaces,
     const scalarField& expansionRatio,
     const indirectPrimitivePatch& pp,
+    const labelList& sidePatchID,
     const labelList& exposedPatchID,
     const labelList& nFaceLayers,
     const labelList& nPointLayers,
@@ -575,7 +795,7 @@ void Foam::addPatchCellLayer::setRefinement
     if (debug)
     {
         Pout<< "addPatchCellLayer::setRefinement : Adding up to "
-            << max(nPointLayers)
+            << gMax(nPointLayers)
             << " layers of cells to indirectPrimitivePatch with "
             << pp.nPoints() << " points" << endl;
     }
@@ -788,8 +1008,6 @@ void Foam::addPatchCellLayer::setRefinement
 
                 label meshEdgeI = meshEdges[edgeI];
 
-                // Mesh faces using edge
-
                 // Mesh faces using edge
                 const labelList& meshFaces = mesh_.edgeFaces(meshEdgeI, ef);
 
@@ -1213,22 +1431,6 @@ void Foam::addPatchCellLayer::setRefinement
     }
 
 
-    // Global indices engine
-    const globalIndex globalFaces(mesh_.nFaces());
-
-    // Get for all pp edgeFaces a unique faceID
-    labelListList globalEdgeFaces
-    (
-         calcGlobalEdgeFaces
-         (
-            mesh_,
-            globalFaces,
-            pp,
-            meshEdges
-        )
-    );
-
-
     // Mark off which edges have been extruded
     boolList doneEdge(pp.nEdges(), false);
 
@@ -1474,16 +1676,17 @@ void Foam::addPatchCellLayer::setRefinement
                         addSideFace
                         (
                             pp,
-                            patchID,
                             addedCells,
-                            newFace,
+
+                            newFace,                // vertices of new face
+                            sidePatchID[startEdgeI],// -1 or patch for face
+
                             patchFaceI,
                             nbrFaceI,
-                            startEdgeI,     // edge to inflate from
+                            meshEdgeI,          // (mesh) edge to inflate
-                            meshEdgeI,      // corresponding mesh edge
+                            i,                  // layer
-                            i,
+                            numEdgeSideFaces,   // num layers
-                            numEdgeSideFaces,
+                            meshFaces,          // edgeFaces
-                            meshFaces,
                             meshMod
                         );
                     }

src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.H

@@ -181,16 +181,6 @@ class addPatchCellLayer
 
     // Private Member Functions
 
-        //- Per patch edge the pp faces (in global indices) using it. Uses
-        //  uniqueEqOp() to remove duplicates.
-        labelListList calcGlobalEdgeFaces
-        (
-            const polyMesh& mesh,
-            const globalIndex& globalFaces,
-            const indirectPrimitivePatch& pp,
-            const labelList& meshEdges
-        );
-
         //- Get the face on the other side of the edge.
         static label nbrFace
         (
@@ -226,12 +216,13 @@ class addPatchCellLayer
         label addSideFace
         (
             const indirectPrimitivePatch&,
-            const labelList& patchID,
             const labelListList& addedCells,
+
             const face& newFace,
+            const label newPatchID,
+
             const label ownFaceI,
             const label nbrFaceI,
-            const label patchEdgeI,
             const label meshEdgeI,
             const label layerI,
             const label numEdgeFaces,
@@ -239,6 +230,18 @@ class addPatchCellLayer
             polyTopoChange&
         ) const;
 
+        //- Find patch to neighbouring processor
+        static label findProcPatch(const polyMesh&, const label nbrProcID);
+
+        //- Extract properties from mesh face
+        static void setFaceProps
+        (
+            const polyMesh&,
+            const label,
+            label&,
+            label&,
+            bool&
+        );
 
         //- Disallow default bitwise copy construct
         addPatchCellLayer(const addPatchCellLayer&);
@@ -256,7 +259,7 @@ public:
     // Constructors
 
         //- Construct from mesh.
-        addPatchCellLayer(const polyMesh& mesh, const bool addToMesh = true);
+        addPatchCellLayer(const polyMesh&, const bool addToMesh = true);
 
 
     // Member Functions
@@ -291,6 +294,33 @@ public:
 
         // Edit
 
+            //- Per patch edge the pp faces (in global indices) using it. Uses
+            //  uniqueEqOp() to remove duplicates.
+            static labelListList globalEdgeFaces
+            (
+                const polyMesh&,
+                const globalIndex& globalFaces,
+                const indirectPrimitivePatch& pp
+            );
+
+            //- Boundary edges get extruded into boundary faces. Determine patch
+            //  for these faces. This might be a to-be-created processor patch
+            //  (patchI >= mesh.boundaryMesh().size()) in which case the
+            //  nbrProcToPatch, patchToNbrProc give the correspondence. nPatches
+            //  is the new number of patches.
+            static void calcSidePatch
+            (
+                const polyMesh&,
+                const globalIndex& globalFaces,
+                const labelListList& globalEdgeFaces,
+                const indirectPrimitivePatch& pp,
+
+                labelList& sidePatchID,
+                label& nPatches,
+                Map<label>& nbrProcToPatch,
+                Map<label>& patchToNbrProc
+            );
+
             //- Play commands into polyTopoChange to create layers on top
             //  of indirectPrimitivePatch (have to be outside faces).
             //  Gets displacement per patch point.
@@ -313,8 +343,11 @@ public:
             //        (instead of e.g. from patch faces)
             void setRefinement
             (
+                const globalIndex& globalFaces,
+                const labelListList& globalEdgeFaces,
                 const scalarField& expansionRatio,
                 const indirectPrimitivePatch& pp,
+                const labelList& sidePatchID,
                 const labelList& exposedPatchID,
                 const labelList& nFaceLayers,
                 const labelList& nPointLayers,
@@ -326,20 +359,26 @@ public:
             //- Add with constant expansion ratio and same nLayers everywhere
             void setRefinement
             (
+                const globalIndex& globalFaces,
+                const labelListList& globalEdgeFaces,
                 const label nLayers,
                 const indirectPrimitivePatch& pp,
+                const labelList& sidePatchID,
                 const vectorField& overallDisplacement,
                 polyTopoChange& meshMod
             )
             {
                 setRefinement
                 (
+                    globalFaces,
+                    globalEdgeFaces,
                     scalarField(pp.nPoints(), 1.0),     // expansion ration
                     pp,
+                    sidePatchID,
                     labelList(0),
-                    labelList(pp.size(), nLayers),
+                    labelList(pp.size(), nLayers),      // nFaceLayers
-                    labelList(pp.nPoints(), nLayers),
+                    labelList(pp.nPoints(), nLayers),   // nPointLayers
-                    overallDisplacement / nLayers,
+                    overallDisplacement / nLayers,      // firstLayerDisp
                     meshMod
                 );
             }

src/edgeMesh/Make/files

@@ -20,7 +20,10 @@ $(edgeFormats)/starcd/STARCDedgeFormatRunTime.C
 $(edgeFormats)/vtk/VTKedgeFormat.C
 $(edgeFormats)/vtk/VTKedgeFormatRunTime.C
 
+extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C
+
 featureEdgeMesh/featureEdgeMesh.C
 
 
+
 LIB = $(FOAM_LIBBIN)/libedgeMesh

src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.H

@@ -0,0 +1,377 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2004-2011 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 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::extendedFeatureEdgeMesh
+
+Description
+
+    Description of feature edges and points.
+
+    Feature points are a sorted subset at the start of the overall points list:
+        0 .. concaveStart_-1                : convex points (w.r.t normals)
+        concaveStart_-1 .. mixedStart_-1    : concave points
+        mixedStart_ .. nonFeatureStart_     : mixed internal/external points
+        nonFeatureStart_ .. size-1          : non-feature points
+
+    Feature edges are the edgeList of the edgeMesh and are sorted:
+        0 .. internalStart_-1           : external edges (convex w.r.t normals)
+        internalStart_ .. flatStart_-1  : internal edges (concave)
+        flatStart_ .. openStart_-1      : flat edges (neither concave or convex)
+                                          can arise from region interfaces on
+                                          flat surfaces
+        openStart_ .. multipleStart_-1  : open edges (e.g. from baffle surfaces)
+        multipleStart_ .. size-1        : multiply connected edges
+
+    The edge direction and feature edge and feature point adjacent normals
+    are stored.
+
+SourceFiles
+    extendedFeatureEdgeMeshI.H
+    extendedFeatureEdgeMesh.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef extendedFeatureEdgeMesh_H
+#define extendedFeatureEdgeMesh_H
+
+#include "edgeMesh.H"
+#include "surfaceFeatures.H"
+#include "objectRegistry.H"
+#include "IOdictionary.H"
+#include "indexedOctree.H"
+#include "treeDataEdge.H"
+#include "pointIndexHit.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class extendedFeatureEdgeMesh Declaration
+\*---------------------------------------------------------------------------*/
+
+class extendedFeatureEdgeMesh
+:
+    public regIOobject,
+    public edgeMesh
+{
+
+public:
+
+    //- Runtime type information
+    TypeName("extendedFeatureEdgeMesh");
+
+    enum pointStatus
+    {
+        CONVEX,     // Fully convex point (w.r.t normals)
+        CONCAVE,    // Fully concave point
+        MIXED,      // A point surrounded by both convex and concave edges
+        NONFEATURE  // Not a feature point
+    };
+
+    enum edgeStatus
+    {
+        EXTERNAL,   // "Convex" edge
+        INTERNAL,   // "Concave" edge
+        FLAT,       // Neither concave or convex, on a flat surface
+        OPEN,       // i.e. only connected to one face
+        MULTIPLE,   // Multiply connected (connected to more than two faces)
+        NONE        // Not a classified feature edge (consistency with
+                    // surfaceFeatures)
+    };
+
+private:
+
+    // Static data
+
+        //- Angular closeness tolerance for treating normals as the same
+        static scalar cosNormalAngleTol_;
+
+        //- Index of the start of the convex feature points - static as 0
+        static label convexStart_;
+
+        //- Index of the start of the external feature edges - static as 0
+        static label externalStart_;
+
+
+    // Private data
+
+        //- Index of the start of the concave feature points
+        label concaveStart_;
+
+        //- Index of the start of the mixed type feature points
+        label mixedStart_;
+
+        //- Index of the start of the non-feature points
+        label nonFeatureStart_;
+
+        //- Index of the start of the internal feature edges
+        label internalStart_;
+
+        //- Index of the start of the flat feature edges
+        label flatStart_;
+
+        //- Index of the start of the open feature edges
+        label openStart_;
+
+        //- Index of the start of the multiply-connected feature edges
+        label multipleStart_;
+
+        //- Normals of the features, to be referred to by index by both feature
+        //  points and edges, unsorted
+        vectorField normals_;
+
+        //- Flat and open edges require the direction of the edge
+        vectorField edgeDirections_;
+
+        //- Indices of the normals that are adjacent to the feature edges
+        labelListList edgeNormals_;
+
+        //- Indices of the normals that are adjacent to the feature points
+        labelListList featurePointNormals_;
+
+        //- Feature edges which are on the boundary between regions
+        labelList regionEdges_;
+
+        //- Search tree for all edges
+        mutable autoPtr<indexedOctree<treeDataEdge> > edgeTree_;
+
+        //- Individual search trees for each type of edge
+        mutable PtrList<indexedOctree<treeDataEdge> > edgeTreesByType_;
+
+
+    // Private Member Functions
+
+        //- Classify the type of feature point.  Requires valid stored member
+        //  data for edges and normals.
+        pointStatus classifyFeaturePoint(label ptI) const;
+
+        //- Classify the type of feature edge.  Requires face centre 0 to face
+        //  centre 1 vector to distinguish internal from external
+        edgeStatus classifyEdge
+        (
+            const List<vector>& norms,
+            const labelList& edNorms,
+            const vector& fC0tofC1
+        ) const;
+
+
+public:
+
+    // Static data
+
+        //- Number of possible point types (i.e. number of slices)
+        static label nPointTypes;
+
+        //- Number of possible feature edge types (i.e. number of slices)
+        static label nEdgeTypes;
+
+    // Constructors
+
+        //- Construct (read) given an IOobject
+        extendedFeatureEdgeMesh(const IOobject&);
+
+        //- Construct as copy
+        explicit extendedFeatureEdgeMesh
+        (
+            const IOobject&,
+            const extendedFeatureEdgeMesh&
+        );
+
+        //- Construct by transferring components (points, edges)
+        extendedFeatureEdgeMesh
+        (
+            const IOobject&,
+            const Xfer<pointField>&,
+            const Xfer<edgeList>&
+        );
+
+        //- Construct (read) given surfaceFeatures, an objectRegistry and a
+        //  fileName to write to.  Extracts, classifies and reorders the data
+        //  from surfaceFeatures.
+        extendedFeatureEdgeMesh
+        (
+            const surfaceFeatures& sFeat,
+            const objectRegistry& obr,
+            const fileName& sFeatFileName
+        );
+
+        //- Construct from all components
+        extendedFeatureEdgeMesh
+        (
+            const IOobject& io,
+            const pointField& pts,
+            const edgeList& eds,
+            label concaveStart,
+            label mixedStart,
+            label nonFeatureStart,
+            label internalStart,
+            label flatStart,
+            label openStart,
+            label multipleStart,
+            const vectorField& normals,
+            const vectorField& edgeDirections,
+            const labelListList& edgeNormals,
+            const labelListList& featurePointNormals,
+            const labelList& regionEdges
+        );
+
+
+    //- Destructor
+    ~extendedFeatureEdgeMesh();
+
+
+    // Member Functions
+
+        // Find
+
+            //- Find nearest surface edge for the sample point.
+            void nearestFeatureEdge
+            (
+                const point& sample,
+                scalar searchDistSqr,
+                pointIndexHit& info
+            ) const;
+
+            //- Find nearest surface edge for each sample point.
+            void nearestFeatureEdge
+            (
+                const pointField& samples,
+                const scalarField& searchDistSqr,
+                List<pointIndexHit>& info
+            ) const;
+
+            //- Find the nearest point on each type of feature edge
+            void nearestFeatureEdgeByType
+            (
+                const point& sample,
+                const scalarField& searchDistSqr,
+                List<pointIndexHit>& info
+            ) const;
+
+        // Access
+
+            //- Return the index of the start of the convex feature points
+            inline label convexStart() const;
+
+            //- Return the index of the start of the concave feature points
+            inline label concaveStart() const;
+
+            //- Return the index of the start of the mixed type feature points
+            inline label mixedStart() const;
+
+            //- Return the index of the start of the non-feature points
+            inline label nonFeatureStart() const;
+
+            //- Return the index of the start of the external feature edges
+            inline label externalStart() const;
+
+            //- Return the index of the start of the internal feature edges
+            inline label internalStart() const;
+
+            //- Return the index of the start of the flat feature edges
+            inline label flatStart() const;
+
+            //- Return the index of the start of the open feature edges
+            inline label openStart() const;
+
+            //- Return the index of the start of the multiply-connected feature
+            //  edges
+            inline label multipleStart() const;
+
+            //- Return whether or not the point index is a feature point
+            inline bool featurePoint(label ptI) const;
+
+            //- Return the normals of the surfaces adjacent to the feature edges
+            //  and points
+            inline const vectorField& normals() const;
+
+            //- Return the edgeDirection vectors
+            inline const vectorField& edgeDirections() const;
+
+            //- Return the direction of edgeI, pointing away from ptI
+            inline vector edgeDirection(label edgeI, label ptI) const;
+
+            //- Return the indices of the normals that are adjacent to the
+            //  feature edges
+            inline const labelListList& edgeNormals() const;
+
+            //- Return the normal vectors for a given set of normal indices
+            inline vectorField edgeNormals(const labelList& edgeNormIs) const;
+
+            //- Return the normal vectors for a given edge
+            inline vectorField edgeNormals(label edgeI) const;
+
+            //- Return the indices of the normals that are adjacent to the
+            //  feature points
+            inline const labelListList& featurePointNormals() const;
+
+            //- Return the normal vectors for a given feature point
+            inline vectorField featurePointNormals(label ptI) const;
+
+            //- Return the feature edges which are on the boundary between
+            //  regions
+            inline const labelList& regionEdges() const;
+
+            //- Return the pointStatus of a specified point
+            inline pointStatus getPointStatus(label ptI) const;
+
+            //- Return the edgeStatus of a specified edge
+            inline edgeStatus getEdgeStatus(label edgeI) const;
+
+            //- Demand driven construction of octree for boundary edges
+            const indexedOctree<treeDataEdge>& edgeTree() const;
+
+            //- Demand driven construction of octree for boundary edges by type
+            const PtrList<indexedOctree<treeDataEdge> >&
+            edgeTreesByType() const;
+
+
+        // Write
+
+            //- Write all components of the extendedFeatureEdgeMesh as obj files
+            void writeObj(const fileName& prefix) const;
+
+            //- Give precedence to the regIOobject write
+            using regIOobject::write;
+
+            //- WriteData function required for regIOobject write operation
+            virtual bool writeData(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "extendedFeatureEdgeMeshI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMeshI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2009-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2009-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,78 +25,79 @@ License
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline Foam::label Foam::featureEdgeMesh::convexStart() const
+inline Foam::label Foam::extendedFeatureEdgeMesh::convexStart() const
 {
     return convexStart_;
 }
 
 
-inline Foam::label Foam::featureEdgeMesh::concaveStart() const
+inline Foam::label Foam::extendedFeatureEdgeMesh::concaveStart() const
 {
     return concaveStart_;
 }
 
 
-inline Foam::label Foam::featureEdgeMesh::mixedStart() const
+inline Foam::label Foam::extendedFeatureEdgeMesh::mixedStart() const
 {
     return mixedStart_;
 }
 
 
-inline Foam::label Foam::featureEdgeMesh::nonFeatureStart() const
+inline Foam::label Foam::extendedFeatureEdgeMesh::nonFeatureStart() const
 {
     return nonFeatureStart_;
 }
 
 
-inline Foam::label Foam::featureEdgeMesh::externalStart() const
+inline Foam::label Foam::extendedFeatureEdgeMesh::externalStart() const
 {
     return externalStart_;
 }
 
 
-inline Foam::label Foam::featureEdgeMesh::internalStart() const
+inline Foam::label Foam::extendedFeatureEdgeMesh::internalStart() const
 {
     return internalStart_;
 }
 
 
-inline Foam::label Foam::featureEdgeMesh::flatStart() const
+inline Foam::label Foam::extendedFeatureEdgeMesh::flatStart() const
 {
     return flatStart_;
 }
 
 
-inline Foam::label Foam::featureEdgeMesh::openStart() const
+inline Foam::label Foam::extendedFeatureEdgeMesh::openStart() const
 {
     return openStart_;
 }
 
 
-inline Foam::label Foam::featureEdgeMesh::multipleStart() const
+inline Foam::label Foam::extendedFeatureEdgeMesh::multipleStart() const
 {
     return multipleStart_;
 }
 
 
-inline bool Foam::featureEdgeMesh::featurePoint(label ptI) const
+inline bool Foam::extendedFeatureEdgeMesh::featurePoint(label ptI) const
 {
     return ptI < nonFeatureStart_;
 }
 
 
-inline const Foam::vectorField& Foam::featureEdgeMesh::normals() const
+inline const Foam::vectorField& Foam::extendedFeatureEdgeMesh::normals() const
 {
     return normals_;
 }
 
-inline const Foam::vectorField& Foam::featureEdgeMesh::edgeDirections() const
+inline const Foam::vectorField& Foam::extendedFeatureEdgeMesh::edgeDirections()
+const
 {
     return edgeDirections_;
 }
 
 
-inline Foam::vector Foam::featureEdgeMesh::edgeDirection
+inline Foam::vector Foam::extendedFeatureEdgeMesh::edgeDirection
 (
     label edgeI,
     label ptI
@@ -114,7 +115,7 @@ inline Foam::vector Foam::featureEdgeMesh::edgeDirection
     }
     else
     {
-        FatalErrorIn("Foam::featureEdgeMesh::edgedirection")
+        FatalErrorIn("Foam::extendedFeatureEdgeMesh::edgedirection")
             << "Requested ptI " << ptI << " is not a point on the requested "
             << "edgeI " << edgeI << ". edgeI start and end: "
             << e.start() << " " << e.end()
@@ -125,13 +126,14 @@ inline Foam::vector Foam::featureEdgeMesh::edgeDirection
 }
 
 
-inline const Foam::labelListList& Foam::featureEdgeMesh::edgeNormals() const
+inline const Foam::labelListList& Foam::extendedFeatureEdgeMesh::edgeNormals()
+const
 {
     return edgeNormals_;
 }
 
 
-inline Foam::vectorField Foam::featureEdgeMesh::edgeNormals
+inline Foam::vectorField Foam::extendedFeatureEdgeMesh::edgeNormals
 (
     const labelList& edgeNormIs
 ) const
@@ -147,27 +149,28 @@ inline Foam::vectorField Foam::featureEdgeMesh::edgeNormals
 }
 
 
-inline Foam::vectorField Foam::featureEdgeMesh::edgeNormals(label edgeI) const
+inline Foam::vectorField Foam::extendedFeatureEdgeMesh::edgeNormals(label edgeI)
+const
 {
     return edgeNormals(edgeNormals_[edgeI]);
 }
 
 
 inline const Foam::labelListList&
-Foam::featureEdgeMesh::featurePointNormals() const
+Foam::extendedFeatureEdgeMesh::featurePointNormals() const
 {
     return featurePointNormals_;
 }
 
 
-inline Foam::vectorField Foam::featureEdgeMesh::featurePointNormals
+inline Foam::vectorField Foam::extendedFeatureEdgeMesh::featurePointNormals
 (
     label ptI
 ) const
 {
     if (!featurePoint(ptI))
     {
-        WarningIn("vectorField featureEdgeMesh::featurePointNormals")
+        WarningIn("vectorField extendedFeatureEdgeMesh::featurePointNormals")
             << "Requesting the normals of a non-feature point. "
             << "Returned zero length vectorField."
             << endl;
@@ -188,13 +191,14 @@ inline Foam::vectorField Foam::featureEdgeMesh::featurePointNormals
 }
 
 
-inline const Foam::labelList& Foam::featureEdgeMesh::regionEdges() const
+inline const Foam::labelList& Foam::extendedFeatureEdgeMesh::regionEdges() const
 {
     return regionEdges_;
 }
 
 
-inline Foam::featureEdgeMesh::pointStatus Foam::featureEdgeMesh::getPointStatus
+inline Foam::extendedFeatureEdgeMesh::pointStatus
+Foam::extendedFeatureEdgeMesh::getPointStatus
 (
     label ptI
 ) const
@@ -218,7 +222,8 @@ inline Foam::featureEdgeMesh::pointStatus Foam::featureEdgeMesh::getPointStatus
 }
 
 
-inline Foam::featureEdgeMesh::edgeStatus Foam::featureEdgeMesh::getEdgeStatus
+inline Foam::extendedFeatureEdgeMesh::edgeStatus
+Foam::extendedFeatureEdgeMesh::getEdgeStatus
 (
     label edgeI
 ) const

src/edgeMesh/featureEdgeMesh/featureEdgeMesh.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,27 +25,12 @@ Class
     Foam::featureEdgeMesh
 
 Description
+    edgeMesh + IO.
 
-    Feature points are a sorted subset at the start of the overall points list:
+    See also extendedFeatureEdgeMesh type which stores additional classification
-        0 .. concaveStart_-1                : convex points (w.r.t normals)
+    of features.
-        concaveStart_-1 .. mixedStart_-1    : concave points
-        mixedStart_ .. nonFeatureStart_     : mixed internal/external points
-        nonFeatureStart_ .. size-1          : non-feature points
-
-    Feature edges are the edgeList of the edgeMesh and are sorted:
-        0 .. internalStart_-1           : external edges (convex w.r.t normals)
-        internalStart_ .. flatStart_-1  : internal edges (concave)
-        flatStart_ .. openStart_-1      : flat edges (neither concave or convex)
-                                          can arise from region interfaces on
-                                          flat surfaces
-        openStart_ .. multipleStart_-1  : open edges (e.g. from baffle surfaces)
-        multipleStart_ .. size-1        : multiply connected edges
-
-    The edge direction and feature edge and feature point adjacent normals
-    are stored.
 
 SourceFiles
-    featureEdgeMeshI.H
     featureEdgeMesh.C
 
 \*---------------------------------------------------------------------------*/
@@ -54,12 +39,7 @@ SourceFiles
 #define featureEdgeMesh_H
 
 #include "edgeMesh.H"
-#include "surfaceFeatures.H"
+#include "regIOobject.H"
-#include "objectRegistry.H"
-#include "IOdictionary.H"
-#include "indexedOctree.H"
-#include "treeDataEdge.H"
-#include "pointIndexHit.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -78,278 +58,28 @@ class featureEdgeMesh
 
 public:
 
-    //- Runtime type information
     TypeName("featureEdgeMesh");
 
-    enum pointStatus
-    {
-        CONVEX,     // Fully convex point (w.r.t normals)
-        CONCAVE,    // Fully concave point
-        MIXED,      // A point surrounded by both convex and concave edges
-        NONFEATURE  // Not a feature point
-    };
-
-    enum edgeStatus
-    {
-        EXTERNAL,   // "Convex" edge
-        INTERNAL,   // "Concave" edge
-        FLAT,       // Neither concave or convex, on a flat surface
-        OPEN,       // i.e. only connected to one face
-        MULTIPLE,   // Multiply connected (connected to more than two faces)
-        NONE        // Not a classified feature edge (consistency with
-                    // surfaceFeatures)
-    };
-
-private:
-
-    // Static data
-
-        //- Angular closeness tolerance for treating normals as the same
-        static scalar cosNormalAngleTol_;
-
-        //- Index of the start of the convex feature points - static as 0
-        static label convexStart_;
-
-        //- Index of the start of the external feature edges - static as 0
-        static label externalStart_;
-
-
-    // Private data
-
-        //- Index of the start of the concave feature points
-        label concaveStart_;
-
-        //- Index of the start of the mixed type feature points
-        label mixedStart_;
-
-        //- Index of the start of the non-feature points
-        label nonFeatureStart_;
-
-        //- Index of the start of the internal feature edges
-        label internalStart_;
-
-        //- Index of the start of the flat feature edges
-        label flatStart_;
-
-        //- Index of the start of the open feature edges
-        label openStart_;
-
-        //- Index of the start of the multiply-connected feature edges
-        label multipleStart_;
-
-        //- Normals of the features, to be referred to by index by both feature
-        //  points and edges, unsorted
-        vectorField normals_;
-
-        //- Flat and open edges require the direction of the edge
-        vectorField edgeDirections_;
-
-        //- Indices of the normals that are adjacent to the feature edges
-        labelListList edgeNormals_;
-
-        //- Indices of the normals that are adjacent to the feature points
-        labelListList featurePointNormals_;
-
-        //- Feature edges which are on the boundary between regions
-        labelList regionEdges_;
-
-        //- Search tree for all edges
-        mutable autoPtr<indexedOctree<treeDataEdge> > edgeTree_;
-
-        //- Individual search trees for each type of edge
-        mutable PtrList<indexedOctree<treeDataEdge> > edgeTreesByType_;
-
-
-    // Private Member Functions
-
-        //- Classify the type of feature point.  Requires valid stored member
-        //  data for edges and normals.
-        pointStatus classifyFeaturePoint(label ptI) const;
-
-        //- Classify the type of feature edge.  Requires face centre 0 to face
-        //  centre 1 vector to distinguish internal from external
-        edgeStatus classifyEdge
-        (
-            const List<vector>& norms,
-            const labelList& edNorms,
-            const vector& fC0tofC1
-        ) const;
-
-
-public:
-
-    // Static data
-
-        //- Number of possible point types (i.e. number of slices)
-        static label nPointTypes;
-
-        //- Number of possible feature edge types (i.e. number of slices)
-        static label nEdgeTypes;
 
     // Constructors
 
         //- Construct (read) given an IOobject
         featureEdgeMesh(const IOobject&);
 
-        //- Construct as copy
+        //- Construct from featureEdgeMesh data
-        explicit featureEdgeMesh(const IOobject&, const featureEdgeMesh&);
-
-        //- Construct by transferring components (points, edges)
         featureEdgeMesh
         (
             const IOobject&,
-            const Xfer<pointField>&,
+            const pointField&,
-            const Xfer<edgeList>&
+            const edgeList&
         );
 
-        //- Construct (read) given surfaceFeatures, an objectRegistry and a
+        //- Construct as copy
-        //  fileName to write to.  Extracts, classifies and reorders the data
+        featureEdgeMesh(const IOobject&, const featureEdgeMesh&);
-        //  from surfaceFeatures.
-        featureEdgeMesh
-        (
-            const surfaceFeatures& sFeat,
-            const objectRegistry& obr,
-            const fileName& sFeatFileName
-        );
-
-        //- Construct from all components
-        featureEdgeMesh
-        (
-            const IOobject& io,
-            const pointField& pts,
-            const edgeList& eds,
-            label concaveStart,
-            label mixedStart,
-            label nonFeatureStart,
-            label internalStart,
-            label flatStart,
-            label openStart,
-            label multipleStart,
-            const vectorField& normals,
-            const vectorField& edgeDirections,
-            const labelListList& edgeNormals,
-            const labelListList& featurePointNormals,
-            const labelList& regionEdges
-        );
 
 
-    //- Destructor
+        //- ReadData function required for regIOobject read operation
-    ~featureEdgeMesh();
+        virtual bool readData(Istream&);
-
-
-    // Member Functions
-
-        // Find
-
-            //- Find nearest surface edge for the sample point.
-            void nearestFeatureEdge
-            (
-                const point& sample,
-                scalar searchDistSqr,
-                pointIndexHit& info
-            ) const;
-
-            //- Find nearest surface edge for each sample point.
-            void nearestFeatureEdge
-            (
-                const pointField& samples,
-                const scalarField& searchDistSqr,
-                List<pointIndexHit>& info
-            ) const;
-
-            //- Find the nearest point on each type of feature edge
-            void nearestFeatureEdgeByType
-            (
-                const point& sample,
-                const scalarField& searchDistSqr,
-                List<pointIndexHit>& info
-            ) const;
-
-        // Access
-
-            //- Return the index of the start of the convex feature points
-            inline label convexStart() const;
-
-            //- Return the index of the start of the concave feature points
-            inline label concaveStart() const;
-
-            //- Return the index of the start of the mixed type feature points
-            inline label mixedStart() const;
-
-            //- Return the index of the start of the non-feature points
-            inline label nonFeatureStart() const;
-
-            //- Return the index of the start of the external feature edges
-            inline label externalStart() const;
-
-            //- Return the index of the start of the internal feature edges
-            inline label internalStart() const;
-
-            //- Return the index of the start of the flat feature edges
-            inline label flatStart() const;
-
-            //- Return the index of the start of the open feature edges
-            inline label openStart() const;
-
-            //- Return the index of the start of the multiply-connected feature
-            //  edges
-            inline label multipleStart() const;
-
-            //- Return whether or not the point index is a feature point
-            inline bool featurePoint(label ptI) const;
-
-            //- Return the normals of the surfaces adjacent to the feature edges
-            //  and points
-            inline const vectorField& normals() const;
-
-            //- Return the edgeDirection vectors
-            inline const vectorField& edgeDirections() const;
-
-            //- Return the direction of edgeI, pointing away from ptI
-            inline vector edgeDirection(label edgeI, label ptI) const;
-
-            //- Return the indices of the normals that are adjacent to the
-            //  feature edges
-            inline const labelListList& edgeNormals() const;
-
-            //- Return the normal vectors for a given set of normal indices
-            inline vectorField edgeNormals(const labelList& edgeNormIs) const;
-
-            //- Return the normal vectors for a given edge
-            inline vectorField edgeNormals(label edgeI) const;
-
-            //- Return the indices of the normals that are adjacent to the
-            //  feature points
-            inline const labelListList& featurePointNormals() const;
-
-            //- Return the normal vectors for a given feature point
-            inline vectorField featurePointNormals(label ptI) const;
-
-            //- Return the feature edges which are on the boundary between
-            //  regions
-            inline const labelList& regionEdges() const;
-
-            //- Return the pointStatus of a specified point
-            inline pointStatus getPointStatus(label ptI) const;
-
-            //- Return the edgeStatus of a specified edge
-            inline edgeStatus getEdgeStatus(label edgeI) const;
-
-            //- Demand driven construction of octree for boundary edges
-            const indexedOctree<treeDataEdge>& edgeTree() const;
-
-            //- Demand driven construction of octree for boundary edges by type
-            const PtrList<indexedOctree<treeDataEdge> >&
-            edgeTreesByType() const;
-
-
-        // Write
-
-            //- Write all components of the featureEdgeMesh as obj files
-            void writeObj(const fileName& prefix) const;
-
-            //- Give precedence to the regIOobject write
-            using regIOobject::write;
 
         //- WriteData function required for regIOobject write operation
         virtual bool writeData(Ostream&) const;
@@ -362,10 +92,6 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#include "featureEdgeMeshI.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 #endif
 
 // ************************************************************************* //

src/finiteVolume/cfdTools/general/adjustPhi/adjustPhi.C

@@ -40,7 +40,8 @@ bool Foam::adjustPhi
 {
     if (p.needReference())
     {
-        p.boundaryField().updateCoeffs();
+        // p coefficients should not be updated here
+        // p.boundaryField().updateCoeffs();
 
         scalar massIn = 0.0;
         scalar fixedMassOut = 0.0;

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -45,6 +45,7 @@ Description
 #include "layerParameters.H"
 #include "combineFaces.H"
 #include "IOmanip.H"
+#include "globalIndex.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -888,68 +889,6 @@ void Foam::autoLayerDriver::handleNonManifolds
     Info<< "Outside of local patch is multiply connected across edges or"
         << " points at " << nNonManif << " points." << endl;
 
-
-    const labelList& meshPoints = pp.meshPoints();
-
-
-    // 2. Parallel check
-    // For now disable extrusion at any shared edges.
-    const labelHashSet sharedEdgeSet(mesh.globalData().sharedEdgeLabels());
-
-    forAll(pp.edges(), edgeI)
-    {
-        if (sharedEdgeSet.found(meshEdges[edgeI]))
-        {
-            const edge& e = mesh.edges()[meshEdges[edgeI]];
-
-            Pout<< "Disabling extrusion at edge "
-                << mesh.points()[e[0]] << mesh.points()[e[1]]
-                << " since it is non-manifold across coupled patches."
-                << endl;
-
-            nonManifoldPoints.insert(e[0]);
-            nonManifoldPoints.insert(e[1]);
-        }
-    }
-
-    // 3b. extrusion can produce multiple faces between 2 cells
-    // across processor boundary
-    // This occurs when a coupled face shares more than 1 edge with a
-    // non-coupled boundary face.
-    // This is now correctly handled by addPatchCellLayer in that it
-    // extrudes a single face from the stringed up edges.
-
-
-    nNonManif = returnReduce(nonManifoldPoints.size(), sumOp<label>());
-
-    if (nNonManif > 0)
-    {
-        Info<< "Outside of patches is multiply connected across edges or"
-            << " points at " << nNonManif << " points." << nl
-            << "Writing " << nNonManif
-            << " points where this happens to pointSet "
-            << nonManifoldPoints.name() << nl
-            << "and setting layer thickness to zero on these points."
-            << endl;
-
-        nonManifoldPoints.instance() = mesh.time().timeName();
-        nonManifoldPoints.write();
-
-        forAll(meshPoints, patchPointI)
-        {
-            if (nonManifoldPoints.found(meshPoints[patchPointI]))
-            {
-                unmarkExtrusion
-                (
-                    patchPointI,
-                    patchDisp,
-                    patchNLayers,
-                    extrudeStatus
-                );
-            }
-        }
-    }
-
     Info<< "Set displacement to zero for all " << nNonManif
         << " non-manifold points" << endl;
 }
@@ -1348,7 +1287,6 @@ void Foam::autoLayerDriver::setNumLayers
 }
 
 
-// Grow no-extrusion layer
 void Foam::autoLayerDriver::growNoExtrusion
 (
     const indirectPrimitivePatch& pp,
@@ -1439,6 +1377,103 @@ void Foam::autoLayerDriver::growNoExtrusion
 }
 
 
+void Foam::autoLayerDriver::determineSidePatches
+(
+    const globalIndex& globalFaces,
+    const labelListList& edgeGlobalFaces,
+    const indirectPrimitivePatch& pp,
+
+    labelList& sidePatchID
+)
+{
+    // Sometimes edges-to-be-extruded are on more than 2 processors.
+    // Work out which 2 hold the faces to be extruded and thus which procpatch
+    // the side-face should be in. As an additional complication this might
+    // mean that 2 procesors that were only edge-connected now suddenly need
+    // to become face-connected i.e. have a processor patch between them.
+
+    fvMesh& mesh = meshRefiner_.mesh();
+
+    // Determine sidePatchID. Any additional processor boundary gets added to
+    // patchToNbrProc,nbrProcToPatch and nPatches gets set to the new number
+    // of patches.
+    label nPatches;
+    Map<label> nbrProcToPatch;
+    Map<label> patchToNbrProc;
+    addPatchCellLayer::calcSidePatch
+    (
+        mesh,
+        globalFaces,
+        edgeGlobalFaces,
+        pp,
+
+        sidePatchID,
+        nPatches,
+        nbrProcToPatch,
+        patchToNbrProc
+    );
+
+    label nOldPatches = mesh.boundaryMesh().size();
+    label nAdded = returnReduce(nPatches-nOldPatches, sumOp<label>());
+    Info<< nl << "Adding in total " << nAdded/2 << " inter-processor patches to"
+        << " handle extrusion of non-manifold processor boundaries."
+        << endl;
+
+    if (nAdded > 0)
+    {
+        // We might not add patches in same order as in patchToNbrProc
+        // so prepare to renumber sidePatchID
+        Map<label> wantedToAddedPatch;
+
+        for (label patchI = nOldPatches; patchI < nPatches; patchI++)
+        {
+            label nbrProcI = patchToNbrProc[patchI];
+            word name =
+                    "procBoundary"
+                  + Foam::name(Pstream::myProcNo())
+                  + "to"
+                  + Foam::name(nbrProcI);
+
+            dictionary patchDict;
+            patchDict.add("type", processorPolyPatch::typeName);
+            patchDict.add("myProcNo", Pstream::myProcNo());
+            patchDict.add("neighbProcNo", nbrProcI);
+            patchDict.add("nFaces", 0);
+            patchDict.add("startFace", mesh.nFaces());
+
+            Pout<< "Adding patch " << patchI
+                << " name:" << name
+                << " between " << Pstream::myProcNo()
+                << " and " << nbrProcI
+                << endl;
+
+            label procPatchI = meshRefiner_.appendPatch
+            (
+                mesh,
+                mesh.boundaryMesh().size(), // new patch index
+                name,
+                patchDict
+            );
+            wantedToAddedPatch.insert(patchI, procPatchI);
+        }
+
+        // Renumber sidePatchID
+        forAll(sidePatchID, i)
+        {
+            label patchI = sidePatchID[i];
+            Map<label>::const_iterator fnd = wantedToAddedPatch.find(patchI);
+            if (fnd != wantedToAddedPatch.end())
+            {
+                sidePatchID[i] = fnd();
+            }
+        }
+
+        mesh.clearOut();
+        const_cast<polyBoundaryMesh&>(mesh.boundaryMesh()).updateMesh();
+    }
+}
+
+
 void Foam::autoLayerDriver::calculateLayerThickness
 (
     const indirectPrimitivePatch& pp,
@@ -2553,6 +2588,37 @@ void Foam::autoLayerDriver::addLayers
         )
     );
 
+
+    // Global face indices engine
+    const globalIndex globalFaces(mesh.nFaces());
+
+    // Determine extrudePatch.edgeFaces in global numbering (so across
+    // coupled patches). This is used only to string up edges between coupled
+    // faces (all edges between same (global)face indices get extruded).
+    labelListList edgeGlobalFaces
+    (
+        addPatchCellLayer::globalEdgeFaces
+        (
+            mesh,
+            globalFaces,
+            pp
+        )
+    );
+
+    // Determine patches for extruded boundary edges. Calculates if any
+    // additional processor patches need to be constructed.
+
+    labelList sidePatchID;
+    determineSidePatches
+    (
+        globalFaces,
+        edgeGlobalFaces,
+        pp,
+
+        sidePatchID
+    );
+
+
     // Construct iterative mesh mover.
     Info<< "Constructing mesh displacer ..." << endl;
 
@@ -3038,8 +3104,12 @@ void Foam::autoLayerDriver::addLayers
         // Not add layer if patchDisp is zero.
         addLayer.setRefinement
         (
+            globalFaces,
+            edgeGlobalFaces,
+
             invExpansionRatio,
             pp(),
+            sidePatchID,        // boundary patch for extruded boundary edges
             labelList(0),       // exposed patchIDs, not used for adding layers
             nPatchFaceLayers,   // layers per face
             nPatchPointLayers,  // layers per point

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -243,6 +243,16 @@ class autoLayerDriver
                     List<extrudeMode>& extrudeStatus
                 ) const;
 
+                //- See what patches boundaryedges should be extruded into
+                void determineSidePatches
+                (
+                    const globalIndex& globalFaces,
+                    const labelListList& edgeGlobalFaces,
+                    const indirectPrimitivePatch& pp,
+
+                    labelList& sidePatchID
+                );
+
                 //- Calculate pointwise wanted and minimum thickness.
                 //  thickness: wanted thickness
                 //  minthickness: when to give up and not extrude

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

@@ -1549,60 +1549,28 @@ void Foam::meshRefinement::checkCoupledFaceZones(const polyMesh& mesh)
 }
 
 
-// Adds patch if not yet there. Returns patchID.
+Foam::label Foam::meshRefinement::appendPatch
-Foam::label Foam::meshRefinement::addPatch
 (
     fvMesh& mesh,
+    const label insertPatchI,
     const word& patchName,
-    const dictionary& patchInfo
+    const dictionary& patchDict
 )
 {
-    polyBoundaryMesh& polyPatches =
-        const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
-
-    const label patchI = polyPatches.findPatchID(patchName);
-    if (patchI != -1)
-    {
-        // Already there
-        return patchI;
-    }
-
-
-    label insertPatchI = polyPatches.size();
-    label startFaceI = mesh.nFaces();
-
-    forAll(polyPatches, patchI)
-    {
-        const polyPatch& pp = polyPatches[patchI];
-
-        if (isA<processorPolyPatch>(pp))
-        {
-            insertPatchI = patchI;
-            startFaceI = pp.start();
-            break;
-        }
-    }
-
-
-    // Below is all quite a hack. Feel free to change once there is a better
-    // mechanism to insert and reorder patches.
-
     // Clear local fields and e.g. polyMesh parallelInfo.
     mesh.clearOut();
 
-    label sz = polyPatches.size();
+    polyBoundaryMesh& polyPatches =
-
+        const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
     fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary());
 
-    dictionary patchDict(patchInfo);
+    label patchI = polyPatches.size();
-    patchDict.set("nFaces", 0);
-    patchDict.set("startFace", startFaceI);
 
     // Add polyPatch at the end
-    polyPatches.setSize(sz+1);
+    polyPatches.setSize(patchI+1);
     polyPatches.set
     (
-        sz,
+        patchI,
         polyPatch::New
         (
             patchName,
@@ -1611,13 +1579,13 @@ Foam::label Foam::meshRefinement::addPatch
             polyPatches
         )
     );
-    fvPatches.setSize(sz+1);
+    fvPatches.setSize(patchI+1);
     fvPatches.set
     (
-        sz,
+        patchI,
         fvPatch::New
         (
-            polyPatches[sz],  // point to newly added polyPatch
+            polyPatches[patchI],  // point to newly added polyPatch
             mesh.boundary()
         )
     );
@@ -1675,21 +1643,69 @@ Foam::label Foam::meshRefinement::addPatch
         mesh,
         calculatedFvPatchField<tensor>::typeName
     );
+    return patchI;
+}
+
+
+// Adds patch if not yet there. Returns patchID.
+Foam::label Foam::meshRefinement::addPatch
+(
+    fvMesh& mesh,
+    const word& patchName,
+    const dictionary& patchInfo
+)
+{
+    polyBoundaryMesh& polyPatches =
+        const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
+    fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary());
+
+    const label patchI = polyPatches.findPatchID(patchName);
+    if (patchI != -1)
+    {
+        // Already there
+        return patchI;
+    }
+
+
+    label insertPatchI = polyPatches.size();
+    label startFaceI = mesh.nFaces();
+
+    forAll(polyPatches, patchI)
+    {
+        const polyPatch& pp = polyPatches[patchI];
+
+        if (isA<processorPolyPatch>(pp))
+        {
+            insertPatchI = patchI;
+            startFaceI = pp.start();
+            break;
+        }
+    }
+
+    dictionary patchDict(patchInfo);
+    patchDict.set("nFaces", 0);
+    patchDict.set("startFace", startFaceI);
+
+    // Below is all quite a hack. Feel free to change once there is a better
+    // mechanism to insert and reorder patches.
+
+    label addedPatchI = appendPatch(mesh, insertPatchI, patchName, patchDict);
+
 
     // Create reordering list
     // patches before insert position stay as is
-    labelList oldToNew(sz+1);
+    labelList oldToNew(addedPatchI+1);
     for (label i = 0; i < insertPatchI; i++)
     {
         oldToNew[i] = i;
     }
     // patches after insert position move one up
-    for (label i = insertPatchI; i < sz; i++)
+    for (label i = insertPatchI; i < addedPatchI; i++)
     {
         oldToNew[i] = i+1;
     }
     // appended patch gets moved to insert position
-    oldToNew[sz] = insertPatchI;
+    oldToNew[addedPatchI] = insertPatchI;
 
     // Shuffle into place
     polyPatches.reorder(oldToNew);

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -747,8 +747,17 @@ public:
 
         // Other topo changes
 
+            //- Helper:append patch to end of mesh.
+            static label appendPatch
+            (
+                fvMesh&,
+                const label insertPatchI,
+                const word&,
+                const dictionary&
+            );
+
             //- Helper:add patch to mesh. Update all registered fields.
-            //  Use addMeshedPatch to add patches originating from surfaces.
+            //  Used by addMeshedPatch to add patches originating from surfaces.
             static label addPatch(fvMesh&, const word& name, const dictionary&);
 
             //- Add patch originating from meshing. Update meshedPatches_.

src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -1007,6 +1007,82 @@ void Foam::refinementSurfaces::findNearestRegion
 }
 
 
+void Foam::refinementSurfaces::findNearestRegion
+(
+    const labelList& surfacesToTest,
+    const pointField& samples,
+    const scalarField& nearestDistSqr,
+    labelList& hitSurface,
+    List<pointIndexHit>& hitInfo,
+    labelList& hitRegion,
+    vectorField& hitNormal
+) const
+{
+    labelList geometries(UIndirectList<label>(surfaces_, surfacesToTest));
+
+    // Do the tests. Note that findNearest returns index in geometries.
+    searchableSurfacesQueries::findNearest
+    (
+        allGeometry_,
+        geometries,
+        samples,
+        nearestDistSqr,
+        hitSurface,
+        hitInfo
+    );
+
+    // Rework the hitSurface to be surface (i.e. index into surfaces_)
+    forAll(hitSurface, pointI)
+    {
+        if (hitSurface[pointI] != -1)
+        {
+            hitSurface[pointI] = surfacesToTest[hitSurface[pointI]];
+        }
+    }
+
+    // Collect the region
+    hitRegion.setSize(hitSurface.size());
+    hitRegion = -1;
+    hitNormal.setSize(hitSurface.size());
+    hitNormal = vector::zero;
+
+    forAll(surfacesToTest, i)
+    {
+        label surfI = surfacesToTest[i];
+
+        // Collect hits for surfI
+        const labelList localIndices(findIndices(hitSurface, surfI));
+
+        List<pointIndexHit> localHits
+        (
+            UIndirectList<pointIndexHit>
+            (
+                hitInfo,
+                localIndices
+            )
+        );
+
+        // Region
+        labelList localRegion;
+        allGeometry_[surfaces_[surfI]].getRegion(localHits, localRegion);
+
+        forAll(localIndices, i)
+        {
+            hitRegion[localIndices[i]] = localRegion[i];
+        }
+
+        // Normal
+        vectorField localNormal;
+        allGeometry_[surfaces_[surfI]].getNormal(localHits, localNormal);
+
+        forAll(localIndices, i)
+        {
+            hitNormal[localIndices[i]] = localNormal[i];
+        }
+    }
+}
+
+
 //// Find intersection with max of edge. Return -1 or the surface
 //// with the highest maxLevel above currentLevel
 //Foam::label Foam::refinementSurfaces::findHighestIntersection

src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -327,6 +327,19 @@ public:
                 labelList& hitRegion
             ) const;
 
+            //- Find nearest point on surfaces. Return surface, region and
+            //  normal on surface (so not global surface)
+            void findNearestRegion
+            (
+                const labelList& surfacesToTest,
+                const pointField& samples,
+                const scalarField& nearestDistSqr,
+                labelList& hitSurface,
+                List<pointIndexHit>& hitInfo,
+                labelList& hitRegion,
+                vectorField& hitNormal
+            ) const;
+
             //- Detect if a point is 'inside' (closed) surfaces.
             //  Returns -1 if not, returns first surface it is.
             void findInside

src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -398,15 +398,24 @@ void Foam::decompositionMethod::calcCellCells
     // Create global cell numbers
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    globalIndex globalCells(mesh.nCells());
+    label nAgglom;
+    if (agglom.size())
+    {
+        nAgglom = max(agglom)+1;
+    }
+    else
+    {
+        nAgglom = 0;
+    }
+    globalIndex globalAgglom(nAgglom);
 
     const labelList& faceOwner = mesh.faceOwner();
     const labelList& faceNeighbour = mesh.faceNeighbour();
     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 
 
-    // Get renumbered owner on other side of coupled faces
+    // Get agglomerate owner on other side of coupled faces
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
     labelList globalNeighbour(mesh.nFaces()-mesh.nInternalFaces());
 
@@ -421,10 +430,13 @@ void Foam::decompositionMethod::calcCellCells
 
             forAll(pp, i)
             {
-                globalNeighbour[bFaceI++] = globalCells.toGlobal
+                globalNeighbour[bFaceI] = globalAgglom.toGlobal
                 (
-                    faceOwner[faceI++]
+                    agglom[faceOwner[faceI]]
                 );
+
+                bFaceI++;
+                faceI++;
             }
         }
     }
@@ -448,8 +460,13 @@ void Foam::decompositionMethod::calcCellCells
         nFacesPerCell[nei]++;
     }
 
-    // Handle coupled faces
+    // Handle coupled faces. In case of agglomeration you might end up
-    HashSet<edge, Hash<edge> > cellPair(mesh.nFaces()-mesh.nInternalFaces());
+    // with multiple connections ('faces') between the same two agglomerations.
+    // This is illegal so mark.
+    HashSet<labelPair, Hash<labelPair> > cellPair
+    (
+        mesh.nFaces()-mesh.nInternalFaces()
+    );
 
     forAll(patches, patchI)
     {
@@ -464,7 +481,7 @@ void Foam::decompositionMethod::calcCellCells
             {
                 label own = agglom[faceOwner[faceI]];
                 label globalNei = globalNeighbour[bFaceI];
-                if (cellPair.insert(edge(own, globalNei)))
+                if (cellPair.insert(labelPair(own, globalNei)))
                 {
                     nFacesPerCell[own]++;
                 }
@@ -491,8 +508,8 @@ void Foam::decompositionMethod::calcCellCells
         label own = agglom[faceOwner[faceI]];
         label nei = agglom[faceNeighbour[faceI]];
 
-        m[offsets[own] + nFacesPerCell[own]++] = globalCells.toGlobal(nei);
+        m[offsets[own] + nFacesPerCell[own]++] = globalAgglom.toGlobal(nei);
-        m[offsets[nei] + nFacesPerCell[nei]++] = globalCells.toGlobal(own);
+        m[offsets[nei] + nFacesPerCell[nei]++] = globalAgglom.toGlobal(own);
     }
 
     // For boundary faces is offsetted coupled neighbour
@@ -510,7 +527,7 @@ void Foam::decompositionMethod::calcCellCells
             {
                 label own = agglom[faceOwner[faceI]];
                 label globalNei = globalNeighbour[bFaceI];
-                if (cellPair.insert(edge(own, globalNei)))
+                if (cellPair.insert(labelPair(own, globalNei)))
                 {
                     m[offsets[own] + nFacesPerCell[own]++] = globalNei;
                 }

src/parallel/decompose/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,6 +26,7 @@ License
 #include "simpleGeomDecomp.H"
 #include "addToRunTimeSelectionTable.H"
 #include "SortableList.H"
+#include "globalIndex.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -53,7 +54,7 @@ void Foam::simpleGeomDecomp::assignToProcessorGroup
 (
     labelList& processorGroup,
     const label nProcGroup
-)
+) const
 {
     label jump = processorGroup.size()/nProcGroup;
     label jumpb = jump + 1;
@@ -90,7 +91,7 @@ void Foam::simpleGeomDecomp::assignToProcessorGroup
     const labelList& indices,
     const scalarField& weights,
     const scalar summedWeights
-)
+) const
 {
     // This routine gets the sorted points.
     // Easiest to explain with an example.
@@ -126,7 +127,10 @@ void Foam::simpleGeomDecomp::assignToProcessorGroup
 }
 
 
-Foam::labelList Foam::simpleGeomDecomp::decompose(const pointField& points)
+Foam::labelList Foam::simpleGeomDecomp::decomposeOneProc
+(
+    const pointField& points
+) const
 {
     // construct a list for the final result
     labelList finalDecomp(points.size());
@@ -195,11 +199,11 @@ Foam::labelList Foam::simpleGeomDecomp::decompose(const pointField& points)
 }
 
 
-Foam::labelList Foam::simpleGeomDecomp::decompose
+Foam::labelList Foam::simpleGeomDecomp::decomposeOneProc
 (
     const pointField& points,
     const scalarField& weights
-)
+) const
 {
     // construct a list for the final result
     labelList finalDecomp(points.size());
@@ -300,4 +304,162 @@ Foam::simpleGeomDecomp::simpleGeomDecomp(const dictionary& decompositionDict)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+Foam::labelList Foam::simpleGeomDecomp::decompose
+(
+    const pointField& points
+)
+{
+    if (!Pstream::parRun())
+    {
+        return decomposeOneProc(points);
+    }
+    else
+    {
+        globalIndex globalNumbers(points.size());
+
+        // Collect all points on master
+        if (Pstream::master())
+        {
+            pointField allPoints(globalNumbers.size());
+
+            label nTotalPoints = 0;
+            // Master first
+            SubField<point>(allPoints, points.size()).assign(points);
+            nTotalPoints += points.size();
+
+            // Add slaves
+            for (int slave=1; slave<Pstream::nProcs(); slave++)
+            {
+                IPstream fromSlave(Pstream::scheduled, slave);
+                pointField nbrPoints(fromSlave);
+                SubField<point>
+                (
+                    allPoints,
+                    nbrPoints.size(),
+                    nTotalPoints
+                ).assign(nbrPoints);
+                nTotalPoints += nbrPoints.size();
+            }
+
+            // Decompose
+            labelList finalDecomp(decomposeOneProc(allPoints));
+
+            // Send back
+            for (int slave=1; slave<Pstream::nProcs(); slave++)
+            {
+                OPstream toSlave(Pstream::scheduled, slave);
+                toSlave << SubField<label>
+                (
+                    finalDecomp,
+                    globalNumbers.localSize(slave),
+                    globalNumbers.offset(slave)
+                );
+            }
+            // Get my own part
+            finalDecomp.setSize(points.size());
+
+            return finalDecomp;
+        }
+        else
+        {
+            // Send my points
+            {
+                OPstream toMaster(Pstream::scheduled, Pstream::masterNo());
+                toMaster<< points;
+            }
+
+            // Receive back decomposition
+            IPstream fromMaster(Pstream::scheduled, Pstream::masterNo());
+            labelList finalDecomp(fromMaster);
+
+            return finalDecomp;
+        }
+    }
+}
+
+
+Foam::labelList Foam::simpleGeomDecomp::decompose
+(
+    const pointField& points,
+    const scalarField& weights
+)
+{
+    if (!Pstream::parRun())
+    {
+        return decomposeOneProc(points, weights);
+    }
+    else
+    {
+        globalIndex globalNumbers(points.size());
+
+        // Collect all points on master
+        if (Pstream::master())
+        {
+            pointField allPoints(globalNumbers.size());
+            scalarField allWeights(allPoints.size());
+
+            label nTotalPoints = 0;
+            // Master first
+            SubField<point>(allPoints, points.size()).assign(points);
+            SubField<scalar>(allWeights, points.size()).assign(weights);
+            nTotalPoints += points.size();
+
+            // Add slaves
+            for (int slave=1; slave<Pstream::nProcs(); slave++)
+            {
+                IPstream fromSlave(Pstream::scheduled, slave);
+                pointField nbrPoints(fromSlave);
+                scalarField nbrWeights(fromSlave);
+                SubField<point>
+                (
+                    allPoints,
+                    nbrPoints.size(),
+                    nTotalPoints
+                ).assign(nbrPoints);
+                SubField<scalar>
+                (
+                    allWeights,
+                    nbrWeights.size(),
+                    nTotalPoints
+                ).assign(nbrWeights);
+                nTotalPoints += nbrPoints.size();
+            }
+
+            // Decompose
+            labelList finalDecomp(decomposeOneProc(allPoints, allWeights));
+
+            // Send back
+            for (int slave=1; slave<Pstream::nProcs(); slave++)
+            {
+                OPstream toSlave(Pstream::scheduled, slave);
+                toSlave << SubField<label>
+                (
+                    finalDecomp,
+                    globalNumbers.localSize(slave),
+                    globalNumbers.offset(slave)
+                );
+            }
+            // Get my own part
+            finalDecomp.setSize(points.size());
+
+            return finalDecomp;
+        }
+        else
+        {
+            // Send my points
+            {
+                OPstream toMaster(Pstream::scheduled, Pstream::masterNo());
+                toMaster<< points << weights;
+            }
+
+            // Receive back decomposition
+            IPstream fromMaster(Pstream::scheduled, Pstream::masterNo());
+            labelList finalDecomp(fromMaster);
+
+            return finalDecomp;
+        }
+    }
+}
+
+
 // ************************************************************************* //

src/parallel/decompose/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -49,7 +49,7 @@ class simpleGeomDecomp
 {
     // Private Member Functions
 
-        void assignToProcessorGroup(labelList& processorGroup, const label);
+        void assignToProcessorGroup(labelList&, const label) const;
 
         void assignToProcessorGroup
         (
@@ -58,7 +58,15 @@ class simpleGeomDecomp
             const labelList& indices,
             const scalarField& weights,
             const scalar summedWeights
-        );
+        ) const;
+
+        labelList decomposeOneProc(const pointField& points) const;
+
+        labelList decomposeOneProc
+        (
+            const pointField& points,
+            const scalarField& weights
+        ) const;
 
         //- Disallow default bitwise copy construct and assignment
         void operator=(const simpleGeomDecomp&);
@@ -86,9 +94,9 @@ public:
 
         virtual bool parallelAware() const
         {
-            // simpleDecomp does not attempt to do anything across proc
+            // simpleDecomp sends all points to the master which does
-            // boundaries
+            // the decomposition.
-            return false;
+            return true;
         }
 
         virtual labelList decompose(const pointField&);

tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/Allrun

@@ -16,11 +16,16 @@ setSet -batch system/sets.setSet > log.setSet1 2>&1
 setsToZones -noFlipMap > log.setsToZones 2>&1
 
 # create the first cyclic - lhs of porous zone
-unset FOAM_SETNAN
+# Note that we don't know what value to give these patches-out-of-nothing so
+# - use binary writing to avoid 'nan'
+# - use setFields to set values
+unset FOAM_SIGFPE
 createBaffles cycLeft '(cycLeft_half0 cycLeft_half1)' -overwrite > log.createBaffles1 2>&1
 
 # create the second cyclic - rhs of porous zone
 createBaffles cycRight '(cycRight_half0 cycRight_half1)' -overwrite > log.createBaffles2 2>&1
+# Initialise newly created patchFields to 0
+changeDictionary
 
 runApplication $application
 

tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/system/changeDictionaryDict

@@ -0,0 +1,197 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.7.1                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      changeDictionaryDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dictionaryReplacement
+{
+    alphat
+    {
+        boundaryField
+        {
+            cycLeft
+            {
+                type            cyclic;
+                value           uniform 0;
+            }
+            cycRight
+            {
+                type            cyclic;
+                value           uniform 0;
+            }
+        }
+    }
+    epsilon
+    {
+        boundaryField
+        {
+            cycLeft
+            {
+                type            cyclic;
+                value           uniform 0;
+            }
+            cycRight
+            {
+                type            cyclic;
+                value           uniform 0;
+            }
+        }
+    }
+    G
+    {
+        boundaryField
+        {
+            cycLeft
+            {
+                type            cyclic;
+                value           uniform 0;
+            }
+            cycRight
+            {
+                type            cyclic;
+                value           uniform 0;
+            }
+        }
+    }
+    H2O
+    {
+        boundaryField
+        {
+            cycLeft
+            {
+                type            cyclic;
+                value           uniform 0;
+            }
+            cycRight
+            {
+                type            cyclic;
+                value           uniform 0;
+            }
+        }
+    }
+    k
+    {
+        boundaryField
+        {
+            cycLeft
+            {
+                type            cyclic;
+                value           uniform 0;
+            }
+            cycRight
+            {
+                type            cyclic;
+                value           uniform 0;
+            }
+        }
+    }
+    mut
+    {
+        boundaryField
+        {
+            cycLeft
+            {
+                type            cyclic;
+                value           uniform 0;
+            }
+            cycRight
+            {
+                type            cyclic;
+                value           uniform 0;
+            }
+        }
+    }
+    N2
+    {
+        boundaryField
+        {
+            cycLeft
+            {
+                type            cyclic;
+                value           uniform 0;
+            }
+            cycRight
+            {
+                type            cyclic;
+                value           uniform 0;
+            }
+        }
+    }
+    O2
+    {
+        boundaryField
+        {
+            cycLeft
+            {
+                type            cyclic;
+                value           uniform 0;
+            }
+            cycRight
+            {
+                type            cyclic;
+                value           uniform 0;
+            }
+        }
+    }
+    p
+    {
+        boundaryField
+        {
+            cycLeft
+            {
+                type            cyclic;
+                value           uniform 0;
+            }
+            cycRight
+            {
+                type            cyclic;
+                value           uniform 0;
+            }
+        }
+    }
+    T
+    {
+        boundaryField
+        {
+            cycLeft
+            {
+                type            cyclic;
+                value           uniform 0;
+            }
+            cycRight
+            {
+                type            cyclic;
+                value           uniform 0;
+            }
+        }
+    }
+    U
+    {
+        boundaryField
+        {
+            cycLeft
+            {
+                type            cyclic;
+                value           uniform (0 0 0);
+            }
+            cycRight
+            {
+                type            cyclic;
+                value           uniform (0 0 0);
+            }
+        }
+    }
+}
+
+// ************************************************************************* //

tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/system/controlDict

@@ -33,7 +33,7 @@ writeInterval   0.1;
 
 purgeWrite      0;
 
-writeFormat     ascii;
+writeFormat     binary;
 
 writePrecision  10;
 

wmake/wmakePrintBuild

@@ -98,14 +98,20 @@ done
 #
 build="$WM_PROJECT_DIR/.build"
 
+unset oldVersion oldPackage
+
+#
 # retrieve old values from the $WM_PROJECT_DIR/.build cache, stored as
 #     version [packager]
+#
+getOldValues()
+{
     set -- $(tail -1 $build 2>/dev/null)
-
+    oldVersion="$1"
-oldVersion="$1"; shift
+    [ "$#" -gt 0 ] && shift
     oldPackage="$@"
     [ "${oldPackage:-none}" = none ] && unset oldPackage
-
+}
 
 #
 # printTag - output the build tag
@@ -148,6 +154,9 @@ else
 fi
 
 
+# retrieve old values
+getOldValues
+
 #
 # update persistent build tag if possible
 #