From: Martin Kronbichler <kronbichler@lnm.mw.tum.de>
Date: Tue, 5 Mar 2019 17:40:11 +0000 (+0100)
Subject: Fix EllipticalManifold::push_forward_gradient
X-Git-Tag: v9.1.0-rc1~303^2~4
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=dc21bf441614471a2a330a28eb1aab455c3180eb;p=dealii.git

Fix EllipticalManifold::push_forward_gradient
---

diff --git a/source/grid/manifold_lib.cc b/source/grid/manifold_lib.cc
index 4dd1c9a31d..c08bff2f32 100644
--- a/source/grid/manifold_lib.cc
+++ b/source/grid/manifold_lib.cc
@@ -1091,11 +1091,7 @@ CylindricalManifold<dim, spacedim>::push_forward(
   Assert(spacedim == 3,
          ExcMessage("CylindricalManifold can only be used for spacedim==3!"));
 
-  // Rotate the orthogonal direction by the given angle.
-  // Formula from Section 5.2 in
-  // http://inside.mines.edu/fs_home/gmurray/ArbitraryAxisRotation/
-  // simplified assuming normal_direction and direction are orthogonal
-  // and unit vectors.
+  // Rotate the orthogonal direction by the given angle
   const double sine_r           = std::sin(chart_point(1)) * chart_point(0);
   const double cosine_r         = std::cos(chart_point(1)) * chart_point(0);
   const Tensor<1, spacedim> dxn = cross_product_3d(direction, normal_direction);
@@ -1118,11 +1114,7 @@ CylindricalManifold<dim, spacedim>::push_forward_gradient(
 
   Tensor<2, 3> derivatives;
 
-  // Rotate the orthogonal direction by the given angle.
-  // Formula from Section 5.2 in
-  // http://inside.mines.edu/fs_home/gmurray/ArbitraryAxisRotation/
-  // simplified assuming normal_direction and direction are orthogonal
-  // and unit vectors.
+  // Rotate the orthogonal direction by the given angle
   const double              sine   = std::sin(chart_point(1));
   const double              cosine = std::cos(chart_point(1));
   const Tensor<1, spacedim> dxn = cross_product_3d(direction, normal_direction);
@@ -1289,14 +1281,19 @@ DerivativeForm<1, 2, 2>
 EllipticalManifold<2, 2>::push_forward_gradient(
   const Point<2> &chart_point) const
 {
-  const double            cs = std::cos(chart_point[1]);
-  const double            sn = std::sin(chart_point[1]);
-  DerivativeForm<1, 2, 2> dX;
+  const double cs = std::cos(chart_point[1]);
+  const double sn = std::sin(chart_point[1]);
+  Tensor<2, 2> dX;
   dX[0][0] = cosh_u * cs;
   dX[0][1] = -chart_point[0] * cosh_u * sn;
   dX[1][0] = sinh_u * sn;
-  dX[1][1] = chart_point[1] * sinh_u * cs;
-  return dX;
+  dX[1][1] = chart_point[0] * sinh_u * cs;
+
+  // rotate according to the major axis direction
+  Tensor<2, 2, double> rot{
+    {{+direction[0], -direction[1]}, {direction[1], direction[0]}}};
+
+  return rot * dX;
 }