From 3a0dc2fd4a10d2386baf2a868287961b0525fbc4 Mon Sep 17 00:00:00 2001
From: Daniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Date: Thu, 4 Jan 2018 17:08:21 +0100
Subject: [PATCH] Implement CylindricalManifold::push_forward_gradient

---
 include/deal.II/grid/manifold_lib.h |  9 ++++++++
 source/grid/manifold_lib.cc         | 36 +++++++++++++++++++++++++++++
 2 files changed, 45 insertions(+)

diff --git a/include/deal.II/grid/manifold_lib.h b/include/deal.II/grid/manifold_lib.h
index feac819df8..27f44282f1 100644
--- a/include/deal.II/grid/manifold_lib.h
+++ b/include/deal.II/grid/manifold_lib.h
@@ -381,6 +381,15 @@ public:
   virtual Point<spacedim>
   push_forward(const Point<3> &chart_point) const override;
 
+  /**
+   * Compute the derivatives of the mapping from cylindrical coordinates
+   * $(r, \phi, \lambda)$ to cartesian coordinates where $r$ denotes the
+   * distance from the axis, $\phi$ the angle between the given point and the
+   * computed normal direction and $\lambda$ the axial position.
+   */
+  virtual DerivativeForm<1, 3, spacedim>
+  push_forward_gradient(const Point<3> &chart_point) const override;
+
   /**
    * Compute new points on the CylindricalManifold. See the documentation of
    * the base class for a detailed description of what this function does.
diff --git a/source/grid/manifold_lib.cc b/source/grid/manifold_lib.cc
index 72dabc7656..52e8361038 100644
--- a/source/grid/manifold_lib.cc
+++ b/source/grid/manifold_lib.cc
@@ -906,6 +906,42 @@ CylindricalManifold<dim, spacedim>::push_forward(const Point<3> &chart_point) co
 
 
 
+template<int dim, int spacedim>
+DerivativeForm<1, 3, spacedim>
+CylindricalManifold<dim, spacedim>::push_forward_gradient(const Point<3> &chart_point) const
+{
+  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.
+  const double sine = std::sin(chart_point(1));
+  const double cosine = std::cos(chart_point(1));
+  const Tensor<1,3> dxn = cross_product_3d(direction, normal_direction);
+  const Tensor<1,3> intermediate = normal_direction*cosine+dxn*sine;
+
+  // derivative w.r.t the radius
+  derivatives[0][0] = intermediate[0];
+  derivatives[1][0] = intermediate[1];
+  derivatives[2][0] = intermediate[2];
+
+  // derivatives w.r.t the angle
+  derivatives[0][1] = -normal_direction[0]*sine + dxn[0]*cosine;
+  derivatives[1][1] = -normal_direction[1]*sine + dxn[1]*cosine;
+  derivatives[2][1] = -normal_direction[2]*sine + dxn[2]*cosine;
+
+  // derivatives w.r.t the direction of the axis
+  derivatives[0][2] = direction[0];
+  derivatives[1][2] = direction[1];
+  derivatives[2][2] = direction[2];
+
+  return derivatives;
+}
+
+
+
 // ============================================================
 // FunctionManifold
 // ============================================================
-- 
2.39.5