Fix distance threshold bug ScalarInterp

src/core/linalg/src/dense/4C_linalg_utils_scalar_interpolation.cpp

@@ -13,6 +13,7 @@
 #include "4C_linalg_fixedsizematrix_solver.hpp"
 #include "4C_utils_exceptions.hpp"
 
+#include <algorithm>
 #include <vector>
 
 FOUR_C_NAMESPACE_OPEN
@@ -199,7 +200,11 @@ std::vector<double> Core::LinAlg::calculate_normalized_weights(
       {
         double norm = ref_locs_tmp[i].norm2();
         if (norm < interp_params.distance_threshold)
+        {
+          std::fill(weights.begin(), weights.end(), 0.0);  // reset all weights to 0.0
           weights[i] = 1.0;
+          break;  // no need to continue
+        }
         else
           weights[i] = 1.0 / std::pow(norm, exponent);
       }
@@ -212,7 +217,11 @@ std::vector<double> Core::LinAlg::calculate_normalized_weights(
       {
         double norm = ref_locs_tmp[i].norm2();
         if (norm < interp_params.distance_threshold)
+        {
+          std::fill(weights.begin(), weights.end(), 0.0);  // reset all weights to 0.0
           weights[i] = 1.0;
+          break;  // no need to continue
+        }
         else
           weights[i] = std::exp(-1.0 * interp_params.exponential_decay_c * norm * norm);
       }
@@ -231,13 +240,11 @@ std::vector<double> Core::LinAlg::calculate_normalized_weights(
       FOUR_C_THROW("Unknown weighting function.");
   }
 
-  // if one of the weights is 1.0: zero out all other weights and return directly
+  // if one of the weights is 1.0: return directly, since the other ones are 0.0 anyway
   for (unsigned int i = 0; i < weights.size(); ++i)
   {
     if (weights[i] == 1.0)
     {
-      for (auto& w : weights) w = 0.0;
-      weights[i] = 1.0;  // set the weight to 1.0 for the closest point
       return weights;
     }
   }
@@ -492,4 +499,4 @@ template class Core::LinAlg::ScalarInterpolator<2, 2, 8>;
 template class Core::LinAlg::ScalarInterpolator<3, 1, 8>;
 template class Core::LinAlg::ScalarInterpolator<3, 2, 10>;
 
-FOUR_C_NAMESPACE_CLOSE
+FOUR_C_NAMESPACE_CLOSE

src/core/linalg/src/dense/4C_linalg_utils_scalar_interpolation.hpp

@@ -46,7 +46,10 @@ namespace Core::LinAlg
     /// (i.e., 1 / (distance^power))
     std::optional<double> inverse_distance_power;
 
-    double distance_threshold = 1.0e-8;  ///< threshold for distance to avoid singularities
+    /// distance threshold below which a reference point is considered coincident with the
+    /// interpolation point; if any reference point is within this threshold, its weight is set
+    /// to 1.0 and all others to 0.0
+    double distance_threshold = 1.0e-8;
   };
 
   /**
@@ -220,4 +223,4 @@ namespace Core::LinAlg
 
 FOUR_C_NAMESPACE_CLOSE
 
-#endif
+#endif

src/core/linalg/tests/4C_linalg_utils_scalar_interpolation_test.cpp

@@ -28,7 +28,7 @@ namespace
   {
     using namespace Core::LinAlg;
 
-    // 1D
+    // Setup dimension of interpolation locations
     constexpr unsigned int loc_dim = 3;
 
     // Reference locations: -1.0, 1.0
@@ -64,11 +64,63 @@ namespace
     for (size_t i = 0; i < weights.size(); ++i) ASSERT_NEAR(weights[i], weights_ref[i], 1e-8);
   }
 
+
+  TEST(LinalgScalarInterpolationTest, CalculateNormalizedWeights1DInverseDistanceAtRefLocations)
+  {
+    using namespace Core::LinAlg;
+
+    // Setup dimension of interpolation locations
+    constexpr unsigned int loc_dim = 3;
+
+    // Reference locations: 0.0, 1.0
+    std::vector<Core::LinAlg::Matrix<loc_dim, 1>> ref_locs;
+    Core::LinAlg::Matrix<loc_dim, 1> loc1(Core::LinAlg::Initialization::zero);
+    Core::LinAlg::Matrix<loc_dim, 1> loc2(Core::LinAlg::Initialization::zero);
+    loc1(0, 0) = 0.0;
+    loc2(0, 0) = 1.0;
+    ref_locs.push_back(loc1);
+    ref_locs.push_back(loc2);
+
+    // Interpolation location: 0.0
+    Core::LinAlg::Matrix<loc_dim, 1> interp_loc(Core::LinAlg::Initialization::zero);
+    interp_loc(0, 0) = 0.0;
+
+    // Weighting function: inverse distance
+    ScalarInterpolationWeightingFunction weight_func =
+        ScalarInterpolationWeightingFunction::inverse_distance;
+
+    // Interpolation parameters
+    ScalarInterpolationParams interp_params;
+    interp_params.distance_threshold = 1e-12;
+    interp_params.inverse_distance_power = 1.0;  // Use power of 1 for inverse distance
+
+    // Call the function
+    std::vector<double> weights = Core::LinAlg::calculate_normalized_weights(
+        ref_locs, interp_loc, weight_func, interp_params);
+
+    // So weights should be [1.0, 0.0]
+    std::vector<double> weights_ref = {1.0, 0.0};
+
+    ASSERT_EQ(weights.size(), 2);
+    for (size_t i = 0; i < weights.size(); ++i) ASSERT_NEAR(weights[i], weights_ref[i], 1e-8);
+
+
+    // And the same for interpolation location: 1.0
+    interp_loc(0, 0) = 1.0;
+    weights = Core::LinAlg::calculate_normalized_weights(
+        ref_locs, interp_loc, weight_func, interp_params);
+    ASSERT_EQ(weights.size(), 2);
+    weights_ref = {0.0, 1.0};
+    for (size_t i = 0; i < weights.size(); ++i) ASSERT_NEAR(weights[i], weights_ref[i], 1e-8);
+  }
+
+
+
   TEST(LinalgScalarInterpolationTest, CalculateNormalizedWeights1DExponential)
   {
     using namespace Core::LinAlg;
 
-    // 1D
+    // Setup dimension of interpolation locations
     constexpr unsigned int loc_dim = 3;
 
     // Reference locations: -1.0, 1.0
@@ -109,6 +161,57 @@ namespace
     for (size_t i = 0; i < weights.size(); ++i) ASSERT_NEAR(weights[i], weights_ref[i], 1e-8);
   }
 
+
+  TEST(LinalgScalarInterpolationTest, CalculateNormalizedWeights1DExponentialAtRefLocations)
+  {
+    using namespace Core::LinAlg;
+
+    // Setup dimension of interpolation locations
+    constexpr unsigned int loc_dim = 3;
+
+    // Reference locations: 0.0, 1.0
+    std::vector<Core::LinAlg::Matrix<loc_dim, 1>> ref_locs;
+    Core::LinAlg::Matrix<loc_dim, 1> loc1(Core::LinAlg::Initialization::zero);
+    Core::LinAlg::Matrix<loc_dim, 1> loc2(Core::LinAlg::Initialization::zero);
+    loc1(0, 0) = 0.0;
+    loc2(0, 0) = 1.0;
+    ref_locs.push_back(loc1);
+    ref_locs.push_back(loc2);
+
+    // Interpolation location: 0.0
+    Core::LinAlg::Matrix<loc_dim, 1> interp_loc(Core::LinAlg::Initialization::zero);
+    interp_loc(0, 0) = 0.0;
+
+    // Weighting function: exponential
+    ScalarInterpolationWeightingFunction weight_func =
+        ScalarInterpolationWeightingFunction::exponential;
+
+    // Interpolation parameters
+    ScalarInterpolationParams interp_params;
+    interp_params.distance_threshold = 1e-12;
+    interp_params.exponential_decay_c = 1.0;  // decay parameter
+
+    // Call the function
+    std::vector<double> weights = Core::LinAlg::calculate_normalized_weights(
+        ref_locs, interp_loc, weight_func, interp_params);
+
+    // Expected weights: 1.0 and 0.0
+    std::vector<double> weights_ref = {1.0, 0.0};
+
+    ASSERT_EQ(weights.size(), 2);
+    for (size_t i = 0; i < weights.size(); ++i) ASSERT_NEAR(weights[i], weights_ref[i], 1e-8);
+
+
+    // And the same for interpolation location: 1.0
+    interp_loc(0, 0) = 1.0;
+    weights = Core::LinAlg::calculate_normalized_weights(
+        ref_locs, interp_loc, weight_func, interp_params);
+    ASSERT_EQ(weights.size(), 2);
+    weights_ref = {0.0, 1.0};
+    for (size_t i = 0; i < weights.size(); ++i) ASSERT_NEAR(weights[i], weights_ref[i], 1e-8);
+  }
+
+
   TEST(LinalgScalarInterpolationTest, LogInterpolation_ManualData)
   {
     // see Satheesh et al., 2024, 10.1002/nme.7373, Section 2.4.4 Comparison of eigenvalue
@@ -437,4 +540,4 @@ namespace
   }
 
 }  // namespace
-FOUR_C_NAMESPACE_CLOSE
+FOUR_C_NAMESPACE_CLOSE