From: Marc Fehling Date: Tue, 12 Oct 2021 02:06:19 +0000 (-0600) Subject: Generalize `GridTools::rotate` for arbitrary rotation axes. X-Git-Tag: v9.4.0-rc1~920^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=4ebbba75f8fb9891f7ff4998d0ac51c133469e1c;p=dealii.git Generalize `GridTools::rotate` for arbitrary rotation axes. --- diff --git a/doc/news/changes/incompatibilities/20211011Fehling b/doc/news/changes/incompatibilities/20211011Fehling new file mode 100644 index 0000000000..7105a29a87 --- /dev/null +++ b/doc/news/changes/incompatibilities/20211011Fehling @@ -0,0 +1,5 @@ +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. +
+(Marc Fehling, 2021/10/11) diff --git a/doc/news/changes/minor/20211011Fehling b/doc/news/changes/minor/20211011Fehling new file mode 100644 index 0000000000..020002bebf --- /dev/null +++ b/doc/news/changes/minor/20211011Fehling @@ -0,0 +1,4 @@ +Changed: The 3D implementation of GridTools::rotate() now accepts unit +vectors as rotation axes. +
+(Marc Fehling, 2021/10/11) diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h index b639f59453..d43eb0435b 100644 --- a/include/deal.II/grid/grid_tools.h +++ b/include/deal.II/grid/grid_tools.h @@ -14,64 +14,65 @@ // --------------------------------------------------------------------- #ifndef dealii_grid_tools_h -# define dealii_grid_tools_h +#define dealii_grid_tools_h -# include +#include -# include -# include -# include +#include +#include +#include +#include -# include +#include -# include +#include -# include +#include -# include -# include -# include +#include +#include +#include -# include -# include -# include -# include +#include +#include +#include +#include -# include +#include -# include -# include -# include -# include -# include +#include +#include +#include +#include +#include -# include +#include DEAL_II_DISABLE_EXTRA_DIAGNOSTICS -# include -# include -# include -# include -# include -# include - -# ifdef DEAL_II_WITH_ZLIB -# include -# include -# include -# include -# endif +#include +#include +#include +#include +#include +#include + +#ifdef DEAL_II_WITH_ZLIB +# include +# include +# include +# include +#endif DEAL_II_ENABLE_EXTRA_DIAGNOSTICS -# include -# include -# include +#include +#include +#include DEAL_II_NAMESPACE_OPEN // Forward declarations -# ifndef DOXYGEN +#ifndef DOXYGEN namespace parallel { namespace distributed @@ -88,7 +89,7 @@ namespace hp } class SparsityPattern; -# endif +#endif namespace internal { @@ -96,14 +97,14 @@ namespace internal class ActiveCellIterator { public: -# ifndef _MSC_VER +#ifndef _MSC_VER using type = typename MeshType::active_cell_iterator; -# else +#else using type = TriaActiveIterator>; -# endif +#endif }; -# ifdef _MSC_VER +#ifdef _MSC_VER template class ActiveCellIterator> { @@ -111,7 +112,7 @@ namespace internal using type = TriaActiveIterator>; }; -# endif +#endif } // namespace internal /** @@ -555,6 +556,23 @@ namespace GridTools void rotate(const double angle, Triangulation &triangulation); + /** + * Rotate all vertices of the given @p triangulation in counter-clockwise + * direction around the @p axis described by a unit vector. Otherwise like the + * function above. + * + * @param[in] angle Angle in radians to rotate the Triangulation by. + * @param[in] axis A unit vector that defines the axis of rotation. + * @param[in,out] triangulation The Triangulation object to rotate. + * + * @note Implemented for spacedim=3 and dim=1, 2, and 3. + */ + template + void + rotate(const double angle, + const Point<3, double> &axis, + Triangulation & triangulation); + /** * Rotate all vertices of the given @p triangulation in counter-clockwise * direction around the axis with the given index. Otherwise like the @@ -566,9 +584,11 @@ namespace GridTools * @param[in,out] triangulation The Triangulation object to rotate. * * @note Implemented for dim=1, 2, and 3. + * + * @deprecated Use the alternative with the unit vector instead. */ template - void + DEAL_II_DEPRECATED_EARLY void rotate(const double angle, const unsigned int axis, Triangulation &triangulation); @@ -904,14 +924,14 @@ namespace GridTools * the performance if the function is called only once on few points. */ template -# ifndef DOXYGEN +#ifndef DOXYGEN std::tuple< std::vector::active_cell_iterator>, std::vector>>, std::vector>> -# else +#else return_type -# endif +#endif compute_point_locations( const Cache & cache, const std::vector> &points, @@ -953,15 +973,15 @@ namespace GridTools * GridTools::compute_point_locations(). */ template -# ifndef DOXYGEN +#ifndef DOXYGEN std::tuple< std::vector::active_cell_iterator>, std::vector>>, std::vector>, std::vector> -# else +#else return_type -# endif +#endif compute_point_locations_try_all( const Cache & cache, const std::vector> &points, @@ -1039,16 +1059,16 @@ namespace GridTools * of this page. */ template -# ifndef DOXYGEN +#ifndef DOXYGEN std::tuple< std::vector::active_cell_iterator>, std::vector>>, std::vector>, std::vector>>, std::vector>> -# else +#else return_type -# endif +#endif distributed_compute_point_locations( const GridTools::Cache & cache, const std::vector> & local_points, @@ -1285,13 +1305,13 @@ namespace GridTools * for this case. */ template class MeshType, int spacedim> -# ifndef _MSC_VER +#ifndef _MSC_VER std::vector::active_cell_iterator> -# else +#else std::vector< typename dealii::internal:: ActiveCellIterator>::type> -# endif +#endif find_cells_adjacent_to_vertex(const MeshType &container, const unsigned int vertex_index); @@ -1358,13 +1378,13 @@ namespace GridTools * processor will return a locally owned cell and the other one a ghost cell. */ template class MeshType, int spacedim> -# ifndef _MSC_VER +#ifndef _MSC_VER std::pair::active_cell_iterator, Point> -# else +#else std::pair>::type, Point> -# endif +#endif find_active_cell_around_point(const Mapping & mapping, const MeshType &mesh, const Point & p, @@ -1379,12 +1399,12 @@ namespace GridTools * @return An iterator into the mesh that points to the surrounding cell. */ template class MeshType, int spacedim> -# ifndef _MSC_VER +#ifndef _MSC_VER typename MeshType::active_cell_iterator -# else +#else typename dealii::internal:: ActiveCellIterator>::type -# endif +#endif find_active_cell_around_point(const MeshType &mesh, const Point & p, const std::vector &marked_vertices = {}, @@ -1481,13 +1501,13 @@ namespace GridTools * call the function above with argument `cache` in this case. */ template class MeshType, int spacedim> -# ifndef _MSC_VER +#ifndef _MSC_VER std::pair::active_cell_iterator, Point> -# else +#else std::pair>::type, Point> -# endif +#endif find_active_cell_around_point( const Mapping & mapping, const MeshType &mesh, @@ -1529,15 +1549,15 @@ namespace GridTools * @endcode */ template class MeshType, int spacedim> -# ifndef _MSC_VER +#ifndef _MSC_VER std::vector::active_cell_iterator, Point>> -# else +#else std::vector>::type, Point>> -# endif +#endif find_all_active_cells_around_point( const Mapping & mapping, const MeshType &mesh, @@ -1553,15 +1573,15 @@ namespace GridTools * function find_all_active_cells_around_point() above. */ template class MeshType, int spacedim> -# ifndef _MSC_VER +#ifndef _MSC_VER std::vector::active_cell_iterator, Point>> -# else +#else std::vector>::type, Point>> -# endif +#endif find_all_active_cells_around_point( const Mapping & mapping, const MeshType &mesh, @@ -1908,13 +1928,13 @@ namespace GridTools * of this page. */ template -# ifndef DOXYGEN +#ifndef DOXYGEN std::tuple>, std::map, std::map>> -# else +#else return_type -# endif +#endif guess_point_owner( const std::vector>> &global_bboxes, const std::vector> & points); @@ -1955,13 +1975,13 @@ namespace GridTools * of this page. */ template -# ifndef DOXYGEN +#ifndef DOXYGEN std::tuple>, std::map, std::map>> -# else +#else return_type -# endif +#endif guess_point_owner( const RTree, unsigned int>> &covering_rtree, const std::vector> & points); @@ -3254,7 +3274,7 @@ namespace GridTools void load(Archive &ar, const unsigned int version); -# ifdef DOXYGEN +#ifdef DOXYGEN /** * Read or write the data of this object to or from a stream for the * purpose of serialization using the [BOOST serialization @@ -3263,11 +3283,11 @@ namespace GridTools template void serialize(Archive &archive, const unsigned int version); -# else +#else // This macro defines the serialize() method that is compatible with // the templated save() and load() method that have been implemented. BOOST_SERIALIZATION_SPLIT_MEMBER() -# endif +#endif }; @@ -3455,7 +3475,7 @@ DeclExceptionMsg(ExcMeshNotOrientable, /* ----------------- Template function --------------- */ -# ifndef DOXYGEN +#ifndef DOXYGEN namespace GridTools { @@ -4260,7 +4280,7 @@ namespace GridTools MeshType::space_dimension> &)> &compute_ghost_owners) { -# ifndef DEAL_II_WITH_MPI +# ifndef DEAL_II_WITH_MPI (void)mesh; (void)pack; (void)unpack; @@ -4268,7 +4288,7 @@ namespace GridTools (void)process_cells; (void)compute_ghost_owners; Assert(false, ExcNeedsMPI()); -# else +# else constexpr int dim = MeshType::dimension; constexpr int spacedim = MeshType::space_dimension; auto tria = @@ -4490,7 +4510,7 @@ namespace GridTools } -# endif // DEAL_II_WITH_MPI +# endif // DEAL_II_WITH_MPI } } // namespace internal @@ -4506,13 +4526,13 @@ namespace GridTools const std::function &cell_filter) { -# ifndef DEAL_II_WITH_MPI +# ifndef DEAL_II_WITH_MPI (void)mesh; (void)pack; (void)unpack; (void)cell_filter; Assert(false, ExcNeedsMPI()); -# else +# else internal::exchange_cell_data( @@ -4526,7 +4546,7 @@ namespace GridTools process(cell, cell->subdomain_id()); }, [](const auto &tria) { return tria.ghost_owners(); }); -# endif +# endif } @@ -4542,13 +4562,13 @@ namespace GridTools const std::function &cell_filter) { -# ifndef DEAL_II_WITH_MPI +# ifndef DEAL_II_WITH_MPI (void)mesh; (void)pack; (void)unpack; (void)cell_filter; Assert(false, ExcNeedsMPI()); -# else +# else internal::exchange_cell_data( @@ -4564,15 +4584,12 @@ namespace GridTools process(cell, cell->level_subdomain_id()); }, [](const auto &tria) { return tria.level_ghost_owners(); }); -# endif +# endif } } // namespace GridTools -# endif +#endif // DOXYGEN DEAL_II_NAMESPACE_CLOSE -/*---------------------------- grid_tools.h ---------------------------*/ -/* end of #ifndef dealii_grid_tools_h */ #endif -/*---------------------------- grid_tools.h ---------------------------*/ diff --git a/source/grid/grid_generator.cc b/source/grid/grid_generator.cc index 814300b3b1..ff9868ae44 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, 1, triangulation); + GridTools::rotate(numbers::PI_2, Point<3>({0., 1., 0.}), 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 4b37c3eb1f..9abd1b7733 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -56,6 +56,8 @@ #include #include +#include + DEAL_II_DISABLE_EXTRA_DIAGNOSTICS #include #include @@ -1970,33 +1972,23 @@ namespace GridTools class Rotate3d { public: - Rotate3d(const double angle, const unsigned int axis) - : angle(angle) - , axis(axis) + Rotate3d(const double angle, const Point<3, double> &axis) + : rotation_matrix( + Physics::Transformations::Rotations::rotation_matrix_3d(axis, + angle)) {} Point<3> operator()(const Point<3> &p) const { - if (axis == 0) - return {p(0), - std::cos(angle) * p(1) - std::sin(angle) * p(2), - std::sin(angle) * p(1) + std::cos(angle) * p(2)}; - else if (axis == 1) - return {std::cos(angle) * p(0) + std::sin(angle) * p(2), - p(1), - -std::sin(angle) * p(0) + std::cos(angle) * p(2)}; - else - return {std::cos(angle) * p(0) - std::sin(angle) * p(1), - std::sin(angle) * p(0) + std::cos(angle) * p(1), - p(2)}; + return static_cast>(rotation_matrix * p); } private: - const double angle; - const unsigned int axis; + const Tensor<2, 3, double> rotation_matrix; }; + template class Scale { @@ -2025,6 +2017,16 @@ namespace GridTools } + template + void + rotate(const double angle, + const Point<3, double> &axis, + Triangulation & triangulation) + { + transform(internal::Rotate3d(angle, axis), triangulation); + } + + template void rotate(const double angle, @@ -2033,9 +2035,13 @@ namespace GridTools { Assert(axis < 3, ExcMessage("Invalid axis given!")); - transform(internal::Rotate3d(angle, axis), triangulation); + Point<3, double> vector; + vector[axis] = 1.; + + transform(internal::Rotate3d(angle, vector), triangulation); } + template void scale(const double scaling_factor, diff --git a/source/grid/grid_tools.inst.in b/source/grid/grid_tools.inst.in index eb2215c3f8..193bfe7dfc 100644 --- a/source/grid/grid_tools.inst.in +++ b/source/grid/grid_tools.inst.in @@ -271,6 +271,12 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) Triangulation &); # if deal_II_space_dimension == 3 + template void + rotate( + const double, + const Point &, + Triangulation &); + template void rotate( const double, diff --git a/source/grid/grid_tools_nontemplates.cc b/source/grid/grid_tools_nontemplates.cc index 158d76d113..38f9d14cc2 100644 --- a/source/grid/grid_tools_nontemplates.cc +++ b/source/grid/grid_tools_nontemplates.cc @@ -17,6 +17,8 @@ #include +#include + #include // GridTools functions that are template specializations (i.e., only compiled @@ -399,24 +401,23 @@ namespace GridTools namespace { - // the following class is only - // needed in 2d, so avoid trouble - // with compilers warning otherwise + // the following class is only needed in 2d, so avoid trouble with compilers + // warning otherwise class Rotate2d { public: explicit Rotate2d(const double angle) - : angle(angle) + : rotation_matrix( + Physics::Transformations::Rotations::rotation_matrix_2d(angle)) {} Point<2> operator()(const Point<2> &p) const { - return {std::cos(angle) * p(0) - std::sin(angle) * p(1), - std::sin(angle) * p(0) + std::cos(angle) * p(1)}; + return static_cast>(rotation_matrix * p); } private: - const double angle; + const Tensor<2, 2, double> rotation_matrix; }; } // namespace diff --git a/tests/fe/bdm_2.cc b/tests/fe/bdm_2.cc index 46703c566d..352de8c28d 100644 --- a/tests/fe/bdm_2.cc +++ b/tests/fe/bdm_2.cc @@ -19,6 +19,7 @@ // and on each of these cells each time we evaluate the shape // functions. +#include #include #include @@ -69,7 +70,7 @@ transform_grid(Triangulation<2> &tria, const unsigned int transform) // second round: rotate // triangulation case 1: - GridTools::rotate(3.14159265358 / 2, tria); + GridTools::rotate(numbers::PI_2, tria); break; // third round: inflate @@ -83,7 +84,7 @@ transform_grid(Triangulation<2> &tria, const unsigned int transform) // stretch case 3: GridTools::scale(.5, tria); - GridTools::rotate(-3.14159265358 / 2, tria); + GridTools::rotate(-numbers::PI_2, tria); GridTools::transform(&stretch_coordinates, tria); break; diff --git a/tests/fe/nedelec.cc b/tests/fe/nedelec.cc index 7e074536cd..4396c24317 100644 --- a/tests/fe/nedelec.cc +++ b/tests/fe/nedelec.cc @@ -25,6 +25,7 @@ // perl -n -e 'print if s/DEAL:NedelecK-TransformN::value//' output +#include #include #include @@ -75,7 +76,7 @@ transform_grid(Triangulation<2> &tria, const unsigned int transform) // second round: rotate // triangulation case 1: - GridTools::rotate(3.14159265358 / 2, tria); + GridTools::rotate(numbers::PI_2, tria); break; // third round: inflate @@ -89,7 +90,7 @@ transform_grid(Triangulation<2> &tria, const unsigned int transform) // stretch case 3: GridTools::scale(.5, tria); - GridTools::rotate(-3.14159265358 / 2, tria); + GridTools::rotate(-numbers::PI_2, tria); GridTools::transform(&stretch_coordinates, tria); break; diff --git a/tests/fe/rt_2.cc b/tests/fe/rt_2.cc index 92df4009d4..1757083cd7 100644 --- a/tests/fe/rt_2.cc +++ b/tests/fe/rt_2.cc @@ -19,6 +19,7 @@ // and on each of these cells each time we evaluate the shape // functions. +#include #include #include @@ -69,7 +70,7 @@ transform_grid(Triangulation<2> &tria, const unsigned int transform) // second round: rotate // triangulation case 1: - GridTools::rotate(3.14159265358 / 2, tria); + GridTools::rotate(numbers::PI_2, tria); break; // third round: inflate @@ -83,7 +84,7 @@ transform_grid(Triangulation<2> &tria, const unsigned int transform) // stretch case 3: GridTools::scale(.5, tria); - GridTools::rotate(-3.14159265358 / 2, tria); + GridTools::rotate(-numbers::PI_2, tria); GridTools::transform(&stretch_coordinates, tria); break; diff --git a/tests/fe/rt_bubbles_2.cc b/tests/fe/rt_bubbles_2.cc index e6462d54a0..5d44fb2a02 100644 --- a/tests/fe/rt_bubbles_2.cc +++ b/tests/fe/rt_bubbles_2.cc @@ -19,6 +19,7 @@ // and on each of these transformed cells we evaluate the shape // functions. +#include #include #include @@ -67,7 +68,7 @@ transform_grid(Triangulation<2> &tria, const unsigned int transform) // second test: rotate cell case 1: - GridTools::rotate(3.14159265358 / 2, tria); + GridTools::rotate(numbers::PI_2, tria); break; // third test: scale cell by a factor of 2 @@ -79,7 +80,7 @@ transform_grid(Triangulation<2> &tria, const unsigned int transform) // stretch case 3: GridTools::scale(.5, tria); - GridTools::rotate(-3.14159265358 / 2, tria); + GridTools::rotate(-numbers::PI_2, tria); GridTools::transform(&stretch_coordinates, tria); break; diff --git a/tests/grid/grid_tools.cc b/tests/grid/grid_tools.cc index 0196fefcda..667fb32566 100644 --- a/tests/grid/grid_tools.cc +++ b/tests/grid/grid_tools.cc @@ -15,6 +15,8 @@ +#include + #include #include #include @@ -77,7 +79,7 @@ test2() GridOut().write_gnuplot(tria, logfile); logfile << "Rotated grid:" << std::endl; - GridTools::rotate(3.14159265258 / 4, tria); + GridTools::rotate(numbers::PI_4, tria); GridOut().write_gnuplot(tria, logfile); } diff --git a/tests/grid/rotate_01.cc b/tests/grid/rotate_01.cc index f3f75a3f7f..a49faef88c 100644 --- a/tests/grid/rotate_01.cc +++ b/tests/grid/rotate_01.cc @@ -15,6 +15,8 @@ // test GridTools::rotate(), output is checked for correctness visually +#include + #include #include #include @@ -63,11 +65,11 @@ test() GridGenerator::subdivided_parallelepiped(tria, origin, edges); // GridOut().write_gnuplot (tria, deallog.get_file_stream()); - GridTools::rotate(numbers::PI / 3.0, 0, tria); + GridTools::rotate(numbers::PI / 3.0, Point<3>({1., 0., 0.}), tria); // GridOut().write_gnuplot (tria, deallog.get_file_stream()); - GridTools::rotate(numbers::PI / 10.0, 1, tria); + GridTools::rotate(numbers::PI / 10.0, Point<3>({0., 1., 0.}), tria); // GridOut().write_gnuplot (tria, deallog.get_file_stream()); - GridTools::rotate(-numbers::PI / 5.0, 2, tria); + GridTools::rotate(-numbers::PI / 5.0, Point<3>({0., 0., 1.}), tria); GridOut().write_gnuplot(tria, deallog.get_file_stream()); }