PISM, A Parallel Ice Sheet Model  stable v2.0.4 committed by Constantine Khrulev on 2022-05-25 12:02:27 -0800
interpolation.hh
Go to the documentation of this file.
1 /* Copyright (C) 2015, 2017, 2018, 2019, 2021 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 #include <map>
25 
26 namespace pism {
27 
29 
30 /**
31  * Class encapsulating linear and piece-wise constant interpolation indexes and weights.
32  *
33  * Most library interpolation routines (GSL's, for example) assume that we are working with a fixed
34  * function given on an input grid; these libraries support evaluating this function at arbitrary
35  * points. In this case an interpolation routine may use values on the input grid to create an
36  * interpolant, but has no knowledge of an "output grid."
37  *
38  * In PISM we have a slightly different situation: we need to transfer several gridded functions
39  * between two *fixed* grids, so we can use both grids to compute interpolation weights, but cannot
40  * use values on the input grid.
41  *
42  * If a point on the output grid is outside the interval defined by the input grid this code uses
43  * *constant extrapolation*.
44  *
45  * This class isolates the code computing interpolation weights so that we can use these 1D weights
46  * in a 2D or 3D regridding code. This avoids code duplication: in bilinear (trilinear)
47  * interpolation X and Y (and Z) grid dimensions are treated the same. This also makes it possible
48  * to test the code computing interpolation weights in isolation.
49  *
50  * We provide getters `left()`, `right()` and `alpha()`.
51  *
52  * Here `left[i]` is the index in the input grid (`x_input`) such that
53  * ~~~ c++
54  * input_x[left[i]] < output_x[i]
55  * ~~~
56  * similar for `right[i]`:
57  * ~~~ c++
58  * input_x[right[i]] >= output_x[i]
59  * ~~~
60  * When `output_x[i]` is outside the input interval `left[i]` and `right[i]` are adjusted so that
61  * the interpolation formula below corresponds to constant extrapolation.
62  *
63  * `alpha[i]` is the linear interpolation weight used like this:
64  * ~~~ c++
65  * const int
66  * L = m_left[k],
67  * R = m_right[k];
68  * const double Alpha = m_alpha[k];
69  * result[k] = input_values[L] + Alpha * (input_values[R] - input_values[L]);
70  * ~~~
71  *
72  * Piecewise constant 1D interpolation used for time-dependent forcing (fluxes that should
73  * be interpreted as piecewise-constant to simplify accounting of mass conservation).
74  *
75  * Here `input_x` defines left end-points of intervals, For example, [0, 1] defines two
76  * intervals: [0, 1) and [1, x_e). Here the value x_e is irrelevant because we use
77  * constant extrapolation for points outside the interval covered by input data.
78  *
79  */
81 public:
82  Interpolation(InterpolationType type, const std::vector<double> &input_x,
83  const std::vector<double> &output_x);
84  Interpolation(InterpolationType type, const double *input_x, unsigned int input_x_size,
85  const double *output_x, unsigned int output_x_size);
86 
87  const std::vector<int>& left() const;
88  const std::vector<int>& right() const;
89  const std::vector<double>& alpha() const;
90 
91  int left(size_t j) const;
92  int right(size_t j) const;
93  double alpha(size_t j) const;
94 
95  //! Return interpolated values (on the output grid) given `input_values` on the input grid.
96  /** This is used for testing. (Regular code calls left(), right(), and alpha().)
97  */
98  std::vector<double> interpolate(const std::vector<double> &input_values) const;
99 
100  /*!
101  * Compute interpolated values on the output grid given values on the input grid.
102  */
103  void interpolate(const double *input, double *output) const;
104 
105 private:
106  //! Interpolation indexes
107  std::vector<int> m_left, m_right;
108  //! Interpolation weights
109  std::vector<double> m_alpha;
110 
111  void init_linear(const double *input_x, unsigned int input_x_size,
112  const double *output_x, unsigned int output_x_size);
113  void init_nearest(const double *input_x, unsigned int input_x_size,
114  const double *output_x, unsigned int output_x_size);
115  void init_piecewise_constant(const double *input_x, unsigned int input_x_size,
116  const double *output_x, unsigned int output_x_size);
117 };
118 
119 std::map<size_t, double> integration_weights(const double *x,
120  size_t x_size,
121  InterpolationType type,
122  double a,
123  double b);
124 
125 std::map<size_t, double> integration_weights(const std::vector<double> &x,
126  InterpolationType type,
127  double a,
128  double b);
129 
130 } // end of namespace pism
131 
132 #endif /* _INTERPOLATION_H_ */
void init_linear(const double *input_x, unsigned int input_x_size, const double *output_x, unsigned int output_x_size)
std::vector< double > m_alpha
Interpolation weights.
Interpolation(InterpolationType type, const std::vector< double > &input_x, const std::vector< double > &output_x)
void init_piecewise_constant(const double *input_x, unsigned int input_x_size, const double *output_x, unsigned int output_x_size)
std::vector< int > m_left
Interpolation indexes.
void init_nearest(const double *input_x, unsigned int input_x_size, const double *output_x, unsigned int output_x_size)
const std::vector< double > & alpha() const
const std::vector< int > & right() const
const std::vector< int > & left() const
std::vector< double > interpolate(const std::vector< double > &input_values) const
Return interpolated values (on the output grid) given input_values on the input grid.
std::vector< int > m_right
std::map< size_t, double > integration_weights(const double *x, size_t x_size, InterpolationType type, double a, double b)
InterpolationType
@ PIECEWISE_CONSTANT