tinterpolation.hh - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
HTML git clone git://src.adamsgaard.dk/pism
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tinterpolation.hh (5717B)
---
1 /* Copyright (C) 2015, 2017, 2018, 2019 PISM Authors
2 *
3 * This file is part of PISM.
4 *
5 * PISM is free software; you can redistribute it and/or modify it under the
6 * terms of the GNU General Public License as published by the Free Software
7 * Foundation; either version 3 of the License, or (at your option) any later
8 * version.
9 *
10 * PISM is distributed in the hope that it will be useful, but WITHOUT ANY
11 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
12 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13 * details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with PISM; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
18 */
19
20 #ifndef _INTERPOLATION_H_
21 #define _INTERPOLATION_H_
22
23 #include <vector>
24
25 namespace pism {
26
27 enum InterpolationType {LINEAR, NEAREST, PIECEWISE_CONSTANT, LINEAR_PERIODIC};
28
29 /**
30 * Class encapsulating linear interpolation indexes and weights.
31 *
32 * Most library interpolation routines (GSL's, for example) assume that we are working with a fixed
33 * function given on an input grid; these libraries support evaluating this function at arbitrary
34 * points. In this case an interpolation routine may use values on the input grid to create an
35 * interpolant, but has no knowledge of an "output grid."
36 *
37 * In PISM we have a slightly different situation: we need to transfer several gridded functions
38 * between two *fixed* grids, so we can use both grids to compute interpolation weights, but cannot
39 * use values on the input grid.
40 *
41 * If a point on the output grid is outside the interval defined by the input grid this code uses
42 * *constant extrapolation*.
43 *
44 * This class isolates the code computing interpolation weights so that we can use these 1D weights
45 * in a 2D or 3D regridding code. This avoids code duplication: in bilinear (trilinear)
46 * interpolation X and Y (and Z) grid dimensions are treated the same. This also makes it possible
47 * to test the code computing interpolation weights in isolation.
48 *
49 * We provide getters `left()`, `right()` and `alpha()`.
50 *
51 * Here `left[i]` is the index in the input grid (`x_input`) such that
52 * ~~~ c++
53 * input_x[left[i]] < output_x[i]
54 * ~~~
55 * similar for `right[i]`:
56 * ~~~ c++
57 * input_x[right[i]] >= output_x[i]
58 * ~~~
59 * When `output_x[i]` is outside the input interval `left[i]` and `right[i]` are adjusted so that
60 * the interpolation formula below corresponds to constant extrapolation.
61 *
62 * `alpha[i]` is the linear interpolation weight used like this:
63 * ~~~ c++
64 * const int
65 * L = m_left[k],
66 * R = m_right[k];
67 * const double Alpha = m_alpha[k];
68 * result[k] = input_values[L] + Alpha * (input_values[R] - input_values[L]);
69 * ~~~
70 *
71 * Piecewise constant 1D interpolation used for time-dependent forcing (fluxes that should
72 * be interpreted as piecewise-constant to simplify accounting of mass conservation).
73 *
74 * Here `input_x` defines left end-points of intervals, For example, [0, 1] defines two
75 * intervals: [0, 1) and [1, x_e). Here the value x_e is irrelevant because we use
76 * constant extrapolation for points outside the interval covered by input data.
77 *
78 *
79 * Linear interpolation for periodic data (annual temperature cycles, etc).
80 *
81 * Input data are assumed to be periodic on the interval from 0 to `period`.
82 *
83 */
84 class Interpolation {
85 public:
86 Interpolation(InterpolationType type, const std::vector<double> &input_x,
87 const std::vector<double> &output_x, double period = 0.0);
88 Interpolation(InterpolationType type, const double *input_x, unsigned int input_x_size,
89 const double *output_x, unsigned int output_x_size, double period = 0.0);
90
91 const std::vector<int>& left() const;
92 const std::vector<int>& right() const;
93 const std::vector<double>& alpha() const;
94
95 int left(size_t i) const;
96 int right(size_t i) const;
97 double alpha(size_t i) const;
98
99 //! Return interpolated values (on the output grid) given `input_values` on the input grid.
100 /** This is used for testing. (Regular code calls left(), right(), and alpha().)
101 */
102 std::vector<double> interpolate(const std::vector<double> &input_values) const;
103
104 /*!
105 * Compute interpolated values on the output grid given values on the input grid.
106 */
107 void interpolate(const double *input, double *output) const;
108
109 double integral(const double *input) const;
110 double integral(const std::vector<double> &input) const;
111 double interval_length() const;
112 private:
113 //! Interpolation indexes
114 std::vector<int> m_left, m_right;
115 //! Interpolation weights
116 std::vector<double> m_alpha;
117
118 //! Integration weights
119 std::vector<double> m_w;
120 //! Length of the interval defined by `output_x`.
121 double m_interval_length;
122
123 void init_linear(const double *input_x, unsigned int input_x_size,
124 const double *output_x, unsigned int output_x_size);
125 void init_nearest(const double *input_x, unsigned int input_x_size,
126 const double *output_x, unsigned int output_x_size);
127 void init_piecewise_constant(const double *input_x, unsigned int input_x_size,
128 const double *output_x, unsigned int output_x_size);
129 void init_linear_periodic(const double *input_x, unsigned int input_x_size,
130 const double *output_x, unsigned int output_x_size,
131 double period);
132
133 void init_weights_linear(const double *x,
134 unsigned int x_size,
135 const double *output_x,
136 unsigned int output_x_size);
137 };
138
139 } // end of namespace pism
140
141 #endif /* _INTERPOLATION_H_ */