1 | /*
|
---|
2 | * Project: MoleCuilder
|
---|
3 | * Description: creates and alters molecular systems
|
---|
4 | * Copyright (C) 2012 University of Bonn. All rights reserved.
|
---|
5 | * Copyright (C) 2013 Frederik Heber. All rights reserved.
|
---|
6 | *
|
---|
7 | *
|
---|
8 | * This file is part of MoleCuilder.
|
---|
9 | *
|
---|
10 | * MoleCuilder is free software: you can redistribute it and/or modify
|
---|
11 | * it under the terms of the GNU General Public License as published by
|
---|
12 | * the Free Software Foundation, either version 2 of the License, or
|
---|
13 | * (at your option) any later version.
|
---|
14 | *
|
---|
15 | * MoleCuilder is distributed in the hope that it will be useful,
|
---|
16 | * but WITHOUT ANY WARRANTY; without even the implied warranty of
|
---|
17 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
---|
18 | * GNU General Public License for more details.
|
---|
19 | *
|
---|
20 | * You should have received a copy of the GNU General Public License
|
---|
21 | * along with MoleCuilder. If not, see <http://www.gnu.org/licenses/>.
|
---|
22 | */
|
---|
23 |
|
---|
24 | /*
|
---|
25 | * SamplingGrid.cpp
|
---|
26 | *
|
---|
27 | * Created on: 25.07.2012
|
---|
28 | * Author: heber
|
---|
29 | */
|
---|
30 |
|
---|
31 | // include config.h
|
---|
32 | #ifdef HAVE_CONFIG_H
|
---|
33 | #include <config.h>
|
---|
34 | #endif
|
---|
35 |
|
---|
36 | // include headers that implement a archive in simple text format
|
---|
37 | // otherwise BOOST_CLASS_EXPORT_IMPLEMENT has no effect
|
---|
38 | #include <boost/archive/text_oarchive.hpp>
|
---|
39 | #include <boost/archive/text_iarchive.hpp>
|
---|
40 |
|
---|
41 | //#include "CodePatterns/MemDebug.hpp"
|
---|
42 |
|
---|
43 | #include "Fragmentation/Summation/SetValues/SamplingGrid.hpp"
|
---|
44 |
|
---|
45 | #include <boost/assign.hpp>
|
---|
46 | #include <boost/bind.hpp>
|
---|
47 | #include <boost/shared_ptr.hpp>
|
---|
48 |
|
---|
49 | #include <algorithm>
|
---|
50 | #include <limits>
|
---|
51 |
|
---|
52 | #include "CodePatterns/Assert.hpp"
|
---|
53 | #include "CodePatterns/Log.hpp"
|
---|
54 |
|
---|
55 | #include "LinearAlgebra/Vector.hpp"
|
---|
56 |
|
---|
57 | // static instances
|
---|
58 | const double SamplingGrid::zeroOffset[NDIM] = { 0., 0., 0. };
|
---|
59 |
|
---|
60 | SamplingGrid::SamplingGrid() :
|
---|
61 | SamplingGridProperties()
|
---|
62 | {
|
---|
63 | setWindowSize(zeroOffset, zeroOffset);
|
---|
64 | ASSERT( getWindowGridPoints() == (size_t)0,
|
---|
65 | "SamplingGrid::SamplingGrid() - incorrect number of samples given for the window.");
|
---|
66 | }
|
---|
67 |
|
---|
68 | SamplingGrid::SamplingGrid(const double _begin[NDIM],
|
---|
69 | const double _end[NDIM],
|
---|
70 | const int _level) :
|
---|
71 | SamplingGridProperties(_begin, _end, _level)
|
---|
72 | {
|
---|
73 | setWindowSize(zeroOffset, zeroOffset);
|
---|
74 | ASSERT( getWindowGridPoints() == (size_t)0,
|
---|
75 | "SamplingGrid::SamplingGrid() - incorrect number of samples given for the window.");
|
---|
76 | }
|
---|
77 |
|
---|
78 | SamplingGrid::SamplingGrid(const double _begin[NDIM],
|
---|
79 | const double _end[NDIM],
|
---|
80 | const int _level,
|
---|
81 | const sampledvalues_t &_sampled_grid) :
|
---|
82 | SamplingGridProperties(_begin, _end, _level),
|
---|
83 | sampled_grid(_sampled_grid)
|
---|
84 | {
|
---|
85 | setWindowSize(_begin, _end);
|
---|
86 | ASSERT( getWindowGridPoints() == (size_t)_sampled_grid.size(),
|
---|
87 | "SamplingGrid::SamplingGrid() - incorrect number of samples given for the window.");
|
---|
88 | }
|
---|
89 |
|
---|
90 | SamplingGrid::SamplingGrid(const SamplingGrid &_grid) :
|
---|
91 | SamplingGridProperties(_grid),
|
---|
92 | sampled_grid(_grid.sampled_grid)
|
---|
93 | {
|
---|
94 | setWindowSize(_grid.begin_window, _grid.end_window);
|
---|
95 | ASSERT( getWindowGridPoints() == _grid.getWindowGridPoints(),
|
---|
96 | "SamplingGrid::SamplingGrid() - incorrect number of samples given for the window.");
|
---|
97 | }
|
---|
98 |
|
---|
99 | SamplingGrid::SamplingGrid(const SamplingGridProperties &_props) :
|
---|
100 | SamplingGridProperties(_props)
|
---|
101 | {
|
---|
102 | setWindowSize(zeroOffset, zeroOffset);
|
---|
103 | ASSERT( getWindowGridPoints() == (size_t)0,
|
---|
104 | "SamplingGrid::SamplingGrid() - incorrect number of samples given for the window.");
|
---|
105 | }
|
---|
106 |
|
---|
107 | SamplingGrid::SamplingGrid(
|
---|
108 | const SamplingGridProperties &_props,
|
---|
109 | const sampledvalues_t &_sampled_grid) :
|
---|
110 | SamplingGridProperties(_props),
|
---|
111 | sampled_grid(_sampled_grid)
|
---|
112 | {
|
---|
113 | setWindowSize(_props.begin, _props.end);
|
---|
114 | ASSERT( getWindowGridPoints() == (size_t)_sampled_grid.size(),
|
---|
115 | "SamplingGrid::SamplingGrid() - incorrect number of samples given for the window.");
|
---|
116 | }
|
---|
117 |
|
---|
118 | SamplingGrid::~SamplingGrid()
|
---|
119 | {}
|
---|
120 |
|
---|
121 | bool SamplingGrid::isCongruent(const SamplingGrid &_props) const
|
---|
122 | {
|
---|
123 | bool status = true;
|
---|
124 | status &= (static_cast<const SamplingGridProperties &>(*this) ==
|
---|
125 | static_cast<const SamplingGridProperties &>(_props));
|
---|
126 | for(size_t i = 0; i<NDIM; ++i) {
|
---|
127 | status &= begin_window[i] == _props.begin_window[i];
|
---|
128 | status &= end_window[i] == _props.end_window[i];
|
---|
129 | }
|
---|
130 | return status;
|
---|
131 | }
|
---|
132 |
|
---|
133 | SamplingGrid& SamplingGrid::operator=(const SamplingGrid& other)
|
---|
134 | {
|
---|
135 | // check for self-assignment
|
---|
136 | if (this != &other) {
|
---|
137 | static_cast<SamplingGridProperties &>(*this) = other;
|
---|
138 | setWindowSize(other.begin_window, other.end_window);
|
---|
139 | sampled_grid = other.sampled_grid;
|
---|
140 | }
|
---|
141 | return *this;
|
---|
142 | }
|
---|
143 |
|
---|
144 | static void multiplyElements(
|
---|
145 | double &dest,
|
---|
146 | const double &source,
|
---|
147 | const double prefactor)
|
---|
148 | {
|
---|
149 | dest *= prefactor*(source);
|
---|
150 | }
|
---|
151 |
|
---|
152 | SamplingGrid& SamplingGrid::operator*=(const double _value)
|
---|
153 | {
|
---|
154 | std::transform(
|
---|
155 | sampled_grid.begin(), sampled_grid.end(),
|
---|
156 | sampled_grid.begin(),
|
---|
157 | boost::bind(std::multiplies<double>(), _1, _value));
|
---|
158 |
|
---|
159 | return *this;
|
---|
160 | }
|
---|
161 |
|
---|
162 | static void getDownsampledGrid(
|
---|
163 | const SamplingGrid &_reference_grid,
|
---|
164 | const SamplingGrid &_grid,
|
---|
165 | boost::shared_ptr<SamplingGrid> &_weight_downsampled)
|
---|
166 | {
|
---|
167 | static const double round_offset(
|
---|
168 | (std::numeric_limits<size_t>::round_style == std::round_toward_zero) ?
|
---|
169 | 0.5 : 0.); // need offset to get to round_toward_nearest behavior
|
---|
170 | const int surplus_level = _reference_grid.getSurplusLevel(_grid)+round_offset;
|
---|
171 | // need to downsample first
|
---|
172 | _weight_downsampled.reset(new SamplingGrid); // may use default cstor as downsample does all settings
|
---|
173 | SamplingGrid::downsample(*_weight_downsampled, _grid, _reference_grid.level-surplus_level);
|
---|
174 | }
|
---|
175 |
|
---|
176 | SamplingGrid& SamplingGrid::operator*=(const SamplingGrid& other)
|
---|
177 | {
|
---|
178 | // check that grids are compatible
|
---|
179 | ASSERT(isCompatible(other),
|
---|
180 | "SamplingGrid::operator*=() - multiplying incomatible grids is so far not in the cards.");
|
---|
181 | const SamplingGrid *other_grid = &other;
|
---|
182 | boost::shared_ptr<SamplingGrid> other_downsampled;
|
---|
183 | if (!isEquivalent(other)) {
|
---|
184 | getDownsampledGrid(*this, other, other_downsampled);
|
---|
185 | other_grid = other_downsampled.get();
|
---|
186 | }
|
---|
187 | /// get minimum of window
|
---|
188 | double min_begin_window[NDIM];
|
---|
189 | double min_end_window[NDIM];
|
---|
190 | bool doShrink = false;
|
---|
191 | for (size_t index=0; index<NDIM;++index) {
|
---|
192 | if (begin_window[index] <= other.begin_window[index]) {
|
---|
193 | min_begin_window[index] = other.begin_window[index];
|
---|
194 | doShrink = true;
|
---|
195 | } else {
|
---|
196 | min_begin_window[index] = begin_window[index];
|
---|
197 | }
|
---|
198 | if (end_window[index] <= other.end_window[index]) {
|
---|
199 | min_end_window[index] = end_window[index];
|
---|
200 | } else {
|
---|
201 | min_end_window[index] = other.end_window[index];
|
---|
202 | doShrink = true;
|
---|
203 | }
|
---|
204 | }
|
---|
205 | LOG(4, "DEBUG: min begin is " << min_begin_window[0] << "," << min_begin_window[1] << "," << min_begin_window[2] << ".");
|
---|
206 | LOG(4, "DEBUG: min end is " << min_end_window[0] << "," << min_end_window[1] << "," << min_end_window[2] << ".");
|
---|
207 | if (doShrink)
|
---|
208 | shrinkWindow(min_begin_window, min_end_window);
|
---|
209 | addWindowOntoWindow(
|
---|
210 | other_grid->begin_window,
|
---|
211 | other_grid->end_window,
|
---|
212 | begin_window,
|
---|
213 | end_window,
|
---|
214 | sampled_grid,
|
---|
215 | other_grid->sampled_grid,
|
---|
216 | boost::bind(multiplyElements, _1, _2, 1.),
|
---|
217 | sourcewindow);
|
---|
218 |
|
---|
219 | return *this;
|
---|
220 | }
|
---|
221 |
|
---|
222 | void SamplingGrid::superposeOtherGrids(const SamplingGrid &other, const double prefactor)
|
---|
223 | {
|
---|
224 | // check that grids are compatible
|
---|
225 | ASSERT(isCompatible(other),
|
---|
226 | "SamplingGrid::superposeOtherGrids() - superposing incompatible grids is so far not in the cards.");
|
---|
227 | const SamplingGrid *other_grid = &other;
|
---|
228 | boost::shared_ptr<SamplingGrid> other_downsampled;
|
---|
229 | if (!isEquivalent(other)) {
|
---|
230 | getDownsampledGrid(*this, other, other_downsampled);
|
---|
231 | other_grid = other_downsampled.get();
|
---|
232 | }
|
---|
233 | /// get maximum of window
|
---|
234 | double max_begin_window[NDIM];
|
---|
235 | double max_end_window[NDIM];
|
---|
236 | bool doExtend = false;
|
---|
237 | for (size_t index=0; index<NDIM;++index) {
|
---|
238 | if (begin_window[index] >= other.begin_window[index]) {
|
---|
239 | max_begin_window[index] = other.begin_window[index];
|
---|
240 | doExtend = true;
|
---|
241 | } else {
|
---|
242 | max_begin_window[index] = begin_window[index];
|
---|
243 | }
|
---|
244 | if (end_window[index] >= other.end_window[index]) {
|
---|
245 | max_end_window[index] = end_window[index];
|
---|
246 | } else {
|
---|
247 | max_end_window[index] = other.end_window[index];
|
---|
248 | doExtend = true;
|
---|
249 | }
|
---|
250 | }
|
---|
251 | LOG(4, "DEBUG: max begin is " << max_begin_window[0] << "," << max_begin_window[1] << "," << max_begin_window[2] << ".");
|
---|
252 | LOG(4, "DEBUG: max end is " << max_end_window[0] << "," << max_end_window[1] << "," << max_end_window[2] << ".");
|
---|
253 | if (doExtend)
|
---|
254 | extendWindow(max_begin_window, max_end_window);
|
---|
255 | /// and copy other into larger window, too
|
---|
256 | addOntoWindow(other_grid->begin_window, other_grid->end_window, other_grid->sampled_grid, prefactor);
|
---|
257 | }
|
---|
258 |
|
---|
259 | const size_t SamplingGrid::getWindowGridPointsPerAxis(const size_t axis) const
|
---|
260 | {
|
---|
261 | static const double round_offset(
|
---|
262 | (std::numeric_limits<size_t>::round_style == std::round_toward_zero) ?
|
---|
263 | 0.5 : 0.); // need offset to get to round_toward_nearest behavior
|
---|
264 | // const double total = getTotalLengthPerAxis(axis);
|
---|
265 | const double delta = getDeltaPerAxis(axis);
|
---|
266 | if (delta == 0)
|
---|
267 | return 0;
|
---|
268 | const double length = getWindowLengthPerAxis(axis);
|
---|
269 | if (length == 0)
|
---|
270 | return 0;
|
---|
271 | return (size_t)(length/delta+round_offset);
|
---|
272 | }
|
---|
273 |
|
---|
274 | double SamplingGrid::integral() const
|
---|
275 | {
|
---|
276 | const double volume_element = getVolume()/(double)getTotalGridPoints();
|
---|
277 | double int_value = 0.;
|
---|
278 | for (sampledvalues_t::const_iterator iter = sampled_grid.begin();
|
---|
279 | iter != sampled_grid.end();
|
---|
280 | ++iter)
|
---|
281 | int_value += *iter;
|
---|
282 | int_value *= volume_element;
|
---|
283 | LOG(2, "DEBUG: SamplingGrid::integral() is " << scientific << setprecision(13) << int_value << ".");
|
---|
284 | return int_value;
|
---|
285 | }
|
---|
286 |
|
---|
287 | double SamplingGrid::integral(const SamplingGrid &weight) const
|
---|
288 | {
|
---|
289 | // check that grids are compatible
|
---|
290 | ASSERT(isCompatible(weight),
|
---|
291 | "SamplingGrid::integral() - integrating with weights from incompatible grids is so far not in the cards.");
|
---|
292 | const SamplingGrid *weight_grid = &weight;
|
---|
293 | boost::shared_ptr<SamplingGrid> weight_downsampled;
|
---|
294 | if (!isEquivalent(weight)) {
|
---|
295 | getDownsampledGrid(*this, weight, weight_downsampled);
|
---|
296 | weight_grid = weight_downsampled.get();
|
---|
297 | }
|
---|
298 | const double volume_element = getVolume()/(double)getTotalGridPoints();
|
---|
299 | double int_value = 0.;
|
---|
300 | sampledvalues_t::const_iterator iter = sampled_grid.begin();
|
---|
301 | sampledvalues_t::const_iterator weightiter = weight_grid->sampled_grid.begin();
|
---|
302 | for (;iter != sampled_grid.end();++iter,++weightiter)
|
---|
303 | int_value += (*weightiter) * (*iter);
|
---|
304 | int_value *= volume_element;
|
---|
305 | //LOG(2, "DEBUG: SamplingGrid::integral() is " << scientific << setprecision(13) << int_value << ".");
|
---|
306 | return int_value;
|
---|
307 | }
|
---|
308 |
|
---|
309 | void SamplingGrid::setWindowSize(
|
---|
310 | const double _begin_window[NDIM],
|
---|
311 | const double _end_window[NDIM])
|
---|
312 | {
|
---|
313 | for (size_t index=0;index<NDIM;++index) {
|
---|
314 | begin_window[index] = getNearestLowerGridPoint(_begin_window[index], index);
|
---|
315 | ASSERT( begin_window[index] >= begin[index],
|
---|
316 | "SamplingGrid::setWindowSize() - window starts earlier than domain on "
|
---|
317 | +toString(index)+"th component.");
|
---|
318 | end_window[index] = getNearestHigherGridPoint(_end_window[index], index);
|
---|
319 | ASSERT( end_window[index] <= end[index],
|
---|
320 | "SamplingGrid::setWindowSize() - window ends later than domain on "
|
---|
321 | +toString(index)+"th component.");
|
---|
322 | }
|
---|
323 | }
|
---|
324 |
|
---|
325 | void SamplingGrid::setWindow(
|
---|
326 | const double _begin_window[NDIM],
|
---|
327 | const double _end_window[NDIM])
|
---|
328 | {
|
---|
329 | setWindowSize(_begin_window, _end_window);
|
---|
330 | const size_t gridpoints_window = getWindowGridPoints();
|
---|
331 | sampled_grid.clear();
|
---|
332 | sampled_grid.resize(gridpoints_window, 0.);
|
---|
333 | }
|
---|
334 |
|
---|
335 | void SamplingGrid::setDomain(
|
---|
336 | const double _begin[NDIM],
|
---|
337 | const double _end[NDIM])
|
---|
338 | {
|
---|
339 | setDomainSize(_begin, _end);
|
---|
340 | setWindowSize(_begin, _end);
|
---|
341 | const size_t gridpoints = getTotalGridPoints();
|
---|
342 | sampled_grid.resize(gridpoints, 0.);
|
---|
343 | }
|
---|
344 |
|
---|
345 | void SamplingGrid::extendWindow(
|
---|
346 | const double _begin_window[NDIM],
|
---|
347 | const double _end_window[NDIM])
|
---|
348 | {
|
---|
349 | #ifndef NDEBUG
|
---|
350 | for(size_t index=0;index < NDIM; ++index) {
|
---|
351 | // check that we truly have to extend the window
|
---|
352 | ASSERT ( begin_window[index] >= _begin_window[index],
|
---|
353 | "SamplingGrid::extendWindow() - component "+toString(index)+
|
---|
354 | " of window start is greater than old value.");
|
---|
355 | ASSERT ( end_window[index] <= _end_window[index],
|
---|
356 | "SamplingGrid::extendWindow() - component "+toString(index)+
|
---|
357 | " of window end is less than old value.");
|
---|
358 |
|
---|
359 | // check that we are still less than domain
|
---|
360 | ASSERT ( _begin_window[index] >= begin[index],
|
---|
361 | "SamplingGrid::extendWindow() - component "+toString(index)+
|
---|
362 | " of window start is less than domain start.");
|
---|
363 | ASSERT ( _end_window[index] <= end[index],
|
---|
364 | "SamplingGrid::extendWindow() - component "+toString(index)+
|
---|
365 | " of window end is greater than domain end.");
|
---|
366 | }
|
---|
367 | #endif
|
---|
368 | // copy old window size and values
|
---|
369 | double old_begin_window[NDIM];
|
---|
370 | double old_end_window[NDIM];
|
---|
371 | for(size_t index=0;index<NDIM;++index) {
|
---|
372 | old_begin_window[index] = begin_window[index];
|
---|
373 | old_end_window[index] = end_window[index];
|
---|
374 | }
|
---|
375 | sampledvalues_t old_values(sampled_grid);
|
---|
376 | // set new window
|
---|
377 | setWindow(_begin_window,_end_window);
|
---|
378 | // now extend it ...
|
---|
379 | addOntoWindow(old_begin_window, old_end_window, old_values, +1.);
|
---|
380 | LOG(6, "DEBUG: Grid after extension is " << sampled_grid << ".");
|
---|
381 | }
|
---|
382 |
|
---|
383 | void SamplingGrid::shrinkWindow(
|
---|
384 | const double _begin_window[NDIM],
|
---|
385 | const double _end_window[NDIM])
|
---|
386 | {
|
---|
387 | #ifndef NDEBUG
|
---|
388 | for(size_t index=0;index < NDIM; ++index) {
|
---|
389 | // check that we truly have to shrink the window
|
---|
390 | ASSERT ( begin_window[index] <= _begin_window[index],
|
---|
391 | "SamplingGrid::shrinkWindow() - component "+toString(index)+
|
---|
392 | " of window start is less than old value.");
|
---|
393 | ASSERT ( end_window[index] >= _end_window[index],
|
---|
394 | "SamplingGrid::shrinkWindow() - component "+toString(index)+
|
---|
395 | " of window end is greater than old value.");
|
---|
396 |
|
---|
397 | // check that we are still less than domain
|
---|
398 | ASSERT ( _begin_window[index] >= begin[index],
|
---|
399 | "SamplingGrid::shrinkWindow() - component "+toString(index)+
|
---|
400 | " of window start is less than domain start.");
|
---|
401 | ASSERT ( _end_window[index] <= end[index],
|
---|
402 | "SamplingGrid::shrinkWindow() - component "+toString(index)+
|
---|
403 | " of window end is greater than domain end.");
|
---|
404 | }
|
---|
405 | #endif
|
---|
406 | // copy old window size and values
|
---|
407 | double old_begin_window[NDIM];
|
---|
408 | double old_end_window[NDIM];
|
---|
409 | for(size_t index=0;index<NDIM;++index) {
|
---|
410 | old_begin_window[index] = begin_window[index];
|
---|
411 | old_end_window[index] = end_window[index];
|
---|
412 | }
|
---|
413 | sampledvalues_t old_values(sampled_grid);
|
---|
414 | // set new window
|
---|
415 | setWindow(_begin_window,_end_window);
|
---|
416 | // now extend it ...
|
---|
417 | addIntoWindow(old_begin_window, old_end_window, old_values, +1.);
|
---|
418 | LOG(6, "DEBUG: Grid after extension is " << sampled_grid << ".");
|
---|
419 | }
|
---|
420 |
|
---|
421 | static void addElements(
|
---|
422 | double &dest,
|
---|
423 | const double &source,
|
---|
424 | const double prefactor)
|
---|
425 | {
|
---|
426 | dest += prefactor*(source);
|
---|
427 | }
|
---|
428 |
|
---|
429 | void SamplingGrid::addOntoWindow(
|
---|
430 | const double _begin_window[NDIM],
|
---|
431 | const double _end_window[NDIM],
|
---|
432 | const sampledvalues_t &_sampled_grid,
|
---|
433 | const double prefactor)
|
---|
434 | {
|
---|
435 | addWindowOntoWindow(
|
---|
436 | begin_window,
|
---|
437 | end_window,
|
---|
438 | _begin_window,
|
---|
439 | _end_window,
|
---|
440 | sampled_grid,
|
---|
441 | _sampled_grid,
|
---|
442 | boost::bind(addElements, _1, _2, boost::cref(prefactor)),
|
---|
443 | destwindow);
|
---|
444 | }
|
---|
445 |
|
---|
446 | void SamplingGrid::addIntoWindow(
|
---|
447 | const double _begin_window[NDIM],
|
---|
448 | const double _end_window[NDIM],
|
---|
449 | const sampledvalues_t &_sampled_grid,
|
---|
450 | const double prefactor)
|
---|
451 | {
|
---|
452 | addWindowOntoWindow(
|
---|
453 | _begin_window,
|
---|
454 | _end_window,
|
---|
455 | begin_window,
|
---|
456 | end_window,
|
---|
457 | sampled_grid,
|
---|
458 | _sampled_grid,
|
---|
459 | boost::bind(addElements, _1, _2, boost::cref(prefactor)),
|
---|
460 | sourcewindow);
|
---|
461 | }
|
---|
462 |
|
---|
463 | void SamplingGrid::getDiscreteWindowIndices(
|
---|
464 | size_t _wbegin[NDIM],
|
---|
465 | size_t _wlength[NDIM],
|
---|
466 | size_t _wend[NDIM]) const
|
---|
467 | {
|
---|
468 | const double round_offset =
|
---|
469 | (std::numeric_limits<size_t>::round_style == std::round_toward_zero) ?
|
---|
470 | 0.5 : 0.; // need offset to get to round_toward_nearest behavior
|
---|
471 | for(size_t index=0;index<NDIM;++index) {
|
---|
472 | if (fabs(end[index] - begin[index]) > std::numeric_limits<double>::epsilon()*1e4) {
|
---|
473 | // we refrain from using floor/ceil as the window's starts and ends,
|
---|
474 | // the grids have to be compatible (equal level), should always be on
|
---|
475 | // discrete grid point locations.
|
---|
476 | const double delta = getDeltaPerAxis(index);
|
---|
477 | // delta is conversion factor from box length to discrete length, i.e. number of points
|
---|
478 | _wbegin[index] = (begin_window[index] - begin[index])/delta+round_offset;
|
---|
479 | _wlength[index] = (end_window[index] - begin_window[index])/delta+round_offset;
|
---|
480 | _wend[index] = (end_window[index] - begin[index])/delta+round_offset;
|
---|
481 | } else {
|
---|
482 | _wbegin[index] = 0;
|
---|
483 | _wlength[index] = 0;
|
---|
484 | _wend[index] = 0;
|
---|
485 | }
|
---|
486 | // total is used as safe-guard against loss due to discrete conversion
|
---|
487 | ASSERT( (_wend[index] - _wbegin[index]) == _wlength[index],
|
---|
488 | "SamplingGrid::getDiscreteWindowCopyIndices() - end - begin is not equal to length for "
|
---|
489 | +toString(index)+"th component.");
|
---|
490 | }
|
---|
491 | }
|
---|
492 |
|
---|
493 | void SamplingGrid::getDiscreteWindowOffsets(
|
---|
494 | size_t _pre_offset[NDIM],
|
---|
495 | size_t _post_offset[NDIM],
|
---|
496 | size_t _length[NDIM],
|
---|
497 | size_t _total[NDIM]) const
|
---|
498 | {
|
---|
499 | const double round_offset =
|
---|
500 | (std::numeric_limits<size_t>::round_style == std::round_toward_zero) ?
|
---|
501 | 0.5 : 0.; // need offset to get to round_toward_nearest behavior
|
---|
502 | for(size_t index=0;index<NDIM;++index) {
|
---|
503 | if (fabs(end[index] - begin[index]) > std::numeric_limits<double>::epsilon()*1e4) {
|
---|
504 | // we refrain from using floor/ceil as the window's starts and ends,
|
---|
505 | // the grids have to be compatible (equal level), should always be on
|
---|
506 | // discrete grid point locations.
|
---|
507 | const double delta = getDeltaPerAxis(index);
|
---|
508 | // delta is conversion factor from box length to discrete length, i.e. number of points
|
---|
509 | _pre_offset[index] = (begin_window[index] - begin[index])/delta+round_offset;
|
---|
510 | _post_offset[index] = (end[index] - end_window[index])/delta+round_offset;
|
---|
511 | _length[index] = (end_window[index] - begin_window[index])/delta+round_offset;
|
---|
512 | _total[index] = (end[index] - begin[index])/delta+round_offset;
|
---|
513 | } else {
|
---|
514 | _pre_offset[index] = 0;
|
---|
515 | _post_offset[index] = 0;
|
---|
516 | _length[index] = 0;
|
---|
517 | _total[index] = 0;
|
---|
518 | }
|
---|
519 | // total is used as safe-guard against loss due to discrete conversion
|
---|
520 | ASSERT( (_pre_offset[index] + _post_offset[index]) + _length[index] == _total[index],
|
---|
521 | "SamplingGrid::getDiscreteWindowCopyIndices() - pre, length, post are not equal to total for "
|
---|
522 | +toString(index)+"th component.");
|
---|
523 | }
|
---|
524 | }
|
---|
525 |
|
---|
526 | void SamplingGrid::getDiscreteWindowCopyIndices(
|
---|
527 | const double *larger_wbegin,
|
---|
528 | const double *larger_wend,
|
---|
529 | const double *smaller_wbegin,
|
---|
530 | const double *smaller_wend,
|
---|
531 | size_t *pre_offset,
|
---|
532 | size_t *post_offset,
|
---|
533 | size_t *length,
|
---|
534 | size_t *total) const
|
---|
535 | {
|
---|
536 | const double round_offset =
|
---|
537 | (std::numeric_limits<size_t>::round_style == std::round_toward_zero) ?
|
---|
538 | 0.5 : 0.; // need offset to get to round_toward_nearest behavior
|
---|
539 | for(size_t index=0;index<NDIM;++index) {
|
---|
540 | if (fabs(end[index] - begin[index]) > std::numeric_limits<double>::epsilon()*1e4) {
|
---|
541 | // we refrain from using floor/ceil as the window's starts and ends,
|
---|
542 | // the grids have to be compatible (equal level), should always be on
|
---|
543 | // discrete grid point locations.
|
---|
544 | const double delta = getDeltaPerAxis(index);
|
---|
545 | // delta is conversion factor from box length to discrete length, i.e. number of points
|
---|
546 | pre_offset[index] = (smaller_wbegin[index] - larger_wbegin[index])/delta+round_offset;
|
---|
547 | length[index] = (smaller_wend[index] - smaller_wbegin[index])/delta+round_offset;
|
---|
548 | post_offset[index] = (larger_wend[index] - smaller_wend[index])/delta+round_offset;
|
---|
549 | total[index] = (larger_wend[index] - larger_wbegin[index])/delta+round_offset;
|
---|
550 | } else {
|
---|
551 | pre_offset[index] = 0;
|
---|
552 | length[index] = 0;
|
---|
553 | post_offset[index] = 0;
|
---|
554 | total[index] = 0;
|
---|
555 | }
|
---|
556 | // total is used as safe-guard against loss due to discrete conversion
|
---|
557 | ASSERT( pre_offset[index]+post_offset[index]+length[index] == total[index],
|
---|
558 | "SamplingGrid::getDiscreteWindowCopyIndices() - pre, post, and length don't sum up to total for "
|
---|
559 | +toString(index)+"th component.");
|
---|
560 | }
|
---|
561 | }
|
---|
562 |
|
---|
563 | void SamplingGrid::addWindowOntoWindow(
|
---|
564 | const double larger_wbegin[NDIM],
|
---|
565 | const double larger_wend[NDIM],
|
---|
566 | const double smaller_wbegin[NDIM],
|
---|
567 | const double smaller_wend[NDIM],
|
---|
568 | sampledvalues_t &dest_sampled_grid,
|
---|
569 | const sampledvalues_t &source_sampled_grid,
|
---|
570 | boost::function<void (double &, const double &)> op,
|
---|
571 | enum eLargerWindow larger_window)
|
---|
572 | {
|
---|
573 | #ifndef NDEBUG
|
---|
574 | for(size_t index=0;index<NDIM;++index) {
|
---|
575 | ASSERT( smaller_wbegin[index] >= larger_wbegin[index],
|
---|
576 | "SamplingGrid::addWindowOntoWindow() - given smaller window starts earlier than larger window in component "
|
---|
577 | +toString(index)+".");
|
---|
578 | ASSERT( smaller_wend[index] <= larger_wend[index],
|
---|
579 | "SamplingGrid::addWindowOntoWindow() - given smaller window ends later than larger window in component "
|
---|
580 | +toString(index)+".");
|
---|
581 | }
|
---|
582 | #endif
|
---|
583 | // the only issue are indices
|
---|
584 | size_t pre_offset[NDIM];
|
---|
585 | size_t post_offset[NDIM];
|
---|
586 | size_t length[NDIM];
|
---|
587 | size_t total[NDIM];
|
---|
588 | getDiscreteWindowCopyIndices(
|
---|
589 | larger_wbegin, larger_wend,
|
---|
590 | smaller_wbegin, smaller_wend,
|
---|
591 | pre_offset,
|
---|
592 | post_offset,
|
---|
593 | length,
|
---|
594 | total
|
---|
595 | );
|
---|
596 | // assert that calculated lengths match with given vector sizes
|
---|
597 | #ifndef NDEBUG
|
---|
598 | const size_t calculated_size = length[0]*length[1]*length[2];
|
---|
599 | if (larger_window == destwindow) {
|
---|
600 | ASSERT( calculated_size == source_sampled_grid.size(),
|
---|
601 | "SamplingGrid::addWindowOntoWindow() - not enough source sampled values given: "
|
---|
602 | +toString(calculated_size)+" != "+toString(source_sampled_grid.size())+".");
|
---|
603 | ASSERT( calculated_size <= dest_sampled_grid.size(),
|
---|
604 | "SamplingGrid::addWindowOntoWindow() - not enough sampled values available: "
|
---|
605 | +toString(calculated_size)+" <= "+toString(dest_sampled_grid.size())+".");
|
---|
606 | } else {
|
---|
607 | ASSERT( calculated_size == dest_sampled_grid.size(),
|
---|
608 | "SamplingGrid::addWindowOntoWindow() - not enough dest sampled values given: "
|
---|
609 | +toString(calculated_size)+" != "+toString(dest_sampled_grid.size())+".");
|
---|
610 | ASSERT( calculated_size <= source_sampled_grid.size(),
|
---|
611 | "SamplingGrid::addWindowOntoWindow() - not enough source sampled values available: "
|
---|
612 | +toString(calculated_size)+" <= "+toString(source_sampled_grid.size())+".");
|
---|
613 | }
|
---|
614 | const size_t total_size = total[0]*total[1]*total[2];
|
---|
615 | if (larger_window == destwindow) {
|
---|
616 | ASSERT( total_size == dest_sampled_grid.size(),
|
---|
617 | "SamplingGrid::addWindowOntoWindow() - total size is not equal to number of present dest points: "
|
---|
618 | +toString(total_size)+" != "+toString(dest_sampled_grid.size())+".");
|
---|
619 | } else {
|
---|
620 | ASSERT( total_size == source_sampled_grid.size(),
|
---|
621 | "SamplingGrid::addWindowOntoWindow() - total size is not equal to number of present source points: "
|
---|
622 | +toString(total_size)+" != "+toString(source_sampled_grid.size())+".");
|
---|
623 | }
|
---|
624 | #endif
|
---|
625 | size_t N[NDIM];
|
---|
626 | // size_t counter = 0;
|
---|
627 | sampledvalues_t::iterator destiter = dest_sampled_grid.begin();
|
---|
628 | sampledvalues_t::const_iterator sourceiter = source_sampled_grid.begin();
|
---|
629 | if (larger_window == destwindow)
|
---|
630 | std::advance(destiter, pre_offset[0]*total[1]*total[2]);
|
---|
631 | else
|
---|
632 | std::advance(sourceiter, pre_offset[0]*total[1]*total[2]);
|
---|
633 | for(N[0]=0; N[0] < length[0]; ++N[0]) {
|
---|
634 | if (larger_window == destwindow)
|
---|
635 | std::advance(destiter, pre_offset[1]*total[2]);
|
---|
636 | else
|
---|
637 | std::advance(sourceiter, pre_offset[1]*total[2]);
|
---|
638 | for(N[1]=0; N[1] < length[1]; ++N[1]) {
|
---|
639 | if (larger_window == destwindow)
|
---|
640 | std::advance(destiter, pre_offset[2]);
|
---|
641 | else
|
---|
642 | std::advance(sourceiter, pre_offset[2]);
|
---|
643 | for(N[2]=0; N[2] < length[2]; ++N[2]) {
|
---|
644 | ASSERT( destiter != dest_sampled_grid.end(),
|
---|
645 | "SamplingGrid::addWindowOntoWindow() - destiter is already at end of window.");
|
---|
646 | ASSERT( sourceiter != source_sampled_grid.end(),
|
---|
647 | "SamplingGrid::addWindowOntoWindow() - sourceiter is already at end of window.");
|
---|
648 | op(*destiter, *sourceiter);
|
---|
649 | ++destiter;
|
---|
650 | ++sourceiter;
|
---|
651 | }
|
---|
652 | if (larger_window == destwindow)
|
---|
653 | std::advance(destiter, post_offset[2]);
|
---|
654 | else
|
---|
655 | std::advance(sourceiter, post_offset[2]);
|
---|
656 | }
|
---|
657 | if (larger_window == destwindow)
|
---|
658 | std::advance(destiter, post_offset[1]*total[2]);
|
---|
659 | else
|
---|
660 | std::advance(sourceiter, post_offset[1]*total[2]);
|
---|
661 | }
|
---|
662 | #ifndef NDEBUG
|
---|
663 | if (larger_window == destwindow)
|
---|
664 | std::advance(destiter, post_offset[0]*total[1]*total[2]);
|
---|
665 | else
|
---|
666 | std::advance(sourceiter, post_offset[0]*total[1]*total[2]);
|
---|
667 | ASSERT( destiter == dest_sampled_grid.end(),
|
---|
668 | "SamplingGrid::addWindowOntoWindow() - destiter is not at end of window.");
|
---|
669 | ASSERT( sourceiter == source_sampled_grid.end(),
|
---|
670 | "SamplingGrid::addWindowOntoWindow() - sourceiter is not at end of window.");
|
---|
671 | #endif
|
---|
672 | LOG(8, "DEBUG: Grid after adding other is " << dest_sampled_grid << ".");
|
---|
673 | }
|
---|
674 |
|
---|
675 |
|
---|
676 | bool SamplingGrid::operator==(const SamplingGrid &other) const
|
---|
677 | {
|
---|
678 | bool status =
|
---|
679 | static_cast<const SamplingGridProperties &>(*this)
|
---|
680 | == static_cast<const SamplingGridProperties &>(other);
|
---|
681 | // compare general properties
|
---|
682 | if (status) {
|
---|
683 | // compare windows
|
---|
684 | for (size_t i=0; i<NDIM; ++i) {
|
---|
685 | status &= begin_window[i] == other.begin_window[i];
|
---|
686 | status &= end_window[i] == other.end_window[i];
|
---|
687 | }
|
---|
688 | // compare grids
|
---|
689 | if (status)
|
---|
690 | status &= sampled_grid == other.sampled_grid;
|
---|
691 | }
|
---|
692 | return status;
|
---|
693 | }
|
---|
694 |
|
---|
695 | /** Struct contains a single point with displacements from the
|
---|
696 | * central point and the weight in the restriction.
|
---|
697 | */
|
---|
698 | struct PointWeight_t {
|
---|
699 | PointWeight_t(const int &d1, const int &d2, const int &d3, const double &_weight) :
|
---|
700 | displacement(NDIM),
|
---|
701 | weight(_weight)
|
---|
702 | {
|
---|
703 | displacement[0] = d1; displacement[1] = d2; displacement[2] = d3;
|
---|
704 | }
|
---|
705 | typedef std::vector<int> displacement_t;
|
---|
706 | displacement_t displacement;
|
---|
707 | double weight;
|
---|
708 | };
|
---|
709 |
|
---|
710 | static void getLengthsOfWindow(
|
---|
711 | int _total[NDIM],
|
---|
712 | const SamplingGrid &_grid)
|
---|
713 | {
|
---|
714 | const size_t gridpoints_axis = _grid.getGridPointsPerAxis();
|
---|
715 | static const double round_offset =
|
---|
716 | (std::numeric_limits<size_t>::round_style == std::round_toward_zero) ?
|
---|
717 | 0.5 : 0.; // need offset to get to round_toward_nearest behavior
|
---|
718 | for (size_t index=0; index<NDIM; ++index) {
|
---|
719 | if (fabs(_grid.end[index] - _grid.begin[index]) > std::numeric_limits<double>::epsilon()*1e4) {
|
---|
720 | const double delta = (double)gridpoints_axis/(_grid.end[index] - _grid.begin[index]);
|
---|
721 | _total[index] = delta*(_grid.end_window[index] - _grid.begin_window[index])+round_offset;
|
---|
722 | } else
|
---|
723 | _total[index] = 0;
|
---|
724 | // we can only assert that its atmost the maximum number of grid points
|
---|
725 | ASSERT (_total[index] <= ::pow(2, _grid.level),
|
---|
726 | "SamplingGrid::downsample() - total "+toString(_total[index])
|
---|
727 | +" is not equal or less than 2^level: "+toString(_grid.level));
|
---|
728 | }
|
---|
729 | }
|
---|
730 |
|
---|
731 | //!> stencil for full weight restriction, see vmg's stencils.hpp
|
---|
732 | static const std::vector< PointWeight_t > FullWeightNearestNeighbor =
|
---|
733 | boost::assign::list_of
|
---|
734 | ( PointWeight_t( 0, 0, 0, 0.125) )
|
---|
735 | ( PointWeight_t( 1, 0, 0, 0.0625) )
|
---|
736 | ( PointWeight_t(-1, 0, 0, 0.0625) )
|
---|
737 | ( PointWeight_t( 0, 1, 0, 0.0625) )
|
---|
738 | ( PointWeight_t( 0, -1, 0, 0.0625) )
|
---|
739 | ( PointWeight_t( 0, 0, 1, 0.0625) )
|
---|
740 | ( PointWeight_t( 0, 0, -1, 0.0625) )
|
---|
741 | ( PointWeight_t( 1, 1, 0, 0.03125) )
|
---|
742 | ( PointWeight_t( 1, -1, 0, 0.03125) )
|
---|
743 | ( PointWeight_t(-1, 1, 0, 0.03125) )
|
---|
744 | ( PointWeight_t(-1, -1, 0, 0.03125) )
|
---|
745 | ( PointWeight_t( 0, 1, 1, 0.03125) )
|
---|
746 | ( PointWeight_t( 0, 1, -1, 0.03125) )
|
---|
747 | ( PointWeight_t( 0, -1, 1, 0.03125) )
|
---|
748 | ( PointWeight_t( 0, -1, -1, 0.03125) )
|
---|
749 | ( PointWeight_t( 1, 0, 1, 0.03125) )
|
---|
750 | ( PointWeight_t( 1, 0, -1, 0.03125) )
|
---|
751 | ( PointWeight_t(-1, 0, 1, 0.03125) )
|
---|
752 | ( PointWeight_t(-1, 0, -1, 0.03125) )
|
---|
753 | ( PointWeight_t( 1, 1, 1, 0.015625) )
|
---|
754 | ( PointWeight_t( 1, 1, -1, 0.015625) )
|
---|
755 | ( PointWeight_t( 1, -1, 1, 0.015625) )
|
---|
756 | ( PointWeight_t(-1, 1, 1, 0.015625) )
|
---|
757 | ( PointWeight_t( 1, -1, -1, 0.015625) )
|
---|
758 | ( PointWeight_t(-1, 1, -1, 0.015625) )
|
---|
759 | ( PointWeight_t(-1, -1, 1, 0.015625) )
|
---|
760 | ( PointWeight_t(-1, -1, -1, 0.015625) )
|
---|
761 | ;
|
---|
762 |
|
---|
763 | int getValidIndex(
|
---|
764 | const PointWeight_t::displacement_t &_disp,
|
---|
765 | const int N[NDIM],
|
---|
766 | const int length[NDIM])
|
---|
767 | {
|
---|
768 | int index = 0;
|
---|
769 | // we simply truncate in case of out of bounds access
|
---|
770 | if ((N[2]+_disp[2] >= 0) && (N[2]+_disp[2] < length[2]))
|
---|
771 | index += _disp[2];
|
---|
772 | if ((N[1]+_disp[1] >= 0) && (N[1]+_disp[1] < length[1]))
|
---|
773 | index += _disp[1]*length[2];
|
---|
774 | if ((N[0]+_disp[0] >= 0) && (N[0]+_disp[0] < length[0]))
|
---|
775 | index += _disp[0]*length[1]*length[2];
|
---|
776 | return index;
|
---|
777 | }
|
---|
778 |
|
---|
779 | void restrictFullWeight(
|
---|
780 | SamplingGrid::sampledvalues_t &_coarse_level,
|
---|
781 | const int length_c[NDIM],
|
---|
782 | const SamplingGrid::sampledvalues_t &_fine_level,
|
---|
783 | const int length_f[NDIM])
|
---|
784 | {
|
---|
785 | int N_c[NDIM];
|
---|
786 | int N_f[NDIM];
|
---|
787 | SamplingGrid::sampledvalues_t::iterator coarseiter = _coarse_level.begin();
|
---|
788 | for(N_c[0]=0, N_f[0]=0; (N_c[0] < length_c[0]) && (N_f[0] < length_f[0]); ++N_c[0], N_f[0] +=2) {
|
---|
789 | for(N_c[1]=0, N_f[1]=0; (N_c[1] < length_c[1]) && (N_f[1] < length_f[1]); ++N_c[1], N_f[1] +=2) {
|
---|
790 | for(N_c[2]=0, N_f[2]=0; (N_c[2] < length_c[2]) && (N_f[2] < length_f[2]); ++N_c[2], N_f[2] +=2) {
|
---|
791 | const int index_base = N_f[2] + (N_f[1] + N_f[0]*length_f[1])*length_f[2];
|
---|
792 | // go through stencil and add each point relative to displacement with weight
|
---|
793 | for (std::vector< PointWeight_t >::const_iterator weightiter = FullWeightNearestNeighbor.begin();
|
---|
794 | weightiter != FullWeightNearestNeighbor.end(); ++weightiter) {
|
---|
795 | const PointWeight_t::displacement_t disp = weightiter->displacement;
|
---|
796 | const int index_disp = getValidIndex(disp, N_f, length_f);
|
---|
797 | *coarseiter += _fine_level[index_base+index_disp]*weightiter->weight;
|
---|
798 | }
|
---|
799 | ++coarseiter;
|
---|
800 | }
|
---|
801 | ASSERT ( (N_c[2] == length_c[2]) && (N_f[2] == length_f[2]),
|
---|
802 | "restrictFullWeight() - N_c "+toString(N_c[2])+" != length_c "+toString(length_c[2])
|
---|
803 | +" or N_f "+toString(N_f[2])+" != length_f "+toString(length_f[2]));
|
---|
804 | }
|
---|
805 | ASSERT ( (N_c[1] == length_c[1]) && (N_f[1] == length_f[1]),
|
---|
806 | "restrictFullWeight() - N_c "+toString(N_c[1])+" != length_c "+toString(length_c[1])
|
---|
807 | +" or N_f "+toString(N_f[1])+" != length_f "+toString(length_f[1]));
|
---|
808 | }
|
---|
809 | ASSERT ( (N_c[0] == length_c[0]) && (N_f[0] == length_f[0]),
|
---|
810 | "restrictFullWeight() - N_c "+toString(N_c[0])+" != length_c "+toString(length_c[0])
|
---|
811 | +" or N_f "+toString(N_f[0])+" != length_f "+toString(length_f[0]));
|
---|
812 | ASSERT( coarseiter == _coarse_level.end(),
|
---|
813 | "restrictFullWeight() - coarseiter is not at end of values.");
|
---|
814 | }
|
---|
815 |
|
---|
816 | void SamplingGrid::padWithZerosForEvenNumberedSamples()
|
---|
817 | {
|
---|
818 | size_t wbegin_index[NDIM];
|
---|
819 | size_t wend_index[NDIM];
|
---|
820 | size_t wlength_index[NDIM];
|
---|
821 | getDiscreteWindowIndices(wbegin_index, wlength_index, wend_index);
|
---|
822 |
|
---|
823 | // calculate new window (always extend it such that both indices are even)
|
---|
824 | bool changed = false;
|
---|
825 | size_t wnewbegin_index[NDIM];
|
---|
826 | size_t wnewend_index[NDIM];
|
---|
827 | for(size_t i=0;i<NDIM;++i) {
|
---|
828 | LOG(2, "INFO: window[" << i << "] starts at " << wbegin_index[i] << " and ends at "
|
---|
829 | << wend_index[i] << " with length " << wlength_index[i]);
|
---|
830 | // We require begin and end of window on even indices (and inside grid).
|
---|
831 | wnewbegin_index[i] = wbegin_index[i];
|
---|
832 | if ((wnewbegin_index[i] % (size_t)2) != 0) {
|
---|
833 | if (wnewbegin_index[i] > 0)
|
---|
834 | --wnewbegin_index[i];
|
---|
835 | else
|
---|
836 | wnewbegin_index[i] = 0;
|
---|
837 | changed = true;
|
---|
838 | }
|
---|
839 | wnewend_index[i] = wend_index[i];
|
---|
840 | if ((wnewend_index[i] % (size_t)2) != 0) {
|
---|
841 | if (wnewend_index[i] < getGridPointsPerAxis())
|
---|
842 | ++wnewend_index[i];
|
---|
843 | else
|
---|
844 | wnewend_index[i] = getGridPointsPerAxis();
|
---|
845 | changed = true;
|
---|
846 | }
|
---|
847 | ASSERT( (wbegin_index[i] >= 0) && (wend_index[i] <= getGridPointsPerAxis()),
|
---|
848 | "SamplingGrid::padWithZerosForEvenNumberedSamples() - indices "
|
---|
849 | +toString(wbegin_index[i])+" and "+toString(wend_index[i])+" larger than grid "
|
---|
850 | +toString(getGridPointsPerAxis())+".");
|
---|
851 | }
|
---|
852 | if (changed) {
|
---|
853 | double begin_newwindow[NDIM];
|
---|
854 | double end_newwindow[NDIM];
|
---|
855 | for(size_t i=0;i<NDIM;++i) {
|
---|
856 | const double delta = getDeltaPerAxis(i);
|
---|
857 | begin_newwindow[i] = begin[i]+delta*wnewbegin_index[i];
|
---|
858 | end_newwindow[i] = begin[i]+delta*wnewend_index[i];
|
---|
859 | }
|
---|
860 | LOG(2, "INFO: Original window is " << Vector(begin_window) << " <-> " << Vector(end_window));
|
---|
861 | LOG(2, "INFO: Padded window is " << Vector(begin_newwindow) << " <-> " << Vector(end_newwindow));
|
---|
862 | // extend window
|
---|
863 | extendWindow(begin_newwindow, end_newwindow);
|
---|
864 | }
|
---|
865 | ASSERT( getWindowGridPoints() % (size_t)8 == 0,
|
---|
866 | "SamplingGrid::padWithZerosForEvenNumberedSamples() - new size "
|
---|
867 | +toString(getWindowGridPoints())+" is still not divisible by 8.");
|
---|
868 | }
|
---|
869 |
|
---|
870 | void SamplingGrid::downsample(
|
---|
871 | SamplingGrid& instance,
|
---|
872 | const SamplingGrid& other,
|
---|
873 | const int _level)
|
---|
874 | {
|
---|
875 | if (&instance != &other) {
|
---|
876 | LOG(2, "INFO: other's window is " << Vector(other.begin_window)
|
---|
877 | << " <-> " << Vector(other.end_window));
|
---|
878 | // take over properties
|
---|
879 | static_cast<SamplingGridProperties &>(instance) = other;
|
---|
880 | instance.setWindowSize(other.begin_window, other.end_window);
|
---|
881 | LOG(2, "INFO: Set instance's window to " << Vector(instance.begin_window)
|
---|
882 | << " <-> " << Vector(instance.end_window));
|
---|
883 | ASSERT( _level <= other.level,
|
---|
884 | "SamplingGrid::downsample() - desired level "+toString(_level)
|
---|
885 | +" larger than level "+toString(other.level)+" of the given values.");
|
---|
886 | if (_level == other.level) {
|
---|
887 | instance.sampled_grid = other.sampled_grid;
|
---|
888 | } else {
|
---|
889 | // if desired level is smaller we need to downsample
|
---|
890 | // we do this similarly to vmg::RestrictionFullWeight (i.e. a full nearest
|
---|
891 | // neighbor interpolation) and always one grid level at a time till we
|
---|
892 | // have reached the desired one
|
---|
893 |
|
---|
894 | // we need to copy the other grid because we might need to pad it with zeros anyway
|
---|
895 | SamplingGrid FinerGrid(other);
|
---|
896 | int length_d[NDIM];
|
---|
897 | int length_s[NDIM];
|
---|
898 | for (instance.level = other.level-1; instance.level >= _level; --instance.level) {
|
---|
899 | // pad with zeros for even indices and get length of fine grid window
|
---|
900 | FinerGrid.padWithZerosForEvenNumberedSamples();
|
---|
901 | getLengthsOfWindow(length_s, FinerGrid);
|
---|
902 | ASSERT( FinerGrid.getWindowGridPoints() % 8 == 0,
|
---|
903 | "SamplingGrid::downsample() - at level "+toString( instance.level)
|
---|
904 | +" given grid points "+toString(FinerGrid.getWindowGridPoints())+" are not even numbered per axis anymore.");
|
---|
905 | // re-adjust the window (length), get the corresponding window length and downsample
|
---|
906 | instance.setWindow(FinerGrid.begin_window, FinerGrid.end_window);
|
---|
907 | getLengthsOfWindow(length_d, instance);
|
---|
908 | ASSERT( instance.sampled_grid.size() == FinerGrid.getWindowGridPoints()/(size_t)8,
|
---|
909 | "SamplingGrid::downsample() - at level "+toString( instance.level)
|
---|
910 | +" points on coarser grid "+toString(instance.sampled_grid.size())
|
---|
911 | +" and downsampled number on finer grid "
|
---|
912 | +toString(FinerGrid.getWindowGridPoints()/(size_t)8)+" do not match.");
|
---|
913 | restrictFullWeight(instance.sampled_grid, length_d, FinerGrid.sampled_grid, length_s);
|
---|
914 | // then use as new finer grid for next downsampling (if it's not the last level)
|
---|
915 | if (instance.level > _level)
|
---|
916 | FinerGrid = instance;
|
---|
917 | }
|
---|
918 | // loop stops at _level-1
|
---|
919 | instance.level = _level;
|
---|
920 |
|
---|
921 | // and finally, renormalize downsampled grid to old value
|
---|
922 | // instance *= other.integral()/instance.integral();
|
---|
923 | }
|
---|
924 | }
|
---|
925 | }
|
---|
926 |
|
---|
927 | std::ostream & operator<<(std::ostream &ost, const SamplingGrid& other)
|
---|
928 | {
|
---|
929 | ost << "SamplingGrid";
|
---|
930 | ost << " starting at " << other.begin[0] << "," << other.begin[1] << "," << other.begin[2];
|
---|
931 | ost << " ending at " << other.end[0] << "," << other.end[1] << "," << other.end[2];
|
---|
932 | ost << ", window starting at " << other.begin_window[0] << "," << other.begin_window[1] << "," << other.begin_window[2];
|
---|
933 | ost << ", window ending at " << other.end_window[0] << "," << other.end_window[1] << "," << other.end_window[2];
|
---|
934 | ost << ", level of " << other.level;
|
---|
935 | ost << " and integrated value of " << other.integral();
|
---|
936 | return ost;
|
---|
937 | }
|
---|
938 |
|
---|
939 | template<> SamplingGrid ZeroInstance<SamplingGrid>()
|
---|
940 | {
|
---|
941 | SamplingGrid returnvalue;
|
---|
942 | return returnvalue;
|
---|
943 | }
|
---|
944 |
|
---|
945 | // we need to explicitly instantiate the serialization functions as
|
---|
946 | // its is only serialized through its base class FragmentJob
|
---|
947 | BOOST_CLASS_EXPORT_IMPLEMENT(SamplingGrid)
|
---|