diff --git a/applications/utilities/mesh/conversion/gmshToFoam/gmshToFoam.C b/applications/utilities/mesh/conversion/gmshToFoam/gmshToFoam.C
index d0dd7c4a4e490fd26fb0b40159579e4347d9281e..51e92215775b25f51df82fe73a02cd3224549e8f 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 809f996afc3b1e1e645c34da30f0a13c76d6a113..48b55bd14b36fec4c8e3c42c6100fd65c192d3a9 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 fb9899438a7c766ac5d604e4e3b46a7850501f4b..00bad93cd888b6c2ae008d252f06b23fe1028e2c 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 1f42b3e429d3fc0ad7473e352055c29e3f61025d..2a37049388999a8024defbb3a211357de420c2ff 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 0afaae03f389e225b883b53389a161e519752be6..37db9dab40674b23eb686328c935856dbb3c05e1 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 cef8d2695b18f85299456b28a35549ad8b795bf4..ede9932cd367460046f8e8d18e386653bc2476bb 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 83c5d147a0949fb0e3f393286acaa14cd1cf9c84..3d5e9fd906a2dfe559a074d4f2b11e9f2eb56e7d 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 f66b48555f3dc72127b7cb0ad549e4d6b1952852..c25c9fb4a8ece2470900d2c6b4de19bcaba4d20a 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 4257571d72a5668a8585f01f98c093f939154e40..845977c910fd12a804053bfc368659cc9546d7ff 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 80714b0bb471e8e7ace2b5a72da25eed90406fce..d04552e9343bb0386faccdcd9b44d2becf6f22c8 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 50e4222be0a54fee5b911ae385cacf592d34ce6f..a513934b28c35bb0ad1e064ece8a6a074226bbfe 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 b0c1f2eaead38c3ef52c1e89c6c23aac8d0b7fe9..9b325c94157f8c1672cb026b57e1b86a9e0e6d8e 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 1dbbc324efb4cba52e89f8d117b19090efe908c6..edb8209d6f5f4e9a43a1cb7fa98b0f342d040d43 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 39807dc780e23b3f7b334f6775e0078a7e5e71be..c5419bfc10b3b3e816470dd6686fab0eb7b458fd 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 426866008b5094dd8e70de8ae1b7feae598c1669..44c3ca0afce090a18ee6eaa42a242f34622edee4 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 e58aa2ae8ea6e69ff5d9104bb3aa5760648efcfb..8d1b96c318499668457b4aa963d92125ced6d1a9 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 4484c7f85241d3b795d77ebb57819d91fa1e4a54..48879430216c50a727e655d9aeec95fe125c7f23 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 541bec22bfe0bc723ae0664fd1ef9ed6d8edb2be..70529371b48d3d35c59acb293b2530f6d16eed91 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 c20f6276678a8a918d29fc553156d15c038240c0..3094f203e35de4490b82798a3d9446ba73a3aa3a 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 25552744506af79abd9d86229a66b741b581c14f..e9cc68d9cc27cf9044d2b3e5d7dc200727bd1b33 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 b0c96720a158370161e3aadb142dfe6a6392e979..8c64de4c8a6012a7dffff4e4cf62ffdcab0997c9 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 41b6eb06bfc5a6a2df4870c81e14d07647b18bd2..ed0c304a218900a07c1d4206b0b6aa2abee618ed 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 dd145c62e1c2c966ee43f07801aaa9803ed9fdd0..ca2bc100418c05ef57c49fdfd0252b6c65b6ed3d 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 21c2cd73a8b56b775b3807fe33fa5113c3182c56..4fe4f0649d798a393d9de5c0d84e9b9ff2da5786 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 a3d182d07b91067f9adfd6f4d2f9fe0c8864929f..2d9a00df536d4c5b01e3cd08d4f66bdc15d061a5 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 4d6928e2b1842dbd8290bbb42c07b28da96e1613..f7e9f445327e2f0e8a73f602a429c099f7acc269 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 771eeec0575bcd9da12ef9d1e213ec1ce92890e3..3698b17f2232b08aae8f030ec1b4e76afa388ca7 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 a43c7763972a286767974e6d8d8c44a02816d5a2..fffeacd3f0bac13864d530481a75aa62db4d071b 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 1a8255ae7d615056a3be9057c3adac92d0b371c6..0adce43d200f13c200ee368c42308591d7491f55 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 7d68828701da8dd440f70f354ac2d6650df8fd50..96407263374762f04cb4ed0c794aebc40055adf6 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 e819c7fc9e9f4191d0bc81acae178eb5f7f809a7..7adb7b5806bd20a6e5caa026a6c6f13b24f71aa0 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 d93a6985c92dba86afc63c268a37f3ebced629f6..52e3a0a8b9ef2b22a8f50813932935b1da1d5d7d 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 568c1b4c7d66721048ca45761f81cc39a162a6e3..dfc1426e760655ee56d70359855b3cade3225650 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 7593b3763e16bbe71e1a4663924269d1cd55fc7b..175ecf96888184028224cedf90e0a71afef86f04 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 fc0c08b8623443bfcd277132f852d8755a958a00..25f1ddfdf9f2019e0b6f7ee06f0836232d8fa1a5 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 714aad7bc93f768f8eb51bab853c37c278bc05ab..f229cfa39139ac730c20f454aa3e72b5d909428d 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 a00e7804014ca33059eaa883e19d22dd20972b4b..8d2d28579b1ef2b707d44962fb4c949ea8759a6b 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 7f7a3ac556a71fd8e2f5cab95f3eed5b0194afca..24e02f5ebb0d2555d9f23fd453ca11b023640bd4 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 0e503da082850722d4316e8e2b876c11a2f8d0e7..9b1659e4b9d70f4dacfb61d0a2b947ed5ebd3e20 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 26ef1811a2f0891b8d8a3b7c8aaac12596de7177..795397fc2910430bbdb43872db34e63537718f37 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 04d67382580e86cddce90113e639314eb4dba78e..3b3a8bb2f3c03cf417d71993f55ac1c058e362fd 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 05e7f05fa2ec0f733cd2f54fe4d140ce96ffba1a..5ab9d5ea3f644b6e5bb4299706b32647a62d1f00 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 752f245196606eaddb198c7633b3c6efd5f080a2..110f57feeb9b75ccdf069d64022a7dc4c689b01b 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 8ab9e1ac9494dcbdd73f21a5a8b2eac50796c31b..581214ba5489f8cb1873cb742d21a8aff93f0ddc 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 391ba268250fbd0d06f9f8f5334c900e70fc48eb..3b4e6ed5b7a7978beca40adcbbcf66936279b8b8 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 34eded1ff313478509b37fc627b996297d358213..0852aa7a8d8f18596d77019544dcfd9772e3f638 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 a142a37ee6be4e4ad58b7d09d8e6b35c5a275e7e..dd6857c0adfa638d238e3966cae93a96ef6aabb6 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 c083f60146787e2f8879ec4aa068542444acfdf6..8bf1e6332f0d15768db6e1f6ee03cc4187489474 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 f64cb0eea2bf2d992f70eac5a851dd9cb3c7babd..e179edd35b7be16b27e4806e8ee3cb301714c639 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 07e3fb0ed412a1f28229decf434cb08ddaf5e7ef..64766a1df7f018cb315b40c304be0f02ca7b46e7 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 d6c0090be4898ac567ce3be758cdc9cca8b3a97f..6cecc030c2180a00fc5820411c3e65ec68fefa62 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 fa421aa73dbc4100baaa973158156c6d4ea5714c..a6cf610c7466b1e90fd824b14b21c728696b1c78 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