From 5c97cd1d39182c4be7478137785a3cd32c6a1511 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs@hunt.opencfd.co.uk>
Date: Thu, 19 Feb 2009 13:08:07 +0000
Subject: [PATCH 01/11] added comment

---
 .../iglooWithFridges/system/snappyHexMeshDict | 192 +++++++++++++++++-
 .../motorBike/system/snappyHexMeshDict        | 191 ++++++++++++++++-
 2 files changed, 377 insertions(+), 6 deletions(-)

diff --git a/tutorials/mesh/snappyHexMesh/iglooWithFridges/system/snappyHexMeshDict b/tutorials/mesh/snappyHexMesh/iglooWithFridges/system/snappyHexMeshDict
index 2efe6901aa..d185a59ce9 100644
--- a/tutorials/mesh/snappyHexMesh/iglooWithFridges/system/snappyHexMeshDict
+++ b/tutorials/mesh/snappyHexMesh/iglooWithFridges/system/snappyHexMeshDict
@@ -15,12 +15,18 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+// Which of the steps to run
 castellatedMesh true;
-
 snap            true;
-
 addLayers       true;
 
+
+// Geometry. Definition of all surfaces. All surfaces are of class
+// searchableSurface.
+// Surfaces are used
+// - to specify refinement for any mesh cell intersecting it
+// - to specify refinement for any mesh cell inside/outside/near
+// - to 'snap' the mesh boundary to the surface
 geometry
 {
     fridgeA
@@ -45,17 +51,68 @@ geometry
     }
 }
 
+
+
+// Settings for the castellatedMesh generation.
 castellatedMeshControls
 {
+
+    // Refinement parameters
+    // ~~~~~~~~~~~~~~~~~~~~~
+
+    // While refining maximum number of cells per processor. This is basically
+    // the number of cells that fit on a processor. If you choose this too small
+    // it will do just more refinement iterations to obtain a similar mesh.
     maxLocalCells   1000000;
+
+    // Overall cell limit (approximately). Refinement will stop immediately
+    // upon reaching this number so a refinement level might not complete.
+    // Note that this is the number of cells before removing the part which
+    // is not 'visible' from the keepPoint. The final number of cells might
+    // actually be a lot less.
     maxGlobalCells  2000000;
+
+    // The surface refinement loop might spend lots of iterations refining just a
+    // few cells. This setting will cause refinement to stop if <= minimumRefine
+    // are selected for refinement. Note: it will at least do one iteration
+    // (unless the number of cells to refine is 0)
     minRefinementCells 0;
+
+    // Number of buffer layers between different levels.
+    // 1 means normal 2:1 refinement restriction, larger means slower
+    // refinement.
     nCellsBetweenLevels 1;
-    features        ( );
+
+
+
+    // Explicit feature edge refinement
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Specifies a level for any cell intersected by its edges.
+    // This is a featureEdgeMesh, read from constant/triSurface for now.
+    features
+    (
+        //{
+        //    file "someLine.eMesh";
+        //    level 2;
+        //}
+    );
+
+
+
+    // Surface based refinement
+    // ~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Specifies two levels for every surface. The first is the minimum level,
+    // every cell intersecting a surface gets refined up to the minimum level.
+    // The second level is the maximum level. Cells that 'see' multiple
+    // intersections where the intersections make an
+    // angle > resolveFeatureAngle get refined up to the maximum level.
     refinementSurfaces
     {
         fridgeA
         {
+            // Surface-wise min and max refinement level
             level           ( 2 2 );
         }
 
@@ -71,23 +128,64 @@ castellatedMeshControls
     }
 
     resolveFeatureAngle 60;
+
+    // Region-wise refinement
+    // ~~~~~~~~~~~~~~~~~~~~~~
+
+    // Specifies refinement level for cells in relation to a surface. One of
+    // three modes
+    // - distance. 'levels' specifies per distance to the surface the
+    //   wanted refinement level. The distances need to be specified in
+    //   descending order.
+    // - inside. 'levels' is only one entry and only the level is used. All
+    //   cells inside the surface get refined up to the level. The surface
+    //   needs to be closed for this to be possible.
+    // - outside. Same but cells outside.
+
     refinementRegions
     {
     }
 
+
+    // Mesh selection
+    // ~~~~~~~~~~~~~~
+
+    // After refinement patches get added for all refinementSurfaces and
+    // all cells intersecting the surfaces get put into these patches. The
+    // section reachable from the locationInMesh is kept.
+    // NOTE: This point should never be on a face, always inside a cell, even
+    // after refinement.
     locationInMesh  ( 3 0.28 0.43 );
 }
 
+
+
+// Settings for the snapping.
 snapControls
 {
+    //- Number of patch smoothing iterations before finding correspondence
+    //  to surface
     nSmoothPatch    3;
+
+    //- Relative distance for points to be attracted by surface feature point
+    //  or edge. True distance is this factor times local
+    //  maximum edge length.
     tolerance       4;
+
+    //- Number of mesh displacement relaxation iterations.
     nSolveIter      30;
+
+    //- Maximum number of snapping relaxation iterations. Should stop
+    //  before upon reaching a correct mesh.
     nRelaxIter      5;
 }
 
+
+
+// Settings for the layer addition.
 addLayersControls
 {
+    // Per final patch (so not geometry!) the layer information
     layers
     {
         fridgeA_region0
@@ -106,41 +204,129 @@ addLayersControls
         }
     }
 
+
+    // Expansion factor for layer mesh
     expansionRatio  1;
+
+    //- Wanted thickness of final added cell layer. If multiple layers
+    //  is the
+    //  thickness of the layer furthest away from the wall.
+    //  Relative to undistorted size of cell outside layer.
     finalLayerRatio 0.5;
+
+    //- Minimum thickness of cell layer. If for any reason layer
+    //  cannot be above minThickness do not add layer.
+    //  Relative to undistorted size of cell outside layer.
     minThickness    0.25;
+
+    //- If points get not extruded do nGrow layers of connected faces that are
+    //  also not grown. This helps convergence of the layer addition process
+    //  close to features.
     nGrow           0;
+
+
+    // Advanced settings
+
+    //- When not to extrude surface. 0 is flat surface, 90 is when two faces
+    //  make straight angle.
     featureAngle    60;
+
+    //- Maximum number of snapping relaxation iterations. Should stop
+    //  before upon reaching a correct mesh.
     nRelaxIter      5;
+
+    // Number of smoothing iterations of surface normals
     nSmoothSurfaceNormals 1;
+
+    // Number of smoothing iterations of interior mesh movement direction
     nSmoothNormals  3;
+
+    // Smooth layer thickness over surface patches
     nSmoothThickness 10;
+
+    // Stop layer growth on highly warped cells
     maxFaceThicknessRatio 0.5;
+
+    // Reduce layer growth where ratio thickness to medial
+    // distance is large
     maxThicknessToMedialRatio 0.3;
+
+    // Angle used to pick up medial axis points
     minMedianAxisAngle 130;
+
+    // Create buffer region for new layer terminations
     nBufferCellsNoExtrude 0;
 }
 
+
+
+// Generic mesh quality settings. At any undoable phase these determine
+// where to undo.
 meshQualityControls
 {
+    //- Maximum non-orthogonality allowed. Set to 180 to disable.
     maxNonOrtho     65;
+
+    //- Max skewness allowed. Set to <0 to disable.
     maxBoundarySkewness 20;
     maxInternalSkewness 4;
+
+    //- Max concaveness allowed. Is angle (in degrees) below which concavity
+    //  is allowed. 0 is straight face, <0 would be convex face.
+    //  Set to 180 to disable.
     maxConcave      80;
+
+    //- Minimum projected area v.s. actual area. Set to -1 to disable.
     minFlatness     0.5;
+
+    //- Minimum pyramid volume. Is absolute volume of cell pyramid.
+    //  Set to a sensible fraction of the smallest cell volume expected.
+    //  Set to very negative number (e.g. -1E30) to disable.
     minVol          1e-13;
+
+    //- Minimum face area. Set to <0 to disable.
     minArea         -1;
+
+    //- Minimum face twist. Set to <-1 to disable. dot product of face normal
+    //- and face centre triangles normal
     minTwist        0.05;
+
+    //- minimum normalised cell determinant
+    //- 1 = hex, <= 0 = folded or flattened illegal cell
     minDeterminant  0.001;
+
+    //- minFaceWeight (0 -> 0.5)
     minFaceWeight   0.05;
+
+    //- minVolRatio (0 -> 1)
     minVolRatio     0.01;
+
+    //must be >0 for Fluent compatibility
     minTriangleTwist -1;
+
+
+    // Advanced
+
+    //- Number of error distribution iterations
     nSmoothScale    4;
+    //- amount to scale back displacement at error points
     errorReduction  0.75;
 }
 
+
+
+// Advanced
+
+// Flags for optional output
+// 0 : only write final meshes
+// 1 : write intermediate meshes
+// 2 : write volScalarField with cellLevel for postprocessing
+// 4 : write current intersections as .obj files
 debug           0;
 
+
+// Merge tolerance. Is fraction of overall bounding box of initial mesh.
+// Note: the write tolerance needs to be higher than this.
 mergeTolerance  1e-06;
 
 
diff --git a/tutorials/mesh/snappyHexMesh/motorBike/system/snappyHexMeshDict b/tutorials/mesh/snappyHexMesh/motorBike/system/snappyHexMeshDict
index fa094e1b03..3f0ff0417e 100644
--- a/tutorials/mesh/snappyHexMesh/motorBike/system/snappyHexMeshDict
+++ b/tutorials/mesh/snappyHexMesh/motorBike/system/snappyHexMeshDict
@@ -13,14 +13,21 @@ FoamFile
     location    "system";
     object      snappyHexMeshDict;
 }
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+// Which of the steps to run
 castellatedMesh true;
-
 snap            true;
-
 addLayers       true;
 
+
+// Geometry. Definition of all surfaces. All surfaces are of class
+// searchableSurface.
+// Surfaces are used
+// - to specify refinement for any mesh cell intersecting it
+// - to specify refinement for any mesh cell inside/outside/near
+// - to 'snap' the mesh boundary to the surface
 geometry
 {
     motorBike.stl
@@ -37,22 +44,89 @@ geometry
     }
 }
 
+
+// Settings for the castellatedMesh generation.
 castellatedMeshControls
 {
+
+    // Refinement parameters
+    // ~~~~~~~~~~~~~~~~~~~~~
+
+    // While refining maximum number of cells per processor. This is basically
+    // the number of cells that fit on a processor. If you choose this too small
+    // it will do just more refinement iterations to obtain a similar mesh.
     maxLocalCells   1000000;
+
+
+    // Overall cell limit (approximately). Refinement will stop immediately
+    // upon reaching this number so a refinement level might not complete.
+    // Note that this is the number of cells before removing the part which
+    // is not 'visible' from the keepPoint. The final number of cells might
+    // actually be a lot less.
     maxGlobalCells  2000000;
+
+    // The surface refinement loop might spend lots of iterations refining just a
+    // few cells. This setting will cause refinement to stop if <= minimumRefine
+    // are selected for refinement. Note: it will at least do one iteration
+    // (unless the number of cells to refine is 0)
     minRefinementCells 10;
+
+    // Number of buffer layers between different levels.
+    // 1 means normal 2:1 refinement restriction, larger means slower
+    // refinement.
     nCellsBetweenLevels 2;
-    features        ( );
+
+
+
+    // Explicit feature edge refinement
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Specifies a level for any cell intersected by its edges.
+    // This is a featureEdgeMesh, read from constant/triSurface for now.
+    features
+    (
+        //{
+        //    file "someLine.eMesh";
+        //    level 2;
+        //}
+    );
+
+
+
+    // Surface based refinement
+    // ~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Specifies two levels for every surface. The first is the minimum level,
+    // every cell intersecting a surface gets refined up to the minimum level.
+    // The second level is the maximum level. Cells that 'see' multiple
+    // intersections where the intersections make an
+    // angle > resolveFeatureAngle get refined up to the maximum level.
+
     refinementSurfaces
     {
         motorBike
         {
+            // Surface-wise min and max refinement level
             level           ( 5 6 );
         }
     }
 
+    // Resolve sharp angles
     resolveFeatureAngle 30;
+
+
+    // Region-wise refinement
+    // ~~~~~~~~~~~~~~~~~~~~~~
+
+    // Specifies refinement level for cells in relation to a surface. One of
+    // three modes
+    // - distance. 'levels' specifies per distance to the surface the
+    //   wanted refinement level. The distances need to be specified in
+    //   descending order.
+    // - inside. 'levels' is only one entry and only the level is used. All
+    //   cells inside the surface get refined up to the level. The surface
+    //   needs to be closed for this to be possible.
+    // - outside. Same but cells outside.
     refinementRegions
     {
         refinementBox
@@ -62,19 +136,45 @@ castellatedMeshControls
         }
     }
 
+
+
+    // Mesh selection
+    // ~~~~~~~~~~~~~~
+
+    // After refinement patches get added for all refinementSurfaces and
+    // all cells intersecting the surfaces get put into these patches. The
+    // section reachable from the locationInMesh is kept.
+    // NOTE: This point should never be on a face, always inside a cell, even
+    // after refinement.
     locationInMesh  ( 3 3 0.43 );
 }
 
+
+// Settings for the snapping.
 snapControls
 {
+    //- Number of patch smoothing iterations before finding correspondence
+    //  to surface
     nSmoothPatch    3;
+
+    //- Relative distance for points to be attracted by surface feature point
+    //  or edge. True distance is this factor times local
+    //  maximum edge length.
     tolerance       4;
+
+    //- Number of mesh displacement relaxation iterations.
     nSolveIter      30;
+
+    //- Maximum number of snapping relaxation iterations. Should stop
+    //  before upon reaching a correct mesh.
     nRelaxIter      5;
 }
 
+
+// Settings for the layer addition.
 addLayersControls
 {
+    // Per final patch (so not geometry!) the layer information
     layers
     {
         minZ
@@ -418,41 +518,126 @@ addLayersControls
         }
     }
 
+    // Expansion factor for layer mesh
     expansionRatio  1;
+
+    //- Wanted thickness of final added cell layer. If multiple layers
+    //  is the
+    //  thickness of the layer furthest away from the wall.
+    //  Relative to undistorted size of cell outside layer.
     finalLayerRatio 0.3;
+
+    //- Minimum thickness of cell layer. If for any reason layer
+    //  cannot be above minThickness do not add layer.
+    //  Relative to undistorted size of cell outside layer.
     minThickness    0.1;
+
+    //- If points get not extruded do nGrow layers of connected faces that are
+    //  also not grown. This helps convergence of the layer addition process
+    //  close to features.
     nGrow           1;
+
+
+    // Advanced settings
+
+    //- When not to extrude surface. 0 is flat surface, 90 is when two faces
+    //  make straight angle.
     featureAngle    30;
+
+    //- Maximum number of snapping relaxation iterations. Should stop
+    //  before upon reaching a correct mesh.
     nRelaxIter      3;
+
+    // Number of smoothing iterations of surface normals
     nSmoothSurfaceNormals 1;
+
+    // Number of smoothing iterations of interior mesh movement direction
     nSmoothNormals  3;
+
+    // Smooth layer thickness over surface patches
     nSmoothThickness 10;
+
+    // Stop layer growth on highly warped cells
     maxFaceThicknessRatio 0.5;
+
+    // Reduce layer growth where ratio thickness to medial
+    // distance is large
     maxThicknessToMedialRatio 0.3;
+
+    // Angle used to pick up medial axis points
     minMedianAxisAngle 130;
+
+    // Create buffer region for new layer terminations
     nBufferCellsNoExtrude 0;
 }
 
+
+// Generic mesh quality settings. At any undoable phase these determine
+// where to undo.
 meshQualityControls
 {
+    //- Maximum non-orthogonality allowed. Set to 180 to disable.
     maxNonOrtho     65;
+
+    //- Max skewness allowed. Set to <0 to disable.
     maxBoundarySkewness 20;
     maxInternalSkewness 4;
+
+    //- Max concaveness allowed. Is angle (in degrees) below which concavity
+    //  is allowed. 0 is straight face, <0 would be convex face.
+    //  Set to 180 to disable.
     maxConcave      80;
+
+    //- Minimum projected area v.s. actual area. Set to -1 to disable.
     minFlatness     0.5;
+
+    //- Minimum pyramid volume. Is absolute volume of cell pyramid.
+    //  Set to very negative number (e.g. -1E30) to disable.
     minVol          1e-13;
+
+    //- Minimum face area. Set to <0 to disable.
     minArea         -1;
+
+    //- Minimum face twist. Set to <-1 to disable. dot product of face normal
+    //- and face centre triangles normal
     minTwist        0.02;
+
+    //- minimum normalised cell determinant
+    //- 1 = hex, <= 0 = folded or flattened illegal cell
     minDeterminant  0.001;
+
+    //- minFaceWeight (0 -> 0.5)
     minFaceWeight   0.02;
+
+    //- minVolRatio (0 -> 1)
     minVolRatio     0.01;
+
+    //must be >0 for Fluent compatibility
     minTriangleTwist -1;
+
+
+    // Advanced
+
+    //- Number of error distribution iterations
     nSmoothScale    4;
+    //- amount to scale back displacement at error points
     errorReduction  0.75;
 }
 
+
+
+// Advanced
+
+// Flags for optional output
+// 0 : only write final meshes
+// 1 : write intermediate meshes
+// 2 : write volScalarField with cellLevel for postprocessing
+// 4 : write current intersections as .obj files
 debug           0;
 
+
+// Merge tolerance. Is fraction of overall bounding box of initial mesh.
+// Note: the write tolerance needs to be higher than this.
 mergeTolerance  1e-06;
 
 

From ef8ac6ff9e12eb2446aae823936bb1ae9073abf4 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs@hunt.opencfd.co.uk>
Date: Thu, 19 Feb 2009 22:41:23 +0000
Subject: [PATCH 02/11] message change

---
 src/OpenFOAM/db/Time/findInstance.C | 21 +++++++++++----------
 1 file changed, 11 insertions(+), 10 deletions(-)

diff --git a/src/OpenFOAM/db/Time/findInstance.C b/src/OpenFOAM/db/Time/findInstance.C
index f21a524c80..a19d3bda63 100644
--- a/src/OpenFOAM/db/Time/findInstance.C
+++ b/src/OpenFOAM/db/Time/findInstance.C
@@ -58,8 +58,9 @@ Foam::word Foam::Time::findInstance
     {
         if (debug)
         {
-            Info<< "Time::findInstance(const fileName&, const word&) : "
-                << "found \"" << name
+            Info<< "Time::findInstance"
+                "(const fileName&, const word&, const IOobject::readOption)"
+                << " : found \"" << name
                 << "\" in " << timeName()/dir
                 << endl;
         }
@@ -98,8 +99,8 @@ Foam::word Foam::Time::findInstance
             if (debug)
             {
                 Info<< "Time::findInstance"
-                    "(const fileName&,const word&) : "
-                    << "found \"" << name
+                    "(const fileName&, const word&, const IOobject::readOption)"
+                    << " : found \"" << name
                     << "\" in " << ts[instanceI].name()/dir
                     << endl;
             }
@@ -129,8 +130,8 @@ Foam::word Foam::Time::findInstance
         if (debug)
         {
             Info<< "Time::findInstance"
-                   "(const fileName&,const word&) : "
-                << "found \"" << name
+                "(const fileName&, const word&, const IOobject::readOption)"
+                << " : found \"" << name
                 << "\" in " << constant()/dir
                 << endl;
         }
@@ -141,10 +142,10 @@ Foam::word Foam::Time::findInstance
     if (rOpt == IOobject::MUST_READ)
     {
         FatalErrorIn
-            (
-                "Time::findInstance(const fileName&,const word&)"
-            )
-            << "Cannot find file \"" << name << "\" in directory "
+        (
+            "Time::findInstance"
+            "(const fileName&, const word&, const IOobject::readOption)"
+        )   << "Cannot find file \"" << name << "\" in directory "
             << constant()/dir
             << exit(FatalError);
     }

From a1525f016abe367e8bde56f4d756031775be056e Mon Sep 17 00:00:00 2001
From: mattijs <mattijs@hunt.opencfd.co.uk>
Date: Fri, 20 Feb 2009 16:47:09 +0000
Subject: [PATCH 03/11] extrapolate to any non-constraint patch

---
 .../PV3FoamReader/PV3FoamReader_SM.xml        |  8 ++--
 .../PV3FoamReader/vtkPV3FoamReader.cxx        |  2 +-
 .../PV3FoamReader/vtkPV3FoamReader.h          |  8 ++--
 .../PV3FoamReader/vtkPV3Foam/vtkPV3Foam.C     | 37 ++++++++++++++++---
 .../vtkPV3Foam/vtkPV3FoamVolFields.H          |  4 +-
 5 files changed, 43 insertions(+), 16 deletions(-)

diff --git a/applications/utilities/postProcessing/graphics/PV3FoamReader/PV3FoamReader/PV3FoamReader_SM.xml b/applications/utilities/postProcessing/graphics/PV3FoamReader/PV3FoamReader/PV3FoamReader_SM.xml
index 644b794a44..ef25ac0378 100644
--- a/applications/utilities/postProcessing/graphics/PV3FoamReader/PV3FoamReader/PV3FoamReader_SM.xml
+++ b/applications/utilities/postProcessing/graphics/PV3FoamReader/PV3FoamReader/PV3FoamReader_SM.xml
@@ -30,16 +30,16 @@
 
     <!-- Global settings -->
 
-    <!-- Extrapolate Walls check-box -->
+    <!-- Extrapolate Patches check-box -->
     <IntVectorProperty
-      name="ExtrapolateWalls"
-      command="SetExtrapolateWalls"
+      name="ExtrapolatePatches"
+      command="SetExtrapolatePatches"
       number_of_elements="1"
       default_values="0"
       animateable="0">
       <BooleanDomain name="bool"/>
       <Documentation>
-        Extrapolate internalField to wall and empty patches
+        Extrapolate internalField to non-constraint patches
       </Documentation>
     </IntVectorProperty>
 
diff --git a/applications/utilities/postProcessing/graphics/PV3FoamReader/PV3FoamReader/vtkPV3FoamReader.cxx b/applications/utilities/postProcessing/graphics/PV3FoamReader/PV3FoamReader/vtkPV3FoamReader.cxx
index 8fc518d745..c16a43deae 100644
--- a/applications/utilities/postProcessing/graphics/PV3FoamReader/PV3FoamReader/vtkPV3FoamReader.cxx
+++ b/applications/utilities/postProcessing/graphics/PV3FoamReader/PV3FoamReader/vtkPV3FoamReader.cxx
@@ -64,7 +64,7 @@ vtkPV3FoamReader::vtkPV3FoamReader()
 
     CacheMesh = 1;
 
-    ExtrapolateWalls = 0;
+    ExtrapolatePatches = 0;
     IncludeSets = 0;
     IncludeZones = 0;
     ShowPatchNames = 0;
diff --git a/applications/utilities/postProcessing/graphics/PV3FoamReader/PV3FoamReader/vtkPV3FoamReader.h b/applications/utilities/postProcessing/graphics/PV3FoamReader/PV3FoamReader/vtkPV3FoamReader.h
index fde87527de..bc21cd8ce8 100644
--- a/applications/utilities/postProcessing/graphics/PV3FoamReader/PV3FoamReader/vtkPV3FoamReader.h
+++ b/applications/utilities/postProcessing/graphics/PV3FoamReader/PV3FoamReader/vtkPV3FoamReader.h
@@ -65,9 +65,9 @@ public:
     vtkGetMacro(CacheMesh, int);
 
     // Description:
-    // FOAM extrapolate internal values onto the walls
-    vtkSetMacro(ExtrapolateWalls, int);
-    vtkGetMacro(ExtrapolateWalls, int);
+    // FOAM extrapolate internal values onto the patches
+    vtkSetMacro(ExtrapolatePatches, int);
+    vtkGetMacro(ExtrapolatePatches, int);
 
     // FOAM read sets control
     vtkSetMacro(IncludeSets, int);
@@ -183,7 +183,7 @@ private:
     int TimeStepRange[2];
     int CacheMesh;
 
-    int ExtrapolateWalls;
+    int ExtrapolatePatches;
     int IncludeSets;
     int IncludeZones;
     int ShowPatchNames;
diff --git a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.C b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.C
index 7a8fd89e9e..0ed3d53185 100644
--- a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.C
+++ b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3Foam.C
@@ -659,29 +659,55 @@ void Foam::vtkPV3Foam::addPatchNames(vtkRenderer* renderer)
         }
     }
 
+    // Count number of zones we're actually going to display. This is truncated
+    // to a max per patch
+
+    const label MAXPATCHZONES = 20;
+
+    label displayZoneI = 0;
+
+    forAll(pbMesh, patchI)
+    {
+        displayZoneI += min(MAXPATCHZONES, nZones[patchI]);
+    }
+
+
     zoneCentre.shrink();
 
     if (debug)
     {
         Info<< "patch zone centres = " << zoneCentre << nl
+            << "displayed zone centres = " << displayZoneI << nl
             << "zones per patch = " << nZones << endl;
     }
 
     // Set the size of the patch labels to max number of zones
-    patchTextActorsPtrs_.setSize(zoneCentre.size());
+    patchTextActorsPtrs_.setSize(displayZoneI);
 
     if (debug)
     {
         Info<< "constructing patch labels" << endl;
     }
 
+    // Actor index
+    displayZoneI = 0;
+
+    // Index in zone centres
     label globalZoneI = 0;
+
     forAll(pbMesh, patchI)
     {
         const polyPatch& pp = pbMesh[patchI];
 
         // Only selected patches will have a non-zero number of zones
-        for (label i=0; i<nZones[patchI]; i++)
+        label nDisplayZones = min(MAXPATCHZONES, nZones[patchI]);
+        label increment = 1;
+        if (nZones[patchI] >= MAXPATCHZONES)
+        {
+            increment = nZones[patchI]/MAXPATCHZONES;
+        }
+
+        for (label i = 0; i < nDisplayZones; i++)
         {
             if (debug)
             {
@@ -719,14 +745,15 @@ void Foam::vtkPV3Foam::addPatchNames(vtkRenderer* renderer)
 
             // Maintain a list of text labels added so that they can be
             // removed later
-            patchTextActorsPtrs_[globalZoneI] = txt;
+            patchTextActorsPtrs_[displayZoneI] = txt;
 
-            globalZoneI++;
+            globalZoneI += increment;
+            displayZoneI++;
         }
     }
 
     // Resize the patch names list to the actual number of patch names added
-    patchTextActorsPtrs_.setSize(globalZoneI);
+    patchTextActorsPtrs_.setSize(displayZoneI);
 
     if (debug)
     {
diff --git a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamVolFields.H b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamVolFields.H
index f217d1a889..6c8d32f470 100644
--- a/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamVolFields.H
+++ b/applications/utilities/postProcessing/graphics/PV3FoamReader/vtkPV3Foam/vtkPV3FoamVolFields.H
@@ -132,8 +132,8 @@ void Foam::vtkPV3Foam::convertVolFields
                 isType<emptyFvPatchField<Type> >(ptf)
              ||
                 (
-                    typeid(patches[patchId]) == typeid(wallPolyPatch)
-                 && reader_->GetExtrapolateWalls()
+                    reader_->GetExtrapolatePatches()
+                && !polyPatch::constraintType(patches[patchId].type())
                 )
             )
             {

From 4b623e594cb4d5bca8b7b7ac3e4c34948b18bcfe Mon Sep 17 00:00:00 2001
From: mattijs <mattijs@hunt.opencfd.co.uk>
Date: Fri, 20 Feb 2009 16:47:33 +0000
Subject: [PATCH 04/11] wildcards in dictionary

---
 .../layerParameters/layerParameters.C         |  41 ++++-
 .../refinementSurfaces/refinementSurfaces.C   | 147 ++++++++++--------
 2 files changed, 117 insertions(+), 71 deletions(-)

diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C
index ecfff62ac3..ddf162dd0f 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C
@@ -29,6 +29,7 @@ License
 #include "mathematicalConstants.H"
 #include "refinementSurfaces.H"
 #include "searchableSurfaces.H"
+#include "regExp.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -300,10 +301,44 @@ Foam::layerParameters::layerParameters
             //    readScalar(layerDict.lookup("minThickness"));
         }
     }
+
+
+    // Check whether layer specification matches any patches
+    const List<keyType> wildCards = layersDict.keys(true);
+
+    forAll(wildCards, i)
+    {
+        regExp re(wildCards[i]);
+
+        bool hasMatch = false;
+        forAll(boundaryMesh, patchI)
+        {
+            if (re.match(boundaryMesh[patchI].name()))
+            {
+                hasMatch = true;
+                break;
+            }
+        }
+        if (!hasMatch)
+        {
+            IOWarningIn("layerParameters::layerParameters(..)", layersDict)
+                << "Wildcard layer specification for " << wildCards[i]
+                << " does not match any patch." << endl;
+        }
+    }
+
+    const List<keyType> nonWildCards = layersDict.keys(false);
+
+    forAll(nonWildCards, i)
+    {
+        if (boundaryMesh.findPatchID(nonWildCards[i]) == -1)
+        {
+            IOWarningIn("layerParameters::layerParameters(..)", layersDict)
+                << "Layer specification for " << nonWildCards[i]
+                << " does not match any patch." << endl;
+        }
+    }
 }
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-
 // ************************************************************************* //
diff --git a/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C b/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
index a626103cd1..6bbb1c8d39 100644
--- a/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
+++ b/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
@@ -265,101 +265,112 @@ Foam::refinementSurfaces::refinementSurfaces
     zoneInside_(surfacesDict.size()),
     regionOffset_(surfacesDict.size())
 {
-    labelList globalMinLevel(surfacesDict.size(), 0);
-    labelList globalMaxLevel(surfacesDict.size(), 0);
-    scalarField globalAngle(surfacesDict.size(), -GREAT);
-    List<Map<label> > regionMinLevel(surfacesDict.size());
-    List<Map<label> > regionMaxLevel(surfacesDict.size());
-    List<Map<scalar> > regionAngle(surfacesDict.size());
+    // Wilcard specification : loop over all surface, all regions
+    // and try to find a match.
 
+    // Count number of surfaces.
     label surfI = 0;
-    forAllConstIter(dictionary, surfacesDict, iter)
+    forAll(allGeometry.names(), geomI)
     {
-        names_[surfI] = iter().keyword();
+        const word& geomName = allGeometry_.names()[geomI];
 
-        surfaces_[surfI] = allGeometry_.findSurfaceID(names_[surfI]);
-
-        if (surfaces_[surfI] == -1)
+        if (surfacesDict.found(geomName))
         {
-            FatalErrorIn
-            (
-                "refinementSurfaces::refinementSurfaces"
-                "(const searchableSurfaces&, const dictionary>&"
-            )   << "No surface called " << iter().keyword() << endl
-                << "Valid surfaces are " << allGeometry_.names()
-                << exit(FatalError);
+            surfI++;
         }
-        const dictionary& dict = surfacesDict.subDict(iter().keyword());
+    }
 
-        const labelPair refLevel(dict.lookup("level"));
-        globalMinLevel[surfI] = refLevel[0];
-        globalMaxLevel[surfI] = refLevel[1];
+    // Size lists
+    surfaces_.setSize(surfI);
+    names_.setSize(surfI);
+    faceZoneNames_.setSize(surfI);
+    cellZoneNames_.setSize(surfI);
+    zoneInside_.setSize(surfI);
+    regionOffset_.setSize(surfI);
 
-        // Global zone names per surface
-        if (dict.found("faceZone"))
+    labelList globalMinLevel(surfI, 0);
+    labelList globalMaxLevel(surfI, 0);
+    scalarField globalAngle(surfI, -GREAT);
+    List<Map<label> > regionMinLevel(surfI);
+    List<Map<label> > regionMaxLevel(surfI);
+    List<Map<scalar> > regionAngle(surfI);
+
+    surfI = 0;
+    forAll(allGeometry.names(), geomI)
+    {
+        const word& geomName = allGeometry_.names()[geomI];
+
+        if (surfacesDict.found(geomName))
         {
-            dict.lookup("faceZone") >> faceZoneNames_[surfI];
-            dict.lookup("cellZone") >> cellZoneNames_[surfI];
-            dict.lookup("zoneInside") >> zoneInside_[surfI];
-        }
+            const dictionary& dict = surfacesDict.subDict(geomName);
 
-        // Global perpendicular angle
-        if (dict.found("perpendicularAngle"))
-        {
-            globalAngle[surfI] = readScalar(dict.lookup("perpendicularAngle"));
-        }
+            names_[surfI] = geomName;
+            surfaces_[surfI] = geomI;
 
-        if (dict.found("regions"))
-        {
-            const dictionary& regionsDict = dict.subDict("regions");
-            const wordList& regionNames =
-                allGeometry_[surfaces_[surfI]].regions();
+            const labelPair refLevel(dict.lookup("level"));
+            globalMinLevel[surfI] = refLevel[0];
+            globalMaxLevel[surfI] = refLevel[1];
 
-            forAllConstIter(dictionary, regionsDict, iter)
+            // Global zone names per surface
+            if (dict.found("faceZone"))
             {
-                const word& key = iter().keyword();
+                dict.lookup("faceZone") >> faceZoneNames_[surfI];
+                dict.lookup("cellZone") >> cellZoneNames_[surfI];
+                dict.lookup("zoneInside") >> zoneInside_[surfI];
+            }
 
-                if (regionsDict.isDict(key))
+            // Global perpendicular angle
+            if (dict.found("perpendicularAngle"))
+            {
+                globalAngle[surfI] = readScalar
+                (
+                    dict.lookup("perpendicularAngle")
+                );
+            }
+
+            if (dict.found("regions"))
+            {
+                const dictionary& regionsDict = dict.subDict("regions");
+                const wordList& regionNames =
+                    allGeometry_[surfaces_[surfI]].regions();
+
+                forAll(regionNames, regionI)
                 {
-                    // Get the dictionary for region iter.keyword()
-                    const dictionary& regionDict = regionsDict.subDict(key);
-
-                    label regionI = findIndex(regionNames, key);
-                    if (regionI == -1)
+                    if (regionsDict.found(regionNames[regionI]))
                     {
-                        FatalErrorIn
+                        // Get the dictionary for region 
+                        const dictionary& regionDict = regionsDict.subDict
                         (
-                            "refinementSurfaces::refinementSurfaces"
-                            "(const searchableSurfaces&, const dictionary>&"
-                        )   << "No region called " << key << " on surface "
-                            << allGeometry_[surfaces_[surfI]].name() << endl
-                            << "Valid regions are " << regionNames
-                            << exit(FatalError);
-                    }
-
-                    const labelPair refLevel(regionDict.lookup("level"));
-
-                    regionMinLevel[surfI].insert(regionI, refLevel[0]);
-                    regionMaxLevel[surfI].insert(regionI, refLevel[1]);
-
-                    if (regionDict.found("perpendicularAngle"))
-                    {
-                        regionAngle[surfI].insert
-                        (
-                            regionI,
-                            readScalar(regionDict.lookup("perpendicularAngle"))
+                            regionNames[regionI]
                         );
+
+                        const labelPair refLevel(regionDict.lookup("level"));
+
+                        regionMinLevel[surfI].insert(regionI, refLevel[0]);
+                        regionMaxLevel[surfI].insert(regionI, refLevel[1]);
+
+                        if (regionDict.found("perpendicularAngle"))
+                        {
+                            regionAngle[surfI].insert
+                            (
+                                regionI,
+                                readScalar
+                                (
+                                    regionDict.lookup("perpendicularAngle")
+                                )
+                            );
+                        }
                     }
                 }
             }
+            surfI++;
         }
-        surfI++;
     }
 
     // Calculate local to global region offset
     label nRegions = 0;
 
-    forAll(surfacesDict, surfI)
+    forAll(surfaces_, surfI)
     {
         regionOffset_[surfI] = nRegions;
         nRegions += allGeometry_[surfaces_[surfI]].regions().size();

From af8a42067c6552c755b89d27606caf8ca20d39d7 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs@hunt.opencfd.co.uk>
Date: Tue, 24 Feb 2009 12:05:35 +0000
Subject: [PATCH 05/11] collapse cell detection

---
 .../mesh/generation/snappyHexMesh/snappyHexMeshDict          | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
index 43d550f406..58f324160e 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
@@ -316,6 +316,11 @@ meshQualityControls
     //must be >0 for Fluent compatibility
     minTriangleTwist -1;
 
+    //- if >0 : preserve single cells with all points on the surface if the
+    //  resulting volume after snapping is larger than minVolFraction times old
+    //  volume. If <0 : delete always.
+    minVolFraction 0.1;
+
 
     // Advanced
 

From 2d81360722bda80b571907be1b4945e91ed203f2 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs@hunt.opencfd.co.uk>
Date: Tue, 24 Feb 2009 12:06:28 +0000
Subject: [PATCH 06/11] removed instance searching

---
 src/meshTools/searchableSurface/searchableSurfaces.C | 10 ----------
 src/meshTools/searchableSurface/searchableSurfaces.H |  2 ++
 2 files changed, 2 insertions(+), 10 deletions(-)

diff --git a/src/meshTools/searchableSurface/searchableSurfaces.C b/src/meshTools/searchableSurface/searchableSurfaces.C
index 46c968f153..cec34bd7e7 100644
--- a/src/meshTools/searchableSurface/searchableSurfaces.C
+++ b/src/meshTools/searchableSurface/searchableSurfaces.C
@@ -39,8 +39,6 @@ defineTypeNameAndDebug(searchableSurfaces, 0);
 }
 
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct with length.
@@ -187,14 +185,6 @@ Foam::searchableSurfaces::searchableSurfaces
         // their object name. Maybe have stlTriSurfaceMesh which appends .stl
         // when reading/writing?
         namedIO().rename(key);  // names_[surfI]
-        if (namedIO().local() != word::null)
-        {
-            namedIO().instance() = namedIO().time().findInstance
-            (
-                namedIO().local(),
-                namedIO().name()
-            );
-        }
 
         // Create and hook surface
         set
diff --git a/src/meshTools/searchableSurface/searchableSurfaces.H b/src/meshTools/searchableSurface/searchableSurfaces.H
index ea91fcdee9..af02ae1883 100644
--- a/src/meshTools/searchableSurface/searchableSurfaces.H
+++ b/src/meshTools/searchableSurface/searchableSurfaces.H
@@ -69,6 +69,8 @@ class searchableSurfaces
         labelList allSurfaces_;
 
 
+    // Private Member Functions
+
         //- Disallow default bitwise copy construct
         searchableSurfaces(const searchableSurfaces&);
 

From 78b10babac17635e1c3126ec4e5584c847987a7a Mon Sep 17 00:00:00 2001
From: mattijs <mattijs@hunt.opencfd.co.uk>
Date: Tue, 24 Feb 2009 12:07:06 +0000
Subject: [PATCH 07/11] collection

---
 src/meshTools/Make/files                      |   1 +
 .../searchableSurfaceCollection.C             | 524 ++++++++++++++++++
 .../searchableSurfaceCollection.H             | 212 +++++++
 3 files changed, 737 insertions(+)
 create mode 100644 src/meshTools/searchableSurface/searchableSurfaceCollection.C
 create mode 100644 src/meshTools/searchableSurface/searchableSurfaceCollection.H

diff --git a/src/meshTools/Make/files b/src/meshTools/Make/files
index a93cd58400..9f849c5894 100644
--- a/src/meshTools/Make/files
+++ b/src/meshTools/Make/files
@@ -60,6 +60,7 @@ $(searchableSurface)/searchablePlane.C
 $(searchableSurface)/searchablePlate.C
 $(searchableSurface)/searchableSphere.C
 $(searchableSurface)/searchableSurface.C
+$(searchableSurface)/searchableSurfaceCollection.C
 $(searchableSurface)/searchableSurfaces.C
 $(searchableSurface)/searchableSurfacesQueries.C
 $(searchableSurface)/searchableSurfaceWithGaps.C
diff --git a/src/meshTools/searchableSurface/searchableSurfaceCollection.C b/src/meshTools/searchableSurface/searchableSurfaceCollection.C
new file mode 100644
index 0000000000..80666ef1f9
--- /dev/null
+++ b/src/meshTools/searchableSurface/searchableSurfaceCollection.C
@@ -0,0 +1,524 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "searchableSurfaceCollection.H"
+#include "addToRunTimeSelectionTable.H"
+#include "SortableList.H"
+#include "Time.H"
+#include "ListOps.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(searchableSurfaceCollection, 0);
+addToRunTimeSelectionTable(searchableSurface, searchableSurfaceCollection, dict);
+
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::searchableSurfaceCollection::findNearest
+(
+    const pointField& samples,
+    scalarField& minDistSqr,
+    List<pointIndexHit>& nearestInfo,
+    labelList& nearestSurf
+) const
+{
+    // Initialise
+    nearestInfo.setSize(samples.size());
+    nearestInfo = pointIndexHit();
+    nearestSurf.setSize(samples.size());
+    nearestSurf = -1;
+
+    List<pointIndexHit> hitInfo(samples.size());
+
+    forAll(subGeom_, surfI)
+    {
+        // Transform then divide
+        tmp<pointField> localSamples = cmptDivide
+        (
+            transform_[surfI].localPosition
+            (
+                samples
+            ),
+            scale_[surfI]
+        );
+        
+        subGeom_[surfI].findNearest(localSamples, minDistSqr, hitInfo);
+
+        forAll(hitInfo, pointI)
+        {
+            if (hitInfo[pointI].hit())
+            {
+                minDistSqr[pointI] = magSqr
+                (   
+                    hitInfo[pointI].hitPoint()
+                  - localSamples()[pointI]
+                );
+
+                // Rework back into global coordinate sys. Multiply then
+                // transform
+                nearestInfo[pointI] = hitInfo[pointI];
+                nearestInfo[pointI].rawPoint() = 
+                    transform_[surfI].globalPosition
+                    (
+                        cmptMultiply
+                        (
+                            nearestInfo[pointI].rawPoint(),
+                            scale_[surfI]
+                        )
+                    );
+                nearestSurf[pointI] = surfI;
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::searchableSurfaceCollection::searchableSurfaceCollection
+(
+    const IOobject& io,
+    const dictionary& dict
+)
+:
+    searchableSurface(io),
+    instance_(dict.size()),
+    scale_(dict.size()),
+    transform_(dict.size()),
+    subGeom_(dict.size())
+{
+    Info<< "SearchableCollection : " << name() << endl;
+
+    label surfI = 0;
+    forAllConstIter(dictionary, dict, iter)
+    {
+        if (dict.isDict(iter().keyword()))
+        {
+            instance_[surfI] = iter().keyword();
+
+            const dictionary& subDict = dict.subDict(instance_[surfI]);
+
+            scale_[surfI] = subDict.lookup("scale");
+            transform_.set
+            (
+                surfI,
+                coordinateSystem::New
+                (
+                    "",
+                    subDict.subDict("transform")
+                )
+            );
+
+            const word subGeomName(subDict.lookup("surface"));
+            //Pout<< "Trying to find " << subGeomName << endl;
+
+            const searchableSurface& s =
+                io.db().lookupObject<searchableSurface>(subGeomName);
+
+            subGeom_.set(surfI, &const_cast<searchableSurface&>(s));
+
+            Info<< "    instance : " << instance_[surfI] << endl;
+            Info<< "    surface  : " << s.name() << endl;
+            Info<< "    scale    : " << scale_[surfI] << endl;
+            Info<< "    coordsys : " << transform_[surfI] << endl;
+
+            surfI++;
+        }
+    }
+    instance_.setSize(surfI);
+    scale_.setSize(surfI);
+    transform_.setSize(surfI);
+    subGeom_.setSize(surfI);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::searchableSurfaceCollection::~searchableSurfaceCollection()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::wordList& Foam::searchableSurfaceCollection::regions() const
+{
+    if (regions_.size() == 0)
+    {
+        regionOffset_.setSize(subGeom_.size());
+
+        DynamicList<word> allRegions;
+        forAll(subGeom_, surfI)
+        {
+            regionOffset_[surfI] = allRegions.size();
+
+            const wordList& subRegions = subGeom_[surfI].regions();
+
+            forAll(subRegions, i)
+            {
+                //allRegions.append(subRegions[i] + "_" + Foam::name(surfI));
+                allRegions.append(instance_[surfI] + "_" + subRegions[i]);
+            }
+        }
+        regions_.transfer(allRegions.shrink());
+    }
+    return regions_;
+}
+
+
+Foam::label Foam::searchableSurfaceCollection::size() const
+{
+    label n = 0;
+    forAll(subGeom_, surfI)
+    {
+        n += subGeom_[surfI].size();
+    }
+    return n;
+}
+
+
+void Foam::searchableSurfaceCollection::findNearest
+(
+    const pointField& samples,
+    const scalarField& nearestDistSqr,
+    List<pointIndexHit>& nearestInfo
+) const
+{
+    // How to scale distance?
+    scalarField minDistSqr(nearestDistSqr);
+
+    labelList nearestSurf;
+    findNearest
+    (
+        samples,
+        minDistSqr,
+        nearestInfo,
+        nearestSurf
+    );
+}
+
+
+void Foam::searchableSurfaceCollection::findLine
+(
+    const pointField& start,
+    const pointField& end,
+    List<pointIndexHit>& info
+) const
+{
+    info.setSize(start.size());
+    info = pointIndexHit();
+
+    // Current nearest (to start) intersection
+    pointField nearest(end);
+
+    List<pointIndexHit> hitInfo(start.size());
+
+    forAll(subGeom_, surfI)
+    {
+        // Starting point
+        tmp<pointField> e0 = cmptDivide
+        (
+            transform_[surfI].localPosition
+            (
+                start
+            ),
+            scale_[surfI]
+        );
+
+        // Current best end point
+        tmp<pointField> e1 = cmptDivide
+        (
+            transform_[surfI].localPosition
+            (
+                nearest
+            ),
+            scale_[surfI]
+        );
+
+        subGeom_[surfI].findLine(e0, e1, hitInfo);
+
+        forAll(hitInfo, pointI)
+        {
+            if (hitInfo[pointI].hit())
+            {
+                // Transform back to global coordinate sys.
+                nearest[pointI] = transform_[surfI].globalPosition
+                (
+                    cmptMultiply
+                    (
+                        hitInfo[pointI].rawPoint(),
+                        scale_[surfI]
+                    )
+                );
+                info[pointI] = hitInfo[pointI];
+                info[pointI].rawPoint() = nearest[pointI];
+            }
+        }
+    }
+
+
+    // Debug check
+    if (false)
+    {
+        forAll(info, pointI)
+        {
+            if (info[pointI].hit())
+            {
+                vector n(end[pointI] - start[pointI]);
+                scalar magN = mag(n);
+
+                if (magN > SMALL)
+                {
+                    n /= mag(n);
+
+                    scalar s = ((info[pointI].rawPoint()-start[pointI])&n);
+
+                    if (s < 0 || s > 1)
+                    {
+                        FatalErrorIn
+                        (
+                            "searchableSurfaceCollection::findLine(..)"
+                        )   << "point:" << info[pointI]
+                            << " s:" << s
+                            << " outside vector "
+                            << " start:" << start[pointI]
+                            << " end:" << end[pointI]
+                            << abort(FatalError);
+                    }
+                }
+            }
+        }
+    }
+}
+
+
+void Foam::searchableSurfaceCollection::findLineAny
+(
+    const pointField& start,
+    const pointField& end,
+    List<pointIndexHit>& info
+) const
+{
+    // To be done ...
+    findLine(start, end, info);
+}
+
+
+void Foam::searchableSurfaceCollection::findLineAll
+(
+    const pointField& start,
+    const pointField& end,
+    List<List<pointIndexHit> >& info
+) const
+{
+    // To be done. Assume for now only one intersection.
+    List<pointIndexHit> nearestInfo;
+    findLine(start, end, nearestInfo);
+
+    info.setSize(start.size());
+    forAll(info, pointI)
+    {
+        if (nearestInfo[pointI].hit())
+        {
+            info[pointI].setSize(1);
+            info[pointI][0] = nearestInfo[pointI];
+        }
+        else
+        {
+            info[pointI].clear();
+        }
+    }
+}
+
+
+void Foam::searchableSurfaceCollection::getRegion
+(
+    const List<pointIndexHit>& info,
+    labelList& region
+) const
+{
+    if (subGeom_.size() == 0)
+    {}
+    else if (subGeom_.size() == 1)
+    {
+        subGeom_[0].getRegion(info, region);
+    }
+    else
+    {
+        region.setSize(info.size());
+        region = -1;
+
+        // Which region did point come from. Retest for now to see which
+        // surface it originates from - crap solution! Should use global indices
+        // in index inside pointIndexHit to do this better.
+
+        pointField samples(info.size());
+        forAll(info, pointI)
+        {
+            if (info[pointI].hit())
+            {
+                samples[pointI] = info[pointI].hitPoint();
+            }
+            else
+            {
+                samples[pointI] = vector::zero;
+            }
+        }
+        //scalarField minDistSqr(info.size(), SMALL);
+        scalarField minDistSqr(info.size(), GREAT);
+
+        labelList nearestSurf;
+        List<pointIndexHit> nearestInfo;
+        findNearest
+        (
+            samples,
+            minDistSqr,
+            nearestInfo,
+            nearestSurf
+        );
+
+        // Check
+        {
+            forAll(info, pointI)
+            {
+                if (info[pointI].hit() && nearestSurf[pointI] == -1)
+                {
+                    FatalErrorIn
+                    (
+                        "searchableSurfaceCollection::getRegion(..)"
+                    )   << "pointI:" << pointI
+                        << " sample:" << samples[pointI]
+                        << " nearest:" << nearestInfo[pointI]
+                        << " nearestsurf:" << nearestSurf[pointI]
+                        << abort(FatalError);
+                }
+            }
+        }
+
+        forAll(subGeom_, surfI)
+        {
+            // Collect points from my surface
+            labelList indices(findIndices(nearestSurf, surfI));
+
+            labelList surfRegion;
+            subGeom_[surfI].getRegion
+            (
+                IndirectList<pointIndexHit>(info, indices),
+                surfRegion
+            );
+            forAll(indices, i)
+            {
+                region[indices[i]] = regionOffset_[surfI] + surfRegion[i];
+            }
+        }
+    }
+}
+
+
+void Foam::searchableSurfaceCollection::getNormal
+(
+    const List<pointIndexHit>& info,
+    vectorField& normal
+) const
+{
+    if (subGeom_.size() == 0)
+    {}
+    else if (subGeom_.size() == 1)
+    {
+        subGeom_[0].getNormal(info, normal);
+    }
+    else
+    {
+        normal.setSize(info.size());
+
+        // See above - crap retest to find surface point originates from.
+        pointField samples(info.size());
+        forAll(info, pointI)
+        {
+            if (info[pointI].hit())
+            {
+                samples[pointI] = info[pointI].hitPoint();
+            }
+            else
+            {
+                samples[pointI] = vector::zero;
+            }
+        }
+        //scalarField minDistSqr(info.size(), SMALL);
+        scalarField minDistSqr(info.size(), GREAT);
+
+        labelList nearestSurf;
+        List<pointIndexHit> nearestInfo;
+        findNearest
+        (
+            samples,
+            minDistSqr,
+            nearestInfo,
+            nearestSurf
+        );
+
+
+        forAll(subGeom_, surfI)
+        {
+            // Collect points from my surface
+            labelList indices(findIndices(nearestSurf, surfI));
+
+            vectorField surfNormal;
+            subGeom_[surfI].getNormal
+            (
+                IndirectList<pointIndexHit>(info, indices),
+                surfNormal
+            );
+            forAll(indices, i)
+            {
+                normal[indices[i]] = surfNormal[i];
+            }
+        }
+    }
+}
+
+
+void Foam::searchableSurfaceCollection::getVolumeType
+(
+    const pointField& points,
+    List<volumeType>& volType
+) const
+{
+    FatalErrorIn
+    (
+        "searchableSurfaceCollection::getVolumeType(const pointField&"
+        ", List<volumeType>&) const"
+    )   << "Volume type not supported for collection."
+        << exit(FatalError);
+}
+
+
+// ************************************************************************* //
diff --git a/src/meshTools/searchableSurface/searchableSurfaceCollection.H b/src/meshTools/searchableSurface/searchableSurfaceCollection.H
new file mode 100644
index 0000000000..c85eaa1026
--- /dev/null
+++ b/src/meshTools/searchableSurface/searchableSurfaceCollection.H
@@ -0,0 +1,212 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::searchableSurfaceCollection
+
+Description
+    Union of transformed searchableSurfaces
+
+SourceFiles
+    searchableSurfaceCollection.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef searchableSurfaceCollection_H
+#define searchableSurfaceCollection_H
+
+#include "searchableSurface.H"
+#include "treeBoundBox.H"
+#include "coordinateSystem.H"
+#include "UPtrList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+
+/*---------------------------------------------------------------------------*\
+                           Class searchableSurfaceCollection Declaration
+\*---------------------------------------------------------------------------*/
+
+class searchableSurfaceCollection
+:
+    public searchableSurface
+{
+private:
+
+    // Private Member Data
+
+        // Per instance data
+
+            //- instance name
+            wordList instance_;
+
+            //- scaling vector
+            vectorField scale_;
+
+            //- transformation
+            PtrList<coordinateSystem> transform_;
+
+            UPtrList<searchableSurface> subGeom_;
+
+        //- Region names
+        mutable wordList regions_;
+        //- From individual regions to collection regions
+        mutable labelList regionOffset_;
+
+
+    // Private Member Functions
+
+        //- Find point nearest to sample. Updates minDistSqr. Sets nearestInfo
+        //  and surface index
+        void findNearest
+        (
+            const pointField& samples,
+            scalarField& minDistSqr,
+            List<pointIndexHit>& nearestInfo,
+            labelList& nearestSurf
+        ) const;
+
+
+        //- Disallow default bitwise copy construct
+        searchableSurfaceCollection(const searchableSurfaceCollection&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const searchableSurfaceCollection&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("searchableSurfaceCollection");
+
+
+    // Constructors
+
+        //- Construct from dictionary (used by searchableSurface)
+        searchableSurfaceCollection
+        (
+            const IOobject& io,
+            const dictionary& dict
+        );
+
+    // Destructor
+
+        virtual ~searchableSurfaceCollection();
+
+
+    // Member Functions
+
+        virtual const wordList& regions() const;
+
+        //- Whether supports volume type below
+        virtual bool hasVolumeType() const
+        {
+            return false;
+        }
+
+        //- Range of local indices that can be returned.
+        virtual label size() const;
+
+
+        // Multiple point queries.
+
+            virtual void findNearest
+            (
+                const pointField& sample,
+                const scalarField& nearestDistSqr,
+                List<pointIndexHit>&
+            ) const;
+
+            virtual void findLine
+            (
+                const pointField& start,
+                const pointField& end,
+                List<pointIndexHit>&
+            ) const;
+
+            virtual void findLineAny
+            (
+                const pointField& start,
+                const pointField& end,
+                List<pointIndexHit>&
+            ) const;
+
+            //- Get all intersections in order from start to end.
+            virtual void findLineAll
+            (
+                const pointField& start,
+                const pointField& end,
+                List<List<pointIndexHit> >&
+            ) const;
+
+            //- From a set of points and indices get the region
+            virtual void getRegion
+            (
+                const List<pointIndexHit>&,
+                labelList& region
+            ) const;
+
+            //- From a set of points and indices get the normal
+            virtual void getNormal
+            (
+                const List<pointIndexHit>&,
+                vectorField& normal
+            ) const;
+
+            //- Determine type (inside/outside/mixed) for point. unknown if
+            //  cannot be determined (e.g. non-manifold surface)
+            virtual void getVolumeType
+            (
+                const pointField&,
+                List<volumeType>&
+            ) const;
+
+
+        // regIOobject implementation
+
+            bool writeData(Ostream&) const
+            {
+                notImplemented
+                (
+                    "searchableSurfaceCollection::writeData(Ostream&) const"
+                );
+                return false;
+            }
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

From 364e8e001ddfb1e84a57dac327ff01f7beab68c8 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs@hunt.opencfd.co.uk>
Date: Tue, 24 Feb 2009 12:07:50 +0000
Subject: [PATCH 08/11] reintroduce comments

---
 .../channel395/constant/postChannelDict       |  4 +++
 .../channelFoam/channel395/system/controlDict | 30 ++++++++++++++++++-
 2 files changed, 33 insertions(+), 1 deletion(-)

diff --git a/tutorials/incompressible/channelFoam/channel395/constant/postChannelDict b/tutorials/incompressible/channelFoam/channel395/constant/postChannelDict
index 89fcd50809..96bbe30d76 100644
--- a/tutorials/incompressible/channelFoam/channel395/constant/postChannelDict
+++ b/tutorials/incompressible/channelFoam/channel395/constant/postChannelDict
@@ -15,10 +15,14 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+// Seed patches to start layering from
 patches         ( bottomWall );
 
+// Direction in which the layers are
 component       y;
 
+// Is the mesh symmetric? If so average(symmetric fields) or
+// subtract(asymmetric) contributions from both halves
 symmetric       true;
 
 
diff --git a/tutorials/incompressible/channelFoam/channel395/system/controlDict b/tutorials/incompressible/channelFoam/channel395/system/controlDict
index 2ca28c29a1..218a76d861 100644
--- a/tutorials/incompressible/channelFoam/channel395/system/controlDict
+++ b/tutorials/incompressible/channelFoam/channel395/system/controlDict
@@ -43,7 +43,35 @@ timePrecision   6;
 
 runTimeModifiable yes;
 
-functions       ( fieldAverage1 { type fieldAverage ; functionObjectLibs ( "libfieldFunctionObjects.so" ) ; enabled true ; outputControl outputTime ; fields ( U { mean on ; prime2Mean on ; base time ; } p { mean on ; prime2Mean on ; base time ; } ) ; } );
+functions
+(
+    fieldAverage1
+    {
+        type fieldAverage;
+
+        functionObjectLibs ( "libfieldFunctionObjects.so" );
+
+        enabled true;
+
+        outputControl outputTime;
+
+        fields
+        (
+            U
+            {
+                mean on;
+                prime2Mean on;
+                base time;
+            }
+            p
+            {
+                mean on;
+                prime2Mean on;
+                base time;
+            }
+        );
+    }
+);
 
 
 // ************************************************************************* //

From ab49c60394698cf42cb579f728b96e328a7c83e7 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs@hunt.opencfd.co.uk>
Date: Tue, 24 Feb 2009 12:08:54 +0000
Subject: [PATCH 09/11] instance searching

---
 .../distributedTriSurfaceMesh.C               |  27 ----
 .../searchableSurface/triSurfaceMesh.C        | 131 +++++++++++++++++-
 .../searchableSurface/triSurfaceMesh.H        |  13 +-
 3 files changed, 137 insertions(+), 34 deletions(-)

diff --git a/src/meshTools/searchableSurface/distributedTriSurfaceMesh.C b/src/meshTools/searchableSurface/distributedTriSurfaceMesh.C
index 0602046703..870f7e8324 100644
--- a/src/meshTools/searchableSurface/distributedTriSurfaceMesh.C
+++ b/src/meshTools/searchableSurface/distributedTriSurfaceMesh.C
@@ -1347,19 +1347,6 @@ Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh
 Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh(const IOobject& io)
 :
     triSurfaceMesh(io),
-//    triSurfaceMesh
-//    (
-//        IOobject
-//        (
-//            io.name(),
-//            io.db().time().findInstanceDir(io.local()),
-//            io.local(),
-//            io.db(),
-//            io.readOpt(),
-//            io.writeOpt(),
-//            io.registerObject()
-//        )
-//    ),
     dict_
     (
         IOobject
@@ -1385,20 +1372,6 @@ Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh
 )
 :
     triSurfaceMesh(io, dict),
-//    triSurfaceMesh
-//    (
-//        IOobject
-//        (
-//            io.name(),
-//            io.db().time().findInstanceDir(io.local()),
-//            io.local(),
-//            io.db(),
-//            io.readOpt(),
-//            io.writeOpt(),
-//            io.registerObject()
-//        ),
-//        dict
-//    ),
     dict_
     (
         IOobject
diff --git a/src/meshTools/searchableSurface/triSurfaceMesh.C b/src/meshTools/searchableSurface/triSurfaceMesh.C
index 7b9abc1aec..121762d7c4 100644
--- a/src/meshTools/searchableSurface/triSurfaceMesh.C
+++ b/src/meshTools/searchableSurface/triSurfaceMesh.C
@@ -43,6 +43,64 @@ addToRunTimeSelectionTable(searchableSurface, triSurfaceMesh, dict);
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+//// Special version of Time::findInstance that does not check headerOk
+//// to search for instances of raw files
+//Foam::word Foam::triSurfaceMesh::findRawInstance
+//(
+//    const Time& runTime,
+//    const fileName& dir,
+//    const word& name
+//)
+//{
+//    // Check current time first
+//    if (isFile(runTime.path()/runTime.timeName()/dir/name))
+//    {
+//        return runTime.timeName();
+//    }
+//    instantList ts = runTime.times();
+//    label instanceI;
+//
+//    for (instanceI = ts.size()-1; instanceI >= 0; --instanceI)
+//    {
+//        if (ts[instanceI].value() <= runTime.timeOutputValue())
+//        {
+//            break;
+//        }
+//    }
+//
+//    // continue searching from here
+//    for (; instanceI >= 0; --instanceI)
+//    {
+//        if (isFile(runTime.path()/ts[instanceI].name()/dir/name))
+//        {
+//            return ts[instanceI].name();
+//        }
+//    }
+//
+//
+//    // not in any of the time directories, try constant
+//
+//    // Note. This needs to be a hard-coded constant, rather than the
+//    // constant function of the time, because the latter points to
+//    // the case constant directory in parallel cases
+//
+//    if (isFile(runTime.path()/runTime.constant()/dir/name))
+//    {
+//        return runTime.constant();
+//    }
+//
+//    FatalErrorIn
+//    (
+//        "searchableSurfaces::findRawInstance"
+//        "(const Time&, const fileName&, const word&)"
+//    )   << "Cannot find file \"" << name << "\" in directory "
+//        << runTime.constant()/dir
+//        << exit(FatalError);
+//
+//    return runTime.constant();
+//}
+
+
 //- Check file existence
 const Foam::fileName& Foam::triSurfaceMesh::checkFile
 (
@@ -149,7 +207,19 @@ bool Foam::triSurfaceMesh::isSurfaceClosed() const
 Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io, const triSurface& s)
 :
     searchableSurface(io),
-    objectRegistry(io),
+    objectRegistry
+    (
+        IOobject
+        (
+            io.name(),
+            io.instance(),
+            io.local(),
+            io.db(),
+            io.readOpt(),
+            io.writeOpt(),
+            false       // searchableSurface already registered under name
+        )
+    ),
     triSurface(s),
     surfaceClosed_(-1)
 {}
@@ -157,8 +227,34 @@ Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io, const triSurface& s)
 
 Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io)
 :
-    searchableSurface(io),
-    objectRegistry(io),
+    // Find instance for triSurfaceMesh
+    searchableSurface
+    (
+        IOobject
+        (
+            io.name(),
+            io.time().findInstance(io.local(), word::null),
+            io.local(),
+            io.db(),
+            io.readOpt(),
+            io.writeOpt(),
+            io.registerObject()
+        )
+    ),
+    // Reused found instance in objectRegistry
+    objectRegistry
+    (
+        IOobject
+        (
+            io.name(),
+            static_cast<const searchableSurface&>(*this).instance(),
+            io.local(),
+            io.db(),
+            io.readOpt(),
+            io.writeOpt(),
+            false       // searchableSurface already registered under name
+        )
+    ),
     triSurface
     (
         checkFile
@@ -177,8 +273,33 @@ Foam::triSurfaceMesh::triSurfaceMesh
     const dictionary& dict
 )
 :
-    searchableSurface(io),
-    objectRegistry(io),
+    searchableSurface
+    (
+        IOobject
+        (
+            io.name(),
+            io.time().findInstance(io.local(), word::null),
+            io.local(),
+            io.db(),
+            io.readOpt(),
+            io.writeOpt(),
+            io.registerObject()
+        )
+    ),
+    // Reused found instance in objectRegistry
+    objectRegistry
+    (
+        IOobject
+        (
+            io.name(),
+            static_cast<const searchableSurface&>(*this).instance(),
+            io.local(),
+            io.db(),
+            io.readOpt(),
+            io.writeOpt(),
+            false       // searchableSurface already registered under name
+        )
+    ),
     triSurface
     (
         checkFile
diff --git a/src/meshTools/searchableSurface/triSurfaceMesh.H b/src/meshTools/searchableSurface/triSurfaceMesh.H
index aa0af5979c..d7aa33cfca 100644
--- a/src/meshTools/searchableSurface/triSurfaceMesh.H
+++ b/src/meshTools/searchableSurface/triSurfaceMesh.H
@@ -76,6 +76,14 @@ private:
 
     // Private Member Functions
 
+        ////- Helper: find instance of files without header
+        //static word findRawInstance
+        //(
+        //    const Time&,
+        //    const fileName&,
+        //    const word&
+        //);
+
         //- Check file existence
         static const fileName& checkFile
         (
@@ -105,10 +113,11 @@ public:
         //- Construct from triSurface
         triSurfaceMesh(const IOobject&, const triSurface&);
 
-        //- Construct read
+        //- Construct read. Does timeInstance search.
         triSurfaceMesh(const IOobject& io);
 
-        //- Construct from dictionary (used by searchableSurface)
+        //- Construct from IO and dictionary (used by searchableSurface).
+        //  Does timeInstance search.
         //  Dictionary may contain a 'scale' entry (eg, 0.001: mm -> m)
         triSurfaceMesh
         (

From 5a30dd1b018a26cacf0a1af7a26f2adc515be473 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs@hunt.opencfd.co.uk>
Date: Tue, 24 Feb 2009 12:09:17 +0000
Subject: [PATCH 10/11] iso surface

---
 .../sampledSurface/isoSurface/isoSurface.C    | 85 +++++++++++++++----
 1 file changed, 68 insertions(+), 17 deletions(-)

diff --git a/src/sampling/sampledSurface/isoSurface/isoSurface.C b/src/sampling/sampledSurface/isoSurface/isoSurface.C
index c69a51f407..a2c98f4f93 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurface.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurface.C
@@ -30,6 +30,7 @@ License
 #include "syncTools.H"
 #include "addToRunTimeSelectionTable.H"
 #include "slicedVolFields.H"
+#include "surfaceFields.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -343,6 +344,7 @@ void Foam::isoSurface::getNeighbour
 {
     const labelList& own = mesh_.faceOwner();
     const labelList& nei = mesh_.faceNeighbour();
+    const surfaceScalarField& weights = mesh_.weights();
 
     if (mesh_.isInternalFace(faceI))
     {
@@ -363,16 +365,26 @@ void Foam::isoSurface::getNeighbour
             nbrValue = cVals[own[faceI]];
             nbrPoint = mesh_.faceCentres()[faceI];
         }
-        else if
-        (
-            pp.coupled()
-        && !refCast<const coupledPolyPatch>(pp).separated()
-        )
+        else if (pp.coupled())
         {
-            // other side value
-            nbrValue = cVals.boundaryField()[patchI][patchFaceI];
-            // other side cell centre
-            nbrPoint = mesh_.C().boundaryField()[patchI][patchFaceI];
+            if (!refCast<const coupledPolyPatch>(pp).separated())
+            {
+                // Behave as internal face:
+                // other side value
+                nbrValue = cVals.boundaryField()[patchI][patchFaceI];
+                // other side cell centre
+                nbrPoint = mesh_.C().boundaryField()[patchI][patchFaceI];
+            }
+            else
+            {
+                // Do some interpolation for now
+                const scalarField& w = weights.boundaryField()[patchI];
+
+                nbrPoint = mesh_.faceCentres()[faceI];
+                nbrValue =
+                    (1-w[patchFaceI])*cVals[own[faceI]]
+                  + w[patchFaceI]*cVals.boundaryField()[patchI][patchFaceI];
+            }
         }
         else
         {
@@ -846,22 +858,26 @@ Foam::triSurface Foam::isoSurface::stitchTriPoints
         forAll(tris, triI)
         {
             const labelledTri& tri = tris[triI];
-
             const labelList& pFaces = pointFaces[tri[0]];
 
-            // Find the minimum of any duplicates
+            // Find the maximum of any duplicates. Maximum since the tris
+            // below triI
+            // get overwritten so we cannot use them in a comparison.
             label dupTriI = -1;
             forAll(pFaces, i)
             {
-                if (pFaces[i] < triI && tris[pFaces[i]] == tri)
+                label nbrTriI = pFaces[i];
+
+                if (nbrTriI > triI && (tris[nbrTriI] == tri))
                 {
-                    dupTriI = pFaces[i];
+                    dupTriI = nbrTriI;
+                    break;
                 }
             }
 
             if (dupTriI == -1)
             {
-                // There is no lower triangle
+                // There is no (higher numbered) duplicate triangle
                 label newTriI = newToOldTri.size();
                 newToOldTri.append(triI);
                 tris[newTriI] = tris[triI];
@@ -876,6 +892,43 @@ Foam::triSurface Foam::isoSurface::stitchTriPoints
             Pout<< "isoSurface : merged from " << nTris
                 << " down to " << tris.size() << " unique triangles." << endl;
         }
+
+
+        if (debug)
+        {
+            triSurface surf(tris, geometricSurfacePatchList(0), newPoints);
+
+            forAll(surf, faceI)
+            {
+                const labelledTri& f = surf[faceI];
+                const labelList& fFaces = surf.faceFaces()[faceI];
+
+                forAll(fFaces, i)
+                {
+                    label nbrFaceI = fFaces[i];
+
+                    if (nbrFaceI <= faceI)
+                    {
+                        // lower numbered faces already checked
+                        continue;
+                    }
+
+                    const labelledTri& nbrF = surf[nbrFaceI];
+
+                    if (f == nbrF)
+                    {
+                        FatalErrorIn("validTri(const triSurface&, const label)")
+                            << "Check : "
+                            << " triangle " << faceI << " vertices " << f
+                            << " fc:" << f.centre(surf.points())
+                            << " has the same vertices as triangle " << nbrFaceI
+                            << " vertices " << nbrF
+                            << " fc:" << nbrF.centre(surf.points())
+                            << abort(FatalError);
+                    }
+                }
+            }
+        }
     }
 
     return triSurface(tris, geometricSurfacePatchList(0), newPoints, true);
@@ -1526,7 +1579,7 @@ Foam::isoSurface::isoSurface
             faceI++;
         }
 
-        // Mark all points that are not physically coupled (so anything
+        // Mark all boundary points that are not physically coupled (so anything
         // but collocated coupled patches)
         if
         (
@@ -1538,8 +1591,6 @@ Foam::isoSurface::isoSurface
 
             forAll(pp, i)
             {
-                boundaryRegion[faceI-mesh_.nInternalFaces()] = patchI;
-
                 const face& f = mesh_.faces()[faceI];
 
                 forAll(f, fp)

From c8944ce200ed844345c1897e53535cd99710ec89 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs@hunt.opencfd.co.uk>
Date: Tue, 24 Feb 2009 13:45:57 +0000
Subject: [PATCH 11/11] iso surface correction

---
 .../sampledSurface/isoSurface/isoSurface.C    | 277 ++++++++++-------
 .../sampledSurface/isoSurface/isoSurface.H    |  46 +--
 .../isoSurface/isoSurfaceTemplates.C          | 286 ++++++++++++------
 3 files changed, 391 insertions(+), 218 deletions(-)

diff --git a/src/sampling/sampledSurface/isoSurface/isoSurface.C b/src/sampling/sampledSurface/isoSurface/isoSurface.C
index a2c98f4f93..939d77cf1b 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurface.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurface.C
@@ -85,9 +85,74 @@ bool Foam::isoSurface::isEdgeOfFaceCut
 }
 
 
+// Get neighbour value and position.
+void Foam::isoSurface::getNeighbour
+(
+    const labelList& boundaryRegion,
+    const volScalarField& cVals,
+    const label cellI,
+    const label faceI,
+    scalar& nbrValue,
+    point& nbrPoint
+) const
+{
+    const labelList& own = mesh_.faceOwner();
+    const labelList& nei = mesh_.faceNeighbour();
+    const surfaceScalarField& weights = mesh_.weights();
+
+    if (mesh_.isInternalFace(faceI))
+    {
+        label nbr = (own[faceI] == cellI ? nei[faceI] : own[faceI]);
+        nbrValue = cVals[nbr];
+        nbrPoint = mesh_.cellCentres()[nbr];
+    }
+    else
+    {
+        label bFaceI = faceI-mesh_.nInternalFaces();
+        label patchI = boundaryRegion[bFaceI];
+        const polyPatch& pp = mesh_.boundaryMesh()[patchI];
+        label patchFaceI = faceI-pp.start();
+
+        if (isA<emptyPolyPatch>(pp))
+        {
+            // Assume zero gradient
+            nbrValue = cVals[own[faceI]];
+            nbrPoint = mesh_.faceCentres()[faceI];
+        }
+        else if (pp.coupled())
+        {
+            if (!refCast<const coupledPolyPatch>(pp).separated())
+            {
+                // Behave as internal face:
+                // other side value
+                nbrValue = cVals.boundaryField()[patchI][patchFaceI];
+                // other side cell centre
+                nbrPoint = mesh_.C().boundaryField()[patchI][patchFaceI];
+            }
+            else
+            {
+                // Do some interpolation for now
+                const scalarField& w = weights.boundaryField()[patchI];
+
+                nbrPoint = mesh_.faceCentres()[faceI];
+                nbrValue =
+                    (1-w[patchFaceI])*cVals[own[faceI]]
+                  + w[patchFaceI]*cVals.boundaryField()[patchI][patchFaceI];
+            }
+        }
+        else
+        {
+            nbrValue = cVals.boundaryField()[patchI][patchFaceI];
+            nbrPoint = mesh_.faceCentres()[faceI];
+        }
+    }
+}
+
+
 // Determine for every face/cell whether it (possibly) generates triangles.
 void Foam::isoSurface::calcCutTypes
 (
+    const labelList& boundaryRegion,
     const volScalarField& cVals,
     const scalarField& pVals
 )
@@ -103,7 +168,20 @@ void Foam::isoSurface::calcCutTypes
     {
         // CC edge.
         bool ownLower = (cVals[own[faceI]] < iso_);
-        bool neiLower = (cVals[nei[faceI]] < iso_);
+
+        scalar nbrValue;
+        point nbrPoint;
+        getNeighbour
+        (
+            boundaryRegion,
+            cVals,
+            own[faceI],
+            faceI,
+            nbrValue,
+            nbrPoint
+        );
+
+        bool neiLower = (nbrValue < iso_);
 
         if (ownLower != neiLower)
         {
@@ -129,15 +207,29 @@ void Foam::isoSurface::calcCutTypes
 
         if (isA<emptyPolyPatch>(pp))
         {
-            // Assume zero gradient so owner and neighbour/boundary value equal
+            // Still needs special treatment?
 
             forAll(pp, i)
             {
                 bool ownLower = (cVals[own[faceI]] < iso_);
 
+                scalar nbrValue;
+                point nbrPoint;
+                getNeighbour
+                (
+                    boundaryRegion,
+                    cVals,
+                    own[faceI],
+                    faceI,
+                    nbrValue,
+                    nbrPoint
+                );
+
+                bool neiLower = (nbrValue < iso_);
+
                 const face f = mesh_.faces()[faceI];
 
-                if (isEdgeOfFaceCut(pVals, f, ownLower, ownLower))
+                if (isEdgeOfFaceCut(pVals, f, ownLower, neiLower))
                 {
                     faceCutType_[faceI] = CUT;
                 }
@@ -150,7 +242,20 @@ void Foam::isoSurface::calcCutTypes
             forAll(pp, i)
             {
                 bool ownLower = (cVals[own[faceI]] < iso_);
-                bool neiLower = (cVals.boundaryField()[patchI][i] < iso_);
+
+                scalar nbrValue;
+                point nbrPoint;
+                getNeighbour
+                (
+                    boundaryRegion,
+                    cVals,
+                    own[faceI],
+                    faceI,
+                    nbrValue,
+                    nbrPoint
+                );
+
+                bool neiLower = (nbrValue < iso_);
 
                 if (ownLower != neiLower)
                 {
@@ -166,6 +271,7 @@ void Foam::isoSurface::calcCutTypes
                         faceCutType_[faceI] = CUT;
                     }
                 }
+
                 faceI++;
             }
         }
@@ -331,70 +437,6 @@ Foam::pointIndexHit Foam::isoSurface::collapseSurface
 }
 
 
-// Get neighbour value and position.
-void Foam::isoSurface::getNeighbour
-(
-    const labelList& boundaryRegion,
-    const volScalarField& cVals,
-    const label cellI,
-    const label faceI,
-    scalar& nbrValue,
-    point& nbrPoint
-) const
-{
-    const labelList& own = mesh_.faceOwner();
-    const labelList& nei = mesh_.faceNeighbour();
-    const surfaceScalarField& weights = mesh_.weights();
-
-    if (mesh_.isInternalFace(faceI))
-    {
-        label nbr = (own[faceI] == cellI ? nei[faceI] : own[faceI]);
-        nbrValue = cVals[nbr];
-        nbrPoint = mesh_.cellCentres()[nbr];
-    }
-    else
-    {
-        label bFaceI = faceI-mesh_.nInternalFaces();
-        label patchI = boundaryRegion[bFaceI];
-        const polyPatch& pp = mesh_.boundaryMesh()[patchI];
-        label patchFaceI = faceI-pp.start();
-
-        if (isA<emptyPolyPatch>(pp))
-        {
-            // Assume zero gradient
-            nbrValue = cVals[own[faceI]];
-            nbrPoint = mesh_.faceCentres()[faceI];
-        }
-        else if (pp.coupled())
-        {
-            if (!refCast<const coupledPolyPatch>(pp).separated())
-            {
-                // Behave as internal face:
-                // other side value
-                nbrValue = cVals.boundaryField()[patchI][patchFaceI];
-                // other side cell centre
-                nbrPoint = mesh_.C().boundaryField()[patchI][patchFaceI];
-            }
-            else
-            {
-                // Do some interpolation for now
-                const scalarField& w = weights.boundaryField()[patchI];
-
-                nbrPoint = mesh_.faceCentres()[faceI];
-                nbrValue =
-                    (1-w[patchFaceI])*cVals[own[faceI]]
-                  + w[patchFaceI]*cVals.boundaryField()[patchI][patchFaceI];
-            }
-        }
-        else
-        {
-            nbrValue = cVals.boundaryField()[patchI][patchFaceI];
-            nbrPoint = mesh_.faceCentres()[faceI];
-        }
-    }
-}
-
-
 // Determine per cell centre whether all the intersections get collapsed
 // to a single point
 void Foam::isoSurface::calcSnappedCc
@@ -618,13 +660,14 @@ void Foam::isoSurface::calcSnappedPoint
 
         forAll(pFaces, pFaceI)
         {
+            // Create points for all intersections close to point
+            // (i.e. from pyramid edges)
+
             label faceI = pFaces[pFaceI];
             const face& f = mesh_.faces()[faceI];
             label own = mesh_.faceOwner()[faceI];
 
-            // Create points for all intersections close to point
-            // (i.e. from pyramid edges)
-
+            // Get neighbour value
             scalar nbrValue;
             point nbrPoint;
             getNeighbour
@@ -843,6 +886,7 @@ Foam::triSurface Foam::isoSurface::stitchTriPoints
     }
 
 
+
     // Determine 'flat hole' situation (see RMT paper).
     // Two unconnected triangles get connected because (some of) the edges
     // separating them get collapsed. Below only checks for duplicate triangles,
@@ -870,6 +914,9 @@ Foam::triSurface Foam::isoSurface::stitchTriPoints
 
                 if (nbrTriI > triI && (tris[nbrTriI] == tri))
                 {
+                    //Pout<< "Duplicate : " << triI << " verts:" << tri
+                    //    << " to " << nbrTriI << " verts:" << tris[nbrTriI]
+                    //    << endl;
                     dupTriI = nbrTriI;
                     break;
                 }
@@ -1526,11 +1573,17 @@ Foam::isoSurface::isoSurface
 {
     if (debug)
     {
-        Pout<< "isoSurface     :" << nl
-            << "    isoField   : " << cVals.name() << nl
-            << "    isoValue   : " << iso << nl
-            << "    regularise : " << regularise << nl
-            << "    mergeTol   : " << mergeTol << nl
+        Pout<< "isoSurface:" << nl
+            << "    isoField      : " << cVals.name() << nl
+            << "    cell min/max  : "
+            << min(cVals.internalField()) << " / "
+            << max(cVals.internalField()) << nl
+            << "    point min/max : "
+            << min(pVals) << " / "
+            << max(pVals) << nl
+            << "    isoValue      : " << iso << nl
+            << "    regularise    : " << regularise << nl
+            << "    mergeTol      : " << mergeTol << nl
             << endl;
     }
 
@@ -1556,15 +1609,7 @@ Foam::isoSurface::isoSurface
         }
     }
 
-
-    // Determine if any cut through face/cell
-    calcCutTypes(cVals, pVals);
-
-
-    // Determine if point is on boundary. Points on boundaries are never
-    // snapped. Coupled boundaries are handled explicitly so not marked here.
-    PackedBoolList isBoundaryPoint(mesh_.nPoints());
-
+    // Pre-calculate patch-per-face to avoid whichPatch call.
     labelList boundaryRegion(mesh_.nFaces()-mesh_.nInternalFaces());
 
     forAll(patches, patchI)
@@ -1578,31 +1623,13 @@ Foam::isoSurface::isoSurface
             boundaryRegion[faceI-mesh_.nInternalFaces()] = patchI;
             faceI++;
         }
-
-        // Mark all boundary points that are not physically coupled (so anything
-        // but collocated coupled patches)
-        if
-        (
-           !pp.coupled()
-        || refCast<const coupledPolyPatch>(pp).separated()
-        )
-        {
-            label faceI = pp.start();
-
-            forAll(pp, i)
-            {
-                const face& f = mesh_.faces()[faceI];
-
-                forAll(f, fp)
-                {
-                    isBoundaryPoint.set(f[fp], 1);
-                }
-                faceI++;
-            }
-        }
     }
 
 
+    // Determine if any cut through face/cell
+    calcCutTypes(boundaryRegion, cVals, pVals);
+
+
     DynamicList<point> snappedPoints(nCutCells_);
 
     // Per cc -1 or a point inside snappedPoints.
@@ -1635,6 +1662,39 @@ Foam::isoSurface::isoSurface
 
     label nCellSnaps = snappedPoints.size();
 
+
+    // Determine if point is on boundary. Points on boundaries are never
+    // snapped. Coupled boundaries are handled explicitly so not marked here.
+    PackedBoolList isBoundaryPoint(mesh_.nPoints());
+
+    forAll(patches, patchI)
+    {
+        // Mark all boundary points that are not physically coupled (so anything
+        // but collocated coupled patches)
+        const polyPatch& pp = patches[patchI];
+
+        if
+        (
+           !pp.coupled()
+        || refCast<const coupledPolyPatch>(pp).separated()
+        )
+        {
+            label faceI = pp.start();
+
+            forAll(pp, i)
+            {
+                const face& f = mesh_.faces()[faceI];
+
+                forAll(f, fp)
+                {
+                    isBoundaryPoint.set(f[fp], 1);
+                }
+                faceI++;
+            }
+        }
+    }
+
+
     // Per point -1 or a point inside snappedPoints.
     labelList snappedPoint;
     if (regularise)
@@ -1727,7 +1787,8 @@ Foam::isoSurface::isoSurface
     if (debug)
     {
         Pout<< "isoSurface : generated " << triMeshCells.size()
-            << " unmerged triangles." << endl;
+            << " unmerged triangles from " << triPoints.size()
+            << " unmerged points." << endl;
     }
 
 
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurface.H b/src/sampling/sampledSurface/isoSurface/isoSurface.H
index 7dc0588a1c..f6c022c036 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurface.H
+++ b/src/sampling/sampledSurface/isoSurface/isoSurface.H
@@ -135,9 +135,20 @@ class isoSurface
             const bool neiLower
         ) const;
 
+        void getNeighbour
+        (
+            const labelList& boundaryRegion,
+            const volScalarField& cVals,
+            const label cellI,
+            const label faceI,
+            scalar& nbrValue,
+            point& nbrPoint
+        ) const;
+
         //- Set faceCutType,cellCutType.
         void calcCutTypes
         (
+            const labelList& boundaryRegion,
             const volScalarField& cVals,
             const scalarField& pVals
         );
@@ -156,16 +167,6 @@ class isoSurface
             DynamicList<labelledTri, 64>& localTris
         );
 
-        void getNeighbour
-        (
-            const labelList& boundaryRegion,
-            const volScalarField& cVals,
-            const label cellI,
-            const label faceI,
-            scalar& nbrValue,
-            point& nbrPoint
-        ) const;
-
         //- Determine per cc whether all near cuts can be snapped to single
         //  point.
         void calcSnappedCc
@@ -193,37 +194,39 @@ class isoSurface
         template<class Type>
         Type generatePoint
         (
-            const DynamicList<Type>& snappedPoints,
-
             const scalar s0,
             const Type& p0,
-            const label p0Index,
+            const bool hasSnap0,
+            const Type& snapP0,
 
             const scalar s1,
             const Type& p1,
-            const label p1Index
+            const bool hasSnap1,
+            const Type& snapP1
         ) const;
 
         template<class Type>
         void generateTriPoints
         (
-            const DynamicList<Type>& snapped,
-
             const scalar s0,
             const Type& p0,
-            const label p0Index,
+            const bool hasSnap0,
+            const Type& snapP0,
 
             const scalar s1,
             const Type& p1,
-            const label p1Index,
+            const bool hasSnap1,
+            const Type& snapP1,
 
             const scalar s2,
             const Type& p2,
-            const label p2Index,
+            const bool hasSnap2,
+            const Type& snapP2,
 
             const scalar s3,
             const Type& p3,
-            const label p3Index,
+            const bool hasSnap3,
+            const Type& snapP3,
 
             DynamicList<Type>& points
         ) const;
@@ -244,7 +247,8 @@ class isoSurface
 
             const scalar neiVal,
             const Type& neiPt,
-            const label neiSnap,
+            const bool hasNeiSnap,
+            const Type& neiSnapPt,
 
             DynamicList<Type>& triPoints,
             DynamicList<label>& triMeshCells
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C b/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
index e2251aad2a..7e7206edcc 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
@@ -33,15 +33,15 @@ License
 template<class Type>
 Type Foam::isoSurface::generatePoint
 (
-    const DynamicList<Type>& snappedPoints,
-
     const scalar s0,
     const Type& p0,
-    const label p0Index,
+    const bool hasSnap0,
+    const Type& snapP0,
 
     const scalar s1,
     const Type& p1,
-    const label p1Index
+    const bool hasSnap1,
+    const Type& snapP1
 ) const
 {
     scalar d = s1-s0;
@@ -50,13 +50,13 @@ Type Foam::isoSurface::generatePoint
     {
         scalar s = (iso_-s0)/d;
 
-        if (s >= 0.5 && s <= 1 && p1Index != -1)
+        if (hasSnap1 && s >= 0.5 && s <= 1)
         {
-            return snappedPoints[p1Index];
+            return snapP1;
         }
-        else if (s >= 0.0 && s <= 0.5 && p0Index != -1)
+        else if (hasSnap0 && s >= 0.0 && s <= 0.5)
         {
-            return snappedPoints[p0Index];
+            return snapP0;
         }
         else
         {
@@ -75,23 +75,25 @@ Type Foam::isoSurface::generatePoint
 template<class Type>
 void Foam::isoSurface::generateTriPoints
 (
-    const DynamicList<Type>& snapped,
-
     const scalar s0,
     const Type& p0,
-    const label p0Index,
+    const bool hasSnap0,
+    const Type& snapP0,
 
     const scalar s1,
     const Type& p1,
-    const label p1Index,
+    const bool hasSnap1,
+    const Type& snapP1,
 
     const scalar s2,
     const Type& p2,
-    const label p2Index,
+    const bool hasSnap2,
+    const Type& snapP2,
 
     const scalar s3,
     const Type& p3,
-    const label p3Index,
+    const bool hasSnap3,
+    const Type& snapP3,
 
     DynamicList<Type>& points
 ) const
@@ -123,29 +125,55 @@ void Foam::isoSurface::generateTriPoints
 
         case 0x0E:
         case 0x01:
-            points.append(generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index));
-            points.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
-            points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
+            points.append
+            (
+                generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1)
+            );
+            points.append
+            (
+                generatePoint(s0,p0,hasSnap0,snapP0,s2,p2,hasSnap2,snapP2)
+            );
+            points.append
+            (
+                generatePoint(s0,p0,hasSnap0,snapP0,s3,p3,hasSnap3,snapP3)
+            );
         break;
 
         case 0x0D:
         case 0x02:
-            points.append(generatePoint(snapped,s1,p1,p1Index,s0,p0,p0Index));
-            points.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
-            points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
+            points.append
+            (
+                generatePoint(s1,p1,hasSnap1,snapP1,s0,p0,hasSnap0,snapP0)
+            );
+            points.append
+            (
+                generatePoint(s1,p1,hasSnap1,snapP1,s3,p3,hasSnap3,snapP3)
+            );
+            points.append
+            (
+                generatePoint(s1,p1,hasSnap1,snapP1,s2,p2,hasSnap2,snapP2)
+            );
         break;
 
         case 0x0C:
         case 0x03:
         {
-            Type tp1 = generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index);
-            Type tp2 = generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index);
+            Type tp1 =
+                generatePoint(s0,p0,hasSnap0,snapP0,s2,p2,hasSnap2,snapP2);
+            Type tp2 =
+                generatePoint(s1,p1,hasSnap1,snapP1,s3,p3,hasSnap3,snapP3);
 
-            points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
+            points.append
+            (
+                generatePoint(s0,p0,hasSnap0,snapP0,s3,p3,hasSnap3,snapP3)
+            );
             points.append(tp1);
             points.append(tp2);
             points.append(tp2);
-            points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
+            points.append
+            (
+                generatePoint(s1,p1,hasSnap1,snapP1,s2,p2,hasSnap2,snapP2)
+            );
             points.append(tp1);
         }
         break;
@@ -153,23 +181,40 @@ void Foam::isoSurface::generateTriPoints
         case 0x0B:
         case 0x04:
         {
-            points.append(generatePoint(snapped,s2,p2,p2Index,s0,p0,p0Index));
-            points.append(generatePoint(snapped,s2,p2,p2Index,s1,p1,p1Index));
-            points.append(generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index));
+            points.append
+            (
+                generatePoint(s2,p2,hasSnap2,snapP2,s0,p0,hasSnap0,snapP0)
+            );
+            points.append
+            (
+                generatePoint(s2,p2,hasSnap2,snapP2,s1,p1,hasSnap1,snapP1)
+            );
+            points.append
+            (
+                generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3)
+            );
         }
         break;
 
         case 0x0A:
         case 0x05:
         {
-            Type tp0 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
-            Type tp1 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
+            Type tp0 =
+                generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1);
+            Type tp1 =
+                generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3);
 
             points.append(tp0);
             points.append(tp1);
-            points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
+            points.append
+            (
+                generatePoint(s0,p0,hasSnap0,snapP0,s3,p3,hasSnap3,snapP3)
+            );
             points.append(tp0);
-            points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
+            points.append
+            (
+                generatePoint(s1,p1,hasSnap1,snapP1,s2,p2,hasSnap2,snapP2)
+            );
             points.append(tp1);
         }
         break;
@@ -177,23 +222,40 @@ void Foam::isoSurface::generateTriPoints
         case 0x09:
         case 0x06:
         {
-            Type tp0 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
-            Type tp1 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
+            Type tp0 =
+                generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1);
+            Type tp1 =
+                generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3);
 
             points.append(tp0);
-            points.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
+            points.append
+            (
+                generatePoint(s1,p1,hasSnap1,snapP1,s3,p3,hasSnap3,snapP3)
+            );
             points.append(tp1);
             points.append(tp0);
-            points.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
+            points.append
+            (
+                generatePoint(s0,p0,hasSnap0,snapP0,s2,p2,hasSnap2,snapP2)
+            );
             points.append(tp1);
         }
         break;
 
         case 0x07:
         case 0x08:
-            points.append(generatePoint(snapped,s3,p3,p3Index,s0,p0,p0Index));
-            points.append(generatePoint(snapped,s3,p3,p3Index,s2,p2,p2Index));
-            points.append(generatePoint(snapped,s3,p3,p3Index,s1,p1,p1Index));
+            points.append
+            (
+                generatePoint(s3,p3,hasSnap3,snapP3,s0,p0,hasSnap0,snapP0)
+            );
+            points.append
+            (
+                generatePoint(s3,p3,hasSnap3,snapP3,s2,p2,hasSnap2,snapP2)
+            );
+            points.append
+            (
+                generatePoint(s3,p3,hasSnap3,snapP3,s1,p1,hasSnap1,snapP1)
+            );
         break;
     }
 }
@@ -215,7 +277,8 @@ Foam::label Foam::isoSurface::generateFaceTriPoints
 
     const scalar neiVal,
     const Type& neiPt,
-    const label neiSnap,
+    const bool hasNeiSnap,
+    const Type& neiSnapPt,
 
     DynamicList<Type>& triPoints,
     DynamicList<label>& triMeshCells
@@ -234,23 +297,37 @@ Foam::label Foam::isoSurface::generateFaceTriPoints
 
         generateTriPoints
         (
-            snappedPoints,
-
             pVals[pointI],
             pCoords[pointI],
-            snappedPoint[pointI],
+            snappedPoint[pointI] != -1,
+            (
+                snappedPoint[pointI] != -1
+              ? snappedPoints[snappedPoint[pointI]]
+              : pTraits<Type>::zero
+            ),
 
             pVals[nextPointI],
             pCoords[nextPointI],
-            snappedPoint[nextPointI],
+            snappedPoint[nextPointI] != -1,
+            (
+                snappedPoint[nextPointI] != -1
+              ? snappedPoints[snappedPoint[nextPointI]]
+              : pTraits<Type>::zero
+            ),
 
             cVals[own],
             cCoords[own],
-            snappedCc[own],
+            snappedCc[own] != -1,
+            (
+                snappedCc[own] != -1
+              ? snappedPoints[snappedCc[own]]
+              : pTraits<Type>::zero
+            ),
 
             neiVal,
             neiPt,
-            neiSnap,
+            hasNeiSnap,
+            neiSnapPt,
 
             triPoints
         );
@@ -311,25 +388,6 @@ void Foam::isoSurface::generateTriPoints
             << abort(FatalError);
     }
 
-    // Determine neighbouring snap status
-    labelList neiSnappedCc(mesh_.nFaces()-mesh_.nInternalFaces(), -1);
-    forAll(patches, patchI)
-    {
-        const polyPatch& pp = patches[patchI];
-
-        if (pp.coupled())
-        {
-            label faceI = pp.start();
-            forAll(pp, i)
-            {
-                neiSnappedCc[faceI-mesh_.nInternalFaces()] =
-                    snappedCc[own[faceI]];
-                faceI++;
-            }
-        }
-    }
-    syncTools::swapBoundaryFaceList(mesh_, neiSnappedCc, false);
-
 
 
     // Generate triangle points
@@ -356,7 +414,12 @@ void Foam::isoSurface::generateTriPoints
 
                 cVals[nei[faceI]],
                 cCoords[nei[faceI]],
-                snappedCc[nei[faceI]],
+                snappedCc[nei[faceI]] != -1,
+                (
+                    snappedCc[nei[faceI]] != -1
+                  ? snappedPoints[snappedCc[nei[faceI]]]
+                  : pTraits<Type>::zero
+                ),
 
                 triPoints,
                 triMeshCells
@@ -365,6 +428,34 @@ void Foam::isoSurface::generateTriPoints
     }
 
 
+    // Determine neighbouring snap status
+    boolList neiSnapped(mesh_.nFaces()-mesh_.nInternalFaces(), false);
+    List<Type> neiSnappedPoint(neiSnapped.size(), pTraits<Type>::zero);
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        if (pp.coupled())
+        {
+            label faceI = pp.start();
+            forAll(pp, i)
+            {
+                label bFaceI = faceI-mesh_.nInternalFaces();
+                label snappedIndex = snappedCc[own[faceI]];
+
+                if (snappedIndex != -1)
+                {
+                    neiSnapped[bFaceI] = true;
+                    neiSnappedPoint[bFaceI] = snappedPoints[snappedIndex];
+                }
+                faceI++;
+            }
+        }
+    }
+    syncTools::swapBoundaryFaceList(mesh_, neiSnapped, false);
+    syncTools::swapBoundaryFaceList(mesh_, neiSnappedPoint, false);
+
+
     forAll(patches, patchI)
     {
         const polyPatch& pp = patches[patchI];
@@ -372,38 +463,53 @@ void Foam::isoSurface::generateTriPoints
         if
         (
             isA<processorPolyPatch>(pp)
-         && refCast<const processorPolyPatch>(pp).owner()
          && !refCast<const processorPolyPatch>(pp).separated()
         )
         {
-            label faceI = pp.start();
-
-            forAll(pp, i)
+            //if (refCast<const processorPolyPatch>(pp).owner())
             {
-                if (faceCutType_[faceI] != NOTCUT)
+                label faceI = pp.start();
+
+                forAll(pp, i)
                 {
-                    generateFaceTriPoints
-                    (
-                        cVals,
-                        pVals,
+                    if (faceCutType_[faceI] != NOTCUT)
+                    {
+                        label bFaceI = faceI-mesh_.nInternalFaces();
+                        if
+                        (
+                            neiSnapped[bFaceI]
+                         && (neiSnappedPoint[bFaceI]==pTraits<Type>::zero)
+                        )
+                        {
+                            FatalErrorIn("isoSurface::generateTriPoints(..)")
+                                << "problem:" << abort(FatalError);
+                        }
 
-                        cCoords,
-                        pCoords,
 
-                        snappedPoints,
-                        snappedCc,
-                        snappedPoint,
-                        faceI,
+                        generateFaceTriPoints
+                        (
+                            cVals,
+                            pVals,
 
-                        cVals.boundaryField()[patchI][i],
-                        cCoords.boundaryField()[patchI][i],
-                        neiSnappedCc[faceI-mesh_.nInternalFaces()],
+                            cCoords,
+                            pCoords,
 
-                        triPoints,
-                        triMeshCells
-                    );
+                            snappedPoints,
+                            snappedCc,
+                            snappedPoint,
+                            faceI,
+
+                            cVals.boundaryField()[patchI][i],
+                            cCoords.boundaryField()[patchI][i],
+                            neiSnapped[faceI-mesh_.nInternalFaces()],
+                            neiSnappedPoint[faceI-mesh_.nInternalFaces()],
+
+                            triPoints,
+                            triMeshCells
+                        );
+                    }
+                    faceI++;
                 }
-                faceI++;
             }
         }
         else if (isA<emptyPolyPatch>(pp))
@@ -430,7 +536,8 @@ void Foam::isoSurface::generateTriPoints
 
                         cVals[own[faceI]],
                         cCoords.boundaryField()[patchI][i],
-                        -1, // fc not snapped
+                        false,  // fc not snapped
+                        pTraits<Type>::zero,
 
                         triPoints,
                         triMeshCells
@@ -462,7 +569,8 @@ void Foam::isoSurface::generateTriPoints
 
                         cVals.boundaryField()[patchI][i],
                         cCoords.boundaryField()[patchI][i],
-                        -1, // fc not snapped
+                        false,  // fc not snapped
+                        pTraits<Type>::zero,
 
                         triPoints,
                         triMeshCells