From d3973743b91c132f1f5962e72ad84aa0aaada02e Mon Sep 17 00:00:00 2001
From: Wolfgang Bangerth <bangerth@math.tamu.edu>
Date: Sun, 28 Feb 2016 12:59:35 -0600
Subject: [PATCH] Provide a default implementation of
 Manifold::get_tangent_vector().

---
 include/deal.II/grid/manifold.h | 15 ++++++++-------
 source/grid/manifold.cc         | 18 ++++++++++++++----
 2 files changed, 22 insertions(+), 11 deletions(-)

diff --git a/include/deal.II/grid/manifold.h b/include/deal.II/grid/manifold.h
index c8390d49ac..e664042ed5 100644
--- a/include/deal.II/grid/manifold.h
+++ b/include/deal.II/grid/manifold.h
@@ -199,9 +199,11 @@ namespace Manifolds
  * \frac{\mathbf s(w)-\mathbf s(0)}{w}$ where all we need to do
  * is compute the new point $\mathbf s(w)$ with weights $w$ and
  * $1-w$ along the geodesic connecting $\mathbf x_1$ and $\mathbf x_2$.
+ * The default implementation of the function does this, by evaluating
+ * the quotient for a small but finite weight $w$.
  * In practice, however, it is almost always possible to explicitly
  * compute the direction vector, i.e., without the need to numerically
- * approximate the limit process.
+ * approximate the limit process, and derived classes should do so.
  *
  *
  * @ingroup manifold
@@ -365,12 +367,11 @@ public:
    * $\mathbf s(t)$ must move "faster" if the two points it connects between
    * arguments $t=0$ and $t=1$ are farther apart.
    *
-   * This function is used, among other cases, in computing normal vectors to
-   * faces or, more generally, surfaces such as the boundary. Since not all
-   * programs need this functionality, this function has a default
-   * implementation that just throws an exception. Consequently, derived
-   * classes only have to implement this function if the program that uses
-   * them does in fact call it directly or indirectly.
+   * The default implementation of this function approximates
+   * $\mathbf s'(0) \approx \frac{$\mathbf s(\epsilon)-\mathbf x_1}{\epsilon}$
+   * for a small value of $\epsilon$, and the evaluation of $\mathbf s(\epsilon)$
+   * is done by calling get_new_point(). If possible, derived classes should
+   * override this function by an implement of the exact derivative.
    *
    * @param x1 The first point that describes the geodesic, and the one
    *   at which the "direction" is to be evaluated.
diff --git a/source/grid/manifold.cc b/source/grid/manifold.cc
index 0bf9d3dcbc..6dbeaec9fc 100644
--- a/source/grid/manifold.cc
+++ b/source/grid/manifold.cc
@@ -214,11 +214,21 @@ get_new_point_on_hex (const Triangulation<3, 3>::hex_iterator &hex) const
 
 template <int dim, int spacedim>
 Tensor<1,spacedim>
-Manifold<dim,spacedim>::get_tangent_vector(const Point<spacedim> &,
-                                           const Point<spacedim> &) const
+Manifold<dim,spacedim>::get_tangent_vector(const Point<spacedim> &x1,
+                                           const Point<spacedim> &x2) const
 {
-  Assert (false, ExcPureFunctionCalled());
-  return Tensor<1,spacedim>();
+  const double epsilon = 1e-8;
+
+  std::vector<Point<spacedim> > q;
+  q.push_back(x1);
+  q.push_back(x2);
+
+  std::vector<double> w;
+  w.push_back(epsilon);
+  w.push_back(1.0-epsilon);
+
+  const Tensor<1,spacedim> neighbor_point = get_new_point (Quadrature<spacedim>(q, w));
+  return (neighbor_point-x1)/epsilon;
 }
 
 /* -------------------------- FlatManifold --------------------- */
-- 
2.39.5