From d28f7bd39f5cf816674f33db85e4498826ef6f74 Mon Sep 17 00:00:00 2001
From: Pietro Incardona <incardon@mpi-cbg.de>
Date: Wed, 21 Sep 2016 22:45:37 +0200
Subject: [PATCH] Fixing example
---
.../3_molecular_dynamic/main_expr_vl.cpp | 151 +++++
.../Vector/3_molecular_dynamic/main_vl.cpp | 482 +++++++++++++++
.../3_molecular_dynamic/main_vl_sym.cpp | 552 ++++++++++++++++++
example/Vector/4_complex_prop/main_ser.cpp | 504 ++++++++++++++++
example/Vector/4_reorder/main_expr.cpp | 139 +++++
.../vector_dist_complex_prp_unit_test.hpp | 231 ++++++++
6 files changed, 2059 insertions(+)
create mode 100644 example/Vector/3_molecular_dynamic/main_expr_vl.cpp
create mode 100644 example/Vector/3_molecular_dynamic/main_vl.cpp
create mode 100644 example/Vector/3_molecular_dynamic/main_vl_sym.cpp
create mode 100644 example/Vector/4_complex_prop/main_ser.cpp
create mode 100644 example/Vector/4_reorder/main_expr.cpp
create mode 100644 src/Vector/vector_dist_complex_prp_unit_test.hpp
diff --git a/example/Vector/3_molecular_dynamic/main_expr_vl.cpp b/example/Vector/3_molecular_dynamic/main_expr_vl.cpp
new file mode 100644
index 00000000..5e214636
--- /dev/null
+++ b/example/Vector/3_molecular_dynamic/main_expr_vl.cpp
@@ -0,0 +1,151 @@
+#include "Vector/vector_dist.hpp"
+#include "Plot/GoogleChart.hpp"
+#include "Operators/Vector/vector_dist_operators.hpp"
+#include "timer.hpp"
+
+constexpr int velocity = 0;
+constexpr int force = 1;
+
+struct ln_potential
+{
+ double sigma12,sigma6,r_cut2,shift;
+
+ ln_potential(double sigma12_, double sigma6_, double r_cut2_, double shift_) {sigma12 = sigma12_; sigma6 = sigma6_; r_cut2 = r_cut2_;shift = shift_;}
+
+ Point<2,double> value(const Point<3,double> & xp, const Point<3,double> xq)
+ {
+ double rn = norm2(xp - xq);
+ if (rn >= r_cut2) return 0.0;
+
+ Point<2,double> E({2.0 * ( sigma12 / (rn*rn*rn*rn*rn*rn) - sigma6 / ( rn*rn*rn) ) - shift,0.0});
+
+ return E;
+ }
+};
+
+struct ln_force
+{
+ double sigma12,sigma6,r_cut2;
+
+ ln_force(double sigma12_, double sigma6_, double r_cut2_) {sigma12 = sigma12_; sigma6 = sigma6_;r_cut2 = r_cut2_;}
+
+ Point<3,double> value(const Point<3,double> & xp, const Point<3,double> xq)
+ {
+ Point<3,double> r = xp - xq;
+ double rn = norm2(r);
+
+ if (rn > r_cut2) return 0.0;
+
+ return 24.0*(2.0 * sigma12 / (rn*rn*rn*rn*rn*rn*rn) - sigma6 / (rn*rn*rn*rn)) * r;
+ }
+};
+
+int main(int argc, char* argv[])
+{
+ double dt = 0.0005, sigma = 0.1, r_cut = 3.0*sigma;
+
+ double sigma6 = pow(sigma,6), sigma12 = pow(sigma,12);
+ double rc2 = r_cut * r_cut;
+ double shift = 2.0 * ( sigma12 / (rc2*rc2*rc2*rc2*rc2*rc2) - sigma6 / ( rc2*rc2*rc2) );
+
+ openfpm::vector<double> x;
+ openfpm::vector<openfpm::vector<double>> y;
+
+
+ openfpm_init(&argc,&argv);
+ Vcluster & vcl = create_vcluster();
+
+ size_t sz[3] = {10,10,10};
+ Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
+ size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
+ Ghost<3,float> ghost(r_cut);
+
+ ln_force lf(sigma12,sigma6,r_cut*r_cut);
+ ln_potential lp(sigma12,sigma6,r_cut*r_cut,shift);
+
+ vector_dist<3,double, aggregate<Point<3,double>,Point<3,double>> > vd(0,box,bc,ghost);
+
+ auto v_force = getV<force>(vd);
+ auto v_velocity = getV<velocity>(vd);
+ auto v_pos = getV<PROP_POS>(vd);
+
+ auto it = vd.getGridIterator(sz);
+
+ while (it.isNext())
+ {
+ vd.add();
+
+ auto key = it.get();
+
+ vd.getLastPos()[0] = key.get(0) * it.getSpacing(0);
+ vd.getLastPos()[1] = key.get(1) * it.getSpacing(1);
+ vd.getLastPos()[2] = key.get(2) * it.getSpacing(2);
+
+ ++it;
+ }
+
+ v_force = 0;
+ v_velocity = 0;
+
+ timer tsim;
+ tsim.start();
+
+ auto NN = vd.getCellList(r_cut);
+
+ vd.updateCellList(NN);
+ v_force = applyKernel_in_sim(vd,NN,lf);
+ unsigned long int f = 0;
+
+ // MD time stepping
+ for (size_t i = 0; i < 10000 ; i++)
+ {
+ assign(v_velocity, v_velocity + 0.5*dt*v_force,
+ v_pos, v_pos + v_velocity*dt);
+
+ vd.map();
+ vd.template ghost_get<>();
+
+ // calculate forces or a(tn + 1) Step 2
+ vd.updateCellList(NN);
+ v_force = applyKernel_in_sim(vd,NN,lf);
+
+ v_velocity = v_velocity + 0.5*dt*v_force;
+
+ if (i % 100 == 0)
+ {
+ vd.deleteGhost();
+ vd.write("particles_",f);
+ vd.ghost_get<>();
+
+ vd.updateCellList(NN);
+ Point<2,double> E = rsum(applyKernel_in_sim(vd,NN,lp) + (v_velocity * v_velocity)/2.0,vd).get();
+
+ vcl.sum(E.get(0));vcl.sum(E.get(1));
+ vcl.execute();
+
+ // we save the energy calculated at time step i c contain the time-step y contain the energy
+ x.add(i);
+ y.add({E.get(0),E.get(1),E.get(0) - E.get(1)});
+
+ if (vcl.getProcessUnitID() == 0)
+ std::cout << "Energy Total: " << E.get(0) << " Kinetic: " << E.get(1) << " Potential: " << E.get(0) - E.get(1) << std::endl;
+
+ f++;
+ }
+ }
+
+ tsim.stop();
+ std::cout << "Time: " << tsim.getwct() << std::endl;
+
+ GCoptions options;
+ options.title = std::string("Energy with time");
+ options.yAxis = std::string("Energy");
+ options.xAxis = std::string("iteration");
+ options.lineWidth = 1.0;
+
+ GoogleChart cg;
+ cg.AddLinesGraph(x,y,options);
+ cg.write("gc_plot2_out.html");
+
+ openfpm_finalize();
+}
diff --git a/example/Vector/3_molecular_dynamic/main_vl.cpp b/example/Vector/3_molecular_dynamic/main_vl.cpp
new file mode 100644
index 00000000..c039c991
--- /dev/null
+++ b/example/Vector/3_molecular_dynamic/main_vl.cpp
@@ -0,0 +1,482 @@
+
+#include "Vector/vector_dist.hpp"
+#include "Decomposition/CartDecomposition.hpp"
+#include "data_type/aggregate.hpp"
+#include "Plot/GoogleChart.hpp"
+#include "Plot/util.hpp"
+#include "timer.hpp"
+
+/*!
+ * \page Vector_3_md_vl Vector 3 molecular dynamic with Verlet list
+ *
+ *
+ * [TOC]
+ *
+ * # Molecular Dynamic with Lennard-Jones potential with verlet list # {#e3_md_vl}
+ *
+ * This example show a simple Lennard-Jones molecular dynamic simulation in a stable regime.
+ * We will use Verlet-list in order to get a speed-up from force calculation
+ *
+ * ## Constants ##
+ *
+ * Here we define some useful constants
+ *
+ * \snippet Vector/3_molecular_dynamic/main_vl.cpp constants
+ *
+ */
+
+//! \cond [constants] \endcond
+
+constexpr int velocity = 0;
+constexpr int force = 1;
+
+//! \cond [constants] \endcond
+
+/*!
+ *
+ * \page Vector_3_md_vl Vector 3 molecular dynamic with Verlet list
+ *
+ * ## Calculate forces ## {#e3_md_vl_cf}
+ *
+ * In this function we calculate the forces between particles. It require the vector of particles,
+ * the Verlet-list and sigma for the Lennard-Jhones potential. The function is exactly the same
+ * as the original with the following changes
+ *
+ * \see \ref e3_md_cf
+ *
+ * * The function accept a VerletList instead of a CellList
+ * \snippet main_vl.cpp arg diff
+ *
+ * * There is no call to updateCellList
+ *
+ * * How to get an iterator over neighborhood of a particle
+ * \snippet main_vl.cpp NN iterator
+ *
+ * Teh rest remain the same
+ *
+ * \snippet Vector/3_molecular_dynamic/main_vl.cpp calc forces vl
+ *
+ */
+
+//! \cond [calc forces vl] \endcond
+
+//! \cond [arg diff] \endcond
+
+void calc_forces(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, VerletList<3, double, FAST, shift<3, double> > & NN, double sigma12, double sigma6, double r_cut)
+{
+ //! \cond [arg diff] \endcond
+
+ // 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;
+
+ //! \cond [NN iterator] \endcond
+
+ // Get an iterator over the neighborhood particles of p
+ auto Np = NN.template getNNIterator<NO_CHECK>(p.getKey());
+
+ //! \cond [NN iterator] \endcond
+
+ // 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);
+
+ if (rn > r_cut * r_cut) {++Np;continue;}
+
+ // 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;
+ }
+}
+
+//! \cond [calc forces vl] \endcond
+
+/*!
+ * \page Vector_3_md_vl Vector 3 molecular dynamic with verlet list
+ *
+ * ## Calculate energy ## {#e3_md_vl_ce}
+ *
+ * We also need a function to calculate energy, this function require the same parameter as calculate forces
+ *
+ * \see \ref e3_md_ce
+ *
+ * The following changes has been made
+ *
+ * * The function accept a VerletList instead of a cell-List
+ * * There is no call to updateCellList
+ * * How to get an iterator over neigborhood particles
+ *
+ * \snippet Vector/3_molecular_dynamic/main_vl.cpp calc energy vl
+ *
+ */
+
+//! \cond [calc energy vl] \endcond
+
+double calc_energy(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, VerletList<3, double, FAST, shift<3, double> > & NN, double sigma12, double sigma6, double r_cut)
+{
+ double E = 0.0;
+
+ double rc = r_cut*r_cut;
+ double shift = 2.0 * ( sigma12 / (rc*rc*rc*rc*rc*rc) - sigma6 / ( rc*rc*rc) );
+
+ // 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);
+
+ // Get an iterator over the neighborhood of the particle p
+ auto Np = NN.template getNNIterator<NO_CHECK>(p.getKey());
+
+ // 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);
+
+ if (rn >= r_cut*r_cut)
+ {++Np;continue;}
+
+ // potential energy (using pow is slower)
+ E += 2.0 * ( sigma12 / (rn*rn*rn*rn*rn*rn) - sigma6 / ( rn*rn*rn) ) - shift;
+
+ // 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;
+ }
+
+ // Calculated energy
+ return E;
+}
+
+//! \cond [calc energy vl] \endcond
+
+int main(int argc, char* argv[])
+{
+ /*!
+ * \page Vector_3_md_vl Vector 3 molecular dynamic with Verlet list
+ *
+ * ## Simulation ## {#e3_md_vl_sim}
+ *
+ * The simulation is equal to the simulation explained in the example molecular dynamic
+ *
+ * \see \ref e3_md
+ *
+ * The differences are that:
+ *
+ * * The Ghost must be r_cut+skin
+ * \snippet
+ *
+ * * We create a Verlet list with skin instead of a Cell list
+ * \snippet Vector/3_molecular_dynamic/main_vl.cpp verlet skin
+ *
+ * * every 10 steps we do a map and update the verlet-list, in all the other case we just do skip labelling
+ * \snippet Vector/3_molecular_dynamic/main_vl.cpp update verlet
+ *
+ * **Explanation**
+ *
+ * Updating the verlet list is extremely expensive. For this reason we create a Verlet list
+ * that contain r_cut + skin particles. Using the fact that during the full simulation each
+ * particle does not move more than 0.0015 in one iteration, if the skin is 0.03
+ * we can update the Verlet list every \f$ \frac{0.03}{2 \cdot 0.0015} = 10 \f$. The 2 factor if given by the fact
+ * that in the worst case where one particle is going left and one on the right from the prospective of
+ * one particle the particle moove \f$ 2 \cdot 0.0015 \f$.
+ *
+ * Because the Verlet lists are constructed based on the local-id of the particles a map or a ghost_get
+ * would invalidate the verlet. For this reason the map is called every 10 time-step (when we
+ * update the verlet), and a particular ghost_get with SKIP_LABELLING is used during every iteration.
+ *
+ * The function ghost_get with skip labeling does not recompute the particle to send but use the
+ * the ids of the old particles updating the positions (and properties if needed) and keeping the old
+ * indexes without invalidating the Verlet-list. Doing this we can avoid to send particles that are
+ * entering the ghost area r_cut+skin. Because we know that no particle in 10 iteration can travel for a
+ * distance bigger than the skin, we are sure that in 10 iteration no-new particle that were not in the
+ * r_cut+skin ghost area can enter the ghost area r_cut.
+ *
+ *
+ * \snippet Vector/3_molecular_dynamic/main_vl.cpp simulation
+ *
+ */
+
+ //! \cond [simulation] \endcond
+
+ double dt = 0.00025;
+ double sigma = 0.1;
+ double r_cut = 3.0*sigma;
+ double r_gskin = 1.3*r_cut;
+ double sigma12 = pow(sigma,12);
+ double sigma6 = pow(sigma,6);
+
+ openfpm::vector<double> x;
+ openfpm::vector<openfpm::vector<double>> y;
+
+ openfpm_init(&argc,&argv);
+ Vcluster & v_cl = create_vcluster();
+
+ // we will use it do place particles on a 10x10x10 Grid like
+ size_t sz[3] = {10,10,10};
+
+ // domain
+ Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
+
+ // Boundary conditions
+ size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
+
+ // ghost, big enough to contain the interaction radius
+ Ghost<3,float> ghost(r_gskin);
+
+ vector_dist<3,double, aggregate<double[3],double[3]> > vd(0,box,bc,ghost);
+
+ auto it = vd.getGridIterator(sz);
+
+ while (it.isNext())
+ {
+ vd.add();
+
+ auto key = it.get();
+
+ vd.getLastPos()[0] = key.get(0) * it.getSpacing(0);
+ vd.getLastPos()[1] = key.get(1) * it.getSpacing(1);
+ vd.getLastPos()[2] = key.get(2) * it.getSpacing(2);
+
+ vd.template getLastProp<velocity>()[0] = 0.0;
+ vd.template getLastProp<velocity>()[1] = 0.0;
+ vd.template getLastProp<velocity>()[2] = 0.0;
+
+ vd.template getLastProp<force>()[0] = 0.0;
+ vd.template getLastProp<force>()[1] = 0.0;
+ vd.template getLastProp<force>()[2] = 0.0;
+
+ ++it;
+ }
+
+ timer tsim;
+ tsim.start();
+
+ //! \cond [verlet skin] \endcond
+
+ // Get the Cell list structure
+ auto NN = vd.getVerlet(r_gskin);
+
+ //! \cond [verlet skin] \endcond
+
+ // calculate forces
+ calc_forces(vd,NN,sigma12,sigma6,r_cut);
+ unsigned long int f = 0;
+
+ int cnt = 0;
+ double max_disp = 0.0;
+
+ // MD time stepping
+ for (size_t i = 0; i < 10000 ; i++)
+ {
+ // Get the iterator
+ auto it3 = vd.getDomainIterator();
+
+ double max_displ = 0.0;
+
+ // integrate velicity and space based on the calculated forces (Step1)
+ while (it3.isNext())
+ {
+ auto p = it3.get();
+
+ // here we calculate v(tn + 0.5)
+ vd.template getProp<velocity>(p)[0] += 0.5*dt*vd.template getProp<force>(p)[0];
+ vd.template getProp<velocity>(p)[1] += 0.5*dt*vd.template getProp<force>(p)[1];
+ vd.template getProp<velocity>(p)[2] += 0.5*dt*vd.template getProp<force>(p)[2];
+
+ Point<3,double> disp({vd.template getProp<velocity>(p)[0]*dt,vd.template getProp<velocity>(p)[1]*dt,vd.template getProp<velocity>(p)[2]*dt});
+
+ // here we calculate x(tn + 1)
+ vd.getPos(p)[0] += disp.get(0);
+ vd.getPos(p)[1] += disp.get(1);
+ vd.getPos(p)[2] += disp.get(2);
+
+ if (disp.norm() > max_displ)
+ max_displ = disp.norm();
+
+ ++it3;
+ }
+
+ if (max_disp < max_displ)
+ max_disp = max_displ;
+
+ //! \cond [update verlet] \endcond
+
+ // Because we moved the particles in space we have to map them and re-sync the ghost
+ if (cnt % 10 == 0)
+ {
+ vd.map();
+ vd.template ghost_get<>();
+ // Get the Cell list structure
+ vd.updateVerlet(NN,r_gskin);
+ }
+ else
+ {
+ vd.template ghost_get<>(SKIP_LABELLING);
+ }
+
+ //! \cond [update verlet] \endcond
+
+ cnt++;
+
+ // calculate forces or a(tn + 1) Step 2
+ calc_forces(vd,NN,sigma12,sigma6,r_cut);
+
+ // Integrate the velocity Step 3
+ auto it4 = vd.getDomainIterator();
+
+ while (it4.isNext())
+ {
+ auto p = it4.get();
+
+ // here we calculate v(tn + 1)
+ vd.template getProp<velocity>(p)[0] += 0.5*dt*vd.template getProp<force>(p)[0];
+ vd.template getProp<velocity>(p)[1] += 0.5*dt*vd.template getProp<force>(p)[1];
+ vd.template getProp<velocity>(p)[2] += 0.5*dt*vd.template getProp<force>(p)[2];
+
+ ++it4;
+ }
+
+ // After every iteration collect some statistic about the confoguration
+ if (i % 100 == 0)
+ {
+ // We write the particle position for visualization (Without ghost)
+ vd.deleteGhost();
+ vd.write("particles_",f);
+
+ // we resync the ghost
+ vd.ghost_get<>(SKIP_LABELLING);
+
+ // We calculate the energy
+ double energy = calc_energy(vd,NN,sigma12,sigma6,r_cut);
+ auto & vcl = create_vcluster();
+ vcl.sum(energy);
+ vcl.max(max_disp);
+ vcl.execute();
+
+ // we save the energy calculated at time step i c contain the time-step y contain the energy
+ x.add(i);
+ y.add({energy});
+
+ // We also print on terminal the value of the energy
+ // only one processor (master) write on terminal
+ if (vcl.getProcessUnitID() == 0)
+ std::cout << "Energy: " << energy << " " << max_disp << " " << std::endl;
+
+ max_disp = 0.0;
+
+ f++;
+ }
+ }
+
+ tsim.stop();
+ std::cout << "Time: " << tsim.getwct() << std::endl;
+
+ //! \cond [simulation] \endcond
+
+ // Google charts options, it store the options to draw the X Y graph
+ GCoptions options;
+
+ // Title of the graph
+ options.title = std::string("Energy with time");
+
+ // Y axis name
+ options.yAxis = std::string("Energy");
+
+ // X axis name
+ options.xAxis = std::string("iteration");
+
+ // width of the line
+ options.lineWidth = 1.0;
+
+ // Object that draw the X Y graph
+ GoogleChart cg;
+
+ // Add the graph
+ // The graph that it produce is in svg format that can be opened on browser
+ cg.AddLinesGraph(x,y,options);
+
+ // Write into html format
+ cg.write("gc_plot2_out.html");
+
+ //! \cond [google chart] \endcond
+
+ /*!
+ * \page Vector_3_md_vl Vector 3 molecular dynamic with Verlet list
+ *
+ * ## Finalize ## {#finalize_v_e3_md_vl}
+ *
+ * At the very end of the program we have always to de-initialize the library
+ *
+ * \snippet Vector/3_molecular_dynamic/main_vl.cpp finalize
+ *
+ */
+
+ //! \cond [finalize] \endcond
+
+ openfpm_finalize();
+
+ //! \cond [finalize] \endcond
+}
diff --git a/example/Vector/3_molecular_dynamic/main_vl_sym.cpp b/example/Vector/3_molecular_dynamic/main_vl_sym.cpp
new file mode 100644
index 00000000..07136e19
--- /dev/null
+++ b/example/Vector/3_molecular_dynamic/main_vl_sym.cpp
@@ -0,0 +1,552 @@
+#include "Vector/vector_dist.hpp"
+#include "Decomposition/CartDecomposition.hpp"
+#include "data_type/aggregate.hpp"
+#include "Plot/GoogleChart.hpp"
+#include "Plot/util.hpp"
+#include "timer.hpp"
+
+/*!
+ * \page Vector_3_md_vl Vector 3 molecular dynamic with Verlet list
+ *
+ *
+ * ## Variations for symmetric interaction ##
+ *
+ * In this example we show the variation in case of symmetric interations.
+ *
+ */
+
+//! \cond [constants] \endcond
+
+constexpr int velocity = 0;
+constexpr int force = 1;
+
+//! \cond [constants] \endcond
+
+/*!
+ *
+ * \page Vector_3_md_vl Vector 3 molecular dynamic with Verlet list symmetric
+ *
+ * ## Calculate forces ##
+ *
+ * In this function we calculate the forces between particles. It require the vector of particles
+ * Cell list and sigma factor for the Lennard-Jhones potential. The function is exactly the same
+ * as the original with the following changes
+ *
+ * \see \ref e3_md_vl_cf
+ *
+ * * We have to reset for force counter for each particles
+ * \snippet Vector/3_molecular_dynamic/main_vl_sym.cpp reset
+ *
+ * * The iterator go through real AND ghost.
+ * \snippet Vector/3_molecular_dynamic/main_vl_sym.cpp real and ghost
+ *
+ * * If we calculate the force for p-q we are also adding this force to q-p
+ * \snippet Vector/3_molecular_dynamic/main_vl_sym.cpp add to q
+ *
+ * ** Explanation**
+ *
+ * The reason why symmetric Verlet exist is to reduce the computation. In a symmetric interaction
+ * we have that the interaction p-q is is the same as q-p. This mean that we are calculating 2 times
+ * the same quantity one when we are on particle p and one when we are on particle q. Symmetric verlet
+ * avoid this. Every time we calculate a force p-q we also increment the force q-p.
+ *
+ * The reason instead why we have to go through ghost is the following
+ *
+ *
+ \verbatim
+
+ | |
+ | Ghost|
+ | |
+ | |
+ | |
+ | |
+ | |
+ | +----+----+----+
+ | | | | | |
+ | | 1| | 2 | 3 |
+ | | | | | |
+ | +--------------+
+ | | | * | 4 |
+ | +--------------------
+ | X | 5 | |
+ | +---------+
+ +---------------------------
+
+ \endverbatim
+ *
+ * Where * is a particle and X is another particle on ghost
+ *
+ * The symmetric verlet list is constructed from the cell-List 1,2,3,4,5.
+ * As you can see the interaction X-* is not present in the neighborhood of *
+ * (cells 1,2,3,4,5). But is present if we draw the neighborhood of X. This
+ * mean we have to iterate also across X if we want the interation X-*
+ *
+ *
+ * \snippet Vector/3_molecular_dynamic/main_vl_sym.cpp calc forces vl
+ *
+ * The symmetric force calculation is faster than the normal one. but is not 2
+ * time faster. The reason can be founded in 2 factors.
+ *
+ * * We have to reset the force
+ *
+ * * Every particle in the ghost
+ * that interact with one in the domain is still calculated twice by the same
+ * processor or some other.
+ *
+ * * The access pattern is worst we alternate read and write stressing more the Memory
+ *
+ * * We have to iterate across domain + ghost
+ *
+ * Overall the function calc_force is faster
+ *
+ * \snippet Vector/3_molecular_dynamic/main_vl_sym.cpp calc forces vl
+ *
+ */
+
+//! \cond [calc forces vl] \endcond
+
+void calc_forces(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, VerletList<3, double, FAST, shift<3, double> > & NN, double sigma12, double sigma6, double r_cut)
+{
+ //! \cond [reset] \endcond
+
+ // Reset the force for all particles
+ auto it3 = vd.getDomainIterator();
+
+ // For each particle p ...
+ while (it3.isNext())
+ {
+ // ... get the particle p
+ auto p = it3.get();
+
+ // 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;
+
+ // Next particle p
+ ++it3;
+ }
+
+ //! \cond [reset] \endcond
+
+ //! \cond [real and ghost] \endcond
+
+ // Get an iterator over particles
+ auto it2 = vd.getDomainAndGhostIterator();
+
+ //! \cond [real and ghost] \endcond
+
+ // 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);
+
+ // Get an iterator over the neighborhood particles of p
+ // Note that in case of symmetric
+ auto Np = NN.template getNNIterator<NO_CHECK>(p.getKey());
+
+ // For each neighborhood particle ...
+ while (Np.isNext())
+ {
+ // ... q
+ auto q = Np.get();
+
+ // if (p == q) skip this particle
+ if (q == p.getKey() || (p.getKey() >= vd.size_local() && q >= vd.size_local())) {++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);
+
+ if (rn > r_cut * r_cut) {++Np;continue;}
+
+ // 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);
+
+ //! \cond [add to q] \endcond
+
+ // we sum the force produced by p on q
+ vd.template getProp<force>(q)[0] -= f.get(0);
+ vd.template getProp<force>(q)[1] -= f.get(1);
+ vd.template getProp<force>(q)[2] -= f.get(2);
+
+ //! \cond [add to q] \endcond
+
+ // Next neighborhood
+ ++Np;
+ }
+
+ // Next particle
+ ++it2;
+ }
+}
+
+//! \cond [calc forces vl] \endcond
+
+
+/*!
+ * \page Vector_3_md_vl Vector 3 molecular dynamic with Verlet list symmetric
+ *
+ * ## Calculate energy ##
+ *
+ * No changes that been made to this code compared to the molecular dynamic example
+ * with Verlet-list
+ *
+ * \see \ref e3_md_vl
+ *
+ *
+ * \snippet Vector/3_molecular_dynamic/main_vl_sym.cpp calc energy vl
+ *
+ */
+
+//! \cond [calc energy vl] \endcond
+
+double calc_energy(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, VerletList<3, double, FAST, shift<3, double> > & NN, double sigma12, double sigma6, double r_cut)
+{
+ double E = 0.0;
+
+ double rc = r_cut*r_cut;
+ double shift = 2.0 * ( sigma12 / (rc*rc*rc*rc*rc*rc) - sigma6 / ( rc*rc*rc) );
+
+ // 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);
+
+ // Get an iterator over the neighborhood of the particle p
+ auto Np = NN.template getNNIterator<NO_CHECK>(p.getKey());
+
+ double Ep = E;
+
+ // 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() && q < vd.size_local()) {++Np; continue;};
+
+ // Get position of the particle q
+ Point<3,double> xq = vd.getPos(q);
+
+ // take the normalized direction
+ double rn = norm2(xp - xq);
+
+ if (rn >= r_cut*r_cut)
+ {++Np;continue;}
+
+ // potential energy (using pow is slower)
+ E += 2.0 * ( sigma12 / (rn*rn*rn*rn*rn*rn) - sigma6 / ( rn*rn*rn) ) - shift;
+
+ // 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;
+ }
+
+ // Calculated energy
+ return E;
+}
+
+//! \cond [calc energy vl] \endcond
+
+int main(int argc, char* argv[])
+{
+ /*!
+ * \page Vector_3_md_vl Vector 3 molecular dynamic with Verlet list symmetric
+ *
+ * ## Simulation ##
+ *
+ * The simulation is equal to the simulation explained in the example molecular dynamic
+ *
+ * \see \ref e3_md_vl
+ *
+ * The differences are that:
+ *
+ * * Create symmetric Verlet-list
+ * \snippet Vector/3_molecular_dynamic/main_vl_sym.cpp sim verlet
+ *
+ * * For Energy calculation we have to create another verlet in this case a normal Verlet-list
+ * \snippet Vector/3_molecular_dynamic/main_vl_sym.cpp norm verlet
+ *
+ * We need 2 verlet-list with one normal because we cannot use
+ * symmetric Verlet to calculate total reductions like total energy
+ *
+ * \snippet Vector/3_molecular_dynamic/main_vl_sym.cpp simulation
+ *
+ */
+
+ //! \cond [simulation] \endcond
+
+ double dt = 0.00025;
+ double sigma = 0.1;
+ double r_cut = 3.0*sigma;
+ double r_gskin = 1.3*r_cut;
+ double sigma12 = pow(sigma,12);
+ double sigma6 = pow(sigma,6);
+
+ openfpm::vector<double> x;
+ openfpm::vector<openfpm::vector<double>> y;
+
+ openfpm_init(&argc,&argv);
+ Vcluster & v_cl = create_vcluster();
+
+ // we will use it do place particles on a 10x10x10 Grid like
+ size_t sz[3] = {10,10,10};
+
+ // domain
+ Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
+
+ // Boundary conditions
+ size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
+
+ // ghost, big enough to contain the interaction radius
+ Ghost<3,float> ghost(r_gskin);
+
+ vector_dist<3,double, aggregate<double[3],double[3]> > vd(0,box,bc,ghost);
+
+ auto it = vd.getGridIterator(sz);
+
+ while (it.isNext())
+ {
+ vd.add();
+
+ auto key = it.get();
+
+ vd.getLastPos()[0] = key.get(0) * it.getSpacing(0);
+ vd.getLastPos()[1] = key.get(1) * it.getSpacing(1);
+ vd.getLastPos()[2] = key.get(2) * it.getSpacing(2);
+
+ vd.template getLastProp<velocity>()[0] = 0.0;
+ vd.template getLastProp<velocity>()[1] = 0.0;
+ vd.template getLastProp<velocity>()[2] = 0.0;
+
+ vd.template getLastProp<force>()[0] = 0.0;
+ vd.template getLastProp<force>()[1] = 0.0;
+ vd.template getLastProp<force>()[2] = 0.0;
+
+ ++it;
+ }
+
+ timer tsim;
+ tsim.start();
+
+ //! \cond [sim verlet] \endcond
+
+ // Get the Cell list structure
+ auto NN = vd.getVerlet(r_gskin,VL_SYMMETRIC);
+
+ //! \cond [sim verlet] \endcond
+
+ //! \cond [norm verlet] \endcond
+
+ auto NNE = vd.getVerlet(r_gskin);
+
+ //! \cond [norm verlet] \endcond
+
+ // calculate forces
+ calc_forces(vd,NN,sigma12,sigma6,r_cut);
+ unsigned long int f = 0;
+
+ int cnt = 0;
+ double max_disp = 0.0;
+
+ // MD time stepping
+ for (size_t i = 0; i < 10000 ; i++)
+ {
+ // Get the iterator
+ auto it3 = vd.getDomainIterator();
+
+ double max_displ = 0.0;
+
+ // integrate velicity and space based on the calculated forces (Step1)
+ while (it3.isNext())
+ {
+ auto p = it3.get();
+
+ // here we calculate v(tn + 0.5)
+ vd.template getProp<velocity>(p)[0] += 0.5*dt*vd.template getProp<force>(p)[0];
+ vd.template getProp<velocity>(p)[1] += 0.5*dt*vd.template getProp<force>(p)[1];
+ vd.template getProp<velocity>(p)[2] += 0.5*dt*vd.template getProp<force>(p)[2];
+
+ Point<3,double> disp({vd.template getProp<velocity>(p)[0]*dt,vd.template getProp<velocity>(p)[1]*dt,vd.template getProp<velocity>(p)[2]*dt});
+
+ // here we calculate x(tn + 1)
+ vd.getPos(p)[0] += disp.get(0);
+ vd.getPos(p)[1] += disp.get(1);
+ vd.getPos(p)[2] += disp.get(2);
+
+ if (disp.norm() > max_displ)
+ max_displ = disp.norm();
+
+ ++it3;
+ }
+
+ if (max_disp < max_displ)
+ max_disp = max_displ;
+
+ // Because we moved the particles in space we have to map them and re-sync the ghost
+ if (cnt % 10 == 0)
+ {
+ vd.map();
+ vd.template ghost_get<>();
+ // Get the Cell list structure
+ vd.updateVerlet(NN,r_gskin,VL_SYMMETRIC);
+ }
+ else
+ {
+ vd.template ghost_get<>(SKIP_LABELLING);
+ }
+
+ cnt++;
+
+ // calculate forces or a(tn + 1) Step 2
+ calc_forces(vd,NN,sigma12,sigma6,r_cut);
+
+ // Integrate the velocity Step 3
+ auto it4 = vd.getDomainIterator();
+
+ while (it4.isNext())
+ {
+ auto p = it4.get();
+
+ // here we calculate v(tn + 1)
+ vd.template getProp<velocity>(p)[0] += 0.5*dt*vd.template getProp<force>(p)[0];
+ vd.template getProp<velocity>(p)[1] += 0.5*dt*vd.template getProp<force>(p)[1];
+ vd.template getProp<velocity>(p)[2] += 0.5*dt*vd.template getProp<force>(p)[2];
+
+ ++it4;
+ }
+
+ // After every iteration collect some statistic about the confoguration
+ if (i % 100 == 0)
+ {
+ // We write the particle position for visualization (Without ghost)
+ vd.deleteGhost();
+ vd.write("particles_",f);
+
+ // we resync the ghost
+ vd.ghost_get<>(SKIP_LABELLING);
+
+ // update the verlet for energy calculation
+ vd.updateVerlet(NNE,r_gskin);
+
+ // We calculate the energy
+ double energy = calc_energy(vd,NNE,sigma12,sigma6,r_cut);
+ auto & vcl = create_vcluster();
+ vcl.sum(energy);
+ vcl.max(max_disp);
+ vcl.execute();
+
+ // we save the energy calculated at time step i c contain the time-step y contain the energy
+ x.add(i);
+ y.add({energy});
+
+ // We also print on terminal the value of the energy
+ // only one processor (master) write on terminal
+ if (vcl.getProcessUnitID() == 0)
+ std::cout << "Energy: " << energy << " " << max_disp << " " << std::endl;
+
+ max_disp = 0.0;
+
+ f++;
+ }
+ }
+
+ tsim.stop();
+ std::cout << "Time: " << tsim.getwct() << std::endl;
+
+ //! \cond [simulation] \endcond
+
+ // Google charts options, it store the options to draw the X Y graph
+ GCoptions options;
+
+ // Title of the graph
+ options.title = std::string("Energy with time");
+
+ // Y axis name
+ options.yAxis = std::string("Energy");
+
+ // X axis name
+ options.xAxis = std::string("iteration");
+
+ // width of the line
+ options.lineWidth = 1.0;
+
+ // Object that draw the X Y graph
+ GoogleChart cg;
+
+ // Add the graph
+ // The graph that it produce is in svg format that can be opened on browser
+ cg.AddLinesGraph(x,y,options);
+
+ // Write into html format
+ cg.write("gc_plot2_out.html");
+
+ //! \cond [google chart] \endcond
+
+ /*!
+ * \page Vector_3_md_vl Vector 3 molecular dynamic with Verlet list symmetric
+ *
+ * ## Finalize ## {#finalize_v_e3_md_sym}
+ *
+ * At the very end of the program we have always to de-initialize the library
+ *
+ * \snippet Vector/3_molecular_dynamic/main_vl_sym.cpp finalize
+ *
+ */
+
+ //! \cond [finalize] \endcond
+
+ openfpm_finalize();
+
+ //! \cond [finalize] \endcond
+
+ /*!
+ * \page Vector_3_md_vl Vector 3 molecular dynamic
+ *
+ * ## Full code ## {#code_v_e3_md_vl}
+ *
+ * \include Vector/3_molecular_dynamic/main_vl.cpp
+ *
+ */
+
+ /*!
+ * \page Vector_3_md_vl Vector 3 molecular dynamic with Verlet list symmetric
+ *
+ * ## Full code symmetric ## {#code_v_e3_md_sym}
+ *
+ * \include Vector/3_molecular_dynamic/main_vl_sym.cpp
+ *
+ */
+}
diff --git a/example/Vector/4_complex_prop/main_ser.cpp b/example/Vector/4_complex_prop/main_ser.cpp
new file mode 100644
index 00000000..3a7e91db
--- /dev/null
+++ b/example/Vector/4_complex_prop/main_ser.cpp
@@ -0,0 +1,504 @@
+/*!
+ *
+ * \page Vector_4_complex_prop_ser Vector 4 property serialization
+ *
+ *
+ * [TOC]
+ *
+ *
+ * # Vector 4 property serialization # {#vector_example_cp_ser}
+ *
+ *
+ * This example show how we can use complex properties in a vector when we use a custom structure with
+ * a pointer inside. In particular we will show how to construct the serialization and de-serialization
+ * methods for the structure my_struct defined here
+ *
+ * \snippet Vector/4_complex_prop/main_ser.cpp my structure
+ *
+ * \snippet Vector/4_complex_prop/main_ser.cpp my structure end
+ *
+ */
+
+#include "Vector/vector_dist.hpp"
+
+//! \cond [my structure] \endcond
+
+struct my_struct
+{
+ //! C string size
+ size_t size;
+
+ //! C string pointer
+ char * ptr;
+
+ //! C++ string
+ std::string str;
+
+ //! vector
+ openfpm::vector<int> v;
+
+//! \cond [my structure] \endcond
+
+public:
+
+ //! Functions to check if the packing object is complex
+ static bool pack() {return true;}
+
+
+ //! \cond [con and dest] \endcond
+
+ my_struct()
+ {
+ size = 0;
+ ptr = NULL;
+ }
+
+ ~my_struct()
+ {
+ if (ptr != NULL)
+ delete [] ptr;
+ }
+
+ //! \cond [con and dest] \endcond
+
+ //! \cond [pack request] \endcond
+
+ //! Serialization request
+ template<int ... prp> inline void packRequest(size_t & req) const
+ {
+ req += 2*sizeof(size_t) + size + str.size()+1;
+ v.packRequest(req);
+ }
+
+ //! \cond [pack request] \endcond
+
+ //! \cond [pack serialize] \endcond
+
+ //! Serialize the data structure
+ template<typename Memory, int ... prp> inline void pack(ExtPreAlloc<Memory> & mem, Pack_stat & sts) const
+ {
+ //! \cond [packer ser pack] \endcond
+
+ //Serialize the number that determine the C-string size
+ Packer<size_t, Memory>::pack(mem,size,sts);
+
+ //! \cond [packer ser pack] \endcond
+
+ //! \cond [pack ser copy] \endcond
+
+ // Allocate and copy the string
+ mem.allocate(size);
+ void * t_ptr = mem.getPointer();
+ memcpy(t_ptr,ptr,size);
+
+ //! \cond [pack ser copy] \endcond
+
+ //! \cond [pack ser other] \endcond
+
+ //Serialize the number that determine the C++-string str.size()+1 given by the null terminator
+ Packer<size_t, Memory>::pack(mem,str.size()+1,sts);
+
+ // Allocate and copy the string
+ mem.allocate(str.size()+1);
+ char * t_ptr2 = (char *)mem.getPointer();
+ memcpy(t_ptr2,str.c_str(),str.size());
+
+ // Add null terminator
+ t_ptr2[str.size()] = 0;
+
+ Packer<decltype(v), Memory>::pack(mem,v,sts);
+
+ //! \cond [pack ser other] \endcond
+ }
+
+ //! \cond [pack serialize] \endcond
+
+ //! \cond [unpack de-serialize] \endcond
+
+ //! De-serialize the data structure
+ template<typename Memory, int ... prp> inline void unpack(ExtPreAlloc<Memory> & mem, Unpack_stat & ps)
+ {
+ //! \cond [unpacker ser pack] \endcond
+
+ // unpack the size of the C string
+ Unpacker<size_t, Memory>::unpack(mem,size,ps);
+
+ //! \cond [unpacker ser pack] \endcond
+
+ //! \cond [unpacker ser create] \endcond
+
+ // Allocate the C string
+ ptr = new char[size];
+
+ // source pointer
+ char * ptr_src = (char *)mem.getPointerOffset(ps.getOffset());
+
+ // copy from the buffer to the destination
+ memcpy(ptr,ptr_src,size);
+
+ // update the pointer
+ ps.addOffset(size);
+
+ //! \cond [unpacker ser create] \endcond
+
+ //! \cond [unpacker ser other] \endcond
+
+ // get the the C++ string size
+ size_t cpp_size;
+ Unpacker<size_t, Memory>::unpack(mem,cpp_size,ps);
+
+ // Create the string from the serialized data (de-serialize)
+ char * ptr_src2 = (char *)mem.getPointerOffset(ps.getOffset());
+ str = std::string(ptr_src2);
+ ps.addOffset(cpp_size);
+
+ // Unpack the vector
+ Unpacker<decltype(v), Memory>::unpack(mem,v,ps);
+
+ //! \cond [unpacker ser other] \endcond
+ }
+
+ //! \cond [unpack de-serialize] \endcond
+
+//! \cond [my structure end] \endcond
+
+};
+
+//! \cond [my structure end] \endcond
+
+/*!
+ *
+ * \page Vector_4_complex_prop_ser Vector 4 property serialization
+ *
+ * In order to make my_struct serializable we need 3 methods
+ *
+ * * **packRequest** This method indicate how many byte are needed to serialize this structure
+ * * **pack** This method serialize the data-structure
+ * * **unpack** This method de-serialize the data structure
+ *
+ * ### packRequest ###
+ *
+ * \snippet Vector/4_complex_prop/main_ser.cpp pack request
+ *
+ * This function calculate the size in byte to serialize this structure in this case we serialize
+ * 2 numbers so 2*sizeof(size_t). One C string of size "size" and one C++ string of size str.size(). Because
+ * std::string has a constructor from null terminating string we add the null terminator to easy construct an
+ * std::string. This mean that
+ * the C++ string will be serialized in str.size()+1 bytes. The result must be summed to counter req that is used
+ * to allocate a buffer big enough for serialization. The function is template and has a variadic argument
+ * "int ... prp" this can be ignored
+ *
+ * ### pack ###
+ *
+ * Here we serialize the object. The function give as argument
+ *
+ * * **ExtPreAlloc<Memory>** mem The memory where we have to serialize the information
+ * * **Pack_stat** unused
+ * * **int ... prp** unused
+ *
+ * The only important parameter is the object mem where we will serialize the information
+ *
+ * We first pack the number that indicate the size of the C string. A convenient way
+ * is use the function
+ *
+ * \snippet Vector/4_complex_prop/main_ser.cpp packer ser pack
+ *
+ * The function **Packer<decltype(object),Memory>::pack(object,sts)** work if **object** is a
+ * fundamental C++ type, a struct that does not contain pointers, any object that
+ * implement the interface pack,unpack,packRequest (like my_struct). Most openfpm
+ * data-structure like openfpm::vector and grid implement such interface and
+ * can be used directly with **Packer<...>::pack(...)**.
+ * After packed the size of the string we have to serialize the content of the pointer string.
+ * unfortunately there is no Packer<...>::pack(...) function to do this and must be
+ * manually done. This can be done with the following steps
+ *
+ * * Allocate the memory to copy the value of the string
+ * * Get the pointer to the memory allocated
+ * * copy the content of the string into the allocated memory
+ *
+ * \snippet Vector/4_complex_prop/main_ser.cpp pack ser copy
+ *
+ * For the C++ string the concept is the same, the difference is that we put a null terminator at the end
+ * of the serialized string
+ *
+ * \snippet Vector/4_complex_prop/main_ser.cpp pack ser other
+ *
+ * ### unpack ###
+ *
+ * Here we de-serialize the object. The function take a reference to object
+ *
+ * * **ExtPreAlloc<Memory>** mem, contain the serialized information that must be de-serialized
+ * * **Pack_stat** contain the offset to the memory to de-serialize
+ * * **int ... prp** unused
+ *
+ * De-serialization must be in general done in the same order as we serialized
+ * We first unpack the number that indicate the size of the C string. A convenient way
+ * is to use the function
+ *
+ * \snippet Vector/4_complex_prop/main_ser.cpp unpacker ser pack
+ *
+ * the variable **size** contain now the size of the packed string we can now
+ *
+ * * create memory with **new** that store the C string
+ * * Get the pointer to the serialized information
+ * * copy the string from mem into the created memory
+ * * update the offset pointer
+ *
+ * \snippet Vector/4_complex_prop/main_ser.cpp unpacker ser create
+ *
+ * For the C++ string is just the same
+ *
+ * * We get the size of the string
+ * * we create an std::string out of the null terminating serialized string
+ * * we assign the created std::string to str
+ *
+ * \snippet Vector/4_complex_prop/main_ser.cpp unpacker ser other
+ *
+ * ### Constructor and destructor ###
+ *
+ * Constructor and destructor are not releated to serialization and de-serialization concept.
+ * But on how my_struct is constructed and destructed in order to avoid memory
+ * leak/corruption. In the constructor we set ptr to NULL in the destructor we
+ * destroy the pointer (if different from NULL) to avoid memory leak
+ *
+ * \snippet Vector/4_complex_prop/main_ser.cpp con and dest
+ *
+ */
+
+int main(int argc, char* argv[])
+{
+ /*!
+ *
+ * \page Vector_4_complex_prop_ser Vector 4 property serialization
+ *
+ *
+ * ## Initialization and vector creation ##
+ *
+ * After we initialize the library we can create a vector with complex properties
+ * with the following line
+ *
+ * \snippet Vector/4_complex_prop/main.cpp vect create
+ *
+ * In this this particular case every particle carry two my_struct object
+ *
+ */
+
+ // initialize the library
+ openfpm_init(&argc,&argv);
+
+ // Here we define our domain a 2D box with internals from 0 to 1.0 for x and y
+ Box<2,float> domain({0.0,0.0},{1.0,1.0});
+
+ // Here we define the boundary conditions of our problem
+ size_t bc[2]={PERIODIC,PERIODIC};
+
+ // extended boundary around the domain, and the processor domain
+ Ghost<2,float> g(0.01);
+
+ // my_struct at position 0 in the aggregate
+ constexpr int my_s1 = 0;
+
+ // my_struct at position 1 in the aggregate
+ constexpr int my_s2 = 1;
+
+ //! \cond [vect create] \endcond
+
+ vector_dist<2,float, aggregate<my_struct,my_struct>>
+ vd(4096,domain,bc,g);
+
+ std::cout << "HAS PACK: " << has_pack_agg<aggregate<my_struct,my_struct>>::result::value << std::endl;
+
+ //! \cond [vect create] \endcond
+
+ /*!
+ *
+ * \page Vector_4_complex_prop_ser Vector 4 property serialization
+ *
+ *
+ * ## Assign values to properties ##
+ *
+ * In this loop we assign position to particles and we fill the two my_struct
+ * that each particle contain. As demostration the first my_struct is filled
+ * with the string representation of the particle coordinates. The second my struct
+ * is filled with the string representation of the particle position multiplied by 2.0.
+ * The the vectors of the two my_struct are filled respectively with the sequence
+ * 1,2,3 and 1,2,3,4
+ *
+ *
+ *
+ * \snippet Vector/4_complex_prop/main_ser.cpp vect assign
+ *
+ */
+
+ //! \cond [vect assign] \endcond
+
+ auto it = vd.getDomainIterator();
+
+ while (it.isNext())
+ {
+ auto p = it.get();
+
+ // we define x, assign a random position between 0.0 and 1.0
+ vd.getPos(p)[0] = (float)rand() / RAND_MAX;
+
+ // we define y, assign a random position between 0.0 and 1.0
+ vd.getPos(p)[1] = (float)rand() / RAND_MAX;
+
+ // Get the particle position as point
+ Point<2,float> pt = vd.getPos(p);
+
+ // create a C string from the particle coordinates
+ // and copy into my struct
+ vd.getProp<my_s1>(p).size = 32;
+ vd.getProp<my_s1>(p).ptr = new char[32];
+ strcpy(vd.getProp<my_s1>(p).ptr,pt.toString().c_str());
+
+ // create a C++ string from the particle coordinates
+ vd.getProp<my_s1>(p).str = std::string(pt.toString());
+
+ vd.getProp<my_s1>(p).v.add(1);
+ vd.getProp<my_s1>(p).v.add(2);
+ vd.getProp<my_s1>(p).v.add(3);
+
+ pt = pt * 2.0;
+
+ // create a C string from the particle coordinates multiplied by 2.0
+ // and copy into my struct
+ vd.getProp<my_s2>(p).size = 32;
+ vd.getProp<my_s2>(p).ptr = new char[32];
+ strcpy(vd.getProp<my_s2>(p).ptr,pt.toString().c_str());
+
+ // create a C++ string from the particle coordinates
+ vd.getProp<my_s2>(p).str = std::string(pt.toString());
+
+ vd.getProp<my_s2>(p).v.add(1);
+ vd.getProp<my_s2>(p).v.add(2);
+ vd.getProp<my_s2>(p).v.add(3);
+ vd.getProp<my_s2>(p).v.add(4);
+
+ // next particle
+ ++it;
+ }
+
+ //! \cond [vect assign] \endcond
+
+ /*!
+ *
+ * \page Vector_4_complex_prop_ser Vector 4 property serialization
+ *
+ *
+ * ## Mapping and ghost_get ##
+ *
+ * Particles are redistributed across processors and we also synchronize the ghost
+ *
+ * \see \ref e0_s_map
+ *
+ * \see \ref e1_part_ghost
+ *
+ * \snippet Vector/4_complex_prop/main_ser.cpp vect map ghost
+ *
+ */
+
+ //! \cond [vect map ghost] \endcond
+
+ // Particles are redistribued across the processors
+ vd.map();
+
+ // Synchronize the ghost
+ vd.ghost_get<my_s1,my_s2>();
+
+ //! \cond [vect map ghost] \endcond
+
+ /*!
+ *
+ * \page Vector_4_complex_prop_ser Vector 4 property serialization
+ *
+ *
+ * ## Output and VTK visualization ##
+ *
+ * Vector with complex properties can be still be visualized, because unknown properties are
+ * automatically excluded
+ *
+ * \see \ref e0_s_vis_vtk
+ *
+ * \snippet Vector/4_complex_prop/main.cpp vtk
+ *
+ */
+
+ //! \cond [vtk] \endcond
+
+ vd.write("particles");
+
+ //! \cond [vtk] \endcond
+
+ /*!
+ *
+ * \page Vector_4_complex_prop_ser Vector 4 property serialization
+ *
+ * ## Print 4 particles in the ghost area ##
+ *
+ * Here we print that the first 4 particles to show that the two my_struct contain the
+ * right information
+ *
+ * \snippet Vector/4_complex_prop/main_ser.cpp print ghost info
+ *
+ */
+
+ //! \cond [print ghost info] \endcond
+
+ size_t fg = vd.size_local();
+
+ Vcluster & v_cl = create_vcluster();
+
+ // Only the master processor print
+ if (v_cl.getProcessUnitID() == 0)
+ {
+ // Print 4 particles
+ for ( ; fg < vd.size_local()+4 ; fg++)
+ {
+ // Print my struct1 information
+ std::cout << "my_struct1:" << std::endl;
+ std::cout << "C-string: " << vd.getProp<my_s1>(fg).ptr << std::endl;
+ std::cout << "Cpp-string: " << vd.getProp<my_s1>(fg).str << std::endl;
+
+ for (size_t i = 0 ; i < vd.getProp<my_s1>(fg).v.size() ; i++)
+ std::cout << "Element: " << i << " " << vd.getProp<my_s1>(fg).v.get(i) << std::endl;
+
+ // Print my struct 2 information
+ std::cout << "my_struct2" << std::endl;
+ std::cout << "C-string: " << vd.getProp<my_s2>(fg).ptr << std::endl;
+ std::cout << "Cpp-string: " << vd.getProp<my_s2>(fg).str << std::endl;
+
+ for (size_t i = 0 ; i < vd.getProp<my_s2>(fg).v.size() ; i++)
+ std::cout << "Element: " << i << " " << vd.getProp<my_s2>(fg).v.get(i) << std::endl;
+ }
+ }
+
+ //! \cond [print ghost info] \endcond
+
+ /*!
+ * \page Vector_4_complex_prop_ser Vector 4 property serialization
+ *
+ * ## Finalize ## {#finalize}
+ *
+ * At the very end of the program we have always to de-initialize the library
+ *
+ * \snippet Vector/4_complex_prop/main_ser.cpp finalize
+ *
+ */
+
+ //! \cond [finalize] \endcond
+
+ openfpm_finalize();
+
+ //! \cond [finalize] \endcond
+
+ /*!
+ * \page Vector_4_complex_prop_ser Vector 4 property serialization
+ *
+ * # Full code # {#code}
+ *
+ * \include Vector/4_complex_prop/main_ser.cpp
+ *
+ */
+}
diff --git a/example/Vector/4_reorder/main_expr.cpp b/example/Vector/4_reorder/main_expr.cpp
new file mode 100644
index 00000000..86c7e31b
--- /dev/null
+++ b/example/Vector/4_reorder/main_expr.cpp
@@ -0,0 +1,139 @@
+#include "Vector/vector_dist.hpp"
+#include "Plot/GoogleChart.hpp"
+#include "Operators/Vector/vector_dist_operators.hpp"
+
+constexpr int velocity = 0;
+constexpr int force = 1;
+
+struct ln_potential
+{
+ double sigma12,sigma6;
+
+ ln_potential(double sigma12_, double sigma6_) {sigma12 = sigma12_, sigma6 = sigma6_;}
+
+ Point<2,double> value(const Point<3,double> & xp, const Point<3,double> xq)
+ {
+ double rn = norm2(xp - xq);
+
+ Point<2,double> E({4.0 * ( sigma12 / (rn*rn*rn*rn*rn*rn) - sigma6 / ( rn*rn*rn) ),
+ 0.0});
+
+ return E;
+ }
+};
+
+struct ln_force
+{
+ double sigma12,sigma6;
+
+ ln_force(double sigma12_, double sigma6_) {sigma12 = sigma12_; sigma6 = sigma6_;}
+
+ Point<3,double> value(const Point<3,double> & xp, const Point<3,double> xq)
+ {
+ Point<3,double> r = xp - xq;
+ double rn = norm2(r);
+
+ return 24.0*(2.0 * sigma12 / (rn*rn*rn*rn*rn*rn*rn) - sigma6 / (rn*rn*rn*rn)) * r;
+ }
+};
+
+int main(int argc, char* argv[])
+{
+ double dt = 0.0005, sigma = 0.1, r_cut = 0.3;
+
+ double sigma6 = pow(sigma,6), sigma12 = pow(sigma,12);
+
+ openfpm::vector<double> x;
+ openfpm::vector<openfpm::vector<double>> y;
+
+
+ openfpm_init(&argc,&argv);
+ Vcluster & vcl = create_vcluster();
+
+ size_t sz[3] = {10,10,10};
+ Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
+ size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
+ Ghost<3,float> ghost(r_cut);
+
+ ln_force lf(sigma12,sigma6);
+ ln_potential lp(sigma12,sigma6);
+
+ vector_dist<3,double, aggregate<Point<3,double>,Point<3,double>> > vd(0,box,bc,ghost);
+
+ auto v_force = getV<force>(vd);
+ auto v_velocity = getV<velocity>(vd);
+ auto v_pos = getV<PROP_POS>(vd);
+
+ auto it = vd.getGridIterator(sz);
+
+ while (it.isNext())
+ {
+ vd.add();
+
+ auto key = it.get();
+
+ vd.getLastPos()[0] = key.get(0) * it.getSpacing(0);
+ vd.getLastPos()[1] = key.get(1) * it.getSpacing(1);
+ vd.getLastPos()[2] = key.get(2) * it.getSpacing(2);
+
+ ++it;
+ }
+
+ v_force = 0;
+ v_velocity = 0;
+
+ auto NN = vd.getCellList(r_cut);
+
+ vd.updateCellList(NN);
+ v_force = applyKernel_in_sim(vd,NN,lf);
+ unsigned long int f = 0;
+
+ // MD time stepping
+ for (size_t i = 0; i < 10000 ; i++)
+ {
+ assign(v_velocity, v_velocity + 0.5*dt*v_force,
+ v_pos, v_pos + v_velocity*dt);
+
+ vd.map();
+ vd.template ghost_get<>();
+
+ // calculate forces or a(tn + 1) Step 2
+ vd.updateCellList(NN);
+ v_force = applyKernel_in_sim(vd,NN,lf);
+
+ v_velocity = v_velocity + 0.5*dt*v_force;
+
+ if (i % 100 == 0)
+ {
+ vd.deleteGhost();
+ vd.write("particles_",f);
+ vd.ghost_get<>();
+
+ Point<2,double> E = rsum(applyKernel_in_sim(vd,NN,lp) + (v_velocity * v_velocity)/2.0,vd).get();
+
+ vcl.sum(E.get(0));vcl.sum(E.get(1));
+ vcl.execute();
+
+ // we save the energy calculated at time step i c contain the time-step y contain the energy
+ x.add(i);
+ y.add({E.get(0),E.get(1),E.get(0) - E.get(1)});
+
+ if (vcl.getProcessUnitID() == 0)
+ std::cout << "Energy Total: " << E.get(0) << " Kinetic: " << E.get(1) << " Potential: " << E.get(0) - E.get(1) << std::endl;
+
+ f++;
+ }
+ }
+
+ GCoptions options;
+ options.title = std::string("Energy with time");
+ options.yAxis = std::string("Energy");
+ options.xAxis = std::string("iteration");
+ options.lineWidth = 1.0;
+
+ GoogleChart cg;
+ cg.AddLinesGraph(x,y,options);
+ cg.write("gc_plot2_out.html");
+
+ openfpm_finalize();
+}
diff --git a/src/Vector/vector_dist_complex_prp_unit_test.hpp b/src/Vector/vector_dist_complex_prp_unit_test.hpp
new file mode 100644
index 00000000..253259a8
--- /dev/null
+++ b/src/Vector/vector_dist_complex_prp_unit_test.hpp
@@ -0,0 +1,231 @@
+/*
+ * vector_dist_complex_prp_unit_test.hpp
+ *
+ * Created on: Sep 18, 2016
+ * Author: i-bird
+ */
+
+#ifndef SRC_VECTOR_VECTOR_DIST_COMPLEX_PRP_UNIT_TEST_HPP_
+#define SRC_VECTOR_VECTOR_DIST_COMPLEX_PRP_UNIT_TEST_HPP_
+
+
+BOOST_AUTO_TEST_CASE( vector_dist_periodic_complex_prp_test_use_3d )
+{
+ Vcluster & v_cl = create_vcluster();
+
+ if (v_cl.getProcessingUnits() > 48)
+ return;
+
+ // set the seed
+ // create the random generator engine
+ std::srand(v_cl.getProcessUnitID());
+ std::default_random_engine eg;
+ std::uniform_real_distribution<float> ud(0.0f, 1.0f);
+
+ long int k = 124288 * v_cl.getProcessingUnits();
+
+ long int big_step = k / 4;
+ big_step = (big_step == 0)?1:big_step;
+
+ print_test_v( "Testing 3D vector full complex prp vector k<=",k);
+
+ // 3D test
+ for ( ; k >= 2 ; k-= decrement(k,big_step) )
+ {
+ BOOST_TEST_CHECKPOINT( "Testing 3D full complex prp vector k=" << k );
+
+ Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
+
+ // Boundary conditions
+ size_t bc[3]={NON_PERIODIC,NON_PERIODIC,NON_PERIODIC};
+
+ // factor
+ float factor = pow(create_vcluster().getProcessingUnits()/2.0f,1.0f/3.0f);
+
+ // ghost
+ Ghost<3,float> ghost(0.05 / factor);
+
+ // ghost2 (a little bigger because of round off error)
+ Ghost<3,float> ghost2(0.05001 / factor);
+
+ // Distributed vector
+ vector_dist<3,float, aggregate<openfpm::vector<float>,grid_cpu<3,aggregate<double,double[3]>> > > vd(k,box,bc,ghost);
+
+ auto it = vd.getIterator();
+
+ while (it.isNext())
+ {
+ auto key = it.get();
+
+ vd.getPos(key)[0] = ud(eg);
+ vd.getPos(key)[1] = ud(eg);
+ vd.getPos(key)[2] = ud(eg);
+
+ // initalize vector property and grid
+
+ vd.template getProp<0>(key).add(vd.getPos(key)[0]);
+ vd.template getProp<0>(key).add(vd.getPos(key)[1]);
+ vd.template getProp<0>(key).add(vd.getPos(key)[2]);
+
+ vd.template getProp<0>(key).add(vd.getPos(key)[0]+vd.getPos(key)[1]);
+ vd.template getProp<0>(key).add(vd.getPos(key)[1]+vd.getPos(key)[2]);
+ vd.template getProp<0>(key).add(vd.getPos(key)[0]+vd.getPos(key)[2]);
+
+ // Grid
+ size_t sz[] = {3,3,3};
+ vd.template getProp<1>(key).resize(sz);
+
+ vd.template getProp<1>(key).template get<0>({0,0,0}) = vd.getPos(key)[0];
+ vd.template getProp<1>(key).template get<0>({0,0,1}) = vd.getPos(key)[1];
+ vd.template getProp<1>(key).template get<0>({0,0,2}) = vd.getPos(key)[2];
+
+ vd.template getProp<1>(key).template get<0>({0,1,0}) = 2.0*vd.getPos(key)[0];
+ vd.template getProp<1>(key).template get<0>({0,1,1}) = 2.0*vd.getPos(key)[1];
+ vd.template getProp<1>(key).template get<0>({0,1,2}) = 2.0*vd.getPos(key)[2];
+
+ vd.template getProp<1>(key).template get<0>({0,2,0}) = 3.0*vd.getPos(key)[0];
+ vd.template getProp<1>(key).template get<0>({0,2,1}) = 3.0*vd.getPos(key)[1];
+ vd.template getProp<1>(key).template get<0>({0,2,2}) = 3.0*vd.getPos(key)[2];
+
+ vd.template getProp<1>(key).template get<0>({1,0,0}) = 4.0*vd.getPos(key)[0];
+ vd.template getProp<1>(key).template get<0>({1,0,1}) = 4.0*vd.getPos(key)[1];
+ vd.template getProp<1>(key).template get<0>({1,0,2}) = 4.0*vd.getPos(key)[2];
+
+ vd.template getProp<1>(key).template get<0>({1,1,0}) = 5.0*vd.getPos(key)[0];
+ vd.template getProp<1>(key).template get<0>({1,1,1}) = 5.0*vd.getPos(key)[1];
+ vd.template getProp<1>(key).template get<0>({1,1,2}) = 5.0*vd.getPos(key)[2];
+
+ vd.template getProp<1>(key).template get<0>({1,2,0}) = 6.0*vd.getPos(key)[0];
+ vd.template getProp<1>(key).template get<0>({1,2,1}) = 6.0*vd.getPos(key)[1];
+ vd.template getProp<1>(key).template get<0>({1,2,2}) = 6.0*vd.getPos(key)[2];
+
+ vd.template getProp<1>(key).template get<0>({2,0,0}) = 7.0*vd.getPos(key)[0];
+ vd.template getProp<1>(key).template get<0>({2,0,1}) = 7.0*vd.getPos(key)[1];
+ vd.template getProp<1>(key).template get<0>({2,0,2}) = 7.0*vd.getPos(key)[2];
+
+ vd.template getProp<1>(key).template get<0>({2,1,0}) = 8.0*vd.getPos(key)[0];
+ vd.template getProp<1>(key).template get<0>({2,1,1}) = 8.0*vd.getPos(key)[1];
+ vd.template getProp<1>(key).template get<0>({2,1,2}) = 8.0*vd.getPos(key)[2];
+
+ vd.template getProp<1>(key).template get<0>({2,2,0}) = 9.0*vd.getPos(key)[0];
+ vd.template getProp<1>(key).template get<0>({2,2,1}) = 9.0*vd.getPos(key)[1];
+ vd.template getProp<1>(key).template get<0>({2,2,2}) = 9.0*vd.getPos(key)[2];
+
+ ++it;
+ }
+
+ vd.map();
+
+ // sync the ghost
+ vd.ghost_get<0,1>();
+
+ // Domain + ghost
+ Box<3,float> dom_ext = box;
+ dom_ext.enlarge(ghost2);
+
+ // Iterate on all particles domain + ghost
+ size_t l_cnt = 0;
+ size_t nl_cnt = 0;
+ size_t n_out = 0;
+
+ auto it2 = vd.getIterator();
+ count_local_n_local<3,vector_dist<3,float, aggregate<openfpm::vector<float>,grid_cpu<3,aggregate<double,double[3]>>>> >(vd,it2,bc,box,dom_ext,l_cnt,nl_cnt,n_out);
+
+ // No particles should be out of domain + ghost
+ BOOST_REQUIRE_EQUAL(n_out,0ul);
+
+ // Ghost must populated because we synchronized them
+ if (k > 524288)
+ {
+ BOOST_REQUIRE(nl_cnt != 0);
+ BOOST_REQUIRE(l_cnt > nl_cnt);
+ }
+
+ // Sum all the particles inside the domain
+ v_cl.sum(l_cnt);
+ v_cl.execute();
+ BOOST_REQUIRE_EQUAL(l_cnt,(size_t)k);
+
+ l_cnt = 0;
+ nl_cnt = 0;
+
+ // Iterate only on the ghost particles
+ auto itg = vd.getGhostIterator();
+ count_local_n_local<3, vector_dist<3,float, aggregate<openfpm::vector<float>,grid_cpu<3,aggregate<double,double[3]>>>> >(vd,itg,bc,box,dom_ext,l_cnt,nl_cnt,n_out);
+
+ // No particle on the ghost must be inside the domain
+ BOOST_REQUIRE_EQUAL(l_cnt,0ul);
+
+ // Ghost must be populated
+ if (k > 524288)
+ {
+ BOOST_REQUIRE(nl_cnt != 0);
+ }
+
+ // check that the particles contain the correct information
+
+ bool ret = true;
+ auto it3 = vd.getDomainAndGhostIterator();
+
+ while (it3.isNext())
+ {
+ auto key = it3.get();
+
+ ret &= vd.template getProp<0>(key).size() == 6;
+
+ ret &= vd.template getProp<0>(key).get(0) == vd.getPos(key)[0];
+ ret &= vd.template getProp<0>(key).get(1) == vd.getPos(key)[1];
+ ret &= vd.template getProp<0>(key).get(2) == vd.getPos(key)[2];
+
+ ret &= vd.template getProp<0>(key).get(3) == vd.getPos(key)[0]+vd.getPos(key)[1];
+ ret &= vd.template getProp<0>(key).get(4) == vd.getPos(key)[1]+vd.getPos(key)[2];
+ ret &= vd.template getProp<0>(key).get(5) == vd.getPos(key)[0]+vd.getPos(key)[2];
+
+// BOOST_REQUIRE_EQUAL(vd.template getProp<0>(key).get(0),vd.getPos(key)[0]);
+
+ ret &= vd.template getProp<1>(key).size() == 3*3*3;
+ ret &= vd.template getProp<1>(key).template get<0>({0,0,0}) == vd.getPos(key)[0];
+ ret &= vd.template getProp<1>(key).template get<0>({0,0,1}) == vd.getPos(key)[1];
+ ret &= vd.template getProp<1>(key).template get<0>({0,0,2}) == vd.getPos(key)[2];
+
+ ret &= vd.template getProp<1>(key).template get<0>({0,1,0}) == 2.0*vd.getPos(key)[0];
+ ret &= vd.template getProp<1>(key).template get<0>({0,1,1}) == 2.0*vd.getPos(key)[1];
+ ret &= vd.template getProp<1>(key).template get<0>({0,1,2}) == 2.0*vd.getPos(key)[2];
+
+ ret &= vd.template getProp<1>(key).template get<0>({0,2,0}) == 3.0*vd.getPos(key)[0];
+ ret &= vd.template getProp<1>(key).template get<0>({0,2,1}) == 3.0*vd.getPos(key)[1];
+ ret &= vd.template getProp<1>(key).template get<0>({0,2,2}) == 3.0*vd.getPos(key)[2];
+
+ ret &= vd.template getProp<1>(key).template get<0>({1,0,0}) == 4.0*vd.getPos(key)[0];
+ ret &= vd.template getProp<1>(key).template get<0>({1,0,1}) == 4.0*vd.getPos(key)[1];
+ ret &= vd.template getProp<1>(key).template get<0>({1,0,2}) == 4.0*vd.getPos(key)[2];
+
+ ret &= vd.template getProp<1>(key).template get<0>({1,1,0}) == 5.0*vd.getPos(key)[0];
+ ret &= vd.template getProp<1>(key).template get<0>({1,1,1}) == 5.0*vd.getPos(key)[1];
+ ret &= vd.template getProp<1>(key).template get<0>({1,1,2}) == 5.0*vd.getPos(key)[2];
+
+ ret &= vd.template getProp<1>(key).template get<0>({1,2,0}) == 6.0*vd.getPos(key)[0];
+ ret &= vd.template getProp<1>(key).template get<0>({1,2,1}) == 6.0*vd.getPos(key)[1];
+ ret &= vd.template getProp<1>(key).template get<0>({1,2,2}) == 6.0*vd.getPos(key)[2];
+
+ ret &= vd.template getProp<1>(key).template get<0>({2,0,0}) == 7.0*vd.getPos(key)[0];
+ ret &= vd.template getProp<1>(key).template get<0>({2,0,1}) == 7.0*vd.getPos(key)[1];
+ ret &= vd.template getProp<1>(key).template get<0>({2,0,2}) == 7.0*vd.getPos(key)[2];
+
+ ret &= vd.template getProp<1>(key).template get<0>({2,1,0}) == 8.0*vd.getPos(key)[0];
+ ret &= vd.template getProp<1>(key).template get<0>({2,1,1}) == 8.0*vd.getPos(key)[1];
+ ret &= vd.template getProp<1>(key).template get<0>({2,1,2}) == 8.0*vd.getPos(key)[2];
+
+ ret &= vd.template getProp<1>(key).template get<0>({2,2,0}) == 9.0*vd.getPos(key)[0];
+ ret &= vd.template getProp<1>(key).template get<0>({2,2,1}) == 9.0*vd.getPos(key)[1];
+ ret &= vd.template getProp<1>(key).template get<0>({2,2,2}) == 9.0*vd.getPos(key)[2];
+
+ ++it3;
+ }
+
+ BOOST_REQUIRE_EQUAL(ret,true);
+ }
+}
+
+
+#endif /* SRC_VECTOR_VECTOR_DIST_COMPLEX_PRP_UNIT_TEST_HPP_ */
--
GitLab