Adding missing example files

example/Vector/8_DEM/Makefile

+include ../../example.mk
+
+CC=mpic++
+
+LDIR =
+
+OBJ = main.o
+
+OPT=
+
+all: dem
+
+dem_test: OPT += -DTEST_RUN
+dem_test: all
+
+%.o: %.cpp
+	$(CC) -O3 $(OPT) -g -c --std=c++11 -o $@ $< $(INCLUDE_PATH)
+
+dem: $(OBJ)
+	$(CC) -o $@ $^ $(OPT) $(CFLAGS) $(LIBS_PATH) $(LIBS)
+
+
+run: dem_test
+	mpirun -np 2 ./dem
+
+.PHONY: clean all run
+
+clean:
+	rm -f *.o *~ core dem
+

example/Vector/8_DEM/config.cfg

+[pack]
+files = main.cpp Makefile

example/Vector/8_DEM/main.cpp

+#define CHECKFOR_POSNAN
+#define CHECKFOR_PROPNAN
+
+/*!
+ * \page Vector_8_DEM Vector 8 Discrete element method
+ *
+ *
+ * [TOC]
+ *
+ * ## Introduction {#Dem_introduction}
+ *
+ * In this example we show how to implement a simulation Discrete element simulation (DEM)
+ * Using the Lorenz-force contact model. We will also use dynamic load balancing to keep the
+ * simulation balanced during all simulation time. The simulation result is show in the
+ * figures and video below
+ *
+ * \htmlonly
+ * <a href="#" onclick="hide_show('vector-video-3')" >Simulation video 1</a><br>
+ * <div style="display:none" id="vector-video-3">
+ * <video id="vid3" width="1200" height="576" controls> <source src="http://openfpm.mpi-cbg.de/web/images/examples/8_DEM/DEM_velocity.mp4" type="video/mp4"></video>
+ * </div>
+ * \endhtmlonly
+ *
+ * \htmlonly
+ * <img src="http://ppmcore.mpi-cbg.de/web/images/examples/8_DEM/DEM_60.png"/>
+ * \endhtmlonly
+ *
+ * A classical model for DEM simulations of spherical granular flows is the Silbert model,
+ * it includes a Herzian contact force and an elastic deformation of the grains. Each particles
+ *  has radius \f$R\f$, mass \f$m\f$, polar momentum \f$I\f$ and is represented by the location of its center of mass \f$r_{i}\f$.
+ * When two particles i and j collide or are in contact, the elastic contact deformation is given by:
+ * \f$\delta_{ij} = 2R-r_{ij}\f$
+ *
+ * where \f$\vec{r}_{ij}=\vec{r}_i - \vec{r}_j\f$ is the distance vector
+ * connecting particle centers and \f$r_{ij} = \|\vec{r}_{ij}\|_2\f$ its module.
+ * The normal and tangential components of the relative velocity at the point of
+ * contact is given by
+ *
+ * \f$\vec{v}_{n_{ij}} = \left(\vec{v}_{ij}\cdot\vec{n}_{ij}\right)\vec{n}_{ij}\f$
+ * \f$\vec{v}_{t_{ij}} = \vec{v}_{ij}-\vec{v}_{n_{ij}}-R\left(\vec{\omega}_i + \vec{\omega}_j \right)\times \vec{n}_{ij} \f$
+ *
+ * whith \f$\vec{n}_{ij}=\vec{r}_{ij}/r_{ij}\f$ is the normal unit vector in direction of the distance vector,
+ * \f$\vec{\omega}_i\f$ is the angular velocity of a particle and
+ * \f$\vec{v}_{ij}=\vec{v}_i-\vec{v}_j\f$ the relative velocity between the two
+ * particles. The evolution of the elastic tangential displacement \f$\vec{u}_{t_{ij}}\f$
+ * is integrated when two particles are in contact using.
+ *
+ * \f$\vec{u'}_{t_{ij}} = \vec{u}_{t_{ij}} + \vec{v}_{t_{ij}} \delta t\f$
+ *
+ * Where $\delta t$ is the time step size. The deformation of the contacts points is
+ *  stored for each particle and for each new contact point the elastic tangential displacement
+ *  is initialized with \f$\vec{u}_{t_{ij}} = 0\f$. Thus for each pair of particle interacting
+ *  the normal and tangential forces become:
+ *
+ * \f$\vec{F}_{n_{ij}}=\sqrt{\frac{\delta_{ij}}{2R}}\,\,\left(k_n\delta_{ij}\vec{n}_{ij}-\gamma_n
+ * m_{\text{eff}}\vec{v}_{n_{ij}}\right) \, ,
+ * \f$
+ * \f$\vec{F}_{t_{ij}}=\sqrt{\frac{\delta_{ij}}{2R}}\,\,\left(-k_t\vec{u}_{t_{ij}}-\gamma_t
+ * m_{\text{eff}}\vec{v}_{t_{ij}}\right) \, ,\f$
+ *
+ * where \f$k_{n,t}\f$ are the elastic constants in normal and tangential direction,
+ * respectively, and \f$\gamma _{n,t}\f$ the corresponding viscoelastic constants. The
+ * effective collision mass is given by \f$m_{\text{eff}}=\frac{m}{2}\f$. For each
+ * contact point in order to enforce Coulomb's law \f$\|\vec{F}_{t_{ij}}\|_2 < \|
+ * \mu\vec{F}_{n_{ij}}\|_2\f$ the tangential force is bounded
+ * by the normal component force. In particular the elastic tangential displacement \f$\vec{u}_{t_{ij}}\f$ is
+ * adjusted with \f$\vec{F}_{t_{ij}} \leftarrow
+ * \vec{F}_{t_{ij}}\frac{\|\mu\vec{F}_{n_{ij}}\|_2}{\|\vec{F}_{t_{ij}}\|_2}
+ * \f$
+ *
+ * This adjustment induce a truncation of the elastic displacement. The Coulomb condition
+ * is equivalent to the case where two spheres slip against each other without
+ * inducing additional deformations. Thus the deformation is truncated using:
+ *
+ * \f$\vec{u}_{t_{ij}} =
+ * -\frac{1}{k_t}\left(\vec{F}_{t_{ij}} \sqrt{\frac{2R}{\delta_{ij}}} + \gamma_t
+ * m_{\text{eff}}\vec{v}_{t_{ij}}\right) \, .\f$
+ *
+ * Considering that each particle \f$i\f$ interact with all the particles $j$ is in touch with , the total resultant force on particle $i$ is then computed by summing the contributions of all pair particles $(i,j)$. Considering that the grains are also under the effect of the gravitational field we obtain that the total force is given by
+ *
+ * \f$\vec{F}_i^{\text{tot}}=m \vec{g} + \sum_j
+ * \left(\vec{F}_{n_{ij}}+\vec{F}_{t_{ij}}\right) \, ,
+ * \f$
+ *
+ * where \f$\vec{g}\f$ is the acceleration due to gravity.
+ * Because particles has also rotational degree of freedoms, the total torque on particle $i$ is calculated using
+ *
+ * \f$\vec{T}_i^{\text{tot}} = -R \sum_j \vec{n}_{ij}\times\vec{F}_{t_{ij}} \, .\f$
+ *
+ * \f$\vec{r}_i\f$ and angular velocities \f$\vec{\omega}_i\f$ for each particle \f$i\f$ at
+ * time step \f$n+1\f$,
+ * We integrate in time the equations
+ * using leap frog scheme
+ * with time step given by $\delta t = 10^{-6}\,$s
+ *
+ * \f$\vec{v}_i^{n+1} = \vec{v}_i^n + \frac{\delta t}{m}\vec{F}_i^{\text{tot}} \, ,
+ * \qquad
+ * \vec{r}_i^{n+1} = \vec{r}_i^n + \delta t \vec{v}_i^{n+1}
+ * \f$
+ * \f$\vec{\omega}_i^{n+1} = \vec{\omega}_i^n + \frac{\delta t}{I_i}\vec{T}_i^{\text{tot}} \, ,\f$
+ *
+ * where \f$\vec{r}_i^{n},\vec{v}_i^{n},\vec{\omega}_i^{n}\f$ denotes respectively
+ *  the position, the speed and the rotational speed of the particle $i$ at time step
+ *  \f$n\f$, and \f$\delta t\f$ the time step size.
+ *
+ * We simulate an avalanche down an inclined plane
+ *
+ * ## Constants {#dem_constants}
+ *
+ * This simulation require several constants to work, here we define mass, radius, delta time
+ * elastic constants and stop time, simulation domain and ghost.
+ *
+ * \snippet Vector/8_DEM/main.cpp constants_def
+ *
+ * ## Sand particles {#sand_part}
+ *
+ * The definition of sand particles require:
+ *
+ *  * **velocity**: particle velocity
+ *  * **force**: particle force
+ *  * **omega**: angular velocity
+ *  * **tau**: moment of the force
+ *  * **cpd**: contact point deformation
+ *  * **cpi**: particle index of the contact point
+ *  * **ncp**: number of contact points
+ *  * **tt**: sand particle or recipient particle
+ *  * **gid**: global id unique for each particles
+ *
+ *
+ * ## Initialization
+ *
+ * Initialize the particle in partice we create particles to form a base
+ * two walls
+ *
+ *  \snippet Vector/8_DEM/main.cpp init sand part
+ *
+ * ## DEM Loop
+ *
+ * The DEM Loop is composition of the previous steps
+ *
+ * 1) integration of the force and calculated force momenta
+ *
+ * \snippet Vector/8_DEM/main.cpp integration pos and angvel
+ *
+ * 2) Update neighborhood data-structure if needed and dynamic load balance
+ *
+ * \snippet Vector/8_DEM/main.cpp dynamic load balancing
+ *
+ * 3) Calculate forces and momenta
+ *
+ * \snippet Vector/8_DEM/main.cpp calculate force
+ *
+ * The discrete element loop follow the the step described above.
+ *
+ */
+//#define SE_CLASS3
+//#define CHECKFOR_POSNAN
+//#define STOP_ON_ERROR
+//#define PRINT_STACKTRACE
+
+#ifdef TEST_RUN
+
+#else
+
+#endif
+
+#include "Vector/vector_dist.hpp"
+#include "Draw/DrawParticles.hpp"
+
+int main(int argc, char* argv[])
+{
+    // initialize the library
+	openfpm_init(&argc,&argv);
+
+	//! \cond [constants_def] \endcond
+
+	double m = 1;
+	double R = 0.06;
+	double cutoff = 2 * R;
+	double skin = 0.1 * cutoff;
+	constexpr int max_contacts = 36;
+	constexpr int max_contacts_def = max_contacts * 3;
+	double dt = 0.0001;
+	double Im = 2.0 * R * R * m / 5.0;
+	double k_n = 7849;
+	double k_t = 7849;
+	double gamma_n = 3401;
+	double gamma_t = 3401;
+	double m_eff = 0.5;
+#ifdef TEST_RUN
+	double t_stop = 0.001;
+#else
+	double t_stop = 16.0000;
+#endif
+	double mu = 0.5;
+	double max_vel;
+	Vcluster & v_cl = create_vcluster();
+
+	Box<3,double> domain({0.0,0.0,0.0},{17.04,6.0,6.6});
+	Ghost<3,double> g(cutoff + skin);
+
+	//! \cond [constants_def] \endcond
+
+	constexpr int velocity = 0;
+	constexpr int force = 1;
+	constexpr int omega = 2;
+	constexpr int tau = 3;
+	constexpr int cpd = 4;
+	constexpr int cpi = 5;
+	constexpr int ncp = 6;
+	constexpr int tt = 7;
+	constexpr int gid = 8;
+
+	size_t bc[3] = {NON_PERIODIC,PERIODIC,NON_PERIODIC};
+
+	vector_dist<3,double,aggregate<Point<3,double>,Point<3,double>,Point<3,double>,Point<3,double>,
+	                                  double[max_contacts_def],int[max_contacts],
+									  int,int,int>> parts(0,domain,bc,g);
+
+	//! \cond [init sand part] \endcond
+
+	// virtual grid
+	size_t sz[3] = {143,51,56};
+
+	Box<3,double> sand_box({1.8,0.0,0.18},{8.58,5.9999,2.7});
+
+    // we draw the initialization
+	auto sand_it = DrawParticles::DrawBox(parts,sz,domain,sand_box);
+
+	while (sand_it.isNext())
+	{
+		// ... add a particle ...
+		parts.add();
+
+		// ... and set it position ...
+		parts.getLastPos()[0] = sand_it.get().get(0);
+		parts.getLastPos()[1] = sand_it.get().get(1);
+		parts.getLastPos()[2] = sand_it.get().get(2);
+
+		parts.getLastProp<velocity>()[0] = 0.0;
+		parts.getLastProp<velocity>()[1] = 0.0;
+		parts.getLastProp<velocity>()[2] = 0.0;
+
+		parts.getLastProp<omega>()[0] = 0.0;
+		parts.getLastProp<omega>()[1] = 0.0;
+		parts.getLastProp<omega>()[2] = 0.0;
+
+		parts.getLastProp<tau>()[0] = 0.0;
+		parts.getLastProp<tau>()[1] = 0.0;
+		parts.getLastProp<tau>()[2] = 0.0;
+
+		parts.getLastProp<force>()[0] = 0.0;
+		parts.getLastProp<force>()[1] = 0.0;
+		parts.getLastProp<force>()[2] = 0.0;
+
+		parts.getLastProp<ncp>() = 0;
+
+		parts.getLastProp<tt>() = 0;
+
+		++sand_it;
+	}
+
+	Box<3,double> base_box({0.06,0.0,0.06},{16.98,5.9999,0.18});
+
+    // we draw the initialization
+	auto base_it = DrawParticles::DrawBox(parts,sz,domain,base_box);
+
+	while (base_it.isNext())
+	{
+		// ... add a particle ...
+		parts.add();
+
+		// ... and set it position ...
+		parts.getLastPos()[0] = base_it.get().get(0);
+		parts.getLastPos()[1] = base_it.get().get(1);
+		parts.getLastPos()[2] = base_it.get().get(2);
+
+		parts.getLastProp<tt>() = 1;
+
+		++base_it;
+	}
+
+	Box<3,double> wall_front({16.86,0.0,0.06},{16.98,5.9999,6.54});
+
+    // we draw the initialization
+	auto wall_f_it = DrawParticles::DrawBox(parts,sz,domain,wall_front);
+
+	while (wall_f_it.isNext())
+	{
+		// ... add a particle ...
+		parts.add();
+
+		// ... and set it position ...
+		parts.getLastPos()[0] = wall_f_it.get().get(0);
+		parts.getLastPos()[1] = wall_f_it.get().get(1);
+		parts.getLastPos()[2] = wall_f_it.get().get(2);
+
+		parts.getLastProp<tt>() = 1;
+
+		++wall_f_it;
+	}
+
+	Box<3,double> wall_back({0.06,0.0,0.06},{0.18,5.9999,6.54});
+
+    // we draw the initialization
+	auto wall_b_it = DrawParticles::DrawBox(parts,sz,domain,wall_back);
+
+	while (wall_b_it.isNext())
+	{
+		// ... add a particle ...
+		parts.add();
+
+		// ... and set it position ...
+		parts.getLastPos()[0] = wall_b_it.get().get(0);
+		parts.getLastPos()[1] = wall_b_it.get().get(1);
+		parts.getLastPos()[2] = wall_b_it.get().get(2);
+
+		parts.getLastProp<tt>() = 1;
+
+		++wall_b_it;
+	}
+
+	//! \cond [init sand part] \endcond
+
+	parts.map();
+	parts.addComputationCosts(ModelSquare());
+	parts.getDecomposition().decompose();
+	parts.map();
+
+	// Fill the gid
+
+	auto it_p = parts.getDomainIterator();
+	size_t accu = parts.accum();
+
+	while (it_p.isNext())
+	{
+		auto p = it_p.get();
+
+		parts.getProp<gid>(p) = accu;
+
+		++accu;
+		++it_p;
+	}
+
+	parts.ghost_get<>();
+
+
+	size_t cnt = 0;
+	size_t cnt_reb = 0;
+	auto nlist = parts.getVerlet<VERLETLIST_BAL(3,double)>(cutoff + skin);
+
+	double tot_sp = 0.0;
+
+	double t = 0.0;
+	while (t < t_stop)
+	{
+		auto pit = parts.getDomainIterator();
+
+		max_vel = 0.0;
+
+		//! \cond [integration pos and angvel] \endcond
+
+		// Update
+		while (pit.isNext())
+		{
+			auto p = pit.get();
+
+			if (parts.getProp<tt>(p) == 1)	{++pit;continue;}
+
+			parts.getProp<velocity>(p)[0] = parts.getProp<velocity>(p)[0] + parts.getProp<force>(p)[0]*dt;
+			parts.getProp<velocity>(p)[1] = parts.getProp<velocity>(p)[1] + parts.getProp<force>(p)[1]*dt;
+			parts.getProp<velocity>(p)[2] = parts.getProp<velocity>(p)[2] + parts.getProp<force>(p)[2]*dt;
+
+			double norm2 = parts.getProp<velocity>(p)[0]*parts.getProp<velocity>(p)[0] +
+						   parts.getProp<velocity>(p)[1]*parts.getProp<velocity>(p)[1] +
+						   parts.getProp<velocity>(p)[2]*parts.getProp<velocity>(p)[2];
+			if (max_vel < norm2)
+			{max_vel = norm2;}
+
+			parts.getPos(p)[0] = parts.getPos(p)[0] + parts.getProp<velocity>(p)[0]*dt;
+			parts.getPos(p)[1] = parts.getPos(p)[1] + parts.getProp<velocity>(p)[1]*dt;
+			parts.getPos(p)[2] = parts.getPos(p)[2] + parts.getProp<velocity>(p)[2]*dt;
+
+		    if (parts.getPos(p)[0] < domain.getLow(0) || parts.getPos(p)[0] > domain.getHigh(0) ||
+				parts.getPos(p)[2] < domain.getLow(2) || parts.getPos(p)[2] > domain.getHigh(2) )
+		    {parts.getProp<tt>(p) = 1;}
+
+		    parts.getProp<omega>(p)[0] = parts.getProp<omega>(p)[0] + parts.getProp<tau>(p)[0]/Im*dt;
+		    parts.getProp<omega>(p)[1] = parts.getProp<omega>(p)[1] + parts.getProp<tau>(p)[1]/Im*dt;
+		    parts.getProp<omega>(p)[2] = parts.getProp<omega>(p)[2] + parts.getProp<tau>(p)[2]/Im*dt;
+
+			++pit;
+		}
+		tot_sp += sqrt(max_vel)*dt;
+		v_cl.max(tot_sp);
+		v_cl.execute();
+
+		//! \cond [integration pos and angvel] \endcond
+
+		//! \cond [dynamic load balancing] \endcond
+
+		if (tot_sp >= skin / 2.0)
+		{
+			parts.map();
+
+			// Check if it is time to rebalance
+
+			if (cnt_reb >= 200)
+			{
+				if (v_cl.rank() == 0)
+				{std::cout << "Redecomposing" << std::endl;}
+				cnt_reb = 0;
+				parts.addComputationCosts(ModelSquare());
+				parts.getDecomposition().redecompose(200);
+				parts.map();
+			}
+
+			if (v_cl.rank() == 0)
+			{std::cout << "Reconstruct Verlet" << std::endl;}
+
+			parts.ghost_get<velocity,omega,gid>();
+			parts.updateVerlet(nlist,cutoff+skin);
+
+			tot_sp = 0.0;
+		}
+		else
+		{
+			parts.ghost_get<velocity,omega,gid>();
+		}
+
+		//! \cond [dynamic load balancing] \endcond
+
+		//! \cond [calculate force] \endcond
+
+		auto pit2 = parts.getDomainIterator();
+
+		while (pit2.isNext())
+		{
+			auto p = pit2.get();
+			Point<3,double> xp = parts.getPos(p);
+			Point<3,double> v_p = parts.getProp<velocity>(p);
+			Point<3,double> omega_p = parts.getProp<omega>(p);
+
+			if (parts.getProp<tt>(p) == 1)	{++pit2;continue;}
+
+			Point<3,double> dF_n({0.0,0.0,0.0});
+			Point<3,double> dF_t({0.0,0.0,0.0});
+			Point<3,double> dTau({0.0,0.0,0.0});
+
+			int contact_ok[max_contacts];
+			for (size_t i = 0; i < max_contacts ; i++)	{contact_ok[i] = 0;}
+
+			auto NN = nlist.getNNIterator(p.getKey());
+
+			while (NN.isNext())
+			{
+				auto q = NN.get();
+
+				if (q == p.getKey())	{++NN;continue;}
+
+				Point<3,double> xq = parts.getPos(q);
+				Point<3,double> v_q = parts.getProp<velocity>(q);
+				Point<3,double> omega_q = parts.getProp<omega>(q);
+
+				Point<3,double> r_pq = xp - xq;
+				double r_s_pq2 = norm2(r_pq);
+
+				// Norm is not defined, next particle
+				if (r_s_pq2 == 0)
+				{continue;}
+
+				double delta_ij = 2.0*R - sqrt(r_s_pq2);
+
+				if (delta_ij < 0.0)	{++NN;continue;}
+
+				size_t cnt_end = parts.getProp<ncp>(p);
+				int this_contact = cnt_end;
+
+				for (size_t k = 0 ; k < cnt_end ; k++)
+				{
+					if (parts.getProp<cpi>(p)[k] == parts.getProp<gid>(q))
+					{
+						this_contact = k;
+					}
+				}
+
+				int cidx;
+				if (this_contact == cnt_end)
+				{
+					parts.getProp<ncp>(p) += 1;
+					this_contact = parts.getProp<ncp>(p) - 1;
+
+					cidx = 3 * this_contact;
+
+					if (this_contact >= max_contacts)
+					{std::cout << "Error reached maximum nunber of contacts points" << std::endl;}
+
+					parts.getProp<cpi>(p)[this_contact] = parts.getProp<gid>(q);
+
+					parts.getProp<cpd>(p)[cidx] = 0.0;
+					parts.getProp<cpd>(p)[cidx + 1] = 0.0;
+					parts.getProp<cpd>(p)[cidx + 2] = 0.0;
+
+				}
+				else
+				{
+					cidx = 3 * this_contact;
+				}
+
+				Point<3,double> n_ij = r_pq / sqrt(r_s_pq2);
+				Point<3,double> v_rel = v_p - v_q;
+				Point<3,double> v_nij = (v_rel * n_ij) * n_ij;
+				Point<3,double> v_omega = (omega_p + omega_q)*R;
+				Point<3,double> v_cross;
+
+				v_cross.get(0) = v_omega.get(1) * n_ij.get(2) - v_omega.get(2) * n_ij.get(1);
+				v_cross.get(1) = v_omega.get(2) * n_ij.get(0) - v_omega.get(0) * n_ij.get(2);
+				v_cross.get(2) = v_omega.get(0) * n_ij.get(1) - v_omega.get(1) * n_ij.get(0);
+
+				Point<3,double> v_tij = v_rel - v_nij - v_cross;
+				Point<3,double> v_dtij = dt * v_tij;
+
+				parts.getProp<cpd>(p)[cidx] += v_dtij.get(0);
+				parts.getProp<cpd>(p)[cidx + 1] += v_dtij.get(1);
+				parts.getProp<cpd>(p)[cidx + 2] += v_dtij.get(2);
+
+				Point<3,double> u_ij;
+
+				u_ij.get(0) = parts.getProp<cpd>(p)[cidx];
+				u_ij.get(1) = parts.getProp<cpd>(p)[cidx + 1];
+				u_ij.get(2) = parts.getProp<cpd>(p)[cidx + 2];
+
+				Point<3,double> F_nij = sqrt(delta_ij/2/R) * (k_n*delta_ij*n_ij - gamma_t*m_eff*v_nij);
+				dF_n = dF_n + F_nij;
+
+				Point<3,double> F_tij = sqrt(delta_ij/2/R) * (-k_t*u_ij - gamma_t*m_eff*v_tij);
+				double F_tij_sq = norm2(F_tij);
+				double F_nij_sq = mu * mu * norm2(F_nij);
+				if (F_tij_sq > F_nij_sq)
+				{
+					F_tij = F_tij * (F_nij_sq / F_tij_sq);
+
+					parts.getProp<cpd>(p)[cidx] = -1.0 / k_t * (F_tij.get(0) * sqrt(2*R/delta_ij) + gamma_t*m_eff*v_tij.get(0));
+					parts.getProp<cpd>(p)[cidx+1] = -1.0 / k_t * (F_tij.get(1) * sqrt(2*R/delta_ij) + gamma_t*m_eff*v_tij.get(1));
+					parts.getProp<cpd>(p)[cidx+2] = -1.0 / k_t * (F_tij.get(2) * sqrt(2*R/delta_ij) + gamma_t*m_eff*v_tij.get(2));
+				}
+
+
+		        dF_t = dF_t + F_tij;
+		        dTau.get(0) = dTau.get(0) - R * (n_ij.get(1) * F_tij.get(2) - n_ij.get(2) * F_tij.get(1));
+		        dTau.get(1) = dTau.get(1) - R * (n_ij.get(2) * F_tij.get(0) - n_ij.get(0) * F_tij.get(2));
+		        dTau.get(2) = dTau.get(2) - R * (n_ij.get(0) * F_tij.get(1) - n_ij.get(1) * F_tij.get(0));
+
+		        contact_ok[this_contact] = 1;
+
+				++NN;
+			}
+
+		    int cnt_end = parts.getProp<ncp>(p);
+		    int i = 0;
+		    for (int iread = 0; iread < cnt_end ; iread++)
+		    {
+		    	if (contact_ok[iread] == 1)
+		    	{
+		    		i = i + 1;
+		    		int j = 3*(i - 1);
+					int k = 3*iread;
+
+					parts.getProp<cpd>(p)[j] = parts.getProp<cpd>(p)[k];
+					parts.getProp<cpd>(p)[j+1] = parts.getProp<cpd>(p)[k+1];
+					parts.getProp<cpd>(p)[j+2] = parts.getProp<cpd>(p)[k+2];
+		    	}
+		    }
+
+		    parts.getProp<ncp>(p) = i;
+
+		    if (parts.getProp<tt>(p) == 0)
+		    {
+		        parts.getProp<force>(p).get(0) = m * 4.905 + dF_n.get(0) + dF_t.get(0);
+		        parts.getProp<force>(p).get(1) = 0.0 + dF_n.get(1) + dF_t.get(1);
+		        parts.getProp<force>(p).get(2) = m * -8.49570921 + dF_n.get(2) + dF_t.get(2);
+
+		        parts.getProp<tau>(p) = dTau;
+		    }
+
+		    if (parts.getProp<tt>(p) == 1)
+		    {
+		    	parts.getProp<force>(p) = 0;
+		    	parts.getProp<tau>(p) = 0;
+		    }
+
+			++pit2;
+		}
+
+		//! \cond [calculate force] \endcond
+
+		if (v_cl.rank() == 0)
+		{std::cout << "Time step" << std::endl;}
+
+		if (cnt % 300 == 0)
+		{
+			std::cout << "Write " << cnt << std::endl;
+			parts.write_frame("output",cnt);
+		}
+
+		cnt_reb++;
+		cnt++;
+		t += dt;
+	}
+
+	openfpm_finalize();
+}