yade-users team mailing list archive

[Yor1 <question270055@xxxxxxxxxxxxxxxxxxxxx>]

New question #270055 on Yade:
https://answers.launchpad.net/yade/+question/270055

Hi everyone,

I try to model the 3 bending point test.
I put 1 million iteration for the test and when i see the results 
i want to add the iteration and not remake again.
This is my code.

Jabrane.



from yade import ymport, utils, pack, export
from yade import plot
from pylab import *
import math

#---------------- SIMULATIONS DEFINED HERE (assembly, material, boundary conditions)

#### packing (previously constructed)
PACKING='poutre240x30x80_1k'
OUT=PACKING+'_flexion_3_points'

#### Simulation Control
DAMP=0.4 # numerical damping
saveData=100 # data record interval
iterMax=1000000 # maximum number of iteration (to be adjusted)
saveVTK=10000 # Vtk files record interval


#### Material microproperties -> Lac du Bonnet granite (cf. A DEM model for soft and hard rock, Scholtes & Donze, JMPS 2013)
intR=1.65# allows near neighbour interaction and coordination number K=13 (determined with coordinationNumber.py -> to be adjusted for every packing)
DENS=4000 # could be adapted to match material density: dens_DEM=dens_rock*(V_rock/V_particles)=dens_rock*1/(1-poro_DEM) -> poro?
YOUNG=65e9 
FRICT=10
ALPHA=0.4
TENS=8e6 
COH=160e6
#### material definition
sphereMat = JCFpmMat(type=1,density=DENS,young=YOUNG,poisson = ALPHA,frictionAngle=radians(FRICT),tensileStrength=TENS,cohesion=COH)
wallMat = JCFpmMat(type=0,density=DENS,young=YOUNG,frictionAngle=radians(0))

for mat in (sphereMat,wallMat):
   O.materials.append(mat) # then wallMat will be used if material is not specified

#### preprocessing to get dimensions
O.bodies.append(ymport.text(PACKING+'.spheres',scale=1.,shift=Vector3(0,0,0)))

dim=utils.aabbExtrema()
xinf=dim[0][0]
xsup=dim[1][0]
X=xsup-xinf
yinf=dim[0][1]
ysup=dim[1][1]
Y=ysup-yinf
zinf=dim[0][2]
zsup=dim[1][2]
Z=zsup-zinf

r=X/15.
#### preprocessing to get dimensions
R=0
Rmax=0
numSpheres=0.
for o in O.bodies:
 if isinstance(o.shape,Sphere):
   numSpheres+=1
   R+=o.shape.radius
   if o.shape.radius>Rmax:
     Rmax=o.shape.radius
Rmean=R/numSpheres

O.reset() # all previous lines were for getting dimensions of the packing to create walls at the right positions (below) because walls have to be genrated after spheres for FlowEngine


O.bodies.append(geom.facetCylinder(center=(xinf+X/4.,yinf-0.9*r,Z/2.),radius=r,height=Z,orientation=Quaternion((1, 0, 0), 0),segmentsNumber=20,wire=False,material=wallMat)) # bas gauche
O.bodies.append(geom.facetCylinder(center=(xinf+0.75*X,yinf-0.9*r,Z/2.),radius=r,height=Z,orientation=Quaternion((1, 0, 0), 0),segmentsNumber=20,wire=False,material=wallMat)) # bas droite
piston=O.bodies.append(geom.facetCylinder(center=(0.5*X,Y+0.9*r,Z/2.),radius=r,height=Z,dynamic=False,orientation=Quaternion((1, 0, 0), 0),segmentsNumber=20,wire=False,material=wallMat)) # haut

p0=O.bodies[piston[0]].state.pos[1]

### packing
beam=O.bodies.append(ymport.text(PACKING+'.spheres',scale=1.,shift=Vector3(0,0,0),material=sphereMat))

# Definition of the facets for joint's geometry
alpha = 90 #angle of the crack with horizontal (pi/4.)
c = 0.375*Y # length of the crack

# 4 points pour l'entaille
ptA =  gts.Vertex( X/2., c, 8./7.*Z)
ptB =  gts.Vertex( X/2., 0, 8./7.*Z)
ptApr = gts.Vertex(X/2., c, -1./7.*Z)
ptBpr = gts.Vertex(X/2., 0, -1./7.*Z)

e1 = gts.Edge(ptA,ptB)
e2 = gts.Edge(ptA,ptApr)
e3 = gts.Edge(ptApr,ptB)
f1 = gts.Face(e1,e2,e3)

e4 = gts.Edge(ptB,ptBpr)
e5 = gts.Edge(ptBpr,ptApr)
f2 = gts.Face(e4,e5,e3)

s1 = gts.Surface()
s1.add(f1)
s1.add(f2)
facet = gtsSurface2Facets(s1,wire = False,material=wallMat)
O.bodies.append(facet)


### set a color to the spheres
for o in O.bodies:
 if isinstance(o.shape,Sphere):
   o.shape.color=(0.7,0.5,0.3)


#---------------- ENGINES DEFINED HERE

#### simulation is defined here (DEM loop, interaction law, servo control, recording, etc...)
##### simulation piston's motion

for i in range(0,len(piston)):
	O.bodies[piston[i]].state.vel[1]=-0.0001

##### simulation beam's grains motion
O.engines=[
        ForceResetter(),
        InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb(aabbEnlargeFactor=intR,label='Saabb'),Bo1_Facet_Aabb()]),
	InteractionLoop(
		[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intR,label='SSgeom'),Ig2_Facet_Sphere_ScGeom()],
		[Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,label='interactionPhys')],
		[Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=True,Key=OUT,label='interactionLaw')]
	),
	
        GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.4, defaultDt=0.1*utils.PWaveTimeStep()),
        NewtonIntegrator(damping=DAMP,label="newton"),
        PyRunner(iterPeriod=int(saveData),initRun=True,command='recorder()',label='data'),
        VTKRecorder(iterPeriod=int(saveVTK),initRun=True,fileName=OUT+'-',recorders=['spheres','bstresses','jcfpm','cracks'],Key=OUT,label='vtk')
]


#### custom recording functions
tensCks=shearCks=cks=cks0=0
def recorder():
    global tensCks, shearCks
    tensCks=0
    shearCks=0
    for o in O.bodies:
	if isinstance(o.shape,Sphere):
		tensCks+=o.state.tensBreak
		shearCks+=o.state.shearBreak
    yade.plot.addData({ 't':O.time
		      ,'i':O.iter
		      ,'f':utils.sumFacetNormalForces(ids=piston,axis=1)
		      ,'eps':p0-O.bodies[piston[0]].state.pos[1]
		      ,'tc':0.5*tensCks,'sc':0.5*shearCks,'unbF':utils.unbalancedForce()}
    )
    yade.plot.saveDataTxt(OUT)

# if you want to plot during simulation
plot.plots={'i':('f')}
#plot.plot()

#---------------- SIMULATION STARTS HERE

#### manage interaction detection factor during the first timestep and then set default interaction range ((cf. A DEM model for soft and hard rock, Scholtes & Donze, JMPS 2013))
O.step();
### initializes the interaction detection factor
SSgeom.interactionDetectionFactor=-1.
Saabb.aabbEnlargeFactor=-1.

#### coordination number verification and reinforcement of boundary particles
numSSlinks=0
numCohesivelinks=0
numFrictionalLinks=0
for i in O.interactions:
    if not i.isReal : continue
    if isinstance(O.bodies[i.id1].shape,Sphere) and isinstance(O.bodies[i.id2].shape,Sphere):
     numSSlinks+=1
     if i.phys.isCohesive :
      numCohesivelinks+=1
     else :
      numFrictionalLinks+=1

print "nbSpheres=", numSpheres," | coordination number =", 2.0*numCohesivelinks/numSpheres


O.run(iterMax)




-- 
You received this question notification because you are a member of
yade-users, which is an answer contact for Yade.