diff --git a/applications/solvers/solidMechanics/elasticIncrSolidFoam/calculateRelativeResidual.H b/applications/solvers/solidMechanics/elasticIncrSolidFoam/calculateRelativeResidual.H
index d5eccd36b..36430dbbb 100644
--- a/applications/solvers/solidMechanics/elasticIncrSolidFoam/calculateRelativeResidual.H
+++ b/applications/solvers/solidMechanics/elasticIncrSolidFoam/calculateRelativeResidual.H
@@ -16,7 +16,6 @@
             (
                 DU.internalField()
               - DU.prevIter().internalField()
-            )
-           /magDU
+            )/magDU
         );
 }
diff --git a/applications/solvers/solidMechanics/elasticIncrSolidFoam/createFields.H b/applications/solvers/solidMechanics/elasticIncrSolidFoam/createFields.H
index 659025e96..c0ec0e32f 100644
--- a/applications/solvers/solidMechanics/elasticIncrSolidFoam/createFields.H
+++ b/applications/solvers/solidMechanics/elasticIncrSolidFoam/createFields.H
@@ -24,7 +24,7 @@
             IOobject::AUTO_WRITE
         ),
         mesh,
-        dimensionedVector("zero", dimLength, vector::zero)	
+        dimensionedVector("zero", dimLength, vector::zero)
     );
 
     volSymmTensorField DEpsilon
diff --git a/applications/solvers/solidMechanics/elasticIncrSolidFoam/elasticIncrSolidFoam.C b/applications/solvers/solidMechanics/elasticIncrSolidFoam/elasticIncrSolidFoam.C
index cc059029b..4f6a0497f 100644
--- a/applications/solvers/solidMechanics/elasticIncrSolidFoam/elasticIncrSolidFoam.C
+++ b/applications/solvers/solidMechanics/elasticIncrSolidFoam/elasticIncrSolidFoam.C
@@ -79,7 +79,7 @@ int main(int argc, char *argv[])
 
 #           include "calculateDivDSigmaExp.H"
 
-            //- linear momentum equation
+            // Linear momentum equation
             fvVectorMatrix DUEqn
             (
                 fvm::d2dt2(rho, DU)
@@ -127,7 +127,8 @@ int main(int argc, char *argv[])
          && ++iCorr < nCorr
         );
 
-        Info<< nl << "Time " << runTime.value() << ", Solving for " << DU.name()
+        Info<< nl << "Time " << runTime.value()
+            << ", Solving for " << DU.name()
             << ", Initial residual = " << initialResidual
             << ", Final residual = " << solverPerf.initialResidual()
             << ", Relative residual = " << relativeResidual
diff --git a/applications/solvers/solidMechanics/elasticIncrSolidFoam/readDivDSigmaExpMethod.H b/applications/solvers/solidMechanics/elasticIncrSolidFoam/readDivDSigmaExpMethod.H
index c75f4c342..de2a549e7 100644
--- a/applications/solvers/solidMechanics/elasticIncrSolidFoam/readDivDSigmaExpMethod.H
+++ b/applications/solvers/solidMechanics/elasticIncrSolidFoam/readDivDSigmaExpMethod.H
@@ -1,9 +1,22 @@
-//- how explicit component of sigma is to be calculated
-word divDSigmaExpMethod(mesh.solutionDict().subDict("solidMechanics").lookup("divSigmaExp"));
-Info << "Selecting divSigmaExp calculation method " << divDSigmaExpMethod <<  endl;
-if(divDSigmaExpMethod != "standard" && divDSigmaExpMethod != "surface" && divDSigmaExpMethod != "decompose" && divDSigmaExpMethod != "laplacian")
-  {
-    FatalError << "divSigmaExp method " << divDSigmaExpMethod << " not found!" << nl
-	       << "valid methods are:\nstandard\nsurface\ndecompose\nlaplacian"
-	       << exit(FatalError);
-  }
+    const dictionary& stressControl =
+        mesh.solutionDict().subDict("solidMechanics");
+
+    // Read how explicit component of sigma is to be calculated
+    word divDSigmaExpMethod(stressControl.lookup("divSigmaExp"));
+
+    Info<< "Selecting divSigmaExp calculation method: "
+        << divDSigmaExpMethod <<  endl;
+
+    if
+    (
+        divDSigmaExpMethod != "standard"
+     && divDSigmaExpMethod != "surface"
+     && divDSigmaExpMethod != "decompose"
+     && divDSigmaExpMethod != "laplacian"
+    )
+    {
+        FatalErrorIn(args.executable())
+            << "divSigmaExp method " << divDSigmaExpMethod << " not found!  "
+	        << "valid methods are:\nstandard\nsurface\ndecompose\nlaplacian"
+	        << exit(FatalError);
+    }
diff --git a/applications/solvers/solidMechanics/elasticIncrSolidFoam/writeFields.H b/applications/solvers/solidMechanics/elasticIncrSolidFoam/writeFields.H
index 574bd09a9..f666ec12b 100644
--- a/applications/solvers/solidMechanics/elasticIncrSolidFoam/writeFields.H
+++ b/applications/solvers/solidMechanics/elasticIncrSolidFoam/writeFields.H
@@ -1,36 +1,34 @@
 if (runTime.outputTime())
-  {
+{
     volScalarField epsilonEq
-      (
-       IOobject
-       (
-	"epsilonEq",
-	runTime.timeName(),
-	mesh,
-	IOobject::NO_READ,
-	IOobject::AUTO_WRITE
-	),
-       sqrt((2.0/3.0)*magSqr(dev(epsilon)))
-       );
+    (
+        IOobject
+        (
+            "epsilonEq",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        sqrt((2.0/3.0)*magSqr(dev(epsilon)))
+    );
 
-    Info<< "Max epsilonEq = " << max(epsilonEq).value()
-	<< endl;
+    Info<< "Max epsilonEq = " << max(epsilonEq).value() << endl;
 
     volScalarField sigmaEq
-      (
-       IOobject
-       (
-	"sigmaEq",
-	runTime.timeName(),
-	mesh,
-	IOobject::NO_READ,
-	IOobject::AUTO_WRITE
-	),
-       sqrt((3.0/2.0)*magSqr(dev(sigma)))
-       );
-    
-    Info<< "Max sigmaEq = " << max(sigmaEq).value()
-	<< endl;
-    
+    (
+        IOobject
+        (
+            "sigmaEq",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        sqrt((3.0/2.0)*magSqr(dev(sigma)))
+    );
+
+    Info<< "Max sigmaEq = " << max(sigmaEq).value() << endl;
+
     runTime.write();
-  }
+}
diff --git a/applications/solvers/solidMechanics/elasticNonLinTLSolidFoam/createFields.H b/applications/solvers/solidMechanics/elasticNonLinTLSolidFoam/createFields.H
index 81a6f557f..0a420cee6 100644
--- a/applications/solvers/solidMechanics/elasticNonLinTLSolidFoam/createFields.H
+++ b/applications/solvers/solidMechanics/elasticNonLinTLSolidFoam/createFields.H
@@ -25,11 +25,11 @@
             IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
-	mesh,
-	dimensionedSymmTensor("zero", dimless, symmTensor::zero)
+        mesh,
+        dimensionedSymmTensor("zero", dimless, symmTensor::zero)
     );
 
-    //- second Piola-Kirchhoff stress tensor
+    //- Second Piola-Kirchhoff stress tensor
     volSymmTensorField sigma
     (
         IOobject
diff --git a/applications/solvers/solidMechanics/elasticNonLinTLSolidFoam/elasticNonLinTLSolidFoam.C b/applications/solvers/solidMechanics/elasticNonLinTLSolidFoam/elasticNonLinTLSolidFoam.C
index 749a9e238..c6479e3cd 100644
--- a/applications/solvers/solidMechanics/elasticNonLinTLSolidFoam/elasticNonLinTLSolidFoam.C
+++ b/applications/solvers/solidMechanics/elasticNonLinTLSolidFoam/elasticNonLinTLSolidFoam.C
@@ -44,21 +44,21 @@ Author
 
 int main(int argc, char *argv[])
 {
-# include "setRootCase.H"
-# include "createTime.H"
-# include "createMesh.H"
-# include "createFields.H"
-# include "createSolidInterfaceNonLin.H"
+#   include "setRootCase.H"
+#   include "createTime.H"
+#   include "createMesh.H"
+#   include "createFields.H"
+#   include "createSolidInterfaceNonLin.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-  Info<< "\nStarting time loop\n" << endl;
+    Info<< "\nStarting time loop\n" << endl;
 
-  while(runTime.loop())
+    while(runTime.loop())
     {
-      Info<< "Time: " << runTime.timeName() << nl << endl;
+        Info<< "Time: " << runTime.timeName() << nl << endl;
 
-#     include "readSolidMechanicsControls.H"
+#       include "readSolidMechanicsControls.H"
 
       int iCorr = 0;
       scalar initialResidual = 0;
@@ -70,36 +70,40 @@ int main(int argc, char *argv[])
       {
           U.storePrevIter();
 
-          surfaceTensorField shearGradU =
-              ((I - n*n)&fvc::interpolate(gradU));
+          surfaceTensorField shearGradU
+          (
+              "shearGradU",
+              (I - sqr(n)) & fvc::interpolate(gradU)
+          );
 
           fvVectorMatrix UEqn
+          (
+              fvm::d2dt2(rho, U)
+           ==
+              fvm::laplacian(2*muf + lambdaf, U, "laplacian(DU,U)")
+//             + fvc::div
+//               (
+//                   -(mu + lambda)*gradU
+//                   + mu*gradU.T()
+//                   + mu*(gradU & gradU.T())
+//                   + lambda*tr(gradU)*I
+//                   + 0.5*lambda*tr(gradU & gradU.T())*I
+//                   + (sigma & gradU),
+//                   "div(sigma)"
+//               )
+            + fvc::div
               (
-                  fvm::d2dt2(rho, U)
-                  ==
-                  fvm::laplacian(2*muf + lambdaf, U, "laplacian(DU,U)")
-                  // + fvc::div(
-                  //         -( (mu + lambda) * gradU )
-                  //         + ( mu * gradU.T() )
-                  //         + ( mu * (gradU & gradU.T()) )
-                  //         + ( lambda * tr(gradU) * I )
-                  //         + ( 0.5 * lambda * tr(gradU & gradU.T()) * I )
-                  //         + ( sigma & gradU ),
-                  //         "div(sigma)"
-                  //         )
-                  + fvc::div(
-                      mesh.magSf()
-                      *(
-                          - (muf + lambdaf)*(fvc::snGrad(U)&(I - n*n))
-                          + lambdaf*tr(shearGradU&(I - n*n))*n
-                          + muf*(shearGradU&n)
-                          + muf * (n & fvc::interpolate(gradU & gradU.T()))
-                          + 0.5*lambdaf
-                          *(n * tr(fvc::interpolate(gradU & gradU.T())))
-                          + (n & fvc::interpolate( sigma & gradU ))
-                          )
-                      )
-                  );
+                  mesh.magSf()*
+                  (
+                      - (muf + lambdaf)*(fvc::snGrad(U) & (I - n*n))
+                      + lambdaf*tr(shearGradU & (I - n*n))*n
+                      + muf*(shearGradU & n)
+                      + muf*(n & fvc::interpolate(gradU & gradU.T()))
+                      + 0.5*lambdaf*(n*tr(fvc::interpolate(gradU & gradU.T())))
+                      + (n & fvc::interpolate(sigma & gradU))
+                  )
+              )
+          );
 
           if (solidInterfaceCorr)
           {
@@ -153,9 +157,9 @@ int main(int argc, char *argv[])
           << " s\n\n" << endl;
     }
 
-  Info<< "End\n" << endl;
+    Info<< "End\n" << endl;
 
-  return(0);
+    return(0);
 }
 
 
diff --git a/applications/solvers/solidMechanics/elasticNonLinULSolidFoam/moveMesh.H b/applications/solvers/solidMechanics/elasticNonLinULSolidFoam/moveMesh.H
index 807c713fd..fe91884e1 100644
--- a/applications/solvers/solidMechanics/elasticNonLinULSolidFoam/moveMesh.H
+++ b/applications/solvers/solidMechanics/elasticNonLinULSolidFoam/moveMesh.H
@@ -1,15 +1,16 @@
 if(moveMeshMethod == "inverseDistance")
-  {
+{
 #   include "moveMeshInverseDistance.H"
-  }
- else if(moveMeshMethod == "leastSquares")
-   {
+}
+else if(moveMeshMethod == "leastSquares")
+{
 #   include "moveMeshLeastSquares.H"
-   }
- else
-   {
-     FatalError << "move mesh method " << moveMeshMethod << " not recognised" << nl
+}
+else
+{
+    FatalErrorIn(args.executable())
+        << "move mesh method " << moveMeshMethod << " not recognised" << nl
 		<< "available methods are:" << nl
 		<< "inverseDistance" << nl
 		<< "leastSquares" << exit(FatalError);
-   }
+}
diff --git a/applications/solvers/solidMechanics/elasticSolidFoam/aitkenRelaxation.H b/applications/solvers/solidMechanics/elasticSolidFoam/aitkenRelaxation.H
index e82330f12..779839046 100644
--- a/applications/solvers/solidMechanics/elasticSolidFoam/aitkenRelaxation.H
+++ b/applications/solvers/solidMechanics/elasticSolidFoam/aitkenRelaxation.H
@@ -1,4 +1,4 @@
-// aitken acceleration 
+// aitken acceleration
 
 aitkenDelta.storePrevIter();
 
@@ -7,26 +7,26 @@ aitkenDelta = (U - U.prevIter()) / aitkenInitialRes;
 
 // update relaxation factor
 if(iCorr == 0)
-  {
+{
     aitkenTheta = 0.1;
-  }
- else
-   {
-     vectorField b = aitkenDelta.internalField() - aitkenDelta.prevIter().internalField();
-     // scalar sumMagB = gSum(mag(b));
-     scalar sumMagB = gSum(magSqr(b));
-     if(sumMagB < SMALL)
-       {
-	 // Warning << "Aitken under-relaxation: denominator less then SMALL"
-	 // << endl;
-	 sumMagB += SMALL;
-       }
-     
-     aitkenTheta = -aitkenTheta*
-       gSum(aitkenDelta.prevIter().internalField() & b)
-       /
-       sumMagB;
-   }
+}
+else
+{
+    vectorField b = aitkenDelta.internalField()
+        - aitkenDelta.prevIter().internalField();
+
+    scalar sumMagB = gSum(magSqr(b));
+    if(sumMagB < SMALL)
+    {
+        // Warning << "Aitken under-relaxation: denominator less then SMALL"
+        // << endl;
+        sumMagB += SMALL;
+    }
+
+    aitkenTheta = -aitkenTheta*
+        gSum(aitkenDelta.prevIter().internalField() & b)/
+        sumMagB;
+}
 
 // correction to the latest U
 U += aitkenTheta*aitkenDelta*aitkenInitialRes;
diff --git a/applications/solvers/solidMechanics/elasticSolidFoam/createFields.H b/applications/solvers/solidMechanics/elasticSolidFoam/createFields.H
index 9e916fb44..c26827227 100644
--- a/applications/solvers/solidMechanics/elasticSolidFoam/createFields.H
+++ b/applications/solvers/solidMechanics/elasticSolidFoam/createFields.H
@@ -12,21 +12,20 @@
         mesh
     );
 
-    //     volTensorField gradU = fvc::grad(U);
     volTensorField gradU
     (
         IOobject
         (
-	 "grad(U)",
-	 runTime.timeName(),
-	 mesh,
-	 IOobject::NO_READ,
-	 IOobject::NO_WRITE
-	 ),
+            "grad(U)",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
         mesh,
         dimensionedTensor("zero", dimless, tensor::zero)
-     );
-    //surfaceVectorField snGradU = fvc::snGrad(U);
+    );
+
     surfaceVectorField snGradU
     (
         IOobject
@@ -38,7 +37,7 @@
             IOobject::NO_WRITE
         ),
         mesh,
-	dimensionedVector("zero", dimless, vector::zero)
+        dimensionedVector("zero", dimless, vector::zero)
     );
 
     volVectorField V
@@ -85,19 +84,19 @@
         dimensionedSymmTensor("zero", dimForce/dimArea, symmTensor::zero)
     );
 
-      volVectorField divSigmaExp
-      (
+    volVectorField divSigmaExp
+    (
         IOobject
         (
-	 "divSigmaExp",
-	 runTime.timeName(),
-	 mesh,
-	 IOobject::NO_READ,
-	 IOobject::NO_WRITE
-	 ),
+            "divSigmaExp",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
         mesh,
         dimensionedVector("zero", dimForce/dimVolume, vector::zero)
-       );
+    );
 
     constitutiveModel rheology(sigma, U);
 
@@ -127,18 +126,20 @@
 
    // for aitken relaxation
    volVectorField aitkenDelta
-    (
-        IOobject
-        (
-            "aitkenDelta",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh,
-	dimensionedVector("zero", dimLength, vector::zero)
-    );
+   (
+       IOobject
+       (
+           "aitkenDelta",
+           runTime.timeName(),
+           mesh,
+           IOobject::NO_READ,
+           IOobject::NO_WRITE
+       ),
+       mesh,
+       dimensionedVector("zero", dimLength, vector::zero)
+   );
+
+
 // aitken relaxation factor
 scalar aitkenInitialRes = 1.0;
 scalar aitkenTheta = 0.1;
diff --git a/applications/solvers/solidMechanics/elasticSolidFoam/elasticSolidFoam.C b/applications/solvers/solidMechanics/elasticSolidFoam/elasticSolidFoam.C
index d5977e294..905dc47b0 100644
--- a/applications/solvers/solidMechanics/elasticSolidFoam/elasticSolidFoam.C
+++ b/applications/solvers/solidMechanics/elasticSolidFoam/elasticSolidFoam.C
@@ -62,7 +62,7 @@ int main(int argc, char *argv[])
   Info<< "\nStarting time loop\n" << endl;
 
   while(runTime.loop())
-    {
+  {
       Info<< "Time: " << runTime.timeName() << nl << endl;
 
 #     include "readSolidMechanicsControls.H"
@@ -74,84 +74,80 @@ int main(int argc, char *argv[])
       lduMatrix::debug = 0;
 
       if (predictor)
-        {
-            Info << "\nPredicting U, gradU and snGradU based on V,"
-                 << "gradV and snGradV\n" << endl;
-            U += V*runTime.deltaT();
-            gradU += gradV*runTime.deltaT();
-            snGradU += snGradV*runTime.deltaT();
-        }
+      {
+          Info<< "\nPredicting U, gradU and snGradU based on V,"
+              << "gradV and snGradV\n" << endl;
+          U += V*runTime.deltaT();
+          gradU += gradV*runTime.deltaT();
+          snGradU += snGradV*runTime.deltaT();
+      }
 
       do
-        {
-            U.storePrevIter();
+      {
+          U.storePrevIter();
 
 #         include "calculateDivSigmaExp.H"
 
-            // linear momentum equation
-            fvVectorMatrix UEqn
-                (
-                    rho*fvm::d2dt2(U)
-                    ==
-                    fvm::laplacian(2*muf + lambdaf, U, "laplacian(DU,U)")
-                    + divSigmaExp
-                    );
+          // linear momentum equation
+          fvVectorMatrix UEqn
+          (
+              rho*fvm::d2dt2(U)
+           ==
+              fvm::laplacian(2*muf + lambdaf, U, "laplacian(DU,U)")
+            + divSigmaExp
+          );
 
-            if (solidInterfaceCorr)
-            {
-                solidInterfacePtr->correct(UEqn);
-            }
+          if (solidInterfaceCorr)
+          {
+              solidInterfacePtr->correct(UEqn);
+          }
 
-//        if (relaxEqn)
-//          {
-//            UEqn.relax();
-//          }
+          solverPerf = UEqn.solve();
 
-            solverPerf = UEqn.solve();
+          if (iCorr == 0)
+          {
+              initialResidual = solverPerf.initialResidual();
+              aitkenInitialRes = gMax(mag(U.internalField()));
+          }
 
-            if (iCorr == 0)
-            {
-                initialResidual = solverPerf.initialResidual();
-                aitkenInitialRes = gMax(mag(U.internalField()));
-            }
-
-            if (aitkenRelax)
-            {
+          if (aitkenRelax)
+          {
 #             include "aitkenRelaxation.H"
-            }
-            else
-            {
-                U.relax();
-            }
+          }
+          else
+          {
+              U.relax();
+          }
 
-            gradU = fvc::grad(U);
+          gradU = fvc::grad(U);
 
 #         include "calculateRelativeResidual.H"
 
-            if (iCorr % infoFrequency == 0)
-            {
-                Info<< "\tTime " << runTime.value()
-                    << ", Corrector " << iCorr
-                    << ", Solving for " << U.name()
-                    << " using " << solverPerf.solverName()
-                    << ", res = " << solverPerf.initialResidual()
-                    << ", rel res = " << relativeResidual;
-                if (aitkenRelax)
-                {
-                    Info<< ", aitken = " << aitkenTheta;
-                }
-                Info<< ", inner iters = " << solverPerf.nIterations() << endl;
-            }
-        }
+          if (iCorr % infoFrequency == 0)
+          {
+              Info<< "\tTime " << runTime.value()
+                  << ", Corrector " << iCorr
+                  << ", Solving for " << U.name()
+                  << " using " << solverPerf.solverName()
+                  << ", res = " << solverPerf.initialResidual()
+                  << ", rel res = " << relativeResidual;
+
+              if (aitkenRelax)
+              {
+                  Info<< ", aitken = " << aitkenTheta;
+              }
+              Info<< ", inner iters = " << solverPerf.nIterations() << endl;
+          }
+      }
       while
-          (
-              iCorr++ == 0
-              ||
-              (solverPerf.initialResidual() > convergenceTolerance
-               //relativeResidual > convergenceTolerance
-               &&
-               iCorr < nCorr)
-              );
+      (
+          iCorr++ == 0
+       || (
+              solverPerf.initialResidual() > convergenceTolerance
+              //relativeResidual > convergenceTolerance
+           && iCorr < nCorr
+          )
+      );
 
       Info<< nl << "Time " << runTime.value() << ", Solving for " << U.name()
           << ", Initial residual = " << initialResidual
@@ -178,9 +174,9 @@ int main(int argc, char *argv[])
           << " s\n\n" << endl;
     }
 
-  Info<< "End\n" << endl;
+    Info<< "End\n" << endl;
 
-  return(0);
+    return(0);
 }