Changeset b4e1be

Timestamp:

Mar 9, 2015, 4:49:20 PM (11 years ago)

Author:

Frederik Heber <heber@…>

Children:

5c22be

Parents:

fb3b4c

Message:

tempcommit: Added functions to bspline.

File:

• src/units/particle/bspline.hpp (modified) (1 diff)

src/units/particle/bspline.hpp

@@ -128 +128 @@
   }
 
+  void changeGridBySelfInducedPotential(Grid& grid, Particle& p, const vmg_float &sign) const
+  {
+        assert(p.Pos().X() >= grid.Extent().Begin().X() && p.Pos().X() < grid.Extent().End().X());
+        assert(p.Pos().Y() >= grid.Extent().Begin().Y() && p.Pos().Y() < grid.Extent().End().Y());
+        assert(p.Pos().Z() >= grid.Extent().Begin().Z() && p.Pos().Z() < grid.Extent().End().Z());
+
+        const int index_global_x = grid.Global().GlobalBegin().X() + std::floor((p.Pos().X() - grid.Extent().Begin().X()) / grid.Extent().MeshWidth().X());
+        const int index_global_y = grid.Global().GlobalBegin().Y() + std::floor((p.Pos().Y() - grid.Extent().Begin().Y()) / grid.Extent().MeshWidth().Y());
+        const int index_global_z = grid.Global().GlobalBegin().Z() + std::floor((p.Pos().Z() - grid.Extent().Begin().Z()) / grid.Extent().MeshWidth().Z());
+
+        assert(index_global_x >= grid.Global().LocalBegin().X() && index_global_x < grid.Global().LocalEnd().X());
+        assert(index_global_y >= grid.Global().LocalBegin().Y() && index_global_y < grid.Global().LocalEnd().Y());
+        assert(index_global_z >= grid.Global().LocalBegin().Z() && index_global_z < grid.Global().LocalEnd().Z());
+
+        const int index_local_x = index_global_x - grid.Global().LocalBegin().X() + grid.Local().HaloSize1().X();
+        const int index_local_y = index_global_y - grid.Global().LocalBegin().Y() + grid.Local().HaloSize1().Y();
+        const int index_local_z = index_global_z - grid.Global().LocalBegin().Z() + grid.Local().HaloSize1().Z();
+
+        assert(index_local_x >= grid.Local().Begin().X() && index_local_x < grid.Local().End().X());
+        assert(index_local_y >= grid.Local().Begin().Y() && index_local_y < grid.Local().End().Y());
+        assert(index_local_z >= grid.Local().Begin().Z() && index_local_z < grid.Local().End().Z());
+
+        const vmg_float pos_beg_x =  grid.Extent().Begin().X() + grid.Extent().MeshWidth().X() * (index_global_x - grid.Global().GlobalBegin().X() - near_field_cells) - p.Pos().X();
+        const vmg_float pos_beg_y =  grid.Extent().Begin().Y() + grid.Extent().MeshWidth().Y() * (index_global_y - grid.Global().GlobalBegin().Y() - near_field_cells) - p.Pos().Y();
+        const vmg_float pos_beg_z =  grid.Extent().Begin().Z() + grid.Extent().MeshWidth().Z() * (index_global_z - grid.Global().GlobalBegin().Z() - near_field_cells) - p.Pos().Z();
+
+        const vmg_float& h_x = grid.Extent().MeshWidth().X();
+        const vmg_float& h_y = grid.Extent().MeshWidth().Y();
+        const vmg_float& h_z = grid.Extent().MeshWidth().Z();
+
+        vmg_float dir_x = pos_beg_x;
+    for (int i=-1*near_field_cells; i<=near_field_cells; ++i) {
+          vmg_float dir_y = pos_beg_y;
+      for (int j=-1*near_field_cells; j<=near_field_cells; ++j) {
+                vmg_float dir_z = pos_beg_z;
+        for (int k=-1*near_field_cells; k<=near_field_cells; ++k) {
+                  const double length_sq = dir_x*dir_x+dir_y*dir_y+dir_z*dir_z;
+                  const double length = std::sqrt(length_sq);
+          grid(index_local_x + i,
+              index_local_y + j,
+              index_local_z + k) += sign * p.Charge() * EvaluatePotential(length);
+
+          dir_z += h_z;
+                }
+            dir_y += h_y;
+          }
+          dir_x += h_x;
+        }
+  }
+
+  void SubtractSelfInducedForces(const Grid& grid, Particle& p) const
+  {
+        assert(p.Pos().X() >= grid.Extent().Begin().X() && p.Pos().X() < grid.Extent().End().X());
+        assert(p.Pos().Y() >= grid.Extent().Begin().Y() && p.Pos().Y() < grid.Extent().End().Y());
+        assert(p.Pos().Z() >= grid.Extent().Begin().Z() && p.Pos().Z() < grid.Extent().End().Z());
+
+        vmg_float temp_val = 0.;
+
+        const int index_global_x = grid.Global().GlobalBegin().X() + std::floor((p.Pos().X() - grid.Extent().Begin().X()) / grid.Extent().MeshWidth().X());
+        const int index_global_y = grid.Global().GlobalBegin().Y() + std::floor((p.Pos().Y() - grid.Extent().Begin().Y()) / grid.Extent().MeshWidth().Y());
+        const int index_global_z = grid.Global().GlobalBegin().Z() + std::floor((p.Pos().Z() - grid.Extent().Begin().Z()) / grid.Extent().MeshWidth().Z());
+
+        assert(index_global_x >= grid.Global().LocalBegin().X() && index_global_x < grid.Global().LocalEnd().X());
+        assert(index_global_y >= grid.Global().LocalBegin().Y() && index_global_y < grid.Global().LocalEnd().Y());
+        assert(index_global_z >= grid.Global().LocalBegin().Z() && index_global_z < grid.Global().LocalEnd().Z());
+
+        const int index_local_x = index_global_x - grid.Global().LocalBegin().X() + grid.Local().HaloSize1().X();
+        const int index_local_y = index_global_y - grid.Global().LocalBegin().Y() + grid.Local().HaloSize1().Y();
+        const int index_local_z = index_global_z - grid.Global().LocalBegin().Z() + grid.Local().HaloSize1().Z();
+
+        assert(index_local_x >= grid.Local().Begin().X() && index_local_x < grid.Local().End().X());
+        assert(index_local_y >= grid.Local().Begin().Y() && index_local_y < grid.Local().End().Y());
+        assert(index_local_z >= grid.Local().Begin().Z() && index_local_z < grid.Local().End().Z());
+
+        const vmg_float pos_beg_x =  grid.Extent().Begin().X() + grid.Extent().MeshWidth().X() * (index_global_x - grid.Global().GlobalBegin().X() - near_field_cells) - p.Pos().X();
+        const vmg_float pos_beg_y =  grid.Extent().Begin().Y() + grid.Extent().MeshWidth().Y() * (index_global_y - grid.Global().GlobalBegin().Y() - near_field_cells) - p.Pos().Y();
+        const vmg_float pos_beg_z =  grid.Extent().Begin().Z() + grid.Extent().MeshWidth().Z() * (index_global_z - grid.Global().GlobalBegin().Z() - near_field_cells) - p.Pos().Z();
+
+        const vmg_float& h_x = grid.Extent().MeshWidth().X();
+        const vmg_float& h_y = grid.Extent().MeshWidth().Y();
+        const vmg_float& h_z = grid.Extent().MeshWidth().Z();
+
+        int int_val = 0.;
+    vmg_float dir_x = pos_beg_x;
+    for (int i=-1*near_field_cells; i<=near_field_cells; ++i) {
+      vmg_float dir_y = pos_beg_y;
+      for (int j=-1*near_field_cells; j<=near_field_cells; ++j) {
+        vmg_float dir_z = pos_beg_z;
+        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));
+
+          dir_z += h_z;
+        }
+        dir_y += h_y;
+      }
+      dir_x += h_x;
+    }
+    // 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);
+    else
+        int_val = 1.0;
+
+    // Iterate over all grid points which lie in the support of the interpolating B-Spline
+        dir_x = pos_beg_x;
+        for (int i=-1*near_field_cells; i<=near_field_cells; ++i) {
+          vmg_float dir_y = pos_beg_y;
+          for (int j=-1*near_field_cells; j<=near_field_cells; ++j) {
+                vmg_float dir_z = pos_beg_z;
+                for (int k=-1*near_field_cells; k<=near_field_cells; ++k) {
+
+                  // Compute distance from grid point to particle
+                  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()) {
+                    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);
+                  } else {
+                          std::cerr << "Value very small " << length_sq << "=("
+                                          << dir_x << ")^2+(" << dir_y << ")^2+(" << dir_z << ")^2 for particle at ("
+                                          << p.Pos().X() << ","
+                                          << p.Pos().Y() << ","
+                                          << p.Pos().Z() << ")"
+                                          << "\n";
+                  }
+                  dir_z += h_z;
+                }
+                dir_y += h_y;
+          }
+          dir_x += h_x;
+        }
+  }
+
   vmg_float EvaluatePotential(const vmg_float& val) const
   {