yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2847
committer: Chiara Modenese <c.modenese@xxxxxxxxx>
branch nick: yade
timestamp: Mon 2011-05-02 18:49:39 +0000
message:
  - For now, add 2d version of particleSD (discrete distribution) - for other type of distributions (like the one generated by makeCloud), something similar can be done.
modified:
  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 'pkg/dem/SpherePack.cpp'
--- pkg/dem/SpherePack.cpp	2011-02-18 16:04:42 +0000
+++ pkg/dem/SpherePack.cpp	2011-05-02 18:49:39 +0000
@@ -299,7 +299,7 @@
 
 // TODO: header, python wrapper, default params
 
-// New code to include the psd giving few points of it
+// Discrete particle size distribution
 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, int seed){
 	vector<Real> numbers;
 	if(!passingIsNotPercentageButCount){
@@ -352,6 +352,59 @@
 	return pack.size();
 }
 
+// 2d function
+long SpherePack::particleSD_2d(Vector2r mn, Vector2r mx, Real rMean, bool periodic, string name, int numSph, const vector<Real>& radii, const vector<Real>& passing, bool passingIsNotPercentageButCount, int seed){
+	vector<Real> numbers;
+	if(!passingIsNotPercentageButCount){
+		Real Vtot=numSph*4./3.*Mathr::PI*pow(rMean,3.); // total volume of the packing (computed with rMean)
+		
+		// calculate number of spheres necessary per each radius to match the wanted psd
+		// passing has to contain increasing values
+		for (size_t i=0; i<radii.size(); i++){
+			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<<"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);
+	}
+
+	static boost::minstd_rand randGen(seed!=0?seed:(int)TimingInfo::getNow(true));
+	static boost::variate_generator<boost::minstd_rand&, boost::uniform_real<> > rnd(randGen, boost::uniform_real<>(0,1));
+
+	const int maxTry=1000;
+	Vector2r size=mx-mn; 
+	//if(periodic)(cellSize=size); in this case, it must be defined in py script as cell needs the third dimension
+	for (int ii=(int)radii.size()-1; ii>=0; ii--){
+		Real r=radii[ii]; // select radius
+		for(int i=0; i<numbers[ii]; i++) { // place as many spheres as required by the psd for the selected radius into the free spot
+			int t;
+			for(t=0; t<maxTry; ++t){
+				Vector3r c;
+				if(!periodic) { for(int axis=0; axis<2; axis++) c[axis]=mn[axis]+r+(size[axis]-2*r)*rnd(); }
+				else { for(int axis=0; axis<2; axis++) c[axis]=mn[axis]+size[axis]*rnd(); }
+				size_t packSize=pack.size(); bool overlap=false;
+				if(!periodic){
+					for(size_t j=0; j<packSize; j++){ if(pow(pack[j].r+r,2) >= (pack[j].c-c).squaredNorm()) { overlap=true; break; } }
+				} else {
+					for(size_t j=0; j<packSize; j++){
+						Vector3r dr;
+						for(int axis=0; axis<2; axis++) dr[axis]=min(cellWrapRel(c[axis],pack[j].c[axis],pack[j].c[axis]+size[axis]),cellWrapRel(pack[j].c[axis],c[axis],c[axis]+size[axis]));
+						if(pow(pack[j].r+r,2)>= dr.squaredNorm()){ overlap=true; break; }
+					}
+				}
+				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]);
+				return i;
+			}
+		}
+	}
+	return pack.size();
+}
+
 long SpherePack::makeClumpCloud(const Vector3r& mn, const Vector3r& mx, const vector<shared_ptr<SpherePack> >& _clumps, bool periodic, int num){
 	// recenter given clumps and compute their margins
 	vector<SpherePack> clumps; /* vector<Vector3r> margins; */ Vector3r boxMargins(Vector3r::Zero()); Real maxR=0;

=== modified file 'pkg/dem/SpherePack.hpp'
--- pkg/dem/SpherePack.hpp	2011-02-10 13:03:39 +0000
+++ pkg/dem/SpherePack.hpp	2011-05-02 18:49:39 +0000
@@ -70,10 +70,10 @@
 
 	// 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, int seed=0);
+	long particleSD_2d(Vector2r min, Vector2r 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, int seed=0);
 
 	long particleSD2(const vector<Real>& radii, const vector<Real>& passing, int numSph, bool periodic=false, Real cloudPorosity=.8, int seed=0);
 
-
 	// 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-02-21 14:45:36 +0000
+++ py/pack/_packSpheres.cpp	2011-05-02 18:49:39 +0000
@@ -23,6 +23,7 @@
 		// 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,python::arg("seed")=0),"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("particleSD2",&SpherePack::particleSD2,(python::arg("radii"),python::arg("passing"),python::arg("numSph"),python::arg("periodic")=false,python::arg("cloudPorosity")=.8,python::arg("seed")=0),"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("particleSD_2d",&SpherePack::particleSD_2d,(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,python::arg("seed")=0),"Create random packing on XY plane, defined 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("aabb",&SpherePack::aabb_py,"Get axis-aligned bounding box coordinates, as 2 3-tuples.")