New question #683821 on Yade: https://answers.launchpad.net/yade/+question/683821
Hi, I'm trying to generate a cylinder sample with target porosity, the MWE comes from [1]. The main changes are: 1. changed the sphere material from FrictMat to CohFrictMat, as well as corresponding engines, 2. cut the cubic sample to fulfil a cylinder shape after target porosity is reached. The problem is after the simulation finished (run about 30 seconds), I can not find any particles in the scene, while I use 'len(O.bodies)' to check, I can get that there are 870 bodies. Here is the MWE: ######## from yade import pack num_spheres=1000 key='_triax_base_' targetPorosity = 0.43 compFricDegree = 30 finalFricDegree = 30 rate=-0.02 damp=0.2 stabilityThreshold=0.01 young=5e6 mn,mx=Vector3(0,0,0),Vector3(0.07,0.07,0.14) cement = CohFrictMat(young=30e9,poisson=0.3,frictionAngle=radians(30),density=2650.0,isCohesive=True,normalCohesion=1e8, shearCohesion=1e8,label='cement') walls = FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls') O.materials.append(cement) O.materials.append(walls) walls=aabbWalls([mn,mx],thickness=0,material='walls') wallIds=O.bodies.append(walls) sp=pack.SpherePack() x=[5.0,6.0005088055510125,6.807147320174279,7.200507128170075,7.399996645238125,7.604595185078019,7.809172757656196,8.008676252898724,8.20306373628217,8.41275033163219,9.005096442414814] y=[0,9.972588799847912,20.08747541588565,29.9397121334748,39.93424666725376,50.065627529175956,59.78605006143409,70.0545568149891,80.04923113051282,89.9064999909142,100.0] #### change the size units psdSizes=x for i in range(0,len(x)): psdSizes[i]=x[i]*0.001 ## change the unit to percentage psdCumm=y for i in range(0,len(y)): psdCumm[i]=y[i]*0.01 sp.makeCloud(mn,mx,psdSizes=psdSizes,psdCumm=psdCumm,distributeMass=True,seed=1) sp.toSimulation(material='cement') Gl1_Sphere.quality=3 triax=TriaxialStressController( maxMultiplier=1.+2e4/young, finalMaxMultiplier=1.+2e3/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_CohFrictMat_CohFrictMat_CohFrictPhys(),Ip2_FrictMat_FrictMat_FrictPhys()], [Law2_ScGeom6D_CohFrictPhys_CohesionMoment(),Law2_ScGeom_FrictPhys_CundallStrack()] ), GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8), triax, TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+key), newton ] triax.goal1=triax.goal2=triax.goal3=-10000 while 1: O.run(1000, True) unb=unbalancedForce() print 'unbalanced force:',unb,' mean stress: ',triax.meanStress if unb<stabilityThreshold and abs(-10000-triax.meanStress)/10000<0.001: break print "### Isotropic state saved ###" import sys while triax.porosity>targetPorosity: compFricDegree = 0.95*compFricDegree setContactFriction(radians(compFricDegree)) print "\r compFrictionDegree: ",compFricDegree," porosity:",triax.porosity, sys.stdout.flush() O.run(500,1) print "### Compacted state saved ###" aabb=utils.aabbExtrema() predX=(aabb[0][0]+aabb[1][0])/2 predY=(aabb[0][1]+aabb[1][1])/2 predR=min((aabb[1][0]-aabb[0][0])/2,(aabb[1][1]-aabb[0][1])/2) pred=pack.inCylinder((predX,predY,0),(predX,predY,0.14),predR) ## erase Walls and spheres which is not in cylinder def cylinderShape(): for i in O.bodies: if not isinstance(i.shape,Sphere): O.bodies.erase(i.id) if isinstance(i.shape,Sphere): if not pred(i.state.pos, i.shape.radius): # if sphere with c,r is in pred(cylinder), then add this sphere into bodies. O.bodies.erase(i.id) cylinderShape() ######## Could you please help me with this problem? Thanks, Leonard [1]https://github.com/yade/trunk/blob/master/examples/triax-tutorial/script-session1.py -- 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