yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2954
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
branch nick: trunk
timestamp: Wed 2011-11-02 16:09:47 +0100
message:
  - yet more fixes. The auto-setting of cell.size was not correct (typo in components)
modified:
  pkg/dem/SpherePack.cpp
  pkg/dem/SpherePack.hpp
  py/pack/_packSpheres.cpp
  py/pack/pack.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-11-02 12:42:29 +0000
+++ pkg/dem/SpherePack.cpp	2011-11-02 15:09:47 +0000
@@ -83,6 +83,7 @@
 }
 
 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){
+	isPeriodic = periodic;
 	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
@@ -92,10 +93,12 @@
 	size=mx-mn;
 	if (!hSizeFound) {hSize=size.asDiagonal();}
 	if (hSizeFound && !periodic) LOG_WARN("hSize can be defined only for periodic cells.");
+	Real volume=hSize.determinant();
 	Matrix3r invHsize =hSize.inverse();
-	Real volume=hSize.determinant();
 	Real area=abs(size[0]*size[2]+size[0]*size[1]+size[1]*size[2]);//2 terms will be null if one coordinate is 0, the other is the area
-	if (!volume) LOG_WARN("The volume of the box is null, we will assume that the packing is 2D. If it is not what you want then you defined wrong input values; check that min and max corners are defined correctly.");
+	if (!volume) {
+		if (hSizeFound || !periodic) throw invalid_argument("The box as null volume, this is not supported. One exception is for for flat box defined by min-max in periodic conditions, if hSize argument defines a non-null volume (or if the hSize argument is left undefined).");
+		else LOG_WARN("The volume of the min-max box is null, we will assume that the packing is 2D. If it is not what you want then you defined wrong input values; check that min and max corners are defined correctly.");}
 	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;}
@@ -126,16 +129,15 @@
 			// 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.");
 		}
-		// check the consistency between sizes, num, and poro. If not consistent, scale the sizes down
+		// check the consistency between sizes, num, and poro. If target number will not fit in (1-poro)*volume, scale down particles sizes
 		if (num>1){
 			appliedPsdScaling=1;
 			if(distributeMass) {
-				//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.);
 				//Normalize psdCumm2 so it's between 0 and 1
 				for(size_t i=1; i<psdSizes.size(); i++) psdCumm2[i]/=psdCumm2[psdSizes.size()-1];
 			} else {
-				//if target number will not fit in (1-poro)*volume, scale down particles size
 				double totVol=0;
 				for(size_t i=1; i<psdSizes.size(); i++) totVol+= 4/3*Mathr::PI*(psdCumm[i]-psdCumm[i-1])*num*
 					pow(0.5*(psdSizes[i]+psdSizes[i-1]),3)*(1+pow(0.5*(psdSizes[i]-psdSizes[i-1]),2));
@@ -143,13 +145,12 @@
 				if (volumeRatio>1) appliedPsdScaling=pow(volumeRatio,-1./3.);
 			}
 			if (appliedPsdScaling<1) for(size_t i=0; i<psdSizes.size(); i++) psdRadii[i]*=appliedPsdScaling;
-		}
-		
+		}		
 	}
 	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;
-	if(periodic)(cellSize=size);
+	if(periodic && volume && !hSizeFound)(cellSize=size);
 	Real r=0;
 	for(int i=0; (i<num) || (num<0); i++) {
 		Real norm, rand;

=== modified file 'pkg/dem/SpherePack.hpp'
--- pkg/dem/SpherePack.hpp	2011-05-02 18:49:39 +0000
+++ pkg/dem/SpherePack.hpp	2011-11-02 15:09:47 +0000
@@ -47,9 +47,9 @@
 	};
 	std::vector<Sph> pack;
 	Vector3r cellSize;
-	Real psdScaleExponent;
 	Real appliedPsdScaling;//a scaling factor that can be applied on size distribution
-	SpherePack(): cellSize(Vector3r::Zero()), psdScaleExponent(2.5), appliedPsdScaling(1.){};
+	bool isPeriodic;
+	SpherePack(): cellSize(Vector3r::Zero()), appliedPsdScaling(1.), isPeriodic(0) {};
 	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-10-31 17:43:05 +0000
+++ py/pack/_packSpheres.cpp	2011-11-02 15:09:47 +0000
@@ -30,6 +30,7 @@
 		.def("dim",&SpherePack::dim,"Return dimensions of the packing in terms of aabb(), as a 3-tuple.")
 		.def("center",&SpherePack::midPt,"Return coordinates of the bounding box center.")
 		.def_readwrite("cellSize",&SpherePack::cellSize,"Size of periodic cell; is Vector3(0,0,0) if not periodic. (Change this property only if you know what you're doing).")
+		.def_readwrite("isPeriodic",&SpherePack::isPeriodic,"was the packing generated in periodic boundaries?")
 		.def("cellFill",&SpherePack::cellFill,"Repeat the packing (if periodic) so that the results has dim() >= given size. The packing retains periodicity, but changes cellSize. Raises exception for non-periodic packing.")
 		.def("cellRepeat",&SpherePack::cellRepeat,"Repeat the packing given number of times in each dimension. Periodicity is retained, cellSize changes. Raises exception for non-periodic packing.")
 		.def("relDensity",&SpherePack::relDensity,"Relative packing density, measured as sum of spheres' volumes / aabb volume.\n(Sphere overlaps are ignored.)")
@@ -41,7 +42,6 @@
 		.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,"[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&>())

=== modified file 'py/pack/pack.py'
--- py/pack/pack.py	2011-08-29 18:23:23 +0000
+++ py/pack/pack.py	2011-11-02 15:09:47 +0000
@@ -75,7 +75,7 @@
 	>>> sp.toSimulation(rot=Quaternion((0,0,1),pi/4),color=(0,0,1))
 	[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19]
 
-Periodic properties are transferred to the simulation correctly, including rotation:
+Periodic properties are transferred to the simulation correctly, including rotation (this could be avoided by explicitely passing "hSize=O.cell.hSize" as an argument):
 
 	>>> O.periodic
 	True
@@ -95,9 +95,9 @@
 """
 	if isinstance(rot,Quaternion): rot=rot.toRotationMatrix()
 	assert(isinstance(rot,Matrix3))
+	if self.isPeriodic: O.periodic=True
 	if self.cellSize!=Vector3.Zero:
-		O.periodic=True
-		O.cell.hSize=rot*Matrix3(self.cellSize[0],0,0, 0,self.cellSize[0],0, 0,0,self.cellSize[1])
+		O.cell.hSize=rot*Matrix3(self.cellSize[0],0,0, 0,self.cellSize[1],0, 0,0,self.cellSize[2])
 		O.cell.trsf=Matrix3.Identity
 	if not self.hasClumps():
 		return O.bodies.append([utils.sphere(rot*c,r,**kw) for c,r in self])