How to pause

Asked by 内山康太郎

Thanks to Jean, I was able to come this far. First of all, thank you Jean.

I am doing an analysis on a direct(simple) shear test.
I have created a program like below.

Next, I want idsCyl2 to move not only in the X direction, but also in the Y direction. In the future, I would like to input seismic waveforms in the XY direction as well.

Now, after moving idsCyl2 (O.bodies[idsCyl2[0]].state.displ()[0]) by 0.020 in the X direction, I want to move idsCyl2 by 0.020 in the Y direction. We also want to repeat in the negative direction.
For example, the locus for dsCyl2 is (0.020, 0) ⇒ (0.020, 0.020) ⇒ (-0.020, 0.020) ⇒ (-0.020, -0.020) ⇒ (0.020, -0.020).

At that time, to run idsCyl2= on Yade's terminal, redefine the engine and enter it. So is there a way to pause the analysis when (0.020 , 0) is reached?

Of course I know how to do it with O.pause() and an input file
def checkDistortion():
          If O.bodies[idsCyl2[0]].state.displ()[0] > .02:
                          O.Pause()
I tried to input, but the analysis itself stopped and I could not move to (0.020, 0.020).

Also, in Yade's terminal
          If O.bodies[idsCyl2[0]].state.displ()[0] > .02:
                          O.Pause()
But it just gets cylDisplacementx at that point and doesn't pause when it reaches (0.020, 0).

Is there a way to move idsCyl2 in the X direction (0.020,0) and pause it?

####################################################################
from __future__ import print_function
from yade import pack,plot,polyhedra_utils,geom
from yade import export,qt
from yade.gridpfacet import *
import math
import numpy as np
import os

frictangle = np.radians(25)#45do
density = 3400.0
young = 3e8
gravity = (0.0, 0.0, 0)
velocity_cyl=1
velocity_cyl_idsCyl2=1
velocity_cyl_Box=1

# create cloud of spheres and insert them into the simulation
# we give corners, mean radius, radius variation

###############################################
#Create sphere
###############################################

mat_sp = CohFrictMat(alphaKr=.2,alphaKtw=.2,young = young, poisson = 0.35,frictionAngle=(frictangle), density=density)
O.materials.append(mat_sp)

pred=pack.inCylinder(centerBottom=(0,0,-0.01),centerTop=(0,0,.02),radius=.06)
sp0=pack.randomDensePack(pred,spheresInCell=3000,radius=.003)
O.bodies.append(sp0)

###############################################
#Boundary creation
###############################################

idsCyl1 = O.bodies.append(geom.facetCylinder((0, 0, 0.00), radius=.06, height=.02, segmentsNumber=30, wallMask=6))

idsCyl2 = O.bodies.append(geom.facetCylinder((0, 0, 0.015), radius=.06, height=.01, segmentsNumber=30, wallMask=5))

#idsCyl3 = O.bodies.append(geom.facetCylinder((0, 0, .0225), radius=0.10, height=.005, segmentsNumber=100, wallMask=2))

radius1 = .06
radius2 = .15

fixBoxIds=[]
moveBoxIds=[]

i=0
for i in range(0,365,5):
 r = math.radians(i)
 x1 = radius1 * math.cos(r)
 y1 = radius1 * math.sin(r)
 x2 = radius2 * math.cos(r-math.radians(5))
 y2 = radius2 * math.sin(r-math.radians(5))
 x3 = radius2 * math.cos(r+math.radians(5))
 y3 = radius2 * math.sin(r+math.radians(5))

 x4 = radius1 * math.cos(r-math.radians(5))
 y4 = radius1 * math.sin(r-math.radians(5))
 x5 = radius1 * math.cos(r+math.radians(5))
 y5 = radius1 * math.sin(r+math.radians(5))
 x6 = radius2 * math.cos(r)
 y6 = radius2 * math.sin(r)

 fixBoxIds.append(O.bodies.append(facet([(x1,y1,.01),(x2,y2,.01),(x3,y3,.01)])))
 fixBoxIds.append(O.bodies.append(facet([(x4,y4,.01),(x5,y5,.01),(x6,y6,.01)])))
 moveBoxIds.append(O.bodies.append(facet([(x1,y1,.01),(x2,y2,.01),(x3,y3,.01)])))
 moveBoxIds.append(O.bodies.append(facet([(x4,y4,.01),(x5,y5,.01),(x6,y6,.01)])))

upper_plate=box((0,0,-0.014),(0.06,0.06,0.001))
upper_plateID=O.bodies.append(upper_plate)
O.bodies[upper_plateID].state.blockedDOFs = 'xyXY'

#O.forces.addF(upper_plateID,(0,0,50e3),permanent=True)
O.forces.setPermF(upper_plateID,(0,0,50e3))
###############################################
#engines
###############################################

O.engines = [
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Box_Aabb()]),
        InteractionLoop(
                # interaction loop
                [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
                [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
                [Law2_ScGeom6D_CohFrictPhys_CohesionMoment(always_use_moment_law=True),Law2_ScGeom_FrictPhys_CundallStrack()]
        ),
        NewtonIntegrator(gravity=gravity,damping=0.2),
        ##TranslationEngine(translationAxis=[0,0,1],velocity=-1,ids=idsCyl3),
        TranslationEngine(translationAxis=[1,0,0],velocity=velocity_cyl_idsCyl2,ids=idsCyl2),
        TranslationEngine(translationAxis=[1,0,0],velocity=velocity_cyl_Box,ids=moveBoxIds),
        PyRunner(command='checkStress1()', iterPeriod=10),
 #PyRunner(command='checkStress2()', iterPeriod=10),
        PyRunner(command='plotPlotData()',iterPeriod=1),
]

def checkStress1():
 cylDisplacement_x = O.bodies[idsCyl2[0]].state.displ()[0]
 cylDisplacement_y = O.bodies[idsCyl2[0]].state.displ()[1]
 plot.addData(
 cylDisplacement_x=cylDisplacement_x,
 cylDisplacement_y=cylDisplacement_y,
 )
 print('{:1.04f}'.format(cylDisplacement_x)," ", '{:1.04f}'.format(cylDisplacement_y))

###############################################
#output
###############################################
sp_num = len([b for b in O.bodies if isinstance(b.shape, Sphere)]) # 粒子の総数をカウント
print("number of spheres = ",sp_num)

O.dt=0.5*PWaveTimeStep()

def plotPlotData():
 pipe_force_x = sum([O.forces.f(i)[0] for i in idsCyl2])
 pipe_force_y = sum([O.forces.f(i)[1] for i in idsCyl2])
 cylDisplacementx = O.bodies[idsCyl2[0]].state.displ()[0]
 cylDisplacementy = O.bodies[idsCyl2[0]].state.displ()[1]
 cylDisplacementz = O.bodies[idsCyl2[0]].state.displ()[2]
 Dz=upper_plate.state.displ()[2]
 plot.addData(
 i=O.iter,
 disp=O.time*velocity_cyl,
 force_x=pipe_force_x,
 force_y=pipe_force_y,
 Dz=Dz,
 cylDisplacementx=cylDisplacementx,
 cylDisplacementy=cylDisplacementy,
 cylDisplacementz=cylDisplacementz,
 )
 plot.saveDataTxt('pipe_force.out',vars=('disp','force_x','force_y',
      'Dz',
      'cylDisplacementx','cylDisplacementx','cylDisplacementx'
      ))

"""
def checkStress2():
 if O.bodies[idsCyl2[0]].state.displ()[0] > 0.0020:
  O.pause()
"""

plot.plots={
 'disp':('force_x','force_y'),
 'i':('Dz'),
 'cylDisplacementx':('cylDisplacementy',)
 }
plot.plot()

#plot.saveDataTxt(O.tags['id'] + '.txt')

####################################################################

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
内山康太郎
Solved:
Last query:
Last reply:
Revision history for this message
内山康太郎 (kenkoutaro) said :
#1

I was able to resolve the issue.

Revision history for this message
Jan Stránský (honzik) said :
#2

Hello,

> I was able to resolve the issue.

please mark the question as solved then.

Optionally, please provide the solution (for others coming here searching same / similar problem)

Cheers
Jan

Revision history for this message
内山康太郎 (kenkoutaro) said :
#3

Run program 1 below.
Once executed, re-enter the engine in the YADE terminal. The program to re-enter is program 2.
At that time, enter the number you defined in Program 1 in NO of PyRunner(command='checkStress_NO()).
By changing the NO defined in Program 1, the PyRunner will stop O.Pause() when it reaches the desired displacement.

If you have any questions, please ask.

program 1↓
############################################################
#program 1
from __future__ import print_function
from yade import pack,plot,polyhedra_utils,geom
from yade import export,qt
from yade.gridpfacet import *
import math
import numpy as np
import os

frictangle = np.radians(25)#45do
density = 3400.0
young = 3e8
gravity = (0.0, 0.0, 0)
velocity_cyl=5
velocity_cyl_idsCyl2=5
velocity_cyl_Box=5

# create cloud of spheres and insert them into the simulation
# we give corners, mean radius, radius variation

###############################################
#Create sphere
###############################################

mat_sp = CohFrictMat(alphaKr=.2,alphaKtw=.2,young = young, poisson = 0.35,frictionAngle=(frictangle), density=density)
O.materials.append(mat_sp)

pred=pack.inCylinder(centerBottom=(0,0,0.0),centerTop=(0,0,.02),radius=.03)
sp0=pack.randomDensePack(pred,spheresInCell=3000,radius=.0004)
O.bodies.append(sp0)

###############################################
#Boundary creation
###############################################

idsCyl1 = O.bodies.append(geom.facetCylinder((0, 0, 0.005), radius=.03, height=.01, segmentsNumber=30, wallMask=6))

idsCyl2 = O.bodies.append(geom.facetCylinder((0, 0, .015), radius=.03, height=.01, segmentsNumber=30, wallMask=5))

#idsCyl3 = O.bodies.append(geom.facetCylinder((0, 0, .0225), radius=0.10, height=.005, segmentsNumber=100, wallMask=2))

radius1 = .03
radius2 = .15

fixBoxIds=[]
moveBoxIds=[]

i=0

for i in range(0,365,5):
 r = math.radians(i)
 x1 = radius1 * math.cos(r)
 y1 = radius1 * math.sin(r)
 x2 = radius2 * math.cos(r-math.radians(5))
 y2 = radius2 * math.sin(r-math.radians(5))
 x3 = radius2 * math.cos(r+math.radians(5))
 y3 = radius2 * math.sin(r+math.radians(5))

 x4 = radius1 * math.cos(r-math.radians(5))
 y4 = radius1 * math.sin(r-math.radians(5))
 x5 = radius1 * math.cos(r+math.radians(5))
 y5 = radius1 * math.sin(r+math.radians(5))
 x6 = radius2 * math.cos(r)
 y6 = radius2 * math.sin(r)

 fixBoxIds.append(O.bodies.append(facet([(x1,y1,.01),(x2,y2,.01),(x3,y3,.01)])))
 fixBoxIds.append(O.bodies.append(facet([(x4,y4,.01),(x5,y5,.01),(x6,y6,.01)])))
 moveBoxIds.append(O.bodies.append(facet([(x1,y1,.01),(x2,y2,.01),(x3,y3,.01)])))
 moveBoxIds.append(O.bodies.append(facet([(x4,y4,.01),(x5,y5,.01),(x6,y6,.01)])))

upper_plate=box((0,0,-0.001),(0.06,0.06,0.001))
upper_plateID=O.bodies.append(upper_plate)
O.bodies[upper_plateID].state.blockedDOFs = 'xyXY'

#O.forces.addF(upper_plateID,(0,0,50e3),permanent=True)
O.forces.setPermF(upper_plateID,(0,0,100e3))
###############################################
#engines
###############################################

O.engines = [
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Box_Aabb()]),
        InteractionLoop(
                # interaction loop
                [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
                [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
                [Law2_ScGeom6D_CohFrictPhys_CohesionMoment(always_use_moment_law=True),Law2_ScGeom_FrictPhys_CundallStrack()]
        ),
        NewtonIntegrator(gravity=gravity,damping=0.2),
        ##TranslationEngine(translationAxis=[0,0,1],velocity=-1,ids=idsCyl3),
        TranslationEngine(translationAxis=[1,0,0],velocity=velocity_cyl_idsCyl2,ids=idsCyl2),
        TranslationEngine(translationAxis=[1,0,0],velocity=velocity_cyl_Box,ids=moveBoxIds),
        PyRunner(command='checkStress()', iterPeriod=10),
        #PyRunner(command='checkStress_NO()', iterPeriod=10),
        PyRunner(command='plotPlotData()',iterPeriod=1),
]

def checkStress():
 cylDisplacement_x = O.bodies[idsCyl2[0]].state.displ()[0]
 cylDisplacement_y = O.bodies[idsCyl2[0]].state.displ()[1]
 plot.addData(
 cylDisplacement_x=cylDisplacement_x,
 cylDisplacement_y=cylDisplacement_y,
 )
 print('{:1.04f}'.format(cylDisplacement_x)," ", '{:1.04f}'.format(cylDisplacement_y))

###############################################
#Definition of distortion
###############################################
def checkStress_1():
 if O.bodies[idsCyl2[0]].state.displ()[0] > 0.020:
  O.pause()

def checkStress_2():
 if O.bodies[idsCyl2[0]].state.displ()[1] > 0.020:
  O.pause()

def checkStress_3():
 if O.bodies[idsCyl2[0]].state.displ()[0] < -0.020:
  O.pause()

def checkStress_4():
 if O.bodies[idsCyl2[0]].state.displ()[1] < -0.020:
  O.pause()

def checkStress_5():
 if O.bodies[idsCyl2[0]].state.displ()[0] > 0.0020:
  O.pause()

###############################################
#output
###############################################
sp_num = len([b for b in O.bodies if isinstance(b.shape, Sphere)]) # 粒子の総数をカウント
print("number of spheres = ",sp_num)

O.dt=0.5*PWaveTimeStep()
def plotPlotData():
 pipe_force_x = sum([O.forces.f(i)[0] for i in idsCyl2])
 pipe_force_y = sum([O.forces.f(i)[1] for i in idsCyl2])
 cylDisplacementx = O.bodies[idsCyl2[0]].state.displ()[0]
 cylDisplacementy = O.bodies[idsCyl2[0]].state.displ()[1]
 cylDisplacementz = O.bodies[idsCyl2[0]].state.displ()[2]
 Dz=upper_plate.state.displ()[2]
 plot.addData(
 i=O.iter,
 disp=O.time*velocity_cyl,
 force_x=pipe_force_x,
 force_y=pipe_force_y,
 Dz=Dz,
 cylDisplacementx=cylDisplacementx,
 cylDisplacementy=cylDisplacementy,
 cylDisplacementz=cylDisplacementz,
 )
 plot.saveDataTxt('pipe_force.out',vars=('disp','force_x','force_y',
      'Dz',
      'cylDisplacementx','cylDisplacementx','cylDisplacementx'
      ))

plot.plots={
 'disp':('force_x','force_y'),
 'i':('Dz'),
 'cylDisplacementx':('cylDisplacementy',)
 }
plot.plot()

#plot.saveDataTxt(O.tags['id'] + '.txt')
############################################################

#program 2↓
############################################################
#program 2
O.engines = [
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Box_Aabb()]),
        InteractionLoop(
                # interaction loop
                [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
                [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
                [Law2_ScGeom6D_CohFrictPhys_CohesionMoment(always_use_moment_law=True),Law2_ScGeom_FrictPhys_CundallStrack()]
        ),
        NewtonIntegrator(gravity=gravity,damping=0.2),
        ##TranslationEngine(translationAxis=[0,0,1],velocity=-1,ids=idsCyl3),
        TranslationEngine(translationAxis=[1,0,0],velocity=velocity_cyl_idsCyl2,ids=idsCyl2),
        TranslationEngine(translationAxis=[1,0,0],velocity=velocity_cyl_Box,ids=moveBoxIds),
        PyRunner(command='checkStress()', iterPeriod=10),
        PyRunner(command='checkStress_NO()', iterPeriod=10),
        PyRunner(command='plotPlotData()',iterPeriod=1),
]
############################################################

Revision history for this message
内山康太郎 (kenkoutaro) said :
#4

Thank you.