VTK

Asked by shiv

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()

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Robert Caulk
Solved:
Last query:
Last reply:

This question was reopened

  • by shiv
Revision history for this message
Robert Caulk (rcaulk) said :
#1

Hello,

I guess this is mostly a paraview question:

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

As I read it, you would like to plot the particle velocities during the slope collapse? Velocities are easily plotted in paraview, so I get the feeling this might be a VTK iterperiod problem, but tough to speculate any further without more information from you. Can you describe exactly what you are seeing in paraview in addition to exactly what you expect to see?

Regarding the statement: "!joint cohesion total degradation! | iteration= 2000"

This is simply a statement that you are printing in your script when the iteration number exceeds a predefined number of iterations and a reference point Y abs(velocity) exceeds (in a Y direction sense) some predefined Y velocity:

if degrade and O.iter>=stableIter and abs(O.bodies[refPoint].state.vel[1])<stableVel :
    print '!joint cohesion total degradation!', ' | iteration=', O.iter

which appears to be followed by an explicit command to remove all joint interaction cohesion and strength:

 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.

So in other words, as soon as your predefined iteration threshold and velocity value at your reference point are met, your script automatically removes joint cohesion and strength.

Cheers,

Robert

Revision history for this message
shiv (shivpreet.ce15) said :
#2

I have used the above script which is similar to given link "https://github.com/yade/trunk/blob/master/examples/jointedCohesiveFrictionalPM/gravityLoading.py"
After the execution of the script, I am getting two files as output whose names are jointDip30_jointFrict20 with format VTP and VTU.
I am not getting an idea about how after loading this files into paraview how can I get a variation of the velocity of spheres while slope is collapsing similar to this image:
https://drive.google.com/file/d/1WXax8_8J9iC_7846oIzRdTiAZ_69AVyi/view?usp=sharing

I have read also the instruction given in YADE Documentation.

Revision history for this message
Robert Caulk (rcaulk) said :
#3

>>
I have read also the instruction given in YADE Documentation.
>>

VTKRecorder saves the simulation at the frequency that you tell it to. This is highlighted in the YADE documentation [1]. As described in the doc, set:

iterPeriod = some number less than 5000.

Cheers,

Robert

[1]https://yade-dem.org/doc/user.html#saving-data-during-the-simulation

Revision history for this message
Robert Caulk (rcaulk) said :
#4

Actually the direct quote from the doc:

"If the period is too high (and data are saved only few times), the video will have few frames."

:-)

Revision history for this message
shiv (shivpreet.ce15) said :
#5

Ok thanks for the help, Robert. But the thing is I am not getting the procedure to generate the scale of variation of the velocity of spheres with different colour variations similar to the attached link:
https://drive.google.com/file/d/1WXax8_8J9iC_7846oIzRdTiAZ_69AVyi/view?usp=sharing

You might see there on the right Y-axis there is a scale which is showing is the variation of velocity and this is the thing that I am not getting in Paraview

Revision history for this message
Robert Caulk (rcaulk) said :
#6
Revision history for this message
shiv (shivpreet.ce15) said :
#7

Thanks Robert Caulk, that solved my question.

Revision history for this message
shiv (shivpreet.ce15) said :
#8

One more thing, can you suggest me some proper procedure to load VTU files in paraview such I can get the proper visualisation of the sphere in paraview. I have followed the instruction given in YADE documentation but that isn't generating a define diagram such as in the given link:

https://drive.google.com/file/d/1FYMnoAD1dpukmm5RS-jUgd2W0beU14Tr/view?usp=sharing

It is generating like this:
https://drive.google.com/file/d/1LF_roHT65-6k2EAnYCqOkm4QgFAmN9hZ/view?usp=sharing

Revision history for this message
Best Robert Caulk (rcaulk) said :
#9

For the sake of clear documentation, please limit each thread to one question [1].

"5. Do not merge different questions in one single post, even if - in your mind - they are connected to one single project. If your question received satisfying answer(s), mark it answered immediately. Do not use the same thread to jump to another question, open a new question instead."

[1]https://yade-dem.org/wiki/Howtoask

Revision history for this message
shiv (shivpreet.ce15) said :
#10

Thanks Robert Caulk, that solved my question.