Changeset 6369bc for src/Fragmentation/Summation/SetValues/SamplingGrid.cpp

Timestamp:

Jun 8, 2016, 10:45:28 PM (9 years ago)

Author:

Frederik Heber <heber@…>

Children:

6c0b3f

Parents:

c739b3

git-author:

Frederik Heber <heber@…> (06/08/16 13:22:11)

git-committer:

Frederik Heber <heber@…> (06/08/16 22:45:28)

Message:

SamplingGrid has downsample() function.

• we need this when small fragment grids are calculated at a higher resolution than full_sample grid is supposed to be calculated at.
• extended SamplingGridUnitTest to downsample() function, needed special attention on boundary of grid.

File:

• src/Fragmentation/Summation/SetValues/SamplingGrid.cpp (modified) (3 diffs)

src/Fragmentation/Summation/SetValues/SamplingGrid.cpp

@@ -43 +43 @@
 #include "Fragmentation/Summation/SetValues/SamplingGrid.hpp"
 
+#include <boost/assign.hpp>
 #include <boost/bind.hpp>
 #include <algorithm>
@@ -610 +611 @@
             "SamplingGrid::addWindowOntoWindow() - destiter is already at end of window.");
         ASSERT( sourceiter != source_sampled_grid.end(),
-            "SamplingGrid::addWindowOntoWindow() - destiter is already at end of window.");
+            "SamplingGrid::addWindowOntoWindow() - sourceiter is already at end of window.");
         op(*destiter, *sourceiter);
         ++destiter;
@@ -658 +659 @@
 }
 
+/** Struct contains a single point with displacements from the
+ * central point and the weight in the restriction.
+ */
+struct PointWeight_t {
+  PointWeight_t(const int &d1, const int &d2, const int &d3, const double &_weight) :
+    displacement(NDIM),
+    weight(_weight)
+  {
+    displacement[0] = d1; displacement[1] = d2; displacement[2] = d3;
+  }
+  typedef std::vector<int>  displacement_t;
+  displacement_t displacement;
+  double weight;
+};
+
+static void getLengthsOfGrid(
+    int _total[NDIM],
+    const SamplingGrid &_grid)
+{
+        const size_t gridpoints_axis = _grid.getGridPointsPerAxis();
+        static 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<NDIM; ++index) {
+          if (fabs(_grid.end[index] - _grid.begin[index]) > std::numeric_limits<double>::epsilon()*1e4) {
+            const double delta = (double)gridpoints_axis/(_grid.end[index] - _grid.begin[index]);
+            _total[index] = delta*(_grid.end_window[index] - _grid.begin_window[index])+round_offset;
+          } else
+            _total[index] = 0;
+          ASSERT (_total[index] == ::pow(2, _grid.level),
+              "SamplingGrid::downsample() - total "+toString(_total[index])
+              +" does not match 2^level: "+toString(_grid.level));
+        }
+}
+
+//!> stencil for full weight restriction, see vmg's stencils.hpp
+static const std::vector< PointWeight_t > FullWeightNearestNeighbor =
+    boost::assign::list_of
+    ( PointWeight_t( 0,  0,  0, 0.125) )
+    ( PointWeight_t( 1,  0,  0, 0.0625) )
+    ( PointWeight_t(-1,  0,  0, 0.0625) )
+    ( PointWeight_t( 0,  1,  0, 0.0625) )
+    ( PointWeight_t( 0, -1,  0, 0.0625) )
+    ( PointWeight_t( 0,  0,  1, 0.0625) )
+    ( PointWeight_t( 0,  0, -1, 0.0625) )
+    ( PointWeight_t( 1,  1,  0, 0.03125) )
+    ( PointWeight_t( 1, -1,  0, 0.03125) )
+    ( PointWeight_t(-1,  1,  0, 0.03125) )
+    ( PointWeight_t(-1, -1,  0, 0.03125) )
+    ( PointWeight_t( 0,  1,  1, 0.03125) )
+    ( PointWeight_t( 0,  1, -1, 0.03125) )
+    ( PointWeight_t( 0, -1,  1, 0.03125) )
+    ( PointWeight_t( 0, -1, -1, 0.03125) )
+    ( PointWeight_t( 1,  0,  1, 0.03125) )
+    ( PointWeight_t( 1,  0, -1, 0.03125) )
+    ( PointWeight_t(-1,  0,  1, 0.03125) )
+    ( PointWeight_t(-1,  0, -1, 0.03125) )
+    ( PointWeight_t( 1,  1,  1, 0.015625) )
+    ( PointWeight_t( 1,  1, -1, 0.015625) )
+    ( PointWeight_t( 1, -1,  1, 0.015625) )
+    ( PointWeight_t(-1,  1,  1, 0.015625) )
+    ( PointWeight_t( 1, -1, -1, 0.015625) )
+    ( PointWeight_t(-1,  1, -1, 0.015625) )
+    ( PointWeight_t(-1, -1,  1, 0.015625) )
+    ( PointWeight_t(-1, -1, -1, 0.015625) )
+;
+
+int getValidIndex(
+    const PointWeight_t::displacement_t &_disp,
+    const int N[NDIM],
+    const int length[NDIM])
+{
+  int index = 0;
+  // we simply truncate in case of out of bounds access
+  if ((N[2]+_disp[2] >= 0) && (N[2]+_disp[2] < length[2]))
+    index += _disp[2];
+  if ((N[1]+_disp[1] >= 0) && (N[1]+_disp[1] < length[1]))
+    index += _disp[1]*length[2];
+  if ((N[0]+_disp[0] >= 0) && (N[0]+_disp[0] < length[0]))
+    index += _disp[0]*length[1]*length[2];
+  return index;
+}
+
+void restrictFullWeight(
+    SamplingGrid::sampledvalues_t &_coarse_level,
+    const int length_c[NDIM],
+    const SamplingGrid::sampledvalues_t &_fine_level,
+    const int length_f[NDIM])
+{
+  int N_c[NDIM];
+  int N_f[NDIM];
+  SamplingGrid::sampledvalues_t::iterator coarseiter = _coarse_level.begin();
+  for(N_c[0]=0, N_f[0]=0; (N_c[0] < length_c[0]) && (N_f[0] < length_f[0]); ++N_c[0], N_f[0] +=2) {
+    for(N_c[1]=0, N_f[1]=0; (N_c[1] < length_c[1]) && (N_f[1] < length_f[1]); ++N_c[1], N_f[1] +=2) {
+      for(N_c[2]=0, N_f[2]=0; (N_c[2] < length_c[2]) && (N_f[2] < length_f[2]); ++N_c[2], N_f[2] +=2) {
+        const int index_base = N_f[2] + (N_f[1] + N_f[0]*length_f[1])*length_f[2];
+        // go through stencil and add each point relative to displacement with weight
+        for (std::vector< PointWeight_t >::const_iterator weightiter = FullWeightNearestNeighbor.begin();
+            weightiter != FullWeightNearestNeighbor.end(); ++weightiter) {
+          const PointWeight_t::displacement_t disp = weightiter->displacement;
+          const int index_disp = getValidIndex(disp, N_f, length_f);
+          *coarseiter += _fine_level[index_base+index_disp]*weightiter->weight;
+        }
+        ++coarseiter;
+      }
+      ASSERT ( (N_c[2] == length_c[2]) && (N_f[2] == length_f[2]),
+          "restrictFullWeight() - N_c "+toString(N_c[2])+" != length_c "+toString(length_c[2])
+          +" or N_f "+toString(N_f[2])+" != length_f "+toString(length_f[2]));
+    }
+    ASSERT ( (N_c[1] == length_c[1]) && (N_f[1] == length_f[1]),
+        "restrictFullWeight() - N_c "+toString(N_c[1])+" != length_c "+toString(length_c[1])
+        +" or N_f "+toString(N_f[1])+" != length_f "+toString(length_f[1]));
+  }
+  ASSERT ( (N_c[0] == length_c[0]) && (N_f[0] == length_f[0]),
+      "restrictFullWeight() - N_c "+toString(N_c[0])+" != length_c "+toString(length_c[0])
+      +" or N_f "+toString(N_f[0])+" != length_f "+toString(length_f[0]));
+  ASSERT( coarseiter ==  _coarse_level.end(),
+      "restrictFullWeight() - coarseiter is not at end of values.");
+}
+
+void SamplingGrid::downsample(
+    SamplingGrid& instance,
+    const SamplingGrid& other,
+    const int _level)
+{
+  if (&instance != &other) {
+    // take over properties
+    static_cast<SamplingGridProperties &>(instance) = other;
+    instance.setWindowSize(other.begin_window, other.end_window);
+    if (_level >= other.level) {
+      instance.sampled_grid = other.sampled_grid;
+    } else {
+      // if desired level is smaller we need to downsample
+      // we do this similarly to vmg::RestrictionFullWeight (i.e. a full nearest
+      // neighbor interpolation) and always one grid level at a time till we
+      // have reached the desired one
+
+      // the reference such that we never have to copy the full grid but only
+      // downsampled ones
+      const sampledvalues_t * sourcevalues = &other.sampled_grid;
+      int length_d[3];
+      int length_s[3];
+      getLengthsOfGrid(length_s, other);
+      for (instance.level = other.level-1; instance.level >= _level; --instance.level) {
+        getLengthsOfGrid(length_d, instance);
+        // we always have an eighth of the number of sample points as we stop
+        sampledvalues_t downsampled(sourcevalues->size()/(size_t)8, 0.);
+        restrictFullWeight(downsampled, length_d, *sourcevalues, length_s);
+        // then copy the downsampled values
+        instance.sampled_grid = downsampled;
+        sourcevalues = &instance.sampled_grid;
+        // and exchange lengths
+        for (size_t i=0;i<3;++i) {
+          length_s[i] = length_d[i];
+        }
+      }
+      instance.level = _level;
+
+      // and finally, renormalize downsampled grid to old value
+//      instance *= other.integral()/instance.integral();
+    }
+  }
+}
+
 std::ostream & operator<<(std::ostream &ost, const SamplingGrid& other)
 {