From c1964d7807b9995647989aad2101ad926d2e809a Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Fri, 10 Aug 2018 15:43:06 +0200
Subject: [PATCH] ENH: distinguish between face areaNormal/unitNormal in the
 code

---
 .../mesh/conversion/gmshToFoam/gmshToFoam.C   | 10 ++--
 .../netgenNeutralToFoam/netgenNeutralToFoam.C |  5 +-
 .../conformalVoronoiMeshCalcDualMesh.C        |  6 +--
 .../conformalVoronoiMeshZones.C               |  5 +-
 .../conformationSurfaces.C                    | 12 ++---
 .../manipulation/polyDualMesh/meshDualiser.C  | 12 ++---
 .../surface/surfaceInflate/surfaceInflate.C   | 19 ++++----
 src/OpenFOAM/meshes/polyMesh/polyMesh.C       |  4 +-
 .../constraint/cyclic/cyclicPolyPatch.C       | 33 ++++++-------
 .../constraint/oldCyclic/oldCyclicPolyPatch.C | 10 ++--
 .../PrimitivePatch/PrimitivePatchMeshData.C   |  9 ++--
 src/conversion/polyDualMesh/polyDualMesh.C    | 10 ++--
 src/dynamicMesh/meshCut/cellCuts/cellCuts.C   | 16 +++----
 .../polyMeshGeometry/polyMeshGeometry.C       |  8 ++--
 .../polyTopoChange/combineFaces.C             | 15 +++---
 .../polyTopoChange/hexRef8/hexRef8.C          | 20 ++++----
 .../enrichedPatch/enrichedPatchCutFaces.C     |  3 +-
 src/fileFormats/stl/STLtriangleI.H            |  5 +-
 .../vtk/read/vtkUnstructuredReader.C          |  4 +-
 .../faMesh/faMeshDemandDrivenData.C           |  4 +-
 .../faMesh/faPatches/faPatch/faPatch.C        |  7 ++-
 .../wallBoundedParticleTemplates.C            |  2 +-
 .../MaxwellianThermal/MaxwellianThermal.C     |  3 +-
 .../MixedDiffuseSpecular.C                    |  3 +-
 .../SpecularReflection/SpecularReflection.C   |  3 +-
 src/lagrangian/basic/particle/particle.H      |  3 +-
 src/lagrangian/basic/particle/particleI.H     |  8 ++--
 .../basic/particle/particleTemplates.C        |  3 +-
 .../ParticleCollector/ParticleCollector.C     |  7 +--
 .../PairCollision/PairCollision.C             | 17 +++----
 .../MPPIC/PackingModels/Implicit/Implicit.C   |  2 +-
 .../molecule/molecule/molecule.C              |  5 +-
 src/lagrangian/solidParticle/solidParticle.C  |  3 +-
 .../blockDescriptor/blockDescriptor.C         |  3 +-
 src/mesh/blockMesh/blockMesh/blockMeshCheck.C |  6 +--
 .../meshRefinementProblemCells.C              |  2 +-
 .../snappyHexMeshDriver/snappySnapDriver.C    |  9 ++--
 .../snappySnapDriverFeature.C                 |  4 +-
 .../AMIMethod/directAMI/directAMI.C           |  6 +--
 .../cyclicAMIPolyPatch/cyclicAMIPolyPatch.C   |  8 ++--
 .../indexedOctree/treeDataPrimitivePatch.C    |  2 +-
 .../primitiveMeshGeometry.C                   | 14 +++---
 .../triSurfaceMesh/triSurfaceMesh.C           | 13 ++----
 .../rotatedBoxToCell/rotatedBoxToCell.C       |  8 +---
 .../setAndNormalToFaceZone.C                  |  2 +-
 .../tetOverlapVolumeTemplates.C               |  4 +-
 .../intersectedSurface/intersectedSurface.C   |  2 +-
 .../faceTriangulation/faceTriangulation.C     | 32 +++++--------
 .../triSurfaceTools/triSurfaceCurvature.C     | 46 +++++++++----------
 .../distributedTriSurfaceMesh.C               |  3 +-
 .../surface/isoSurface/isoSurfaceCell.C       |  4 +-
 .../surfaceFormats/stl/STLsurfaceFormat.C     | 12 ++---
 52 files changed, 200 insertions(+), 256 deletions(-)

diff --git a/applications/utilities/mesh/conversion/gmshToFoam/gmshToFoam.C b/applications/utilities/mesh/conversion/gmshToFoam/gmshToFoam.C
index d0dd7c4a4e..51e9221577 100644
--- a/applications/utilities/mesh/conversion/gmshToFoam/gmshToFoam.C
+++ b/applications/utilities/mesh/conversion/gmshToFoam/gmshToFoam.C
@@ -188,19 +188,17 @@ label findInternalFace(const primitiveMesh& mesh, const labelList& meshF)
 bool correctOrientation(const pointField& points, const cellShape& shape)
 {
     // Get centre of shape.
-    point cc(shape.centre(points));
+    const point cc(shape.centre(points));
 
     // Get outwards pointing faces.
     faceList faces(shape.faces());
 
-    forAll(faces, i)
+    for (const face& f : faces)
     {
-        const face& f = faces[i];
-
-        vector n(f.normal(points));
+        const vector areaNorm(f.areaNormal(points));
 
         // Check if vector from any point on face to cc points outwards
-        if (((points[f[0]] - cc) & n) < 0)
+        if (((points[f[0]] - cc) & areaNorm) < 0)
         {
             // Incorrectly oriented
             return false;
diff --git a/applications/utilities/mesh/conversion/netgenNeutralToFoam/netgenNeutralToFoam.C b/applications/utilities/mesh/conversion/netgenNeutralToFoam/netgenNeutralToFoam.C
index 809f996afc..48b55bd14b 100644
--- a/applications/utilities/mesh/conversion/netgenNeutralToFoam/netgenNeutralToFoam.C
+++ b/applications/utilities/mesh/conversion/netgenNeutralToFoam/netgenNeutralToFoam.C
@@ -225,9 +225,10 @@ int main(int argc, char *argv[])
                 // Determine orientation of tri v.s. cell centre.
                 point cc(cll.centre(points));
                 point fc(tri.centre(points));
-                vector fn(tri.normal(points));
 
-                if (((fc - cc) & fn) < 0)
+                const vector areaNorm(tri.areaNormal(points));
+
+                if (((fc - cc) & areaNorm) < 0)
                 {
                     // Boundary face points inwards. Flip.
                     boundaryFaces[facei].flip();
diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
index fb9899438a..00bad93cd8 100644
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
@@ -1832,11 +1832,11 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
                         }
 
                         vector correctNormal = calcSharedPatchNormal(vc1, vc2);
-                        correctNormal /= mag(correctNormal);
+                        correctNormal.normalise();
 
                         Info<< "    cN " << correctNormal << endl;
 
-                        vector fN = f.normal(pts);
+                        vector fN = f.areaNormal(pts);
 
                         if (mag(fN) < SMALL)
                         {
@@ -1844,7 +1844,7 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
                             continue;
                         }
 
-                        fN /= mag(fN);
+                        fN.normalise();
                         Info<< "    fN " << fN << endl;
 
                         if ((fN & correctNormal) > 0)
diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshZones.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshZones.C
index 1f42b3e429..2a37049388 100644
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshZones.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshZones.C
@@ -481,10 +481,9 @@ void Foam::conformalVoronoiMesh::calcFaceZones
                     norm
                 );
 
-                vector fN = faces[facei].normal(mesh.points());
-                fN /= mag(fN) + SMALL;
+                const vector areaNorm = faces[facei].areaNormal(mesh.points());
 
-                if ((norm[0] & fN) < 0)
+                if ((norm[0] & areaNorm) < 0)
                 {
                     flipMap[facei] = true;
                 }
diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
index 0afaae03f3..37db9dab40 100644
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
@@ -89,11 +89,9 @@ void Foam::conformationSurfaces::hasBoundedVolume
             Info<< "        Index = " << surfaces_[s] << endl;
             Info<< "        Offset = " << regionOffset_[s] << endl;
 
-            forAll(triSurf, sI)
+            for (const labelledTri& f : triSurf)
             {
-                const label patchID =
-                    triSurf[sI].region()
-                  + regionOffset_[s];
+                const label patchID = f.region() + regionOffset_[s];
 
                 // Don't include baffle surfaces in the calculation
                 if
@@ -102,15 +100,15 @@ void Foam::conformationSurfaces::hasBoundedVolume
                  != extendedFeatureEdgeMesh::BOTH
                 )
                 {
-                    sum += triSurf[sI].normal(surfPts);
+                    sum += f.areaNormal(surfPts);
                 }
                 else
                 {
-                    nBaffles++;
+                    ++nBaffles;
                 }
             }
             Info<< "        has " << nBaffles << " baffles out of "
-                << triSurf.size() << " triangles" << endl;
+                << triSurf.size() << " triangles" << nl;
 
             totalTriangles += triSurf.size();
         }
diff --git a/applications/utilities/mesh/manipulation/polyDualMesh/meshDualiser.C b/applications/utilities/mesh/manipulation/polyDualMesh/meshDualiser.C
index cef8d2695b..ede9932cd3 100644
--- a/applications/utilities/mesh/manipulation/polyDualMesh/meshDualiser.C
+++ b/applications/utilities/mesh/manipulation/polyDualMesh/meshDualiser.C
@@ -271,8 +271,8 @@ Foam::label Foam::meshDualiser::addInternalFace
         );
 
         //pointField dualPoints(meshMod.points());
-        //vector n(newFace.normal(dualPoints));
-        //n /= mag(n);
+        //const vector n(newFace.unitNormal(dualPoints));
+        //
         //Pout<< "Generated internal dualFace:" << dualFacei
         //    << " verts:" << newFace
         //    << " points:" << UIndirectList<point>(meshMod.points(), newFace)
@@ -298,8 +298,8 @@ Foam::label Foam::meshDualiser::addInternalFace
         );
 
         //pointField dualPoints(meshMod.points());
-        //vector n(newFace.normal(dualPoints));
-        //n /= mag(n);
+        //const vector n(newFace.unitNormal(dualPoints));
+        //
         //Pout<< "Generated internal dualFace:" << dualFacei
         //    << " verts:" << newFace
         //    << " points:" << UIndirectList<point>(meshMod.points(), newFace)
@@ -355,8 +355,8 @@ Foam::label Foam::meshDualiser::addBoundaryFace
     );
 
     //pointField dualPoints(meshMod.points());
-    //vector n(newFace.normal(dualPoints));
-    //n /= mag(n);
+    //const vector n(newFace.unitNormal(dualPoints));
+    //
     //Pout<< "Generated boundary dualFace:" << dualFacei
     //    << " verts:" << newFace
     //    << " points:" << UIndirectList<point>(meshMod.points(), newFace)
diff --git a/applications/utilities/surface/surfaceInflate/surfaceInflate.C b/applications/utilities/surface/surfaceInflate/surfaceInflate.C
index 83c5d147a0..3d5e9fd906 100644
--- a/applications/utilities/surface/surfaceInflate/surfaceInflate.C
+++ b/applications/utilities/surface/surfaceInflate/surfaceInflate.C
@@ -89,11 +89,8 @@ tmp<vectorField> calcVertexNormals(const triSurface& surf)
 {
     // Weighted average of normals of faces attached to the vertex
     // Weight = fA / (mag(e0)^2 * mag(e1)^2);
-    tmp<vectorField> tpointNormals
-    (
-        new pointField(surf.nPoints(), Zero)
-    );
-    vectorField& pointNormals = tpointNormals.ref();
+    auto tpointNormals = tmp<vectorField>::New(surf.nPoints(), Zero);
+    auto& pointNormals = tpointNormals.ref();
 
     const pointField& points = surf.points();
     const labelListList& pointFaces = surf.pointFaces();
@@ -108,20 +105,20 @@ tmp<vectorField> calcVertexNormals(const triSurface& surf)
             const label faceI = pFaces[fI];
             const triFace& f = surf[faceI];
 
-            vector fN = f.normal(points);
+            vector areaNorm = f.areaNormal(points);
 
             scalar weight = calcVertexNormalWeight
             (
                 f,
                 meshPoints[pI],
-                fN,
+                areaNorm,
                 points
             );
 
-            pointNormals[pI] += weight*fN;
+            pointNormals[pI] += weight * areaNorm;
         }
 
-        pointNormals[pI] /= mag(pointNormals[pI]) + VSMALL;
+        pointNormals[pI].normalise();
     }
 
     return tpointNormals;
@@ -168,9 +165,9 @@ tmp<vectorField> calcPointNormals
 
                 // Get average edge normal
                 vector n = Zero;
-                forAll(eFaces, i)
+                for (const label facei : eFaces)
                 {
-                    n += s.faceNormals()[eFaces[i]];
+                    n += s.faceNormals()[facei];
                 }
                 n /= eFaces.size();
 
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMesh.C b/src/OpenFOAM/meshes/polyMesh/polyMesh.C
index f66b48555f..c25c9fb4a8 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMesh.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMesh.C
@@ -1375,7 +1375,7 @@ bool Foam::polyMesh::pointInCell
 
                     vector proj = p - faceTri.centre();
 
-                    if ((faceTri.normal() & proj) > 0)
+                    if ((faceTri.areaNormal() & proj) > 0)
                     {
                         return false;
                     }
@@ -1405,7 +1405,7 @@ bool Foam::polyMesh::pointInCell
 
                     vector proj = p - faceTri.centre();
 
-                    if ((faceTri.normal() & proj) > 0)
+                    if ((faceTri.areaNormal() & proj) > 0)
                     {
                         return false;
                     }
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C
index 4257571d72..845977c910 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C
@@ -60,9 +60,9 @@ Foam::label Foam::cyclicPolyPatch::findMaxArea
 
     forAll(faces, facei)
     {
-        scalar areaSqr = magSqr(faces[facei].normal(points));
+        scalar areaSqr = magSqr(faces[facei].areaNormal(points));
 
-        if (areaSqr > maxAreaSqr)
+        if (maxAreaSqr < areaSqr)
         {
             maxAreaSqr = areaSqr;
             maxI = facei;
@@ -81,7 +81,7 @@ void Foam::cyclicPolyPatch::calcTransforms()
         vectorField half0Areas(half0.size());
         forAll(half0, facei)
         {
-            half0Areas[facei] = half0[facei].normal(half0.points());
+            half0Areas[facei] = half0[facei].areaNormal(half0.points());
         }
 
         // Half1
@@ -89,7 +89,7 @@ void Foam::cyclicPolyPatch::calcTransforms()
         vectorField half1Areas(half1.size());
         forAll(half1, facei)
         {
-            half1Areas[facei] = half1[facei].normal(half1.points());
+            half1Areas[facei] = half1[facei].areaNormal(half1.points());
         }
 
         calcTransforms
@@ -260,19 +260,17 @@ void Foam::cyclicPolyPatch::calcTransforms
             // use calculated normals.
             vector n0 = findFaceMaxRadius(half0Ctrs);
             vector n1 = -findFaceMaxRadius(half1Ctrs);
-            n0 /= mag(n0) + VSMALL;
-            n1 /= mag(n1) + VSMALL;
+            n0.normalise();
+            n1.normalise();
 
             if (debug)
             {
-                scalar theta = radToDeg(acos(n0 & n1));
-
                 Pout<< "cyclicPolyPatch::calcTransforms :"
                     << " patch:" << name()
                     << " Specified rotation :"
                     << " n0:" << n0 << " n1:" << n1
-                    << " swept angle: " << theta << " [deg]"
-                    << endl;
+                    << " swept angle: " << radToDeg(acos(n0 & n1))
+                    << " [deg]" << endl;
             }
 
             // Extended tensor from two local coordinate systems calculated
@@ -420,18 +418,17 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors
             {
                 vector n0 = findFaceMaxRadius(half0Ctrs);
                 vector n1 = -findFaceMaxRadius(half1Ctrs);
-                n0 /= mag(n0) + VSMALL;
-                n1 /= mag(n1) + VSMALL;
+                n0.normalise();
+                n1.normalise();
 
                 if (debug)
                 {
-                    scalar theta = radToDeg(acos(n0 & n1));
-
                     Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
                         << " patch:" << name()
                         << " Specified rotation :"
                         << " n0:" << n0 << " n1:" << n1
-                        << " swept angle: " << theta << " [deg]"
+                        << " swept angle: "
+                        << radToDeg(acos(n0 & n1)) << " [deg]"
                         << endl;
                 }
 
@@ -499,12 +496,10 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors
                 // Determine the face with max area on both halves. These
                 // two faces are used to determine the transformation tensors
                 label max0I = findMaxArea(pp0.points(), pp0);
-                vector n0 = pp0[max0I].normal(pp0.points());
-                n0 /= mag(n0) + VSMALL;
+                const vector n0 = pp0[max0I].unitNormal(pp0.points());
 
                 label max1I = findMaxArea(pp1.points(), pp1);
-                vector n1 = pp1[max1I].normal(pp1.points());
-                n1 /= mag(n1) + VSMALL;
+                const vector n1 = pp1[max1I].unitNormal(pp1.points());
 
                 if (mag(n0 & n1) < 1-matchTolerance())
                 {
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/oldCyclic/oldCyclicPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/oldCyclic/oldCyclicPolyPatch.C
index 80714b0bb4..d04552e934 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/oldCyclic/oldCyclicPolyPatch.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/oldCyclic/oldCyclicPolyPatch.C
@@ -92,9 +92,9 @@ Foam::label Foam::oldCyclicPolyPatch::findMaxArea
 
     forAll(faces, facei)
     {
-        scalar areaSqr = magSqr(faces[facei].normal(points));
+        scalar areaSqr = magSqr(faces[facei].areaNormal(points));
 
-        if (areaSqr > maxAreaSqr)
+        if (maxAreaSqr < areaSqr)
         {
             maxAreaSqr = areaSqr;
             maxI = facei;
@@ -359,12 +359,10 @@ void Foam::oldCyclicPolyPatch::getCentresAndAnchors
             // Determine the face with max area on both halves. These
             // two faces are used to determine the transformation tensors
             label max0I = findMaxArea(pp.points(), half0Faces);
-            vector n0 = half0Faces[max0I].normal(pp.points());
-            n0 /= mag(n0) + VSMALL;
+            const vector n0 = half0Faces[max0I].unitNormal(pp.points());
 
             label max1I = findMaxArea(pp.points(), half1Faces);
-            vector n1 = half1Faces[max1I].normal(pp.points());
-            n1 /= mag(n1) + VSMALL;
+            const vector n1 = half1Faces[max1I].unitNormal(pp.points());
 
             if (mag(n0 & n1) < 1-matchTolerance())
             {
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchMeshData.C b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchMeshData.C
index 50e4222be0..a513934b28 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchMeshData.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchMeshData.C
@@ -281,9 +281,9 @@ calcPointNormals() const
 
         const labelList& curFaces = pf[pointi];
 
-        forAll(curFaces, facei)
+        for (const label facei : curFaces)
         {
-            curNormal += faceUnitNormals[curFaces[facei]];
+            curNormal += faceUnitNormals[facei];
         }
 
         curNormal /= mag(curNormal) + VSMALL;
@@ -425,7 +425,7 @@ calcFaceAreas() const
 
     forAll(n, facei)
     {
-        n[facei] = this->operator[](facei).normal(points_);
+        n[facei] = this->operator[](facei).areaNormal(points_);
     }
 
     if (debug)
@@ -472,8 +472,7 @@ calcFaceNormals() const
 
     forAll(n, facei)
     {
-        n[facei] = this->operator[](facei).normal(points_);
-        n[facei] /= mag(n[facei]) + VSMALL;
+        n[facei] = this->operator[](facei).unitNormal(points_);
     }
 
     if (debug)
diff --git a/src/conversion/polyDualMesh/polyDualMesh.C b/src/conversion/polyDualMesh/polyDualMesh.C
index b0c1f2eaea..9b325c9415 100644
--- a/src/conversion/polyDualMesh/polyDualMesh.C
+++ b/src/conversion/polyDualMesh/polyDualMesh.C
@@ -1006,8 +1006,9 @@ void Foam::polyDualMesh::calcDual
         {
             // Check orientation.
             const face& f = dynDualFaces.last();
-            vector n = f.normal(dualPoints);
-            if (((mesh.points()[owner] - dualPoints[f[0]]) & n) > 0)
+            const vector areaNorm = f.areaNormal(dualPoints);
+
+            if (((mesh.points()[owner] - dualPoints[f[0]]) & areaNorm) > 0)
             {
                 WarningInFunction
                     << " on boundary edge:" << edgeI
@@ -1121,8 +1122,9 @@ void Foam::polyDualMesh::calcDual
             {
                 // Check orientation.
                 const face& f = dynDualFaces.last();
-                vector n = f.normal(dualPoints);
-                if (((mesh.points()[owner] - dualPoints[f[0]]) & n) > 0)
+                const vector areaNorm = f.areaNormal(dualPoints);
+
+                if (((mesh.points()[owner] - dualPoints[f[0]]) & areaNorm) > 0)
                 {
                     WarningInFunction
                         << " on internal edge:" << edgeI
diff --git a/src/dynamicMesh/meshCut/cellCuts/cellCuts.C b/src/dynamicMesh/meshCut/cellCuts/cellCuts.C
index 1dbbc324ef..edb8209d6f 100644
--- a/src/dynamicMesh/meshCut/cellCuts/cellCuts.C
+++ b/src/dynamicMesh/meshCut/cellCuts/cellCuts.C
@@ -1363,28 +1363,26 @@ bool Foam::cellCuts::loopAnchorConsistent
     // Create identity face for ease of calculation of normal etc.
     face f(identity(loopPts.size()));
 
-    vector normal = f.normal(loopPts);
-    point ctr = f.centre(loopPts);
+    const vector areaNorm = f.areaNormal(loopPts);
+    const point ctr = f.centre(loopPts);
 
 
     // Get average position of anchor points.
     vector avg(Zero);
 
-    forAll(anchorPoints, ptI)
+    for (const label pointi : anchorPoints)
     {
-        avg += mesh().points()[anchorPoints[ptI]];
+        avg += mesh().points()[pointi];
     }
     avg /= anchorPoints.size();
 
 
-    if (((avg - ctr) & normal) > 0)
+    if (((avg - ctr) & areaNorm) > 0)
     {
         return true;
     }
-    else
-    {
-        return false;
-    }
+
+    return false;
 }
 
 
diff --git a/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C b/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C
index 39807dc780..c5419bfc10 100644
--- a/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C
+++ b/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C
@@ -1646,7 +1646,7 @@ bool Foam::polyMeshGeometry::checkFaceTwist
 //                        p[f[fpI]],
 //                        p[f.nextLabel(fpI)],
 //                        fc
-//                    ).normal()
+//                    ).areaNormal()
 //                );
 //
 //                scalar magTri = mag(triArea);
@@ -1721,7 +1721,7 @@ bool Foam::polyMeshGeometry::checkFaceTwist
                             p[f[fpI]],
                             p[f.nextLabel(fpI)],
                             fc
-                        ).normal()
+                        ).areaNormal()
                     );
 
                     scalar magTri = mag(triArea);
@@ -1826,7 +1826,7 @@ bool Foam::polyMeshGeometry::checkTriangleTwist
                     p[f[fp]],
                     p[f.nextLabel(fp)],
                     fc
-                ).normal();
+                ).areaNormal();
 
                 scalar magTri = mag(prevN);
 
@@ -1853,7 +1853,7 @@ bool Foam::polyMeshGeometry::checkTriangleTwist
                             p[f[fp]],
                             p[f.nextLabel(fp)],
                             fc
-                        ).normal()
+                        ).areaNormal()
                     );
                     scalar magTri = mag(triN);
 
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.C
index 426866008b..44c3ca0afc 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.C
@@ -55,9 +55,8 @@ bool Foam::combineFaces::convexFace
     const face& f
 )
 {
-    // Get outwards pointing normal of f.
-    vector n = f.normal(points);
-    n /= mag(n);
+    // Get outwards pointing normal of f, only the sign matters.
+    const vector areaNorm = f.areaNormal(points);
 
     // Get edge from f[0] to f[size-1];
     vector ePrev(points[f.first()] - points[f.last()]);
@@ -78,7 +77,7 @@ bool Foam::combineFaces::convexFace
         {
             vector edgeNormal = ePrev ^ e10;
 
-            if ((edgeNormal & n) < 0)
+            if ((edgeNormal & areaNorm) < 0)
             {
                 // Concave. Check angle.
                 if ((ePrev & e10) < minConcaveCos)
@@ -150,12 +149,10 @@ void Foam::combineFaces::regioniseFaces
         // - small angle
         if (p0 != -1 && p0 == p1 && !patches[p0].coupled())
         {
-            vector f0Normal = mesh_.faceAreas()[f0];
-            f0Normal /= mag(f0Normal);
-            vector f1Normal = mesh_.faceAreas()[f1];
-            f1Normal /= mag(f1Normal);
+            vector f0Normal = normalised(mesh_.faceAreas()[f0]);
+            vector f1Normal = normalised(mesh_.faceAreas()[f1]);
 
-            if ((f0Normal&f1Normal) > minCos)
+            if ((f0Normal & f1Normal) > minCos)
             {
                 Map<label>::const_iterator f0Fnd = faceRegion.find(f0);
 
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.C
index e58aa2ae8e..8d1b96c318 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.C
@@ -835,11 +835,11 @@ void Foam::hexRef8::checkInternalOrientation
     face compactFace(identity(newFace.size()));
     pointField compactPoints(meshMod.points(), newFace);
 
-    vector n(compactFace.normal(compactPoints));
+    const vector areaNorm(compactFace.areaNormal(compactPoints));
 
-    vector dir(neiPt - ownPt);
+    const vector dir(neiPt - ownPt);
 
-    if ((dir & n) < 0)
+    if ((dir & areaNorm) < 0)
     {
         FatalErrorInFunction
             << "cell:" << celli << " old face:" << facei
@@ -850,9 +850,9 @@ void Foam::hexRef8::checkInternalOrientation
             << abort(FatalError);
     }
 
-    vector fcToOwn(compactFace.centre(compactPoints) - ownPt);
+    const vector fcToOwn(compactFace.centre(compactPoints) - ownPt);
 
-    scalar s = (fcToOwn&n) / (dir&n);
+    const scalar s = (fcToOwn & areaNorm) / (dir & areaNorm);
 
     if (s < 0.1 || s > 0.9)
     {
@@ -881,11 +881,11 @@ void Foam::hexRef8::checkBoundaryOrientation
     face compactFace(identity(newFace.size()));
     pointField compactPoints(meshMod.points(), newFace);
 
-    vector n(compactFace.normal(compactPoints));
+    const vector areaNorm(compactFace.areaNormal(compactPoints));
 
-    vector dir(boundaryPt - ownPt);
+    const vector dir(boundaryPt - ownPt);
 
-    if ((dir & n) < 0)
+    if ((dir & areaNorm) < 0)
     {
         FatalErrorInFunction
             << "cell:" << celli << " old face:" << facei
@@ -896,9 +896,9 @@ void Foam::hexRef8::checkBoundaryOrientation
             << abort(FatalError);
     }
 
-    vector fcToOwn(compactFace.centre(compactPoints) - ownPt);
+    const vector fcToOwn(compactFace.centre(compactPoints) - ownPt);
 
-    scalar s = (fcToOwn&dir) / magSqr(dir);
+    const scalar s = (fcToOwn & dir) / magSqr(dir);
 
     if (s < 0.7 || s > 1.3)
     {
diff --git a/src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatchCutFaces.C b/src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatchCutFaces.C
index 4484c7f852..4887943021 100644
--- a/src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatchCutFaces.C
+++ b/src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatchCutFaces.C
@@ -150,8 +150,7 @@ void Foam::enrichedPatch::calcCutFaces() const
         }
 
         // Grab face normal
-        vector normal = curLocalFace.normal(lp);
-        normal /= mag(normal);
+        const vector normal = curLocalFace.unitNormal(lp);
 
         while (edgeSeeds.size())
         {
diff --git a/src/fileFormats/stl/STLtriangleI.H b/src/fileFormats/stl/STLtriangleI.H
index 541bec22bf..70529371b4 100644
--- a/src/fileFormats/stl/STLtriangleI.H
+++ b/src/fileFormats/stl/STLtriangleI.H
@@ -143,9 +143,8 @@ inline void Foam::STLtriangle::write
     const point& pt2
 )
 {
-    // calculate the normal ourselves
-    vector norm = triPointRef(pt0, pt1, pt2).normal();
-    norm /= mag(norm) + VSMALL;
+    // Calculate the normal ourselves
+    const vector norm = triPointRef(pt0, pt1, pt2).unitNormal();
 
     write(os, norm, pt0, pt1, pt2);
 }
diff --git a/src/fileFormats/vtk/read/vtkUnstructuredReader.C b/src/fileFormats/vtk/read/vtkUnstructuredReader.C
index c20f627667..3094f203e3 100644
--- a/src/fileFormats/vtk/read/vtkUnstructuredReader.C
+++ b/src/fileFormats/vtk/read/vtkUnstructuredReader.C
@@ -979,10 +979,10 @@ void Foam::vtkUnstructuredReader::read(ISstream& inFile)
                 );
 
                 const point bottomCc(bottom.centre());
-                const vector bottomNormal(bottom.normal());
+                const vector bottomNormal(bottom.areaNormal());
                 const point topCc(top.centre());
 
-                if (((topCc-bottomCc)&bottomNormal) < 0)
+                if (((topCc - bottomCc) & bottomNormal) < 0)
                 {
                     // Flip top and bottom
                     Swap(shape[0], shape[3]);
diff --git a/src/finiteArea/faMesh/faMeshDemandDrivenData.C b/src/finiteArea/faMesh/faMeshDemandDrivenData.C
index 2555274450..e9cc68d9cc 100644
--- a/src/finiteArea/faMesh/faMeshDemandDrivenData.C
+++ b/src/finiteArea/faMesh/faMeshDemandDrivenData.C
@@ -425,9 +425,7 @@ void Foam::faMesh::calcFaceAreaNormals() const
     vectorField& nInternal = faceAreaNormals.ref();
     forAll(localFaces, faceI)
     {
-        nInternal[faceI] =
-            localFaces[faceI].normal(localPoints)/
-            localFaces[faceI].mag(localPoints);
+        nInternal[faceI] = localFaces[faceI].unitNormal(localPoints);
     }
 
     forAll(boundary(), patchI)
diff --git a/src/finiteArea/faMesh/faPatches/faPatch/faPatch.C b/src/finiteArea/faMesh/faPatches/faPatch/faPatch.C
index b0c96720a1..8c64de4c8a 100644
--- a/src/finiteArea/faMesh/faPatches/faPatch/faPatch.C
+++ b/src/finiteArea/faMesh/faPatches/faPatch/faPatch.C
@@ -306,8 +306,8 @@ Foam::labelList Foam::faPatch::ngbPolyPatchFaces() const
 
 Foam::tmp<Foam::vectorField> Foam::faPatch::ngbPolyPatchFaceNormals() const
 {
-    tmp<vectorField> tfN(new vectorField());
-    vectorField& fN = tfN.ref();
+    auto tfN = tmp<vectorField>::New();
+    auto& fN = tfN.ref();
 
     if (ngbPolyPatchIndex() == -1)
     {
@@ -325,8 +325,7 @@ Foam::tmp<Foam::vectorField> Foam::faPatch::ngbPolyPatchFaceNormals() const
 
     forAll(fN, faceI)
     {
-        fN[faceI] = faces[ngbFaces[faceI]].normal(points)
-            /faces[ngbFaces[faceI]].mag(points);
+        fN[faceI] = faces[ngbFaces[faceI]].unitNormal(points);
     }
 
     return tfN;
diff --git a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticleTemplates.C b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticleTemplates.C
index 41b6eb06bf..ed0c304a21 100644
--- a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticleTemplates.C
+++ b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticleTemplates.C
@@ -127,7 +127,7 @@ Foam::scalar Foam::wallBoundedParticle::trackToEdge
             const bool posVol = (nbrTi.tet(mesh()).mag() > 0);
             const vector path(endPosition - localPosition_);
 
-            if (posVol == ((nbrTi.faceTri(mesh()).normal() & path) < 0))
+            if (posVol == ((nbrTi.faceTri(mesh()).areaNormal() & path) < 0))
             {
                 // Change into nbrCell. No need to change tetFace, tetPt.
                 //Pout<< "    crossed from cell:" << celli_
diff --git a/src/lagrangian/DSMC/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C b/src/lagrangian/DSMC/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C
index dd145c62e1..ca2bc10041 100644
--- a/src/lagrangian/DSMC/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C
+++ b/src/lagrangian/DSMC/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C
@@ -68,8 +68,7 @@ void Foam::MaxwellianThermal<CloudType>::correct
 
     label wppLocalFace = wpp.whichFace(p.face());
 
-    vector nw = p.normal();
-    nw /= mag(nw);
+    const vector nw = p.normal();
 
     // Normal velocity magnitude
     scalar U_dot_nw = U & nw;
diff --git a/src/lagrangian/DSMC/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C b/src/lagrangian/DSMC/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C
index 21c2cd73a8..4fe4f0649d 100644
--- a/src/lagrangian/DSMC/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C
+++ b/src/lagrangian/DSMC/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C
@@ -66,8 +66,7 @@ void Foam::MixedDiffuseSpecular<CloudType>::correct
 
     label wppLocalFace = wpp.whichFace(p.face());
 
-    vector nw = p.normal();
-    nw /= mag(nw);
+    const vector nw = p.normal();
 
     // Normal velocity magnitude
     scalar U_dot_nw = U & nw;
diff --git a/src/lagrangian/DSMC/submodels/WallInteractionModel/SpecularReflection/SpecularReflection.C b/src/lagrangian/DSMC/submodels/WallInteractionModel/SpecularReflection/SpecularReflection.C
index a3d182d07b..2d9a00df53 100644
--- a/src/lagrangian/DSMC/submodels/WallInteractionModel/SpecularReflection/SpecularReflection.C
+++ b/src/lagrangian/DSMC/submodels/WallInteractionModel/SpecularReflection/SpecularReflection.C
@@ -57,8 +57,7 @@ void Foam::SpecularReflection<CloudType>::correct
 {
     vector& U = p.U();
 
-    vector nw = p.normal();
-    nw /= mag(nw);
+    const vector nw = p.normal();
 
     scalar U_dot_nw = U & nw;
 
diff --git a/src/lagrangian/basic/particle/particle.H b/src/lagrangian/basic/particle/particle.H
index 4d6928e2b1..f7e9f44532 100644
--- a/src/lagrangian/basic/particle/particle.H
+++ b/src/lagrangian/basic/particle/particle.H
@@ -522,8 +522,7 @@ public:
             //- Return the current tet transformation tensor
             inline barycentricTensor currentTetTransform() const;
 
-            //- Return the normal of the tri on tetFacei_ for the
-            //  current tet.
+            //- The (unit) normal of the tri on tetFacei_ for the current tet.
             inline vector normal() const;
 
             //- Is the particle on a face?
diff --git a/src/lagrangian/basic/particle/particleI.H b/src/lagrangian/basic/particle/particleI.H
index 771eeec057..3698b17f22 100644
--- a/src/lagrangian/basic/particle/particleI.H
+++ b/src/lagrangian/basic/particle/particleI.H
@@ -276,7 +276,7 @@ inline Foam::barycentricTensor Foam::particle::currentTetTransform() const
 
 inline Foam::vector Foam::particle::normal() const
 {
-    return currentTetIndices().faceTri(mesh_).normal();
+    return currentTetIndices().faceTri(mesh_).unitNormal();
 }
 
 
@@ -324,8 +324,7 @@ void Foam::particle::patchData(vector& n, vector& U) const
         Pair<vector> centre, base, vertex1, vertex2;
         movingTetGeometry(1, centre, base, vertex1, vertex2);
 
-        n = triPointRef(base[0], vertex1[0], vertex2[0]).normal();
-        n /= mag(n);
+        n = triPointRef(base[0], vertex1[0], vertex2[0]).unitNormal();
 
         // Interpolate the motion of the three face vertices to the current
         // coordinates
@@ -343,8 +342,7 @@ void Foam::particle::patchData(vector& n, vector& U) const
         vector centre, base, vertex1, vertex2;
         stationaryTetGeometry(centre, base, vertex1, vertex2);
 
-        n = triPointRef(base, vertex1, vertex2).normal();
-        n /= mag(n);
+        n = triPointRef(base, vertex1, vertex2).unitNormal();
 
         U = Zero;
     }
diff --git a/src/lagrangian/basic/particle/particleTemplates.C b/src/lagrangian/basic/particle/particleTemplates.C
index a43c776397..fffeacd3f0 100644
--- a/src/lagrangian/basic/particle/particleTemplates.C
+++ b/src/lagrangian/basic/particle/particleTemplates.C
@@ -250,8 +250,7 @@ void Foam::particle::hitSymmetryPlanePatch
 template<class TrackCloudType>
 void Foam::particle::hitSymmetryPatch(TrackCloudType&, trackingData&)
 {
-    vector nf = normal();
-    nf /= mag(nf);
+    const vector nf = normal();
 
     transformProperties(I - 2.0*nf*nf);
 }
diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.C b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.C
index 1a8255ae7d..0adce43d20 100644
--- a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.C
+++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.C
@@ -315,13 +315,14 @@ void Foam::ParticleCollector<CloudType>::collectParcelPolygon
         // the face's decomposed triangles does not work due to ambiguity along
         // the diagonals.
         const face& f = faces_[facei];
-        const vector n = f.normal(points_);
+        const vector areaNorm = f.areaNormal(points_);
         bool inside = true;
-        for (label i = 0; i < f.size(); ++ i)
+        for (label i = 0; i < f.size(); ++i)
         {
             const label j = f.fcIndex(i);
             const triPointRef t(pIntersect, points_[f[i]], points_[f[j]]);
-            if ((n & t.normal()) < 0)
+
+            if ((areaNorm & t.areaNormal()) < 0)
             {
                 inside = false;
                 break;
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairCollision.C b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairCollision.C
index 7d68828701..9640726337 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairCollision.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairCollision.C
@@ -238,15 +238,14 @@ void Foam::PairCollision<CloudType>::wallInteraction()
 
                 if (nearest.distance() < r)
                 {
-                    vector normal = mesh.faceAreas()[realFacei];
-
-                    normal /= mag(normal);
+                    const vector normal =
+                        normalised(mesh.faceAreas()[realFacei]);
 
                     const vector& nearPt = nearest.rawPoint();
 
-                    vector pW = nearPt - pos;
+                    const vector pW = normalised(nearPt - pos);
 
-                    scalar normalAlignment = normal & pW/(mag(pW) + SMALL);
+                    const scalar normalAlignment = normal & pW;
 
                     // Find the patchIndex and wallData for WallSiteData object
                     label patchi = patchID[realFacei - mesh.nInternalFaces()];
@@ -315,15 +314,13 @@ void Foam::PairCollision<CloudType>::wallInteraction()
 
                 if (nearest.distance() < r)
                 {
-                    vector normal = rwf.normal(pts);
-
-                    normal /= mag(normal);
+                    const vector normal = rwf.unitNormal(pts);
 
                     const vector& nearPt = nearest.rawPoint();
 
-                    vector pW = nearPt - pos;
+                    const vector pW = normalised(nearPt - pos);
 
-                    scalar normalAlignment = normal & pW/mag(pW);
+                    const scalar normalAlignment = normal & pW;
 
                     // Find the patchIndex and wallData for WallSiteData object
 
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C
index e819c7fc9e..7adb7b5806 100644
--- a/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C
+++ b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C
@@ -344,7 +344,7 @@ Foam::vector Foam::PackingModels::Implicit<CloudType>::velocityCorrection
     const vector U = uCorrect_()[celli];
 
     // face geometry
-    vector nHat = mesh.faces()[facei].normal(mesh.points());
+    vector nHat = mesh.faces()[facei].areaNormal(mesh.points());
     const scalar nMag = mag(nHat);
     nHat /= nMag;
 
diff --git a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C
index d93a6985c9..52e3a0a8b9 100644
--- a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C
+++ b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C
@@ -249,10 +249,9 @@ void Foam::molecule::hitProcessorPatch(moleculeCloud&, trackingData& td)
 
 void Foam::molecule::hitWallPatch(moleculeCloud&, trackingData&)
 {
-    vector nw = normal();
-    nw /= mag(nw);
+    const vector nw = normal();
 
-    scalar vn = v_ & nw;
+    const scalar vn = v_ & nw;
 
     // Specular reflection
     if (vn > 0)
diff --git a/src/lagrangian/solidParticle/solidParticle.C b/src/lagrangian/solidParticle/solidParticle.C
index 568c1b4c7d..dfc1426e76 100644
--- a/src/lagrangian/solidParticle/solidParticle.C
+++ b/src/lagrangian/solidParticle/solidParticle.C
@@ -104,8 +104,7 @@ void Foam::solidParticle::hitProcessorPatch
 
 void Foam::solidParticle::hitWallPatch(solidParticleCloud& cloud, trackingData&)
 {
-    vector nw = normal();
-    nw /= mag(nw);
+    const vector nw = normal();
 
     scalar Un = U_ & nw;
     vector Ut = U_ - Un*nw;
diff --git a/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C
index 7593b3763e..175ecf9688 100644
--- a/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C
+++ b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C
@@ -59,7 +59,8 @@ void Foam::blockDescriptor::check(const Istream& is)
     forAll(faces, i)
     {
         point faceCentre(faces[i].centre(vertices_));
-        vector faceNormal(faces[i].normal(vertices_));
+        vector faceNormal(faces[i].areaNormal(vertices_));
+
         if (mag(faceNormal) > SMALL)
         {
             if (((faceCentre - blockCentre) & faceNormal) > 0)
diff --git a/src/mesh/blockMesh/blockMesh/blockMeshCheck.C b/src/mesh/blockMesh/blockMesh/blockMeshCheck.C
index fc0c08b862..25f1ddfdf9 100644
--- a/src/mesh/blockMesh/blockMesh/blockMeshCheck.C
+++ b/src/mesh/blockMesh/blockMesh/blockMeshCheck.C
@@ -187,9 +187,9 @@ void Foam::blockMesh::check(const polyMesh& bm, const dictionary& dict) const
                         if
                         (
                             (
-                                patchFace.normal(points)
-                              & faces[cellFaces[cellFacei]].normal(points)
-                            ) < 0.0
+                                patchFace.areaNormal(points)
+                              & faces[cellFaces[cellFacei]].areaNormal(points)
+                            ) < 0
                         )
                         {
                             Info<< tab << tab
diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementProblemCells.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementProblemCells.C
index 714aad7bc9..f229cfa391 100644
--- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementProblemCells.C
+++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementProblemCells.C
@@ -296,7 +296,7 @@ bool Foam::meshRefinement::isCollapsedFace
     const scalar severeNonorthogonalityThreshold =
         ::cos(degToRad(maxNonOrtho));
 
-    vector s = mesh_.faces()[facei].normal(points);
+    vector s = mesh_.faces()[facei].areaNormal(points);
     scalar magS = mag(s);
 
     // Check face area
diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.C b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.C
index a00e780401..8d2d28579b 100644
--- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.C
+++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.C
@@ -2458,10 +2458,11 @@ void Foam::snappySnapDriver::detectWarpedFaces
                     //Info<< "Splitting face:" << f << " into f0:" << f0
                     //    << " f1:" << f1 << endl;
 
-                    vector n0 = f0.normal(localPoints);
-                    scalar n0Mag = mag(n0);
-                    vector n1 = f1.normal(localPoints);
-                    scalar n1Mag = mag(n1);
+                    const vector n0 = f0.areaNormal(localPoints);
+                    const scalar n0Mag = mag(n0);
+
+                    const vector n1 = f1.areaNormal(localPoints);
+                    const scalar n1Mag = mag(n1);
 
                     if (n0Mag > ROOTVSMALL && n1Mag > ROOTVSMALL)
                     {
diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriverFeature.C b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriverFeature.C
index 7f7a3ac556..24e02f5ebb 100644
--- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriverFeature.C
+++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriverFeature.C
@@ -1860,9 +1860,9 @@ Foam::labelPair Foam::snappySnapDriver::findDiagonalAttraction
                         isConcave
                         (
                             compact0.centre(points0),
-                            compact0.normal(points0),
+                            compact0.areaNormal(points0),
                             compact1.centre(points1),
-                            compact1.normal(points1),
+                            compact1.areaNormal(points1),
                             concaveCos
                         )
                     )
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.C
index 0e503da082..9b1659e4b9 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.C
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.C
@@ -87,7 +87,7 @@ void Foam::directAMI<SourcePatch, TargetPatch>::appendToDirectSeeds
             // second attempt: match by shooting a ray into the tgt face
             if (!found)
             {
-                const vector srcN = srcF.normal(srcPoints);
+                const vector srcN = srcF.areaNormal(srcPoints);
 
                 for (const label tgtI : tgtNbr)
                 {
@@ -117,7 +117,7 @@ void Foam::directAMI<SourcePatch, TargetPatch>::appendToDirectSeeds
                 {
                     Pout<< "source face not found: id=" << srcI
                         << " centre=" << srcCf[srcI]
-                        << " normal=" << srcF.normal(srcPoints)
+                        << " normal=" << srcF.areaNormal(srcPoints)
                         << " points=" << srcF.points(srcPoints)
                         << endl;
 
@@ -128,7 +128,7 @@ void Foam::directAMI<SourcePatch, TargetPatch>::appendToDirectSeeds
 
                         Pout<< "face id: " << tgtI
                             << " centre=" << tgtF.centre(tgtPoints)
-                            << " normal=" << tgtF.normal(tgtPoints)
+                            << " normal=" << tgtF.areaNormal(tgtPoints)
                             << " points=" << tgtF.points(tgtPoints)
                             << endl;
                     }
diff --git a/src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C b/src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C
index 26ef1811a2..795397fc29 100644
--- a/src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C
+++ b/src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C
@@ -197,8 +197,8 @@ void Foam::cyclicAMIPolyPatch::calcTransforms
                 reduce(n0, maxMagSqrOp<point>());
                 reduce(n1, maxMagSqrOp<point>());
 
-                n0 /= mag(n0) + VSMALL;
-                n1 /= mag(n1) + VSMALL;
+                n0.normalise();
+                n1.normalise();
 
                 // Extended tensor from two local coordinate systems calculated
                 // using normal and rotation axis
@@ -367,14 +367,14 @@ void Foam::cyclicAMIPolyPatch::calcTransforms()
     vectorField half0Areas(half0.size());
     forAll(half0, facei)
     {
-        half0Areas[facei] = half0[facei].normal(half0.points());
+        half0Areas[facei] = half0[facei].areaNormal(half0.points());
     }
 
     const cyclicAMIPolyPatch& half1 = neighbPatch();
     vectorField half1Areas(half1.size());
     forAll(half1, facei)
     {
-        half1Areas[facei] = half1[facei].normal(half1.points());
+        half1Areas[facei] = half1[facei].areaNormal(half1.points());
     }
 
     calcTransforms
diff --git a/src/meshTools/indexedOctree/treeDataPrimitivePatch.C b/src/meshTools/indexedOctree/treeDataPrimitivePatch.C
index 04d6738258..3b3a8bb2f3 100644
--- a/src/meshTools/indexedOctree/treeDataPrimitivePatch.C
+++ b/src/meshTools/indexedOctree/treeDataPrimitivePatch.C
@@ -172,7 +172,7 @@ Foam::volumeType Foam::treeDataPrimitivePatch<PatchType>::getVolumeType
     // already have point.
 
     pointHit curHit = f.nearestPoint(sample, points);
-    const vector area = f.normal(points);
+    const vector area = f.areaNormal(points);
     const point& curPt = curHit.rawPoint();
 
     //
diff --git a/src/meshTools/primitiveMeshGeometry/primitiveMeshGeometry.C b/src/meshTools/primitiveMeshGeometry/primitiveMeshGeometry.C
index 05e7f05fa2..5ab9d5ea3f 100644
--- a/src/meshTools/primitiveMeshGeometry/primitiveMeshGeometry.C
+++ b/src/meshTools/primitiveMeshGeometry/primitiveMeshGeometry.C
@@ -1023,13 +1023,11 @@ bool Foam::primitiveMeshGeometry::checkFaceTwist
 
     label nWarped = 0;
 
-    forAll(checkFaces, i)
+    for (const label facei : checkFaces)
     {
-        label facei = checkFaces[i];
-
         const face& f = fcs[facei];
 
-        scalar magArea = mag(faceAreas[facei]);
+        const scalar magArea = mag(faceAreas[facei]);
 
         if (f.size() > 3 && magArea > VSMALL)
         {
@@ -1039,21 +1037,21 @@ bool Foam::primitiveMeshGeometry::checkFaceTwist
 
             forAll(f, fpI)
             {
-                vector triArea
+                const vector triArea
                 (
                     triPointRef
                     (
                         p[f[fpI]],
                         p[f.nextLabel(fpI)],
                         fc
-                    ).normal()
+                    ).areaNormal()
                 );
 
-                scalar magTri = mag(triArea);
+                const scalar magTri = mag(triArea);
 
                 if (magTri > VSMALL && ((nf & triArea/magTri) < minTwist))
                 {
-                    nWarped++;
+                    ++nWarped;
 
                     if (setPtr)
                     {
diff --git a/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.C b/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.C
index 752f245196..110f57feeb 100644
--- a/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.C
+++ b/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.C
@@ -752,8 +752,8 @@ void Foam::triSurfaceMesh::getNormal
         {
             if (info[i].hit())
             {
-                label facei = info[i].index();
-                normal[i] = s[facei].normal(pts);
+                const label facei = info[i].index();
+                normal[i] = s[facei].unitNormal(pts);
 
                 scalar qual = s[facei].tri(pts).quality();
 
@@ -768,12 +768,10 @@ void Foam::triSurfaceMesh::getNormal
                         if (nbrQual > qual)
                         {
                             qual = nbrQual;
-                            normal[i] = s[nbri].normal(pts);
+                            normal[i] = s[nbri].unitNormal(pts);
                         }
                     }
                 }
-
-                normal[i] /= mag(normal[i]) + VSMALL;
             }
             else
             {
@@ -789,12 +787,9 @@ void Foam::triSurfaceMesh::getNormal
             if (info[i].hit())
             {
                 const label facei = info[i].index();
-                // Cached:
-                //normal[i] = faceNormals()[facei];
 
                 // Uncached
-                normal[i] = s[facei].normal(pts);
-                normal[i] /= mag(normal[i]) + VSMALL;
+                normal[i] = s[facei].unitNormal(pts);
             }
             else
             {
diff --git a/src/meshTools/sets/cellSources/rotatedBoxToCell/rotatedBoxToCell.C b/src/meshTools/sets/cellSources/rotatedBoxToCell/rotatedBoxToCell.C
index 8ab9e1ac94..581214ba54 100644
--- a/src/meshTools/sets/cellSources/rotatedBoxToCell/rotatedBoxToCell.C
+++ b/src/meshTools/sets/cellSources/rotatedBoxToCell/rotatedBoxToCell.C
@@ -62,11 +62,7 @@ void Foam::rotatedBoxToCell::combine(topoSet& set, const bool add) const
     boxPoints[6] = origin_ + k_ + i_ + j_;
     boxPoints[7] = origin_ + k_ + j_;
 
-    labelList boxVerts(8);
-    forAll(boxVerts, i)
-    {
-        boxVerts[i] = i;
-    }
+    labelList boxVerts(identity(8));
 
     const cellModel& hex = cellModel::ref(cellModel::HEX);
 
@@ -77,7 +73,7 @@ void Foam::rotatedBoxToCell::combine(topoSet& set, const bool add) const
     vectorField boxFaceNormals(boxFaces.size());
     forAll(boxFaces, i)
     {
-        boxFaceNormals[i] = boxFaces[i].normal(boxPoints);
+        boxFaceNormals[i] = boxFaces[i].areaNormal(boxPoints);
 
         //Pout<< "Face:" << i << " position:" << boxFaces[i].centre(boxPoints)
         //    << " normal:" << boxFaceNormals[i] << endl;
diff --git a/src/meshTools/sets/faceZoneSources/setAndNormalToFaceZone/setAndNormalToFaceZone.C b/src/meshTools/sets/faceZoneSources/setAndNormalToFaceZone/setAndNormalToFaceZone.C
index 391ba26825..3b4e6ed5b7 100644
--- a/src/meshTools/sets/faceZoneSources/setAndNormalToFaceZone/setAndNormalToFaceZone.C
+++ b/src/meshTools/sets/faceZoneSources/setAndNormalToFaceZone/setAndNormalToFaceZone.C
@@ -124,7 +124,7 @@ void Foam::setAndNormalToFaceZone::applyToSet
                 {
                     newAddressing.append(facei);
 
-                    vector n = faces[facei].normal(points);
+                    const vector n = faces[facei].areaNormal(points);
                     if ((n & normal_) > 0)
                     {
                         newFlipMap.append(false);
diff --git a/src/meshTools/tetOverlapVolume/tetOverlapVolumeTemplates.C b/src/meshTools/tetOverlapVolume/tetOverlapVolumeTemplates.C
index 34eded1ff3..0852aa7a8d 100644
--- a/src/meshTools/tetOverlapVolume/tetOverlapVolumeTemplates.C
+++ b/src/meshTools/tetOverlapVolume/tetOverlapVolumeTemplates.C
@@ -52,7 +52,7 @@ void Foam::tetOverlapVolume::tetTetOverlap
     tetPointRef tetATet = tetA.tet();
     for (label facei = 0; facei < 4; ++facei)
     {
-        tetAFaceAreas[facei] = -tetATet.tri(facei).normal();
+        tetAFaceAreas[facei] = -tetATet.tri(facei).areaNormal();
         tetAMag2FaceAreas[facei] = magSqr(tetAFaceAreas[facei]);
         if (tetAMag2FaceAreas[facei] < ROOTVSMALL)
         {
@@ -65,7 +65,7 @@ void Foam::tetOverlapVolume::tetTetOverlap
     tetPointRef tetBTet = tetB.tet();
     for (label facei = 0; facei < 4; ++facei)
     {
-        tetBFaceAreas[facei] = -tetBTet.tri(facei).normal();
+        tetBFaceAreas[facei] = -tetBTet.tri(facei).areaNormal();
         tetBMag2FaceAreas[facei] = magSqr(tetBFaceAreas[facei]);
         if (tetBMag2FaceAreas[facei] < ROOTVSMALL)
         {
diff --git a/src/meshTools/triSurface/booleanOps/intersectedSurface/intersectedSurface.C b/src/meshTools/triSurface/booleanOps/intersectedSurface/intersectedSurface.C
index a142a37ee6..dd6857c0ad 100644
--- a/src/meshTools/triSurface/booleanOps/intersectedSurface/intersectedSurface.C
+++ b/src/meshTools/triSurface/booleanOps/intersectedSurface/intersectedSurface.C
@@ -1021,7 +1021,7 @@ Foam::faceList Foam::intersectedSurface::splitFace
     // See if normal needs flipping.
     faces.shrink();
 
-    vector n = faces[0].normal(eSurf.points());
+    const vector n = faces[0].areaNormal(eSurf.points());
 
     if ((n & surf.faceNormals()[facei]) < 0)
     {
diff --git a/src/meshTools/triSurface/faceTriangulation/faceTriangulation.C b/src/meshTools/triSurface/faceTriangulation/faceTriangulation.C
index c083f60146..8bf1e6332f 100644
--- a/src/meshTools/triSurface/faceTriangulation/faceTriangulation.C
+++ b/src/meshTools/triSurface/faceTriangulation/faceTriangulation.C
@@ -54,8 +54,8 @@ Foam::tmp<Foam::vectorField> Foam::faceTriangulation::calcEdges
     const pointField& points
 )
 {
-    tmp<vectorField> tedges(new vectorField(f.size()));
-    vectorField& edges = tedges.ref();
+    auto tedges = tmp<vectorField>::New(f.size());
+    auto& edges = tedges.ref();
 
     forAll(f, i)
     {
@@ -63,9 +63,8 @@ Foam::tmp<Foam::vectorField> Foam::faceTriangulation::calcEdges
         point nextPt = points[f[f.fcIndex(i)]];
 
         vector vec(nextPt - thisPt);
-        vec /= mag(vec) + VSMALL;
 
-        edges[i] = vec;
+        edges[i] = vec.normalise();
     }
 
     return tedges;
@@ -159,9 +158,9 @@ bool Foam::faceTriangulation::triangleContainsPoint
     const point& pt
 )
 {
-    scalar area01Pt = triPointRef(p0, p1, pt).normal() & n;
-    scalar area12Pt = triPointRef(p1, p2, pt).normal() & n;
-    scalar area20Pt = triPointRef(p2, p0, pt).normal() & n;
+    const scalar area01Pt = triPointRef(p0, p1, pt).areaNormal() & n;
+    const scalar area12Pt = triPointRef(p1, p2, pt).areaNormal() & n;
+    const scalar area20Pt = triPointRef(p2, p0, pt).areaNormal() & n;
 
     if ((area01Pt > 0) && (area12Pt > 0) && (area20Pt > 0))
     {
@@ -172,10 +171,8 @@ bool Foam::faceTriangulation::triangleContainsPoint
         FatalErrorInFunction << abort(FatalError);
         return false;
     }
-    else
-    {
-        return false;
-    }
+
+    return false;
 }
 
 
@@ -209,7 +206,7 @@ void Foam::faceTriangulation::findDiagonal
     );
     // rayDir should be normalized already but is not due to rounding errors
     // so normalize.
-    rayDir /= mag(rayDir) + VSMALL;
+    rayDir.normalise();
 
 
     //
@@ -309,8 +306,8 @@ void Foam::faceTriangulation::findDiagonal
         {
             // pt inside triangle (so perhaps visible)
             // Select based on minimal angle (so guaranteed visible).
-            vector edgePt0 = pt - startPt;
-            edgePt0 /= mag(edgePt0);
+            vector edgePt0 = (pt - startPt);
+            edgePt0.normalise();
 
             scalar cos = rayDir & edgePt0;
             if (cos > maxCos)
@@ -598,14 +595,12 @@ bool Foam::faceTriangulation::split
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Null constructor
 Foam::faceTriangulation::faceTriangulation()
 :
     triFaceList()
 {}
 
 
-// Construct from components
 Foam::faceTriangulation::faceTriangulation
 (
     const pointField& points,
@@ -615,8 +610,7 @@ Foam::faceTriangulation::faceTriangulation
 :
     triFaceList(f.size()-2)
 {
-    vector avgNormal = f.normal(points);
-    avgNormal /= mag(avgNormal) + VSMALL;
+    const vector avgNormal = f.unitNormal(points);
 
     label triI = 0;
 
@@ -629,7 +623,6 @@ Foam::faceTriangulation::faceTriangulation
 }
 
 
-// Construct from components
 Foam::faceTriangulation::faceTriangulation
 (
     const pointField& points,
@@ -651,7 +644,6 @@ Foam::faceTriangulation::faceTriangulation
 }
 
 
-// Construct from Istream
 Foam::faceTriangulation::faceTriangulation(Istream& is)
 :
     triFaceList(is)
diff --git a/src/meshTools/triSurface/triSurfaceTools/triSurfaceCurvature.C b/src/meshTools/triSurface/triSurfaceTools/triSurfaceCurvature.C
index f64cb0eea2..e179edd35b 100644
--- a/src/meshTools/triSurface/triSurfaceTools/triSurfaceCurvature.C
+++ b/src/meshTools/triSurface/triSurfaceTools/triSurfaceCurvature.C
@@ -68,8 +68,8 @@ Foam::triSurfaceTools::vertexNormals(const triSurface& surf)
 
     Info<< "Calculating vertex normals" << endl;
 
-    tmp<vectorField> tfld(new vectorField(surf.nPoints(), Zero));
-    vectorField& pointNormals = tfld.ref();
+    auto tpointNormals = tmp<vectorField>::New(surf.nPoints(), Zero);
+    auto& pointNormals = tpointNormals.ref();
 
     const pointField& points = surf.points();
     const labelListList& pointFaces = surf.pointFaces();
@@ -79,28 +79,27 @@ Foam::triSurfaceTools::vertexNormals(const triSurface& surf)
     {
         const labelList& pFaces = pointFaces[pI];
 
-        forAll(pFaces, fI)
+        for (const label facei : pFaces)
         {
-            const label facei = pFaces[fI];
             const triFace& f = surf[facei];
 
-            vector fN = f.normal(points);
+            const vector areaNorm = f.areaNormal(points);
 
             const scalar weight = vertexNormalWeight
             (
                 f,
                 meshPoints[pI],
-                fN,
+                areaNorm,
                 points
             );
 
-            pointNormals[pI] += weight*fN;
+            pointNormals[pI] += weight * areaNorm;
         }
 
-        pointNormals[pI] /= mag(pointNormals[pI]) + VSMALL;
+        pointNormals[pI].normalise();
     }
 
-    return tfld;
+    return tpointNormals;
 }
 
 
@@ -114,8 +113,8 @@ Foam::triSurfaceTools::vertexTriads
     const pointField& points = surf.points();
     const Map<label>& meshPointMap = surf.meshPointMap();
 
-    tmp<triadField> tfld(new triadField(points.size()));
-    triadField& pointTriads = tfld.ref();
+    auto tpointTriads = tmp<triadField>::New(points.size());
+    auto& pointTriads = tpointTriads.ref();
 
     forAll(points, pI)
     {
@@ -131,16 +130,13 @@ Foam::triSurfaceTools::vertexTriads
         plane p(pt, normal);
 
         // Pick arbitrary point in plane
-        vector dir1 = pt - p.somePointInPlane(1e-3);
-        dir1 /= mag(dir1);
-
-        vector dir2 = dir1 ^ normal;
-        dir2 /= mag(dir2);
+        vector dir1 = normalised(pt - p.somePointInPlane(1e-3));
+        vector dir2 = normalised(dir1 ^ normal);
 
         pointTriads[meshPointMap[pI]] = triad(dir1, dir2, normal);
     }
 
-    return tfld;
+    return tpointTriads;
 }
 
 
@@ -187,7 +183,7 @@ Foam::triSurfaceTools::curvatures
 
         // Set up a local coordinate system for the face
         const vector& e0 = edgeVectors[0];
-        const vector eN = f.normal(points);
+        const vector eN = f.areaNormal(points);
         const vector e1 = (e0 ^ eN);
 
         if (magSqr(eN) < ROOTVSMALL)
@@ -276,7 +272,7 @@ Foam::triSurfaceTools::curvatures
             (
                 f,
                 meshPoints[patchPointIndex],
-                f.normal(points),
+                f.areaNormal(points),
                 points
             );
 
@@ -304,8 +300,8 @@ Foam::triSurfaceTools::curvatures
         }
     }
 
-    tmp<scalarField> tfld(new scalarField(points.size(), Zero));
-    scalarField& curvatureAtPoints = tfld.ref();
+    auto tcurvatureAtPoints = tmp<scalarField>::New(points.size(), Zero);
+    scalarField& curvatureAtPoints = tcurvatureAtPoints.ref();
 
     forAll(curvatureAtPoints, pI)
     {
@@ -325,7 +321,7 @@ Foam::triSurfaceTools::curvatures
         curvatureAtPoints[meshPoints[pI]] = curvature;
     }
 
-    return tfld;
+    return tcurvatureAtPoints;
 }
 
 
@@ -356,8 +352,8 @@ Foam::triSurfaceTools::writeCurvature
 {
     Info<< "Extracting curvature of surface at the points." << endl;
 
-    tmp<scalarField> tfld = triSurfaceTools::curvatures(surf);
-    scalarField& curv = tfld.ref();
+    tmp<scalarField> tcurv = triSurfaceTools::curvatures(surf);
+    scalarField& curv = tcurv.ref();
 
     triSurfacePointScalarField outputField
     (
@@ -379,7 +375,7 @@ Foam::triSurfaceTools::writeCurvature
     outputField.write();
     outputField.swap(curv);
 
-    return tfld;
+    return tcurv;
 }
 
 
diff --git a/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C b/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C
index 07e3fb0ed4..64766a1df7 100644
--- a/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C
+++ b/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C
@@ -1891,8 +1891,7 @@ void Foam::distributedTriSurfaceMesh::getNormal
     forAll(triangleIndex, i)
     {
         label trii = triangleIndex[i];
-        normal[i] = s[trii].normal(s.points());
-        normal[i] /= mag(normal[i]) + VSMALL;
+        normal[i] = s[trii].unitNormal(s.points());
     }
 
 
diff --git a/src/sampling/surface/isoSurface/isoSurfaceCell.C b/src/sampling/surface/isoSurface/isoSurfaceCell.C
index d6c0090be4..6cecc030c2 100644
--- a/src/sampling/surface/isoSurface/isoSurfaceCell.C
+++ b/src/sampling/surface/isoSurface/isoSurfaceCell.C
@@ -295,8 +295,8 @@ Foam::pointIndexHit Foam::isoSurfaceCell::collapseSurface
 
         if (shared[0] != -1)
         {
-            const vector n0 = tri0.normal(localPoints);
-            const vector n1 = tri1.normal(localPoints);
+            const vector n0 = tri0.areaNormal(localPoints);
+            const vector n1 = tri1.areaNormal(localPoints);
 
             // Merge any zero-sized triangles,
             // or if they point in the same direction.
diff --git a/src/surfMesh/surfaceFormats/stl/STLsurfaceFormat.C b/src/surfMesh/surfaceFormats/stl/STLsurfaceFormat.C
index fa421aa73d..a6cf610c74 100644
--- a/src/surfMesh/surfaceFormats/stl/STLsurfaceFormat.C
+++ b/src/surfMesh/surfaceFormats/stl/STLsurfaceFormat.C
@@ -37,9 +37,9 @@ inline void Foam::fileFormats::STLsurfaceFormat<Face>::writeShell
     const Face& f
 )
 {
-    // calculate the normal ourselves, for flexibility and speed
-    vector norm = triPointRef(pts[f[0]], pts[f[1]], pts[f[2]]).normal();
-    norm /= mag(norm) + VSMALL;
+    // Calculate the normal ourselves, for flexibility and speed
+    const vector norm =
+        triPointRef(pts[f[0]], pts[f[1]], pts[f[2]]).unitNormal();
 
     // simple triangulation about f[0].
     // better triangulation should have been done before
@@ -71,9 +71,9 @@ inline void Foam::fileFormats::STLsurfaceFormat<Face>::writeShell
     const label zoneI
 )
 {
-    // calculate the normal ourselves, for flexibility and speed
-    vector norm = triPointRef(pts[f[0]], pts[f[1]], pts[f[2]]).normal();
-    norm /= mag(norm) + VSMALL;
+    // Calculate the normal ourselves, for flexibility and speed
+    const vector norm =
+        triPointRef(pts[f[0]], pts[f[1]], pts[f[2]]).unitNormal();
 
     // simple triangulation about f[0].
     // better triangulation should have been done before