yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 3745
committer: Jerome Duriez <jerome.duriez@xxxxxxxxxxxxxxx>
timestamp: Mon 2013-11-04 16:32:10 +0100
message:
  A script example for JCFpm : two rock parts with one joint in the middle, like in laboratory experiments
added:
  examples/jointedCohesiveFrictionalPM/testingJoint.py


--
lp:yade
https://code.launchpad.net/~yade-pkg/yade/git-trunk

Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription

=== added file 'examples/jointedCohesiveFrictionalPM/testingJoint.py'
--- examples/jointedCohesiveFrictionalPM/testingJoint.py	1970-01-01 00:00:00 +0000
+++ examples/jointedCohesiveFrictionalPM/testingJoint.py	2013-11-04 15:32:10 +0000
@@ -0,0 +1,135 @@
+# encoding: utf-8
+
+
+# Abstract : this script defines a "rock joint" sample : two rectangular blocks separated by an horizontal joint surface. Imposing relative movements of the blocks, that are clumps, allows to test directly the behaviour of the joint, described by JCFpm model.
+# jerome.duriez@xxxxxxxxxxxxxxx
+
+
+# Mechanical properties of rock matrix and rock joint :
+def mat(): return JCFpmMat(type=1,young=15.e9,frictionAngle=radians(35),density=3000,poisson=0.35,tensileStrength=4.5e6,cohesion=45.e6,jointNormalStiffness=5.e7,jointShearStiffness=2.5e7,jointCohesion=0.,jointTensileStrength=0.,jointFrictionAngle=radians(35.),jointDilationAngle=0.0)
+
+
+# --- Creating a sample of spheres
+
+# definition of a predicate 
+from yade import pack
+Lx = 10
+Ly = 10
+Lz = 6
+pred = pack.inAlignedBox((0,0,0),(Lx,Ly,Lz))
+# use of randomDensePack() function
+nSpheres = 1500.0
+poros=0.13 # apparently the value of porosity of samples generated by pack.randomDensePack
+rMeanSpheres = pow(Lx*Ly*Lz*3.0/4.0*(1-poros)/(pi*nSpheres),1.0/3.0)
+print '\nGenerating sphere sample, be patient'
+sp = pack.randomDensePack(pred,radius=rMeanSpheres,rRelFuzz=0.3,memoizeDb='/tmp/gts-triax-packings.sqlite',returnSpherePack=True)
+sp.toSimulation(color=(0.9,0.8,0.6),wire=False,material=mat)
+print 'Sphere sample generated !'
+
+
+# --- The joint surface : half of the height
+import gts
+v1 = gts.Vertex(0 , 0 , Lz/2.0)
+v2 = gts.Vertex(Lx, 0 , Lz/2.0)
+v3 = gts.Vertex(Lx, Ly, Lz/2.0)
+v4 = gts.Vertex(0 , Ly, Lz/2.0)
+
+e1 = gts.Edge(v1,v2)
+e2 = gts.Edge(v2,v4)
+e3 = gts.Edge(v4,v1)
+f1 = gts.Face(e1,e2,e3)
+
+e4 = gts.Edge(v4,v3)
+e5 = gts.Edge(v3,v2)
+f2 = gts.Face(e2,e4,e5)
+
+s1 = gts.Surface()
+s1.add(f1)
+s1.add(f2)
+
+facet = gtsSurface2Facets(s1,wire = False,material=mat)
+O.bodies.append(facet)
+
+
+# --- Identification of spheres onJoint, and so on:
+execfile('identifBis.py')
+
+
+# --- Engines definition
+O.engines=[
+	ForceResetter(),
+	InsertionSortCollider([Bo1_Sphere_Aabb()]),
+	InteractionLoop(
+		[Ig2_Sphere_Sphere_ScGeom()],
+		[Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1)],
+		[Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(smoothJoint=True)]),
+	GlobalStiffnessTimeStepper(timestepSafetyCoefficient=0.8),
+        NewtonIntegrator(damping=0.2),
+        PyRunner(command='afficheIt()',initRun=True,iterPeriod=1000),
+]
+def afficheIt():
+	print 'It', O.iter
+
+O.step()
+
+# --- Clumping the blocks
+upperBlock=[]
+lowerBlock=[]
+for inte in O.interactions:
+    if not inte.phys.isOnJoint:
+       bod1 = O.bodies[inte.id1]
+       bod2 = O.bodies[inte.id2]
+       if bod1.state.pos[2]<Lz/2.0:
+          if not (bod1.id in lowerBlock):
+		  lowerBlock.append(bod1.id)
+		  bod1.shape.color=Vector3(1,0,0)
+	  if bod2.state.pos[2]>Lz/2.0:
+             print '\n **** ERROR !!!! ******* \n\n'
+	  else:
+             if not (bod2.id in lowerBlock):
+		     lowerBlock.append(bod2.id)
+		     bod2.shape.color=Vector3(1,0,0)
+       else:
+          if not (bod1.id in upperBlock):
+		  upperBlock.append(bod1.id)
+		  bod1.shape.color=Vector3(0,0,1)
+	  if bod2.state.pos[2]<Lz/2.0:
+             print '\n **** ERROR !!!! ******* \n\n'
+	  else:
+		  if not (bod2.id in upperBlock):
+			  upperBlock.append(bod2.id)
+			  bod2.shape.color=Vector3(0,0,1)
+print '\n Clumping upper block, be patient'
+idUpperClump=O.bodies.clump(upperBlock)
+print 'Clumped !'
+
+print '\n Clumping lower block, be patient'
+idLowerClump=O.bodies.clump(lowerBlock)
+print 'Clumped !'
+
+upperClump = O.bodies[idUpperClump]
+lowerClump = O.bodies[idLowerClump]
+z0 = upperClump.state.pos[2]
+upperClump.dynamic=False
+lowerClump.dynamic=False
+
+
+# --- Saving data
+O.engines =O.engines+[PyRunner(command='dataCollector()',initRun=True,iterPeriod=500)]
+def dataCollector():
+	plot.addData(iterations=O.iter,un=-(upperClump.state.pos[2]-z0),Fz=O.forces.f(idUpperClump,sync=True)[2],Ft=sqrt(pow(O.forces.f(idUpperClump,sync=True)[1],2.)+pow(O.forces.f(idUpperClump,sync=True)[0],2.)),Fzinf=O.forces.f(idLowerClump,sync=True)[2],uf=unbalancedForce(),ec=kineticEnergy())
+
+from yade import plot
+plot.plots={'un':('Fz')}
+
+# --- Simulation !
+compSpeed = 0.0075 # put here a positive value, for compression
+upperClump.state.vel=Vector3(0,0,-compSpeed)
+nIt = 10000
+print '\nComputation begins now, for', nIt,'iterations'
+O.run(nIt,wait=True)
+print 'Computation just finished ! \n'
+
+yade.qt.View()
+yade.qt.Controller()
+plot.plot()