← Back to team overview

yade-dev team mailing list archive

[Branch ~yade-dev/yade/trunk] Rev 2718: - Add a script demonstrating use of inclined initial cell and reseting trsf

 

------------------------------------------------------------
revno: 2718
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
branch nick: yade
timestamp: Wed 2011-02-09 10:46:15 +0100
message:
  - Add a script demonstrating use of inclined initial cell and reseting trsf
  - Small doc or formatting fixes here and there
added:
  scripts/test/periodic-triax-settingHsize.py
modified:
  pkg/dem/PeriIsoCompressor.cpp
  pkg/dem/PeriIsoCompressor.hpp
  pkg/dem/SpherePack.cpp


--
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/PeriIsoCompressor.cpp'
--- pkg/dem/PeriIsoCompressor.cpp	2011-01-31 15:29:44 +0000
+++ pkg/dem/PeriIsoCompressor.cpp	2011-02-09 09:46:15 +0000
@@ -29,7 +29,6 @@
 			if(!b || !b->bound) continue;
 			for(int i=0; i<3; i++) maxSpan=max(maxSpan,b->bound->max[i]-b->bound->min[i]);
 		}
-
 	}
 	if(maxDisplPerStep<0) maxDisplPerStep=1e-2*charLen; // this should be tuned somehow…
 	const long& step=scene->iter;
@@ -175,7 +174,6 @@
 				assert( mass>0 );//user set
 				Real dampFactor = 1 - growDamping*Mathr::Sign ( strain_rate * ( goal[axis]-stress[axis] ) );
 				strain_rate+=dampFactor*scene->dt* ( goal[axis]-stress[axis] ) /mass;
-				//if ((scene->iter%5000)==0){cerr << axis<<": stress="<<stress[axis]<<", goal="<<goal[axis]<<", velGrad="<<strain_rate<<endl;}
 				LOG_TRACE ( axis<<": stress="<<stress[axis]<<", goal="<<goal[axis]<<", velGrad="<<strain_rate );}
 
 		} else {    // control strain, see "true strain" definition here http://en.wikipedia.org/wiki/Finite_strain_theory

=== modified file 'pkg/dem/PeriIsoCompressor.hpp'
--- pkg/dem/PeriIsoCompressor.hpp	2011-01-02 10:12:16 +0000
+++ pkg/dem/PeriIsoCompressor.hpp	2011-02-09 09:46:15 +0000
@@ -47,7 +47,7 @@
 		virtual void action();
 		void strainStressStiffUpdate();
 	YADE_CLASS_BASE_DOC_ATTRS(PeriTriaxController,BoundaryController,"Engine for independently controlling stress or strain in periodic simulations.\n\n``strainStress`` contains absolute values for the controlled quantity, and ``stressMask`` determines meaning of those values (0 for strain, 1 for stress): e.g. ``( 1<<0 | 1<<2 ) = 1 | 4 = 5`` means that ``strainStress[0]`` and ``strainStress[2]`` are stress values, and ``strainStress[1]`` is strain. \n\nSee scripts/test/periodic-triax.py for a simple example.",
-		((bool,reversedForces,false,,"For broken constitutive laws, normalForce and shearForce on interactions are in the reverse sense. see `bugreport <https://bugs.launchpad.net/yade/+bug/493102>`_"))
+		((bool,reversedForces,false,,"For some constitutive laws (practicaly all laws based on ScGeom), normalForce and shearForce on interactions are in the reverse sense and this flag must be true (mandatory). see `bugreport <https://bugs.launchpad.net/yade/+bug/493102>`_"))
 		((bool,dynCell,false,,"Imposed stress can be controlled using the packing stiffness or by applying the laws of dynamic (dynCell=true). Don't forget to assign a :yref:`mass<PeriTriaxController.mass>` to the cell."))
 		((Vector3r,goal,Vector3r::Zero(),,"Desired stress or strain values (depending on stressMask), strains defined as ``strain(i)=log(Fii)``.\n\n.. warning:: Strains are relative to the :yref:`O.cell.refSize<Cell.refSize>` (reference cell size), not the current one (e.g. at the moment when the new strain value is set)."))
 		((int,stressMask,((void)"all strains",0),,"mask determining strain/stress (0/1) meaning for goal components"))

=== modified file 'pkg/dem/SpherePack.cpp'
--- pkg/dem/SpherePack.cpp	2011-01-31 11:38:54 +0000
+++ pkg/dem/SpherePack.cpp	2011-02-09 09:46:15 +0000
@@ -90,7 +90,7 @@
 	Vector3r size=mx-mn; Real volume=size.x()*size.y()*size.z();
 	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 num<=0."); porosity=0.5;}
+	if(porosity<=0) {LOG_WARN("porosity must be >0, changing it for you. It will be ineffective if rMean>0."); porosity=0.5;}
 	//If rMean is not defined, then in will be defined in RDIST_NUM
 	if(rMean>0) mode=RDIST_RMEAN;
 	else if(num>0 && psdSizes.size()==0) {

=== added file 'scripts/test/periodic-triax-settingHsize.py'
--- scripts/test/periodic-triax-settingHsize.py	1970-01-01 00:00:00 +0000
+++ scripts/test/periodic-triax-settingHsize.py	2011-02-09 09:46:15 +0000
@@ -0,0 +1,49 @@
+# coding: utf-8
+# 2011 ©Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
+"Demonstrate the compression of a periodic cell with non-trivial initial geometry."
+from yade import *
+from yade import pack,log,qt
+log.setLevel('PeriTriaxController',log.TRACE)
+O.periodic=True
+
+O.cell.hSize=Matrix3(0.1,0.1,0, 0,0.2,0, 0,0,0.1)
+sp=pack.SpherePack()
+num=sp.makeCloud(minCorner=Vector3().Zero, maxCorner=(0.1,0.2,0.1), rMean=-0.01,rRelFuzz=.2, num=500,periodic=True, porosity=0.52,distributeMass=False)
+O.bodies.append([utils.sphere(s[0],s[1]) for s in sp])
+
+
+O.engines=[
+	ForceResetter(),
+	InsertionSortCollider([Bo1_Sphere_Aabb()]),
+	InteractionLoop(
+		[Ig2_Sphere_Sphere_ScGeom()],
+		[Ip2_FrictMat_FrictMat_FrictPhys()],
+		[Law2_ScGeom_FrictPhys_CundallStrack()]
+	),
+	#We compress the packing isotropicaly first
+	PeriTriaxController( dynCell=True,mass=0.2,maxUnbalanced=0.01, relStressTol=0.01,goal=(-1e4,-1e4,-1e4), stressMask=7,globUpdate=5, maxStrainRate=(1.,1.,1.), doneHook='triaxDone()', reversedForces=True ,label='triax'),
+	NewtonIntegrator(damping=.2),
+	#PyRunner(iterPeriod=500,command='print "strain: ",triax.strain," stress: ",triax.stress')
+]
+O.dt=utils.PWaveTimeStep()
+qt.View()
+
+phase=0
+def triaxDone():
+	global phase
+	if phase==0:
+		print 'Here we are: stress',triax.stress,'strain',triax.strain
+		#Here we reset the transformation, the compressed packing corresponds to null strain
+		O.cell.trsf=Matrix3.Identity
+		print 'Now εzz will go from 0 to .4 while σxx and σyy will be kept the same.'
+		triax.stressMask=3
+		triax.goal=(-1e4,-1e4,-0.4)
+
+		phase+=1
+	elif phase==1:
+		print 'Here we are: stress',triax.stress,'strain',triax.strain
+		print 'Done, pausing now.'
+		O.pause()
+		
+	
+