Commit 87b24ba9 authored by Lukas Riedel's avatar Lukas Riedel

Merge branch 'master' into 127-add-data-assimilation-interface-to-richardssimulation

parents da420aed 272e78ff
......@@ -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
......
......@@ -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;
......
......@@ -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));
}
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment