Implement CylindricalManifold::push_forward_gradient

diff --git a/include/deal.II/grid/manifold_lib.h b/include/deal.II/grid/manifold_lib.h
index feac819df874639e5e0f63c0fa302bf4731ca1ac..27f44282f15b226375d23a98a336137bbbbd2984 100644 (file)
--- 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 72dabc76567c1768832bf334ce71fe45cd9cb6f7..52e836103812ad3558e8bd3db987e5d17b446f7a 100644 (file)
--- 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
 // ============================================================