Ive made a model with spheres to represent pircas in Peru ( very similar to the function of a dry stone retaining walls) and now for comparing results and making my model more representative to the form of the rocks, Im using clumps . I use templates and when running it , what only appears is one sphere and not all clumps with the packing in the volume I want. At my first attempt using templates,geometrically it showed the way i wanted but once i run it , it appeared the "fatal error" . Then I change it to JCFpmState following recommendations ( https://answers.launchpad.net/yade/+question/686318 ) but it only appears one sphere as i mentioned. Maybe there is something with the code what its wrong. Im truly dont know how to make it work. Im copying the script in case i can have a hand with it. Thank you so much in advance. from yade import pack, qt, plot from numpy import arange import numpy as np import random #DATOS DE ENTRADA r=0.07 H=1 #altura del muro B=.45 #base mayor b=.25 #base menor L=4 #longitud del muro #EN SERVICIO EmpujeTop=2400#Pa2400 EmpujeBottom=9100#Pa9100 damp=0.7 aceloutofplane=0 ii=0 # GEOMETRIA ******************************************************************** a=(B-b)*.5 Vol1=pack.inParallelepiped(o=(0,0,0), a=(L,0,0), b=(0,B,0), c=(0,0,H)) Vol2=pack.inParallelepiped(o=(L,B+a,0), a=(0,B+a,0), b=(0,a+b,H), c=(L,a+b,H)) Vol3=pack.inParallelepiped(o=(L,-a,0), a=(0,-a,0), b=(0,a,H), c=(L,a,H)) Vol4=Vol1-Vol2-Vol3-Vol3 Esferas=pack.randomDensePack(Vol4,radius=r,rRelFuzz=0) O.bodies.append(Esferas) relRadList1 = [.04,.03] relPosList1 = [[0.035,0,0],[.083,0,0]] #peanut: relRadList2 = [.05,.05] relPosList2 = [[0,0,0],[.073,0,0]] relRadList4 = [0.025,0.025,0.025,0.025] relPosList4 = [[0,0,0],[.043,0,0],[0.02,0.045,0],[0.02,0,045]] #stick: relRadList3 = [0.05,0.025,0.05,0.025] relPosList3 = [[.04,0,0],[.123,0,0.022],[0.14,0,0],[0.123,0,-0.022]] templates= [] templates.append(clumpTemplate(relRadii=relRadList1,relPositions=relPosList1)) templates.append(clumpTemplate(relRadii=relRadList2,relPositions=relPosList2)) templates.append(clumpTemplate(relRadii=relRadList3,relPositions=relPosList3)) templates.append(clumpTemplate(relRadii=relRadList4,relPositions=relPosList4)) O.bodies.replaceByClumps(templates,[.5,.3,.1,.1]) # DEFINICION DEL MATERIAL ****************************************************** mat1 = JCFpmMat() mat1.cohesion =0 #Pa32.37e3 mat1.density = 2821 #kg/m3 mat1.frictionAngle = radians(58.74) #rad mat1.poisson = 0.18 mat1.young = 37.04e9 #Pa # DEFINICION DE LA JUNTA ****************************************************** mat1.jointCohesion = 262000 #Pa mat1.jointFrictionAngle = radians(40) mat1.jointNormalStiffness= 1.75e8 #Pa/m mat1.jointShearStiffness=0.7e8 #Pa/m O.materials.append(mat1) # APLICACION DE FUERZAS ******************************************************** listSphereA=[] listSphereB=[] listSphereC=[] listSphereD=[] for x1 in O.bodies: if isinstance(x1.shape,Sphere) and x1.state.pos[1]<3*r and x1.state.pos[2]<=H and x1.state.pos[2]>3*H/4.000: listSphereA.append(x1.id) for x2 in O.bodies: if isinstance(x2.shape,Sphere) and x2.state.pos[1]<3*r and x2.state.pos[2]<=3*H/4.000 and x2.state.pos[2]>H/2.000: listSphereB.append(x2.id) for x3 in O.bodies: if isinstance(x3.shape,Sphere) and x3.state.pos[1]<3.5*r and x3.state.pos[2]<=H/2.000 and x3.state.pos[2]>H/4.000: listSphereC.append(x3.id) for x4 in O.bodies: if isinstance(x4.shape,Sphere) and x4.state.pos[1]<3*r and x4.state.pos[2]<=H/4.000 and x4.state.pos[2]>0.000: listSphereD.append(x4.id) ha=1*H/8.000 hb=3*H/8.000 hc=5*H/8.000 hd=7*H/8.000 m=EmpujeBottom-EmpujeTop Ea=m*ha+EmpujeTop Eb=m*hb+EmpujeTop Ec=m*hc+EmpujeTop Ed=m*hd+EmpujeTop Fa=Ea*(H/4.000)*L/len(listSphereA) Fb=Eb*(H/4.000)*L/len(listSphereB) Fc=Ec*(H/4.000)*L/len(listSphereC) Fd=Ed*(H/4.000)*L/len(listSphereD) def aplicarFuerzaA(): for i in listSphereA: O.forces.setPermF(i,(0,Fa,0)) def aplicarFuerzaB(): for i in listSphereB: O.forces.setPermF(i,(0,Fb,0)) def aplicarFuerzaC(): for i in listSphereC: O.forces.setPermF(i,(0,Fc,0)) def aplicarFuerzaD(): for i in listSphereD: O.forces.setPermF(i,(0,Fd,0.98)) # CONDICIONES DE BORDE ********************************************************* Borde1=utils.wall(0.43,axis=0, sense=0) Borde2=utils.wall(L-0.43,axis=0, sense=0) Borde3=utils.wall(0, axis=1, sense=0) Borde4=utils.wall(0.05, axis=2, sense=0) O.bodies.append(Borde1) O.bodies.append(Borde2) O.bodies.append(Borde3) O.bodies.append(Borde4) listSphereBase=[] for i in listSphereBase: O.bodies[i].state.blockedDOFs='xyz' # DEFINICION DE FUNCIONES ****************************************************** for m in O.bodies: if isinstance(m.shape,Sphere) and m.state.pos[2]<2*H: listSphereBase.append(m.id) listSphereTodo=[] for i in listSphereTodo: O.bodies[i].state.blockedDOFs='XYZ' def checkUnbalanced(): print ('iter %d, unbalanced forces = %f'%(O.iter, utils.unbalancedForce())) print 'Fzas en c/ esfera (N)', Fa, Fb, Fc, Fd print ( 'Desplazamiento superior',O.bodies[w].state.pos[1],O.bodies[w].state.pos[2]) listTopBodies=[] for z1 in O.bodies: if isinstance(z1.shape,Sphere) and z1.state.pos[2]<=H and z1.state.pos[2]>H-3.000*r and z1.state.pos[0]>(L/2.000)-1.1*r and z1.state.pos[0]<(L/2.000)+1.1*r and z1.state.pos[1]>B/2.000: listTopBodies.append(z1.id) w=max(listTopBodies) def addPlotData(): plot.addData(Ypos=O.bodies[w].state.pos[1],Zpos=O.bodies[w].state.pos[2],iteracion=O.iter,unbal=unbalancedForce()) plot.plots={'Ypos':('Zpos')} # ENGINE *********************************************************************** O.engines=[ ForceResetter(), InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Wall_Aabb()]), InteractionLoop( [Ig2_Sphere_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()], [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys()], [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM()] ), NewtonIntegrator(gravity=(0,aceloutofplane,-9.81),damping=damp), PyRunner(firstIterRun=10000,command='aplicarFuerzaA()'), PyRunner(firstIterRun=10000,command='aplicarFuerzaB()'), PyRunner(firstIterRun=10000,command='aplicarFuerzaC()'), PyRunner(firstIterRun=10000,command='aplicarFuerzaD()'), PyRunner(iterPeriod=10000,command='checkUnbalanced()'), PyRunner(iterPeriod=1,command='addPlotData()'), VTKRecorder(iterPeriod=100,recorders=['spheres','jcfpm','facets','colors'],fileName='/tmp/p1-') ] for b in O.bodies: b.mat=mat1 b.state=JCFpmState() # DETALLES FINALES ************************************************************* O.dt=0.5*PWaveTimeStep() # establece el delta de tiempo como una fraccion del tiempo critico plot.plot(subPlots=True) # llama a los sub plots plot.live=True # ploteo en tiempo real qt.Controller() # abre la ventana del controlador V=qt.View() # abre la ventada de la vista R=yade.qt.Renderer() # llama a la renderizacion R.bgColor=(1.,1.,1.) # definel color blan V.screenSize = (900,900) # tamano de pantalla

Yade

History of question #687565