Question #706766 on Yade changed: https://answers.launchpad.net/yade/+question/706766
Karol Brzezinski proposed the following answer: Hi Leonard, Try this (please note that I am using python 3, so the call pf print function is modified): ########## ################### 2D MWE ############# from __future__ import division from yade import pack, plot, export, ymport num_spheres=1000 rate=-0.01 damp=0.6 stabilityThreshold=0.001 young=5e6 confinement=6.7e3 mn,mx=Vector3(0,0,0.02),Vector3(1,2,0.02) O.materials.append(FrictMat(young=young, poisson=0.5, frictionAngle=radians(30), density=2600, label='spheres')) O.materials.append(FrictMat(young=young, poisson=0.5, frictionAngle=0, density=0, label='walls')) walls=aabbWalls([Vector3(0,0,0),Vector3(1,2,0.04)],thickness=0,material='walls') wallIds=O.bodies.append(walls) sp=pack.SpherePack() sp.makeCloud(mn,mx,-1,0.3333,num_spheres,False, 0.75,seed=1) O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp]) Gl1_Sphere.quality=3 for b in O.bodies: if isinstance(b.shape,Sphere): b.state.blockedDOFs='zXY' b.shape.color=[1,1,1] triax=TriaxialStressController( thickness = 0, stressMask = 7, internalCompaction=False, ) newton=NewtonIntegrator(damping=damp) O.engines = [ ForceResetter(), InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Box_Aabb()]), InteractionLoop([Ig2_Sphere_Sphere_ScGeom(), Ig2_Box_Sphere_ScGeom()], [Ip2_FrictMat_FrictMat_FrictPhys()], [Law2_ScGeom_FrictPhys_CundallStrack()]), GlobalStiffnessTimeStepper(active=1, timeStepUpdateInterval=100, timestepSafetyCoefficient=0.8), triax, newton ] triax.goal1=triax.goal2=triax.goal3=-confinement while 1: O.run(1000, True) unb=unbalancedForce() print('unbF:',unb,' meanStress: ',-triax.meanStress,'top:',-triax.stress(triax.wall_top_id)[1],'left:',-triax.stress(triax.wall_left_id)[0],'front:',-triax.stress(triax.wall_front_id)[2]) if unb<stabilityThreshold and abs(-confinement-triax.meanStress)/confinement<0.0001: break print("### state 1 completed ###") ######## export files and run simulation export.text('state_0.txt')## export files state 0 O.run(500,True) export.text('state_1.txt')## export files state 1 ###################### layerThickness = 0.015*2 # this is more or less diameter of your particle O.bodies.clear() for i in range(3): bb = ymport.text('state_0.txt', shift = Vector3(0,0,i*layerThickness)) O.bodies.append(bb) TW=TesselationWrapper() TW.triangulate() TW.computeVolumes() TW.setState(0) O.bodies.clear() for i in range(3): bb = ymport.text('state_1.txt', shift = Vector3(0,0,i*layerThickness)) O.bodies.append(bb) TW.setState(1) TW.defToVtk("strain_pseudo2D.vtk") ####### Cheers, Karol -- You received this question notification because your team yade-users is an answer contact for Yade. _______________________________________________ Mailing list: https://launchpad.net/~yade-users Post to : yade-users@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp