yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2687
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
branch nick: yade
timestamp: Thu 2011-01-27 18:05:23 +0100
message:
  makecloud update
  - Fix psd transformation from mass distrib to number distrib.
  - Make sure the function will always return the number of particles required by the user as soon as it is psecified.
  - If num is specified, generate sizes the deterministic way by descending order.
  - make porosity parameter optional and give it a default value.
  - update doc and script. 
modified:
  pkg/dem/SpherePack.cpp
  pkg/dem/SpherePack.hpp
  py/pack/_packSpheres.cpp
  scripts/test/psd.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/SpherePack.cpp'
--- pkg/dem/SpherePack.cpp	2011-01-21 11:28:22 +0000
+++ pkg/dem/SpherePack.cpp	2011-01-27 17:05:23 +0000
@@ -85,18 +85,19 @@
 long SpherePack::makeCloud(Vector3r mn, Vector3r mx, Real rMean, Real rRelFuzz, int num, bool periodic, Real porosity, const vector<Real>& psdSizes, const vector<Real>& psdCumm, bool distributeMass){
 	static boost::minstd_rand randGen(TimingInfo::getNow(/* get the number even if timing is disabled globally */ true));
 	static boost::variate_generator<boost::minstd_rand&, boost::uniform_real<Real> > rnd(randGen, boost::uniform_real<Real>(0,1));
-	vector<Real> psdRadii; // hold psdSizes as volume if distributeMass is in use, as plain radii (rather than diameters) if not
+	vector<Real> psdRadii; // holds plain radii (rather than diameters), scaled down in some situations to get the target number
 	vector<Real> psdCumm2; // psdCumm but dimensionally transformed to match mass distribution
+	Vector3r size=mx-mn; Real volume=size.x()*size.y()*size.z();
 	int mode=-1; bool err=false;
 	// determine the way we generate radii
+	if(porosity<=0) {
+		LOG_WARN("porosity must be >0, changing it for you. It will be ineffective if num<=0.");
+		porosity=0.5;}
 	if(rMean>0) mode=RDIST_RMEAN;
-	if(porosity>0){
-		err=(mode>=0); mode=RDIST_POROSITY;
-		if(num<0) throw invalid_argument("SpherePack.makeCloud: when specifying porosity, num must be given as well.");
-		Vector3r dimensions=mx-mn; Real volume=dimensions.x()*dimensions.y()*dimensions.z();
+	if(num>0 && psdSizes.size()==0){
+		err=(mode>=0); mode=RDIST_NUM;
 		// the term (1+rRelFuzz²) comes from the mean volume for uniform distribution : Vmean = 4/3*pi*Rmean*(1+rRelFuzz²) 
-		rMean=pow(volume*(1-porosity)/(Mathr::PI*(4/3.)*(1+rRelFuzz*rRelFuzz)*num),1/3.);
-	}
+		rMean=pow(volume*(1-porosity)/(Mathr::PI*(4/3.)*(1+rRelFuzz*rRelFuzz)*num),1/3.);}
 	if(psdSizes.size()>0){
 		err=(mode>=0); mode=RDIST_PSD;
 		if(psdSizes.size()!=psdCumm.size()) throw invalid_argument(("SpherePack.makeCloud: psdSizes and psdCumm must have same dimensions ("+lexical_cast<string>(psdSizes.size())+"!="+lexical_cast<string>(psdCumm.size())).c_str());
@@ -105,36 +106,51 @@
 		psdRadii.reserve(psdSizes.size());
 		for(size_t i=0; i<psdSizes.size(); i++) {
 			psdRadii.push_back(/* radius, not diameter */ .5*psdSizes[i]);
-			if(distributeMass) psdCumm2.push_back(1-pow(1-psdCumm[i],psdScaleExponent)); // 5/2. is an approximate (fractal) dimension, empirically determined
+			if(distributeMass) {
+				//psdCumm2 is first obtained by integrating the number of particles over the volumic PSD (dN/dSize = totV*(dPassing/dSize)*1/(4/3πr³)). I suspect similar reasoning behind pow3Interp, since the expressions are a bit similar. The total cumulated number will be the number of spheres in volume*(1-porosity), it is used to decide if the PSD will be scaled down. psdCumm2 is normalized below in order to fit in [0,1]. (B.C.)
+				if (i==0) psdCumm2.push_back(0);
+				else psdCumm2.push_back(psdCumm2[i-1] + 3.0*volume*(1-porosity)/Mathr::PI*(psdCumm[i]-psdCumm[i-1])/(psdSizes[i]-psdSizes[i-1])*(pow(psdSizes[i-1],-2)-pow(psdSizes[i],-2)));
+			}
 			LOG_DEBUG(i<<". "<<psdRadii[i]<<", cdf="<<psdCumm[i]<<", cdf2="<<(distributeMass?lexical_cast<string>(psdCumm2[i]):string("--")));
 			// check monotonicity
 			if(i>0 && (psdSizes[i-1]>psdSizes[i] || psdCumm[i-1]>psdCumm[i])) throw invalid_argument("SpherePack:makeCloud: psdSizes and psdCumm must be both non-decreasing.");
 		}
+		if(distributeMass) {
+			if (num>0) {
+				//if target number will not fit in (1-poro)*volume, scale down particles size
+				if (psdCumm2[psdSizes.size()-1]<num){
+					appliedPsdScaling=pow(psdCumm2[psdSizes.size()-1]/num,1./3.);
+					for(size_t i=0; i<psdSizes.size(); i++) psdRadii[i]*=appliedPsdScaling;}
+			}
+			//Normalize psdCumm2 so it's between 0 and 1
+			for(size_t i=1; i<psdSizes.size(); i++) psdCumm2[i]/=psdCumm2[psdSizes.size()-1];
+		}
 	}
-	if(err || mode<0) throw invalid_argument("SpherePack.makeCloud: exactly one of rMean, porosity, psdSizes & psdCumm arguments must be specified.");
+	if(err || mode<0) throw invalid_argument("SpherePack.makeCloud: at least one of rMean, porosity, psdSizes & psdCumm arguments must be specified. rMean can't be combined with psdSizes.");
 	// adjust uniform distribution parameters with distributeMass; rMean has the meaning (dimensionally) of _volume_
 	const int maxTry=1000;
-	Vector3r size=mx-mn;
 	if(periodic)(cellSize=size);
 	Real r=0;
 	for(int i=0; (i<num) || (num<0); i++) {
-		// determine radius of the next sphere we will attempt to place in space
+		Real norm, rand;
+		//Determine radius of the next sphere we will attempt to place in space. If (num>0), generate radii the deterministic way, in decreasing order, else radii are stochastic since we don't know what the final number will be
+		if (num>0) rand = ((Real)num-(Real)i+0.5)/((Real)num+1.);
+		else rand = rnd();
 		int t;
 		switch(mode){
 			case RDIST_RMEAN:
-			case RDIST_POROSITY: 
-				if(distributeMass) r=pow3Interp(rnd(),rMean*(1-rRelFuzz),rMean*(1+rRelFuzz));
-				else r=rMean*(2*(rnd()-.5)*rRelFuzz+1); // uniform distribution in rMean*(1±rRelFuzz)
+				//FIXME : r is never defined, it will be zero at first iteration, but it will have values in the next ones. Some magic?
+			case RDIST_NUM:
+				if(distributeMass) r=pow3Interp(rand,rMean*(1-rRelFuzz),rMean*(1+rRelFuzz));
+				else r=rMean*(2*(rand-.5)*rRelFuzz+1); // uniform distribution in rMean*(1±rRelFuzz)
 				break;
 			case RDIST_PSD:
 				if(distributeMass){
-					Real norm;
-					int piece=psdGetPiece(rnd(),psdCumm2,norm);
+					int piece=psdGetPiece(rand,psdCumm2,norm);
 					r=pow3Interp(norm,psdRadii[piece],psdRadii[piece+1]);
 				} else {
-					Real norm; int piece=psdGetPiece(rnd(),psdCumm,norm);
-					r=psdRadii[piece]+norm*(psdRadii[piece+1]-psdRadii[piece]);
-				}
+					int piece=psdGetPiece(rand,psdCumm,norm);
+					r=psdRadii[piece]+norm*(psdRadii[piece+1]-psdRadii[piece]);}
 		}
 		// try to put the sphere into a free spot
 		for(t=0; t<maxTry; ++t){
@@ -154,10 +170,16 @@
 			if(!overlap) { pack.push_back(Sph(c,r)); break; }
 		}
 		if (t==maxTry) {
-			if(num>0) LOG_WARN("Exceeded "<<maxTry<<" tries to insert non-overlapping sphere to packing. Only "<<i<<" spheres was added, although you requested "<<num);
-			return i;
-		}
+			if(num>0) {
+				if (mode==RDIST_PSD || mode==RDIST_NUM) {
+					Real nextPoro = porosity+(1-porosity)/10.;
+					if (mode==RDIST_PSD) LOG_WARN("Exceeded "<<maxTry<<" tries to insert non-overlapping sphere to packing. Only "<<i<<" spheres was added, although you requested "<<num<<". Trying again with porosity "<<nextPoro<<". The size distribution is being scaled down");
+					pack.clear();
+					return makeCloud(mn, mx, -1., rRelFuzz, num, periodic, nextPoro, psdSizes, psdCumm, distributeMass);}
+			}
+			return i;}
 	}
+	if (appliedPsdScaling<1) LOG_WARN("The size distribution has been scaled down by a factor pack.appliedPsdScaling="<<appliedPsdScaling);
 	return pack.size();
 }
 

=== modified file 'pkg/dem/SpherePack.hpp'
--- pkg/dem/SpherePack.hpp	2011-01-21 11:28:22 +0000
+++ pkg/dem/SpherePack.hpp	2011-01-27 17:05:23 +0000
@@ -35,7 +35,7 @@
 	struct ClumpInfo{ int clumpId; Vector3r center; Real rad; int minId, maxId; };
 
 public:
-	enum {RDIST_RMEAN, RDIST_POROSITY, RDIST_PSD};
+	enum {RDIST_RMEAN, RDIST_NUM, RDIST_PSD};
 	struct Sph{
 		Vector3r c; Real r; int clumpId;
 		Sph(const Vector3r& _c, Real _r, int _clumpId=-1): c(_c), r(_r), clumpId(_clumpId) {};
@@ -48,7 +48,8 @@
 	std::vector<Sph> pack;
 	Vector3r cellSize;
 	Real psdScaleExponent;
-	SpherePack(): cellSize(Vector3r::Zero()), psdScaleExponent(2.5){};
+	Real appliedPsdScaling;//a scaling factor that can be applied on size distribution
+	SpherePack(): cellSize(Vector3r::Zero()), psdScaleExponent(2.5), appliedPsdScaling(1.){};
 	SpherePack(const python::list& l):cellSize(Vector3r::Zero()){ fromList(l); }
 	// add single sphere
 	void add(const Vector3r& c, Real r){ pack.push_back(Sph(c,r)); }

=== modified file 'py/pack/_packSpheres.cpp'
--- py/pack/_packSpheres.cpp	2011-01-21 11:28:22 +0000
+++ py/pack/_packSpheres.cpp	2011-01-27 17:05:23 +0000
@@ -14,9 +14,10 @@
 		.def("load",&SpherePack::fromFile,(python::arg("fileName")),"Load packing from external text file (current data will be discarded).")
 		.def("save",&SpherePack::toFile,(python::arg("fileName")),"Save packing to external text file (will be overwritten).")
 		.def("fromSimulation",&SpherePack::fromSimulation,"Make packing corresponding to the current simulation. Discards current data.")
-		.def("makeCloud",&SpherePack::makeCloud,(python::arg("minCorner"),python::arg("maxCorner"),python::arg("rMean")=-1,python::arg("rRelFuzz")=0,python::arg("num")=-1,python::arg("periodic")=false,python::arg("porosity")=-1,python::arg("psdSizes")=vector<Real>(),python::arg("psdCumm")=vector<Real>(),python::arg("distributeMass")=false),"Create random loose packing enclosed in box."
-		"\nSphere radius distribution can be specified using one of the following ways (they are mutually exclusive):\n\n#. *rMean* and *rRelFuzz* gives uniform radius distribution between *rMean*(1 ± *rRelFuzz*).\n#. *porosity*, *num* and *rRelFuzz* which estimates mean radius so that *porosity* is attained at the end\n#. *psdSizes* and *psdCumm*, two arrays specifying points of the `particle size distribution <http://en.wikipedia.org/wiki/Particle_size_distribution>`__ function.\n\nBy default (with ``distributeMass==False``), the distribution is applied to particle radii. The usual sense of \"particle size distribution\" is the distribution of *mass fraction* (rather than particle count); this can be achieved with ``distributeMass=True``."
-		"\n\n:param Vector3 minCorner: lower corner of the box\n:param Vector3 maxCorner: upper corner of the box\n:param float rMean: mean radius or spheres\n:param float rRelFuzz: dispersion of radius relative to rMean\n:param int num: number of spheres to be generated (if negavite, generate as many as possible, ending after a fixed number of tries to place the sphere in space)\n:param bool periodic: whether the packing to be generated should be periodic\n:param float porosity: if specified, estimate mean radius $r_m$ (*rMean* must not be given) using *rRelFuzz* ($z$) and *num* ($N$) so that the porosity given ($\\rho$) is approximately achieved at the end of generation, $r_m=\\sqrt[3]{\\frac{V(1-\\rho)}{\\frac{4}{3}\\pi(1+z^2)N}}$. The value of $\\rho$=0.65 is recommended.\n:param psdSizes: sieve sizes (particle diameters) when particle size distribution (PSD) is specified\n:param psdCumm: cummulative fractions of particle sizes given by *psdSizes*; must be the same length as *psdSizes* and should be non-decreasing\n:param bool distributeMass: if ``True``, given distribution will be used to distribute sphere's mass rather than radius of them.\n:returns: number of created spheres, which can be lower than *num* is the packing is too tight.\n")
+		.def("makeCloud",&SpherePack::makeCloud,(python::arg("minCorner"),python::arg("maxCorner"),python::arg("rMean")=-1,python::arg("rRelFuzz")=0,python::arg("num")=-1,python::arg("periodic")=false,python::arg("porosity")=0.5,python::arg("psdSizes")=vector<Real>(),python::arg("psdCumm")=vector<Real>(),python::arg("distributeMass")=false),"Create random loose packing enclosed in box."
+		"\nSphere radius distribution can be specified using one of the following ways:\n\n#. *rMean* and *rRelFuzz* gives uniform radius distribution between *rMean*(1 ± *rRelFuzz*).\n#. *num*,*rRelFuzz*, and (optional) *porosity*, which estimates mean radius so that *porosity* is attained at the end. *porosity* is only an initial guess for the generation algorithm, which will retry with higher porosity until the prescibed *num* is obtained.\n#. *psdSizes* and *psdCumm*, two arrays specifying points of the `particle size distribution <http://en.wikipedia.org/wiki/Particle_size_distribution>`__ function. As many spheres as possible are generated.\n#. *psdSizes*, *psdCumm*, *num*, and (optional) *porosity*, like above but if *num* is not obtained, *psdSizes* will be scaled down uniformly, until *num* is obtained (see :yref:`appliedPsdScaling<yade._packSpheres.SpherePack.appliedPsdScaling>`).\n\nBy default (with ``distributeMass==False``), the distribution is applied to particle radii. The usual sense of \"particle size distribution\" is the distribution of *mass fraction* (rather than particle count); this can be achieved with ``distributeMass=True``."
+		"\n\nIf *num* is defined, then sizes generation is deterministic, giving the best fit of target distribution. It enables spheres placement in descending size order, thus giving lower porosity than the random generation."
+		"\n\n:param Vector3 minCorner: lower corner of the box\n:param Vector3 maxCorner: upper corner of the box\n:param float rMean: mean radius or spheres\n:param float rRelFuzz: dispersion of radius relative to rMean\n:param int num: number of spheres to be generated. If negavite (default), generate as many as possible with stochastic sizes, ending after a fixed number of tries to place the sphere in space, else generate exactly *num* spheres with deterministic size distribution.\n:param bool periodic: whether the packing to be generated should be periodic\n:param float porosity: initial guess for the iterative generation procedure (if *num*>1). The algorithm will be retrying until the number of generated spheres is *num*. The first iteration tries with the provided porosity, but next iterations increase it if necessary (hence an initialy high porosity can speed-up the algorithm). If *psdSizes* is not defined, *rRelFuzz* ($z$) and *num* ($N$) are used so that the porosity given ($\\rho$) is approximately achieved at the end of generation, $r_m=\\sqrt[3]{\\frac{V(1-\\rho)}{\\frac{4}{3}\\pi(1+z^2)N}}$. The default is $\\rho$=0.5. The optimal value depends on *rRelFuzz* or  *psdSizes*.\n:param psdSizes: sieve sizes (particle diameters) when particle size distribution (PSD) is specified\n:param psdCumm: cummulative fractions of particle sizes given by *psdSizes*; must be the same length as *psdSizes* and should be non-decreasing\n:param bool distributeMass: if ``True``, given distribution will be used to distribute sphere's mass rather than radius of them.\n:returns: number of created spheres, which can be lower than *num* is the packing is too tight.\n")
 		.def("psd",&SpherePack::psd,(python::arg("bins")=50,python::arg("mass")=true),"Return `particle size distribution <http://en.wikipedia.org/wiki/Particle_size_distribution>`__ of the packing.\n:param int bins: number of bins between minimum and maximum diameter\n:param mass: Compute relative mass rather than relative particle count for each bin. Corresponds to :yref:`distributeMass parameter for makeCloud<yade.pack.SpherePack.makeCloud>`.\n:returns: tuple of ``(cumm,edges)``, where ``cumm`` are cummulative fractions for respective diameters  and ``edges`` are those diameter values. Dimension of both arrays is equal to ``bins+1``.")
 		// new psd
 		.def("particleSD",&SpherePack::particleSD,(python::arg("minCorner"),python::arg("maxCorner"),python::arg("rMean"),python::arg("periodic")=false,python::arg("name"),python::arg("numSph"),python::arg("radii")=vector<Real>(),python::arg("passing")=vector<Real>(),python::arg("passingIsNotPercentageButCount")=false),"Create random packing enclosed in box given by minCorner and maxCorner, containing numSph spheres. Returns number of created spheres, which can be < num if the packing is too tight. The computation is done according to the given psd.")
@@ -38,7 +39,8 @@
 		.def("__len__",&SpherePack::len,"Get number of spheres in the packing")
 		.def("__getitem__",&SpherePack::getitem,"Get entry at given index, as tuple of center and radius.")
 		.def("__iter__",&SpherePack::getIterator,"Return iterator over spheres.")
-		.def_readwrite("psdScaleExponent",&SpherePack::psdScaleExponent,"Exponent used to scale cummulative distribution function values so that standard uniform distribution is mapped to uniform mass distribution. Theoretically, this value is 3, but that usually overfavors small particles. The default value is 2.5.")
+		.def_readwrite("psdScaleExponent",&SpherePack::psdScaleExponent,"[Deprecated] Defined for compatibility, no effect.")
+		.def_readonly("appliedPsdScaling",&SpherePack::appliedPsdScaling,"A factor between 0 and 1, uniformly applied on all sizes of of the PSD.")
 		;
 	python::class_<SpherePack::_iterator>("SpherePackIterator",python::init<SpherePack::_iterator&>())
 		.def("__iter__",&SpherePack::_iterator::iter)

=== modified file 'scripts/test/psd.py'
--- scripts/test/psd.py	2010-09-09 09:51:23 +0000
+++ scripts/test/psd.py	2011-01-27 17:05:23 +0000
@@ -7,26 +7,30 @@
 from yade import pack
 import pylab
 # PSD given as points of piecewise-linear function
-psdSizes,psdCumm=[.02,.04,.06,.08,.1],[0.,.1,.7,.9,1.]
+psdSizes,psdCumm=[.02,0.04,0.045,.05,.06,.08,.12],[0.,0.1,0.3,0.3,.3,.7,1.]
 pylab.plot(psdSizes,psdCumm,label='precribed mass PSD')
-sp1=pack.SpherePack(); sp1.psdScaleExponent=2 # this exponent is somewhat empirical, see docs of SpherePack
-sp1.makeCloud((0,0,0),(1,1,1),psdSizes=psdSizes,psdCumm=psdCumm,distributeMass=True)
-pylab.semilogx(*sp1.psd(bins=30,mass=True),label='Mass PSD of %d random spheres'%len(sp1))
-pylab.semilogx(*sp1.psd(bins=30,mass=False),label='Size PSD of %d andom spheres'%len(sp1))
+sp0=pack.SpherePack();
+sp0.makeCloud((0,0,0),(1,1,1),psdSizes=psdSizes,psdCumm=psdCumm,distributeMass=True)
+sp1=pack.SpherePack();
+sp1.makeCloud((0,0,0),(1,1,1),psdSizes=psdSizes,psdCumm=psdCumm,distributeMass=True,num=5000)
+sp2=pack.SpherePack();
+sp2.makeCloud((0,0,0),(1,1,1),psdSizes=psdSizes,psdCumm=psdCumm,distributeMass=True,num=20000)
+pylab.semilogx(*sp0.psd(bins=30,mass=True),label='Mass PSD of (free) %d random spheres'%len(sp0))
+pylab.semilogx(*sp1.psd(bins=30,mass=True),label='Mass PSD of (imposed) %d random spheres'%len(sp1))
+pylab.semilogx(*sp2.psd(bins=30,mass=True),label='Mass PSD of (imposed) %d random spheres (scaled down)'%len(sp2))
+
 pylab.legend()
 
-# uniform distribution of size (sp2) and of mass (sp3)
-sp2=pack.SpherePack(); sp2.makeCloud((0,0,0),(1,1,1),rMean=0.03,rRelFuzz=2/3.);
-sp3=pack.SpherePack(); sp3.makeCloud((0,0,0),(1,1,1),rMean=0.03,rRelFuzz=2/3.,distributeMass=True);
+# uniform distribution of size (sp3) and of mass (sp4)
+sp3=pack.SpherePack(); sp3.makeCloud((0,0,0),(1,1,1),rMean=0.03,rRelFuzz=2/3.,distributeMass=False);
+sp4=pack.SpherePack(); sp4.makeCloud((0,0,0),(1,1,1),rMean=0.03,rRelFuzz=2/3.,distributeMass=True);
 pylab.figure()
-pylab.plot(*(sp2.psd(mass=True)+('g',)+sp3.psd(mass=True)+('r',)))
+pylab.plot(*(sp3.psd(mass=True)+('g',)+sp4.psd(mass=True)+('r',)))
 pylab.legend(['Mass PSD of size-uniform distribution','Mass PSD of mass-uniform distribution'])
 
 pylab.figure()
-pylab.plot(*(sp2.psd()+('g',)+sp3.psd()+('r',)))
+pylab.plot(*(sp3.psd(mass=False)+('g',)+sp4.psd(mass=False)+('r',)))
 pylab.legend(['Size PSD of size-uniform distribution','Size PSD of mass-uniform distribution'])
 pylab.show()
 
-pylab.show()
-
-
+pylab.show()
\ No newline at end of file