From: Rene Gassmoeller <rene.gassmoeller@mailbox.org>
Date: Wed, 5 Apr 2017 00:11:05 +0000 (-0600)
Subject: Handle special case SphericalManifold::get_tangent_vector
X-Git-Tag: v9.0.0-rc1~1744^2~1
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=6ec0fe7336a6dcc0ff3bfc89391a7df00af75f36;p=dealii.git

Handle special case SphericalManifold::get_tangent_vector
---

diff --git a/source/grid/manifold_lib.cc b/source/grid/manifold_lib.cc
index 02a9147bf3..3579acfd3b 100644
--- a/source/grid/manifold_lib.cc
+++ b/source/grid/manifold_lib.cc
@@ -239,22 +239,28 @@ get_tangent_vector (const Point<spacedim> &p1,
   const double r1 = (p1 - center).norm();
   const double r2 = (p2 - center).norm();
 
-  Assert(r1 > 1e-10,
+  const double tolerance = 1e-10;
+
+  Assert(r1 > tolerance,
          ExcMessage("p1 cannot coincide with the center."));
 
-  Assert(r2 > 1e-10,
+  Assert(r2 > tolerance,
          ExcMessage("p2 cannot coincide with the center."));
 
   const Tensor<1,spacedim> e1 = (p1 - center)/r1;
   const Tensor<1,spacedim> e2 = (p2 - center)/r2;
 
-  Assert(e1*e2 + 1.0 > 1e-10,
+  Assert(e1*e2 + 1.0 > tolerance,
          ExcMessage("p1 and p2 cannot lie on the same diameter and be opposite "
                     "respect to the center."));
 
   // Tangent vector to the unit sphere along the geodesic given by e1 and e2.
   Tensor<1,spacedim> tg = (e2-(e2*e1)*e1);
-  tg = tg / tg.norm();
+
+  // There is a special case if e2*e1, in which case tg=0
+  const double tg_norm = tg.norm();
+  if (tg_norm > tolerance)
+    tg /= tg_norm;
 
   const double gamma = std::acos(e1*e2);