From 45193dbfab25c131757c14dc2c1f009f15a6d0cf Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Tue, 7 May 2019 19:13:49 +0200 Subject: [PATCH] extend HSFC to handle degenerate cases --- source/base/utilities.cc | 93 +++++++++++++++++++++++++++++++--------- 1 file changed, 73 insertions(+), 20 deletions(-) diff --git a/source/base/utilities.cc b/source/base/utilities.cc index 17641f7175..ed21748cc0 100644 --- a/source/base/utilities.cc +++ b/source/base/utilities.cc @@ -109,21 +109,23 @@ namespace Utilities { std::vector> 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(points[i][d]) - - static_cast(bl[d])) / - extents[d]; - Assert(v >= 0. && v <= 1., ExcInternalError()); - int_points[i][d] = - static_cast(v * static_cast(max_int)); - } + if (valid_extents[d]) + { + Assert(extents[d] > 0, ExcInternalError()); + const LongDouble v = (static_cast(points[i][d]) - + static_cast(bl[d])) / + extents[d]; + Assert(v >= 0. && v <= 1., ExcInternalError()); + AssertIndexRange(eff_d, effective_dim); + int_points[i][eff_d] = + static_cast(v * static_cast(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::max() : (Integer(1) << min_bits) - 1); - return inverse_Hilbert_space_filling_curve_effective( - 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( + points, bl, extents, valid_extents, min_bits, max_int); + } + + // various degenerate cases + std::array zero_ind; + for (unsigned int d = 0; d < dim; ++d) + zero_ind[d] = 0; + + std::vector> 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( + 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( + 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; } -- 2.39.5