Use correct coordinates in linear interpolator

Fetch the enclosing multi indices with the coordinates in mesh units,
not physical units. Improve the unit test to cover more complicated
cases, thus ensuring that the interpolator works as intended.

CHANGELOG.md

@@ -40,7 +40,7 @@
 * Coupling between transient models for water flow and solute transport !96
 * Initial conditions generated from H5 input data !130
 * Generic Python VTK file reader !143
-* Linear interpolator for initial conditions and scaling fields !145
+* Linear interpolator for initial conditions and scaling fields !145, !156
 * Parameterizations for hydrodynamic dispersion in solute transport !141
 * Generic Python VTK file reader !143, !150
 

dune/dorie/common/interpolator.hh

@@ -267,7 +267,7 @@ public:
 
         // compute indices and position vectors in mesh units
         const auto pos_unit = get_pos_unit(pos);
-        const auto [idx_lower, idx_upper] = get_knot_indices(pos);
+        const auto [idx_lower, idx_upper] = get_knot_indices(pos_unit);
         const auto pos_unit_diff = get_pos_unit_diff(pos_unit, idx_lower);
         const auto& data = this->_data;
 

dune/dorie/test/test-interpolators.cc

@@ -138,10 +138,12 @@ template<>
 void test_linear<2> ()
 {
     constexpr int dim = 2;
-    std::vector<double> data({0.0, 1.0, 1.0, 2.0});
+    std::vector<double> data({0.0, 0.0, 1.0,
+                              1.0, 1.0, 2.0,
+                              3.0, 3.0, 4.0});
 
     std::vector<size_t> shape(dim);
-    std::fill(begin(shape), end(shape), 2);
+    std::fill(begin(shape), end(shape), 3);
 
     Dune::FieldVector<double, dim> extensions(1.0);
     Dune::FieldVector<double, dim> offset(0.0);
@@ -156,22 +158,25 @@ void test_linear<2> ()
     
     // check without offset
     using Dune::FloatCmp::eq; // floating-point comparison
-    std::vector<Dune::FieldVector<double, dim>> points(5);
+    std::vector<Dune::FieldVector<double, dim>> points(7);
 
     points[0] = {0.0, 0.0};
     points[1] = {1.0, 0.0};
     points[2] = {0.0, 1.0};
     points[3] = {1.0, 1.0};
-    points[4] = {0.75, 0.75};
+    points[4] = {0.5, 0.5};
+    points[5] = {0.25, 0.5};
+    points[6] = {1.0, 0.75};
     assert(eq(interp->evaluate(points[0]), 0.0));
     assert(eq(interp->evaluate(points[1]), 1.0));
-    assert(eq(interp->evaluate(points[2]), 1.0));
-    assert(eq(interp->evaluate(points[3]), 2.0));
-    assert(eq(interp->evaluate(points[4]), 1.5));
-
+    assert(eq(interp->evaluate(points[2]), 3.0));
+    assert(eq(interp->evaluate(points[3]), 4.0));
+    assert(eq(interp->evaluate(points[4]), 1.0));
+    assert(eq(interp->evaluate(points[5]), 1.0));
+    assert(eq(interp->evaluate(points[6]), 3.0));
 
     // check with offset
-    std::fill(offset.begin(), offset.end(), -1.0);
+    std::fill(offset.begin(), offset.end(), -0.5);
     interp = Dune::Dorie::InterpolatorFactory<InterpolatorTraits<dim>>
         ::create("linear",
                  data,
@@ -179,11 +184,8 @@ void test_linear<2> ()
                  extensions,
                  offset);
 
-    points.resize(2);
-    points[0] = {0.0, 0.0};
-    points[1] = {-0.25, -0.25};
-    assert(eq(interp->evaluate(points[0]), 2.0));
-    assert(eq(interp->evaluate(points[1]), 1.5));
+    assert(eq(interp->evaluate(points[0]), 1.0));
+    assert(eq(interp->evaluate(points[5]), 3.5));
 }
 
 /// Test the linear interpolator in 3D
@@ -197,7 +199,7 @@ void test_linear<3> ()
     std::vector<size_t> shape(dim);
     std::fill(begin(shape), end(shape), 2);
 
-    Dune::FieldVector<double, dim> extensions(1.0);
+    Dune::FieldVector<double, dim> extensions(2.0);
     Dune::FieldVector<double, dim> offset(0.0);
 
     // build interpolator
@@ -213,10 +215,10 @@ void test_linear<3> ()
     std::vector<Dune::FieldVector<double, dim>> points(5);
 
     points[0] = {0.0, 0.0, 0.0};
-    points[1] = {1.0, 0.0, 0.0};
-    points[2] = {0.0, 0.0, 1.0};
-    points[3] = {1.0, 1.0, 1.0};
-    points[4] = {0.5, 0.5, 1.0};
+    points[1] = {2.0, 0.0, 0.0};
+    points[2] = {0.0, 0.0, 2.0};
+    points[3] = {2.0, 2.0, 2.0};
+    points[4] = {1.0, 1.0, 2.0};
     assert(eq(interp->evaluate(points[0]), 0.0));
     assert(eq(interp->evaluate(points[1]), 0.0));
     assert(eq(interp->evaluate(points[2]), 1.0));
@@ -224,7 +226,7 @@ void test_linear<3> ()
     assert(eq(interp->evaluate(points[4]), 1.5));
 
     // check with offset
-    std::fill(offset.begin(), offset.end(), -1.0);
+    std::fill(offset.begin(), offset.end(), -2.0);
     interp = Dune::Dorie::InterpolatorFactory<InterpolatorTraits<dim>>
         ::create("linear",
                  data,
@@ -234,7 +236,7 @@ void test_linear<3> ()
 
     points.resize(2);
     points[0] = {0.0, 0.0, 0.0};
-    points[1] = {-0.5, -0.5, 0.0};
+    points[1] = {-1.0, -1.0, 0.0};
     assert(eq(interp->evaluate(points[0]), 2.0));
     assert(eq(interp->evaluate(points[1]), 1.5));
 }