yade-users team mailing list archive

[Shivpreet sharma <question671372@xxxxxxxxxxxxxxxxxxxxx>]

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

hi, I m new on YADE. I have problems with vtk recorder. I m able to export the particles positions of my deposition and the facets that I have, and I am able to generate vtp and vtu files but I am not able to set the conditions in paraview such that I can get a variation of velocities of sphere during collapsing of the slope. Also while executing this script I am getting a statement like 
"!joint cohesion total degradation!  | iteration= 2000"
So does  it have some reference why all the sphere are not falling down during collapsing of the slope above the persistent plane



from yade import plot, pack,utils,ymport

#### controling parameters
packing='iit_10_plane'
smoothContact=True
jointFrict=radians(20)
jointDil=radians(0)
output='jointDip30_jointFrict20'
maxIter=10000

#### Import of the sphere assembly
def sphereMat(): return JCFpmMat(type=1,young=1e8,frictionAngle=radians(30),density=3000,poisson=0.3,tensileStrength=1e6,cohesion=1e6,jointNormalStiffness=1e7,jointShearStiffness=1e7,jointCohesion=1e6,jointFrictionAngle=jointFrict,jointDilationAngle=jointDil) ## Rq: density needs to be adapted as porosity of real rock is different to granular assembly due to difference in porosity (utils.sumForces(baseBodies,(0,1,0))/(Z*X) should be equal to Gamma*g*h with h=Y, g=9.82 and Gamma=2700 kg/m3
O.bodies.append(ymport.text(packing+'.spheres',scale=1.,shift=Vector3(5,-32.4,-4),material=sphereMat))

## preprocessing to get dimensions of the packing
dim=utils.aabbExtrema()
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

## preprocessing to get spheres 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

#### Identification of the spheres on joint (some DIY here!) -> work to do on import function textExt to directly load material properties from the ascii file
inFile=open(packing+'_jointedPM.spheres','r')
for line in inFile:
    if '#' in line : continue
    id = int(line.split()[0])
    onJ = int(line.split()[1])
    nj = int(line.split()[2])
    j11 = float(line.split()[3])
    j12 = float(line.split()[4])
    j13 = float(line.split()[5])
    j21 = float(line.split()[6])
    j22 = float(line.split()[7])
    j23 = float(line.split()[8])
    j31 = float(line.split()[9])
    j32 = float(line.split()[10])
    j33 = float(line.split()[11])
    O.bodies[id].state.onJoint=onJ
    O.bodies[id].state.joint=nj
    O.bodies[id].state.jointNormal1=(j11,j12,j13)
    O.bodies[id].state.jointNormal2=(j21,j22,j23)
    O.bodies[id].state.jointNormal3=(j31,j32,j33)
inFile.close

#### Boundary conditions
e=2*Rmean
Xmax=0
Ymax=0
baseBodies=[]

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[1]<(yinf+2*e) :
	 o.state.blockedDOFs+='xyz'
	 baseBodies.append(o.id)
	 o.shape.color=(1,1,1)

      ## to identify indicator on top
      if o.state.pos[1]>(yinf-e) and o.state.pos[0]>(18) and o.state.pos[2]>(zinf+(Z-e)/2) and o.state.pos[2]<(zsup-(Z-e)/2) : 
	refPoint=o.id
	p0=o.state.pos[1]

baseBodies=tuple(baseBodies)
O.bodies[refPoint].shape.color=(1,0,0)

#### Engines definition
interactionRadius=1. # to set initial contacts to larger neighbours
O.engines=[

	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=interactionRadius,label='is2aabb'),]),
	InteractionLoop(
		[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=interactionRadius,label='ss2d3dg')],
		[Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,label='interactionPhys')],
		[Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(smoothJoint=smoothContact,label='interactionLaw')]
	),
	GlobalStiffnessTimeStepper(timestepSafetyCoefficient=0.8),
	PyRunner(iterPeriod=1000,initRun=False,command='jointStrengthDegradation()'),
	VTKRecorder(iterPeriod=5000,initRun=True,fileName=(output+'-'),recorders=['spheres','velocity','intr']),
	PyRunner(iterPeriod=10,initRun=True,command='dataCollector()'),
	NewtonIntegrator(damping=0.7,gravity=(0.,-9.82,0.)),

]

#### dataCollector
plot.plots={'iterations':('p',)}
def dataCollector():
	R=O.bodies[refPoint]
	plot.addData(v=R.state.vel[1],p=R.state.pos[1]-p0,iterations=O.iter,t=O.realtime)
	plot.saveDataTxt(output)

#### 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.


		
#### YADE windows
from yade import qt
v=qt.Controller()
v=qt.View()

O.dt=0.1*utils.PWaveTimeStep()

#### set cohesive links with interaction radius>=1
O.step();
#### initializes now the interaction detection factor to strictly 1
ss2d3dg.interactionDetectionFactor=-1.
is2aabb.aabbEnlargeFactor=-1.
#### if you want to avoid contact detection (Lattice like)
#O.engines=O.engines[:1]+O.engines[3:]

#### RUN!!!
O.dt=-0.1*utils.PWaveTimeStep()

O.run(maxIter)
plot.plot()




-- 
You received this question notification because your team yade-users is
an answer contact for Yade.