From de79c6377bc0f75e9936ec98f7807944383b02d5 Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Mon, 7 Jan 2019 09:20:51 +0100
Subject: [PATCH] ENH: for-range, forAllIters() ... in meshTools/

- reduced clutter when iterating over containers
---
 .../cyclicACMIGAMGInterface.C                 |  50 +--
 .../cyclicAMIGAMGInterface.C                  |  50 +--
 .../PatchFunction1/makePatchFunction1s.C      |   2 +-
 .../PatchEdgeFaceWave/patchPatchDist.C        |  12 +-
 .../PatchEdgeFaceWave/patchPatchDist.H        |   8 +-
 src/meshTools/meshStructure/meshStructure.C   |  11 +-
 src/meshTools/meshTools/meshTools.C           |   9 +-
 .../regionCoupledBaseGAMGInterface.C          |  42 +--
 src/meshTools/regionSplit/localPointRegion.C  |  30 +-
 .../cellSources/regionToCell/regionToCell.C   |   6 +-
 .../setsToFaceZone/setsToFaceZone.C           |   4 +-
 .../intersectedSurface/intersectedSurface.C   | 161 +++------
 .../surfaceIntersectionFuncs.C                |  23 +-
 .../surfaceFeatures/surfaceFeatures.C         |  28 +-
 .../triSurfaceTools/triSurfaceTools.C         | 342 +++++++-----------
 .../twoDPointCorrector/twoDPointCorrector.C   |  47 ++-
 16 files changed, 313 insertions(+), 512 deletions(-)

diff --git a/src/meshTools/AMIInterpolation/GAMG/interfaces/cyclicACMIGAMGInterface/cyclicACMIGAMGInterface.C b/src/meshTools/AMIInterpolation/GAMG/interfaces/cyclicACMIGAMGInterface/cyclicACMIGAMGInterface.C
index efbfef120a..7b8b086fd4 100644
--- a/src/meshTools/AMIInterpolation/GAMG/interfaces/cyclicACMIGAMGInterface/cyclicACMIGAMGInterface.C
+++ b/src/meshTools/AMIInterpolation/GAMG/interfaces/cyclicACMIGAMGInterface/cyclicACMIGAMGInterface.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2013-2016 OpenFOAM Foundation
@@ -80,27 +80,22 @@ Foam::cyclicACMIGAMGInterface::cyclicACMIGAMGInterface
 
         Map<label> masterToCoarseFace(localRestrictAddressing.size());
 
-        forAll(localRestrictAddressing, ffi)
+        for (const label curMaster : localRestrictAddressing)
         {
-            label curMaster = localRestrictAddressing[ffi];
+            const auto iter = masterToCoarseFace.cfind(curMaster);
 
-            Map<label>::const_iterator fnd = masterToCoarseFace.find
-            (
-                curMaster
-            );
-
-            if (fnd == masterToCoarseFace.end())
+            if (iter.found())
             {
-                // New coarse face
-                label coarseI = dynFaceCells.size();
-                dynFaceRestrictAddressing.append(coarseI);
-                dynFaceCells.append(curMaster);
-                masterToCoarseFace.insert(curMaster, coarseI);
+                // Already have coarse face
+                dynFaceRestrictAddressing.append(iter.val());
             }
             else
             {
-                // Already have coarse face
-                dynFaceRestrictAddressing.append(fnd());
+                // New coarse face
+                const label coarseI = dynFaceCells.size();
+                dynFaceRestrictAddressing.append(coarseI);
+                dynFaceCells.append(curMaster);
+                masterToCoarseFace.insert(curMaster, coarseI);
             }
         }
 
@@ -127,26 +122,21 @@ Foam::cyclicACMIGAMGInterface::cyclicACMIGAMGInterface
 
             Map<label> masterToCoarseFace(neighbourRestrictAddressing.size());
 
-            forAll(neighbourRestrictAddressing, ffi)
+            for (const label curMaster : neighbourRestrictAddressing)
             {
-                label curMaster = neighbourRestrictAddressing[ffi];
+                const auto iter = masterToCoarseFace.cfind(curMaster);
 
-                Map<label>::const_iterator fnd = masterToCoarseFace.find
-                (
-                    curMaster
-                );
-
-                if (fnd == masterToCoarseFace.end())
+                if (iter.found())
                 {
-                    // New coarse face
-                    label coarseI = masterToCoarseFace.size();
-                    dynNbrFaceRestrictAddressing.append(coarseI);
-                    masterToCoarseFace.insert(curMaster, coarseI);
+                    // Already have coarse face
+                    dynNbrFaceRestrictAddressing.append(iter.val());
                 }
                 else
                 {
-                    // Already have coarse face
-                    dynNbrFaceRestrictAddressing.append(fnd());
+                    // New coarse face
+                    const label coarseI = masterToCoarseFace.size();
+                    dynNbrFaceRestrictAddressing.append(coarseI);
+                    masterToCoarseFace.insert(curMaster, coarseI);
                 }
             }
 
diff --git a/src/meshTools/AMIInterpolation/GAMG/interfaces/cyclicAMIGAMGInterface/cyclicAMIGAMGInterface.C b/src/meshTools/AMIInterpolation/GAMG/interfaces/cyclicAMIGAMGInterface/cyclicAMIGAMGInterface.C
index fff8ac187d..e8bae00bb2 100644
--- a/src/meshTools/AMIInterpolation/GAMG/interfaces/cyclicAMIGAMGInterface/cyclicAMIGAMGInterface.C
+++ b/src/meshTools/AMIInterpolation/GAMG/interfaces/cyclicAMIGAMGInterface/cyclicAMIGAMGInterface.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2016 OpenFOAM Foundation
@@ -80,27 +80,22 @@ Foam::cyclicAMIGAMGInterface::cyclicAMIGAMGInterface
 
         Map<label> masterToCoarseFace(localRestrictAddressing.size());
 
-        forAll(localRestrictAddressing, ffi)
+        for (const label curMaster : localRestrictAddressing)
         {
-            label curMaster = localRestrictAddressing[ffi];
+            const auto iter = masterToCoarseFace.cfind(curMaster);
 
-            Map<label>::const_iterator fnd = masterToCoarseFace.find
-            (
-                curMaster
-            );
-
-            if (fnd == masterToCoarseFace.end())
+            if (iter.found())
             {
-                // New coarse face
-                label coarseI = dynFaceCells.size();
-                dynFaceRestrictAddressing.append(coarseI);
-                dynFaceCells.append(curMaster);
-                masterToCoarseFace.insert(curMaster, coarseI);
+                // Already have coarse face
+                dynFaceRestrictAddressing.append(iter.val());
             }
             else
             {
-                // Already have coarse face
-                dynFaceRestrictAddressing.append(fnd());
+                // New coarse face
+                const label coarseI = dynFaceCells.size();
+                dynFaceRestrictAddressing.append(coarseI);
+                dynFaceCells.append(curMaster);
+                masterToCoarseFace.insert(curMaster, coarseI);
             }
         }
 
@@ -127,26 +122,21 @@ Foam::cyclicAMIGAMGInterface::cyclicAMIGAMGInterface
 
             Map<label> masterToCoarseFace(neighbourRestrictAddressing.size());
 
-            forAll(neighbourRestrictAddressing, ffi)
+            for (const label curMaster : neighbourRestrictAddressing)
             {
-                label curMaster = neighbourRestrictAddressing[ffi];
+                const auto iter = masterToCoarseFace.cfind(curMaster);
 
-                Map<label>::const_iterator fnd = masterToCoarseFace.find
-                (
-                    curMaster
-                );
-
-                if (fnd == masterToCoarseFace.end())
+                if (iter.found())
                 {
-                    // New coarse face
-                    label coarseI = masterToCoarseFace.size();
-                    dynNbrFaceRestrictAddressing.append(coarseI);
-                    masterToCoarseFace.insert(curMaster, coarseI);
+                    // Already have coarse face
+                    dynNbrFaceRestrictAddressing.append(iter.val());
                 }
                 else
                 {
-                    // Already have coarse face
-                    dynNbrFaceRestrictAddressing.append(fnd());
+                    // New coarse face
+                    const label coarseI = masterToCoarseFace.size();
+                    dynNbrFaceRestrictAddressing.append(coarseI);
+                    masterToCoarseFace.insert(curMaster, coarseI);
                 }
             }
 
diff --git a/src/meshTools/PatchFunction1/makePatchFunction1s.C b/src/meshTools/PatchFunction1/makePatchFunction1s.C
index 1b59d549fe..5e78f4bcc9 100644
--- a/src/meshTools/PatchFunction1/makePatchFunction1s.C
+++ b/src/meshTools/PatchFunction1/makePatchFunction1s.C
@@ -144,7 +144,7 @@ namespace Foam
     //            );
     //
     //        // Add the UniformValueField under the Function1 name
-    //        forAllConstIter(typename F1Type, F1Table, iter)
+    //        forAllConstIters(F1Table, iter)
     //        {
     //            //bool ok =
     //            PF1Table.insert(iter.key(), cstrIter());
diff --git a/src/meshTools/algorithms/PatchEdgeFaceWave/patchPatchDist.C b/src/meshTools/algorithms/PatchEdgeFaceWave/patchPatchDist.C
index 431d0e2487..dcaae88483 100644
--- a/src/meshTools/algorithms/PatchEdgeFaceWave/patchPatchDist.C
+++ b/src/meshTools/algorithms/PatchEdgeFaceWave/patchPatchDist.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2018-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2016 OpenFOAM Foundation
@@ -47,12 +47,6 @@ Foam::patchPatchDist::patchPatchDist
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::patchPatchDist::~patchPatchDist()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 void Foam::patchPatchDist::correct()
@@ -120,8 +114,8 @@ void Foam::patchPatchDist::correct()
     {
         const edge& e = patch_.edges()[edgeI];
         const edge meshE = edge(mp[e[0]], mp[e[1]]);
-        EdgeMap<label>::const_iterator edgeFnd = nbrEdges.find(meshE);
-        if (edgeFnd != nbrEdges.end())
+
+        if (nbrEdges.found(meshE))
         {
             initialEdges.append(edgeI);
             initialEdgesInfo.append
diff --git a/src/meshTools/algorithms/PatchEdgeFaceWave/patchPatchDist.H b/src/meshTools/algorithms/PatchEdgeFaceWave/patchPatchDist.H
index bf60139369..2d56fa3307 100644
--- a/src/meshTools/algorithms/PatchEdgeFaceWave/patchPatchDist.H
+++ b/src/meshTools/algorithms/PatchEdgeFaceWave/patchPatchDist.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2013 OpenFOAM Foundation
@@ -57,10 +57,6 @@ class patchPatchDist
 :
     public scalarField
 {
-
-
-private:
-
     // Private Member Data
 
         //- Patch to operate on
@@ -85,7 +81,7 @@ public:
 
 
     //- Destructor
-    virtual ~patchPatchDist();
+    virtual ~patchPatchDist() = default;
 
 
     // Member Functions
diff --git a/src/meshTools/meshStructure/meshStructure.C b/src/meshTools/meshStructure/meshStructure.C
index dc66a60572..6617a41000 100644
--- a/src/meshTools/meshStructure/meshStructure.C
+++ b/src/meshTools/meshStructure/meshStructure.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2013-2016 OpenFOAM Foundation
@@ -36,7 +36,7 @@ License
 
 namespace Foam
 {
-defineTypeNameAndDebug(meshStructure, 0);
+    defineTypeNameAndDebug(meshStructure, 0);
 }
 
 
@@ -327,7 +327,7 @@ void Foam::meshStructure::correct
                         label pointi = f[fp];
                         label nextPointi = f.nextLabel(fp);
 
-                        EdgeMap<label>::const_iterator fnd = pointsToEdge.find
+                        const auto fnd = pointsToEdge.cfind
                         (
                             edge
                             (
@@ -335,9 +335,10 @@ void Foam::meshStructure::correct
                                 pointData[nextPointi].data()
                             )
                         );
-                        if (fnd != pointsToEdge.end())
+
+                        if (fnd.found())
                         {
-                            faceToPatchEdgeAddressing_[facei] = fnd();
+                            faceToPatchEdgeAddressing_[facei] = fnd.val();
                             faceToPatchFaceAddressing_[facei] = 0;
                             label own = mesh.faceOwner()[facei];
                             faceLayer_[facei] = cellData[own].distance();
diff --git a/src/meshTools/meshTools/meshTools.C b/src/meshTools/meshTools/meshTools.C
index 0a421cdfcc..fd04cfecd8 100644
--- a/src/meshTools/meshTools/meshTools.C
+++ b/src/meshTools/meshTools/meshTools.C
@@ -276,10 +276,8 @@ void Foam::meshTools::writeOBJ
 {
     writeOBJ(os, bb.points());
 
-    forAll(treeBoundBox::edges, edgei)
+    for (const edge& e : treeBoundBox::edges)
     {
-        const edge& e = treeBoundBox::edges[edgei];
-
         os << "l " << (e[0] + 1) <<  ' ' << (e[1] + 1) << nl;
     }
 }
@@ -296,10 +294,9 @@ void Foam::meshTools::writeOBJ
 {
     labelHashSet usedFaces(4*cellLabels.size());
 
-    forAll(cellLabels, i)
+    for (const label celli : cellLabels)
     {
-        const cell& cFaces = cells[cellLabels[i]];
-        usedFaces.insert(cFaces);
+        usedFaces.insert(cells[celli]);
     }
 
     writeOBJ(os, faces, points, usedFaces.toc());
diff --git a/src/meshTools/regionCoupled/GAMG/interfaces/regionCoupledGAMGInterface/regionCoupledBaseGAMGInterface.C b/src/meshTools/regionCoupled/GAMG/interfaces/regionCoupledGAMGInterface/regionCoupledBaseGAMGInterface.C
index 8e64313d01..e9feb65ef3 100644
--- a/src/meshTools/regionCoupled/GAMG/interfaces/regionCoupledGAMGInterface/regionCoupledBaseGAMGInterface.C
+++ b/src/meshTools/regionCoupled/GAMG/interfaces/regionCoupledGAMGInterface/regionCoupledBaseGAMGInterface.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C) 2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2015 OpenFOAM Foundation
@@ -71,23 +71,20 @@ Foam::regionCoupledBaseGAMGInterface::regionCoupledBaseGAMGInterface
         {
             label curMaster = localRestrictAddressing[ffi];
 
-            Map<label>::const_iterator fnd = masterToCoarseFace.find
-            (
-                curMaster
-            );
+            const auto fnd = masterToCoarseFace.cfind(curMaster);
 
-            if (fnd == masterToCoarseFace.end())
+            if (fnd.found())
             {
-                // New coarse face
-                label coarseI = dynFaceCells.size();
-                dynFaceRestrictAddressing.append(coarseI);
-                dynFaceCells.append(curMaster);
-                masterToCoarseFace.insert(curMaster, coarseI);
+                // Already have coarse face
+                dynFaceRestrictAddressing.append(fnd.val());
             }
             else
             {
-                // Already have coarse face
-                dynFaceRestrictAddressing.append(fnd());
+                // New coarse face
+                const label coarseI = dynFaceCells.size();
+                dynFaceRestrictAddressing.append(coarseI);
+                dynFaceCells.append(curMaster);
+                masterToCoarseFace.insert(curMaster, coarseI);
             }
         }
 
@@ -153,22 +150,19 @@ Foam::regionCoupledBaseGAMGInterface::regionCoupledBaseGAMGInterface
                 {
                     label curMaster = nbrPatchRestrictMap[ffi];
 
-                    Map<label>::const_iterator fnd = masterToCoarseFace.find
-                    (
-                        curMaster
-                    );
+                    const auto fnd = masterToCoarseFace.cfind(curMaster);
 
-                    if (fnd == masterToCoarseFace.end())
+                    if (fnd.found())
                     {
-                        // New coarse face
-                        label coarseI = masterToCoarseFace.size();
-                        dynNbrFaceRestrictAddressing.append(coarseI);
-                        masterToCoarseFace.insert(curMaster, coarseI);
+                        // Already have coarse face
+                        dynNbrFaceRestrictAddressing.append(fnd.val());
                     }
                     else
                     {
-                        // Already have coarse face
-                        dynNbrFaceRestrictAddressing.append(fnd());
+                        // New coarse face
+                        const label coarseI = masterToCoarseFace.size();
+                        dynNbrFaceRestrictAddressing.append(coarseI);
+                        masterToCoarseFace.insert(curMaster, coarseI);
                     }
                 }
 
diff --git a/src/meshTools/regionSplit/localPointRegion.C b/src/meshTools/regionSplit/localPointRegion.C
index edbd7db15c..59c62fe507 100644
--- a/src/meshTools/regionSplit/localPointRegion.C
+++ b/src/meshTools/regionSplit/localPointRegion.C
@@ -124,9 +124,9 @@ void Foam::localPointRegion::countPointRegions
     // From faces with any duplicated point on it to local face
     meshFaceMap_.resize(meshPointMap_.size());
 
-    forAllConstIter(Map<label>, candidateFace, iter)
+    forAllConstIters(candidateFace, iter)
     {
-        label facei = iter.key();
+        const label facei = iter.key();
 
         if (!mesh.isInternalFace(facei))
         {
@@ -142,14 +142,14 @@ void Foam::localPointRegion::countPointRegions
 
             forAll(f, fp)
             {
-                label pointi = f[fp];
+                const label pointi = f[fp];
 
                 // Even points which were not candidates for splitting might
                 // be on multiple baffles that are being split so check.
 
                 if (candidatePoint[pointi])
                 {
-                    label region = minRegion[facei][fp];
+                    const label region = minRegion[facei][fp];
 
                     if (minPointRegion[pointi] == -1)
                     {
@@ -158,10 +158,10 @@ void Foam::localPointRegion::countPointRegions
                     else if (minPointRegion[pointi] != region)
                     {
                         // Multiple regions for this point. Add.
-                        Map<label>::iterator iter = meshPointMap_.find(pointi);
-                        if (iter != meshPointMap_.end())
+                        const auto iter = meshPointMap_.cfind(pointi);
+                        if (iter.found())
                         {
-                            labelList& regions = pointRegions[iter()];
+                            labelList& regions = pointRegions[iter.val()];
                             if (!regions.found(region))
                             {
                                 label sz = regions.size();
@@ -190,9 +190,9 @@ void Foam::localPointRegion::countPointRegions
 
     // Add internal faces that use any duplicated point. Can only have one
     // region!
-    forAllConstIter(Map<label>, candidateFace, iter)
+    forAllConstIters(candidateFace, iter)
     {
-        label facei = iter.key();
+        const label facei = iter.key();
 
         if (mesh.isInternalFace(facei))
         {
@@ -204,7 +204,7 @@ void Foam::localPointRegion::countPointRegions
                 // speeds up rejection.
                 if (candidatePoint[f[fp]] && meshPointMap_.found(f[fp]))
                 {
-                    label meshFaceMapI = meshFaceMap_.size();
+                    const label meshFaceMapI = meshFaceMap_.size();
                     meshFaceMap_.insert(facei, meshFaceMapI);
                 }
             }
@@ -222,7 +222,7 @@ void Foam::localPointRegion::countPointRegions
 
     // Compact minRegion
     faceRegions_.setSize(meshFaceMap_.size());
-    forAllConstIter(Map<label>, meshFaceMap_, iter)
+    forAllConstIters(meshFaceMap_, iter)
     {
         faceRegions_[iter()].labelList::transfer(minRegion[iter.key()]);
 
@@ -321,7 +321,7 @@ void Foam::localPointRegion::calcPointRegions
     // only ones using a
     // candidate point so the only ones that can be affected)
     faceList minRegion(mesh.nFaces());
-    forAllConstIter(Map<label>, candidateFace, iter)
+    forAllConstIters(candidateFace, iter)
     {
         label facei = iter.key();
         const face& f = mesh.faces()[facei];
@@ -350,7 +350,7 @@ void Foam::localPointRegion::calcPointRegions
 
         Map<label> minPointValue(128);
         label nChanged = 0;
-        forAllConstIter(Map<label>, candidateCell, iter)
+        forAllConstIters(candidateCell, iter)
         {
             minPointValue.clear();
 
@@ -368,8 +368,8 @@ void Foam::localPointRegion::calcPointRegions
 
                     forAll(f, fp)
                     {
-                        label pointi = f[fp];
-                        Map<label>::iterator iter = minPointValue.find(pointi);
+                        const label pointi = f[fp];
+                        auto iter = minPointValue.find(pointi);
 
                         if (iter == minPointValue.end())
                         {
diff --git a/src/meshTools/sets/cellSources/regionToCell/regionToCell.C b/src/meshTools/sets/cellSources/regionToCell/regionToCell.C
index e4a304e282..f6803c4367 100644
--- a/src/meshTools/sets/cellSources/regionToCell/regionToCell.C
+++ b/src/meshTools/sets/cellSources/regionToCell/regionToCell.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2011, 2015-2018 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011, 2015-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2016 OpenFOAM Foundation
@@ -355,9 +355,9 @@ void Foam::regionToCell::combine(topoSet& set, const bool add) const
         cellSet subSet(mesh_, setName_);
 
         selectedCell = false;
-        forAllConstIter(cellSet, subSet, iter)
+        for (const label celli : subSet)
         {
-            selectedCell[iter.key()] = true;
+            selectedCell[celli] = true;
         }
     }
 
diff --git a/src/meshTools/sets/faceZoneSources/setsToFaceZone/setsToFaceZone.C b/src/meshTools/sets/faceZoneSources/setsToFaceZone/setsToFaceZone.C
index c354e2b541..0275cde080 100644
--- a/src/meshTools/sets/faceZoneSources/setsToFaceZone/setsToFaceZone.C
+++ b/src/meshTools/sets/faceZoneSources/setsToFaceZone/setsToFaceZone.C
@@ -136,10 +136,8 @@ void Foam::setsToFaceZone::applyToSet
             DynamicList<label> newAddressing(zoneSet.addressing());
             DynamicList<bool> newFlipMap(zoneSet.flipMap());
 
-            forAllConstIter(faceSet, fSet, iter)
+            for (const label facei : fSet)
             {
-                label facei = iter.key();
-
                 if (!zoneSet.found(facei))
                 {
                     bool flipFace = false;
diff --git a/src/meshTools/triSurface/booleanOps/intersectedSurface/intersectedSurface.C b/src/meshTools/triSurface/booleanOps/intersectedSurface/intersectedSurface.C
index 77aca92211..31b2042897 100644
--- a/src/meshTools/triSurface/booleanOps/intersectedSurface/intersectedSurface.C
+++ b/src/meshTools/triSurface/booleanOps/intersectedSurface/intersectedSurface.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010, 2015 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2010, 2015-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2016 OpenFOAM Foundation
@@ -60,16 +60,13 @@ void Foam::intersectedSurface::writeOBJ
     Ostream& os
 )
 {
-    forAll(points, pointi)
+    for (const point& pt : points)
     {
-        const point& pt = points[pointi];
-
         os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << nl;
     }
-    forAll(edges, edgei)
-    {
-        const edge& e = edges[edgei];
 
+    for (const edge& e : edges)
+    {
         os << "l " << e.start()+1 << ' ' << e.end()+1 << nl;
     }
 }
@@ -83,15 +80,13 @@ void Foam::intersectedSurface::writeOBJ
     Ostream& os
 )
 {
-    forAll(points, pointi)
+    for (const point& pt : points)
     {
-        const point& pt = points[pointi];
-
         os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << nl;
     }
-    forAll(faceEdges, i)
+    for (const label edgei : faceEdges)
     {
-        const edge& e = edges[faceEdges[i]];
+        const edge& e = edges[edgei];
 
         os << "l " << e.start()+1 << ' ' << e.end()+1 << nl;
     }
@@ -112,13 +107,13 @@ void Foam::intersectedSurface::writeLocalOBJ
 
     label maxVerti = 0;
 
-    forAll(faceEdges, i)
+    for (const label edgei : faceEdges)
     {
-        const edge& e = edges[faceEdges[i]];
+        const edge& e = edges[edgei];
 
         forAll(e, i)
         {
-            label pointi = e[i];
+            const label pointi = e[i];
 
             if (pointMap[pointi] == -1)
             {
@@ -131,11 +126,11 @@ void Foam::intersectedSurface::writeLocalOBJ
         }
     }
 
-    forAll(faceEdges, i)
+    for (const label edgei : faceEdges)
     {
-        const edge& e = edges[faceEdges[i]];
+        const edge& e = edges[edgei];
 
-        os << "l " << pointMap[e.start()]+1 << ' ' << pointMap[e.end()]+1
+        os  << "l " << pointMap[e.start()]+1 << ' ' << pointMap[e.end()]+1
             << nl;
     }
 }
@@ -148,18 +143,16 @@ void Foam::intersectedSurface::writeOBJ
     Ostream& os
 )
 {
-    forAll(points, pointi)
+    for (const point& pt : points)
     {
-        const point& pt = points[pointi];
-
         os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << nl;
     }
 
     os << 'f';
 
-    forAll(f, fp)
+    for (const label pointi : f)
     {
-        os << ' ' << f[fp]+1;
+        os << ' ' << pointi+1;
     }
     os << nl;
 }
@@ -173,10 +166,8 @@ void Foam::intersectedSurface::printVisit
 )
 {
     Pout<< "Visited:" << nl;
-    forAll(edgeLabels, i)
+    for (const label edgei : edgeLabels)
     {
-        label edgei = edgeLabels[i];
-
         const edge& e = edges[edgei];
 
         label stat = visited[edgei];
@@ -261,16 +252,7 @@ void Foam::intersectedSurface::incCount
     const label offset
 )
 {
-    Map<label>::iterator iter = visited.find(key);
-
-    if (iter == visited.end())
-    {
-        visited.insert(key, offset);
-    }
-    else
-    {
-        iter() += offset;
-    }
+    visited(key, 0) += offset;
 }
 
 
@@ -288,77 +270,50 @@ Foam::intersectedSurface::calcPointEdgeAddressing
 
     Map<DynamicList<label>> facePointEdges(4*fEdges.size());
 
-    forAll(fEdges, i)
+    for (const label edgei : fEdges)
     {
-        label edgei = fEdges[i];
-
         const edge& e = edges[edgei];
 
         // Add e.start to point-edges
-        Map<DynamicList<label>>::iterator iter =
-            facePointEdges.find(e.start());
-
-        if (iter == facePointEdges.end())
-        {
-            DynamicList<label> oneEdge;
-            oneEdge.append(edgei);
-            facePointEdges.insert(e.start(), oneEdge);
-        }
-        else
-        {
-            iter().append(edgei);
-        }
+        facePointEdges(e.start()).append(edgei);
 
         // Add e.end to point-edges
-        Map<DynamicList<label>>::iterator iter2 =
-            facePointEdges.find(e.end());
-
-        if (iter2 == facePointEdges.end())
-        {
-            DynamicList<label> oneEdge;
-            oneEdge.append(edgei);
-            facePointEdges.insert(e.end(), oneEdge);
-        }
-        else
-        {
-            iter2().append(edgei);
-        }
+        facePointEdges(e.end()).append(edgei);
     }
 
     // Shrink it
-    forAllIter(Map<DynamicList<label>>, facePointEdges, iter)
+    forAllIters(facePointEdges, iter)
     {
-        iter().shrink();
+        iter.val().shrink();
 
         // Check on dangling points.
-        if (iter().empty())
+        if (iter.val().empty())
         {
             FatalErrorInFunction
                 << "Point:" << iter.key() << " used by too few edges:"
-                << iter() << abort(FatalError);
+                << iter.val() << abort(FatalError);
         }
     }
 
     if (debug & 2)
     {
         // Print facePointEdges
-        Pout<< "calcPointEdgeAddressing: face consisting of edges:" << endl;
-        forAll(fEdges, i)
+        Pout<< "calcPointEdgeAddressing: face consisting of edges:" << nl;
+        for (const label edgei : fEdges)
         {
-            label edgei = fEdges[i];
             const edge& e = edges[edgei];
             Pout<< "    " << edgei << ' ' << e
                 << points[e.start()]
-                << points[e.end()] << endl;
+                << points[e.end()] << nl;
         }
 
-        Pout<< "    Constructed point-edge adressing:" << endl;
-        forAllConstIter(Map<DynamicList<label>>, facePointEdges, iter)
+        Pout<< "    Constructed point-edge adressing:" << nl;
+        forAllConstIters(facePointEdges, iter)
         {
             Pout<< "    vertex " << iter.key() << " is connected to edges "
-                << iter() << endl;
+                << iter.val() << nl;
         }
-        Pout<<endl;
+        Pout<< endl;
     }
 
     return facePointEdges;
@@ -469,10 +424,8 @@ Foam::label Foam::intersectedSurface::nextEdge
     scalar maxAngle = -GREAT;
     label maxEdgei = -1;
 
-    forAll(connectedEdges, conni)
+    for (const label edgei : connectedEdges)
     {
-        label edgei = connectedEdges[conni];
-
         if (edgei != prevEdgei)
         {
             label stat = visited[edgei];
@@ -638,14 +591,13 @@ void Foam::intersectedSurface::findNearestVisited
     minVerti = -1;
     minDist = GREAT;
 
-    forAllConstIter(Map<label>, pointVisited, iter)
+    forAllConstIters(pointVisited, iter)
     {
-        label pointi = iter.key();
+        const label pointi = iter.key();
+        const label nVisits = iter.val();
 
         if (pointi != excludePointi)
         {
-            label nVisits = iter();
-
             if (nVisits == 2*facePointEdges[pointi].size())
             {
                 // Fully visited (i.e. both sides of all edges)
@@ -689,11 +641,12 @@ Foam::faceList Foam::intersectedSurface::resplitFace
     // can compare number to facePointEdges.
     Map<label> pointVisited(2*facePointEdges.size());
 
-    forAllConstIter(Map<label>, visited, iter)
+    forAllConstIters(visited, iter)
     {
-        label edgei = iter.key();
+        const label edgei = iter.key();
+        const label stat = iter.val();
+
         const edge& e = eSurf.edges()[edgei];
-        label stat = iter();
 
         if (stat == STARTTOEND || stat == ENDTOSTART)
         {
@@ -715,13 +668,13 @@ Foam::faceList Foam::intersectedSurface::resplitFace
 
     if (debug)
     {
-        forAllConstIter(Map<label>, pointVisited, iter)
+        forAllConstIters(pointVisited, iter)
         {
-            label pointi = iter.key();
+            const label pointi = iter.key();
+            const label nVisits = iter.val();
 
-            label nVisits = iter();
-
-            Pout<< "point:" << pointi << "  nVisited:" << nVisits
+            Pout<< "point:" << pointi
+                << "  nVisited:" << nVisits
                 << "  pointEdges:" << facePointEdges[pointi].size() << endl;
         }
     }
@@ -736,13 +689,13 @@ Foam::faceList Foam::intersectedSurface::resplitFace
     {
         scalar minDist = GREAT;
 
-        forAllConstIter(Map<DynamicList<label>>, facePointEdges, iter)
+        forAllConstIters(facePointEdges, iter)
         {
-            label pointi = iter.key();
+            const label pointi = iter.key();
 
-            label nVisits = pointVisited[pointi];
+            const DynamicList<label>& pEdges = iter.val();
 
-            const DynamicList<label>& pEdges = iter();
+            const label nVisits = pointVisited[pointi];
 
             if (nVisits < 2*pEdges.size())
             {
@@ -779,15 +732,15 @@ Foam::faceList Foam::intersectedSurface::resplitFace
     {
         scalar minDist = GREAT;
 
-        forAllConstIter(Map<DynamicList<label>>, facePointEdges, iter)
+        forAllConstIters(facePointEdges, iter)
         {
-            label pointi = iter.key();
+            const label pointi = iter.key();
+
+            const DynamicList<label>& pEdges = iter.val();
 
             if (pointi != unvisitedVert0)
             {
-                label nVisits = pointVisited[pointi];
-
-                const DynamicList<label>& pEdges = iter();
+                const label nVisits = pointVisited[pointi];
 
                 if (nVisits < 2*pEdges.size())
                 {
@@ -1032,9 +985,9 @@ Foam::faceList Foam::intersectedSurface::splitFace
 
     if ((n & surf.faceNormals()[facei]) < 0)
     {
-        forAll(faces, i)
+        for (face& f : faces)
         {
-            reverse(faces[i]);
+            reverse(f);
         }
     }
 
@@ -1224,7 +1177,7 @@ Foam::intersectedSurface::intersectedSurface
 
             // Edge already connected to surfStart for sure. See if also
             // connects to surfEnd
-            if (surfE.start() == surfEndi || surfE.end() == surfEndi)
+            if (surfE.found(surfEndi))
             {
                 surfEdgei = pEdges[i];
 
diff --git a/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersectionFuncs.C b/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersectionFuncs.C
index f36e22134f..f25989f719 100644
--- a/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersectionFuncs.C
+++ b/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersectionFuncs.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2010, 2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2016 OpenFOAM Foundation
@@ -41,9 +41,8 @@ namespace Foam
 // Write points in obj format
 static void writeObjPoints(const UList<point>& pts, Ostream& os)
 {
-    forAll(pts, i)
+    for (const point& pt : pts)
     {
-        const point& pt = pts[i];
         os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << nl;
     }
 }
@@ -177,14 +176,9 @@ Foam::edgeList Foam::surfaceIntersection::filterEdges
     {
         const edge& e = edges[edgeI];
 
-        if
-        (
-            (e.start() != e.end())
-         && (uniqueEdges.find(e) == uniqueEdges.end())
-        )
+        if ((e.start() != e.end()) && uniqueEdges.insert(e))
         {
-            // Edge is -non degenerate and -not yet seen.
-            uniqueEdges.insert(e);
+            // Edge is non-degenerate and not yet seen.
 
             map[edgeI] = newEdgeI;
 
@@ -218,10 +212,9 @@ Foam::labelList Foam::surfaceIntersection::filterLabels
     {
         label elem = elems[elemI];
 
-        if (uniqueElems.find(elem) == uniqueElems.end())
+        if (uniqueElems.insert(elem))
         {
             // First time elem is seen
-            uniqueElems.insert(elem);
 
             map[elemI] = newElemI;
 
@@ -291,10 +284,8 @@ Foam::label Foam::surfaceIntersection::classify
     {
         return 1;
     }
-    else
-    {
-        return -1;
-    }
+
+    return -1;
 }
 
 
diff --git a/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C b/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C
index 27e7ea919a..aa7d459e53 100644
--- a/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C
+++ b/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2011, 2017 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011, 2017-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2016 OpenFOAM Foundation
@@ -803,11 +803,11 @@ Foam::surfaceFeatures::surfaceFeatures
 
     forAll(allEdgeStat, eI)
     {
-        EdgeMap<label>::const_iterator iter = dynFeatEdges.find(surfEdges[eI]);
+        const auto iter = dynFeatEdges.cfind(surfEdges[eI]);
 
-        if (iter != dynFeatEdges.end())
+        if (iter.found())
         {
-            allEdgeStat[eI] = edgeStat[iter()];
+            allEdgeStat[eI] = edgeStat[iter.val()];
         }
     }
 
@@ -1215,10 +1215,8 @@ void Foam::surfaceFeatures::writeObj(const fileName& prefix) const
     OFstream pointStr(prefix + "_points.obj");
     Pout<< "Writing feature points to " << pointStr.name() << endl;
 
-    forAll(featurePoints_, i)
+    for (const label pointi : featurePoints_)
     {
-        label pointi = featurePoints_[i];
-
         meshTools::writeOBJ(pointStr, surf_.localPoints()[pointi]);
     }
 }
@@ -1265,7 +1263,7 @@ Foam::Map<Foam::label> Foam::surfaceFeatures::nearestSamples
 
     forAll(pointLabels, i)
     {
-        label surfPointi = pointLabels[i];
+        const label surfPointi = pointLabels[i];
 
         const point& surfPt = surfPoints[surfPointi];
 
@@ -1304,10 +1302,10 @@ Foam::Map<Foam::label> Foam::surfaceFeatures::nearestSamples
         OFstream objStream("nearestSamples.obj");
 
         label vertI = 0;
-        forAllConstIter(Map<label>, nearest, iter)
+        forAllConstIters(nearest, iter)
         {
             meshTools::writeOBJ(objStream, samples[iter.key()]); vertI++;
-            meshTools::writeOBJ(objStream, surfPoints[iter()]); vertI++;
+            meshTools::writeOBJ(objStream, surfPoints[iter.val()]); vertI++;
             objStream<< "l " << vertI-1 << ' ' << vertI << endl;
         }
     }
@@ -1428,13 +1426,13 @@ Foam::Map<Foam::label> Foam::surfaceFeatures::nearestSamples
         OFstream objStream("nearestEdges.obj");
 
         label vertI = 0;
-        forAllConstIter(Map<label>, nearest, iter)
+        forAllConstIters(nearest, iter)
         {
             const label sampleI = iter.key();
 
-            meshTools::writeOBJ(objStream, samples[sampleI]); vertI++;
+            const edge& e = surfEdges[iter.val()];
 
-            const edge& e = surfEdges[iter()];
+            meshTools::writeOBJ(objStream, samples[sampleI]); vertI++;
 
             point nearPt =
                 e.line(surfPoints).nearestDist(samples[sampleI]).rawPoint();
@@ -1582,7 +1580,7 @@ Foam::Map<Foam::pointIndexHit> Foam::surfaceFeatures::nearestEdges
         OFstream objStream("nearestEdges.obj");
 
         label vertI = 0;
-        forAllConstIter(Map<pointIndexHit>, nearest, iter)
+        forAllConstIters(nearest, iter)
         {
             const label sampleEdgeI = iter.key();
 
@@ -1592,7 +1590,7 @@ Foam::Map<Foam::pointIndexHit> Foam::surfaceFeatures::nearestEdges
             meshTools::writeOBJ(objStream, sampleEdge.centre(samplePoints));
             vertI++;
 
-            meshTools::writeOBJ(objStream, iter().rawPoint());
+            meshTools::writeOBJ(objStream, iter.val().rawPoint());
             vertI++;
 
             objStream<< "l " << vertI-1 << ' ' << vertI << endl;
diff --git a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C
index 06811ab4a1..2dd7c980ae 100644
--- a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C
+++ b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010, 2016-2018 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2010, 2016-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2016 OpenFOAM Foundation
@@ -70,10 +70,8 @@ void Foam::triSurfaceTools::calcRefineStatus
 
         const labelList& myNeighbours = surf.faceFaces()[facei];
 
-        forAll(myNeighbours, myNeighbourI)
+        for (const label neighbourFacei : myNeighbours)
         {
-            label neighbourFacei = myNeighbours[myNeighbourI];
-
             if (refine[neighbourFacei] == GREEN)
             {
                 // Change to red refinement and propagate
@@ -165,13 +163,11 @@ Foam::triSurface Foam::triSurfaceTools::doRefine
 )
 {
     // Storage for new points. (start after old points)
-    DynamicList<point> newPoints(surf.nPoints());
-    forAll(surf.localPoints(), pointi)
-    {
-        newPoints.append(surf.localPoints()[pointi]);
-    }
     label newVertI = surf.nPoints();
 
+    DynamicList<point> newPoints(newVertI);
+    newPoints.append(surf.localPoints());
+
     // Storage for new faces
     DynamicList<labelledTri> newFaces(surf.size());
 
@@ -186,24 +182,15 @@ Foam::triSurface Foam::triSurfaceTools::doRefine
             // Create new vertices on all edges to be refined.
             const labelList& fEdges = surf.faceEdges()[facei];
 
-            forAll(fEdges, i)
+            for (const label edgei : fEdges)
             {
-                label edgeI = fEdges[i];
-
-                if (edgeMid[edgeI] == -1)
+                if (edgeMid[edgei] == -1)
                 {
-                    const edge& e = surf.edges()[edgeI];
+                    const edge& e = surf.edges()[edgei];
 
                     // Create new point on mid of edge
-                    newPoints.append
-                    (
-                        0.5
-                      * (
-                            surf.localPoints()[e.start()]
-                          + surf.localPoints()[e.end()]
-                        )
-                    );
-                    edgeMid[edgeI] = newVertI++;
+                    newPoints.append(e.centre(surf.localPoints()));
+                    edgeMid[edgei] = newVertI++;
                 }
             }
 
@@ -260,11 +247,9 @@ Foam::triSurface Foam::triSurfaceTools::doRefine
 
 
             // Create triangles for GREEN refinement.
-            forAll(fEdges, i)
+            for (const label edgei : fEdges)
             {
-                const label edgeI = fEdges[i];
-
-                label otherFacei = otherFace(surf, facei, edgeI);
+                label otherFacei = otherFace(surf, facei, edgei);
 
                 if ((otherFacei != -1) && (refineStatus[otherFacei] == GREEN))
                 {
@@ -272,8 +257,8 @@ Foam::triSurface Foam::triSurfaceTools::doRefine
                     (
                         surf,
                         otherFacei,
-                        edgeI,
-                        edgeMid[edgeI],
+                        edgei,
+                        edgeMid[edgei],
                         newFaces
                     );
                 }
@@ -335,26 +320,18 @@ void Foam::triSurfaceTools::protectNeighbours
     labelList& faceStatus
 )
 {
-//    const labelList& myFaces = surf.pointFaces()[vertI];
-//    forAll(myFaces, i)
+//    for (const label facei : surf.pointFaces()[vertI])
 //    {
-//        label facei = myFaces[i];
-//
 //        if ((faceStatus[facei] == ANYEDGE) || (faceStatus[facei] >= 0))
 //        {
 //            faceStatus[facei] = NOEDGE;
 //        }
 //    }
 
-    const labelList& myEdges = surf.pointEdges()[vertI];
-    forAll(myEdges, i)
+    for (const label edgei : surf.pointEdges()[vertI])
     {
-        const labelList& myFaces = surf.edgeFaces()[myEdges[i]];
-
-        forAll(myFaces, myFacei)
+        for (const label facei : surf.edgeFaces()[edgei])
         {
-            label facei = myFaces[myFacei];
-
             if ((faceStatus[facei] == ANYEDGE) || (faceStatus[facei] >= 0))
             {
                 faceStatus[facei] = NOEDGE;
@@ -376,8 +353,8 @@ Foam::labelHashSet Foam::triSurfaceTools::getCollapsedFaces
 )
 {
     const edge& e = surf.edges()[edgeI];
-    label v1 = e.start();
-    label v2 = e.end();
+    const label v1 = e.start();
+    const label v2 = e.end();
 
     // Faces using edge will certainly get collapsed.
     const labelList& myFaces = surf.edgeFaces()[edgeI];
@@ -391,10 +368,8 @@ Foam::labelHashSet Foam::triSurfaceTools::getCollapsedFaces
     //    collapses
     const labelList& v1Faces = surf.pointFaces()[v1];
 
-    forAll(v1Faces, v1Facei)
+    for (const label face1I : v1Faces)
     {
-        label face1I = v1Faces[v1Facei];
-
         label otherEdgeI = oppositeEdge(surf, face1I, v1);
 
         // Step across edge to other face
@@ -425,18 +400,13 @@ Foam::label Foam::triSurfaceTools::vertexUsesFace
     const label vertI
 )
 {
-    const labelList& myFaces = surf.pointFaces()[vertI];
-
-    forAll(myFaces, myFacei)
+    for (const label face1I : surf.pointFaces()[vertI])
     {
-        label face1I = myFaces[myFacei];
-
         if (faceUsed.found(face1I))
         {
             return face1I;
         }
     }
-
     return -1;
 }
 
@@ -452,8 +422,8 @@ void Foam::triSurfaceTools::getMergedEdges
 )
 {
     const edge& e = surf.edges()[edgeI];
-    label v1 = e.start();
-    label v2 = e.end();
+    const label v1 = e.start();
+    const label v2 = e.end();
 
     const labelList& v1Faces = surf.pointFaces()[v1];
     const labelList& v2Faces = surf.pointFaces()[v2];
@@ -461,19 +431,17 @@ void Foam::triSurfaceTools::getMergedEdges
     // Mark all (non collapsed) faces using v2
     labelHashSet v2FacesHash(v2Faces.size());
 
-    forAll(v2Faces, v2Facei)
+    for (const label facei : v2Faces)
     {
-        if (!collapsedFaces.found(v2Faces[v2Facei]))
+        if (!collapsedFaces.found(facei))
         {
-            v2FacesHash.insert(v2Faces[v2Facei]);
+            v2FacesHash.insert(facei);
         }
     }
 
 
-    forAll(v1Faces, v1Facei)
+    for (const label face1I: v1Faces)
     {
-        label face1I = v1Faces[v1Facei];
-
         if (collapsedFaces.found(face1I))
         {
             continue;
@@ -616,7 +584,6 @@ Foam::scalar Foam::triSurfaceTools::edgeCosAngle
                 << v1 << " of collapsed edge" << abort(FatalError);
         }
     }
-
     return cosAngle;
 }
 
@@ -635,21 +602,15 @@ Foam::scalar Foam::triSurfaceTools::collapseMinCosAngle
 
     scalar minCos = 1;
 
-    forAll(v1Faces, v1Facei)
+    for (const label facei : v1Faces)
     {
-        label facei = v1Faces[v1Facei];
-
         if (collapsedFaces.found(facei))
         {
             continue;
         }
 
-        const labelList& myEdges = surf.faceEdges()[facei];
-
-        forAll(myEdges, myEdgeI)
+        for (const label edgeI : surf.faceEdges()[facei])
         {
-            label edgeI = myEdges[myEdgeI];
-
             minCos =
                 min
                 (
@@ -688,10 +649,8 @@ bool Foam::triSurfaceTools::collapseCreatesFold
 {
     const labelList& v1Faces = surf.pointFaces()[v1];
 
-    forAll(v1Faces, v1Facei)
+    for (const label facei : v1Faces)
     {
-        label facei = v1Faces[v1Facei];
-
         if (collapsedFaces.found(facei))
         {
             continue;
@@ -699,10 +658,8 @@ bool Foam::triSurfaceTools::collapseCreatesFold
 
         const labelList& myEdges = surf.faceEdges()[facei];
 
-        forAll(myEdges, myEdgeI)
+        for (const label edgeI : myEdges)
         {
-            label edgeI = myEdges[myEdgeI];
-
             if
             (
                 edgeCosAngle
@@ -750,10 +707,8 @@ bool Foam::triSurfaceTools::collapseCreatesFold
 //
 //    labelList collapsed = collapsedFaces.toc();
 //
-//    forAll(collapsed, collapseI)
+//    for (const label facei : collapsed)
 //    {
-//        const label facei = collapsed[collapseI];
-//
 //        const labelList& myEdges = surf.faceEdges()[facei];
 //
 //        Pout<< "collapsing facei:" << facei << " uses edges:" << myEdges
@@ -1022,7 +977,6 @@ Foam::surfaceLocation Foam::triSurfaceTools::cutEdge
             }
         }
     }
-
     return cut;
 }
 
@@ -1173,10 +1127,8 @@ Foam::surfaceLocation Foam::triSurfaceTools::visitFaces
 
     scalar minDistSqr = Foam::sqr(GREAT);
 
-    forAll(eFaces, i)
+    for (const label triI : eFaces)
     {
-        label triI = eFaces[i];
-
         // Make sure we don't revisit previous face
         if (triI != start.triangle())
         {
@@ -1253,10 +1205,8 @@ void Foam::triSurfaceTools::writeOBJ
 {
     OFstream outFile(fName);
 
-    forAll(pts, pointi)
+    for (const point& pt : pts)
     {
-        const point& pt = pts[pointi];
-
         outFile<< "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
     }
     Pout<< "Written " << pts.size() << " vertices to file " << fName << endl;
@@ -1319,15 +1269,13 @@ void Foam::triSurfaceTools::getVertexTriangles
     edgeTris.setSize(startFaces.size() + endFaces.size() - myFaces.size());
 
     label nTris = 0;
-    forAll(startFaces, startFacei)
+    for (const label facei : startFaces)
     {
-        edgeTris[nTris++] = startFaces[startFacei];
+        edgeTris[nTris++] = facei;
     }
 
-    forAll(endFaces, endFacei)
+    for (const label facei : endFaces)
     {
-        label facei = endFaces[endFacei];
-
         if ((facei != face1I) && (facei != face2I))
         {
             edgeTris[nTris++] = facei;
@@ -1344,30 +1292,24 @@ Foam::labelList Foam::triSurfaceTools::getVertexVertices
 )
 {
     const edgeList& edges = surf.edges();
-
-    label v1 = e.start();
-    label v2 = e.end();
+    const label v1 = e.start();
+    const label v2 = e.end();
 
     // Get all vertices connected to v1 or v2 through an edge
     labelHashSet vertexNeighbours;
 
     const labelList& v1Edges = surf.pointEdges()[v1];
 
-    forAll(v1Edges, v1EdgeI)
+    for (const label edgei : v1Edges)
     {
-        const edge& e = edges[v1Edges[v1EdgeI]];
-        vertexNeighbours.insert(e.otherVertex(v1));
+        vertexNeighbours.insert(edges[edgei].otherVertex(v1));
     }
 
     const labelList& v2Edges = surf.pointEdges()[v2];
 
-    forAll(v2Edges, v2EdgeI)
+    for (const label edgei : v2Edges)
     {
-        const edge& e = edges[v2Edges[v2EdgeI]];
-
-        label vertI = e.otherVertex(v2);
-
-        vertexNeighbours.insert(vertI);
+        vertexNeighbours.insert(edges[edgei].otherVertex(v2));
     }
     return vertexNeighbours.toc();
 }
@@ -1387,16 +1329,13 @@ Foam::label Foam::triSurfaceTools::otherFace
     {
         return -1;
     }
+    else if (facei == myFaces[0])
+    {
+        return myFaces[1];
+    }
     else
     {
-        if (facei == myFaces[0])
-        {
-            return myFaces[1];
-        }
-        else
-        {
-            return myFaces[0];
-        }
+        return myFaces[0];
     }
 }
 
@@ -1460,7 +1399,8 @@ void Foam::triSurfaceTools::otherVertices
     else
     {
         FatalErrorInFunction
-            << "Vertex " << vertI << " not in face " << f << abort(FatalError);
+            << "Vertex " << vertI << " not in face " << f << nl
+            << abort(FatalError);
     }
 }
 
@@ -1475,20 +1415,19 @@ Foam::label Foam::triSurfaceTools::oppositeEdge
 {
     const labelList& myEdges = surf.faceEdges()[facei];
 
-    forAll(myEdges, myEdgeI)
+    for (const label edgei : myEdges)
     {
-        label edgeI = myEdges[myEdgeI];
+        const edge& e = surf.edges()[edgei];
 
-        const edge& e = surf.edges()[edgeI];
-
-        if ((e.start() != vertI) && (e.end() != vertI))
+        if (!e.found(vertI))
         {
-            return edgeI;
+            return edgei;
         }
     }
 
     FatalErrorInFunction
         << "Cannot find vertex " << vertI << " in edges of face " << facei
+        << nl
         << abort(FatalError);
 
     return -1;
@@ -1506,13 +1445,11 @@ Foam::label Foam::triSurfaceTools::oppositeVertex
     const triSurface::FaceType& f = surf.localFaces()[facei];
     const edge& e = surf.edges()[edgeI];
 
-    forAll(f, fp)
+    for (const label pointi : f)
     {
-        label vertI = f[fp];
-
-        if (vertI != e.start() && vertI != e.end())
+        if (!e.found(pointi))
         {
-            return vertI;
+            return pointi;
         }
     }
 
@@ -1534,14 +1471,13 @@ Foam::label Foam::triSurfaceTools::getEdge
 {
     const labelList& v1Edges = surf.pointEdges()[v1];
 
-    forAll(v1Edges, v1EdgeI)
+    for (const label edgei : v1Edges)
     {
-        label edgeI = v1Edges[v1EdgeI];
-        const edge& e = surf.edges()[edgeI];
+        const edge& e = surf.edges()[edgei];
 
-        if ((e.start() == v2) || (e.end() == v2))
+        if (e.found(v2))
         {
-            return edgeI;
+            return edgei;
         }
     }
     return -1;
@@ -1567,10 +1503,8 @@ Foam::label Foam::triSurfaceTools::getTriangle
 
     const labelList& eFaces = surf.edgeFaces()[e0I];
 
-    forAll(eFaces, eFacei)
+    for (const label facei : eFaces)
     {
-        label facei = eFaces[eFacei];
-
         const labelList& myEdges = surf.faceEdges()[facei];
 
         if
@@ -1607,13 +1541,7 @@ Foam::triSurface Foam::triSurfaceTools::collapseEdges
     forAll(edgeMids, edgeI)
     {
         const edge& e = surf.edges()[edgeI];
-
-        edgeMids[edgeI] =
-            0.5
-          * (
-                surf.localPoints()[e.start()]
-              + surf.localPoints()[e.end()]
-            );
+        edgeMids[edgeI] = e.centre(surf.localPoints());
     }
 
 
@@ -1674,20 +1602,18 @@ Foam::triSurface Foam::triSurfaceTools::collapseEdges
 
     // Do actual 'collapsing' of edges
 
-    forAll(collapseEdgeLabels, collapseEdgeI)
+    for (const label edgei : collapseEdgeLabels)
     {
-        const label edgeI = collapseEdgeLabels[collapseEdgeI];
-
-        if ((edgeI < 0) || (edgeI >= surf.nEdges()))
+        if (edgei < 0 || edgei >= surf.nEdges())
         {
             FatalErrorInFunction
                 << "Edge label outside valid range." << endl
-                << "edge label:" << edgeI << endl
+                << "edge label:" << edgei << endl
                 << "total number of edges:" << surf.nEdges() << endl
                 << abort(FatalError);
         }
 
-        const labelList& neighbours = edgeFaces[edgeI];
+        const labelList& neighbours = edgeFaces[edgei];
 
         if (neighbours.size() == 2)
         {
@@ -1697,11 +1623,11 @@ Foam::triSurface Foam::triSurfaceTools::collapseEdges
             // Check faceStatus to make sure this one can be collapsed
             if
             (
-                ((stat0 == ANYEDGE) || (stat0 == edgeI))
-             && ((stat1 == ANYEDGE) || (stat1 == edgeI))
+                ((stat0 == ANYEDGE) || (stat0 == edgei))
+             && ((stat1 == ANYEDGE) || (stat1 == edgei))
             )
             {
-                const edge& e = edges[edgeI];
+                const edge& e = edges[edgei];
 
                 // Set up mapping to 'collapse' points of edge
                 if
@@ -1712,7 +1638,7 @@ Foam::triSurface Foam::triSurfaceTools::collapseEdges
                 {
                     FatalErrorInFunction
                         << "points already mapped. Double collapse." << endl
-                        << "edgeI:" << edgeI
+                        << "edgei:" << edgei
                         << "  start:" << e.start()
                         << "  end:" << e.end()
                         << "  pointMap[start]:" << pointMap[e.start()]
@@ -1725,7 +1651,7 @@ Foam::triSurface Foam::triSurfaceTools::collapseEdges
                 pointMap[e.end()] = minVert;
 
                 // Move shared vertex to mid of edge
-                newPoints[minVert] = edgeMids[edgeI];
+                newPoints[minVert] = edgeMids[edgei];
 
                 // Protect neighbouring faces
                 protectNeighbours(surf, e.start(), faceStatus);
@@ -1733,13 +1659,13 @@ Foam::triSurface Foam::triSurfaceTools::collapseEdges
                 protectNeighbours
                 (
                     surf,
-                    oppositeVertex(surf, neighbours[0], edgeI),
+                    oppositeVertex(surf, neighbours[0], edgei),
                     faceStatus
                 );
                 protectNeighbours
                 (
                     surf,
-                    oppositeVertex(surf, neighbours[1], edgeI),
+                    oppositeVertex(surf, neighbours[1], edgei),
                     faceStatus
                 );
 
@@ -1748,7 +1674,7 @@ Foam::triSurface Foam::triSurfaceTools::collapseEdges
                     getCollapsedFaces
                     (
                         surf,
-                        edgeI
+                        edgei
                     ).toc();
 
                 forAll(collapseFaces, collapseI)
@@ -1848,17 +1774,15 @@ Foam::triSurface Foam::triSurfaceTools::greenRefine
 
 
     // Refine edges
-    forAll(refineEdges, refineEdgeI)
+    for (const label edgei : refineEdges)
     {
-        label edgeI = refineEdges[refineEdgeI];
-
-        const labelList& myFaces = surf.edgeFaces()[edgeI];
+        const labelList& myFaces = surf.edgeFaces()[edgei];
 
         bool neighbourIsRefined= false;
 
-        forAll(myFaces, myFacei)
+        for (const label facei : myFaces)
         {
-            if (refineStatus[myFaces[myFacei]] != NONE)
+            if (refineStatus[facei] != NONE)
             {
                 neighbourIsRefined =  true;
             }
@@ -1868,35 +1792,28 @@ Foam::triSurface Foam::triSurfaceTools::greenRefine
         if (!neighbourIsRefined)
         {
             // Refine edge
-            const edge& e = surf.edges()[edgeI];
+            const edge& e = surf.edges()[edgei];
 
-            point mid =
-                0.5
-              * (
-                    surf.localPoints()[e.start()]
-                  + surf.localPoints()[e.end()]
-                );
-
-            newPoints[newPointi] = mid;
+            newPoints[newPointi] = e.centre(surf.localPoints());
 
             // Refine faces using edge
-            forAll(myFaces, myFacei)
+            for (const label facei : myFaces)
             {
                 // Add faces to newFaces
                 greenRefine
                 (
                     surf,
-                    myFaces[myFacei],
-                    edgeI,
+                    facei,
+                    edgei,
                     newPointi,
                     newFaces
                 );
 
                 // Mark as refined
-                refineStatus[myFaces[myFacei]] = GREEN;
+                refineStatus[facei] = GREEN;
             }
 
-            newPointi++;
+            ++newPointi;
         }
     }
 
@@ -1928,26 +1845,23 @@ Foam::label Foam::triSurfaceTools::minEdge
     const labelList& edgeIndices
 )
 {
-    scalar minLength = GREAT;
-    label minIndex = -1;
-    forAll(edgeIndices, i)
+    scalar minLen = GREAT;
+    label minEdge = -1;
+
+    for (const label edgei : edgeIndices)
     {
-        const edge& e = surf.edges()[edgeIndices[i]];
+        const edge& e = surf.edges()[edgei];
 
-        scalar length =
-            mag
-            (
-                surf.localPoints()[e.end()]
-              - surf.localPoints()[e.start()]
-            );
+        const scalar length = e.mag(surf.localPoints());
 
-        if (length < minLength)
+        if (length < minLen)
         {
-            minLength = length;
-            minIndex = i;
+            minLen = length;
+            minEdge = edgei;
         }
     }
-    return edgeIndices[minIndex];
+
+    return minEdge;
 }
 
 
@@ -1958,26 +1872,23 @@ Foam::label Foam::triSurfaceTools::maxEdge
     const labelList& edgeIndices
 )
 {
-    scalar maxLength = -GREAT;
-    label maxIndex = -1;
-    forAll(edgeIndices, i)
+    scalar maxLen = -GREAT;
+    label maxEdge = -1;
+
+    for (const label edgei : edgeIndices)
     {
-        const edge& e = surf.edges()[edgeIndices[i]];
+        const edge& e = surf.edges()[edgei];
 
-        scalar length =
-            mag
-            (
-                surf.localPoints()[e.end()]
-              - surf.localPoints()[e.start()]
-            );
+        const scalar length = e.mag(surf.localPoints());
 
-        if (length > maxLength)
+        if (length > maxLen)
         {
-            maxLength = length;
-            maxIndex = i;
+            maxLen = length;
+            maxEdge = edgei;
         }
     }
-    return edgeIndices[maxIndex];
+
+    return maxEdge;
 }
 
 
@@ -2368,9 +2279,8 @@ Foam::triSurface Foam::triSurfaceTools::triangulate
 
     label newPatchi = 0;
 
-    forAllConstIter(labelHashSet, includePatches, iter)
+    for (const label patchi : includePatches)
     {
-        const label patchi = iter.key();
         const polyPatch& patch = bMesh[patchi];
         const pointField& points = patch.points();
 
@@ -2425,15 +2335,14 @@ Foam::triSurface Foam::triSurfaceTools::triangulate
 
     newPatchi = 0;
 
-    forAllConstIter(labelHashSet, includePatches, iter)
+    for (const label patchi : includePatches)
     {
-        const label patchi = iter.key();
         const polyPatch& patch = bMesh[patchi];
 
         surface.patches()[newPatchi].name() = patch.name();
         surface.patches()[newPatchi].geometricType() = patch.type();
 
-        newPatchi++;
+        ++newPatchi;
     }
 
     return surface;
@@ -2544,10 +2453,10 @@ Foam::triSurface Foam::triSurfaceTools::delaunay2D(const List<vector2D>& pts)
     // pts.begin() if double precision.
     List<doubleScalar> geompackVertices(2*pts.size());
     label doubleI = 0;
-    forAll(pts, i)
+    for (const vector2D& pt : pts)
     {
-        geompackVertices[doubleI++] = pts[i][0];
-        geompackVertices[doubleI++] = pts[i][1];
+        geompackVertices[doubleI++] = pt[0];
+        geompackVertices[doubleI++] = pt[1];
     }
 
     // Storage for triangles
@@ -2656,10 +2565,8 @@ void Foam::triSurfaceTools::calcInterpolationWeights
 
         scalar minDistance = GREAT;
 
-        forAll(s, facei)
+        for (const labelledTri& f : s)
         {
-            const labelledTri& f = s[facei];
-
             triPointRef tri(f.tri(points));
 
             label nearType, nearLabel;
@@ -2774,9 +2681,9 @@ bool Foam::triSurfaceTools::validTri
     const FaceType& f = surf[facei];
 
     // Simple check on indices ok.
-    forAll(f, fp)
+    for (const label pointi : f)
     {
-        if (f[fp] < 0 || f[fp] >= surf.points().size())
+        if (pointi < 0 || pointi >= surf.points().size())
         {
             if (verbose)
             {
@@ -2807,10 +2714,8 @@ bool Foam::triSurfaceTools::validTri
 
     // Check if faceNeighbours use same points as this face.
     // Note: discards normal information - sides of baffle are merged.
-    forAll(fFaces, i)
+    for (const label nbrFacei : fFaces)
     {
-        label nbrFacei = fFaces[i];
-
         if (nbrFacei <= facei)
         {
             // lower numbered faces already checked
@@ -2867,9 +2772,9 @@ bool Foam::triSurfaceTools::validTri
     }
 
     // Simple check on indices ok.
-    forAll(f, fp)
+    for (const label pointi : f)
     {
-        if (f[fp] < 0 || f[fp] >= surf.points().size())
+        if (pointi < 0 || pointi >= surf.points().size())
         {
             if (verbose)
             {
@@ -2900,10 +2805,8 @@ bool Foam::triSurfaceTools::validTri
 
     // Check if faceNeighbours use same points as this face.
     // Note: discards normal information - sides of baffle are merged.
-    forAll(fFaces, i)
+    for (const label nbrFacei : fFaces)
     {
-        label nbrFacei = fFaces[i];
-
         if (nbrFacei <= facei)
         {
             // lower numbered faces already checked
@@ -3058,7 +2961,6 @@ Foam::surfaceLocation Foam::triSurfaceTools::trackToEdge
         }
         snapToEnd(s, end, nearest);
     }
-
     return nearest;
 }
 
diff --git a/src/meshTools/twoDPointCorrector/twoDPointCorrector.C b/src/meshTools/twoDPointCorrector/twoDPointCorrector.C
index 53aac0ac96..5a13901b24 100644
--- a/src/meshTools/twoDPointCorrector/twoDPointCorrector.C
+++ b/src/meshTools/twoDPointCorrector/twoDPointCorrector.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011, 2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2016 OpenFOAM Foundation
@@ -59,14 +59,13 @@ void Foam::twoDPointCorrector::calcAddressing() const
     // Try and find a wedge patch
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
-    forAll(patches, patchi)
+    for (const polyPatch& p : patches)
     {
-        if (isA<wedgePolyPatch>(patches[patchi]))
+        if (isA<wedgePolyPatch>(p))
         {
             isWedge_ = true;
 
-            const wedgePolyPatch& wp =
-                refCast<const wedgePolyPatch>(patches[patchi]);
+            const wedgePolyPatch& wp = refCast<const wedgePolyPatch>(p);
 
             pn = wp.centreNormal();
 
@@ -75,7 +74,7 @@ void Foam::twoDPointCorrector::calcAddressing() const
 
             if (polyMesh::debug)
             {
-                Pout<< "Found normal from wedge patch " << patchi;
+                Pout<< "Found normal from wedge patch " << p.index() << nl;
             }
 
             break;
@@ -85,15 +84,15 @@ void Foam::twoDPointCorrector::calcAddressing() const
     // Try to find an empty patch with faces
     if (!isWedge_)
     {
-        forAll(patches, patchi)
+        for (const polyPatch& p : patches)
         {
-            if (isA<emptyPolyPatch>(patches[patchi]) && patches[patchi].size())
+            if (isA<emptyPolyPatch>(p) && p.size())
             {
-                pn = patches[patchi].faceAreas()[0];
+                pn = p.faceAreas()[0];
 
                 if (polyMesh::debug)
                 {
-                    Pout<< "Found normal from empty patch " << patchi;
+                    Pout<< "Found normal from empty patch " << p.index() << nl;
                 }
 
                 break;
@@ -115,7 +114,7 @@ void Foam::twoDPointCorrector::calcAddressing() const
 
     if (polyMesh::debug)
     {
-        Pout<< " twoDPointCorrector normal: " << pn << endl;
+        Pout<< " twoDPointCorrector normal: " << pn << nl;
     }
 
     // Select edges to be included in check.
@@ -132,7 +131,7 @@ void Foam::twoDPointCorrector::calcAddressing() const
     {
         const edge& e = meshEdges[edgeI];
 
-        vector edgeVector = e.vec(meshPoints)/(e.mag(meshPoints) + VSMALL);
+        vector edgeVector = e.unitVec(meshPoints);
 
         if (mag(edgeVector & pn) > edgeOrthogonalityTol)
         {
@@ -235,15 +234,13 @@ Foam::direction Foam::twoDPointCorrector::normalDir() const
     {
         return vector::Z;
     }
-    else
-    {
-        FatalErrorInFunction
-            << "Plane normal not aligned with the coordinate system" << nl
-            << "    pn = " << pn
-            << abort(FatalError);
 
-        return vector::Z;
-    }
+    FatalErrorInFunction
+        << "Plane normal not aligned with the coordinate system" << nl
+        << "    pn = " << pn
+        << abort(FatalError);
+
+    return vector::Z;
 }
 
 
@@ -284,11 +281,11 @@ void Foam::twoDPointCorrector::correctPoints(pointField& p) const
     const labelList& neIndices = normalEdgeIndices();
     const vector& pn = planeNormal();
 
-    forAll(neIndices, edgeI)
+    for (const label edgei : neIndices)
     {
-        point& pStart = p[meshEdges[neIndices[edgeI]].start()];
+        point& pStart = p[meshEdges[edgei].start()];
 
-        point& pEnd = p[meshEdges[neIndices[edgeI]].end()];
+        point& pEnd = p[meshEdges[edgei].end()];
 
         // calculate average point position
         point A = 0.5*(pStart + pEnd);
@@ -328,9 +325,9 @@ void Foam::twoDPointCorrector::correctDisplacement
     const labelList& neIndices = normalEdgeIndices();
     const vector& pn = planeNormal();
 
-    forAll(neIndices, edgeI)
+    for (const label edgei : neIndices)
     {
-        const edge& e = meshEdges[neIndices[edgeI]];
+        const edge& e = meshEdges[edgei];
 
         label startPointi = e.start();
         point pStart = p[startPointi] + disp[startPointi];