/*
* Project: MoleCuilder
* Description: creates and alters molecular systems
* Copyright (C) 2012 University of Bonn. All rights reserved.
* Copyright (C) 2013 Frederik Heber. All rights reserved.
*
*
* This file is part of MoleCuilder.
*
* MoleCuilder is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 2 of the License, or
* (at your option) any later version.
*
* MoleCuilder is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with MoleCuilder. If not, see .
*/
/*
* SamplingGrid.cpp
*
* Created on: 25.07.2012
* Author: heber
*/
// include config.h
#ifdef HAVE_CONFIG_H
#include
#endif
// include headers that implement a archive in simple text format
// otherwise BOOST_CLASS_EXPORT_IMPLEMENT has no effect
#include
#include
#include "CodePatterns/MemDebug.hpp"
#include "Fragmentation/Summation/SetValues/SamplingGrid.hpp"
#include
#include
#include
#include
#include "CodePatterns/Assert.hpp"
#include "CodePatterns/Log.hpp"
// static instances
const double SamplingGrid::zeroOffset[NDIM] = { 0., 0., 0. };
SamplingGrid::SamplingGrid() :
SamplingGridProperties()
{
setWindowSize(zeroOffset, zeroOffset);
ASSERT( getWindowGridPoints() == (size_t)0,
"SamplingGrid::SamplingGrid() - incorrect number of samples given for the window.");
}
SamplingGrid::SamplingGrid(const double _begin[NDIM],
const double _end[NDIM],
const int _level) :
SamplingGridProperties(_begin, _end, _level)
{
setWindowSize(zeroOffset, zeroOffset);
ASSERT( getWindowGridPoints() == (size_t)0,
"SamplingGrid::SamplingGrid() - incorrect number of samples given for the window.");
}
SamplingGrid::SamplingGrid(const double _begin[NDIM],
const double _end[NDIM],
const int _level,
const sampledvalues_t &_sampled_grid) :
SamplingGridProperties(_begin, _end, _level),
sampled_grid(_sampled_grid)
{
setWindowSize(_begin, _end);
ASSERT( getWindowGridPoints() == (size_t)_sampled_grid.size(),
"SamplingGrid::SamplingGrid() - incorrect number of samples given for the window.");
}
SamplingGrid::SamplingGrid(const SamplingGrid &_grid) :
SamplingGridProperties(_grid),
sampled_grid(_grid.sampled_grid)
{
setWindowSize(_grid.begin_window, _grid.end_window);
ASSERT( getWindowGridPoints() == _grid.getWindowGridPoints(),
"SamplingGrid::SamplingGrid() - incorrect number of samples given for the window.");
}
SamplingGrid::SamplingGrid(const SamplingGridProperties &_props) :
SamplingGridProperties(_props)
{
setWindowSize(zeroOffset, zeroOffset);
ASSERT( getWindowGridPoints() == (size_t)0,
"SamplingGrid::SamplingGrid() - incorrect number of samples given for the window.");
}
SamplingGrid::SamplingGrid(
const SamplingGridProperties &_props,
const sampledvalues_t &_sampled_grid) :
SamplingGridProperties(_props),
sampled_grid(_sampled_grid)
{
setWindowSize(_props.begin, _props.end);
ASSERT( getWindowGridPoints() == (size_t)_sampled_grid.size(),
"SamplingGrid::SamplingGrid() - incorrect number of samples given for the window.");
}
SamplingGrid::~SamplingGrid()
{}
bool SamplingGrid::isCongruent(const SamplingGrid &_props) const
{
bool status = true;
status &= (static_cast(*this) ==
static_cast(_props));
for(size_t i = 0; i(*this) = other;
setWindowSize(other.begin_window, other.end_window);
sampled_grid = other.sampled_grid;
}
return *this;
}
static void multiplyElements(
double &dest,
const double &source,
const double prefactor)
{
dest *= prefactor*(source);
}
SamplingGrid& SamplingGrid::operator*=(const double _value)
{
std::transform(
sampled_grid.begin(), sampled_grid.end(),
sampled_grid.begin(),
boost::bind(std::multiplies(), _1, _value));
return *this;
}
SamplingGrid& SamplingGrid::operator*=(const SamplingGrid& other)
{
// check that grids are compatible
if (isEquivalent(other)) {
/// get minimum of window
double min_begin_window[NDIM];
double min_end_window[NDIM];
bool doShrink = false;
for (size_t index=0; index= other.begin_window[index]) {
max_begin_window[index] = other.begin_window[index];
doExtend = true;
} else {
max_begin_window[index] = begin_window[index];
}
if (end_window[index] >= other.end_window[index]) {
max_end_window[index] = end_window[index];
} else {
max_end_window[index] = other.end_window[index];
doExtend = true;
}
}
LOG(4, "DEBUG: max begin is " << max_begin_window[0] << "," << max_begin_window[1] << "," << max_begin_window[2] << ".");
LOG(4, "DEBUG: max end is " << max_end_window[0] << "," << max_end_window[1] << "," << max_end_window[2] << ".");
if (doExtend)
extendWindow(max_begin_window, max_end_window);
/// and copy other into larger window, too
addOntoWindow(other.begin_window, other.end_window, other.sampled_grid, prefactor);
} else {
ASSERT(0, "SamplingGrid::superposeOtherGrids() - superposing incompatible grids is so far not in the cards.");
}
}
const size_t SamplingGrid::getWindowGridPointsPerAxis(const size_t axis) const
{
static const double round_offset(
(std::numeric_limits::round_style == std::round_toward_zero) ?
0.5 : 0.); // need offset to get to round_toward_nearest behavior
// const double total = getTotalLengthPerAxis(axis);
const double delta = getDeltaPerAxis(axis);
if (delta == 0)
return 0;
const double length = getWindowLengthPerAxis(axis);
if (length == 0)
return 0;
return (size_t)(length/delta+round_offset);
}
double SamplingGrid::integral() const
{
const double volume_element = getVolume()/(double)getTotalGridPoints();
double int_value = 0.;
for (sampledvalues_t::const_iterator iter = sampled_grid.begin();
iter != sampled_grid.end();
++iter)
int_value += *iter;
int_value *= volume_element;
LOG(2, "DEBUG: SamplingGrid::integral() is " << scientific << setprecision(13) << int_value << ".");
return int_value;
}
double SamplingGrid::integral(const SamplingGrid &weight) const
{
if (isEquivalent(weight)) {
const double volume_element = getVolume()/(double)getTotalGridPoints();
double int_value = 0.;
sampledvalues_t::const_iterator iter = sampled_grid.begin();
sampledvalues_t::const_iterator weightiter = weight.sampled_grid.begin();
for (;iter != sampled_grid.end();++iter,++weightiter)
int_value += (*weightiter) * (*iter);
int_value *= volume_element;
//LOG(2, "DEBUG: SamplingGrid::integral() is " << scientific << setprecision(13) << int_value << ".");
return int_value;
} else
return 0.;
}
void SamplingGrid::setWindowSize(
const double _begin_window[NDIM],
const double _end_window[NDIM])
{
for (size_t index=0;index= begin[index],
"SamplingGrid::setWindowSize() - window starts earlier than domain on "
+toString(index)+"th component.");
end_window[index] = getNearestHigherGridPoint(_end_window[index], index);
ASSERT( end_window[index] <= end[index],
"SamplingGrid::setWindowSize() - window ends later than domain on "
+toString(index)+"th component.");
}
}
void SamplingGrid::setWindow(
const double _begin_window[NDIM],
const double _end_window[NDIM])
{
setWindowSize(_begin_window, _end_window);
const size_t gridpoints_window = getWindowGridPoints();
sampled_grid.clear();
sampled_grid.resize(gridpoints_window, 0.);
}
void SamplingGrid::setDomain(
const double _begin[NDIM],
const double _end[NDIM])
{
setDomainSize(_begin, _end);
setWindowSize(_begin, _end);
const size_t gridpoints = getTotalGridPoints();
sampled_grid.resize(gridpoints, 0.);
}
void SamplingGrid::extendWindow(
const double _begin_window[NDIM],
const double _end_window[NDIM])
{
#ifndef NDEBUG
for(size_t index=0;index < NDIM; ++index) {
// check that we truly have to extend the window
ASSERT ( begin_window[index] >= _begin_window[index],
"SamplingGrid::extendWindow() - component "+toString(index)+
" of window start is greater than old value.");
ASSERT ( end_window[index] <= _end_window[index],
"SamplingGrid::extendWindow() - component "+toString(index)+
" of window end is less than old value.");
// check that we are still less than domain
ASSERT ( _begin_window[index] >= begin[index],
"SamplingGrid::extendWindow() - component "+toString(index)+
" of window start is less than domain start.");
ASSERT ( _end_window[index] <= end[index],
"SamplingGrid::extendWindow() - component "+toString(index)+
" of window end is greater than domain end.");
}
#endif
// copy old window size and values
double old_begin_window[NDIM];
double old_end_window[NDIM];
for(size_t index=0;index= _end_window[index],
"SamplingGrid::shrinkWindow() - component "+toString(index)+
" of window end is greater than old value.");
// check that we are still less than domain
ASSERT ( _begin_window[index] >= begin[index],
"SamplingGrid::shrinkWindow() - component "+toString(index)+
" of window start is less than domain start.");
ASSERT ( _end_window[index] <= end[index],
"SamplingGrid::shrinkWindow() - component "+toString(index)+
" of window end is greater than domain end.");
}
#endif
// copy old window size and values
double old_begin_window[NDIM];
double old_end_window[NDIM];
for(size_t index=0;index::round_style == std::round_toward_zero) ?
0.5 : 0.; // need offset to get to round_toward_nearest behavior
for(size_t index=0;index std::numeric_limits::epsilon()*1e4) {
// we refrain from using floor/ceil as the window's starts and ends,
// the grids have to be compatible (equal level), should always be on
// discrete grid point locations.
const double delta = getDeltaPerAxis(index);
// delta is conversion factor from box length to discrete length, i.e. number of points
_wbegin[index] = (begin_window[index] - begin[index])/delta+round_offset;
_wlength[index] = (end_window[index] - begin_window[index])/delta+round_offset;
_wend[index] = (end_window[index] - begin[index])/delta+round_offset;
} else {
_wbegin[index] = 0;
_wlength[index] = 0;
_wend[index] = 0;
}
// total is used as safe-guard against loss due to discrete conversion
ASSERT( (_wend[index] - _wbegin[index]) == _wlength[index],
"SamplingGrid::getDiscreteWindowCopyIndices() - end - begin is not equal to length for "
+toString(index)+"th component.");
}
}
void SamplingGrid::getDiscreteWindowOffsets(
size_t _pre_offset[NDIM],
size_t _post_offset[NDIM],
size_t _length[NDIM],
size_t _total[NDIM]) const
{
const double round_offset =
(std::numeric_limits::round_style == std::round_toward_zero) ?
0.5 : 0.; // need offset to get to round_toward_nearest behavior
for(size_t index=0;index std::numeric_limits::epsilon()*1e4) {
// we refrain from using floor/ceil as the window's starts and ends,
// the grids have to be compatible (equal level), should always be on
// discrete grid point locations.
const double delta = getDeltaPerAxis(index);
// delta is conversion factor from box length to discrete length, i.e. number of points
_pre_offset[index] = (begin_window[index] - begin[index])/delta+round_offset;
_post_offset[index] = (end[index] - end_window[index])/delta+round_offset;
_length[index] = (end_window[index] - begin_window[index])/delta+round_offset;
_total[index] = (end[index] - begin[index])/delta+round_offset;
} else {
_pre_offset[index] = 0;
_post_offset[index] = 0;
_length[index] = 0;
_total[index] = 0;
}
// total is used as safe-guard against loss due to discrete conversion
ASSERT( (_pre_offset[index] + _post_offset[index]) + _length[index] == _total[index],
"SamplingGrid::getDiscreteWindowCopyIndices() - pre, length, post are not equal to total for "
+toString(index)+"th component.");
}
}
void SamplingGrid::getDiscreteWindowCopyIndices(
const double *larger_wbegin,
const double *larger_wend,
const double *smaller_wbegin,
const double *smaller_wend,
size_t *pre_offset,
size_t *post_offset,
size_t *length,
size_t *total) const
{
const double round_offset =
(std::numeric_limits::round_style == std::round_toward_zero) ?
0.5 : 0.; // need offset to get to round_toward_nearest behavior
for(size_t index=0;index std::numeric_limits::epsilon()*1e4) {
// we refrain from using floor/ceil as the window's starts and ends,
// the grids have to be compatible (equal level), should always be on
// discrete grid point locations.
const double delta = getDeltaPerAxis(index);
// delta is conversion factor from box length to discrete length, i.e. number of points
pre_offset[index] = (smaller_wbegin[index] - larger_wbegin[index])/delta+round_offset;
length[index] = (smaller_wend[index] - smaller_wbegin[index])/delta+round_offset;
post_offset[index] = (larger_wend[index] - smaller_wend[index])/delta+round_offset;
total[index] = (larger_wend[index] - larger_wbegin[index])/delta+round_offset;
} else {
pre_offset[index] = 0;
length[index] = 0;
post_offset[index] = 0;
total[index] = 0;
}
// total is used as safe-guard against loss due to discrete conversion
ASSERT( pre_offset[index]+post_offset[index]+length[index] == total[index],
"SamplingGrid::getDiscreteWindowCopyIndices() - pre, post, and length don't sum up to total for "
+toString(index)+"th component.");
}
}
void SamplingGrid::addWindowOntoWindow(
const double larger_wbegin[NDIM],
const double larger_wend[NDIM],
const double smaller_wbegin[NDIM],
const double smaller_wend[NDIM],
sampledvalues_t &dest_sampled_grid,
const sampledvalues_t &source_sampled_grid,
boost::function op,
enum eLargerWindow larger_window)
{
#ifndef NDEBUG
for(size_t index=0;index= larger_wbegin[index],
"SamplingGrid::addWindowOntoWindow() - given smaller window starts earlier than larger window in component "
+toString(index)+".");
ASSERT( smaller_wend[index] <= larger_wend[index],
"SamplingGrid::addWindowOntoWindow() - given smaller window ends later than larger window in component "
+toString(index)+".");
}
#endif
// the only issue are indices
size_t pre_offset[NDIM];
size_t post_offset[NDIM];
size_t length[NDIM];
size_t total[NDIM];
getDiscreteWindowCopyIndices(
larger_wbegin, larger_wend,
smaller_wbegin, smaller_wend,
pre_offset,
post_offset,
length,
total
);
// assert that calculated lengths match with given vector sizes
#ifndef NDEBUG
const size_t calculated_size = length[0]*length[1]*length[2];
if (larger_window == destwindow) {
ASSERT( calculated_size == source_sampled_grid.size(),
"SamplingGrid::addWindowOntoWindow() - not enough source sampled values given: "
+toString(calculated_size)+" != "+toString(source_sampled_grid.size())+".");
ASSERT( calculated_size <= dest_sampled_grid.size(),
"SamplingGrid::addWindowOntoWindow() - not enough sampled values available: "
+toString(calculated_size)+" <= "+toString(dest_sampled_grid.size())+".");
} else {
ASSERT( calculated_size == dest_sampled_grid.size(),
"SamplingGrid::addWindowOntoWindow() - not enough dest sampled values given: "
+toString(calculated_size)+" != "+toString(dest_sampled_grid.size())+".");
ASSERT( calculated_size <= source_sampled_grid.size(),
"SamplingGrid::addWindowOntoWindow() - not enough source sampled values available: "
+toString(calculated_size)+" <= "+toString(source_sampled_grid.size())+".");
}
const size_t total_size = total[0]*total[1]*total[2];
if (larger_window == destwindow) {
ASSERT( total_size == dest_sampled_grid.size(),
"SamplingGrid::addWindowOntoWindow() - total size is not equal to number of present dest points: "
+toString(total_size)+" != "+toString(dest_sampled_grid.size())+".");
} else {
ASSERT( total_size == source_sampled_grid.size(),
"SamplingGrid::addWindowOntoWindow() - total size is not equal to number of present source points: "
+toString(total_size)+" != "+toString(source_sampled_grid.size())+".");
}
#endif
size_t N[NDIM];
// size_t counter = 0;
sampledvalues_t::iterator destiter = dest_sampled_grid.begin();
sampledvalues_t::const_iterator sourceiter = source_sampled_grid.begin();
if (larger_window == destwindow)
std::advance(destiter, pre_offset[0]*total[1]*total[2]);
else
std::advance(sourceiter, pre_offset[0]*total[1]*total[2]);
for(N[0]=0; N[0] < length[0]; ++N[0]) {
if (larger_window == destwindow)
std::advance(destiter, pre_offset[1]*total[2]);
else
std::advance(sourceiter, pre_offset[1]*total[2]);
for(N[1]=0; N[1] < length[1]; ++N[1]) {
if (larger_window == destwindow)
std::advance(destiter, pre_offset[2]);
else
std::advance(sourceiter, pre_offset[2]);
for(N[2]=0; N[2] < length[2]; ++N[2]) {
ASSERT( destiter != dest_sampled_grid.end(),
"SamplingGrid::addWindowOntoWindow() - destiter is already at end of window.");
ASSERT( sourceiter != source_sampled_grid.end(),
"SamplingGrid::addWindowOntoWindow() - sourceiter is already at end of window.");
op(*destiter, *sourceiter);
++destiter;
++sourceiter;
}
if (larger_window == destwindow)
std::advance(destiter, post_offset[2]);
else
std::advance(sourceiter, post_offset[2]);
}
if (larger_window == destwindow)
std::advance(destiter, post_offset[1]*total[2]);
else
std::advance(sourceiter, post_offset[1]*total[2]);
}
#ifndef NDEBUG
if (larger_window == destwindow)
std::advance(destiter, post_offset[0]*total[1]*total[2]);
else
std::advance(sourceiter, post_offset[0]*total[1]*total[2]);
ASSERT( destiter == dest_sampled_grid.end(),
"SamplingGrid::addWindowOntoWindow() - destiter is not at end of window.");
ASSERT( sourceiter == source_sampled_grid.end(),
"SamplingGrid::addWindowOntoWindow() - sourceiter is not at end of window.");
#endif
LOG(8, "DEBUG: Grid after adding other is " << dest_sampled_grid << ".");
}
bool SamplingGrid::operator==(const SamplingGrid &other) const
{
bool status =
static_cast(*this)
== static_cast(other);
// compare general properties
if (status) {
// compare windows
for (size_t i=0; i displacement_t;
displacement_t displacement;
double weight;
};
static void getLengthsOfWindow(
int _total[NDIM],
const SamplingGrid &_grid)
{
const size_t gridpoints_axis = _grid.getGridPointsPerAxis();
static const double round_offset =
(std::numeric_limits::round_style == std::round_toward_zero) ?
0.5 : 0.; // need offset to get to round_toward_nearest behavior
for (size_t index=0; index std::numeric_limits::epsilon()*1e4) {
const double delta = (double)gridpoints_axis/(_grid.end[index] - _grid.begin[index]);
_total[index] = delta*(_grid.end_window[index] - _grid.begin_window[index])+round_offset;
} else
_total[index] = 0;
// we can only assert that its atmost the maximum number of grid points
ASSERT (_total[index] <= ::pow(2, _grid.level),
"SamplingGrid::downsample() - total "+toString(_total[index])
+" is not equal or less than 2^level: "+toString(_grid.level));
}
}
//!> stencil for full weight restriction, see vmg's stencils.hpp
static const std::vector< PointWeight_t > FullWeightNearestNeighbor =
boost::assign::list_of
( PointWeight_t( 0, 0, 0, 0.125) )
( PointWeight_t( 1, 0, 0, 0.0625) )
( PointWeight_t(-1, 0, 0, 0.0625) )
( PointWeight_t( 0, 1, 0, 0.0625) )
( PointWeight_t( 0, -1, 0, 0.0625) )
( PointWeight_t( 0, 0, 1, 0.0625) )
( PointWeight_t( 0, 0, -1, 0.0625) )
( PointWeight_t( 1, 1, 0, 0.03125) )
( PointWeight_t( 1, -1, 0, 0.03125) )
( PointWeight_t(-1, 1, 0, 0.03125) )
( PointWeight_t(-1, -1, 0, 0.03125) )
( PointWeight_t( 0, 1, 1, 0.03125) )
( PointWeight_t( 0, 1, -1, 0.03125) )
( PointWeight_t( 0, -1, 1, 0.03125) )
( PointWeight_t( 0, -1, -1, 0.03125) )
( PointWeight_t( 1, 0, 1, 0.03125) )
( PointWeight_t( 1, 0, -1, 0.03125) )
( PointWeight_t(-1, 0, 1, 0.03125) )
( PointWeight_t(-1, 0, -1, 0.03125) )
( PointWeight_t( 1, 1, 1, 0.015625) )
( PointWeight_t( 1, 1, -1, 0.015625) )
( PointWeight_t( 1, -1, 1, 0.015625) )
( PointWeight_t(-1, 1, 1, 0.015625) )
( PointWeight_t( 1, -1, -1, 0.015625) )
( PointWeight_t(-1, 1, -1, 0.015625) )
( PointWeight_t(-1, -1, 1, 0.015625) )
( PointWeight_t(-1, -1, -1, 0.015625) )
;
int getValidIndex(
const PointWeight_t::displacement_t &_disp,
const int N[NDIM],
const int length[NDIM])
{
int index = 0;
// we simply truncate in case of out of bounds access
if ((N[2]+_disp[2] >= 0) && (N[2]+_disp[2] < length[2]))
index += _disp[2];
if ((N[1]+_disp[1] >= 0) && (N[1]+_disp[1] < length[1]))
index += _disp[1]*length[2];
if ((N[0]+_disp[0] >= 0) && (N[0]+_disp[0] < length[0]))
index += _disp[0]*length[1]*length[2];
return index;
}
void restrictFullWeight(
SamplingGrid::sampledvalues_t &_coarse_level,
const int length_c[NDIM],
const SamplingGrid::sampledvalues_t &_fine_level,
const int length_f[NDIM])
{
int N_c[NDIM];
int N_f[NDIM];
SamplingGrid::sampledvalues_t::iterator coarseiter = _coarse_level.begin();
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) {
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) {
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) {
const int index_base = N_f[2] + (N_f[1] + N_f[0]*length_f[1])*length_f[2];
// go through stencil and add each point relative to displacement with weight
for (std::vector< PointWeight_t >::const_iterator weightiter = FullWeightNearestNeighbor.begin();
weightiter != FullWeightNearestNeighbor.end(); ++weightiter) {
const PointWeight_t::displacement_t disp = weightiter->displacement;
const int index_disp = getValidIndex(disp, N_f, length_f);
*coarseiter += _fine_level[index_base+index_disp]*weightiter->weight;
}
++coarseiter;
}
ASSERT ( (N_c[2] == length_c[2]) && (N_f[2] == length_f[2]),
"restrictFullWeight() - N_c "+toString(N_c[2])+" != length_c "+toString(length_c[2])
+" or N_f "+toString(N_f[2])+" != length_f "+toString(length_f[2]));
}
ASSERT ( (N_c[1] == length_c[1]) && (N_f[1] == length_f[1]),
"restrictFullWeight() - N_c "+toString(N_c[1])+" != length_c "+toString(length_c[1])
+" or N_f "+toString(N_f[1])+" != length_f "+toString(length_f[1]));
}
ASSERT ( (N_c[0] == length_c[0]) && (N_f[0] == length_f[0]),
"restrictFullWeight() - N_c "+toString(N_c[0])+" != length_c "+toString(length_c[0])
+" or N_f "+toString(N_f[0])+" != length_f "+toString(length_f[0]));
ASSERT( coarseiter == _coarse_level.end(),
"restrictFullWeight() - coarseiter is not at end of values.");
}
void SamplingGrid::padWithZerosForEvenNumberedSamples()
{
size_t wbegin_index[NDIM];
size_t wend_index[NDIM];
size_t wlength_index[NDIM];
getDiscreteWindowIndices(wbegin_index, wlength_index, wend_index);
// calculate new window (always extend it such that both indices are even)
bool changed = false;
size_t wnewbegin_index[NDIM];
size_t wnewend_index[NDIM];
for(size_t i=0;i 0)
--wnewbegin_index[i];
else
wnewbegin_index[i] = 0;
changed = true;
}
wnewend_index[i] = wend_index[i];
if ((wnewend_index[i] % (size_t)2) != 0) {
if (wnewend_index[i] < getGridPointsPerAxis())
++wnewend_index[i];
else
wnewend_index[i] = getGridPointsPerAxis();
changed = true;
}
ASSERT( (wbegin_index[i] >= 0) && (wend_index[i] <= getGridPointsPerAxis()),
"SamplingGrid::padWithZerosForEvenNumberedSamples() - indices "
+toString(wbegin_index[i])+" and "+toString(wend_index[i])+" larger than grid "
+toString(getGridPointsPerAxis())+".");
}
if (changed) {
double begin_newwindow[NDIM];
double end_newwindow[NDIM];
for(size_t i=0;i(instance) = other;
instance.setWindowSize(other.begin_window, other.end_window);
ASSERT( _level <= other.level,
"SamplingGrid::downsample() - desired level "+toString(_level)
+" larger than level "+toString(other.level)+" of the given values.");
if (_level == other.level) {
instance.sampled_grid = other.sampled_grid;
} else {
// if desired level is smaller we need to downsample
// we do this similarly to vmg::RestrictionFullWeight (i.e. a full nearest
// neighbor interpolation) and always one grid level at a time till we
// have reached the desired one
// the reference such that we never have to copy the full grid but only
// downsampled ones
const sampledvalues_t * sourcevalues = &other.sampled_grid;
int length_d[3];
int length_s[3];
getLengthsOfWindow(length_s, other);
for (instance.level = other.level-1; instance.level >= _level; --instance.level) {
getLengthsOfWindow(length_d, instance);
// we always have an eighth of the number of sample points as we stop
ASSERT( sourcevalues->size() % 8 == 0,
"SamplingGrid::downsample() - at level "+toString( instance.level)
+" given grid points "+toString(sourcevalues->size())+" are not even numbered per axis anymore.");
sampledvalues_t downsampled(sourcevalues->size()/(size_t)8, 0.);
restrictFullWeight(downsampled, length_d, *sourcevalues, length_s);
// then copy the downsampled values
instance.sampled_grid = downsampled;
sourcevalues = &instance.sampled_grid;
// and exchange lengths
for (size_t i=0;i<3;++i) {
length_s[i] = length_d[i];
}
}
instance.level = _level;
// and finally, renormalize downsampled grid to old value
// instance *= other.integral()/instance.integral();
}
}
}
std::ostream & operator<<(std::ostream &ost, const SamplingGrid& other)
{
ost << "SamplingGrid";
ost << " starting at " << other.begin[0] << "," << other.begin[1] << "," << other.begin[2];
ost << " ending at " << other.end[0] << "," << other.end[1] << "," << other.end[2];
ost << ", window starting at " << other.begin_window[0] << "," << other.begin_window[1] << "," << other.begin_window[2];
ost << ", window ending at " << other.end_window[0] << "," << other.end_window[1] << "," << other.end_window[2];
ost << ", level of " << other.level;
ost << " and integrated value of " << other.integral();
return ost;
}
template<> SamplingGrid ZeroInstance()
{
SamplingGrid returnvalue;
return returnvalue;
}
// we need to explicitly instantiate the serialization functions as
// its is only serialized through its base class FragmentJob
BOOST_CLASS_EXPORT_IMPLEMENT(SamplingGrid)