From 5022ca913fbb5703a5c44e1cd45b9ca18154510f Mon Sep 17 00:00:00 2001
From: Martin Kronbichler <kronbichler@lnm.mw.tum.de>
Date: Sun, 1 Nov 2020 10:09:04 +0100
Subject: [PATCH] Transfinite interpolation: Use inverse quadratic
 approximation for pull_back

---
 include/deal.II/fe/mapping_q_internal.h |  8 ++++++-
 include/deal.II/grid/manifold_lib.h     | 19 +++++++++++++++
 source/grid/manifold_lib.cc             | 31 ++++++++++++++++---------
 3 files changed, 46 insertions(+), 12 deletions(-)

diff --git a/include/deal.II/fe/mapping_q_internal.h b/include/deal.II/fe/mapping_q_internal.h
index bb3c28ca59..d4d86f9080 100644
--- a/include/deal.II/fe/mapping_q_internal.h
+++ b/include/deal.II/fe/mapping_q_internal.h
@@ -996,12 +996,18 @@ namespace internal
             }
       }
 
+      /**
+       * Copy constructor.
+       */
+      InverseQuadraticApproximation(const InverseQuadraticApproximation &) =
+        default;
+
       /**
        * Evaluate the quadratic approximation.
        */
       template <typename Number>
       Point<dim, Number>
-      compute(const Point<spacedim, Number> &p)
+      compute(const Point<spacedim, Number> &p) const
       {
         Point<dim, Number> result;
         for (unsigned int d = 0; d < dim; ++d)
diff --git a/include/deal.II/grid/manifold_lib.h b/include/deal.II/grid/manifold_lib.h
index e8b4d86339..ba1fe83e7b 100644
--- a/include/deal.II/grid/manifold_lib.h
+++ b/include/deal.II/grid/manifold_lib.h
@@ -26,6 +26,17 @@
 
 DEAL_II_NAMESPACE_OPEN
 
+// forward declaration
+namespace internal
+{
+  namespace MappingQGenericImplementation
+  {
+    template <int, int>
+    class InverseQuadraticApproximation;
+  }
+} // namespace internal
+
+
 /**
  * Manifold description for a polar coordinate system.
  *
@@ -1133,6 +1144,14 @@ private:
    */
   FlatManifold<dim> chart_manifold;
 
+  /**
+   * A vector of quadratic approximations to the inverse map from real points
+   * to chart points for each of the coarse mesh cells.
+   */
+  std::vector<internal::MappingQGenericImplementation::
+                InverseQuadraticApproximation<dim, spacedim>>
+    quadratic_approximation;
+
   /**
    * The connection to Triangulation::signals::clear that must be reset once
    * this class goes out of scope.
diff --git a/source/grid/manifold_lib.cc b/source/grid/manifold_lib.cc
index a2110f226b..91f03bed62 100644
--- a/source/grid/manifold_lib.cc
+++ b/source/grid/manifold_lib.cc
@@ -17,6 +17,7 @@
 #include <deal.II/base/tensor.h>
 
 #include <deal.II/fe/mapping.h>
+#include <deal.II/fe/mapping_q_internal.h>
 
 #include <deal.II/grid/manifold_lib.h>
 #include <deal.II/grid/tria.h>
@@ -1682,10 +1683,13 @@ TransfiniteInterpolationManifold<dim, spacedim>::initialize(
   });
   level_coarse = triangulation.last()->level();
   coarse_cell_is_flat.resize(triangulation.n_cells(level_coarse), false);
-  typename Triangulation<dim, spacedim>::active_cell_iterator
-    cell = triangulation.begin(level_coarse),
-    endc = triangulation.end(level_coarse);
-  for (; cell != endc; ++cell)
+  quadratic_approximation.clear();
+
+  std::vector<Point<dim>> unit_points =
+    QIterated<dim>(QTrapez<1>(), 2).get_points();
+  std::vector<Point<spacedim>> real_points(unit_points.size());
+
+  for (const auto &cell : triangulation.active_cell_iterators())
     {
       bool cell_is_flat = true;
       for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
@@ -1700,6 +1704,11 @@ TransfiniteInterpolationManifold<dim, spacedim>::initialize(
       AssertIndexRange(static_cast<unsigned int>(cell->index()),
                        coarse_cell_is_flat.size());
       coarse_cell_is_flat[cell->index()] = cell_is_flat;
+
+      // build quadratic interpolation
+      for (unsigned int i = 0; i < unit_points.size(); ++i)
+        real_points[i] = push_forward(cell, unit_points[i]);
+      quadratic_approximation.emplace_back(real_points, unit_points);
     }
 }
 
@@ -2281,7 +2290,7 @@ TransfiniteInterpolationManifold<dim, spacedim>::
       for (unsigned int i = 0; i < points.size(); ++i)
         {
           Point<dim> point =
-            cell->real_to_unit_cell_affine_approximation(points[i]);
+            quadratic_approximation[cell->index()].compute(points[i]);
           current_distance += GeometryInfo<dim>::distance_to_unit_cell(point);
         }
       distances_and_cells.push_back(
@@ -2434,11 +2443,11 @@ TransfiniteInterpolationManifold<dim, spacedim>::compute_chart_points(
     [&](const typename Triangulation<dim, spacedim>::cell_iterator &cell,
         const unsigned int point_index) {
       Point<dim> guess;
-      // an optimization: keep track of whether or not we used the affine
+      // an optimization: keep track of whether or not we used the quadratic
       // approximation so that we don't call pull_back with the same
       // initial guess twice (i.e., if pull_back fails the first time,
       // don't try again with the same function arguments).
-      bool used_affine_approximation = false;
+      bool used_quadratic_approximation = false;
       // if we have already computed three points, we can guess the fourth
       // to be the missing corner point of a rectangle
       if (point_index == 3 && surrounding_points.size() >= 8)
@@ -2468,9 +2477,9 @@ TransfiniteInterpolationManifold<dim, spacedim>::compute_chart_points(
         }
       else
         {
-          guess = cell->real_to_unit_cell_affine_approximation(
+          guess = quadratic_approximation[cell->index()].compute(
             surrounding_points[point_index]);
-          used_affine_approximation = true;
+          used_quadratic_approximation = true;
         }
       chart_points[point_index] =
         pull_back(cell, surrounding_points[point_index], guess);
@@ -2480,9 +2489,9 @@ TransfiniteInterpolationManifold<dim, spacedim>::compute_chart_points(
       // than the cheap methods used above)
       if (chart_points[point_index][0] ==
             internal::invalid_pull_back_coordinate &&
-          !used_affine_approximation)
+          !used_quadratic_approximation)
         {
-          guess = cell->real_to_unit_cell_affine_approximation(
+          guess = quadratic_approximation[cell->index()].compute(
             surrounding_points[point_index]);
           chart_points[point_index] =
             pull_back(cell, surrounding_points[point_index], guess);
-- 
2.39.5