yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 1720
committer: Václav Šmilauer <vaclav@flux>
branch nick: trunk
timestamp: Sat 2009-08-22 21:26:41 +0200
message:
  PeriIsoCompressor: fix absolute displacements changing proportions, operate on averages instead
  pack.randomDensePack: fix proportions in periodic packs; change residual stress to 1e8 instead of 1e9
modified:
  extra/PeriodicInsertionSortCollider.cpp
  pkg/dem/DataClass/SpherePack.hpp
  py/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 'extra/PeriodicInsertionSortCollider.cpp'
--- extra/PeriodicInsertionSortCollider.cpp	2009-08-22 18:21:25 +0000
+++ extra/PeriodicInsertionSortCollider.cpp	2009-08-22 19:26:41 +0000
@@ -380,7 +380,7 @@
 	const long& step=rb->currentIteration;
 	Vector3r cellSize=rb->cellMax-rb->cellMin; //unused: Real cellVolume=cellSize[0]*cellSize[1]*cellSize[2];
 	Vector3r cellArea=Vector3r(cellSize[1]*cellSize[2],cellSize[0]*cellSize[2],cellSize[0]*cellSize[1]);
-	Real minSize=min(cellSize[0],min(cellSize[1],cellSize[2]));
+	Real minSize=min(cellSize[0],min(cellSize[1],cellSize[2])), maxSize=max(cellSize[0],max(cellSize[1],cellSize[2]));
 	if(minSize<2.1*maxSpan){ throw runtime_error("Minimum cell size is smaller than 2.1*span_of_the_biggest_body! (periodic collider requirement)"); }
 	if(((step%globalUpdateInt)==0) || avgStiffness<0 || sigma[0]<0 || sigma[1]<0 || sigma[2]<0){
 		Vector3r sumForces=Shop::totalForceInVolume(avgStiffness,rb);
@@ -393,13 +393,14 @@
 	// is the stress condition satisfied in all directions?
 	bool allStressesOK=true;
 	if(keepProportions){ // the same algo as below, but operating on quantitites averaged over all dimensions
-		Real sigAvg=(sigma[0]+sigma[1]+sigma[2])/3., avgArea=(cellArea[0]+cellArea[1]+cellArea[2])/3.;
-		Real grow=1e-4*(sigAvg-sigmaGoal)*avgArea/(avgStiffness>0?avgStiffness:1);
-		if(abs(grow)>maxDisplPerStep) grow=Mathr::Sign(grow)*maxDisplPerStep;
-		grow=max(grow,-(minSize-2.1*maxSpan));
-		if(avgStiffness>0) { sigma-=(grow*avgStiffness)*Vector3r::ONE; sigAvg-=grow*avgStiffness; }
+		Real sigAvg=(sigma[0]+sigma[1]+sigma[2])/3., avgArea=(cellArea[0]+cellArea[1]+cellArea[2])/3., avgSize=(cellSize[0]+cellSize[1]+cellSize[2])/3.;
+		Real avgGrow=1e-4*(sigAvg-sigmaGoal)*avgArea/(avgStiffness>0?avgStiffness:1);
+		Real maxToAvg=maxSize/avgSize;
+		if(abs(maxToAvg*avgGrow)>maxDisplPerStep) avgGrow=Mathr::Sign(avgGrow)*maxDisplPerStep/maxToAvg;
+		avgGrow=max(avgGrow,-(minSize-2.1*maxSpan)/maxToAvg);
+		if(avgStiffness>0) { sigma-=(avgGrow*avgStiffness)*Vector3r::ONE; sigAvg-=avgGrow*avgStiffness; }
 		if(abs((sigAvg-sigmaGoal)/sigmaGoal)>5e-3) allStressesOK=false;
-		cellGrow=Vector3r(grow,grow,grow);
+		cellGrow=(avgGrow/avgSize)*cellSize;
 	}
 	else{ // handle each dimension separately
 		for(int axis=0; axis<3; axis++){

=== modified file 'pkg/dem/DataClass/SpherePack.hpp'
--- pkg/dem/DataClass/SpherePack.hpp	2009-08-22 15:05:20 +0000
+++ pkg/dem/DataClass/SpherePack.hpp	2009-08-22 19:26:41 +0000
@@ -84,7 +84,7 @@
 	// transformations
 	void translate(const Vector3r& shift){ FOREACH(Sph& s, pack) s.c+=shift; }
 	void rotate(const Vector3r& axis, Real angle){
-		if(cellSize!=Vector3r::ZERO) { LOG_WARN("Periodicity reset when rotating periodic packing."); cellSize=Vector3r::ZERO; }
+		if(cellSize!=Vector3r::ZERO) { LOG_WARN("Periodicity reset when rotating periodic packing (non-zero cellSize="<<cellSize<<")"); cellSize=Vector3r::ZERO; }
 		Vector3r mid=midPt(); Quaternionr q(axis,angle); FOREACH(Sph& s, pack) s.c=q*(s.c-mid)+mid;
 	}
 	void scale(Real scale){ Vector3r mid=midPt(); cellSize*=scale; FOREACH(Sph& s, pack) {s.c=scale*(s.c-mid)+mid; s.r*=abs(scale); } }

=== modified file 'py/pack.py'
--- py/pack.py	2009-08-22 15:05:20 +0000
+++ py/pack.py	2009-08-22 19:26:41 +0000
@@ -272,7 +272,7 @@
 			if (rRelFuzz==0 and rDev!=0) or (rRelFuzz==0 and rDev!=0) or (rRelFuzz!=0 and abs((rDev-rRelFuzz)/rRelFuzz)>1e-2): memoDbgMsg("REJECT: radius fuzz differs too much (%g, %g desired)"%(rDev,rRelFuzz)); continue # radius fuzz differs too much
 			if isPeri and wantPeri:
 				if spheresInCell>NN: memoDbgMsg("REJECT: Number of spheres in the packing too small"); continue
-				if abs((fullDim[1]/fullDim[0])/(Y/X)-1)>0.3 or abs((fullDim[2]/fullDim[0])/(Z/X)-1)>0.3: memoDbgMsg("REJECT: proportions differ too much from what is desired."); continue
+				if abs((y1/x1)/(Y/X)-1)>0.3 or abs((z1/x1)/(Z/X)-1)>0.3: memoDbgMsg("REJECT: proportions (y/x=%g, z/x=%g) differ too much from what is desired (%g, %g)."%(Y/X,Z/X,y1/x1,z1/x1)); continue
 			else:
 				if (X<fullDim[0] or Y<fullDim[1] or Z<fullDim[2]): memoDbgMsg("REJECT: not large enough"); continue # not large enough
 			memoDbgMsg("ACCEPTED");
@@ -294,11 +294,12 @@
 		#print cloudPorosity,beta,gamma,N100,x1,y1,z1,O.periodicCell
 		#print x1,y1,z1,radius,rRelFuzz
 		num=sp.makeCloud(O.periodicCell[0],O.periodicCell[1],radius,rRelFuzz,spheresInCell,True)
-		O.engines=[BexResetter(),BoundingVolumeMetaEngine([InteractingSphere2AABB()]),PeriodicInsertionSortCollider(),InteractionDispatchers([ef2_Sphere_Sphere_Dem3DofGeom()],[SimpleElasticRelationships()],[Law2_Dem3Dof_Elastic_Elastic()]),PeriIsoCompressor(charLen=radius/5.,stresses=[100e9,1e9],maxUnbalanced=1e-2,doneHook='O.pause();',globalUpdateInt=5),NewtonsDampedLaw(damping=.6)]
+		O.engines=[BexResetter(),BoundingVolumeMetaEngine([InteractingSphere2AABB()]),PeriodicInsertionSortCollider(),InteractionDispatchers([ef2_Sphere_Sphere_Dem3DofGeom()],[SimpleElasticRelationships()],[Law2_Dem3Dof_Elastic_Elastic()]),PeriIsoCompressor(charLen=radius/5.,stresses=[100e9,1e8],maxUnbalanced=1e-2,doneHook='O.pause();',globalUpdateInt=5,keepProportions=True),NewtonsDampedLaw(damping=.6)]
 		for s in sp: O.bodies.append(utils.sphere(s[0],s[1],density=1000))
 		O.dt=utils.PWaveTimeStep()
 		O.run(); O.wait()
 		sp=SpherePack(); sp.fromSimulation()
+		#print 'Resulting cellSize',sp.cellSize,'proportions',sp.cellSize[1]/sp.cellSize[0],sp.cellSize[2]/sp.cellSize[0]
 		# repetition to the required cell size will be done below, after memoizing the result
 	else:
 		assumedFinalDensity=0.6
@@ -328,7 +329,9 @@
 		conn.commit()
 		print "Packing saved to the database",memoizeDb
 	if wantPeri: sp.cellFill(Vector3(fullDim[0],fullDim[1],fullDim[2]))
-	if orientation: sp.rotate(*orientation.ToAxisAngle())
+	if orientation:
+		sp.cellSize=(0,0,0); # reset periodicity to avoid warning when rotating periodic packing
+		sp.rotate(*orientation.ToAxisAngle())
 	return filterSpherePack(predicate,sp,**kw)
 
 # compatibility with the deprecated name, can be removed in the future