extend HSFC to handle degenerate cases

diff --git a/source/base/utilities.cc b/source/base/utilities.cc
index 17641f71758a038c5509bdac8ce8514fc1f5af10..ed21748cc0a75417888e3220cd67b9f0d815279e 100644 (file)
--- a/source/base/utilities.cc
+++ b/source/base/utilities.cc
@@ -109,21 +109,23 @@ namespace Utilities
     {
       std::vector<std::array<Integer, effective_dim>> int_points(points.size());
 
-      static_assert(dim == effective_dim, "Not implemented yet");
-      (void)valid_extents;
-
       for (unsigned int i = 0; i < points.size(); ++i)
         {
           // convert into integers:
+          unsigned int eff_d = 0;
           for (unsigned int d = 0; d < dim; ++d)
-            {
-              const LongDouble v = (static_cast<LongDouble>(points[i][d]) -
-                                    static_cast<LongDouble>(bl[d])) /
-                                   extents[d];
-              Assert(v >= 0. && v <= 1., ExcInternalError());
-              int_points[i][d] =
-                static_cast<Integer>(v * static_cast<LongDouble>(max_int));
-            }
+            if (valid_extents[d])
+              {
+                Assert(extents[d] > 0, ExcInternalError());
+                const LongDouble v = (static_cast<LongDouble>(points[i][d]) -
+                                      static_cast<LongDouble>(bl[d])) /
+                                     extents[d];
+                Assert(v >= 0. && v <= 1., ExcInternalError());
+                AssertIndexRange(eff_d, effective_dim);
+                int_points[i][eff_d] =
+                  static_cast<Integer>(v * static_cast<LongDouble>(max_int));
+                ++eff_d;
+              }
         }
 
       // note that we call this with "min_bits"
@@ -165,9 +167,6 @@ namespace Utilities
         valid_extents[i] = (extents[i] > 0.);
       }
 
-    Assert(valid_extents.any(),
-           ExcMessage("Bounding box is degenerate in all dimensions."));
-
     // make sure our conversion from fractional coordinates to
     // Integers work as expected, namely our cast (LongDouble)max_int
     const int min_bits =
@@ -180,12 +179,66 @@ namespace Utilities
                                std::numeric_limits<Integer>::max() :
                                (Integer(1) << min_bits) - 1);
 
-    return inverse_Hilbert_space_filling_curve_effective<dim,
-                                                         Number,
-                                                         dim,
-                                                         LongDouble,
-                                                         Integer>(
-      points, bl, extents, valid_extents, min_bits, max_int);
+    const unsigned int effective_dim = valid_extents.count();
+    if (effective_dim == dim)
+      {
+        return inverse_Hilbert_space_filling_curve_effective<dim,
+                                                             Number,
+                                                             dim,
+                                                             LongDouble,
+                                                             Integer>(
+          points, bl, extents, valid_extents, min_bits, max_int);
+      }
+
+    // various degenerate cases
+    std::array<std::uint64_t, dim> zero_ind;
+    for (unsigned int d = 0; d < dim; ++d)
+      zero_ind[d] = 0;
+
+    std::vector<std::array<std::uint64_t, dim>> ind(points.size(), zero_ind);
+    // manually check effective_dim == 1 and effective_dim == 2
+    if (dim == 3 && effective_dim == 2)
+      {
+        const auto ind2 =
+          inverse_Hilbert_space_filling_curve_effective<dim,
+                                                        Number,
+                                                        2,
+                                                        LongDouble,
+                                                        Integer>(
+            points, bl, extents, valid_extents, min_bits, max_int);
+
+        for (unsigned int i = 0; i < ind.size(); ++i)
+          for (unsigned int d = 0; d < 2; ++d)
+            ind[i][d + 1] = ind2[i][d];
+
+        return ind;
+      }
+    else if (effective_dim == 1)
+      {
+        const auto ind1 =
+          inverse_Hilbert_space_filling_curve_effective<dim,
+                                                        Number,
+                                                        1,
+                                                        LongDouble,
+                                                        Integer>(
+            points, bl, extents, valid_extents, min_bits, max_int);
+
+        for (unsigned int i = 0; i < ind.size(); ++i)
+          ind[i][dim - 1] = ind1[i][0];
+
+        return ind;
+      }
+
+    // we should get here only if effective_dim == 0
+    Assert(effective_dim == 0, ExcInternalError());
+
+    // if the bounding box is degenerate in all dimensions,
+    // can't do much but exit gracefully by setting index according
+    // to the index of each point so that there is no re-ordering
+    for (unsigned int i = 0; i < points.size(); ++i)
+      ind[i][dim - 1] = i;
+
+    return ind;
   }