Update step-18 to use rotation matrix defined in physics module.

doc/news/changes/minor/20170109Jean-PaulPelteret_1

@@ -0,0 +1,6 @@
+Fixed: A bug in the rotation matrix used in step-18 has now been corrected.
+This tutorial now uses the function
+Physics::Transformations::Rotations::rotation_matrix_3d to compute the 3d 
+rotation matrix.
+<br>
+(Jean-Paul Pelteret, Paul Kuberry, 2017/01/09)

examples/step-18/step-18.cc

@@ -1,6 +1,6 @@
 /* ---------------------------------------------------------------------
  *
- * Copyright (C) 2000 - 2016 by the deal.II authors
+ * Copyright (C) 2000 - 2017 by the deal.II authors
  *
  * This file is part of the deal.II library.
  *
@@ -56,7 +56,7 @@
 #include <deal.II/numerics/data_out.h>
 #include <deal.II/numerics/error_estimator.h>
 
-// And here the only two new things among the header files: an include file in
+// And here the only three new things among the header files: an include file in
 // which symmetric tensors of rank 2 and 4 are implemented, as introduced in
 // the introduction:
 #include <deal.II/base/symmetric_tensor.h>
@@ -66,6 +66,11 @@
 // owned by the present process in a %parallel program:
 #include <deal.II/grid/filtered_iterator.h>
 
+// And lastly a header that contains some functions that will help us compute
+// rotaton matrices of the local coordinate systems at specific points in the
+// domain.
+#include <deal.II/physics/transformations.h>
+
 // This is then simply C++ again:
 #include <fstream>
 #include <iostream>
@@ -280,11 +285,12 @@ namespace Step18
     // From this, compute the angle of rotation:
     const double angle = std::atan (curl);
 
-    // And from this, build the antisymmetric rotation matrix:
-    const double t[2][2] = {{ cos(angle), sin(angle) },
-      {-sin(angle), cos(angle) }
-    };
-    return Tensor<2,2>(t);
+    // And from this, build the antisymmetric rotation matrix. We want this
+    // rotation matrix to represent the rotation of the local coordinate system
+    // with respect to the global Cartesian basis, to we construct it with a
+    // negative angle. The rotation matrix therefore represents the rotation
+    // required to move from the local to the global coordinate system.
+    return Physics::Transformations::Rotations::rotation_matrix_2d(-angle);
   }
 
 
@@ -299,7 +305,8 @@ namespace Step18
                          grad_u[1][0] - grad_u[0][1]);
 
     // From this vector, using its magnitude, compute the tangent of the angle
-    // of rotation, and from it the actual angle:
+    // of rotation, and from it the actual angle of rotation with respect to
+    // the Cartesian basis:
     const double tan_angle = std::sqrt(curl*curl);
     const double angle = std::atan (tan_angle);
 
@@ -313,44 +320,19 @@ namespace Step18
     // into trouble when dividing doing so. Therefore, let's shortcut this and
     // simply return the identity matrix if the angle of rotation is really
     // small:
-    if (angle < 1e-9)
+    if (std::abs(angle) < 1e-9)
       {
         static const double rotation[3][3]
         = {{ 1, 0, 0}, { 0, 1, 0 }, { 0, 0, 1 } };
         static const Tensor<2,3> rot(rotation);
         return rot;
       }
 
-    // Otherwise compute the real rotation matrix. The algorithm for this is
-    // not exactly obvious, but can be found in a number of books,
-    // particularly on computer games where rotation is a very frequent
-    // operation. Online, you can find a description at
-    // http://www.makegames.com/3drotation/ and (this particular form, with
-    // the signs as here) at
-    // http://www.gamedev.net/reference/articles/article1199.asp:
-    const double c = std::cos(angle);
-    const double s = std::sin(angle);
-    const double t = 1-c;
-
+    // Otherwise compute the real rotation matrix. For this, again we rely on
+    // a predefined function to compute the rotation matrix of the local
+    // coordinate system.
     const Point<3> axis = curl/tan_angle;
-    const double rotation[3][3]
-    = {{
-        t *axis[0] *axis[0]+c,
-        t *axis[0] *axis[1]+s *axis[2],
-        t *axis[0] *axis[2]-s *axis[1]
-      },
-      {
-        t *axis[0] *axis[1]-s *axis[2],
-        t *axis[1] *axis[1]+c,
-        t *axis[1] *axis[2]+s *axis[0]
-      },
-      {
-        t *axis[0] *axis[2]+s *axis[1],
-        t *axis[1] *axis[1]-s *axis[0],
-        t *axis[2] *axis[2]+c
-      }
-    };
-    return Tensor<2,3>(rotation);
+    return Physics::Transformations::Rotations::rotation_matrix_3d(axis, -angle);
   }