yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2975
committer: Jan Stransky <_honzik@xxxxxxxxxx>
branch nick: yade
timestamp: Thu 2011-12-01 13:52:20 +0100
message:
  1) modification of CpmState class
  2) some changes in py/pack/pack.py
modified:
  pkg/dem/ConcretePM.cpp
  pkg/dem/ConcretePM.hpp
  pkg/dem/VTKRecorder.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/ConcretePM.cpp'
--- pkg/dem/ConcretePM.cpp	2011-08-30 18:52:11 +0000
+++ pkg/dem/ConcretePM.cpp	2011-12-01 12:52:20 +0000
@@ -159,7 +159,7 @@
 		BC->crossSection=Mathr::PI*pow(minRad,2);
 		BC->kn=BC->crossSection*BC->E/contGeom->refLength;
 		BC->ks=BC->crossSection*BC->G/contGeom->refLength;
-		BC->epsFracture = isnan(BC->crackOpening)? BC->epsCrackOnset*BC->relDuctility : BC->crackOpening/contGeom->refLength;
+		BC->epsFracture = isnan(BC->crackOpening)? BC->epsCrackOnset*BC->relDuctility : BC->crackOpening/(2*minRad);//*contGeom->refLength;
 	}
 
 	// shorthands
@@ -175,8 +175,8 @@
 		const Real& omegaThreshold(Law2_Dem3DofGeom_CpmPhys_Cpm::omegaThreshold);
 		const Real& epsCrackOnset(BC->epsCrackOnset);
 		Real& relResidualStrength(BC->relResidualStrength);
-		const Real& crackOpening(BC->crackOpening);
-		const Real& relDuctility(BC->relDuctility);
+		//const Real& crackOpening(BC->crackOpening);
+		//const Real& relDuctility(BC->relDuctility);
 		const Real& epsFracture(BC->epsFracture);
 		const int& damLaw(BC->damLaw);
 		const bool& neverDamage(BC->neverDamage);
@@ -442,16 +442,15 @@
 		const Body::id_t id1=I->getId1(), id2=I->getId2();
 		Dem3DofGeom* geom=YADE_CAST<Dem3DofGeom*>(I->geom.get());
 
-		Matrix3r stressTimesV = Matrix3r::Zero();
 		const Vector3r& n= geom->normal;
 		const Real& Fn = phys->Fn;
 		const Vector3r& Fs = phys->Fs;
-		//stressTimesV[i,j] += geom->refLength*(Fn*n[i]*n[j]+0.5*(Fs[i]*n[j]+Fs[j]*n[i]));
-		//stressTimesV += geom->refLength*(Fn*outer(n,n)+.5*(outer(Fs,n)+outer(n,Fs)));
-		stressTimesV += geom->refLength*(Fn*n*n.transpose()+.5*(Fs*n.transpose()+n*Fs.transpose()));
+		//stress[i,j] += geom->refLength*(Fn*n[i]*n[j]+0.5*(Fs[i]*n[j]+Fs[j]*n[i]));
+		//stress += geom->refLength*(Fn*outer(n,n)+.5*(outer(Fs,n)+outer(n,Fs)));
+		Matrix3r stress = geom->refLength*(Fn*n*n.transpose()+.5*(Fs*n.transpose()+n*Fs.transpose()));
 		
-		bodyStats[id1].stressTimesV += stressTimesV;
-		bodyStats[id2].stressTimesV += stressTimesV;
+		bodyStats[id1].stress += stress;
+		bodyStats[id2].stress += stress;
 		bodyStats[id1].nLinks++; bodyStats[id2].nLinks++;
 		
 		if(!phys->isCohesive) continue;
@@ -467,8 +466,7 @@
 		// add damaged contacts that have already been deleted
 		CpmState* state=dynamic_cast<CpmState*>(B->state.get());
 		if(!state) continue;
-		state->stressTimesV=bodyStats[id].stressTimesV;
-		//state->tau=bodyStats[id].tau;
+		state->stress=bodyStats[id].stress;
 		int cohLinksWhenever=bodyStats[id].nCohLinks+state->numBrokenCohesive;
 		if(cohLinksWhenever>0){
 			state->normDmg=(bodyStats[id].dmgSum+state->numBrokenCohesive)/cohLinksWhenever;
@@ -480,6 +478,10 @@
 		else { state->normDmg=0; state->normEpsPl=0;}
 		B->shape->color=Vector3r(state->normDmg,1-state->normDmg,B->state->blockedDOFs==State::DOF_ALL?0:1);
 		nAvgRelResidual+=0.5*state->numBrokenCohesive; // add half or broken interactions, other body has the other half
+		Sphere* sphere=dynamic_cast<Sphere*>(B->shape.get());
+		if(!sphere) continue;
+		Real& r = sphere->radius;
+		state->stress=bodyStats[id].stress/(4/3.*Mathr::PI*r*r*r)*.5;
 	}
 	avgRelResidual/=nAvgRelResidual;
 }

=== modified file 'pkg/dem/ConcretePM.hpp'
--- pkg/dem/ConcretePM.hpp	2011-08-30 18:52:11 +0000
+++ pkg/dem/ConcretePM.hpp	2011-12-01 12:52:20 +0000
@@ -66,9 +66,7 @@
 		((Real,normDmg,0,,"Average damage including already deleted contacts (it is really not damage, but 1-relResidualStrength now)"))
 		((Real,epsPlBroken,0,,"Plastic strain on contacts already deleted (bogus values)"))
 		((Real,normEpsPl,0,,"Sum of plastic strains normalized by number of contacts (bogus values)"))
-		//((Vector3r,sigma,Vector3r::Zero(),,"Normal stresses on the particle"))
-		//((Vector3r,tau,Vector3r::Zero(),,"Shear stresses on the particle."))
-		((Matrix3r,stressTimesV,Matrix3r::Zero(),,"Stress tensor on the particle (multiplied by volume surrounded it). To get actual stress, divide this matrix with the volume (for dense packing something like (4/3.*pi*r*r*r/0.62)")),
+		((Matrix3r,stress,Matrix3r::Zero(),,"Stress tensor of the spherical particle (under assumption that particle volume = pi*r*r*r*4/3.). To get actual stress, multiply this value by packing fraction (for random dense packing something like 0.63)")),
 		/*ctor*/ createIndex();
 	);
 	REGISTER_CLASS_INDEX(CpmState,State);
@@ -87,7 +85,7 @@
 		((bool,neverDamage,false,,"If true, no damage will occur (for testing only)."))
 		((Real,epsCrackOnset,NaN,,"Limit elastic strain [-]"))
 		((Real,crackOpening,NaN,,"Crack opening when the crack is fully broken in tension. [m]"))
-		((Real,relDuctility,NaN,,"Deprecated"))
+		((Real,relDuctility,NaN,,"relative ductility of bonds in normal direction"))
 		((int,damLaw,1,,"Law for gamage evolution in uniaxial tension. 0 for linear stress-strain softening branch, 1 for exponential damage evolution law"))
 		((Real,dmgTau,((void)"deactivated if negative",-1),,"Characteristic time for normal viscosity. [s]"))
 		((Real,dmgRateExp,0,,"Exponent for normal viscosity function. [-]"))
@@ -127,7 +125,7 @@
 			((Real,crossSection,NaN,,"equivalent cross-section associated with this contact [m²]"))
 			((Real,epsCrackOnset,NaN,,"strain at which the material starts to behave non-linearly"))
 			((Real,crackOpening,NaN,,"Crack opening (extansion of the bond) when the bond is fully broken in tension. [m]"))
-			((Real,relDuctility,NaN,,"Deprecated, use :yref:`CpmMat::crackOpening` instead"))
+			((Real,relDuctility,NaN,,"Relative ductility of bonds in normal direction"))
 			((Real,epsFracture,NaN,,"strain at which the bond is fully broken [-]"))
 			((Real,dmgTau,-1,,"characteristic time for damage (if non-positive, the law without rate-dependence is used)"))
 			((Real,dmgRateExp,0,,"exponent in the rate-dependent damage evolution"))
@@ -249,7 +247,7 @@
 #endif
 
 class CpmStateUpdater: public PeriodicEngine {
-	struct BodyStats{ int nCohLinks; int nLinks; Real dmgSum, epsPlSum; Matrix3r stressTimesV; BodyStats(): nCohLinks(0), nLinks(0), dmgSum(0.), epsPlSum(0.), stressTimesV(Matrix3r::Zero()) {} };
+	struct BodyStats{ int nCohLinks; int nLinks; Real dmgSum, epsPlSum; Matrix3r stress; BodyStats(): nCohLinks(0), nLinks(0), dmgSum(0.), epsPlSum(0.), stress(Matrix3r::Zero()) {} };
 	public:
 		virtual void action(){ update(scene); }
 		void update(Scene* rb=NULL);

=== modified file 'pkg/dem/VTKRecorder.cpp'
--- pkg/dem/VTKRecorder.cpp	2011-08-29 18:23:23 +0000
+++ pkg/dem/VTKRecorder.cpp	2011-12-01 12:52:20 +0000
@@ -306,7 +306,7 @@
 				
 				if (recActive[REC_CPM]){
 					cpmDamage->InsertNextValue(YADE_PTR_CAST<CpmState>(b->state)->normDmg);
-					const Matrix3r& ss=YADE_PTR_CAST<CpmState>(b->state)->stressTimesV;
+					const Matrix3r& ss=YADE_PTR_CAST<CpmState>(b->state)->stress;
 					//float s[3]={ss[0],ss[1],ss[2]};
 					float s[9]={ss[0,0],ss[0,1],ss[0,2],ss[1,0],ss[1,1],ss[1,2],ss[2,0],ss[2,1],ss[2,2]};
 					cpmStress->InsertNextTupleValue(s);

=== modified file 'py/pack/pack.py'
--- py/pack/pack.py	2011-11-02 15:09:47 +0000
+++ py/pack/pack.py	2011-12-01 12:52:20 +0000
@@ -302,7 +302,7 @@
 			if predicate(s[0],s[1]): ret+=[utils.sphere(s[0],radius=s[1],**kw)]
 		return ret
 
-def _memoizePacking(memoizeDb,sp,radius,rRelFuzz,wantPeri,fullDim):
+def _memoizePacking(memoizeDb,sp,radius,rRelFuzz,wantPeri,fullDim,noPrint=False):
 	if not memoizeDb: return
 	import cPickle,sqlite3,time,os
 	if os.path.exists(memoizeDb):
@@ -317,16 +317,16 @@
 	c.execute('insert into packings values (?,?,?,?,?,?,?,?,?)',(radius,rRelFuzz,packDim[0],packDim[1],packDim[2],len(sp),time.time(),wantPeri,packBlob,))
 	c.close()
 	conn.commit()
-	print "Packing saved to the database",memoizeDb
+	if not noPrint: print "Packing saved to the database",memoizeDb
 
-def _getMemoizedPacking(memoizeDb,radius,rRelFuzz,x1,y1,z1,fullDim,wantPeri,fillPeriodic,spheresInCell,memoDbg=False):
+def _getMemoizedPacking(memoizeDb,radius,rRelFuzz,x1,y1,z1,fullDim,wantPeri,fillPeriodic,spheresInCell,memoDbg=False,noPrint=False):
 	"""Return suitable SpherePack read from *memoizeDb* if found, None otherwise.
 		
 		:param fillPeriodic: whether to fill fullDim by repeating periodic packing
 		:param wantPeri: only consider periodic packings
 	"""
 	import os,os.path,sqlite3,time,cPickle,sys
-	if memoDbg:
+	if memoDbg and not noPrint:
 		def memoDbgMsg(s): print s
 	else:
 		def memoDbgMsg(s): pass
@@ -352,7 +352,7 @@
 		else:
 			if (X<fullDim[0] or Y<fullDim[1] or Z<fullDim[2]): memoDbgMsg("REJECT: not large enough"); continue # not large enough
 		memoDbgMsg("ACCEPTED");
-		print "Found suitable packing in %s (radius=%g±%g,N=%g,dim=%g×%g×%g,%s,scale=%g), created %s"%(memoizeDb,R,rDev,NN,X,Y,Z,"periodic" if isPeri else "non-periodic",scale,time.asctime(time.gmtime(timestamp)))
+		if not noPrint: print "Found suitable packing in %s (radius=%g±%g,N=%g,dim=%g×%g×%g,%s,scale=%g), created %s"%(memoizeDb,R,rDev,NN,X,Y,Z,"periodic" if isPeri else "non-periodic",scale,time.asctime(time.gmtime(timestamp)))
 		c.execute('select pack from packings where timestamp=?',(timestamp,))
 		sp=SpherePack(cPickle.loads(str(c.fetchone()[0])))
 		sp.scale(scale);
@@ -435,7 +435,8 @@
 		# x1,y1,z1 already computed above
 		sp=SpherePack()
 		O.periodic=True
-		O.cell.refSize=(x1,y1,z1)
+		#O.cell.refSize=(x1,y1,z1)
+		O.cell.setBox((x1,y1,z1))
 		#print cloudPorosity,beta,gamma,N100,x1,y1,z1,O.cell.refSize
 		#print x1,y1,z1,radius,rRelFuzz
 		O.materials.append(FrictMat(young=3e10,density=2400))
@@ -468,7 +469,7 @@
 		sp.rotate(*orientation.toAxisAngle())
 	return filterSpherePack(predicate,sp,material=material,color=color,returnSpherePack=returnSpherePack)
 
-def randomPeriPack(radius,initSize,rRelFuzz=0.0,memoizeDb=None):
+def randomPeriPack(radius,initSize,rRelFuzz=0.0,memoizeDb=None,noPrint=False):
 	"""Generate periodic dense packing.
 
 	A cell of initSize is stuffed with as many spheres as possible, then we run periodic compression with PeriIsoCompressor, just like with
@@ -481,12 +482,13 @@
 	:return: SpherePack object, which also contains periodicity information.
 	"""
 	from math import pi
-	sp=_getMemoizedPacking(memoizeDb,radius,rRelFuzz,initSize[0],initSize[1],initSize[2],fullDim=Vector3(0,0,0),wantPeri=True,fillPeriodic=False,spheresInCell=-1,memoDbg=True)
+	sp=_getMemoizedPacking(memoizeDb,radius,rRelFuzz,initSize[0],initSize[1],initSize[2],fullDim=Vector3(0,0,0),wantPeri=True,fillPeriodic=False,spheresInCell=-1,memoDbg=True,noPrint=noPrint)
 	if sp: return sp
 	O.switchScene(); O.resetThisScene()
 	sp=SpherePack()
 	O.periodic=True
-	O.cell.refSize=initSize
+	#O.cell.refSize=initSize
+	O.cell.setBox(initSize)
 	sp.makeCloud(Vector3().Zero,O.cell.refSize,radius,rRelFuzz,-1,True)
 	O.engines=[ForceResetter(),InsertionSortCollider([Bo1_Sphere_Aabb()],nBins=2,verletDist=.05*radius),InteractionLoop([Ig2_Sphere_Sphere_Dem3DofGeom()],[Ip2_FrictMat_FrictMat_FrictPhys()],[Law2_Dem3DofGeom_FrictPhys_CundallStrack()]),PeriIsoCompressor(charLen=2*radius,stresses=[-100e9,-1e8],maxUnbalanced=1e-2,doneHook='O.pause();',globalUpdateInt=20,keepProportions=True),NewtonIntegrator(damping=.8)]
 	O.materials.append(FrictMat(young=30e9,frictionAngle=.1,poisson=.3,density=1e3))
@@ -496,7 +498,7 @@
 	O.run(); O.wait()
 	ret=SpherePack()
 	ret.fromSimulation()
-	_memoizePacking(memoizeDb,ret,radius,rRelFuzz,wantPeri=True,fullDim=Vector3(0,0,0)) # fullDim unused
+	_memoizePacking(memoizeDb,ret,radius,rRelFuzz,wantPeri=True,fullDim=Vector3(0,0,0),noPrint=noPrint) # fullDim unused
 	O.switchScene()
 	return ret