yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2724
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
branch nick: yade
timestamp: Thu 2011-02-10 14:03:39 +0100
message:
  - add hSize parameter to makeCloud so that it can generate arbitrary parallelepipeds (previously working only by luck on a few special geometries)
  - documentation
  - usage in script, which uses now more funny shape (try it!) 
modified:
  pkg/dem/SpherePack.cpp
  pkg/dem/SpherePack.hpp
  py/pack/_packSpheres.cpp
  scripts/test/periodic-triax-settingHsize.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-02-09 09:46:15 +0000
+++ pkg/dem/SpherePack.cpp	2011-02-10 13:03:39 +0000
@@ -82,12 +82,18 @@
 	if(scene->isPeriodic) { cellSize=scene->cell->getSize(); }
 }
 
-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, int seed){
+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, int seed, Matrix3r hSize){
 	static boost::minstd_rand randGen(seed!=0?seed:(int)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; // 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();
+	vector<Real> psdCumm2; // psdCumm but dimensionally transformed to match mass distribution	
+	Vector3r size;
+	bool hSizeFound =(hSize!=Matrix3r::Zero());//is hSize passed to the function?
+	if (!hSizeFound) {size=mx-mn; hSize=size.asDiagonal();}
+	if (hSizeFound && !periodic) throw invalid_argument("hSize can be defined only for periodic cells.");
+	Matrix3r invHsize =hSize.inverse();
+	Real volume=hSize.determinant();
+	if (!volume) throw invalid_argument("The box defined as null volume. Define at least maxCorner of the box, or hSize if periodic.");
 	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 rMean>0."); porosity=0.5;}
@@ -155,14 +161,22 @@
 		for(t=0; t<maxTry; ++t){
 			Vector3r c;
 			if(!periodic) { for(int axis=0; axis<3; axis++) c[axis]=mn[axis]+r+(size[axis]-2*r)*rnd(); }
-			else { for(int axis=0; axis<3; axis++) c[axis]=mn[axis]+size[axis]*rnd(); }
+			else { 	for(int axis=0; axis<3; axis++) c[axis]=rnd();//coordinates in [0,1]
+				c=mn+hSize*c;}//coordinates in reference frame (inside the base cell)
 			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 {
+			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<3; 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]));
+					Vector3r dr=Vector3r::Zero();
+					if (!hSizeFound) {//The box is axis-aligned, use the wrap methods
+						for(int axis=0; axis<3; 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]));
+					} else {//not aligned, find closest neighbor in a cube of size 1, then transform distance to cartesian coordinates
+						Vector3r c1c2=invHsize*(pack[j].c-c);
+						for(int axis=0; axis<3; axis++){
+							if (abs(c1c2[axis])<abs(c1c2[axis] - Mathr::Sign(c1c2[axis]))) dr[axis]=c1c2[axis];
+							else dr[axis] = c1c2[axis] - Mathr::Sign(c1c2[axis]);}
+						dr=hSize*dr;//now in cartesian coordinates
+					}
 					if(pow(pack[j].r+r,2)>= dr.squaredNorm()){ overlap=true; break; }
 				}
 			}
@@ -172,9 +186,9 @@
 			if(num>0) {
 				if (mode!=RDIST_RMEAN) {
 					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");
+					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 makeCloud(mn, mx, -1., rRelFuzz, num, periodic, nextPoro, psdSizes, psdCumm, distributeMass,seed,hSize);}
 				else LOG_WARN("Exceeded "<<maxTry<<" tries to insert non-overlapping sphere to packing. Only "<<i<<" spheres was added, although you requested "<<num<<".");
 			}
 			return i;}

=== modified file 'pkg/dem/SpherePack.hpp'
--- pkg/dem/SpherePack.hpp	2011-01-31 11:38:54 +0000
+++ pkg/dem/SpherePack.hpp	2011-02-10 13:03:39 +0000
@@ -63,7 +63,7 @@
 	void fromSimulation();
 
 	// random generation; if num<0, insert as many spheres as possible; if porosity>0, recompute meanRadius (porosity>0.65 recommended) and try generating this porosity with num spheres.
-	long makeCloud(Vector3r min, Vector3r max, Real rMean=-1, Real rFuzz=0, int num=-1, bool periodic=false, Real porosity=-1, const vector<Real>& psdSizes=vector<Real>(), const vector<Real>& psdCumm=vector<Real>(), bool distributeMass=false, int seed=0);
+	long makeCloud(Vector3r min, Vector3r max, Real rMean=-1, Real rFuzz=0, int num=-1, bool periodic=false, Real porosity=-1, const vector<Real>& psdSizes=vector<Real>(), const vector<Real>& psdCumm=vector<Real>(), bool distributeMass=false, int seed=0, Matrix3r hSize=Matrix3r::Zero());
 	// return number of piece for x in piecewise function defined by cumm with non-decreasing elements ∈(0,1)
 	// norm holds normalized coordinate withing the piece
 	int psdGetPiece(Real x, const vector<Real>& cumm, Real& norm);

=== modified file 'py/pack/_packSpheres.cpp'
--- py/pack/_packSpheres.cpp	2011-01-31 11:38:54 +0000
+++ py/pack/_packSpheres.cpp	2011-02-10 13:03:39 +0000
@@ -2,6 +2,7 @@
 
 #include<yade/pkg/dem/SpherePack.hpp>
 #include<yade/lib/pyutil/doc_opts.hpp>
+#include<yade/lib/base/Math.hpp>
 
 BOOST_PYTHON_MODULE(_packSpheres){
 	python::scope().attr("__doc__")="Creation, manipulation, IO for generic sphere packings.";
@@ -14,10 +15,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")=0.5,python::arg("psdSizes")=vector<Real>(),python::arg("psdCumm")=vector<Real>(),python::arg("distributeMass")=false,python::arg("seed")=0),"Create random loose packing enclosed in box."
+		.def("makeCloud",&SpherePack::makeCloud,(python::arg("minCorner")=Vector3r(Vector3r::Zero()),python::arg("maxCorner")=Vector3r(Vector3r::Zero()),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,python::arg("seed")=0,python::arg("hSize")=Matrix3r(Matrix3r::Zero())),"Create random loose packing enclosed in a parallelepiped."
 		"\nSphere radius distribution can be specified using one of the following ways:\n\n#. *rMean*, *rRelFuzz* and *num* gives uniform radius distribution in *rMean* (1 ± *rRelFuzz* ). Less than *num* spheres can be generated if it is too high.\n#. *rRelFuzz*, *num* and (optional) *porosity*, which estimates mean radius so that *porosity* is attained at the end.  *rMean* must be less than 0 (default). *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:param seed: number used to initialize the random number generator.\n:returns: number of created spheres, which can be lower than *num* depending on the method used.\n")
+		"\n\n:param Vector3 minCorner: lower corner of an axis-aligned box\n:param Vector3 maxCorner: upper corner of an axis-aligned box\n:param Matrix3 hSize: base vectors of a generalized box (arbitrary parallelepiped, typically :yref:`Cell::hSize`), superseeds minCorner and maxCorner if defined. For periodic boundaries only.\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:param seed: number used to initialize the random number generator.\n:returns: number of created spheres, which can be lower than *num* depending on the method used.\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,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.")

=== modified file 'scripts/test/periodic-triax-settingHsize.py'
--- scripts/test/periodic-triax-settingHsize.py	2011-02-09 09:46:15 +0000
+++ scripts/test/periodic-triax-settingHsize.py	2011-02-10 13:03:39 +0000
@@ -3,12 +3,13 @@
 "Demonstrate the compression of a periodic cell with non-trivial initial geometry."
 from yade import *
 from yade import pack,log,qt
-log.setLevel('PeriTriaxController',log.TRACE)
 O.periodic=True
 
-O.cell.hSize=Matrix3(0.1,0.1,0, 0,0.2,0, 0,0,0.1)
+O.cell.hSize=Matrix3(1.0, -0.15, -0.10,
+		     -0.2 ,1.5, 0.3,
+		    0.3, -0.3, 1.0)
 sp=pack.SpherePack()
-num=sp.makeCloud(minCorner=Vector3().Zero, maxCorner=(0.1,0.2,0.1), rMean=-0.01,rRelFuzz=.2, num=500,periodic=True, porosity=0.52,distributeMass=False)
+num=sp.makeCloud(hSize=O.cell.hSize, rMean=-0.01,rRelFuzz=.2, num=500,periodic=True, porosity=0.52,distributeMass=False)
 O.bodies.append([utils.sphere(s[0],s[1]) for s in sp])
 
 
@@ -25,7 +26,7 @@
 	NewtonIntegrator(damping=.2),
 	#PyRunner(iterPeriod=500,command='print "strain: ",triax.strain," stress: ",triax.stress')
 ]
-O.dt=utils.PWaveTimeStep()
+O.dt=0.5*utils.PWaveTimeStep()
 qt.View()
 
 phase=0
@@ -35,15 +36,12 @@
 		print 'Here we are: stress',triax.stress,'strain',triax.strain
 		#Here we reset the transformation, the compressed packing corresponds to null strain
 		O.cell.trsf=Matrix3.Identity
-		print 'Now εzz will go from 0 to .4 while σxx and σyy will be kept the same.'
-		triax.stressMask=3
-		triax.goal=(-1e4,-1e4,-0.4)
+		print 'Now εxx will go from 0 to .4 while σzz and σyy will be kept the same.'
+		triax.stressMask=6
+		triax.goal=(-0.4,-1e4,-1e4)
 
 		phase+=1
 	elif phase==1:
 		print 'Here we are: stress',triax.stress,'strain',triax.strain
 		print 'Done, pausing now.'
 		O.pause()
-		
-	
-