ENH: use vector::normalise and VectorSpace::normalised for clarity

applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/interpolatedFaces.H

@@ -239,11 +239,9 @@ forAll(isNeiInterpolatedFace, faceI)
                 n1 = vector(-n.z(), 0, n.x());
             }
         }
+        n1.normalise();
 
-        n1 /= mag(n1);
+        const vector n2 = normalised(n ^ n1);
-
-        vector n2 = n ^ n1;
-        n2 /= mag(n2);
 
         tensor rot =
             tensor

applications/solvers/multiphase/icoReactingMultiphaseInterFoam/laserDTRM/laserDTRM.C

@@ -165,8 +165,7 @@ void Foam::radiation::laserDTRM::initialise()
 
     const scalar t = mesh_.time().value();
     const vector lPosition = focalLaserPosition_->value(t);
-    vector lDir = laserDirection_->value(t);
+    const vector lDir = normalised(laserDirection_->value(t));
-    lDir /= mag(lDir);
 
     if (debug)
     {
@@ -175,7 +174,7 @@ void Foam::radiation::laserDTRM::initialise()
     }
 
     // Find a vector on the area plane. Normal to laser direction
-    vector rArea = vector::zero;
+    vector rArea = Zero;
     scalar magr = 0.0;
 
     {
@@ -188,7 +187,7 @@ void Foam::radiation::laserDTRM::initialise()
             magr = mag(rArea);
         }
     }
-    rArea /= mag(rArea);
+    rArea.normalise();
 
     scalar dr =  focalLaserRadius_/ndr_;
     scalar dTheta =  mathematical::twoPi/ndTheta_;

applications/utilities/mesh/advanced/splitCells/splitCells.C

@@ -135,8 +135,11 @@ bool largerAngle
     // Get cos between faceCentre and normal vector to determine in
     // which quadrant angle is. (Is correct for unwarped faces only!)
     // Correct for non-outwards pointing normal.
-    vector c1c0(mesh.faceCentres()[f1] - mesh.faceCentres()[f0]);
+    const vector c1c0 =
-    c1c0 /= mag(c1c0) + VSMALL;
+        normalised
+        (
+            mesh.faceCentres()[f1] - mesh.faceCentres()[f0]
+        );
 
     scalar fcCosAngle = n0 & c1c0;
 
@@ -327,8 +330,7 @@ bool splitCell
     // halfway on fully regular meshes (since we want cuts
     // to be snapped to vertices)
     planeN += 0.01*halfNorm;
-
+    planeN.normalise();
-    planeN /= mag(planeN);
 
     // Define plane through edge
     plane cutPlane(mesh.points()[e.start()], planeN);
@@ -409,11 +411,8 @@ void collectCuts
             label f0, f1;
             meshTools::getEdgeFaces(mesh, celli, edgeI, f0, f1);
 
-            vector n0 = faceAreas[f0];
+            const vector n0 = normalised(faceAreas[f0]);
-            n0 /= mag(n0);
+            const vector n1 = normalised(faceAreas[f1]);
-
-            vector n1 = faceAreas[f1];
-            n1 /= mag(n1);
 
             if
             (

applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/cellAspectRatioControl/cellAspectRatioControl.C

@@ -45,7 +45,7 @@ Foam::cellAspectRatioControl::cellAspectRatioControl
     )
 {
     // Normalise the direction
-    aspectRatioDirection_ /= mag(aspectRatioDirection_) + SMALL;
+    aspectRatioDirection_.normalise();
 
     Info<< nl
         << "Cell Aspect Ratio Control" << nl

applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/cellSizeAndAlignmentControl/searchableSurfaceControl/searchableSurfaceControl.C

@@ -428,7 +428,7 @@ void Foam::searchableSurfaceControl::initialVertices
 
                 if (mag(normals[0]) < SMALL)
                 {
-                    normals[0] = vector(1, 1, 1);
+                    normals[0] = vector::one;
                 }
 
                 pointAlignment.reset(new triad(normals[0]));

applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C

@@ -1137,7 +1137,7 @@ void Foam::conformalVoronoiMesh::move()
 
                         alignmentDirs[aA] = a + sign(dotProduct)*b;
 
-                        alignmentDirs[aA] /= mag(alignmentDirs[aA]);
+                        alignmentDirs[aA].normalise();
                     }
                 }
             }

applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C

@@ -220,8 +220,7 @@ void Foam::conformalVoronoiMesh::createEdgePointGroupByCirculating
 //            << endl;
 
         // Calculate master point
-        vector masterPtVec(normalDir + nextNormalDir);
+        const vector masterPtVec = normalised(normalDir + nextNormalDir);
-        masterPtVec /= mag(masterPtVec) + SMALL;
 
         if
         (

applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C

@@ -738,8 +738,11 @@ Foam::Field<bool> Foam::conformationSurfaces::wellInOutSide
 
                 surface.findNearest(sample, nearestDistSqr, info);
 
-                vector hitDir = info[0].rawPoint() - samplePts[i];
+                const vector hitDir =
-                hitDir /= mag(hitDir) + SMALL;
+                    normalised
+                    (
+                        info[0].rawPoint() - samplePts[i]
+                    );
 
                 pointIndexHit surfHit;
                 label hitSurface;

applications/utilities/mesh/generation/foamyMesh/foamyQuadMesh/CV2D.C

@@ -558,7 +558,7 @@ void Foam::CV2D::newPoints()
 
                     alignmentDirs[aA] = a + sign(dotProduct)*b;
 
-                    alignmentDirs[aA] /= mag(alignmentDirs[aA]);
+                    alignmentDirs[aA].normalise();
                 }
             }
         }
@@ -845,7 +845,7 @@ void Foam::CV2D::newPoints()
             cd0 = vector2D(cd0.x() + cd0.y(), cd0.y() - cd0.x());
 
             // Normalise the primary coordinate direction
-            cd0 /= mag(cd0);
+            cd0.normalise();
 
             // Calculate the orthogonal coordinate direction
             vector2D cd1(-cd0.y(), cd0.x());

applications/utilities/mesh/generation/foamyMesh/foamyQuadMesh/insertFeaturePoints.C

@@ -123,7 +123,7 @@ void Foam::CV2D::insertFeaturePoints()
                 vector2DField fpn = toPoint2D(feMesh.edgeNormals(edgeI));
 
                 vector2D cornerNormal = sum(fpn);
-                cornerNormal /= mag(cornerNormal);
+                cornerNormal.normalise();
 
                 if (debug)
                 {

applications/utilities/mesh/manipulation/transformPoints/transformPoints.C

@@ -280,8 +280,8 @@ int main(int argc, char *argv[])
         (
             args.lookup("rotate")()
         );
-        n1n2[0] /= mag(n1n2[0]);
+        n1n2[0].normalise();
-        n1n2[1] /= mag(n1n2[1]);
+        n1n2[1].normalise();
 
         const tensor rotT = rotationTensor(n1n2[0], n1n2[1]);
 

applications/utilities/preProcessing/engineSwirl/createFields.H

@@ -74,6 +74,6 @@ if (mag(yT) < SMALL)
 }
 
 //swirlAxis doesn't have to be of unit length.
-xT /= mag(xT);
+xT.normalise();
-yT /= mag(yT);
+yT.normalise();
-zT /= mag(zT);
+zT.normalise();

applications/utilities/surface/surfaceBooleanFeatures/surfaceBooleanFeatures.C

@@ -317,11 +317,9 @@ label calcNormalDirection
     const vector& pointOnEdge
 )
 {
-    vector cross = (normal ^ edgeDir);
+    const vector cross = normalised(normal ^ edgeDir);
-    cross /= mag(cross);
 
-    vector fC0tofE0 = faceCentre - pointOnEdge;
+    const vector fC0tofE0 = normalised(faceCentre - pointOnEdge);
-    fC0tofE0 /= mag(fC0tofE0);
 
     label nDir = ((cross & fC0tofE0) > 0.0 ? 1 : -1);
 

applications/utilities/surface/surfaceInflate/surfaceInflate.C

@@ -189,7 +189,7 @@ tmp<vectorField> calcPointNormals
         {
             if (nNormals[pointI] > 0)
             {
-                pointNormals[pointI] /= mag(pointNormals[pointI]);
+                pointNormals[pointI].normalise();
             }
         }
     }

applications/utilities/surface/surfaceSplitNonManifolds/surfaceSplitNonManifolds.C

@@ -565,16 +565,24 @@ void calcPointVecs
             if (face0I != -1)
             {
                 label v0 = triSurfaceTools::oppositeVertex(surf, face0I, edgeI);
-                vector e0 = (points[v0] - points[e.start()]) ^ eVec;
+                const vector e0 =
-                e0 /= mag(e0);
+                    normalised
+                    (
+                        (points[v0] - points[e.start()]) ^ eVec
+                    );
+
                 midVec = e0;
             }
 
             if (face1I != -1)
             {
                 label v1 = triSurfaceTools::oppositeVertex(surf, face1I, edgeI);
-                vector e1 = (points[e.start()] - points[v1]) ^ eVec;
+                const vector e1 =
-                e1 /= mag(e1);
+                    normalised
+                    (
+                        (points[e.start()] - points[v1]) ^ eVec
+                    );
+
                 midVec += e1;
             }
 
@@ -895,8 +903,7 @@ int main(int argc, char *argv[])
             {
                 scalar minLen = minEdgeLen(surf, pointi);
 
-                vector n = borderPointVec[pointi];
+                const vector n = normalised(borderPointVec[pointi]);
-                n /= mag(n);
 
                 newPoints[newPointi] = newPoints[pointi] + 0.1 * minLen * n;
             }

applications/utilities/surface/surfaceTransformPoints/surfaceTransformPoints.C

@@ -175,8 +175,8 @@ int main(int argc, char *argv[])
         (
             args.lookup("rotate")()
         );
-        n1n2[0] /= mag(n1n2[0]);
+        n1n2[0].normalise();
-        n1n2[1] /= mag(n1n2[1]);
+        n1n2[1].normalise();
 
         const tensor rotT = rotationTensor(n1n2[0], n1n2[1]);
 

src/OpenFOAM/meshes/polyMesh/polyMesh.C

@@ -94,7 +94,7 @@ void Foam::polyMesh::calcDirections() const
     {
         reduce(emptyDirVec, sumOp<vector>());
 
-        emptyDirVec /= mag(emptyDirVec);
+        emptyDirVec.normalise();
 
         for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
         {
@@ -118,7 +118,7 @@ void Foam::polyMesh::calcDirections() const
     {
         reduce(wedgeDirVec, sumOp<vector>());
 
-        wedgeDirVec /= mag(wedgeDirVec);
+        wedgeDirVec.normalise();
 
         for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
         {

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/oldCyclic/oldCyclicPolyPatch.C

@@ -308,10 +308,17 @@ void Foam::oldCyclicPolyPatch::getCentresAndAnchors
             label face0 = getConsistentRotationFace(half0Ctrs);
             label face1 = getConsistentRotationFace(half1Ctrs);
 
-            vector n0 = ((half0Ctrs[face0] - rotationCentre_) ^ rotationAxis_);
+            const vector n0 =
-            vector n1 = ((half1Ctrs[face1] - rotationCentre_) ^ -rotationAxis_);
+                normalised
-            n0 /= mag(n0) + VSMALL;
+                (
-            n1 /= mag(n1) + VSMALL;
+                    (half0Ctrs[face0] - rotationCentre_) ^ rotationAxis_
+                );
+
+            const vector n1 =
+                normalised
+                (
+                    (half1Ctrs[face1] - rotationCentre_) ^ -rotationAxis_
+                );
 
             if (debug)
             {

src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.C

@@ -86,7 +86,7 @@ void Foam::wedgePolyPatch::calcGeometry(PstreamBuffers&)
                 sign(n_.y())*(max(mag(n_.y()), 0.5) - 0.5),
                 sign(n_.z())*(max(mag(n_.z()), 0.5) - 0.5)
             );
-        centreNormal_ /= mag(centreNormal_);
+        centreNormal_.normalise();
 
         cosAngle_ = centreNormal_ & n_;
 

src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsSortEdges.C

@@ -76,8 +76,11 @@ Foam::PatchTools::sortedEdgeFaces
             {
                 if (f0[fpI] != e.start())
                 {
-                    vector faceEdgeDir = localPoints[f0[fpI]] - edgePt;
+                    const vector faceEdgeDir =
-                    faceEdgeDir /= mag(faceEdgeDir) + VSMALL;
+                        normalised
+                        (
+                            localPoints[f0[fpI]] - edgePt
+                        );
 
                     const scalar angle = e2 & faceEdgeDir;
 
@@ -91,11 +94,10 @@ Foam::PatchTools::sortedEdgeFaces
 
             // Get vector normal both to e2 and to edge from opposite vertex
             // to edge (will be x-axis of our coordinate system)
-            vector e0 = e2 ^ maxAngleEdgeDir;
+            const vector e0 = normalised(e2 ^ maxAngleEdgeDir);
-            e0 /= mag(e0) + VSMALL;
 
             // Get y-axis of coordinate system
-            vector e1 = e2 ^ e0;
+            const vector e1 = e2 ^ e0;
 
             SortableList<scalar> faceAngles(faceNbs.size());
 
@@ -115,8 +117,11 @@ Foam::PatchTools::sortedEdgeFaces
                 {
                     if (f[fpI] != e.start())
                     {
-                        vector faceEdgeDir = localPoints[f[fpI]] - edgePt;
+                        const vector faceEdgeDir =
-                        faceEdgeDir /= mag(faceEdgeDir) + VSMALL;
+                            normalised
+                            (
+                                localPoints[f[fpI]] - edgePt
+                            );
 
                         const scalar angle = e2 & faceEdgeDir;
 
@@ -128,8 +133,7 @@ Foam::PatchTools::sortedEdgeFaces
                     }
                 }
 
-                vector vec = e2 ^ maxAngleEdgeDir;
+                const vector vec = normalised(e2 ^ maxAngleEdgeDir);
-                vec /= mag(vec) + VSMALL;
 
                 faceAngles[nbI] = pseudoAngle
                 (

src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchMeshData.C

@@ -286,7 +286,7 @@ calcPointNormals() const
             curNormal += faceUnitNormals[facei];
         }
 
-        curNormal /= mag(curNormal) + VSMALL;
+        curNormal.normalise();
     }
 
     if (debug)

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C

@@ -78,8 +78,7 @@ Foam::scalar Foam::primitiveMeshTools::boundaryFaceSkewness
 {
     vector Cpf = fCtrs[facei] - ownCc;
 
-    vector normal = fAreas[facei];
+    vector normal = normalised(fAreas[facei]);
-    normal /= mag(normal) + ROOTVSMALL;
     vector d = normal*(normal & Cpf);
 
 

src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C

@@ -389,8 +389,7 @@ Foam::tmp<Foam::surfaceScalarField> Foam::fvMeshDistribute::generateTestField
     const fvMesh& mesh
 )
 {
-    vector testNormal(1, 1, 1);
+    const vector testNormal = normalised(vector::one);
-    testNormal /= mag(testNormal);
 
     tmp<surfaceScalarField> tfld
     (
@@ -440,8 +439,7 @@ void Foam::fvMeshDistribute::testField(const surfaceScalarField& fld)
 {
     const fvMesh& mesh = fld.mesh();
 
-    vector testNormal(1, 1, 1);
+    const vector testNormal = normalised(vector::one);
-    testNormal /= mag(testNormal);
 
     const surfaceVectorField n(mesh.Sf()/mesh.magSf());
 

src/dynamicMesh/meshCut/cellLooper/geomCellLooper.C

@@ -159,9 +159,7 @@ void Foam::geomCellLooper::getBase(const vector& n, vector& e0, vector& e1)
 
 
     // Use component normal to n as base vector.
-    e0 = base - nComp*n;
+    e0 = normalised(base - nComp * n);
-
-    e0 /= mag(e0) + VSMALL;
 
     e1 = n ^ e0;
 
@@ -404,8 +402,7 @@ bool Foam::geomCellLooper::cut
 
     forAll(sortedAngles, i)
     {
-        vector toCtr(loopPoints[i] - ctr);
+        const vector toCtr = normalised(loopPoints[i] - ctr);
-        toCtr /= mag(toCtr);
 
         sortedAngles[i] = pseudoAngle(e0, e1, toCtr);
     }

src/dynamicMesh/meshCut/cellLooper/topoCellLooper.C

@@ -262,8 +262,7 @@ Foam::label Foam::topoCellLooper::getAlignedNonFeatureEdge
             vector e0 = mesh().points()[e.start()] - ctr;
             vector e1 = mesh().points()[e.end()] - ctr;
 
-            vector n = e0 ^ e1;
+            const vector n = normalised(e0 ^ e1);
-            n /= mag(n);
 
             scalar cosAngle = mag(refDir & n);
 

src/dynamicMesh/meshCut/directions/directions.C

@@ -315,8 +315,7 @@ Foam::directions::directions
         vector tan2(globalDict.lookup("tan2"));
         check2D(correct2DPtr, tan2);
 
-        vector normal = tan1 ^ tan2;
+        const vector normal = normalised(tan1 ^ tan2);
-        normal /= mag(normal);
 
         Info<< "Global Coordinate system:" << endl
             << "     normal : " << normal << endl

src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C

@@ -1498,8 +1498,7 @@ bool Foam::polyMeshGeometry::checkFaceAngles
 
         const face& f = fcs[facei];
 
-        vector faceNormal = faceAreas[facei];
+        const vector faceNormal = normalised(faceAreas[facei]);
-        faceNormal /= mag(faceNormal) + VSMALL;
 
         // Get edge from f[0] to f[size-1];
         vector ePrev(p[f.first()] - p[f.last()]);
@@ -1692,20 +1691,30 @@ bool Foam::polyMeshGeometry::checkFaceTwist
 
             if (mesh.isInternalFace(facei))
             {
-                nf = cellCentres[nei[facei]] - cellCentres[own[facei]];
+                nf =
-                nf /= mag(nf) + VSMALL;
+                    normalised
+                    (
+                        cellCentres[nei[facei]]
+                      - cellCentres[own[facei]]
+                    );
             }
             else if (patches[patches.whichPatch(facei)].coupled())
             {
                 nf =
-                    neiCc[facei-mesh.nInternalFaces()]
+                    normalised
-                  - cellCentres[own[facei]];
+                    (
-                nf /= mag(nf) + VSMALL;
+                        neiCc[facei-mesh.nInternalFaces()]
+                      - cellCentres[own[facei]]
+                    );
             }
             else
             {
-                nf = faceCentres[facei] - cellCentres[own[facei]];
+                nf =
-                nf /= mag(nf) + VSMALL;
+                    normalised
+                    (
+                        faceCentres[facei]
+                      - cellCentres[own[facei]]
+                    );
             }
 
             if (nf != vector::zero)

src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C

@@ -458,8 +458,11 @@ Foam::label Foam::faceCoupleInfo::mostAlignedCutEdge
 
             eVec /= magEVec;
 
-            vector eToEndPoint(localPoints[edgeEnd] - localPoints[otherPointi]);
+            const vector eToEndPoint =
-            eToEndPoint /= mag(eToEndPoint);
+                normalised
+                (
+                    localPoints[edgeEnd] - localPoints[otherPointi]
+                );
 
             scalar cosAngle = eVec & eToEndPoint;
 

src/dynamicMesh/polyTopoChange/polyTopoChange/removePoints.C

@@ -213,11 +213,8 @@ Foam::label Foam::removePoints::countPointUsage
             label vLeft = e0.otherVertex(common);
             label vRight = e1.otherVertex(common);
 
-            vector e0Vec = points[common] - points[vLeft];
+            const vector e0Vec = normalised(points[common] - points[vLeft]);
-            e0Vec /= mag(e0Vec) + VSMALL;
+            const vector e1Vec = normalised(points[vRight] - points[common]);
-
-            vector e1Vec = points[vRight] - points[common];
-            e1Vec /= mag(e1Vec) + VSMALL;
 
             if ((e0Vec & e1Vec) > minCos)
             {

src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatchCutFaces.C

@@ -220,10 +220,9 @@ void Foam::enrichedPatch::calcCutFaces() const
                 // Get the vector along the edge and the right vector
                 vector ahead = curPoint - lp[prevPointLabel];
                 ahead -= normal*(normal & ahead);
-                ahead /= mag(ahead);
+                ahead.normalise();
 
-                vector right = normal ^ ahead;
+                const vector right = normalised(normal ^ ahead);
-                right /= mag(right);
 
                 // Pout<< "normal: " << normal
                 //     << " ahead: " << ahead

src/dynamicMesh/slidingInterface/slidingInterfaceProjectPoints.C

@@ -920,15 +920,16 @@ bool Foam::slidingInterface::projectPoints() const
             // projected-master-to-edge distance is used, to avoid
             // problems with badly defined master planes.  HJ,
             // 17/Oct/2004
-            vector edgeNormalInPlane =
+            const vector edgeNormalInPlane =
-                edgeVec
+                normalised
-              ^ (
+                (
-                    slavePointNormals[curEdge.start()]
+                    edgeVec
-                  + slavePointNormals[curEdge.end()]
+                  ^ (
+                        slavePointNormals[curEdge.start()]
+                      + slavePointNormals[curEdge.end()]
+                    )
                 );
 
-            edgeNormalInPlane /= mag(edgeNormalInPlane);
-
             for (const label cmp : curMasterPoints)
             {
                 // Skip the current point if the edge start or end has

src/finiteArea/faMesh/faMeshDemandDrivenData.C

@@ -1281,10 +1281,11 @@ void Foam::faMesh::calcPointAreaNormalsByQuadricsFit() const
 
         // Transform points
         const vector& origin = points[curPoint];
-        vector axis(result[curPoint]/mag(result[curPoint]));
+        const vector axis = normalised(result[curPoint]);
         vector dir(allPoints[0] - points[curPoint]);
         dir -= axis*(axis&dir);
-        dir /= mag(dir);
+        dir.normalise();
+
         coordinateSystem cs("cs", origin, axis, dir);
 
         forAll(allPoints, pI)
@@ -1563,10 +1564,11 @@ void Foam::faMesh::calcPointAreaNormalsByQuadricsFit() const
 
                 // Transform points
                 const vector& origin = points[curPoint];
-                vector axis(result[curPoint]/mag(result[curPoint]));
+                const vector axis = normalised(result[curPoint]);
                 vector dir(allPoints[0] - points[curPoint]);
                 dir -= axis*(axis&dir);
-                dir /= mag(dir);
+                dir.normalise();
+
                 const coordinateSystem cs("cs", origin, axis, dir);
 
                 scalarField W(allPoints.size(), 1.0);
@@ -1734,10 +1736,11 @@ void Foam::faMesh::calcPointAreaNormalsByQuadricsFit() const
 
                 // Transform points
                 const vector& origin = points[curPoint];
-                vector axis(result[curPoint]/mag(result[curPoint]));
+                const vector axis = normalised(result[curPoint]);
                 vector dir(allPoints[0] - points[curPoint]);
                 dir -= axis*(axis&dir);
-                dir /= mag(dir);
+                dir.normalise();
+
                 coordinateSystem cs("cs", origin, axis, dir);
 
                 scalarField W(allPoints.size(), 1.0);

src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C

@@ -404,19 +404,16 @@ void Foam::edgeInterpolation::makeDeltaCoeffs() const
     forAll(owner, edgeI)
     {
         // Edge normal - area normal
-        vector edgeNormal = lengths[edgeI]^edges[edgeI].vec(points);
+        vector edgeNormal =
-        edgeNormal /= mag(edgeNormal);
+            normalised(lengths[edgeI] ^ edges[edgeI].vec(points));
-
 
         // Unit delta vector
         vector unitDelta =
             faceCentres[neighbour[edgeI]]
           - faceCentres[owner[edgeI]];
 
-        unitDelta -=
+        unitDelta -= edgeNormal*(edgeNormal & unitDelta);
-            edgeNormal*(edgeNormal&unitDelta);
+        unitDelta.normalise();
-
-        unitDelta /= mag(unitDelta);
 
 
         // Calc PN arc length
@@ -447,7 +444,7 @@ void Foam::edgeInterpolation::makeDeltaCoeffs() const
 
 
         // Edge normal - area tangent
-        edgeNormal = lengths[edgeI]/mag(lengths[edgeI]);
+        edgeNormal = normalised(lengths[edgeI]);
 
         // Delta coefficient as per definition
 //         dc[edgeI] = 1.0/(lPN*(unitDelta & edgeNormal));
@@ -496,7 +493,6 @@ void Foam::edgeInterpolation::makeCorrectionVectors() const
     const labelUList& neighbour = mesh().neighbour();
 
     const edgeVectorField& lengths = mesh().Le();
-    const edgeScalarField& magLengths = mesh().magLe();
 
     const edgeList& edges = mesh().edges();
     const pointField& points = mesh().points();
@@ -508,8 +504,8 @@ void Foam::edgeInterpolation::makeCorrectionVectors() const
     forAll(owner, edgeI)
     {
         // Edge normal - area normal
-        vector edgeNormal = lengths[edgeI] ^ edges[edgeI].vec(points);
+        vector edgeNormal =
-        edgeNormal /= mag(edgeNormal);
+            normalised(lengths[edgeI] ^ edges[edgeI].vec(points));
 
         // Unit delta vector
         vector unitDelta =
@@ -517,11 +513,10 @@ void Foam::edgeInterpolation::makeCorrectionVectors() const
           - faceCentres[owner[edgeI]];
 
         unitDelta -= edgeNormal*(edgeNormal & unitDelta);
-
+        unitDelta.normalise();
-        unitDelta /= mag(unitDelta);
 
         // Edge normal - area tangent
-        edgeNormal = lengths[edgeI]/magLengths[edgeI];
+        edgeNormal = normalised(lengths[edgeI]);
 
         // Delta coeffs
         deltaCoeffs[edgeI] = 1.0/(unitDelta & edgeNormal);

src/finiteVolume/cfdTools/general/SRF/SRFModel/SRFModel/SRFModel.C

@@ -59,14 +59,11 @@ Foam::SRF::SRFModel::SRFModel
     ),
     Urel_(Urel),
     mesh_(Urel_.mesh()),
-    origin_("origin", dimLength, lookup("origin")),
+    origin_("origin", dimLength, get<vector>("origin")),
-    axis_(lookup("axis")),
+    axis_(normalised(get<vector>("axis"))),
     SRFModelCoeffs_(optionalSubDict(type + "Coeffs")),
     omega_(dimensionedVector("omega", dimless/dimTime, Zero))
-{
+{}
-    // Normalise the axis
-    axis_ /= mag(axis_);
-}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
@@ -82,21 +79,19 @@ bool Foam::SRF::SRFModel::read()
     if (regIOobject::read())
     {
         // Re-read origin
-        lookup("origin") >> origin_;
+        readEntry("origin", origin_);
 
         // Re-read axis
-        lookup("axis") >> axis_;
+        readEntry("axis", axis_);
-        axis_ /= mag(axis_);
+        axis_.normalise();
 
         // Re-read sub-model coeffs
         SRFModelCoeffs_ = optionalSubDict(type() + "Coeffs");
 
         return true;
     }
-    else
+
-    {
+    return false;
-        return false;
-    }
 }
 
 

src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutCell/isoCutCell.C

@@ -267,12 +267,11 @@ void Foam::isoCutCell::calcIsoFacePointsFromEdges()
 
     // Defining local coordinates with zhat along isoface normal and xhat from
     // isoface centre to first point in isoFaceEdges_
-    const vector zhat = isoFaceArea_/mag(isoFaceArea_);
+    const vector zhat = normalised(isoFaceArea_);
     vector xhat = isoFaceEdges_[0][0] - isoFaceCentre_;
     xhat = (xhat - (xhat & zhat)*zhat);
-    xhat /= mag(xhat);
+    xhat.normalise();
-    vector yhat = zhat^xhat;
+    vector yhat = normalised(zhat ^ xhat);
-    yhat /= mag(yhat);
 
     DebugPout << "Calculated local coordinates" << endl;
 

src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutFace/isoCutFace.C

@@ -735,12 +735,12 @@ void Foam::isoCutFace::quadAreaCoeffs
         if (Bx > 10*SMALL)
         {
             // If |AB| > 0 ABCD we use AB to define xhat
-            xhat = (B - A)/mag(B - A);
+            xhat = normalised(B - A);
         }
         else if (mag(C - D) > 10*SMALL)
         {
             // If |AB| ~ 0 ABCD is a triangle ACD and we use CD for xhat
-            xhat = (C - D)/mag(C - D);
+            xhat = normalised(C - D);
         }
         else
         {

src/finiteVolume/interpolation/interpolation/interpolationCellPointFace/findCellPointFaceTet.H

@@ -81,10 +81,11 @@ bool interpolationCellPointFace<Type>::findTet
             referencePoint = tetPoints[p1];
 
             faceNormal =
-                (tetPoints[p3] - tetPoints[p1])
+                normalised
-              ^ (tetPoints[p2] - tetPoints[p1]);
+                (
-
+                    (tetPoints[p3] - tetPoints[p1])
-            faceNormal /= mag(faceNormal);
+                  ^ (tetPoints[p2] - tetPoints[p1])
+                );
 
             // correct normal to point into the tet
             vector v0 = tetPoints[n] - referencePoint;

src/finiteVolume/interpolation/interpolation/interpolationCellPointFace/findCellPointFaceTriangle.H

@@ -57,8 +57,7 @@ bool interpolationCellPointFace<Type>::findTriangle
         // calculate edge normal (pointing inwards)
         for (label i=0; i<3; i++)
         {
-            normal[i] = triangleFaceNormal ^ edge[i];
+            normal[i] = normalised(triangleFaceNormal ^ edge[i]);
-            normal[i] /= mag(normal[i]) + VSMALL;
         }
 
         // check if position is inside triangle

src/finiteVolume/interpolation/interpolation/interpolationCellPointFace/interpolationCellPointFace.C

@@ -88,12 +88,9 @@ Type Foam::interpolationCellPointFace<Type>::interpolate
         label closestFace = -1;
         scalar minDistance = GREAT;
 
-        forAll(cellFaces, facei)
+        for (const label nFace : cellFaces)
         {
-            label nFace = cellFaces[facei];
+            const vector normal = normalised(this->pMeshFaceAreas_[nFace]);
-
-            vector normal = this->pMeshFaceAreas_[nFace];
-            normal /= mag(normal);
 
             const vector& faceCentreTmp = this->pMeshFaceCentres_[nFace];
 

src/finiteVolume/interpolation/interpolation/interpolationPointMVC/pointMVCWeight.C

@@ -156,10 +156,8 @@ void Foam::pointMVCWeight::calcWeights
             label jPlus1 = f.fcIndex(j);
             //Pout<< "    uj:" << u[j] << " ujPlus1:" << u[jPlus1] << endl;
 
-            vector n0 = u[j]^v;
+            const vector n0 = normalised(u[j] ^ v);
-            n0 /= mag(n0);
+            const vector n1 = normalised(u[jPlus1] ^ v);
-            vector n1 = u[jPlus1]^v;
-            n1 /= mag(n1);
 
             scalar l = min(mag(n0 - n1), 2.0);
             //Pout<< "    l:" << l << endl;

src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.C

@@ -76,8 +76,7 @@ void Foam::FitData<FitDataType, ExtendedStencil, Polynomial>::findFaceDirs
 {
     const fvMesh& mesh = this->mesh();
 
-    idir = mesh.faceAreas()[facei];
+    idir = normalised(mesh.faceAreas()[facei]);
-    idir /= mag(idir);
 
     #ifndef SPHERICAL_GEOMETRY
     if (mesh.nGeometricD() <= 2) // find the normal direction

src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.C

@@ -116,8 +116,7 @@ void Foam::DSMCParcel<ParcelType>::hitWallPatch
 
     scalar m = constProps.mass();
 
-    vector nw = wpp.faceAreas()[wppLocalFace];
+    const vector nw = normalised(wpp.faceAreas()[wppLocalFace]);
-    nw /= mag(nw);
 
     scalar U_dot_nw = U_ & nw;
 

src/lagrangian/DSMC/submodels/InflowBoundaryModel/FreeStream/FreeStream.C

@@ -265,18 +265,15 @@ void Foam::FreeStream<CloudType>::inflow()
 
             // Normal unit vector *negative* so normal is pointing into the
             // domain
-            vector n = patch.faceAreas()[pFI];
+            const vector n = -normalised(patch.faceAreas()[pFI]);
-            n /= -mag(n);
 
             // Wall tangential unit vector. Use the direction between the
             // face centre and the first vertex in the list
-            vector t1 = fC - (mesh.points()[f[0]]);
+            const vector t1 = normalised(fC - (mesh.points()[f[0]]));
-            t1 /= mag(t1);
 
             // Other tangential unit vector.  Rescaling in case face is not
             // flat and n and t1 aren't perfectly orthogonal
-            vector t2 = n^t1;
+            const vector t2 = normalised(n ^ t1);
-            t2 /= mag(t2);
 
             scalar faceTemperature = boundaryT[patchi][pFI];
 

src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.C

@@ -167,7 +167,7 @@ void Foam::ParticleCollector<CloudType>::initConcentricCircles()
     {
         refDir = this->coeffDict().lookup("refDir");
         refDir -= normal_[0]*(normal_[0] & refDir);
-        refDir /= mag(refDir);
+        refDir.normalise();
     }
     else
     {
@@ -581,8 +581,7 @@ Foam::ParticleCollector<CloudType>::ParticleCollector
         forAll(polygons, polyI)
         {
             polygons[polyI] = polygonAndNormal[polyI].first();
-            normal_[polyI] = polygonAndNormal[polyI].second();
+            normal_[polyI] = normalised(polygonAndNormal[polyI].second());
-            normal_[polyI] /= mag(normal_[polyI]) + ROOTVSMALL;
         }
 
         initPolygons(polygons);

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ConeInjection/ConeInjection.C

@@ -104,8 +104,7 @@ Foam::ConeInjection<CloudType>::ConeInjection
     forAll(positionAxis_, i)
     {
         vector& axis = positionAxis_[i].second();
-
+        axis.normalise();
-        axis /= mag(axis);
 
         vector tangent = Zero;
         scalar magTangent = 0.0;
@@ -277,7 +276,7 @@ void Foam::ConeInjection<CloudType>::setProperties
     vector normal = alpha*(tanVec1_[i]*cos(beta) + tanVec2_[i]*sin(beta));
     vector dirVec = dcorr*positionAxis_[i].second();
     dirVec += normal;
-    dirVec /= mag(dirVec);
+    dirVec.normalise();
 
     parcel.U() = Umag_.value(t)*dirVec;
 

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ConeNozzleInjection/ConeNozzleInjection.C

@@ -206,7 +206,7 @@ Foam::ConeNozzleInjection<CloudType>::ConeNozzleInjection
     Random& rndGen = this->owner().rndGen();
 
     // Normalise direction vector
-    direction_ /= mag(direction_);
+    direction_.normalise();
 
     // Determine direction vectors tangential to direction
     vector tangent = Zero;
@@ -436,7 +436,7 @@ void Foam::ConeNozzleInjection<CloudType>::setProperties
     vector normal = alpha*normal_;
     vector dirVec = dcorr*direction_;
     dirVec += normal;
-    dirVec /= mag(dirVec);
+    dirVec.normalise();
 
     switch (flowType_)
     {

src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H

@@ -90,9 +90,11 @@ inline Foam::molecule::constantProperties::constantProperties
 
         Info<< nl << "Linear molecule." << endl;
 
-        vector dir = siteReferencePositions_[1] - siteReferencePositions_[0];
+        const vector dir =
-
+            normalised
-        dir /= mag(dir);
+            (
+                siteReferencePositions_[1] - siteReferencePositions_[0]
+            );
 
         tensor Q = rotationTensor(dir, vector(1,0,0));
 
@@ -334,9 +336,11 @@ inline bool Foam::molecule::constantProperties::linearMoleculeTest() const
         return true;
     }
 
-    vector refDir = siteReferencePositions_[1] - siteReferencePositions_[0];
+    const vector refDir =
-
+        normalised
-    refDir /= mag(refDir);
+        (
+            siteReferencePositions_[1] - siteReferencePositions_[0]
+        );
 
     for
     (
@@ -345,9 +349,11 @@ inline bool Foam::molecule::constantProperties::linearMoleculeTest() const
         i++
     )
     {
-        vector dir = siteReferencePositions_[i] - siteReferencePositions_[i-1];
+        const vector dir =
-
+            normalised
-        dir /= mag(dir);
+            (
+                siteReferencePositions_[i] - siteReferencePositions_[i-1]
+            );
 
         if (mag(refDir & dir) < 1 - SMALL)
         {

src/lagrangian/spray/submodels/AtomizationModel/LISAAtomization/LISAAtomization.C

@@ -43,7 +43,7 @@ Foam::LISAAtomization<CloudType>::LISAAtomization
     SMDCalcMethod_(this->coeffDict().lookup("SMDCalculationMethod"))
 {
     // Note: Would be good if this could be picked up from the injector
-    injectorDirection_ /= mag(injectorDirection_);
+    injectorDirection_.normalise();
 
     if (SMDCalcMethod_ == "method1")
     {

src/mesh/extrudeModel/linearDirection/linearDirection.C

@@ -45,11 +45,9 @@ addToRunTimeSelectionTable(extrudeModel, linearDirection, dictionary);
 linearDirection::linearDirection(const dictionary& dict)
 :
     extrudeModel(typeName, dict),
-    direction_(coeffDict_.lookup("direction")),
+    direction_(coeffDict_.get<vector>("direction").normalise()),
-    thickness_(readScalar(coeffDict_.lookup("thickness")))
+    thickness_(coeffDict_.get<scalar>("thickness"))
 {
-    direction_ /= mag(direction_);
-
     if (thickness_ <= 0)
     {
         FatalErrorInFunction

src/mesh/snappyHexMesh/externalDisplacementMeshMover/fieldSmoother/fieldSmoother.C

@@ -129,8 +129,7 @@ void Foam::fieldSmoother::smoothNormals
             {
                 //full smoothing neighbours + point value
                 average[pointI] = 0.5*(normals[pointI]+average[pointI]);
-                normals[pointI] = average[pointI];
+                normals[pointI] = normalised(average[pointI]);
-                normals[pointI] /= mag(normals[pointI]) + VSMALL;
             }
         }
     }
@@ -196,8 +195,7 @@ void Foam::fieldSmoother::smoothPatchNormals
         {
             // full smoothing neighbours + point value
             average[pointI] = 0.5*(normals[pointI]+average[pointI]);
-            normals[pointI] = average[pointI];
+            normals[pointI] = normalised(average[pointI]);
-            normals[pointI] /= mag(normals[pointI]) + VSMALL;
         }
     }
 }

src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C

@@ -1441,12 +1441,14 @@ void Foam::medialAxisMeshMover::calculateDisplacement
 
             //- Option 2: Look at component in the direction
             //  of nearest (medial axis or static) point.
-            vector n =
+            const vector n = normalised(patchDisp[patchPointI]);
-                patchDisp[patchPointI]
+            const vector mVec =
-              / (mag(patchDisp[patchPointI]) + VSMALL);
+                normalised
-            vector mVec = medialVec_[pointI]-mesh().points()[pointI];
+                (
-            mVec /= mag(mVec)+VSMALL;
+                    medialVec_[pointI] - mesh().points()[pointI]
-            thicknessRatio *= (n&mVec);
+                );
+
+            thicknessRatio *= (n & mVec);
 
             if (thicknessRatio > maxThicknessToMedialRatio)
             {

src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C

@@ -172,8 +172,7 @@ void Foam::meshRefinement::calcNeighbourData
             forAll(faceCells, i)
             {
                 // Extrapolate the face centre.
-                vector fn = faceAreas[i];
+                const vector fn = normalised(faceAreas[i]);
-                fn /= mag(fn)+VSMALL;
 
                 label own = faceCells[i];
                 label ownLevel = cellLevel[own];
@@ -2344,7 +2343,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMeshRegions
     label nRemove = findRegions
     (
         mesh_,
-        mergeDistance_*vector(1,1,1),   // perturbVec
+        mergeDistance_ * vector::one,   // perturbVec
         locationsInMesh,
         locationsOutsideMesh,
         cellRegion.nRegions(),

src/mesh/snappyHexMesh/meshRefinement/meshRefinementBaffles.C

@@ -1649,7 +1649,7 @@ void Foam::meshRefinement::findCellZoneInsideWalk
         (
             mesh_,
             cellRegion,
-            mergeDistance_*vector(1,1,1),
+            mergeDistance_ * vector::one,
             insidePoint
         );
 
@@ -1916,7 +1916,7 @@ void Foam::meshRefinement::findCellZoneTopo
         (
             mesh_,
             cellRegion,
-            mergeDistance_*vector(1,1,1),
+            mergeDistance_ * vector::one,
             keepPoint
         );
 
@@ -3875,7 +3875,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMesh
     findRegions
     (
         mesh_,
-        mergeDistance_*vector(1,1,1),   // perturbVec
+        mergeDistance_ * vector::one,   // perturbVec
         locationsInMesh,
         locationsOutsideMesh,
         cellRegion.nRegions(),

src/mesh/snappyHexMesh/meshRefinement/meshRefinementProblemCells.C

@@ -243,8 +243,7 @@ Foam::Map<Foam::label> Foam::meshRefinement::findEdgeConnectedProblemCells
     {
         label facei = candidateFaces[i];
 
-        vector n = mesh_.faceAreas()[facei];
+        const vector n = normalised(mesh_.faceAreas()[facei]);
-        n /= mag(n);
 
         label region = surfaces_.globalRegion
         (

src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.C

@@ -1101,9 +1101,7 @@ void Foam::snappySnapDriver::detectNearSurfaces
     //// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     //
     //{
-    //    const scalar cos45 = Foam::cos(degToRad(45.0));
+    //    const vector n = normalised(vector::one);
-    //    vector n(cos45, cos45, cos45);
-    //    n /= mag(n);
     //
     //    pointField start(14*pp.nPoints());
     //    pointField end(start.size());

src/meshTools/cellFeatures/cellFeatures.C

@@ -124,15 +124,11 @@ bool Foam::cellFeatures::isCellFeatureEdge
 
     // Check the angle between them by comparing the face normals.
 
-    vector n0 = mesh_.faceAreas()[face0];
+    const vector n0 = normalised(mesh_.faceAreas()[face0]);
-    n0 /= mag(n0);
+    const vector n1 = normalised(mesh_.faceAreas()[face1]);
-
-    vector n1 = mesh_.faceAreas()[face1];
-    n1 /= mag(n1);
 
     scalar cosAngle = n0 & n1;
 
-
     const edge& e = mesh_.edges()[edgeI];
 
     const face& f0 = mesh_.faces()[face0];
@@ -167,10 +163,8 @@ bool Foam::cellFeatures::isCellFeatureEdge
     {
         return true;
     }
-    else
+
-    {
+    return false;
-        return false;
-    }
 }
 
 

src/meshTools/indexedOctree/treeDataFace.C

@@ -250,14 +250,11 @@ Foam::volumeType Foam::treeDataFace::getVolumeType
 
             vector pointNormal(Zero);
 
-            forAll(pFaces, i)
+            for (const label facei : pFaces)
             {
-                if (isTreeFace_.test(pFaces[i]))
+                if (isTreeFace_.test(facei))
                 {
-                    vector n = mesh_.faceAreas()[pFaces[i]];
+                    pointNormal += normalised(mesh_.faceAreas()[facei]);
-                    n /= mag(n) + VSMALL;
-
-                    pointNormal += n;
                 }
             }
 
@@ -319,14 +316,11 @@ Foam::volumeType Foam::treeDataFace::getVolumeType
 
             vector edgeNormal(Zero);
 
-            forAll(eFaces, i)
+            for (const label facei : eFaces)
             {
-                if (isTreeFace_.test(eFaces[i]))
+                if (isTreeFace_.test(facei))
                 {
-                    vector n = mesh_.faceAreas()[eFaces[i]];
+                    edgeNormal += normalised(mesh_.faceAreas()[facei]);
-                    n /= mag(n) + VSMALL;
-
-                    edgeNormal += n;
                 }
             }
 
@@ -369,11 +363,8 @@ Foam::volumeType Foam::treeDataFace::getVolumeType
             vector ePrev = points[f[f.rcIndex(fp)]] - fc;
             vector eNext = points[f[f.fcIndex(fp)]] - fc;
 
-            vector nLeft = ePrev ^ e;
+            vector nLeft = normalised(ePrev ^ e);
-            nLeft /= mag(nLeft) + VSMALL;
+            vector nRight = normalised(e ^ eNext);
-
-            vector nRight = e ^ eNext;
-            nRight /= mag(nRight) + VSMALL;
 
             if (debug & 2)
             {

src/meshTools/indexedOctree/treeDataPrimitivePatch.C

@@ -310,11 +310,8 @@ Foam::volumeType Foam::treeDataPrimitivePatch<PatchType>::getVolumeType
             vector ePrev = points[f[f.rcIndex(fp)]] - fc;
             vector eNext = points[f[f.fcIndex(fp)]] - fc;
 
-            vector nLeft = ePrev ^ e;
+            vector nLeft = normalised(ePrev ^ e);
-            nLeft /= mag(nLeft) + VSMALL;
+            vector nRight = normalised(e ^ eNext);
-
-            vector nRight = e ^ eNext;
-            nRight /= mag(nRight) + VSMALL;
 
             if (debug & 2)
             {

src/meshTools/meshSearch/meshSearch.C

@@ -810,8 +810,7 @@ Foam::List<Foam::pointIndexHit> Foam::meshSearch::intersections
 {
     DynamicList<pointIndexHit> hits;
 
-    vector edgeVec = pEnd - pStart;
+    const vector edgeVec = normalised(pEnd - pStart);
-    edgeVec /= mag(edgeVec);
 
     point pt = pStart;
 

src/meshTools/meshTools/meshTools.C

@@ -800,7 +800,7 @@ Foam::vector Foam::meshTools::edgeToCutDir
         edgeI = meshTools::walkFace(mesh, facei, edgeI, vertI, 2);
     }
 
-    avgVec /= mag(avgVec) + VSMALL;
+    avgVec.normalise();
 
     return avgVec;
 }

src/meshTools/primitiveMeshGeometry/primitiveMeshGeometry.C

@@ -579,8 +579,7 @@ bool Foam::primitiveMeshGeometry::checkFaceSkewness
             // Boundary faces: consider them to have only skewness error.
             // (i.e. treat as if mirror cell on other side)
 
-            vector faceNormal = faceAreas[facei];
+            const vector faceNormal = normalised(faceAreas[facei]);
-            faceNormal /= mag(faceNormal) + VSMALL;
 
             vector dOwn = faceCentres[facei] - cellCentres[own[facei]];
 
@@ -767,8 +766,7 @@ bool Foam::primitiveMeshGeometry::checkFaceAngles
 
         const face& f = fcs[facei];
 
-        vector faceNormal = faceAreas[facei];
+        const vector faceNormal = normalised(faceAreas[facei]);
-        faceNormal /= mag(faceNormal) + VSMALL;
 
         // Get edge from f[0] to f[size-1];
         vector ePrev(p[f.first()] - p[f.last()]);

src/meshTools/searchableSurfaces/searchableCone/searchableCone.C

@@ -107,15 +107,8 @@ void Foam::searchableCone::findNearestAndNormal
     // Remove the parallel component and normalise
     v -= parallel*unitDir_;
 
-    scalar magV = mag(v);
+    const scalar magV = mag(v);
-    if (magV < ROOTVSMALL)
+    v.normalise();
-    {
-       v = Zero;
-    }
-    else
-    {
-       v /= magV;
-    }
 
     // Nearest and normal on disk at point1
     point disk1Point(point1_ + min(max(magV, innerRadius1_), radius1_)*v);
@@ -143,31 +136,30 @@ void Foam::searchableCone::findNearestAndNormal
         p1 /= mag(p1);
 
         // Find vector along the two end of cone
-        vector b(projPt2 - projPt1);
+        const vector b = normalised(projPt2 - projPt1);
-        scalar magS = mag(b);
-        b /= magS;
 
         // Find the vector along sample pt and pt at one end of cone
-        vector a(sample - projPt1);
+        vector a = (sample - projPt1);
 
         if (mag(a) <= ROOTVSMALL)
         {
             // Exception: sample on disk1. Redo with projPt2.
-            vector a(sample - projPt2);
+            a = (sample - projPt2);
+
             // Find normal unitvector
-            nearCone = (a & b)*b+projPt2;
+            nearCone = (a & b)*b + projPt2;
+
             vector b1 = (p1 & b)*b;
-            normalCone = p1 - b1;
+            normalCone = normalised(p1 - b1);
-            normalCone /= mag(normalCone);
         }
         else
         {
-            // Find neartest point on cone surface
+            // Find nearest point on cone surface
-            nearCone = (a & b)*b+projPt1;
+            nearCone = (a & b)*b + projPt1;
+
             // Find projection along surface of cone
             vector b1 = (p1 & b)*b;
-            normalCone = p1 - b1;
+            normalCone = normalised(p1 - b1);
-            normalCone /= mag(normalCone);
         }
 
         if (innerRadius1_ > 0 || innerRadius2_ > 0)
@@ -176,13 +168,10 @@ void Foam::searchableCone::findNearestAndNormal
             point iCprojPt1 = point1_+ innerRadius1_*v;
             point iCprojPt2 = point2_+ innerRadius2_*v;
 
-            vector iCp1 = (iCprojPt1 - point1_);
+            const vector iCp1 = normalised(iCprojPt1 - point1_);
-            iCp1 /= mag(iCp1);
 
             // Find vector along the two end of cone
-            vector iCb(iCprojPt2 - iCprojPt1);
+            const vector iCb = normalised(iCprojPt2 - iCprojPt1);
-            magS = mag(iCb);
-            iCb /= magS;
 
 
             // Find the vector along sample pt and pt at one end of conde
@@ -190,22 +179,22 @@ void Foam::searchableCone::findNearestAndNormal
 
             if (mag(iCa) <= ROOTVSMALL)
             {
-                vector iCa(sample - iCprojPt2);
+                iCa = (sample - iCprojPt2);
 
                 // Find normal unitvector
                 iCnearCone = (iCa & iCb)*iCb+iCprojPt2;
+
                 vector b1 = (iCp1 & iCb)*iCb;
-                iCnormalCone = iCp1 - b1;
+                iCnormalCone = normalised(iCp1 - b1);
-                iCnormalCone /= mag(iCnormalCone);
             }
             else
             {
                 // Find nearest point on cone surface
                 iCnearCone = (iCa & iCb)*iCb+iCprojPt1;
+
                 // Find projection along surface of cone
                 vector b1 = (iCp1 & iCb)*iCb;
-                iCnormalCone = iCp1 - b1;
+                iCnormalCone = normalised(iCp1 - b1);
-                iCnormalCone /= mag(iCnormalCone);
             }
         }
     }
@@ -234,8 +223,9 @@ void Foam::searchableCone::findNearestAndNormal
             scalar para = (v1 & unitDir_);
             // Remove the parallel component and normalise
             v1 -= para*unitDir_;
-            scalar magV1 = mag(v1);
+            const scalar magV1 = mag(v1);
             v1 = v1/magV1;
+
             if (para < 0.0 && magV1 >= radius1_)
             {
                 // Near point 1. Set point to intersection of disk and cone.
@@ -281,7 +271,8 @@ void Foam::searchableCone::findNearestAndNormal
             scalar para = (v1 & unitDir_);
             // Remove the parallel component and normalise
             v1 -= para*unitDir_;
-            scalar magV1 = mag(v1);
+
+            const scalar magV1 = mag(v1);
             v1 = v1/magV1;
 
             if (para < 0.0 && magV1 >= innerRadius1_)
@@ -382,7 +373,7 @@ void Foam::searchableCone::findLineAll
 
     {
         // Find dot product: mag(s)>VSMALL suggest that it is greater
-        scalar s = (V&unitDir_);
+        scalar s = (V & unitDir_);
         if (mag(s) > VSMALL)
         {
             tPoint1 = -s1/s;
@@ -471,8 +462,7 @@ void Foam::searchableCone::findLineAll
     else
     {
         vector va = cone.unitDir_;
-        vector v1(end-start);
+        vector v1 = normalised(end-start);
-        v1 = v1/mag(v1);
         scalar p  = (va&v1);
         vector a1 = (v1-p*va);
 

src/meshTools/searchableSurfaces/searchableExtrudedCircle/searchableExtrudedCircle.C

@@ -89,8 +89,8 @@ Foam::searchableExtrudedCircle::searchableExtrudedCircle
     vector halfSpan(0.5*bounds().span());
     point ctr(bounds().midpoint());
 
-    bounds().min() = ctr - mag(halfSpan)*vector(1, 1, 1);
+    bounds().min() = ctr - mag(halfSpan) * vector::one;
-    bounds().max() = ctr + mag(halfSpan)*vector(1, 1, 1);
+    bounds().max() = ctr + mag(halfSpan) * vector::one;
 
     // Calculate bb of all points
     treeBoundBox bb(bounds());
@@ -182,8 +182,9 @@ void Foam::searchableExtrudedCircle::findNearest
 
         if (info[i].hit())
         {
-            vector d(samples[i]-info[i].hitPoint());
+            const vector d = normalised(samples[i] - info[i].hitPoint());
-            info[i].setPoint(info[i].hitPoint() + d/mag(d)*radius_);
+
+            info[i].setPoint(info[i].hitPoint() + d*radius_);
         }
     }
 }
@@ -360,7 +361,8 @@ void Foam::searchableExtrudedCircle::findParametricNearest
     {
         radialStart = start-curvePoints[0];
         radialStart -= (radialStart&axialVecs[0])*axialVecs[0];
-        radialStart /= mag(radialStart);
+        radialStart.normalise();
+
         qStart = quaternion(radialStart, 0.0);
 
         info[0] = pointIndexHit(true, start, 0);
@@ -370,11 +372,12 @@ void Foam::searchableExtrudedCircle::findParametricNearest
     {
         vector radialEnd(end-curvePoints.last());
         radialEnd -= (radialEnd&axialVecs.last())*axialVecs.last();
-        radialEnd /= mag(radialEnd);
+        radialEnd.normalise();
 
         vector projectedEnd = radialEnd;
         projectedEnd -= (projectedEnd&axialVecs[0])*axialVecs[0];
-        projectedEnd /= mag(projectedEnd);
+        projectedEnd.normalise();
+
         qProjectedEnd = quaternion(projectedEnd, 0.0);
 
         info.last() = pointIndexHit(true, end, 0);
@@ -385,8 +388,8 @@ void Foam::searchableExtrudedCircle::findParametricNearest
         quaternion q(slerp(qStart, qProjectedEnd, lambdas[i]));
         vector radialDir(q.transform(radialStart));
 
-        radialDir -= (radialDir&axialVecs[i])*axialVecs.last();
+        radialDir -= (radialDir & axialVecs[i]) * axialVecs.last();
-        radialDir /= mag(radialDir);
+        radialDir.normalise();
 
         info[i] = pointIndexHit(true, curvePoints[i]+radius_*radialDir, 0);
     }
@@ -434,8 +437,8 @@ void Foam::searchableExtrudedCircle::getNormal
             // Subtract axial direction
             const vector axialVec = edges[curvePt.index()].unitVec(points);
 
-            normal[i] -= (normal[i]&axialVec)*axialVec;
+            normal[i] -= (normal[i] & axialVec) * axialVec;
-            normal[i] /= mag(normal[i]);
+            normal[i].normalise();
         }
     }
 }

src/meshTools/searchableSurfaces/searchableSphere/searchableSphere.C

@@ -336,9 +336,7 @@ void Foam::searchableSphere::getNormal
     {
         if (info[i].hit())
         {
-            normal[i] = info[i].hitPoint() - centre_;
+            normal[i] = normalised(info[i].hitPoint() - centre_);
-
-            normal[i] /= mag(normal[i])+VSMALL;
         }
         else
         {

src/meshTools/sets/faceSources/normalToFace/normalToFace.C

@@ -51,7 +51,7 @@ Foam::topoSetSource::addToUsageTable Foam::normalToFace::usage_
 
 void Foam::normalToFace::setNormal()
 {
-    normal_ /= mag(normal_) + VSMALL;
+    normal_.normalise();
 
     Info<< "    normalToFace : Normalized vector to " << normal_ << endl;
 
@@ -116,8 +116,7 @@ void Foam::normalToFace::applyToSet
 
         forAll(mesh_.faceAreas(), facei)
         {
-            vector n = mesh_.faceAreas()[facei];
+            const vector n = normalised(mesh_.faceAreas()[facei]);
-            n /= mag(n) + VSMALL;
 
             if (mag(1 - (n & normal_)) < tol_)
             {
@@ -136,8 +135,7 @@ void Foam::normalToFace::applyToSet
         {
             const label facei = iter.key();
 
-            vector n = mesh_.faceAreas()[facei];
+            const vector n = normalised(mesh_.faceAreas()[facei]);
-            n /= mag(n) + VSMALL;
 
             if (mag(1 - (n & normal_)) < tol_)
             {

src/meshTools/triSurface/booleanOps/intersectedSurface/intersectedSurface.C

@@ -433,11 +433,14 @@ Foam::label Foam::intersectedSurface::nextEdge
     const edge& prevE = edges[prevEdgei];
 
     // x-axis of coordinate system
-    vector e0 = n ^ (points[prevE.otherVertex(prevVerti)] - points[prevVerti]);
+    const vector e0 =
-    e0 /= mag(e0) + VSMALL;
+        normalised
+        (
+            n ^ (points[prevE.otherVertex(prevVerti)] - points[prevVerti])
+        );
 
-    // Get y-axis of coordinate system
+    // y-axis of coordinate system
-    vector e1 = n ^ e0;
+    const vector e1 = n ^ e0;
 
     if (mag(mag(e1) - 1) > SMALL)
     {
@@ -489,10 +492,12 @@ Foam::label Foam::intersectedSurface::nextEdge
             )
             {
                 // Calculate angle of edge with respect to base e0, e1
-                vector vec =
+                const vector vec =
-                    n ^ (points[e.otherVertex(prevVerti)] - points[prevVerti]);
+                    normalised
-
+                    (
-                vec /= mag(vec) + VSMALL;
+                        n
+                      ^ (points[e.otherVertex(prevVerti)] - points[prevVerti])
+                    );
 
                 scalar angle = pseudoAngle(e0, e1, vec);
 

src/meshTools/triSurface/booleanOps/surfaceIntersection/edgeIntersections.C

@@ -118,8 +118,7 @@ void Foam::edgeIntersections::intersectEdges
         const point& pStart = points1[meshPoints[e.start()]];
         const point& pEnd = points1[meshPoints[e.end()]];
 
-        const vector eVec(pEnd - pStart);
+        const vector n = normalised(pEnd - pStart);
-        const vector n(eVec/(mag(eVec) + VSMALL));
 
         // Start tracking somewhat before pStart and up to somewhat after p1.
         // Note that tolerances here are smaller than those used to classify
@@ -252,8 +251,7 @@ bool Foam::edgeIntersections::inlinePerturb
             label v0 = surf1.meshPoints()[e[0]];
             label v1 = surf1.meshPoints()[e[1]];
 
-            vector eVec(points1[v1] - points1[v0]);
+            const vector n = normalised(points1[v1] - points1[v0]);
-            vector n = eVec/mag(eVec);
 
             if (perturbStart)
             {
@@ -326,9 +324,7 @@ bool Foam::edgeIntersections::rotatePerturb
             n /= magN;
 
             rndVec -= n*(n & rndVec);
-
+            rndVec.normalise();
-            // Normalize
-            rndVec /= mag(rndVec) + VSMALL;
 
             // Scale to be moved by tolerance.
             rndVec *= 0.01*magN;

src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.C

@@ -410,8 +410,7 @@ void Foam::triSurfaceSearch::findLineAll
 
             if (inter.hit())
             {
-                vector lineVec = end[pointi] - start[pointi];
+                const vector lineVec = normalised(end[pointi] - start[pointi]);
-                lineVec /= mag(lineVec) + VSMALL;
 
                 if
                 (

src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.C

@@ -107,8 +107,7 @@ Foam::pointToPointPlanarInterpolation::calcCoordinateSystem
             << exit(FatalError);
     }
 
-    vector n = e1^(points[index2]-p0);
+    const vector n = normalised(e1 ^ (points[index2]-p0));
-    n /= mag(n);
 
     if (debug)
     {

src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C

@@ -315,11 +315,8 @@ Foam::scalar Foam::triSurfaceTools::faceCosAngle
     const vector base0(pLeft - pStart);
     const vector base1(pRight - pStart);
 
-    vector n0(common ^ base0);
+    const vector n0 = normalised(common ^ base0);
-    n0 /= Foam::mag(n0);
+    const vector n1 = normalised(base1 ^ common);
-
-    vector n1(base1 ^ common);
-    n1 /= Foam::mag(n1);
 
     return n0 & n1;
 }
@@ -2068,11 +2065,11 @@ Foam::vector Foam::triSurfaceTools::surfaceNormal
 
         vector edgeNormal(Zero);
 
-        forAll(eFaces, i)
+        for (const label facei : eFaces)
         {
-            edgeNormal += surf.faceNormals()[eFaces[i]];
+            edgeNormal += surf.faceNormals()[facei];
         }
-        return edgeNormal/(mag(edgeNormal) + VSMALL);
+        return normalised(edgeNormal);
     }
     else
     {
@@ -2626,14 +2623,13 @@ void Foam::triSurfaceTools::calcInterpolationWeights
     edge[1] = tri.a()-tri.c();
     edge[2] = tri.b()-tri.a();
 
-    vector triangleFaceNormal = edge[1] ^ edge[2];
+    const vector triangleFaceNormal = edge[1] ^ edge[2];
 
     // calculate edge normal (pointing inwards)
     FixedList<vector, 3> normal;
     for (label i=0; i<3; i++)
     {
-        normal[i] = triangleFaceNormal ^ edge[i];
+        normal[i] = normalised(triangleFaceNormal ^ edge[i]);
-        normal[i] /= mag(normal[i]) + VSMALL;
     }
 
     weights[0] = ((p-tri.b()) & normal[0]) / max(VSMALL, normal[0] & edge[1]);

src/sampling/sampledSet/circle/circleSet.C

@@ -94,8 +94,7 @@ void Foam::circleSet::calcSamples
     label nPoint = 1;
     while (theta < 360)
     {
-        axis1 = axis1*cosAlpha + (axis1^circleAxis_)*sinAlpha;
+        axis1 = normalised(axis1*cosAlpha + (axis1^circleAxis_)*sinAlpha);
-        axis1 /= mag(axis1);
         point pt = origin_ + radius*axis1;
 
         label celli = searchEngine().findCell(pt);

src/sampling/sampledSet/sampledSet/sampledSet.C

@@ -172,9 +172,7 @@ Foam::scalar Foam::sampledSet::calcSign
 
     vec /= magVec;
 
-    vector n = mesh().faceAreas()[facei];
+    const vector n = normalised(mesh().faceAreas()[facei]);
-
-    n /= mag(n) + VSMALL;
 
     return n & vec;
 }

src/thermophysicalModels/radiation/radiationModels/solarLoad/solarLoad.C

@@ -294,10 +294,9 @@ void Foam::radiation::solarLoad::updateSkyDiffusiveRadiation
 
 void Foam::radiation::solarLoad::initialise(const dictionary& coeffs)
 {
-    if (coeffs.found("gridUp"))
+    if (coeffs.readIfPresent("gridUp", verticalDir_))
     {
-         coeffs.lookup("gridUp") >> verticalDir_;
+         verticalDir_.normalise();
-         verticalDir_ /= mag(verticalDir_);
     }
     else if (mesh_.foundObject<uniformDimensionedVectorField>("g"))
     {

src/thermophysicalModels/radiation/submodels/solarCalculator/solarCalculator.C

@@ -123,12 +123,8 @@ void Foam::solarCalculator::calculateBetaTetha()
 
 void Foam::solarCalculator::calculateSunDirection()
 {
-
+    gridUp_  = normalised(dict_.get<vector>("gridUp"));
-    dict_.lookup("gridUp") >> gridUp_;
+    eastDir_ = normalised(dict_.get<vector>("gridEast"));
-    gridUp_ /= mag(gridUp_);
-
-    dict_.lookup("gridEast") >> eastDir_;
-    eastDir_ /= mag(eastDir_);
 
     coord_.reset
     (
@@ -140,7 +136,7 @@ void Foam::solarCalculator::calculateSunDirection()
     direction_.y() =  cos(beta_)*cos(tetha_); // South axis
     direction_.x() =  cos(beta_)*sin(tetha_); // West axis
 
-    direction_ /= mag(direction_);
+    direction_.normalise();
 
     if (debug)
     {
@@ -163,10 +159,9 @@ void Foam::solarCalculator::init()
     {
         case mSunDirConstant:
         {
-            if (dict_.found("sunDirection"))
+            if (dict_.readIfPresent("sunDirection", direction_))
             {
-                dict_.lookup("sunDirection") >> direction_;
+                direction_.normalise();
-                direction_ /= mag(direction_);
             }
             else
             {

src/waveModels/waveModel/waveModel.C

@@ -57,8 +57,7 @@ void Foam::waveModel::initialiseGeometry()
     // - X: streamwise: patch normal
     // - Y: spanwise: Z^X
     // - Z: up: (negative) gravity direction
-    vector x(-gAverage(patch_.faceAreas()));
+    vector x = normalised(-gAverage(patch_.faceAreas()));
-    x /= mag(x) + ROOTVSMALL;
     vector z = -g_/mag(g_);
     vector y = z^x;