source: src/Jobs/WindowGrid_converter.cpp@ 30e5dc

/*
 * 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.
 * Please see the COPYING file or "Copyright notice" in builder.cpp for details.
 * 
 *
 *   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/>. 
 */

/*
 * WindowGrid_converter.cpp
 *
 *  Created on: Dec 20, 2012
 *      Author: heber
 */


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

#include "grid/grid.hpp"

#include "WindowGrid_converter.hpp"

#include "CodePatterns/MemDebug.hpp"

#include <iostream>
#include <iterator>
#include <limits>

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

#include "Fragmentation/Summation/SetValues/SamplingGrid.hpp"
#include "Jobs/ChargeSmearer.hpp"

void WindowGrid_converter::addGridOntoWindow(
    VMG::Grid &grid,
    SamplingGrid &window,
    const double prefactor,
    const bool OpenBoundaryConditions)
{
#ifndef NDEBUG
  for(size_t index=0;index<3;++index) {
    ASSERT( window.begin[index] >= window.begin_window[index],
        "InterfaceVMGJob::addGridOntoWindow() - given grid starts earlier than window in component "
        +toString(index)+".");
    ASSERT( window.end[index] <= window.end_window[index],
        "InterfaceVMGJob::addGridOntoWindow() - given grid ends later than window in component "
        +toString(index)+".");
  }
#endif
  // the only issue are indices
  const size_t gridpoints_axis = window.getGridPointsPerAxis();
  size_t pre_offset[3];
  size_t post_offset[3];
  size_t length[3];
  size_t total[3];
  const double round_offset =
      (std::numeric_limits<size_t>::round_style == std::round_toward_zero) ?
          0.5 : 0.; // need offset to get to round_toward_nearest behavior
  for(size_t index=0;index<3;++index) {
    const double delta = (double)gridpoints_axis/(window.end[index] - window.begin[index]);
    // delta is conversion factor from box length to discrete length, i.e. number of points
    pre_offset[index] = delta*(window.begin[index] - window.begin_window[index])+round_offset;
    length[index] = delta*(window.end[index] - window.begin[index])+round_offset;
    post_offset[index] = delta*(window.end_window[index] - window.end[index])+round_offset;
    total[index] = delta*(window.end_window[index] - window.begin_window[index])+round_offset;
    // total is used as safe-guard against loss due to discrete conversion
    ASSERT( pre_offset[index]+post_offset[index]+length[index] == total[index],
        "InterfaceVMGJob::addGridOntoWindow() - pre, post, and length don't sum up to total for "
        +toString(index)+"th component.");
  }
#ifndef NDEBUG
  const size_t calculated_size = length[0]*length[1]*length[2];
//  const size_t given_size = std::distance(grid_begin, grid_end);
//  ASSERT( calculated_size == given_size,
//      "InterfaceVMGJob::addGridOntoWindow() - not enough sampled values given: "
//      +toString(calculated_size)+" != "+toString(given_size)+".");
  ASSERT( calculated_size <=  window.sampled_grid.size(),
      "InterfaceVMGJob::addGridOntoWindow() - not enough sampled values available: "
      +toString(calculated_size)+" <= "+toString(window.sampled_grid.size())+".");
  const size_t total_size = total[0]*total[1]*total[2];
  ASSERT( total_size == window.sampled_grid.size(),
      "InterfaceVMGJob::addGridOntoWindow() - total size is not equal to number of present points: "
      +toString(total_size)+" != "+toString(window.sampled_grid.size())+".");
#endif
  size_t N[3];
  SamplingGrid::sampledvalues_t::iterator griditer = window.sampled_grid.begin();
  std::advance(griditer, pre_offset[0]*total[1]*total[2]);

  VMG::Grid::iterator copyiter = grid.Iterators().Local().Begin();
  const VMG::Index size = grid.Local().Size();
  if (OpenBoundaryConditions)
    for(N[0]=(size_t)0; N[0] < (size_t)size[0]/(size_t)4; ++N[0])
      for(N[1]=(size_t)0; N[1] < (size_t)size[1]; ++N[1])
        for(N[2]=(size_t)0; N[2] < (size_t)size[2]; ++N[2]) {
          ASSERT( copyiter != grid.Iterators().Local().End(),
              "InterfaceVMGJob::addGridOntoWindow() - copyiter is already at end of window.");
          ++copyiter;
        }
  for(N[0]=(size_t)0; N[0] < length[0]; ++N[0]) {
    std::advance(griditer, pre_offset[1]*total[2]);
    if (OpenBoundaryConditions)
      for(N[1]=(size_t)0; N[1] < (size_t)size[1]/(size_t)4; ++N[1])
        for(N[2]=(size_t)0; N[2] < (size_t)size[2]; ++N[2]) {
          ASSERT( copyiter != grid.Iterators().Local().End(),
              "InterfaceVMGJob::addGridOntoWindow() - copyiter is already at end of window.");
          ++copyiter;
        }
    for(N[1]=(size_t)0; N[1] < length[1]; ++N[1]) {
      std::advance(griditer, pre_offset[2]);
      if (OpenBoundaryConditions)
        for(N[2]=(size_t)0; N[2] < (size_t)size[2]/(size_t)4; ++N[2]) {
          ASSERT( copyiter != grid.Iterators().Local().End(),
              "InterfaceVMGJob::addGridOntoWindow() - copyiter is already at end of window.");
          ++copyiter;
        }
      for(N[2]=(size_t)0; N[2] < length[2]; ++N[2]) {
        ASSERT( griditer != window.sampled_grid.end(),
            "InterfaceVMGJob::addGridOntoWindow() - griditer is already at end of window.");
        ASSERT( copyiter != grid.Iterators().Local().End(),
            "InterfaceVMGJob::addGridOntoWindow() - copyiter is already at end of window.");
        *griditer++ += prefactor*grid(*copyiter++);
      }
      std::advance(griditer, post_offset[2]);
      if (OpenBoundaryConditions)
        for(N[2]=(size_t)0; N[2] < (size_t)size[2] - (size_t)size[2]/(size_t)4 - length[2]; ++N[2]) {
          ASSERT( copyiter != grid.Iterators().Local().End(),
              "InterfaceVMGJob::addGridOntoWindow() - copyiter is already at end of window.");
          ++copyiter;
        }
    }
    std::advance(griditer, post_offset[1]*total[2]);
    if (OpenBoundaryConditions)
      for(N[1]=(size_t)0; N[1] < (size_t)size[1] - (size_t)size[1]/(size_t)4 - length[1]; ++N[1])
        for(N[2]=(size_t)0; N[2] < (size_t)size[2]; ++N[2]) {
          ASSERT( copyiter != grid.Iterators().Local().End(),
              "InterfaceVMGJob::addGridOntoWindow() - copyiter is already at end of window.");
          ++copyiter;
        }
  }
  if (OpenBoundaryConditions)
    for(N[0]=(size_t)0; N[0] < (size_t)size[0] - (size_t)size[0]/(size_t)4 - length[0]; ++N[0])
      for(N[1]=(size_t)0; N[1] < (size_t)size[1]; ++N[1])
        for(N[2]=(size_t)0; N[2] < (size_t)size[2]; ++N[2]) {
          ASSERT( copyiter != grid.Iterators().Local().End(),
              "InterfaceVMGJob::addGridOntoWindow() - copyiter is already at end of window.");
          ++copyiter;
        }
#ifndef NDEBUG
  std::advance(griditer, post_offset[0]*total[1]*total[2]);
  ASSERT( griditer == window.sampled_grid.end(),
      "InterfaceVMGJob::addGridOntoWindow() - griditer is not at end of window.");
  ASSERT( copyiter == grid.Iterators().Local().End(),
      "InterfaceVMGJob::addGridOntoWindow() - copyiter is not at end of window.");
#endif
}

void WindowGrid_converter::addWindowOntoGrid(
    VMG::Grid& window,
    const SamplingGrid &grid,
    const double prefactor,
    const bool OpenBoundaryConditions,
    const bool DoSmearCharges)
{
#ifndef NDEBUG
  for(size_t index=0;index<3;++index) {
    ASSERT( grid.begin_window[index] >= grid.begin[index],
        "InterfaceVMGJob::addWindowOntoGrid() - given window starts earlier than grid in component "
        +toString(index)+".");
    ASSERT( grid.end_window[index] <= grid.end[index],
        "InterfaceVMGJob::addWindowOntoGrid() - given window ends later than grid in component "
        +toString(index)+".");
  }
#endif
  // the only issue are indices
  const size_t gridpoints_axis = grid.getGridPointsPerAxis();
  size_t pre_offset[3];
  size_t post_offset[3];
  size_t length[3];
  size_t total[3];
  const double round_offset =
      (std::numeric_limits<size_t>::round_style == std::round_toward_zero) ?
          0.5 : 0.; // need offset to get to round_toward_nearest behavior
  for(size_t index=0;index<3;++index) {
    const double delta = (double)gridpoints_axis/(grid.end[index] - grid.begin[index]);
    // delta is conversion factor from box length to discrete length, i.e. number of points
    pre_offset[index] = delta*(grid.begin_window[index] - grid.begin[index])+round_offset;
    length[index] = delta*(grid.end_window[index] - grid.begin_window[index])+round_offset;
    post_offset[index] = delta*(grid.end[index] - grid.end_window[index])+round_offset;
    total[index] = delta*(grid.end[index] - grid.begin[index])+round_offset;
    // total is used as safe-guard against loss due to discrete conversion
    ASSERT( pre_offset[index]+post_offset[index]+length[index] == total[index],
        "InterfaceVMGJob::addWindowOntoGrid() - pre, post, and length don't sum up to total for "
        +toString(index)+"th component.");
  }
#ifndef NDEBUG
  const size_t calculated_size = length[0]*length[1]*length[2];
  ASSERT( calculated_size == grid.sampled_grid.size(),
      "InterfaceVMGJob::addWindowOntoGrid() - not enough sampled values given: "
      +toString(calculated_size)+" != "+toString(grid.sampled_grid.size())+".");
//  ASSERT( calculated_size <=  given_size,
//      "InterfaceVMGJob::addWindowOntoGrid() - not enough sampled values available: "
//      +toString(calculated_size)+" <= "+toString(given_size)+".");
//  const size_t total_size = total[0]*total[1]*total[2];
//  ASSERT( total_size == given_size,
//      "InterfaceVMGJob::addWindowOntoGrid() - total size is not equal to number of present points: "
//      +toString(total_size)+" != "+toString(given_size)+".");
#endif
  size_t N[3];
  // in open boundary case grid.Local() contains more than just the inner area. The grid
  // is actually twice as large to allow for the FV discretization to work. Hence, we first
  // have to seek our starting point ... see VMG::Grid::GetSpatialPos()
  if (OpenBoundaryConditions) {
    const VMG::Index size = window.Local().Size();
//    const VMG::Index boundary1size = window.Local().BoundarySize1();
//    const VMG::Index boundary2size = window.Local().BoundarySize2();
//    const VMG::Index halo1size = window.Local().HaloSize1();
//    const VMG::Index halo2size = window.Local().HaloSize2();
    // this mimicks VMG::GridIndexTranslations::EndOffset()
    const size_t off = OpenBoundaryConditions ? 1 : 0;
    for (size_t i=0;i<3;++i)
      pre_offset[i] += ((size_t)size[i] - off) / 4;
    for (size_t i=0;i<3;++i)
      total[i] = ((size_t)size[i]);
    for (size_t i=0;i<3;++i)
      post_offset[i] = ((size_t)size[i]) - pre_offset[i] - length[i];
  }

  /// reset grid to zero everywhere

  VMG::Grid::iterator griditer = window.Iterators().Local().Begin();
//  griditer.advance(pre_offset[0]*total[1]*total[2]);
  for(N[0]=0; N[0] < total[0]; ++N[0])
    for(N[1]=0; N[1] < total[1]; ++N[1])
      for(N[2]=0; N[2] < total[2]; ++N[2]) {
        ASSERT( griditer != window.Iterators().Local().End(),
            "InterfaceVMGJob::addWindowOntoGrid() - griditer is already at end of window.");
        window(*griditer++) = 0.;
      }

#ifndef NDEBUG
  ASSERT( griditer == window.Iterators().Local().End(),
      "InterfaceVMGJob::addWindowOntoGrid() - griditer is not at end of window.");
#endif

  /// then apply charges onto grid

  griditer = window.Iterators().Local().Begin();
  SamplingGrid::sampledvalues_t::const_iterator copyiter = grid.sampled_grid.begin();
  const ChargeSmearer &smearer = ChargeSmearer::getInstance();
  // NOTE: GridIterator::operator+=()'s implementation in vmg is broken. DON'T USE!
//  griditer += pre_offset[0] * total[1] * total[2];
  for(N[0]=0; N[0] < pre_offset[0]; ++N[0])
    for(N[1]=0; N[1] < total[1]; ++N[1])
      for(N[2]=0; N[2] < total[2]; ++N[2]) {
        ASSERT( griditer != window.Iterators().Local().End(),
            "InterfaceVMGJob::addWindowOntoGrid() - griditer is already at end of window.");
        griditer++;
      }
  for(N[0]=0; N[0] < length[0]; ++N[0]) {
//    griditer += pre_offset[1] * total[2];
      for(N[1]=0; N[1] < pre_offset[1]; ++N[1])
        for(N[2]=0; N[2] < total[2]; ++N[2]) {
          ASSERT( griditer != window.Iterators().Local().End(),
              "InterfaceVMGJob::addWindowOntoGrid() - griditer is already at end of window.");
          griditer++;
        }
    for(N[1]=0; N[1] < length[1]; ++N[1]) {
//      griditer += pre_offset[2];
      for(N[2]=0; N[2] < pre_offset[2]; ++N[2]) {
        griditer++;
      }

      for(N[2]=0; N[2] < length[2]; ++N[2]) {
        ASSERT( griditer != window.Iterators().Local().End(),
            "InterfaceVMGJob::addWindowOntoGrid() - griditer is already at end of window.");
        ASSERT( copyiter != grid.sampled_grid.end(),
            "InterfaceVMGJob::addWindowOntoGrid() - griditer is already at end of window.");
        // either smear out charges or not.
        if (DoSmearCharges) {
          smearer(window, griditer++, prefactor*(*copyiter++));
        } else {
          window(*griditer++) += prefactor*(*copyiter++);
        }
      }
//      griditer += post_offset[2];
      for(N[2]=0; N[2] < post_offset[2]; ++N[2]) {
        ASSERT( griditer != window.Iterators().Local().End(),
            "InterfaceVMGJob::addWindowOntoGrid() - griditer is already at end of window.");
        griditer++;
      }
    }
//    griditer += post_offset[1] * total[2];
    for(N[1]=0; N[1] < post_offset[1]; ++N[1])
      for(N[2]=0; N[2] < total[2]; ++N[2]) {
        ASSERT( griditer != window.Iterators().Local().End(),
            "InterfaceVMGJob::addWindowOntoGrid() - griditer is already at end of window.");
        griditer++;
      }
  }
//  griditer += post_offset[0] * total[1] * total[2];
  for(N[0]=0; N[0] < post_offset[0]; ++N[0])
    for(N[1]=0; N[1] < total[1]; ++N[1])
      for(N[2]=0; N[2] < total[2]; ++N[2]) {
        ASSERT( griditer != window.Iterators().Local().End(),
            "InterfaceVMGJob::addWindowOntoGrid() - griditer is already at end of window.");
        griditer++;
      }
#ifndef NDEBUG
  ASSERT( griditer == window.Iterators().Local().End(),
      "InterfaceVMGJob::addWindowOntoGrid() - griditer is not at end of window.");
  ASSERT( copyiter == grid.sampled_grid.end(),
      "InterfaceVMGJob::addWindowOntoGrid() - copyiter is not at end of window.");
#endif

  //  LOG(2, "DEBUG: Grid after adding other is " << grid << ".");
}