]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Fix SphericalManifold::get_intermediate_point for particular point combination.
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Mon, 20 Nov 2017 11:59:45 +0000 (12:59 +0100)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Mon, 20 Nov 2017 12:21:59 +0000 (13:21 +0100)
source/grid/manifold_lib.cc

index e4d0b77e7c5cd4a167f5dac4d1b85a643b61af1b..7ea56756c2a1a0dbcb06a3688b1189470bc2b921 100644 (file)
@@ -260,9 +260,9 @@ get_intermediate_point (const Point<spacedim> &p1,
 {
   const double tol = 1e-10;
 
-  if ( p1.distance(p2) < tol || w < tol)
+  if ( (p1-p2).norm_square() < tol*tol || std::abs(w) < tol)
     return p1;
-  else if (w > 1.0 - tol)
+  else if (std::abs(w-1.0) < tol)
     return p2;
 
   // If the points are one dimensional then there is no need for anything but
@@ -281,18 +281,21 @@ get_intermediate_point (const Point<spacedim> &p1,
   const Tensor<1,spacedim> e1 = v1/r1;
   const Tensor<1,spacedim> e2 = v2/r2;
 
-  // Treat points that are collinear with the center special.
-  if ((e1 + e2).norm_square() == 0.)
+  // Find the cosine of the angle gamma described by v1 and v2.
+  const double cosgamma = e1*e2;
+
+  // Points are collinear with the center (allow for 8*eps as a tolerance)
+  if (cosgamma < -1 + 8.*std::numeric_limits<double>::epsilon())
     return center;
 
-  if ((e1 - e2).norm_square() < tol*tol)
+  // Points are along a line, in which case e1 and e2 are essentially the same.
+  if (cosgamma > 1 - 8.*std::numeric_limits<double>::epsilon())
     return Point<spacedim>(center + w*v2 + (1-w)*v1);
 
-  // Find the angle gamma described by v1 and v2:
-  const double gamma = std::acos(e1*e2);
-
-  // Find the angle sigma that corresponds to arclength equal to w
-  const double sigma = w * gamma;
+  // Find the angle sigma that corresponds to arclength equal to w. acos
+  // should never be undefined because we have ruled out the two special cases
+  // above.
+  const double sigma = w * std::acos(cosgamma);
 
   // Normal to v1 in the plane described by v1,v2,and the origin.
   // Since p1 and p2 do not coincide n is not zero and well defined.

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.