Index: src/Fragmentation/SetValues/SamplingGrid.hpp
===================================================================
--- src/Fragmentation/SetValues/SamplingGrid.hpp	(revision 7b049409753f488e726192ddd4e7610aed7e12b5)
+++ src/Fragmentation/SetValues/SamplingGrid.hpp	(revision 7b049409753f488e726192ddd4e7610aed7e12b5)
@@ -0,0 +1,416 @@
+/*
+ * SamplingGrid.hpp
+ *
+ *  Created on: 25.07.2012
+ *      Author: heber
+ */
+
+#ifndef SAMPLINGGRID_HPP_
+#define SAMPLINGGRID_HPP_
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <boost/function.hpp>
+#include <iosfwd>
+#include <vector>
+
+#include "boost/serialization/export.hpp"
+#include "boost/serialization/vector.hpp"
+
+#include "Fragmentation/SetValues/SamplingGridProperties.hpp"
+
+class MPQCData;
+class SamplingGridTest;
+
+/** This class stores a sample function on a three-dimensional grid.
+ *
+ *  \note We do not use boost::multi_array because it is not trivial to serialize.
+ *
+ */
+class SamplingGrid : public SamplingGridProperties
+{
+  //!> grant unit test access to private parts
+  friend class SamplingGridTest;
+  //!> grant output operator access
+  friend std::ostream & operator<<(std::ostream &ost, const SamplingGrid& other);
+public:
+  //!> typedef for sampled values
+  typedef std::vector< double > sampledvalues_t;
+
+  /** Constructor for class SamplingGrid for full window.
+   *
+   * Here, the window of sampled values spans the given domain.
+   *
+   * \param _begin offset for grid per axis
+   * \param _end edge length of grid per axis
+   * \param _level number of grid points in \f$2^{\mathrm{level}}\f$
+   * \param _sampled_grid sample points
+   */
+  SamplingGrid(const double _begin[3],
+      const double _end[3],
+      const int _level,
+      const sampledvalues_t &_sampled_grid);
+
+  /** Constructor for class SamplingGrid for empty window.
+   *
+   * Here, the window is initially of size zero.
+   *
+   * \param _begin offset for grid per axis
+   * \param _end edge length of grid per axis
+   * \param _level number of grid points in \f$2^{\mathrm{level}}\f$
+   */
+  SamplingGrid(const double _begin[3],
+      const double _end[3],
+      const int _level);
+
+  /** Copy constructor for class SamplingGrid with full window from SamplingGridProperties.
+   *
+   * Here, the window is initially empty.
+   *
+   * \param _props properties to copy
+   */
+  SamplingGrid(const SamplingGridProperties &_props);
+
+  /** Copy constructor for class SamplingGrid with empty window from SamplingGridProperties.
+   *
+   * Here, the window must span the whole domain
+   *
+   * \param _props properties to copy
+   * \param _sampled_grid sample points
+   */
+  SamplingGrid(
+      const SamplingGridProperties &_props,
+      const sampledvalues_t &_sampled_grid);
+
+  /** Copy constructor for class SamplingGrid.
+   *
+   * The window of sampled values corresponds to the one on \a _grid.
+   *
+   * \param _grid grid to copy
+   */
+  SamplingGrid(const SamplingGrid &_grid);
+
+  /** default cstor.
+   */
+  SamplingGrid();
+
+  virtual ~SamplingGrid();
+
+  /** Checks whether another instance is consistent with this one.
+   *
+   *  \note Conistency is stronger as grids must have the same window.
+   *
+   * \param _props other properties to check against
+   * \return true - are consistent, false - else
+   */
+  bool isCongruent(const SamplingGrid &_props) const;
+
+  /** Assignment operator.
+   *
+   * \param other other instance to assign ourselves to
+   */
+  SamplingGrid& operator=(const SamplingGrid& other);
+
+  /** Addition operator with another SamplingGrid instance \a other.
+   *
+   * \param other other instance to sum onto this one.
+   * \return ref to this instance
+   */
+  SamplingGrid& operator+=(const SamplingGrid& other)
+  {
+    superposeOtherGrids(other, +1.);
+    return *this;
+  }
+
+  /** Element-wise multiplication operator with another SamplingGrid instance \a other.
+   *
+   * With non-zero windows we have to pay some more attention here.
+   * Now, the windows may not be congruent but we have to find the intersection
+   * of the two windows and then construct the new window only of this size and
+   * multiply. The trick then is not to copy&change the other grid but to
+   * access it properly.
+   *
+   * \param other other instance to sum onto this one.
+   * \return ref to this instance
+   */
+  SamplingGrid& operator*=(const SamplingGrid& other);
+
+  /** Subtraction operator with another SamplingGrid instance \a other.
+   *
+   * \param other other instance to subtract from this one.
+   * \return ref to this instance
+   */
+  SamplingGrid& operator-=(const SamplingGrid& other)
+  {
+    superposeOtherGrids(other, -1.);
+    return *this;
+  }
+
+  /** Returns the numeric integral over the grid.
+   *
+   * @return sum of grid values times volume element
+   */
+  double integral() const;
+
+  /** Returns the numeric integral over the grid where the grid is element-wise multiplied with \a weight.
+   *
+   * @param weight grid of weights
+   * @return sum of grid values weighted by respective element in weight times volume element
+   */
+  double integral(const SamplingGrid &weight) const;
+
+  /** Returns the volume of the domain for this sampled function.
+   *
+   * @return volume
+   */
+  const double getVolume() const;
+
+  /** Returns the total number of gridpoints of the discrete mesh covering the (window) volume.
+   *
+   * @return number of gridpoints sampled_values should have
+   */
+  const size_t getWindowGridPoints() const;
+
+  /** Returns the number of gridpoints of the discrete mesh for the current
+   * window size for given axis \axis.
+   *
+   * \param axis axis to calculate number of gridpoints for
+   * \return number of gridpoints along this axis
+   */
+  const size_t getWindowGridPointsPerAxis(const size_t axis) const;
+
+  /** Returns the total number of gridpoints of the discrete mesh covering the volume.
+   *
+   * @return number of gridpoints sampled_values should have
+   */
+  const size_t getTotalGridPoints() const;
+
+  /** Returns the number of grid points per axis.
+   *
+   * @return number of grid points per unit length
+   */
+  const size_t getGridPointsPerAxis() const;
+
+  /** Returns the length of the domain for the given \a axis.
+   *
+   * \param axis axis for which to get step length
+   * \return domain length for the given axis, i.e. end - begin
+   */
+  const double getTotalLengthPerAxis(const size_t axis) const;
+
+  /** Returns the length of the window for the given \a axis.
+   *
+   * \param axis axis for which to get step length
+   * \return window length for the given axis, i.e. end - begin
+   */
+  const double getWindowLengthPerAxis(const size_t axis) const;
+
+  /** Returns the real step length from one discrete grid point to the next.
+   *
+   * \param axis axis for which to get step length
+   * \return step length for the given axis, as domain length over getGridPointsPerAxis()
+   */
+  const double getDeltaPerAxis(const size_t axis) const;
+
+  /** Returns the volume of the domain covered by the current window.
+   *
+   * @return volume of window
+   */
+  const double getWindowVolume() const;
+
+  /** Sets the size of the window.
+   *
+   * \note also resets the sampled points so far.
+   *
+   * \param _begin_window start of new window
+   * \param _end_window end of window
+   */
+  void setWindow(const double _begin_window[3], const double _end_window[3]);
+
+private:
+  /** Sets the size of the domain.
+   *
+   * \note also resets the sampled points so far and the window.
+   *
+   * \param _begin start of new window
+   * \param _end end of window
+   */
+  void setDomain(const double _begin[3], const double _end[3]);
+
+  /** Sets the size of the domain.
+   *
+   * \note this is just internally used for easing the array setting.
+   *
+   * \param _begin start of domain
+   * \param _end end of domain
+   */
+  void setDomainSize(const double _begin[3], const double _end[3]);
+
+  /** Extends the window while keeping the values.
+   *
+   * \param _begin_window new start of window
+   * \param _end_window new end of window
+   */
+  void extendWindow(const double _begin_window[3], const double _end_window[3]);
+
+  /** Shrinks the window while keeping the values.
+   *
+   * \param _begin_window new start of window
+   * \param _end_window new end of window
+   */
+  void shrinkWindow(const double _begin_window[3], const double _end_window[3]);
+
+  /** Adds another (smaller) window onto the one in this instance.
+   *
+   * \note We assume here that the given window fits on the this one.
+   *
+   * \param _begin_window start of other window
+   * \param _end_window end of other window
+   * \param _sampled_grid other set of sampled values
+   * @param prefactor +1. is then addition, -1. is subtraction.
+   */
+  void addOntoWindow(
+      const double _begin_window[3],
+      const double _end_window[3],
+      const sampledvalues_t &_sampled_grid,
+      const double prefactor);
+
+  /** Adds another (larger) window into the one in this instance.
+   *
+   * \note We assume here that the given window is larger than this one.
+   *
+   * \param _begin_window start of other window
+   * \param _end_window end of other window
+   * \param _sampled_grid other set of sampled values
+   * @param prefactor +1. is then addition, -1. is subtraction.
+   */
+  void addIntoWindow(
+      const double _begin_window[3],
+      const double _end_window[3],
+      const sampledvalues_t &_sampled_grid,
+      const double prefactor);
+
+  /** Enum to help in addWindowOntoWindow() decide which iterator needs to be
+   * advanced.
+   */
+  enum eLargerWindow {
+    destwindow,
+    sourcewindow
+  };
+
+  /** Helper function to copy one (larger) window into a (smaller) window.
+   *
+   * \note Why do we need the extra \a choice? We need to know which window
+   * tuples is associated with which sampled values that are constrained by
+   * one of them being constant, hence the source values
+   *
+   * \param larger_wbegin start of larger window
+   * \param larger_wend end of larger window
+   * \param smaller_wbegin start of smaller window
+   * \param smaller_wend end of smaller window
+   * \param dest_sampled_grid larger set of sampled values
+   * \param source_sampled_grid smaller set of sampled values
+   * \param op operation to perform with the two elements
+   * \param larger_window indicates which is the larger window
+   */
+  void addWindowOntoWindow(
+      const double larger_wbegin[3],
+      const double larger_wend[3],
+      const double smaller_wbegin[3],
+      const double smaller_wend[3],
+      sampledvalues_t &dest_sampled_grid,
+      const sampledvalues_t &source_sampled_grid,
+      boost::function<void (double &, const double &)> op,
+      enum eLargerWindow larger_window);
+
+  /** Helper function that contains all the logic of how to superpose two
+   * grids.
+   *
+   * Is called by SamplingGrid::operator+=() and SamplingGrid::operator-=()
+   *
+   * @param other other histogram
+   * @param prefactor +1. is then addition, -1. is subtraction.
+   */
+  void superposeOtherGrids(const SamplingGrid &other, const double prefactor);
+
+  /** Sets the size of the window.
+   *
+   * \note also resets the sampled points so far.
+   *
+   * \param _begin_window start of new window
+   * \param _end_window end of window
+   */
+  void setWindowSize(const double _begin_window[3], const double _end_window[3]);
+
+  /** Helper function to get point at grid point for given \a axis and less than value.
+   *
+   * @param value value to find nearest grid point to
+   * @param axis axis of the value
+   * @return nearest lower grid point
+   */
+  double getNearestLowerGridPoint(
+      const double value, const size_t axis) const;
+
+  /** Helper function to get point at grid point for given \a axis and greater than value.
+   *
+   * @param value value to find nearest grid point to
+   * @param axis axis of the value
+   * @return nearest lower grid point
+   */
+  double getNearestHigherGridPoint(
+      const double value, const size_t axis) const;
+
+public:
+  /// We do not store the whole grid if many entries are actually zero
+  /// but only a window wherein the sampled function is non-zero.
+
+  //!> sample points of the window
+  sampledvalues_t sampled_grid;
+
+  //!> start of the window relative to SamplingGridProperties::begin and SamplingGridProperties::size
+  double begin_window[3];
+  //!> end of the window relative to SamplingGridProperties::begin and SamplingGridProperties::size
+  double end_window[3];
+
+private:
+  friend class MPQCData;
+
+  friend class boost::serialization::access;
+  // serialization
+  template <typename Archive>
+  void serialize(Archive& ar, const unsigned int version)
+  {
+    ar & boost::serialization::base_object<SamplingGridProperties>(*this);
+    ar & const_cast< sampledvalues_t &>(sampled_grid);
+    for(size_t i=0;i<3;++i) {
+      ar & begin_window[i];
+      ar & end_window[i];
+    }
+  }
+
+  //!> static typedef to use in cstor when no initial values are given
+  static const double zeroOffset[3];
+};
+
+/** Output operator for class SamplingGrid.
+ *
+ * \param ost output stream to print to
+ * \param other instance to print
+ * \return ref to stream for concatenation
+ */
+std::ostream & operator<<(std::ostream &ost, const SamplingGrid& other);
+
+template<typename T> T ZeroInstance();
+template<> SamplingGrid ZeroInstance<SamplingGrid>();
+
+// we need to give this class a unique key for serialization
+// its is only serialized through its base class FragmentJob
+BOOST_CLASS_EXPORT_KEY(SamplingGrid)
+
+// define inline functions
+#include "SamplingGrid_inline.hpp"
+
+#endif /* SAMPLINGGRID_HPP_ */
