diff --git a/tutorials/solidMechanics/elasticThermalSolidFoam/hotCylinder/analyticalHotCylinder/analyticalHotCylinder.C b/tutorials/solidMechanics/elasticThermalSolidFoam/hotCylinder/analyticalHotCylinder/analyticalHotCylinder.C
index 39a3dd7a4..c89c7c5f3 100644
--- a/tutorials/solidMechanics/elasticThermalSolidFoam/hotCylinder/analyticalHotCylinder/analyticalHotCylinder.C
+++ b/tutorials/solidMechanics/elasticThermalSolidFoam/hotCylinder/analyticalHotCylinder/analyticalHotCylinder.C
@@ -26,7 +26,7 @@ Description
     Generate analytical solution for a thick-walled cylinder with a
     temperature gradient.
     Temperature field T and stress field sigma and generated.
-    Based on solution outlined in Timoshenko, Theory of Elasticity. 
+    Based on solution outlined in Timoshenko, Theory of Elasticity.
 
 Author
     philip.cardiff@ucd.ie
@@ -43,351 +43,275 @@ int main(int argc, char *argv[])
 #   include "createTime.H"
 #   include "createMesh.H"
 
-  runTime++;
+    runTime++;
 
-  Info << "Writing analytical solution for a plain strain cylinder with concentric hole,\nwhere"
-       << "\n\tinner radius = 0.5"
-       << "\n\touter radius = 0.7"
-       << "\n\tinner temperature = 100"
-       << "\n\touter temperature = 0"
-       << "\n\tinner pressure = 0"
-       << "\n\touter pressure = 0"
-       << "\n\tE = 200e9"
-       << "\n\tu = 0.3"
-       << "\n\talpha = 1e-5"
-       << nl << endl;
- 
-  //- inner and outer radii and temperatures
-  scalar a = 0.5;
-  scalar b = 0.7;
-  scalar Ti = 100;
-  scalar To = 0;
+    Info<< "Writing analytical solution for a plain strain cylinder "
+        << "with concentric hole,\nwhere"
+        << "\n\tinner radius = 0.5"
+        << "\n\touter radius = 0.7"
+        << "\n\tinner temperature = 100"
+        << "\n\touter temperature = 0"
+        << "\n\tinner pressure = 0"
+        << "\n\touter pressure = 0"
+        << "\n\tE = 200e9"
+        << "\n\tu = 0.3"
+        << "\n\talpha = 1e-5"
+        << nl << endl;
 
-  //- mechanical and thermal properties
-  scalar E = 200e9;
-  scalar nu = 0.3;
-  scalar alpha = 1e-5;
+    //- inner and outer radii and temperatures
+    scalar a = 0.5;
+    scalar b = 0.7;
+    scalar Ti = 100;
+    scalar To = 0;
 
-  //- create T field
-  volScalarField T
+    //- mechanical and thermal properties
+    scalar E = 200e9;
+    scalar nu = 0.3;
+    scalar alpha = 1e-5;
+
+    const volVectorField& C = mesh.C();
+
+    //- radial coordinate
+    volScalarField radii
     (
-     IOobject
-     (
-      "analyticalT",
-      runTime.timeName(),
-      mesh,
-      IOobject::NO_READ,
-      IOobject::AUTO_WRITE
-      ),
-     mesh,
-     dimensionedScalar("zero", dimTemperature, 0.0)
-     );
-  
-  const volVectorField& C = mesh.C();
+        sqrt
+        (
+            sqr(C.component(vector::X))
+          + sqr(C.component(vector::Y))
+        )/dimensionedScalar("one", dimLength, 1)
+    );
 
-  //- radial coordinate
-  volScalarField radii =
-    C.component(vector::X)*C.component(vector::X) + C.component(vector::Y)*C.component(vector::Y);
-  forAll(radii.internalField(), celli)
-    {
-      radii.internalField()[celli] = ::sqrt(radii.internalField()[celli]);
-    }
-  forAll(radii.boundaryField(), patchi)
-    {
-      forAll(radii.boundaryField()[patchi], facei)
-	{
-	  radii.boundaryField()[patchi][facei] = ::sqrt(radii.boundaryField()[patchi][facei]);
-	}
-    }
+    const scalarField& rIn = radii.internalField();
 
-  forAll(T.internalField(), celli)
-    {
-      const scalar& r = radii[celli];
-
-      T.internalField()[celli] =
-	( (Ti-To)/Foam::log(b/a) ) * Foam::log(b/r);
-    }
-
-  forAll(T.boundaryField(), patchi)
-    {
-      forAll(T.boundaryField()[patchi], facei)
-	{
-	  const scalar& r = radii.boundaryField()[patchi][facei];
-
-	  T.boundaryField()[patchi][facei] =
-	    ( (Ti-To)/Foam::log(b/a) ) * Foam::log(b/r);
-	}
-    }
-
-  //- write temperature file
-  Info << "Writing analytical termpature field" << endl;
-  T.write();
-
-
-  //- create sigma field
-  volScalarField sigmaR
+    Info << "Writing analytical termpature field" << endl;
+    //- create T field
+    volScalarField T
     (
-     IOobject
-     (
-      "sigmaR",
-      runTime.timeName(),
-      mesh,
-      IOobject::NO_READ,
-      IOobject::AUTO_WRITE
-      ),
-     mesh,
-     dimensionedScalar("zero", dimForce/dimArea, 0.0)
-     );
+        IOobject
+        (
+            "analyticalT",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        ((Ti - To)/Foam::log(b/a))*Foam::log(b/radii)
+    );
+    T.write();
 
-  forAll(sigmaR.internalField(), celli)
-    {
-      const scalar& r = radii.internalField()[celli];
-
-      sigmaR.internalField()[celli] =
-	( (alpha*E*(Ti-To))/(2*(1-nu)*Foam::log(b/a)) ) *
-	(-Foam::log(b/r) -( a*a/(b*b - a*a))*(1 - (b*b)/(r*r))*Foam::log(b/a));
-    }
-
-  forAll(sigmaR.boundaryField(), patchi)
-    {
-      forAll(sigmaR.boundaryField()[patchi], facei)
-	{
-	  const scalar& r = radii.boundaryField()[patchi][facei];
-	  
-	  sigmaR.boundaryField()[patchi][facei] =
-	    ( (alpha*E*(Ti-To))/(2*(1-nu)*Foam::log(b/a)) ) *
-	    ( -Foam::log(b/r) - ( a*a/(b*b - a*a))*(1 - (b*b)/(r*r))*Foam::log(b/a) );
-	}
-    }
-
-  //- write temperature file
-  Info << "\nWriting analytical sigmaR field" << endl;
-  sigmaR.write();
-
-
-  volScalarField sigmaTheta
+    //- create sigma field
+    Info << "\nWriting analytical sigmaR field" << endl;
+    volScalarField sigmaR
     (
-     IOobject
-     (
-      "sigmaTheta",
-      runTime.timeName(),
-      mesh,
-      IOobject::NO_READ,
-      IOobject::AUTO_WRITE
-      ),
-     mesh,
-     dimensionedScalar("zero", dimForce/dimArea, 0.0)
-     );
-
-  forAll(sigmaTheta.internalField(), celli)
-    {
-      const scalar& r = radii.internalField()[celli];
-
-      sigmaTheta.internalField()[celli] =
-	( (alpha*E*(Ti-To))/(2*(1-nu)*Foam::log(b/a)) ) *
-	(1 -Foam::log(b/r) - ( a*a/(b*b - a*a))*(1 + (b*b)/(r*r))*Foam::log(b/a) );
-    }
-
-  forAll(sigmaTheta.boundaryField(), patchi)
-    {
-      forAll(sigmaTheta.boundaryField()[patchi], facei)
-	{
-	  const scalar& r = radii.boundaryField()[patchi][facei];
-	  
-	  sigmaTheta.boundaryField()[patchi][facei] =
-	    ( (alpha*E*(Ti-To))/(2*(1-nu)*Foam::log(b/a)) ) *
-	    (1 -Foam::log(b/r) - ( a*a/(b*b - a*a))*(1 + (b*b)/(r*r))*Foam::log(b/a) );
-	}
-    }
+        IOobject
+        (
+            "sigmaR",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        ((alpha*E*(Ti - To))/(2*(1 - nu)*Foam::log(b/a)))*
+        (
+            -Foam::log(b/radii)
+          - (sqr(a)/(sqr(b) - sqr(a)))*(1 - sqr(b)/sqr(radii))*Foam::log(b/a)
+        )
+    );
+    sigmaR.write();
 
 
-  //- write temperature file
-  Info << "\nWriting analytical sigmaTheta field" << endl;
-  sigmaTheta.write();
-
-  volScalarField sigmaZ
+    Info << "\nWriting analytical sigmaTheta field" << endl;
+    volScalarField sigmaTheta
     (
-     IOobject
-     (
-      "sigmaZ",
-      runTime.timeName(),
-      mesh,
-      IOobject::NO_READ,
-      IOobject::AUTO_WRITE
-      ),
-     mesh,
-     dimensionedScalar("zero", dimForce/dimArea, 0.0)
-     );
+        IOobject
+        (
+            "sigmaTheta",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        ((alpha*E*(Ti - To))/(2*(1 - nu)*Foam::log(b/a)))*
+        (
+            1 - Foam::log(b/radii)
+          - (sqr(a)/(sqr(b) - sqr(a)))*(1 + sqr(b)/sqr(radii))*Foam::log(b/a)
+        )
+    );
+    sigmaTheta.write();
 
-  forAll(sigmaZ.internalField(), celli)
-    {
-      //- Timoshenko says this but I am not sure I am not sure the BCs in
-      //- the z direction
-      // sigmaZ.internalField()[celli] =
-      // 	( (alpha*E*(Ti-To))/(2*(1-nu)*Foam::log(b/a)) ) *
-      // 	(1 - 2*Foam::log(b/r) - ( 2*a*a/(b*b - a*a))*Foam::log(b/a));
-
-      sigmaZ.internalField()[celli] = 
-	0.3*(sigmaR.internalField()[celli] + sigmaTheta.internalField()[celli])
-	- E*alpha*(T.internalField()[celli]);
-    }
-
-  forAll(sigmaZ.boundaryField(), patchi)
-    {
-      forAll(sigmaZ.boundaryField()[patchi], facei)
-	{
-	  //- Timoshenko says this but I am not sure I am not sure the BCs in
-	  //- the z direction
-	  //sigmaZ.boundaryField()[patchi][facei] =
-	  //( (alpha*E*(Ti-To))/(2*(1-nu)*Foam::log(b/a)) ) *
-	  //(1 - 2*Foam::log(b/r) - ( 2*a*a/(b*b - a*a))*Foam::log(b/a));
-
-	  //-for general 2-D plain strain problems, the axial stress is given by this:
-	  sigmaZ.boundaryField()[patchi][facei] = 
-	    nu*(sigmaR.boundaryField()[patchi][facei] + sigmaTheta.boundaryField()[patchi][facei])
-	    - E*alpha*(T.boundaryField()[patchi][facei]);
-	}
-    }
-
-
-  //- write temperature file
-  Info << "\nWriting analytical sigmaZ field" << endl;
-  sigmaZ.write();
-
-
-  //- create analytical sigma tensor
-
-  //- create theta field
-  volScalarField theta
+    Info << "\nWriting analytical sigmaZ field" << endl;
+    volScalarField sigmaZ
     (
-     IOobject
-     (
-      "theta",
-      runTime.timeName(),
-      mesh,
-      IOobject::NO_READ,
-      IOobject::NO_WRITE
-      ),
-     mesh,
-     dimensionedScalar("zero", dimless, 0.0)
-     );
-  forAll(theta.internalField(), celli)
-    {
-      const scalar& x = mesh.C().internalField()[celli][vector::X];
-      const scalar& y = mesh.C().internalField()[celli][vector::Y];
+        IOobject
+        (
+            "sigmaZ",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        // Timoshenko says this but I am not sure I am not sure the BCs in
+        // the z direction
+        // ((alpha*E*(Ti - To))/(2*(1 - nu)*Foam::log(b/a)))*
+        // (1 - 2*Foam::log(b/radii) - ( 2*sqr(a)/(sqr(b) - sqr(a)))*Foam::log(b/a));
+        0.3*(sigmaR + sigmaTheta) - E*alpha*(T)
+    );
+    sigmaZ.write();
 
-      theta.internalField()[celli] = Foam::atan(y/x);
-    }
-
-  forAll(theta.boundaryField(), patchi)
-    {
-      forAll(theta.boundaryField()[patchi], facei)
-	{
-	  const scalar& x = mesh.C().boundaryField()[patchi][facei][vector::X];
-	  const scalar& y = mesh.C().boundaryField()[patchi][facei][vector::Y];
-	  
-	  theta.boundaryField()[patchi][facei] = Foam::atan(y/x);
-	}
-    }
-
-  //- rotation matrix to convert polar stresses to cartesian
-  volTensorField rotMat
+    //- create theta field
+    volScalarField yOverX
     (
-     IOobject
-     (
-      "rotMat",
-      runTime.timeName(),
-      mesh,
-      IOobject::NO_READ,
-      IOobject::NO_WRITE
-      ),
-     mesh,
-     dimensionedTensor("zero", dimless, tensor::zero)
-     );
+        "yOverX",
+        Foam::max
+        (
+            scalar(-1),
+            Foam::min
+            (
+                scalar(1),
+                mesh.C().component(vector::Y)/
+                stabilise
+                (
+                    mesh.C().component(vector::X),
+                    dimensionedScalar("small", dimLength, SMALL)
+                )
+            )
+        )
+    );
 
-  forAll(rotMat.internalField(), celli)
-    {
-      const scalar& t = theta.internalField()[celli];
-
-      rotMat.internalField()[celli] = tensor(::cos(t), ::sin(t), 0,
-					     -::sin(t), ::cos(t), 0,
-					     0, 0, 1);
-    }
-
-  forAll(rotMat.boundaryField(), patchi)
-    {
-      forAll(rotMat.boundaryField()[patchi], facei)
-	{
-	  const scalar& t = theta.boundaryField()[patchi][facei];
-	  
-	  rotMat.boundaryField()[patchi][facei] =  tensor(::cos(t), ::sin(t), 0,
-					     -::sin(t), ::cos(t), 0,
-					     0, 0, 1);
-	}
-    }
-
-  volSymmTensorField sigma
+    volScalarField theta
     (
-     IOobject
-     (
-      "analyticalSigma",
-      runTime.timeName(),
-      mesh,
-      IOobject::NO_READ,
-      IOobject::AUTO_WRITE
-      ),
-     mesh,
-     dimensionedSymmTensor("zero", dimForce/dimArea, symmTensor::zero)
-     );
+        IOobject
+        (
+            "theta",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        Foam::atan(yOverX)
+    );
 
-  forAll(sigma.internalField(), celli)
+    //- rotation matrix to convert polar stresses to cartesian
+    volTensorField rotMat
+    (
+        IOobject
+        (
+            "rotMat",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedTensor("zero", dimless, tensor::zero)
+    );
+
+    tensorField& rotMatIn = rotMat.internalField();
+    const scalarField tIn = theta.internalField();
+
+    forAll (rotMatIn, celli)
     {
-      const scalar& r = sigmaR.internalField()[celli];
-      const scalar& t = sigmaTheta.internalField()[celli];
-      const scalar& z = sigmaZ.internalField()[celli];
+        const scalar& t = tIn[celli];
 
-      const tensor& rot = rotMat.internalField()[celli];
-
-      symmTensor sigmaCart(r, 0, 0,
-			      t, 0,
-			         z);
-
-      sigma.internalField()[celli] =
-	    symm(rot.T() & sigmaCart & rot);
-
-      //-for general 2-D plain strain problems, the axial stress is given by this:
-      //- (which is not equal to the solution by Timoshenko... hmmmnn)
-//       sigma.internalField()[celli][symmTensor::ZZ] = 
-// 	0.3*(sigma.internalField()[celli][symmTensor::XX] + sigma.internalField()[celli][symmTensor::YY])
-// 	- E*alpha*(T.internalField()[celli]);
-    }
-
-  forAll(sigma.boundaryField(), patchi)
-    {
-      forAll(sigma.boundaryField()[patchi], facei)
-	{
-	  const scalar& r = sigmaR.boundaryField()[patchi][facei];
-	  const scalar& t = sigmaTheta.boundaryField()[patchi][facei];
-	  const scalar& z = sigmaZ.boundaryField()[patchi][facei];
-
-	  const tensor& rot = rotMat.boundaryField()[patchi][facei];
-
-	  symmTensor sigmaCart(r, 0, 0,
-		                  t, 0,
-			             z);
-	  sigma.boundaryField()[patchi][facei] =
-	    symm(rot.T() & sigmaCart & rot);
-	}
+        rotMatIn[celli] =
+            tensor
+            (
+                Foam::cos(t),  Foam::sin(t), 0,
+               -Foam::sin(t),  Foam::cos(t), 0,
+                0, 0, 1
+            );
     }
 
 
+    forAll (rotMat.boundaryField(), patchi)
+    {
+        forAll (rotMat.boundaryField()[patchi], facei)
+        {
+            const scalar& t = theta.boundaryField()[patchi][facei];
 
-  Info << "\nWriting analytical sigma tensor" << endl;
-  sigma.write();
+            rotMat.boundaryField()[patchi][facei] =
+                tensor
+                (
+                    Foam::cos(t),  Foam::sin(t), 0,
+                   -Foam::sin(t),  Foam::cos(t), 0,
+                    0, 0, 1
+                );
+        }
+    }
 
-  Info << nl << "End" << endl;
-      
-  return 0;
+    volSymmTensorField sigma
+    (
+        IOobject
+        (
+            "analyticalSigma",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedSymmTensor("zero", dimForce/dimArea, symmTensor::zero)
+    );
+
+    {
+        symmTensorField& sigmaIn = sigma.internalField();
+
+        const scalarField& rIn = sigmaR.internalField();
+        const scalarField& tIn = sigmaTheta.internalField();
+        const scalarField& zIn = sigmaZ.internalField();
+
+        forAll (sigmaIn, celli)
+        {
+            symmTensor sigmaCart
+            (
+                rIn[celli], 0, 0,
+                tIn[celli], 0,
+                zIn[celli]
+            );
+
+            const tensor& rot = rotMatIn[celli];
+
+            sigmaIn[celli] = symm(rot.T() & sigmaCart & rot);
+
+            // for general 2-D plain strain problems, the axial stress is:
+            // (which is not equal to the solution by Timoshenko... hmmmnn)
+            //       sigmaIn[celli][symmTensor::ZZ] =
+            //           0.3*(sigmaIn[celli][symmTensor::XX]
+            //         + sigmaIn[celli][symmTensor::YY])
+            //         - E*alpha*(T.internalField()[celli]);
+        }
+    }
+
+    forAll (sigma.boundaryField(), patchi)
+    {
+        symmTensorField& pSigma = sigma.boundaryField()[patchi];
+        const scalarField& pR = sigmaR.boundaryField()[patchi];
+        const scalarField& pT = sigmaTheta.boundaryField()[patchi];
+        const scalarField& pZ = sigmaZ.boundaryField()[patchi];
+
+        const tensorField pRot = rotMat.boundaryField()[patchi];
+
+        forAll (pSigma, facei)
+        {
+            const tensor& rot = pRot[facei];
+
+            symmTensor sigmaCart
+            (
+                pR[facei], 0, 0,
+                pT[facei], 0,
+                pZ[facei]
+            );
+
+            pSigma[facei] = symm(rot.T() & sigmaCart & rot);
+        }
+    }
+
+    Info << "\nWriting analytical sigma tensor" << endl;
+    sigma.write();
+
+    Info << nl << "End" << endl;
+
+    return 0;
 }