source: src/Fragmentation/Summation/SetValues/SamplingGridProperties.cpp@ 3d444a

/*
 * 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 <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 "CodePatterns/Assert.hpp"

#include "Fragmentation/Summation/SetValues/SamplingGridProperties.hpp"

SamplingGridProperties::SamplingGridProperties(
    const double _begin[NDIM],
    const double _end[NDIM],
    const int _level) :
  level(_level)
{
  for(size_t i=0; i<NDIM; ++i) {
    begin[i] = _begin[i];
    end[i] = _end[i];
  }
}

SamplingGridProperties::SamplingGridProperties(const SamplingGridProperties &_props) :
  level(_props.level)
{
  for(size_t i=0; i<NDIM; ++i) {
    begin[i] = _props.begin[i];
    end[i] = _props.end[i];
  }
}

SamplingGridProperties::SamplingGridProperties() :
  level(0)
{
  for(size_t i=0; i<NDIM; ++i) {
    begin[i] = 0.;
    end[i] = 0.;
  }
}

SamplingGridProperties::~SamplingGridProperties()
{}

SamplingGridProperties& SamplingGridProperties::operator=(const SamplingGridProperties& other)
{
  // check for self-assignment
  if (this != &other) {
    for(size_t index=0; index<NDIM; ++index) {
      begin[index] = other.begin[index];
      end[index] = other.end[index];
    }
    level = other.level;
  }
  return *this;
}

bool SamplingGridProperties::operator==(const SamplingGridProperties &_props) const
{
  bool status = true;
  for (size_t i=0; i<NDIM; ++i) {
    status &= begin[i] == _props.begin[i];
    status &= end[i] == _props.end[i];
  }
  status &= level == _props.level;
  return status;
}


double SamplingGridProperties::getNearestLowerGridPoint(
    const double value, const size_t axis) const
{
  const double length = getTotalLengthPerAxis(axis);
  if ((fabs(length) < std::numeric_limits<double>::epsilon()) || (getGridPointsPerAxis() == 0))
    return begin[axis];
  if ((value - end[axis]) > -std::numeric_limits<double>::epsilon()*1.e4)
    return end[axis];
  const double offset = value - begin[axis];
  if (offset < std::numeric_limits<double>::epsilon()*1.e4)
    return begin[axis];
  // modify value a little if value actually has been on a grid point but
  // perturbed by numerical rounding: here up, as we always go lower
  const double factor =
      floor((double)getGridPointsPerAxis()*(offset/length)
          +std::numeric_limits<double>::epsilon()*1.e4);
  return begin[axis]+factor*getDeltaPerAxis(axis);
}

double SamplingGridProperties::getNearestHigherGridPoint(
    const double value, const size_t axis) const
{
  const double length = getTotalLengthPerAxis(axis);
  if ((fabs(length) < std::numeric_limits<double>::epsilon()) || (getGridPointsPerAxis() == 0))
    return end[axis];
  if ((value - end[axis]) > -std::numeric_limits<double>::epsilon()*1.e4)
    return end[axis];
  const double offset = value - begin[axis];
  if (offset < std::numeric_limits<double>::epsilon()*1.e4)
    return begin[axis];
  // modify value a little if value actually has been on a grid point but
  // perturbed by numerical rounding: here down, as we always go higher
  const double factor =
      ceil((double)getGridPointsPerAxis()*(offset/length)
          -std::numeric_limits<double>::epsilon()*1.e4);
  return begin[axis]+factor*getDeltaPerAxis(axis);
}

double SamplingGridProperties::getSurplusLevel(const SamplingGridProperties &_props) const
{
  static const double log_two = log(2.);
  const double domain_extent =
      std::max(getTotalLengthPerAxis(0),
          std::max(getTotalLengthPerAxis(1),
              getTotalLengthPerAxis(2)));
  const double props_extent =
      std::max(_props.getTotalLengthPerAxis(0),
          std::max(_props.getTotalLengthPerAxis(1),
              _props.getTotalLengthPerAxis(2)));
  //!> states how many more levels we have because of smaller grid (assuming equal levels)
  const double surplus_level = log(domain_extent/props_extent)/log_two;
  return surplus_level;
}

bool SamplingGridProperties::isCompatible(const SamplingGridProperties &_props) const
{
  // only equally sized or finer grids are compatible: we only downsample
  const double surplus_level = getSurplusLevel(_props);
  if (fabs(surplus_level - round(surplus_level)) > std::numeric_limits<double>::epsilon()*1e4)
    return false;
  if (level > (_props.level+surplus_level))
    return false;
  // check whether grid point corners coincide
  for (size_t i=0;i<NDIM;++i) {
    const double lowercorner = getNearestLowerGridPoint(_props.begin[i], i);
    if (fabs(lowercorner - _props.begin[i]) > std::numeric_limits<double>::epsilon()*1e4)
      return false;
    const double uppercorner = getNearestHigherGridPoint(_props.end[i], i);
    if (fabs(uppercorner - _props.end[i]) > std::numeric_limits<double>::epsilon()*1e4)
      return false;
  }
  return true;
}

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