test-interpolators.cc 8.27 KB
Newer Older
Lukas Riedel's avatar
Lukas Riedel committed
1 2 3 4 5 6 7 8 9 10 11 12
#ifdef HAVE_CONFIG_H
    #include "config.h"
#endif

#include <cassert>
#include <exception>
#include <vector>
#include <cmath>
#include <numeric>
#include <iostream>

#include <dune/common/exceptions.hh>
13
#include <dune/common/float_cmp.hh>
Lukas Riedel's avatar
Lukas Riedel committed
14
#include <dune/common/parallel/mpihelper.hh>
Lukas Riedel's avatar
Lukas Riedel committed
15
#include <dune/common/fvector.hh>
Lukas Riedel's avatar
Lukas Riedel committed
16

17
#include <dune/dorie/common/interpolator.hh>
Lukas Riedel's avatar
Lukas Riedel committed
18 19 20 21 22

template<int dimension>
struct InterpolatorTraits
{
    static constexpr int dim = dimension;
Lukas Riedel's avatar
Lukas Riedel committed
23
    using Domain = Dune::FieldVector<double, dim>;
Lukas Riedel's avatar
Lukas Riedel committed
24
    using DF = double;
25
    using RF = double;
Lukas Riedel's avatar
Lukas Riedel committed
26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42
};

/// Test the nearest neighbor interpolator.
/** 
 *  * Create a 3x3 grid with extensions 3x3.
 *  * Offset by -1x-1.
 *  * Check values.
 */
template<int dim=2>
void test_nearest_neighbor ()
{
    std::vector<double> data(std::pow(3, dim));
    std::iota(begin(data), end(data), 0.0);

    std::vector<size_t> shape(dim);
    std::fill(begin(shape), end(shape), 3);

Lukas Riedel's avatar
Lukas Riedel committed
43 44
    Dune::FieldVector<double, dim> extensions(3.0);
    // std::fill(begin(extensions), end(extensions), 3.0);
Lukas Riedel's avatar
Lukas Riedel committed
45

Lukas Riedel's avatar
Lukas Riedel committed
46 47
    Dune::FieldVector<double, dim> offset(0.0);
    // std::fill(begin(offset), end(offset), 0.0);
Lukas Riedel's avatar
Lukas Riedel committed
48 49 50

    // build interpolator
    auto interp = Dune::Dorie::InterpolatorFactory<InterpolatorTraits<dim>>
51 52 53 54 55
        ::create("nearest",
                 data,
                 shape,
                 extensions,
                 offset);
Lukas Riedel's avatar
Lukas Riedel committed
56 57
    
    // check without offset
58
    using Dune::FloatCmp::eq; // floating-point comparison
Lukas Riedel's avatar
Lukas Riedel committed
59
    std::vector<Dune::FieldVector<double, dim>> corners;
Lukas Riedel's avatar
Lukas Riedel committed
60 61 62 63 64 65
    if constexpr (dim == 2) {
        corners.resize(4);
        corners[0] = {0.0, 0.0};
        corners[1] = {1.0, 0.0};
        corners[2] = {0.0, 1.0};
        corners[3] = {1.0, 1.0};
66 67 68 69
        assert(eq(interp->evaluate(corners[0]), 0.0));
        assert(eq(interp->evaluate(corners[1]), 1.0));
        assert(eq(interp->evaluate(corners[2]), 3.0));
        assert(eq(interp->evaluate(corners[3]), 4.0));
Lukas Riedel's avatar
Lukas Riedel committed
70 71 72 73 74 75 76 77 78 79 80
    }
    else if (dim == 3) {
        corners.resize(8);
        corners[0] = {0.0, 0.0, 0.0};
        corners[1] = {1.0, 0.0, 0.0};
        corners[2] = {0.0, 1.0, 0.0};
        corners[3] = {1.0, 1.0, 0.0};
        corners[4] = {0.0, 0.0, 1.0};
        corners[5] = {1.0, 0.0, 1.0};
        corners[6] = {0.0, 1.0, 1.0};
        corners[7] = {1.0, 1.0, 1.0};
81 82 83 84 85 86 87 88
        assert(eq(interp->evaluate(corners[0]), 0.0));
        assert(eq(interp->evaluate(corners[1]), 1.0));
        assert(eq(interp->evaluate(corners[2]), 3.0));
        assert(eq(interp->evaluate(corners[3]), 4.0));
        assert(eq(interp->evaluate(corners[4]), 9.0));
        assert(eq(interp->evaluate(corners[5]), 10.0));
        assert(eq(interp->evaluate(corners[6]), 12.0));
        assert(eq(interp->evaluate(corners[7]), 13.0));
Lukas Riedel's avatar
Lukas Riedel committed
89 90 91
    }

    // check with offset
Lukas Riedel's avatar
Lukas Riedel committed
92
    std::fill(offset.begin(), offset.end(), -1.0);
Lukas Riedel's avatar
Lukas Riedel committed
93
    interp = Dune::Dorie::InterpolatorFactory<InterpolatorTraits<dim>>
94 95 96 97 98
        ::create("nearest",
                 data,
                 shape,
                 extensions,
                 offset);
Lukas Riedel's avatar
Lukas Riedel committed
99 100 101 102 103 104 105

    if constexpr (dim == 2) {
        corners.resize(4);
        corners[0] = {0.0, 0.0};
        corners[1] = {1.0, 0.0};
        corners[2] = {0.0, 1.0};
        corners[3] = {1.0, 1.0};
106 107 108 109
        assert(eq(interp->evaluate(corners[0]), 4.0));
        assert(eq(interp->evaluate(corners[1]), 5.0));
        assert(eq(interp->evaluate(corners[2]), 7.0));
        assert(eq(interp->evaluate(corners[3]), 8.0));
Lukas Riedel's avatar
Lukas Riedel committed
110 111 112 113 114 115 116 117 118 119 120
    }
    else if (dim == 3) {
        corners.resize(8);
        corners[0] = {0.0, 0.0, 0.0};
        corners[1] = {1.0, 0.0, 0.0};
        corners[2] = {0.0, 1.0, 0.0};
        corners[3] = {1.0, 1.0, 0.0};
        corners[4] = {0.0, 0.0, 1.0};
        corners[5] = {1.0, 0.0, 1.0};
        corners[6] = {0.0, 1.0, 1.0};
        corners[7] = {1.0, 1.0, 1.0};
121 122 123 124 125 126 127 128
        assert((interp->evaluate(corners[0]), 13.0));
        assert((interp->evaluate(corners[1]), 14.0));
        assert((interp->evaluate(corners[2]), 16.0));
        assert((interp->evaluate(corners[3]), 17.0));
        assert((interp->evaluate(corners[4]), 22.0));
        assert((interp->evaluate(corners[5]), 23.0));
        assert((interp->evaluate(corners[6]), 25.0));
        assert((interp->evaluate(corners[7]), 26.0));
Lukas Riedel's avatar
Lukas Riedel committed
129 130 131
    }
}

Lukas Riedel's avatar
Lukas Riedel committed
132 133 134 135 136 137 138 139 140
/// Test the linear interpolator
template<int dim>
void test_linear ();

/// Test the linear interpolator in 2D
template<>
void test_linear<2> ()
{
    constexpr int dim = 2;
141 142 143
    std::vector<double> data({0.0, 0.0, 1.0,
                              1.0, 1.0, 2.0,
                              3.0, 3.0, 4.0});
Lukas Riedel's avatar
Lukas Riedel committed
144 145

    std::vector<size_t> shape(dim);
146
    std::fill(begin(shape), end(shape), 3);
Lukas Riedel's avatar
Lukas Riedel committed
147 148 149 150 151 152 153 154 155 156 157 158 159 160

    Dune::FieldVector<double, dim> extensions(1.0);
    Dune::FieldVector<double, dim> offset(0.0);

    // build interpolator
    auto interp = Dune::Dorie::InterpolatorFactory<InterpolatorTraits<dim>>
        ::create("linear",
                 data,
                 shape,
                 extensions,
                 offset);
    
    // check without offset
    using Dune::FloatCmp::eq; // floating-point comparison
161
    std::vector<Dune::FieldVector<double, dim>> points(7);
Lukas Riedel's avatar
Lukas Riedel committed
162 163 164 165 166

    points[0] = {0.0, 0.0};
    points[1] = {1.0, 0.0};
    points[2] = {0.0, 1.0};
    points[3] = {1.0, 1.0};
167 168 169
    points[4] = {0.5, 0.5};
    points[5] = {0.25, 0.5};
    points[6] = {1.0, 0.75};
Lukas Riedel's avatar
Lukas Riedel committed
170 171
    assert(eq(interp->evaluate(points[0]), 0.0));
    assert(eq(interp->evaluate(points[1]), 1.0));
172 173 174 175 176
    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));
Lukas Riedel's avatar
Lukas Riedel committed
177 178

    // check with offset
179
    std::fill(offset.begin(), offset.end(), -0.5);
Lukas Riedel's avatar
Lukas Riedel committed
180 181 182 183 184 185 186
    interp = Dune::Dorie::InterpolatorFactory<InterpolatorTraits<dim>>
        ::create("linear",
                 data,
                 shape,
                 extensions,
                 offset);

187 188
    assert(eq(interp->evaluate(points[0]), 1.0));
    assert(eq(interp->evaluate(points[5]), 3.5));
Lukas Riedel's avatar
Lukas Riedel committed
189 190 191 192 193 194 195 196 197 198 199 200 201
}

/// Test the linear interpolator in 3D
template<>
void test_linear<3> ()
{
    constexpr int dim = 3;
    std::vector<double> data({0.0, 0.0, 0.0, 0.0,
                              1.0, 1.0, 2.0, 2.0});

    std::vector<size_t> shape(dim);
    std::fill(begin(shape), end(shape), 2);

202
    Dune::FieldVector<double, dim> extensions(2.0);
Lukas Riedel's avatar
Lukas Riedel committed
203 204 205 206 207 208 209 210 211 212 213 214 215 216 217
    Dune::FieldVector<double, dim> offset(0.0);

    // build interpolator
    auto interp = Dune::Dorie::InterpolatorFactory<InterpolatorTraits<dim>>
        ::create("linear",
                 data,
                 shape,
                 extensions,
                 offset);

    // check without offset
    using Dune::FloatCmp::eq; // floating-point comparison
    std::vector<Dune::FieldVector<double, dim>> points(5);

    points[0] = {0.0, 0.0, 0.0};
218 219 220 221
    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};
Lukas Riedel's avatar
Lukas Riedel committed
222 223 224 225 226 227 228
    assert(eq(interp->evaluate(points[0]), 0.0));
    assert(eq(interp->evaluate(points[1]), 0.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));

    // check with offset
229
    std::fill(offset.begin(), offset.end(), -2.0);
Lukas Riedel's avatar
Lukas Riedel committed
230 231 232 233 234 235 236 237 238
    interp = Dune::Dorie::InterpolatorFactory<InterpolatorTraits<dim>>
        ::create("linear",
                 data,
                 shape,
                 extensions,
                 offset);

    points.resize(2);
    points[0] = {0.0, 0.0, 0.0};
239
    points[1] = {-1.0, -1.0, 0.0};
Lukas Riedel's avatar
Lukas Riedel committed
240 241 242 243
    assert(eq(interp->evaluate(points[0]), 2.0));
    assert(eq(interp->evaluate(points[1]), 1.5));
}

Lukas Riedel's avatar
Lukas Riedel committed
244 245 246 247
int main (int argc, char** argv)
{
    try{
        // initialize MPI if needed
248
        auto& helper = Dune::MPIHelper::instance(argc, argv);
Lukas Riedel's avatar
Lukas Riedel committed
249 250
        auto log = Dune::Dorie::create_base_logger(helper);
        log->set_level(spdlog::level::trace);
Lukas Riedel's avatar
Lukas Riedel committed
251 252 253 254

        // test the NearestNeighbor interpolator
        test_nearest_neighbor<2>();
        test_nearest_neighbor<3>();
Lukas Riedel's avatar
Lukas Riedel committed
255 256 257

        test_linear<2>();
        test_linear<3>();
Lukas Riedel's avatar
Lukas Riedel committed
258 259 260 261 262 263 264 265 266 267 268 269
    }

    catch (Dune::Exception &e) {
        std::cerr << "Dune reported error: " << e << std::endl;
        return 1;
    }
    catch (std::exception &e) {
        std::cerr << "Exception thrown: " << e.what() << std::endl;
        return 1;
    }
    catch (...) {
        std::cerr << "Unknown exception!" << std::endl;
270
        return 1;
Lukas Riedel's avatar
Lukas Riedel committed
271
    }
272
}