From bddd1966f8fc2ddac7dce358c5c1bfd947d1b605 Mon Sep 17 00:00:00 2001
From: Matthias Maier <tamiko@43-1.org>
Date: Mon, 16 May 2022 12:36:16 -0500
Subject: [PATCH] TransfiniteInterpolationManifold: fix codimension one case

The codimension one case currently fails with an assertion:
```
An error occurred in line <290> of file <source/grid/grid_tools.cc> in function
    std::pair<DerivativeForm<1, dim, spacedim>, Tensor<1, spacedim>> dealii::GridTools::affine_cell_approximation(const ArrayView<const Point<spacedim>> &) [dim = 2, spacedim = 3]
The violated condition was:
    ::dealii::deal_II_exceptions::internals::compare_for_equality(vertices.size(), GeometryInfo<dim>::vertices_per_cell)
Additional information:
    Dimension 9 not equal to 4.
```
This is triggered because we take a compatibility branch in case of
`dim < spacedim` (in `include/deal.II/fe/mapping_q_internal.h`):
```
908         if (real_support_points.size() ==
909               GeometryInfo<dim>::vertices_per_cell ||
910             dim < spacedim)
911           {
```
but still initialize the vector `unit_points` by a subdivided quadrature
(which is needed for the default quadratic approximation; in
`source/grid/manifold_lib.cc`):
```
1686   std::vector<Point<dim>> unit_points =
1687     QIterated<dim>(QTrapez<1>(), 2).get_points();
1688   std::vector<Point<spacedim>> real_points(unit_points.size());
```

Fix this by simply initializing `unit_points` with the right number of
interpolation points for the `dim < spacedim` variant.
---
 include/deal.II/fe/mapping_q_internal.h | 18 +++++++++++-------
 source/grid/manifold_lib.cc             | 12 ++++++++++--
 2 files changed, 21 insertions(+), 9 deletions(-)

diff --git a/include/deal.II/fe/mapping_q_internal.h b/include/deal.II/fe/mapping_q_internal.h
index ade514c553..2174e1480b 100644
--- a/include/deal.II/fe/mapping_q_internal.h
+++ b/include/deal.II/fe/mapping_q_internal.h
@@ -901,13 +901,17 @@ namespace internal
         AssertDimension(real_support_points.size(), unit_support_points.size());
 
         // For the bi-/trilinear approximation, we cannot build a quadratic
-        // polynomial due to a lack of points (interpolation matrix would get
-        // singular), so pick the affine approximation. Similarly, it is not
-        // entirely clear how to gather enough information for the case dim <
-        // spacedim
-        if (real_support_points.size() ==
-              GeometryInfo<dim>::vertices_per_cell ||
-            dim < spacedim)
+        // polynomial due to a lack of points (interpolation matrix would
+        // get singular). Similarly, it is not entirely clear how to gather
+        // enough information for the case dim < spacedim.
+        //
+        // In both cases we require the vector real_support_points to
+        // contain the vertex positions and fall back to an affine
+        // approximation:
+        Assert(dim == spacedim || real_support_points.size() ==
+                                    GeometryInfo<dim>::vertices_per_cell,
+               ExcInternalError());
+        if (real_support_points.size() == GeometryInfo<dim>::vertices_per_cell)
           {
             const auto affine = GridTools::affine_cell_approximation<dim>(
               make_array_view(real_support_points));
diff --git a/source/grid/manifold_lib.cc b/source/grid/manifold_lib.cc
index 1b9ad3730c..085f8e8248 100644
--- a/source/grid/manifold_lib.cc
+++ b/source/grid/manifold_lib.cc
@@ -1673,7 +1673,7 @@ TransfiniteInterpolationManifold<dim, spacedim>::initialize(
   const Triangulation<dim, spacedim> &triangulation)
 {
   this->triangulation = &triangulation;
-  // in case the triangulation is cleared, remove the pointers by a signal
+  // In case the triangulation is cleared, remove the pointers by a signal:
   clear_signal.disconnect();
   clear_signal = triangulation.signals.clear.connect([&]() -> void {
     this->triangulation = nullptr;
@@ -1683,8 +1683,16 @@ TransfiniteInterpolationManifold<dim, spacedim>::initialize(
   coarse_cell_is_flat.resize(triangulation.n_cells(level_coarse), false);
   quadratic_approximation.clear();
 
+  // In case of dim == spacedim we perform a quadratic approximation in
+  // InverseQuadraticApproximation(), thus initialize the unit_points
+  // vector with one subdivision to get 3^dim unit_points.
+  //
+  // In the co-dimension one case (meaning  dim < spacedim) we have to fall
+  // back to a simple GridTools::affine_cell_approximation<dim>() which
+  // requires 2^dim points, instead. Thus, initialize the QIteraded
+  // quadrature with no subdivisions.
   std::vector<Point<dim>> unit_points =
-    QIterated<dim>(QTrapez<1>(), 2).get_points();
+    QIterated<dim>(QTrapez<1>(), (dim == spacedim ? 2 : 1)).get_points();
   std::vector<Point<spacedim>> real_points(unit_points.size());
 
   for (const auto &cell : triangulation.active_cell_iterators())
-- 
2.39.5