Changeset 1a92cf

Timestamp:

Apr 29, 2012, 10:50:41 AM (14 years ago)

Author:

Julian Iseringhausen <isering@…>

Children:

39a6d9

Parents:

4571da

Message:

vmg: Prepared force calculation.

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

Location:

src

Files:

• base/interpolate_polynomial.cpp (modified) (8 diffs)
• base/interpolate_polynomial.hpp (modified) (4 diffs)
• interface/interface_particles.cpp (modified) (3 diffs)
• samples/bspline.cpp (modified) (3 diffs)
• samples/bspline.hpp (modified) (2 diffs)

src/base/interpolate_polynomial.cpp

@@ -14 +14 @@
 
 InterpolatePolynomial::InterpolatePolynomial(const unsigned int& degree) :
+  deg(degree),
   deg_1(degree+1),
-  pos_begin(0.0)
+  buffer(degree+1),
+  buffer_diff(degree)
 {
   coeff = new vmg_float[Helper::intpow(deg_1, 3)];
+  coeff_buffer = new vmg_float[deg_1];
+  for (int i=0; i<degree; ++i)
+    buffer_diff[i].resize(i+1);
 }
 
@@ -23 +28 @@
 {
   delete [] coeff;
+  delete [] coeff_buffer;
 }
 
@@ -59 +65 @@
 void InterpolatePolynomial::_compute_coefficients_1d(const Index& index, const unsigned int& direction)
 {
-  vmg_float temp_array[deg_1];
   vmg_float power = 1.0;
   unsigned long faculty = 1;
@@ -66 +71 @@
 
   for (i=index, c=0; c<deg_1; ++i[direction], ++c)
-    temp_array[c] = _access_coeff(i);
+    coeff_buffer[c] = _access_coeff(i);
 
   i=index;
@@ -72 +77 @@
   for (c=1; c<deg_1; ++i[direction], ++c) {
     for (unsigned int j=0; j<deg_1-c; ++j)
-      temp_array[j] = temp_array[j+1] - temp_array[j];
+      coeff_buffer[j] = coeff_buffer[j+1] - coeff_buffer[j];
     faculty *= c;
     power *= h[direction];
-    _access_coeff(i) = temp_array[0] / (faculty*power);
+    _access_coeff(i) = coeff_buffer[0] / (faculty*power);
   }
 }
@@ -86 +91 @@
 
   prod[0] = 1.0;
-  offset[0] = pos_begin[0];
+  offset[0] = pos[0] - pos_begin[0];
   for (i[0]=0; i[0]<deg_1; ++i[0]) {
     prod[1] = prod[0];
-    offset[1] = pos_begin[1];
+    offset[1] = pos[1] - pos_begin[1];
     for (i[1]=0; i[1]<deg_1; ++i[1]) {
       prod[2] = prod[1];
-      offset[2] = pos_begin[2];
+      offset[2] = pos[2] - pos_begin[2];
       for (i[2]=0; i[2]<deg_1; ++i[2]) {
         result += _access_coeff(i) * prod[2];
-        prod[2] *= pos[2] - offset[2];
-        offset[2] += h[2];
+        prod[2] *= offset[2];
+        offset[2] -= h[2];
       }
-      prod[1] *= pos[1] - offset[1];
-      offset[1] += h[1];
+      prod[1] *= offset[1];
+      offset[1] -= h[1];
     }
-    prod[0] *= pos[0] - offset[0];
-    offset[0] += h[0];
+    prod[0] *= offset[0];
+    offset[0] -= h[0];
   }
 
@@ -108 +113 @@
 }
 
-Vector InterpolatePolynomial::EvaluateGradient(const Vector& pos)
+void InterpolatePolynomial::Evaluate(const Vector& pos, vmg_float& pot, Vector& field)
 {
-  Vector buffer[deg_1];
-  std::vector<Vector> buffer_diff[deg_1-1];
+ Vector offset;
+
+ pot = 0.0;
+ field = 0.0;
+
+  buffer[0] = 1.0;
+  offset = pos - pos_begin;
+  for (int i=1; i<deg_1; ++i) {
+    buffer[i] = buffer[i-1] * offset;
+    offset -= h;
+  }
+
+  offset = pos - pos_begin;
+  for (int i=0; i<deg; ++i) {
+    for (int j=0; j<i; ++j)
+      buffer_diff[i][j] = buffer_diff[i-1][j] * offset;
+    buffer_diff[i][i] = buffer[i];
+    offset -= h;
+  }
+
+  for (int i=1; i<deg; ++i)
+    for (int j=1; j<=i; ++j)
+      buffer_diff[i][0] += buffer_diff[i][j];
+
+  for (int i=0; i<deg_1; ++i)
+    for (int j=0; j<deg_1; ++j)
+      for (int k=0; k<deg_1; ++k)
+        pot += _access_coeff(i,j,k) * buffer[i][0] * buffer[j][1] * buffer[k][2];
+
+  for (int i=0; i<deg_1; ++i)
+    for (int j=0; j<deg_1; ++j)
+      for (int k=0; k<deg; ++k) {
+        field[0] -= _access_coeff(k+1, i, j) * buffer[i][1] * buffer[j][2] * buffer_diff[k][0][0];
+        field[1] -= _access_coeff(i, k+1, j) * buffer[i][0] * buffer[j][2] * buffer_diff[k][0][1];
+        field[2] -= _access_coeff(i, j, k+1) * buffer[i][0] * buffer[j][1] * buffer_diff[k][0][2];
+      }
+
+}
+
+Vector InterpolatePolynomial::EvaluateNegGradient(const Vector& pos)
+{
   Vector offset, result;
 
@@ -118 +162 @@
   for (int i=1; i<deg_1; ++i) {
     buffer[i] = buffer[i-1] * offset;
-    offset += h;
+    offset -= h;
   }
 
   offset = pos - pos_begin;
-  for (int i=0; i<deg_1-1; ++i) {
-    buffer_diff[i].resize(i+1);
-    offset += h;
+  for (int i=0; i<deg; ++i) {
     for (int j=0; j<i; ++j)
       buffer_diff[i][j] = buffer_diff[i-1][j] * offset;
-    buffer_diff[i][i] = buffer[i+1];
+    buffer_diff[i][i] = buffer[i];
+    offset -= h;
   }
+
+  for (int i=1; i<deg; ++i)
+    for (int j=1; j<=i; ++j)
+      buffer_diff[i][0] += buffer_diff[i][j];
 
   for (int i=0; i<deg_1; ++i)
     for (int j=0; j<deg_1; ++j)
-      for (int k=0; k<deg_1-1; ++k) {
-        Vector sum;
-        for (int l=0; l<=k; ++l)
-          sum += buffer_diff[k][l];
-        result[0] += _access_coeff(k+1, i, j) * buffer[i][1] * buffer[j][2] * sum[0];
-        result[1] += _access_coeff(i, k+1, j) * buffer[i][0] * buffer[j][2] * sum[1];
-        result[2] += _access_coeff(i, j, k+1) * buffer[i][0] * buffer[j][1] * sum[2];
+      for (int k=0; k<deg; ++k) {
+        result[0] -= _access_coeff(k+1, i, j) * buffer[i][1] * buffer[j][2] * buffer_diff[k][0][0];
+        result[1] -= _access_coeff(i, k+1, j) * buffer[i][0] * buffer[j][2] * buffer_diff[k][0][1];
+        result[2] -= _access_coeff(i, j, k+1) * buffer[i][0] * buffer[j][1] * buffer_diff[k][0][2];
       }
 

src/base/interpolate_polynomial.hpp

@@ -2 +2 @@
 #define INTERPOLATE_POLYNOMIAL_HPP_
 
-#include <base/vector.hpp>
+#include <vector>
+
+#include "base/index.hpp"
+#include "base/vector.hpp"
 
 namespace VMG
@@ -8 +11 @@
 
 class Grid;
-class Index;
 
 class InterpolatePolynomial
@@ -19 +21 @@
 
   vmg_float Evaluate(const Vector& pos);
-  Vector EvaluateGradient(const Vector& pos);
+  void Evaluate(const Vector& pos, vmg_float& pot, Vector& field);
+  Vector EvaluateNegGradient(const Vector& pos);
 
 private:
@@ -35 +38 @@
 
   vmg_float* coeff;
-  unsigned int deg_1;
+  vmg_float* coeff_buffer;
+  unsigned int deg, deg_1;
   Vector pos_begin;
   Vector h;
+  std::vector<Vector> buffer;
+  std::vector< std::vector<Vector> > buffer_diff;
 };
 

src/interface/interface_particles.cpp

@@ -101 +101 @@
 {
   Index i;
-  vmg_float length;
 
 #ifdef DEBUG_OUTPUT
@@ -161 +160 @@
     for (p1=lc_1.begin(); p1!=lc_1.end(); ++p1) {
 
-      (*p1)->Pot() = ip.Evaluate((*p1)->Pos()) - (*p1)->Charge() * spl.GetAntiDerivativeAtZero();
+      ip.Evaluate((*p1)->Pos(), (*p1)->Pot(), (*p1)->Field());
+      (*p1)->Pot() -= (*p1)->Charge() * spl.GetAntiDerivativeAtZero();
 
 #ifdef DEBUG_OUTPUT
@@ -177 +177 @@
               if (*p1 != *p2) {
 
-                length = ((*p2)->Pos() - (*p1)->Pos()).Length();
+                const Vector dir = (*p1)->Pos() - (*p2)->Pos();
+                const vmg_float length = dir.Length();
 
                 if (length < r_cut) {
+
                   (*p1)->Pot() += (*p2)->Charge() / length * (1.0 + spl.EvaluatePotential(length));
+                  (*p1)->Field() += (*p2)->Charge() * dir * spl.EvaluateField(length);
 
 #ifdef DEBUG_OUTPUT

src/samples/bspline.cpp

@@ -33 +33 @@
   potential_nom((BSPLINE_DEGREE+1)/2+1),
   potential_denom((BSPLINE_DEGREE+1)/2+1),
+  field_nom((BSPLINE_DEGREE+1)/2),
+  field_denom((BSPLINE_DEGREE+1)/2),
   intervals((BSPLINE_DEGREE+1)/2),
   R(width)
@@ -68 +70 @@
   potential_denom[1] = Polynomial(0, 160.0*POW(R,5));
   potential_denom[2] = Polynomial(0, 1.0);
+
+  field_nom[0] = Polynomial(5,
+                            20.0 * POW(R,5),
+                            0.0,
+                            0.0,
+                            -135.0 * POW(R,2),
+                            0.0,
+                            243.0);
+
+  field_nom[1] = Polynomial(5,
+                            81.0 * POW(R,5),
+                            0.0,
+                            0.0,
+                            -810.0 * POW(R,2),
+                            1215.0 * R,
+                            -486.0);
+
+  field_denom[0] = Polynomial(3,
+                              0.0,
+                              0.0,
+                              0.0,
+                              20.0 * POW(R,5));
+
+  field_denom[1] = Polynomial(3,
+                              0.0,
+                              0.0,
+                              0.0,
+                              80.0 * POW(R,5));
 
   antid = 39.0 / (16.0 * R);
@@ -141 +171 @@
   potential_denom[3] = Polynomial(0, 0.0);
 
+  field_nom[0] = Polynomial(8,
+                            280.0    * POW(R,8),
+                            0.0,
+                            0.0,
+                            -5544.0  * POW(R,5),
+                            0.0,
+                            27216.0  * POW(R,3),
+                            0.0,
+                            -87480.0 * R,
+                            76545.0);
+
+  field_nom[1] = Polynomial(8,
+                            1675.0 * POW(R,8),
+                            0.0,
+                            0.0,
+                            -25704.0 * POW(R,5),
+                            -85050.0 * POW(R,4),
+                            571536.0 * POW(R,3),
+                            -1020600.0 * POW(R,2),
+                            787320.0 * R,
+                            -229635.0);
+
+  field_nom[2] = Polynomial(8,
+                            729.0 * POW(R,8),
+                            0.0,
+                            0.0,
+                            -40824.0 * POW(R,5),
+                            153090.0 * POW(R,4),
+                            -244944.0 * POW(R,3),
+                            204120.0 * POW(R,2),
+                            -87480.0 * R,
+                            15309.0);
+
+  field_denom[0] = Polynomial(3,
+                              0.0,
+                              0.0,
+                              0.0,
+                              280.0 * POW(R,8));
+  field_denom[1] = Polynomial(3,
+                              0.0,
+                              0.0,
+                              0.0,
+                              1680.0 * POW(R,8));
+
+  field_denom[2] = Polynomial(3,
+                              0.0,
+                              0.0,
+                              0.0,
+                              560.0 * POW(R,8));
+
   antid = 239.0 / (70.0 * R);
 

src/samples/bspline.hpp

@@ -91 +91 @@
   }
 
+  vmg_float EvaluateField(const vmg_float& val) const
+  {
+    for (unsigned int i=0; i<intervals.size(); ++i)
+      if (val < intervals[i])
+        return field_nom[i](val) / field_denom[i](val);
+    return 0.0;
+  }
+
   const vmg_float& GetAntiDerivativeAtZero() const
   {
@@ -99 +107 @@
   std::vector<Polynomial> spline_nom, spline_denom;
   std::vector<Polynomial> potential_nom, potential_denom;
+  std::vector<Polynomial> field_nom, field_denom;
   vmg_float antid;
   std::vector<vmg_float> intervals;