diff --git a/applications/solvers/coupled/conjugateHeatTransfer/fvPatchFields/externalRadiation/ersPointSource/ersPointSource.C b/applications/solvers/coupled/conjugateHeatTransfer/fvPatchFields/externalRadiation/ersPointSource/ersPointSource.C
index f8b206c57..b66a430df 100644
--- a/applications/solvers/coupled/conjugateHeatTransfer/fvPatchFields/externalRadiation/ersPointSource/ersPointSource.C
+++ b/applications/solvers/coupled/conjugateHeatTransfer/fvPatchFields/externalRadiation/ersPointSource/ersPointSource.C
@@ -46,7 +46,7 @@ ersPointSource::ersPointSource
     alpha_(readScalar(dict.lookup("alpha"))),
     direction_(dict.lookup("direction"))
 {
-    q() = -alpha_*qmax_*min(direction_ & p.Sf()/p.magSf(), 0.0);
+    q() = -alpha_*qmax_*min(direction_ & p.Sf()/p.magSf(), scalar(0));
 }
 
 
diff --git a/applications/solvers/solidMechanics/deprecatedSolvers/contactStressFoam/contactPatchPair.C b/applications/solvers/solidMechanics/deprecatedSolvers/contactStressFoam/contactPatchPair.C
index 4804383a1..2af8a4675 100644
--- a/applications/solvers/solidMechanics/deprecatedSolvers/contactStressFoam/contactPatchPair.C
+++ b/applications/solvers/solidMechanics/deprecatedSolvers/contactStressFoam/contactPatchPair.C
@@ -30,18 +30,10 @@ Description
 #include "contactPatchPair.H"
 #include "volFields.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct from components
-contactPatchPair::contactPatchPair
+    Foam::contactPatchPair::contactPatchPair
 (
     const fvMesh& m,
     const label master,
@@ -75,6 +67,4 @@ contactPatchPair::contactPatchPair
 {}
 
 
-} // End namespace Foam
-
 // ************************************************************************* //
diff --git a/applications/solvers/solidMechanics/deprecatedSolvers/newContactStressFoam/contactPatchPair.C b/applications/solvers/solidMechanics/deprecatedSolvers/newContactStressFoam/contactPatchPair.C
index 9e5cf7992..59a70f0e4 100644
--- a/applications/solvers/solidMechanics/deprecatedSolvers/newContactStressFoam/contactPatchPair.C
+++ b/applications/solvers/solidMechanics/deprecatedSolvers/newContactStressFoam/contactPatchPair.C
@@ -302,7 +302,7 @@ void Foam::contactPatchPair::correct
             (
                 slavePressure
             ),
-            0.0
+            scalar(0)
         );
 
     // Calculate master traction, using the positive part of
diff --git a/applications/solvers/solidMechanics/elasticAcpSolidFoam/updateCrack.H b/applications/solvers/solidMechanics/elasticAcpSolidFoam/updateCrack.H
index 50bb49ddd..a9f4c0b83 100644
--- a/applications/solvers/solidMechanics/elasticAcpSolidFoam/updateCrack.H
+++ b/applications/solvers/solidMechanics/elasticAcpSolidFoam/updateCrack.H
@@ -6,13 +6,15 @@ nCoupledFacesToBreak = 0;
     // scalarField effTraction =
     //   cohesiveZone.internalField() *
     //   mag(traction.internalField());
-  scalarField normalTraction =
-    cohesiveZone.internalField() *
-    ( n.internalField() & traction.internalField() );
-  normalTraction = max(normalTraction, 0.0); // only consider tensile tractions
-  scalarField shearTraction =
-    cohesiveZone.internalField() *
-    mag( (I - Foam::sqr(n.internalField())) & traction.internalField() );
+    scalarField normalTraction =
+        cohesiveZone.internalField()*
+        ( n.internalField() & traction.internalField() );
+
+    // only consider tensile tractions
+    normalTraction = max(normalTraction, scalar(0));
+    scalarField shearTraction =
+        cohesiveZone.internalField() *
+        mag( (I - Foam::sqr(n.internalField())) & traction.internalField() );
 
     // the traction fraction is monitored to decide which faces to break:
     // ie (tN/tNC)^2 + (tS/tSC)^2 >1 to crack a face
@@ -27,16 +29,18 @@ nCoupledFacesToBreak = 0;
     scalarField effTractionFraction(normalTraction.size(), 0.0);
 
     if(cohesivePatchUPtr)
-      {
-	effTractionFraction =
-	  (normalTraction/sigmaMaxI)*(normalTraction/sigmaMaxI) + (shearTraction/tauMaxI)*(shearTraction/tauMaxI);
-      }
+    {
+        effTractionFraction =
+            (normalTraction/sigmaMaxI)*(normalTraction/sigmaMaxI)
+          + (shearTraction/tauMaxI)*(shearTraction/tauMaxI);
+    }
     else
-      {
-	// solidCohesiveFixedModeMix only uses sigmaMax
-	effTractionFraction =
-	  (normalTraction/sigmaMaxI)*(normalTraction/sigmaMaxI) + (shearTraction/sigmaMaxI)*(shearTraction/sigmaMaxI);
-      }
+    {
+        // solidCohesiveFixedModeMix only uses sigmaMax
+        effTractionFraction =
+            (normalTraction/sigmaMaxI)*(normalTraction/sigmaMaxI)
+          + (shearTraction/sigmaMaxI)*(shearTraction/sigmaMaxI);
+    }
     maxEffTractionFraction = gMax(effTractionFraction);
 
     SLList<label> facesToBreakList;
@@ -44,15 +48,21 @@ nCoupledFacesToBreak = 0;
 
     forAll(effTractionFraction, faceI)
     {
-      if (effTractionFraction[faceI] > 1.0)
+        if (effTractionFraction[faceI] > 1.0)
         {
             facesToBreakList.insert(faceI);
-            facesToBreakEffTractionFractionList.insert(effTractionFraction[faceI]);
+            facesToBreakEffTractionFractionList.insert
+            (
+                effTractionFraction[faceI]
+            );
         }
     }
 
     labelList facesToBreak(facesToBreakList);
-    List<scalar> facesToBreakEffTractionFraction(facesToBreakEffTractionFractionList);
+    List<scalar> facesToBreakEffTractionFraction
+    (
+        facesToBreakEffTractionFractionList
+    );
 
     nFacesToBreak = facesToBreak.size();
 
@@ -122,7 +132,7 @@ nCoupledFacesToBreak = 0;
 	  scalarField pNormalTraction =
 	    cohesiveZone.boundaryField()[patchI] *
 	    ( n.boundaryField()[patchI] & traction.boundaryField()[patchI] );
-	  pNormalTraction = max(pNormalTraction, 0.0); // only consider tensile tractions
+	  pNormalTraction = max(pNormalTraction, scalar(0)); // only consider tensile tractions
 	  scalarField pShearTraction =
 	    cohesiveZone.boundaryField()[patchI] *
 	    mag( (I - Foam::sqr(n.boundaryField()[patchI])) & traction.boundaryField()[patchI] );
diff --git a/applications/solvers/solidMechanics/elasticIncrAcpSolidFoam/updateCrack.H b/applications/solvers/solidMechanics/elasticIncrAcpSolidFoam/updateCrack.H
index c105bdedd..649ba6e4c 100644
--- a/applications/solvers/solidMechanics/elasticIncrAcpSolidFoam/updateCrack.H
+++ b/applications/solvers/solidMechanics/elasticIncrAcpSolidFoam/updateCrack.H
@@ -6,13 +6,16 @@ nCoupledFacesToBreak = 0;
     // scalarField effTraction =
     //   cohesiveZone.internalField() *
     //   mag(traction.internalField());
-  scalarField normalTraction =
-    cohesiveZone.internalField() *
-    ( n.internalField() & traction.internalField() );
-  normalTraction = max(normalTraction, 0.0); // only consider tensile tractions
-  scalarField shearTraction =
-    cohesiveZone.internalField() *
-    mag( (I - Foam::sqr(n.internalField())) & traction.internalField() );
+    scalarField normalTraction =
+        cohesiveZone.internalField()*
+        ( n.internalField() & traction.internalField() );
+
+     // only consider tensile tractions
+    normalTraction = max(normalTraction, scalar(0));
+    
+    scalarField shearTraction =
+        cohesiveZone.internalField() *
+        mag( (I - Foam::sqr(n.internalField())) & traction.internalField() );
 
     // the traction fraction is monitored to decide which faces to break:
     // ie (tN/tNC)^2 + (tS/tSC)^2 >1 to crack a face
@@ -28,16 +31,18 @@ nCoupledFacesToBreak = 0;
     scalarField effTractionFraction(normalTraction.size(), 0.0);
 
     if(cohesivePatchDUPtr)
-      {
-	effTractionFraction =
-	  (normalTraction/sigmaMaxI)*(normalTraction/sigmaMaxI) + (shearTraction/tauMaxI)*(shearTraction/tauMaxI);
-      }
+    {
+        effTractionFraction =
+            (normalTraction/sigmaMaxI)*(normalTraction/sigmaMaxI)
+          + (shearTraction/tauMaxI)*(shearTraction/tauMaxI);
+    }
     else
-      {
-	// solidCohesiveFixedModeMix only uses sigmaMax
-	effTractionFraction =
-	  (normalTraction/sigmaMaxI)*(normalTraction/sigmaMaxI) + (shearTraction/sigmaMaxI)*(shearTraction/sigmaMaxI);
-      }
+    {
+        // solidCohesiveFixedModeMix only uses sigmaMax
+        effTractionFraction =
+            (normalTraction/sigmaMaxI)*(normalTraction/sigmaMaxI)
+          + (shearTraction/sigmaMaxI)*(shearTraction/sigmaMaxI);
+    }
 
     maxEffTractionFraction = gMax(effTractionFraction);
 
@@ -46,15 +51,21 @@ nCoupledFacesToBreak = 0;
 
     forAll(effTractionFraction, faceI)
     {
-      if (effTractionFraction[faceI] > 1.0)
+        if (effTractionFraction[faceI] > 1.0)
         {
             facesToBreakList.insert(faceI);
-            facesToBreakEffTractionFractionList.insert(effTractionFraction[faceI]);
+            facesToBreakEffTractionFractionList.insert
+            (
+                effTractionFraction[faceI]
+            );
         }
     }
 
     labelList facesToBreak(facesToBreakList);
-    List<scalar> facesToBreakEffTractionFraction(facesToBreakEffTractionFractionList);
+    List<scalar> facesToBreakEffTractionFraction
+    (
+        facesToBreakEffTractionFractionList
+    );
 
     nFacesToBreak = facesToBreak.size();
 
@@ -69,9 +80,14 @@ nCoupledFacesToBreak = 0;
     scalar faceToBreakEffTractionFraction = 0;
     forAll(facesToBreakEffTractionFraction, faceI)
     {
-        if (facesToBreakEffTractionFraction[faceI] > faceToBreakEffTractionFraction)
+        if
+        (
+            facesToBreakEffTractionFraction[faceI]
+          > faceToBreakEffTractionFraction
+        )
         {
-            faceToBreakEffTractionFraction = facesToBreakEffTractionFraction[faceI];
+            faceToBreakEffTractionFraction =
+                facesToBreakEffTractionFraction[faceI];
             faceToBreakIndex = facesToBreak[faceI];
         }
     }
@@ -83,7 +99,11 @@ nCoupledFacesToBreak = 0;
         bool procHasFaceToBreak = false;
         if (nFacesToBreak > 0)
         {
-            if ( mag(gMaxEffTractionFraction - faceToBreakEffTractionFraction) < SMALL )
+            if
+            (
+                mag(gMaxEffTractionFraction - faceToBreakEffTractionFraction)
+              < SMALL
+            )
             {
                 // philipc - Maximum traction fraction is on this processor
                 procHasFaceToBreak = true;
@@ -123,9 +143,12 @@ nCoupledFacesToBreak = 0;
 	    // scalarField pEffTractionFraction = pEffTraction/sigmaMax.boundaryField()[patchI];
 
 	  scalarField pNormalTraction =
-	    cohesiveZone.boundaryField()[patchI] *
-	    ( n.boundaryField()[patchI] & traction.boundaryField()[patchI] );
-	  pNormalTraction = max(pNormalTraction, 0.0); // only consider tensile tractions
+          cohesiveZone.boundaryField()[patchI]*
+          ( n.boundaryField()[patchI] & traction.boundaryField()[patchI] );
+
+      // only consider tensile tractions
+	  pNormalTraction = max(pNormalTraction, scalar(0));
+
 	  scalarField pShearTraction =
 	    cohesiveZone.boundaryField()[patchI] *
 	    mag( (I - Foam::sqr(n.boundaryField()[patchI])) & traction.boundaryField()[patchI] );
@@ -139,24 +162,24 @@ nCoupledFacesToBreak = 0;
 	  //   (pNormalTraction/pSigmaMax)*(pNormalTraction/pSigmaMax) + (pShearTraction/pTauMax)*(pShearTraction/pTauMax);
 	  scalarField pEffTractionFraction(pNormalTraction.size(), 0.0);
 	  if(cohesivePatchDUPtr)
-	    {
+      {
 	      pEffTractionFraction =
 		(pNormalTraction/pSigmaMax)*(pNormalTraction/pSigmaMax) + (pShearTraction/pTauMax)*(pShearTraction/pTauMax);
-	    }
+      }
 	  else
-	    {
+      {
 	      // solidCohesiveFixedModeMix only uses sigmaMax
 	      pEffTractionFraction =
 		(pNormalTraction/pSigmaMax)*(pNormalTraction/pSigmaMax) + (pShearTraction/pSigmaMax)*(pShearTraction/pSigmaMax);
-	    }
+      }
 
             label start = mesh.boundaryMesh()[patchI].start();
 
             forAll(pEffTractionFraction, faceI)
             {
-	      if (pEffTractionFraction[faceI] > maxEffTractionFraction)
-		{
-		  maxEffTractionFraction = pEffTractionFraction[faceI];
+                if (pEffTractionFraction[faceI] > maxEffTractionFraction)
+                {
+                    maxEffTractionFraction = pEffTractionFraction[faceI];
                 }
 
 	      if (pEffTractionFraction[faceI] > 1.0)
@@ -368,30 +391,41 @@ nCoupledFacesToBreak = 0;
 	faceToBreakSigmaMax = sigmaMax.boundaryField()[patchID][localIndex];
 	faceToBreakTauMax = tauMax.boundaryField()[patchID][localIndex];
 	scalar normalTrac = faceToBreakNormal & faceToBreakTraction;
-	normalTrac = max(normalTrac, 0.0);
+	normalTrac = max(normalTrac, scalar(0));
 	scalar shearTrac = mag( (I - sqr(faceToBreakNormal)) & faceToBreakTraction );
 	scalar scaleFactor = 1;
 	if(cohesivePatchDUPtr)
-	  {
+    {
 	    scaleFactor =
-	      ::sqrt(1 / (
-			  (normalTrac/faceToBreakSigmaMax)*(normalTrac/faceToBreakSigmaMax)
-			  + (shearTrac/faceToBreakTauMax)*(shearTrac/faceToBreakTauMax)
-			  ) );
-	  }
+            Foam::sqrt
+            (
+                1/
+                (
+                    (normalTrac/faceToBreakSigmaMax)*
+                    (normalTrac/faceToBreakSigmaMax)
+                  + (shearTrac/faceToBreakTauMax)*(shearTrac/faceToBreakTauMax)
+                )
+            );
+    }
 	else
-	  {
+    {
 	      // solidCohesiveFixedModeMix only uses sigmaMax
 	    scaleFactor =
-	      ::sqrt(1 / (
-			  (normalTrac/faceToBreakSigmaMax)*(normalTrac/faceToBreakSigmaMax)
-			  + (shearTrac/faceToBreakSigmaMax)*(shearTrac/faceToBreakSigmaMax)
-			  ) );
-	  }
+            Foam::sqrt
+            (
+                1/
+                (
+                    (normalTrac/faceToBreakSigmaMax)*
+                    (normalTrac/faceToBreakSigmaMax)
+                  + (shearTrac/faceToBreakSigmaMax)*
+                    (shearTrac/faceToBreakSigmaMax)
+                )
+            );
+    }
 
-        faceToBreakTraction *= scaleFactor;
-
-        topoChange = true;
+    faceToBreakTraction *= scaleFactor;
+    
+    topoChange = true;
     }
 
     reduce(topoChange, orOp<bool>());
diff --git a/applications/solvers/solidMechanics/elasticOrthoAcpSolidFoam/updateCrack.H b/applications/solvers/solidMechanics/elasticOrthoAcpSolidFoam/updateCrack.H
index 40b129eec..8d6a68324 100644
--- a/applications/solvers/solidMechanics/elasticOrthoAcpSolidFoam/updateCrack.H
+++ b/applications/solvers/solidMechanics/elasticOrthoAcpSolidFoam/updateCrack.H
@@ -6,13 +6,15 @@ nCoupledFacesToBreak = 0;
     // scalarField effTraction =
     //   cohesiveZone.internalField() *
     //   mag(traction.internalField());
-  scalarField normalTraction =
-    cohesiveZone.internalField() *
-    ( n.internalField() & traction.internalField() );
-  normalTraction = max(normalTraction, 0.0); // only consider tensile tractions
-  scalarField shearTraction =
-    cohesiveZone.internalField() *
-    mag( (I - Foam::sqr(n.internalField())) & traction.internalField() );
+    scalarField normalTraction =
+        cohesiveZone.internalField()*
+        ( n.internalField() & traction.internalField() );
+
+    // only consider tensile tractions
+    normalTraction = max(normalTraction, scalar(0));
+    scalarField shearTraction =
+        cohesiveZone.internalField()*
+        mag( (I - Foam::sqr(n.internalField())) & traction.internalField() );
 
     // the traction fraction is monitored to decide which faces to break:
     // ie (tN/tNC)^2 + (tS/tSC)^2 >1 to crack a face
@@ -26,17 +28,19 @@ nCoupledFacesToBreak = 0;
     //scalarField effTractionFraction = effTraction/sigmaMax;
     scalarField effTractionFraction(normalTraction.size(), 0.0);
 
-    if(cohesivePatchUPtr)
-      {
-	effTractionFraction =
-	  (normalTraction/sigmaMaxI)*(normalTraction/sigmaMaxI) + (shearTraction/tauMaxI)*(shearTraction/tauMaxI);
-      }
+    if (cohesivePatchUPtr)
+    {
+        effTractionFraction =
+            (normalTraction/sigmaMaxI)*(normalTraction/sigmaMaxI)
+          + (shearTraction/tauMaxI)*(shearTraction/tauMaxI);
+    }
     else
-      {
-	// solidCohesiveFixedModeMix only uses sigmaMax
-	effTractionFraction =
-	  (normalTraction/sigmaMaxI)*(normalTraction/sigmaMaxI) + (shearTraction/sigmaMaxI)*(shearTraction/sigmaMaxI);
-      }
+    {
+        // solidCohesiveFixedModeMix only uses sigmaMax
+        effTractionFraction =
+            (normalTraction/sigmaMaxI)*(normalTraction/sigmaMaxI)
+          + (shearTraction/sigmaMaxI)*(shearTraction/sigmaMaxI);
+    }
     maxEffTractionFraction = gMax(effTractionFraction);
 
     SLList<label> facesToBreakList;
@@ -44,15 +48,21 @@ nCoupledFacesToBreak = 0;
 
     forAll(effTractionFraction, faceI)
     {
-      if (effTractionFraction[faceI] > 1.0)
+        if (effTractionFraction[faceI] > 1.0)
         {
             facesToBreakList.insert(faceI);
-            facesToBreakEffTractionFractionList.insert(effTractionFraction[faceI]);
+            facesToBreakEffTractionFractionList.insert
+            (
+                effTractionFraction[faceI]
+            );
         }
     }
 
     labelList facesToBreak(facesToBreakList);
-    List<scalar> facesToBreakEffTractionFraction(facesToBreakEffTractionFractionList);
+    List<scalar> facesToBreakEffTractionFraction
+    (
+        facesToBreakEffTractionFractionList
+    );
 
     nFacesToBreak = facesToBreak.size();
 
@@ -67,9 +77,15 @@ nCoupledFacesToBreak = 0;
     scalar faceToBreakEffTractionFraction = 0;
     forAll(facesToBreakEffTractionFraction, faceI)
     {
-        if (facesToBreakEffTractionFraction[faceI] > faceToBreakEffTractionFraction)
+        if
+        (
+            facesToBreakEffTractionFraction[faceI]
+          > faceToBreakEffTractionFraction
+        )
         {
-            faceToBreakEffTractionFraction = facesToBreakEffTractionFraction[faceI];
+            faceToBreakEffTractionFraction =
+                facesToBreakEffTractionFraction[faceI];
+
             faceToBreakIndex = facesToBreak[faceI];
         }
     }
@@ -82,7 +98,11 @@ nCoupledFacesToBreak = 0;
         bool procHasFaceToBreak = false;
         if (nFacesToBreak > 0)
         {
-            if ( mag(gMaxEffTractionFraction - faceToBreakEffTractionFraction) < SMALL )
+            if
+            (
+                mag(gMaxEffTractionFraction - faceToBreakEffTractionFraction)
+              < SMALL
+            )
             {
                 // philipc - Maximum traction fraction is on this processor
                 procHasFaceToBreak = true;
@@ -90,7 +110,7 @@ nCoupledFacesToBreak = 0;
         }
 
         // Check if maximum is present on more then one processors
-	label procID = Pstream::nProcs();
+        label procID = Pstream::nProcs();
         if (procHasFaceToBreak)
         {
             procID = Pstream::myProcNo();
@@ -115,54 +135,59 @@ nCoupledFacesToBreak = 0;
         if (mesh.boundary()[patchI].coupled())
         {
             // scalarField pEffTraction =
-	    //   cohesiveZone.boundaryField()[patchI] *
-	    //   mag(traction.boundaryField()[patchI]);
-	    // scalarField pEffTractionFraction = pEffTraction/sigmaMax.boundaryField()[patchI];
+        //   cohesiveZone.boundaryField()[patchI] *
+        //   mag(traction.boundaryField()[patchI]);
+        // scalarField pEffTractionFraction = pEffTraction/sigmaMax.boundaryField()[patchI];
 
-	  scalarField pNormalTraction =
-	    cohesiveZone.boundaryField()[patchI] *
-	    ( n.boundaryField()[patchI] & traction.boundaryField()[patchI] );
-	  pNormalTraction = max(pNormalTraction, 0.0); // only consider tensile tractions
-	  scalarField pShearTraction =
-	    cohesiveZone.boundaryField()[patchI] *
-	    mag( (I - Foam::sqr(n.boundaryField()[patchI])) & traction.boundaryField()[patchI] );
+      scalarField pNormalTraction =
+        cohesiveZone.boundaryField()[patchI] *
+        ( n.boundaryField()[patchI] & traction.boundaryField()[patchI] );
 
-	  // the traction fraction is monitored to decide which faces to break:
-	  // ie (tN/tNC)^2 + (tS/tSC)^2 >1 to crack a face
-	  const scalarField& pSigmaMax = sigmaMax.boundaryField()[patchI];
-	  const scalarField& pTauMax = tauMax.boundaryField()[patchI];
+      // only consider tensile tractions
+      pNormalTraction = max(pNormalTraction, scalar(0));
+      scalarField pShearTraction =
+        cohesiveZone.boundaryField()[patchI] *
+        mag( (I - Foam::sqr(n.boundaryField()[patchI])) & traction.boundaryField()[patchI] );
 
-	  scalarField pEffTractionFraction(pNormalTraction.size(), 0.0);
-	  if(cohesivePatchUPtr)
-	    {
-	      pEffTractionFraction =
-		(pNormalTraction/pSigmaMax)*(pNormalTraction/pSigmaMax) + (pShearTraction/pTauMax)*(pShearTraction/pTauMax);
-	    }
-	  else
-	    {
-	      // solidCohesiveFixedModeMix only uses sigmaMax
-	      pEffTractionFraction =
-		(pNormalTraction/pSigmaMax)*(pNormalTraction/pSigmaMax) + (pShearTraction/pSigmaMax)*(pShearTraction/pSigmaMax);
-	    }
+      // the traction fraction is monitored to decide which faces to break:
+      // ie (tN/tNC)^2 + (tS/tSC)^2 >1 to crack a face
+      const scalarField& pSigmaMax = sigmaMax.boundaryField()[patchI];
+      const scalarField& pTauMax = tauMax.boundaryField()[patchI];
 
-            label start = mesh.boundaryMesh()[patchI].start();
+      scalarField pEffTractionFraction(pNormalTraction.size(), 0);
 
-            forAll(pEffTractionFraction, faceI)
-            {
-	      if (pEffTractionFraction[faceI] > maxEffTractionFraction)
-		{
-		  maxEffTractionFraction = pEffTractionFraction[faceI];
-                }
+      if(cohesivePatchUPtr)
+      {
+          pEffTractionFraction =
+              (pNormalTraction/pSigmaMax)*(pNormalTraction/pSigmaMax)
+            + (pShearTraction/pTauMax)*(pShearTraction/pTauMax);
+      }
+      else
+      {
+          // solidCohesiveFixedModeMix only uses sigmaMax
+          pEffTractionFraction =
+              (pNormalTraction/pSigmaMax)*(pNormalTraction/pSigmaMax)
+            + (pShearTraction/pSigmaMax)*(pShearTraction/pSigmaMax);
+      }
 
-	      if (pEffTractionFraction[faceI] > 1.0)
-                {
-                    coupledFacesToBreakList.insert(start + faceI);
-                    coupledFacesToBreakEffTractionFractionList.insert
-                    (
-                        pEffTractionFraction[faceI]
-                    );
-                }
-            }
+      label start = mesh.boundaryMesh()[patchI].start();
+
+      forAll(pEffTractionFraction, faceI)
+      {
+          if (pEffTractionFraction[faceI] > maxEffTractionFraction)
+          {
+              maxEffTractionFraction = pEffTractionFraction[faceI];
+          }
+
+          if (pEffTractionFraction[faceI] > 1.0)
+          {
+              coupledFacesToBreakList.insert(start + faceI);
+              coupledFacesToBreakEffTractionFractionList.insert
+                  (
+                      pEffTractionFraction[faceI]
+                  );
+          }
+      }
         }
     }
 
@@ -250,7 +275,7 @@ nCoupledFacesToBreak = 0;
     if (nCoupledFacesToBreak)
     {
       label patchID =
-	mesh.boundaryMesh().whichPatch(coupledFaceToBreakIndex);
+    mesh.boundaryMesh().whichPatch(coupledFaceToBreakIndex);
 
         label start = mesh.boundaryMesh()[patchID].start();
         label localIndex = coupledFaceToBreakIndex - start;
@@ -315,31 +340,31 @@ nCoupledFacesToBreak = 0;
         faceToBreakNormal = n.internalField()[faceToBreakIndex];
 
         // Scale broken face traction
-	faceToBreakSigmaMax = sigmaMaxI[faceToBreakIndex];
-	faceToBreakTauMax = tauMaxI[faceToBreakIndex];
-	scalar normalTrac = faceToBreakNormal & faceToBreakTraction;
-	normalTrac = max(normalTrac, 0.0);
-	scalar shearTrac = mag( (I - sqr(faceToBreakNormal)) & faceToBreakTraction );
-	scalar scaleFactor = 1;
-	if(cohesivePatchUPtr)
-	  {
-	    scaleFactor =
-	      ::sqrt(1 / (
-			  (normalTrac/faceToBreakSigmaMax)*(normalTrac/faceToBreakSigmaMax)
-			  + (shearTrac/faceToBreakTauMax)*(shearTrac/faceToBreakTauMax)
-			  ) );
-	    }
-	  else
-	    {
-	      // solidCohesiveFixedModeMix only uses sigmaMax
-	    scaleFactor =
-	      ::sqrt(1 / (
-			  (normalTrac/faceToBreakSigmaMax)*(normalTrac/faceToBreakSigmaMax)
-			  + (shearTrac/faceToBreakSigmaMax)*(shearTrac/faceToBreakSigmaMax)
-			  ) );
-	    }
+    faceToBreakSigmaMax = sigmaMaxI[faceToBreakIndex];
+    faceToBreakTauMax = tauMaxI[faceToBreakIndex];
+    scalar normalTrac = faceToBreakNormal & faceToBreakTraction;
+    normalTrac = max(normalTrac, 0.0);
+    scalar shearTrac = mag( (I - sqr(faceToBreakNormal)) & faceToBreakTraction );
+    scalar scaleFactor = 1;
+    if(cohesivePatchUPtr)
+      {
+        scaleFactor =
+          ::sqrt(1 / (
+              (normalTrac/faceToBreakSigmaMax)*(normalTrac/faceToBreakSigmaMax)
+              + (shearTrac/faceToBreakTauMax)*(shearTrac/faceToBreakTauMax)
+              ) );
+        }
+      else
+        {
+          // solidCohesiveFixedModeMix only uses sigmaMax
+        scaleFactor =
+          ::sqrt(1 / (
+              (normalTrac/faceToBreakSigmaMax)*(normalTrac/faceToBreakSigmaMax)
+              + (shearTrac/faceToBreakSigmaMax)*(shearTrac/faceToBreakSigmaMax)
+              ) );
+        }
 
-	faceToBreakTraction *= scaleFactor;
+    faceToBreakTraction *= scaleFactor;
 
         topoChange = true;
     }
@@ -354,29 +379,29 @@ nCoupledFacesToBreak = 0;
         faceToBreakNormal = n.boundaryField()[patchID][localIndex];
 
         // Scale broken face traction
-	faceToBreakSigmaMax = sigmaMax.boundaryField()[patchID][localIndex];
-	faceToBreakTauMax = tauMax.boundaryField()[patchID][localIndex];
-	scalar normalTrac = faceToBreakNormal & faceToBreakTraction;
-	normalTrac = max(normalTrac, 0.0);
-	scalar shearTrac = mag( (I - sqr(faceToBreakNormal)) & faceToBreakTraction );
-	scalar scaleFactor = 1;
-	if(cohesivePatchUPtr)
-	  {
-	    scaleFactor =
-	      ::sqrt(1 / (
-			  (normalTrac/faceToBreakSigmaMax)*(normalTrac/faceToBreakSigmaMax)
-			  + (shearTrac/faceToBreakTauMax)*(shearTrac/faceToBreakTauMax)
-			  ) );
-	  }
-	else
-	  {
-	      // solidCohesiveFixedModeMix only uses sigmaMax
-	    scaleFactor =
-	      ::sqrt(1 / (
-			  (normalTrac/faceToBreakSigmaMax)*(normalTrac/faceToBreakSigmaMax)
-			  + (shearTrac/faceToBreakSigmaMax)*(shearTrac/faceToBreakSigmaMax)
-			  ) );
-	  }
+    faceToBreakSigmaMax = sigmaMax.boundaryField()[patchID][localIndex];
+    faceToBreakTauMax = tauMax.boundaryField()[patchID][localIndex];
+    scalar normalTrac = faceToBreakNormal & faceToBreakTraction;
+    normalTrac = max(normalTrac, 0.0);
+    scalar shearTrac = mag( (I - sqr(faceToBreakNormal)) & faceToBreakTraction );
+    scalar scaleFactor = 1;
+    if(cohesivePatchUPtr)
+      {
+        scaleFactor =
+          ::sqrt(1 / (
+              (normalTrac/faceToBreakSigmaMax)*(normalTrac/faceToBreakSigmaMax)
+              + (shearTrac/faceToBreakTauMax)*(shearTrac/faceToBreakTauMax)
+              ) );
+      }
+    else
+      {
+          // solidCohesiveFixedModeMix only uses sigmaMax
+        scaleFactor =
+          ::sqrt(1 / (
+              (normalTrac/faceToBreakSigmaMax)*(normalTrac/faceToBreakSigmaMax)
+              + (shearTrac/faceToBreakSigmaMax)*(shearTrac/faceToBreakSigmaMax)
+              ) );
+      }
 
         faceToBreakTraction *= scaleFactor;
 
@@ -416,20 +441,20 @@ nCoupledFacesToBreak = 0;
         Cf = fvc::interpolate(C);
         Kf = fvc::interpolate(K);
 
-	// we need to modify propertiess after cracking otherwise momentum equation is wrong
-	// but solidInterface seems to hold some information about old mesh
-	// so we will delete it and make another
-	// we could probably add a public clearout function
-	// create new solidInterface
-	//Pout << "Creating new solidInterface" << endl;
-	//delete solidInterfacePtr;
-	//solidInterfacePtr = new solidInterface(mesh, rheology);
-	// delete demand driven data as the mesh has changed
-	if(rheology.solidInterfaceActive())
-	  {
-	    rheology.solInterface().clearOut();
-	    solidInterfacePtr->modifyProperties(Cf, Kf);
-	  }
+    // we need to modify propertiess after cracking otherwise momentum equation is wrong
+    // but solidInterface seems to hold some information about old mesh
+    // so we will delete it and make another
+    // we could probably add a public clearout function
+    // create new solidInterface
+    //Pout << "Creating new solidInterface" << endl;
+    //delete solidInterfacePtr;
+    //solidInterfacePtr = new solidInterface(mesh, rheology);
+    // delete demand driven data as the mesh has changed
+    if(rheology.solidInterfaceActive())
+      {
+        rheology.solInterface().clearOut();
+        solidInterfacePtr->modifyProperties(Cf, Kf);
+      }
 
         // Local crack displacement
         vectorField UpI =
@@ -441,21 +466,21 @@ nCoupledFacesToBreak = 0;
         vectorField globalUpI = mesh.globalCrackField(UpI);
         vectorField globalOldUpI = mesh.globalCrackField(oldUpI);
 
-	// C and K field on new crack faces must be updated
+    // C and K field on new crack faces must be updated
         symmTensor4thOrderField CPI = C.boundaryField()[cohesivePatchID].patchInternalField();
         diagTensorField KPI = K.boundaryField()[cohesivePatchID].patchInternalField();
         symmTensor4thOrderField globalCPI = mesh.globalCrackField(CPI);
         diagTensorField globalKPI = mesh.globalCrackField(KPI);
 
-	// cohesivePatchU.size()
-	int cohesivePatchSize(cohesivePatchUPtr ? cohesivePatchUPtr->size() : cohesivePatchUFixedModePtr->size());
+    // cohesivePatchU.size()
+    int cohesivePatchSize(cohesivePatchUPtr ? cohesivePatchUPtr->size() : cohesivePatchUFixedModePtr->size());
 
         // Initialise U for new cohesive face
         const labelList& gcfa = mesh.globalCrackFaceAddressing();
         label globalIndex = mesh.localCrackStart();
         // for (label i=0; i<cohesivePatchU.size(); i++)
         for (label i=0; i<cohesivePatchSize; i++)
-	  {
+      {
             label oldFaceIndex = faceMap[start+i];
 
             // If new face
@@ -474,10 +499,10 @@ nCoupledFacesToBreak = 0;
                      + globalOldUpI[gcfa[globalIndex]]
                     );
 
-		// initialise C and K on new faces
-		// set new face value to value of internal cell
-		Cf.boundaryField()[cohesivePatchID][i] = globalCPI[globalIndex];
-		Kf.boundaryField()[cohesivePatchID][i] = globalKPI[globalIndex];
+        // initialise C and K on new faces
+        // set new face value to value of internal cell
+        Cf.boundaryField()[cohesivePatchID][i] = globalCPI[globalIndex];
+        Kf.boundaryField()[cohesivePatchID][i] = globalKPI[globalIndex];
 
                 globalIndex++;
             }
@@ -488,81 +513,87 @@ nCoupledFacesToBreak = 0;
         }
 
         // we must calculate grad using interface
-	// U at the interface has not been calculated yet as interface.correct()
-	// has not been called yet
-	// not really a problem as gradU is correct in second outer iteration
-	// as long as this does not cause convergence problems for the first iterations.
-	// we should be able to calculate the interface displacements without
-	// having to call interface.correct()
-	// todo: add calculateInterfaceU() function
-	// interface grad uses Gauss, we need least squares
-	//gradU = solidInterfacePtr->grad(U);
+    // U at the interface has not been calculated yet as interface.correct()
+    // has not been called yet
+    // not really a problem as gradU is correct in second outer iteration
+    // as long as this does not cause convergence problems for the first iterations.
+    // we should be able to calculate the interface displacements without
+    // having to call interface.correct()
+    // todo: add calculateInterfaceU() function
+    // interface grad uses Gauss, we need least squares
+    //gradU = solidInterfacePtr->grad(U);
         gradU = fvc::grad(U); // leastSquaresSolidInterface grad scheme
         //snGradU = fvc::snGrad(U);
 
 #       include "calculateTraction.H"
-	//if (nFacesToBreak || nCoupledFacesToBreak) mesh.write(); traction.write();
+    //if (nFacesToBreak || nCoupledFacesToBreak) mesh.write(); traction.write();
 
         // Initialise initiation traction for new cohesive patch face
         // for (label i=0; i<cohesivePatchU.size(); i++)
-	for (label i=0; i<cohesivePatchSize; i++)
+    for (label i=0; i<cohesivePatchSize; i++)
+    {
+        label oldFaceIndex = faceMap[start+i];
+
+        // If new face
+        if
+        (
+            (oldFaceIndex == faceToBreakIndex)
+         || (oldFaceIndex == coupledFaceToBreakIndex)
+        )
         {
-            label oldFaceIndex = faceMap[start+i];
+            vector n0 =
+                mesh.Sf().boundaryField()[cohesivePatchID][i]
+                /mesh.magSf().boundaryField()[cohesivePatchID][i];
+            //vector n1 = -n0;
 
-            // If new face
-            if
-            (
-                (oldFaceIndex == faceToBreakIndex)
-             || (oldFaceIndex == coupledFaceToBreakIndex)
-            )
+            if ((n0&faceToBreakNormal) > SMALL)
             {
-                vector n0 =
-                    mesh.Sf().boundaryField()[cohesivePatchID][i]
-                   /mesh.magSf().boundaryField()[cohesivePatchID][i];
-                //vector n1 = -n0;
-
-                if ((n0&faceToBreakNormal) > SMALL)
-                {
-		  traction.boundaryField()[cohesivePatchID][i] =
+                traction.boundaryField()[cohesivePatchID][i] =
                     faceToBreakTraction;
 
-		  traction.oldTime().boundaryField()[cohesivePatchID][i] =
-		    faceToBreakTraction;
+                traction.oldTime().boundaryField()[cohesivePatchID][i] =
+                    faceToBreakTraction;
 
-		  if(cohesivePatchUPtr)
-		    {
-		      cohesivePatchUPtr->traction()[i] = faceToBreakTraction;
-		    }
-		  else
-		    {
-		      cohesivePatchUFixedModePtr->traction()[i] = faceToBreakTraction;
-		      cohesivePatchUFixedModePtr->initiationTraction()[i] = faceToBreakTraction;
-		    }
+                if(cohesivePatchUPtr)
+                {
+                    cohesivePatchUPtr->traction()[i] = faceToBreakTraction;
                 }
                 else
                 {
-		  traction.boundaryField()[cohesivePatchID][i] =
+                    cohesivePatchUFixedModePtr->traction()[i] =
+                        faceToBreakTraction;
+
+                    cohesivePatchUFixedModePtr->initiationTraction()[i] =
+                        faceToBreakTraction;
+                }
+            }
+            else
+            {
+                traction.boundaryField()[cohesivePatchID][i] =
                     -faceToBreakTraction;
 
-		  traction.oldTime().boundaryField()[cohesivePatchID][i] =
-		    -faceToBreakTraction;
+                traction.oldTime().boundaryField()[cohesivePatchID][i] =
+                    -faceToBreakTraction;
 
-		  //cohesivePatchU.traction()[i] = -faceToBreakTraction;
-		  if(cohesivePatchUPtr)
-		    {
-		      cohesivePatchUPtr->traction()[i] = -faceToBreakTraction;
-		    }
-		  else
-		    {
-		      cohesivePatchUFixedModePtr->traction()[i] = -faceToBreakTraction;
-		      cohesivePatchUFixedModePtr->initiationTraction()[i] = -faceToBreakTraction;
-		    }
+                //cohesivePatchU.traction()[i] = -faceToBreakTraction;
+                if(cohesivePatchUPtr)
+                {
+                    cohesivePatchUPtr->traction()[i] = -faceToBreakTraction;
+                }
+                else
+                {
+                    cohesivePatchUFixedModePtr->traction()[i] =
+                        -faceToBreakTraction;
+
+                    cohesivePatchUFixedModePtr->initiationTraction()[i] =
+                        -faceToBreakTraction;
                 }
             }
         }
+    }
 
-	// hmmnn we only need a reference for very small groups of cells
-	// turn off for now
-	//#       include "updateReference.H"
+    // hmmnn we only need a reference for very small groups of cells
+    // turn off for now
+    //#       include "updateReference.H"
     }
 }
diff --git a/applications/utilities/miscellaneous/autoDerivative/testFadField/testFadField.C b/applications/utilities/miscellaneous/autoDerivative/testFadField/testFadField.C
index d07a6942f..547e0198d 100644
--- a/applications/utilities/miscellaneous/autoDerivative/testFadField/testFadField.C
+++ b/applications/utilities/miscellaneous/autoDerivative/testFadField/testFadField.C
@@ -43,7 +43,7 @@ int main(int argc, char *argv[])
     Info << "tf value: " << FadOneValue(tf) << endl;
     Info << "tf deriv: " << FadOneDeriv(tf, 0) << endl;
 
-    FadOneSetDeriv(tf, 0, Field<double>(3, 1));
+    FadOneSetDeriv(tf, 0, scalarField(3, 1));
 
     fadScalar ten(10.0);
     fadScalar bs = tf[0] - ten;
diff --git a/etc/bashrc b/etc/bashrc
index 479f0a88e..4a73e2504 100755
--- a/etc/bashrc
+++ b/etc/bashrc
@@ -138,7 +138,7 @@ export WM_COMPILER_LIB_ARCH=
 # WM_ARCH_OPTION = 32 | 64
 : ${WM_ARCH_OPTION:=64}; export WM_ARCH_OPTION
 
-# WM_PRECISION_OPTION = DP | SP
+# WM_PRECISION_OPTION = LDP | DP | SP
 : ${WM_PRECISION_OPTION:=DP}; export WM_PRECISION_OPTION
 
 # WM_COMPILE_OPTION = Opt | Debug | Prof
diff --git a/src/Pstream/mpi/Pstream.C b/src/Pstream/mpi/Pstream.C
index b7cddd6ec..b21f097c9 100644
--- a/src/Pstream/mpi/Pstream.C
+++ b/src/Pstream/mpi/Pstream.C
@@ -37,6 +37,8 @@ License
 #   define MPI_SCALAR MPI_FLOAT
 #elif defined(WM_DP)
 #   define MPI_SCALAR MPI_DOUBLE
+#elif defined(WM_LDP)
+#   define MPI_SCALAR MPI_LONG_DOUBLE
 #endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/cudaSolvers/cudaSolver/cudaTypes.H b/src/cudaSolvers/cudaSolver/cudaTypes.H
index 64b6a9e29..f50f06733 100644
--- a/src/cudaSolvers/cudaSolver/cudaTypes.H
+++ b/src/cudaSolvers/cudaSolver/cudaTypes.H
@@ -38,15 +38,19 @@ Author
 #ifndef cudaTypes_H
 #define cudaTypes_H
 
-#if defined(WM_DP)
+#if defined(WM_SP)
 
-    #define MPI_SCALAR MPI_DOUBLE
+#   define MPI_SCALAR MPI_FLOAT
+    typedef float ValueType;
+
+#elif defined(WM_DP)
+
+#   define MPI_SCALAR MPI_DOUBLE
     typedef double ValueType;
 
-#elif defined(WM_SP)
+#elif defined(WM_LDP)
 
-    #define MPI_SCALAR MPI_FLOAT
-    typedef float ValueType;
+#   error Long double not supported in CUDA
 
 #endif
 
diff --git a/src/dbns/multigrid/mgMeshLevel/mgMeshLevel.C b/src/dbns/multigrid/mgMeshLevel/mgMeshLevel.C
index 247a05154..9f94a70cd 100644
--- a/src/dbns/multigrid/mgMeshLevel/mgMeshLevel.C
+++ b/src/dbns/multigrid/mgMeshLevel/mgMeshLevel.C
@@ -163,6 +163,8 @@ void Foam::mgMeshLevel::makeChild() const
     child_.setSize(nCells());
     int nMoves = -1;
 
+#   ifdef WM_DP
+
     MGridGen
     (
         nCells(),
@@ -175,10 +177,51 @@ void Foam::mgMeshLevel::makeChild() const
         mgMaxClusterSize_,
         options.begin(),
         &nMoves,
+        &nCoarseCells,
+        child_.begin()
+    );
+
+#   else
+
+    // Conversion of type for MGridGen interface
+    const scalarField& vols = cellVolumes();
+
+    Field<double> dblVols(vols.size());
+    forAll (dblVols, i)
+    {
+        dblVols[i] = vols[i];
+    }
+
+    Field<double> dblAreas(boundaryAreas.size());
+    forAll (dblAreas, i)
+    {
+        dblAreas[i] = boundaryAreas[i];
+    }
+
+    Field<double> dblFaceWeights(faceWeights.size());
+    forAll (dblFaceWeights, i)
+    {
+        dblFaceWeights[i] = faceWeights[i];
+    }
+
+    MGridGen
+    (
+        nCells(),
+        cellCellOffsets.begin(),
+        dblVols.begin(),
+        dblAreas.begin(),
+        cellCells.begin(),
+        dblFaceWeights.begin(),
+        mgMinClusterSize_,
+        mgMaxClusterSize_,
+        options.begin(),
+        &nMoves,
         &nCoarseCells_,
         child_.begin()
     );
 
+#   endif
+
     Info<< "Number of coarse cells = " << nCoarseCells_ << endl;
 }
 
diff --git a/src/foam/Make/files b/src/foam/Make/files
index 7a39691bd..47a6daa96 100644
--- a/src/foam/Make/files
+++ b/src/foam/Make/files
@@ -26,6 +26,7 @@ $(ints)/label/label.C
 $(ints)/uLabel/uLabel.C
 
 primitives/Scalar/doubleScalar/doubleScalar.C
+primitives/Scalar/longDoubleScalar/longDoubleScalar.C
 primitives/Scalar/floatScalar/floatScalar.C
 primitives/Scalar/scalar/scalar.C
 primitives/DiagTensor/diagTensor/diagTensor.C
diff --git a/src/foam/db/IOstreams/Fstreams/IFstream.H b/src/foam/db/IOstreams/Fstreams/IFstream.H
index bfb8e910e..303e99635 100644
--- a/src/foam/db/IOstreams/Fstreams/IFstream.H
+++ b/src/foam/db/IOstreams/Fstreams/IFstream.H
@@ -109,8 +109,8 @@ public:
         IFstream
         (
             const fileName& pathname,
-            streamFormat format=ASCII,
-            versionNumber version=currentVersion
+            streamFormat format = ASCII,
+            versionNumber version = currentVersion
         );
 
 
diff --git a/src/foam/db/IOstreams/IOstreams/Istream.H b/src/foam/db/IOstreams/IOstreams/Istream.H
index dbafacb17..9de625a98 100644
--- a/src/foam/db/IOstreams/IOstreams/Istream.H
+++ b/src/foam/db/IOstreams/IOstreams/Istream.H
@@ -75,8 +75,8 @@ public:
         //- Set stream status
         Istream
         (
-            streamFormat format=ASCII,
-            versionNumber version=currentVersion,
+            streamFormat format = ASCII,
+            versionNumber version = currentVersion,
             compressionType compression=UNCOMPRESSED
         )
         :
@@ -122,6 +122,9 @@ public:
             //- Read a doubleScalar
             virtual Istream& read(doubleScalar&) = 0;
 
+            //- Read a longDoubleScalar
+            virtual Istream& read(longDoubleScalar&) = 0;
+
             //- Read binary block
             virtual Istream& read(char*, std::streamsize) = 0;
 
diff --git a/src/foam/db/IOstreams/IOstreams/Ostream.H b/src/foam/db/IOstreams/IOstreams/Ostream.H
index d9d1f8506..5c160da95 100644
--- a/src/foam/db/IOstreams/IOstreams/Ostream.H
+++ b/src/foam/db/IOstreams/IOstreams/Ostream.H
@@ -77,8 +77,8 @@ public:
         //- Set stream status
         Ostream
         (
-            streamFormat format=ASCII,
-            versionNumber version=currentVersion,
+            streamFormat format = ASCII,
+            versionNumber version = currentVersion,
             compressionType compression=UNCOMPRESSED
         )
         :
@@ -131,6 +131,9 @@ public:
             //- Write doubleScalar
             virtual Ostream& write(const doubleScalar) = 0;
 
+            //- Write doubleScalar
+            virtual Ostream& write(const longDoubleScalar) = 0;
+
             //- Write binary block
             virtual Ostream& write(const char*, std::streamsize) = 0;
 
diff --git a/src/foam/db/IOstreams/Pstreams/IPstream.C b/src/foam/db/IOstreams/Pstreams/IPstream.C
index 99280eeb7..4b6e57e3f 100644
--- a/src/foam/db/IOstreams/Pstreams/IPstream.C
+++ b/src/foam/db/IOstreams/Pstreams/IPstream.C
@@ -211,6 +211,21 @@ Foam::Istream& Foam::IPstream::read(token& t)
             return *this;
         }
 
+        // longDoubleScalar
+        case token::LONG_DOUBLE_SCALAR :
+        {
+            longDoubleScalar val;
+            if (read(val))
+            {
+                t = val;
+            }
+            else
+            {
+                t.setBad();
+            }
+            return *this;
+        }
+
         // Character (returned as a single character word) or error
         default:
         {
@@ -281,6 +296,13 @@ Foam::Istream& Foam::IPstream::read(doubleScalar& val)
 }
 
 
+Foam::Istream& Foam::IPstream::read(longDoubleScalar& val)
+{
+    readFromBuffer(val);
+    return *this;
+}
+
+
 Foam::Istream& Foam::IPstream::read(char* data, std::streamsize count)
 {
     if (format() != BINARY)
diff --git a/src/foam/db/IOstreams/Pstreams/IPstream.H b/src/foam/db/IOstreams/Pstreams/IPstream.H
index 5774d7c11..d19d466f9 100644
--- a/src/foam/db/IOstreams/Pstreams/IPstream.H
+++ b/src/foam/db/IOstreams/Pstreams/IPstream.H
@@ -84,7 +84,7 @@ public:
             const int fromProcNo,
             const label bufSize = 0,
             streamFormat format=BINARY,
-            versionNumber version=currentVersion
+            versionNumber version = currentVersion
         );
 
 
@@ -143,6 +143,9 @@ public:
             //- Read a doubleScalar
             Istream& read(doubleScalar&);
 
+            //- Read a longDoubleScalar
+            Istream& read(longDoubleScalar&);
+
             //- Read binary block
             Istream& read(char*, std::streamsize);
 
diff --git a/src/foam/db/IOstreams/Pstreams/OPstream.C b/src/foam/db/IOstreams/Pstreams/OPstream.C
index bbccc8fb3..05cad70b2 100644
--- a/src/foam/db/IOstreams/Pstreams/OPstream.C
+++ b/src/foam/db/IOstreams/Pstreams/OPstream.C
@@ -210,6 +210,14 @@ Foam::Ostream& Foam::OPstream::write(const doubleScalar val)
 }
 
 
+Foam::Ostream& Foam::OPstream::write(const longDoubleScalar val)
+{
+    write(char(token::LONG_DOUBLE_SCALAR));
+    writeToBuffer(val);
+    return *this;
+}
+
+
 Foam::Ostream& Foam::OPstream::write(const char* data, std::streamsize count)
 {
     if (format() != BINARY)
diff --git a/src/foam/db/IOstreams/Pstreams/OPstream.H b/src/foam/db/IOstreams/Pstreams/OPstream.H
index 5a0249223..22835337a 100644
--- a/src/foam/db/IOstreams/Pstreams/OPstream.H
+++ b/src/foam/db/IOstreams/Pstreams/OPstream.H
@@ -83,7 +83,7 @@ public:
             const int toProcNo,
             const label bufSize = 0,
             streamFormat format=BINARY,
-            versionNumber version=currentVersion
+            versionNumber version = currentVersion
         );
 
 
@@ -152,6 +152,9 @@ public:
             //- Write doubleScalar
             Ostream& write(const doubleScalar);
 
+            //- Write longDoubleScalar
+            Ostream& write(const longDoubleScalar);
+
             //- Write binary block
             Ostream& write(const char*, std::streamsize);
 
diff --git a/src/foam/db/IOstreams/Sstreams/ISstream.C b/src/foam/db/IOstreams/Sstreams/ISstream.C
index 3f3a4ae1f..d2c7d2054 100644
--- a/src/foam/db/IOstreams/Sstreams/ISstream.C
+++ b/src/foam/db/IOstreams/Sstreams/ISstream.C
@@ -489,6 +489,14 @@ Foam::Istream& Foam::ISstream::read(doubleScalar& val)
 }
 
 
+Foam::Istream& Foam::ISstream::read(longDoubleScalar& val)
+{
+    is_ >> val;
+    setState(is_.rdstate());
+    return *this;
+}
+
+
 // read binary block
 Foam::Istream& Foam::ISstream::read(char* buf, std::streamsize count)
 {
diff --git a/src/foam/db/IOstreams/Sstreams/ISstream.H b/src/foam/db/IOstreams/Sstreams/ISstream.H
index d4b20fad6..3bc92e61e 100644
--- a/src/foam/db/IOstreams/Sstreams/ISstream.H
+++ b/src/foam/db/IOstreams/Sstreams/ISstream.H
@@ -97,9 +97,9 @@ public:
         (
             istream& is,
             const string& name,
-            streamFormat format=ASCII,
-            versionNumber version=currentVersion,
-            compressionType compression=UNCOMPRESSED
+            streamFormat format = ASCII,
+            versionNumber version = currentVersion,
+            compressionType compression = UNCOMPRESSED
         );
 
 
@@ -165,6 +165,9 @@ public:
             //- Read a doubleScalar
             virtual Istream& read(doubleScalar&);
 
+            //- Read a longDoubleScalar
+            virtual Istream& read(longDoubleScalar&);
+
             //- Read binary block
             virtual Istream& read(char*, std::streamsize);
 
diff --git a/src/foam/db/IOstreams/Sstreams/OSstream.C b/src/foam/db/IOstreams/Sstreams/OSstream.C
index 11b15e000..4184de482 100644
--- a/src/foam/db/IOstreams/Sstreams/OSstream.C
+++ b/src/foam/db/IOstreams/Sstreams/OSstream.C
@@ -195,6 +195,14 @@ Foam::Ostream& Foam::OSstream::write(const doubleScalar val)
 }
 
 
+Foam::Ostream& Foam::OSstream::write(const longDoubleScalar val)
+{
+    os_ << val;
+    setState(os_.rdstate());
+    return *this;
+}
+
+
 Foam::Ostream& Foam::OSstream::write(const char* buf, std::streamsize count)
 {
     if (format() != BINARY)
diff --git a/src/foam/db/IOstreams/Sstreams/OSstream.H b/src/foam/db/IOstreams/Sstreams/OSstream.H
index d2adfb115..2f10c327e 100644
--- a/src/foam/db/IOstreams/Sstreams/OSstream.H
+++ b/src/foam/db/IOstreams/Sstreams/OSstream.H
@@ -90,8 +90,8 @@ public:
         (
             ostream& os,
             const string& name,
-            streamFormat format=ASCII,
-            versionNumber version=currentVersion,
+            streamFormat format = ASCII,
+            versionNumber version = currentVersion,
             compressionType compression=UNCOMPRESSED
         );
 
@@ -155,6 +155,9 @@ public:
             //- Write doubleScalar
             virtual Ostream& write(const doubleScalar);
 
+            //- Write longDoubleScalar
+            virtual Ostream& write(const longDoubleScalar);
+
             //- Write binary block
             virtual Ostream& write(const char*, std::streamsize);
 
diff --git a/src/foam/db/IOstreams/Sstreams/prefixOSstream.C b/src/foam/db/IOstreams/Sstreams/prefixOSstream.C
index e59639890..a29be1a7a 100644
--- a/src/foam/db/IOstreams/Sstreams/prefixOSstream.C
+++ b/src/foam/db/IOstreams/Sstreams/prefixOSstream.C
@@ -146,6 +146,13 @@ Foam::Ostream& Foam::prefixOSstream::write(const doubleScalar val)
 }
 
 
+Foam::Ostream& Foam::prefixOSstream::write(const longDoubleScalar val)
+{
+    checkWritePrefix();
+    return OSstream::write(val);
+}
+
+
 Foam::Ostream& Foam::prefixOSstream::write
 (
     const char* buf,
diff --git a/src/foam/db/IOstreams/Sstreams/prefixOSstream.H b/src/foam/db/IOstreams/Sstreams/prefixOSstream.H
index 21b395dda..3ea631726 100644
--- a/src/foam/db/IOstreams/Sstreams/prefixOSstream.H
+++ b/src/foam/db/IOstreams/Sstreams/prefixOSstream.H
@@ -73,8 +73,8 @@ public:
         (
             ostream& os,
             const string& name,
-            streamFormat format=ASCII,
-            versionNumber version=currentVersion,
+            streamFormat format = ASCII,
+            versionNumber version = currentVersion,
             compressionType compression=UNCOMPRESSED
         );
 
@@ -130,6 +130,9 @@ public:
             //- Write doubleScalar
             virtual Ostream& write(const doubleScalar);
 
+            //- Write doubleScalar
+            virtual Ostream& write(const longDoubleScalar);
+
             //- Write binary block
             virtual Ostream& write(const char*, std::streamsize);
 
diff --git a/src/foam/db/IOstreams/StringStreams/IStringStream.H b/src/foam/db/IOstreams/StringStreams/IStringStream.H
index 32c626016..5aada9c35 100644
--- a/src/foam/db/IOstreams/StringStreams/IStringStream.H
+++ b/src/foam/db/IOstreams/StringStreams/IStringStream.H
@@ -60,8 +60,8 @@ public:
         IStringStream
         (
             const string& buffer,
-            streamFormat format=ASCII,
-            versionNumber version=currentVersion
+            streamFormat format = ASCII,
+            versionNumber version = currentVersion
         )
         :
             ISstream
@@ -78,8 +78,8 @@ public:
         IStringStream
         (
             const char* buffer,
-            streamFormat format=ASCII,
-            versionNumber version=currentVersion
+            streamFormat format = ASCII,
+            versionNumber version = currentVersion
         )
         :
             ISstream
diff --git a/src/foam/db/IOstreams/StringStreams/OStringStream.H b/src/foam/db/IOstreams/StringStreams/OStringStream.H
index cb85b528e..a94e34eb1 100644
--- a/src/foam/db/IOstreams/StringStreams/OStringStream.H
+++ b/src/foam/db/IOstreams/StringStreams/OStringStream.H
@@ -59,8 +59,8 @@ public:
         //- Construct and set stream status
         OStringStream
         (
-            streamFormat format=ASCII,
-            versionNumber version=currentVersion
+            streamFormat format = ASCII,
+            versionNumber version = currentVersion
         )
         :
             OSstream
diff --git a/src/foam/db/IOstreams/Tstreams/ITstream.C b/src/foam/db/IOstreams/Tstreams/ITstream.C
index ba3d6192d..4a4183199 100644
--- a/src/foam/db/IOstreams/Tstreams/ITstream.C
+++ b/src/foam/db/IOstreams/Tstreams/ITstream.C
@@ -149,6 +149,13 @@ Foam::Istream& Foam::ITstream::read(doubleScalar&)
 }
 
 
+Foam::Istream& Foam::ITstream::read(longDoubleScalar&)
+{
+    notImplemented("Istream& ITstream::read(longDoubleScalar&)");
+    return *this;
+}
+
+
 Foam::Istream& Foam::ITstream::read(char*, std::streamsize)
 {
     notImplemented("Istream& ITstream::read(char*, std::streamsize)");
diff --git a/src/foam/db/IOstreams/Tstreams/ITstream.H b/src/foam/db/IOstreams/Tstreams/ITstream.H
index 43c163cce..f5aba41cb 100644
--- a/src/foam/db/IOstreams/Tstreams/ITstream.H
+++ b/src/foam/db/IOstreams/Tstreams/ITstream.H
@@ -70,8 +70,8 @@ public:
         (
             const string& name,
             const tokenList& tokens,
-            streamFormat format=ASCII,
-            versionNumber version=currentVersion
+            streamFormat format = ASCII,
+            versionNumber version = currentVersion
         )
         :
             Istream(format, version),
@@ -167,6 +167,9 @@ public:
             //- Read a doubleScalar
             virtual Istream& read(doubleScalar&);
 
+            //- Read a longDoubleScalar
+            virtual Istream& read(longDoubleScalar&);
+
             //- Read binary block
             virtual Istream& read(char*, std::streamsize);
 
diff --git a/src/foam/db/IOstreams/hashes/OSHA1stream.H b/src/foam/db/IOstreams/hashes/OSHA1stream.H
index 9e17fe9f8..8c3ab7d31 100644
--- a/src/foam/db/IOstreams/hashes/OSHA1stream.H
+++ b/src/foam/db/IOstreams/hashes/OSHA1stream.H
@@ -151,8 +151,8 @@ public:
         //- Construct and set stream status
         OSHA1stream
         (
-            streamFormat format=ASCII,
-            versionNumber version=currentVersion
+            streamFormat format = ASCII,
+            versionNumber version = currentVersion
         )
         :
             OSstream
diff --git a/src/foam/db/IOstreams/token/token.H b/src/foam/db/IOstreams/token/token.H
index efe7a2030..f0f045be4 100644
--- a/src/foam/db/IOstreams/token/token.H
+++ b/src/foam/db/IOstreams/token/token.H
@@ -81,6 +81,7 @@ public:
         LABEL,
         FLOAT_SCALAR,
         DOUBLE_SCALAR,
+        LONG_DOUBLE_SCALAR,
         COMPOUND,
 
         ERROR
@@ -258,6 +259,7 @@ private:
             label labelToken_;
             floatScalar floatScalarToken_;
             doubleScalar doubleScalarToken_;
+            longDoubleScalar longDoubleScalarToken_;
             mutable compound* compoundTokenPtr_;
         };
 
@@ -307,6 +309,9 @@ public:
         //- Construct doubleScalar token
         inline token(const doubleScalar, label lineNumber=0);
 
+        //- Construct longDoubleScalar token
+        inline token(const longDoubleScalar, label lineNumber=0);
+
         //- Construct from Istream
         token(Istream&);
 
@@ -344,6 +349,9 @@ public:
             inline bool isDoubleScalar() const;
             inline doubleScalar doubleScalarToken() const;
 
+            inline bool isLongDoubleScalar() const;
+            inline longDoubleScalar longDoubleScalarToken() const;
+
             inline bool isScalar() const;
             inline scalar scalarToken() const;
 
@@ -391,6 +399,7 @@ public:
             inline void operator=(const label);
             inline void operator=(const floatScalar);
             inline void operator=(const doubleScalar);
+            inline void operator=(const longDoubleScalar);
 
             inline void operator=(compound*);
 
@@ -404,6 +413,7 @@ public:
             inline bool operator==(const label) const;
             inline bool operator==(const floatScalar) const;
             inline bool operator==(const doubleScalar) const;
+            inline bool operator==(const longDoubleScalar) const;
 
 
         // Inequality
@@ -415,6 +425,7 @@ public:
             inline bool operator!=(const label) const;
             inline bool operator!=(const floatScalar) const;
             inline bool operator!=(const doubleScalar) const;
+            inline bool operator!=(const longDoubleScalar) const;
 
 
     // IOstream operators
diff --git a/src/foam/db/IOstreams/token/tokenI.H b/src/foam/db/IOstreams/token/tokenI.H
index 7e16adb83..425b7c338 100644
--- a/src/foam/db/IOstreams/token/tokenI.H
+++ b/src/foam/db/IOstreams/token/tokenI.H
@@ -103,6 +103,10 @@ inline token::token(const token& t)
             doubleScalarToken_ = t.doubleScalarToken_;
         break;
 
+        case LONG_DOUBLE_SCALAR:
+            longDoubleScalarToken_ = t.longDoubleScalarToken_;
+        break;
+
         case COMPOUND:
             compoundTokenPtr_ = t.compoundTokenPtr_;
             compoundTokenPtr_->refCount::operator++();
@@ -162,6 +166,15 @@ inline token::token(const doubleScalar s, label lineNumber)
 {}
 
 
+// Construct longDoubleScalar token
+inline token::token(const longDoubleScalar s, label lineNumber)
+:
+    type_(LONG_DOUBLE_SCALAR),
+    longDoubleScalarToken_(s),
+    lineNumber_(lineNumber)
+{}
+
+
 // * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
 
 // Delete token clearing the storage used by word or string
@@ -289,6 +302,7 @@ inline bool token::isDoubleScalar() const
     return (type_ == DOUBLE_SCALAR);
 }
 
+
 inline doubleScalar token::doubleScalarToken() const
 {
     if (type_ == DOUBLE_SCALAR)
@@ -303,9 +317,34 @@ inline doubleScalar token::doubleScalarToken() const
 }
 
 
+inline bool token::isLongDoubleScalar() const
+{
+    return (type_ == LONG_DOUBLE_SCALAR);
+}
+
+
+inline longDoubleScalar token::longDoubleScalarToken() const
+{
+    if (type_ == LONG_DOUBLE_SCALAR)
+    {
+        return longDoubleScalarToken_;
+    }
+    else
+    {
+        parseError("longDoubleScalar");
+        return 0.0;
+    }
+}
+
+
 inline bool token::isScalar() const
 {
-    return (type_ == FLOAT_SCALAR || type_ == DOUBLE_SCALAR);
+    return
+    (
+        type_ == FLOAT_SCALAR
+     || type_ == DOUBLE_SCALAR
+     || type_ == LONG_DOUBLE_SCALAR
+    );
 }
 
 inline scalar token::scalarToken() const
@@ -318,6 +357,10 @@ inline scalar token::scalarToken() const
     {
         return doubleScalarToken_;
     }
+    else if (type_ == LONG_DOUBLE_SCALAR)
+    {
+        return longDoubleScalarToken_;
+    }
     else
     {
         parseError("scalar");
@@ -420,6 +463,10 @@ inline void token::operator=(const token& t)
             doubleScalarToken_ = t.doubleScalarToken_;
         break;
 
+        case LONG_DOUBLE_SCALAR:
+            longDoubleScalarToken_ = t.longDoubleScalarToken_;
+        break;
+
         case COMPOUND:
             compoundTokenPtr_ = t.compoundTokenPtr_;
             compoundTokenPtr_->refCount::operator++();
@@ -432,6 +479,7 @@ inline void token::operator=(const token& t)
     lineNumber_ = t.lineNumber_;
 }
 
+
 inline void token::operator=(const punctuationToken p)
 {
     clear();
@@ -484,6 +532,13 @@ inline void token::operator=(const doubleScalar s)
     doubleScalarToken_ = s;
 }
 
+inline void token::operator=(const longDoubleScalar s)
+{
+    clear();
+    type_ = LONG_DOUBLE_SCALAR;
+    longDoubleScalarToken_ = s;
+}
+
 inline void token::operator=(token::compound* cPtr)
 {
     clear();
@@ -522,6 +577,9 @@ inline bool token::operator==(const token& t) const
         case DOUBLE_SCALAR:
             return equal(doubleScalarToken_, t.doubleScalarToken_);
 
+        case LONG_DOUBLE_SCALAR:
+            return equal(longDoubleScalarToken_, t.longDoubleScalarToken_);
+
         case COMPOUND:
             return compoundTokenPtr_ == t.compoundTokenPtr_;
 
@@ -562,6 +620,11 @@ inline bool token::operator==(const doubleScalar s) const
     return (type_ == DOUBLE_SCALAR && equal(doubleScalarToken_, s));
 }
 
+inline bool token::operator==(const longDoubleScalar s) const
+{
+    return (type_ == LONG_DOUBLE_SCALAR && equal(longDoubleScalarToken_, s));
+}
+
 inline bool token::operator!=(const token& t) const
 {
     return !operator==(t);
@@ -592,6 +655,11 @@ inline bool token::operator!=(const doubleScalar s) const
     return !operator==(s);
 }
 
+inline bool token::operator!=(const longDoubleScalar s) const
+{
+    return !operator==(s);
+}
+
 inline bool token::operator!=(const label l) const
 {
     return !operator==(l);
diff --git a/src/foam/db/IOstreams/token/tokenIO.C b/src/foam/db/IOstreams/token/tokenIO.C
index 9036b965d..58dcafcfd 100644
--- a/src/foam/db/IOstreams/token/tokenIO.C
+++ b/src/foam/db/IOstreams/token/tokenIO.C
@@ -85,6 +85,10 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const token& t)
             os << t.doubleScalarToken_;
         break;
 
+        case token::LONG_DOUBLE_SCALAR:
+            os << t.longDoubleScalarToken_;
+        break;
+
         case token::COMPOUND:
             os << *t.compoundTokenPtr_;
         break;
@@ -168,6 +172,10 @@ ostream& Foam::operator<<(ostream& os, const InfoProxy<token>& ip)
             os  << " the doubleScalar " << t.doubleScalarToken();
         break;
 
+        case token::LONG_DOUBLE_SCALAR:
+            os  << " the longDoubleScalar " << t.doubleScalarToken();
+        break;
+
         case token::COMPOUND:
         {
             if (t.compoundToken().empty())
@@ -238,6 +246,10 @@ Ostream& operator<<(Ostream& os, const InfoProxy<token>& ip)
             os  << " the doubleScalar " << t.doubleScalarToken();
         break;
 
+        case token::LONG_DOUBLE_SCALAR:
+            os  << " the longDoubleScalar " << t.longDoubleScalarToken();
+        break;
+
         case token::COMPOUND:
         {
             if (t.compoundToken().empty())
diff --git a/src/foam/db/Time/Time.C b/src/foam/db/Time/Time.C
index 758bd11ae..28a9d60b5 100644
--- a/src/foam/db/Time/Time.C
+++ b/src/foam/db/Time/Time.C
@@ -75,7 +75,7 @@ void Foam::Time::adjustDeltaT()
     {
         scalar timeToNextWrite = max
         (
-            0.0,
+            scalar(0),
             (outputTimeIndex_ + 1)*writeInterval_ - (value() - startTime_)
         );
 
diff --git a/src/foam/fields/PointPatchFields/constraint/processor/ProcessorPointPatchField.C b/src/foam/fields/PointPatchFields/constraint/processor/ProcessorPointPatchField.C
index 7178f5eb7..382dbbfc7 100644
--- a/src/foam/fields/PointPatchFields/constraint/processor/ProcessorPointPatchField.C
+++ b/src/foam/fields/PointPatchFields/constraint/processor/ProcessorPointPatchField.C
@@ -516,8 +516,8 @@ evaluate
 
             // ZT, 22-07-2014 - point ordering is not same
             // at master and slave side after topology change
-            const labelList& neiPoints =
-                procPatch_.procPolyPatch().neighbPoints();
+//             const labelList& neiPoints =
+//                 procPatch_.procPolyPatch().neighbPoints();
 
             // VS, 2015-04-12 - This doesn't work in parallel!
 //            tpn =
diff --git a/src/foam/meshes/meshShapes/face/face.C b/src/foam/meshes/meshShapes/face/face.C
index d674e00ab..71c6b5efc 100644
--- a/src/foam/meshes/meshShapes/face/face.C
+++ b/src/foam/meshes/meshShapes/face/face.C
@@ -225,7 +225,9 @@ Foam::label Foam::face::split
             splitEdge /= Foam::mag(splitEdge) + VSMALL;
 
             const scalar splitCos = splitEdge & rightEdge;
-            const scalar splitAngle = acos(max(-1.0, min(1.0, splitCos)));
+            const scalar splitAngle =
+                acos(max(scalar(-1), min(scalar(1), splitCos)));
+
             const scalar angleDiff = fabs(splitAngle - bisectAngle);
 
             if (angleDiff < minDiff)
diff --git a/src/foam/meshes/primitiveShapes/triangle/triangleI.H b/src/foam/meshes/primitiveShapes/triangle/triangleI.H
index b417df5e5..2c5fec51a 100644
--- a/src/foam/meshes/primitiveShapes/triangle/triangleI.H
+++ b/src/foam/meshes/primitiveShapes/triangle/triangleI.H
@@ -400,7 +400,7 @@ inline scalar triangle<Point, PointRef>::circumRadius() const
     {
         scalar a = (d1 + d2)*(d2 + d3)*(d3 + d1) / denom;
 
-        return 0.5*Foam::sqrt(min(GREAT, max(0, a)));
+        return 0.5*Foam::sqrt(min(GREAT, max(scalar(0), a)));
     }
 }
 
@@ -411,9 +411,7 @@ inline scalar triangle<Point, PointRef>::quality() const
     return
         mag()
        /(
-            mathematicalConstant::pi
-          *Foam::sqr(circumRadius())
-          *0.413497
+            0.413497*mathematicalConstant::pi*sqr(circumRadius())
           + VSMALL
         );
 }
@@ -775,8 +773,8 @@ inline bool triangle<Point, PointRef>::classify
     if (hit)
     {
         // alpha,beta might get negative due to precision errors
-        alpha = max(0.0, min(1.0, alpha));
-        beta = max(0.0, min(1.0, beta));
+        alpha = max(scalar(0), min(scalar(1), alpha));
+        beta = max(scalar(0), min(scalar(1), beta));
     }
 
     nearType = NONE;
diff --git a/src/foam/primitives/Scalar/Scalar.H b/src/foam/primitives/Scalar/Scalar.H
index 7fdef931a..0f3f4e676 100644
--- a/src/foam/primitives/Scalar/Scalar.H
+++ b/src/foam/primitives/Scalar/Scalar.H
@@ -161,11 +161,11 @@ inline Scalar sqrtSumSqr(const Scalar a, const Scalar b)
 
     if (maga > magb)
     {
-        return maga*sqrt(1.0 + sqr(magb/maga));
+        return maga*::sqrt(1.0 + sqr(magb/maga));
     }
     else
     {
-        return magb < ScalarVSMALL ? 0.0 : magb*sqrt(1.0 + sqr(maga/magb));
+        return magb < ScalarVSMALL ? 0.0 : magb*::sqrt(1.0 + sqr(maga/magb));
     }
 }
 
diff --git a/src/foam/primitives/Scalar/doubleFloat.H b/src/foam/primitives/Scalar/doubleFloat.H
index c34274233..0f2e420ca 100644
--- a/src/foam/primitives/Scalar/doubleFloat.H
+++ b/src/foam/primitives/Scalar/doubleFloat.H
@@ -84,6 +84,13 @@ inline double pow(const type1 s, const type2 e)   \
     return ::pow(double(s), double(e));           \
 }
 
+MAXMINPOW(long double, long double, long double)
+MAXMINPOW(long double, long double, double)
+MAXMINPOW(long double, long double, float)
+MAXMINPOW(long double, long double, int)
+MAXMINPOW(long double, double, long double)
+MAXMINPOW(long double, float, long double)
+MAXMINPOW(long double, int, long double)
 
 MAXMINPOW(double, double, double)
 MAXMINPOW(double, double, float)
diff --git a/src/foam/primitives/Scalar/longDoubleScalar/longDoubleScalar.C b/src/foam/primitives/Scalar/longDoubleScalar/longDoubleScalar.C
new file mode 100644
index 000000000..5cc7567be
--- /dev/null
+++ b/src/foam/primitives/Scalar/longDoubleScalar/longDoubleScalar.C
@@ -0,0 +1,43 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     |
+    \\  /    A nd           | For copyright notice see file Copyright
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation, either version 3 of the License, or (at your
+    option) any later version.
+
+    foam-extend is distributed in the hope that it will be useful, but
+    WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+    General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "longDoubleScalar.H"
+#include "IOstreams.H"
+
+#include <sstream>
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#define Scalar longDoubleScalar
+#define ScalarVGREAT longDoubleScalarVGREAT
+#define ScalarVSMALL longDoubleScalarVSMALL
+#define readScalar readLongDoubleScalar
+#include "Scalar.C"
+#undef Scalar
+#undef ScalarVGREAT
+#undef ScalarVSMALL
+#undef readScalar
+
+// ************************************************************************* //
diff --git a/src/foam/primitives/Scalar/longDoubleScalar/longDoubleScalar.H b/src/foam/primitives/Scalar/longDoubleScalar/longDoubleScalar.H
new file mode 100644
index 000000000..184be0933
--- /dev/null
+++ b/src/foam/primitives/Scalar/longDoubleScalar/longDoubleScalar.H
@@ -0,0 +1,115 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     |
+    \\  /    A nd           | For copyright notice see file Copyright
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation, either version 3 of the License, or (at your
+    option) any later version.
+
+    foam-extend is distributed in the hope that it will be useful, but
+    WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+    General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Typedef
+    Foam::doubleScalar
+
+Description
+    Long double precision floating point scalar type.
+
+SourceFiles
+    longDoubleScalar.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef longDoubleScalar_H
+#define longDoubleScalar_H
+
+#include "doubleFloat.H"
+#include "direction.H"
+#include "word.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+typedef long double longDoubleScalar;
+
+// Largest and smallest scalar values allowed in certain parts of the code.
+// (15 is the number of significant figures in an
+//  IEEE double precision number.  See limits.h or float.h)
+static const longDoubleScalar longDoubleScalarGREAT = 1.0e+15;
+static const longDoubleScalar longDoubleScalarVGREAT = 1.0e+300;
+static const longDoubleScalar longDoubleScalarROOTVGREAT = 1.0e+150;
+static const longDoubleScalar longDoubleScalarSMALL = 1.0e-15;
+static const longDoubleScalar longDoubleScalarVSMALL = 1.0e-300;
+static const longDoubleScalar longDoubleScalarROOTVSMALL = 1.0e-150;
+
+
+#define Scalar longDoubleScalar
+#define ScalarVGREAT longDoubleScalarVGREAT
+#define ScalarVSMALL longDoubleScalarVSMALL
+#define readScalar readLongDoubleScalar
+
+inline Scalar mag(const Scalar s)
+{
+    return ::fabs(s);
+}
+
+#define transFunc(func)            \
+inline Scalar func(const Scalar s) \
+{                                  \
+    return ::func(s);              \
+}
+
+#include "Scalar.H"
+
+inline Scalar hypot(const Scalar x, const Scalar y)
+{
+    return ::hypot(x, y);
+}
+
+inline Scalar atan2(const Scalar y, const Scalar x)
+{
+    return ::atan2(y, x);
+}
+
+inline Scalar jn(const int n, const Scalar s)
+{
+    return ::jn(n, s);
+}
+
+inline Scalar yn(const int n, const Scalar s)
+{
+    return ::yn(n, s);
+}
+
+#undef Scalar
+#undef ScalarVGREAT
+#undef ScalarVSMALL
+#undef readScalar
+#undef transFunc
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/foam/primitives/Scalar/scalar/scalar.H b/src/foam/primitives/Scalar/scalar/scalar.H
index 929d8018f..db0ac4920 100644
--- a/src/foam/primitives/Scalar/scalar/scalar.H
+++ b/src/foam/primitives/Scalar/scalar/scalar.H
@@ -38,6 +38,7 @@ SourceFiles
 
 #include "floatScalar.H"
 #include "doubleScalar.H"
+#include "longDoubleScalar.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -77,6 +78,24 @@ namespace Foam
     scalar readScalar(Istream& is);
 }
 
+#elif defined(WM_LDP)
+
+// Define scalar as a long double
+
+namespace Foam
+{
+    typedef longDoubleScalar scalar;
+
+    static const scalar GREAT = longDoubleScalarGREAT;
+    static const scalar VGREAT = longDoubleScalarVGREAT;
+    static const scalar ROOTVGREAT = longDoubleScalarROOTVGREAT;
+    static const scalar SMALL = longDoubleScalarSMALL;
+    static const scalar VSMALL = longDoubleScalarVSMALL;
+    static const scalar ROOTVSMALL = longDoubleScalarROOTVSMALL;
+
+    scalar readScalar(Istream& is);
+}
+
 #endif
 
 
diff --git a/src/foam/primitives/Tensor/tensor/tensor.C b/src/foam/primitives/Tensor/tensor/tensor.C
index be1e70201..f036c8288 100644
--- a/src/foam/primitives/Tensor/tensor/tensor.C
+++ b/src/foam/primitives/Tensor/tensor/tensor.C
@@ -118,7 +118,7 @@ vector eigenValues(const tensor& t)
 
             if (disc >= -SMALL)
             {
-                scalar q = -0.5*sqrt(max(0.0, disc));
+                scalar q = -0.5*sqrt(max(scalar(0), disc));
 
                 i = 0;
                 ii = -0.5*a + q;
@@ -381,7 +381,7 @@ vector eigenValues(const symmTensor& t)
         // If there is a zero root
         if (mag(c) < SMALL)
         {
-            const scalar disc = Foam::max(sqr(a) - 4*b, 0.0);
+            const scalar disc = Foam::max(sqr(a) - 4*b, scalar(0));
 
             scalar q = -0.5*sqrt(max(scalar(0), disc));
 
diff --git a/src/foam/primitives/transform/symmTransform.H b/src/foam/primitives/transform/symmTransform.H
index c69862e6d..11fe83afd 100644
--- a/src/foam/primitives/transform/symmTransform.H
+++ b/src/foam/primitives/transform/symmTransform.H
@@ -117,163 +117,163 @@ inline SymmTensor4thOrder<Cmpt> transform
  const SymmTensor4thOrder<Cmpt>& C
  )
 {
-  //- represent fourth order tensors in 6x6 notation then rotation is given by
-  //- C_rotated_af = A_ba * C_cd * A_ef
-  //- where A is a function of stt
+    //- represent fourth order tensors in 6x6 notation.  Rotation is given by
+    //- C_rotated_af = A_ba * C_cd * A_ef
+    //- where A is a function of stt
 
-  const double s = ::sqrt(2);
-  const double A[6][6] =
+    const scalar s = ::sqrt(2);
+    const scalar A[6][6] =
     {
-      { stt.xx()*stt.xx(), stt.xy()*stt.xy(), stt.xz()*stt.xz(), s*stt.xx()*stt.xy(), s*stt.xy()*stt.xz(), s*stt.xz()*stt.xx() },
-      { stt.xy()*stt.xy(), stt.yy()*stt.yy(), stt.yz()*stt.yz(), s*stt.xy()*stt.yy(), s*stt.yy()*stt.yz(), s*stt.yz()*stt.xy() },
-      { stt.xz()*stt.xz(), stt.yz()*stt.yz(), stt.zz()*stt.zz(), s*stt.xz()*stt.yz(), s*stt.yz()*stt.zz(), s*stt.zz()*stt.xz() },
-      { s*stt.xx()*stt.xy(), s*stt.xy()*stt.yy(), s*stt.xz()*stt.yz(),
-	(stt.xy()*stt.xy()+stt.xx()*stt.yy()), (stt.xz()*stt.yy()+stt.xy()*stt.yz()), (stt.xx()*stt.yz()+stt.xz()*stt.xy()) },
-      { s*stt.xy()*stt.xz(), s*stt.yy()*stt.yz(), s*stt.yz()*stt.zz(),
-	(stt.yy()*stt.xz()+stt.xy()*stt.yz()), (stt.yz()*stt.yz()+stt.yy()*stt.zz()), (stt.xy()*stt.zz()+stt.yz()*stt.xz()) },
-      { s*stt.xz()*stt.xx(), s*stt.yz()*stt.xy(), s*stt.zz()*stt.xz(),
-	(stt.yz()*stt.xx()+stt.xz()*stt.xy()), (stt.zz()*stt.xy()+stt.yz()*stt.xz()), (stt.xz()*stt.xz()+stt.zz()*stt.xx()) }
+        { stt.xx()*stt.xx(), stt.xy()*stt.xy(), stt.xz()*stt.xz(), s*stt.xx()*stt.xy(), s*stt.xy()*stt.xz(), s*stt.xz()*stt.xx() },
+        { stt.xy()*stt.xy(), stt.yy()*stt.yy(), stt.yz()*stt.yz(), s*stt.xy()*stt.yy(), s*stt.yy()*stt.yz(), s*stt.yz()*stt.xy() },
+        { stt.xz()*stt.xz(), stt.yz()*stt.yz(), stt.zz()*stt.zz(), s*stt.xz()*stt.yz(), s*stt.yz()*stt.zz(), s*stt.zz()*stt.xz() },
+        { s*stt.xx()*stt.xy(), s*stt.xy()*stt.yy(), s*stt.xz()*stt.yz(),
+          (stt.xy()*stt.xy()+stt.xx()*stt.yy()), (stt.xz()*stt.yy()+stt.xy()*stt.yz()), (stt.xx()*stt.yz()+stt.xz()*stt.xy()) },
+        { s*stt.xy()*stt.xz(), s*stt.yy()*stt.yz(), s*stt.yz()*stt.zz(),
+          (stt.yy()*stt.xz()+stt.xy()*stt.yz()), (stt.yz()*stt.yz()+stt.yy()*stt.zz()), (stt.xy()*stt.zz()+stt.yz()*stt.xz()) },
+        { s*stt.xz()*stt.xx(), s*stt.yz()*stt.xy(), s*stt.zz()*stt.xz(),
+          (stt.yz()*stt.xx()+stt.xz()*stt.xy()), (stt.zz()*stt.xy()+stt.yz()*stt.xz()), (stt.xz()*stt.xz()+stt.zz()*stt.xx()) }
     };
 
-  return symmTensor4thOrder
+    return symmTensor4thOrder
     (
-     // xxxx
-     A[0][0] * C.xxxx() * A[0][0] +
-     A[1][0] * C.xxyy() * A[0][0] +
-     A[2][0] * C.xxzz() * A[0][0] +
-     A[0][0] * C.xxyy() * A[1][0] +
-     A[1][0] * C.yyyy() * A[1][0] +
-     A[2][0] * C.yyzz() * A[1][0] +
-     A[0][0] * C.xxzz() * A[2][0] +
-     A[1][0] * C.yyzz() * A[2][0] +
-     A[2][0] * C.zzzz() * A[2][0] +
-     A[3][0] * C.xyxy() * A[3][0] +
-     A[4][0] * C.yzyz() * A[4][0] +
-     A[5][0] * C.zxzx() * A[5][0],
+        // xxxx
+        A[0][0] * C.xxxx() * A[0][0] +
+        A[1][0] * C.xxyy() * A[0][0] +
+        A[2][0] * C.xxzz() * A[0][0] +
+        A[0][0] * C.xxyy() * A[1][0] +
+        A[1][0] * C.yyyy() * A[1][0] +
+        A[2][0] * C.yyzz() * A[1][0] +
+        A[0][0] * C.xxzz() * A[2][0] +
+        A[1][0] * C.yyzz() * A[2][0] +
+        A[2][0] * C.zzzz() * A[2][0] +
+        A[3][0] * C.xyxy() * A[3][0] +
+        A[4][0] * C.yzyz() * A[4][0] +
+        A[5][0] * C.zxzx() * A[5][0],
 
-     // xxyy
-     A[0][0] * C.xxxx() * A[0][1] +
-     A[1][0] * C.xxyy() * A[0][1] +
-     A[2][0] * C.xxzz() * A[0][1] +
-     A[0][0] * C.xxyy() * A[1][1] +
-     A[1][0] * C.yyyy() * A[1][1] +
-     A[2][0] * C.yyzz() * A[1][1] +
-     A[0][0] * C.xxzz() * A[2][1] +
-     A[1][0] * C.yyzz() * A[2][1] +
-     A[2][0] * C.zzzz() * A[2][1] +
-     A[3][0] * C.xyxy() * A[3][1] +
-     A[4][0] * C.yzyz() * A[4][1] +
-     A[5][0] * C.zxzx() * A[5][1],
+        // xxyy
+        A[0][0] * C.xxxx() * A[0][1] +
+        A[1][0] * C.xxyy() * A[0][1] +
+        A[2][0] * C.xxzz() * A[0][1] +
+        A[0][0] * C.xxyy() * A[1][1] +
+        A[1][0] * C.yyyy() * A[1][1] +
+        A[2][0] * C.yyzz() * A[1][1] +
+        A[0][0] * C.xxzz() * A[2][1] +
+        A[1][0] * C.yyzz() * A[2][1] +
+        A[2][0] * C.zzzz() * A[2][1] +
+        A[3][0] * C.xyxy() * A[3][1] +
+        A[4][0] * C.yzyz() * A[4][1] +
+        A[5][0] * C.zxzx() * A[5][1],
 
-     // xzz
-     A[0][0] * C.xxxx() * A[0][2] +
-     A[1][0] * C.xxyy() * A[0][2] +
-     A[2][0] * C.xxzz() * A[0][2] +
-     A[0][0] * C.xxyy() * A[1][2] +
-     A[1][0] * C.yyyy() * A[1][2] +
-     A[2][0] * C.yyzz() * A[1][2] +
-     A[0][0] * C.xxzz() * A[2][2] +
-     A[1][0] * C.yyzz() * A[2][2] +
-     A[2][0] * C.zzzz() * A[2][2] +
-     A[3][0] * C.xyxy() * A[3][2] +
-     A[4][0] * C.yzyz() * A[4][2] +
-     A[5][0] * C.zxzx() * A[5][2],
+        // xzz
+        A[0][0] * C.xxxx() * A[0][2] +
+        A[1][0] * C.xxyy() * A[0][2] +
+        A[2][0] * C.xxzz() * A[0][2] +
+        A[0][0] * C.xxyy() * A[1][2] +
+        A[1][0] * C.yyyy() * A[1][2] +
+        A[2][0] * C.yyzz() * A[1][2] +
+        A[0][0] * C.xxzz() * A[2][2] +
+        A[1][0] * C.yyzz() * A[2][2] +
+        A[2][0] * C.zzzz() * A[2][2] +
+        A[3][0] * C.xyxy() * A[3][2] +
+        A[4][0] * C.yzyz() * A[4][2] +
+        A[5][0] * C.zxzx() * A[5][2],
 
-     // yyyy
-     A[0][1] * C.xxxx() * A[0][1] +
-     A[1][1] * C.xxyy() * A[0][1] +
-     A[2][1] * C.xxzz() * A[0][1] +
-     A[0][1] * C.xxyy() * A[1][1] +
-     A[1][1] * C.yyyy() * A[1][1] +
-     A[2][1] * C.yyzz() * A[1][1] +
-     A[0][1] * C.xxzz() * A[2][1] +
-     A[1][1] * C.yyzz() * A[2][1] +
-     A[2][1] * C.zzzz() * A[2][1] +
-     A[3][1] * C.xyxy() * A[3][1] +
-     A[4][1] * C.yzyz() * A[4][1] +
-     A[5][1] * C.zxzx() * A[5][1],
+        // yyyy
+        A[0][1] * C.xxxx() * A[0][1] +
+        A[1][1] * C.xxyy() * A[0][1] +
+        A[2][1] * C.xxzz() * A[0][1] +
+        A[0][1] * C.xxyy() * A[1][1] +
+        A[1][1] * C.yyyy() * A[1][1] +
+        A[2][1] * C.yyzz() * A[1][1] +
+        A[0][1] * C.xxzz() * A[2][1] +
+        A[1][1] * C.yyzz() * A[2][1] +
+        A[2][1] * C.zzzz() * A[2][1] +
+        A[3][1] * C.xyxy() * A[3][1] +
+        A[4][1] * C.yzyz() * A[4][1] +
+        A[5][1] * C.zxzx() * A[5][1],
 
-     // yyzz
-     A[0][1] * C.xxxx() * A[0][2] +
-     A[1][1] * C.xxyy() * A[0][2] +
-     A[2][1] * C.xxzz() * A[0][2] +
-     A[0][1] * C.xxyy() * A[1][2] +
-     A[1][1] * C.yyyy() * A[1][2] +
-     A[2][1] * C.yyzz() * A[1][2] +
-     A[0][1] * C.xxzz() * A[2][2] +
-     A[1][1] * C.yyzz() * A[2][2] +
-     A[2][1] * C.zzzz() * A[2][2] +
-     A[3][1] * C.xyxy() * A[3][2] +
-     A[4][1] * C.yzyz() * A[4][2] +
-     A[5][1] * C.zxzx() * A[5][2],
+        // yyzz
+        A[0][1] * C.xxxx() * A[0][2] +
+        A[1][1] * C.xxyy() * A[0][2] +
+        A[2][1] * C.xxzz() * A[0][2] +
+        A[0][1] * C.xxyy() * A[1][2] +
+        A[1][1] * C.yyyy() * A[1][2] +
+        A[2][1] * C.yyzz() * A[1][2] +
+        A[0][1] * C.xxzz() * A[2][2] +
+        A[1][1] * C.yyzz() * A[2][2] +
+        A[2][1] * C.zzzz() * A[2][2] +
+        A[3][1] * C.xyxy() * A[3][2] +
+        A[4][1] * C.yzyz() * A[4][2] +
+        A[5][1] * C.zxzx() * A[5][2],
 
-     // zzzz
-     A[0][2] * C.xxxx() * A[0][2] +
-     A[1][2] * C.xxyy() * A[0][2] +
-     A[2][2] * C.xxzz() * A[0][2] +
-     A[0][2] * C.xxyy() * A[1][2] +
-     A[1][2] * C.yyyy() * A[1][2] +
-     A[2][2] * C.yyzz() * A[1][2] +
-     A[0][2] * C.xxzz() * A[2][2] +
-     A[1][2] * C.yyzz() * A[2][2] +
-     A[2][2] * C.zzzz() * A[2][2] +
-     A[3][2] * C.xyxy() * A[3][2] +
-     A[4][2] * C.yzyz() * A[4][2] +
-     A[5][2] * C.zxzx() * A[5][2],
+        // zzzz
+        A[0][2] * C.xxxx() * A[0][2] +
+        A[1][2] * C.xxyy() * A[0][2] +
+        A[2][2] * C.xxzz() * A[0][2] +
+        A[0][2] * C.xxyy() * A[1][2] +
+        A[1][2] * C.yyyy() * A[1][2] +
+        A[2][2] * C.yyzz() * A[1][2] +
+        A[0][2] * C.xxzz() * A[2][2] +
+        A[1][2] * C.yyzz() * A[2][2] +
+        A[2][2] * C.zzzz() * A[2][2] +
+        A[3][2] * C.xyxy() * A[3][2] +
+        A[4][2] * C.yzyz() * A[4][2] +
+        A[5][2] * C.zxzx() * A[5][2],
 
-     // xyxy
-     A[0][3] * C.xxxx() * A[0][3] +
-     A[1][3] * C.xxyy() * A[0][3] +
-     A[2][3] * C.xxzz() * A[0][3] +
-     A[0][3] * C.xxyy() * A[1][3] +
-     A[1][3] * C.yyyy() * A[1][3] +
-     A[2][3] * C.yyzz() * A[1][3] +
-     A[0][3] * C.xxzz() * A[2][3] +
-     A[1][3] * C.yyzz() * A[2][3] +
-     A[2][3] * C.zzzz() * A[2][3] +
-     A[3][3] * C.xyxy() * A[3][3] +
-     A[4][3] * C.yzyz() * A[4][3] +
-     A[5][3] * C.zxzx() * A[5][3],
+        // xyxy
+        A[0][3] * C.xxxx() * A[0][3] +
+        A[1][3] * C.xxyy() * A[0][3] +
+        A[2][3] * C.xxzz() * A[0][3] +
+        A[0][3] * C.xxyy() * A[1][3] +
+        A[1][3] * C.yyyy() * A[1][3] +
+        A[2][3] * C.yyzz() * A[1][3] +
+        A[0][3] * C.xxzz() * A[2][3] +
+        A[1][3] * C.yyzz() * A[2][3] +
+        A[2][3] * C.zzzz() * A[2][3] +
+        A[3][3] * C.xyxy() * A[3][3] +
+        A[4][3] * C.yzyz() * A[4][3] +
+        A[5][3] * C.zxzx() * A[5][3],
 
-     // yzyz
-     A[0][4] * C.xxxx() * A[0][4] +
-     A[1][4] * C.xxyy() * A[0][4] +
-     A[2][4] * C.xxzz() * A[0][4] +
-     A[0][4] * C.xxyy() * A[1][4] +
-     A[1][4] * C.yyyy() * A[1][4] +
-     A[2][4] * C.yyzz() * A[1][4] +
-     A[0][4] * C.xxzz() * A[2][4] +
-     A[1][4] * C.yyzz() * A[2][4] +
-     A[2][4] * C.zzzz() * A[2][4] +
-     A[3][4] * C.xyxy() * A[3][4] +
-     A[4][4] * C.yzyz() * A[4][4] +
-     A[5][4] * C.zxzx() * A[5][4],
+        // yzyz
+        A[0][4] * C.xxxx() * A[0][4] +
+        A[1][4] * C.xxyy() * A[0][4] +
+        A[2][4] * C.xxzz() * A[0][4] +
+        A[0][4] * C.xxyy() * A[1][4] +
+        A[1][4] * C.yyyy() * A[1][4] +
+        A[2][4] * C.yyzz() * A[1][4] +
+        A[0][4] * C.xxzz() * A[2][4] +
+        A[1][4] * C.yyzz() * A[2][4] +
+        A[2][4] * C.zzzz() * A[2][4] +
+        A[3][4] * C.xyxy() * A[3][4] +
+        A[4][4] * C.yzyz() * A[4][4] +
+        A[5][4] * C.zxzx() * A[5][4],
 
-     // zxzx
-     A[0][5] * C.xxxx() * A[0][5] +
-     A[1][5] * C.xxyy() * A[0][5] +
-     A[2][5] * C.xxzz() * A[0][5] +
-     A[0][5] * C.xxyy() * A[1][5] +
-     A[1][5] * C.yyyy() * A[1][5] +
-     A[2][5] * C.yyzz() * A[1][5] +
-     A[0][5] * C.xxzz() * A[2][5] +
-     A[1][5] * C.yyzz() * A[2][5] +
-     A[2][5] * C.zzzz() * A[2][5] +
-     A[3][5] * C.xyxy() * A[3][5] +
-     A[4][5] * C.yzyz() * A[4][5] +
-     A[5][5] * C.zxzx() * A[5][5]
-     );
+        // zxzx
+        A[0][5] * C.xxxx() * A[0][5] +
+        A[1][5] * C.xxyy() * A[0][5] +
+        A[2][5] * C.xxzz() * A[0][5] +
+        A[0][5] * C.xxyy() * A[1][5] +
+        A[1][5] * C.yyyy() * A[1][5] +
+        A[2][5] * C.yyzz() * A[1][5] +
+        A[0][5] * C.xxzz() * A[2][5] +
+        A[1][5] * C.yyzz() * A[2][5] +
+        A[2][5] * C.zzzz() * A[2][5] +
+        A[3][5] * C.xyxy() * A[3][5] +
+        A[4][5] * C.yzyz() * A[4][5] +
+        A[5][5] * C.zxzx() * A[5][5]
+    );
 }
 
 
 template<class Cmpt>
 inline DiagTensor<Cmpt> transform
 (
- const symmTensor& stt,
- const DiagTensor<Cmpt>& dt
+    const symmTensor& stt,
+    const DiagTensor<Cmpt>& dt
  )
 {
-  return dt;
+    return dt;
 }
 
 
@@ -328,21 +328,26 @@ inline symmTensor transformMask<symmTensor>(const symmTensor& st)
 
 
 template<>
-inline symmTensor4thOrder transformMask<symmTensor4thOrder>(const symmTensor& st)
+inline symmTensor4thOrder transformMask<symmTensor4thOrder>
+(
+    const symmTensor& st
+)
 {
-  notImplemented("symmTransform.H\n"
-		 "template<>\n"
-		 "inline symmTensor4thOrder transformMask<symmTensor4thOrder>(const symmTensor& st)\n"
-		 "not implemented");
+    notImplemented
+    (
+        "template<>\n"
+        "inline symmTensor4thOrder transformMask<symmTensor4thOrder>"
+        "(const symmTensor& st)"
+    );
 
-  return symmTensor4thOrder::zero;
+    return symmTensor4thOrder::zero;
 }
 
 
 template<>
 inline diagTensor transformMask<diagTensor>(const symmTensor& st)
 {
-  return diagTensor(st.xx(), st.yy(), st.zz());
+    return diagTensor(st.xx(), st.yy(), st.zz());
 }
 
 
diff --git a/src/foam/primitives/transform/transform.H b/src/foam/primitives/transform/transform.H
index a0d1e083e..5e87abbf1 100644
--- a/src/foam/primitives/transform/transform.H
+++ b/src/foam/primitives/transform/transform.H
@@ -146,151 +146,156 @@ inline SphericalTensor<Cmpt> transform
 
 
 template<class Cmpt>
-inline SymmTensor4thOrder<Cmpt> transform(const tensor& L, const SymmTensor4thOrder<Cmpt>& C)
+inline SymmTensor4thOrder<Cmpt> transform
+(
+    const tensor& L,
+    const SymmTensor4thOrder<Cmpt>& C
+)
 {
-  //- represent fourth order tensors in 6x6 notation then rotation is given by
-  //- C_rotated_af = A_ba * C_cd * A_ef
-  //- where A is a function of L
+    //- represent fourth order tensors in 6x6 notation.  Rotation is given by
+    //- C_rotated_af = A_ba * C_cd * A_ef
+    //- where A is a function of L
 
-  const double s = ::sqrt(2);
-  const double A[6][6] =
+    const scalar s = ::sqrt(2);
+    const scalar A[6][6] =
     {
-      { L.xx()*L.xx(), L.xy()*L.xy(), L.xz()*L.xz(), s*L.xx()*L.xy(), s*L.xy()*L.xz(), s*L.xz()*L.xx() },
-      { L.yx()*L.yx(), L.yy()*L.yy(), L.yz()*L.yz(), s*L.yx()*L.yy(), s*L.yy()*L.yz(), s*L.yz()*L.yx() },
-      { L.zx()*L.zx(), L.zy()*L.zy(), L.zz()*L.zz(), s*L.zx()*L.zy(), s*L.zy()*L.zz(), s*L.zz()*L.zx() },
-      { s*L.xx()*L.yx(), s*L.xy()*L.yy(), s*L.xz()*L.yz(), (L.xy()*L.yx()+L.xx()*L.yy()), (L.xz()*L.yy()+L.xy()*L.yz()), (L.xx()*L.yz()+L.xz()*L.yx()) },
-      { s*L.yx()*L.zx(), s*L.yy()*L.zy(), s*L.yz()*L.zz(), (L.yy()*L.zx()+L.yx()*L.zy()), (L.yz()*L.zy()+L.yy()*L.zz()), (L.yx()*L.zz()+L.yz()*L.zx()) },
-      { s*L.zx()*L.xx(), s*L.zy()*L.xy(), s*L.zz()*L.xz(), (L.zy()*L.xx()+L.zx()*L.xy()), (L.zz()*L.xy()+L.zy()*L.xz()), (L.zx()*L.xz()+L.zz()*L.xx()) }
+        {
+            L.xx()*L.xx(), L.xy()*L.xy(), L.xz()*L.xz(), s*L.xx()*L.xy(), s*L.xy()*L.xz(), s*L.xz()*L.xx() },
+        { L.yx()*L.yx(), L.yy()*L.yy(), L.yz()*L.yz(), s*L.yx()*L.yy(), s*L.yy()*L.yz(), s*L.yz()*L.yx() },
+        { L.zx()*L.zx(), L.zy()*L.zy(), L.zz()*L.zz(), s*L.zx()*L.zy(), s*L.zy()*L.zz(), s*L.zz()*L.zx() },
+        { s*L.xx()*L.yx(), s*L.xy()*L.yy(), s*L.xz()*L.yz(), (L.xy()*L.yx()+L.xx()*L.yy()), (L.xz()*L.yy()+L.xy()*L.yz()), (L.xx()*L.yz()+L.xz()*L.yx()) },
+        { s*L.yx()*L.zx(), s*L.yy()*L.zy(), s*L.yz()*L.zz(), (L.yy()*L.zx()+L.yx()*L.zy()), (L.yz()*L.zy()+L.yy()*L.zz()), (L.yx()*L.zz()+L.yz()*L.zx()) },
+        { s*L.zx()*L.xx(), s*L.zy()*L.xy(), s*L.zz()*L.xz(), (L.zy()*L.xx()+L.zx()*L.xy()), (L.zz()*L.xy()+L.zy()*L.xz()), (L.zx()*L.xz()+L.zz()*L.xx()) }
     };
 
-  return symmTensor4thOrder
+    return symmTensor4thOrder
     (
-     // xxxx
-     A[0][0] * C.xxxx() * A[0][0] +
-     A[1][0] * C.xxyy() * A[0][0] +
-     A[2][0] * C.xxzz() * A[0][0] +
-     A[0][0] * C.xxyy() * A[1][0] +
-     A[1][0] * C.yyyy() * A[1][0] +
-     A[2][0] * C.yyzz() * A[1][0] +
-     A[0][0] * C.xxzz() * A[2][0] +
-     A[1][0] * C.yyzz() * A[2][0] +
-     A[2][0] * C.zzzz() * A[2][0] +
-     A[3][0] * C.xyxy() * A[3][0] +
-     A[4][0] * C.yzyz() * A[4][0] +
-     A[5][0] * C.zxzx() * A[5][0],
+        // xxxx
+        A[0][0] * C.xxxx() * A[0][0] +
+        A[1][0] * C.xxyy() * A[0][0] +
+        A[2][0] * C.xxzz() * A[0][0] +
+        A[0][0] * C.xxyy() * A[1][0] +
+        A[1][0] * C.yyyy() * A[1][0] +
+        A[2][0] * C.yyzz() * A[1][0] +
+        A[0][0] * C.xxzz() * A[2][0] +
+        A[1][0] * C.yyzz() * A[2][0] +
+        A[2][0] * C.zzzz() * A[2][0] +
+        A[3][0] * C.xyxy() * A[3][0] +
+        A[4][0] * C.yzyz() * A[4][0] +
+        A[5][0] * C.zxzx() * A[5][0],
 
-     // xxyy
-     A[0][0] * C.xxxx() * A[0][1] +
-     A[1][0] * C.xxyy() * A[0][1] +
-     A[2][0] * C.xxzz() * A[0][1] +
-     A[0][0] * C.xxyy() * A[1][1] +
-     A[1][0] * C.yyyy() * A[1][1] +
-     A[2][0] * C.yyzz() * A[1][1] +
-     A[0][0] * C.xxzz() * A[2][1] +
-     A[1][0] * C.yyzz() * A[2][1] +
-     A[2][0] * C.zzzz() * A[2][1] +
-     A[3][0] * C.xyxy() * A[3][1] +
-     A[4][0] * C.yzyz() * A[4][1] +
-     A[5][0] * C.zxzx() * A[5][1],
+        // xxyy
+        A[0][0] * C.xxxx() * A[0][1] +
+        A[1][0] * C.xxyy() * A[0][1] +
+        A[2][0] * C.xxzz() * A[0][1] +
+        A[0][0] * C.xxyy() * A[1][1] +
+        A[1][0] * C.yyyy() * A[1][1] +
+        A[2][0] * C.yyzz() * A[1][1] +
+        A[0][0] * C.xxzz() * A[2][1] +
+        A[1][0] * C.yyzz() * A[2][1] +
+        A[2][0] * C.zzzz() * A[2][1] +
+        A[3][0] * C.xyxy() * A[3][1] +
+        A[4][0] * C.yzyz() * A[4][1] +
+        A[5][0] * C.zxzx() * A[5][1],
 
-     // xzz
-     A[0][0] * C.xxxx() * A[0][2] +
-     A[1][0] * C.xxyy() * A[0][2] +
-     A[2][0] * C.xxzz() * A[0][2] +
-     A[0][0] * C.xxyy() * A[1][2] +
-     A[1][0] * C.yyyy() * A[1][2] +
-     A[2][0] * C.yyzz() * A[1][2] +
-     A[0][0] * C.xxzz() * A[2][2] +
-     A[1][0] * C.yyzz() * A[2][2] +
-     A[2][0] * C.zzzz() * A[2][2] +
-     A[3][0] * C.xyxy() * A[3][2] +
-     A[4][0] * C.yzyz() * A[4][2] +
-     A[5][0] * C.zxzx() * A[5][2],
+        // xzz
+        A[0][0] * C.xxxx() * A[0][2] +
+        A[1][0] * C.xxyy() * A[0][2] +
+        A[2][0] * C.xxzz() * A[0][2] +
+        A[0][0] * C.xxyy() * A[1][2] +
+        A[1][0] * C.yyyy() * A[1][2] +
+        A[2][0] * C.yyzz() * A[1][2] +
+        A[0][0] * C.xxzz() * A[2][2] +
+        A[1][0] * C.yyzz() * A[2][2] +
+        A[2][0] * C.zzzz() * A[2][2] +
+        A[3][0] * C.xyxy() * A[3][2] +
+        A[4][0] * C.yzyz() * A[4][2] +
+        A[5][0] * C.zxzx() * A[5][2],
 
-     // yyyy
-     A[0][1] * C.xxxx() * A[0][1] +
-     A[1][1] * C.xxyy() * A[0][1] +
-     A[2][1] * C.xxzz() * A[0][1] +
-     A[0][1] * C.xxyy() * A[1][1] +
-     A[1][1] * C.yyyy() * A[1][1] +
-     A[2][1] * C.yyzz() * A[1][1] +
-     A[0][1] * C.xxzz() * A[2][1] +
-     A[1][1] * C.yyzz() * A[2][1] +
-     A[2][1] * C.zzzz() * A[2][1] +
-     A[3][1] * C.xyxy() * A[3][1] +
-     A[4][1] * C.yzyz() * A[4][1] +
-     A[5][1] * C.zxzx() * A[5][1],
+        // yyyy
+        A[0][1] * C.xxxx() * A[0][1] +
+        A[1][1] * C.xxyy() * A[0][1] +
+        A[2][1] * C.xxzz() * A[0][1] +
+        A[0][1] * C.xxyy() * A[1][1] +
+        A[1][1] * C.yyyy() * A[1][1] +
+        A[2][1] * C.yyzz() * A[1][1] +
+        A[0][1] * C.xxzz() * A[2][1] +
+        A[1][1] * C.yyzz() * A[2][1] +
+        A[2][1] * C.zzzz() * A[2][1] +
+        A[3][1] * C.xyxy() * A[3][1] +
+        A[4][1] * C.yzyz() * A[4][1] +
+        A[5][1] * C.zxzx() * A[5][1],
 
-     // yyzz
-     A[0][1] * C.xxxx() * A[0][2] +
-     A[1][1] * C.xxyy() * A[0][2] +
-     A[2][1] * C.xxzz() * A[0][2] +
-     A[0][1] * C.xxyy() * A[1][2] +
-     A[1][1] * C.yyyy() * A[1][2] +
-     A[2][1] * C.yyzz() * A[1][2] +
-     A[0][1] * C.xxzz() * A[2][2] +
-     A[1][1] * C.yyzz() * A[2][2] +
-     A[2][1] * C.zzzz() * A[2][2] +
-     A[3][1] * C.xyxy() * A[3][2] +
-     A[4][1] * C.yzyz() * A[4][2] +
-     A[5][1] * C.zxzx() * A[5][2],
+        // yyzz
+        A[0][1] * C.xxxx() * A[0][2] +
+        A[1][1] * C.xxyy() * A[0][2] +
+        A[2][1] * C.xxzz() * A[0][2] +
+        A[0][1] * C.xxyy() * A[1][2] +
+        A[1][1] * C.yyyy() * A[1][2] +
+        A[2][1] * C.yyzz() * A[1][2] +
+        A[0][1] * C.xxzz() * A[2][2] +
+        A[1][1] * C.yyzz() * A[2][2] +
+        A[2][1] * C.zzzz() * A[2][2] +
+        A[3][1] * C.xyxy() * A[3][2] +
+        A[4][1] * C.yzyz() * A[4][2] +
+        A[5][1] * C.zxzx() * A[5][2],
 
-     // zzzz
-     A[0][2] * C.xxxx() * A[0][2] +
-     A[1][2] * C.xxyy() * A[0][2] +
-     A[2][2] * C.xxzz() * A[0][2] +
-     A[0][2] * C.xxyy() * A[1][2] +
-     A[1][2] * C.yyyy() * A[1][2] +
-     A[2][2] * C.yyzz() * A[1][2] +
-     A[0][2] * C.xxzz() * A[2][2] +
-     A[1][2] * C.yyzz() * A[2][2] +
-     A[2][2] * C.zzzz() * A[2][2] +
-     A[3][2] * C.xyxy() * A[3][2] +
-     A[4][2] * C.yzyz() * A[4][2] +
-     A[5][2] * C.zxzx() * A[5][2],
+        // zzzz
+        A[0][2] * C.xxxx() * A[0][2] +
+        A[1][2] * C.xxyy() * A[0][2] +
+        A[2][2] * C.xxzz() * A[0][2] +
+        A[0][2] * C.xxyy() * A[1][2] +
+        A[1][2] * C.yyyy() * A[1][2] +
+        A[2][2] * C.yyzz() * A[1][2] +
+        A[0][2] * C.xxzz() * A[2][2] +
+        A[1][2] * C.yyzz() * A[2][2] +
+        A[2][2] * C.zzzz() * A[2][2] +
+        A[3][2] * C.xyxy() * A[3][2] +
+        A[4][2] * C.yzyz() * A[4][2] +
+        A[5][2] * C.zxzx() * A[5][2],
 
-     // xyxy
-     A[0][3] * C.xxxx() * A[0][3] +
-     A[1][3] * C.xxyy() * A[0][3] +
-     A[2][3] * C.xxzz() * A[0][3] +
-     A[0][3] * C.xxyy() * A[1][3] +
-     A[1][3] * C.yyyy() * A[1][3] +
-     A[2][3] * C.yyzz() * A[1][3] +
-     A[0][3] * C.xxzz() * A[2][3] +
-     A[1][3] * C.yyzz() * A[2][3] +
-     A[2][3] * C.zzzz() * A[2][3] +
-     A[3][3] * C.xyxy() * A[3][3] +
-     A[4][3] * C.yzyz() * A[4][3] +
-     A[5][3] * C.zxzx() * A[5][3],
+        // xyxy
+        A[0][3] * C.xxxx() * A[0][3] +
+        A[1][3] * C.xxyy() * A[0][3] +
+        A[2][3] * C.xxzz() * A[0][3] +
+        A[0][3] * C.xxyy() * A[1][3] +
+        A[1][3] * C.yyyy() * A[1][3] +
+        A[2][3] * C.yyzz() * A[1][3] +
+        A[0][3] * C.xxzz() * A[2][3] +
+        A[1][3] * C.yyzz() * A[2][3] +
+        A[2][3] * C.zzzz() * A[2][3] +
+        A[3][3] * C.xyxy() * A[3][3] +
+        A[4][3] * C.yzyz() * A[4][3] +
+        A[5][3] * C.zxzx() * A[5][3],
 
-     // yzyz
-     A[0][4] * C.xxxx() * A[0][4] +
-     A[1][4] * C.xxyy() * A[0][4] +
-     A[2][4] * C.xxzz() * A[0][4] +
-     A[0][4] * C.xxyy() * A[1][4] +
-     A[1][4] * C.yyyy() * A[1][4] +
-     A[2][4] * C.yyzz() * A[1][4] +
-     A[0][4] * C.xxzz() * A[2][4] +
-     A[1][4] * C.yyzz() * A[2][4] +
-     A[2][4] * C.zzzz() * A[2][4] +
-     A[3][4] * C.xyxy() * A[3][4] +
-     A[4][4] * C.yzyz() * A[4][4] +
-     A[5][4] * C.zxzx() * A[5][4],
+        // yzyz
+        A[0][4] * C.xxxx() * A[0][4] +
+        A[1][4] * C.xxyy() * A[0][4] +
+        A[2][4] * C.xxzz() * A[0][4] +
+        A[0][4] * C.xxyy() * A[1][4] +
+        A[1][4] * C.yyyy() * A[1][4] +
+        A[2][4] * C.yyzz() * A[1][4] +
+        A[0][4] * C.xxzz() * A[2][4] +
+        A[1][4] * C.yyzz() * A[2][4] +
+        A[2][4] * C.zzzz() * A[2][4] +
+        A[3][4] * C.xyxy() * A[3][4] +
+        A[4][4] * C.yzyz() * A[4][4] +
+        A[5][4] * C.zxzx() * A[5][4],
 
-     // zxzx
-     A[0][5] * C.xxxx() * A[0][5] +
-     A[1][5] * C.xxyy() * A[0][5] +
-     A[2][5] * C.xxzz() * A[0][5] +
-     A[0][5] * C.xxyy() * A[1][5] +
-     A[1][5] * C.yyyy() * A[1][5] +
-     A[2][5] * C.yyzz() * A[1][5] +
-     A[0][5] * C.xxzz() * A[2][5] +
-     A[1][5] * C.yyzz() * A[2][5] +
-     A[2][5] * C.zzzz() * A[2][5] +
-     A[3][5] * C.xyxy() * A[3][5] +
-     A[4][5] * C.yzyz() * A[4][5] +
-     A[5][5] * C.zxzx() * A[5][5]
-     );
+        // zxzx
+        A[0][5] * C.xxxx() * A[0][5] +
+        A[1][5] * C.xxyy() * A[0][5] +
+        A[2][5] * C.xxzz() * A[0][5] +
+        A[0][5] * C.xxyy() * A[1][5] +
+        A[1][5] * C.yyyy() * A[1][5] +
+        A[2][5] * C.yyzz() * A[1][5] +
+        A[0][5] * C.xxzz() * A[2][5] +
+        A[1][5] * C.yyzz() * A[2][5] +
+        A[2][5] * C.zzzz() * A[2][5] +
+        A[3][5] * C.xyxy() * A[3][5] +
+        A[4][5] * C.yzyz() * A[4][5] +
+        A[5][5] * C.zxzx() * A[5][5]
+    );
 }
 
 
diff --git a/src/fvAgglomerationMethods/MGridGenGamgAgglomeration/MGridGenGAMGAgglomerate.C b/src/fvAgglomerationMethods/MGridGenGamgAgglomeration/MGridGenGAMGAgglomerate.C
index 744c02d4f..c0f1d78ab 100644
--- a/src/fvAgglomerationMethods/MGridGenGamgAgglomeration/MGridGenGAMGAgglomerate.C
+++ b/src/fvAgglomerationMethods/MGridGenGamgAgglomeration/MGridGenGAMGAgglomerate.C
@@ -146,6 +146,8 @@ Foam::tmp<Foam::labelField> Foam::MGridGenGAMGAgglomeration::agglomerate
     List<int> finalAgglom(nFineCells);
     int nMoves = -1;
 
+#   ifdef WM_DP
+
     MGridGen
     (
         nFineCells,
@@ -162,6 +164,46 @@ Foam::tmp<Foam::labelField> Foam::MGridGenGAMGAgglomeration::agglomerate
         finalAgglom.begin()
     );
 
+#   else
+
+    // Conversion of type for MGridGen interface
+
+    Field<realtype> dblVols(V.size());
+    forAll (dblVols, i)
+    {
+        dblVols[i] = V[i];
+    }
+
+    Field<realtype> dblAreas(Sb.size());
+    forAll (dblAreas, i)
+    {
+        dblAreas[i] = Sb[i];
+    }
+
+    Field<realtype> dblFaceWeights(faceWeights.size());
+    forAll (dblFaceWeights, i)
+    {
+        dblFaceWeights[i] = faceWeights[i];
+    }
+
+    MGridGen
+    (
+        nFineCells,
+        cellCellOffsets.begin(),
+        dblVols.begin(),
+        dblAreas.begin(),
+        cellCells.begin(),
+        dblFaceWeights.begin(),
+        minSize,
+        maxSize,
+        options.begin(),
+        &nMoves,
+        &nCoarseCells,
+        finalAgglom.begin()
+    );
+
+#   endif
+
     return tmp<labelField>(new labelField(finalAgglom));
 }
 
diff --git a/src/solidModels/arbitraryCrack/solidCohesive/solidCohesiveFvPatchVectorField.C b/src/solidModels/arbitraryCrack/solidCohesive/solidCohesiveFvPatchVectorField.C
index a7a78cb44..6bf305cc9 100644
--- a/src/solidModels/arbitraryCrack/solidCohesive/solidCohesiveFvPatchVectorField.C
+++ b/src/solidModels/arbitraryCrack/solidCohesive/solidCohesiveFvPatchVectorField.C
@@ -729,7 +729,7 @@ void solidCohesiveFvPatchVectorField::updateCoeffs()
                 unloadingDeltaEff_,
                 Foam::sqrt
                 (
-                    max(deltaN_, 0.0)*max(deltaN_, 0.0) + deltaS_*deltaS_
+                    sqr(max(deltaN_, scalar(0))) + sqr(deltaS_)
                 )
             );
 
diff --git a/src/solidModels/contactModels/normalContactModels/dirichletNeumann/dirichletNeumann.C b/src/solidModels/contactModels/normalContactModels/dirichletNeumann/dirichletNeumann.C
index 1703a046e..fc4e852fb 100644
--- a/src/solidModels/contactModels/normalContactModels/dirichletNeumann/dirichletNeumann.C
+++ b/src/solidModels/contactModels/normalContactModels/dirichletNeumann/dirichletNeumann.C
@@ -261,8 +261,8 @@ void dirichletNeumann::correct
        slaveDispMag =
      localSlaveInterpolator.pointToFaceInterpolate<scalar>
      (
-      min(slavePointPenetration, 0.0)
-      );
+         min(slavePointPenetration, scalar(0))
+     );
 
        // for visualisation
        slaveContactPointGap() = slavePointPenetration;
@@ -277,11 +277,14 @@ void dirichletNeumann::correct
      = mesh.boundaryMesh()[slavePatchIndex].start();
        forAll(slaveDispMag, facei)
      {
-       slaveDispMag[facei] =
-         globalSlavePenetration
-         [
-          mesh.faceZones()[slaveFaceZoneID()].whichFace(slavePatchStart + facei)
-          ];
+         slaveDispMag[facei] =
+             globalSlavePenetration
+             [
+                 mesh.faceZones()[slaveFaceZoneID()].whichFace
+                 (
+                     slavePatchStart + facei
+                 )
+             ];
 
        //- when the master surface surrounds the slave (like the pelvis
        // and femur head) then
@@ -317,7 +320,7 @@ void dirichletNeumann::correct
        slaveContactPointGap() = slavePointPenetration;
 
        // set all positive penetrations to zero
-       slaveDispMag = min(slaveDispMag,0.0);
+       slaveDispMag = min(slaveDispMag, scalar(0));
 
        minSlavePointPenetration = min(slavePointPenetration);
        reduce(minSlavePointPenetration, minOp<scalar>());