From: Marc Fehling Date: Tue, 12 Oct 2021 20:10:40 +0000 (-0600) Subject: Update and unify interface for `Physics` and `GridTools` rotations. X-Git-Tag: v9.4.0-rc1~920^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=8542fabe3aaa857f8b156a6d8674b2762ce78879;p=dealii.git Update and unify interface for `Physics` and `GridTools` rotations. --- diff --git a/doc/news/changes/incompatibilities/20211011Fehling b/doc/news/changes/incompatibilities/20211011Fehling index 7105a29a87..6f18ea6fe4 100644 --- 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.
(Marc Fehling, 2021/10/11) diff --git a/examples/step-18/step-18.cc b/examples/step-18/step-18.cc index f0906082ac..de93a76b2c 100644 --- 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 bc2d57376c..80bb1252f4 100644 --- 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 d43eb0435b..e2d62c7977 100644 --- a/include/deal.II/grid/grid_tools.h +++ b/include/deal.II/grid/grid_tools.h @@ -569,9 +569,9 @@ namespace GridTools */ template void - rotate(const double angle, - const Point<3, double> &axis, - Triangulation & triangulation); + rotate(const Tensor<1, 3, double> &axis, + const double angle, + Triangulation & 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 e8bcb8a695..268a162f97 100644 --- a/include/deal.II/physics/transformations.h +++ b/include/deal.II/physics/transformations.h @@ -18,7 +18,6 @@ #include -#include #include #include @@ -90,6 +89,15 @@ namespace Physics */ template 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 + 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 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 +Tensor<2, 3, Number> +Physics::Transformations::Rotations::rotation_matrix_3d( + const Point<3, Number> &axis, + const Number & angle) +{ + return rotation_matrix_3d(static_cast>(axis), angle); +} + + + template 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 ff9868ae44..f7d6f104a2 100644 --- 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 9abd1b7733..ba6e7f2a0b 100644 --- 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 void - rotate(const double angle, - const Point<3, double> &axis, - Triangulation & triangulation) + rotate(const Tensor<1, 3, double> &axis, + const double angle, + Triangulation & 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 193bfe7dfc..170164eae3 100644 --- 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( + const Tensor<1, deal_II_space_dimension, double> &, const double, - const Point &, Triangulation &); template void diff --git a/tests/grid/rotate_01.cc b/tests/grid/rotate_01.cc index a49faef88c..9f2e08c0b0 100644 --- a/tests/grid/rotate_01.cc +++ b/tests/grid/rotate_01.cc @@ -65,11 +65,11 @@ test() GridGenerator::subdivided_parallelepiped(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 a6e79f794d..0e14c7ca95 100644 --- 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> &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 diff --git a/tests/physics/transformations-rotations_01.cc b/tests/physics/transformations-rotations_01.cc index 606422e1f3..84554d4541 100644 --- 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 269024ec61..c715661a7f 100644 --- 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> &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); }