yade-users team mailing list archive

[procrastinator <question684213@xxxxxxxxxxxxxxxxxxxxx>]

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

Hi all,
 I am trying to get the microstrain field. I did as the document said[1]. I want to use PyRunner command to export a series of .vtk files. It seems that I can just get one VTK file.  is there someone can help me to check the ### get the strain field part ### in the MWE.


from yade import pack, plot
from yade import export
from yade import utils
compFricDegree=45.0
finalFricDegree=19.5
rate=-15
damp=0.7
target_strain = 0.15
stabilityThreshold=0.01
young=4e8
mn,mx=Vector3(0,0,0),Vector3(4e-3,8e-3,4e-3)
O.materials.append(CohFrictMat(young=young,poisson=0.3,frictionAngle=radians(compFricDegree),isCohesive=False,alphaKr=0.2,alphaKtw=0,etaRoll=0.5,momentRotationLaw=True,density=2648,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=0,density=0,label='walls'))
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)
sp=pack.SpherePack()
psdSizes=[0.12e-3,0.14e-3,0.16e-3,0.18e-3,0.19e-3,0.20e-3,0.22e-3,0.24e-3,0.29e-3,0.38e-3]
psdCumm=[0,11.1,22.2,33.3,44.4,55.5,66.6,77.7,88.8,100]
sp.makeCloud(mn,mx,psdSizes=psdSizes,psdCumm=psdCumm,num=500,distributeMass=1,seed=1)
sp.toSimulation(material='spheres')
O.dt=.5*PWaveTimeStep()
triax=TriaxialStressController(
	maxMultiplier=1.0+2e5/young,
	finalMaxMultiplier=1.0+2e4/young,
	thickness=0,
	stressMask=7,
	internalCompaction=False
)
newton=NewtonIntegrator(damping=damp)
O.engines=[
	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
	InteractionLoop(
			[Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
			[Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
			[Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm=True,always_use_moment_law=True,label='cohesiveLaw')]
	),
	GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.8),
	triax,
	newton,
]
triax.goal1=triax.goal2=triax.goal3=-100000
while 1:
	O.run(3000,True)
	unb=unbalancedForce()
	print 'unbalanced force: ',unb,'mean stress: ',triax.meanStress
	if unb<stabilityThreshold and abs(-100000-triax.meanStress)/-100000<0.01:
		break
triax.wall_bottom_activated=True
triax.wall_top_activated=True
triax.wall_left_activated=True
triax.wall_right_activated=True
triax.wall_back_activated=True
triax.wall_front_activated=True
setContactFriction(radians(finalFricDegree))
re11=-triax.strain[0]
re22=-triax.strain[1]
re33=-triax.strain[2]
triax.stressMask = 5
newton=NewtonIntegrator(damping=0.0)
triax.goal2=rate
triax.strainRate2=15
triax.strainRate1=triax.strainRate3=1000.0
def stop_loading():
	if -triax.strain[1]-re22> target_strain:
		O.pause()
###################### get the strain field #####################################################
TW=TesselationWrapper()
TW.triangulate()
TW.computeVolumes()
TW.setState(0)
#step = 0
#while -triax.strain[1]-re22< target_strain:
	#step = step+500
	#O.run(400,True)
	#TW.setState(1)
	#TW.defToVtk("strain_"+str(step)+".vtk")
O.run(100,True)
TW.setState(1)
def strain_export():
	#O.run(100,True)
	TW.setState(1)
	TW.defToVtk("strain.vtk")
#################################################################################################
O.engines=[
	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
	InteractionLoop(
			[Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
			[Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
			[Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm=True,always_use_moment_law=True,label='cohesiveLaw')]
	),
	GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.8),
	triax,
	newton,
	PyRunner(command='stop_loading()',iterPeriod=10000),
	PyRunner(command='strain_export()',iterPeriod=100),
]	
O.run(50000,True)


I found that there is a question is almost the same as my question[2], they use a while function to solve this problem. But I want to know how can we solve this problem by using pyrunner function. thanks.

[1]https://github.com/yade/trunk/blob/master/examples/tesselationwrapper/tesselationWrapper.py
[2]https://answers.launchpad.net/yade/+question/675267

best,
yong

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