ENH: simplify assign zero for Dimensioned and Geometric fields (#3402)

- regular or forced assignment from `zero` does not need dimension
  checking:

    Old:  U = dimensionedVector(U.dimensions(), Zero);
    New:  U = Zero;

  this eliminates a fair bit of clutter and simplifies new coding
  without losing general dimension checking capabilities.

applications/solvers/basic/potentialFoam/createFields.H

@@ -15,11 +15,10 @@ volVectorField U
 // Initialise the velocity internal field to zero
 // Note: explicitly bypass evaluation of contraint patch overrides
 //       (e.g. swirlFanVelocity might lookup phi,rho)
-//U = dimensionedVector(U.dimensions(), Zero);
+//U = Zero;
 {
-    const dimensionedVector dt(U.dimensions(), Zero);
-    U.internalFieldRef() = dt;
-    U.boundaryFieldRef() = dt.value();
+    U.internalFieldRef() = Zero;
+    U.boundaryFieldRef() = Zero;
 }
 
 surfaceScalarField phi

applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C

@@ -1086,7 +1086,7 @@ void Foam::multiphaseMixtureThermo::solveAlphas
 
     MULES::limitSum(alphaPhiCorrs);
 
-    rhoPhi_ = dimensionedScalar(dimensionSet(1, 0, -1, 0, 0), Zero);
+    rhoPhi_ = Zero;
 
     volScalarField sumAlpha
     (

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

@@ -557,7 +557,7 @@ void Foam::radiation::laserDTRM::calculate()
 
 
     // Reset the field
-    Q_ == dimensionedScalar(Q_.dimensions(), Zero);
+    Q_ == Zero;
 
     a_ = absorptionEmission_->a();
     e_ = absorptionEmission_->e();

applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H

@@ -232,7 +232,7 @@
             surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAUf);
 
             phasei = 0;
-            phi = dimensionedScalar("phi", phi.dimensions(), Zero);
+            phi = Zero;
 
             for (phaseModel& phase : fluid.phases())
             {
@@ -261,7 +261,7 @@
 
             mSfGradp = pEqnIncomp.flux()/rAUf;
 
-            U = dimensionedVector("U", dimVelocity, Zero);
+            U = Zero;
 
             phasei = 0;
             for (phaseModel& phase : fluid.phases())

applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C

@@ -626,7 +626,7 @@ void Foam::multiphaseMixture::solveAlphas
 
     MULES::limitSum(alphaPhiCorrs);
 
-    rhoPhi_ = dimensionedScalar(dimMass/dimTime, Zero);
+    rhoPhi_ = Zero;
 
     volScalarField sumAlpha
     (

applications/solvers/stressAnalysis/solidEquilibriumDisplacementFoam/kineticEnergyLimiter.H

@@ -14,6 +14,6 @@ if (!(runTime.timeIndex() % 5))
     if (smi < -SMALL)
     {
         Info<< "Resetting Dcorr to 0" << endl;
-        Dcorr == dimensionedVector(Dcorr.dimensions(), Zero);
+        Dcorr == Zero;
     }
 }

src/OpenFOAM/dimensionSet/dimensionSet.C

@@ -70,7 +70,7 @@ static inline bool checkDims
 
 Foam::dimensionSet::dimensionSet()
 :
-    exponents_(Zero)
+    exponents_(Foam::zero{})
 {}
 
 
@@ -142,7 +142,7 @@ Foam::dimensionSet::values() noexcept
 
 void Foam::dimensionSet::clear()
 {
-    exponents_ = Zero;
+    exponents_ = Foam::zero{};
 }
 
 
@@ -166,13 +166,13 @@ Foam::scalar& Foam::dimensionSet::operator[](const dimensionType type)
 }
 
 
-Foam::scalar Foam::dimensionSet::operator[](const label type) const
+Foam::scalar Foam::dimensionSet::operator[](const int type) const
 {
     return exponents_[type];
 }
 
 
-Foam::scalar& Foam::dimensionSet::operator[](const label type)
+Foam::scalar& Foam::dimensionSet::operator[](const int type)
 {
     return exponents_[type];
 }

src/OpenFOAM/dimensionSet/dimensionSet.H

@@ -329,17 +329,26 @@ public:
         scalar operator[](const dimensionType) const;
         scalar& operator[](const dimensionType);
 
-        scalar operator[](const label) const;
-        scalar& operator[](const label);
+        scalar operator[](const int) const;
+        scalar& operator[](const int);
 
         bool operator==(const dimensionSet&) const;
         bool operator!=(const dimensionSet&) const;
 
+        //- Assignment operation, checks for identical dimensions.
+        //- Use reset() to change the dimensions.
         bool operator=(const dimensionSet&) const;
 
+        //- Addition operation, checks for identical dimensions.
         bool operator+=(const dimensionSet&) const;
+
+        //- Subtraction operation, checks for identical dimensions.
         bool operator-=(const dimensionSet&) const;
+
+        //- Multiplication, modifies the exponents.
         bool operator*=(const dimensionSet&);
+
+        //- Division, modifies the exponents.
         bool operator/=(const dimensionSet&);
 };
 
@@ -414,7 +423,7 @@ dimensionSet invTransform(const dimensionSet& ds);
 //- The dimensionSet inverted
 dimensionSet operator~(const dimensionSet& ds);
 
-//- Unary negation.
+//- Unary negation; does not change the dimensions
 dimensionSet operator-(const dimensionSet& ds);
 
 dimensionSet operator+(const dimensionSet& ds1, const dimensionSet& ds2);

src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedField.C

@@ -560,6 +560,13 @@ void Foam::DimensionedField<Type, GeoMesh>::operator=
 }
 
 
+template<class Type, class GeoMesh>
+void Foam::DimensionedField<Type, GeoMesh>::operator=(Foam::zero)
+{
+    Field<Type>::operator=(Foam::zero{});
+}
+
+
 #define COMPUTED_ASSIGNMENT(TYPE, op)                                          \
                                                                                \
 template<class Type, class GeoMesh>                                            \

src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedField.H

@@ -631,9 +631,12 @@ public:
         void operator=(const DimensionedField<Type, GeoMesh>& df);
         void operator=(const tmp<DimensionedField<Type, GeoMesh>>& tdf);
 
-        //- Assign dimensions and value.
+        //- Assign value (with dimension check)
         void operator=(const dimensioned<Type>& dt);
 
+        //- Assign value zero (no dimension checks)
+        void operator=(Foam::zero);
+
         void operator+=(const DimensionedField<Type, GeoMesh>& df);
         void operator+=(const tmp<DimensionedField<Type, GeoMesh>>& tdf);
 

src/OpenFOAM/fields/FieldFields/FieldField/FieldField.C

@@ -450,9 +450,9 @@ void FieldField<Field, Type>::operator=(const tmp<FieldField>& tf)
 template<template<class> class Field, class Type>
 void FieldField<Field, Type>::operator=(const Type& val)
 {
-    forAll(*this, i)
+    for (auto& pfld : *this)
     {
-        this->operator[](i) = val;
+        pfld = val;
     }
 }
 
@@ -460,9 +460,9 @@ void FieldField<Field, Type>::operator=(const Type& val)
 template<template<class> class Field, class Type>
 void FieldField<Field, Type>::operator=(Foam::zero)
 {
-    forAll(*this, i)
+    for (auto& pfld : *this)
     {
-        this->operator[](i) = Foam::zero{};
+        pfld = Foam::zero{};
     }
 }
 
@@ -489,11 +489,11 @@ void FieldField<Field, Type>::operator op                                      \
 }                                                                              \
                                                                                \
 template<template<class> class Field, class Type>                              \
-void FieldField<Field, Type>::operator op(const TYPE& t)                       \
+void FieldField<Field, Type>::operator op(const TYPE& val)                     \
 {                                                                              \
-    forAll(*this, i)                                                           \
+    for (auto& pfld : *this)                                                   \
     {                                                                          \
-        this->operator[](i) op t;                                              \
+        pfld op val;                                                           \
     }                                                                          \
 }
 

src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricBoundaryField.C

@@ -951,6 +951,16 @@ void Foam::GeometricBoundaryField<Type, PatchField, GeoMesh>::operator=
 }
 
 
+template<class Type, template<class> class PatchField, class GeoMesh>
+void Foam::GeometricBoundaryField<Type, PatchField, GeoMesh>::operator=
+(
+    Foam::zero
+)
+{
+    FieldField<PatchField, Type>::operator=(Foam::zero{});
+}
+
+
 template<class Type, template<class> class PatchField, class GeoMesh>
 void Foam::GeometricBoundaryField<Type, PatchField, GeoMesh>::operator==
 (
@@ -983,9 +993,9 @@ void Foam::GeometricBoundaryField<Type, PatchField, GeoMesh>::operator==
     const Type& val
 )
 {
-    forAll(*this, patchi)
+    for (auto& pfld : *this)
     {
-        this->operator[](patchi) == val;
+        pfld == val;
     }
 }
 

src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricBoundaryField.H

@@ -260,6 +260,9 @@ public:
         //- Assignment to uniform value
         void operator=(const Type& val);
 
+        //- Assignment to zero
+        void operator=(Foam::zero);
+
         //- Forced assignment from GeometricBoundaryField
         void operator==(const GeometricBoundaryField& bf);
 
@@ -269,10 +272,14 @@ public:
         //- Forced assignment to uniform value
         void operator==(const Type& val);
 
+        //- Forced assignment to zero
+        void operator==(Foam::zero) { *this == Type(Foam::zero{}); }
+
         // Prevent automatic comparison rewriting (c++20)
         bool operator!=(const GeometricBoundaryField&) = delete;
         bool operator!=(const FieldField<PatchField, Type>&) = delete;
         bool operator!=(const Type&) = delete;
+        bool operator!=(Foam::zero) = delete;
 };
 
 

src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C

@@ -1478,6 +1478,21 @@ void Foam::GeometricField<Type, PatchField, GeoMesh>::operator=
 }
 
 
+template<class Type, template<class> class PatchField, class GeoMesh>
+void Foam::GeometricField<Type, PatchField, GeoMesh>::operator=(Foam::zero)
+{
+    // No dimension checking
+    primitiveFieldRef() = Foam::zero{};
+    boundaryFieldRef() = Foam::zero{};
+
+    // Make sure any e.g. jump-cyclic are updated.
+    boundaryFieldRef().evaluate_if
+    (
+        [](const auto& pfld) { return pfld.constraintOverride(); }
+    );
+}
+
+
 template<class Type, template<class> class PatchField, class GeoMesh>
 void Foam::GeometricField<Type, PatchField, GeoMesh>::operator==
 (
@@ -1520,6 +1535,21 @@ void Foam::GeometricField<Type, PatchField, GeoMesh>::operator==
 }
 
 
+template<class Type, template<class> class PatchField, class GeoMesh>
+void Foam::GeometricField<Type, PatchField, GeoMesh>::operator==(Foam::zero)
+{
+    // No dimension checking
+    primitiveFieldRef() = Foam::zero{};
+    boundaryFieldRef() == Foam::zero{};
+
+    // Make sure any e.g. jump-cyclic are updated.
+    boundaryFieldRef().evaluate_if
+    (
+        [](const auto& pfld) { return pfld.constraintOverride(); }
+    );
+}
+
+
 #define COMPUTED_ASSIGNMENT(TYPE, op)                                          \
                                                                                \
 template<class Type, template<class> class PatchField, class GeoMesh>          \

src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H

@@ -1028,11 +1028,21 @@ public:
 
         void operator=(const GeometricField<Type, PatchField, GeoMesh>&);
         void operator=(const tmp<GeometricField<Type, PatchField, GeoMesh>>&);
-        void operator=(const dimensioned<Type>&);
+
+        //- Assign value (with dimension check)
+        void operator=(const dimensioned<Type>& dt);
+
+        //- Assign value zero (no dimension checks)
+        void operator=(Foam::zero);
 
         void operator==(const tmp<GeometricField<Type, PatchField, GeoMesh>>&);
+
+        //- Assign value (with dimension check). Force assign for boundary.
         void operator==(const dimensioned<Type>&);
 
+        //- Assign value zero (no dimension checks). Force assign for boundary.
+        void operator==(Foam::zero);
+
         void operator+=(const GeometricField<Type, PatchField, GeoMesh>&);
         void operator+=(const tmp<GeometricField<Type, PatchField, GeoMesh>>&);
 
@@ -1061,6 +1071,7 @@ public:
             const tmp<GeometricField<Type, PatchField, GeoMesh>>&
         ) = delete;
         bool operator!=(const dimensioned<Type>&) = delete;
+        bool operator!=(Foam::zero) = delete;
 
 
     // Write

src/TurbulenceModels/turbulenceModels/LES/LESdeltas/SLADelta/SLADelta.C

@@ -50,7 +50,7 @@ tmp<scalarField> sumNeighbours
 {
     const fvMesh& mesh = field.mesh();
 
-    result == dimensioned<Type>(field.dimensions(), Zero);
+    result == Zero;
 
     auto tscaling = tmp<scalarField>::New(mesh.nCells(), Zero);
     auto& scaling = tscaling.ref();

src/combustionModels/FSD/FSD.C

@@ -307,7 +307,7 @@ void FSD<ReactionThermo, ThermoType>::calculateSourceNorm()
 template<class ReactionThermo, class ThermoType>
 void FSD<ReactionThermo, ThermoType>::correct()
 {
-    this->wFuel_ == dimensionedScalar(dimMass/dimTime/dimVolume, Zero);
+    this->wFuel_ == Zero;
 
     if (this->active())
     {

src/combustionModels/diffusion/diffusion.C

@@ -69,7 +69,7 @@ diffusion<ReactionThermo, ThermoType>::~diffusion()
 template<class ReactionThermo, class ThermoType>
 void diffusion<ReactionThermo, ThermoType>::correct()
 {
-    this->wFuel_ == dimensionedScalar(dimMass/dimVolume/dimTime, Zero);
+    this->wFuel_ == Zero;
 
     if (this->active())
     {

src/combustionModels/diffusionMulticomponent/diffusionMulticomponent.C

@@ -216,15 +216,13 @@ diffusionMulticomponent<ReactionThermo, ThermoType>::correct()
             forAll(lhs, l)
             {
                 const label lIndex = lhs[l].index;
-                this->chemistryPtr_->RR(lIndex) =
-                    dimensionedScalar(dimMass/dimTime/dimVolume, Zero);
+                this->chemistryPtr_->RR(lIndex) = Zero;
             }
 
             forAll(rhs, l)
             {
                 const label rIndex = rhs[l].index;
-                this->chemistryPtr_->RR(rIndex) =
-                    dimensionedScalar(dimMass/dimTime/dimVolume, Zero);
+                this->chemistryPtr_->RR(rIndex) = Zero;
             }
         }
 

src/combustionModels/eddyDissipationModelBase/eddyDissipationModelBase.C

@@ -81,7 +81,7 @@ eddyDissipationModelBase<ReactionThermo, ThermoType>::rtTurb() const
 template<class ReactionThermo, class ThermoType>
 void eddyDissipationModelBase<ReactionThermo, ThermoType>::correct()
 {
-    this->wFuel_ == dimensionedScalar(dimMass/dimVolume/dimTime, Zero);
+    this->wFuel_ == Zero;
 
     if (this->active())
     {

src/combustionModels/infinitelyFastChemistry/infinitelyFastChemistry.C

@@ -67,7 +67,7 @@ infinitelyFastChemistry<ReactionThermo, ThermoType>::~infinitelyFastChemistry()
 template<class ReactionThermo, class ThermoType>
 void infinitelyFastChemistry<ReactionThermo, ThermoType>::correct()
 {
-    this->wFuel_ == dimensionedScalar(dimMass/dimVolume/dimTime, Zero);
+    this->wFuel_ == Zero;
 
     if (this->active())
     {

src/finiteArea/faMatrices/faMatrix/faMatrix.H

@@ -510,8 +510,8 @@ public:
         void operator+=(const dimensioned<Type>&);
         void operator-=(const dimensioned<Type>&);
 
-        void operator+=(const Foam::zero) {}
-        void operator-=(const Foam::zero) {}
+        void operator+=(Foam::zero) {}
+        void operator-=(Foam::zero) {}
 
         void operator*=(const areaScalarField::Internal&);
         void operator*=(const tmp<areaScalarField::Internal>&);

src/finiteVolume/cfdTools/compressible/createDpdt.H

@@ -27,6 +27,6 @@ volScalarField dpdt(dpdtHeader, fvc::ddt(p));
 
 if (!thermo.dpdt())
 {
-    dpdt == dimensionedScalar(dpdt.dimensions(), Zero);
+    dpdt == Zero;
     dpdt.writeOpt(IOobject::NO_WRITE);
 }

src/finiteVolume/finiteVolume/snGradSchemes/skewCorrectedSnGrad/skewCorrectedSnGrad.C

@@ -60,7 +60,7 @@ Foam::fv::skewCorrectedSnGrad<Type>::fullGradCorrection
     auto& ssf = tssf.ref();
 
     ssf.setOriented();
-    ssf = dimensioned<Type>(ssf.dimensions(), Zero);
+    ssf = Zero;
 
 
     typedef typename

src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H

@@ -693,8 +693,8 @@ public:
         void operator+=(const dimensioned<Type>&);
         void operator-=(const dimensioned<Type>&);
 
-        void operator+=(const Foam::zero) {}
-        void operator-=(const Foam::zero) {}
+        void operator+=(Foam::zero) {}
+        void operator-=(Foam::zero) {}
 
         void operator*=(const volScalarField::Internal&);
         void operator*=(const tmp<volScalarField::Internal>&);

src/finiteVolume/fvMesh/fvMeshGeometry.C

@@ -427,7 +427,7 @@ const Foam::surfaceScalarField& Foam::fvMesh::phi() const
     // mesh motion fluxes if the time has been incremented
     if (!time().subCycling() && phiPtr_->timeIndex() != time().timeIndex())
     {
-        (*phiPtr_) = dimensionedScalar(dimVolume/dimTime, Foam::zero{});
+        (*phiPtr_) = Zero;
     }
 
     phiPtr_->setOriented();

src/finiteVolume/interpolation/surfaceInterpolation/schemes/zoneBlended/zoneBlended.H

@@ -310,7 +310,7 @@ public:
             auto tcorr =
                 SurfaceField::New("correction", vf.mesh(), vf.dimensions());
             auto& corr = tcorr.ref();
-            corr = dimensioned<Type>(vf.dimensions(), Zero);
+            corr = Zero;
 
             // Set default scheme correction
             const auto& scheme0 = this->scheme(0);

src/functionObjects/forces/forceCoeffs/forceCoeffs.C

@@ -122,8 +122,8 @@ void Foam::functionObjects::forceCoeffs::reset()
     Cf_.reset();
     Cm_.reset();
 
-    forceCoeff() == dimensionedVector(dimless, Zero);
-    momentCoeff() == dimensionedVector(dimless, Zero);
+    forceCoeff() == Zero;
+    momentCoeff() == Zero;
 }
 
 

src/functionObjects/forces/forces/forces.C

@@ -200,8 +200,8 @@ void Foam::functionObjects::forces::reset()
 
     if (porosity_)
     {
-        force == dimensionedVector(force.dimensions(), Zero);
-        moment == dimensionedVector(moment.dimensions(), Zero);
+        force == Zero;
+        moment == Zero;
     }
     else
     {

src/functionObjects/solvers/scalarTransport/scalarTransport.C

@@ -195,7 +195,7 @@ Foam::functionObjects::scalarTransport::scalarTransport
 
     if (resetOnStartUp_)
     {
-        s == dimensionedScalar(dimless, Zero);
+        s == Zero;
     }
 }
 

src/fused/finiteVolume/fusedGaussLaplacianScheme.C

@@ -188,7 +188,7 @@ void fusedGaussLaplacianScheme<Type, GType>::gradComponent
     GeometricField<Type, fvPatchField, volMesh>& gGrad
 )
 {
-    gGrad = dimensioned<Type>(vf.dimensions()/dimLength, Zero);
+    gGrad = Zero;
 
     // Calculate grad of vf.component(cmpt)
     fvc::GaussOp

src/fvMotionSolver/motionInterpolation/patchCorrected/patchCorrectedInterpolationTemplates.C

@@ -109,7 +109,7 @@ void Foam::patchCorrectedInterpolation::interpolateDataFromPatchGroups
         dimensionedScalar(dimless, Zero),
         pointPatchFieldBase::zeroGradientType()
     );
-    data = dimensioned<Type>(data.dimensions(), Zero);
+    data = Zero;
 
     forAll(patchGroups_, patchGroupI)
     {

src/optimisation/adjointOptimisation/adjoint/objectives/incompressible/objectiveIncompressible/objectiveIncompressible.C

@@ -240,29 +240,27 @@ void objectiveIncompressible::nullify()
     {
         if (hasdJdv())
         {
-            dJdvPtr_() == dimensionedVector(dJdvPtr_().dimensions(), Zero);
+            dJdvPtr_() == Zero;
         }
         if (hasdJdp())
         {
-            dJdpPtr_() == dimensionedScalar(dJdpPtr_().dimensions(), Zero);
+            dJdpPtr_() == Zero;
         }
         if (hasdJdT())
         {
-            dJdTPtr_() == dimensionedScalar(dJdTPtr_().dimensions(), Zero);
+            dJdTPtr_() == Zero;
         }
         if (hasdJdTMVar1())
         {
-            dJdTMvar1Ptr_() ==
-                dimensionedScalar(dJdTMvar1Ptr_().dimensions(), Zero);
+            dJdTMvar1Ptr_() == Zero;
         }
         if (hasdJdTMVar2())
         {
-            dJdTMvar2Ptr_() ==
-                dimensionedScalar(dJdTMvar2Ptr_().dimensions(), Zero);
+            dJdTMvar2Ptr_() == Zero;
         }
         if (hasBoundarydJdv())
         {
-            bdJdvPtr_() == vector::zero;
+            bdJdvPtr_() == Zero;
         }
         if (hasBoundarydJdvn())
         {
@@ -270,11 +268,11 @@ void objectiveIncompressible::nullify()
         }
         if (hasBoundarydJdvt())
         {
-            bdJdvtPtr_() == vector::zero;
+            bdJdvtPtr_() == Zero;
         }
         if (hasBoundarydJdp())
         {
-            bdJdpPtr_() == vector::zero;
+            bdJdpPtr_() == Zero;
         }
         if (hasBoundarydJdT())
         {
@@ -294,7 +292,7 @@ void objectiveIncompressible::nullify()
         }
         if (hasBoundarydJdGradU())
         {
-            bdJdGradUPtr_() == tensor::zero;
+            bdJdGradUPtr_() == Zero;
         }
 
         // Nullify geometric fields and sets nullified_ to true

src/optimisation/adjointOptimisation/adjoint/objectives/objective/objective.C

@@ -484,7 +484,7 @@ void objective::nullify()
     {
         if (hasdJdb())
         {
-            dJdbPtr_() == dimensionedScalar(dJdbPtr_().dimensions(), Zero);
+            dJdbPtr_() == Zero;
         }
         if (hasdJdbField())
         {
@@ -492,40 +492,38 @@ void objective::nullify()
         }
         if (hasBoundarydJdb())
         {
-            bdJdbPtr_() == vector::zero;
+            bdJdbPtr_() == Zero;
         }
         if (hasdSdbMult())
         {
-            bdSdbMultPtr_() == vector::zero;
+            bdSdbMultPtr_() == Zero;
         }
         if (hasdndbMult())
         {
-            bdndbMultPtr_() == vector::zero;
+            bdndbMultPtr_() == Zero;
         }
         if (hasdxdbMult())
         {
-            bdxdbMultPtr_() == vector::zero;
+            bdxdbMultPtr_() == Zero;
         }
         if (hasdxdbDirectMult())
         {
-            bdxdbDirectMultPtr_() == vector::zero;
+            bdxdbDirectMultPtr_() == Zero;
         }
         if (hasBoundaryEdgeContribution())
         {
-            for (Field<vectorField>& field : bEdgeContribution_())
+            for (auto& field : bEdgeContribution_())
             {
-                field = vector::zero;
+                field = Zero;
             }
         }
         if (hasDivDxDbMult())
         {
-            divDxDbMultPtr_() ==
-                dimensionedScalar(divDxDbMultPtr_().dimensions(), Zero);
+            divDxDbMultPtr_() == Zero;
         }
         if (hasGradDxDbMult())
         {
-            gradDxDbMultPtr_() ==
-                dimensionedTensor(gradDxDbMultPtr_().dimensions(), Zero);
+            gradDxDbMultPtr_() == Zero;
         }
 
         nullified_ = true;

src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/adjointSensitivity/adjointEikonalSolver/adjointEikonalSolver.C

@@ -266,8 +266,8 @@ void adjointEikonalSolver::solve()
 
 void adjointEikonalSolver::reset()
 {
-    source_ == dimensionedScalar(source_.dimensions(), Zero);
-    distanceSensPtr_() = vector::zero;
+    source_ == Zero;
+    distanceSensPtr_() = Zero;
 }
 
 

src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/adjointSensitivity/shape/adjointMeshMovementSolver/adjointMeshMovementSolver.C

@@ -187,7 +187,7 @@ void adjointMeshMovementSolver::solve()
 
 void adjointMeshMovementSolver::reset()
 {
-    source_ == dimensionedVector(source_.dimensions(), Zero);
+    source_ == Zero;
     meshMovementSensPtr_() = Zero;
 }
 

src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/adjointSensitivity/shape/shapeSensitivityBase/ShapeSensitivitiesBase.C

@@ -237,7 +237,7 @@ void Foam::ShapeSensitivitiesBase::allocateMultipliers()
 
 void Foam::ShapeSensitivitiesBase::clearMultipliers()
 {
-    gradDxDbMult_() = dimensionedTensor(gradDxDbMult_().dimensions(), Zero);
+    gradDxDbMult_() = Zero;
     if (divDxDbMult_)
     {
         divDxDbMult_() = Zero;

src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressible/incompressibleVars.C

@@ -402,9 +402,9 @@ void incompressibleVars::resetMeanFields()
         Info<< "Resetting mean fields to zero" << endl;
 
         // Reset fields to zero
-        pMeanPtr_() == dimensionedScalar(pInst().dimensions(), Zero);
-        UMeanPtr_() == dimensionedVector(UInst().dimensions(), Zero);
-        phiMeanPtr_() == dimensionedScalar(phiInst().dimensions(), Zero);
+        pMeanPtr_() == Zero;
+        UMeanPtr_() == Zero;
+        phiMeanPtr_() == Zero;
         RASModelVariables_().resetMeanFields();
 
         // Reset averaging iteration index to 0

src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/incompressibleAdjoint/incompressibleAdjointVars.C

@@ -72,9 +72,9 @@ void incompressibleAdjointVars::restoreInitValues()
     if (solverControl_.storeInitValues())
     {
         Info<< "Restoring adjoint field to initial ones" << endl;
-        paInst() == dimensionedScalar(paInst().dimensions(), Zero);
-        UaInst() == dimensionedVector(UaInst().dimensions(), Zero);
-        phiaInst() == dimensionedScalar(phiaInst().dimensions(), Zero);
+        paInst() == Zero;
+        UaInst() == Zero;
+        phiaInst() == Zero;
         adjointTurbulence_().restoreInitValues();
     }
 }
@@ -87,9 +87,9 @@ void incompressibleAdjointVars::resetMeanFields()
         Info<< "Resetting adjoint mean fields to zero" << endl;
 
         // Reset fields to zero
-        paMeanPtr_() == dimensionedScalar(paPtr_().dimensions(), Zero);
-        UaMeanPtr_() == dimensionedVector(UaPtr_().dimensions(), Zero);
-        phiaMeanPtr_() == dimensionedScalar(phiaPtr_().dimensions(), Zero);
+        paMeanPtr_() == Zero;
+        UaMeanPtr_() == Zero;
+        phiaMeanPtr_() == Zero;
         adjointTurbulence_().resetMeanFields();
 
         // Reset averaging iteration index to 0

src/optimisation/adjointOptimisation/adjoint/solvers/variablesSet/variablesSet/variablesSetTemplates.C

@@ -348,11 +348,10 @@ void variablesSet::nullifyField
     GeometricField<Type, PatchField, GeoMesh>& field
 )
 {
-    typedef GeometricField<Type, PatchField, GeoMesh> fieldType;
-    field == dimensioned<Type>(field.dimensions(), Zero);
+    field == Zero;
     if (field.nOldTimes())
     {
-        fieldType& oldTime = field.oldTime();
+        auto& oldTime = field.oldTime();
         variablesSet::nullifyField(oldTime);
     }
 }

src/optimisation/adjointOptimisation/adjoint/turbulenceModels/incompressibleAdjoint/adjointRAS/adjointRASModel/adjointRASModel.C

@@ -68,14 +68,12 @@ void adjointRASModel::restoreInitValues()
     {
         if (adjointTMVariable1Ptr_)
         {
-            volScalarField& var1 = adjointTMVariable1Ptr_.ref();
-            var1 == dimensionedScalar(var1.dimensions(), Zero);
+            adjointTMVariable1Ptr_.ref() == Zero;
         }
 
         if (adjointTMVariable2Ptr_)
         {
-            volScalarField& var2 = adjointTMVariable2Ptr_.ref();
-            var2 == dimensionedScalar(var2.dimensions(), Zero);
+            adjointTMVariable2Ptr_.ref() == Zero;
         }
     }
 }
@@ -465,13 +463,11 @@ void adjointRASModel::resetMeanFields()
     {
         if (adjointTMVariable1MeanPtr_)
         {
-            adjointTMVariable1MeanPtr_() ==
-                dimensionedScalar(adjointTMVariable1Ptr_().dimensions(), Zero);
+            adjointTMVariable1MeanPtr_.ref() == Zero;
         }
         if (adjointTMVariable2MeanPtr_)
         {
-            adjointTMVariable2MeanPtr_() ==
-                dimensionedScalar(adjointTMVariable2Ptr_().dimensions(), Zero);
+            adjointTMVariable2MeanPtr_.ref() == Zero;
         }
     }
 }

src/optimisation/adjointOptimisation/adjoint/turbulenceModels/turbulenceModelVariables/RAS/RASModelVariables/RASModelVariables.C

@@ -363,18 +363,15 @@ void RASModelVariables::resetMeanFields()
         // Reset fields to zero
         if (TMVar1MeanPtr_)
         {
-            TMVar1MeanPtr_.ref() ==
-                dimensionedScalar(TMVar1Inst().dimensions(), Zero);
+            TMVar1MeanPtr_.ref() == Zero;
         }
         if (TMVar2MeanPtr_)
         {
-            TMVar2MeanPtr_.ref() ==
-                dimensionedScalar(TMVar2Inst().dimensions(), Zero);
+            TMVar2MeanPtr_.ref() == Zero;
         }
         if (nutMeanPtr_)
         {
-            nutMeanPtr_.ref() ==
-                dimensionedScalar(nutRefInst().dimensions(), Zero);
+            nutMeanPtr_.ref() == Zero;
         }
     }
 }

src/phaseSystemModels/multiphaseInter/phasesSystem/multiphaseSystem/multiphaseSystem.C

@@ -262,8 +262,8 @@ void Foam::multiphaseInter::multiphaseSystem::solveAlphas()
         // Set Su and Sp to zero
         for (const multiphaseInter::phaseModel& phase : phases_)
         {
-            Su_[phase.name()] = dimensionedScalar("Su", dimless/dimTime, Zero);
-            Sp_[phase.name()] = dimensionedScalar("Sp", dimless/dimTime, Zero);
+            Su_[phase.name()] = Zero;
+            Sp_[phase.name()] = Zero;
 
             // Add alpha*div(U)
             //const volScalarField& alpha = phase;
@@ -368,7 +368,7 @@ void Foam::multiphaseInter::multiphaseSystem::solveAlphas()
             );
 
             // Reset rhoPhi
-            rhoPhi_ = dimensionedScalar("rhoPhi", dimMass/dimTime, Zero);
+            rhoPhi_ = Zero;
 
             for (multiphaseInter::phaseModel& phase : phases_)
             {

src/phaseSystemModels/multiphaseInter/phasesSystem/phaseModel/StaticPhaseModel/StaticPhaseModel.C

@@ -101,7 +101,7 @@ template<class BasePhaseModel>
 const Foam::surfaceScalarField&
 Foam::StaticPhaseModel<BasePhaseModel>::phi()
 {
-    phi_ = dimensionedScalar(dimensionSet(0, 3, -1, 0, 0), Zero);
+    phi_ = Zero;
     return phi_;
 }
 
@@ -128,7 +128,7 @@ template<class BasePhaseModel>
 Foam::surfaceScalarField&
 Foam::StaticPhaseModel<BasePhaseModel>::alphaPhi()
 {
-    alphaPhi_ = dimensionedScalar(dimensionSet(0, 3, -1, 0, 0), Zero);
+    alphaPhi_ = Zero;
     return alphaPhi_;
 }
 

src/phaseSystemModels/reactingEuler/multiphaseSystem/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C

@@ -418,7 +418,7 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::correctInterfaceThermo()
         }
         else
         {
-            iDmdtNew == dimensionedScalar(iDmdt.dimensions());
+            iDmdtNew == Zero;
         }
 
         volScalarField H1(heatTransferModelIter().first()->K());

src/phaseSystemModels/reactingEuler/twoPhaseCompressibleTurbulenceModels/phasePressureModel/phasePressureModel.C

@@ -63,7 +63,7 @@ Foam::RASModels::phasePressureModel::phasePressureModel
     expMax_(coeffDict_.get<scalar>("expMax")),
     g0_("g0", dimPressure, coeffDict_)
 {
-    nut_ == dimensionedScalar(nut_.dimensions());
+    nut_ == Zero;
 
     if (type == typeName)
     {
@@ -205,11 +205,8 @@ Foam::RASModels::phasePressureModel::devRhoReff
         IOobject::groupName("devRhoReff", U.group()),
         IOobject::NO_REGISTER,
         mesh_,
-        dimensioned<symmTensor>
-        (
-            rho_.dimensions()*dimensionSet(0, 2, -2, 0, 0),
-            Zero
-        )
+        Zero,
+        rho_.dimensions()*dimensionSet(0, 2, -2, 0, 0)
     );
 }
 

src/phaseSystemModels/reactingEuler/twoPhaseSystem/diameterModels/IATE/IATEsources/dummy/dummy.C

@@ -61,7 +61,7 @@ Foam::diameterModels::IATEsources::dummy::R
             iate_.phase().mesh()
         ),
         iate_.phase().mesh(),
-        dimensionedScalar(kappai.dimensions()/dimTime, 0)
+        dimensionedScalar(kappai.dimensions()/dimTime, Zero)
     );
 
     return fvm::Su(R, kappai);

src/phaseSystemModels/twoPhaseEuler/phaseCompressibleTurbulenceModels/phasePressureModel/phasePressureModel.C

@@ -66,7 +66,7 @@ Foam::RASModels::phasePressureModel::phasePressureModel
     expMax_(coeffDict_.get<scalar>("expMax")),
     g0_("g0", dimPressure, coeffDict_)
 {
-    nut_ == dimensionedScalar(nut_.dimensions(), Zero);
+    nut_ == Zero;
 
     if (type == typeName)
     {
@@ -213,10 +213,8 @@ Foam::RASModels::phasePressureModel::devRhoReff
         IOobject::groupName("devRhoReff", U.group()),
         IOobject::NO_REGISTER,
         mesh_,
-        dimensioned<symmTensor>
-        (
-            rho_.dimensions()*dimensionSet(0, 2, -2, 0, 0), Zero
-        )
+        Zero,
+        rho_.dimensions()*dimensionSet(0, 2, -2, 0, 0)
     );
 }
 

src/regionFaModels/liquidFilm/liquidFilmModel/liquidFilmModel.C

@@ -266,8 +266,8 @@ void liquidFilmModel::preEvolveRegion()
     liquidFilmBase::preEvolveRegion();
 
 
-    cloudMassTrans_ == dimensionedScalar(dimMass, Zero);
-    cloudDiameterTrans_ == dimensionedScalar(dimLength, Zero);
+    cloudMassTrans_ == Zero;
+    cloudDiameterTrans_ == Zero;
 
     const scalar deltaT = primaryMesh().time().deltaTValue();
     const scalarField rAreaDeltaT(scalar(1)/deltaT/regionMesh().S().field());

src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C

@@ -92,7 +92,7 @@ bool reactingOneDim::read(const dictionary& dict)
 void reactingOneDim::updateqr()
 {
     // Update local qr from coupled qr field
-    qr_ == dimensionedScalar(qr_.dimensions(), Zero);
+    qr_ == Zero;
 
     // Retrieve field from coupled region using mapped boundary conditions
     qr_.correctBoundaryConditions();
@@ -143,8 +143,8 @@ void reactingOneDim::updateqr()
 
 void reactingOneDim::updatePhiGas()
 {
-    phiHsGas_ == dimensionedScalar(phiHsGas_.dimensions(), Zero);
-    phiGas_ == dimensionedScalar(phiGas_.dimensions(), Zero);
+    phiHsGas_ == Zero;
+    phiGas_ == Zero;
 
     const speciesTable& gasTable = solidChemistry_->gasTable();
 

src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C

@@ -85,9 +85,9 @@ void kinematicSingleLayer::resetPrimaryRegionSourceTerms()
 {
     DebugInFunction << endl;
 
-    rhoSpPrimary_ == dimensionedScalar(rhoSp_.dimensions(), Zero);
-    USpPrimary_ == dimensionedVector(USp_.dimensions(), Zero);
-    pSpPrimary_ == dimensionedScalar(pSp_.dimensions(), Zero);
+    rhoSpPrimary_ == Zero;
+    USpPrimary_ == Zero;
+    pSpPrimary_ == Zero;
 }
 
 
@@ -849,9 +849,9 @@ void kinematicSingleLayer::preEvolveRegion()
 
     // Reset transfer fields
     availableMass_ = mass();
-    cloudMassTrans_ == dimensionedScalar(dimMass, Zero);
-    cloudDiameterTrans_ == dimensionedScalar(dimLength, Zero);
-    primaryMassTrans_ == dimensionedScalar(dimMass, Zero);
+    cloudMassTrans_ == Zero;
+    cloudDiameterTrans_ == Zero;
+    primaryMassTrans_ == Zero;
 }
 
 

src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.C

@@ -98,7 +98,7 @@ void thermoSingleLayer::resetPrimaryRegionSourceTerms()
 
     kinematicSingleLayer::resetPrimaryRegionSourceTerms();
 
-    hsSpPrimary_ == dimensionedScalar(hsSp_.dimensions(), Zero);
+    hsSpPrimary_ == Zero;
 }
 
 
@@ -607,7 +607,7 @@ void thermoSingleLayer::preEvolveRegion()
     DebugInFunction << endl;
 
     kinematicSingleLayer::preEvolveRegion();
-    primaryEnergyTrans_ == dimensionedScalar(dimEnergy, Zero);
+    primaryEnergyTrans_ == Zero;
 }
 
 

src/thermophysicalModels/radiation/radiationModels/fvDOM/fvDOM/fvDOM.C

@@ -725,10 +725,10 @@ void Foam::radiation::fvDOM::updateBlackBodyEmission()
 
 void Foam::radiation::fvDOM::updateG()
 {
-    G_ = dimensionedScalar(dimMass/pow3(dimTime), Zero);
-    qr_ = dimensionedScalar(dimMass/pow3(dimTime), Zero);
-    qem_ = dimensionedScalar(dimMass/pow3(dimTime), Zero);
-    qin_ = dimensionedScalar(dimMass/pow3(dimTime), Zero);
+    G_ = Zero;
+    qr_ = Zero;
+    qem_ = Zero;
+    qin_ = Zero;
 
     forAll(IRay_, rayI)
     {

src/thermophysicalModels/radiation/radiationModels/fvDOM/radiativeIntensityRay/radiativeIntensityRay.C

@@ -300,7 +300,7 @@ Foam::scalar Foam::radiation::radiativeIntensityRay::correct()
 
 void Foam::radiation::radiativeIntensityRay::addIntensity()
 {
-    I_ = dimensionedScalar(dimMass/pow3(dimTime), Zero);
+    I_ = Zero;
 
     forAll(ILambda_, lambdaI)
     {

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

@@ -868,7 +868,7 @@ void Foam::radiation::solarLoad::calculate()
     if (firstIter_ || facesChanged || timeDependentLoad)
     {
         // Reset Ru
-        Ru_ = dimensionedScalar("Ru", dimMass/dimLength/pow3(dimTime), Zero);
+        Ru_ = Zero;
 
         solarCalc_.correctDirectSolarRad();
         solarCalc_.correctDiffuseSolarRad();

src/thermophysicalModels/radiation/submodels/absorptionEmissionModel/wideBandAbsorptionEmission/wideBandAbsorptionEmission.C

@@ -328,7 +328,7 @@ void Foam::radiation::wideBandAbsorptionEmission::correct
     PtrList<volScalarField>& aLambda
 ) const
 {
-    a = dimensionedScalar(dimless/dimLength, Zero);
+    a = Zero;
 
     for (label j=0; j<nBands_; j++)
     {

src/thermophysicalModels/solidThermo/solidThermo/heSolidThermo.C

@@ -151,8 +151,8 @@ heSolidThermo
     heThermo<BasicSolidThermo, MixtureType>(mesh, phaseName)
 {
     calculate();
-    this->mu_ == dimensionedScalar(this->mu_.dimensions(), Zero);
-    this->psi_ == dimensionedScalar(this->psi_.dimensions(), Zero);
+    this->mu_ == Zero;
+    this->psi_ == Zero;
 }
 
 
@@ -168,8 +168,8 @@ heSolidThermo
     heThermo<BasicSolidThermo, MixtureType>(mesh, dict, phaseName)
 {
     calculate();
-    this->mu_ == dimensionedScalar(this->mu_.dimensions(), Zero);
-    this->psi_ == dimensionedScalar(this->psi_.dimensions(), Zero);
+    this->mu_ == Zero;
+    this->psi_ == Zero;
 }
 
 
@@ -188,8 +188,8 @@ heSolidThermo
 
     // TBD. initialise psi, mu (at heThermo level) since these do not
     // get initialised. Move to heThermo constructor?
-    this->mu_ == dimensionedScalar(this->mu_.dimensions(), Zero);
-    this->psi_ == dimensionedScalar(this->psi_.dimensions(), Zero);
+    this->mu_ == Zero;
+    this->psi_ == Zero;
 }
 
 

src/transportModels/geometricVoF/advectionSchemes/isoAdvection/isoAdvectionTemplates.C

@@ -150,7 +150,7 @@ void Foam::isoAdvection::limitFluxes
             DebugInfo << "boundAlpha... " << endl;
 
             DynamicList<label> correctedFaces(3*nOvershoots);
-            dVfcorrectionValues = dimensionedScalar("0",dimVolume,0.0);
+            dVfcorrectionValues = Zero;
             boundFlux(needBounding, dVfcorrectionValues, correctedFaces,Sp,Su);
 
             correctedFaces.append

src/transportModels/geometricVoF/reconstructionSchemes/plicSchemes/gradAlpha/gradAlpha.C

@@ -151,8 +151,8 @@ void Foam::reconstruction::gradAlpha::reconstruct(bool forceUpdate)
         }
     }
     interfaceNormal_.resize(interfaceLabels_.size());
-    centre_ = dimensionedVector("centre", dimLength, Zero);
-    normal_ = dimensionedVector("normal", dimArea, Zero);
+    centre_ = Zero;
+    normal_ = Zero;
 
     gradSurf(alpha1_);
 

src/transportModels/geometricVoF/reconstructionSchemes/plicSchemes/plicRDF/plicRDF.C

@@ -372,8 +372,8 @@ Foam::reconstruction::plicRDF::plicRDF
 {
     setInitNormals(false);
 
-    centre_ = dimensionedVector("centre", dimLength, Zero);
-    normal_ = dimensionedVector("normal", dimArea, Zero);
+    centre_ = Zero;
+    normal_ = Zero;
 
     forAll(interfaceLabels_, i)
     {
@@ -436,8 +436,8 @@ void Foam::reconstruction::plicRDF::reconstruct(bool forceUpdate)
     // Sets interfaceCell_ and interfaceNormal
     setInitNormals(interpolateNormal_);
 
-    centre_ = dimensionedVector("centre", dimLength, Zero);
-    normal_ = dimensionedVector("normal", dimArea, Zero);
+    centre_ = Zero;
+    normal_ = Zero;
 
     // nextToInterface is update on setInitNormals
     const boolList& nextToInterface_ = RDF_.nextToInterface();