From 60044566173b47475dfe72dc90584f43a641d050 Mon Sep 17 00:00:00 2001
From: Lukas Riedel <mail@lukasriedel.com>
Date: Thu, 16 May 2019 10:34:26 +0000
Subject: [PATCH] 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                          |  2 +-
 dune/dorie/common/interpolator.hh     |  2 +-
 dune/dorie/test/test-interpolators.cc | 44 ++++++++++++++-------------
 3 files changed, 25 insertions(+), 23 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 9311c973..82376fed 100644
--- a/CHANGELOG.md
+++ b/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
 
diff --git a/dune/dorie/common/interpolator.hh b/dune/dorie/common/interpolator.hh
index 6ff3aa10..7441108f 100644
--- a/dune/dorie/common/interpolator.hh
+++ b/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;
 
diff --git a/dune/dorie/test/test-interpolators.cc b/dune/dorie/test/test-interpolators.cc
index c34e8c38..82724780 100644
--- a/dune/dorie/test/test-interpolators.cc
+++ b/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));
 }
-- 
GitLab