Added tangent vector to SphericalManifold.

diff --git a/include/deal.II/grid/manifold_lib.h b/include/deal.II/grid/manifold_lib.h
index 78694794e8a9a60e97d310f428431172043850dc..21de7b05552a1eefca372bd3609ec65b09f8826a 100644 (file)
--- a/include/deal.II/grid/manifold_lib.h
+++ b/include/deal.II/grid/manifold_lib.h
@@ -84,6 +84,24 @@ public:
   virtual Point<spacedim>
   push_forward(const Point<spacedim> &chart_point) const;
 
+
+  /**
+   * Given a point in the spacedim dimensional Euclidean space, this
+   * method returns the derivatives of the function $F$ that maps from
+   * the polar coordinate system to the Euclidean coordinate
+   * system. In other words, it is a matrix of size
+   * $\text{spacedim}\times\text{spacedim}$.
+   *
+   * This function is used in the computations required by the
+   * get_tangent_vector() function.
+   *
+   * Refer to the general documentation of this class for more information.
+   */
+  virtual
+  DerivativeForm<1,spacedim,spacedim>
+  push_forward_gradient(const Point<spacedim> &chart_point) const;
+
+
   /**
    * Let the new point be the average sum of surrounding vertices.
    *

diff --git a/source/grid/manifold_lib.cc b/source/grid/manifold_lib.cc
index cfe05206f9c3877bad85ce1f5e77b0ae304ba3aa..bb5b91c8610d4330c11b5360990248c0db32b251 100644 (file)
--- a/source/grid/manifold_lib.cc
+++ b/source/grid/manifold_lib.cc
@@ -131,6 +131,47 @@ SphericalManifold<dim,spacedim>::pull_back(const Point<spacedim> &space_point) c
 }
 
 
+template <int dim, int spacedim>
+DerivativeForm<1,spacedim,spacedim>
+SphericalManifold<dim,spacedim>::push_forward_gradient(const Point<spacedim> &spherical_point) const
+{
+  Assert(spherical_point[0] >=0.0,
+         ExcMessage("Negative radius for given point."));
+  const double rho = spherical_point[0];
+  const double theta = spherical_point[1];
+
+  DerivativeForm<1,spacedim,spacedim> DX;
+  if (rho > 1e-10)
+    switch (spacedim)
+      {
+      case 2:
+        DX[0][0] = cos(theta);
+        DX[0][1] = -rho*sin(theta);
+        DX[1][0] = sin(theta);
+        DX[1][1] = rho*cos(theta);
+        break;
+      case 3:
+      {
+        const double &phi= spherical_point[2];
+        DX[0][0] =      sin(theta)*cos(phi);
+        DX[0][1] =  rho*cos(theta)*cos(phi);
+        DX[0][2] = -rho*sin(theta)*sin(phi);
+
+        DX[1][0] =     sin(theta)*sin(phi);
+        DX[1][1] = rho*cos(theta)*sin(phi);
+        DX[1][2] = rho*sin(theta)*cos(phi);
+
+        DX[2][0] =      cos(theta);
+        DX[2][1] = -rho*sin(theta);
+        DX[2][2] = 0;
+      }
+      break;
+      default:
+        Assert(false, ExcInternalError());
+      }
+  return DX;
+}
+
 // ============================================================
 // CylindricalManifold
 // ============================================================

diff --git a/tests/manifold/spherical_manifold_06.cc b/tests/manifold/spherical_manifold_06.cc
new file mode 100644 (file)
index 0000000..c2800ae
--- /dev/null
+++ b/tests/manifold/spherical_manifold_06.cc
@@ -0,0 +1,68 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// Check get_tangent_vector for spherical manifold, on simple points.
+
+#include "../tests.h"
+
+#include <deal.II/base/utilities.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/fe/mapping_q_generic.h>
+#include <deal.II/fe/mapping_manifold.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/grid/manifold_lib.h>
+
+#include <numeric>
+
+template<int dim, int spacedim>
+void test()
+{
+  deallog << "dim=" << dim << ", spacedim=" << spacedim << std::endl;
+
+  Point<spacedim> center;
+  static  const SphericalManifold<dim,spacedim> manifold(center);
+
+  // Go from 0,1 to 1,0
+  Point<spacedim> p0,p1;
+  p0[1] = 1.0;
+  p1[0] = 1.0;
+
+  Tensor<1,spacedim> T = manifold.get_tangent_vector(p0, p1);
+
+  deallog << "P0      : " << p0 << std::endl;
+  deallog << "P1      : " << p1 << std::endl;
+  deallog << "T(P0-P1): " << T << std::endl;
+  deallog << "Error   : " << T.norm() - numbers::PI/2 << std::endl;
+}
+
+
+int
+main()
+{
+  std::ofstream logfile ("output");
+  deallog.attach(logfile);
+  deallog.threshold_double(1.e-10);
+
+  test<2,2>();
+  test<2,3>();
+
+  test<3,3>();
+
+  return 0;
+}
+
+
+

diff --git a/tests/manifold/spherical_manifold_06.output b/tests/manifold/spherical_manifold_06.output
new file mode 100644 (file)
index 0000000..d906a07
--- /dev/null
+++ b/tests/manifold/spherical_manifold_06.output
@@ -0,0 +1,16 @@
+
+DEAL::dim=2, spacedim=2
+DEAL::P0      : 0.00000 1.00000
+DEAL::P1      : 1.00000 0.00000
+DEAL::T(P0-P1): 1.57080 -9.61835e-17
+DEAL::Error   : 0
+DEAL::dim=2, spacedim=3
+DEAL::P0      : 0.00000 1.00000 0.00000
+DEAL::P1      : 1.00000 0.00000 0.00000
+DEAL::T(P0-P1): 1.57080 -9.61835e-17 0.00000
+DEAL::Error   : 0
+DEAL::dim=3, spacedim=3
+DEAL::P0      : 0.00000 1.00000 0.00000
+DEAL::P1      : 1.00000 0.00000 0.00000
+DEAL::T(P0-P1): 1.57080 -9.61835e-17 0.00000
+DEAL::Error   : 0