TransfiniteInterpolationManifold: fix codimension one case

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));

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())

tests/manifold/transfinite_manifold_12.cc

@@ -0,0 +1,44 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 - 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+// This test verifies that the transfinite interpolation can be initialized
+// in the codimension one case.
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/tria.h>
+
+#include <memory>
+
+#include "../tests.h"
+
+int
+main()
+{
+  initlog();
+
+  using namespace dealii;
+
+  Triangulation<2, 3> tria;
+
+  GridGenerator::hyper_cube(tria, -1, 1);
+  TransfiniteInterpolationManifold<2, 3> transfinite;
+  transfinite.initialize(tria);
+
+  deallog << "OK" << std::endl;
+
+  return 0;
+}

tests/manifold/transfinite_manifold_12.debug.output

@@ -0,0 +1,2 @@
+
+DEAL::OK