From 5f6e325d4fce8c497244c6603cbb117a1f9cd6a2 Mon Sep 17 00:00:00 2001 From: Jean-Paul Pelteret Date: Mon, 9 Jan 2017 19:15:08 +0100 Subject: [PATCH] Update step-18 to use rotation matrix defined in physics module. --- .../changes/minor/20170109Jean-PaulPelteret_1 | 6 ++ examples/step-18/step-18.cc | 58 +++++++------------ 2 files changed, 26 insertions(+), 38 deletions(-) create mode 100644 doc/news/changes/minor/20170109Jean-PaulPelteret_1 diff --git a/doc/news/changes/minor/20170109Jean-PaulPelteret_1 b/doc/news/changes/minor/20170109Jean-PaulPelteret_1 new file mode 100644 index 0000000000..b5d8e58aa3 --- /dev/null +++ b/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. +
+(Jean-Paul Pelteret, Paul Kuberry, 2017/01/09) diff --git a/examples/step-18/step-18.cc b/examples/step-18/step-18.cc index a3fae45e1b..4f2ba2dbc4 100644 --- a/examples/step-18/step-18.cc +++ b/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 #include -// 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 @@ -66,6 +66,11 @@ // owned by the present process in a %parallel program: #include +// 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 + // This is then simply C++ again: #include #include @@ -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,7 +320,7 @@ 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 } }; @@ -321,36 +328,11 @@ namespace Step18 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); } -- 2.39.5