source: src/Jobs/Grid/SamplingGrid.cpp@ 98f8fe

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2012 University of Bonn. 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 <http://www.gnu.org/licenses/>.
 */

/*
 * SamplingGrid.cpp
 *
 *  Created on: 25.07.2012
 *      Author: heber
 */

// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

// include headers that implement a archive in simple text format
// otherwise BOOST_CLASS_EXPORT_IMPLEMENT has no effect
#include <boost/archive/text_oarchive.hpp>
#include <boost/archive/text_iarchive.hpp>

#include "CodePatterns/MemDebug.hpp"

#include "Jobs/Grid/SamplingGrid.hpp"

#include <limits>

#include "CodePatterns/Assert.hpp"
#include "CodePatterns/Log.hpp"

SamplingGrid::SamplingGrid(const double _begin[3],
    const double _end[3],
    const int _level,
    const sampledvalues_t &_sampled_grid) :
  SamplingGridProperties(_begin, _end, _level),
  sampled_grid(_sampled_grid)
{}

SamplingGrid::SamplingGrid(const SamplingGrid &_grid) :
  SamplingGridProperties(_grid),
  sampled_grid(_grid.sampled_grid)
{}

SamplingGrid::SamplingGrid(const SamplingGridProperties &_props) :
  SamplingGridProperties(_props)
{}

SamplingGrid::SamplingGrid(
    const SamplingGridProperties &_props,
    const sampledvalues_t &_sampled_grid) :
  SamplingGridProperties(_props),
  sampled_grid(_sampled_grid)
{}

SamplingGrid::~SamplingGrid()
{}

SamplingGrid& SamplingGrid::operator=(const SamplingGrid& other)
{
  // check for self-assignment
  if (this != &other) {
    static_cast<SamplingGridProperties &>(*this) = other;
    sampled_grid = other.sampled_grid;
  }
  return *this;
}

SamplingGrid& SamplingGrid::operator*=(const SamplingGrid& other)
{
  // check that grids are compatible
  if (isCompatible(other)) {
    sampledvalues_t::iterator iter = sampled_grid.begin();
    sampledvalues_t::const_iterator otheriter = other.sampled_grid.begin();
    for(; iter != sampled_grid.end(); ++iter, ++otheriter )
      *iter *= (*otheriter);
  } else {
    ASSERT(0, "SamplingGrid::operator*=() - superposing incompatible grids is so far not in the cards.");
  }
  return *this;
}

void SamplingGrid::superposeOtherGrids(const SamplingGrid &other, const double prefactor)
{
  // check that grids are compatible
  if (isCompatible(other)) {
    // then add other's grid onto this one
    sampledvalues_t::iterator iter = sampled_grid.begin();
    sampledvalues_t::const_iterator otheriter = other.sampled_grid.begin();
    for(; iter != sampled_grid.end(); ++iter, ++otheriter )
      *iter += prefactor * (*otheriter);
  } else {
    ASSERT(0, "SamplingGrid::superposeOtherGrids() - superposing incompatible grids is so far not in the cards.");
  }
}

double SamplingGrid::getVolume() const
{
  double volume = 1.;
  for (size_t i=0;i<3;++i)
    volume *= end[i]-begin[i];
  return volume;
}

const size_t SamplingGrid::getTotalGridPoints() const
{
  return pow(getGridPointsPerAxis(),3);
}

const size_t SamplingGrid::getGridPointsPerAxis() const
{
  return pow(2, level);
}

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 (isCompatible(weight)) {
    const size_t gridpoints = pow(2, level);
    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.;
}

std::ostream & operator<<(std::ostream &ost, const SamplingGrid& other)
{
  ost << "SamplingGrid starting at " << other.begin[0] << "," << other.begin[1] << "," << other.begin[2];
  ost << ", ending at " << other.end[0] << "," << other.end[1] << "," << other.end[2];
  ost << ", level of " << other.level;
  ost << " and integrated value of " << other.integral();
  return ost;
}

template<> SamplingGrid ZeroInstance<SamplingGrid>()
{
  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)