yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2894
committer: Anton Gladky <gladky.anton@xxxxxxxxx>
branch nick: yade
timestamp: Wed 2011-07-20 17:30:57 +0200
message:
  1. Add PSD-support to SpheresFactory (not finished, but already working)
  2. Add exception in utils.psd(), if all particles have the same diameter
modified:
  pkg/dem/SpheresFactory.cpp
  pkg/dem/SpheresFactory.hpp
  py/utils.py


--
lp:yade
https://code.launchpad.net/~yade-dev/yade/trunk

Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-dev/yade/trunk/+edit-subscription

=== modified file 'pkg/dem/SpheresFactory.cpp'
--- pkg/dem/SpheresFactory.cpp	2011-07-19 07:06:08 +0000
+++ pkg/dem/SpheresFactory.cpp	2011-07-20 15:30:57 +0000
@@ -28,14 +28,49 @@
 		if(!collider) throw runtime_error("SpheresFactory: No Collider instance found in engines (needed for collision detection).");
 	}
 	goalMass+=massFlowRate*scene->dt; // totalMass that we want to attain in the current step
-
+	if ((PSDcum.size()>0) and (!PSDuse)) {			//Defined, that we will use PSD
+		if (PSDcum.size() != PSDsizes.size()) {
+			LOG_ERROR("PSDcum and PSDsizes should have an equal number of elements.");
+			throw std::logic_error("PSDcum and PSDsizes should have an equal number of elements.");
+		}
+		PSDuse = true;
+		
+		//Prepare main vectors
+		for (unsigned int i=0; i<PSDcum.size(); i++) {
+			if (i==0) {
+				PSDNeedProc.push_back(PSDcum[i]);
+			} else {
+				PSDNeedProc.push_back(PSDcum[i]-PSDcum[i-1]);
+			}
+			PSDCurMean.push_back(0);
+			PSDCurProc.push_back(0);
+		}
+	}
+		
 	normal.normalize();
 
 	LOG_TRACE("totalMass="<<totalMass<<", goalMass="<<goalMass);
 
 	while(totalMass<goalMass && (maxParticles<0 || numParticles<maxParticles)){
-		// pick random radius (later: based on the particle size distribution)
-		Real r=rMin+randomUnit()*(rMax-rMin);
+		Real r=0.0;
+		
+		Real maxdiff=0.0;		//This and next variable are for PSD-distribution
+		int maxdiffID=0;
+		
+		if (PSDuse) {
+			//find in what "bin" we have maximal difference between required number of material and current:
+			for (unsigned int k=0; k<PSDcum.size(); k++) {
+				if ((maxdiff < (PSDNeedProc[k]-PSDCurProc[k]))) {
+					maxdiff = PSDNeedProc[k]-PSDCurProc[k];
+					maxdiffID = k;
+				}
+			}
+			r=PSDsizes[maxdiffID]/2.0;
+		} else {
+			// pick random radius
+			r=rMin+randomUnit()*(rMax-rMin);
+		}
+		
 		LOG_TRACE("Radius "<<r);
 		Vector3r c=Vector3r::Zero();
 		// until there is no overlap, pick a random position for the new particle
@@ -83,10 +118,22 @@
 		// insert particle in the simulation
 		scene->bodies->insert(b);
 		ids.push_back(b->getId());
-		// increment total mass we've spit out
+		// increment total mass, volume and numparticles we've spit out
 		totalMass+=state->mass;
+		totalVolume+= vol;
 		numParticles++;
 		
+		if (PSDuse) {		//Add newly created "material" into the bin
+			Real summMaterial = 0.0;
+			if (PSDcalculate=="mass") { PSDCurMean[maxdiffID]=PSDCurMean[maxdiffID]+state->mass; summMaterial = totalMass;}
+			else if (PSDcalculate=="volume") { PSDCurMean[maxdiffID]=PSDCurMean[maxdiffID]+vol; summMaterial = totalVolume;}
+			else if (PSDcalculate=="number") { PSDCurMean[maxdiffID]=PSDCurMean[maxdiffID]+1; summMaterial = numParticles;}
+			
+			for (unsigned int k=0; k<PSDcum.size(); k++) {			//Update  relationships in bins
+				PSDCurProc[k] = PSDCurMean[k]/summMaterial;
+			}
+		}	
+		
 	} 
 	//std::cout<<"mass flow rate: "<<totalMass<<endl;totalMass =0.0;
 };

=== modified file 'pkg/dem/SpheresFactory.hpp'
--- pkg/dem/SpheresFactory.hpp	2011-07-20 07:49:38 +0000
+++ pkg/dem/SpheresFactory.hpp	2011-07-20 15:30:57 +0000
@@ -12,10 +12,11 @@
 		vector<Real> PSDCurMean;  //Current value of material in each bin
 		vector<Real> PSDCurProc;  //Current value of material in each bin, in procents
 		vector<Real> PSDNeedProc; //Need value of procent in each bin
+		bool PSDuse;        //PSD or not
 	public:
 		virtual void action();
 	DECLARE_LOGGER;
-	YADE_CLASS_BASE_DOC_ATTRS(SpheresFactory,GlobalEngine,"Engine for spitting spheres based on mass flow rate, particle size distribution etc. Initial velocity of particles is given by *vMin*, *vMax*, the *massFlowRate* determines how many particles to generate at each step. When *goalMass* is attained or positive *maxParticles* is reached, the engine does not produce particles anymore. Geometry of the region should be defined in a derived engine by overridden SpheresFactory::pickRandomPosition(). \n\nA sample script for this engine is in :ysrc:`scripts/spheresFactory.py`.",
+	YADE_CLASS_BASE_DOC_ATTRS_CTOR(SpheresFactory,GlobalEngine,"Engine for spitting spheres based on mass flow rate, particle size distribution etc. Initial velocity of particles is given by *vMin*, *vMax*, the *massFlowRate* determines how many particles to generate at each step. When *goalMass* is attained or positive *maxParticles* is reached, the engine does not produce particles anymore. Geometry of the region should be defined in a derived engine by overridden SpheresFactory::pickRandomPosition(). \n\nA sample script for this engine is in :ysrc:`scripts/spheresFactory.py`.",
 		((Real,massFlowRate,NaN,,"Mass flow rate [kg/s]"))
 		((Real,rMin,NaN,,"Minimum radius of generated spheres (uniform distribution)"))
 		((Real,rMax,NaN,,"Maximum radius of generated spheres (uniform distribution)"))
@@ -27,15 +28,17 @@
 		((int,mask,-1,,"groupMask to apply for newly created spheres "))
 		((vector<int>,ids,,,"ids of created bodies"))
 		((Real,totalMass,0,,"Mass of spheres that was produced so far. |yupdate|"))
+		((Real,totalVolume,0,,"Volume of spheres that was produced so far. |yupdate|"))
 		((Real,goalMass,0,,"Total mass that should be attained at the end of the current step. |yupdate|"))
 		((int,maxParticles,100,,"The number of particles at which to stop generating new ones (regardless of massFlowRate"))
 		((int,numParticles,0,,"Cummulative number of particles produces so far |yupdate|"))
 		((int,maxAttempt,5000 ,,"Maximum number of attempts to position a new sphere randomly."))
 		((bool,silent,false ,,"If true no complain about excessing maxAttempt but disable the factory (by set massFlowRate=0)."))
 		((std::string,blockedDOFs,"" ,,"Blocked degress of freedom"))
-		((vector<Real>,PSDsizes,,,"PSD-dispersion, sizes of cells [m]"))
+		((vector<Real>,PSDsizes,,,"PSD-dispersion, sizes of cells, Diameter [m]"))
 		((vector<Real>,PSDcum,,,"PSD-dispersion, cumulative procent meanings [-]"))
-		((std::string,PSDcalculate,"mass",,"How the PSD will be calculated: mass, number or volume of particles"))
+		((std::string,PSDcalculate,"mass",,"How the PSD will be calculated: mass, number or volume of particles")),
+		;
 	);
 };
 REGISTER_SERIALIZABLE(SpheresFactory);

=== modified file 'py/utils.py'
--- py/utils.py	2011-07-19 14:52:47 +0000
+++ py/utils.py	2011-07-20 15:30:57 +0000
@@ -785,6 +785,8 @@
 			if ((2*b.shape.radius)	> maxD) : maxD = 2*b.shape.radius
 			if (((2*b.shape.radius)	< minD) or (minD==0.0)): minD = 2*b.shape.radius
 
+	if (minD==maxD): return false       #All particles are having the same size
+  
 	binsSizes = numpy.linspace(minD, maxD, bins+1)
 	
 	deltaBinD = (maxD-minD)/bins