From: David Wells <wellsd2@rpi.edu>
Date: Thu, 18 Jan 2018 21:42:00 +0000 (-0500)
Subject: Make the TFI chart point calculation more robust.
X-Git-Tag: v9.0.0-rc1~544^2~2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b24f0fdd723a717446409c1aed6c69767da4936d;p=dealii.git

Make the TFI chart point calculation more robust.

If the guess computed by the shortcut is not accurate enough then we
should try to use the affine approximation instead.
---

diff --git a/source/grid/manifold_lib.cc b/source/grid/manifold_lib.cc
index 43989a311a..7a0bf39e92 100644
--- a/source/grid/manifold_lib.cc
+++ b/source/grid/manifold_lib.cc
@@ -28,6 +28,12 @@ DEAL_II_NAMESPACE_OPEN
 
 namespace internal
 {
+  // The pull_back function fails regularly in the compute_chart_points
+  // method, and, instead of throwing an exception, returns a point outside
+  // the unit cell. The individual coordinates of that point are given by the
+  // value below.
+  static constexpr double invalid_pull_back_coordinate = 20.0;
+
   // Rotate a given unit vector u around the axis dir
   // where the angle is given by the length of dir.
   // This is the exponential map for a sphere.
@@ -1600,7 +1606,7 @@ TransfiniteInterpolationManifold<dim,spacedim>
 {
   Point<dim> outside;
   for (unsigned int d=0; d<dim; ++d)
-    outside[d] = 20.;
+    outside[d] = internal::invalid_pull_back_coordinate;
 
   // project the user-given input to unit cell
   Point<dim> chart_point = GeometryInfo<dim>::project_to_unit_cell(initial_guess);
@@ -1887,7 +1893,15 @@ TransfiniteInterpolationManifold<dim, spacedim>
                   Assert(false, ExcInternalError());
                 }
 
+              // This guess should be pretty good, but if the pull_back fails
+              // then try again with affine approximation (which is more
+              // expensive)
               chart_points[i] = pull_back(cell, surrounding_points[i], guess);
+              if (chart_points[i][0] == internal::invalid_pull_back_coordinate)
+                {
+                  chart_points[i] = pull_back(cell, surrounding_points[i],
+                                              cell->real_to_unit_cell_affine_approximation(surrounding_points[i]));
+                }
             }
           else
             chart_points[i] = pull_back(cell, surrounding_points[i],