yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #10074
[Branch ~yade-pkg/yade/git-trunk] Rev 3730: Proposal of new (similar but more compact) examples using JCFpm classes. A choice can be made one...
------------------------------------------------------------
revno: 3730
committer: Jerome Duriez <jerome.duriez@xxxxxxxxxxxxxxx>
timestamp: Fri 2013-10-04 11:47:21 +0200
message:
Proposal of new (similar but more compact) examples using JCFpm classes. A choice can be made one day after some feedback ?
added:
examples/jointedCohesiveFrictionalPM/gravityBis.py
examples/jointedCohesiveFrictionalPM/identifBis.py
modified:
examples/jointedCohesiveFrictionalPM/README
--
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
=== modified file 'examples/jointedCohesiveFrictionalPM/README'
--- examples/jointedCohesiveFrictionalPM/README 2013-08-01 11:54:13 +0000
+++ examples/jointedCohesiveFrictionalPM/README 2013-10-04 09:47:21 +0000
@@ -11,3 +11,6 @@
2-identificationSpheresOnJoint.py is a script used to identify the spheres belonging to a predefined packing (parallellepiped_10.spheres) interacting along pre-existing discontinuity surfaces defined by a meshed surface (persistentPlane30Deg.stl). Executing this script produces 2 text files containing respectively spheres from packing (parallellepiped_10_persistentPlane30Deg.spheres) and spheres attributes regarding jcfPM variables (parallellepiped_10_persistentPlane30Deg_jointedPM.spheres).
3-gravityLoading.py is a simple example showing how to use Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM. It simulates the loading of a pre-fractured packing (parallellepiped_10_persistentPlane30Deg.spheres) by gravity to emphasize how the pre-existing discontinuity affects the behavior. User can play along with joint surface properties (smooth contact logic or not, joint friction angle,...) to see the effect on the simulated behavior.
+
+
+A more compact way of use is proposed in gravityBis.py. This script is stand-alone, and calls itself the scrit identifBis.py. This later script can be called from any script in JCFpm modelling to detect spheres that should interact according to smooth contact logic.
\ No newline at end of file
=== added file 'examples/jointedCohesiveFrictionalPM/gravityBis.py'
--- examples/jointedCohesiveFrictionalPM/gravityBis.py 1970-01-01 00:00:00 +0000
+++ examples/jointedCohesiveFrictionalPM/gravityBis.py 2013-10-04 09:47:21 +0000
@@ -0,0 +1,125 @@
+# encoding: utf-8
+
+# example of use JCFpm classes : Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM...
+# a cubic rock massif, containing a rock joint with ~ 31 deg. dip angle. At one time, the mechanical properties of the joint are degraded, triggering a smooth sliding
+
+
+# definition of a predicate for use of randomDensePack() function
+
+from yade import pack
+dimModele = 10.0
+pred = pack.inAlignedBox((0,0,0),(dimModele,dimModele,dimModele))
+
+# the packing of spheres :
+
+def mat(): return JCFpmMat(type=1,young=1e8,frictionAngle=radians(30),density=3000)
+nSpheres = 3000.0
+poros=0.13
+rMeanSpheres = dimModele * pow(3.0/4.0*(1-poros)/(pi*nSpheres),1.0/3.0)
+print ''
+print 'Creating a cubic sample of spheres (may take some time)'
+print ''
+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)
+
+# Definition of the facets for joint's geometry
+
+import gts
+# joint with ~ 31 deg. dip angle
+v1 = gts.Vertex(0 , 0 , 0.8*dimModele)
+v2 = gts.Vertex(0 , dimModele , 0.8*dimModele )
+v3 = gts.Vertex(dimModele , dimModele , 0.2*dimModele)
+v4 = gts.Vertex(dimModele , 0 , 0.2*dimModele)
+
+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)
+
+yade.qt.View()
+yade.qt.Controller()
+
+O.saveTmp()
+# identification of spheres onJoint, and so on:
+execfile('identifBis.py')
+dim=utils.aabbExtrema()
+xsup=dim[1][0]
+yinf=dim[0][1]
+ysup=dim[1][1]
+zinf=dim[0][2]
+zsup=dim[1][2]
+
+e=2*rMeanSpheres
+
+for o in O.bodies:
+ if isinstance(o.shape,Sphere):
+ o.shape.color=(0.9,0.8,0.6)
+ ## to fix boundary particles on ground
+ if o.state.pos[2]<(zinf+2*e) :
+ o.state.blockedDOFs+='xyz'
+ o.shape.color=(1,1,1)
+
+ ## to identify indicator on top
+ if o.state.pos[2]>(zsup-e) and o.state.pos[0]>(xsup-e) and o.state.pos[1]>((yinf+ysup-e)/2.0) and o.state.pos[1]<((yinf+ysup+e)/2) :
+ refPoint=o.id
+ p0=o.state.pos[2] # its initial position
+
+O.bodies[refPoint].shape.highlight=True
+
+#### Engines definition
+interactionRadius=1. # to set initial contacts to larger neighbours
+O.engines=[
+ ForceResetter(),
+ InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=interactionRadius),]),
+ InteractionLoop(
+ [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=interactionRadius)],
+ [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,alpha=0.3,tensileStrength=1e6,cohesion=1e6,jointNormalStiffness=1e7,jointShearStiffness=1e7,jointCohesion=1e6,jointFrictionAngle=radians(20),jointDilationAngle=0.0)],
+ [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(smoothJoint=True)]
+ ),
+ GlobalStiffnessTimeStepper(timestepSafetyCoefficient=0.8),
+ PyRunner(iterPeriod=1000,initRun=False,command='jointStrengthDegradation()'),
+ PyRunner(iterPeriod=10,initRun=True,command='dataCollector()'),
+ NewtonIntegrator(damping=0.7,gravity=(0.,0.,-10.)),
+
+]
+
+#### dataCollector
+from yade import plot
+plot.plots={'iterations':('v',)}
+def dataCollector():
+ R=O.bodies[refPoint]
+ plot.addData(v=R.state.vel[2],p=R.state.pos[2]-p0,iterations=O.iter,t=O.realtime)
+
+#### joint strength degradation
+stableIter=2000
+stableVel=0.001
+degrade=True
+def jointStrengthDegradation():
+ global degrade
+ if degrade and O.iter>=stableIter and abs(O.bodies[refPoint].state.vel[1])<stableVel :
+ print '!joint cohesion total degradation!', ' | iteration=', O.iter
+ degrade=False
+ for i in O.interactions:
+ if i.phys.isOnJoint :
+ if i.phys.isCohesive:
+ i.phys.isCohesive=False
+ i.phys.FnMax=0.
+ i.phys.FsMax=0.
+
+O.step()
+print ''
+print 'Seeking first after an initial equilibrium state'
+print ''
+O.run(10000)
+plot.plot()
=== added file 'examples/jointedCohesiveFrictionalPM/identifBis.py'
--- examples/jointedCohesiveFrictionalPM/identifBis.py 1970-01-01 00:00:00 +0000
+++ examples/jointedCohesiveFrictionalPM/identifBis.py 2013-10-04 09:47:21 +0000
@@ -0,0 +1,161 @@
+# -*- coding: utf-8 -*-
+
+# ---- Script to detect spheres which are "onJoint", according to JCFpm. -----
+# To be called directly within an other script, for example, with execfile('identifBis.py')
+# The sample (spheres + facets) has to exist already, with their JCFpmMat
+
+
+############################ engines definition
+interactionRadius=1.;
+O.engines=[
+
+ ForceResetter(),
+ InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=interactionRadius,label='is2aabb'),Bo1_Facet_Aabb()]),
+ InteractionLoop(
+ [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=interactionRadius,label='ss2d3dg'),Ig2_Facet_Sphere_ScGeom()],
+ [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,alpha=1,tensileStrength=1e6,cohesion=1e6,jointNormalStiffness=1,jointShearStiffness=1,jointTensileStrength=1e6,jointCohesion=1e6,jointFrictionAngle=1,label='interactionPhys')],
+ [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(smoothJoint=True,label='interactionLaw')]
+ ),
+ NewtonIntegrator(damping=1)
+
+]
+
+############################ timestep + opening yade windows
+O.dt=0.001*utils.PWaveTimeStep()
+
+# from yade import qt
+# v=qt.Controller()
+# v=qt.View()
+
+############################ Identification spheres on joint
+#### color set for particles on joint
+jointcolor1=(0,1,0)
+jointcolor2=(1,0,0)
+jointcolor3=(0,0,1)
+jointcolor4=(1,1,1)
+jointcolor5=(0,0,0)
+
+#### first step-> find spheres on facet
+O.step();
+
+for i in O.interactions:
+ ##if not i.isReal : continue
+ ### Rk: facet are only stored in id1
+ if isinstance(O.bodies[i.id1].shape,Facet) and isinstance(O.bodies[i.id2].shape,Sphere):
+ vertices=O.bodies[i.id1].shape.vertices
+ normalRef=vertices[0].cross(vertices[1]) # defines the normal to the facet normalRef
+ nRef=normalRef/(normalRef.norm()) ## normalizes normalRef
+ normalFacetSphere=i.geom.normal # geom.normal is oriented from id1 to id2 -> normalFacetSphere from facet (i.id1) to sphere (i.id2)
+
+ if O.bodies[i.id2].mat.onJoint==False : ## particles has not yet been identified as belonging to a joint plane
+ O.bodies[i.id2].mat.onJoint=True
+ O.bodies[i.id2].mat.joint=1
+ O.bodies[i.id2].shape.color=jointcolor1
+ if nRef.dot(normalFacetSphere)>=0 :
+ O.bodies[i.id2].mat.jointNormal1=nRef
+ elif nRef.dot(normalFacetSphere)<0 :
+ O.bodies[i.id2].mat.jointNormal1=-nRef
+ elif O.bodies[i.id2].mat.onJoint==True : ## particles has already been identified as belonging to, at least, 1 facet
+ if O.bodies[i.id2].mat.joint==1 and ((O.bodies[i.id2].mat.jointNormal1.cross(nRef)).norm()>0.05) : ## particles has already been identified as belonging to only 1 facet
+ O.bodies[i.id2].mat.joint=2
+ O.bodies[i.id2].shape.color=jointcolor2
+ if nRef.dot(normalFacetSphere)>=0 :
+ O.bodies[i.id2].mat.jointNormal2=nRef
+ elif nRef.dot(normalFacetSphere)<0 :
+ O.bodies[i.id2].mat.jointNormal2=-nRef
+ elif O.bodies[i.id2].mat.joint==2 and ((O.bodies[i.id2].mat.jointNormal1.cross(nRef)).norm()>0.05) and ((O.bodies[i.id2].mat.jointNormal2.cross(nRef)).norm()>0.05): ## particles has already been identified as belonging to more than 1 facet
+ O.bodies[i.id2].mat.joint=3
+ O.bodies[i.id2].shape.color=jointcolor3
+ if nRef.dot(normalFacetSphere)>=0 :
+ O.bodies[i.id2].mat.jointNormal3=nRef
+ elif nRef.dot(normalFacetSphere)<0 :
+ O.bodies[i.id2].mat.jointNormal3=-nRef
+ elif O.bodies[i.id2].mat.joint==3 and ((O.bodies[i.id2].mat.jointNormal1.cross(nRef)).norm()>0.05) and ((O.bodies[i.id2].mat.jointNormal2.cross(nRef)).norm()>0.05) and ((O.bodies[i.id2].mat.jointNormal3.cross(nRef)).norm()>0.05):
+ O.bodies[i.id2].mat.joint=4
+ O.bodies[i.id2].shape.color=jointcolor5
+
+#### second step -> find spheres interacting with spheres on facet (could be executed in the same timestep as step 1?)
+for j in O.interactions:
+ #if not i.isReal : continue
+ ## Rk: facet are only stored in id1
+ if isinstance(O.bodies[j.id1].shape,Facet) and isinstance(O.bodies[j.id2].shape,Sphere):
+ vertices=O.bodies[j.id1].shape.vertices
+ normalRef=vertices[0].cross(vertices[1]) # defines the normal to the facet normalRef
+ nRef=normalRef/(normalRef.norm()) ## normalizes normalRef
+ if ((O.bodies[j.id2].mat.jointNormal1.cross(nRef)).norm()<0.05) :
+ jointNormalRef=O.bodies[j.id2].mat.jointNormal1
+ elif ((O.bodies[j.id2].mat.jointNormal2.cross(nRef)).norm()<0.05) :
+ jointNormalRef=O.bodies[j.id2].mat.jointNormal2
+ elif ((O.bodies[j.id2].mat.jointNormal3.cross(nRef)).norm()<0.05) :
+ jointNormalRef=O.bodies[j.id2].mat.jointNormal3
+ else : continue
+ facetCenter=O.bodies[j.id1].state.pos
+ #### seek for each sphere interacting with the identified sphere j.id2
+ for n in O.interactions.withBody(j.id2) :
+ if isinstance(O.bodies[n.id1].shape,Sphere) and isinstance(O.bodies[n.id2].shape,Sphere):
+ if j.id2==n.id1: # the sphere that was detected on facet (that is, j.id2) is id1 of interaction n
+ sphOnF=n.id1
+ othSph=n.id2
+ elif j.id2==n.id2: # here, this sphere that was detected on facet (that is, j.id2) is id2 of interaction n
+ sphOnF=n.id2
+ othSph=n.id1
+ facetSphereDir=(O.bodies[othSph].state.pos-facetCenter)
+ if O.bodies[othSph].mat.onJoint==True :
+ if O.bodies[othSph].mat.joint==3 and ((O.bodies[othSph].mat.jointNormal1.cross(jointNormalRef)).norm()>0.05) and ((O.bodies[othSph].mat.jointNormal2.cross(jointNormalRef)).norm()>0.05) and ((O.bodies[othSph].mat.jointNormal3.cross(jointNormalRef)).norm()>0.05):
+ O.bodies[othSph].mat.joint=4
+ O.bodies[othSph].shape.color=jointcolor5
+ elif O.bodies[othSph].mat.joint==2 and ((O.bodies[othSph].mat.jointNormal1.cross(jointNormalRef)).norm()>0.05) and ((O.bodies[othSph].mat.jointNormal2.cross(jointNormalRef)).norm()>0.05):
+ O.bodies[othSph].mat.joint=3
+ if facetSphereDir.dot(jointNormalRef)>=0:
+ O.bodies[othSph].mat.jointNormal3=jointNormalRef
+ elif facetSphereDir.dot(jointNormalRef)<0:
+ O.bodies[othSph].mat.jointNormal3=-jointNormalRef
+ elif O.bodies[othSph].mat.joint==1 and ((O.bodies[othSph].mat.jointNormal1.cross(jointNormalRef)).norm()>0.05) :
+ O.bodies[othSph].mat.joint=2
+ if facetSphereDir.dot(jointNormalRef)>=0:
+ O.bodies[othSph].mat.jointNormal2=jointNormalRef
+ elif facetSphereDir.dot(jointNormalRef)<0:
+ O.bodies[othSph].mat.jointNormal2=-jointNormalRef
+ elif O.bodies[othSph].mat.onJoint==False :
+ O.bodies[othSph].mat.onJoint=True
+ O.bodies[othSph].mat.joint=1
+ O.bodies[othSph].shape.color=jointcolor4
+ if facetSphereDir.dot(jointNormalRef)>=0:
+ O.bodies[othSph].mat.jointNormal1=jointNormalRef
+ elif facetSphereDir.dot(jointNormalRef)<0:
+ O.bodies[othSph].mat.jointNormal1=-jointNormalRef
+
+#### for visualization:
+#bj=0
+#vert=(0.,1.,0.)
+#hor=(0.,1.,1.)
+#for o in O.bodies:
+ #if o.mat.onJoint==True : # or o.shape.name=='Facet':
+ ##if o.shape.name=='Facet':
+ ##o.shape.wire=True
+ ##o.state.pos+=(0,50,0)
+ ##bj+=1
+ #if o.mat.jointNormal1.dot(hor)>0 :
+ ##o.state.pos+=(0,50,0)
+ #o.shape.color=jointcolor1
+ #elif o.mat.jointNormal1.dot(hor)<0 :
+ ##o.state.pos+=(0,55,0)
+ #o.shape.color=jointcolor2
+ #if o.mat.type>2 :
+ #bj+=1
+ #o.shape.color=jointcolor5
+ ##print o.mat.jointNormal.dot(vert)
+
+
+##### to delete facets (should be OK to start the simulation after that!
+for b in O.bodies:
+ if isinstance(b.shape,Facet):
+ O.bodies.erase(b.id)
+
+O.resetTime()
+O.interactions.clear()
+print ''
+print 'IdentificationSpheresOnJoint executed ! Spheres onJoint (and so on...) detected, facets deleted'
+print ''
+
+
