diff --git a/example/Vector/0_simple/main.cpp b/example/Vector/0_simple/main.cpp
index 54b4fccc64e8b3cb1f00dae622b9f67ac3d363d9..badf7e47200278ed28f959528a860c87f22d8f34 100644
--- a/example/Vector/0_simple/main.cpp
+++ b/example/Vector/0_simple/main.cpp
@@ -7,6 +7,8 @@
* \subpage Vector_4_reo_root
* \subpage Vector_4_cp
* \subpage Vector_4_mp_cl
+ * \subpage Vector_5_md_vl_sym
+ * \subpage Vector_6_complex_usage
*
*/
diff --git a/example/Vector/3_molecular_dynamic/Makefile b/example/Vector/3_molecular_dynamic/Makefile
index 6d55806ec7a5643f4e95e499a145b277d5639604..d535a223f7b46c91e7f335459f98299debb699bb 100644
--- a/example/Vector/3_molecular_dynamic/Makefile
+++ b/example/Vector/3_molecular_dynamic/Makefile
@@ -7,9 +7,8 @@ LDIR =
OBJ = main.o
OBJ_EXPR = main_expr.o
OBJ_VL = main_vl.o
-OBJ_VL_SYM = main_vl_sym.o
-all: md_dyn md_dyn_expr md_dyn_vl md_dyn_vl_sym
+all: md_dyn md_dyn_expr md_dyn_vl
%.o: %.cpp
$(CC) -fext-numeric-literals -O3 -g -c --std=c++11 -o $@ $< $(INCLUDE_PATH)
@@ -23,14 +22,11 @@ md_dyn_expr: $(OBJ_EXPR)
md_dyn_vl: $(OBJ_VL)
$(CC) -o $@ $^ $(CFLAGS) $(LIBS_PATH) $(LIBS)
-md_dyn_vl_sym: $(OBJ_VL_SYM)
- $(CC) -o $@ $^ $(CFLAGS) $(LIBS_PATH) $(LIBS)
-
run: all
- source $$HOME/openfpm_vars; mpirun -np 3 ./md_dyn; mpirun -np 3 ./md_dyn_expr; mpirun -np 3 ./md_dyn_vl; mpirun -np 3 ./md_dyn_vl_sym;
+ source $$HOME/openfpm_vars; mpirun -np 3 ./md_dyn; mpirun -np 3 ./md_dyn_expr; mpirun -np 3 ./md_dyn_vl;
.PHONY: clean all run
clean:
- rm -f *.o *~ core md_dyn md_dyn_expr md_dyn_vl md_dyn_vl_sym
+ rm -f *.o *~ core md_dyn md_dyn_expr md_dyn_vl
diff --git a/example/Vector/3_molecular_dynamic/main_vl_sym.cpp b/example/Vector/3_molecular_dynamic/main_vl_sym.cpp
deleted file mode 100644
index 07136e19c202c479ad2559567ce58b1ea5fb419f..0000000000000000000000000000000000000000
--- a/example/Vector/3_molecular_dynamic/main_vl_sym.cpp
+++ /dev/null
@@ -1,552 +0,0 @@
-#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/6_complex_usage/Makefile b/example/Vector/6_complex_usage/Makefile
new file mode 100644
index 0000000000000000000000000000000000000000..98a3b9c708a9c4f1d4a3ecf84d83ca2caa43c686
--- /dev/null
+++ b/example/Vector/6_complex_usage/Makefile
@@ -0,0 +1,24 @@
+include ../../example.mk
+
+CC=mpic++
+
+LDIR =
+
+OBJ_DORD = main.o
+
+all: complex_use
+
+%.o: %.cpp
+ $(CC) -fext-numeric-literals -O3 -g -c --std=c++11 $(OPT) -o $@ $< $(INCLUDE_PATH)
+
+complex_use: $(OBJ_DORD)
+ $(CC) -o $@ $^ $(CFLAGS) $(LIBS_PATH) $(LIBS)
+
+run: all_test
+ source $$HOME/openfpm_vars; mpirun -np 3 ./complex_use
+
+.PHONY: clean all run all_test on_test
+
+clean:
+ rm -f *.o *~ core complex_use
+
diff --git a/example/Vector/6_complex_usage/main.cpp b/example/Vector/6_complex_usage/main.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ecbd4b9a07f09ea12cab26888bdc1fdba1136255
--- /dev/null
+++ b/example/Vector/6_complex_usage/main.cpp
@@ -0,0 +1,409 @@
+#include "Vector/vector_dist.hpp"
+#include "data_type/aggregate.hpp"
+#include "Plot/GoogleChart.hpp"
+#include "Plot/util.hpp"
+#include "timer.hpp"
+
+/*!
+ * \page Vector_6_complex_usage Vector 6 complex usage for validation and debugging
+ *
+ * [TOC]
+ *
+ * # Complex usage for validation and debugging # {#e6_cp_deb}
+ *
+ * In this example we show how the flexibility of the library can be used to perform complex
+ * tasks for validation and debugging. We will use a lot of feature that has been explained in
+ * the previous examples
+ *
+ * In a previous example we show in case of symmetric interaction between particles how to symmetric
+ * cell list could be used to speed up the calculation.
+ *
+ * \see not present
+ *
+ * In this example we validate that both computation match, in particular because the computation depend
+ * from the neighborhood, we will check and validate that the neighborhood of each
+ * particle is equivalent with symmetric cell list and normal cell-list.
+ *
+ */
+
+int main(int argc, char* argv[])
+{
+ /*!
+ * \page Vector_6_complex_usage Vector 6 complex usage for validation
+ *
+ * ## Initialization ##
+ *
+ * The initialization is classical we define cutoff radius, domain, boundary conditions,
+ * ghost part and some useful constants. We also define a struct that will contain the debug information. This structure
+ * contain an id that is the global id of the particle and a point that is the the position of the
+ * neighborhood.
+ *
+ * \see \ref e0_s_init
+ *
+ * \snippet Vector/6_complex_usage/main.cpp initialization
+ *
+ */
+
+ //! \cond [initialization] \endcond
+
+ openfpm_init(&argc,&argv);
+
+ constexpr int gid = 0;
+ constexpr int nn_norm = 1;
+ constexpr int nn_sym = 2;
+
+ // used to define the domain
+ float L = 1000.0;
+
+ // Domain
+ Box<3,float> box({-L,-L,-L},{L,L,L});
+
+ // Boundary conditions
+ size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
+
+ // cut-off radius
+ float r_cut = 100.0;
+
+ // ghost
+ Ghost<3,float> ghost(r_cut);
+
+ // Point and global id
+ struct point_and_gid
+ {
+ // global id of the particle
+ size_t id;
+
+ // Position of the neighborhood particle
+ Point<3,float> xq;
+
+ // Used to reorder the neighborhood particles by id
+ bool operator<(const struct point_and_gid & pag)
+ {
+ return (id < pag.id);
+ }
+ };
+
+ //! \cond [initialization] \endcond
+
+ /*!
+ * \page Vector_6_complex_usage Vector 6 complex usage for debugging
+ *
+ * ## Particle ##
+ *
+ * For this example we will use a particle with 3 properties
+ *
+ * * The first property is a global index for the particle unique across processors
+ * * The second is a vector that contain for each neighborhood the global id of the particle
+ * and its position
+ * * The last one is again a vector equivalent to the previous but is produced with the symmetric
+ * cell-list
+ *
+ * \snippet Vector/6_complex_usage/main.cpp particle prop
+ *
+ */
+
+ //! \cond [particle prop] \endcond
+
+ // Particle properties list
+ typedef aggregate<size_t,openfpm::vector<point_and_gid>,openfpm::vector<point_and_gid>> part_prop;
+
+ //! \cond [particle prop] \endcond
+
+ /*!
+ * \page Vector_6_complex_usage Vector 6 complex usage for debugging
+ *
+ * ## Distributed vector ##
+ *
+ * Here we create a distributed vector of 4096 particles. We initialize the particles
+ * randomly and we assign a global index to the first property of the particle
+ *
+ * The function accum return the following number
+ *
+ * \f$ \sum_{i < proc} size\_local(i) \f$
+ *
+ * where \f$ proc \f$ is the processor where we are computing the formula and \f$ size\_local(i) \f$
+ * is the number of particles the processor \f$ i \f$ has. This number is used to
+ * produce a global id of the particles
+ *
+ * \snippet Vector/6_complex_usage/main.cpp glob id part
+ *
+ */
+
+ //! \cond [glob id part] \endcond
+
+ // Distributed vector
+ vector_dist<3,float, part_prop > vd(4096,box,bc,ghost);
+
+ // used to calculate a global index for each particle
+ size_t start = vd.accum();
+
+ auto it = vd.getDomainIterator();
+
+ while (it.isNext())
+ {
+ auto key = it.get();
+
+ vd.getPos(key)[0] = 2.0*L*((float)rand()/RAND_MAX) - L;
+ vd.getPos(key)[1] = 2.0*L*((float)rand()/RAND_MAX) - L;
+ vd.getPos(key)[2] = 2.0*L*((float)rand()/RAND_MAX) - L;
+
+ vd.getProp<gid>(key) = key.getKey() + start;
+
+ ++it;
+ }
+
+ //! \cond [glob id part] \endcond
+
+ /*!
+ * \page Vector_6_complex_usage Vector 6 complex usage for debugging
+ *
+ * ## Redistribute ##
+ *
+ * Redistribute the particles and synchronize the ghosts. In this case we
+ * are only interested in synchronizing the global id property.
+ *
+ * \snippet Vector/6_complex_usage/main.cpp map and ghost
+ *
+ */
+
+ //! \cond [map and ghost] \endcond
+
+ vd.map();
+
+ // sync the ghost
+ vd.ghost_get<gid>();
+
+ //! \cond [map and ghost] \endcond
+
+ /*!
+ * \page Vector_6_complex_usage Vector 6 complex usage for debugging
+ *
+ * ## Calculate the neighborhood of each particle ##
+ *
+ * Here we calculate the neighborhood of each particles using a simple **cell list**.
+ * In case the particle q has a distance from the particle p smaller than r_cut.
+ * We add it in the first list of neighborhood
+ *
+ * \snippet Vector/6_complex_usage/main.cpp add nn particles
+ *
+ */
+
+ //! \cond [add nn particles] \endcond
+
+ auto NN = vd.getCellList(r_cut);
+ auto p_it = vd.getDomainIterator();
+
+ while (p_it.isNext())
+ {
+ auto p = p_it.get();
+
+ Point<3,float> xp = vd.getPos(p);
+
+ auto Np = NN.getNNIterator(NN.getCell(vd.getPos(p)));
+
+ while (Np.isNext())
+ {
+ auto q = Np.get();
+
+ // skip self interaction
+ if (p.getKey() == q)
+ {
+ ++Np;
+ continue;
+ }
+
+ Point<3,float> xq = vd.getPos(q);
+ Point<3,float> f = (xp - xq);
+
+ float distance = f.norm();
+
+ // if the distance smalle than the cut-off radius add it to the neighborhood list
+ if (distance < r_cut )
+ {
+ vd.getProp<nn_norm>(p).add();
+ vd.getProp<nn_norm>(p).last().xq = xq;
+ vd.getProp<nn_norm>(p).last().id = vd.getProp<0>(q);
+ }
+
+ // Next neighborhood
+ ++Np;
+ }
+
+ // Next particle
+ ++p_it;
+ }
+
+ //! \cond [add nn particles] \endcond
+
+ /*!
+ * \page Vector_6_complex_usage Vector 6 complex usage for debugging
+ *
+ * ## Calculate the neighborhood of each particle with symmetric cell-list ##
+ *
+ * Here we calculate the neighborhood of each particle using instead a
+ * **symmetric cell list**. In case of symmetric cell-list if we find that a
+ * particle p is neighborhood of particle q we have to add p to the neighborhood
+ * of q and q to the neighborhood of p. Because q can be a ghost particle when
+ * we add the neighborhood p to q we have to transmit such information to the real
+ * owner of the particle. This can be done with the function **ghost_put**. In this
+ * case we use the operation **merge_** that add the already calculated neighborhood
+ * with the transmitted one. More in general it merge the information together.
+ *
+ * \snippet Vector/6_complex_usage/main.cpp add nn particles sym
+ *
+ */
+
+ //! \cond [add nn particles sym] \endcond
+
+ auto NN2 = vd.getCellListSym(r_cut);
+
+ auto p_it2 = vd.getDomainIterator();
+
+ while (p_it2.isNext())
+ {
+ auto p = p_it2.get();
+
+ Point<3,float> xp = vd.getPos(p);
+
+ auto Np = NN2.template getNNIteratorSym<NO_CHECK>(NN2.getCell(vd.getPos(p)),p.getKey(),vd.getPosVector());
+
+ while (Np.isNext())
+ {
+ auto q = Np.get();
+
+ if (p.getKey() == q)
+ {
+ ++Np;
+ continue;
+ }
+
+ // repulsive
+
+ Point<3,float> xq = vd.getPos(q);
+ Point<3,float> f = (xp - xq);
+
+ float distance = f.norm();
+
+ // Particle should be inside r_cut range
+
+ if (distance < r_cut )
+ {
+ vd.getProp<nn_sym>(p).add();
+ vd.getProp<nn_sym>(q).add();
+
+ vd.getProp<nn_sym>(p).last().xq = xq;
+ vd.getProp<nn_sym>(q).last().xq = xp;
+ vd.getProp<nn_sym>(p).last().id = vd.getProp<0>(q);
+ vd.getProp<nn_sym>(q).last().id = vd.getProp<0>(p);
+ }
+
+ ++Np;
+ }
+
+ ++p_it2;
+ }
+
+ vd.ghost_put<merge_,2>();
+
+ //! \cond [add nn particles sym] \endcond
+
+ /*!
+ * \page Vector_6_complex_usage Vector 6 complex usage for debugging
+ *
+ * ## Cheking for validation ##
+ *
+ * Here we check that the two calculated neighborhood match. In particular,
+ * because the order of the particles does not match, we have to first reorder
+ * by global-id, and than check the list. We cannot instead easily compare the
+ * position. The reason can be seen in this figure
+ *
+ \verbatim
+
+ +----------------+
+ | |
+ | 1 |
+ | * 2 |
+ | 3 * |
+ +<----- ghost * +<-------- real
+ |* <--- real |* <------- ghost
+ |4 |4
+ | |
+ | |
+ +----------------+
+
+
+
+ \endverbatim
+ *
+ * In the case we are calculating the neighborhood of the particle \f$ + \f$ in case of normal Cell-list.
+ * The real particle 1,2,3 are added, while the particle 4 is added with the ghost particle coordinates.
+ * In the case of symmetric cell-list the real 1,2,3 are added to \f$ + \f$. The particle 4
+ * instead is added by the ghost put. In particular 4 real add the particle to the \f$ + \f$ ghost particle
+ * and than ghost put merge the information. This mean that the symmetric add the particle 4 with the real coordinates
+ * of the particle 4
+ *
+ * \snippet Vector/6_complex_usage/main.cpp checking
+ *
+ */
+
+ //! \cond [checking] \endcond
+
+ auto p_it3 = vd.getDomainIterator();
+
+ bool ret = true;
+ while (p_it3.isNext())
+ {
+ auto p = p_it3.get();
+
+ vd.getProp<nn_norm>(p).sort();
+ vd.getProp<nn_sym>(p).sort();
+
+ ret &= vd.getProp<nn_norm>(p).size() == vd.getProp<nn_sym>(p).size();
+
+ for (size_t i = 0 ; i < vd.getProp<1>(p).size() ; i++)
+ {
+ ret &= vd.getProp<nn_norm>(p).get(i).id == vd.getProp<nn_sym>(p).get(i).id;
+
+ if (box.isInside(vd.getProp<nn_norm>(p).get(i).xq) == true)
+ {
+ ret &= vd.getProp<nn_norm>(p).get(i).xq == vd.getProp<nn_sym>(p).get(i).xq;
+ }
+ }
+
+ ++p_it3;
+ }
+
+ if (ret != true)
+ {
+ std::cout << "ERROR" << std::endl;
+ exit(1);
+ }
+
+ //! \cond [checking] \endcond
+
+ /*!
+ * \page Vector_6_complex_usage Vector 6 complex usage for validation and debugging
+ *
+ * ## Finalize ##
+ *
+ * At the very end of the program we have always to de-initialize the library
+ *
+ * \snippet Vector/6_complex_usage/main.cpp finalize
+ *
+ */
+
+ //! \cond [finalize] \endcond
+
+ openfpm_finalize();
+
+ //! \cond [finalize] \endcond
+
+ /*!
+ * \page Vector_6_complex_usage Vector 6 complex usage for validation and debugging
+ *
+ * # Full code # {#code}
+ *
+ * \include Vector/6_complex_usage/main.cpp
+ *
+ */
+}
diff --git a/openfpm_data b/openfpm_data
index 21a5f6ee7bf426ad13106a26f581b788b8c7c95d..1a4dc99034168bd02982055185682b3c3a9c2ff4 160000
--- a/openfpm_data
+++ b/openfpm_data
@@ -1 +1 @@
-Subproject commit 21a5f6ee7bf426ad13106a26f581b788b8c7c95d
+Subproject commit 1a4dc99034168bd02982055185682b3c3a9c2ff4
diff --git a/openfpm_io b/openfpm_io
index 62b8316404bf98c5cc5ad299563425656cfd3003..5cf769662d6e57a7a2fd12d1cfc835fd2793bc7b 160000
--- a/openfpm_io
+++ b/openfpm_io
@@ -1 +1 @@
-Subproject commit 62b8316404bf98c5cc5ad299563425656cfd3003
+Subproject commit 5cf769662d6e57a7a2fd12d1cfc835fd2793bc7b
diff --git a/openfpm_pdata.doc b/openfpm_pdata.doc
index 94101cdb7c8a63ebe6e1be42f6608f00bf1a879b..1375042ee7fb74340aa4414415ecafd645f56e6d 100644
--- a/openfpm_pdata.doc
+++ b/openfpm_pdata.doc
@@ -38,7 +38,7 @@ PROJECT_NAME = "OpenFPM_pdata"
# could be handy for archiving the generated documentation or if some version
# control system is used.
-PROJECT_NUMBER = 0.5.1
+PROJECT_NUMBER = 0.6.0
# Using the PROJECT_BRIEF tag one can provide an optional one line description
# for a project that appears at the top of each page and should give viewer a
diff --git a/script/install_PETSC.sh b/script/install_PETSC.sh
index a3bbfb8f436034836b63a5243f7324495c17c842..9341cb8a7ca54689009818169687a7be0396ba5b 100755
--- a/script/install_PETSC.sh
+++ b/script/install_PETSC.sh
@@ -181,6 +181,9 @@ if [ ! -d "$1/MUMPS" ]; then
$sed_command -i "/LIBBLAS\s=\s-lblas/c\LIBBLAS = -lopenblas" Makefile.inc
+ $sed_command -i "/INCPAR\s\+=\s\-I\/usr\/include/c\INCPAR =" Makefile.inc
+ $sed_command -i "/LIBPAR\s\+=\s\$(SCALAP)\s\-L\/usr\/lib\s\-lmpi/c\LIBPAR = \$(SCALAP)" Makefile.inc
+
make -j $2
if [ $? -eq 0 ]; then
diff --git a/src/Vector/vector_dist.hpp b/src/Vector/vector_dist.hpp
index 6606fdf33885cfcc00fc13edc4b542a03a17e36f..e8e0832b1d211d350daf2a9a687f7d00ddc539d6 100644
--- a/src/Vector/vector_dist.hpp
+++ b/src/Vector/vector_dist.hpp
@@ -526,15 +526,29 @@ public:
return cell_list;
}
+ /*! \brief for each particle get the symmetric verlet list
+ *
+ * \param r_cut cut-off radius
+ *
+ */
+ VerletList<dim,St,FAST,shift<dim,St> > getVerletSym(St r_cut)
+ {
+ VerletList<dim,St,FAST,shift<dim,St>> ver;
+ // Processor bounding box
+ Box<dim, St> pbox = getDecomposition().getProcessorBounds();
+
+ ver.InitializeSym(getDecomposition().getDomain(),pbox,getDecomposition().getGhost(),r_cut,v_pos,g_m);
+
+ return ver;
+ }
/*! \brief for each particle get the verlet list
*
* \param r_cut cut-off radius
- * \param opt option like VL_SYMMETRIC and VL_NON_SYMMETRIC
*
*/
- VerletList<dim,St,FAST,shift<dim,St> > getVerlet(St r_cut, size_t opt = VL_NON_SYMMETRIC)
+ VerletList<dim,St,FAST,shift<dim,St> > getVerlet(St r_cut)
{
VerletList<dim,St,FAST,shift<dim,St>> ver;
@@ -548,7 +562,7 @@ public:
// enlarge the box where the Verlet is defined
bt.enlarge(g);
- ver.Initialize(bt,getDecomposition().getProcessorBounds(),r_cut,v_pos,g_m,opt);
+ ver.Initialize(bt,getDecomposition().getProcessorBounds(),r_cut,v_pos,g_m,VL_NON_SYMMETRIC);
return ver;
}
@@ -1126,6 +1140,28 @@ public:
{
return v_pos;
}
+
+ /*! \brief It return the sum of the particles in the previous processors
+ *
+ * \return the particles number
+ *
+ */
+ size_t accum()
+ {
+ openfpm::vector<size_t> accu;
+
+ size_t sz = size_local();
+
+ v_cl.allGather(sz,accu);
+ v_cl.execute();
+
+ sz = 0;
+
+ for (size_t i = 0 ; i < v_cl.getProcessUnitID() ; i++)
+ sz += accu.get(i);
+
+ return sz;
+ }
};
diff --git a/src/Vector/vector_dist_cell_list_tests.hpp b/src/Vector/vector_dist_cell_list_tests.hpp
index 3826bbb9b74c49fb69a3cad570434fe8c682c5de..b6c4265a4729be1f64cf4ae2f776acbab0730d53 100644
--- a/src/Vector/vector_dist_cell_list_tests.hpp
+++ b/src/Vector/vector_dist_cell_list_tests.hpp
@@ -530,173 +530,206 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_cell_list )
BOOST_REQUIRE_EQUAL(ret,true);
}
-BOOST_AUTO_TEST_CASE( vector_dist_sym_verlet_list_test )
+
+BOOST_AUTO_TEST_CASE( vector_dist_symmetric_verlet_list )
{
- if (create_vcluster().getProcessingUnits() > 48)
+ Vcluster & v_cl = create_vcluster();
+
+ if (v_cl.getProcessingUnits() > 24)
return;
- long int k = 4096*create_vcluster().getProcessingUnits();
+ float L = 1000.0;
- long int big_step = k / 30;
+ // set the seed
+ // create the random generator engine
+ std::srand(0);
+ std::default_random_engine eg;
+ std::uniform_real_distribution<float> ud(-L,L);
+
+ long int k = 4096 * v_cl.getProcessingUnits();
+
+ long int big_step = k / 4;
big_step = (big_step == 0)?1:big_step;
- long int small_step = 21;
- print_test( "Testing vector symmetric verlet-list k<=",k);
+ print_test("Testing 3D periodic vector symmetric cell-list k=",k);
+ BOOST_TEST_CHECKPOINT( "Testing 3D periodic vector symmetric cell-list k=" << k );
- // 3D test
- for ( ; k > 8*big_step ; k-= (k > 2*big_step)?big_step:small_step )
- {
- double r_cut = 0.1;
+ Box<3,float> box({-L,-L,-L},{L,L,L});
- // domain
- Box<3,double> box({0.0,0.0,0.0},{1.0,1.0,1.0});
+ // Boundary conditions
+ size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
- // Boundary conditions
- size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
+ float r_cut = 100.0;
- // ghost, big enough to contain the interaction radius
- Ghost<3,float> ghost(r_cut);
+ // ghost
+ Ghost<3,float> ghost(r_cut);
- vector_dist<3,double, aggregate<double> > vd(k,box,bc,ghost);
+ // Point and global id
+ struct point_and_gid
+ {
+ size_t id;
+ Point<3,float> xq;
+ bool operator<(const struct point_and_gid & pag)
{
- auto it = vd.getDomainIterator();
+ return (id < pag.id);
+ }
+ };
- while (it.isNext())
- {
- auto p = it.get();
+ typedef aggregate<size_t,size_t,size_t,openfpm::vector<point_and_gid>,openfpm::vector<point_and_gid>> part_prop;
- vd.getPos(p)[0] = (double)rand()/RAND_MAX;
- vd.getPos(p)[1] = (double)rand()/RAND_MAX;
- vd.getPos(p)[2] = (double)rand()/RAND_MAX;
+ // Distributed vector
+ vector_dist<3,float, part_prop > vd(k,box,bc,ghost);
+ size_t start = vd.init_size_accum(k);
- ++it;
- }
- }
+ auto it = vd.getIterator();
- vd.map();
- vd.ghost_get<>();
+ while (it.isNext())
+ {
+ auto key = it.get();
- // Get the Cell list structure
- auto NN = vd.getVerlet(r_cut,VL_SYMMETRIC);
- auto NNc = vd.getVerlet(r_cut);
+ vd.getPos(key)[0] = ud(eg);
+ vd.getPos(key)[1] = ud(eg);
+ vd.getPos(key)[2] = ud(eg);
- // Get an iterator over particles
- auto it2 = vd.getDomainAndGhostIterator();
+ // Fill some properties randomly
- openfpm::vector<openfpm::vector<size_t>> idx;
- idx.resize(vd.size_local_with_ghost());
+ vd.getProp<0>(key) = 0;
+ vd.getProp<1>(key) = 0;
+ vd.getProp<2>(key) = key.getKey() + start;
- /////////// SYMMETRIC CASE VERLET-LIST ////////
+ ++it;
+ }
- // For each particle ...
- while (it2.isNext())
- {
- // ... p
- auto p = it2.get();
+ vd.map();
- // Get the position of the particle p
- Point<3,double> xp = vd.getPos(p);
+ // sync the ghost
+ vd.ghost_get<0,2>();
- // Get an iterator over the neighborhood of the particle p symmetric
- auto NpSym = NN.template getNNIterator<NO_CHECK>(p.getKey());
+ auto NN = vd.getVerlet(r_cut);
+ auto p_it = vd.getDomainIterator();
- // For each neighborhood of the particle p
- while (NpSym.isNext())
- {
- // Neighborhood particle q
- auto q = NpSym.get();
+ while (p_it.isNext())
+ {
+ auto p = p_it.get();
- // if p == q skip this particle
- if (q == p.getKey()) {++NpSym; continue;};
+ Point<3,float> xp = vd.getPos(p);
- // Get position of the particle q
- Point<3,double> xq = vd.getPos(q);
+ auto Np = NN.getNNIterator(p.getKey());
- // take the normalized direction
- double rn = norm2(xp - xq);
+ while (Np.isNext())
+ {
+ auto q = Np.get();
- // potential energy (using pow is slower)
- vd.getProp<0>(p) += rn;
- vd.getProp<0>(q) += rn;
+ if (p.getKey() == q)
+ {
+ ++Np;
+ continue;
+ }
- idx.get(p.getKey()).add(q);
- idx.get(q).add(p.getKey());
+ // repulsive
+ Point<3,float> xq = vd.getPos(q);
+ Point<3,float> f = (xp - xq);
+
+ float distance = f.norm();
+
+ // Particle should be inside 2 * r_cut range
- // Next neighborhood
- ++NpSym;
+ if (distance < r_cut )
+ {
+ vd.getProp<0>(p)++;
+ vd.getProp<3>(p).add();
+ vd.getProp<3>(p).last().xq = xq;
+ vd.getProp<3>(p).last().id = vd.getProp<2>(q);
}
- // Next Particle
- ++it2;
+ ++Np;
}
- /////////////// NON SYMMETRIC CASE VERLET-LIST ////////////////////////
+ ++p_it;
+ }
- openfpm::vector<openfpm::vector<size_t>> idx2;
- idx2.resize(vd.size_local());
+ // We now try symmetric Cell-list
- auto it = vd.getDomainIterator();
+ auto NN2 = vd.getVerletSym(r_cut);
- // For each particle ...
- while (it.isNext())
- {
- // ... p
- auto p = it.get();
+ auto p_it2 = vd.getDomainIterator();
+
+ while (p_it2.isNext())
+ {
+ auto p = p_it2.get();
- // Get the position of the particle p
- Point<3,double> xp = vd.getPos(p);
+ Point<3,float> xp = vd.getPos(p);
- // Get an iterator over the neighborhood of the particle p
- auto Np = NNc.template getNNIterator<NO_CHECK>(p.getKey());
+ auto Np = NN2.template getNNIterator<NO_CHECK>(p.getKey());
- double Ep = 0.0;
+ while (Np.isNext())
+ {
+ auto q = Np.get();
- // For each neighborhood of the particle p
- while (Np.isNext())
+ if (p.getKey() == q)
{
- // Neighborhood particle q
- auto q = Np.get();
+ ++Np;
+ continue;
+ }
- // if p == q skip this particle
- if (q == p.getKey()) {++Np; continue;};
+ // repulsive
- // Get position of the particle q
- Point<3,double> xq = vd.getPos(q);
+ Point<3,float> xq = vd.getPos(q);
+ Point<3,float> f = (xp - xq);
- // take the normalized direction
- double rn = norm2(xp - xq);
+ float distance = f.norm();
- idx2.get(p.getKey()).add(q);
+ // Particle should be inside r_cut range
- // potential energy (using pow is slower)
- Ep += rn;
+ if (distance < r_cut )
+ {
+ vd.getProp<1>(p)++;
+ vd.getProp<1>(q)++;
- // Next neighborhood
- ++Np;
+ vd.getProp<4>(p).add();
+ vd.getProp<4>(q).add();
+
+ vd.getProp<4>(p).last().xq = xq;
+ vd.getProp<4>(q).last().xq = xp;
+ vd.getProp<4>(p).last().id = vd.getProp<2>(q);
+ vd.getProp<4>(q).last().id = vd.getProp<2>(p);
}
- idx.get(p.getKey()).sort();
- idx2.get(p.getKey()).sort();
+ ++Np;
+ }
- bool ret = true;
+ ++p_it2;
+ }
- BOOST_REQUIRE_EQUAL(idx.get(p.getKey()).size(),idx2.get(p.getKey()).size());
+ vd.ghost_put<add_,1>();
+ vd.ghost_put<merge_,4>();
- for (size_t i = 0 ; i < idx.get(p.getKey()).size() ; i++)
- {
- ret &= idx.get(p.getKey()).get(i) == idx2.get(p.getKey()).get(i);
- }
+ auto p_it3 = vd.getDomainIterator();
- BOOST_REQUIRE_EQUAL(ret,true);
+ bool ret = true;
+ while (p_it3.isNext())
+ {
+ auto p = p_it3.get();
- BOOST_REQUIRE_CLOSE(Ep,vd.getProp<0>(p),0.01);
+ ret &= vd.getProp<1>(p) == vd.getProp<0>(p);
- // Next Particle
- ++it;
- }
+ vd.getProp<3>(p).sort();
+ vd.getProp<4>(p).sort();
+
+ ret &= vd.getProp<3>(p).size() == vd.getProp<4>(p).size();
+
+ for (size_t i = 0 ; i < vd.getProp<3>(p).size() ; i++)
+ ret &= vd.getProp<3>(p).get(i).id == vd.getProp<4>(p).get(i).id;
+
+ if (ret == false)
+ break;
+
+ ++p_it3;
}
+
+ BOOST_REQUIRE_EQUAL(ret,true);
}
diff --git a/src/Vector/vector_dist_comm.hpp b/src/Vector/vector_dist_comm.hpp
index 6a42f4c4658a0dd4ddd987c1d74744e54a1764be..8283f30c53add880daf3656a38eb42ea86a73f89 100644
--- a/src/Vector/vector_dist_comm.hpp
+++ b/src/Vector/vector_dist_comm.hpp
@@ -1052,12 +1052,13 @@ public:
size_t i2 = 0;
- #ifdef SE_CLASS1
if (v_prp.size() - lg_m != o_part_loc.size())
- std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " v_prp.size() - lg_m = " << v_prp.size() - lg_m << " != " << o_part_loc.size() << std::endl;
+ {
+ std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " Local ghost particles = " << v_prp.size() - lg_m << " != " << o_part_loc.size() << std::endl;
+ std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << "Check that you did a ghost_get before a ghost_put" << std::endl;
+ }
- #endif
for (size_t i = lg_m ; i < v_prp.size() ; i++)
{
diff --git a/src/main.cpp b/src/main.cpp
index 38a88d937269940f14d9282b2a8c655dda788382..d82bfd852c95a5a8526c338fe9e33cc3130a572a 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -22,6 +22,7 @@ int main(int argc, char* argv[])
return boost::unit_test::unit_test_main( &init_unit_test, argc, argv );
}
+#include "debug.hpp"
#include "Grid/grid_dist_id.hpp"
#include "Point_test.hpp"
#include "Decomposition/CartDecomposition.hpp"