------------------------------------------------------------ revno: 3745 committer: Jerome Duriez <[email protected]> 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. +# [email protected] + + +# 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()
_______________________________________________ Mailing list: https://launchpad.net/~yade-dev Post to : [email protected] Unsubscribe : https://launchpad.net/~yade-dev More help : https://help.launchpad.net/ListHelp
