Changeset b303d0

Timestamp:

Sep 6, 2015, 10:25:02 PM (10 years ago)

Author:

Frederik Heber <heber@…>

Children:

3f6604

Parents:

ebc866

git-author:

Frederik Heber <heber@…> (03/13/15 15:02:44)

git-committer:

Frederik Heber <heber@…> (09/06/15 22:25:02)

Message:

tempcommit: Fixes to SubtractSelfInducedForces().

File:

• src/units/particle/bspline.hpp (modified) (5 diffs)

src/units/particle/bspline.hpp

@@ -210 +210 @@
         const vmg_float& h_z = grid.Extent().MeshWidth().Z();
 
-        int int_val = 0.;
+        vmg_float int_val = 0.;
     vmg_float dir_x = pos_beg_x;
     for (int i=-1*near_field_cells; i<=near_field_cells; ++i) {
@@ -218 +218 @@
         for (int k=-1*near_field_cells; k<=near_field_cells; ++k) {
 
-          // Compute distance from grid point to particle
-          int_val += EvaluateSpline(std::sqrt(dir_x*dir_x+dir_y*dir_y+dir_z*dir_z));
+                  const double length_sq = dir_x*dir_x+dir_y*dir_y+dir_z*dir_z;
+                  if (fabs(length_sq) > std::numeric_limits<double>::epsilon()) {
+                          // Compute distance from grid point to particle
+                          int_val += EvaluateSpline(std::sqrt(length_sq));
+                  }
 
           dir_z += h_z;
@@ -229 +232 @@
     // Reciprocal value of the numerically integrated spline
     if (fabs(int_val) > std::numeric_limits<double>::epsilon())
-        int_val = 1.0 / (int_val * h_x * h_y * h_z);
+        int_val = 1. / (int_val * h_x * h_y * h_z);
     else
-        int_val = 1.0;
+        int_val = 1. / (h_x * h_y * h_z);
 
     // Iterate over all grid points which lie in the support of the interpolating B-Spline
+    vmg_float test_int_val = 0.;
         dir_x = pos_beg_x;
         for (int i=-1*near_field_cells; i<=near_field_cells; ++i) {
@@ -245 +249 @@
                   if (fabs(length_sq) > std::numeric_limits<double>::epsilon()) {
                     const double length = std::sqrt(length_sq);
-//                  temp_val = EvaluateField(length);
-//                  p.Field()[0] -= p.Charge() * dir_x * temp_val;
-//                  p.Field()[1] -= p.Charge() * dir_y * temp_val;
-//                  p.Field()[2] -= p.Charge() * dir_z * temp_val;
-                    temp_val = int_val * EvaluateSpline(length);
-                    p.Field()[0] -= p.Charge()*p.Charge() * temp_val * temp_val * dir_x / (length_sq*length);
-                    p.Field()[1] -= p.Charge()*p.Charge() * temp_val * temp_val * dir_y / (length_sq*length);
-                    p.Field()[2] -= p.Charge()*p.Charge() * temp_val * temp_val * dir_z / (length_sq*length);
+
+                    temp_val = h_x * h_y * h_z * int_val * EvaluateSpline(length);
+                    test_int_val += temp_val;
+                    p.Field()[0] -= p.Charge()* temp_val * dir_x / (length_sq*length);
+                    p.Field()[1] -= p.Charge()* temp_val * dir_y / (length_sq*length);
+                    p.Field()[2] -= p.Charge()* temp_val * dir_z / (length_sq*length);
                   } else {
                           std::cerr << "Value very small " << length_sq << "=("
@@ -267 +269 @@
           dir_x += h_x;
         }
+        if ( fabs(test_int_val -1.) > std::numeric_limits<double>::epsilon()*1e4 ) {
+                std::cerr << "Integrated spline value should be 1 but is " << test_int_val << std::endl;
+                assert(0);
+        }
   }