Changeset ab63b6

Timestamp:

Dec 8, 2011, 12:53:48 PM (14 years ago)

Author:

Julian Iseringhausen <isering@…>

Children:

138f86

Parents:

facba0

Message:

Ensure that interpolating splines are normalized wrt the grid.

git-svn-id: ​https://svn.version.fz-juelich.de/scafacos/trunk@1243 5161e1c8-67bf-11de-9fd5-51895aff932f

Location:

src

Files:

• interface/interface_particles.cpp (modified) (1 diff)
• samples/bspline.hpp (modified) (2 diffs)

src/interface/interface_particles.cpp

@@ -100 +100 @@
 #endif
 
-  for (iter=particles.begin(); iter!=particles.end(); ++iter) {
-
-    // Compute global and local grid index of the lower left corner of the grid cell containing the particle
-    index_global = static_cast<Index>((iter->Pos() - temp_grid->Extent().Begin()) / temp_grid->Extent().MeshWidth());
-    index_local = index_global - temp_grid->Global().BeginLocal() + temp_grid->Local().Begin();
-
-    // Iterate over all grid points which lie in the support of the interpolating B-Spline
-    for (index.X()=-1*near_field_cells-1; index.X()<=near_field_cells+1; ++index.X())
-      for (index.Y()=-1*near_field_cells-1; index.Y()<=near_field_cells+1; ++index.Y())
-        for (index.Z()=-1*near_field_cells-1; index.Z()<=near_field_cells+1; ++index.Z()) {
-
-          // Compute distance from grid point to particle
-          length = ( iter->Pos() - grid.Extent().Begin() - grid.Extent().MeshWidth() * (index_global+index) ).Length();
-
-          // If grid point lies in the support of the B-Spline, do charge interpolation
-          if (length < r_cut)
-            (*temp_grid)(index_local+index) += spl.EvaluateSpline(length) * iter->Charge();
-
-        }
-
-  }
+  for (iter=particles.begin(); iter!=particles.end(); ++iter)
+    spl.SetSpline(*temp_grid, *iter, near_field_cells);
 
   // Communicate charges over halo

src/samples/bspline.hpp

@@ -12 +12 @@
 
 #include "base/helper.hpp"
+#include "base/index.hpp"
+#include "base/particle.hpp"
 #include "base/polynomial.hpp"
 #include "base/vector.hpp"
+#include "grid/grid.hpp"
 
 namespace VMG
@@ -29 +32 @@
         return spline_nom[i](val) / spline_denom[i](val);
     return 0.0;
+  }
+
+  void SetSpline(Grid& grid, const Particle::Particle p, const int& near_field_cells)
+  {
+    Index i, index_global, index_local;
+    vmg_float temp_val, length;
+    std::vector<vmg_float> vals;
+    std::vector<vmg_float>::const_iterator iter;
+
+    vmg_float int_val = 0.0;
+
+    index_global = static_cast<Index>((p.Pos() - grid.Extent().Begin()) / grid.Extent().MeshWidth());
+    index_local = index_global - grid.Global().BeginLocal() + grid.Local().Begin();
+
+    // Iterate over all grid points which lie in the support of the interpolating B-Spline
+    for (i.X()=-1*near_field_cells-1; i.X()<=near_field_cells+1; ++i.X())
+      for (i.Y()=-1*near_field_cells-1; i.Y()<=near_field_cells+1; ++i.Y())
+        for (i.Z()=-1*near_field_cells-1; i.Z()<=near_field_cells+1; ++i.Z()) {
+
+          // Compute distance from grid point to particle
+          length = ( p.Pos() - grid.Extent().Begin() - grid.Extent().MeshWidth() * (index_global+i) ).Length();
+
+          temp_val = EvaluateSpline(length);
+          vals.push_back(temp_val * p.Charge());
+          int_val += temp_val;
+
+        }
+
+    int_val = 1.0 / (Helper::pow(grid.MeshWidth(), 3) * int_val);
+
+    iter = vals.begin();
+
+    for (i.X()=-1*near_field_cells-1; i.X()<=near_field_cells+1; ++i.X())
+      for (i.Y()=-1*near_field_cells-1; i.Y()<=near_field_cells+1; ++i.Y())
+        for (i.Z()=-1*near_field_cells-1; i.Z()<=near_field_cells+1; ++i.Z())
+          grid(index_local+i) += *iter++ * int_val;
+
   }