Context Navigation

-              r6a5eb2
+              r57fda0
 #include <boost/assign.hpp>
 #include <boost/bind.hpp>
+#include <boost/shared_ptr.hpp>
 #include <algorithm>
 #include <limits>
 …
+}
+static void getDownsampledGrid(
+    const SamplingGrid &_reference_grid,
+    const SamplingGrid &_grid,
+    boost::shared_ptr<SamplingGrid> &_weight_downsampled)
+{
+  static const double round_offset(
+      (std::numeric_limits<size_t>::round_style == std::round_toward_zero) ?
+.5 : 0.); // need offset to get to round_toward_nearest behavior
+  const int surplus_level = _reference_grid.getSurplusLevel(_grid)+round_offset;
+  // need to downsample first
+  _weight_downsampled.reset(new SamplingGrid); // may use default cstor as downsample does all settings
+  SamplingGrid::downsample(*_weight_downsampled, _grid, _reference_grid.level-surplus_level);
+}
 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<NDIM;++index) {
+      if (begin_window[index] <= other.begin_window[index]) {
+        min_begin_window[index] = other.begin_window[index];
+        doShrink = true;
+      } else {
+        min_begin_window[index] = begin_window[index];
+      }
+      if (end_window[index] <= other.end_window[index]) {
+        min_end_window[index] = end_window[index];
+      } else {
+        min_end_window[index] = other.end_window[index];
+        doShrink = true;
+      }
+    }
+    LOG(4, "DEBUG: min begin is " << min_begin_window[0] << "," << min_begin_window[1] << "," << min_begin_window[2] << ".");
+    LOG(4, "DEBUG: min end is " << min_end_window[0] << "," << min_end_window[1] << "," << min_end_window[2] << ".");
+    if (doShrink)
+      shrinkWindow(min_begin_window, min_end_window);
+    addWindowOntoWindow(
+        other.begin_window,
+        other.end_window,
+        begin_window,
+        end_window,
+        sampled_grid,
+        other.sampled_grid,
+        boost::bind(multiplyElements, _1, _2, 1.),
+        sourcewindow);
+  } else {
+    ASSERT(0, "SamplingGrid::operator*=() - multiplying incongruent grids is so far not in the cards.");
+  }
+  ASSERT(isCompatible(other),
+      "SamplingGrid::operator*=() - multiplying incomatible grids is so far not in the cards.");
+  const SamplingGrid *other_grid = &other;
+  boost::shared_ptr<SamplingGrid> other_downsampled;
+  if (!isEquivalent(other)) {
+    getDownsampledGrid(*this, other, other_downsampled);
+    other_grid = other_downsampled.get();
+  }
+  /// get minimum of window
+  double min_begin_window[NDIM];
+  double min_end_window[NDIM];
+  bool doShrink = false;
+  for (size_t index=0; index<NDIM;++index) {
+    if (begin_window[index] <= other.begin_window[index]) {
+      min_begin_window[index] = other.begin_window[index];
+      doShrink = true;
+    } else {
+      min_begin_window[index] = begin_window[index];
+    }
+    if (end_window[index] <= other.end_window[index]) {
+      min_end_window[index] = end_window[index];
+    } else {
+      min_end_window[index] = other.end_window[index];
+      doShrink = true;
+    }
+  }
+  LOG(4, "DEBUG: min begin is " << min_begin_window[0] << "," << min_begin_window[1] << "," << min_begin_window[2] << ".");
+  LOG(4, "DEBUG: min end is " << min_end_window[0] << "," << min_end_window[1] << "," << min_end_window[2] << ".");
+  if (doShrink)
+    shrinkWindow(min_begin_window, min_end_window);
+  addWindowOntoWindow(
+      other_grid->begin_window,
+      other_grid->end_window,
+      begin_window,
+      end_window,
+      sampled_grid,
+      other_grid->sampled_grid,
+      boost::bind(multiplyElements, _1, _2, 1.),
+      sourcewindow);
   return *this;
+}
 …
 void SamplingGrid::superposeOtherGrids(const SamplingGrid &other, const double prefactor)
+{
+  /// check that grids are compatible
+  if (isEquivalent(other)) {
+    /// get maximum of window
+    double max_begin_window[NDIM];
+    double max_end_window[NDIM];
+    bool doExtend = false;
+    for (size_t index=0; index<NDIM;++index) {
+      if (begin_window[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.");
+  }
+  // check that grids are compatible
+  ASSERT(isCompatible(other),
+      "SamplingGrid::superposeOtherGrids() - superposing incompatible grids is so far not in the cards.");
+  const SamplingGrid *other_grid = &other;
+  boost::shared_ptr<SamplingGrid> other_downsampled;
+  if (!isEquivalent(other)) {
+    getDownsampledGrid(*this, other, other_downsampled);
+    other_grid = other_downsampled.get();
+  }
+  /// get maximum of window
+  double max_begin_window[NDIM];
+  double max_end_window[NDIM];
+  bool doExtend = false;
+  for (size_t index=0; index<NDIM;++index) {
+    if (begin_window[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_grid->begin_window, other_grid->end_window, other_grid->sampled_grid, prefactor);
+}
 …
 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.;
+}
+  // check that grids are compatible
+  ASSERT(isCompatible(weight),
+      "SamplingGrid::integral() - integrating with weights from incompatible grids is so far not in the cards.");
+  const SamplingGrid *weight_grid = &weight;
+  boost::shared_ptr<SamplingGrid> weight_downsampled;
+  if (!isEquivalent(weight)) {
+    getDownsampledGrid(*this, weight, weight_downsampled);
+    weight_grid = weight_downsampled.get();
+  }
+  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_grid->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;
+}
 void SamplingGrid::setWindowSize(
     const double _begin_window[NDIM],

Note: See TracChangeset for help on using the changeset viewer.

Context Navigation

Changeset 57fda0 for src/Fragmentation/Summation

Legend:

src/Fragmentation/Summation/SetValues/SamplingGrid.cpp

Download in other formats: