Commit cecff1f7 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts

fix number of particle seeding by element roundings

parent 847be77e
Pipeline #1700 failed with stage
in 30 minutes and 20 seconds
......@@ -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;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment