code seems to compile

bfps/DNS.py

@@ -865,6 +865,7 @@ class DNS(_code):
                     rseed = opt.particle_rand_seed)
             if not os.path.exists(self.get_particle_file_name()):
                 with h5py.File(self.get_particle_file_name(), 'w') as particle_file:
+                    particle_file.create_group('tracers0/position')
                     particle_file.create_group('tracers0/velocity')
                     particle_file.create_group('tracers0/acceleration')
         return None

bfps/cpp/full_code/NSVEparticles.cpp

@@ -77,6 +77,14 @@ int NSVEparticles<rnumber>::do_stats()
     if (!(this->iteration % this->niter_part == 0))
         return EXIT_SUCCESS;
 
+    /// sample position
+    sample_particles_system_position(
+                                 this->ps,
+                                 (this->simname + "_particles.h5"), // filename
+                                 "tracers0",                        // hdf5 parent group
+                                 "position"                         // dataset basename TODO
+                                 );
+
     /// sample velocity
     sample_from_particles_system(*this->tmp_vec_field,              // field to save
                                  this->ps,

bfps/cpp/particles/particles_sampling.hpp

@@ -48,5 +48,37 @@ void sample_from_particles_system(const field<rnumber, be, fc>& in_field, // a p
                      ps->get_step_idx());
 }
 
-#endif
+template <class partsize_t, class particles_rnumber>
+void sample_particles_system_position(
+        std::unique_ptr<abstract_particles_system<partsize_t, particles_rnumber>>& ps, // a pointer to an particles_system<double>
+                                  const std::string& filename,
+                                  const std::string& parent_groupname,
+                                  const std::string& fname){
+    const std::string datasetname = fname + std::string("/") + std::to_string(ps->get_step_idx());
+
+    // Stop here if already exists
+    if(particles_output_sampling_hdf5<partsize_t, particles_rnumber, 3, 3>::DatasetExistsCol(MPI_COMM_WORLD,
+                                                                                             filename,
+                                                                                             parent_groupname,
+                                                                                             datasetname)){
+        return;
+    }
+
+    const partsize_t nb_particles = ps->getLocalNbParticles();
+    std::unique_ptr<particles_rnumber[]> sample_rhs(new particles_rnumber[3*nb_particles]);
+    std::copy(ps->getParticlesPositions(), ps->getParticlesPositions() + 3*nb_particles, sample_rhs.get());
+
+    particles_output_sampling_hdf5<partsize_t, particles_rnumber, 3, 3> outputclass(MPI_COMM_WORLD,
+                                                                                    ps->getGlobalNbParticles(),
+                                                                                    filename,
+                                                                                    parent_groupname,
+                                                                                    datasetname);
+    outputclass.save(ps->getParticlesPositions(),
+                     &sample_rhs,
+                     ps->getParticlesIndexes(),
+                     ps->getLocalNbParticles(),
+                     ps->get_step_idx());
+}
+
+#endif//PARTICLES_SAMPLING_HPP