Merge branch 'olesenm'

applications/solvers/DNS/dnsFoam/dnsFoam.C

@@ -85,7 +85,7 @@ int main(int argc, char *argv[])
 
         for (int corr=1; corr<=1; corr++)
         {
-            volScalarField rAU = 1.0/UEqn.A();
+            volScalarField rAU(1.0/UEqn.A());
 
             U = rAU*UEqn.H();
             phi = (fvc::interpolate(U) & mesh.Sf())

applications/solvers/DNS/dnsFoam/readTurbulenceProperties.H

@@ -12,8 +12,10 @@
         )
     );
 
-    volVectorField force =
+    volVectorField force
-        U/dimensionedScalar("dt", dimTime, runTime.deltaTValue());
+    (
+        U/dimensionedScalar("dt", dimTime, runTime.deltaTValue())
+    );
 
     Kmesh K(mesh);
     UOprocess forceGen(K, runTime.deltaTValue(), turbulenceProperties);

applications/solvers/basic/laplacianFoam/write.H

@@ -1,6 +1,6 @@
     if (runTime.outputTime())
     {
-        volVectorField gradT = fvc::grad(T);
+        volVectorField gradT(fvc::grad(T));
 
         volScalarField gradTx
         (

applications/solvers/combustion/PDRFoam/PDRModels/XiEqModels/basicXiSubXiEq/basicXiSubXiEq.C

@@ -34,8 +34,8 @@ namespace XiEqModels
 {
     defineTypeNameAndDebug(basicSubGrid, 0);
     addToRunTimeSelectionTable(XiEqModel, basicSubGrid, dictionary);
-};
+}
-};
+}
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -123,20 +123,24 @@ Foam::tmp<Foam::volScalarField> Foam::XiEqModels::basicSubGrid::XiEq() const
     const objectRegistry& db = Su_.db();
     const volVectorField& U = db.lookupObject<volVectorField>("U");
 
-    volScalarField magU = mag(U);
+    volScalarField magU(mag(U));
-    volVectorField Uhat =
+    volVectorField Uhat
-    U/(mag(U) + dimensionedScalar("Usmall", U.dimensions(), 1e-4));
+    (
+        U/(mag(U) + dimensionedScalar("Usmall", U.dimensions(), 1e-4))
+    );
 
-    volScalarField n = max(N_ - (Uhat & ns_ & Uhat), scalar(1e-4));
+    volScalarField n(max(N_ - (Uhat & ns_ & Uhat), scalar(1e-4)));
 
-    volScalarField b = (Uhat & B_ & Uhat)/n;
+    volScalarField b((Uhat & B_ & Uhat)/n);
 
-    volScalarField up = sqrt((2.0/3.0)*turbulence_.k());
+    volScalarField up(sqrt((2.0/3.0)*turbulence_.k()));
 
-    volScalarField XiSubEq =
+    volScalarField XiSubEq
+    (
         scalar(1)
       + max(2.2*sqrt(b), min(0.34*magU/up, scalar(1.6)))
-       *min(0.25*n, scalar(1));
+       *min(0.25*n, scalar(1))
+    );
 
     return XiSubEq*XiEqModel_->XiEq();
 }

applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.C

@@ -78,7 +78,7 @@ void PDRkEpsilon::correct()
 
     RASModel::correct();
 
-    volScalarField divU = fvc::div(phi_/fvc::interpolate(rho_));
+    volScalarField divU(fvc::div(phi_/fvc::interpolate(rho_)));
 
     if (mesh_.moving())
     {
@@ -99,7 +99,7 @@ void PDRkEpsilon::correct()
     const PDRDragModel& drag =
         U_.db().lookupObject<PDRDragModel>("PDRDragModel");
 
-    volScalarField GR = drag.Gk();
+    volScalarField GR(drag.Gk());
 
     // Dissipation equation
     tmp<fvScalarMatrix> epsEqn

applications/solvers/combustion/PDRFoam/StCourantNo.H

@@ -34,9 +34,11 @@ Description
 
     if (mesh.nInternalFaces())
     {
-        scalarField sumPhi =
+        scalarField sumPhi
+        (
             fvc::surfaceSum(mag(phiSt))().internalField()
-           /rho.internalField();
+          / rho.internalField()
+        );
 
         StCoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
 

applications/solvers/combustion/PDRFoam/UEqn.H

@@ -7,7 +7,7 @@
         betav*rho*g
     );
 
-    volSymmTensorField invA = inv(I*UEqn.A() + drag->Dcu());
+    volSymmTensorField invA(inv(I*UEqn.A() + drag->Dcu()));
 
     if (momentumPredictor)
     {

applications/solvers/combustion/PDRFoam/XiEqns

@@ -1,91 +0,0 @@
-    // Calculate Xi according to the selected flame wrinkling model
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    // Calculate coefficients for Gulder's flame speed correlation
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    volScalarField up = uPrimeCoef*sqrt((2.0/3.0)*k);
-    volScalarField epsilonPlus = pow(uPrimeCoef, 3)*epsilon;
-
-    volScalarField tauEta = sqrt(mag(thermo->muu()/(rhou*epsilonPlus)));
-    volScalarField Reta = up/
-    (
-        sqrt(epsilonPlus*tauEta)
-      + dimensionedScalar("1e-8", up.dimensions(), 1e-8)
-    );
-
-    else if (XiModel == "algebraic")
-    {
-        // Simple algebraic model for Xi based on Gulders correlation
-        // with a linear correction function to give a plausible profile for Xi
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-        volScalarField GEta = GEtaCoef/tauEta;
-        volScalarField XiEqEta = 1.0 + XiCoef*sqrt(up/(Su + SuMin))*Reta;
-
-        volScalarField R =
-            GEta*XiEqEta/(XiEqEta - 0.999) + GIn*XiIn/(XiIn - 0.999);
-
-        volScalarField XiEqStar = R/(R - GEta - GIn);
-
-        //- Calculate the unweighted b
-        //volScalarField Rrho = thermo->rhou()/thermo->rhob();
-        //volScalarField bbar = b/(b + (Rrho*(1.0 - b)));
-
-        Xi == 1.0 + (1.0 + (2*XiShapeCoef)*(0.5 - b))*(XiEqStar - 1.0);
-    }
-    else if (XiModel == "transport")
-    {
-        // Calculate Xi transport coefficients based on Gulders correlation
-        // and DNS data for the rate of generation
-        // with a linear correction function to give a plausible profile for Xi
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-        volScalarField GEta = GEtaCoef/tauEta;
-        volScalarField XiEqEta = 1.0 + XiCoef*sqrt(up/(Su + SuMin))*Reta;
-
-        volScalarField R =
-            GEta*XiEqEta/(XiEqEta - 0.999) + GIn*XiIn/(XiIn - 0.999);
-
-        volScalarField XiEqStar = R/(R - GEta - GIn);
-
-        volScalarField XiEq =
-            1.0 + (1.0 + (2*XiShapeCoef)*(0.5 - b))*(XiEqStar - 1.0);
-
-        volScalarField G = R*(XiEq - 1.0)/XiEq;
-
-    // Calculate Xi flux
-    // ~~~~~~~~~~~~~~~~~
-    surfaceScalarField phiXi =
-        phiSt
-      + (
-          - fvc::interpolate(fvc::laplacian(Dbf, b)/mgb)*nf
-          + fvc::interpolate(rho)*fvc::interpolate(Su*(1.0/Xi - Xi))*nf
-        );
-
-
-        // Solve for the flame wrinkling
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-        solve
-        (
-            betav*fvm::ddt(rho, Xi)
-          + mvConvection->fvmDiv(phi, Xi)
-          + fvm::div(phiXi, Xi, "div(phiXi,Xi)")
-          - fvm::Sp(fvc::div(phiXi), Xi)
-         ==
-            betav*rho*R
-          - fvm::Sp(betav*rho*(R - G), Xi)
-        );
-
-
-        // Correct boundedness of Xi
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~
-        Xi.max(1.0);
-        Xi = min(Xi, 2.0*XiEq);
-        Info<< "max(Xi) = " << max(Xi).value() << endl;
-        Info<< "max(XiEq) = " << max(XiEq).value() << endl;
-    }
-    else
-    {
-        FatalError
-            << args.executable() << " : Unknown Xi model " << XiModel
-            << abort(FatalError);
-    }

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/Gulder/Gulder.C

@@ -34,8 +34,8 @@ namespace XiEqModels
 {
     defineTypeNameAndDebug(Gulder, 0);
     addToRunTimeSelectionTable(XiEqModel, Gulder, dictionary);
-};
+}
-};
+}
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -64,15 +64,18 @@ Foam::XiEqModels::Gulder::~Gulder()
 
 Foam::tmp<Foam::volScalarField> Foam::XiEqModels::Gulder::XiEq() const
 {
-    volScalarField up = sqrt((2.0/3.0)*turbulence_.k());
+    volScalarField up(sqrt((2.0/3.0)*turbulence_.k()));
     const volScalarField& epsilon = turbulence_.epsilon();
 
-    volScalarField tauEta = sqrt(mag(thermo_.muu()/(thermo_.rhou()*epsilon)));
+    volScalarField tauEta(sqrt(mag(thermo_.muu()/(thermo_.rhou()*epsilon))));
 
-    volScalarField Reta = up/
+    volScalarField Reta
     (
+        up
+      / (
             sqrt(epsilon*tauEta)
           + dimensionedScalar("1e-8", up.dimensions(), 1e-8)
+        )
     );
 
     return 1.0 + XiEqCoef*sqrt(up/(Su_ + SuMin))*Reta;

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/SCOPEXiEq/SCOPEXiEq.C

@@ -34,8 +34,8 @@ namespace XiEqModels
 {
     defineTypeNameAndDebug(SCOPEXiEq, 0);
     addToRunTimeSelectionTable(XiEqModel, SCOPEXiEq, dictionary);
-};
+}
-};
+}
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -83,13 +83,13 @@ Foam::tmp<Foam::volScalarField> Foam::XiEqModels::SCOPEXiEq::XiEq() const
     const volScalarField& k = turbulence_.k();
     const volScalarField& epsilon = turbulence_.epsilon();
 
-    volScalarField up = sqrt((2.0/3.0)*k);
+    volScalarField up(sqrt((2.0/3.0)*k));
-    volScalarField l = (lCoef*sqrt(3.0/2.0))*up*k/epsilon;
+    volScalarField l((lCoef*sqrt(3.0/2.0))*up*k/epsilon);
-    volScalarField Rl = up*l*thermo_.rhou()/thermo_.muu();
+    volScalarField Rl(up*l*thermo_.rhou()/thermo_.muu());
 
-    volScalarField upBySu = up/(Su_ + SuMin);
+    volScalarField upBySu(up/(Su_ + SuMin));
-    volScalarField K = 0.157*upBySu/sqrt(Rl);
+    volScalarField K(0.157*upBySu/sqrt(Rl));
-    volScalarField Ma = MaModel.Ma();
+    volScalarField Ma(MaModel.Ma());
 
     tmp<volScalarField> tXiEq
     (

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/instabilityXiEq/instabilityXiEq.C

@@ -34,8 +34,8 @@ namespace XiEqModels
 {
     defineTypeNameAndDebug(instability, 0);
     addToRunTimeSelectionTable(XiEqModel, instability, dictionary);
-};
+}
-};
+}
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -64,7 +64,7 @@ Foam::XiEqModels::instability::~instability()
 
 Foam::tmp<Foam::volScalarField> Foam::XiEqModels::instability::XiEq() const
 {
-    volScalarField turbXiEq = XiEqModel_->XiEq();
+    volScalarField turbXiEq(XiEqModel_->XiEq());
     return XiEqIn/turbXiEq + turbXiEq;
 }
 

applications/solvers/combustion/PDRFoam/XiModels/XiGModels/KTS/KTS.C

@@ -34,8 +34,8 @@ namespace XiGModels
 {
     defineTypeNameAndDebug(KTS, 0);
     addToRunTimeSelectionTable(XiGModel, KTS, dictionary);
-};
+}
-};
+}
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -63,10 +63,11 @@ Foam::XiGModels::KTS::~KTS()
 
 Foam::tmp<Foam::volScalarField> Foam::XiGModels::KTS::G() const
 {
-    volScalarField up = sqrt((2.0/3.0)*turbulence_.k());
+    // volScalarField up(sqrt((2.0/3.0)*turbulence_.k()));
     const volScalarField& epsilon = turbulence_.epsilon();
 
-    volScalarField tauEta = sqrt(mag(thermo_.muu()/(thermo_.rhou()*epsilon)));
+    tmp<volScalarField> tauEta =
+        sqrt(mag(thermo_.muu()/(thermo_.rhou()*epsilon)));
 
     return GEtaCoef/tauEta;
 }

applications/solvers/combustion/PDRFoam/XiModels/XiGModels/instabilityG/instabilityG.C

@@ -65,7 +65,7 @@ Foam::XiGModels::instabilityG::~instabilityG()
 
 Foam::tmp<Foam::volScalarField> Foam::XiGModels::instabilityG::G() const
 {
-    volScalarField turbXiG = XiGModel_->G();
+    volScalarField turbXiG(XiGModel_->G());
     return GIn*GIn/(GIn + turbXiG) + turbXiG;
 }
 

applications/solvers/combustion/PDRFoam/XiModels/algebraic/algebraic.C

@@ -34,8 +34,8 @@ namespace XiModels
 {
     defineTypeNameAndDebug(algebraic, 0);
     addToRunTimeSelectionTable(XiModel, algebraic, dictionary);
-};
+}
-};
+}
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -74,12 +74,12 @@ Foam::tmp<Foam::volScalarField> Foam::XiModels::algebraic::Db() const
 
 void Foam::XiModels::algebraic::correct()
 {
-    volScalarField XiEqEta = XiEqModel_->XiEq();
+    volScalarField XiEqEta(XiEqModel_->XiEq());
-    volScalarField GEta = XiGModel_->G();
+    volScalarField GEta(XiGModel_->G());
 
-    volScalarField R = GEta*XiEqEta/(XiEqEta - 0.999);
+    volScalarField R(GEta*XiEqEta/(XiEqEta - 0.999));
 
-    volScalarField XiEqStar = R/(R - GEta);
+    volScalarField XiEqStar(R/(R - GEta));
 
     Xi_ == 1.0 + (1.0 + (2*XiShapeCoef)*(0.5 - b_))*(XiEqStar - 1.0);
 }

applications/solvers/combustion/PDRFoam/XiModels/transport/transport.C

@@ -34,8 +34,8 @@ namespace XiModels
 {
     defineTypeNameAndDebug(transport, 0);
     addToRunTimeSelectionTable(XiModel, transport, dictionary);
-};
+}
-};
+}
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -77,17 +77,19 @@ void Foam::XiModels::transport::correct
     const fv::convectionScheme<scalar>& mvConvection
 )
 {
-    volScalarField XiEqEta = XiEqModel_->XiEq();
+    volScalarField XiEqEta(XiEqModel_->XiEq());
-    volScalarField GEta = XiGModel_->G();
+    volScalarField GEta(XiGModel_->G());
 
-    volScalarField R = GEta*XiEqEta/(XiEqEta - 0.999);
+    volScalarField R(GEta*XiEqEta/(XiEqEta - 0.999));
 
-    volScalarField XiEqStar = R/(R - GEta);
+    volScalarField XiEqStar(R/(R - GEta));
 
-    volScalarField XiEq =
+    volScalarField XiEq
-        1.0 + (1.0 + (2*XiShapeCoef)*(0.5 - b_))*(XiEqStar - 1.0);
+    (
+        1.0 + (1.0 + (2*XiShapeCoef)*(0.5 - b_))*(XiEqStar - 1.0)
+    );
 
-    volScalarField G = R*(XiEq - 1.0)/XiEq;
+    volScalarField G(R*(XiEq - 1.0)/XiEq);
 
     const objectRegistry& db = b_.db();
     const volScalarField& betav = db.lookupObject<volScalarField>("betav");

applications/solvers/combustion/PDRFoam/bEqn.H

@@ -25,20 +25,20 @@ if (ign.ignited())
 
     // Unburnt gas density
     // ~~~~~~~~~~~~~~~~~~~
-    volScalarField rhou = thermo.rhou();
+    volScalarField rhou(thermo.rhou());
 
     // Calculate flame normal etc.
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    //volVectorField n = fvc::grad(b);
+    // volVectorField n(fvc::grad(b));
-    volVectorField n = fvc::reconstruct(fvc::snGrad(b)*mesh.magSf());
+    volVectorField n(fvc::reconstruct(fvc::snGrad(b)*mesh.magSf()));
 
     volScalarField mgb("mgb", mag(n));
 
     dimensionedScalar dMgb("dMgb", mgb.dimensions(), SMALL);
 
     {
-        volScalarField bc = b*c;
+        volScalarField bc(b*c);
 
         dMgb += 1.0e-3*
             (bc*mgb)().weightedAverage(mesh.V())
@@ -47,8 +47,8 @@ if (ign.ignited())
 
     mgb += dMgb;
 
-    surfaceVectorField Sfhat = mesh.Sf()/mesh.magSf();
+    surfaceVectorField Sfhat(mesh.Sf()/mesh.magSf());
-    surfaceVectorField nfVec = fvc::interpolate(n);
+    surfaceVectorField nfVec(fvc::interpolate(n));
     nfVec += Sfhat*(fvc::snGrad(b) - (Sfhat & nfVec));
     nfVec /= (mag(nfVec) + dMgb);
     surfaceScalarField nf("nf", mesh.Sf() & nfVec);

applications/solvers/combustion/PDRFoam/pEqn.H

@@ -1,6 +1,6 @@
 rho = thermo.rho();
 
-volScalarField rAU = 1.0/UEqn.A();
+volScalarField rAU(1.0/UEqn.A());
 U = invA & UEqn.H();
 
 if (transonic)

applications/solvers/combustion/XiFoam/bEqn.H

@@ -2,18 +2,18 @@ if (ign.ignited())
 {
     // progress variable
     // ~~~~~~~~~~~~~~~~~
-    volScalarField c = scalar(1) - b;
+    volScalarField c(scalar(1) - b);
 
     // Unburnt gas density
     // ~~~~~~~~~~~~~~~~~~~
-    volScalarField rhou = thermo.rhou();
+    volScalarField rhou(thermo.rhou());
 
     // Calculate flame normal etc.
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    volVectorField n = fvc::grad(b);
+    volVectorField n(fvc::grad(b));
 
-    volScalarField mgb = mag(n);
+    volScalarField mgb(mag(n));
 
     dimensionedScalar dMgb = 1.0e-3*
         (b*c*mgb)().weightedAverage(mesh.V())
@@ -22,11 +22,11 @@ if (ign.ignited())
 
     mgb += dMgb;
 
-    surfaceVectorField SfHat = mesh.Sf()/mesh.magSf();
+    surfaceVectorField SfHat(mesh.Sf()/mesh.magSf());
-    surfaceVectorField nfVec = fvc::interpolate(n);
+    surfaceVectorField nfVec(fvc::interpolate(n));
     nfVec += SfHat*(fvc::snGrad(b) - (SfHat & nfVec));
     nfVec /= (mag(nfVec) + dMgb);
-    surfaceScalarField nf = (mesh.Sf() & nfVec);
+    surfaceScalarField nf((mesh.Sf() & nfVec));
     n /= mgb;
 
 
@@ -34,7 +34,7 @@ if (ign.ignited())
 
     // Calculate turbulent flame speed flux
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    surfaceScalarField phiSt = fvc::interpolate(rhou*StCorr*Su*Xi)*nf;
+    surfaceScalarField phiSt(fvc::interpolate(rhou*StCorr*Su*Xi)*nf);
 
     scalar StCoNum = max
     (
@@ -71,16 +71,20 @@ if (ign.ignited())
     // Calculate coefficients for Gulder's flame speed correlation
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    volScalarField up = uPrimeCoef*sqrt((2.0/3.0)*turbulence->k());
+    volScalarField up(uPrimeCoef*sqrt((2.0/3.0)*turbulence->k()));
-  //volScalarField up = sqrt(mag(diag(n * n) & diag(turbulence->r())));
+  //volScalarField up(sqrt(mag(diag(n * n) & diag(turbulence->r()))));
 
-    volScalarField epsilon = pow(uPrimeCoef, 3)*turbulence->epsilon();
+    volScalarField epsilon(pow(uPrimeCoef, 3)*turbulence->epsilon());
 
-    volScalarField tauEta = sqrt(thermo.muu()/(rhou*epsilon));
+    volScalarField tauEta(sqrt(thermo.muu()/(rhou*epsilon)));
 
-    volScalarField Reta = up/
+    volScalarField Reta
     (
-        sqrt(epsilon*tauEta) + dimensionedScalar("1e-8", up.dimensions(), 1e-8)
+        up
+      / (
+            sqrt(epsilon*tauEta)
+          + dimensionedScalar("1e-8", up.dimensions(), 1e-8)
+        )
     );
 
   //volScalarField l = 0.337*k*sqrt(k)/epsilon;
@@ -88,34 +92,38 @@ if (ign.ignited())
 
     // Calculate Xi flux
     // ~~~~~~~~~~~~~~~~~
-    surfaceScalarField phiXi =
+    surfaceScalarField phiXi
+    (
         phiSt
       - fvc::interpolate(fvc::laplacian(turbulence->alphaEff(), b)/mgb)*nf
-      + fvc::interpolate(rho)*fvc::interpolate(Su*(1.0/Xi - Xi))*nf;
+      + fvc::interpolate(rho)*fvc::interpolate(Su*(1.0/Xi - Xi))*nf
+    );
 
 
     // Calculate mean and turbulent strain rates
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    volVectorField Ut = U + Su*Xi*n;
+    volVectorField Ut(U + Su*Xi*n);
-    volScalarField sigmat = (n & n)*fvc::div(Ut) - (n & fvc::grad(Ut) & n);
+    volScalarField sigmat((n & n)*fvc::div(Ut) - (n & fvc::grad(Ut) & n));
 
-    volScalarField sigmas =
+    volScalarField sigmas
+    (
         ((n & n)*fvc::div(U) - (n & fvc::grad(U) & n))/Xi
       + (
             (n & n)*fvc::div(Su*n)
           - (n & fvc::grad(Su*n) & n)
-        )*(Xi + scalar(1))/(2*Xi);
+        )*(Xi + scalar(1))/(2*Xi)
+    );
 
 
     // Calculate the unstrained laminar flame speed
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    volScalarField Su0 = unstrainedLaminarFlameSpeed()();
+    volScalarField Su0(unstrainedLaminarFlameSpeed()());
 
 
     // Calculate the laminar flame speed in equilibrium with the applied strain
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    volScalarField SuInf = Su0*max(scalar(1) - sigmas/sigmaExt, scalar(0.01));
+    volScalarField SuInf(Su0*max(scalar(1) - sigmas/sigmaExt, scalar(0.01)));
 
     if (SuModel == "unstrained")
     {
@@ -130,9 +138,11 @@ if (ign.ignited())
         // Solve for the strained laminar flame speed
         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-        volScalarField Rc =
+        volScalarField Rc
+        (
             (sigmas*SuInf*(Su0 - SuInf) + sqr(SuMin)*sigmaExt)
-           /(sqr(Su0 - SuInf) + sqr(SuMin));
+            /(sqr(Su0 - SuInf) + sqr(SuMin))
+        );
 
         fvScalarMatrix SuEqn
         (
@@ -183,17 +193,21 @@ if (ign.ignited())
         // with a linear correction function to give a plausible profile for Xi
         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-        volScalarField XiEqStar =
+        volScalarField XiEqStar
-            scalar(1.001) + XiCoef*sqrt(up/(Su + SuMin))*Reta;
+        (
+            scalar(1.001) + XiCoef*sqrt(up/(Su + SuMin))*Reta
+        );
 
-        volScalarField XiEq =
+        volScalarField XiEq
+        (
             scalar(1.001)
           + (scalar(1) + (2*XiShapeCoef)*(scalar(0.5) - b))
-           *(XiEqStar - scalar(1.001));
+            *(XiEqStar - scalar(1.001))
+        );
 
-        volScalarField Gstar = 0.28/tauEta;
+        volScalarField Gstar(0.28/tauEta);
-        volScalarField R = Gstar*XiEqStar/(XiEqStar - scalar(1));
+        volScalarField R(Gstar*XiEqStar/(XiEqStar - scalar(1)));
-        volScalarField G = R*(XiEq - scalar(1.001))/XiEq;
+        volScalarField G(R*(XiEq - scalar(1.001))/XiEq);
 
         //R *= (Gstar + 2*mag(dev(symm(fvc::grad(U)))))/Gstar;
 

applications/solvers/combustion/XiFoam/createFields.H

@@ -61,8 +61,10 @@
     );
 
     Info<< "Creating field DpDt\n" << endl;
-    volScalarField DpDt =
+    volScalarField DpDt
-        fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+    (
+        fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p)
+    );
 
 
     Info<< "Creating field Xi\n" << endl;

applications/solvers/combustion/XiFoam/pEqn.H

@@ -1,6 +1,6 @@
 rho = thermo.rho();
 
-volScalarField rAU = 1.0/UEqn.A();
+volScalarField rAU(1.0/UEqn.A());
 U = rAU*UEqn.H();
 
 if (transonic)

applications/solvers/combustion/coldEngineFoam/createFields.H

@@ -55,5 +55,7 @@
     );
 
     Info<< "Creating field DpDt\n" << endl;
-    volScalarField DpDt =
+    volScalarField DpDt
-        fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+    (
+        fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p)
+    );

applications/solvers/combustion/dieselEngineFoam/YEqn.H

@@ -12,7 +12,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
 {
 
     label inertIndex = -1;
-    volScalarField Yt = 0.0*Y[0];
+    volScalarField Yt(0.0*Y[0]);
 
     forAll(Y, i)
     {

applications/solvers/combustion/dieselEngineFoam/createFields.H

@@ -82,8 +82,10 @@ autoPtr<compressible::turbulenceModel> turbulence
 );
 
 Info<< "Creating field DpDt\n" << endl;
-volScalarField DpDt =
+volScalarField DpDt
-    fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+(
+    fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p)
+);
 
 
 multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;

applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.C

@@ -94,11 +94,13 @@ int main(int argc, char *argv[])
 
         // turbulent time scale
         {
-            volScalarField tk =
+            volScalarField tk
-                Cmix*sqrt(turbulence->muEff()/rho/turbulence->epsilon());
+            (
-            volScalarField tc = chemistry.tc();
+                Cmix*sqrt(turbulence->muEff()/rho/turbulence->epsilon())
+            );
+            volScalarField tc(chemistry.tc());
 
-            //Chalmers PaSR model
+            // Chalmers PaSR model
             kappa = (runTime.deltaT() + tc)/(runTime.deltaT() + tc + tk);
         }
 

applications/solvers/combustion/dieselEngineFoam/pEqn.H

@@ -1,6 +1,6 @@
 rho = thermo.rho();
 
-volScalarField A = UEqn.A();
+volScalarField A(UEqn.A());
 U = UEqn.H()/A;
 
 if (transonic)

applications/solvers/combustion/dieselFoam/dieselFoam.C

@@ -85,9 +85,11 @@ int main(int argc, char *argv[])
 
         // turbulent time scale
         {
-            volScalarField tk =
+            volScalarField tk
-                Cmix*sqrt(turbulence->muEff()/rho/turbulence->epsilon());
+            (
-            volScalarField tc = chemistry.tc();
+                Cmix*sqrt(turbulence->muEff()/rho/turbulence->epsilon())
+            );
+            volScalarField tc(chemistry.tc());
 
             // Chalmers PaSR model
             kappa = (runTime.deltaT() + tc)/(runTime.deltaT()+tc+tk);

applications/solvers/combustion/dieselFoam/pEqn.H

@@ -1,6 +1,6 @@
 rho = thermo.rho();
 
-volScalarField rAU = 1.0/UEqn.A();
+volScalarField rAU(1.0/UEqn.A());
 U = rAU*UEqn.H();
 
 if (transonic)

applications/solvers/combustion/engineFoam/pEqn.H

@@ -1,6 +1,6 @@
 rho = thermo.rho();
 
-volScalarField rAU = 1.0/UEqn.A();
+volScalarField rAU(1.0/UEqn.A());
 U = rAU*UEqn.H();
 
 if (transonic)

applications/solvers/combustion/fireFoam/combustionModels/combustionModel/combustionModel.C

@@ -32,7 +32,7 @@ namespace Foam
 {
     defineTypeNameAndDebug(combustionModel, 0);
     defineRunTimeSelectionTable(combustionModel, dictionary);
-};
+}
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -78,8 +78,8 @@ Foam::combustionModel::combustionModel::R(volScalarField& fu) const
 {
     const basicMultiComponentMixture& composition = thermo_.composition();
     const volScalarField& ft = composition.Y("ft");
-    volScalarField fres = composition.fres(ft, stoicRatio_.value());
+    volScalarField fres(composition.fres(ft, stoicRatio_.value()));
-    volScalarField wFuelNorm = this->wFuelNorm()*pos(fu - fres);
+    volScalarField wFuelNorm(this->wFuelNorm()*pos(fu - fres));
 
     return wFuelNorm*fres - fvm::Sp(wFuelNorm, fu);
 }

applications/solvers/combustion/fireFoam/createFields.H

@@ -123,8 +123,10 @@ volScalarField dQ
 
 
 Info<< "Creating field DpDt\n" << endl;
-volScalarField DpDt =
+volScalarField DpDt
-    fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+(
+    fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p)
+);
 
 
 dimensionedScalar initialMass = fvc::domainIntegrate(rho);

applications/solvers/combustion/fireFoam/fuhsEqn.H

@@ -1,7 +1,7 @@
 {
     // Solve fuel equation
     // ~~~~~~~~~~~~~~~~~~~
-    fvScalarMatrix R = combustion->R(fu);
+    fvScalarMatrix R(combustion->R(fu));
 
     {
         fvScalarMatrix fuEqn

applications/solvers/combustion/fireFoam/pEqn.H

@@ -1,6 +1,6 @@
 rho = thermo.rho();
 
-volScalarField rAU = 1.0/UEqn.A();
+volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("(rho*(1|A(U)))", fvc::interpolate(rho*rAU));
 U = rAU*UEqn.H();
 
@@ -17,7 +17,7 @@ phi = phiU - rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf();
 
 for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 {
-    surfaceScalarField rhorAUf = fvc::interpolate(rho*rAU);
+    surfaceScalarField rhorAUf(fvc::interpolate(rho*rAU));
 
     fvScalarMatrix p_rghEqn
     (

applications/solvers/combustion/reactingFoam/YEqn.H

@@ -11,7 +11,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
 
 {
     label inertIndex = -1;
-    volScalarField Yt = 0.0*Y[0];
+    volScalarField Yt(0.0*Y[0]);
 
     forAll(Y, i)
     {

applications/solvers/combustion/reactingFoam/chemistry.H

@@ -10,9 +10,14 @@
     // turbulent time scale
     if (turbulentReaction)
     {
-        volScalarField tk =
+        volScalarField tk
-            Cmix*sqrt(turbulence->muEff()/rho/turbulence->epsilon());
+        (
-        volScalarField tc = chemistry.tc();
+            Cmix*sqrt(turbulence->muEff()/rho/turbulence->epsilon())
+        );
+        volScalarField tc
+        (
+            chemistry.tc()
+        );
 
         // Chalmers PaSR model
         kappa = (runTime.deltaT() + tc)/(runTime.deltaT() + tc + tk);

applications/solvers/combustion/reactingFoam/createFields.H

@@ -72,8 +72,10 @@ autoPtr<compressible::turbulenceModel> turbulence
 );
 
 Info<< "Creating field DpDt\n" << endl;
-volScalarField DpDt =
+volScalarField DpDt
-    fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+(
+    fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p)
+);
 
 multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 

applications/solvers/combustion/reactingFoam/pEqn.H

@@ -1,6 +1,6 @@
 rho = thermo.rho();
 
-volScalarField rAU = 1.0/UEqn.A();
+volScalarField rAU(1.0/UEqn.A());
 U = rAU*UEqn.H();
 
 if (transonic)

applications/solvers/combustion/rhoReactingFoam/YEqn.H

@@ -11,7 +11,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
 
 {
     label inertIndex = -1;
-    volScalarField Yt = 0.0*Y[0];
+    volScalarField Yt(0.0*Y[0]);
 
     forAll(Y, i)
     {

applications/solvers/combustion/rhoReactingFoam/chemistry.H

@@ -10,9 +10,11 @@
     // turbulent time scale
     if (turbulentReaction)
     {
-        volScalarField tk =
+        volScalarField tk
-                Cmix*sqrt(turbulence->muEff()/rho/turbulence->epsilon());
+        (
-        volScalarField tc = chemistry.tc();
+            Cmix*sqrt(turbulence->muEff()/rho/turbulence->epsilon())
+        );
+        volScalarField tc(chemistry.tc());
 
         // Chalmers PaSR model
         kappa = (runTime.deltaT() + tc)/(runTime.deltaT() + tc + tk);

applications/solvers/combustion/rhoReactingFoam/createFields.H

@@ -73,8 +73,10 @@ autoPtr<compressible::turbulenceModel> turbulence
 );
 
 Info<< "Creating field DpDt\n" << endl;
-volScalarField DpDt =
+volScalarField DpDt
-    fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+(
+    fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p)
+);
 
 multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 

applications/solvers/combustion/rhoReactingFoam/pEqn.H

@@ -5,14 +5,16 @@
     // pressure solution - done in 2 parts. Part 1:
     thermo.rho() -= psi*p;
 
-    volScalarField rAU = 1.0/UEqn.A();
+    volScalarField rAU(1.0/UEqn.A());
     U = rAU*UEqn.H();
 
     if (transonic)
     {
-        surfaceScalarField phiv =
+        surfaceScalarField phiv
+        (
             (fvc::interpolate(U) & mesh.Sf())
-          + fvc::ddtPhiCorr(rAU, rho, U, phi);
+          + fvc::ddtPhiCorr(rAU, rho, U, phi)
+        );
 
         phi = fvc::interpolate(rho)*phiv;
 

applications/solvers/compressible/rhoCentralFoam/BCs/T/smoluchowskiJumpTFvPatchScalarField.C

@@ -181,13 +181,16 @@ void Foam::smoluchowskiJumpTFvPatchScalarField::updateCoeffs()
         )
     );
 
-    Field<scalar> C2 = pmu/prho
+    Field<scalar> C2
+    (
+        pmu/prho
         *sqrt(ppsi*constant::mathematical::piByTwo)
         *2.0*gamma_/Pr.value()/(gamma_ + 1.0)
-        *(2.0 - accommodationCoeff_)/accommodationCoeff_;
+        *(2.0 - accommodationCoeff_)/accommodationCoeff_
+    );
 
-    Field<scalar> aCoeff = prho.snGrad() - prho/C2;
+    Field<scalar> aCoeff(prho.snGrad() - prho/C2);
-    Field<scalar> KEbyRho = 0.5*magSqr(pU);
+    Field<scalar> KEbyRho(0.5*magSqr(pU));
 
     valueFraction() = (1.0/(1.0 + patch().deltaCoeffs()*C2));
     refValue() = Twall_;
@@ -214,9 +217,11 @@ void Foam::smoluchowskiJumpTFvPatchScalarField::write(Ostream& os) const
 
 namespace Foam
 {
-
+    makeNonTemplatedPatchTypeField
-makePatchTypeField(fvPatchScalarField, smoluchowskiJumpTFvPatchScalarField);
+    (
-
+        fvPatchScalarField,
+        smoluchowskiJumpTFvPatchScalarField
+    );
 }
 
 

applications/solvers/compressible/rhoCentralFoam/BCs/U/maxwellSlipUFvPatchVectorField.C

@@ -83,8 +83,7 @@ maxwellSlipUFvPatchVectorField::maxwellSlipUFvPatchVectorField
     if
     (
         mag(accommodationCoeff_) < SMALL
-        ||
+     || mag(accommodationCoeff_) > 2.0
-        mag(accommodationCoeff_) > 2.0
     )
     {
         FatalIOErrorIn
@@ -142,10 +141,13 @@ void maxwellSlipUFvPatchVectorField::updateCoeffs()
     const fvPatchField<scalar>& ppsi =
         patch().lookupPatchField<volScalarField, scalar>("psi");
 
-    Field<scalar> C1 = sqrt(ppsi*constant::mathematical::piByTwo)
+    Field<scalar> C1
-        *(2.0 - accommodationCoeff_)/accommodationCoeff_;
+    (
+        sqrt(ppsi*constant::mathematical::piByTwo)
+      * (2.0 - accommodationCoeff_)/accommodationCoeff_
+    );
 
-    Field<scalar> pnu = pmu/prho;
+    Field<scalar> pnu(pmu/prho);
     valueFraction() = (1.0/(1.0 + patch().deltaCoeffs()*C1*pnu));
 
     refValue() = Uwall_;
@@ -156,8 +158,8 @@ void maxwellSlipUFvPatchVectorField::updateCoeffs()
             this->db().objectRegistry::lookupObject<volScalarField>("T");
         label patchi = this->patch().index();
         const fvPatchScalarField& pT = vsfT.boundaryField()[patchi];
-        Field<vector> gradpT = fvc::grad(vsfT)().boundaryField()[patchi];
+        Field<vector> gradpT(fvc::grad(vsfT)().boundaryField()[patchi]);
-        vectorField n = patch().nf();
+        vectorField n(patch().nf());
 
         refValue() -= 3.0*pnu/(4.0*pT)*transform(I - n*n, gradpT);
     }
@@ -166,7 +168,7 @@ void maxwellSlipUFvPatchVectorField::updateCoeffs()
     {
         const fvPatchTensorField& ptauMC =
             patch().lookupPatchField<volTensorField, tensor>("tauMC");
-        vectorField n = patch().nf();
+        vectorField n(patch().nf());
 
         refValue() -= C1/prho*transform(I - n*n, (n & ptauMC));
     }
@@ -196,7 +198,11 @@ void maxwellSlipUFvPatchVectorField::write(Ostream& os) const
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-makePatchTypeField(fvPatchVectorField, maxwellSlipUFvPatchVectorField);
+makeNonTemplatedPatchTypeField
+(
+    fvPatchVectorField,
+    maxwellSlipUFvPatchVectorField
+);
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/compressible/rhoCentralFoam/BCs/mixedFixedValueSlip/mixedFixedValueSlipFvPatchField.C

@@ -135,8 +135,8 @@ void mixedFixedValueSlipFvPatchField<Type>::rmap
 template<class Type>
 tmp<Field<Type> > mixedFixedValueSlipFvPatchField<Type>::snGrad() const
 {
-    vectorField nHat = this->patch().nf();
+    tmp<vectorField> nHat = this->patch().nf();
-    Field<Type> pif = this->patchInternalField();
+    Field<Type> pif(this->patchInternalField());
 
     return
     (
@@ -155,7 +155,7 @@ void mixedFixedValueSlipFvPatchField<Type>::evaluate(const Pstream::commsTypes)
         this->updateCoeffs();
     }
 
-    vectorField nHat = this->patch().nf();
+    vectorField nHat(this->patch().nf());
 
     Field<Type>::operator=
     (
@@ -174,7 +174,7 @@ template<class Type>
 tmp<Field<Type> >
 mixedFixedValueSlipFvPatchField<Type>::snGradTransformDiag() const
 {
-    vectorField nHat = this->patch().nf();
+    vectorField nHat(this->patch().nf());
     vectorField diag(nHat.size());
 
     diag.replace(vector::X, mag(nHat.component(vector::X)));

applications/solvers/compressible/rhoCentralFoam/BCs/mixedFixedValueSlip/mixedFixedValueSlipFvPatchFields.H

@@ -36,7 +36,7 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-makePatchTypeFieldTypedefs(mixedFixedValueSlip)
+makePatchTypeFieldTypedefs(mixedFixedValueSlip);
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/compressible/rhoCentralFoam/BCs/mixedFixedValueSlip/mixedFixedValueSlipFvPatchFieldsFwd.H

@@ -37,7 +37,7 @@ namespace Foam
 
 template<class Type> class mixedFixedValueSlipFvPatchField;
 
-makePatchTypeFieldTypedefs(mixedFixedValueSlip)
+makePatchTypeFieldTypedefs(mixedFixedValueSlip);
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/compressible/rhoCentralFoam/BCs/rho/fixedRhoFvPatchScalarField.C

@@ -112,7 +112,11 @@ void fixedRhoFvPatchScalarField::updateCoeffs()
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-makePatchTypeField(fvPatchScalarField, fixedRhoFvPatchScalarField);
+makeNonTemplatedPatchTypeField
+(
+    fvPatchScalarField,
+    fixedRhoFvPatchScalarField
+);
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/compressible/rhoCentralFoam/compressibleCourantNo.H

@@ -34,8 +34,10 @@ scalar meanCoNum = 0.0;
 
 if (mesh.nInternalFaces())
 {
-    surfaceScalarField amaxSfbyDelta =
+    surfaceScalarField amaxSfbyDelta
-        mesh.surfaceInterpolation::deltaCoeffs()*amaxSf;
+    (
+        mesh.surfaceInterpolation::deltaCoeffs()*amaxSf
+    );
 
     CoNum = max(amaxSfbyDelta/mesh.magSf()).value()*runTime.deltaTValue();
 

applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.C

@@ -62,52 +62,76 @@ int main(int argc, char *argv[])
     {
         // --- upwind interpolation of primitive fields on faces
 
-        surfaceScalarField rho_pos =
+        surfaceScalarField rho_pos
-            fvc::interpolate(rho, pos, "reconstruct(rho)");
+        (
-        surfaceScalarField rho_neg =
+            fvc::interpolate(rho, pos, "reconstruct(rho)")
-            fvc::interpolate(rho, neg, "reconstruct(rho)");
+        );
+        surfaceScalarField rho_neg
+        (
+            fvc::interpolate(rho, neg, "reconstruct(rho)")
+        );
 
-        surfaceVectorField rhoU_pos =
+        surfaceVectorField rhoU_pos
-            fvc::interpolate(rhoU, pos, "reconstruct(U)");
+        (
-        surfaceVectorField rhoU_neg =
+            fvc::interpolate(rhoU, pos, "reconstruct(U)")
-            fvc::interpolate(rhoU, neg, "reconstruct(U)");
+        );
+        surfaceVectorField rhoU_neg
+        (
+            fvc::interpolate(rhoU, neg, "reconstruct(U)")
+        );
 
-        volScalarField rPsi = 1.0/psi;
+        volScalarField rPsi(1.0/psi);
-        surfaceScalarField rPsi_pos =
+        surfaceScalarField rPsi_pos
-            fvc::interpolate(rPsi, pos, "reconstruct(T)");
+        (
-        surfaceScalarField rPsi_neg =
+            fvc::interpolate(rPsi, pos, "reconstruct(T)")
-            fvc::interpolate(rPsi, neg, "reconstruct(T)");
+        );
+        surfaceScalarField rPsi_neg
+        (
+            fvc::interpolate(rPsi, neg, "reconstruct(T)")
+        );
 
-        surfaceScalarField e_pos =
+        surfaceScalarField e_pos
-            fvc::interpolate(e, pos, "reconstruct(T)");
+        (
-        surfaceScalarField e_neg =
+            fvc::interpolate(e, pos, "reconstruct(T)")
-            fvc::interpolate(e, neg, "reconstruct(T)");
+        );
+        surfaceScalarField e_neg
+        (
+            fvc::interpolate(e, neg, "reconstruct(T)")
+        );
 
-        surfaceVectorField U_pos = rhoU_pos/rho_pos;
+        surfaceVectorField U_pos(rhoU_pos/rho_pos);
-        surfaceVectorField U_neg = rhoU_neg/rho_neg;
+        surfaceVectorField U_neg(rhoU_neg/rho_neg);
 
-        surfaceScalarField p_pos = rho_pos*rPsi_pos;
+        surfaceScalarField p_pos(rho_pos*rPsi_pos);
-        surfaceScalarField p_neg = rho_neg*rPsi_neg;
+        surfaceScalarField p_neg(rho_neg*rPsi_neg);
 
-        surfaceScalarField phiv_pos = U_pos & mesh.Sf();
+        surfaceScalarField phiv_pos(U_pos & mesh.Sf());
-        surfaceScalarField phiv_neg = U_neg & mesh.Sf();
+        surfaceScalarField phiv_neg(U_neg & mesh.Sf());
 
-        volScalarField c = sqrt(thermo.Cp()/thermo.Cv()*rPsi);
+        volScalarField c(sqrt(thermo.Cp()/thermo.Cv()*rPsi));
-        surfaceScalarField cSf_pos =
+        surfaceScalarField cSf_pos
-            fvc::interpolate(c, pos, "reconstruct(T)")*mesh.magSf();
+        (
-        surfaceScalarField cSf_neg =
+            fvc::interpolate(c, pos, "reconstruct(T)")*mesh.magSf()
-            fvc::interpolate(c, neg, "reconstruct(T)")*mesh.magSf();
+        );
+        surfaceScalarField cSf_neg
+        (
+            fvc::interpolate(c, neg, "reconstruct(T)")*mesh.magSf()
+        );
 
-        surfaceScalarField ap =
+        surfaceScalarField ap
-            max(max(phiv_pos + cSf_pos, phiv_neg + cSf_neg), v_zero);
+        (
-        surfaceScalarField am =
+            max(max(phiv_pos + cSf_pos, phiv_neg + cSf_neg), v_zero)
-            min(min(phiv_pos - cSf_pos, phiv_neg - cSf_neg), v_zero);
+        );
+        surfaceScalarField am
+        (
+            min(min(phiv_pos - cSf_pos, phiv_neg - cSf_neg), v_zero)
+        );
 
-        surfaceScalarField a_pos = ap/(ap - am);
+        surfaceScalarField a_pos(ap/(ap - am));
 
         surfaceScalarField amaxSf("amaxSf", max(mag(am), mag(ap)));
 
-        surfaceScalarField aSf = am*a_pos;
+        surfaceScalarField aSf(am*a_pos);
 
         if (fluxScheme == "Tadmor")
         {
@@ -115,13 +139,13 @@ int main(int argc, char *argv[])
             a_pos = 0.5;
         }
 
-        surfaceScalarField a_neg = (1.0 - a_pos);
+        surfaceScalarField a_neg(1.0 - a_pos);
 
         phiv_pos *= a_pos;
         phiv_neg *= a_neg;
 
-        surfaceScalarField aphiv_pos = phiv_pos - aSf;
+        surfaceScalarField aphiv_pos(phiv_pos - aSf);
-        surfaceScalarField aphiv_neg = phiv_neg + aSf;
+        surfaceScalarField aphiv_neg(phiv_neg + aSf);
 
         // Reuse amaxSf for the maximum positive and negative fluxes
         // estimated by the central scheme
@@ -148,14 +172,18 @@ int main(int argc, char *argv[])
 
         surfaceScalarField phi("phi", aphiv_pos*rho_pos + aphiv_neg*rho_neg);
 
-        surfaceVectorField phiUp =
+        surfaceVectorField phiUp
+        (
             (aphiv_pos*rhoU_pos + aphiv_neg*rhoU_neg)
-          + (a_pos*p_pos + a_neg*p_neg)*mesh.Sf();
+          + (a_pos*p_pos + a_neg*p_neg)*mesh.Sf()
+        );
 
-        surfaceScalarField phiEp =
+        surfaceScalarField phiEp
+        (
             aphiv_pos*(rho_pos*(e_pos + 0.5*magSqr(U_pos)) + p_pos)
           + aphiv_neg*(rho_neg*(e_neg + 0.5*magSqr(U_neg)) + p_neg)
-          + aSf*p_pos - aSf*p_neg;
+          + aSf*p_pos - aSf*p_neg
+        );
 
         volTensorField tauMC("tauMC", mu*dev2(Foam::T(fvc::grad(U))));
 
@@ -185,7 +213,7 @@ int main(int argc, char *argv[])
         }
 
         // --- Solve energy
-        surfaceScalarField sigmaDotU =
+        surfaceScalarField sigmaDotU
         (
             (
                 fvc::interpolate(mu)*mesh.magSf()*fvc::snGrad(U)

applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C

@@ -59,48 +59,72 @@ int main(int argc, char *argv[])
     {
         // --- upwind interpolation of primitive fields on faces
 
-        surfaceScalarField rho_pos =
+        surfaceScalarField rho_pos
-            fvc::interpolate(rho, pos, "reconstruct(rho)");
+        (
-        surfaceScalarField rho_neg =
+            fvc::interpolate(rho, pos, "reconstruct(rho)")
-            fvc::interpolate(rho, neg, "reconstruct(rho)");
+        );
+        surfaceScalarField rho_neg
+        (
+            fvc::interpolate(rho, neg, "reconstruct(rho)")
+        );
 
-        surfaceVectorField rhoU_pos =
+        surfaceVectorField rhoU_pos
-            fvc::interpolate(rhoU, pos, "reconstruct(U)");
+        (
-        surfaceVectorField rhoU_neg =
+            fvc::interpolate(rhoU, pos, "reconstruct(U)")
-            fvc::interpolate(rhoU, neg, "reconstruct(U)");
+        );
+        surfaceVectorField rhoU_neg
+        (
+            fvc::interpolate(rhoU, neg, "reconstruct(U)")
+        );
 
-        volScalarField rPsi = 1.0/psi;
+        volScalarField rPsi(1.0/psi);
-        surfaceScalarField rPsi_pos =
+        surfaceScalarField rPsi_pos
-            fvc::interpolate(rPsi, pos, "reconstruct(T)");
+        (
-        surfaceScalarField rPsi_neg =
+            fvc::interpolate(rPsi, pos, "reconstruct(T)")
-            fvc::interpolate(rPsi, neg, "reconstruct(T)");
+        );
+        surfaceScalarField rPsi_neg
+        (
+            fvc::interpolate(rPsi, neg, "reconstruct(T)")
+        );
 
-        surfaceScalarField e_pos =
+        surfaceScalarField e_pos
-            fvc::interpolate(e, pos, "reconstruct(T)");
+        (
-        surfaceScalarField e_neg =
+            fvc::interpolate(e, pos, "reconstruct(T)")
-            fvc::interpolate(e, neg, "reconstruct(T)");
+        );
+        surfaceScalarField e_neg
+        (
+            fvc::interpolate(e, neg, "reconstruct(T)")
+        );
 
-        surfaceVectorField U_pos = rhoU_pos/rho_pos;
+        surfaceVectorField U_pos(rhoU_pos/rho_pos);
-        surfaceVectorField U_neg = rhoU_neg/rho_neg;
+        surfaceVectorField U_neg(rhoU_neg/rho_neg);
 
-        surfaceScalarField p_pos = rho_pos*rPsi_pos;
+        surfaceScalarField p_pos(rho_pos*rPsi_pos);
-        surfaceScalarField p_neg = rho_neg*rPsi_neg;
+        surfaceScalarField p_neg(rho_neg*rPsi_neg);
 
-        surfaceScalarField phiv_pos = U_pos & mesh.Sf();
+        surfaceScalarField phiv_pos(U_pos & mesh.Sf());
-        surfaceScalarField phiv_neg = U_neg & mesh.Sf();
+        surfaceScalarField phiv_neg(U_neg & mesh.Sf());
 
-        volScalarField c = sqrt(thermo.Cp()/thermo.Cv()*rPsi);
+        volScalarField c(sqrt(thermo.Cp()/thermo.Cv()*rPsi));
-        surfaceScalarField cSf_pos =
+        surfaceScalarField cSf_pos
-            fvc::interpolate(c, pos, "reconstruct(T)")*mesh.magSf();
+        (
-        surfaceScalarField cSf_neg =
+            fvc::interpolate(c, pos, "reconstruct(T)")*mesh.magSf()
-            fvc::interpolate(c, neg, "reconstruct(T)")*mesh.magSf();
+        );
+        surfaceScalarField cSf_neg
+        (
+            fvc::interpolate(c, neg, "reconstruct(T)")*mesh.magSf()
+        );
 
-        surfaceScalarField ap =
+        surfaceScalarField ap
-            max(max(phiv_pos + cSf_pos, phiv_neg + cSf_neg), v_zero);
+        (
-        surfaceScalarField am =
+            max(max(phiv_pos + cSf_pos, phiv_neg + cSf_neg), v_zero)
-            min(min(phiv_pos - cSf_pos, phiv_neg - cSf_neg), v_zero);
+        );
+        surfaceScalarField am
+        (
+            min(min(phiv_pos - cSf_pos, phiv_neg - cSf_neg), v_zero)
+        );
 
-        surfaceScalarField a_pos = ap/(ap - am);
+        surfaceScalarField a_pos(ap/(ap - am));
 
         surfaceScalarField amaxSf("amaxSf", max(mag(am), mag(ap)));
 
@@ -112,7 +136,7 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        surfaceScalarField aSf = am*a_pos;
+        surfaceScalarField aSf(am*a_pos);
 
         if (fluxScheme == "Tadmor")
         {
@@ -120,24 +144,28 @@ int main(int argc, char *argv[])
             a_pos = 0.5;
         }
 
-        surfaceScalarField a_neg = (1.0 - a_pos);
+        surfaceScalarField a_neg(1.0 - a_pos);
 
         phiv_pos *= a_pos;
         phiv_neg *= a_neg;
 
-        surfaceScalarField aphiv_pos = phiv_pos - aSf;
+        surfaceScalarField aphiv_pos(phiv_pos - aSf);
-        surfaceScalarField aphiv_neg = phiv_neg + aSf;
+        surfaceScalarField aphiv_neg(phiv_neg + aSf);
 
         surfaceScalarField phi("phi", aphiv_pos*rho_pos + aphiv_neg*rho_neg);
 
-        surfaceVectorField phiUp =
+        surfaceVectorField phiUp
+        (
             (aphiv_pos*rhoU_pos + aphiv_neg*rhoU_neg)
-          + (a_pos*p_pos + a_neg*p_neg)*mesh.Sf();
+          + (a_pos*p_pos + a_neg*p_neg)*mesh.Sf()
+        );
 
-        surfaceScalarField phiEp =
+        surfaceScalarField phiEp
+        (
             aphiv_pos*(rho_pos*(e_pos + 0.5*magSqr(U_pos)) + p_pos)
           + aphiv_neg*(rho_neg*(e_neg + 0.5*magSqr(U_neg)) + p_neg)
-          + aSf*p_pos - aSf*p_neg;
+          + aSf*p_pos - aSf*p_neg
+        );
 
         volTensorField tauMC("tauMC", mu*dev2(Foam::T(fvc::grad(U))));
 
@@ -167,7 +195,7 @@ int main(int argc, char *argv[])
         }
 
         // --- Solve energy
-        surfaceScalarField sigmaDotU =
+        surfaceScalarField sigmaDotU
         (
             (
                 fvc::interpolate(mu)*mesh.magSf()*fvc::snGrad(U)

applications/solvers/compressible/rhoPimpleFoam/UEqn.H

@@ -9,7 +9,7 @@ tmp<fvVectorMatrix> UEqn
 
 UEqn().relax();
 
-volScalarField rAU = 1.0/UEqn().A();
+volScalarField rAU(1.0/UEqn().A());
 
 if (momentumPredictor)
 {

applications/solvers/compressible/rhoPimpleFoam/createFields.H

@@ -52,5 +52,7 @@
     );
 
     Info<< "Creating field DpDt\n" << endl;
-    volScalarField DpDt =
+    volScalarField DpDt
-        fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+    (
+        fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p)
+    );

applications/solvers/compressible/rhoPorousMRFPimpleFoam/UEqn.H

@@ -12,7 +12,7 @@ UEqn().relax();
 mrfZones.addCoriolis(rho, UEqn());
 pZones.addResistance(UEqn());
 
-volScalarField rAU = 1.0/UEqn().A();
+volScalarField rAU(1.0/UEqn().A());
 
 if (momentumPredictor)
 {

applications/solvers/compressible/rhoPorousMRFPimpleFoam/createFields.H

@@ -65,8 +65,10 @@
 
 
     Info<< "Creating field DpDt\n" << endl;
-    volScalarField DpDt =
+    volScalarField DpDt
-        fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+    (
+        fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p)
+    );
 
     MRFZones mrfZones(mesh);
     mrfZones.correctBoundaryVelocity(U);

applications/solvers/compressible/rhoPorousMRFPimpleFoam/pEqn.H

@@ -1,6 +1,6 @@
 rho = thermo.rho();
 
-volScalarField rAU = 1.0/UEqn().A();
+volScalarField rAU(1.0/UEqn().A());
 U = rAU*UEqn().H();
 
 if (nCorr <= 1)

applications/solvers/compressible/rhoPorousSimpleFoam/UEqn.H

@@ -22,7 +22,7 @@
         trTU = inv(tTU());
         trTU().rename("rAU");
 
-        volVectorField gradp = fvc::grad(p);
+        volVectorField gradp(fvc::grad(p));
 
         for (int UCorr=0; UCorr<nUCorr; UCorr++)
         {

applications/solvers/compressible/rhoSimpleFoam/pEqn.H

@@ -3,7 +3,7 @@ rho = max(rho, rhoMin);
 rho = min(rho, rhoMax);
 rho.relax();
 
-volScalarField rAU = 1.0/UEqn().A();
+volScalarField rAU(1.0/UEqn().A());
 U = rAU*UEqn().H();
 UEqn.clear();
 

applications/solvers/compressible/rhoSimplecFoam/pEqn.H

@@ -3,10 +3,10 @@ rho = max(rho, rhoMin);
 rho = min(rho, rhoMax);
 rho.relax();
 
-volScalarField p0 = p;
+volScalarField p0(p);
 
-volScalarField AU = UEqn().A();
+volScalarField AU(UEqn().A());
-volScalarField AtU = AU - UEqn().H1();
+volScalarField AtU(AU - UEqn().H1());
 U = UEqn().H()/AU;
 UEqn.clear();
 

applications/solvers/compressible/sonicDyMFoam/createFields.H

@@ -51,5 +51,7 @@
     );
 
     Info<< "Creating field DpDt\n" << endl;
-    volScalarField DpDt =
+    volScalarField DpDt
-        fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+    (
+        fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p)
+    );

applications/solvers/compressible/sonicFoam/pEqn.H

@@ -1,6 +1,6 @@
 rho = thermo.rho();
 
-volScalarField rAU = 1.0/UEqn.A();
+volScalarField rAU(1.0/UEqn.A());
 U = rAU*UEqn.H();
 
 surfaceScalarField phid

applications/solvers/compressible/sonicLiquidFoam/sonicLiquidFoam.C

@@ -71,7 +71,7 @@ int main(int argc, char *argv[])
 
         for (int corr=0; corr<nCorr; corr++)
         {
-            volScalarField rAU = 1.0/UEqn.A();
+            volScalarField rAU(1.0/UEqn.A());
             U = rAU*UEqn.H();
 
             surfaceScalarField phid

applications/solvers/electromagnetics/mhdFoam/mhdFoam.C

@@ -92,7 +92,7 @@ int main(int argc, char *argv[])
 
             for (int corr=0; corr<nCorr; corr++)
             {
-                volScalarField rAU = 1.0/UEqn.A();
+                volScalarField rAU(1.0/UEqn.A());
 
                 U = rAU*UEqn.H();
 
@@ -136,7 +136,7 @@ int main(int argc, char *argv[])
 
             BEqn.solve();
 
-            volScalarField rBA = 1.0/BEqn.A();
+            volScalarField rBA(1.0/BEqn.A());
 
             phiB = (fvc::interpolate(B) & mesh.Sf())
                 + fvc::ddtPhiCorr(rBA, B, phiB);

applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H

@@ -7,7 +7,7 @@
     phi = (fvc::interpolate(U) & mesh.Sf())
         + fvc::ddtPhiCorr(rAU, U, phi);
 
-    surfaceScalarField buoyancyPhi = rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf();
+    surfaceScalarField buoyancyPhi(rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf());
     phi -= buoyancyPhi;
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H

@@ -8,7 +8,7 @@
     phi = fvc::interpolate(U) & mesh.Sf();
     adjustPhi(phi, U, p_rgh);
 
-    surfaceScalarField buoyancyPhi = rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf();
+    surfaceScalarField buoyancyPhi(rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf());
     phi -= buoyancyPhi;
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)

applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H

@@ -5,7 +5,7 @@
     // pressure solution - done in 2 parts. Part 1:
     thermo.rho() -= psi*p_rgh;
 
-    volScalarField rAU = 1.0/UEqn.A();
+    volScalarField rAU(1.0/UEqn.A());
     surfaceScalarField rhorAUf("(rho*(1|A(U)))", fvc::interpolate(rho*rAU));
 
     U = rAU*UEqn.H();
@@ -16,7 +16,7 @@
       + fvc::ddtPhiCorr(rAU, rho, U, phi)
     );
 
-    surfaceScalarField buoyancyPhi = -rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf();
+    surfaceScalarField buoyancyPhi(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
     phi += buoyancyPhi;
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)

applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H

@@ -2,7 +2,7 @@
     rho = thermo.rho();
     rho.relax();
 
-    volScalarField rAU = 1.0/UEqn().A();
+    volScalarField rAU(1.0/UEqn().A());
     surfaceScalarField rhorAUf("(rho*(1|A(U)))", fvc::interpolate(rho*rAU));
 
     U = rAU*UEqn().H();
@@ -11,7 +11,7 @@
     phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
     bool closedVolume = adjustPhi(phi, U, p_rgh);
 
-    surfaceScalarField buoyancyPhi = rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf();
+    surfaceScalarField buoyancyPhi(rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
     phi -= buoyancyPhi;
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)

applications/solvers/heatTransfer/buoyantSimpleFoam_old/pEqn.H

@@ -1,7 +1,7 @@
 {
     rho = thermo.rho();
 
-    volScalarField rAU = 1.0/UEqn().A();
+    volScalarField rAU(1.0/UEqn().A());
     surfaceScalarField rhorAUf("(rho*(1|A(U)))", fvc::interpolate(rho*rAU));
 
     U = rAU*UEqn().H();
@@ -10,8 +10,10 @@
     phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
     bool closedVolume = adjustPhi(phi, U, p);
 
-    surfaceScalarField buoyancyPhi =
+    surfaceScalarField buoyancyPhi
-        rhorAUf*fvc::interpolate(rho)*(g & mesh.Sf());
+    (
+        rhorAUf*fvc::interpolate(rho)*(g & mesh.Sf())
+    );
     phi += buoyancyPhi;
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H

@@ -4,7 +4,7 @@
     rho = min(rho, rhoMax[i]);
     rho.relax();
 
-    volScalarField rAU = 1.0/UEqn().A();
+    volScalarField rAU(1.0/UEqn().A());
     surfaceScalarField rhorAUf("(rho*(1|A(U)))", fvc::interpolate(rho*rAU));
 
     U = rAU*UEqn().H();
@@ -15,7 +15,7 @@
     dimensionedScalar compressibility = fvc::domainIntegrate(psi);
     bool compressible = (compressibility.value() > SMALL);
 
-    surfaceScalarField buoyancyPhi = rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf();
+    surfaceScalarField buoyancyPhi(rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
     phi -= buoyancyPhi;
 
     // Solve pressure

applications/solvers/heatTransfer/chtMultiRegionFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.C

@@ -138,7 +138,7 @@ Foam::solidWallHeatFluxTemperatureFvPatchScalarField::K() const
         const symmTensorField& KWall =
             patch().lookupPatchField<volSymmTensorField, scalar>(KName_);
 
-        vectorField n = patch().nf();
+        vectorField n(patch().nf());
 
         return n & KWall & n;
     }
@@ -203,7 +203,7 @@ void Foam::solidWallHeatFluxTemperatureFvPatchScalarField::write
 
 namespace Foam
 {
-    makePatchTypeField
+    makeNonTemplatedPatchTypeField
     (
         fvPatchScalarField,
         solidWallHeatFluxTemperatureFvPatchScalarField

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleCourantNo.C

@@ -34,9 +34,11 @@ Foam::scalar Foam::compressibleCourantNo
     const surfaceScalarField& phi
 )
 {
-    scalarField sumPhi =
+    scalarField sumPhi
+    (
         fvc::surfaceSum(mag(phi))().internalField()
-       /rho.internalField();
+      / rho.internalField()
+    );
 
     scalar CoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
 

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H

@@ -5,7 +5,7 @@
 
     rho = thermo.rho();
 
-    volScalarField rAU = 1.0/UEqn().A();
+    volScalarField rAU(1.0/UEqn().A());
     surfaceScalarField rhorAUf("(rho*(1|A(U)))", fvc::interpolate(rho*rAU));
 
     U = rAU*UEqn().H();

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffNo.C

@@ -39,10 +39,12 @@ Foam::scalar Foam::solidRegionDiffNo
 
     //- Take care: can have fluid domains with 0 cells so do not test for
     //  zero internal faces.
-    surfaceScalarField KrhoCpbyDelta =
+    surfaceScalarField KrhoCpbyDelta
+    (
         mesh.surfaceInterpolation::deltaCoeffs()
       * fvc::interpolate(K)
-      / fvc::interpolate(Cprho);
+      / fvc::interpolate(Cprho)
+    );
 
     DiNum = gMax(KrhoCpbyDelta.internalField())*runTime.deltaT().value();
 
@@ -66,14 +68,16 @@ Foam::scalar Foam::solidRegionDiffNo
     scalar DiNum = 0.0;
     scalar meanDiNum = 0.0;
 
-    volScalarField K = mag(Kdirectional);
+    volScalarField K(mag(Kdirectional));
 
     //- Take care: can have fluid domains with 0 cells so do not test for
     //  zero internal faces.
-    surfaceScalarField KrhoCpbyDelta =
+    surfaceScalarField KrhoCpbyDelta
+    (
         mesh.surfaceInterpolation::deltaCoeffs()
       * fvc::interpolate(K)
-      / fvc::interpolate(Cprho);
+      / fvc::interpolate(Cprho)
+    );
 
     DiNum = gMax(KrhoCpbyDelta.internalField())*runTime.deltaT().value();
 

applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointOutletPressure/adjointOutletPressureFvPatchScalarField.C

@@ -121,7 +121,7 @@ void Foam::adjointOutletPressureFvPatchScalarField::write(Ostream& os) const
 
 namespace Foam
 {
-    makePatchTypeField
+    makeNonTemplatedPatchTypeField
     (
         fvPatchScalarField,
         adjointOutletPressureFvPatchScalarField

applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointOutletVelocity/adjointOutletVelocityFvPatchVectorField.C

@@ -96,10 +96,10 @@ void Foam::adjointOutletVelocityFvPatchVectorField::updateCoeffs()
     const fvPatchField<vector>& Up =
         patch().lookupPatchField<volVectorField, vector>("U");
 
-    scalarField Un = mag(patch().nf() & Up);
+    scalarField Un(mag(patch().nf() & Up));
-    vectorField UtHat = (Up - patch().nf()*Un)/(Un + SMALL);
+    vectorField UtHat((Up - patch().nf()*Un)/(Un + SMALL));
 
-    vectorField Uan = patch().nf()*(patch().nf() & patchInternalField());
+    vectorField Uan(patch().nf()*(patch().nf() & patchInternalField()));
 
     vectorField::operator=(phiap*patch().Sf()/sqr(patch().magSf()) + UtHat);
     //vectorField::operator=(Uan + UtHat);
@@ -119,7 +119,7 @@ void Foam::adjointOutletVelocityFvPatchVectorField::write(Ostream& os) const
 
 namespace Foam
 {
-    makePatchTypeField
+    makeNonTemplatedPatchTypeField
     (
         fvPatchVectorField,
         adjointOutletVelocityFvPatchVectorField

applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C

@@ -114,7 +114,7 @@ int main(int argc, char *argv[])
             solve(UEqn() == -fvc::grad(p));
 
             p.boundaryField().updateCoeffs();
-            volScalarField rAU = 1.0/UEqn().A();
+            volScalarField rAU(1.0/UEqn().A());
             U = rAU*UEqn().H();
             UEqn.clear();
             phi = fvc::interpolate(U) & mesh.Sf();
@@ -153,10 +153,13 @@ int main(int argc, char *argv[])
         {
             // Adjoint Momentum predictor
 
-            volVectorField adjointTransposeConvection = (fvc::grad(Ua) & U);
+            volVectorField adjointTransposeConvection((fvc::grad(Ua) & U));
-            //volVectorField adjointTransposeConvection = fvc::reconstruct
+            //volVectorField adjointTransposeConvection
             //(
+            //    fvc::reconstruct
+            //    (
             //        mesh.magSf()*(fvc::snGrad(Ua) & fvc::interpolate(U))
+            //    )
             //);
 
             zeroCells(adjointTransposeConvection, inletCells);
@@ -174,7 +177,7 @@ int main(int argc, char *argv[])
             solve(UaEqn() == -fvc::grad(pa));
 
             pa.boundaryField().updateCoeffs();
-            volScalarField rAUa = 1.0/UaEqn().A();
+            volScalarField rAUa(1.0/UaEqn().A());
             Ua = rAUa*UaEqn().H();
             UaEqn.clear();
             phia = fvc::interpolate(Ua) & mesh.Sf();

applications/solvers/incompressible/boundaryFoam/boundaryFoam.C

@@ -60,7 +60,7 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        fvVectorMatrix divR = turbulence->divDevReff(U);
+        fvVectorMatrix divR(turbulence->divDevReff(U));
         divR.source() = flowMask & divR.source();
 
         fvVectorMatrix UEqn

applications/solvers/incompressible/boundaryFoam/evaluateNearWall.H

@@ -1,7 +1,7 @@
 {
     // Evaluate near-wall behaviour
 
-    scalar nu = turbulence->nu().boundaryField()[patchId][faceId];
+    scalar nu = turbulence->nu()().boundaryField()[patchId][faceId];
     scalar nut = turbulence->nut()().boundaryField()[patchId][faceId];
     symmTensor R = turbulence->devReff()().boundaryField()[patchId][faceId];
     scalar epsilon = turbulence->epsilon()()[cellId];

applications/solvers/incompressible/boundaryFoam/interrogateWallPatches.H

@@ -13,7 +13,7 @@ forAll(patches, patchi)
 
     if (isA<wallFvPatch>(currPatch))
     {
-        const vectorField nf = currPatch.nf();
+        const vectorField nf(currPatch.nf());
 
         forAll(nf, facei)
         {
@@ -67,8 +67,10 @@ else
 label cellId = patches[patchId].faceCells()[faceId];
 
 // create position array for graph generation
-scalarField y =
+scalarField y
+(
     wallNormal
-  & (mesh.C().internalField() - mesh.C().boundaryField()[patchId][faceId]);
+  & (mesh.C().internalField() - mesh.C().boundaryField()[patchId][faceId])
+);
 
 Info<< "    Height to first cell centre y0 = " << y[cellId] << endl;

applications/solvers/incompressible/channelFoam/channelFoam.C

@@ -77,7 +77,7 @@ int main(int argc, char *argv[])
 
         // --- PISO loop
 
-        volScalarField rAU = 1.0/UEqn.A();
+        volScalarField rAU(1.0/UEqn.A());
 
         for (int corr=0; corr<nCorr; corr++)
         {

applications/solvers/incompressible/icoFoam/icoFoam.C

@@ -66,7 +66,7 @@ int main(int argc, char *argv[])
 
         for (int corr=0; corr<nCorr; corr++)
         {
-            volScalarField rAU = 1.0/UEqn.A();
+            volScalarField rAU(1.0/UEqn.A());
 
             U = rAU*UEqn.H();
             phi = (fvc::interpolate(U) & mesh.Sf())

applications/solvers/incompressible/nonNewtonianIcoFoam/nonNewtonianIcoFoam.C

@@ -69,7 +69,7 @@ int main(int argc, char *argv[])
 
         for (int corr=0; corr<nCorr; corr++)
         {
-            volScalarField rAU = 1.0/UEqn.A();
+            volScalarField rAU(1.0/UEqn.A());
 
             U = rAU*UEqn.H();
             phi = (fvc::interpolate(U) & mesh.Sf())

applications/solvers/incompressible/pimpleFoam/UEqn.H

@@ -9,7 +9,7 @@ tmp<fvVectorMatrix> UEqn
 
 UEqn().relax();
 
-volScalarField rAU = 1.0/UEqn().A();
+volScalarField rAU(1.0/UEqn().A());
 
 if (momentumPredictor)
 {

applications/solvers/incompressible/pisoFoam/pisoFoam.C

@@ -79,7 +79,7 @@ int main(int argc, char *argv[])
 
             for (int corr=0; corr<nCorr; corr++)
             {
-                volScalarField rAU = 1.0/UEqn.A();
+                volScalarField rAU(1.0/UEqn.A());
 
                 U = rAU*UEqn.H();
                 phi = (fvc::interpolate(U) & mesh.Sf())

applications/solvers/incompressible/porousSimpleFoam/UEqn.H

@@ -22,7 +22,7 @@
         trTU = inv(tTU());
         trTU().rename("rAU");
 
-        volVectorField gradp = fvc::grad(p);
+        volVectorField gradp(fvc::grad(p));
 
         for (int UCorr=0; UCorr<nUCorr; UCorr++)
         {

applications/solvers/incompressible/shallowWaterFoam/CourantNo.H

@@ -36,9 +36,11 @@ scalar waveCoNum = 0.0;
 
 if (mesh.nInternalFaces())
 {
-    scalarField sumPhi =
+    scalarField sumPhi
+    (
         fvc::surfaceSum(mag(phi))().internalField()
-       /h.internalField();
+      / h.internalField()
+    );
 
     CoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
 

applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C

@@ -89,10 +89,10 @@ int main(int argc, char *argv[])
             // --- PISO loop
             for (int corr=0; corr<nCorr; corr++)
             {
-                volScalarField rAU = 1.0/hUEqn.A();
+                volScalarField rAU(1.0/hUEqn.A());
-                surfaceScalarField ghrAUf = magg*fvc::interpolate(h*rAU);
+                surfaceScalarField ghrAUf(magg*fvc::interpolate(h*rAU));
 
-                surfaceScalarField phih0 = ghrAUf*mesh.magSf()*fvc::snGrad(h0);
+                surfaceScalarField phih0(ghrAUf*mesh.magSf()*fvc::snGrad(h0));
 
                 if (rotating)
                 {

applications/solvers/incompressible/simpleFoam/pEqn.H

@@ -1,6 +1,6 @@
     p.boundaryField().updateCoeffs();
 
-    volScalarField AU = UEqn().A();
+    volScalarField AU(UEqn().A());
     U = UEqn().H()/AU;
     UEqn.clear();
     phi = fvc::interpolate(U) & mesh.Sf();

applications/solvers/lagrangian/coalChemistryFoam/YEqn.H

@@ -12,7 +12,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
 
 {
     label inertIndex = -1;
-    volScalarField Yt = 0.0*Y[0];
+    volScalarField Yt(0.0*Y[0]);
 
     forAll(Y, i)
     {

applications/solvers/lagrangian/coalChemistryFoam/chemistry.H

@@ -10,10 +10,14 @@
     // turbulent time scale
     if (turbulentReaction)
     {
-        DimensionedField<scalar, volMesh> tk =
+        DimensionedField<scalar, volMesh> tk
-            Cmix*sqrt(turbulence->muEff()/rho/turbulence->epsilon());
+        (
-        DimensionedField<scalar, volMesh> tc =
+            Cmix*sqrt(turbulence->muEff()/rho/turbulence->epsilon())
-            chemistry.tc()().dimensionedInternalField();
+        );
+        DimensionedField<scalar, volMesh> tc
+        (
+            chemistry.tc()().dimensionedInternalField()
+        );
 
         // Chalmers PaSR model
         kappa = (runTime.deltaT() + tc)/(runTime.deltaT() + tc + tk);

applications/solvers/lagrangian/coalChemistryFoam/pEqn.H

@@ -1,6 +1,6 @@
 rho = thermo.rho();
 
-volScalarField rAU = 1.0/UEqn.A();
+volScalarField rAU(1.0/UEqn.A());
 U = rAU*UEqn.H();
 
 if (transonic)

applications/solvers/lagrangian/incompressibleUncoupledKinematicParcelDyMFoam/createFields.H

@@ -51,7 +51,7 @@
 
     singlePhaseTransportModel laminarTransport(U, phi);
 
-    const volScalarField nu = laminarTransport.nu();
+    const volScalarField nu(laminarTransport.nu());
 
     autoPtr<incompressible::turbulenceModel> turbulence
     (

applications/solvers/lagrangian/incompressibleUncoupledKinematicParcelFoam/createFields.H

@@ -51,7 +51,7 @@
 
     singlePhaseTransportModel laminarTransport(U, phi);
 
-    const volScalarField nu = laminarTransport.nu();
+    const volScalarField nu(laminarTransport.nu());
 
     autoPtr<incompressible::turbulenceModel> turbulence
     (

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/YEqn.H

@@ -13,7 +13,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
 
 {
     label inertIndex = -1;
-    volScalarField Yt = 0.0*Y[0];
+    volScalarField Yt(0.0*Y[0]);
 
     forAll(Y, i)
     {

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/chemistry.H

@@ -10,10 +10,14 @@
     // turbulent time scale
     if (turbulentReaction)
     {
-        DimensionedField<scalar, volMesh> tk =
+        DimensionedField<scalar, volMesh> tk
-            Cmix*sqrt(turbulence->muEff()/rho/turbulence->epsilon());
+        (
-        DimensionedField<scalar, volMesh> tc =
+            Cmix*sqrt(turbulence->muEff()/rho/turbulence->epsilon())
-            chemistry.tc()().dimensionedInternalField();
+        );
+        DimensionedField<scalar, volMesh> tc
+        (
+            chemistry.tc()().dimensionedInternalField()
+        );
 
         // Chalmers PaSR model
         kappa = (runTime.deltaT() + tc)/(runTime.deltaT() + tc + tk);

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H

@@ -5,7 +5,7 @@
     // pressure solution - done in 2 parts. Part 1:
     thermo.rho() -= psi*p;
 
-    volScalarField rAU = 1.0/UEqn.A();
+    volScalarField rAU(1.0/UEqn.A());
     U = rAU*UEqn.H();
 
     if (pZones.size() > 0)

applications/solvers/lagrangian/reactingParcelFilmFoam/YEqn.H

@@ -12,7 +12,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
 
 {
     label inertIndex = -1;
-    volScalarField Yt = 0.0*Y[0];
+    volScalarField Yt(0.0*Y[0]);
 
     forAll(Y, i)
     {

applications/solvers/lagrangian/reactingParcelFilmFoam/chemistry.H

@@ -10,10 +10,14 @@
     // turbulent time scale
     if (turbulentReaction)
     {
-        DimensionedField<scalar, volMesh> tk =
+        DimensionedField<scalar, volMesh> tk
-            Cmix*sqrt(turbulence->muEff()/rho/turbulence->epsilon());
+        (
-        DimensionedField<scalar, volMesh> tc =
+            Cmix*sqrt(turbulence->muEff()/rho/turbulence->epsilon())
-            chemistry.tc()().dimensionedInternalField();
+        );
+        DimensionedField<scalar, volMesh> tc
+        (
+            chemistry.tc()().dimensionedInternalField()
+        );
 
         // Chalmers PaSR model
         kappa = (runTime.deltaT() + tc)/(runTime.deltaT() + tc + tk);

applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H

@@ -1,6 +1,6 @@
 rho = thermo.rho();
 
-volScalarField rAU = 1.0/UEqn.A();
+volScalarField rAU(1.0/UEqn.A());
 U = rAU*UEqn.H();
 
 if (transonic)

applications/solvers/lagrangian/reactingParcelFoam/YEqn.H

@@ -12,7 +12,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
 
 {
     label inertIndex = -1;
-    volScalarField Yt = 0.0*Y[0];
+    volScalarField Yt(0.0*Y[0]);
 
     forAll(Y, i)
     {

applications/solvers/lagrangian/reactingParcelFoam/chemistry.H

@@ -10,10 +10,14 @@
     // turbulent time scale
     if (turbulentReaction)
     {
-        DimensionedField<scalar, volMesh> tk =
+        DimensionedField<scalar, volMesh> tk
-            Cmix*sqrt(turbulence->muEff()/rho/turbulence->epsilon());
+        (
-        DimensionedField<scalar, volMesh> tc =
+            Cmix*sqrt(turbulence->muEff()/rho/turbulence->epsilon())
-            chemistry.tc()().dimensionedInternalField();
+        );
+        DimensionedField<scalar, volMesh> tc
+        (
+            chemistry.tc()().dimensionedInternalField()
+        );
 
         // Chalmers PaSR model
         kappa = (runTime.deltaT() + tc)/(runTime.deltaT() + tc + tk);

applications/solvers/lagrangian/reactingParcelFoam/pEqn.H

@@ -1,6 +1,6 @@
 rho = thermo.rho();
 
-volScalarField rAU = 1.0/UEqn.A();
+volScalarField rAU(1.0/UEqn.A());
 U = rAU*UEqn.H();
 
 if (transonic)