Merge branch...

Merge branch '154-linear-interpolator-only-works-for-minimal-shapes-and-trivial-extensions' into 'master'

Resolve "Linear interpolator only works for minimal shapes and trivial extensions"

Closes #154

See merge request !156

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