yade-dev team mailing list archive

[eudoxos at BerliOS <eudoxos@xxxxxxxxxxxxxxx>]

Author: eudoxos
Date: 2009-04-15 10:19:44 +0200 (Wed, 15 Apr 2009)
New Revision: 1755

Modified:
   trunk/extra/Brefcom.cpp
   trunk/extra/Brefcom.hpp
   trunk/gui/py/yade-multi
   trunk/scripts/multi.py
Log:
1. Add approximate viscosity equations to brefcom (not working)
2. Fixes in yade-multi 


Modified: trunk/extra/Brefcom.cpp
===================================================================
--- trunk/extra/Brefcom.cpp	2009-04-10 11:36:17 UTC (rev 1754)
+++ trunk/extra/Brefcom.cpp	2009-04-15 08:19:44 UTC (rev 1755)
@@ -105,6 +105,7 @@
 		contPhys->dmgRateExp=dmgRateExp;
 		contPhys->plTau=plTau;
 		contPhys->plRateExp=plRateExp;
+		contPhys->viscApprox=viscApprox;
 
 		interaction->interactionPhysics=contPhys;
 	}
@@ -156,6 +157,14 @@
 }
 
 Real BrefcomContact::computeDmgOverstress(Real dt){
+	if(viscApprox){
+		Real prevDmgStrain=dmgStrain;
+		dmgStrain=max(0.,epsN*omega-dmgOverstress/E);
+		if(dmgStrain<=prevDmgStrain) dmgOverstress=0; // damage doesn't grow, no viscous response
+		else dmgOverstress=epsCrackOnset*E*pow(dmgTau*(dmgStrain-prevDmgStrain)/dt,dmgRateExp);
+		LOG_TRACE("dmgStrain="<<dmgStrain<<", dmgOverstress="<<dmgOverstress);
+		return dmgOverstress;
+	}
 	if(dmgStrain>=epsN*omega){ // unloading, no viscous stress
 		dmgStrain=epsN*omega;
 		LOG_TRACE("Elastic/unloading, no viscous overstress");
@@ -174,7 +183,7 @@
 	if(sigmaTNorm<sigmaTYield) return 1.;
 	Real c=undamagedCohesion*pow(plTau/(G*dt),plRateExp)*pow(sigmaTNorm-sigmaTYield,plRateExp-1.);
 	Real beta=solveBeta(c,plRateExp);
-	LOG_DEBUG("scaling factor "<<1.-exp(beta)*(1-sigmaTYield/sigmaTNorm));
+	//LOG_DEBUG("scaling factor "<<1.-exp(beta)*(1-sigmaTYield/sigmaTNorm));
 	return 1.-exp(beta)*(1-sigmaTYield/sigmaTNorm);
 }
 

Modified: trunk/extra/Brefcom.hpp
===================================================================
--- trunk/extra/Brefcom.hpp	2009-04-10 11:36:17 UTC (rev 1754)
+++ trunk/extra/Brefcom.hpp	2009-04-15 08:19:44 UTC (rev 1755)
@@ -72,14 +72,16 @@
 			dmgTau,
 			//! exponent in the rate-dependent damage evolution
 			dmgRateExp,
-			//! damage strain
+			//! damage strain (at previous or current step)
 			dmgStrain,
+			//! damage viscous overstress (at previous step or at current step)
+			dmgOverstress,
 			//! characteristic time for viscoplasticity (if non-positive, no rate-dependence for shear)
 			plTau,
 			//! exponent in the rate-dependent viscoplasticity
-			plRateExp,
-			//! coefficient that takes transversal strain into accound when calculating kappaDReduced
-			transStrainCoeff;
+			plRateExp;
+		//! whether to approximate viscosity with difference relation instead of iterating to find solution
+		bool viscApprox;
 		/*! Up to now maximum normal strain (semi-norm), non-decreasing in time. */
 		Real kappaD;
 		/*! Transversal strain (perpendicular to the contact axis) */
@@ -104,7 +106,7 @@
 
 
 
-		BrefcomContact(): NormalShearInteraction(),E(0), G(0), tanFrictionAngle(0), undamagedCohesion(0), crossSection(0), xiShear(0), dmgTau(-1), dmgRateExp(0), dmgStrain(0), plTau(-1), plRateExp(0), epsPlSum(0.) { createIndex(); epsT=Vector3r::ZERO; kappaD=0; isCohesive=false; neverDamage=false; omega=0; Fn=0; Fs=Vector3r::ZERO; epsPlSum=0; }
+		BrefcomContact(): NormalShearInteraction(),E(0), G(0), tanFrictionAngle(0), undamagedCohesion(0), crossSection(0), xiShear(0), dmgTau(-1), dmgRateExp(0), dmgStrain(0), plTau(-1), plRateExp(0), epsPlSum(0.) { createIndex(); epsT=Vector3r::ZERO; kappaD=0; isCohesive=false; neverDamage=false; omega=0; Fn=0; Fs=Vector3r::ZERO; epsPlSum=0; viscApprox=false; dmgOverstress=0; dmgStrain=0; }
 		//	BrefcomContact(Real _E, Real _G, Real _tanFrictionAngle, Real _undamagedCohesion, Real _equilibriumDist, Real _crossSection, Real _epsCrackOnset, Real _epsFracture, Real _expBending, Real _xiShear, Real _tau=0, Real _expDmgRate=1): InteractionPhysics(), E(_E), G(_G), tanFrictionAngle(_tanFrictionAngle), undamagedCohesion(_undamagedCohesion), equilibriumDist(_equilibriumDist), crossSection(_crossSection), epsCrackOnset(_epsCrackOnset), epsFracture(_epsFracture), expBending(_expBending), xiShear(_xiShear), tau(_tau), expDmgRate(_expDmgRate) { epsT=Vector3r::ZERO; kappaD=0; isCohesive=false; neverDamage=false; omega=0; Fn=0; Fs=Vector3r::ZERO; /*TRVAR5(epsCrackOnset,epsFracture,Kn,crossSection,equilibriumDist); */ }
 		virtual ~BrefcomContact();
 
@@ -121,9 +123,10 @@
 			(dmgTau)
 			(dmgRateExp)
 			(dmgStrain)
+			(dmgOverstress)
 			(plTau)
 			(plRateExp)
-			(transStrainCoeff)
+			(viscApprox)
 
 			(cummBetaIter)
 			(cummBetaCount)
@@ -225,6 +228,7 @@
 		
 		*/
 		Real sigmaT, epsCrackOnset, relDuctility, G_over_E, tau, expDmgRate, omegaThreshold, dmgTau, dmgRateExp, plTau, plRateExp;
+		bool viscApprox;
 		//! Should new contacts be cohesive? They will before this iter#, they will not be afterwards. If 0, they will never be. If negative, they will always be created as cohesive.
 		long cohesiveThresholdIter;
 		//! Create contacts that don't receive any damage (BrefcomContact::neverDamage=true); defaults to false
@@ -233,7 +237,7 @@
 		BrefcomMakeContact(){
 			// init to signaling_NaN to force crash if not initialized (better than unknowingly using garbage values)
 			sigmaT=epsCrackOnset=relDuctility=G_over_E=std::numeric_limits<Real>::signaling_NaN();
-			neverDamage=false;
+			neverDamage=viscApprox=false;
 			cohesiveThresholdIter=-1;
 			dmgTau=-1; dmgRateExp=0; plTau=-1; plRateExp=-1;
 			omegaThreshold=0.999;
@@ -251,6 +255,7 @@
 			(dmgRateExp)
 			(plTau)
 			(plRateExp)
+			(viscApprox)
 			(omegaThreshold)
 		);
 

Modified: trunk/gui/py/yade-multi
===================================================================
--- trunk/gui/py/yade-multi	2009-04-10 11:36:17 UTC (rev 1754)
+++ trunk/gui/py/yade-multi	2009-04-15 08:19:44 UTC (rev 1755)
@@ -82,7 +82,7 @@
 
 print "Will run `%s' on `%s' with nice value %d, output redirected to `%s', %d jobs at a time."%(executable,simul,nice,logFormat,maxJobs)
 
-ll=['']+open(table,'r').readlines()
+ll=[re.sub('\s*#.*','',l) for l in ['']+open(table,'r').readlines()] # remove comments
 availableLines=[i for i in range(len(ll)) if not re.match(r'^\s*(#.*)?$',ll[i][:-1]) and i>1]
 
 # read actual data

Modified: trunk/scripts/multi.py
===================================================================
--- trunk/scripts/multi.py	2009-04-10 11:36:17 UTC (rev 1754)
+++ trunk/scripts/multi.py	2009-04-15 08:19:44 UTC (rev 1755)
@@ -2,8 +2,8 @@
 # see http://yade.wikia.com/wiki/ScriptParametricStudy#scripts.2Fmulti.py for explanations
 #
 ## read parameters from table here
-from yade import utils
-utils.readParamsFromTable(gravity=-9.81,density=2400,initialSpeed=20,noTableOk=True)
+from yade import utils, plot
+utils.readParamsFromTable(gravity=-9.81,density=2400,initialSpeed=10,noTableOk=True)
 print gravity,density,initialSpeed
 
 o=Omega()
@@ -22,8 +22,9 @@
 		[InteractingSphere2InteractingSphere4SpheresContactGeometry(),InteractingBox2InteractingSphere4SpheresContactGeometry()],
 		[SimpleElasticRelationships(),],
 		[ef2_Spheres_Elastic_ElasticLaw(),]
-	)
+	),
 	GravityEngine(gravity=(0,0,gravity)), ## here we use the 'gravity' parameter
+	PeriodicPythonRunner(iterPeriod=100,command='myAddPlotData()',label='plotDataCollector'),
 	NewtonsDampedLaw(damping=0.4)
 ]
 o.bodies.append(utils.box([0,50,0],extents=[1,50,1],dynamic=False,color=[1,0,0],young=30e9,poisson=.3,density=density)) ## here we use the density parameter
@@ -33,11 +34,16 @@
 
 o.dt=.8*utils.PWaveTimeStep()
 ## o.saveTmp('initial')
+def myAddPlotData():
+	s=O.bodies[1]
+	plot.addData({'t':O.time,'y_sph':s.phys.pos[1],'z_sph':s.phys.pos[2]})
+plot.plots={'y_sph':('z_sph',)}
 
 # run 30000 iterations
-o.run(30000,True)
+o.run(20000,True)
 
 # write some results to a common file
 # (we rely on the fact that 2 processes will not write results at exactly the same time)
 # 'a' means to open for appending
 file('multi.out','a').write('%s %g %g %g %g\n'%(o.tags['description'],gravity,density,initialSpeed,o.bodies[1].phys.pos[1]))
+print 'gnuplot',plot.saveGnuplot(O.tags['id'])