Line data Source code
1 : /*
2 : * Copyright (c) 2015: G-CSC, Goethe University Frankfurt
3 : * Author: Dmitry Logashenko
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : /*
34 : * Inverse Distance Weighting (IDW) interpolation for data sets
35 : */
36 : #ifndef __H__UG__LIB_DISC__SPATIAL_DISC__USER_DATA__IDW_USER_DATA__
37 : #define __H__UG__LIB_DISC__SPATIAL_DISC__USER_DATA__IDW_USER_DATA__
38 :
39 : #include <vector>
40 :
41 : // ug4 headers
42 : #include "common/common.h"
43 : #include "common/math/ugmath.h"
44 : #include "lib_disc/spatial_disc/user_data/std_glob_pos_data.h"
45 :
46 : namespace ug {
47 :
48 : /// Class for inverse distance weighting based on a general data type.
49 : /**
50 : * The static functions in the class compute the field given by the IDW
51 : * interpolation of a given order for given interpolation points and values.
52 : * Type of the values specified by the template parameter (and may be a scalar
53 : * or a tensor arithmetic type).
54 : *
55 : * Having a set of interpolation points \f$ \mathbf{x}_i \f$ with values
56 : * \f$ u_i \f$ at them, the interpolated value \f$ u \f$ at point
57 : * \f$ \mathbf{x} \not\in \{ \mathbf{x}_i \} \f$ is computed as
58 : * \f{eqnarray*}{
59 : * u = \frac {\sum_i w_i \cdot u_i} {\sum_i w_i},
60 : * \f}
61 : * where
62 : * \f{eqnarray*}{
63 : * w_i := \frac{1}{\| \mathbf{x} - \mathbf{x}_i \|_2^p}
64 : * \f}
65 : * (\f$ p \f$ being the order of the interpolation).
66 : * For \f$ \mathbf{x} = \mathbf{x}_i \f$ (up to some numerical precision), we
67 : * set \f$ u = u_i \f$.
68 : *
69 : * \remark For \f$ p \f$ less or equal the dimension of the geometric space,
70 : * the interpolated value \f$ u \f$ may be essentially influenced by the values
71 : * at \f$ \mathbf{x}_i \f$ located far away from \f$ \mathbf{x} \f$. However,
72 : * this influence may be avoided by specifying the radius \f$ R \f$: Then only
73 : * those \f$ \mathbf{x}_i \f$ (and therefore \f$ u_i \f$) are considered in the
74 : * sums, for which \f$ \| \mathbf{x} - \mathbf{x}_i \|_2 \le R \f$ holds.
75 : *
76 : * \remark Note that if the radius \f$ R \f$ is specified, then it can happen
77 : * that the interpolated value depends only on values at points located on one
78 : * side of \f$ \mathbf{x} \f$, completely ignoring a trend prescribed by points
79 : * on the other sides. Thus, \f$ R \f$ should be large enough.
80 : *
81 : * To loop the interpolation points, iterators of the templated type TPntIterator
82 : * are used. Every element of the reference elements should have two members:
83 : * pos and value:
84 : * TPntIterator ptr;
85 : * ptr->pos is a MathVector<WDim> object with the coordinates of the interpolation
86 : * point;
87 : * ptr->value is a TData object of the value at that point.
88 : *
89 : * \tparam WDim dimensionality of the space
90 : * \tparam TPntIterator interpolation point iterator type
91 : * \tparam TData type of the values to interpolate
92 : */
93 : template <int WDim, typename TPntIterator, typename TData = number>
94 : class IDWInterpolation
95 : {
96 : public:
97 :
98 : /// dimensionality of the space (i.e. of the coordinate vectors)
99 : static const int dim = WDim;
100 :
101 : /// type of the interpolation point iterator
102 : typedef TPntIterator t_pnt_iter;
103 :
104 : /// type of the data to extrapolate
105 : typedef TData data_type;
106 :
107 : public:
108 :
109 : /// computes the interpolation basing on all the interpolation points
110 : static void compute
111 : (
112 : data_type & res, ///< interpolated value
113 : const MathVector<dim> & pos, ///< geometric position where to interpolate
114 : t_pnt_iter pnt_beg, ///< the first interpolation point
115 : t_pnt_iter pnt_end, ///< delimiter of the iterpolation points
116 : number order, ///< order of the interpolation
117 : number small_dist = 1e-7 ///< distance at which we do not distinguish the points
118 : );
119 :
120 : /// computes the interpolation basing on the interpolation points in a given ball
121 : static void compute
122 : (
123 : data_type & res, ///< interpolated value
124 : const MathVector<dim> & pos, ///< geometric position where to interpolate
125 : number R, ///< radius of the ball (if 0 then the whole space)
126 : t_pnt_iter pnt_beg, ///< the first interpolation point
127 : t_pnt_iter pnt_end, ///< delimiter of the iterpolation points
128 : number order, ///< order of the interpolation
129 : number small_dist = 1e-7 ///< distance at which we do not distinguish the points
130 : );
131 : };
132 :
133 : /// UserData interface for the IDW interpolation
134 : /**
135 : * This class implements the UserData interface for the inverse-distance-weighting
136 : * interpolation.
137 : *
138 : * \sa IDWInterpolation
139 : *
140 : * Setting the radius to 0 means the unconstrained version of the IDW interpolation.
141 : *
142 : * \tparam WDim dimensionality of the geometric space (the world dimension)
143 : * \tparam TData type of the data to interpolate
144 : */
145 : template <int WDim, typename TData = number>
146 : class IDWUserData
147 : : public StdGlobPosData<IDWUserData<WDim, TData>, TData, WDim>
148 : {
149 : public:
150 :
151 : /// dimensionality of the space (i.e. of the coordinate vectors)
152 : static const int dim = WDim;
153 :
154 : /// type of the data to extrapolate
155 : typedef TData data_type;
156 :
157 : private:
158 :
159 : /// type of a interpolation point data item
160 : struct data_item
161 : {
162 : MathVector<dim> pos; ///< (global) geometrical coordinates of the point
163 : data_type value; ///< value at that point
164 :
165 0 : data_item (const MathVector<dim> & x, const data_type & v) : pos (x), value (v) {};
166 0 : data_item (const data_item & dat) : pos (dat.pos), value (dat.value) {};
167 : };
168 :
169 : public:
170 :
171 : /// class constructor that creates an empty object with default parameters
172 0 : IDWUserData ()
173 0 : : m_order (dim + 1), m_R (0)
174 0 : {}
175 :
176 : /// class constructor that creates an empty object with given parameters
177 0 : IDWUserData (number order, number R)
178 0 : : m_order (order), m_R (R)
179 0 : {}
180 :
181 : /// virtual destructor
182 0 : ~IDWUserData () {}
183 :
184 : public:
185 :
186 : /// sets the radius of the neighbourhood where the interpolation points are taken from
187 0 : void set_radius (number R) {m_R = R;}
188 :
189 : /// sets the order of the interpolation
190 0 : void set_order (number order) {m_order = order;}
191 :
192 : /// deletes all the interpolation points from the list
193 : void clear () {m_data.clear ();}
194 :
195 : /// loads data from a given stream (and appends the loaded points to the current list)
196 : void load_data_from (std::istream & in);
197 :
198 : /// loads data from a given file (and appends the loaded points to the current list)
199 : void load_data_from (const char * file_name);
200 :
201 : /// appends an interpolation point to the list
202 0 : void append (const MathVector<dim> & x, const data_type & val) {m_data.push_back (data_item (x, val));}
203 :
204 : public:
205 :
206 : /// evaluates the data at a given point
207 : inline void evaluate (data_type & value, const MathVector<dim> & x, number time, int si) const
208 : {
209 : typedef typename std::vector<data_item>::const_iterator pnt_iter_type;
210 0 : IDWInterpolation<dim, pnt_iter_type, data_type>::compute (value, x, m_R,
211 0 : m_data.begin (), m_data.end (), m_order);
212 : }
213 :
214 : private:
215 :
216 : std::vector<data_item> m_data; ///< interpolation points
217 : number m_order; ///< order of the interpolation
218 : number m_R; ///< radius of the neighbourhood to look for the interpolation points in (0 == infinite)
219 : };
220 :
221 : } // end namespace ug
222 :
223 : #include "invdist_user_data_impl.h"
224 :
225 : #endif // __H__UG__LIB_DISC__SPATIAL_DISC__USER_DATA__IDW_USER_DATA__
226 :
227 : /* End of File */
|
