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 : * Implementation of the Inverse Distance Weighting (IDW) interpolation for data sets
35 : */
36 :
37 : namespace ug {
38 :
39 : /**
40 : * Computes the interpolation basing on all the interpolation point.
41 : */
42 : template <int WDim, typename TPntIterator, typename TData>
43 0 : void IDWInterpolation<WDim, TPntIterator, TData>::compute
44 : (
45 : data_type & res, ///< interpolated value
46 : const MathVector<dim> & pos, ///< geometric position where to interpolate
47 : t_pnt_iter pnt_beg, ///< the first interpolation point
48 : t_pnt_iter pnt_end, ///< delimiter of the iterpolation points
49 : number order, ///< order of the interpolation
50 : number small_dist ///< distance at which we do not distinguish the points
51 : )
52 : {
53 0 : if (pnt_beg == pnt_end)
54 0 : UG_THROW ("IDWInterpolation: Cannot interpolate using 0 data points!");
55 :
56 : data_type sum = 0;
57 : number factor = 0;
58 0 : for (t_pnt_iter pnt = pnt_beg; pnt != pnt_end; ++pnt)
59 : {
60 : number dist;
61 0 : if ((dist = VecDistance (pos, pnt->pos)) < small_dist)
62 : { /* We are at a data point: */
63 0 : res = pnt->value;
64 : return;
65 : }
66 0 : dist = pow (dist, order);
67 0 : data_type value = pnt->value;
68 0 : value /= dist;
69 0 : sum += value;
70 0 : factor += 1 / dist;
71 : }
72 0 : sum /= factor;
73 0 : res = sum;
74 : }
75 :
76 : /**
77 : * Computes the interpolation basing on the interpolation points in a given ball.
78 : */
79 : template <int WDim, typename TPntIterator, typename TData>
80 0 : void IDWInterpolation<WDim, TPntIterator, TData>::compute
81 : (
82 : data_type & res, ///< interpolated value
83 : const MathVector<dim> & pos, ///< geometric position where to interpolate
84 : number R, ///< radius of the ball (if 0 then the whole space)
85 : t_pnt_iter pnt_beg, ///< the first interpolation point
86 : t_pnt_iter pnt_end, ///< delimiter of the iterpolation points
87 : number order, ///< order of the interpolation
88 : number small_dist ///< distance at which we do not distinguish the points
89 : )
90 : {
91 0 : if (R == 0.0) // a special case: we do not consider the ball
92 : {
93 : IDWInterpolation<dim, t_pnt_iter, data_type>::compute
94 0 : (res, pos, pnt_beg, pnt_end, order, small_dist);
95 0 : return;
96 : }
97 :
98 0 : if (pnt_beg == pnt_end)
99 0 : UG_THROW ("IDWInterpolation: Cannot interpolate using 0 data points!");
100 :
101 : data_type sum = 0;
102 : number factor = 0;
103 : bool no_data_in_ball = true;
104 0 : for (t_pnt_iter pnt = pnt_beg; pnt != pnt_end; ++pnt)
105 : {
106 : number dist;
107 0 : if ((dist = VecDistance (pos, pnt->pos)) < small_dist)
108 : { /* We are at a data point: */
109 0 : res = pnt->value;
110 : return;
111 : }
112 0 : if (dist > R) continue;
113 : no_data_in_ball = false;
114 0 : dist = pow (dist, order);
115 0 : data_type value = pnt->value;
116 0 : value /= dist;
117 0 : sum += value;
118 0 : factor += 1 / dist;
119 : }
120 0 : if (no_data_in_ball)
121 0 : UG_THROW ("IDWInterpolation: No interpolation points in the ball with R = " << R
122 : << " and center at" << pos << ".");
123 0 : sum /= factor;
124 0 : res = sum;
125 : }
126 :
127 : /**
128 : * Loads interpolation points from a given stream.
129 : */
130 : template <int WDim, typename TData>
131 0 : void IDWUserData<WDim, TData>::load_data_from (std::istream & input)
132 : {
133 : std::string input_line;
134 0 : std::istringstream line_stream;
135 : MathVector<dim, number> pos;
136 : data_type value;
137 :
138 : line_stream.exceptions (std::istream::failbit | std::istream::badbit);
139 :
140 0 : while (! input.eof ())
141 : {
142 0 : if (input.fail ())
143 0 : UG_THROW ("IDWUserData: Could not load the interpolation points from the file!");
144 0 : std::getline (input, input_line);
145 0 : if (input_line.length () == 0)
146 0 : continue;
147 : try
148 : {
149 : line_stream.str (input_line);
150 0 : read_plain_txt (line_stream, pos);
151 : line_stream >> value;
152 : }
153 0 : catch (std::istream::failure& e)
154 : {
155 0 : UG_THROW ("IDWUserData: Failed to parse line '" << input_line << "' in the input file.");
156 : };
157 :
158 : this->append (pos, value);
159 : }
160 0 : }
161 :
162 : /**
163 : * Loads interpolation points from a given file.
164 : */
165 : template <int WDim, typename TData>
166 0 : void IDWUserData<WDim, TData>::load_data_from (const char * file_name)
167 : {
168 0 : std::ifstream input (file_name, std::ifstream::in);
169 0 : if (input.fail ())
170 0 : UG_THROW ("IDWUserData: Cannot open data file '" << file_name << "' for input!");
171 0 : this->load_data_from (input);
172 0 : }
173 :
174 : } // end namespace ug
175 :
176 : /* End of File */
|
