Question #691407 on Yade changed: https://answers.launchpad.net/yade/+question/691407
Chien-Cheng Hung gave more information on the question: Hi Jan, Thanks for your email. I've tested your microstrain vtk file and it works with my Paraview. So I think the crash is not related to the Paraview version (btw my Paraview version is 5.0.1 (64-bit)). Maybe it's related to my code. Could you help me have a look? I'd appreciate that. ### My simplified code (direct shear simulation of a granular layer) from yade import pack,plot,export import numpy as np import math sp1=pack.SpherePack() sp2=pack.SpherePack() sp3=pack.SpherePack() O.periodic=True # dimensions of sample RADIUS=0.25 a=15 b=1 c=1 length=a*(2*RADIUS) height=length/b width=length/c thickness=RADIUS ### Guassian distribution psdSizes=[.456,.5,.544] psdCumm=[0,0.5,1] # friction angles wallFRIC=0 boundaryFRIC=0.5 spFRIC=0.5 # boundary conditions PI=1.e5 SN=5.e6 # normal stress RATE_NS1=1 # velocity of top plate during compaction RATE_NS2=1 # velocity of top plate during shear RATE_shear=1 # shear velocity roll_stiff=0 roll_fric=0 # simulation control DAMPSHEAR=0. O.cell.hSize=Matrix3(length,0,0,0,3*height,0,0,0,width) O.materials.append(FrictMat(density=3000,young=1e8,poisson=0.5,frictionAngle=wallFRIC,label='boxMat')) O.materials.append(FrictMat(density=3000,young=1e8,poisson=0.5,frictionAngle=boundaryFRIC,label='boundaryMat')) O.materials.append(FrictMat(density=3000,young=1e8,poisson=0.5,frictionAngle=spFRIC,label='sphereMat')) upBox = utils.box(center=(length/2,2*height+thickness,1.5*width),orientation=Quaternion(1,0,0,0),extents=(2*length,thickness/2.,width),fixed=1,wire=False,color=(1,0,0),material='boxMat') lowBox = utils.box(center=(length/2,height-thickness,1.5*width),orientation=Quaternion(1,0,0,0),extents=(2*length,thickness/2.,width),fixed=1,wire=False,color=(1,0,0),material='boxMat') O.bodies.append([upBox,lowBox]) sp1.makeCloud((0,height+3*RADIUS,width),(length,2*height-3*RADIUS,2*width), psdSizes =psdSizes, psdCumm =psdCumm, periodic=True, seed =1) sp2.makeCloud((0,height+RADIUS,width),(length,height+RADIUS-1e-10,2*width), rMean=RADIUS, periodic=True, seed =1) sp3.makeCloud((0,2*height-RADIUS,width),(length,2*height-RADIUS-1e-10,2*width), rMean=RADIUS, periodic=True, seed =1) sphere_id = O.bodies.append([utils.sphere(s[0],s[1],color=(0,0,1),material='sphereMat') for s in sp1]) bottomLayer_id = O.bodies.append([utils.sphere(s[0],s[1],color=(1,0,1),material='boundaryMat') for s in sp2]) topLayer_id = O.bodies.append([utils.sphere(s[0],s[1],color=(1,0,1),material='boundaryMat') for s in sp3]) effCellVol=(O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2] volRatio=(O.cell.hSize[0,0]*O.cell.hSize[1,1]*O.cell.hSize[2,2])/effCellVol O.engines=[ ForceResetter() ,InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb()],verletDist=-0.1,allowBiggerThanPeriod=True) ,InteractionLoop( [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D()], [Ip2_FrictMat_FrictMat_MindlinPhys(krot=roll_stiff,eta=roll_fric)], [Law2_ScGeom_MindlinPhys_Mindlin(includeMoment=True)] ) ,GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=1,timestepSafetyCoefficient=0.8,defaultDt=-1) ,PyRunner(command='fixVelocity(RATE_shear)',iterPeriod=1,label='fixVel',dead=True) ,PeriTriaxController(dynCell=True,mass=10,maxUnbalanced=1e-3,relStressTol=1e-4,stressMask=7,goal=(-PI/volRatio,-PI/volRatio,-PI/volRatio),globUpdate=1,maxStrainRate=(1,1,1),doneHook='triaxDone()',label='triax') ,NewtonIntegrator(gravity=(0,0,0),damping=0.3,label='newton') ] def triaxDone(): global phase volRatio=(O.cell.hSize[0,0]*O.cell.hSize[1,1]*O.cell.hSize[2,2])/((O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2]) h=O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1] vol=h*O.cell.hSize[0,0]*O.cell.hSize[2,2] contactStress=getStress(vol) vol_s=Rmean=Rmax=nbSph=0 Rmin=1e6 x_ref=O.bodies[0].state.pos[0] for o in O.bodies: if isinstance(o.shape,Sphere): nbSph+=1 Rmean+=o.shape.radius if o.shape.radius>Rmax: Rmax=o.shape.radius if o.shape.radius<Rmin: Rmin=o.shape.radius vol_s += 4.*pi/3.*(o.shape.radius)**3 Rmean=Rmean/nbSph n = 1-vol_s/vol print ('DONE! iter=',O.iter,'| sample generated: nb spheres',nbSph,', Rmean=',Rmean,', Rratio=',Rmax/Rmin,', porosity=',n,'wall_refPos=',x_ref) tt=TriaxialCompressionEngine() #tt.setContactProperties(shearFRIC) triax.dead=True O.pause() #### If top particle layer contact with wall then give the particle velocity, and vice versa. def identify_topLayer(): top = [] for i in O.interactions.withBody(0): if i.id1 != 0: top.append(i.id1) if i.id2 != 0: top.append(i.id2) return top O.run(10000000,1) stage=0 stiff=fnPlaten=currentSN=0. def servo(): global stage,stiff,fnPlaten,currentSN if stage==0: currentSN=O.forces.f(0)[1]/(O.cell.hSize[0,0]*O.cell.hSize[2,2]) unbF=unbalancedForce() print ('SN=',SN,'| current SN = ',currentSN,' | unbF=',unbF ) boundaryVel=copysign(min(RATE_NS1,abs(0.5*(currentSN-SN))),currentSN-SN) O.bodies[0].state.vel[1]=boundaryVel if ( (abs(currentSN-SN)/SN)<0.001 and unbF<0.005 ): stage+=1 fnPlaten=O.forces.f(0)[1] ## the following computes the stiffness of the plate (used for stress control of the top plate) for i in O.interactions.withBody(O.bodies[0].id): stiff+=i.phys.kn O.pause() if stage==1: fnDesired=SN*(O.cell.hSize[0,0]*O.cell.hSize[2,2]) boundaryVel=copysign(min(RATE_NS2,abs(0.333*(O.forces.f(0)[1]-fnDesired)/stiff/O.dt)),O.forces.f(0)[1]-fnDesired) O.bodies[0].state.vel[1]=boundaryVel O.engines = O.engines[:5]+[PyRunner(command='servo()',iterPeriod=1,label='servo')]+O.engines[5:] O.run(10000000,1) #### fixed bottom bottomLayer_id for i in bottomLayer_id: O.bodies[i].state.blockedDOFs='xyzXYZ' O.bodies[i].state.dynamics = False O.bodies[i].state.vel = Vector3(0,0,0) ### Start shearing newton.damping=DAMPSHEAR fixVel.dead = False refTime = O.time def fixVelocity(RATE_shear): O.bodies[0].state.vel[0] = RATE_shear topLayer_id = identify_topLayer() for i in topLayer_id: O.bodies[i].state.vel[0] = RATE_shear slip = (O.time-refTime)*RATE_shear h=O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1] ss = slip/h if ss > 1: O.pause() O.run(10000,1) TW=TesselationWrapper() TW.triangulate() TW.computeVolumes() TW.volume(10) TW.setState(0) O.run(100,True) TW.setState(1) TW.defToVtk("strain.vtk") ### Cheers, Chien-Cheng -- 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