fix number of particle seeding by element roundings

modules/slimParticles/dgParticleTracker2D.cpp

@@ -1113,7 +1113,10 @@ void particleArray::addParticlesWithDofContainer(dgGroupCollection* groups, dgDo
       fullMatrix<double> ecoord;
       groupOfElements->elementCoordinatesProxy(elementId, ecoord);
       double elSurface = groupOfElements->getFunctionSpace().volume(ecoord);
-      int nbParticlesInElem = static_cast<int>((concentration/1000000.0)*elSurface + 0.5);  //static_cast rounds down
+      double nbParticlesInElemDouble = (concentration/1000000.0)*elSurface;
+      int nbParticlesInElem = static_cast<int>(nbParticlesInElemDouble);  //static_cast rounds down
+      if( rand()*1./RAND_MAX < (nbParticlesInElemDouble - nbParticlesInElem))
+        nbParticlesInElem += 1;
       if (nbParticlesInElem != 0) {
         //If reef does not have enough particles for minNbParticlesPerHabitat, add however many needed on this element
         if (particlesMissing[habitatId] > 0) {
@@ -1129,7 +1132,7 @@ void particleArray::addParticlesWithDofContainer(dgGroupCollection* groups, dgDo
       
       //Now seed particles
       for (int j = 0; j < nbParticlesInElem; j++) {
-//        if (j==0) printf("Number of particles in reef element %d (habitatId %d) is %d\n",nbOfHabitatElements, habitatId, nbParticlesInElem );
+        if (j==0) printf("Number of particles in reef element %d (habitatId %d) is %d\n",nbOfHabitatElements, habitatId, nbParticlesInElem );
         double u=rand()*1./RAND_MAX;
         double v=rand()*1./RAND_MAX;
         if(u+v>1){
@@ -1285,6 +1288,7 @@ void particleArray::addParticlesWithReefMap(dgGroupCollection* groups, dgDofCont
       if (reefId > nbOfShallowReefs && seedOverWhichReefs == 0) continue;
       if (reefId <= nbOfShallowReefs && seedOverWhichReefs == 1) continue;
 
+
       //Seed element in parametric space:
       //First calc. no. of particles to seed in element, and count reef as having been seeded
       fullMatrix<double> ecoord;