diff --git a/applications/utilities/surface/surfaceInflate/surfaceInflate.C b/applications/utilities/surface/surfaceInflate/surfaceInflate.C
index 2990c26a18..e30ff5fc55 100644
--- a/applications/utilities/surface/surfaceInflate/surfaceInflate.C
+++ b/applications/utilities/surface/surfaceInflate/surfaceInflate.C
@@ -88,7 +88,7 @@ tmp<vectorField> calcVertexNormals(const triSurface& surf)
     // Weight = fA / (mag(e0)^2 * mag(e1)^2);
     tmp<vectorField> tpointNormals
     (
-        new pointField(surf.nPoints(), vector::zero)
+        new pointField(surf.nPoints(), Zero)
     );
     vectorField& pointNormals = tpointNormals.ref();
 
@@ -151,7 +151,7 @@ tmp<vectorField> calcPointNormals
                 {
                     if (!isFeaturePoint[e[i]])
                     {
-                        pointNormals[e[i]] = vector::zero;
+                        pointNormals[e[i]] = Zero;
                     }
                 }
             }
@@ -164,7 +164,7 @@ tmp<vectorField> calcPointNormals
                 const labelList& eFaces = edgeFaces[edgeI];
 
                 // Get average edge normal
-                vector n = vector::zero;
+                vector n = Zero;
                 forAll(eFaces, i)
                 {
                     n += s.faceNormals()[eFaces[i]];
@@ -483,7 +483,7 @@ void lloydsSmoothing
             {
                 const labelList& pFaces = pointFaces[pointI];
 
-                point avg = point::zero;
+                point avg(Zero);
                 forAll(pFaces, pFaceI)
                 {
                     avg += faceCentres[pFaces[pFaceI]];
@@ -498,7 +498,7 @@ void lloydsSmoothing
 
         const pointField& points = s.points();
 
-        vectorField pointSum(s.nPoints(), vector::zero);
+        vectorField pointSum(s.nPoints(), Zero);
         labelList nPointSum(s.nPoints(), 0);
 
         forAll(edges, edgeI)
diff --git a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H
index de477d20e0..299ea32661 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H
@@ -230,7 +230,7 @@ public:
             //  point and density specification
             inline tensor inertia
             (
-                PointRef refPt = vector::zero,
+                PointRef refPt = Zero,
                 scalar density = 1.0
             ) const;
 
diff --git a/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C b/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C
index 860fa12489..6cd22195e5 100644
--- a/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C
+++ b/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C
@@ -525,7 +525,7 @@ void Foam::motionSmootherAlgo::setDisplacement
             label pointI = cppMeshPoints[i];
             if (isPatchPoint[pointI])
             {
-                displacement[pointI] = vector::zero;
+                displacement[pointI] = Zero;
             }
         }
     }
diff --git a/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C b/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C
index 54169aaa1e..c90a197ce7 100644
--- a/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C
+++ b/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C
@@ -118,7 +118,7 @@ void Foam::polyMeshGeometry::updateCellCentresAndVols
         const labelList& cFaces(cells[cellI]);
 
         // Estimate the cell centre and bounding box using the face centres
-        vector cEst = vector::zero;
+        vector cEst(Zero);
         boundBox bb(boundBox::invertedBox);
 
         forAll(cFaces, cFaceI)
diff --git a/src/finiteVolume/interpolation/interpolation/interpolationPointMVC/pointMVCWeight.C b/src/finiteVolume/interpolation/interpolation/interpolationPointMVC/pointMVCWeight.C
index d77af00d4e..6660414716 100644
--- a/src/finiteVolume/interpolation/interpolation/interpolationPointMVC/pointMVCWeight.C
+++ b/src/finiteVolume/interpolation/interpolation/interpolationPointMVC/pointMVCWeight.C
@@ -112,7 +112,7 @@ void Foam::pointMVCWeight::calcWeights
             u(j) = uVec[toLocal[f[j]]];
         }
 
-        vector v(point::zero);
+        vector v(Zero);
         forAll(f, j)
         {
             label jPlus1 = f.fcIndex(j);
diff --git a/src/fvMotionSolver/fvMotionSolvers/displacement/SBRStress/displacementSBRStressFvMotionSolver.C b/src/fvMotionSolver/fvMotionSolvers/displacement/SBRStress/displacementSBRStressFvMotionSolver.C
index 5ffef6acab..2b26d8d58a 100644
--- a/src/fvMotionSolver/fvMotionSolvers/displacement/SBRStress/displacementSBRStressFvMotionSolver.C
+++ b/src/fvMotionSolver/fvMotionSolvers/displacement/SBRStress/displacementSBRStressFvMotionSolver.C
@@ -124,7 +124,7 @@ displacementSBRStressFvMotionSolver
         (
             "cellDisplacement",
             displacementMotionSolver::pointDisplacement().dimensions(),
-            vector::zero
+            Zero
         ),
         cellMotionBoundaryTypes<vector>
         (
diff --git a/src/fvMotionSolver/fvMotionSolvers/displacement/laplacian/displacementLaplacianFvMotionSolver.C b/src/fvMotionSolver/fvMotionSolvers/displacement/laplacian/displacementLaplacianFvMotionSolver.C
index 3021e09675..3d5f22d09d 100644
--- a/src/fvMotionSolver/fvMotionSolvers/displacement/laplacian/displacementLaplacianFvMotionSolver.C
+++ b/src/fvMotionSolver/fvMotionSolvers/displacement/laplacian/displacementLaplacianFvMotionSolver.C
@@ -169,7 +169,7 @@ displacementLaplacianFvMotionSolver
         (
             "cellDisplacement",
             pointDisplacement_.dimensions(),
-            vector::zero
+            Zero
         ),
         cellMotionBoundaryTypes<vector>(pointDisplacement_.boundaryField())
     ),
diff --git a/src/fvMotionSolver/fvMotionSolvers/displacement/surfaceAlignedSBRStress/surfaceAlignedSBRStressFvMotionSolver.C b/src/fvMotionSolver/fvMotionSolvers/displacement/surfaceAlignedSBRStress/surfaceAlignedSBRStressFvMotionSolver.C
index 165ba536b8..7a342f1998 100644
--- a/src/fvMotionSolver/fvMotionSolvers/displacement/surfaceAlignedSBRStress/surfaceAlignedSBRStressFvMotionSolver.C
+++ b/src/fvMotionSolver/fvMotionSolvers/displacement/surfaceAlignedSBRStress/surfaceAlignedSBRStressFvMotionSolver.C
@@ -74,7 +74,7 @@ surfaceAlignedSBRStressFvMotionSolver
             IOobject::NO_WRITE
         ),
         fvMesh_,
-        dimensionedVector("zero", dimless, vector::zero)
+        dimensionedVector("zero", dimless, Zero)
     ),
     maxAng_(coeffDict().lookupOrDefault<scalar>("maxAng", 80.0)),
     minAng_(coeffDict().lookupOrDefault<scalar>("minAng", 20.0)),
@@ -93,7 +93,7 @@ surfaceAlignedSBRStressFvMotionSolver
             IOobject::NO_WRITE
         ),
         fvMesh_,
-        dimensionedSymmTensor("zero", dimless, symmTensor::zero)
+        dimensionedSymmTensor("zero", dimless, Zero)
     ),
     minSigmaDiff_(coeffDict().lookupOrDefault<scalar>("minSigmaDiff", 1e-4))
 {
@@ -131,8 +131,8 @@ Foam::surfaceAlignedSBRStressFvMotionSolver::
 
 void Foam::surfaceAlignedSBRStressFvMotionSolver::calculateCellRot()
 {
-    cellRot_.internalField() = vector::zero;
-    pointDisplacement_.internalField() = vector::zero;
+    cellRot_.internalField() = Zero;
+    pointDisplacement_.internalField() = Zero;
 
     // Find intersections
     pointField start(fvMesh_.nInternalFaces());
@@ -316,7 +316,7 @@ void Foam::surfaceAlignedSBRStressFvMotionSolver::solve()
                 IOobject::NO_WRITE
             ),
             fvMesh_,
-            dimensionedVector("zero", dimLength, vector::zero),
+            dimensionedVector("zero", dimLength, Zero),
             cellMotionBoundaryTypes<vector>
             (
                 pointDisplacement().boundaryField()
diff --git a/src/fvMotionSolver/motionInterpolation/patchTransformed/patchTransformedInterpolation.C b/src/fvMotionSolver/motionInterpolation/patchTransformed/patchTransformedInterpolation.C
index 22c9d28c32..331c73c858 100644
--- a/src/fvMotionSolver/motionInterpolation/patchTransformed/patchTransformedInterpolation.C
+++ b/src/fvMotionSolver/motionInterpolation/patchTransformed/patchTransformedInterpolation.C
@@ -123,11 +123,11 @@ void Foam::patchTransformedInterpolation::interpolate
 
     pointDisplacement.correctBoundaryConditions();
 
-    vectorField pointRotation(nPoints, vector::zero);
+    vectorField pointRotation(nPoints, Zero);
     scalarField pointExpansion(nPoints, scalar(0));
 
     labelList pointDisplacementNSum(nPoints, 0);
-    vectorField pointDisplacementSum(nPoints, vector::zero);
+    vectorField pointDisplacementSum(nPoints, Zero);
 
     forAll(patches_, patchI)
     {
diff --git a/src/fvOptions/sources/derived/directionalPressureGradientExplicitSource/directionalPressureGradientExplicitSource.C b/src/fvOptions/sources/derived/directionalPressureGradientExplicitSource/directionalPressureGradientExplicitSource.C
index fec009d7f8..2c170d976b 100644
--- a/src/fvOptions/sources/derived/directionalPressureGradientExplicitSource/directionalPressureGradientExplicitSource.C
+++ b/src/fvOptions/sources/derived/directionalPressureGradientExplicitSource/directionalPressureGradientExplicitSource.C
@@ -179,9 +179,9 @@ directionalPressureGradientExplicitSource
 :
     cellSetOption(sourceName, modelType, dict, mesh),
     model_(PressureDropModelNames_.read(coeffs_.lookup("model"))),
-    gradP0_(cells_.size(), vector::zero),
-    dGradP_(cells_.size(), vector::zero),
-    gradPporous_(cells_.size(), vector::zero),
+    gradP0_(cells_.size(), Zero),
+    dGradP_(cells_.size(), Zero),
+    gradPporous_(cells_.size(), Zero),
     flowDir_(coeffs_.lookup("flowDir")),
     invAPtr_(NULL),
     D_(0),
@@ -386,7 +386,7 @@ void Foam::fv::directionalPressureGradientExplicitSource::correct
     }
 
     // Accumulate 'upstream' velocity into cells
-    vectorField UfCells(cells_.size(), vector::zero);
+    vectorField UfCells(cells_.size(), Zero);
     scalarField UfCellWeights(cells_.size(), 0.0);
 
     const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
@@ -486,7 +486,7 @@ void Foam::fv::directionalPressureGradientExplicitSource::addSup
             IOobject::NO_WRITE
         ),
         mesh_,
-        dimensionedVector("zero", eqn.dimensions()/dimVolume, vector::zero)
+        dimensionedVector("zero", eqn.dimensions()/dimVolume, Zero)
     );
 
     UIndirectList<vector>(Su, cells_) = gradP0_ + dGradP_ + gradPporous_;
@@ -536,7 +536,7 @@ void Foam::fv::directionalPressureGradientExplicitSource::constrain
     }
 
     gradP0_ += dGradP_;
-    dGradP_ = vector::zero;
+    dGradP_ = Zero;
 }
 
 
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/KinematicLookupTableInjection/kinematicParcelInjectionData.C b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/KinematicLookupTableInjection/kinematicParcelInjectionData.C
index 12364b9dfe..a831d77283 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/KinematicLookupTableInjection/kinematicParcelInjectionData.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/KinematicLookupTableInjection/kinematicParcelInjectionData.C
@@ -36,7 +36,7 @@ namespace Foam
 
 Foam::kinematicParcelInjectionData::kinematicParcelInjectionData()
 :
-    x_(point::zero),
+    x_(Zero),
     U_(Zero),
     d_(0.0),
     rho_(0.0),
diff --git a/src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C b/src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C
index 46ae5f4686..f68aa1d09a 100644
--- a/src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C
+++ b/src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C
@@ -630,13 +630,13 @@ bool Foam::medialAxisMeshMover::unmarkExtrusion
     if (extrudeStatus[patchPointI] == snappyLayerDriver::EXTRUDE)
     {
         extrudeStatus[patchPointI] = snappyLayerDriver::NOEXTRUDE;
-        patchDisp[patchPointI] = vector::zero;
+        patchDisp[patchPointI] = Zero;
         return true;
     }
     else if (extrudeStatus[patchPointI] == snappyLayerDriver::EXTRUDEREMOVE)
     {
         extrudeStatus[patchPointI] = snappyLayerDriver::NOEXTRUDE;
-        patchDisp[patchPointI] = vector::zero;
+        patchDisp[patchPointI] = Zero;
         return true;
     }
     else
@@ -1224,7 +1224,7 @@ Foam::medialAxisMeshMover::medialAxisMeshMover
             false
         ),
         pMesh(),
-        dimensionedVector("dispVec", dimLength, vector::zero)
+        dimensionedVector("dispVec", dimLength, Zero)
     ),
     medialRatio_
     (
@@ -1266,7 +1266,7 @@ Foam::medialAxisMeshMover::medialAxisMeshMover
             false
         ),
         pMesh(),
-        dimensionedVector("medialVec", dimLength, vector::zero)
+        dimensionedVector("medialVec", dimLength, Zero)
     )
 {
     update(dict);
@@ -1594,7 +1594,7 @@ void Foam::medialAxisMeshMover::calculateDisplacement
     {
         if (!pointWallDist[pointI].valid(dummyTrackData))
         {
-            displacement[pointI] = vector::zero;
+            displacement[pointI] = Zero;
         }
         else
         {
diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C
index 749989e44f..db59fd4fc3 100644
--- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C
+++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C
@@ -1737,7 +1737,7 @@ Foam::tmp<Foam::pointVectorField> Foam::meshRefinement::makeDisplacementField
                 IOobject::AUTO_WRITE
             ),
             pMesh,
-            dimensionedVector("displacement", dimLength, vector::zero),
+            dimensionedVector("displacement", dimLength, Zero),
             patchFieldTypes
         )
     );
diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementRefine.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementRefine.C
index 4ffa1248f8..c9dba90ec0 100644
--- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementRefine.C
+++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementRefine.C
@@ -1558,7 +1558,7 @@ bool Foam::meshRefinement::isGap
     {
         scalar cosAngle = (normal0 & normal1);
 
-        vector avg = vector::zero;
+        vector avg = Zero;
         if (cosAngle < (-1+planarCos))
         {
             // Opposite normals
@@ -1615,7 +1615,7 @@ bool Foam::meshRefinement::isNormalGap
     {
         scalar cosAngle = (normal0 & normal1);
 
-        vector avg = vector::zero;
+        vector avg = Zero;
         if (cosAngle < (-1+planarCos))
         {
             // Opposite normals
@@ -1782,8 +1782,8 @@ Foam::label Foam::meshRefinement::markProximityRefinement
     // minLevel) and cache per cell the max surface level and the local normal
     // on that surface.
     labelList cellMaxLevel(mesh_.nCells(), -1);
-    vectorField cellMaxNormal(mesh_.nCells(), vector::zero);
-    pointField cellMaxLocation(mesh_.nCells(), vector::zero);
+    vectorField cellMaxNormal(mesh_.nCells(), Zero);
+    pointField cellMaxLocation(mesh_.nCells(), Zero);
 
     {
         // Per segment the normals of the surfaces
diff --git a/src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.C b/src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.C
index c9532641bd..514b1e3b06 100644
--- a/src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.C
+++ b/src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.C
@@ -536,7 +536,7 @@ void Foam::refinementFeatures::findNearestEdge
     nearInfo.setSize(samples.size());
     nearInfo = pointIndexHit();
     nearNormal.setSize(samples.size());
-    nearNormal = vector::zero;
+    nearNormal = Zero;
 
     forAll(edgeTrees_, featI)
     {
@@ -595,7 +595,7 @@ void Foam::refinementFeatures::findNearestRegionEdge
     nearInfo.setSize(samples.size());
     nearInfo = pointIndexHit();
     nearNormal.setSize(samples.size());
-    nearNormal = vector::zero;
+    nearNormal = Zero;
 
 
     const PtrList<indexedOctree<treeDataEdge>>& regionTrees =
diff --git a/src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.C b/src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.C
index f37714e561..b3f56f621a 100644
--- a/src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.C
+++ b/src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.C
@@ -1216,7 +1216,7 @@ void Foam::refinementSurfaces::findNearestIntersection
     region1.setSize(start.size());
     region1 = -1;
     normal1.setSize(start.size());
-    normal1 = vector::zero;
+    normal1 = Zero;
 
     // Current end of segment to test.
     pointField nearest(end);
@@ -1325,7 +1325,7 @@ void Foam::refinementSurfaces::findNearestIntersection
     surface1.setSize(start.size());
     surface1 = -1;
     normal1.setSize(start.size());
-    normal1 = vector::zero;
+    normal1 = Zero;
 
     // Current end of segment to test.
     pointField nearest(end);
@@ -1374,7 +1374,7 @@ void Foam::refinementSurfaces::findNearestIntersection
     hitInfo1.setSize(start.size());
     hitInfo1 = pointIndexHit();
     normal1.setSize(start.size());
-    normal1 = vector::zero;
+    normal1 = Zero;
 
     // Current end of segment to test.
     pointField nearest(end);
@@ -1559,7 +1559,7 @@ void Foam::refinementSurfaces::findNearestRegion
     hitRegion.setSize(hitSurface.size());
     hitRegion = -1;
     hitNormal.setSize(hitSurface.size());
-    hitNormal = vector::zero;
+    hitNormal = Zero;
 
     forAll(surfacesToTest, i)
     {
@@ -1784,7 +1784,7 @@ void Foam::refinementSurfaces::findNearestRegion
     hitRegion.setSize(hitSurface.size());
     hitRegion = -1;
     hitNormal.setSize(hitSurface.size());
-    hitNormal = vector::zero;
+    hitNormal = Zero;
 
     forAll(surfacesToTest, i)
     {
diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.C b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.C
index 2d5d3621a2..cfbb7c3de8 100644
--- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.C
+++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.C
@@ -223,14 +223,14 @@ bool Foam::snappyLayerDriver::unmarkExtrusion
     {
         extrudeStatus[patchPointI] = NOEXTRUDE;
         patchNLayers[patchPointI] = 0;
-        patchDisp[patchPointI] = vector::zero;
+        patchDisp[patchPointI] = Zero;
         return true;
     }
     else if (extrudeStatus[patchPointI] == EXTRUDEREMOVE)
     {
         extrudeStatus[patchPointI] = NOEXTRUDE;
         patchNLayers[patchPointI] = 0;
-        patchDisp[patchPointI] = vector::zero;
+        patchDisp[patchPointI] = Zero;
         return true;
     }
     else
@@ -866,7 +866,7 @@ Foam::snappyLayerDriver::makeLayerDisplacementField
                 IOobject::AUTO_WRITE
             ),
             pMesh,
-            dimensionedVector("displacement", dimLength, vector::zero),
+            dimensionedVector("displacement", dimLength, Zero),
             patchFieldTypes,
             actualPatchTypes
         )
@@ -953,7 +953,7 @@ void Foam::snappyLayerDriver::growNoExtrusion
     {
         if (extrudeStatus[patchPointI] == NOEXTRUDE)
         {
-            patchDisp[patchPointI] = vector::zero;
+            patchDisp[patchPointI] = Zero;
             patchNLayers[patchPointI] = 0;
         }
     }
@@ -1532,7 +1532,7 @@ void Foam::snappyLayerDriver::getPatchDisplacement
         {
             // Do not use unmarkExtrusion; forcibly set to zero extrusion.
             patchNLayers[patchPointI] = 0;
-            patchDisp[patchPointI] = vector::zero;
+            patchDisp[patchPointI] = Zero;
         }
         else
         {
@@ -1561,7 +1561,7 @@ void Foam::snappyLayerDriver::getPatchDisplacement
     {
         if (extrudeStatus[patchPointI] == EXTRUDEREMOVE)
         {
-            point avg(vector::zero);
+            point avg(Zero);
             label nPoints = 0;
 
             const labelList& pEdges = pp.pointEdges()[patchPointI];
@@ -1798,7 +1798,7 @@ Foam::label Foam::snappyLayerDriver::truncateDisplacement
         else if (extrudeStatus[patchPointI] == NOEXTRUDE)
         {
             // Make sure displacement is 0. Should already be so but ...
-            patchDisp[patchPointI] = vector::zero;
+            patchDisp[patchPointI] = Zero;
             patchNLayers[patchPointI] = 0;
         }
     }
@@ -3681,7 +3681,7 @@ void Foam::snappyLayerDriver::addLayers
 
             // Calculate displacement for final layer for addPatchLayer.
             // (layer of cells next to the original mesh)
-            vectorField finalDisp(patchNLayers.size(), vector::zero);
+            vectorField finalDisp(patchNLayers.size(), Zero);
 
             forAll(nPatchPointLayers, i)
             {
diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C
index 99957c9fa4..484229bf93 100644
--- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C
+++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C
@@ -1025,7 +1025,7 @@ Foam::label Foam::snappyRefineDriver::refinementInterfaceRefine
                 //    const scalar rCVol = pow(cellVolumes[cellI], -5.0/3.0);
                 //
                 //    // Determine principal axes of cell
-                //    symmTensor R(symmTensor::zero);
+                //    symmTensor R(Zero);
                 //
                 //    forAll(cFaces, i)
                 //    {
diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.C b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.C
index 0af06a6e75..95178d2881 100644
--- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.C
+++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.C
@@ -250,7 +250,7 @@ Foam::tmp<Foam::pointField> Foam::snappySnapDriver::smoothInternalDisplacement
 
     // Calculate average of connected cells
     labelList nCells(mesh.nPoints(), 0);
-    pointField sumLocation(mesh.nPoints(), vector::zero);
+    pointField sumLocation(mesh.nPoints(), Zero);
 
     forAll(isMovingPoint, pointI)
     {
@@ -276,7 +276,7 @@ Foam::tmp<Foam::pointField> Foam::snappySnapDriver::smoothInternalDisplacement
         vector::zero
     );
 
-    tmp<pointField> tdisplacement(new pointField(mesh.nPoints(), vector::zero));
+    tmp<pointField> tdisplacement(new pointField(mesh.nPoints(), Zero));
     pointField& displacement = tdisplacement.ref();
 
     label nAdapted = 0;
@@ -370,7 +370,7 @@ Foam::tmp<Foam::pointField> Foam::snappySnapDriver::smoothPatchDisplacement
     // Get average position of boundary face centres
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    vectorField avgBoundary(pointFaces.size(), vector::zero);
+    vectorField avgBoundary(pointFaces.size(), Zero);
     labelList nBoundary(pointFaces.size(), 0);
 
     forAll(pointFaces, patchPointI)
@@ -418,7 +418,7 @@ Foam::tmp<Foam::pointField> Foam::snappySnapDriver::smoothPatchDisplacement
     vectorField avgInternal;
     labelList nInternal;
     {
-        vectorField globalSum(mesh.nPoints(), vector::zero);
+        vectorField globalSum(mesh.nPoints(), Zero);
         labelList globalNum(mesh.nPoints(), 0);
 
         // Note: no use of pointFaces
@@ -521,7 +521,7 @@ Foam::tmp<Foam::pointField> Foam::snappySnapDriver::smoothPatchDisplacement
 
 
     // Displacement to calculate.
-    tmp<pointField> tpatchDisp(new pointField(meshPoints.size(), vector::zero));
+    tmp<pointField> tpatchDisp(new pointField(meshPoints.size(), Zero));
     pointField& patchDisp = tpatchDisp.ref();
 
     forAll(pointFaces, i)
@@ -597,7 +597,7 @@ Foam::tmp<Foam::pointField> Foam::snappySnapDriver::smoothPatchDisplacement
 //    const labelListList& pointEdges = pp.pointEdges();
 //    const edgeList& edges = pp.edges();
 //
-//    tmp<pointField> tavg(new pointField(pointEdges.size(), vector::zero));
+//    tmp<pointField> tavg(new pointField(pointEdges.size(), Zero));
 //    pointField& avg = tavg();
 //
 //    forAll(pointEdges, vertI)
@@ -997,7 +997,7 @@ Foam::tmp<Foam::pointField> Foam::snappySnapDriver::avgCellCentres
 
     tmp<pointField> tavgBoundary
     (
-        new pointField(pointFaces.size(), vector::zero)
+        new pointField(pointFaces.size(), Zero)
     );
     pointField& avgBoundary = tavgBoundary.ref();
     labelList nBoundary(pointFaces.size(), 0);
@@ -1805,7 +1805,7 @@ Foam::vectorField Foam::snappySnapDriver::calcNearestSurface
     const fvMesh& mesh = meshRefiner.mesh();
 
     // Displacement per patch point
-    vectorField patchDisp(localPoints.size(), vector::zero);
+    vectorField patchDisp(localPoints.size(), Zero);
 
     if (returnReduce(localPoints.size(), sumOp<label>()) > 0)
     {
@@ -2754,7 +2754,7 @@ void Foam::snappySnapDriver::doSnap
             if (snapParams.detectNearSurfacesSnap())
             {
                 nearestPoint.setSize(pp.nPoints(), vector::max);
-                nearestNormal.setSize(pp.nPoints(), vector::zero);
+                nearestNormal.setSize(pp.nPoints(), Zero);
             }
 
             vectorField disp = calcNearestSurface
diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriverFeature.C b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriverFeature.C
index 80bc890eb2..45535a7f19 100644
--- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriverFeature.C
+++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriverFeature.C
@@ -152,7 +152,7 @@ void Foam::snappySnapDriver::smoothAndConstrain
         // - same for feature points. They are already attracted to the
         //   nearest feature point.
 
-        vectorField dispSum(pp.nPoints(), vector::zero);
+        vectorField dispSum(pp.nPoints(), Zero);
         labelList dispCount(pp.nPoints(), 0);
 
         const labelListList& pointEdges = pp.pointEdges();
@@ -233,9 +233,9 @@ void Foam::snappySnapDriver::calcNearestFace
 
     // Displacement and orientation per pp face.
     faceDisp.setSize(pp.size());
-    faceDisp = vector::zero;
+    faceDisp = Zero;
     faceSurfaceNormal.setSize(pp.size());
-    faceSurfaceNormal = vector::zero;
+    faceSurfaceNormal = Zero;
     faceSurfaceGlobalRegion.setSize(pp.size());
     faceSurfaceGlobalRegion = -1;
 
@@ -408,7 +408,7 @@ void Foam::snappySnapDriver::calcNearestFace
     //
     //// Determine rotation axis
     //faceRotation.setSize(pp.size());
-    //faceRotation = vector::zero;
+    //faceRotation = Zero;
     //
     //forAll(faceRotation, faceI)
     //{
@@ -745,7 +745,7 @@ Foam::pointIndexHit Foam::snappySnapDriver::findMultiPatchPoint
             }
         }
     }
-    return pointIndexHit(false, vector::zero, labelMax);
+    return pointIndexHit(false, Zero, labelMax);
 }
 
 
@@ -933,7 +933,7 @@ void Foam::snappySnapDriver::featureAttractionUsingReconstruction
     pointConstraint& patchConstraint
 ) const
 {
-    patchAttraction = vector::zero;
+    patchAttraction = Zero;
     patchConstraint = pointConstraint();
 
     const List<point>& pfSurfNormals = pointFaceSurfNormals[pointI];
@@ -1155,7 +1155,7 @@ void Foam::snappySnapDriver::featureAttractionUsingReconstruction
 
     forAll(pp.localPoints(), pointI)
     {
-        vector attraction = vector::zero;
+        vector attraction = Zero;
         pointConstraint constraint;
 
         featureAttractionUsingReconstruction
@@ -1510,7 +1510,7 @@ void Foam::snappySnapDriver::releasePointsNextToMultiPatch
                     //Pout<< "Knocking out constraint"
                     //    << " on non-multiPatchPoint:"
                     //    << pp.localPoints()[pointI] << endl;
-                    patchAttraction[pointI] = vector::zero;
+                    patchAttraction[pointI] = Zero;
                     patchConstraints[pointI] = pointConstraint();
                     nChanged++;
 
@@ -2193,14 +2193,14 @@ Foam::snappySnapDriver::findNearFeaturePoint
                 // Current pointI nearer.
                 pointAttractor[featI][featPointI] = pointI;
                 pointConstraints[featI][featPointI].first() = 3;
-                pointConstraints[featI][featPointI].second() = vector::zero;
+                pointConstraints[featI][featPointI].second() = Zero;
 
                 // Store for later use
                 patchAttraction[pointI] = featPt-pt;
                 patchConstraints[pointI] = pointConstraints[featI][featPointI];
 
                 // Reset oldPointI to nearest on feature edge
-                patchAttraction[oldPointI] = vector::zero;
+                patchAttraction[oldPointI] = Zero;
                 patchConstraints[oldPointI] = pointConstraint();
 
                 findNearFeatureEdge
@@ -2224,7 +2224,7 @@ Foam::snappySnapDriver::findNearFeaturePoint
             // Current pointI nearer.
             pointAttractor[featI][featPointI] = pointI;
             pointConstraints[featI][featPointI].first() = 3;
-            pointConstraints[featI][featPointI].second() = vector::zero;
+            pointConstraints[featI][featPointI].second() = Zero;
 
             // Store for later use
             patchAttraction[pointI] = featPt-pt;
@@ -2346,7 +2346,7 @@ void Foam::snappySnapDriver::determineFeatures
         //    2: attract to feature line, constraint is feature line direction
         //    3: attract to feature point, constraint is zero)
 
-        vector attraction = vector::zero;
+        vector attraction = Zero;
         pointConstraint constraint;
 
         featureAttractionUsingReconstruction
@@ -3029,8 +3029,7 @@ void Foam::snappySnapDriver::determineBaffleFeatures
                 {
                     pointAttractor[featI][featPointI] = pointI;
                     pointConstraints[featI][featPointI].first() = 3;
-                    pointConstraints[featI][featPointI].second() =
-                        vector::zero;
+                    pointConstraints[featI][featPointI].second() = Zero;
 
                     // Store for later use
                     patchAttraction[pointI] = featPt-pt;
@@ -3245,7 +3244,7 @@ void Foam::snappySnapDriver::reverseAttractMeshPoints
 
     // Per mesh point the point on nearest feature edge.
     patchAttraction.setSize(pp.nPoints());
-    patchAttraction = vector::zero;
+    patchAttraction = Zero;
     patchConstraints.setSize(pp.nPoints());
     patchConstraints = pointConstraint();
 
@@ -3440,7 +3439,7 @@ void Foam::snappySnapDriver::featureAttractionUsingFeatureEdges
     }
 
     // Reverse: from pp point to nearest feature
-    vectorField rawPatchAttraction(pp.nPoints(), vector::zero);
+    vectorField rawPatchAttraction(pp.nPoints(), Zero);
     List<pointConstraint> rawPatchConstraints(pp.nPoints());
 
     determineFeatures
@@ -3875,11 +3874,11 @@ Foam::vectorField Foam::snappySnapDriver::calcNearestSurfaceFeature
         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
         // vector from point on surface back to face centre
-        vectorField faceDisp(pp.size(), vector::zero);
+        vectorField faceDisp(pp.size(), Zero);
         // normal of surface at point on surface
-        vectorField faceSurfaceNormal(pp.size(), vector::zero);
+        vectorField faceSurfaceNormal(pp.size(), Zero);
         labelList faceSurfaceGlobalRegion(pp.size(), -1);
-        //vectorField faceRotation(pp.size(), vector::zero);
+        //vectorField faceRotation(pp.size(), Zero);
 
         calcNearestFace
         (
@@ -3927,7 +3926,7 @@ Foam::vectorField Foam::snappySnapDriver::calcNearestSurfaceFeature
 
     // Nearest feature
     patchAttraction.setSize(localPoints.size());
-    patchAttraction = vector::zero;
+    patchAttraction = Zero;
     // Constraints at feature
     patchConstraints.setSize(localPoints.size());
     patchConstraints = pointConstraint();
diff --git a/src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C b/src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C
index c9db04a276..71773ffc53 100644
--- a/src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C
+++ b/src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C
@@ -490,7 +490,7 @@ Foam::cyclicAMIPolyPatch::cyclicAMIPolyPatch
     nbrPatchName_(word::null),
     nbrPatchID_(-1),
     rotationAxis_(Zero),
-    rotationCentre_(point::zero),
+    rotationCentre_(Zero),
     rotationAngleDefined_(false),
     rotationAngle_(0.0),
     separationVector_(Zero),
@@ -521,7 +521,7 @@ Foam::cyclicAMIPolyPatch::cyclicAMIPolyPatch
     coupleGroup_(dict),
     nbrPatchID_(-1),
     rotationAxis_(Zero),
-    rotationCentre_(point::zero),
+    rotationCentre_(Zero),
     rotationAngleDefined_(false),
     rotationAngle_(0.0),
     separationVector_(Zero),
diff --git a/src/meshTools/coordinateSystems/coordinateRotation/cylindrical.C b/src/meshTools/coordinateSystems/coordinateRotation/cylindrical.C
index b5e03a3896..4bccbfed5f 100644
--- a/src/meshTools/coordinateSystems/coordinateRotation/cylindrical.C
+++ b/src/meshTools/coordinateSystems/coordinateRotation/cylindrical.C
@@ -99,7 +99,7 @@ Foam::cylindrical::cylindrical
 )
 :
     Rptr_(),
-    origin_(point::zero),
+    origin_(Zero),
     e3_(Zero)
 {
     // If origin is specified in the coordinateSystem
diff --git a/src/meshTools/coordinateSystems/coordinateSystem.C b/src/meshTools/coordinateSystems/coordinateSystem.C
index 3abb25524b..e3150ff8f5 100644
--- a/src/meshTools/coordinateSystems/coordinateSystem.C
+++ b/src/meshTools/coordinateSystems/coordinateSystem.C
@@ -43,7 +43,7 @@ Foam::coordinateSystem::coordinateSystem()
 :
     name_(),
     note_(),
-    origin_(point::zero),
+    origin_(Zero),
     R_(new axesRotation(sphericalTensor::I))
 {}
 
@@ -98,7 +98,7 @@ Foam::coordinateSystem::coordinateSystem
 :
     name_(name),
     note_(),
-    origin_(point::zero),
+    origin_(Zero),
     R_()
 {
     init(dict);
@@ -109,7 +109,7 @@ Foam::coordinateSystem::coordinateSystem(const dictionary& dict)
 :
     name_(),
     note_(),
-    origin_(point::zero),
+    origin_(Zero),
     R_()
 {
     init(dict);
@@ -124,7 +124,7 @@ Foam::coordinateSystem::coordinateSystem
 :
     name_(),
     note_(),
-    origin_(point::zero),
+    origin_(Zero),
     R_()
 {
     const entry* entryPtr = dict.lookupEntryPtr(typeName_(), false, false);
@@ -168,7 +168,7 @@ Foam::coordinateSystem::coordinateSystem(Istream& is)
 :
     name_(is),
     note_(),
-    origin_(point::zero),
+    origin_(Zero),
     R_()
 {
     dictionary dict(is);
diff --git a/src/meshTools/searchableSurface/searchableCone.C b/src/meshTools/searchableSurface/searchableCone.C
index 9cbf6651d3..02d0b1b16c 100644
--- a/src/meshTools/searchableSurface/searchableCone.C
+++ b/src/meshTools/searchableSurface/searchableCone.C
@@ -100,7 +100,7 @@ void Foam::searchableCone::findNearestAndNormal
     scalar magV = mag(v);
     if (magV < ROOTVSMALL)
     {
-       v = vector::zero;
+       v = Zero;
     }
     else
     {
@@ -1045,7 +1045,7 @@ void Foam::searchableCone::getNormal
 ) const
 {
     normal.setSize(info.size());
-    normal = vector::zero;
+    normal = Zero;
 
     forAll(info, i)
     {
diff --git a/src/meshTools/searchableSurface/searchableRotatedBox.C b/src/meshTools/searchableSurface/searchableRotatedBox.C
index dd4dd06227..ffb8dc9c5d 100644
--- a/src/meshTools/searchableSurface/searchableRotatedBox.C
+++ b/src/meshTools/searchableSurface/searchableRotatedBox.C
@@ -57,7 +57,7 @@ Foam::searchableRotatedBox::searchableRotatedBox
             io.writeOpt(),
             false      //io.registerObject(),
         ),
-        treeBoundBox(point::zero, dict.lookup("span"))
+        treeBoundBox(Zero, dict.lookup("span"))
     ),
     transform_
     (
diff --git a/src/meshTools/tetOverlapVolume/tetOverlapVolume.H b/src/meshTools/tetOverlapVolume/tetOverlapVolume.H
index 879268a009..b027599218 100644
--- a/src/meshTools/tetOverlapVolume/tetOverlapVolume.H
+++ b/src/meshTools/tetOverlapVolume/tetOverlapVolume.H
@@ -70,7 +70,7 @@ class tetOverlapVolume
 
             inline sumMomentOp()
             :
-                vol_(0.0, vector::zero)
+                vol_(0.0, Zero)
             {}
 
             inline void operator()(const tetPoints& tet)
diff --git a/src/postProcessing/functionObjects/forces/forceCoeffs/forceCoeffs.C b/src/postProcessing/functionObjects/forces/forceCoeffs/forceCoeffs.C
index 93b9dff88e..5ec4846025 100644
--- a/src/postProcessing/functionObjects/forces/forceCoeffs/forceCoeffs.C
+++ b/src/postProcessing/functionObjects/forces/forceCoeffs/forceCoeffs.C
@@ -273,7 +273,7 @@ void Foam::forceCoeffs::read(const dictionary& dict)
                     IOobject::NO_WRITE
                 ),
                 mesh,
-                dimensionedVector("0", dimless, vector::zero)
+                dimensionedVector("0", dimless, Zero)
             )
         );
 
@@ -292,7 +292,7 @@ void Foam::forceCoeffs::read(const dictionary& dict)
                     IOobject::NO_WRITE
                 ),
                 mesh,
-                dimensionedVector("0", dimless, vector::zero)
+                dimensionedVector("0", dimless, Zero)
             )
         );
 
diff --git a/src/postProcessing/functionObjects/forces/forces/forces.C b/src/postProcessing/functionObjects/forces/forces/forces.C
index 5a3e1ca796..433cdbfcb9 100644
--- a/src/postProcessing/functionObjects/forces/forces/forces.C
+++ b/src/postProcessing/functionObjects/forces/forces/forces.C
@@ -304,7 +304,7 @@ void Foam::forces::resetFields()
                 obr_.lookupObject<volVectorField>(fieldName("force"))
             );
 
-        force == dimensionedVector("0", force.dimensions(), vector::zero);
+        force == dimensionedVector("0", force.dimensions(), Zero);
 
         volVectorField& moment =
             const_cast<volVectorField&>
@@ -312,7 +312,7 @@ void Foam::forces::resetFields()
                 obr_.lookupObject<volVectorField>(fieldName("moment"))
             );
 
-        moment == dimensionedVector("0", moment.dimensions(), vector::zero);
+        moment == dimensionedVector("0", moment.dimensions(), Zero);
     }
 }
 
@@ -968,7 +968,7 @@ void Foam::forces::read(const dictionary& dict)
                     IOobject::NO_WRITE
                 ),
                 mesh,
-                dimensionedVector("0", dimForce, vector::zero)
+                dimensionedVector("0", dimForce, Zero)
             )
         );
 
@@ -987,7 +987,7 @@ void Foam::forces::read(const dictionary& dict)
                     IOobject::NO_WRITE
                 ),
                 mesh,
-                dimensionedVector("0", dimForce*dimLength, vector::zero)
+                dimensionedVector("0", dimForce*dimLength, Zero)
             )
         );
 
diff --git a/src/postProcessing/functionObjects/utilities/fluxSummary/fluxSummary.C b/src/postProcessing/functionObjects/utilities/fluxSummary/fluxSummary.C
index 60097f08af..0ac60ab308 100644
--- a/src/postProcessing/functionObjects/utilities/fluxSummary/fluxSummary.C
+++ b/src/postProcessing/functionObjects/utilities/fluxSummary/fluxSummary.C
@@ -186,7 +186,7 @@ void Foam::fluxSummary::initialiseFaceZoneAndDirection
     const surfaceVectorField& Sf = mesh.Sf();
     const surfaceScalarField& magSf = mesh.magSf();
 
-    vector n = vector::zero;
+    vector n(Zero);
 
     forAll(fZone, i)
     {
diff --git a/src/sampling/meshToMesh/calcMethod/meshToMeshMethod/meshToMeshMethod.C b/src/sampling/meshToMesh/calcMethod/meshToMeshMethod/meshToMeshMethod.C
index 53b13c6da9..04afc2cef8 100644
--- a/src/sampling/meshToMesh/calcMethod/meshToMeshMethod/meshToMeshMethod.C
+++ b/src/sampling/meshToMesh/calcMethod/meshToMeshMethod/meshToMeshMethod.C
@@ -146,7 +146,7 @@ Foam::meshToMeshMethod::interVolAndCentroid
     if (volAndInertia.first() <= ROOTVSMALL)
     {
         volAndInertia.first() = 0.0;
-        volAndInertia.second() = vector::zero;
+        volAndInertia.second() = Zero;
     }
     else
     {
diff --git a/src/sampling/sampledSurface/writers/foamFile/foamFileSurfaceWriter.C b/src/sampling/sampledSurface/writers/foamFile/foamFileSurfaceWriter.C
index e351afd676..1248f42c15 100644
--- a/src/sampling/sampledSurface/writers/foamFile/foamFileSurfaceWriter.C
+++ b/src/sampling/sampledSurface/writers/foamFile/foamFileSurfaceWriter.C
@@ -80,7 +80,7 @@ Foam::fileName Foam::foamFileSurfaceWriter::write
 
     // Face centers. Not really necessary but very handy when reusing as inputs
     // for e.g. timeVaryingMapped bc.
-    pointField faceCentres(faces.size(),point::zero);
+    pointField faceCentres(faces.size(), Zero);
 
     forAll(faces, faceI)
     {
diff --git a/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.C b/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.C
index 1e3e1fe7ad..baad475656 100644
--- a/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.C
+++ b/src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.C
@@ -323,7 +323,7 @@ Foam::radiation::fvDOM::fvDOM(const volScalarField& T)
     solarLoad_(),
     meshOrientation_
     (
-        coeffs_.lookupOrDefault<vector>("meshOrientation", vector::zero)
+        coeffs_.lookupOrDefault<vector>("meshOrientation", Zero)
     )
 {
     initialise();
@@ -418,7 +418,7 @@ Foam::radiation::fvDOM::fvDOM
     solarLoad_(),
     meshOrientation_
     (
-        coeffs_.lookupOrDefault<vector>("meshOrientation", vector::zero)
+        coeffs_.lookupOrDefault<vector>("meshOrientation", Zero)
     )
 {
     initialise();
diff --git a/src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C b/src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C
index d569d454f5..cb6f8417a2 100644
--- a/src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C
+++ b/src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C
@@ -137,7 +137,7 @@ Foam::radiation::radiativeIntensityRay::radiativeIntensityRay
 
     if (mesh_.nSolutionD() == 2)
     {
-        vector meshDir(vector::zero);
+        vector meshDir(Zero);
         if (dom_.meshOrientation() != vector::zero)
         {
             meshDir = dom_.meshOrientation();
@@ -161,7 +161,7 @@ Foam::radiation::radiativeIntensityRay::radiativeIntensityRay
     }
     else if (mesh_.nSolutionD() == 1)
     {
-        vector meshDir(vector::zero);
+        vector meshDir(Zero);
         if (dom_.meshOrientation() != vector::zero)
         {
             meshDir = dom_.meshOrientation();
@@ -195,7 +195,7 @@ Foam::radiation::radiativeIntensityRay::radiativeIntensityRay
             IOobject::AUTO_WRITE
         );
 
-        // check if field exists and can be read
+        // Check if field exists and can be read
         if (IHeader.typeHeaderOk<volScalarField>(true))
         {
             ILambda_.set
diff --git a/src/thermophysicalModels/radiation/radiationModels/solarLoad/solarLoad.C b/src/thermophysicalModels/radiation/radiationModels/solarLoad/solarLoad.C
index 6a38bb86b4..a69cfbd0eb 100644
--- a/src/thermophysicalModels/radiation/radiationModels/solarLoad/solarLoad.C
+++ b/src/thermophysicalModels/radiation/radiationModels/solarLoad/solarLoad.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2015 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -44,8 +44,8 @@ namespace Foam
     }
 }
 
-const Foam::word Foam::radiation::solarLoad::viewFactorWalls
-    = "viewFactorWall";
+const Foam::word Foam::radiation::solarLoad::viewFactorWalls = "viewFactorWall";
+
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -126,9 +126,7 @@ void Foam::radiation::solarLoad::updateDirectHitRadiation
         label patchID = patches.whichPatch(faceI);
         const polyPatch& pp = patches[patchID];
         const label localFaceI = faceI - pp.start();
-        const vector qPrim =
-              solarCalc_.directSolarRad()
-            * solarCalc_.direction();
+        const vector qPrim = solarCalc_.directSolarRad()*solarCalc_.direction();
 
         if (includeMappedPatchBasePatches[patchID])
         {
@@ -220,9 +218,9 @@ void Foam::radiation::solarLoad::updateSkyDiffusiveRadiation
                         // Ground reflected
                         Er =
                             solarCalc_.directSolarRad()
-                         * (solarCalc_.C() + Foam::sin(solarCalc_.beta()))
-                         * solarCalc_.groundReflectivity()
-                         * (1.0 - cosEpsilon)/2.0;
+                          * (solarCalc_.C() + Foam::sin(solarCalc_.beta()))
+                          * solarCalc_.groundReflectivity()
+                          * (1.0 - cosEpsilon)/2.0;
                     }
 
                     const label cellI = cellIds[faceI];
@@ -241,10 +239,10 @@ void Foam::radiation::solarLoad::updateSkyDiffusiveRadiation
                         for (label bandI = 0; bandI < nBands_; bandI++)
                         {
                             Ru_[cellI] +=
-                                  (Ed + Er)
-                                * spectralDistribution_[bandI]
-                                * absorptivity_[patchID][bandI]()[faceI]
-                                * sf[faceI]/V[cellI];
+                                (Ed + Er)
+                              * spectralDistribution_[bandI]
+                              * absorptivity_[patchID][bandI]()[faceI]
+                              * sf[faceI]/V[cellI];
                         }
                     }
                 }
@@ -337,20 +335,9 @@ void Foam::radiation::solarLoad::initialise(const dictionary& coeffs)
         );
     }
 
-    if (coeffs.found("solidCoupled"))
-    {
-         coeffs.lookup("solidCoupled") >> solidCoupled_;
-    }
-
-    if (coeffs.found("wallCoupled"))
-    {
-         coeffs.lookup("wallCoupled") >> wallCoupled_;
-    }
-
-    if (coeffs.found("updateAbsorptivity"))
-    {
-        coeffs.lookup("updateAbsorptivity") >> updateAbsorptivity_;
-    }
+    coeffs.readIfPresent("solidCoupled", solidCoupled_);
+    coeffs.readIfPresent("wallCoupled", wallCoupled_);
+    coeffs.readIfPresent("updateAbsorptivity", updateAbsorptivity_);
 }
 
 
@@ -399,8 +386,7 @@ void Foam::radiation::solarLoad::calculateQdiff
     forAll(includePatches_, i)
     {
         const label patchI = includePatches_[i];
-        nLocalVFCoarseFaces +=
-            coarseMesh_->boundaryMesh()[patchI].size();
+        nLocalVFCoarseFaces += coarseMesh_->boundaryMesh()[patchI].size();
     }
 
     label totalFVNCoarseFaces = nLocalVFCoarseFaces;
@@ -413,7 +399,7 @@ void Foam::radiation::solarLoad::calculateQdiff
 
     scalarField compactCoarseRave(map_->constructSize(), 0.0);
     scalarField compactCoarsePartialArea(map_->constructSize(), 0.0);
-    vectorList compactCoarseNorm(map_->constructSize(), vector::zero);
+    vectorList compactCoarseNorm(map_->constructSize(), Zero);
 
     const boundaryRadiationProperties& boundaryRadiation =
         boundaryRadiationProperties::New(mesh_);
@@ -443,8 +429,8 @@ void Foam::radiation::solarLoad::calculateQdiff
         for (label bandI = 0; bandI < nBands_; bandI++)
         {
             const tmp<scalarField> tr =
-               spectralDistribution_[bandI]
-              *boundaryRadiation.reflectivity(patchID, bandI);
+                spectralDistribution_[bandI]
+               *boundaryRadiation.reflectivity(patchID, bandI);
             r += tr();
         }
 
@@ -459,14 +445,9 @@ void Foam::radiation::solarLoad::calculateQdiff
         forAll(cpp, coarseI)
         {
             const label coarseFaceID = coarsePatchFace[coarseI];
-            const labelList& fineFaces =
-                coarseToFine_[i][coarseFaceID];
+            const labelList& fineFaces = coarseToFine_[i][coarseFaceID];
 
-            UIndirectList<scalar> fineSf
-            (
-                sf,
-                fineFaces
-            );
+            UIndirectList<scalar> fineSf(sf, fineFaces);
             scalar fineArea = sum(fineSf());
 
             scalar fullArea = 0.0;
@@ -503,18 +484,18 @@ void Foam::radiation::solarLoad::calculateQdiff
     }
 
 
-     SubList<scalar>(compactCoarsePartialArea, nLocalVFCoarseFaces) =
-         localCoarsePartialArea;
+    SubList<scalar>(compactCoarsePartialArea, nLocalVFCoarseFaces) =
+        localCoarsePartialArea;
 
-     SubList<scalar>(compactCoarseRave, nLocalVFCoarseFaces) =
-         localCoarseRave;
+    SubList<scalar>(compactCoarseRave, nLocalVFCoarseFaces) =
+        localCoarseRave;
 
-     SubList<vector>(compactCoarseNorm, nLocalVFCoarseFaces) =
-         localCoarseNorm;
+    SubList<vector>(compactCoarseNorm, nLocalVFCoarseFaces) =
+        localCoarseNorm;
 
-     map_->distribute(compactCoarsePartialArea);
-     map_->distribute(compactCoarseRave);
-     map_->distribute(compactCoarseNorm);
+    map_->distribute(compactCoarsePartialArea);
+    map_->distribute(compactCoarseRave);
+    map_->distribute(compactCoarseNorm);
 
 
     // Calculate coarse hitFaces and Sun direct hit heat fluxes
@@ -527,8 +508,7 @@ void Foam::radiation::solarLoad::calculateQdiff
         const polyPatch& pp = coarseMesh_->boundaryMesh()[patchID];
         const polyPatch& ppf = mesh_.boundaryMesh()[patchID];
 
-        const labelList& coarsePatchFace =
-            coarseMesh_->patchFaceMap()[patchID];
+        const labelList& coarsePatchFace = coarseMesh_->patchFaceMap()[patchID];
         const scalarField& sf = mesh_.magSf().boundaryField()[patchID];
 
         scalarField a(ppf.size(), 0.0);
@@ -544,11 +524,7 @@ void Foam::radiation::solarLoad::calculateQdiff
         {
             const label coarseFaceID = coarsePatchFace[coarseI];
             const labelList& fineFaces = coarseToFine_[i][coarseFaceID];
-            UIndirectList<scalar> fineSf
-            (
-                sf,
-                fineFaces
-            );
+            UIndirectList<scalar> fineSf(sf, fineFaces);
             scalar fineArea = sum(fineSf());
 
             scalar aAve = 0.0;
@@ -568,12 +544,12 @@ void Foam::radiation::solarLoad::calculateQdiff
                 label compactI = compactFaces[j];
 
                 localqDiffusive[locaFaceI] +=
-                      compactCoarsePartialArea[compactI]
-                    * aAve
-                    * (solarCalc_.directSolarRad()*solarCalc_.direction())
-                    & compactCoarseNorm[compactI]
-                    * vf[j]
-                    * compactCoarseRave[compactI];
+                    compactCoarsePartialArea[compactI]
+                  * aAve
+                  * (solarCalc_.directSolarRad()*solarCalc_.direction())
+                  & compactCoarseNorm[compactI]
+                  * vf[j]
+                  * compactCoarseRave[compactI];
 
             }
             locaFaceI++;
@@ -693,7 +669,7 @@ Foam::radiation::solarLoad::solarLoad(const volScalarField& T)
         dimensionedScalar("Ru", dimMass/dimLength/pow3(dimTime), 0.0)
     ),
     solarCalc_(this->subDict(typeName + "Coeffs"), mesh_),
-    verticalDir_(vector::zero),
+    verticalDir_(Zero),
     useVFbeamToDiffuse_(false),
     includePatches_(mesh_.boundary().size(), -1),
     coarseToFine_(),
@@ -783,7 +759,7 @@ Foam::radiation::solarLoad::solarLoad
         dimensionedScalar("Ru", dimMass/dimLength/pow3(dimTime), 0.0)
     ),
     solarCalc_(coeffs_, mesh_),
-    verticalDir_(vector::zero),
+    verticalDir_(Zero),
     useVFbeamToDiffuse_(false),
     includePatches_(mesh_.boundary().size(), -1),
     coarseToFine_(),
@@ -875,7 +851,7 @@ Foam::radiation::solarLoad::solarLoad
         dimensionedScalar("Ru", dimMass/dimLength/pow3(dimTime), 0.0)
     ),
     solarCalc_(dict, mesh_),
-    verticalDir_(vector::zero),
+    verticalDir_(Zero),
     useVFbeamToDiffuse_(false),
     includePatches_(mesh_.boundary().size(), -1),
     coarseToFine_(),
@@ -1034,4 +1010,5 @@ Foam::radiation::solarLoad::Ru() const
     return Ru_;
 }
 
+
 // ************************************************************************* //
diff --git a/src/thermophysicalModels/radiation/submodels/solarCalculator/solarCalculator.C b/src/thermophysicalModels/radiation/submodels/solarCalculator/solarCalculator.C
index 08abcf2853..cc712b4cb2 100644
--- a/src/thermophysicalModels/radiation/submodels/solarCalculator/solarCalculator.C
+++ b/src/thermophysicalModels/radiation/submodels/solarCalculator/solarCalculator.C
@@ -123,7 +123,7 @@ void Foam::solarCalculator::calculateSunDirection()
 
     coord_.reset
     (
-        new coordinateSystem("grid", vector::zero, gridUp_, eastDir_)
+        new coordinateSystem("grid", Zero, gridUp_, eastDir_)
     );
 
     direction_.z() = -sin(beta_);
@@ -236,7 +236,7 @@ Foam::solarCalculator::solarCalculator
 :
     mesh_(mesh),
     dict_(dict),
-    direction_(vector::zero),
+    direction_(Zero),
     directSolarRad_(0.0),
     diffuseSolarRad_(0.0),
     groundReflectivity_(0.0),