From 7a2c1d182e9d4e4dc20237acb96d5449cbecc0ce Mon Sep 17 00:00:00 2001
From: David Wells <wellsd2@rpi.edu>
Date: Thu, 18 Jan 2018 13:37:43 -0500
Subject: [PATCH] Implement a 2D-specific chart point estimate.

The current algorithm calculates a reasonable guess for the chart point
in 3D. Using this 3D algorithm in 2D results in relatively bad guesses:

for i = 5:
affine approx chart point: 0.5 3.64977e-14
pulled back chart point:   1 0.5
for i = 6:
affine approx chart point: 0.5 0.5
pulled back chart point:   0.5 8.678e-17
for i = 7:
affine approx chart point: 0.5 1
pulled back chart point:   0.5 1

This patch adds a 2D-specific algorithm that simply averages the two
adjacent vertices (e.g., chart point 5 is assumed to be an average of
chart points 1 and 3).

This issue was reported by Juan Carlos Araujo Cabarcas, who found a test
case where the affine guess was sufficiently inaccurate that
TransfiniteInterpolationManifold::pull_back failed to find the right
point. This change fixes that test case.
---
 source/grid/manifold_lib.cc | 66 ++++++++++++++++++++++++++++++++++---
 1 file changed, 61 insertions(+), 5 deletions(-)

diff --git a/source/grid/manifold_lib.cc b/source/grid/manifold_lib.cc
index 48f2878e12..43989a311a 100644
--- a/source/grid/manifold_lib.cc
+++ b/source/grid/manifold_lib.cc
@@ -1825,12 +1825,68 @@ TransfiniteInterpolationManifold<dim, spacedim>
               chart_points[i] = pull_back(cell, surrounding_points[i], p3);
             }
           // 8 points usually form either a cube or a rectangle with vertices
-          // and line mid points. assume a cube here which gives us some new
-          // initial guesses
-          else if (surrounding_points.size() == 8 && i > 4)
+          // and line mid points. Get the initial guess with line segment
+          // midpoints in 2D and assuming a cube for 3D.
+          else if (surrounding_points.size() == 8 &&
+                   ((dim == 3 && i > 4) || (dim == 2 && i > 3)))
             {
-              const Point<dim> guess = chart_points[i-4] +
-                                       (chart_points[4] - chart_points[0]);
+              Point<dim> guess;
+              switch (dim)
+                {
+                case 2:
+                  // inline the standard numbering
+                  //
+                  // 2 - 7 - 3
+                  // |       |
+                  // 4       5
+                  // |       |
+                  // 0 - 6 - 1
+                  //
+                  // to calculate guesses based on averaging already computed
+                  // chart points.
+                  switch (i)
+                    {
+                    case 4:
+                      guess = 0.5*(chart_points[0] + chart_points[2]);
+                      break;
+                    case 5:
+                      guess = 0.5*(chart_points[1] + chart_points[3]);
+                      break;
+                    case 6:
+                      guess = 0.5*(chart_points[0] + chart_points[1]);
+                      break;
+                    case 7:
+                      guess = 0.5*(chart_points[2] + chart_points[3]);
+                      break;
+                    default:
+                      Assert(false, ExcInternalError());
+                    }
+                  break;
+                case 3:
+                  // Assuming that we are in 3D and have the points around a
+                  // cube numbered as
+                  //
+                  //         6-------7
+                  //        /|      /|
+                  //       /       / |
+                  //      /  |    /  |
+                  //     4-------5   |
+                  //     |   2- -|- -3
+                  //     |  /    |  /
+                  //     |       | /
+                  //     |/      |/
+                  //     0-------1
+                  //
+                  // (where vertex 2 is the back left vertex) we can estimate
+                  // where chart points 5 - 7 are by computing the height (in
+                  // chart coordinates) as c4 - c0 and then adding that onto the
+                  // appropriate bottom vertex.
+                  guess = chart_points[i - 4] + (chart_points[4] - chart_points[0]);
+                  break;
+                default:
+                  Assert(false, ExcInternalError());
+                }
+
               chart_points[i] = pull_back(cell, surrounding_points[i], guess);
             }
           else
-- 
2.39.5