Changeset 8f3cdd

Timestamp:

Feb 15, 2013, 5:07:31 PM (12 years ago)

Author:

Frederik Heber <heber@…>

Branches:

Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable

Children:

b123a5

Parents:

64bafe0

git-author:

Frederik Heber <heber@…> (12/19/12 12:48:13)

git-committer:

Frederik Heber <heber@…> (02/15/13 17:07:31)

Message:

Adapted VMGJob and its Interface to new non-zero window use.

• we required two new functions to copy window contents onto the grid and back again.
• SamplingGrid::setWindow() is now private.

Location:

src/Jobs

Files:

• Grid/SamplingGrid.hpp (modified) (2 diffs)
• InterfaceVMGJob.cpp (modified) (7 diffs)
• InterfaceVMGJob.hpp (modified) (5 diffs)
• VMGJob.cpp (modified) (1 diff)

src/Jobs/Grid/SamplingGrid.hpp

@@ -177 +177 @@
   const double getWindowVolume() const;
 
-private:
   /** Sets the size of the window.
    *
@@ -187 +186 @@
   void setWindow(const double _begin_window[3], const double _end_window[3]);
 
+private:
   /** Sets the size of the domain.
    *

src/Jobs/InterfaceVMGJob.cpp

@@ -52 +52 @@
 #include <cmath>
 #include <iostream>
+#include <limits>
 
 #include "CodePatterns/Log.hpp"
 
+#include "Jobs/Grid/SamplingGrid.hpp"
 
 using namespace VMG;
 using VMGInterfaces::InterfaceVMGJob;
 
-InterfaceVMGJob::InterfaceVMGJob(const std::vector< double > &_sampled_input,
+InterfaceVMGJob::InterfaceVMGJob(const SamplingGrid &_sampled_input,
     VMGData &_returndata,
     const std::vector< std::vector<double> > &_particle_positions,
@@ -66 +68 @@
     int levelMin,
     int levelMax,
-    const VMG::Vector &box_begin,
-    vmg_float box_end,
+    const VMG::Vector &_box_begin,
+    vmg_float _box_end,
     const int& near_field_cells,
     int coarseningSteps,
     double alpha) :
   VMG::Interface(boundary, levelMin, levelMax,
-      box_begin, box_end, coarseningSteps, alpha),
+      _box_begin, _box_end, coarseningSteps, alpha),
   spl(near_field_cells, Extent(MaxLevel()).MeshWidth().Max()),
   sampled_input(_sampled_input),
@@ -78 +80 @@
   level(levelMax)
 {
+  for (size_t i=0;i<3;++i) {
+    box_begin[i] = _box_begin[i];
+    box_end[i] = _box_end;
+  }
   std::vector< std::vector<double> >::const_iterator positer = _particle_positions.begin();
   std::vector<double>::const_iterator chargeiter = _particle_charges.begin();
@@ -94 +100 @@
 {}
 
+void InterfaceVMGJob::addGridOntoWindow(
+    VMG::Grid &grid,
+    SamplingGrid &window,
+    const double prefactor) const
+{
+#ifndef NDEBUG
+  for(size_t index=0;index<3;++index) {
+    ASSERT( box_begin[index] >= window.begin_window[index],
+        "InterfaceVMGJob::addGridOntoWindow() - given grid starts earlier than window in component "
+        +toString(index)+".");
+    ASSERT( box_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*(box_begin[index] - window.begin_window[index])+round_offset;
+    length[index] = delta*(box_end[index] - box_begin[index])+round_offset;
+    post_offset[index] = delta*(window.end_window[index] - box_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]);
+  Grid::iterator copyiter = grid.Iterators().Local().Begin();
+  for(N[0]=0; N[0] < length[0]; ++N[0]) {
+    std::advance(griditer, pre_offset[1]*total[2]);
+    for(N[1]=0; N[1] < length[1]; ++N[1]) {
+      std::advance(griditer, pre_offset[2]);
+      for(N[2]=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() - griditer is already at end of window.");
+        *griditer++ += prefactor*grid(*copyiter++);
+      }
+      std::advance(griditer, post_offset[2]);
+    }
+    std::advance(griditer, post_offset[1]*total[2]);
+  }
+#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 InterfaceVMGJob::addWindowOntoGrid(
+    VMG::Grid& window,
+    const SamplingGrid &grid,
+    const double prefactor) const
+{
+#ifndef NDEBUG
+  for(size_t index=0;index<3;++index) {
+    ASSERT( grid.begin_window[index] >= box_begin[index],
+        "InterfaceVMGJob::addWindowOntoGrid() - given window starts earlier than grid in component "
+        +toString(index)+".");
+    ASSERT( grid.end_window[index] <= box_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] - box_begin[index])+round_offset;
+    length[index] = delta*(grid.end_window[index] - grid.begin_window[index])+round_offset;
+    post_offset[index] = delta*(box_end[index] - grid.end_window[index])+round_offset;
+    total[index] = delta*(box_end[index] - box_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];
+  Grid::iterator griditer = window.Iterators().Local().Begin();
+  griditer += pre_offset[0]*total[1]*total[2];
+  SamplingGrid::sampledvalues_t::const_iterator copyiter = grid.sampled_grid.begin();
+  for(N[0]=0; N[0] < length[0]; ++N[0]) {
+    griditer += pre_offset[1]*total[2];
+    for(N[1]=0; N[1] < length[1]; ++N[1]) {
+      griditer += pre_offset[2];
+      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.");
+        window(*griditer++) += prefactor*(*copyiter++);
+      }
+      griditer += post_offset[2];
+    }
+    griditer += post_offset[1]*total[2];
+  }
+#ifndef NDEBUG
+  griditer += post_offset[0]*total[1]*total[2];
+  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 << ".");
+}
+
 void InterfaceVMGJob::ImportRightHandSide(Multigrid& multigrid)
 {
@@ -142 +304 @@
 
   // print nuclei grid to vtk
-//  comm.PrintGrid(particle_grid, "Sampled Nuclei Density");
+  comm.PrintGrid(particle_grid, "Sampled Nuclei Density");
 
   // add sampled electron charge density onto grid
-  std::vector<double>::const_iterator sample_iter = sampled_input.begin();
-  for (Grid::iterator iter = grid.Iterators().Local().Begin();
-      iter != grid.Iterators().Local().End();
-      ++iter)
-    grid(*iter) = -1. * (*sample_iter++);
-  assert( sample_iter == sampled_input.end() );
+  addWindowOntoGrid(
+      grid,
+      sampled_input,
+      -1.);
 
   // print electron grid to vtk
-//  comm.PrintGrid(grid, "Sampled Electron Density");
+  comm.PrintGrid(grid, "Sampled Electron Density");
 
   // add particle_grid onto grid
@@ -180 +340 @@
 
   // print total grid to vtk
-//  comm.PrintGrid(grid, "Total Charge Density");
+  comm.PrintGrid(grid, "Total Charge Density");
 
 //  delete temp_grid;
@@ -200 +360 @@
 
   // obtain sampled potential from grid
-  returndata.sampled_potential.sampled_grid.clear();
-  for (i.X()=grid.Local().Begin().X(); i.X()<grid.Local().End().X(); ++i.X())
-    for (i.Y()=grid.Local().Begin().Y(); i.Y()<grid.Local().End().Y(); ++i.Y())
-      for (i.Z()=grid.Local().Begin().Z(); i.Z()<grid.Local().End().Z(); ++i.Z()) {
-        returndata.sampled_potential.sampled_grid.push_back(grid(i));
-      }
+  returndata.sampled_potential.setWindow(
+      box_begin,
+      box_end
+      );
+  addGridOntoWindow(
+      grid,
+      returndata.sampled_potential,
+      +1.
+      );
 
   // calculate integral over potential as long-range energy contribution

src/Jobs/InterfaceVMGJob.hpp

@@ -22 +22 @@
 #include "Jobs/VMGData.hpp"
 
+class SamplingGrid;
+
 namespace VMG
 {
+class Grid;
 class MGGrid;
 class MGMultigrid;
@@ -34 +37 @@
 {
 public:
-  InterfaceVMGJob(const std::vector< double > &_sampled_input,
+  InterfaceVMGJob(const SamplingGrid &_sampled_input,
       VMGData &returndata,
       const std::vector< std::vector< double > > &_particle_positions,
@@ -41 +44 @@
       int levelMin,
       int levelMax,
-      const VMG::Vector &box_begin,
-      vmg_float box_end,
+      const VMG::Vector &_box_begin,
+      vmg_float _box_end,
       const int& near_field_cells,
       int coarseningSteps=0,
@@ -56 +59 @@
 
 private:
+  /** Helper function to copy VMG's Grid onto the non-zero window of a SamplingGrid
+   *
+   * The code is very similar to SamplingGrid::addOntoWindow().
+   *
+   * @param grid grid with non-zero window
+   * @param prefactor factor for each added value, e.g. -1. makes a subtraction
+   */
+  void addGridOntoWindow(
+      VMG::Grid& grid,
+      SamplingGrid &window,
+      const double prefactor) const;
+
+  /** Helper function to copy non-zero window in SamplingGrid onto VMG's Grid.
+   *
+   * The code is very similar to SamplingGrid::addOntoWindow().
+   *
+   * @param grid grid with non-zero window
+   * @param prefactor factor for each added value, e.g. -1. makes a subtraction
+   */
+  void addWindowOntoGrid(
+      VMG::Grid& window,
+      const SamplingGrid &grid,
+      const double prefactor) const;
+private:
   //!> sampled density on the grid as input
-  const std::vector< double > sampled_input;
+  const SamplingGrid sampled_input;
   //!> sampled potential on the grid as output
   VMGData &returndata;
@@ -64 +91 @@
   //!> nuclei charges
   std::list<VMG::Particle::Particle> particles;
+
+  //!> contains start of grid intervals
+  double box_begin[3];
+  //!> contains end of grid intervals
+  double box_end[3];
 };
 

src/Jobs/VMGJob.cpp

@@ -207 +207 @@
   new DiscretizationPoissonFD(4);
   new VMGInterfaces::InterfaceVMGJob(
-      density_grid.sampled_grid,
+      density_grid,
       returndata,
       particle_positions,