yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2669
committer: Václav Šmilauer <eu@xxxxxxxx>
branch nick: yade
timestamp: Fri 2011-01-21 12:28:22 +0100
message:
  1. Add a more straightforward SpherePack.ParticleSD2 method, adjust the tutorial
modified:
  doc/sphinx/tutorial/06-periodic-triaxial-test.py
  pkg/dem/SpherePack.cpp
  pkg/dem/SpherePack.hpp
  py/pack/_packSpheres.cpp


--
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 'doc/sphinx/tutorial/06-periodic-triaxial-test.py'
--- doc/sphinx/tutorial/06-periodic-triaxial-test.py	2011-01-20 16:35:46 +0000
+++ doc/sphinx/tutorial/06-periodic-triaxial-test.py	2011-01-21 11:28:22 +0000
@@ -35,12 +35,8 @@
 	sp.makeCloud((0,0,0),(2,2,2),rMean=.1,rRelFuzz=.3,periodic=True)
 elif 0:
 	## per-fraction distribution
-	## rMean & numSph needed by the algorithm (name is useless but required)
 	## passing: cummulative percentage
-	sp.particleSD((0,0,0),(2,2,2),
-		rMean=.12,name='',
-		radii=[.09,.1,.2],passing=[40,80,100],periodic=True,numSph=300
-	)
+	sp.particleSD2(radii=[.09,.1,.2],passing=[40,80,100],periodic=True,numSph=1000)
 else:
 	## create packing from clumps
 	# configuration of one clump
@@ -70,7 +66,7 @@
 		# call this function when goal is reached and the packing is stable
 		doneHook='compactionFinished()'
 	),
-	PyRunner(command='addPlotData()',iterPeriod=200),
+	PyRunner(command='addPlotData()',iterPeriod=100),
 ]
 O.dt=.5*utils.PWaveTimeStep()
 
@@ -89,7 +85,6 @@
 plot.plots={'i':('unbalanced',),'i ':('sxx','syy','szz'),' i':('exx','eyy','ezz'),
 	# energy plot
 	' i ':(O.energy.keys,None,'Etot'),
-	'snapshot':None
 }
 # show the plot
 plot.plot()

=== modified file 'pkg/dem/SpherePack.cpp'
--- pkg/dem/SpherePack.cpp	2011-01-09 16:34:50 +0000
+++ pkg/dem/SpherePack.cpp	2011-01-21 11:28:22 +0000
@@ -217,6 +217,52 @@
 	return py::make_tuple(edges,cumm);
 }
 
+/* possible enhacement: proportions parameter, so that the domain is not cube, but box with sides having given proportions */
+long SpherePack::particleSD2(const vector<Real>& radii, const vector<Real>& passing, int numSph, bool periodic, Real cloudPorosity){
+	typedef Eigen::Matrix<Real,Eigen::Dynamic,Eigen::Dynamic> MatrixXr;
+	typedef Eigen::Matrix<Real,Eigen::Dynamic,1> VectorXr;
+	
+	int dim=radii.size()+1;
+	if(passing.size()!=radii.size()) throw std::invalid_argument("SpherePack.particleSD2: radii and passing must have the same length.");
+	MatrixXr M=MatrixXr::Zero(dim,dim);
+	VectorXr rhs=VectorXr::Zero(dim);
+	/*
+	
+	We know percentages for each fraction (Δpi) and their radii (ri), and want to find
+	the number of sphere for each fraction Ni and total solid volume Vs. For each fraction,
+	we know that the volume is equal to Ni*(4/3*πri³), which must be equal to Vs*Δpi (Δpi is
+	relative solid volume of the i-th fraction).
+
+	The last equation says that total number of particles (sum of fractions) is equal to N,
+	which is the total number of particles requested by the user.
+
+	   N1     N2     N3    Vs       rhs
+
+	4/3πr₁³   0      0     -Δp₁   | 0
+	  0     4/3πr₂³  0     -Δp₂   | 0
+	  0       0    4/3πr₃³ -Δp₃   | 0
+     1       1      1      0     | N
+
+	*/
+	for(int i=0; i<dim-1; i++){
+		M(i,i)=(4/3.)*Mathr::PI*pow(radii[i],3);
+		M(i,dim-1)=-(passing[i]-(i>0?passing[i-1]:0))/100.;
+		M(dim-1,i)=1;
+	}
+	rhs[dim-1]=numSph;
+	// NumsVs=M^-1*rhs: number of spheres and volume of solids
+	VectorXr NumsVs(dim); NumsVs=M.inverse()*rhs;
+	Real Vs=NumsVs[dim-1]; // total volume of solids
+	Real Vtot=Vs/(1-cloudPorosity); // total volume of cell containing the packing
+	Vector3r cellSize=pow(Vtot,1/3.)*Vector3r().Ones(); // for now, assume always cubic sample
+	Real rMean=pow(Vs/(numSph*(4/3.)*Mathr::PI),1/3.); // make rMean such that particleSD will compute the right Vs (called a bit confusingly Vtot anyway) inversely
+	// cerr<<"Vs="<<Vs<<", Vtot="<<Vtot<<", rMean="<<rMean<<endl;
+	// cerr<<"cellSize="<<cellSize<<", rMean="<<rMean<<", numSph="<<numSph<<endl;
+	return particleSD(Vector3r::Zero(),cellSize,rMean,periodic,"",numSph,radii,passing,false);
+};
+
+// TODO: header, python wrapper, default params
+
 // New code to include the psd giving few points of it
 long SpherePack::particleSD(Vector3r mn, Vector3r mx, Real rMean, bool periodic, string name, int numSph, const vector<Real>& radii, const vector<Real>& passing, bool passingIsNotPercentageButCount){
 	vector<Real> numbers;
@@ -229,12 +275,7 @@
 			Real volS=4./3.*Mathr::PI*pow(radii[i],3.);
 			if (i==0) {numbers.push_back(passing[i]/100.*Vtot/volS);}
 			else {numbers.push_back((passing[i]-passing[i-1])/100.*Vtot/volS);} // 
-			
-			cout<<"numbers["<<i<<"] = "<<numbers[i]<<endl;
-			cout<<"radii["<<i<<"] = "<<radii[i]<<endl;
-			cout<<"vol tot = "<<Vtot<<endl;
-			cout<<"v_sphere = "<<volS<<endl;
-			cout<<"passing["<<i<<"] = "<<passing[i]<<endl;
+			cout<<"fraction #"<<i<<" ("<<passing[i]<<"%, r="<<radii[i]<<"): "<<numbers[i]<<" spheres, fraction/cloud volumes "<<volS<<"/"<<Vtot<<endl;
 		}
 	} else {
 		FOREACH(Real p, passing) numbers.push_back(p);
@@ -267,7 +308,7 @@
 				if(!overlap) { pack.push_back(Sph(c,r)); break; }
 			}
 			if (t==maxTry) {
-				if(numbers[ii]>0) LOG_WARN("Exceeded "<<maxTry<<" tries to insert non-overlapping sphere to packing. Only "<<i<<" spheres was added, although you requested "<<numbers[ii]<<" with radius"<<radii[ii]);
+				if(numbers[ii]>0) LOG_WARN("Exceeded "<<maxTry<<" tries to insert non-overlapping sphere to packing. Only "<<i<<" spheres was added, although you requested "<<numbers[ii]<<" with radius "<<radii[ii]);
 				return i;
 			}
 		}

=== modified file 'pkg/dem/SpherePack.hpp'
--- pkg/dem/SpherePack.hpp	2010-11-12 08:03:16 +0000
+++ pkg/dem/SpherePack.hpp	2011-01-21 11:28:22 +0000
@@ -70,6 +70,9 @@
 	// function to generate the packing based on a given psd
 	long particleSD(Vector3r min, Vector3r max, Real rMean, bool periodic=false, string name="", int numSph=400, const vector<Real>& radii=vector<Real>(), const vector<Real>& passing=vector<Real>(),bool passingIsNotPercentageButCount=false);
 
+	long particleSD2(const vector<Real>& radii, const vector<Real>& passing, int numSph, bool periodic=false, Real cloudPorosity=.8);
+
+
 	// interpolate a variable with power distribution (exponent -3) between two margin values, given uniformly distributed x∈(0,1)
 	Real pow3Interp(Real x,Real a,Real b){ return pow(x*(pow(b,-2)-pow(a,-2))+pow(a,-2),-1./2); }
 

=== modified file 'py/pack/_packSpheres.cpp'
--- py/pack/_packSpheres.cpp	2011-01-21 08:14:28 +0000
+++ py/pack/_packSpheres.cpp	2011-01-21 11:28:22 +0000
@@ -17,10 +17,11 @@
 		.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("psd",&SpherePack::psd,(python::arg("bins")=10,python::arg("mass")=false),"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``.")
+		.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.")
-		.def("makeClumpCloud",&SpherePack::makeClumpCloud,(python::arg("minCorner"),python::arg("maxCorner"),python::arg("clumps"),python::arg("periodic")=false,python::arg("num")=-1),"Create random loose packing of clumps within box given by *minCorner* and *maxCorner*. Clumps are selected with equal probability. At most *num* clumps will be positioned if *num* is positive; otherwise, as many clumps as possible will be put in space, until maximum number of attempts to place a new clump randomly is attained.")
+		.def("particleSD2",&SpherePack::particleSD2,(python::arg("radii"),python::arg("passing"),python::arg("numSph"),python::arg("periodic")=false,python::arg("cloudPorosity")=.8),"Create random packing following the given particle size distribution (radii and volume/mass passing for each fraction) and total number of particles *numSph*. The cloud size (periodic or aperiodic) is computed from the PSD and is always cubic.")
+		.def("makeClumpCloud",&SpherePack::makeClumpCloud,(python::arg("minCorner"),python::arg("maxCorner"),python::arg("clumps"),python::arg("periodic")=false,python::arg("num")=-1),"Create random loose packing of clumps within box given by *minCorner* and *maxCorner*. Clumps are selected with equal probability. At most *num* clumps will be positioned if *num* is positive; otherwise, as many clumps as possible will be put in space, until maximum number of attemps to place a new clump randomly is attained.")
 		//
 		.def("aabb",&SpherePack::aabb_py,"Get axis-aligned bounding box coordinates, as 2 3-tuples.")
 		.def("dim",&SpherePack::dim,"Return dimensions of the packing in terms of aabb(), as a 3-tuple.")