Update and unify interface for `Physics` and `GridTools` rotations.

diff --git a/doc/news/changes/incompatibilities/20211011Fehling b/doc/news/changes/incompatibilities/20211011Fehling
index 7105a29a87dc322a3dddf723b98d4bc45f9e266a..6f18ea6fe4902c240d93a5d245eabd227a6b005d 100644 (file)
--- a/doc/news/changes/incompatibilities/20211011Fehling
+++ b/doc/news/changes/incompatibilities/20211011Fehling
@@ -1,5 +1,7 @@
 Changed: The 3D implementation of GridTools::rotate() with an integer
 for a Cartesian coordinate direction has been superseded by the version
-that accepts unit vectors as rotation axes.
+that accepts unit vectors as rotation axes. Further,
+Physics::Transformations::Rotations::rotation_matrix_3d() now requires
+a Tensor<1,3> object instead of a Point<3> as an axis.
 <br>
 (Marc Fehling, 2021/10/11)

diff --git a/examples/step-18/step-18.cc b/examples/step-18/step-18.cc
index f0906082ac0a43457b395b6ae1e0e618eaef2391..de93a76b2cfd454a1e7a00537037172b0989d396 100644 (file)
--- a/examples/step-18/step-18.cc
+++ b/examples/step-18/step-18.cc
@@ -287,9 +287,9 @@ namespace Step18
   {
     // Again first compute the curl of the velocity field. This time, it is a
     // real vector:
-    const Point<3> curl(grad_u[2][1] - grad_u[1][2],
-                        grad_u[0][2] - grad_u[2][0],
-                        grad_u[1][0] - grad_u[0][1]);
+    const Tensor<1, 3> curl({grad_u[2][1] - grad_u[1][2],
+                             grad_u[0][2] - grad_u[2][0],
+                             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 with respect to
@@ -317,7 +317,7 @@ namespace Step18
     // 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 Tensor<1, 3> axis = curl / tan_angle;
     return Physics::Transformations::Rotations::rotation_matrix_3d(axis,
                                                                    -angle);
   }

diff --git a/include/deal.II/base/symmetric_tensor.templates.h b/include/deal.II/base/symmetric_tensor.templates.h
index bc2d57376c5eb6bb688abf1eeb9938d52b917c99..80bb1252f4c9dc7846ef346c111675b946378a5b 100644 (file)
--- a/include/deal.II/base/symmetric_tensor.templates.h
+++ b/include/deal.II/base/symmetric_tensor.templates.h
@@ -759,15 +759,15 @@ namespace internal
         {
           case (0):
             R = dealii::Physics::Transformations::Rotations::rotation_matrix_3d(
-              {1, 0, 0}, rotation_angle);
+              Tensor<1, 3>({1., 0., 0.}), rotation_angle);
             break;
           case (1):
             R = dealii::Physics::Transformations::Rotations::rotation_matrix_3d(
-              {0, 1, 0}, rotation_angle);
+              Tensor<1, 3>({0., 1., 0.}), rotation_angle);
             break;
           case (2):
             R = dealii::Physics::Transformations::Rotations::rotation_matrix_3d(
-              {0, 0, 1}, rotation_angle);
+              Tensor<1, 3>({0., 0., 1.}), rotation_angle);
             break;
           default:
             AssertThrow(false, ExcNotImplemented());

diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h
index d43eb0435b332cc4e72564efc4838204b3fa06d1..e2d62c7977845fa253ac3eb51383d6839fd66307 100644 (file)
--- a/include/deal.II/grid/grid_tools.h
+++ b/include/deal.II/grid/grid_tools.h
@@ -569,9 +569,9 @@ namespace GridTools
    */
   template <int dim>
   void
-  rotate(const double            angle,
-         const Point<3, double> &axis,
-         Triangulation<dim, 3> & triangulation);
+  rotate(const Tensor<1, 3, double> &axis,
+         const double                angle,
+         Triangulation<dim, 3> &     triangulation);
 
   /**
    * Rotate all vertices of the given @p triangulation in counter-clockwise

diff --git a/include/deal.II/physics/transformations.h b/include/deal.II/physics/transformations.h
index e8bcb8a6955a0848f8ea96d2a67be86aaa634f3a..268a162f979008fe2e7253483c9ddf4b13aa4f6f 100644 (file)
--- a/include/deal.II/physics/transformations.h
+++ b/include/deal.II/physics/transformations.h
@@ -18,7 +18,6 @@
 
 #include <deal.II/base/config.h>
 
-#include <deal.II/base/point.h>
 #include <deal.II/base/symmetric_tensor.h>
 #include <deal.II/base/tensor.h>
 
@@ -90,6 +89,15 @@ namespace Physics
        */
       template <typename Number>
       Tensor<2, 3, Number>
+      rotation_matrix_3d(const Tensor<1, 3, Number> &axis, const Number &angle);
+
+      /**
+       * @copydoc Physics::Transformations::Rotations::rotation_matrix_3d()
+       *
+       * @deprecated Use the variant with a Tensor as an axis.
+       */
+      template <typename Number>
+      DEAL_II_DEPRECATED_EARLY Tensor<2, 3, Number>
       rotation_matrix_3d(const Point<3, Number> &axis, const Number &angle);
 
       //@}
@@ -914,8 +922,8 @@ Physics::Transformations::Rotations::rotation_matrix_2d(const Number &angle)
 template <typename Number>
 Tensor<2, 3, Number>
 Physics::Transformations::Rotations::rotation_matrix_3d(
-  const Point<3, Number> &axis,
-  const Number &          angle)
+  const Tensor<1, 3, Number> &axis,
+  const Number &              angle)
 {
   Assert(std::abs(axis.norm() - 1.0) < 1e-9,
          ExcMessage("The supplied axial vector is not a unit vector."));
@@ -936,6 +944,17 @@ Physics::Transformations::Rotations::rotation_matrix_3d(
 
 
 
+template <typename Number>
+Tensor<2, 3, Number>
+Physics::Transformations::Rotations::rotation_matrix_3d(
+  const Point<3, Number> &axis,
+  const Number &          angle)
+{
+  return rotation_matrix_3d(static_cast<Tensor<1, 3, Number>>(axis), angle);
+}
+
+
+
 template <int dim, typename Number>
 inline Tensor<1, dim, Number>
 Physics::Transformations::Contravariant::push_forward(

diff --git a/source/grid/grid_generator.cc b/source/grid/grid_generator.cc
index ff9868ae4411c9c78fe48ab468434470fb778f39..f7d6f104a26436a5cb7edd6a69d4e658722ea430 100644 (file)
--- a/source/grid/grid_generator.cc
+++ b/source/grid/grid_generator.cc
@@ -4940,7 +4940,7 @@ namespace GridGenerator
                                          n_slices,
                                          2 * half_length,
                                          triangulation);
-    GridTools::rotate(numbers::PI_2, Point<3>({0., 1., 0.}), triangulation);
+    GridTools::rotate(Tensor<1, 3>({0., 1., 0.}), numbers::PI_2, triangulation);
     GridTools::shift(Tensor<1, 3>({-half_length, 0.0, 0.0}), triangulation);
     // At this point we have a cylinder. Multiply the y and z coordinates by a
     // factor that scales (with x) linearly between radius_0 and radius_1 to fix

diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc
index 9abd1b773362184e30b2b386fb66220e4003af1a..ba6e7f2a0b9f09f9cb0bf0238ddbdac45fb619b7 100644 (file)
--- a/source/grid/grid_tools.cc
+++ b/source/grid/grid_tools.cc
@@ -1972,7 +1972,7 @@ namespace GridTools
     class Rotate3d
     {
     public:
-      Rotate3d(const double angle, const Point<3, double> &axis)
+      Rotate3d(const Tensor<1, 3, double> &axis, const double angle)
         : rotation_matrix(
             Physics::Transformations::Rotations::rotation_matrix_3d(axis,
                                                                     angle))
@@ -2019,11 +2019,11 @@ namespace GridTools
 
   template <int dim>
   void
-  rotate(const double            angle,
-         const Point<3, double> &axis,
-         Triangulation<dim, 3> & triangulation)
+  rotate(const Tensor<1, 3, double> &axis,
+         const double                angle,
+         Triangulation<dim, 3> &     triangulation)
   {
-    transform(internal::Rotate3d(angle, axis), triangulation);
+    transform(internal::Rotate3d(axis, angle), triangulation);
   }
 
 
@@ -2035,10 +2035,10 @@ namespace GridTools
   {
     Assert(axis < 3, ExcMessage("Invalid axis given!"));
 
-    Point<3, double> vector;
+    Tensor<1, 3, double> vector;
     vector[axis] = 1.;
 
-    transform(internal::Rotate3d(angle, vector), triangulation);
+    transform(internal::Rotate3d(vector, angle), triangulation);
   }
 
 

diff --git a/source/grid/grid_tools.inst.in b/source/grid/grid_tools.inst.in
index 193bfe7dfc705b6a2c594048c83d5b7398128cf8..170164eae3943b7ea3201d9974bf7f2883ae7010 100644 (file)
--- a/source/grid/grid_tools.inst.in
+++ b/source/grid/grid_tools.inst.in
@@ -273,8 +273,8 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
 #  if deal_II_space_dimension == 3
       template void
       rotate<deal_II_dimension>(
+        const Tensor<1, deal_II_space_dimension, double> &,
         const double,
-        const Point<deal_II_space_dimension, double> &,
         Triangulation<deal_II_dimension, deal_II_space_dimension> &);
 
       template void

diff --git a/tests/grid/rotate_01.cc b/tests/grid/rotate_01.cc
index a49faef88c76d2718335dfe381fbfd9101d6252b..9f2e08c0b027a9aee096606497e08f74dd50a297 100644 (file)
--- a/tests/grid/rotate_01.cc
+++ b/tests/grid/rotate_01.cc
@@ -65,11 +65,11 @@ test()
   GridGenerator::subdivided_parallelepiped<dim, spacedim>(tria, origin, edges);
 
   // GridOut().write_gnuplot (tria, deallog.get_file_stream());
-  GridTools::rotate(numbers::PI / 3.0, Point<3>({1., 0., 0.}), tria);
+  GridTools::rotate(Tensor<1, 3>({1., 0., 0.}), numbers::PI / 3.0, tria);
   // GridOut().write_gnuplot (tria, deallog.get_file_stream());
-  GridTools::rotate(numbers::PI / 10.0, Point<3>({0., 1., 0.}), tria);
+  GridTools::rotate(Tensor<1, 3>({0., 1., 0.}), numbers::PI / 10.0, tria);
   // GridOut().write_gnuplot (tria, deallog.get_file_stream());
-  GridTools::rotate(-numbers::PI / 5.0, Point<3>({0., 0., 1.}), tria);
+  GridTools::rotate(Tensor<1, 3>({0., 0., 1.}), -numbers::PI / 5.0, tria);
   GridOut().write_gnuplot(tria, deallog.get_file_stream());
 }
 

diff --git a/tests/physics/step-18-rotation_matrix.cc b/tests/physics/step-18-rotation_matrix.cc
index a6e79f794d37a0584e9c2e316fa4265984f03792..0e14c7ca95f8955c115f78d6f74c66bdc8a26df8 100644 (file)
--- a/tests/physics/step-18-rotation_matrix.cc
+++ b/tests/physics/step-18-rotation_matrix.cc
@@ -128,14 +128,14 @@ namespace Step18
   Tensor<2, 3>
   get_rotation_matrix(const std::vector<Tensor<1, 3>> &grad_u)
   {
-    const Point<3> curl(grad_u[2][1] - grad_u[1][2],
-                        grad_u[0][2] - grad_u[2][0],
-                        grad_u[1][0] - grad_u[0][1]);
-    const double   tan_angle = std::sqrt(curl * curl);
+    const Tensor<1, 3> curl({grad_u[2][1] - grad_u[1][2],
+                             grad_u[0][2] - grad_u[2][0],
+                             grad_u[1][0] - grad_u[0][1]});
+    const double       tan_angle = std::sqrt(curl * curl);
     // Note: Here the negative angle suggests that we're computing the rotation
     // of the coordinate system around a fixed point
-    const double   angle = -std::atan(tan_angle);
-    const Point<3> axis  = curl / tan_angle;
+    const double       angle = -std::atan(tan_angle);
+    const Tensor<1, 3> axis  = curl / tan_angle;
     return Physics::Transformations::Rotations::rotation_matrix_3d(axis, angle);
   }
   template <int dim>

diff --git a/tests/physics/transformations-rotations_01.cc b/tests/physics/transformations-rotations_01.cc
index 606422e1f3d93712608580b1e4a82f66164aa44a..84554d454134099fe2fb0dc3fea1d880fd597881 100644 (file)
--- a/tests/physics/transformations-rotations_01.cc
+++ b/tests/physics/transformations-rotations_01.cc
@@ -42,7 +42,8 @@ test_rotation_matrix_3d_z_axis(const double angle)
   Assert(std::abs(determinant(R_z) - 1.0) < 1e-9,
          ExcMessage("Rodrigues rotation matrix determinant is not unity"));
   const Tensor<2, 3> R =
-    Transformations::Rotations::rotation_matrix_3d(Point<3>({0, 0, 1}), angle);
+    Transformations::Rotations::rotation_matrix_3d(Tensor<1, 3>({0., 0., 1.}),
+                                                   angle);
   Assert(std::abs(determinant(R) - 1.0) < 1e-9,
          ExcMessage("Rotation matrix determinant is not unity"));
 
@@ -53,7 +54,7 @@ test_rotation_matrix_3d_z_axis(const double angle)
 }
 
 void
-test_rotation_matrix_3d(const Point<3> &axis, const double angle)
+test_rotation_matrix_3d(const Tensor<1, 3> &axis, const double angle)
 {
   // http://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula
   // http://en.wikipedia.org/wiki/Rotation_matrix
@@ -84,8 +85,8 @@ test_rotation_matrix_3d(const Point<3> &axis, const double angle)
          ExcMessage("Incorrect computation of R in 3d"));
 }
 
-Point<3>
-normalise(const Point<3> &p)
+Tensor<1, 3>
+normalise(const Tensor<1, 3> &p)
 {
   Assert(p.norm() > 0.0, ExcMessage("Point vector has zero norm"));
   return p / p.norm();
@@ -146,9 +147,12 @@ main()
   test_rotation_matrix_3d_z_axis(45.0 * deg_to_rad);
   test_rotation_matrix_3d_z_axis(60.0 * deg_to_rad);
 
-  test_rotation_matrix_3d(normalise(Point<3>({1, 1, 1})), 90.0 * deg_to_rad);
-  test_rotation_matrix_3d(normalise(Point<3>({0, 2, 1})), 45.0 * deg_to_rad);
-  test_rotation_matrix_3d(normalise(Point<3>({-1, 3, 2})), 60.0 * deg_to_rad);
+  test_rotation_matrix_3d(normalise(Tensor<1, 3>({1., 1., 1.})),
+                          90.0 * deg_to_rad);
+  test_rotation_matrix_3d(normalise(Tensor<1, 3>({0., 2., 1.})),
+                          45.0 * deg_to_rad);
+  test_rotation_matrix_3d(normalise(Tensor<1, 3>({-1., 3., 2.})),
+                          60.0 * deg_to_rad);
 
   deallog << "OK" << std::endl;
 }

diff --git a/tests/simplex/step-18.cc b/tests/simplex/step-18.cc
index 269024ec61d6d9444a7e904c9c72fa15313f17c0..c715661a7f81e67c08ebc246d17f67fc729ca16f 100644 (file)
--- a/tests/simplex/step-18.cc
+++ b/tests/simplex/step-18.cc
@@ -163,9 +163,9 @@ namespace Step18
   Tensor<2, 3>
   get_rotation_matrix(const std::vector<Tensor<1, 3>> &grad_u)
   {
-    const Point<3> curl(grad_u[2][1] - grad_u[1][2],
-                        grad_u[0][2] - grad_u[2][0],
-                        grad_u[1][0] - grad_u[0][1]);
+    const Tensor<1, 3> curl({grad_u[2][1] - grad_u[1][2],
+                             grad_u[0][2] - grad_u[2][0],
+                             grad_u[1][0] - grad_u[0][1]});
 
     const double tan_angle = std::sqrt(curl * curl);
     const double angle     = std::atan(tan_angle);
@@ -177,7 +177,7 @@ namespace Step18
         return rot;
       }
 
-    const Point<3> axis = curl / tan_angle;
+    const Tensor<1, 3> axis = curl / tan_angle;
     return Physics::Transformations::Rotations::rotation_matrix_3d(axis,
                                                                    -angle);
   }