From 3170267f202f5c8983985adb6b6048d19f018719 Mon Sep 17 00:00:00 2001
From: Pietro Incardona <i-bird@linux.fritz.box>
Date: Wed, 10 Aug 2016 00:48:32 +0200
Subject: [PATCH] fixing examples
---
example/Vector/4_reorder/energy_force.hpp | 138 ++++++++++++++++++++++
1 file changed, 138 insertions(+)
create mode 100644 example/Vector/4_reorder/energy_force.hpp
diff --git a/example/Vector/4_reorder/energy_force.hpp b/example/Vector/4_reorder/energy_force.hpp
new file mode 100644
index 000000000..5836b0fc8
--- /dev/null
+++ b/example/Vector/4_reorder/energy_force.hpp
@@ -0,0 +1,138 @@
+/*
+ * energy_force.hpp
+ *
+ * Created on: Aug 6, 2016
+ * Author: i-bird
+ */
+
+#ifndef EXAMPLE_VECTOR_4_REORDER_ENERGY_FORCE_HPP_
+#define EXAMPLE_VECTOR_4_REORDER_ENERGY_FORCE_HPP_
+
+constexpr int velocity = 0;
+constexpr int force = 1;
+
+
+void calc_forces(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, CellList<3, double, FAST, shift<3, double> > & NN, double sigma12, double sigma6)
+{
+ // update the Cell-list
+ vd.updateCellList(NN);
+
+ // Get an iterator over particles
+ auto it2 = vd.getDomainIterator();
+
+ // For each particle p ...
+ while (it2.isNext())
+ {
+ // ... get the particle p
+ auto p = it2.get();
+
+ // Get the position xp of the particle
+ Point<3,double> xp = vd.getPos(p);
+
+ // Reset the forice counter
+ vd.template getProp<force>(p)[0] = 0.0;
+ vd.template getProp<force>(p)[1] = 0.0;
+ vd.template getProp<force>(p)[2] = 0.0;
+
+ // Get an iterator over the neighborhood particles of p
+ auto Np = NN.template getNNIterator<NO_CHECK>(NN.getCell(vd.getPos(p)));
+
+ // For each neighborhood particle ...
+ while (Np.isNext())
+ {
+ // ... q
+ auto q = Np.get();
+
+ // if (p == q) skip this particle
+ if (q == p.getKey()) {++Np; continue;};
+
+ // Get the position of p
+ Point<3,double> xq = vd.getPos(q);
+
+ // Get the distance between p and q
+ Point<3,double> r = xp - xq;
+
+ // take the norm of this vector
+ double rn = norm2(r);
+
+ // Calculate the force, using pow is slower
+ Point<3,double> f = 24.0*(2.0 *sigma12 / (rn*rn*rn*rn*rn*rn*rn) - sigma6 / (rn*rn*rn*rn)) * r;
+
+ // we sum the force produced by q on p
+ vd.template getProp<force>(p)[0] += f.get(0);
+ vd.template getProp<force>(p)[1] += f.get(1);
+ vd.template getProp<force>(p)[2] += f.get(2);
+
+ // Next neighborhood
+ ++Np;
+ }
+
+ // Next particle
+ ++it2;
+ }
+}
+
+
+double calc_energy(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, CellList<3, double, FAST, shift<3, double> > & NN, double sigma12, double sigma6)
+{
+ double E = 0.0;
+
+ // Update the cell-list
+ vd.updateCellList(NN);
+
+ // Get the iterator
+ auto it2 = vd.getDomainIterator();
+
+ // For each particle ...
+ while (it2.isNext())
+ {
+ // ... p
+ auto p = it2.get();
+
+ // Get the position of the particle p
+ Point<3,double> xp = vd.getPos(p);
+
+ // Reset the force
+ vd.template getProp<force>(p)[0] = 0.0;
+ vd.template getProp<force>(p)[1] = 0.0;
+ vd.template getProp<force>(p)[2] = 0.0;
+
+ // Get an iterator over the neighborhood of the particle p
+ auto Np = NN.template getNNIterator<NO_CHECK>(NN.getCell(vd.getPos(p)));
+
+ // For each neighborhood of the particle p
+ while (Np.isNext())
+ {
+ // Neighborhood particle q
+ auto q = Np.get();
+
+ // if p == q skip this particle
+ if (q == p.getKey()) {++Np; continue;};
+
+ // Get position of the particle q
+ Point<3,double> xq = vd.getPos(q);
+
+ // take the normalized direction
+ double rn = norm2(xp - xq);
+
+ // potential energy (using pow is slower)
+ E += 4.0 * ( sigma12 / (rn*rn*rn*rn*rn*rn) - sigma6 / ( rn*rn*rn) );
+
+ // Next neighborhood
+ ++Np;
+ }
+
+ // Kinetic energy of the particle given by its actual speed
+ E += (vd.template getProp<velocity>(p)[0]*vd.template getProp<velocity>(p)[0] +
+ vd.template getProp<velocity>(p)[1]*vd.template getProp<velocity>(p)[1] +
+ vd.template getProp<velocity>(p)[2]*vd.template getProp<velocity>(p)[2]) / 2;
+
+ // Next Particle
+ ++it2;
+ }
+
+ return E;
+}
+
+
+#endif /* EXAMPLE_VECTOR_4_REORDER_ENERGY_FORCE_HPP_ */
--
GitLab