yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2911
committer: Jan Stransky <honzik.stransky@xxxxxxxxx>
branch nick: yade
timestamp: Mon 2011-08-29 20:23:23 +0200
message:
  Change of CpmMat.relDuctility to CpmMat.crackOpening (making Cpm model more particle size independent)
modified:
  examples/concrete/periodic.py
  examples/concrete/uniax.py
  pkg/dem/ConcretePM.cpp
  pkg/dem/ConcretePM.hpp
  pkg/dem/VTKRecorder.cpp
  py/eudoxos.py
  py/export.py
  py/pack/pack.py
  scripts/test/cpm-dem3dof-scgeom.py
  scripts/test/peri3dController_example1.py
  scripts/test/peri3dController_shear.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 'examples/concrete/periodic.py'
--- examples/concrete/periodic.py	2011-01-20 14:55:15 +0000
+++ examples/concrete/periodic.py	2011-08-29 18:23:23 +0000
@@ -41,7 +41,7 @@
 	sigmaT=3.5e6,
 	frictionAngle=atan(0.8),
 	epsCrackOnset=1e-4,
-	relDuctility=30,
+	crackOpening=1e-6,
 
 	intRadius=1.5,
 	dtSafety=.4,
@@ -71,7 +71,7 @@
 # load the packing (again);
 #
 import cPickle as pickle
-concreteId=O.materials.append(CpmMat(young=young,frictionAngle=frictionAngle,poisson=poisson,density=4800,sigmaT=sigmaT,relDuctility=relDuctility,epsCrackOnset=epsCrackOnset,G_over_E=G_over_E,isoPrestress=isoPrestress))
+concreteId=O.materials.append(CpmMat(young=young,frictionAngle=frictionAngle,poisson=poisson,density=4800,sigmaT=sigmaT,crackOpening=crackOpening,epsCrackOnset=epsCrackOnset,G_over_E=G_over_E,isoPrestress=isoPrestress))
 sphDict=pickle.load(open(packingFile))
 from yade import pack
 sp=pack.SpherePack()

=== modified file 'examples/concrete/uniax.py'
--- examples/concrete/uniax.py	2010-12-08 15:20:55 +0000
+++ examples/concrete/uniax.py	2011-08-29 18:23:23 +0000
@@ -47,7 +47,7 @@
 	sigmaT=3.5e6,
 	frictionAngle=atan(0.8),
 	epsCrackOnset=1e-4,
-	relDuctility=30,
+	crackOpening=1e-6,
 
 	intRadius=1.5,
 	dtSafety=.8,
@@ -76,7 +76,7 @@
 # make geom; the dimensions are hard-coded here; could be in param table if desired
 # z-oriented hyperboloid, length 20cm, diameter 10cm, skirt 8cm
 # using spheres 7mm of diameter
-concreteId=O.materials.append(CpmMat(young=young,frictionAngle=frictionAngle,poisson=poisson,density=4800,sigmaT=sigmaT,relDuctility=relDuctility,epsCrackOnset=epsCrackOnset,G_over_E=G_over_E,isoPrestress=isoPrestress))
+concreteId=O.materials.append(CpmMat(young=young,frictionAngle=frictionAngle,poisson=poisson,density=4800,sigmaT=sigmaT,crackOpening=crackOpening,epsCrackOnset=epsCrackOnset,G_over_E=G_over_E,isoPrestress=isoPrestress))
 
 spheres=pack.randomDensePack(pack.inHyperboloid((0,0,-.5*specimenLength),(0,0,.5*specimenLength),.25*specimenLength,.17*specimenLength),spheresInCell=2000,radius=sphereRadius,memoizeDb='/tmp/triaxPackCache.sqlite',material=concreteId)
 #spheres=pack.randomDensePack(pack.inAlignedBox((-.25*specimenLength,-.25*specimenLength,-.5*specimenLength),(.25*specimenLength,.25*specimenLength,.5*specimenLength)),spheresInCell=2000,radius=sphereRadius,memoizeDb='/tmp/triaxPackCache.sqlite')

=== modified file 'pkg/dem/ConcretePM.cpp'
--- pkg/dem/ConcretePM.cpp	2011-08-27 17:40:30 +0000
+++ pkg/dem/ConcretePM.cpp	2011-08-29 18:23:23 +0000
@@ -29,13 +29,13 @@
 	if(!mat1->neverDamage) {
 		assert(!isnan(mat1->sigmaT));
 		assert(!isnan(mat1->epsCrackOnset));
-		assert(!isnan(mat1->crackOpening));
+		assert(!isnan(mat1->crackOpening) || !isnan(mat1->relDuctility));
 		assert(!isnan(mat1->G_over_E));
 	}
 	if(!mat2->neverDamage) {
 		assert(!isnan(mat2->sigmaT));
 		assert(!isnan(mat2->epsCrackOnset));
-		assert(!isnan(mat2->crackOpening));
+		assert(!isnan(mat2->crackOpening) || !isnan(mat2->relDuctility));
 		assert(!isnan(mat2->G_over_E));
 	}
 
@@ -50,6 +50,7 @@
 		#define _CPATTR(a) cpmPhys->a=mat1->a
 			_CPATTR(epsCrackOnset);
 			_CPATTR(crackOpening);
+			_CPATTR(relDuctility);
 			_CPATTR(neverDamage);
 			_CPATTR(dmgTau);
 			_CPATTR(dmgRateExp);
@@ -64,9 +65,10 @@
 			cpmPhys->G=.5*(mat1->G_over_E+mat2->G_over_E)*.5*(mat1->young+mat2->young);
 			cpmPhys->tanFrictionAngle=tan(.5*(mat1->frictionAngle+mat2->frictionAngle));
 			cpmPhys->undamagedCohesion=.5*(mat1->sigmaT+mat2->sigmaT);
-			_AVGATTR(crackOpening);
 			cpmPhys->isCohesive=(cohesiveThresholdIter<0 || scene->iter<cohesiveThresholdIter);
 			_AVGATTR(epsCrackOnset);
+			_AVGATTR(crackOpening);
+			_AVGATTR(relDuctility);
 			cpmPhys->neverDamage=(mat1->neverDamage || mat2->neverDamage);
 			_AVGATTR(dmgTau);
 			_AVGATTR(dmgRateExp);
@@ -173,6 +175,7 @@
 		const Real& epsCrackOnset(BC->epsCrackOnset);
 		Real& relResidualStrength(BC->relResidualStrength);
 		const Real& crackOpening(BC->crackOpening);
+		const Real& epsF(BC->epsCrackOnset*BC->relDuctility);
 		const int& damLaw(BC->damLaw);
 		const bool& neverDamage(BC->neverDamage);
 		Real& omega(BC->omega);
@@ -224,7 +227,8 @@
 		epsN+=BC->isoPrestress/E;
 		// very simplified version of the constitutive law
 		kappaD=max(max(0.,epsN),kappaD); // internal variable, max positive strain (non-decreasing)
-		Real epsFracture = crackOpening/contGeom->refLength;
+		//Real epsFracture = crackOpening/contGeom->refLength;
+		Real epsFracture = isnan(epsF)? crackOpening/contGeom->refLength : epsF;
 		omega=isCohesive?funcG(kappaD,epsCrackOnset,epsFracture,neverDamage,damLaw):1.; // damage variable (non-decreasing, as funcG is also non-decreasing)
 		sigmaN=(1-(epsN>0?omega:0))*E*epsN; // damage taken in account in tension only
 		sigmaT=G*epsT; // trial stress
@@ -437,33 +441,32 @@
 		const Body::id_t id1=I->getId1(), id2=I->getId2();
 		Dem3DofGeom* geom=YADE_CAST<Dem3DofGeom*>(I->geom.get());
 
-		Matrix3r stress = Matrix3r::Zero();
+		Matrix3r stressTimesV = Matrix3r::Zero();
 		const Vector3r& n= geom->normal;
 		const Real& Fn = phys->Fn;
 		const Vector3r& Fs = phys->Fs;
-		//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)));
-		stress += geom->refLength*(Fn*n*n.transpose()+.5*(Fs*n.transpose()+n*Fs.transpose()));
+		//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()));
 		
-		bodyStats[id1].stress += stress;
-		bodyStats[id2].stress += stress;
+		bodyStats[id1].stressTimesV += stressTimesV;
+		bodyStats[id2].stressTimesV += stressTimesV;
 		bodyStats[id1].nLinks++; bodyStats[id2].nLinks++;
 		
 		if(!phys->isCohesive) continue;
 		bodyStats[id1].nCohLinks++; bodyStats[id1].dmgSum+=(1-phys->relResidualStrength); bodyStats[id1].epsPlSum+=phys->epsPlSum;
 		bodyStats[id2].nCohLinks++; bodyStats[id2].dmgSum+=(1-phys->relResidualStrength); bodyStats[id2].epsPlSum+=phys->epsPlSum;
 		maxOmega=max(maxOmega,phys->omega);
-
 		avgRelResidual+=phys->relResidualStrength;
 		nAvgRelResidual+=1;
 	}
 	FOREACH(shared_ptr<Body> B, *scene->bodies){
-		//if (!B) continue;
+		if (!B) continue;
 		const Body::id_t& id=B->getId();
 		// add damaged contacts that have already been deleted
 		CpmState* state=dynamic_cast<CpmState*>(B->state.get());
 		if(!state) continue;
-		//state->sigma=bodyStats[id].sigma;
+		state->stressTimesV=bodyStats[id].stressTimesV;
 		//state->tau=bodyStats[id].tau;
 		int cohLinksWhenever=bodyStats[id].nCohLinks+state->numBrokenCohesive;
 		if(cohLinksWhenever>0){

=== modified file 'pkg/dem/ConcretePM.hpp'
--- pkg/dem/ConcretePM.hpp	2011-08-27 17:40:30 +0000
+++ pkg/dem/ConcretePM.hpp	2011-08-29 18:23:23 +0000
@@ -68,7 +68,7 @@
 		((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,stress,Matrix3r::Zero(),,"Stress tensor 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)")),
 		/*ctor*/ createIndex();
 	);
 	REGISTER_CLASS_INDEX(CpmState,State);
@@ -86,8 +86,8 @@
 		((Real,sigmaT,NaN,,"Initial cohesion [Pa]"))
 		((bool,neverDamage,false,,"If true, no damage will occur (for testing only)."))
 		((Real,epsCrackOnset,NaN,,"Limit elastic strain [-]"))
-		((Real,relDuctility,NaN,,"Relative ductility, for damage evolution law peak right-tangent. [-]"))
 		((Real,crackOpening,NaN,,"Crack opening when the crack is fully broken in tension. [m]"))
+		((Real,relDuctility,NaN,,"Deprecated"))
 		((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,6 +127,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 when the crack is fully broken in tension. [m]"))
+			((Real,relDuctility,NaN,,"Deprecated"))
 			((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"))
 			((Real,dmgStrain,0,,"damage strain (at previous or current step)"))
@@ -247,7 +248,7 @@
 #endif
 
 class CpmStateUpdater: public PeriodicEngine {
-	struct BodyStats{ int nCohLinks; int nLinks; Real dmgSum, epsPlSum; Matrix3r stress; BodyStats(): nCohLinks(0), nLinks(0), dmgSum(0.), epsPlSum(0.), stress(Matrix3r::Zero()) {} };
+	struct BodyStats{ int nCohLinks; int nLinks; Real dmgSum, epsPlSum; Matrix3r stressTimesV; BodyStats(): nCohLinks(0), nLinks(0), dmgSum(0.), epsPlSum(0.), stressTimesV(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-27 17:40:30 +0000
+++ pkg/dem/VTKRecorder.cpp	2011-08-29 18:23:23 +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)->stress;
+					const Matrix3r& ss=YADE_PTR_CAST<CpmState>(b->state)->stressTimesV;
 					//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/eudoxos.py'
--- py/eudoxos.py	2010-09-27 17:47:59 +0000
+++ py/eudoxos.py	2011-08-29 18:23:23 +0000
@@ -127,7 +127,7 @@
 	The format is line-oriented as follows::
 
 		E G                                                 # elastic material parameters
-		epsCrackOnset relDuctility xiShear transStrainCoeff # tensile parameters; epsFr=epsCrackOnset*relDuctility
+		epsCrackOnset crackOpening xiShear transStrainCoeff # tensile parameters; epsFr=crackOpening/len
 		cohesionT tanPhi                                    # shear parameters
 		number_of_spheres number_of_links
 		id x y z r boundary                                 # spheres; boundary: -1 negative, 0 none, 1 positive

=== modified file 'py/export.py'
--- py/export.py	2011-08-27 17:40:30 +0000
+++ py/export.py	2011-08-29 18:23:23 +0000
@@ -216,6 +216,7 @@
 		bodies = []
 		for i in ids:
 			b = O.bodies[i]
+			if not b: continue
 			if b.shape.__class__.__name__!="Sphere":
 				if not allIds: print "Warning: body %d is not Sphere"%(i)
 				continue
@@ -268,6 +269,7 @@
 		bodies = []
 		for i in ids:
 			b = O.bodies[i]
+			if not b: continue
 			if b.shape.__class__.__name__!="Facet":
 				if not allIds: print "Warning: body %d is not Facet"%(i)
 				continue
@@ -306,3 +308,4 @@
 				for b in bodies:
 					outFile.write("%g\n"%(eval(command)))
 		outFile.close()
+		self.facetsSnapCount += 1

=== modified file 'py/pack/pack.py'
--- py/pack/pack.py	2011-07-29 02:38:16 +0000
+++ py/pack/pack.py	2011-08-29 18:23:23 +0000
@@ -368,7 +368,7 @@
 
 
 
-def randomDensePack(predicate,radius,material=-1,dim=None,cropLayers=0,rRelFuzz=0.,spheresInCell=0,memoizeDb=None,useOBB=True,memoDbg=False,color=None):
+def randomDensePack(predicate,radius,material=-1,dim=None,cropLayers=0,rRelFuzz=0.,spheresInCell=0,memoizeDb=None,useOBB=True,memoDbg=False,color=None,returnSpherePack=None):
 	"""Generator of random dense packing with given geometry properties, using TriaxialTest (aperiodic)
 	or PeriIsoCompressor (periodic). The periodicity depens on whether	the spheresInCell parameter is given.
 
@@ -392,6 +392,7 @@
 		computed first; it can reduce substantially number of spheres for the triaxial compression (like 10× depending on
 		how much asymmetric the body is), see scripts/test/gts-triax-pack-obb.py.
 	:param memoDbg: show packigns that are considered and reasons why they are rejected/accepted
+	:param returnSpherePack: see :yref:`filterSpherePack`
 
 	:return: SpherePack object with spheres, filtered by the predicate.
 	"""
@@ -426,7 +427,7 @@
 			if orientation:
 				sp.cellSize=(0,0,0) # resetting cellSize avoids warning when rotating
 				sp.rotate(*orientation.toAxisAngle())
-			return filterSpherePack(predicate,sp,material=material)
+			return filterSpherePack(predicate,sp,material=material,returnSpherePack=returnSpherePack)
 		else: print "No suitable packing in database found, running",'PERIODIC compression' if wantPeri else 'triaxial'
 		sys.stdout.flush()
 	O.switchScene(); O.resetThisScene() ### !!
@@ -465,7 +466,7 @@
 	if orientation:
 		sp.cellSize=(0,0,0); # reset periodicity to avoid warning when rotating periodic packing
 		sp.rotate(*orientation.toAxisAngle())
-	return filterSpherePack(predicate,sp,material=material,color=color)
+	return filterSpherePack(predicate,sp,material=material,color=color,returnSpherePack=returnSpherePack)
 
 def randomPeriPack(radius,initSize,rRelFuzz=0.0,memoizeDb=None):
 	"""Generate periodic dense packing.

=== modified file 'scripts/test/cpm-dem3dof-scgeom.py'
--- scripts/test/cpm-dem3dof-scgeom.py	2010-10-29 10:12:44 +0000
+++ scripts/test/cpm-dem3dof-scgeom.py	2011-08-29 18:23:23 +0000
@@ -13,7 +13,7 @@
 dist=r1+r2
 offset=Vector3(0,0,2*dist)
 p1,p2=Vector3(0,0,0),Vector3(dist,0,0)
-O.materials.append(CpmMat(young=30e9,poisson=.2,frictionAngle=atan(.8),sigmaT=3e6,relDuctility=5,epsCrackOnset=1e-4,G_over_E=.2,neverDamage=False,plTau=-1,plRateExp=0,dmgTau=-1,dmgRateExp=0))
+O.materials.append(CpmMat(young=30e9,poisson=.2,frictionAngle=atan(.8),sigmaT=3e6,crackOpening=1e-6,epsCrackOnset=1e-4,G_over_E=.2,neverDamage=False,plTau=-1,plRateExp=0,dmgTau=-1,dmgRateExp=0))
 # first 2 spheres used for Dem3DofGeom
 # the other 2 used for ScGeom (#3 is dynamic, since ScGeom needs that)
 O.bodies.append([utils.sphere(p1,r1,dynamic=False),utils.sphere(p2,r2,dynamic=False)])

=== modified file 'scripts/test/peri3dController_example1.py'
--- scripts/test/peri3dController_example1.py	2011-01-20 14:55:15 +0000
+++ scripts/test/peri3dController_example1.py	2011-08-29 18:23:23 +0000
@@ -4,7 +4,7 @@
 from yade import pack, plot
 
 # create some material
-O.materials.append(CpmMat(young=25e9,frictionAngle=.7,G_over_E=.2,sigmaT=3e6,epsCrackOnset=1e-4,relDuctility=30))
+O.materials.append(CpmMat(young=25e9,frictionAngle=.7,G_over_E=.2,sigmaT=3e6,epsCrackOnset=1e-4,crackOpening=5e-6))
 
 # create periodic assembly of particles
 initSize=1.2

=== modified file 'scripts/test/peri3dController_shear.py'
--- scripts/test/peri3dController_shear.py	2011-01-20 14:55:15 +0000
+++ scripts/test/peri3dController_shear.py	2011-08-29 18:23:23 +0000
@@ -8,7 +8,7 @@
 from yade import pack,plot,qt
 
 # define material
-O.materials.append(CpmMat(young=25e9,G_over_E=.2,sigmaT=3e6,epsCrackOnset=1e-4,relDuctility=30))
+O.materials.append(CpmMat(young=25e9,G_over_E=.2,sigmaT=3e6,epsCrackOnset=1e-4,crackOpening=1e-6))
 
 # create periodic assembly of particles
 initSize=1.2