diff --git a/meshLibrary/Make/files b/meshLibrary/Make/files
index 1fd86117efa65f40fe29307740687a99968ac524..b41b92a0e530f740e1738d6c3321c36dfd7813e2 100644
--- a/meshLibrary/Make/files
+++ b/meshLibrary/Make/files
@@ -5,8 +5,10 @@ meshSurfacePartitioner = utilities/surfaceTools/meshSurfacePartitioner
 boundaryFacesGenerator = utilities/surfaceTools/boundaryFacesGenerator
 boundaryLayers = utilities/boundaryLayers
 checkMeshDict = utilities/checkMeshDict
+
 meshSurfaceCheckInvertedVertices = utilities/surfaceTools/meshSurfaceCheckInvertedVertices
 meshSurfaceCheckEdgeTypes = utilities/surfaceTools/meshSurfaceCheckEdgeTypes
+meshSurfaceDetectPlanarRegions = utilities/surfaceTools/meshSurfaceDetectPlanarRegions
 meshSurfaceCutter  =  utilities/surfaceTools/meshSurfaceCutter
 meshSurfaceMapper  =  utilities/surfaceTools/meshSurfaceMapper
 meshSurfaceMapper2D = utilities/surfaceTools/meshSurfaceMapper2D
@@ -17,10 +19,12 @@ meshSurfaceEdgeExtractorNonTopo = utilities/surfaceTools/meshSurfaceEdgeExtracto
 meshSurfaceEdgeExtractorFUN = utilities/surfaceTools/meshSurfaceEdgeExtractorFUN
 meshSurfaceEdgeExtractor2D = utilities/surfaceTools/meshSurfaceEdgeExtractor2D
 correctEdgesBetweenPatches = utilities/surfaceTools/correctEdgesBetweenPatches
+
 createFundamentalSheetsBase = utilities/surfaceTools/createFundamentalSheets
 createFundamentalSheets = $(createFundamentalSheetsBase)/createFundamentalSheets
 createFundamentalSheetsJFS = $(createFundamentalSheetsBase)/createFundamentalSheetsJFS
 createFundamentalSheetsFJ = $(createFundamentalSheetsBase)/createFundamentalSheetsFJ
+
 decomposeCellsNearConcaveEdges = utilities/surfaceTools/decomposeCellsNearConcaveEdges
 renameBoundaryPatches = utilities/surfaceTools/renameBoundaryPatches
 
@@ -74,17 +78,18 @@ checkCellConnectionsOverFaces = $(topology)/checkCellConnectionsOverFaces
 checkIrregularSurfaceConnections = $(topology)/checkIrregularSurfaceConnections
 checkNonMappableCellConnections = $(topology)/checkNonMappableCellConnections
 
-triSurfaceCleanupDuplicates = utilities/triSurfaceTools/triSurfaceCleanupDuplicates
-triSurfaceCopyParts = utilities/triSurfaceTools/triSurfaceCopyParts
-triSurfaceCurvatureEstimator = utilities/triSurfaceTools/triSurfaceCurvatureEstimator
-triSurfacePartitioner = utilities/triSurfaceTools/triSurfacePartitioner
-triSurfaceDetectMaterials = utilities/triSurfaceTools/triSurfaceDetectMaterials
-triSurfaceDetectPlanarRegions = utilities/triSurfaceTools/triSurfaceDetectPlanarRegions
-triSurfaceDetectFeatureEdges = utilities/triSurfaceTools/triSurfaceDetectFeatureEdges
-triSurfaceClassifyEdges = utilities/triSurfaceTools/triSurfaceClassifyEdges
-triSurfacePatchManipulator = utilities/triSurfaceTools/triSurfacePatchManipulator
-triSurfaceRemoveFacets = utilities/triSurfaceTools/triSurfaceRemoveFacets
-triSurfaceExtrude2DEdges = utilities/triSurfaceTools/triSurfaceExtrude2DEdges
+triSurfaceTools = utilities/triSurfaceTools
+triSurfaceCleanupDuplicates = $(triSurfaceTools)/triSurfaceCleanupDuplicates
+triSurfaceCopyParts = $(triSurfaceTools)/triSurfaceCopyParts
+triSurfaceCurvatureEstimator = $(triSurfaceTools)/triSurfaceCurvatureEstimator
+triSurfacePartitioner = $(triSurfaceTools)/triSurfacePartitioner
+triSurfaceDetectMaterials = $(triSurfaceTools)/triSurfaceDetectMaterials
+triSurfaceDetectPlanarRegions = $(triSurfaceTools)/triSurfaceDetectPlanarRegions
+triSurfaceDetectFeatureEdges = $(triSurfaceTools)/triSurfaceDetectFeatureEdges
+triSurfaceClassifyEdges = $(triSurfaceTools)/triSurfaceClassifyEdges
+triSurfacePatchManipulator = $(triSurfaceTools)/triSurfacePatchManipulator
+triSurfaceRemoveFacets = $(triSurfaceTools)/triSurfaceRemoveFacets
+triSurfaceExtrude2DEdges = $(triSurfaceTools)/triSurfaceExtrude2DEdges
 
 polyMeshGen = utilities/meshes/polyMeshGen
 boundaryPatch  = utilities/meshes/polyMeshGen/boundaryPatch
@@ -287,6 +292,8 @@ $(meshSurfaceCheckInvertedVertices)/meshSurfaceCheckInvertedVertices.C
 
 $(meshSurfaceCheckEdgeTypes)/meshSurfaceCheckEdgeTypes.C
 
+$(meshSurfaceDetectPlanarRegions)/meshSurfaceDetectPlanarRegions.C
+
 $(meshSurfaceCutter)/meshSurfaceCutter.C
 $(meshSurfaceCutter)/meshSurfaceCutterFaceCutter.C
 $(meshSurfaceCutter)/meshSurfaceCutterInternalFaces.C
@@ -303,6 +310,7 @@ $(meshSurfaceMapper2D)/meshSurfaceMapper2DMapVertices.C
 $(meshSurfaceMapper2D)/meshSurfaceMapper2DPremapVertices.C
 
 $(edgeExtractor)/edgeExtractor.C
+$(edgeExtractor)/edgeExtractorCorners.C
 
 $(meshSurfaceEdgeExtractor)/meshSurfaceEdgeExtractor.C
 $(meshSurfaceEdgeExtractor)/meshSurfaceEdgeCreateEdgeVertices.C
diff --git a/meshLibrary/cartesianMesh/cartesianMeshGenerator/cartesianMeshGenerator.C b/meshLibrary/cartesianMesh/cartesianMeshGenerator/cartesianMeshGenerator.C
index 1baf2e6cf46b8885b2b74000ca0651dc7bb641be..62f69016c1886c13a5b35950c093b2c698e30656 100644
--- a/meshLibrary/cartesianMesh/cartesianMeshGenerator/cartesianMeshGenerator.C
+++ b/meshLibrary/cartesianMesh/cartesianMeshGenerator/cartesianMeshGenerator.C
@@ -163,7 +163,7 @@ void cartesianMeshGenerator::mapMeshToSurface()
     # endif
 
     //- untangle surface faces
-    meshSurfaceOptimizer(*msePtr, *octreePtr_).preOptimizeSurface();
+    meshSurfaceOptimizer(*msePtr, *octreePtr_).untangleSurface();
 
     # ifdef DEBUG
     mesh_.write();
@@ -196,34 +196,8 @@ void cartesianMeshGenerator::optimiseMeshSurface()
 
 void cartesianMeshGenerator::generateBoundaryLayers()
 {
+    //- add boundary layers
     boundaryLayers bl(mesh_);
-
-/*    if( meshDict_.found("boundaryLayers") )
-    {
-        if( )
-        wordList createLayers;
-
-        if( meshDict_.isDict("boundaryLayers") )
-        {
-            const dictionary& dict = meshDict_.subDict("boundaryLayers");
-            createLayers = dict.toc();
-        }
-        else
-        {
-            wordList bndLayers(meshDict_.lookup("boundaryLayers"));
-            createLayers.transfer(bndLayers);
-        }
-
-        forAll(createLayers, patchI)
-            bl.addLayerForPatch(createLayers[patchI]);
-    }
-    else
-    {
-        //bl.createOTopologyLayers();
-        bl.addLayerForAllPatches();
-    }
-    */
-
     bl.addLayerForAllPatches();
 
     # ifdef DEBUG
@@ -251,14 +225,16 @@ void cartesianMeshGenerator::refBoundaryLayers()
 
 void cartesianMeshGenerator::optimiseFinalMesh()
 {
-    //- final optimisation
-    meshOptimizer optimizer(mesh_);
-
-    optimizer.optimizeSurface(*octreePtr_);
-
+    //- untangle the surface if needed
+    meshSurfaceEngine mse(mesh_);
+    meshSurfaceOptimizer(mse, *octreePtr_).optimizeSurface();
     deleteDemandDrivenData(octreePtr_);
 
+    //- final optimisation
+    meshOptimizer optimizer(mesh_);
     optimizer.optimizeMeshFV();
+    optimizer.optimizeLowQualityFaces();
+    optimizer.untangleMeshFV();
 
     # ifdef DEBUG
     mesh_.write();
diff --git a/meshLibrary/hexMesh/hexMeshGenerator/hexMeshGenerator.C b/meshLibrary/hexMesh/hexMeshGenerator/hexMeshGenerator.C
index b3b7d8b1164b6b923c901d2938c96ac8d2c86b07..2000d729056d12b8a436fb67d0464e9bb24a5c91 100644
--- a/meshLibrary/hexMesh/hexMeshGenerator/hexMeshGenerator.C
+++ b/meshLibrary/hexMesh/hexMeshGenerator/hexMeshGenerator.C
@@ -120,7 +120,7 @@ void hexMeshGenerator::mapMeshToSurface()
     # endif
 
     //- untangle surface faces
-    meshSurfaceOptimizer(*msePtr, *octreePtr_).preOptimizeSurface();
+    meshSurfaceOptimizer(*msePtr, *octreePtr_).untangleSurface();
 
     # ifdef DEBUG
     mesh_.write();
diff --git a/meshLibrary/tetMesh/tetMeshGenerator/tetMeshGenerator.C b/meshLibrary/tetMesh/tetMeshGenerator/tetMeshGenerator.C
index 70aadd056a2ac9f6b314c59fd70c3f83d9a4e279..06020bf34539fc0db76fec2f059fe349840945f3 100644
--- a/meshLibrary/tetMesh/tetMeshGenerator/tetMeshGenerator.C
+++ b/meshLibrary/tetMesh/tetMeshGenerator/tetMeshGenerator.C
@@ -90,6 +90,7 @@ void tetMeshGenerator::surfacePreparation()
 
 void tetMeshGenerator::mapMeshToSurface()
 {
+    return;
     //- calculate mesh surface
     meshSurfaceEngine* msePtr = new meshSurfaceEngine(mesh_);
 
@@ -102,7 +103,7 @@ void tetMeshGenerator::mapMeshToSurface()
     # endif
 
     //- untangle surface faces
-    meshSurfaceOptimizer(*msePtr, *octreePtr_).preOptimizeSurface();
+    meshSurfaceOptimizer(*msePtr, *octreePtr_).untangleSurface();
 
     # ifdef DEBUG
     mesh_.write();
diff --git a/meshLibrary/utilities/boundaryLayers/boundaryLayerCells.C b/meshLibrary/utilities/boundaryLayers/boundaryLayerCells.C
index 688fcc3e1804442235b72cb6bb6d8ef39369560e..c8d5bae310b8fff4fe4039ea988fd0d4c4676d4b 100644
--- a/meshLibrary/utilities/boundaryLayers/boundaryLayerCells.C
+++ b/meshLibrary/utilities/boundaryLayers/boundaryLayerCells.C
@@ -761,14 +761,20 @@ void boundaryLayers::createNewFacesFromPointsParallel
                 const label fI = pointFaceCandidates(pointI, i);
                 const DynList<label, 4>& f = faceCandidates[fI];
 
-                Pair<label> pk
+                const labelPair pk
                 (
                     patchKey_[candidatePatches[fI][0]],
                     patchKey_[candidatePatches[fI][1]]
                 );
 
+                const labelPair rpk
+                (
+                    patchKey_[candidatePatches[fI][1]],
+                    patchKey_[candidatePatches[fI][0]]
+                );
+
                 if(
-                    (f[0] == pointI) && ((pk == lp) || (pk.reversePair() == lp))
+                    (f[0] == pointI) && ((pk == lp) || (rpk == lp))
                 )
                 {
                     //- found the processor containing other face
diff --git a/meshLibrary/utilities/helperClasses/parallelHelpers/labelledPair/labelledPair.H b/meshLibrary/utilities/helperClasses/parallelHelpers/labelledPair/labelledPair.H
index 4b2258871b4393cb512669940980c13d4d8ff975..b10e89a53a7b09b05b4cd7c55990690f7df7270d 100644
--- a/meshLibrary/utilities/helperClasses/parallelHelpers/labelledPair/labelledPair.H
+++ b/meshLibrary/utilities/helperClasses/parallelHelpers/labelledPair/labelledPair.H
@@ -47,18 +47,18 @@ namespace Foam
 /*---------------------------------------------------------------------------*\
                 Class labelledPair Declaration
 \*---------------------------------------------------------------------------*/
-    
+
 class labelledPair
 {
     // Private data
         //- label
         label pLabel_;
-    
+
         //- pair data
         labelPair pair_;
-    
+
     public:
-        
+
     // Constructors
         //- Null construct
         labelledPair()
@@ -66,7 +66,7 @@ class labelledPair
             pLabel_(-1),
             pair_()
         {}
-        
+
         //- Construct from label and pair
         labelledPair
         (
@@ -77,32 +77,32 @@ class labelledPair
             pLabel_(pl),
             pair_(lp)
         {}
-            
+
     // Destructor
         ~labelledPair()
         {}
-            
+
     // Member functions
         //- return pair label
         inline label pairLabel() const
         {
             return pLabel_;
         }
-        
+
         //- return the pair
         inline const labelPair& pair() const
         {
             return pair_;
         }
-        
+
     // Member operators
-        
+
         inline void operator=(const labelledPair& lpp)
         {
             pLabel_ = lpp.pLabel_;
             pair_ = lpp.pair_;
         }
-        
+
         inline bool operator==
         (
             const labelledPair& lpp
@@ -110,13 +110,17 @@ class labelledPair
         {
             if( pLabel_ != lpp.pLabel_ )
                 return false;
-            
-            if( (pair_ != lpp.pair_) && (pair_.reversePair() != lpp.pair_) )
-                return false;
-            
-            return true;
+
+            const labelPair& op = lpp.pair_;
+
+            if( (pair_.first() == op.first()) && (pair_.second() == op.second()) )
+                return true;
+            if( (pair_.first() == op.second()) && (pair_.second() == op.first()) )
+                return true;
+
+            return false;
         }
-        
+
         inline bool operator!=
         (
             const labelledPair& lcc
@@ -124,7 +128,7 @@ class labelledPair
         {
             return !this->operator==(lcc);
         }
-        
+
         inline bool operator<(const labelledPair& lpp) const
         {
             if( pLabel_ < lpp.pLabel_ )
@@ -135,22 +139,22 @@ class labelledPair
             {
                 return false;
             }
-            
+
             if(
                 (pair_.first() + pair_.second()) <
                 (lpp.pair().first() + lpp.pair().second())
             )
                 return true;
-            
+
             if(
                 Foam::min(pair_.first(), pair_.second()) <
                 Foam::min(lpp.pair().first(), lpp.pair().second())
             )
                 return true;
-            
+
             return false;
         }
-        
+
     // Friend operators
         friend Ostream& operator<<
         (
@@ -170,7 +174,7 @@ class labelledPair
 
             return os;
         }
-        
+
         friend Istream& operator>>
         (
             Istream& is,
@@ -179,16 +183,16 @@ class labelledPair
         {
             // Read beginning of labelledPair
             is.readBegin("labelledPair");
-        
+
             is >> lpp.pLabel_;
             is >> lpp.pair_;
-            
+
             // Read end of labelledPair
             is.readEnd("labelledPair");
-        
+
             // Check state of Istream
             is.check("operator>>(Istream&, labelledPair");
-            
+
             return is;
         }
 };
diff --git a/meshLibrary/utilities/octrees/meshOctree/meshOctree.H b/meshLibrary/utilities/octrees/meshOctree/meshOctree.H
index 47790d1265ac24eaffb35d3d04281f9a4870ce56..a33fca7cf1b380f4cee6fac0984f0091e08b7440 100644
--- a/meshLibrary/utilities/octrees/meshOctree/meshOctree.H
+++ b/meshLibrary/utilities/octrees/meshOctree/meshOctree.H
@@ -164,6 +164,7 @@ public:
         (
             point& nearest,
             scalar& distSq,
+            label& nearestTriangle,
             label& region,
             const point& p
         ) const;
@@ -173,6 +174,7 @@ public:
         (
             point& nearest,
             scalar& distSq,
+            label& nearestTriangle,
             const label region,
             const point& p
         ) const;
@@ -180,18 +182,20 @@ public:
         //- find nearest feature-edges vertex to a given vertex
         bool findNearestEdgePoint
         (
-            const point& p,
-            const DynList<label>& regions,
             point& edgePoint,
-            scalar& distSq
+            scalar& distSq,
+            label& nearestEdge,
+            const point& p,
+            const DynList<label>& regions
         ) const;
 
-        bool findNearestVertexToTheEdge
+        bool findNearestPointToEdge
         (
-            const FixedList<point, 2>& edgePoints,
-            const FixedList<label, 2>& edgePointRegions,
             point& nearest,
-            scalar& distSq
+            scalar& distSq,
+            label& nearestEdge,
+            const FixedList<point, 2>& edgePoints,
+            const FixedList<label, 2>& edgePointRegions
         ) const;
 
         //- find nearest corner point
@@ -199,12 +203,13 @@ public:
         (
             point& nearest,
             scalar& distSq,
+            label& nearestPoint,
             const point& p,
             const DynList<label>& patches
         ) const;
 
         //- find the nearest vertex to the given patches
-        bool findNearestVertexToPatches
+        bool findNearestPointToPatches
         (
             point& nearest,
             scalar& distSq,
diff --git a/meshLibrary/utilities/octrees/meshOctree/meshOctreeFindNearestSurfacePoint.C b/meshLibrary/utilities/octrees/meshOctree/meshOctreeFindNearestSurfacePoint.C
index adeb2f53adbff2a828c46c57587a4571290afa6e..494a7cb1c283ce837105889d194af4502072ca01 100644
--- a/meshLibrary/utilities/octrees/meshOctree/meshOctreeFindNearestSurfacePoint.C
+++ b/meshLibrary/utilities/octrees/meshOctree/meshOctreeFindNearestSurfacePoint.C
@@ -45,11 +45,13 @@ void meshOctree::findNearestSurfacePoint
 (
     point& nearest,
     scalar& distSq,
+    label& nearestTriangle,
     label& region,
     const point& p
 ) const
 {
     region = -1;
+    nearestTriangle = 1;
 
     const label cLabel = findLeafContainingVertex(p);
     vector sizeVec;
@@ -96,6 +98,7 @@ void meshOctree::findNearestSurfacePoint
                 {
                     distSq = dSq;
                     nearest = p0;
+                    nearestTriangle = el[tI];
                     region = surface_[el[tI]].region();
                     found = true;
                 }
@@ -129,6 +132,7 @@ void meshOctree::findNearestSurfacePointInRegion
 (
     point& nearest,
     scalar& distSq,
+    label& nearestTriangle,
     const label region,
     const point& p
 ) const
@@ -149,12 +153,13 @@ void meshOctree::findNearestSurfacePointInRegion
     bool found(false);
     label iterationI(0);
     DynList<const meshOctreeCube*, 256> neighbours;
+    nearestTriangle = -1;
 
     distSq = VGREAT;
 
     do
     {
-        boundBox bb(p - sizeVec, p + sizeVec);
+        const boundBox bb(p - sizeVec, p + sizeVec);
 
         neighbours.clear();
         findLeavesContainedInBox(bb, neighbours);
@@ -182,6 +187,7 @@ void meshOctree::findNearestSurfacePointInRegion
                 {
                     distSq = dSq;
                     nearest = p0;
+                    nearestTriangle = el[tI];
                     found = true;
                 }
             }
@@ -198,10 +204,11 @@ void meshOctree::findNearestSurfacePointInRegion
 
 bool meshOctree::findNearestEdgePoint
 (
-    const point& p,
-    const DynList<label>& regions,
     point& edgePoint,
-    scalar& distSq
+    scalar& distSq,
+    label& nearestEdge,
+    const point& p,
+    const DynList<label>& regions
 ) const
 {
     //- find the estimate for the searching range
@@ -228,6 +235,7 @@ bool meshOctree::findNearestEdgePoint
     label iterationI(0);
 
     distSq = VGREAT;
+    nearestEdge = -1;
 
     do
     {
@@ -274,6 +282,7 @@ bool meshOctree::findNearestEdgePoint
                 {
                     distSq = dSq;
                     edgePoint = np;
+                    nearestEdge = ce[eI];
                     foundAnEdge = true;
                 }
             }
@@ -287,12 +296,13 @@ bool meshOctree::findNearestEdgePoint
     return foundAnEdge;
 }
 
-bool meshOctree::findNearestVertexToTheEdge
+bool meshOctree::findNearestPointToEdge
 (
-    const FixedList<point, 2>& edgePoints,
-    const FixedList<label, 2>& edgePointRegions,
     point& nearest,
-    scalar& distSq
+    scalar& distSq,
+    label& nearestEdge,
+    const FixedList<point, 2>& edgePoints,
+    const FixedList<label, 2>& edgePointRegions
 ) const
 {
     const point c = 0.5 * (edgePoints[0] + edgePoints[1]);
@@ -309,6 +319,7 @@ bool meshOctree::findNearestVertexToTheEdge
     const edgeLongList& surfaceEdges = surface_.edges();
 
     distSq = VGREAT;
+    nearestEdge = -1;
 
     bool found(false);
 
@@ -357,6 +368,7 @@ bool meshOctree::findNearestVertexToTheEdge
                     if( magSqr(nearestOnEdge - nearestOnLine) < distSq )
                     {
                         nearest = nearestOnEdge;
+                        nearestEdge = edges[eI];
                         distSq = magSqr(nearestOnEdge - nearestOnLine);
                         found = true;
                     }
@@ -372,6 +384,7 @@ bool meshOctree::findNearestCorner
 (
     point& nearest,
     scalar& distSq,
+    label& nearestPoint,
     const point& p,
     const DynList<label>& patches
 ) const
@@ -400,6 +413,7 @@ bool meshOctree::findNearestCorner
     const VRWGraph& eFacets = surface_.edgeFacets();
 
     distSq = VGREAT;
+    nearestPoint = -1;
 
     do
     {
@@ -438,7 +452,7 @@ bool meshOctree::findNearestCorner
                     {
                         const label eI = pEdges(spI, i);
 
-                        if( pEdges.sizeOfRow(eI) < 2 )
+                        if( eFacets.sizeOfRow(eI) != 2 )
                             break;
 
                         if(
@@ -480,6 +494,7 @@ bool meshOctree::findNearestCorner
                                 distSq = dSq;
                                 found = true;
                                 nearest = points[spI];
+                                nearestPoint = spI;
                             }
                         }
                     }
@@ -495,7 +510,7 @@ bool meshOctree::findNearestCorner
     return found;
 }
 
-bool meshOctree::findNearestVertexToPatches
+bool meshOctree::findNearestPointToPatches
 (
     point& nearest,
     scalar& distSq,
@@ -507,6 +522,21 @@ bool meshOctree::findNearestVertexToPatches
     if( patches.size() == 0 )
         return false;
 
+    point mapPoint(p);
+    scalar dSq(VGREAT);
+
+    bool found(false);
+    if( patches.size() == 2 )
+    {
+        label nse;
+        found = findNearestEdgePoint(mapPoint, dSq, nse, p, patches);
+    }
+    else if( patches.size() > 2 )
+    {
+        label nsp;
+        found = findNearestCorner(mapPoint, dSq, nsp, p, patches);
+    }
+
     point mapPointApprox(p);
     scalar distSqApprox;
     label iter(0);
@@ -516,10 +546,12 @@ bool meshOctree::findNearestVertexToPatches
         forAll(patches, patchI)
         {
             point np;
+            label nearestTri;
             this->findNearestSurfacePointInRegion
             (
                 np,
                 distSqApprox,
+                nearestTri,
                 patches[patchI],
                 mapPointApprox
             );
@@ -528,7 +560,7 @@ bool meshOctree::findNearestVertexToPatches
         }
 
         newP /= patches.size();
-        if( Foam::magSqr(newP - mapPointApprox) < tol * searchRange_ )
+        if( Foam::magSqr(newP - mapPointApprox) < tol * dSq )
             break;
 
         mapPointApprox = newP;
@@ -537,6 +569,12 @@ bool meshOctree::findNearestVertexToPatches
     distSq = Foam::magSqr(mapPointApprox - p);
     nearest = mapPointApprox;
 
+    if( found && (dSq < 1.5 * distSq) )
+    {
+        nearest = mapPoint;
+        distSq = dSq;
+    }
+
     return true;
 }
 
diff --git a/meshLibrary/utilities/octrees/surfaceIntersectionsOctree/surfaceIntersectionsOctree.C b/meshLibrary/utilities/octrees/surfaceIntersectionsOctree/surfaceIntersectionsOctree.C
index 1eff2e43ae626ccb2f4034ccea7a38997a850d08..fb1515330b48fc2285feb8b0f44178e4ff0c8d3f 100644
--- a/meshLibrary/utilities/octrees/surfaceIntersectionsOctree/surfaceIntersectionsOctree.C
+++ b/meshLibrary/utilities/octrees/surfaceIntersectionsOctree/surfaceIntersectionsOctree.C
@@ -29,6 +29,7 @@ Description
 #include "surfaceIntersectionsOctree.H"
 #include "triSurf.H"
 #include "boundBox.H"
+#include "SLList.H"
 
 // #define DEBUGSearch
 
@@ -89,15 +90,12 @@ surfaceIntersectionsOctree::surfaceIntersectionsOctree
     //- in a single octree cube
     initialCube_.refineTree(maxN, k);
 
-    SLList<surfaceIntersectionsOctreeCube*> leaves;
+    LongList<surfaceIntersectionsOctreeCube*> leaves;
     findLeavesForCube(&initialCube_, leaves);
 
-    for(SLList<surfaceIntersectionsOctreeCube*>::iterator cIter = leaves.begin();
-        cIter != leaves.end();
-        ++cIter
-    )
+    forAll(leaves, leafI)
     {
-        surfaceIntersectionsOctreeCube& oc = *cIter();
+        surfaceIntersectionsOctreeCube& oc = *leaves[leafI];
 
         if( oc.containedElements().size() > 0 )
         {
@@ -132,7 +130,7 @@ surfaceIntersectionsOctree::~surfaceIntersectionsOctree()
 void surfaceIntersectionsOctree::findLeavesForCube
 (
     surfaceIntersectionsOctreeCube* oc,
-    SLList<surfaceIntersectionsOctreeCube*>& lvs
+    LongList<surfaceIntersectionsOctreeCube*>& lvs
 ) const
 {
     if( oc->isLeaf() )
@@ -151,7 +149,8 @@ void surfaceIntersectionsOctree::findLeavesForCube
     }
 }
 
-surfaceIntersectionsOctreeCube* surfaceIntersectionsOctree::findLeafContainingVertex(const point& p) const
+surfaceIntersectionsOctreeCube*
+surfaceIntersectionsOctree::findLeafContainingVertex(const point& p) const
 {
     # ifdef OCTREE_DEBUG
     Info << "Finding leaf for vertex " << p << endl;
@@ -175,7 +174,8 @@ surfaceIntersectionsOctreeCube* surfaceIntersectionsOctree::findLeafContainingVe
         if( !oc->isLeaf() )
         {
             //- find a subCube containing the vertex;
-            const FixedList<surfaceIntersectionsOctreeCube*, 8>& sc = *oc->subCubes();
+            const FixedList<surfaceIntersectionsOctreeCube*, 8>& sc =
+                *oc->subCubes();
 
             bool found(false);
 
@@ -193,8 +193,8 @@ surfaceIntersectionsOctreeCube* surfaceIntersectionsOctree::findLeafContainingVe
             if( !found )
                 FatalErrorIn
                 (
-                    "surfaceIntersectionsOctreeCube* meshOctree::findLeafContainingVertex"
-                    "(const point& p) const"
+                    "surfaceIntersectionsOctreeCube* meshOctree::"
+                    "findLeafContainingVertex(const point& p) const"
                 ) << "Vertex is not found in any subCubes!!"
                     << abort(FatalError);
         }
diff --git a/meshLibrary/utilities/octrees/surfaceIntersectionsOctree/surfaceIntersectionsOctree.H b/meshLibrary/utilities/octrees/surfaceIntersectionsOctree/surfaceIntersectionsOctree.H
index 5997b586cf233e7eeeb32e46bc099eb4bcc8a512..d5d7aa9d507963a0371a13aeba73ffbd34db748a 100644
--- a/meshLibrary/utilities/octrees/surfaceIntersectionsOctree/surfaceIntersectionsOctree.H
+++ b/meshLibrary/utilities/octrees/surfaceIntersectionsOctree/surfaceIntersectionsOctree.H
@@ -36,7 +36,7 @@ SourceFiles
 #ifndef surfaceIntersectionsOctree_H
 #define surfaceIntersectionsOctree_H
 
-#include "SLList.H"
+#include "LongList.H"
 #include "boolList.H"
 #include "pointIndexHit.H"
 #include "surfaceIntersectionsOctreeCube.H"
@@ -90,7 +90,7 @@ class surfaceIntersectionsOctree
         void findLeavesForCube
         (
             surfaceIntersectionsOctreeCube*,
-            SLList<surfaceIntersectionsOctreeCube*>&
+            LongList<surfaceIntersectionsOctreeCube*>&
         ) const;
 
         //- Disallow default bitwise copy construct
@@ -142,7 +142,10 @@ public:
         ) const;
 
         //- find a cube containing the vertex
-        surfaceIntersectionsOctreeCube* findLeafContainingVertex(const point&) const;
+        surfaceIntersectionsOctreeCube* findLeafContainingVertex
+        (
+            const point&
+        ) const;
 };
 
 
diff --git a/meshLibrary/utilities/octrees/surfaceIntersectionsOctree/surfaceIntersectionsOctreeIntersections.C b/meshLibrary/utilities/octrees/surfaceIntersectionsOctree/surfaceIntersectionsOctreeIntersections.C
index e805a7e6f74d67411beb874762b06a54c968624c..e67631681febc2e51d1ed12cff0eaa0440e7cc38 100644
--- a/meshLibrary/utilities/octrees/surfaceIntersectionsOctree/surfaceIntersectionsOctreeIntersections.C
+++ b/meshLibrary/utilities/octrees/surfaceIntersectionsOctree/surfaceIntersectionsOctreeIntersections.C
@@ -31,6 +31,7 @@ Description
 #include "boundBox.H"
 #include "DynList.H"
 #include "helperFunctions.H"
+#include "SLList.H"
 
 // #define DEBUGSearch
 
diff --git a/meshLibrary/utilities/smoothers/geometry/meshOptimizer/meshOptimizer.C b/meshLibrary/utilities/smoothers/geometry/meshOptimizer/meshOptimizer.C
index a3e291085836d8db1fa1a20aab29131730396a9c..15feeb158c897306dcc6f62705ea7867e1f0de0e 100644
--- a/meshLibrary/utilities/smoothers/geometry/meshOptimizer/meshOptimizer.C
+++ b/meshLibrary/utilities/smoothers/geometry/meshOptimizer/meshOptimizer.C
@@ -39,14 +39,14 @@ Description
 
 namespace Foam
 {
-    
+
 // * * * * * * * * Private member functions * * * * * * * * * * * * * * * * * //
-    
+
 const meshSurfaceEngine& meshOptimizer::meshSurface() const
 {
     if( !msePtr_ )
         msePtr_ = new meshSurfaceEngine(mesh_);
-    
+
     return *msePtr_;
 }
 
@@ -62,7 +62,7 @@ label meshOptimizer::findBadFaces
 ) const
 {
     badFaces.clear();
-    
+
     polyMeshGenChecks::checkFacePyramids
     (
         mesh_,
@@ -71,7 +71,7 @@ label meshOptimizer::findBadFaces
         &badFaces,
         &changedFace
     );
-    
+
     polyMeshGenChecks::checkFaceFlatness
     (
         mesh_,
@@ -80,7 +80,7 @@ label meshOptimizer::findBadFaces
         &badFaces,
         &changedFace
     );
-    
+
     polyMeshGenChecks::checkCellPartTetrahedra
     (
         mesh_,
@@ -89,7 +89,7 @@ label meshOptimizer::findBadFaces
         &badFaces,
         &changedFace
     );
-    
+
     polyMeshGenChecks::checkFaceAreas
     (
         mesh_,
@@ -98,9 +98,9 @@ label meshOptimizer::findBadFaces
         &badFaces,
         &changedFace
     );
-    
+
     const label nBadFaces = returnReduce(badFaces.size(), sumOp<label>());
-    
+
     return nBadFaces;
 }
 
@@ -111,7 +111,7 @@ label meshOptimizer::findLowQualityFaces
 ) const
 {
     badFaces.clear();
-    
+
     polyMeshGenChecks::checkFaceDotProduct
     (
         mesh_,
@@ -119,7 +119,7 @@ label meshOptimizer::findLowQualityFaces
         70.0,
         &badFaces
     );
-    
+
     polyMeshGenChecks::checkFaceSkewness
     (
         mesh_,
@@ -127,9 +127,9 @@ label meshOptimizer::findLowQualityFaces
         2.0,
         &badFaces
     );
-    
+
     const label nBadFaces = returnReduce(badFaces.size(), sumOp<label>());
-    
+
     return nBadFaces;
 }
 
@@ -144,25 +144,25 @@ meshOptimizer::meshOptimizer(polyMeshGen& mesh)
 {
     const meshSurfaceEngine& mse = meshSurface();
     const labelList& bPoints = mse.boundaryPoints();
-    
+
     //- mark boundary vertices
     forAll(bPoints, bpI)
         vertexLocation_[bPoints[bpI]] = BOUNDARY;
-    
+
     //- mark edge vertices
     meshSurfacePartitioner mPart(mse);
-    forAllConstIter(labelHashSet, mPart.edgeNodes(), it)
+    forAllConstIter(labelHashSet, mPart.edgePoints(), it)
         vertexLocation_[bPoints[it.key()]] = EDGE;
-    
+
     //- mark corner vertices
     forAllConstIter(labelHashSet, mPart.corners(), it)
         vertexLocation_[bPoints[it.key()]] = CORNER;
-    
+
     if( Pstream::parRun() )
     {
         const polyMeshGenAddressing& addresing = mesh_.addressingData();
         const VRWGraph& pointAtProcs = addresing.pointAtProcs();
-        
+
         forAll(pointAtProcs, pointI)
             if( pointAtProcs.sizeOfRow(pointI) != 0 )
                 vertexLocation_[pointI] |= PARALLELBOUNDARY;
diff --git a/meshLibrary/utilities/smoothers/geometry/meshOptimizer/optimizeMeshFV.C b/meshLibrary/utilities/smoothers/geometry/meshOptimizer/optimizeMeshFV.C
index fb332646a06c65c8a8f2bb9f3cc4204b83c70c4e..252e84f078560925e859fbb8874787afd2e3c173 100644
--- a/meshLibrary/utilities/smoothers/geometry/meshOptimizer/optimizeMeshFV.C
+++ b/meshLibrary/utilities/smoothers/geometry/meshOptimizer/optimizeMeshFV.C
@@ -235,8 +235,7 @@ void meshOptimizer::optimizeMeshFV()
     Info << "Starting smoothing the mesh" << endl;
 
     laplaceSmoother lps(mesh_, vertexLocation_);
-    lps.optimizeLaplacianPC(2);
-    lps.optimizeLaplacianWPC(3);
+    lps.optimizeLaplacianPC(5);
 
     untangleMeshFV();
 
diff --git a/meshLibrary/utilities/smoothers/geometry/meshSurfaceOptimizer/meshSurfaceOptimizer.C b/meshLibrary/utilities/smoothers/geometry/meshSurfaceOptimizer/meshSurfaceOptimizer.C
index e9d12f51f1319fa103071988b570e52f43b84401..b16986165cb0b7f21ceb36e26ace329395d5887e 100644
--- a/meshLibrary/utilities/smoothers/geometry/meshSurfaceOptimizer/meshSurfaceOptimizer.C
+++ b/meshLibrary/utilities/smoothers/geometry/meshSurfaceOptimizer/meshSurfaceOptimizer.C
@@ -44,7 +44,7 @@ namespace Foam
 void meshSurfaceOptimizer::classifySurfaceVertices()
 {
     const labelHashSet& corners = partitioner_.corners();
-    const labelHashSet& edgePoints = partitioner_.edgeNodes();
+    const labelHashSet& edgePoints = partitioner_.edgePoints();
 
     //- set all vertices to partition
     vertexType_ = PARTITION;
diff --git a/meshLibrary/utilities/smoothers/geometry/meshSurfaceOptimizer/meshSurfaceOptimizer.H b/meshLibrary/utilities/smoothers/geometry/meshSurfaceOptimizer/meshSurfaceOptimizer.H
index 51fadf905a6a6d0e91ad592f3b95ccef590e9da9..73c4fe3955296e4bf4c8203f27f19d0caa798d24 100644
--- a/meshLibrary/utilities/smoothers/geometry/meshSurfaceOptimizer/meshSurfaceOptimizer.H
+++ b/meshLibrary/utilities/smoothers/geometry/meshSurfaceOptimizer/meshSurfaceOptimizer.H
@@ -258,8 +258,19 @@ public:
         ~meshSurfaceOptimizer();
 
     // Member Functions
-        //- optimisation before the edges and corners are created
-        bool preOptimizeSurface(const label nAdditionalLayers = 2);
+        //- runs a surface smoother on the selected boundary points
+        bool untangleSurface
+        (
+            const labelLongList& activeBoundaryPoints,
+            const label nAdditionalLayers = 0
+        );
+
+        //- checks for invertex surface elements and tries to untangle them
+        //- it tries to keep the points on the surface for a couple
+        //- of iteration and gives up the final iterations
+        //- by default, it smooths two additional layers of elements
+        //- around the inverted ones
+        bool untangleSurface(const label nAdditionalLayers = 2);
 
         //- optimize boundary nodes after boundary regions are created
         void optimizeSurface(const label nIterations = 5);
diff --git a/meshLibrary/utilities/smoothers/geometry/meshSurfaceOptimizer/meshSurfaceOptimizerOptimizeSurface.C b/meshLibrary/utilities/smoothers/geometry/meshSurfaceOptimizer/meshSurfaceOptimizerOptimizeSurface.C
index 2cdd72d6a4a18ef01d776daa76a275921a78953e..f2afa44f93b85705a338dc08034ebe5a6cbfbfd4 100644
--- a/meshLibrary/utilities/smoothers/geometry/meshSurfaceOptimizer/meshSurfaceOptimizerOptimizeSurface.C
+++ b/meshLibrary/utilities/smoothers/geometry/meshSurfaceOptimizer/meshSurfaceOptimizerOptimizeSurface.C
@@ -153,9 +153,13 @@ label meshSurfaceOptimizer::findInvertedVertices
     return nInvertedTria;
 }
 
-bool meshSurfaceOptimizer::preOptimizeSurface(const label nLayers)
+bool meshSurfaceOptimizer::untangleSurface
+(
+    const labelLongList& selectedBoundaryPoints,
+    const label nAdditionalLayers
+)
 {
-    Info << "Optimizing positions of surface nodes" << endl;
+    Info << "Starting untangling the surface of the volume mesh" << endl;
 
     bool changed(false);
 
@@ -177,7 +181,12 @@ bool meshSurfaceOptimizer::preOptimizeSurface(const label nLayers)
         surfaceEngine_.bpNeiProcs();
     }
 
-    boolList smoothVertex;
+    boolList smoothVertex(bPoints.size(), false);
+    # ifdef USE_OMP
+    # pragma omp parallel for schedule(dynamic, 50)
+    # endif
+    forAll(selectedBoundaryPoints, i)
+        smoothVertex[selectedBoundaryPoints[i]] = true;
 
     meshSurfaceEngineModifier surfaceModifier(surfaceEngine_);
     meshSurfaceMapper mapper(surfaceEngine_, meshOctree_);
@@ -186,7 +195,8 @@ bool meshSurfaceOptimizer::preOptimizeSurface(const label nLayers)
     label nInvertedTria;
     label nGlobalIter(0);
 
-    labelLongList procBndNodes, movedPoints;
+    labelLongList procBndPoints, movedPoints;
+    labelLongList procEdgePoints, movedEdgePoints;
 
     do
     {
@@ -194,28 +204,81 @@ bool meshSurfaceOptimizer::preOptimizeSurface(const label nLayers)
 
         do
         {
-            nInvertedTria = findInvertedVertices(smoothVertex, nLayers);
+            nInvertedTria = findInvertedVertices(smoothVertex, nAdditionalLayers);
 
             if( nInvertedTria == 0 ) break;
 
             changed = true;
 
-            procBndNodes.clear();
+            procBndPoints.clear();
             movedPoints.clear();
+            procEdgePoints.clear();
+            movedEdgePoints.clear();
+
             forAll(bPoints, bpI)
             {
-                if( smoothVertex[bpI] && (vertexType_[bpI] & PARTITION) )
+                if( !smoothVertex[bpI] )
+                    continue;
+
+                if( vertexType_[bpI] & PARTITION )
                 {
                     movedPoints.append(bpI);
 
                     if( vertexType_[bpI] & PROCBND )
                     {
-                        procBndNodes.append(bpI);
+                        procBndPoints.append(bpI);
+                        continue;
+                    }
+                }
+                else if( vertexType_[bpI] & EDGE )
+                {
+                    movedEdgePoints.append(bpI);
+
+                    if( vertexType_[bpI] & PROCBND )
+                    {
+                        procEdgePoints.append(bpI);
                         continue;
                     }
                 }
             }
 
+            //- smooth edge vertices
+            # ifdef USE_OMP
+            # pragma omp parallel
+            # endif
+            {
+                LongList<labelledPoint> newPos;
+
+                # ifdef USE_OMP
+                # pragma omp for schedule(dynamic, 40)
+                # endif
+                forAll(movedEdgePoints, i)
+                {
+                    const label bpI = movedEdgePoints[i];
+
+                    if( vertexType_[bpI] & PROCBND )
+                        continue;
+
+                    newPos.append
+                    (
+                        labelledPoint(bpI, newEdgePositionLaplacian(bpI))
+                    );
+                }
+
+                forAll(newPos, i)
+                    surfaceModifier.moveBoundaryVertexNoUpdate
+                    (
+                        newPos[i].pointLabel(),
+                        newPos[i].coordinates()
+                    );
+            }
+
+            if( Pstream::parRun() )
+                edgeNodeDisplacementParallel(procEdgePoints);
+
+            //- update normals and other geometric data
+            surfaceModifier.updateGeometry(movedEdgePoints);
+
             //- use laplacian smoothing
             # ifdef USE_OMP
             # pragma omp parallel
@@ -248,7 +311,7 @@ bool meshSurfaceOptimizer::preOptimizeSurface(const label nLayers)
             }
 
             if( Pstream::parRun() )
-                nodeDisplacementLaplacianFCParallel(procBndNodes, true);
+                nodeDisplacementLaplacianFCParallel(procBndPoints, true);
 
             //- update normals and other geometric data
             surfaceModifier.updateGeometry(movedPoints);
@@ -285,7 +348,7 @@ bool meshSurfaceOptimizer::preOptimizeSurface(const label nLayers)
             }
 
             if( Pstream::parRun() )
-                nodeDisplacementSurfaceOptimizerParallel(procBndNodes);
+                nodeDisplacementSurfaceOptimizerParallel(procBndPoints);
 
             if( remapVertex )
                 mapper.mapVerticesOntoSurface(movedPoints);
@@ -300,7 +363,7 @@ bool meshSurfaceOptimizer::preOptimizeSurface(const label nLayers)
             Info << "Smoothing remaining inverted vertices " << endl;
 
             movedPoints.clear();
-            procBndNodes.clear();
+            procBndPoints.clear();
             forAll(smoothVertex, bpI)
                 if( smoothVertex[bpI] )
                 {
@@ -308,7 +371,7 @@ bool meshSurfaceOptimizer::preOptimizeSurface(const label nLayers)
 
                     if( vertexType_[bpI] & PROCBND )
                     {
-                        procBndNodes.append(bpI);
+                        procBndPoints.append(bpI);
                         continue;
                     }
 
@@ -317,7 +380,7 @@ bool meshSurfaceOptimizer::preOptimizeSurface(const label nLayers)
 
             if( Pstream::parRun() )
             {
-                nodeDisplacementLaplacianFCParallel(procBndNodes, false);
+                nodeDisplacementLaplacianFCParallel(procBndPoints, false);
             }
 
             if( remapVertex )
@@ -332,11 +395,20 @@ bool meshSurfaceOptimizer::preOptimizeSurface(const label nLayers)
 
     } while( nInvertedTria && (++nGlobalIter < 10) );
 
-    Info << "Finished optimizing positions of surface nodes" << endl;
+    Info << "Finished untangling the surface of the volume mesh" << endl;
 
     return changed;
 }
 
+bool meshSurfaceOptimizer::untangleSurface(const label nAdditionalLayers)
+{
+    labelLongList selectedNodes(surfaceEngine_.boundaryPoints().size());
+    forAll(selectedNodes, i)
+        selectedNodes[i] = i;
+
+    return untangleSurface(selectedNodes, nAdditionalLayers);
+}
+
 void meshSurfaceOptimizer::optimizeSurface(const label nIterations)
 {
     const labelList& bPoints = surfaceEngine_.boundaryPoints();
@@ -468,7 +540,7 @@ void meshSurfaceOptimizer::optimizeSurface(const label nIterations)
 
     Info << endl;
 
-    preOptimizeSurface(0);
+    untangleSurface(0);
 }
 
 void meshSurfaceOptimizer::optimizeSurface2D(const label nIterations)
@@ -528,7 +600,7 @@ void meshSurfaceOptimizer::optimizeSurface2D(const label nIterations)
 
         meshSurfaceEngineModifier bMod(surfaceEngine_);
         # ifdef USE_OMP
-        # pragma omp parallel if( edgePoints.size() > 1000 )
+        # pragma omp parallel if( edgePoints.size() > 100 )
         # endif
         {
             LongList<labelledPoint> newPos;
diff --git a/meshLibrary/utilities/surfaceTools/correctEdgesBetweenPatches/correctEdgesBetweenPatchesDistributeFaces.C b/meshLibrary/utilities/surfaceTools/correctEdgesBetweenPatches/correctEdgesBetweenPatchesDistributeFaces.C
index 612493b42342ed3cda434928e612394f00765f99..24968eb6b831e42641873f932a60063670b84f3d 100644
--- a/meshLibrary/utilities/surfaceTools/correctEdgesBetweenPatches/correctEdgesBetweenPatchesDistributeFaces.C
+++ b/meshLibrary/utilities/surfaceTools/correctEdgesBetweenPatches/correctEdgesBetweenPatchesDistributeFaces.C
@@ -288,8 +288,8 @@ void correctEdgesBetweenPatches::patchCorrection()
         nodeType[it.key()] |= 1;
 
     //- set flgs to edge vertices
-    const labelHashSet& edgeNodes = surfacePartitioner.edgeNodes();
-    forAllConstIter(labelHashSet, edgeNodes, it)
+    const labelHashSet& edgePoints = surfacePartitioner.edgePoints();
+    forAllConstIter(labelHashSet, edgePoints, it)
         nodeType[it.key()] |= 2;
 
     //- set flags for feature edges
diff --git a/meshLibrary/utilities/surfaceTools/edgeExtraction/edgeExtractor/edgeExtractor.C b/meshLibrary/utilities/surfaceTools/edgeExtraction/edgeExtractor/edgeExtractor.C
index d01e5b6fc0efe73ddf483ea7b749bab6a9091579..cf87f1604e42a4572cad5f9a329b057aa2232be2 100644
--- a/meshLibrary/utilities/surfaceTools/edgeExtraction/edgeExtractor/edgeExtractor.C
+++ b/meshLibrary/utilities/surfaceTools/edgeExtraction/edgeExtractor/edgeExtractor.C
@@ -45,6 +45,8 @@ Description
 #include "triSurfacePartitioner.H"
 #include "triSurfaceClassifyEdges.H"
 #include "meshSurfaceMapper.H"
+#include "meshSurfaceCheckInvertedVertices.H"
+#include "meshSurfaceCheckEdgeTypes.H"
 
 # ifdef USE_OMP
 #include <omp.h>
@@ -122,6 +124,8 @@ void edgeExtractor::calculateValence()
 void edgeExtractor::calculateSingleCellEdge()
 {
     const meshSurfaceEngine& mse = this->surfaceEngine();
+    const edgeList& edges = mse.edges();
+    const VRWGraph& bpEdges = mse.boundaryPointEdges();
     const VRWGraph& edgeFaces = mse.edgeFaces();
     const labelList& faceCells = mse.faceOwners();
 
@@ -140,6 +144,55 @@ void edgeExtractor::calculateSingleCellEdge()
                 edgeType_[eI] |= SINGLECELLEDGE;
         }
     }
+
+    //- calculate the number of cells attache to a boundary edge
+    const labelList& bp = mse.bp();
+    const cellListPMG& cells = mse.mesh().cells();
+    const faceListPMG& faces = mse.faces();
+
+    nCellsAtEdge_.setSize(edgeFaces.size());
+    nCellsAtEdge_ = 0;
+
+    # ifdef USE_OMP
+    # pragma omp parallel for schedule(dynamic, 100)
+    # endif
+    forAll(cells, cellI)
+    {
+        const cell& c = cells[cellI];
+
+        DynList<edge> foundEdge;
+
+        forAll(c, fI)
+        {
+            const face& f = faces[c[fI]];
+
+            forAll(f, eI)
+            {
+                const edge e = f.faceEdge(eI);
+
+                const label bps = bp[e.start()];
+
+                if( bps < 0 )
+                    continue;
+
+                forAllRow(bpEdges, bps, i)
+                {
+                    const label beI = bpEdges(bps, i);
+                    const edge& be = edges[beI];
+
+                    if( (e == be) && !foundEdge.contains(be) )
+                    {
+                        foundEdge.append(be);
+
+                        # ifdef USE_OMP
+                        # pragma omp atomic
+                        # endif
+                        ++nCellsAtEdge_[beI];
+                    }
+                }
+            }
+        }
+    }
 }
 
 void edgeExtractor::findPatchesNearSurfaceFace()
@@ -458,7 +511,7 @@ const triSurfaceClassifyEdges& edgeExtractor::edgeClassifier() const
 void edgeExtractor::findFaceCandidates
 (
     labelLongList& faceCandidates,
-    const labelLongList* facePatchPtr,
+    const labelList* facePatchPtr,
     const Map<label>* otherFacePatchPtr
 ) const
 {
@@ -466,7 +519,7 @@ void edgeExtractor::findFaceCandidates
     if( !facePatchPtr )
         facePatchPtr = &facePatch_;
 
-    const labelLongList& fPatches = *facePatchPtr;
+    const labelList& fPatches = *facePatchPtr;
 
     if( !otherFacePatchPtr )
     {
@@ -535,7 +588,7 @@ void edgeExtractor::findFaceCandidates
 void edgeExtractor::findOtherFacePatchesParallel
 (
     Map<label>& otherFacePatch,
-    const labelLongList* facePatchPtr
+    const labelList* facePatchPtr
 ) const
 {
     otherFacePatch.clear();
@@ -543,7 +596,7 @@ void edgeExtractor::findOtherFacePatchesParallel
     if( !facePatchPtr )
         facePatchPtr = &facePatch_;
 
-    const labelLongList& fPatches = *facePatchPtr;
+    const labelList& fPatches = *facePatchPtr;
 
     if( Pstream::parRun() )
     {
@@ -607,6 +660,7 @@ edgeExtractor::edgeExtractor
     pointValence_(),
     pointPatch_(),
     facePatch_(),
+    nCellsAtEdge_(),
     edgeType_(),
     patchesNearFace_(),
     featureEdgesNearEdge_()
@@ -663,11 +717,12 @@ void edgeExtractor::moveVerticesTowardsDiscontinuities(const label nIterations)
 
                 point newP;
                 scalar distSq;
-                label patchI;
+                label patchI, nearestTri;
                 meshOctree_.findNearestSurfacePoint
                 (
                     newP,
                     distSq,
+                    nearestTri,
                     patchI,
                     centre
                 );
@@ -761,10 +816,10 @@ void edgeExtractor::moveVerticesTowardsDiscontinuities(const label nIterations)
             //Info << "Moved point " << mp << endl;
 
             point newPoint;
-            label patchI;
+            label patchI, nt;
             scalar distSq;
 
-            meshOctree_.findNearestSurfacePoint(newPoint, distSq, patchI, mp);
+            meshOctree_.findNearestSurfacePoint(newPoint, distSq, nt, patchI, mp);
 
             //Info << "Mapped point " << newPoint << nl << endl;
 
@@ -784,6 +839,7 @@ void edgeExtractor::distributeBoundaryFaces()
 {
     const meshSurfaceEngine& mse = this->surfaceEngine();
     const labelList& bPoints = mse.boundaryPoints();
+    const labelList& bp = mse.bp();
     const faceList::subList& bFaces = mse.boundaryFaces();
     const pointFieldPMG& points = mse.points();
 
@@ -808,11 +864,11 @@ void edgeExtractor::distributeBoundaryFaces()
     {
         const point& bp = points[bPoints[bpI]];
 
-        label fPatch;
+        label fPatch, nTri;
         point p;
         scalar distSq;
 
-        meshOctree_.findNearestSurfacePoint(p, distSq, fPatch, bp);
+        meshOctree_.findNearestSurfacePoint(p, distSq, nTri, fPatch, bp);
 
         if( (fPatch > -1) && (fPatch < nPatches) )
         {
@@ -841,11 +897,11 @@ void edgeExtractor::distributeBoundaryFaces()
         const point c = bf.centre(points);
 
         //- find the nearest surface patch to face centre
-        label fPatch;
+        label fPatch, nTri;
         point p;
         scalar distSq;
 
-        meshOctree_.findNearestSurfacePoint(p, distSq, fPatch, c);
+        meshOctree_.findNearestSurfacePoint(p, distSq, nTri, fPatch, c);
 
         if( (fPatch > -1) && (fPatch < nPatches) )
         {
@@ -863,6 +919,160 @@ void edgeExtractor::distributeBoundaryFaces()
     }
 }
 
+bool edgeExtractor::distributeBoundaryFacesNormalAlignment()
+{
+    bool changed(false);
+
+    const pointFieldPMG& points = mesh_.points();
+    const meshSurfaceEngine& mse = this->surfaceEngine();
+    const faceList::subList& bFaces = mse.boundaryFaces();
+    const labelList& bp = mse.bp();
+    const VRWGraph& faceEdges = mse.faceEdges();
+    const VRWGraph& edgeFaces = mse.edgeFaces();
+
+    const triSurf& surf = meshOctree_.surface();
+    const pointField& sPoints = surf.points();
+
+    label nCorrected;
+    Map<label> otherProcNewPatch;
+
+    do
+    {
+        nCorrected = 0;
+
+        //- allocate a copy of boundary patches
+        labelList newBoundaryPatches(facePatch_);
+
+        //- check whether there exist situations where a boundary face
+        //- is surrounded by more faces in different patches than the
+        //- faces in the current patch
+        if( Pstream::parRun() )
+        {
+            findOtherFacePatchesParallel
+            (
+                otherProcNewPatch,
+                &facePatch_
+            );
+        }
+
+        //- find the faces which have neighbouring faces in other patches
+        labelLongList candidates;
+        findFaceCandidates(candidates, &facePatch_, &otherProcNewPatch);
+
+        //- go through the list of faces and check if they shall remain
+        //- in the current patch
+        # ifdef USE_OMP
+        # pragma omp parallel for schedule(dynamic, 40) \
+        reduction(+ : nCorrected)
+        # endif
+        forAll(candidates, i)
+        {
+            const label bfI = candidates[i];
+            const face& bf = bFaces[bfI];
+
+            DynList<label> allNeiPatches;
+            DynList<label> neiPatches;
+            neiPatches.setSize(faceEdges.sizeOfRow(bfI));
+
+            forAllRow(faceEdges, bfI, eI)
+            {
+                const label beI = faceEdges(bfI, eI);
+
+                if( edgeFaces.sizeOfRow(beI) == 2 )
+                {
+                    label fNei = edgeFaces(beI, 0);
+                    if( fNei == bfI )
+                        fNei = edgeFaces(faceEdges(bfI, eI), 1);
+
+                    allNeiPatches.appendIfNotIn(facePatch_[fNei]);
+                    neiPatches[eI] = facePatch_[fNei];
+                }
+                else if( edgeFaces.sizeOfRow(beI) == 1 )
+                {
+                    allNeiPatches.appendIfNotIn(otherProcNewPatch[beI]);
+                    neiPatches[eI] = otherProcNewPatch[beI];
+                }
+            }
+
+            //- do not modify faces with all neighbours in the same patch
+            if
+            (
+                (allNeiPatches.size() == 1) &&
+                (allNeiPatches[0] == facePatch_[bfI])
+            )
+                continue;
+
+            //- check whether there exist edges which are more suitable for
+            //- projection onto feature edges than the currently selected ones
+            label newPatch(-1);
+            DynList<scalar> normalAlignment(allNeiPatches.size());
+            DynList<scalar> distanceSq(allNeiPatches.size());
+            scalar maxDSq(0.0);
+            forAll(allNeiPatches, i)
+            {
+                point pMap;
+                scalar dSq(VGREAT);
+                label nearestTriangle;
+
+                point p = bf.centre(points);
+                meshOctree_.findNearestSurfacePointInRegion
+                (
+                    pMap,
+                    dSq,
+                    nearestTriangle,
+                    allNeiPatches[i],
+                    p
+                );
+
+                maxDSq = Foam::max(dSq, maxDSq);
+
+                //- calculate normal vectors
+                vector tn = surf[nearestTriangle].normal(sPoints);
+                tn /= (mag(tn) + VSMALL);
+                vector fn = bf.normal(points);
+                fn /= (mag(fn) + SMALL);
+
+                //- calculate alignment
+                normalAlignment[i] = mag(tn & fn);
+                distanceSq[i] = dSq;
+            }
+
+            scalar maxAlignment(0.0);
+            forAll(normalAlignment, i)
+            {
+                const scalar metric
+                (
+                    sqrt(maxDSq / (distanceSq[i] + VSMALL)) * normalAlignment[i]
+                );
+
+                if( metric > maxAlignment )
+                {
+                    maxAlignment = metric;
+                    newPatch = allNeiPatches[i];
+                }
+            }
+
+            if( (newPatch >= 0) && (newPatch != facePatch_[bfI]) )
+            {
+                newBoundaryPatches[bfI] = newPatch;
+                ++nCorrected;
+            }
+        }
+
+        reduce(nCorrected, sumOp<label>());
+
+        if( nCorrected )
+        {
+            changed = true;
+
+            //- transfer the new patches back
+            facePatch_.transfer(newBoundaryPatches);
+        }
+    } while( nCorrected != 0 );
+
+    return changed;
+}
+
 void edgeExtractor::findEdgeCandidates()
 {
     const triSurf& surface = meshOctree_.surface();
@@ -1077,7 +1287,7 @@ bool edgeExtractor::checkConcaveEdgeCells()
     const List<direction>& edgeType = edgeClassifier.edgeTypes();
 
     //- create a copy of facePatch array for local modifications
-    labelLongList newBoundaryPatches(facePatch_);
+    labelList newBoundaryPatches(facePatch_);
 
     //- start checking the surface of the mesh
     label nChanged;
@@ -1263,27 +1473,32 @@ bool edgeExtractor::checkConcaveEdgeCells()
     return changed;
 }
 
-bool edgeExtractor::checkFacePatches()
+bool edgeExtractor::checkFacePatchesTopology()
 {
     bool changed(false);
 
-    //const pointFieldPMG& points = mesh_.points();
+    const pointFieldPMG& points = mesh_.points();
     const meshSurfaceEngine& mse = this->surfaceEngine();
     const faceList::subList& bFaces = mse.boundaryFaces();
     const labelList& bp = mse.bp();
     const VRWGraph& faceEdges = mse.faceEdges();
     const VRWGraph& edgeFaces = mse.edgeFaces();
 
-    //- allocate a copy of boundary patches
-    labelLongList newBoundaryPatches(facePatch_);
+    const triSurf& surf = meshOctree_.surface();
+    const pointField& sPoints = surf.points();
 
     label nCorrected;
     Map<label> otherProcNewPatch;
 
+    label nIter(0);
     do
     {
+        Info << "Iteration " << nIter++ << endl;
         nCorrected = 0;
 
+        //- allocate a copy of boundary patches
+        labelList newBoundaryPatches(facePatch_);
+
         //- check whether there exist situations where a boundary face
         //- is surrounded by more faces in different patches than the
         //- faces in the current patch
@@ -1292,13 +1507,13 @@ bool edgeExtractor::checkFacePatches()
             findOtherFacePatchesParallel
             (
                 otherProcNewPatch,
-                &newBoundaryPatches
+                &facePatch_
             );
         }
 
         //- find the faces which have neighbouring faces in other patches
         labelLongList candidates;
-        findFaceCandidates(candidates, &newBoundaryPatches, &otherProcNewPatch);
+        findFaceCandidates(candidates, &facePatch_, &otherProcNewPatch);
 
         //- go through the list of faces and check if they shall remain
         //- in the current patch
@@ -1315,7 +1530,7 @@ bool edgeExtractor::checkFacePatches()
             //- onto the same patch
             bool allInSamePatch(true);
             forAll(bf, pI)
-                if( pointPatch_[bp[bf[pI]]] != newBoundaryPatches[bfI] )
+                if( pointPatch_[bp[bf[pI]]] != facePatch_[bfI] )
                 {
                     allInSamePatch = false;
                     break;
@@ -1338,8 +1553,8 @@ bool edgeExtractor::checkFacePatches()
                     if( fNei == bfI )
                         fNei = edgeFaces(faceEdges(bfI, eI), 1);
 
-                    allNeiPatches.appendIfNotIn(newBoundaryPatches[fNei]);
-                    neiPatches[eI] = newBoundaryPatches[fNei];
+                    allNeiPatches.appendIfNotIn(facePatch_[fNei]);
+                    neiPatches[eI] = facePatch_[fNei];
                 }
                 else if( edgeFaces.sizeOfRow(beI) == 1 )
                 {
@@ -1352,10 +1567,168 @@ bool edgeExtractor::checkFacePatches()
             if
             (
                 (allNeiPatches.size() == 1) &&
-                (allNeiPatches[0] == newBoundaryPatches[bfI])
+                (allNeiPatches[0] == facePatch_[bfI])
             )
                 continue;
 
+            //- check whether there exist edges which are more suitable for
+            //- projection onto feature edges than the currently selected ones
+            label newPatch(-1);
+/*            DynList<scalar> smallestDistanceToPatch(allNeiPatches.size(), VGREAT);
+            DynList<scalar> normalAlignment(allNeiPatches.size());
+            DynList<scalar> alignmentToPatch(allNeiPatches.size(), 0.0);
+            DynList<bool> analysed(allNeiPatches.size(), false);
+            forAll(allNeiPatches, i)
+            {
+                point pMap;
+                scalar minDistSq(VGREAT);
+                label nearestTriangle;
+
+                point p = bf.centre(points);
+                meshOctree_.findNearestSurfacePointInRegion
+                (
+                    pMap,
+                    minDistSq,
+                    nearestTriangle,
+                    allNeiPatches[i],
+                    p
+                );
+                vector tn = surf[nearestTriangle].normal(sPoints);
+                tn /= (mag(tn) + VSMALL);
+                vector fn = bf.normal(points);
+                fn /= (mag(fn) + SMALL);
+                normalAlignment[i] = mag(tn & fn);
+
+                if( allNeiPatches[i] == facePatch_[bfI] )
+                    continue;
+
+                DynList<label> ePatches(2);
+                ePatches[0] = facePatch_[bfI];
+                ePatches[1] = allNeiPatches[i];
+
+                //- calculate distance of face points from the nearest
+                //- point on the feature edge between the selected patches
+                DynList<point> nearestOnEdge(bf.size());
+                DynList<scalar> distFromEdge(bf.size());
+
+                label nearest(-1);
+                minDistSq = VGREAT;
+
+                forAll(bf, pI)
+                {
+                    const point& p = points[bf[pI]];
+
+                    meshOctree_.findNearestPointToPatches
+                    (
+                        nearestOnEdge[pI],
+                        distFromEdge[pI],
+                        p,
+                        ePatches
+                    );
+
+                    if( distFromEdge[pI] < minDistSq )
+                    {
+                        minDistSq = distFromEdge[pI];
+                        nearest = pI;
+                    }
+                }
+
+                //- calculate edge alignment
+                DynList<scalar> edgeAlignment(bf.size());
+                forAll(bf, eI)
+                {
+                    const edge e = bf.faceEdge(eI);
+                    vector ev = e.vec(points);
+                    const scalar magEv = mag(ev);
+                    ev /= (magEv + VSMALL);
+
+                    vector pev = nearestOnEdge[bf.fcIndex(eI)] - nearestOnEdge[eI];
+                    pev /= (mag(pev) + VSMALL);
+
+                    edgeAlignment[eI] = (1.0 + (ev & pev)) / 2.0;
+                }
+
+                if( nearest < 0 )
+                    FatalErrorIn
+                    (
+                        "void edgeExtractor::checkFacePatches()"
+                    ) << "Cannot find nearest point in patch "
+                      << ePatches[1] << abort(FatalError);
+
+                Info << "\nBoundary face " << bfI << " has normal " << bf.normal(points) << endl;
+                forAll(bf, pI)
+                    Info << "Face point " << pI << " has coordinates " << points[bf[pI]] << endl;
+                Info << "ePatches " << ePatches << endl;
+                Info << "neiPatches " << neiPatches << endl;
+                Info << "Nearest point position " << nearest << endl;
+                Info << "nearestOnEdge " << nearestOnEdge << endl;
+                Info << "distFromEdge " << distFromEdge << endl;
+
+                //- find the edge which fits best to the feature edge
+                //- between the selected patches
+                const scalar aPrev = edgeAlignment[bf.rcIndex(nearest)];
+                const scalar aNext = edgeAlignment[nearest];
+
+                Info << "Alignment prev " << aPrev << endl;
+                Info << "Alignment next " << aNext << endl;
+
+                const label currPatch = facePatch_[bfI];
+                const label prevPatch = neiPatches[bf.rcIndex(nearest)];
+                const label nextPatch = neiPatches[nearest];
+
+                if( (aPrev >= aNext) && (prevPatch == currPatch) )
+                {
+                    Info << "1. Changing patch from " << currPatch
+                         << " to " << allNeiPatches[i] << endl;
+
+                    smallestDistanceToPatch[i] = minDistSq;
+                    alignmentToPatch[i] = aPrev;
+                    analysed[i] = true;
+                }
+                else if( (aNext > aPrev) && (nextPatch == currPatch) )
+                {
+                    Info << "2. Changing patch from " << currPatch
+                         << " to " << allNeiPatches[i] << endl;
+
+                    smallestDistanceToPatch[i] = minDistSq;
+                    alignmentToPatch[i] = aNext;
+                    analysed[i] = true;
+                }
+            }
+
+            Info << "allNeiPatches " << allNeiPatches << endl;
+            Info << "Analysed " << analysed << endl;
+            Info << "Smallest distance " << smallestDistanceToPatch << endl;
+            Info << "Best edge alignment " << alignmentToPatch << endl;
+            Info << "Nornal alignment " << normalAlignment << endl;
+
+            bool alreadyAnalysed(false);
+            forAll(analysed, aI)
+            {
+                if( analysed[aI] )
+                {
+                    if( alreadyAnalysed )
+                    {
+                        newPatch = -1;
+                        break;
+                    }
+
+                    newPatch = allNeiPatches[aI];
+                    alreadyAnalysed = true;
+                }
+            }
+
+            if( (newPatch >= 0) && (newPatch != facePatch_[bfI]) )
+            {
+                Info << "All patches " << allNeiPatches
+                     << " current patch " << facePatch_[bfI]
+                     << " to patch " << newPatch << endl;
+
+                newBoundaryPatches[bfI] = newPatch;
+                ++nCorrected;
+                continue;
+            }
+*/
             //- check if some faces have to be distributed to another patch
             //- in order to reduce the number of feature edges
             Map<label> nNeiInPatch(allNeiPatches.size());
@@ -1364,7 +1737,8 @@ bool edgeExtractor::checkFacePatches()
             forAll(neiPatches, eI)
                 ++nNeiInPatch[neiPatches[eI]];
 
-            label newPatch(-1), nNeiEdges(0);
+            newPatch = -1;
+            label nNeiEdges(0);
             forAllConstIter(Map<label>, nNeiInPatch, it)
             {
                 if( it() > nNeiEdges )
@@ -1374,16 +1748,24 @@ bool edgeExtractor::checkFacePatches()
                 }
                 else if
                 (
-                    (it() == nNeiEdges) && (it.key() == newBoundaryPatches[bfI])
+                    (it() == nNeiEdges) && (it.key() == facePatch_[bfI])
                 )
                 {
                     newPatch = it.key();
                 }
             }
 
-            if( (newPatch < 0) || (newPatch == newBoundaryPatches[bfI]) )
+            //- do not swap in case the
+            if( (newPatch < 0) || (newPatch == facePatch_[bfI]) )
                 continue;
 
+            //- do not swap if the number of neighbours in the patch
+            //- does not dominate
+            if( nNeiEdges <= (bf.size() / 2) )
+                continue;
+
+            //- check whether the edges shared ith the neighbour patch form
+            //- a singly linked chain
             DynList<bool> sharedEdge;
             sharedEdge.setSize(bFaces[bfI].size());
             sharedEdge = false;
@@ -1403,71 +1785,718 @@ bool edgeExtractor::checkFacePatches()
         reduce(nCorrected, sumOp<label>());
 
         if( nCorrected )
+        {
             changed = true;
 
-    } while( nCorrected != 0 );
+            //- transfer the new patches back
+            facePatch_.transfer(newBoundaryPatches);
 
-    //- transfer the new patches back
-    facePatch_.transfer(newBoundaryPatches);
+            const triSurf* surfPtr = surfaceWithPatches();
+            surfPtr->writeSurface("checkFacePatches_"+help::scalarToText(nIter)+".stl");
+            deleteDemandDrivenData(surfPtr);
+        }
+
+        //break;
+
+    } while( nCorrected != 0 && (nIter < 3) );
 
     return changed;
 }
 
-bool edgeExtractor::findCornerCandidates()
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *//
+
+namespace featureEdgeHelpers
 {
-    bool changed(false);
 
-    const triSurf& surf = meshOctree_.surface();
-    const pointField& sp = surf.points();
+class featureEdgesNeiOp
+{
+    // Private data
+        //- reference to meshSurfaceEngine
+        const meshSurfaceEngine& mse_;
 
-    //- create the surface partitioner and find the corners on the surface mesh
-    triSurfacePartitioner surfPartitioner(surf);
+        //- refence to a list holding information which edges are feature edges
+        const boolList& isFeatureEdge_;
 
-    const labelList& corners = surfPartitioner.corners();
-    const List<DynList<label> >& cornerPatches =
-        surfPartitioner.cornerPatches();
+        //- number of feature edges at a surface point
+        labelList nFeatureEdgesAtPoint_;
 
-    Map<label> surfacePointToCorner;
-    forAll(corners, cornerI)
-        surfacePointToCorner.insert(corners[cornerI], cornerI);
+    // Private member functions
+        //- calculate the number of feature edges connected to a surface vertex
+        void calculateNumberOfEdgesAtPoint()
+        {
+            const labelList& bp = mse_.bp();
+            const edgeList& edges = mse_.edges();
 
-    List<labelledScalar> nearestPoint
-    (
-        corners.size(),
-        labelledScalar(0, VGREAT)
-    );
+            nFeatureEdgesAtPoint_.setSize(mse_.boundaryPoints().size());
+            nFeatureEdgesAtPoint_ = 0;
 
-    //- calculate the search range
-    const meshSurfaceEngine& mse = this->surfaceEngine();
-    const pointFieldPMG& points = mesh_.points();
-    const labelList& bPoints = mse.boundaryPoints();
-    const labelList& bp = mse.bp();
-    const faceList::subList& bFaces = mse.boundaryFaces();
+            forAll(isFeatureEdge_, edgeI)
+            {
+                if( !isFeatureEdge_[edgeI] )
+                    continue;
 
-    scalarList searchRange(bPoints.size(), 0.0);
-    forAll(bFaces, bfI)
-    {
-        const face& bf = bFaces[bfI];
-        const point c = bf.centre(points);
+                const edge& e = edges[edgeI];
+                ++nFeatureEdgesAtPoint_[bp[e.start()]];
+                ++nFeatureEdgesAtPoint_[bp[e.end()]];
+            }
 
-        forAll(bf, pI)
-        {
-            const scalar d = 2.0 * Foam::mag(c - points[bf[pI]]);
-            const label bpI = bp[bf[pI]];
+            if( Pstream::parRun() )
+            {
+                const Map<label>& globalToLocal =
+                    mse_.globalToLocalBndPointAddressing();
+                const DynList<label>& neiProcs = mse_.bpNeiProcs();
+                const VRWGraph& bpAtProcs = mse_.bpAtProcs();
 
-            searchRange[bpI] = Foam::max(d, searchRange[bpI]);
-        }
-    }
+                std::map<label, DynList<labelPair> > exchangeData;
+                forAll(neiProcs, i)
+                    exchangeData[neiProcs[i]].clear();
 
-    if( Pstream::parRun() )
-    {
-        const VRWGraph& bpAtProcs = mse.beAtProcs();
-        const DynList<label>& bpNeiProcs = mse.bpNeiProcs();
-        const Map<label>& globalToLocal =
-            mse.globalToLocalBndPointAddressing();
+                //- fill the data from sending
+                forAllConstIter(Map<label>, globalToLocal, it)
+                {
+                    const label bpI = it();
 
-        std::map<label, LongList<labelledScalar> > exchangeData;
-        forAll(bpNeiProcs, i)
+                    forAllRow(bpAtProcs, bpI, i)
+                    {
+                        const label neiProc = bpAtProcs(bpI, i);
+
+                        if( neiProc == Pstream::myProcNo() )
+                            continue;
+
+                        exchangeData[neiProc].append
+                        (
+                            labelPair(it.key(), nFeatureEdgesAtPoint_[bpI])
+                        );
+                    }
+                }
+
+                //- exchange the data between the procesors
+                LongList<labelPair> receivedData;
+                help::exchangeMap(exchangeData, receivedData);
+
+                forAll(receivedData, i)
+                {
+                    const labelPair& lp = receivedData[i];
+
+                    nFeatureEdgesAtPoint_[globalToLocal[lp.first()]] +=
+                        lp.second();
+                }
+            }
+        }
+
+public:
+
+    featureEdgesNeiOp
+    (
+        const meshSurfaceEngine& mse,
+        const boolList& isFeatureEdge
+    )
+    :
+        mse_(mse),
+        isFeatureEdge_(isFeatureEdge),
+        nFeatureEdgesAtPoint_()
+    {
+        calculateNumberOfEdgesAtPoint();
+    }
+
+    label size() const
+    {
+        return isFeatureEdge_.size();
+    }
+
+    void operator()(const label beI, DynList<label>& neighbourEdges) const
+    {
+        neighbourEdges.clear();
+
+        const VRWGraph& bpEdges = mse_.boundaryPointEdges();
+        const labelList& bp = mse_.bp();
+        const edgeList& edges = mse_.edges();
+
+        const edge& e = edges[beI];
+
+        const label bps = bp[e.start()];
+        const label bpe = bp[e.end()];
+
+        if( nFeatureEdgesAtPoint_[bps] == 2 )
+        {
+            forAllRow(bpEdges, bps, peI)
+            {
+                const label beJ = bpEdges(bps, peI);
+
+                if( (beJ == beI) || !isFeatureEdge_[beJ] )
+                    continue;
+
+                neighbourEdges.append(beJ);
+            }
+        }
+
+        if( nFeatureEdgesAtPoint_[bpe] == 2 )
+        {
+            forAllRow(bpEdges, bpe, peI)
+            {
+                const label beJ = bpEdges(bpe, peI);
+
+                if( (beJ == beI) || !isFeatureEdge_[beJ] )
+                    continue;
+
+                neighbourEdges.append(beJ);
+            }
+        }
+    }
+
+    template<class labelListType>
+    void collectGroups
+    (
+        std::map<label, DynList<label> >& neiGroups,
+        const labelListType& elementInGroup,
+        const DynList<label>& localGroupLabel
+    ) const
+    {
+        const Map<label>& globalToLocal = mse_.globalToLocalBndPointAddressing();
+        const VRWGraph& bpAtProcs = mse_.bpAtProcs();
+        const VRWGraph& bpEdges = mse_.boundaryPointEdges();
+
+        const DynList<label>& neiProcs = mse_.beNeiProcs();
+
+        std::map<label, DynList<labelPair> > exchangeData;
+        forAll(neiProcs, i)
+            exchangeData[neiProcs[i]].clear();
+
+        forAllConstIter(Map<label>, globalToLocal, it)
+        {
+            const label bpI = it();
+
+            if( nFeatureEdgesAtPoint_[bpI] != 2 )
+                continue;
+
+            forAllRow(bpEdges, bpI, i)
+            {
+                const label beI = bpEdges(bpI, i);
+
+                if( !isFeatureEdge_[beI] )
+                    continue;
+
+                const label groupI = elementInGroup[beI];
+
+                forAllRow(bpAtProcs, bpI, ppI)
+                {
+                    const label neiProc = bpAtProcs(bpI, ppI);
+
+                    if( neiProc == Pstream::myProcNo() )
+                        continue;
+
+                    exchangeData[neiProc].append
+                    (
+                        labelPair(it.key(), localGroupLabel[groupI])
+                    );
+                }
+            }
+        }
+
+        LongList<labelPair> receivedData;
+        help::exchangeMap(exchangeData, receivedData);
+
+        forAll(receivedData, i)
+        {
+            const labelPair& lp = receivedData[i];
+            const label groupI = elementInGroup[globalToLocal[lp.first()]];
+
+            DynList<label>& ng = neiGroups[localGroupLabel[groupI]];
+
+            //- store the connection over the inter-processor boundary
+            ng.appendIfNotIn(lp.second());
+        }
+    }
+};
+
+class featureEdgesSelOp
+{
+    // Private data
+        //- reference to a list holding information which edge is afeture edge
+        const boolList& isFeatureEdge_;
+
+public:
+
+    featureEdgesSelOp(const boolList& isFeatureEdge)
+    :
+        isFeatureEdge_(isFeatureEdge)
+    {}
+
+    bool operator()(const label beI) const
+    {
+        return isFeatureEdge_[beI];
+    }
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *//
+
+} // End namespace featureEdgeHelpers
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+bool edgeExtractor::checkFacePatchesGeometry()
+{
+    bool changed(false);
+
+    const pointFieldPMG& points = mesh_.points();
+    const meshSurfaceEngine& mse = this->surfaceEngine();
+    const labelList& bPoints = mse.boundaryPoints();
+    const faceList::subList& bFaces = mse.boundaryFaces();
+    const labelList& bp = mse.bp();
+    const edgeList& edges = mse.edges();
+    const VRWGraph& faceEdges = mse.faceEdges();
+    const VRWGraph& edgeFaces = mse.edgeFaces();
+
+    //- check if there exist any inverted faces
+    meshSurfaceCheckInvertedVertices surfCheck(mse);
+    const labelHashSet& invertedPoints = surfCheck.invertedVertices();
+
+    if( invertedPoints.size() == 0 )
+        return false;
+
+    WarningIn
+    (
+        "void edgeExtractor::extractEdges()"
+    ) << "Found " << invertedPoints.size()
+      << " points with inverted surface normals. Getting rid of them..."
+      << endl;
+
+    //- untangle the surface
+    labelLongList smoothNodes;
+    forAllConstIter(labelHashSet, invertedPoints, it)
+        smoothNodes.append(bp[it.key()]);
+
+    meshSurfaceOptimizer mso(*surfaceEnginePtr_, meshOctree_);
+    mso.untangleSurface(smoothNodes, 1);
+
+    //- allocate a copy of boundary patches
+    labelList newBoundaryPatches(facePatch_.size());
+
+    meshSurfaceCheckEdgeTypes edgeChecker(mse);
+    const List<direction>& edgeType = edgeChecker.edgeTypes();
+
+    label nCorrected;
+    Map<label> otherProcNewPatch;
+    boolList isFeatureEdge(edgeType.size());
+
+    label iter(0);
+
+    do
+    {
+        nCorrected = 0;
+        newBoundaryPatches = facePatch_;
+
+        //- check whether there exist situations where a boundary face
+        //- is surrounded by more faces in different patches than the
+        //- faces in the current patch
+        if( Pstream::parRun() )
+        {
+            findOtherFacePatchesParallel
+            (
+                otherProcNewPatch,
+                &facePatch_
+            );
+        }
+
+        //- update the information which edges are assigned as feature edges
+        # ifdef USE_OMP
+        # pragma omp parallel for schedule(dynamic, 50)
+        # endif
+        forAll(edgeFaces, beI)
+        {
+            if( edgeFaces.sizeOfRow(beI) == 2 )
+            {
+                const label patch0 = facePatch_[edgeFaces(beI, 0)];
+                const label patch1 = facePatch_[edgeFaces(beI, 1)];
+
+                if( patch0 != patch1 )
+                {
+                    isFeatureEdge[beI] = true;
+                }
+                else
+                {
+                    isFeatureEdge[beI] = false;
+                }
+            }
+            else if( otherProcNewPatch.found(beI) )
+            {
+                const label patch0 = facePatch_[edgeFaces(beI, 0)];
+                const label patch1 = otherProcNewPatch[beI];
+
+                if( patch0 != patch1 )
+                {
+                    isFeatureEdge[beI] = true;
+                }
+                else
+                {
+                    isFeatureEdge[beI] = false;
+                }
+            }
+        }
+
+        //- find patches at each surface point
+        LongList<DynList<label, 4> > pointPatches(bPoints.size());
+        # ifdef USE_OMP
+        # pragma omp parallel for schedule(dynamic, 40)
+        # endif
+        forAll(bFaces, bfI)
+        {
+            const face& bf = bFaces[bfI];
+
+            forAll(bf, pI)
+            {
+                DynList<label, 4>& pp = pointPatches[bp[bf[pI]]];
+
+                # ifdef USE_OMP
+                # pragma omp critical
+                # endif
+                {
+                    pp.appendIfNotIn(facePatch_[bfI]);
+                }
+            }
+        }
+
+        //- find a network of edges marked to become feature edges
+        //- and check which parts of the mesh are concave/convex
+        Info << "Finding Groups of feature edges " << endl;
+        labelList featureEdgeGroup(edgeFaces.size(), -1);
+        label nEdgeGroups =
+            help::groupMarking
+            (
+                featureEdgeGroup,
+                featureEdgeHelpers::featureEdgesNeiOp(mse, isFeatureEdge),
+                featureEdgeHelpers::featureEdgesSelOp(isFeatureEdge)
+            );
+
+        Info << "Found " << nEdgeGroups << " edge groups" << endl;
+
+        //- check whether a group of edges is predominantly convex or concave
+        List<labelPair> groupStatistics(nEdgeGroups, labelPair(0, 0));
+        List<DynList<label, 4> > pointEdgeGroup(pointPatches.size());
+
+        forAll(edgeType, edgeI)
+        {
+            if( !isFeatureEdge[edgeI] )
+                continue;
+
+            const edge& e = edges[edgeI];
+            pointEdgeGroup[bp[e.start()]].appendIfNotIn(featureEdgeGroup[edgeI]);
+            pointEdgeGroup[bp[e.end()]].appendIfNotIn(featureEdgeGroup[edgeI]);
+
+            if( edgeType[edgeI] & meshSurfaceCheckEdgeTypes::CONVEXEDGE )
+            {
+                ++groupStatistics[featureEdgeGroup[edgeI]].first();
+            }
+            else if( edgeType[edgeI] & meshSurfaceCheckEdgeTypes::CONCAVEEDGE )
+            {
+                ++groupStatistics[featureEdgeGroup[edgeI]].second();
+            }
+        }
+
+        List<direction> groupType
+        (
+            nEdgeGroups,
+            meshSurfaceCheckEdgeTypes::UNDETERMINED
+        );
+
+        forAll(groupStatistics, groupI)
+        {
+            const labelPair& gs = groupStatistics[groupI];
+
+            if( gs.first() > 2 * gs.second() )
+            {
+                groupType[groupI] = meshSurfaceCheckEdgeTypes::CONVEXEDGE;
+            }
+            else if( gs.second() > 2 * gs.first() )
+            {
+                groupType[groupI] = meshSurfaceCheckEdgeTypes::CONCAVEEDGE;
+            }
+        }
+
+        forAll(groupType, i)
+            Info << "Group " << i << " is of type " << label(groupType[i])
+                << " num convex " << groupStatistics[i].first()
+                << " num concave " << groupStatistics[i].second() << endl;
+
+        //- find the faces which have neighbouring faces in other patches
+        labelLongList candidates;
+        forAll(bFaces, bfI)
+        {
+            const face& bf = bFaces[bfI];
+
+            forAll(bf, pI)
+                if( invertedPoints.found(bf[pI]) )
+                {
+                    candidates.append(bfI);
+                    break;
+                }
+        }
+
+        Info << "Number of face candidates " << candidates.size() << endl;
+
+        //- check the inverted edge/corner points whether they can be swapped
+        //- with some other point in the vicinity which fits better
+        forAll(candidates, i)
+        {
+            const label bfI = candidates[i];
+            const face& bf = bFaces[bfI];
+
+            //- check which edges of the face are currently marked
+            //- as feature edges, check if there exists corner points
+            DynList<label> allNeiPatches;
+            allNeiPatches.append(facePatch_[bfI]);
+            DynList<label> neiInPatch(bf.size(), facePatch_[bfI]);
+            DynList<label> faceEdgeType(bf.size(), direction(0));
+
+            forAll(bf, pI)
+            {
+                const label beI = faceEdges(bfI, pI);
+
+                if( isFeatureEdge[beI] )
+                {
+                    //- check the neighbour patch over the feature edge
+                    if( edgeFaces.sizeOfRow(beI) == 2 )
+                    {
+                        label neiPatchI = facePatch_[edgeFaces(beI, 0)];
+
+                        if( neiPatchI == neiInPatch[pI] )
+                            neiPatchI = facePatch_[edgeFaces(beI, 1)];
+
+                        neiInPatch[pI] = neiPatchI;
+                        allNeiPatches.appendIfNotIn(neiPatchI);
+                    }
+                    else if( edgeFaces.sizeOfRow(beI) == 1 )
+                    {
+                        neiInPatch[pI] = otherProcNewPatch[beI];
+                        allNeiPatches.appendIfNotIn(neiInPatch[pI]);
+                    }
+                    else
+                    {
+                        FatalErrorIn
+                        (
+                            "bool edgeExtractor::checkFacePatchesGeometry()"
+                        ) << "Invalid surface of the volume mesh at boundary"
+                          << " edge " << beI << abort(FatalError);
+                    }
+
+                    //- mark feature edges in the face
+                    faceEdgeType[pI] = groupType[featureEdgeGroup[beI]];
+                }
+            }
+
+            //- check if there exist better candidates for feature edges
+            scalar minDist(VGREAT);
+            label bestPatch(-1);
+            DynList<scalar> distanceToPatch(allNeiPatches.size());
+            forAll(distanceToPatch, i)
+            {
+                const label searchPatch = allNeiPatches[i];
+
+                distanceToPatch[i] = 0.0;
+
+                forAll(bf, pI)
+                {
+                    const point& p = points[bf[pI]];
+                    point np;
+                    scalar dSq;
+                    label nt;
+
+                    meshOctree_.findNearestSurfacePointInRegion
+                    (
+                        np,
+                        dSq,
+                        nt,
+                        searchPatch,
+                        p
+                    );
+
+                    distanceToPatch[i] += dSq;
+                }
+
+                if( distanceToPatch[i] < minDist )
+                {
+                    minDist = distanceToPatch[i];
+                    bestPatch = allNeiPatches[i];
+                }
+            }
+
+            Info << "allNeiPatches " << allNeiPatches << endl;
+            Info << "Distance to patch " << distanceToPatch << endl;
+            Info << "minDist " << minDist << endl;
+            Info << "best patch " << bestPatch << endl;
+
+            if( bestPatch != facePatch_[bfI] )
+            {
+                newBoundaryPatches[bfI] = bestPatch;
+                ++nCorrected;
+            }
+
+            forAll(faceEdgeType, eI)
+            {
+                const direction eType = faceEdgeType[eI];
+
+                if( !eType )
+                    continue;
+
+
+/*
+                const label beI = faceEdges(bfI, eI);
+                const label pbeI = faceEdges(bfI, bf.rcIndex(eI));
+                const label nbeI = faceEdges(bfI, bf.fcIndex(eI));
+
+                if( eType & meshSurfaceCheckEdgeTypes::CONVEXEDGE )
+                {
+                    //- check if there exist any face edge which less
+                    //- cells attached to it than to the current edge
+                    if
+                    (
+                         !isFeatureEdge[pbeI] &&
+                         (nCellsAtEdge_[pbeI] < nCellsAtEdge_[beI])
+                    )
+                    {
+
+                        //# ifdef DEBUGEdges
+                        Info << "2.1.1 Face " << bf
+                             << " moved from patch " << facePatch_[bfI]
+                             << " to patch  " << neiInPatch[eI] << endl;
+                        //# endif
+
+                        ++nCorrected;
+                        newBoundaryPatches[bfI] = neiInPatch[eI];
+                        break;
+                    }
+                    else if
+                    (
+                        !isFeatureEdge[nbeI] &&
+                        (nCellsAtEdge_[nbeI] < nCellsAtEdge_[beI])
+                    )
+                    {
+                        //# ifdef DEBUGEdges
+                        Info << "2.1.2 Face " << bf
+                             << " moved from patch " << facePatch_[bfI]
+                             << " to patch  " << neiInPatch[eI] << endl;
+                        //# endif
+
+                        ++nCorrected;
+                        newBoundaryPatches[bfI] = neiInPatch[eI];
+                        break;
+                    }
+                }
+                else if( eType & meshSurfaceCheckEdgeTypes::CONCAVEEDGE )
+                {
+                    //- check if there exist any edge in the face
+                    //- with more cells attached to it than to the
+                    //- current feature edge
+                    if
+                    (
+                        !isFeatureEdge[pbeI] &&
+                        (nCellsAtEdge_[pbeI] > nCellsAtEdge_[beI])
+                    )
+                    {
+                        //# ifdef DEBUGEdges
+                        Info << "2.2.1 Face " << bf
+                             << " moved from patch " << facePatch_[bfI]
+                             << " to patch  " << neiInPatch[eI] << endl;
+                        //# endif
+
+                        ++nCorrected;
+                        newBoundaryPatches[bfI] = neiInPatch[eI];
+                        break;
+                    }
+                    else if
+                    (
+                        !isFeatureEdge[nbeI] &&
+                        (nCellsAtEdge_[nbeI] > nCellsAtEdge_[beI])
+                    )
+                    {
+                        //# ifdef DEBUGEdges
+                        Info << "2.2.2 Face " << bf
+                             << " moved from patch " << facePatch_[bfI]
+                             << " to patch  " << neiInPatch[eI] << endl;
+                        //# endif
+
+                        ++nCorrected;
+                        newBoundaryPatches[bfI] = neiInPatch[eI];
+                        break;
+                    }
+                }
+                */
+            }
+        }
+
+        reduce(nCorrected, sumOp<label>());
+
+        //::exit(0);
+
+        if( nCorrected )
+        {
+            changed = true;
+            facePatch_ = newBoundaryPatches;
+        }
+
+        if( ++iter > 0 )
+            break;
+
+    } while( nCorrected != 0 );
+
+    return changed;
+}
+/*
+bool edgeExtractor::findCornerCandidates()
+{
+    bool changed(false);
+
+    const triSurf& surf = meshOctree_.surface();
+    const pointField& sp = surf.points();
+
+    //- create the surface partitioner and find the corners on the surface mesh
+    triSurfacePartitioner surfPartitioner(surf);
+
+    const labelList& corners = surfPartitioner.corners();
+    const List<DynList<label> >& cornerPatches =
+        surfPartitioner.cornerPatches();
+
+    Map<label> surfacePointToCorner;
+    forAll(corners, cornerI)
+        surfacePointToCorner.insert(corners[cornerI], cornerI);
+
+    List<labelledScalar> nearestPoint
+    (
+        corners.size(),
+        labelledScalar(0, VGREAT)
+    );
+
+    //- calculate the search range
+    const meshSurfaceEngine& mse = this->surfaceEngine();
+    const pointFieldPMG& points = mesh_.points();
+    const labelList& bPoints = mse.boundaryPoints();
+    const labelList& bp = mse.bp();
+    const faceList::subList& bFaces = mse.boundaryFaces();
+
+    scalarList searchRange(bPoints.size(), 0.0);
+    forAll(bFaces, bfI)
+    {
+        const face& bf = bFaces[bfI];
+        const point c = bf.centre(points);
+
+        forAll(bf, pI)
+        {
+            const scalar d = 2.0 * Foam::mag(c - points[bf[pI]]);
+            const label bpI = bp[bf[pI]];
+
+            searchRange[bpI] = Foam::max(d, searchRange[bpI]);
+        }
+    }
+
+    if( Pstream::parRun() )
+    {
+        const VRWGraph& bpAtProcs = mse.beAtProcs();
+        const DynList<label>& bpNeiProcs = mse.bpNeiProcs();
+        const Map<label>& globalToLocal =
+            mse.globalToLocalBndPointAddressing();
+
+        std::map<label, LongList<labelledScalar> > exchangeData;
+        forAll(bpNeiProcs, i)
             exchangeData.insert
             (
                 std::make_pair(bpNeiProcs[i], LongList<labelledScalar>())
@@ -1548,7 +2577,7 @@ bool edgeExtractor::findCornerCandidates()
 
     return changed;
 }
-
+*/
 void edgeExtractor::projectDeterminedFeatureVertices()
 {
     List<DynList<label, 5> > pointPatches;
@@ -1559,6 +2588,7 @@ void edgeExtractor::projectDeterminedFeatureVertices()
     const labelList& bPoints = mse.boundaryPoints();
     const labelList& bp = mse.bp();
     const faceList::subList& bFaces = mse.boundaryFaces();
+    meshOctree_.surface().pointEdges();
 
     //- calculate patches for each point
     forAll(bFaces, bfI)
@@ -1583,7 +2613,7 @@ void edgeExtractor::projectDeterminedFeatureVertices()
                 std::make_pair(bpNeiProcs[i], LongList<labelPair>())
             );
 
-        //- collet the data distributed to others
+        //- collect the data distributed to others
         forAllConstIter(Map<label>, globalToLocal, it)
         {
             const label bpI = it();
@@ -1619,46 +2649,74 @@ void edgeExtractor::projectDeterminedFeatureVertices()
 
     meshSurfaceEngineModifier surfMod(mse);
 
-    for(label iterI=0;iterI<2;++iterI)
+    # ifdef USE_OMP
+    # pragma omp parallel for schedule(dynamic, 10)
+    # endif
+    forAll(pointPatches, bpI)
     {
-        # ifdef USE_OMP
-        # pragma omp parallel for schedule(dynamic, 10)
-        # endif
-        forAll(pointPatches, bpI)
-        {
-            const DynList<label, 5>& pPatches = pointPatches[bpI];
+        if( pointPatches[bpI].size() < 2 )
+            continue;
 
-            if( pPatches.size() < 2 )
-                continue;
+        const DynList<label> pPatches = pointPatches[bpI];
 
-            const point& p = points[bPoints[bpI]];
+        const point& p = points[bPoints[bpI]];
+
+        //- find the nearest object on the surface mesh
+        point newP;
+        scalar dSqExact;
+        if( pPatches.size() == 2 )
+        {
+            label nse;
+            meshOctree_.findNearestEdgePoint(newP, dSqExact, nse, p, pPatches);
+        }
+        else
+        {
+            label nsp;
+            meshOctree_.findNearestCorner(newP, dSqExact, nsp, p, pPatches);
+        }
 
-            point newP(vector::zero);
-            label counter(0);
+        //- find the nearest object in an iterative procedure
+        point pp(p);
+        for(label iterI=0;iterI<20;++iterI)
+        {
+            point inp(vector::zero);
 
             forAll(pPatches, i)
             {
-                point pp;
+                point np;
                 scalar dSq;
+                label nt;
+
                 meshOctree_.findNearestSurfacePointInRegion
                 (
-                    pp,
+                    np,
                     dSq,
+                    nt,
                     pPatches[i],
-                    p
+                    pp
                 );
 
-                newP += pp;
-                ++counter;
+                inp += np;
             }
 
-            newP /= counter;
-            surfMod.moveBoundaryVertexNoUpdate(bpI, newP);
+            inp /= pPatches.size();
+            const scalar currDSq = magSqr(inp - pp);
+            pp = inp;
+
+            if( currDSq < 1e-2 * dSqExact )
+                break;
         }
 
-        surfMod.syncVerticesAtParallelBoundaries();
-        surfMod.updateGeometry();
+        //- check if the exact position of the corner is further away
+        //- than the iteratively found object
+        if( dSqExact > 1.1 * magSqr(pp - p) )
+            newP = pp;
+
+        surfMod.moveBoundaryVertexNoUpdate(bpI, newP);
     }
+
+    surfMod.syncVerticesAtParallelBoundaries();
+    surfMod.updateGeometry();
 }
 
 bool edgeExtractor::untangleSurface()
@@ -1668,18 +2726,21 @@ bool edgeExtractor::untangleSurface()
     meshSurfaceEngine& mse =
         const_cast<meshSurfaceEngine&>(this->surfaceEngine());
     meshSurfaceOptimizer optimizer(mse, meshOctree_);
-    changed = optimizer.preOptimizeSurface();
+    changed = optimizer.untangleSurface();
 
     return changed;
 }
 
 void edgeExtractor::extractEdges()
 {
-    bool changed;
     label nIter(0);
 
+    //const meshSurfaceEngine& mse = surfaceEngine();
+
     distributeBoundaryFaces();
 
+    distributeBoundaryFacesNormalAlignment();
+
     //moveVerticesTowardsDiscontinuities(20);
 
 //    updateMeshPatches();
@@ -1688,54 +2749,49 @@ void edgeExtractor::extractEdges()
 //    returnReduce(1, sumOp<label>());
 //    ::exit(0);
 
- //   return;
-
     # ifdef DEBUGEdgeExtractor
     const triSurf* sPtr = surfaceWithPatches();
     sPtr->writeSurface("initialDistributionOfPatches.stl");
     deleteDemandDrivenData(sPtr);
     # endif
 
-    do
+    Info << "Starting topological adjustment of patches" << endl;
+    if( checkFacePatchesTopology() )
     {
-        changed = false;
-
-        Info << "Checking patch in the neighbourhood of each face" << endl;
-        if( checkFacePatches() )
-        {
-            # ifdef DEBUGEdgeExtractor
-            Info << "Changes due to face patches" << endl;
-            fileName sName("checkFacePatches"+help::scalarToText(nIter)+".stl");
-            sPtr = surfaceWithPatches();
-            sPtr->writeSurface(sName);
-            deleteDemandDrivenData(sPtr);
-            # endif
-
-            changed = true;
-        }
+        Info << "Finished topological adjustment of pathes" << endl;
+
+        # ifdef DEBUGEdgeExtractor
+        Info << "Changes due to face patches" << endl;
+        fileName sName("checkFacePatches"+help::scalarToText(nIter)+".stl");
+        sPtr = surfaceWithPatches();
+        sPtr->writeSurface(sName);
+        deleteDemandDrivenData(sPtr);
+        # endif
+    }
+    else
+    {
+        Info << "No topological adjustment was needed" << endl;
+    }
 
-        //findEdgeCandidates();
+    return;
 
-//        if( findCornerCandidates() )
-//        {
-//            # ifdef DEBUGEdgeExtractor
-//            Info << "Changes due to corner candidates" << endl;
-//            fileName sName("findCornerCandidates"+help::scalarToText(nIter)+".stl");
-//            sPtr = surfaceWithPatches();
-//            sPtr->writeSurface(sName);
-//            deleteDemandDrivenData(sPtr);
-//            # endif
+    //- project the selected feature vertices onto the corresponding
+    //- objects on the surface mesh and improve the distribution of patches
+    //- in case the mesh gets tangled
+    bool changed;
 
-//            changed = true;
-//        }
+    do
+    {
+        changed = false;
 
+        //- project feature vertices onto the surface mesh
         projectDeterminedFeatureVertices();
 
-        if( untangleSurface() )
+        if( checkFacePatchesGeometry()  )
         {
             # ifdef DEBUGEdgeExtractor
             Info << "Changes due to untangling" << endl;
-            fileName sName("untangleSurface"+help::scalarToText(nIter)+".stl");
+            fileName sName("untangledSurface"+help::scalarToText(nIter)+".stl");
             sPtr = surfaceWithPatches();
             sPtr->writeSurface(sName);
             deleteDemandDrivenData(sPtr);
@@ -1743,10 +2799,6 @@ void edgeExtractor::extractEdges()
 
             changed = true;
         }
-        else
-        {
-            break;
-        }
 
     } while( changed && ++nIter < 3 );
 
@@ -1816,6 +2868,7 @@ void edgeExtractor::updateMeshPatches()
     wordList patchNames(nPatches);
     VRWGraph newBoundaryFaces;
     labelLongList newBoundaryOwners(bFaces.size());
+    labelLongList newBoundaryPatches(bFaces.size());
 
     //- set patchNames
     forAll(surface.patches(), patchI)
@@ -1826,6 +2879,7 @@ void edgeExtractor::updateMeshPatches()
     {
         newBoundaryFaces.appendList(bFaces[bfI]);
         newBoundaryOwners[bfI] = faceOwner[bfI];
+        newBoundaryPatches[bfI] = facePatch_[bfI];
     }
 
     //- replace the boundary with the new patches
@@ -1834,7 +2888,7 @@ void edgeExtractor::updateMeshPatches()
         patchNames,
         newBoundaryFaces,
         newBoundaryOwners,
-        facePatch_
+        newBoundaryPatches
     );
 }
 
diff --git a/meshLibrary/utilities/surfaceTools/edgeExtraction/edgeExtractor/edgeExtractor.H b/meshLibrary/utilities/surfaceTools/edgeExtraction/edgeExtractor/edgeExtractor.H
index dce70449c2b4e0af6a4470fb222e56f24966edd0..0edf98ccca87c9266d0fd888c4f92445244daa5c 100644
--- a/meshLibrary/utilities/surfaceTools/edgeExtraction/edgeExtractor/edgeExtractor.H
+++ b/meshLibrary/utilities/surfaceTools/edgeExtraction/edgeExtractor/edgeExtractor.H
@@ -67,7 +67,7 @@ class edgeExtractor
         polyMeshGen& mesh_;
 
         //- surface engine
-        mutable const meshSurfaceEngine* surfaceEnginePtr_;
+        mutable meshSurfaceEngine* surfaceEnginePtr_;
 
         //- reference to the octree
         const meshOctree& meshOctree_;
@@ -85,7 +85,10 @@ class edgeExtractor
         labelLongList pointPatch_;
 
         //- boundary face patch
-        labelLongList facePatch_;
+        labelList facePatch_;
+
+        //- number of cells attached to a boundary edge
+        labelLongList nCellsAtEdge_;
 
         //- edge classification
         LongList<direction> edgeType_;
@@ -131,7 +134,7 @@ class edgeExtractor
         void findFaceCandidates
         (
             labelLongList& faceCandidates,
-            const labelLongList* facePatchPtr = NULL,
+            const labelList* facePatchPtr = NULL,
             const Map<label>* otherFacePatchPtr = NULL
         ) const;
 
@@ -139,13 +142,16 @@ class edgeExtractor
         void findOtherFacePatchesParallel
         (
             Map<label>& otherFacePatch,
-            const labelLongList* facePatchPtr = NULL
+            const labelList* facePatchPtr = NULL
         ) const;
 
         //- project face centres on the nearest location at the surface mesh
         //- and assign the patch to the patch of the surface element
         void distributeBoundaryFaces();
 
+        //- move faces into the patch with the best alignment
+        bool distributeBoundaryFacesNormalAlignment();
+
         //- go through the boundary faces which have at least one neighour
         //- assigned to a different patch and check which of its edges
         //- are best candidates to be used as feature edges
@@ -193,16 +199,23 @@ public:
         //- input geometry
         void moveVerticesTowardsDiscontinuities(const label nIterations = 2);
 
-
-
         //- check if there exist cells at concave feature edges which have more
         //- than one face at the boundary and the faces are
         //- distributed into patches at the concave edge
         bool checkConcaveEdgeCells();
 
         //- check and improve the distribution of mesh faces into patches
-        //- in order to minimize the number of decomposed faces and to optimize
-        bool checkFacePatches();
+        //- in order to minimize the number of decomposed faces
+        bool checkFacePatchesTopology();
+
+        //- checks whether there exist corners which do not exist in the surface
+        //- mesh, and checks whether the locations of corners in the volume mesh
+        //- are near the existing counterparts in the surface mesh
+        bool checkCorners();
+
+        //- optimise distribution of mesh faces into patches
+        //- in order to get better geometric quality of the mesh
+        bool checkFacePatchesGeometry();
 
         //- find the nearest points on the surface of the volume mesh
         //- to the corners on the surface mesh
diff --git a/meshLibrary/utilities/surfaceTools/edgeExtraction/edgeExtractor/edgeExtractorCorners.C b/meshLibrary/utilities/surfaceTools/edgeExtraction/edgeExtractor/edgeExtractorCorners.C
index c10ac62047b12c24361e7d560a4f55297b28e3a1..7c99e5292cea14b103912f8d59353055b7ec7178 100644
--- a/meshLibrary/utilities/surfaceTools/edgeExtraction/edgeExtractor/edgeExtractorCorners.C
+++ b/meshLibrary/utilities/surfaceTools/edgeExtraction/edgeExtractor/edgeExtractorCorners.C
@@ -192,6 +192,188 @@ bool edgeExtractor::findCornerCandidates()
     return changed;
 }
 
+bool edgeExtractor::checkCorners()
+{
+    bool changed(false);
+
+    const pointFieldPMG& points = mesh_.points();
+    const meshSurfaceEngine& mse = this->surfaceEngine();
+    const labelList& bPoints = mse.boundaryPoints();
+    const faceList::subList& bFaces = mse.boundaryFaces();
+    const labelList& bp = mse.bp();
+    const edgeList& edges = mse.edges();
+    const VRWGraph& faceEdges = mse.faceEdges();
+    const VRWGraph& edgeFaces = mse.edgeFaces();
+
+    const triSurfacePartitioner& sPart = this->partitioner();
+    const labelList& surfEdgeIndex = sPart.edgePartitions();
+
+    //- allocate a copy of boundary patches
+    labelList newBoundaryPatches(facePatch_.size());
+
+    label nCorrected;
+    Map<label> otherProcNewPatch;
+
+    label iter(0);
+
+    do
+    {
+        nCorrected = 0;
+        newBoundaryPatches = facePatch_;
+
+        //- check whether there exist situations where a boundary face
+        //- is surrounded by more faces in different patches than the
+        //- faces in the current patch
+        if( Pstream::parRun() )
+        {
+            findOtherFacePatchesParallel
+            (
+                otherProcNewPatch,
+                &facePatch_
+            );
+        }
+
+        //- update the information which edges are assigned as feature edges
+        meshSurfacePartitioner mPart(mse, facePatch_);
+
+        //- find corners in the current constelation and their nearest
+        //- counterparts in the suface mesh
+        const std::map<label, DynList<label> >& cornerPatches =
+            mPart.cornerPatches();
+
+        std::map<label, label> cornerIndex;
+        std::map<label, std::pair<point, scalar> > nearestToCorner;
+
+        typedef std::map<label, DynList<label> > mapType;
+        forAllConstIter(mapType, cornerPatches, it)
+        {
+            const label bpI = it->first;
+            const DynList<label>& neiPatches = it->second;
+
+            const point& p = points[bPoints[bpI]];
+            point pMap;
+            scalar dSq;
+            label nsp;
+
+            if( !meshOctree_.findNearestCorner(pMap, dSq, nsp, p, neiPatches) )
+            {
+                nsp = -1;
+                meshOctree_.findNearestPointToPatches(pMap, dSq, p, neiPatches);
+            }
+
+            nearestToCorner[bpI] = std::make_pair(pMap, dSq);
+            cornerIndex[bpI] = nsp;
+        }
+
+        //- for all edge nodes find their nearest counterparts in the surface
+        const labelHashSet& featureEdge = mPart.featureEdges();
+
+        std::map<label, std::pair<point, scalar> > nearestToEdgePoint;
+        std::map<label, label> edgePointIndex;
+
+        forAllConstIter(labelHashSet, featureEdge, it)
+        {
+            const label beI = it.key();
+
+            //- get the patches bounded by this feature edge
+            DynList<label> patches(2);
+
+            if( edgeFaces.sizeOfRow(beI) == 2 )
+            {
+                patches[0] = facePatch_[edgeFaces(beI, 0)];
+                patches[1] = facePatch_[edgeFaces(beI, 1)];
+            }
+            else if( edgeFaces.sizeOfRow(beI) == 1 )
+            {
+                patches[0] = facePatch_[edgeFaces(beI, 0)];
+                patches[1] = otherProcNewPatch[beI];
+            }
+
+            //- check if some points have alreay been checked
+            const edge& e = edges[beI];
+            const label bps = bp[e.start()];
+            const label bpe = bp[e.end()];
+
+            DynList<label> checkPoints;
+            if
+            (
+                (cornerPatches.find(bps) == cornerPatches.end()) &&
+                (nearestToEdgePoint.find(bps) == nearestToEdgePoint.end())
+            )
+                checkPoints.append(bps);
+
+            if
+            (
+                (cornerPatches.find(bpe) == cornerPatches.end()) &&
+                (nearestToEdgePoint.find(bpe) == nearestToEdgePoint.end())
+            )
+                checkPoints.append(bpe);
+
+            forAll(checkPoints, i)
+            {
+                const point& p = points[bPoints[checkPoints[i]]];
+                point pMap;
+                scalar dSq;
+                label nse;
+
+                if( !meshOctree_.findNearestEdgePoint(pMap, dSq, nse, p, patches) )
+                {
+                    nse = -1;
+                    meshOctree_.findNearestPointToPatches(pMap, dSq, p, patches);
+                }
+
+                nearestToEdgePoint[checkPoints[i]] = std::make_pair(pMap, dSq);
+                edgePointIndex[checkPoints[i]] = nse;
+            }
+        }
+
+
+        forAll(bFaces, bfI)
+        {
+            const face& bf = bFaces[bfI];
+
+            DynList<label> neiPatch(bf.size(), facePatch_[bfI]);
+
+            forAll(bf, eI)
+            {
+                const label beI = faceEdges(bfI, eI);
+
+                if( featureEdge.found(beI) )
+                {
+                    if( edgeFaces.sizeOfRow(beI) == 2 )
+                    {
+                        label otherFace = edgeFaces(beI, 0);
+                        if( otherFace == bfI )
+                            otherFace = edgeFaces(beI, 1);
+
+                        neiPatch[eI] = facePatch_[otherFace];
+                    }
+                    else if( edgeFaces.sizeOfRow(beI) == 1 )
+                    {
+                        neiPatch[eI] = otherProcNewPatch[beI];
+                    }
+                }
+            }
+        }
+
+        reduce(nCorrected, sumOp<label>());
+
+        //::exit(0);
+
+        if( nCorrected )
+        {
+            changed = true;
+            facePatch_ = newBoundaryPatches;
+        }
+
+        if( ++iter > 0 )
+            break;
+
+    } while( nCorrected != 0 );
+
+    return changed;
+}
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *//
 
 } // End namespace Foam
diff --git a/meshLibrary/utilities/surfaceTools/meshSurfaceCutter/meshSurfaceCutterInternalFaces.C b/meshLibrary/utilities/surfaceTools/meshSurfaceCutter/meshSurfaceCutterInternalFaces.C
index c67e3af7b810eadfd257192534deccf1ff3b613d..0abf4d8b2c54bb4cacc212c93975487fce3dab8f 100644
--- a/meshLibrary/utilities/surfaceTools/meshSurfaceCutter/meshSurfaceCutterInternalFaces.C
+++ b/meshLibrary/utilities/surfaceTools/meshSurfaceCutter/meshSurfaceCutterInternalFaces.C
@@ -32,12 +32,6 @@ Description
 
 //#define DEBUGCutter
 
-# ifdef DEBUGCutter
-#include "triSurfSearch.H"
-#include "octreeDataTriSurface.H"
-#include "octree.H"
-# endif
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
diff --git a/meshLibrary/utilities/surfaceTools/meshSurfaceDetectPlanarRegions/meshSurfaceDetectPlanarRegions.C b/meshLibrary/utilities/surfaceTools/meshSurfaceDetectPlanarRegions/meshSurfaceDetectPlanarRegions.C
new file mode 100644
index 0000000000000000000000000000000000000000..d83c8203022287fd0277470aa9716b303ff92e70
--- /dev/null
+++ b/meshLibrary/utilities/surfaceTools/meshSurfaceDetectPlanarRegions/meshSurfaceDetectPlanarRegions.C
@@ -0,0 +1,242 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2005-2007 Franjo Juretic
+     \\/     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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+Description
+
+\*---------------------------------------------------------------------------*/
+
+#include "meshSurfaceDetectPlanarRegions.H"
+#include "meshSurfaceEngine.H"
+#include "labelledPoint.H"
+#include "helperFunctions.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace meshSurfaceDetectPlanarRegionsHelper
+{
+
+class meshSurfaceDetectPlanarRegionsNeiOp
+{
+    // Private data
+        //- const reference to meshSurfaceEngine
+        const meshSurfaceEngine& mse_;
+
+        //- angle tolerance
+        const scalar cosTol_;
+
+        //- face normals
+        vectorField faceNormals_;
+
+public:
+
+    // Constructors
+        //- Construct from meshSurfaceEngine
+        meshSurfaceDetectPlanarRegionsNeiOp
+        (
+            const meshSurfaceEngine& mse,
+            const scalar cosTol
+        )
+        :
+            mse_(mse),
+            cosTol_(cosTol),
+            faceNormals_()
+        {
+            //- create face-faces addressing outside of a parallel region
+            mse_.faceFaces();
+
+            //- calculate face normals
+            const pointFieldPMG& pts = mse_.points();
+            const faceList::subList& bFaces = mse_.boundaryFaces();
+            faceNormals_.setSize(bFaces.size());
+
+            # ifdef USE_OMP
+            # pragma omp parallel for schedule(dynamic, 50)
+            # endif
+            forAll(bFaces, bfI)
+            {
+                const face& bf = bFaces[bfI];
+
+                vector& fn = faceNormals_[bfI];
+                fn = bf.normal(pts);
+                fn /= mag(fn);
+            }
+        }
+
+    // Public member functions
+        //- return the size (number of boundary faces)
+        inline label size() const
+        {
+            return mse_.boundaryFaces().size();
+        }
+
+        //- find neighbours of a boundary face
+        void operator()(const label bfI, DynList<label>& neiFaces) const
+        {
+            neiFaces.clear();
+            const VRWGraph& faceFaces = mse_.faceFaces();
+
+            const vector& fn = faceNormals_[bfI];
+            forAllRow(faceFaces, bfI, ffI)
+            {
+                const vector& on = faceNormals_[faceFaces(bfI, ffI)];
+
+                if( (on & fn) > cosTol_ )
+                    neiFaces.append(faceFaces(bfI, ffI));
+            }
+        }
+
+        //- unify region in case of an MPI parallel run
+        template<class labelListType>
+        void collectGroups
+        (
+            std::map<label, DynList<label> >& neiGroups,
+            const labelListType& elementInGroup,
+            const DynList<label>& localGroupLabel
+        ) const
+        {
+        /*
+            const Map<label>& globalToLocal =
+                mse_.globalToLocalBndEdgeAddressing();
+            const Map<label>& otherFaceAtProc = mse_.otherEdgeFaceAtProc();
+
+            const DynList<label>& neiProcs = mse_.beNeiProcs();
+
+            std::map<label, DynList<labelledPoint> > exchangeData;
+            forAll(neiProcs, i)
+                exchangeData[neiProcs[i]].clear();
+
+            forAllConstIter(Map<label>, globalToLocal, it)
+            {
+                const label beI = it();
+
+                forAllRow(bpEdges, bpI, i)
+                {
+                    const label beI = bpEdges(bpI, i);
+
+                    if( !isFeatureEdge_[beI] )
+                        continue;
+
+                    const label groupI = elementInGroup[beI];
+
+                    forAllRow(bpAtProcs, bpI, ppI)
+                    {
+                        const label neiProc = bpAtProcs(bpI, ppI);
+
+                        if( neiProc == Pstream::myProcNo() )
+                            continue;
+
+                        exchangeData[neiProc].append
+                        (
+                            labelPair(it.key(), localGroupLabel[groupI])
+                        );
+                    }
+                }
+            }
+
+            LongList<labelPair> receivedData;
+            help::exchangeMap(exchangeData, receivedData);
+
+            forAll(receivedData, i)
+            {
+                const labelPair& lp = receivedData[i];
+                const label groupI = elementInGroup[globalToLocal[lp.first()]];
+
+                DynList<label>& ng = neiGroups[localGroupLabel[groupI]];
+
+                //- store the connection over the inter-processor boundary
+                ng.appendIfNotIn(lp.second());
+            }
+            */
+        }
+};
+
+class meshSurfaceDetectPlanarRegionsSelOp
+{
+public:
+
+    // Constructors
+        //- Null construct
+        meshSurfaceDetectPlanarRegionsSelOp()
+        {}
+
+    // Public member functions
+        //- a dummy operator
+        bool operator()(const label) const
+        {
+            return true;
+        }
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace meshSurfaceDetectPlanarRegionsHelper
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+void meshSurfaceDetectPlanarRegions::findPlanarRegions()
+{
+    const scalar cosTol = Foam::cos(angleTol_);
+
+    meshSurfaceDetectPlanarRegionsHelper::meshSurfaceDetectPlanarRegionsNeiOp nop
+    (
+        meshSurface_,
+        cosTol
+    );
+
+    meshSurfaceDetectPlanarRegionsHelper::meshSurfaceDetectPlanarRegionsSelOp sop;
+
+    nPlanarRegions_ = help::groupMarking(faceInPlanarRegion_, nop, sop);
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+meshSurfaceDetectPlanarRegions::meshSurfaceDetectPlanarRegions
+(
+    const meshSurfaceEngine& meshSurface,
+    const scalar angleTol
+)
+:
+    meshSurface_(meshSurface),
+    angleTol_(angleTol),
+    faceInPlanarRegion_(),
+    nPlanarRegions_(0)
+{
+    findPlanarRegions();
+}
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+meshSurfaceDetectPlanarRegions::~meshSurfaceDetectPlanarRegions()
+{}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/meshLibrary/utilities/surfaceTools/meshSurfaceDetectPlanarRegions/meshSurfaceDetectPlanarRegions.H b/meshLibrary/utilities/surfaceTools/meshSurfaceDetectPlanarRegions/meshSurfaceDetectPlanarRegions.H
new file mode 100644
index 0000000000000000000000000000000000000000..d0c55f1874d0eff988f2192505b9cf05e83bddad
--- /dev/null
+++ b/meshLibrary/utilities/surfaceTools/meshSurfaceDetectPlanarRegions/meshSurfaceDetectPlanarRegions.H
@@ -0,0 +1,122 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2005-2007 Franjo Juretic
+     \\/     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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+Class
+    meshSurfaceDetectPlanarRegions
+
+Description
+    Detects planar regions in the surface of the volume mesh
+
+SourceFiles
+    meshSurfaceDetectPlanarRegions.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef meshSurfaceDetectPlanarRegions_H
+#define meshSurfaceDetectPlanarRegions_H
+
+#include "polyMeshGenModifier.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declarations
+class meshSurfaceEngine;
+
+/*---------------------------------------------------------------------------*\
+                Class meshSurfaceDetectPlanarRegions Declaration
+\*---------------------------------------------------------------------------*/
+
+class meshSurfaceDetectPlanarRegions
+{
+    // Private data
+        //- const reference to mesh surface
+        const meshSurfaceEngine& meshSurface_;
+
+        //- angle deviation tolerance
+        scalar angleTol_;
+
+        //- planar region index fo each boundary face
+        labelList faceInPlanarRegion_;
+
+        //- number of detected planar regions
+        label nPlanarRegions_;
+
+    // Private member functions
+        //- calculate new boundary and replace the existing one
+        void findPlanarRegions();
+
+        //- Disallow default bitwise copy construct
+        meshSurfaceDetectPlanarRegions(const meshSurfaceDetectPlanarRegions&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const meshSurfaceDetectPlanarRegions&);
+
+public:
+
+    // Constructors
+
+        //- Construct from mesh and IOdictionary
+        meshSurfaceDetectPlanarRegions
+        (
+            const meshSurfaceEngine& meshSurface,
+            const scalar angleTol = (M_PI / 180.0)
+        );
+
+    // Destructor
+
+        ~meshSurfaceDetectPlanarRegions();
+
+    // Member Functions
+        //- return the number of detected planar regions
+        inline label nPlanarRegions() const
+        {
+            return nPlanarRegions_;
+        }
+
+        //- returns a list containing indices of planar regions
+        inline const labelList& planarRegions() const
+        {
+            return faceInPlanarRegion_;
+        }
+
+        //- return a region index for a given face
+        inline label planarRegion(const label bfI) const
+        {
+            return faceInPlanarRegion_[bfI];
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/meshLibrary/utilities/surfaceTools/meshSurfaceEdgeExtractor/meshSurfaceEdgeCreateEdgeVertices.C b/meshLibrary/utilities/surfaceTools/meshSurfaceEdgeExtractor/meshSurfaceEdgeCreateEdgeVertices.C
index 4b235d5f9d31295b348c7cf9d9d418342b0ea972..96cd5c3464aac1fa97e5fb99ac1d97cc2a6e2a1e 100644
--- a/meshLibrary/utilities/surfaceTools/meshSurfaceEdgeExtractor/meshSurfaceEdgeCreateEdgeVertices.C
+++ b/meshLibrary/utilities/surfaceTools/meshSurfaceEdgeExtractor/meshSurfaceEdgeCreateEdgeVertices.C
@@ -81,6 +81,7 @@ void meshSurfaceEdgeExtractor::createEdgeVertices()
             {
                 point newP;
                 scalar distSq;
+                label nse;
 
                 FixedList<point, 2> edgePoints;
                 FixedList<label, 2> patches;
@@ -91,12 +92,13 @@ void meshSurfaceEdgeExtractor::createEdgeVertices()
                 patches[1] = pointRegions_(e, 0);
 
                 const bool found =
-                meshOctree_.findNearestVertexToTheEdge
+                meshOctree_.findNearestPointToEdge
                 (
-                    edgePoints,
-                    patches,
                     newP,
-                    distSq
+                    distSq,
+                    nse,
+                    edgePoints,
+                    patches
                 );
 
                 if( found )
diff --git a/meshLibrary/utilities/surfaceTools/meshSurfaceEdgeExtractor2D/meshSurfaceEdgeExtractor2DDistributeFaces.C b/meshLibrary/utilities/surfaceTools/meshSurfaceEdgeExtractor2D/meshSurfaceEdgeExtractor2DDistributeFaces.C
index 70a6d208a6d1342f8642c67434da2899c3f2fd94..163b4452bbceac16e3cba2dea9ba1ee789a51562 100644
--- a/meshLibrary/utilities/surfaceTools/meshSurfaceEdgeExtractor2D/meshSurfaceEdgeExtractor2DDistributeFaces.C
+++ b/meshLibrary/utilities/surfaceTools/meshSurfaceEdgeExtractor2D/meshSurfaceEdgeExtractor2DDistributeFaces.C
@@ -125,8 +125,8 @@ void meshSurfaceEdgeExtractor2D::distributeBoundaryFaces()
             //- find the patch index of the nearest location on the surface mesh
             point mapPoint;
             scalar distSq;
-            label patchI;
-            meshOctree_.findNearestSurfacePoint(mapPoint, distSq, patchI, c);
+            label patchI, nt;
+            meshOctree_.findNearestSurfacePoint(mapPoint, distSq, nt, patchI, c);
 
             bndFacePatch[bfI] = patchI;
         }
diff --git a/meshLibrary/utilities/surfaceTools/meshSurfaceEdgeExtractorFUN/meshSurfaceEdgeExtractorFUNDistributeFaces.C b/meshLibrary/utilities/surfaceTools/meshSurfaceEdgeExtractorFUN/meshSurfaceEdgeExtractorFUNDistributeFaces.C
index 9f7945dcbfebb8abea4ececbb818f8b4cc510b2e..42626222828b6395adc9652e4694f9119ce627b8 100644
--- a/meshLibrary/utilities/surfaceTools/meshSurfaceEdgeExtractorFUN/meshSurfaceEdgeExtractorFUNDistributeFaces.C
+++ b/meshLibrary/utilities/surfaceTools/meshSurfaceEdgeExtractorFUN/meshSurfaceEdgeExtractorFUNDistributeFaces.C
@@ -39,7 +39,9 @@ Description
 #include "meshSurfaceCheckEdgeTypes.H"
 #include "meshSurfaceEngine.H"
 
+# ifdef USE_OMP
 #include <omp.h>
+# endif
 
 //#define DEBUGMapping
 
@@ -87,11 +89,11 @@ void meshSurfaceEdgeExtractorFUN::distributeBoundaryFaces()
     {
         const point c = bFaces[bfI].centre(points);
 
-        label facePatch;
+        label facePatch, nt;
         point p;
         scalar distSq;
 
-        meshOctree_.findNearestSurfacePoint(p, distSq, facePatch, c);
+        meshOctree_.findNearestSurfacePoint(p, distSq, nt, facePatch, c);
 
         if( (facePatch > -1) && (facePatch < nPatches) )
         {
diff --git a/meshLibrary/utilities/surfaceTools/meshSurfaceMapper/meshSurfaceMapperCornersAndEdges.C b/meshLibrary/utilities/surfaceTools/meshSurfaceMapper/meshSurfaceMapperCornersAndEdges.C
index 14bc754f566e5edcf8f91eda23edf63ae5e038ac..d8949fd1d4d254841783ff90db8884fa5469c56a 100644
--- a/meshLibrary/utilities/surfaceTools/meshSurfaceMapper/meshSurfaceMapperCornersAndEdges.C
+++ b/meshLibrary/utilities/surfaceTools/meshSurfaceMapper/meshSurfaceMapperCornersAndEdges.C
@@ -187,10 +187,12 @@ void meshSurfaceMapper::mapCorners(const labelLongList& nodesToMap)
             forAllRow(pPatches, bpI, patchI)
             {
                 point np;
+                label nt;
                 meshOctree_.findNearestSurfacePointInRegion
                 (
                     np,
                     distSqApprox,
+                    nt,
                     pPatches(bpI, patchI),
                     mapPointApprox
                 );
@@ -292,10 +294,12 @@ void meshSurfaceMapper::mapEdgeNodes(const labelLongList& nodesToMap)
             forAll(patches, patchI)
             {
                 point np;
+                label nt;
                 meshOctree_.findNearestSurfacePointInRegion
                 (
                     np,
                     distSqApprox,
+                    nt,
                     patches[patchI],
                     mapPointApprox
                 );
@@ -314,7 +318,8 @@ void meshSurfaceMapper::mapEdgeNodes(const labelLongList& nodesToMap)
         //- find the nearest vertex on the triSurface feature edge
         point mapPoint;
         scalar distSq;
-        meshOctree_.findNearestEdgePoint(p, patches, mapPoint, distSq);
+        label nse;
+        meshOctree_.findNearestEdgePoint(mapPoint, distSq, nse, p, patches);
 
         //- use the vertex with the smallest mapping distance
         if( distSq > 1.2 * distSqApprox )
@@ -370,7 +375,7 @@ void meshSurfaceMapper::mapCornersAndEdges()
     mapCorners(selectedPoints);
 
     selectedPoints.clear();
-    const labelHashSet& edgePoints = mPart.edgeNodes();
+    const labelHashSet& edgePoints = mPart.edgePoints();
     forAllConstIter(labelHashSet, edgePoints, it)
         selectedPoints.append(it.key());
 
diff --git a/meshLibrary/utilities/surfaceTools/meshSurfaceMapper/meshSurfaceMapperMapVertices.C b/meshLibrary/utilities/surfaceTools/meshSurfaceMapper/meshSurfaceMapperMapVertices.C
index 01bc0636b85864e967e9e253dfba32c54347816e..9cae0a365b4e54b1af2c7a2ba097089057fef731 100644
--- a/meshLibrary/utilities/surfaceTools/meshSurfaceMapper/meshSurfaceMapperMapVertices.C
+++ b/meshLibrary/utilities/surfaceTools/meshSurfaceMapper/meshSurfaceMapperMapVertices.C
@@ -37,7 +37,9 @@ Description
 
 #include <map>
 
+# ifdef USE_OMP
 #include <omp.h>
+# endif
 
 //#define DEBUGMapping
 
@@ -170,7 +172,7 @@ void meshSurfaceMapper::mapToSmallestDistance(LongList<parMapperHelper>& parN)
 
 void meshSurfaceMapper::mapNodeToPatch(const label bpI, const label patchI)
 {
-    label patch;
+    label patch, nt;
     point mapPoint;
     scalar dSq;
 
@@ -180,11 +182,11 @@ void meshSurfaceMapper::mapNodeToPatch(const label bpI, const label patchI)
 
     if( patchI < 0 )
     {
-        meshOctree_.findNearestSurfacePoint(mapPoint, dSq, patch, p);
+        meshOctree_.findNearestSurfacePoint(mapPoint, dSq, nt, patch, p);
     }
     else
     {
-        meshOctree_.findNearestSurfacePointInRegion(mapPoint, dSq, patchI, p);
+        meshOctree_.findNearestSurfacePointInRegion(mapPoint, dSq, nt, patchI, p);
     }
 
     meshSurfaceEngineModifier surfModifier(surfaceEngine_);
@@ -230,7 +232,7 @@ void meshSurfaceMapper::mapVerticesOntoSurface(const labelLongList& nodesToMap)
             << points[boundaryPoints[bpI]] << endl;
         # endif
 
-        label patch;
+        label patch, nt;
         point mapPoint;
         scalar dSq;
 
@@ -238,6 +240,7 @@ void meshSurfaceMapper::mapVerticesOntoSurface(const labelLongList& nodesToMap)
         (
             mapPoint,
             dSq,
+            nt,
             patch,
             points[boundaryPoints[bpI]]
         );
@@ -293,7 +296,7 @@ void meshSurfaceMapper::mapVerticesOntoSurfacePatches
 {
     const meshSurfacePartitioner& mPart = meshPartitioner();
     const labelHashSet& cornerPoints = mPart.corners();
-    const labelHashSet& edgePoints = mPart.edgeNodes();
+    const labelHashSet& edgePoints = mPart.edgePoints();
 
     boolList treatedPoint(surfaceEngine_.boundaryPoints().size(), false);
 
@@ -340,10 +343,12 @@ void meshSurfaceMapper::mapVerticesOntoSurfacePatches
         const point& p = points[bPoints[bpI]];
         point mapPoint;
         scalar dSq;
+        label nt;
         meshOctree_.findNearestSurfacePointInRegion
         (
             mapPoint,
             dSq,
+            nt,
             pointPatches(bpI, 0),
             p
         );
diff --git a/meshLibrary/utilities/surfaceTools/meshSurfaceMapper/meshSurfaceMapperPremapVertices.C b/meshLibrary/utilities/surfaceTools/meshSurfaceMapper/meshSurfaceMapperPremapVertices.C
index 4be686ff5237fd2cccc652d1038ee44b60aba4e7..2a3e5a519bbb0aae1e3ce312c89e54a52edad468 100644
--- a/meshLibrary/utilities/surfaceTools/meshSurfaceMapper/meshSurfaceMapperPremapVertices.C
+++ b/meshLibrary/utilities/surfaceTools/meshSurfaceMapper/meshSurfaceMapperPremapVertices.C
@@ -140,11 +140,14 @@ void meshSurfaceMapper::preMapVertices(const label nIterations)
             }
         }
 
-        //- calculate coordinates of points for searching
+        //- create the surface modifier and move the surface points
+        meshSurfaceEngineModifier surfaceModifier(surfaceEngine_);
+        LongList<parMapperHelper> parallelBndNodes;
+
         # ifdef USE_OMP
-        # pragma omp parallel for
+        # pragma omp parallel for schedule(dynamic, 50)
         # endif
-        forAll(preMapPositions, bpI)
+        forAll(boundaryPoints, bpI)
         {
             labelledPoint& lp = preMapPositions[bpI];
 
@@ -156,20 +159,10 @@ void meshSurfaceMapper::preMapVertices(const label nIterations)
             }
 
             lp.coordinates() /= lp.pointLabel();
-        }
-
-        //- create the surface modifier and move the surface points
-        meshSurfaceEngineModifier surfaceModifier(surfaceEngine_);
-        LongList<parMapperHelper> parallelBndNodes;
 
-        # ifdef USE_OMP
-        # pragma omp parallel for schedule(dynamic, 50)
-        # endif
-        forAll(boundaryPoints, bpI)
-        {
             const point& p = points[boundaryPoints[bpI]];
 
-            label patch;
+            label patch, nearestTri;
             point pMap = p;
             scalar dSq;
 
@@ -177,8 +170,9 @@ void meshSurfaceMapper::preMapVertices(const label nIterations)
             (
                 pMap,
                 dSq,
+                nearestTri,
                 patch,
-                preMapPositions[bpI].coordinates()
+                lp.coordinates()
             );
 
             const point newP = 0.5 * (pMap + p);
diff --git a/meshLibrary/utilities/surfaceTools/meshSurfaceMapper2D/meshSurfaceMapper2DMapVertices.C b/meshLibrary/utilities/surfaceTools/meshSurfaceMapper2D/meshSurfaceMapper2DMapVertices.C
index 471e93c78b3b9cb96efe05919bd6e978d54946de..f0e0a5d4d05aa347b46f576bfdfbccc8b06c62e6 100644
--- a/meshLibrary/utilities/surfaceTools/meshSurfaceMapper2D/meshSurfaceMapper2DMapVertices.C
+++ b/meshLibrary/utilities/surfaceTools/meshSurfaceMapper2D/meshSurfaceMapper2DMapVertices.C
@@ -306,7 +306,7 @@ void meshSurfaceMapper2D::mapVerticesOntoSurface(const labelLongList& edgesToMap
              << e << endl;
         # endif
 
-        label patch;
+        label patch, nt;
         point mapPoint;
         scalar dSq;
 
@@ -314,6 +314,7 @@ void meshSurfaceMapper2D::mapVerticesOntoSurface(const labelLongList& edgesToMap
         (
             mapPoint,
             dSq,
+            nt,
             patch,
             points[e.start()]
         );
@@ -426,10 +427,12 @@ void meshSurfaceMapper2D::mapCorners(const labelLongList& edgesToMap)
             forAll(ePatches, epI)
             {
                 point np;
+                label nt;
                 meshOctree_.findNearestSurfacePointInRegion
                 (
                     np,
                     distSqApprox,
+                    nt,
                     ePatches[epI],
                     mapPointApprox
                 );
@@ -449,7 +452,8 @@ void meshSurfaceMapper2D::mapCorners(const labelLongList& edgesToMap)
         //- find the nearest triSurface corner for the given corner
         scalar distSq;
         point mapPoint;
-        meshOctree_.findNearestEdgePoint(p, ePatches, mapPoint, distSq);
+        label nse;
+        meshOctree_.findNearestEdgePoint(mapPoint, distSq, nse, p, ePatches);
 
         if( distSq > mappingDistance[beI] )
         {
@@ -537,11 +541,13 @@ void meshSurfaceMapper2D::mapVerticesOntoSurfacePatches
         const point& p = points[e.start()];
         point mapPoint;
         scalar dSq;
+        label nt;
 
         meshOctree_.findNearestSurfacePointInRegion
         (
             mapPoint,
             dSq,
+            nt,
             edgePatches(beI, 0),
             p
         );
diff --git a/meshLibrary/utilities/surfaceTools/meshSurfaceMapper2D/meshSurfaceMapper2DPremapVertices.C b/meshLibrary/utilities/surfaceTools/meshSurfaceMapper2D/meshSurfaceMapper2DPremapVertices.C
index fbfa70b464a9f175bd8d80d66b8c83b104672fc3..7eac99c311a7ab5d18b63d3975e057939ff37db2 100644
--- a/meshLibrary/utilities/surfaceTools/meshSurfaceMapper2D/meshSurfaceMapper2DPremapVertices.C
+++ b/meshLibrary/utilities/surfaceTools/meshSurfaceMapper2D/meshSurfaceMapper2DPremapVertices.C
@@ -191,7 +191,7 @@ void meshSurfaceMapper2D::preMapVertices(const label nIterations)
 
             const point& p = points[e.start()];
 
-            label patch;
+            label patch, nt;
             point pMap = p;
             scalar dSq;
 
@@ -199,6 +199,7 @@ void meshSurfaceMapper2D::preMapVertices(const label nIterations)
             (
                 pMap,
                 dSq,
+                nt,
                 patch,
                 preMapPositions[eI].coordinates()
             );
diff --git a/meshLibrary/utilities/surfaceTools/meshSurfacePartitioner/meshSurfacePartitioner.C b/meshLibrary/utilities/surfaceTools/meshSurfacePartitioner/meshSurfacePartitioner.C
index 936c04ca7454ff0ceb4401362579604aec027b9d..b46d49926188c2b38a1d4c9f4a4f21c13c67c5e3 100644
--- a/meshLibrary/utilities/surfaceTools/meshSurfacePartitioner/meshSurfacePartitioner.C
+++ b/meshLibrary/utilities/surfaceTools/meshSurfacePartitioner/meshSurfacePartitioner.C
@@ -35,15 +35,34 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from meshSurfaceEngine. Holds reference!
 meshSurfacePartitioner::meshSurfacePartitioner
 (
     const meshSurfaceEngine& meshSurface
 )
 :
     meshSurface_(meshSurface),
+    facePatch_(meshSurface.boundaryFacePatches()),
     corners_(),
-    edgeNodes_(),
+    cornerPatches_(),
+    edgePoints_(),
+    partitionPartitions_(),
+    nEdgesAtPoint_(),
+    featureEdges_()
+{
+    calculateCornersEdgesAndAddressing();
+}
+
+meshSurfacePartitioner::meshSurfacePartitioner
+(
+    const meshSurfaceEngine& meshSurface,
+    const labelList& facePatch
+)
+:
+    meshSurface_(meshSurface),
+    facePatch_(facePatch),
+    corners_(),
+    cornerPatches_(),
+    edgePoints_(),
     partitionPartitions_(),
     nEdgesAtPoint_(),
     featureEdges_()
diff --git a/meshLibrary/utilities/surfaceTools/meshSurfacePartitioner/meshSurfacePartitioner.H b/meshLibrary/utilities/surfaceTools/meshSurfacePartitioner/meshSurfacePartitioner.H
index cd5a56e31b64cb8f1ff5868deccbb90091dcb4ef..3d73cc23b16483380d02613bf13abb328c84bd1d 100644
--- a/meshLibrary/utilities/surfaceTools/meshSurfacePartitioner/meshSurfacePartitioner.H
+++ b/meshLibrary/utilities/surfaceTools/meshSurfacePartitioner/meshSurfacePartitioner.H
@@ -39,6 +39,8 @@ SourceFiles
 #include "meshSurfaceEngine.H"
 #include "HashSet.H"
 
+#include <map>
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -54,11 +56,17 @@ class meshSurfacePartitioner
         //- reference to mesh surface engine
         const meshSurfaceEngine& meshSurface_;
 
+        //- pointer to face patches
+        const labelList& facePatch_;
+
         //- labels of corner nodes
         labelHashSet corners_;
 
-        //- labels of edge nodes
-        labelHashSet edgeNodes_;
+        //- patches at a corner node
+        std::map<label, DynList<label> > cornerPatches_;
+
+        //- labels of edge points
+        labelHashSet edgePoints_;
 
         //- partition-partitions addressing
         List<labelHashSet> partitionPartitions_;
@@ -83,24 +91,37 @@ public:
 
     // Constructors
 
-        //- Construct from mesh
+        //- Construct from meshSurfaceEngine
         meshSurfacePartitioner(const meshSurfaceEngine&);
 
+        //- Construct from meshSurfaceEngine and face patches
+        meshSurfacePartitioner
+        (
+            const meshSurfaceEngine& meshSurface,
+            const labelList& facePatch
+        );
+
     // Destructor
 
         ~meshSurfacePartitioner();
 
     // Member Functions
-        //- return labels of corner nodes (from the list of boundary points)
+        //- return labels of corner points (from the list of boundary points)
         inline const labelHashSet& corners() const
         {
             return corners_;
         }
 
-        //- return labels of edge nodes (from the list of boundary points)
-        inline const labelHashSet& edgeNodes() const
+        //- return label of patches near corner points
+        inline const std::map<label, DynList<label> >& cornerPatches() const
+        {
+            return cornerPatches_;
+        }
+
+        //- return labels of edge points (from the list of boundary points)
+        inline const labelHashSet& edgePoints() const
         {
-            return edgeNodes_;
+            return edgePoints_;
         }
 
         //- return the number of feature edges attached to a boundary point
diff --git a/meshLibrary/utilities/surfaceTools/meshSurfacePartitioner/meshSurfacePartitionerFunctions.C b/meshLibrary/utilities/surfaceTools/meshSurfacePartitioner/meshSurfacePartitionerFunctions.C
index a3075d1c6d60fcedb1d6edb3026dfd8ff86bef82..4c0dbc4e0fd01e67ee9316ea913fa5634dbbbd22 100644
--- a/meshLibrary/utilities/surfaceTools/meshSurfacePartitioner/meshSurfacePartitionerFunctions.C
+++ b/meshLibrary/utilities/surfaceTools/meshSurfacePartitioner/meshSurfacePartitionerFunctions.C
@@ -42,10 +42,10 @@ void meshSurfacePartitioner::calculateCornersEdgesAndAddressing()
     const labelList& bp = meshSurface_.bp();
     const edgeList& edges = meshSurface_.edges();
     const VRWGraph& edgeFaces = meshSurface_.edgeFaces();
-    const labelList& facePatches = meshSurface_.boundaryFacePatches();
+    const VRWGraph& pointFaces = meshSurface_.pointFaces();
 
     corners_.clear();
-    edgeNodes_.clear();
+    edgePoints_.clear();
     partitionPartitions_.setSize(meshSurface_.mesh().boundaries().size());
 
     nEdgesAtPoint_.clear();
@@ -59,8 +59,8 @@ void meshSurfacePartitioner::calculateCornersEdgesAndAddressing()
         if( edgeFaces.sizeOfRow(edgeI) != 2 )
             continue;
 
-        const label patch0 = facePatches[edgeFaces(edgeI, 0)];
-        const label patch1 = facePatches[edgeFaces(edgeI, 1)];
+        const label patch0 = facePatch_[edgeFaces(edgeI, 0)];
+        const label patch1 = facePatch_[edgeFaces(edgeI, 1)];
 
         if( patch0 != patch1 )
         {
@@ -78,7 +78,38 @@ void meshSurfacePartitioner::calculateCornersEdgesAndAddressing()
     if( Pstream::parRun() )
     {
         const Map<label>& otherFaceAtProc = meshSurface_.otherEdgeFaceAtProc();
-        const Map<label>& otherFacePatch = meshSurface_.otherEdgeFacePatch();
+
+        //- find patches on other procs sharing surface edges
+        Map<label> otherFacePatch;
+
+        const DynList<label>& beNeiProcs = meshSurface_.beNeiProcs();
+        const Map<label>& globalToLocalEdges =
+            meshSurface_.globalToLocalBndEdgeAddressing();
+
+        std::map<label, labelLongList> exchangeData;
+        forAll(beNeiProcs, i)
+            exchangeData[beNeiProcs[i]].clear();
+
+        forAllConstIter(Map<label>, globalToLocalEdges, it)
+        {
+            const label beI = it();
+
+            labelLongList& data = exchangeData[otherFaceAtProc[beI]];
+
+            data.append(it.key());
+            data.append(facePatch_[edgeFaces(beI, 0)]);
+        }
+
+        labelLongList receivedData;
+        help::exchangeMap(exchangeData, receivedData);
+
+        for(label i=0;i<receivedData.size();)
+        {
+            const label geI = receivedData[i++];
+            const label patchI = receivedData[i++];
+
+            otherFacePatch.insert(globalToLocalEdges[geI], patchI);
+        }
 
         //- take into account feature edges at processor boundaries
         forAllConstIter(Map<label>, otherFaceAtProc, it)
@@ -87,7 +118,7 @@ void meshSurfacePartitioner::calculateCornersEdgesAndAddressing()
 
             if( it() <= Pstream::myProcNo() )
                 continue;
-            if( otherFacePatch[beI] != facePatches[edgeFaces(beI, 0)] )
+            if( otherFacePatch[beI] != facePatch_[edgeFaces(beI, 0)] )
             {
                 const edge& e = edges[beI];
                 ++nEdgesAtPoint_[bp[e.start()]];
@@ -96,7 +127,7 @@ void meshSurfacePartitioner::calculateCornersEdgesAndAddressing()
         }
 
         //- gather data on all processors
-        std::map<label, labelLongList> exchangeData;
+        exchangeData.clear();
         const DynList<label>& bpNeiProcs = meshSurface_.bpNeiProcs();
         forAll(bpNeiProcs, i)
             exchangeData.insert
@@ -129,7 +160,7 @@ void meshSurfacePartitioner::calculateCornersEdgesAndAddressing()
         }
 
         //- exchange information
-        labelLongList receivedData;
+        receivedData.clear();
         help::exchangeMap(exchangeData, receivedData);
 
         //- add the edges from other processors to the points
@@ -143,15 +174,79 @@ void meshSurfacePartitioner::calculateCornersEdgesAndAddressing()
         }
     }
 
+    //- mark edges and corners
     forAll(nEdgesAtPoint_, bpI)
     {
         if( nEdgesAtPoint_[bpI] > 2 )
         {
             corners_.insert(bpI);
+            cornerPatches_[bpI].clear();
         }
         else if( nEdgesAtPoint_[bpI] == 2 )
         {
-            edgeNodes_.insert(bpI);
+            edgePoints_.insert(bpI);
+        }
+    }
+
+    //- find patches attached to corners
+    forAllConstIter(labelHashSet, corners_, it)
+    {
+        const label bpI = it.key();
+
+        forAllRow(pointFaces, bpI, pfI)
+        {
+            cornerPatches_[bpI].appendIfNotIn(facePatch_[pointFaces(bpI, pfI)]);
+        }
+    }
+
+    if( Pstream::parRun() )
+    {
+        const Map<label>& globalToLocal =
+            meshSurface_.globalToLocalBndPointAddressing();
+        const DynList<label>& bpNeiProcs = meshSurface_.bpNeiProcs();
+        const VRWGraph& bpAtProcs = meshSurface_.bpAtProcs();
+
+        std::map<label, labelLongList> exchangeData;
+        forAll(bpNeiProcs, i)
+            exchangeData[bpNeiProcs[i]].clear();
+
+        forAllConstIter(Map<label>, globalToLocal, it)
+        {
+            const label bpI = it();
+
+            if( corners_.found(bpI) )
+            {
+                const DynList<label>& cPatches = cornerPatches_[bpI];
+
+                forAllRow(bpAtProcs, bpI, i)
+                {
+                    const label neiProc = bpAtProcs(bpI, i);
+
+                    if( neiProc == Pstream::myProcNo() )
+                        continue;
+
+                    labelLongList& data = exchangeData[neiProc];
+
+                    data.append(it.key());
+                    data.append(cPatches.size());
+                    forAll(cPatches, ppI)
+                        data.append(cPatches[ppI]);
+                }
+            }
+        }
+
+        //- exchange data with other prcessors
+        labelLongList receivedData;
+        help::exchangeMap(exchangeData, receivedData);
+
+        label counter(0);
+        while( counter < receivedData.size() )
+        {
+            const label bpI = globalToLocal[receivedData[counter++]];
+            const label size = receivedData[counter++];
+
+            for(label i=0;i<size;++i)
+                cornerPatches_[bpI].appendIfNotIn(receivedData[counter++]);
         }
     }
 
@@ -159,16 +254,10 @@ void meshSurfacePartitioner::calculateCornersEdgesAndAddressing()
     reduce(counter, sumOp<label>());
     Info << "Found " << counter
         << " corners at the surface of the volume mesh" << endl;
-    counter = edgeNodes_.size();
+    counter = edgePoints_.size();
     reduce(counter, sumOp<label>());
     Info << "Found " << counter
         << " edge points at the surface of the volume mesh" << endl;
-
-//    const pointFieldPMG& points = meshSurface_.points();
-//    forAllConstIter(labelHashSet, corners_, it)
-//        Info << "Corner point " << it.key() << " has coordinates " << points[bPoints[it.key()]] << endl;
-//    forAllConstIter(labelHashSet, edgeNodes_, it)
-//        Info << "Edge point " << it.key() << " has coordinates " << points[bPoints[it.key()]] << endl;
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/meshLibrary/utilities/triSurfaceTools/triSurfaceCleanupDuplicates/triSurfaceCleanupDuplicatesFunctions.C b/meshLibrary/utilities/triSurfaceTools/triSurfaceCleanupDuplicates/triSurfaceCleanupDuplicatesFunctions.C
index 45b112645eefa756850d513658951e08ee6209e0..864fce864b71c54d8e0f528a3edf5f496042f80e 100644
--- a/meshLibrary/utilities/triSurfaceTools/triSurfaceCleanupDuplicates/triSurfaceCleanupDuplicatesFunctions.C
+++ b/meshLibrary/utilities/triSurfaceTools/triSurfaceCleanupDuplicates/triSurfaceCleanupDuplicatesFunctions.C
@@ -31,7 +31,9 @@ Description
 #include "meshOctree.H"
 #include "demandDrivenData.H"
 
+# ifdef USE_OMP
 #include <omp.h>
+# endif
 #include <set>
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/meshLibrary/utilities/triSurfaceTools/triSurfacePartitioner/triSurfacePartitionerCreateAddressing.C b/meshLibrary/utilities/triSurfaceTools/triSurfacePartitioner/triSurfacePartitionerCreateAddressing.C
index aaf4d854df19b4ecbdf421f0deb0e1fddeac1d6a..6e152d2121f7d1299870bef1cef76028238e2f64 100644
--- a/meshLibrary/utilities/triSurfaceTools/triSurfacePartitioner/triSurfacePartitionerCreateAddressing.C
+++ b/meshLibrary/utilities/triSurfaceTools/triSurfacePartitioner/triSurfacePartitionerCreateAddressing.C
@@ -159,7 +159,7 @@ void triSurfacePartitioner::calculateEdgePartitions()
     forAll(corners_, i)
         corners.insert(corners_[i]);
 
-    edgePartitions_.setSize(edgeFaces. size());
+    edgePartitions_.setSize(edgeFaces.size());
     edgePartitions_ = -1;
 
     label nPartitions(0);
diff --git a/meshLibrary/voronoiMesh/voronoiMeshGenerator/voronoiMeshGenerator.C b/meshLibrary/voronoiMesh/voronoiMeshGenerator/voronoiMeshGenerator.C
index 0ee39c49b5f4bb5dd5948b4cd78b8af2e88cd155..3278825ff42864689531f8b7ec3b97de4c5cd7e9 100644
--- a/meshLibrary/voronoiMesh/voronoiMeshGenerator/voronoiMeshGenerator.C
+++ b/meshLibrary/voronoiMesh/voronoiMeshGenerator/voronoiMeshGenerator.C
@@ -128,7 +128,7 @@ void voronoiMeshGenerator::mapMeshToSurface()
     # endif
 
     //- untangle surface faces
-    meshSurfaceOptimizer(*msePtr, *octreePtr_).preOptimizeSurface();
+    meshSurfaceOptimizer(*msePtr, *octreePtr_).untangleSurface();
 
     # ifdef DEBUG
     mesh_.write();
@@ -152,7 +152,9 @@ void voronoiMeshGenerator::mapEdgesAndCorners()
 void voronoiMeshGenerator::optimiseMeshSurface()
 {
     meshSurfaceEngine mse(mesh_);
-    meshSurfaceOptimizer(mse, *octreePtr_).optimizeSurface();
+    meshSurfaceOptimizer surfOptimiser(mse, *octreePtr_);
+    surfOptimiser.optimizeSurface();
+    surfOptimiser.untangleSurface();
 
     # ifdef DEBUG
     mesh_.write();
diff --git a/testingInterfaces/testCartesianMeshExtractor/Make/options b/testingInterfaces/testCartesianMeshExtractor/Make/options
index 0fa7a4adb3c4b15013f80f1b3773f0124d78fb17..5206b3c7c0f2a6eacab0bd746cfd2488e3563452 100644
--- a/testingInterfaces/testCartesianMeshExtractor/Make/options
+++ b/testingInterfaces/testCartesianMeshExtractor/Make/options
@@ -1,2 +1,2 @@
 EXE_INC = -I$(LIB_SRC)/triSurface/lnInclude -I$(LIB_SRC)/meshTools/lnInclude -I../../meshLibrary/lnInclude
-EXE_LIBS = -lMeshLibrary -ltriSurface -lmeshTools
+EXE_LIBS = -lmeshLibrary -ltriSurface -lmeshTools
diff --git a/testingInterfaces/testCartesianMeshExtractor/testCartesianMeshExtractor.C b/testingInterfaces/testCartesianMeshExtractor/testCartesianMeshExtractor.C
index b96561df471722e7019efde79fc4bd811f8ba2ff..5dae2c706fdaefbd2708af8ce4282d65f680a5bf 100644
--- a/testingInterfaces/testCartesianMeshExtractor/testCartesianMeshExtractor.C
+++ b/testingInterfaces/testCartesianMeshExtractor/testCartesianMeshExtractor.C
@@ -71,7 +71,7 @@ int main(int argc, char *argv[])
     triSurf surf(runTime.path()/surfaceFile);
 
     // construct the octree
-    meshOctree mo(surf, true);
+    meshOctree mo(surf);
     meshOctreeCreator moc(mo, meshDict);
     moc.createOctreeBoxes();
 
@@ -88,7 +88,7 @@ int main(int argc, char *argv[])
     cmg.createMesh();
 
     const pointFieldPMG& points = pmg.points();
-    forAll(points, pointI)
+/*    forAll(points, pointI)
     {
         const point p = points[pointI];
         for(label pointJ=pointI+1;pointJ<points.size();++pointJ)
@@ -96,7 +96,7 @@ int main(int argc, char *argv[])
                 Info << "Points " << pointI << " and " << pointJ
                      << " are duplicates " << endl;
     }
-
+*/
     boolList usedPoint(points.size());
     const faceListPMG& faces = pmg.faces();
     forAll(faces, faceI)
diff --git a/testingInterfaces/testSurfaceSmoothing/testSurfaceSmoothing.C b/testingInterfaces/testSurfaceSmoothing/testSurfaceSmoothing.C
index f099f39bf37e024b9892b950ace044d65dde6c7d..db2a659b60eac2758d33088dd7febfda6a5d5e23 100644
--- a/testingInterfaces/testSurfaceSmoothing/testSurfaceSmoothing.C
+++ b/testingInterfaces/testSurfaceSmoothing/testSurfaceSmoothing.C
@@ -70,14 +70,6 @@ int main(int argc, char *argv[])
 
     triSurf surf(registry.path()/surfaceFile);
 
-    if( meshDict.found("subsetFileName") )
-    {
-        fileName subsetFileName = meshDict.lookup("subsetFileName");
-        if( Pstream::parRun() )
-            subsetFileName = ".."/subsetFileName;
-        surf.readFaceSubsets(registry.path()/subsetFileName);
-    }
-
     // construct the octree
     meshOctree moc(surf);
     meshOctreeCreator(moc, meshDict).createOctreeBoxes();
@@ -85,58 +77,9 @@ int main(int argc, char *argv[])
     polyMeshGen pmg(registry);
     pmg.read();
 
-    // make sure that mesh contains all surface patches
-    if(
-        Pstream::parRun() &&
-        (pmg.patchNames().size() != surf.patches().size())
-    )
-    {
-        wordList patchNames(surf.patches().size());
-        labelList patchStart(patchNames.size(), pmg.nInternalFaces());
-        labelList nFacesInPatch(patchNames.size(), 0);
-
-        forAll(patchNames, patchI)
-            patchNames[patchI] = surf.patches()[patchI].name();
-
-        forAll(pmg.patchNames(), patchI)
-        {
-            label pos(-1);
-            forAll(patchNames, nameI)
-                if( patchNames[nameI] == pmg.patchNames()[patchI] )
-                {
-                    pos = nameI;
-                    break;
-                }
-
-            nFacesInPatch[pos] = pmg.numFacesInPatch()[patchI];
-            patchStart[pos] = pmg.patchStart()[patchI];
-
-            for(label i=pos+1;i<patchNames.size();++i)
-                patchStart[i] += nFacesInPatch[pos];
-        }
-
-        polyMeshGenModifier meshModifier(pmg);
-        meshModifier.patchNamesAccess() = patchNames;
-        meshModifier.patchStartAccess() = patchStart;
-        meshModifier.nFacesInPatchAccess() = nFacesInPatch;
-    }
-
-    meshSurfaceEngine ms(pmg);
-    meshSurfaceOptimizer mo(ms, moc);
-
-    /*
-    labelList pointRegion(pmg.points().size(), -1);
-    for(label faceI=mesh.nInternalFaces();faceI<pmg.faces().size();faceI++)
-    {
-        const face& f = pmg.faces()[faceI];
-
-        forAll(f, pI)
-            pointRegion[f[pI]] = 1;
-    }
-
-    mo.preOptimizeSurface(pointRegion);
-    */
+
     mo.optimizeSurface();
+    mo.untangleSurface();
 
     pmg.write();