Skip to content

Commit

Permalink
Update and unify interface.
Browse files Browse the repository at this point in the history
  • Loading branch information
marcfehling committed Oct 12, 2021
1 parent 0cea6f7 commit 8e23d53
Show file tree
Hide file tree
Showing 10 changed files with 40 additions and 37 deletions.
6 changes: 3 additions & 3 deletions include/deal.II/base/symmetric_tensor.templates.h
Original file line number Diff line number Diff line change
Expand Up @@ -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());
Expand Down
6 changes: 3 additions & 3 deletions include/deal.II/grid/grid_tools.h
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
7 changes: 3 additions & 4 deletions include/deal.II/physics/transformations.h
Original file line number Diff line number Diff line change
Expand Up @@ -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>

Expand Down Expand Up @@ -90,7 +89,7 @@ namespace Physics
*/
template <typename Number>
Tensor<2, 3, Number>
rotation_matrix_3d(const Point<3, Number> &axis, const Number &angle);
rotation_matrix_3d(const Tensor<1, 3, Number> &axis, const Number &angle);

//@}

Expand Down Expand Up @@ -914,8 +913,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."));
Expand Down
2 changes: 1 addition & 1 deletion source/grid/grid_generator.cc
Original file line number Diff line number Diff line change
Expand Up @@ -4976,7 +4976,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
Expand Down
14 changes: 7 additions & 7 deletions source/grid/grid_tools.cc
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand Down Expand Up @@ -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);
}


Expand All @@ -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);
}


Expand Down
2 changes: 1 addition & 1 deletion source/grid/grid_tools.inst.in
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
6 changes: 3 additions & 3 deletions tests/grid/rotate_01.cc
Original file line number Diff line number Diff line change
Expand Up @@ -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());
}

Expand Down
12 changes: 6 additions & 6 deletions tests/physics/step-18-rotation_matrix.cc
Original file line number Diff line number Diff line change
Expand Up @@ -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>
Expand Down
14 changes: 9 additions & 5 deletions tests/physics/transformations-rotations_01.cc
Original file line number Diff line number Diff line change
Expand Up @@ -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"));

Expand All @@ -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
Expand Down Expand Up @@ -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;
}
8 changes: 4 additions & 4 deletions tests/simplex/step-18.cc
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand All @@ -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);
}
Expand Down

0 comments on commit 8e23d53

Please sign in to comment.