Rewrite and clean-up

tutorials/solidMechanics/elasticThermalSolidFoam/hotCylinder/analyticalHotCylinder/analyticalHotCylinder.C

@@ -45,7 +45,8 @@ int main(int argc, char *argv[])
 
     runTime++;
 
-  Info << "Writing analytical solution for a plain strain cylinder with concentric hole,\nwhere"
+    Info<< "Writing analytical solution for a plain strain cylinder "
+        << "with concentric hole,\nwhere"
         << "\n\tinner radius = 0.5"
         << "\n\touter radius = 0.7"
         << "\n\tinner temperature = 100"
@@ -68,6 +69,21 @@ int main(int argc, char *argv[])
     scalar nu = 0.3;
     scalar alpha = 1e-5;
 
+    const volVectorField& C = mesh.C();
+
+    //- radial coordinate
+    volScalarField radii
+    (
+        sqrt
+        (
+            sqr(C.component(vector::X))
+          + sqr(C.component(vector::Y))
+        )/dimensionedScalar("one", dimLength, 1)
+    );
+
+    const scalarField& rIn = radii.internalField();
+
+    Info << "Writing analytical termpature field" << endl;
     //- create T field
     volScalarField T
     (
@@ -79,52 +95,12 @@ int main(int argc, char *argv[])
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-     mesh,
+        ((Ti - To)/Foam::log(b/a))*Foam::log(b/radii)
-     dimensionedScalar("zero", dimTemperature, 0.0)
     );
-  
-  const volVectorField& C = mesh.C();
-
-  //- radial coordinate
-  volScalarField radii =
-    C.component(vector::X)*C.component(vector::X) + C.component(vector::Y)*C.component(vector::Y);
-  forAll(radii.internalField(), celli)
-    {
-      radii.internalField()[celli] = ::sqrt(radii.internalField()[celli]);
-    }
-  forAll(radii.boundaryField(), patchi)
-    {
-      forAll(radii.boundaryField()[patchi], facei)
-	{
-	  radii.boundaryField()[patchi][facei] = ::sqrt(radii.boundaryField()[patchi][facei]);
-	}
-    }
-
-  forAll(T.internalField(), celli)
-    {
-      const scalar& r = radii[celli];
-
-      T.internalField()[celli] =
-	( (Ti-To)/Foam::log(b/a) ) * Foam::log(b/r);
-    }
-
-  forAll(T.boundaryField(), patchi)
-    {
-      forAll(T.boundaryField()[patchi], facei)
-	{
-	  const scalar& r = radii.boundaryField()[patchi][facei];
-
-	  T.boundaryField()[patchi][facei] =
-	    ( (Ti-To)/Foam::log(b/a) ) * Foam::log(b/r);
-	}
-    }
-
-  //- write temperature file
-  Info << "Writing analytical termpature field" << endl;
     T.write();
 
-
     //- create sigma field
+    Info << "\nWriting analytical sigmaR field" << endl;
     volScalarField sigmaR
     (
         IOobject
@@ -135,36 +111,16 @@ int main(int argc, char *argv[])
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-     mesh,
+        ((alpha*E*(Ti - To))/(2*(1 - nu)*Foam::log(b/a)))*
-     dimensionedScalar("zero", dimForce/dimArea, 0.0)
+        (
+            -Foam::log(b/radii)
+          - (sqr(a)/(sqr(b) - sqr(a)))*(1 - sqr(b)/sqr(radii))*Foam::log(b/a)
+        )
     );
-
-  forAll(sigmaR.internalField(), celli)
-    {
-      const scalar& r = radii.internalField()[celli];
-
-      sigmaR.internalField()[celli] =
-	( (alpha*E*(Ti-To))/(2*(1-nu)*Foam::log(b/a)) ) *
-	(-Foam::log(b/r) -( a*a/(b*b - a*a))*(1 - (b*b)/(r*r))*Foam::log(b/a));
-    }
-
-  forAll(sigmaR.boundaryField(), patchi)
-    {
-      forAll(sigmaR.boundaryField()[patchi], facei)
-	{
-	  const scalar& r = radii.boundaryField()[patchi][facei];
-	  
-	  sigmaR.boundaryField()[patchi][facei] =
-	    ( (alpha*E*(Ti-To))/(2*(1-nu)*Foam::log(b/a)) ) *
-	    ( -Foam::log(b/r) - ( a*a/(b*b - a*a))*(1 - (b*b)/(r*r))*Foam::log(b/a) );
-	}
-    }
-
-  //- write temperature file
-  Info << "\nWriting analytical sigmaR field" << endl;
     sigmaR.write();
 
 
+    Info << "\nWriting analytical sigmaTheta field" << endl;
     volScalarField sigmaTheta
     (
         IOobject
@@ -175,36 +131,15 @@ int main(int argc, char *argv[])
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-     mesh,
+        ((alpha*E*(Ti - To))/(2*(1 - nu)*Foam::log(b/a)))*
-     dimensionedScalar("zero", dimForce/dimArea, 0.0)
+        (
+            1 - Foam::log(b/radii)
+          - (sqr(a)/(sqr(b) - sqr(a)))*(1 + sqr(b)/sqr(radii))*Foam::log(b/a)
+        )
     );
-
-  forAll(sigmaTheta.internalField(), celli)
-    {
-      const scalar& r = radii.internalField()[celli];
-
-      sigmaTheta.internalField()[celli] =
-	( (alpha*E*(Ti-To))/(2*(1-nu)*Foam::log(b/a)) ) *
-	(1 -Foam::log(b/r) - ( a*a/(b*b - a*a))*(1 + (b*b)/(r*r))*Foam::log(b/a) );
-    }
-
-  forAll(sigmaTheta.boundaryField(), patchi)
-    {
-      forAll(sigmaTheta.boundaryField()[patchi], facei)
-	{
-	  const scalar& r = radii.boundaryField()[patchi][facei];
-	  
-	  sigmaTheta.boundaryField()[patchi][facei] =
-	    ( (alpha*E*(Ti-To))/(2*(1-nu)*Foam::log(b/a)) ) *
-	    (1 -Foam::log(b/r) - ( a*a/(b*b - a*a))*(1 + (b*b)/(r*r))*Foam::log(b/a) );
-	}
-    }
-
-
-  //- write temperature file
-  Info << "\nWriting analytical sigmaTheta field" << endl;
     sigmaTheta.write();
 
+    Info << "\nWriting analytical sigmaZ field" << endl;
     volScalarField sigmaZ
     (
         IOobject
@@ -215,49 +150,34 @@ int main(int argc, char *argv[])
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-     mesh,
+        // Timoshenko says this but I am not sure I am not sure the BCs in
-     dimensionedScalar("zero", dimForce/dimArea, 0.0)
+        // the z direction
+        // ((alpha*E*(Ti - To))/(2*(1 - nu)*Foam::log(b/a)))*
+        // (1 - 2*Foam::log(b/radii) - ( 2*sqr(a)/(sqr(b) - sqr(a)))*Foam::log(b/a));
+        0.3*(sigmaR + sigmaTheta) - E*alpha*(T)
     );
-
-  forAll(sigmaZ.internalField(), celli)
-    {
-      //- Timoshenko says this but I am not sure I am not sure the BCs in
-      //- the z direction
-      // sigmaZ.internalField()[celli] =
-      // 	( (alpha*E*(Ti-To))/(2*(1-nu)*Foam::log(b/a)) ) *
-      // 	(1 - 2*Foam::log(b/r) - ( 2*a*a/(b*b - a*a))*Foam::log(b/a));
-
-      sigmaZ.internalField()[celli] = 
-	0.3*(sigmaR.internalField()[celli] + sigmaTheta.internalField()[celli])
-	- E*alpha*(T.internalField()[celli]);
-    }
-
-  forAll(sigmaZ.boundaryField(), patchi)
-    {
-      forAll(sigmaZ.boundaryField()[patchi], facei)
-	{
-	  //- Timoshenko says this but I am not sure I am not sure the BCs in
-	  //- the z direction
-	  //sigmaZ.boundaryField()[patchi][facei] =
-	  //( (alpha*E*(Ti-To))/(2*(1-nu)*Foam::log(b/a)) ) *
-	  //(1 - 2*Foam::log(b/r) - ( 2*a*a/(b*b - a*a))*Foam::log(b/a));
-
-	  //-for general 2-D plain strain problems, the axial stress is given by this:
-	  sigmaZ.boundaryField()[patchi][facei] = 
-	    nu*(sigmaR.boundaryField()[patchi][facei] + sigmaTheta.boundaryField()[patchi][facei])
-	    - E*alpha*(T.boundaryField()[patchi][facei]);
-	}
-    }
-
-
-  //- write temperature file
-  Info << "\nWriting analytical sigmaZ field" << endl;
     sigmaZ.write();
 
-
-  //- create analytical sigma tensor
-
     //- create theta field
+    volScalarField yOverX
+    (
+        "yOverX",
+        Foam::max
+        (
+            scalar(-1),
+            Foam::min
+            (
+                scalar(1),
+                mesh.C().component(vector::Y)/
+                stabilise
+                (
+                    mesh.C().component(vector::X),
+                    dimensionedScalar("small", dimLength, SMALL)
+                )
+            )
+        )
+    );
+
     volScalarField theta
     (
         IOobject
@@ -268,27 +188,8 @@ int main(int argc, char *argv[])
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
-     mesh,
+        Foam::atan(yOverX)
-     dimensionedScalar("zero", dimless, 0.0)
     );
-  forAll(theta.internalField(), celli)
-    {
-      const scalar& x = mesh.C().internalField()[celli][vector::X];
-      const scalar& y = mesh.C().internalField()[celli][vector::Y];
-
-      theta.internalField()[celli] = Foam::atan(y/x);
-    }
-
-  forAll(theta.boundaryField(), patchi)
-    {
-      forAll(theta.boundaryField()[patchi], facei)
-	{
-	  const scalar& x = mesh.C().boundaryField()[patchi][facei][vector::X];
-	  const scalar& y = mesh.C().boundaryField()[patchi][facei][vector::Y];
-	  
-	  theta.boundaryField()[patchi][facei] = Foam::atan(y/x);
-	}
-    }
 
     //- rotation matrix to convert polar stresses to cartesian
     volTensorField rotMat
@@ -305,24 +206,36 @@ int main(int argc, char *argv[])
         dimensionedTensor("zero", dimless, tensor::zero)
     );
 
-  forAll(rotMat.internalField(), celli)
+    tensorField& rotMatIn = rotMat.internalField();
-    {
+    const scalarField tIn = theta.internalField();
-      const scalar& t = theta.internalField()[celli];
 
-      rotMat.internalField()[celli] = tensor(::cos(t), ::sin(t), 0,
+    forAll (rotMatIn, celli)
-					     -::sin(t), ::cos(t), 0,
+    {
-					     0, 0, 1);
+        const scalar& t = tIn[celli];
+
+        rotMatIn[celli] =
+            tensor
+            (
+                Foam::cos(t),  Foam::sin(t), 0,
+               -Foam::sin(t),  Foam::cos(t), 0,
+                0, 0, 1
+            );
     }
 
-  forAll(rotMat.boundaryField(), patchi)
+
+    forAll (rotMat.boundaryField(), patchi)
     {
-      forAll(rotMat.boundaryField()[patchi], facei)
+        forAll (rotMat.boundaryField()[patchi], facei)
         {
             const scalar& t = theta.boundaryField()[patchi][facei];
 
-	  rotMat.boundaryField()[patchi][facei] =  tensor(::cos(t), ::sin(t), 0,
+            rotMat.boundaryField()[patchi][facei] =
-					     -::sin(t), ::cos(t), 0,
+                tensor
-					     0, 0, 1);
+                (
+                    Foam::cos(t),  Foam::sin(t), 0,
+                   -Foam::sin(t),  Foam::cos(t), 0,
+                    0, 0, 1
+                );
         }
     }
 
@@ -340,47 +253,58 @@ int main(int argc, char *argv[])
         dimensionedSymmTensor("zero", dimForce/dimArea, symmTensor::zero)
     );
 
-  forAll(sigma.internalField(), celli)
     {
-      const scalar& r = sigmaR.internalField()[celli];
+        symmTensorField& sigmaIn = sigma.internalField();
-      const scalar& t = sigmaTheta.internalField()[celli];
-      const scalar& z = sigmaZ.internalField()[celli];
 
-      const tensor& rot = rotMat.internalField()[celli];
+        const scalarField& rIn = sigmaR.internalField();
+        const scalarField& tIn = sigmaTheta.internalField();
+        const scalarField& zIn = sigmaZ.internalField();
 
-      symmTensor sigmaCart(r, 0, 0,
+        forAll (sigmaIn, celli)
-			      t, 0,
-			         z);
-
-      sigma.internalField()[celli] =
-	    symm(rot.T() & sigmaCart & rot);
-
-      //-for general 2-D plain strain problems, the axial stress is given by this:
-      //- (which is not equal to the solution by Timoshenko... hmmmnn)
-//       sigma.internalField()[celli][symmTensor::ZZ] = 
-// 	0.3*(sigma.internalField()[celli][symmTensor::XX] + sigma.internalField()[celli][symmTensor::YY])
-// 	- E*alpha*(T.internalField()[celli]);
-    }
-
-  forAll(sigma.boundaryField(), patchi)
         {
-      forAll(sigma.boundaryField()[patchi], facei)
+            symmTensor sigmaCart
-	{
+            (
-	  const scalar& r = sigmaR.boundaryField()[patchi][facei];
+                rIn[celli], 0, 0,
-	  const scalar& t = sigmaTheta.boundaryField()[patchi][facei];
+                tIn[celli], 0,
-	  const scalar& z = sigmaZ.boundaryField()[patchi][facei];
+                zIn[celli]
+            );
 
-	  const tensor& rot = rotMat.boundaryField()[patchi][facei];
+            const tensor& rot = rotMatIn[celli];
 
-	  symmTensor sigmaCart(r, 0, 0,
+            sigmaIn[celli] = symm(rot.T() & sigmaCart & rot);
-		                  t, 0,
+
-			             z);
+            // for general 2-D plain strain problems, the axial stress is:
-	  sigma.boundaryField()[patchi][facei] =
+            // (which is not equal to the solution by Timoshenko... hmmmnn)
-	    symm(rot.T() & sigmaCart & rot);
+            //       sigmaIn[celli][symmTensor::ZZ] =
+            //           0.3*(sigmaIn[celli][symmTensor::XX]
+            //         + sigmaIn[celli][symmTensor::YY])
+            //         - E*alpha*(T.internalField()[celli]);
         }
     }
 
+    forAll (sigma.boundaryField(), patchi)
+    {
+        symmTensorField& pSigma = sigma.boundaryField()[patchi];
+        const scalarField& pR = sigmaR.boundaryField()[patchi];
+        const scalarField& pT = sigmaTheta.boundaryField()[patchi];
+        const scalarField& pZ = sigmaZ.boundaryField()[patchi];
 
+        const tensorField pRot = rotMat.boundaryField()[patchi];
+
+        forAll (pSigma, facei)
+        {
+            const tensor& rot = pRot[facei];
+
+            symmTensor sigmaCart
+            (
+                pR[facei], 0, 0,
+                pT[facei], 0,
+                pZ[facei]
+            );
+
+            pSigma[facei] = symm(rot.T() & sigmaCart & rot);
+        }
+    }
 
     Info << "\nWriting analytical sigma tensor" << endl;
     sigma.write();