ENH: distinguish between face areaNormal/unitNormal in the code

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];
+        const vector areaNorm(f.areaNormal(points));
-
-        vector n(f.normal(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;

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();

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)

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());
+                const vector areaNorm = faces[facei].areaNormal(mesh.points());
-                fN /= mag(fN) + SMALL;
 
-                if ((norm[0] & fN) < 0)
+                if ((norm[0] & areaNorm) < 0)
                 {
                     flipMap[facei] = true;
                 }

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 =
+                const label patchID = f.region() + regionOffset_[s];
-                    triSurf[sI].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();
         }

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));
+        //const vector n(newFace.unitNormal(dualPoints));
-        //n /= mag(n);
+        //
         //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));
+        //const vector n(newFace.unitNormal(dualPoints));
-        //n /= mag(n);
+        //
         //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));
+    //const vector n(newFace.unitNormal(dualPoints));
-    //n /= mag(n);
+    //
     //Pout<< "Generated boundary dualFace:" << dualFacei
     //    << " verts:" << newFace
     //    << " points:" << UIndirectList<point>(meshMod.points(), newFace)

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
+    auto tpointNormals = tmp<vectorField>::New(surf.nPoints(), Zero);
-    (
+    auto& pointNormals = tpointNormals.ref();
-        new pointField(surf.nPoints(), Zero)
-    );
-    vectorField& 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();
 

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;
                     }

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;
+            n0.normalise();
-            n1 /= mag(n1) + VSMALL;
+            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]"
+                    << " swept angle: " << radToDeg(acos(n0 & n1))
-                    << endl;
+                    << " [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;
+                n0.normalise();
-                n1 /= mag(n1) + VSMALL;
+                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());
+                const vector n0 = pp0[max0I].unitNormal(pp0.points());
-                n0 /= mag(n0) + VSMALL;
 
                 label max1I = findMaxArea(pp1.points(), pp1);
-                vector n1 = pp1[max1I].normal(pp1.points());
+                const vector n1 = pp1[max1I].unitNormal(pp1.points());
-                n1 /= mag(n1) + VSMALL;
 
                 if (mag(n0 & n1) < 1-matchTolerance())
                 {

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());
+            const vector n0 = half0Faces[max0I].unitNormal(pp.points());
-            n0 /= mag(n0) + VSMALL;
 
             label max1I = findMaxArea(pp.points(), half1Faces);
-            vector n1 = half1Faces[max1I].normal(pp.points());
+            const vector n1 = half1Faces[max1I].unitNormal(pp.points());
-            n1 /= mag(n1) + VSMALL;
 
             if (mag(n0 & n1) < 1-matchTolerance())
             {

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] = this->operator[](facei).unitNormal(points_);
-        n[facei] /= mag(n[facei]) + VSMALL;
     }
 
     if (debug)

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);
+            const vector areaNorm = f.areaNormal(dualPoints);
-            if (((mesh.points()[owner] - dualPoints[f[0]]) & n) > 0)
+
+            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);
+                const vector areaNorm = f.areaNormal(dualPoints);
-                if (((mesh.points()[owner] - dualPoints[f[0]]) & n) > 0)
+
+                if (((mesh.points()[owner] - dualPoints[f[0]]) & areaNorm) > 0)
                 {
                     WarningInFunction
                         << " on internal edge:" << edgeI

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);
+    const vector areaNorm = f.areaNormal(loopPts);
-    point ctr = f.centre(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;
-    }
 }
 
 

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);
 

src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.C

@@ -55,9 +55,8 @@ bool Foam::combineFaces::convexFace
     const face& f
 )
 {
-    // Get outwards pointing normal of f.
+    // Get outwards pointing normal of f, only the sign matters.
-    vector n = f.normal(points);
+    const vector areaNorm = f.areaNormal(points);
-    n /= mag(n);
 
     // 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];
+            vector f0Normal = normalised(mesh_.faceAreas()[f0]);
-            f0Normal /= mag(f0Normal);
+            vector f1Normal = normalised(mesh_.faceAreas()[f1]);
-            vector f1Normal = mesh_.faceAreas()[f1];
-            f1Normal /= mag(f1Normal);
 
-            if ((f0Normal&f1Normal) > minCos)
+            if ((f0Normal & f1Normal) > minCos)
             {
                 Map<label>::const_iterator f0Fnd = faceRegion.find(f0);
 

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)
     {

src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatchCutFaces.C

@@ -150,8 +150,7 @@ void Foam::enrichedPatch::calcCutFaces() const
         }
 
         // Grab face normal
-        vector normal = curLocalFace.normal(lp);
+        const vector normal = curLocalFace.unitNormal(lp);
-        normal /= mag(normal);
 
         while (edgeSeeds.size())
         {

src/fileFormats/stl/STLtriangleI.H

@@ -143,9 +143,8 @@ inline void Foam::STLtriangle::write
     const point& pt2
 )
 {
-    // calculate the normal ourselves
+    // Calculate the normal ourselves
-    vector norm = triPointRef(pt0, pt1, pt2).normal();
+    const vector norm = triPointRef(pt0, pt1, pt2).unitNormal();
-    norm /= mag(norm) + VSMALL;
 
     write(os, norm, pt0, pt1, pt2);
 }

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]);

src/finiteArea/faMesh/faMeshDemandDrivenData.C

@@ -425,9 +425,7 @@ void Foam::faMesh::calcFaceAreaNormals() const
     vectorField& nInternal = faceAreaNormals.ref();
     forAll(localFaces, faceI)
     {
-        nInternal[faceI] =
+        nInternal[faceI] = localFaces[faceI].unitNormal(localPoints);
-            localFaces[faceI].normal(localPoints)/
-            localFaces[faceI].mag(localPoints);
     }
 
     forAll(boundary(), patchI)

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());
+    auto tfN = tmp<vectorField>::New();
-    vectorField& fN = tfN.ref();
+    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)
+        fN[faceI] = faces[ngbFaces[faceI]].unitNormal(points);
-            /faces[ngbFaces[faceI]].mag(points);
     }
 
     return tfN;

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_

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();
+    const vector nw = p.normal();
-    nw /= mag(nw);
 
     // Normal velocity magnitude
     scalar U_dot_nw = U & nw;

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();
+    const vector nw = p.normal();
-    nw /= mag(nw);
 
     // Normal velocity magnitude
     scalar U_dot_nw = U & nw;

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();
+    const vector nw = p.normal();
-    nw /= mag(nw);
 
     scalar U_dot_nw = U & nw;
 

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
+            //- The (unit) normal of the tri on tetFacei_ for the current tet.
-            //  current tet.
             inline vector normal() const;
 
             //- Is the particle on a face?

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 = triPointRef(base[0], vertex1[0], vertex2[0]).unitNormal();
-        n /= mag(n);
 
         // 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 = triPointRef(base, vertex1, vertex2).unitNormal();
-        n /= mag(n);
 
         U = Zero;
     }

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();
+    const vector nf = normal();
-    nf /= mag(nf);
 
     transformProperties(I - 2.0*nf*nf);
 }

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;

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];
+                    const vector normal =
-
+                        normalised(mesh.faceAreas()[realFacei]);
-                    normal /= mag(normal);
 
                     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);
+                    const vector normal = rwf.unitNormal(pts);
-
-                    normal /= mag(normal);
 
                     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
 

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;
 

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();
+    const vector nw = normal();
-    nw /= mag(nw);
 
-    scalar vn = v_ & nw;
+    const scalar vn = v_ & nw;
 
     // Specular reflection
     if (vn > 0)

src/lagrangian/solidParticle/solidParticle.C

@@ -104,8 +104,7 @@ void Foam::solidParticle::hitProcessorPatch
 
 void Foam::solidParticle::hitWallPatch(solidParticleCloud& cloud, trackingData&)
 {
-    vector nw = normal();
+    const vector nw = normal();
-    nw /= mag(nw);
 
     scalar Un = U_ & nw;
     vector Ut = U_ - Un*nw;

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)

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)
+                                patchFace.areaNormal(points)
-                              & faces[cellFaces[cellFacei]].normal(points)
+                              & faces[cellFaces[cellFacei]].areaNormal(points)
-                            ) < 0.0
+                            ) < 0
                         )
                         {
                             Info<< tab << tab

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

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);
+                    const vector n0 = f0.areaNormal(localPoints);
-                    scalar n0Mag = mag(n0);
+                    const scalar n0Mag = mag(n0);
-                    vector n1 = f1.normal(localPoints);
+
-                    scalar n1Mag = mag(n1);
+                    const vector n1 = f1.areaNormal(localPoints);
+                    const scalar n1Mag = mag(n1);
 
                     if (n0Mag > ROOTVSMALL && n1Mag > ROOTVSMALL)
                     {

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
                         )
                     )

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;
                     }

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;
+                n0.normalise();
-                n1 /= mag(n1) + VSMALL;
+                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

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();
 
     //

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)
                     {

src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.C

@@ -752,8 +752,8 @@ void Foam::triSurfaceMesh::getNormal
         {
             if (info[i].hit())
             {
-                label facei = info[i].index();
+                const label facei = info[i].index();
-                normal[i] = s[facei].normal(pts);
+                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] = s[facei].unitNormal(pts);
-                normal[i] /= mag(normal[i]) + VSMALL;
             }
             else
             {

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);
+    labelList boxVerts(identity(8));
-    forAll(boxVerts, i)
-    {
-        boxVerts[i] = i;
-    }
 
     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;

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);

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)
         {

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)
     {

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()));
+    auto tedges = tmp<vectorField>::New(f.size());
-    vectorField& edges = tedges.ref();
+    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;
+    const scalar area01Pt = triPointRef(p0, p1, pt).areaNormal() & n;
-    scalar area12Pt = triPointRef(p1, p2, pt).normal() & n;
+    const scalar area12Pt = triPointRef(p1, p2, pt).areaNormal() & n;
-    scalar area20Pt = triPointRef(p2, p0, pt).normal() & 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;
-    }
 }
 
 
@@ -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;
+            vector edgePt0 = (pt - startPt);
-            edgePt0 /= mag(edgePt0);
+            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);
+    const vector avgNormal = f.unitNormal(points);
-    avgNormal /= mag(avgNormal) + VSMALL;
 
     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)

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));
+    auto tpointNormals = tmp<vectorField>::New(surf.nPoints(), Zero);
-    vectorField& pointNormals = tfld.ref();
+    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()));
+    auto tpointTriads = tmp<triadField>::New(points.size());
-    triadField& pointTriads = tfld.ref();
+    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);
+        vector dir1 = normalised(pt - p.somePointInPlane(1e-3));
-        dir1 /= mag(dir1);
+        vector dir2 = normalised(dir1 ^ normal);
-
-        vector dir2 = dir1 ^ normal;
-        dir2 /= mag(dir2);
 
         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));
+    auto tcurvatureAtPoints = tmp<scalarField>::New(points.size(), Zero);
-    scalarField& curvatureAtPoints = tfld.ref();
+    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);
+    tmp<scalarField> tcurv = triSurfaceTools::curvatures(surf);
-    scalarField& curv = tfld.ref();
+    scalarField& curv = tcurv.ref();
 
     triSurfacePointScalarField outputField
     (
@@ -379,7 +375,7 @@ Foam::triSurfaceTools::writeCurvature
     outputField.write();
     outputField.swap(curv);
 
-    return tfld;
+    return tcurv;
 }
 
 

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] = s[trii].unitNormal(s.points());
-        normal[i] /= mag(normal[i]) + VSMALL;
     }
 
 

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 n0 = tri0.areaNormal(localPoints);
-            const vector n1 = tri1.normal(localPoints);
+            const vector n1 = tri1.areaNormal(localPoints);
 
             // Merge any zero-sized triangles,
             // or if they point in the same direction.

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
+    // Calculate the normal ourselves, for flexibility and speed
-    vector norm = triPointRef(pts[f[0]], pts[f[1]], pts[f[2]]).normal();
+    const vector norm =
-    norm /= mag(norm) + VSMALL;
+        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
+    // Calculate the normal ourselves, for flexibility and speed
-    vector norm = triPointRef(pts[f[0]], pts[f[1]], pts[f[2]]).normal();
+    const vector norm =
-    norm /= mag(norm) + VSMALL;
+        triPointRef(pts[f[0]], pts[f[1]], pts[f[2]]).unitNormal();
 
     // simple triangulation about f[0].
     // better triangulation should have been done before