Good evening, I'm from Peru and I'm part of an investigation that tries to represent the pirca performance at seismic action. The pirca is a dry stone retaining very used here. For making the model, I was trying to use spherepack with clumps to represent the blocks of rock form. First, If i use a volume like a regular prism using vector,it works perfectly but what i want is a trapezoidal form (cross section) of the dry stone retaining wall. I tried many times but i couldnt find a way to make it work,thats why im sharing the code i was working on and hopefully i can have a hand with this problem. Is there a way to put the clumps in a volume with a trapezoidal form as cross section? Second,I also want to recognize the clump that has the greatest displacement thats why i wrote that code . I would like to know if there is another way to do it. listTopBodies=[] for z1 in O.bodies: if isinstance(z1.shape,Sphere) and z1.state.pos[2]<=H+0.1 and z1.state.pos[2]>H-3.000*r and z1.state.pos[0]>(L/2.000)-2.000*r and z1.state.pos[0]<(L/2.000)+2.000*r and z1.state.pos[1]>B/2.000: listTopBodies.append(z1.id) i=max(listTopBodies) def addPlotData(): plot.addData(Ypos=O.bodies[i].state.pos[1],Zpos=O.bodies[i].state.pos[2],iteracion=O.iter,unbal=unbalancedForce()) plot.plots={'Ypos':('Zpos'),'iteracion':('unbal')} -------------------------------------------------------------------------------------------------------------------------------------------------------------------- Here i share the complete code for the model: # MODELO NUMERICO CON CLUMPS ************************* 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=2 #longitud del muro #EN SERVICIO EmpujeTop=2400 #Pa EmpujeBottom=9100 #Pa damp=0.7 aceloutofplane=0 ii=0 # Definicion del material mat1 para la pirca ********** mat1 = JCFpmMat() mat1.cohesion = 0 #32.37e3 #Pa mat1.density = 2233.2 #kg/m3 mat1.frictionAngle = radians(58.74) #rad mat1.poisson = 0.18 mat1.young = 37.04e9 #Pa # DEFINICION DE LA JUNTA ****************************************************** mat1.jointCohesion = 232500 mat1.jointFrictionAngle = radians(40) mat1.jointNormalStiffness= 1.75e8 mat1.jointShearStiffness=0.875e8 #mat1.jointDilationAngle = 0 #mat1.jointTensileStrength = 0 O.materials.append(mat1) # 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 # Conjunto de esferas representadas como centros y radios sp=pack.SpherePack() biesfera=pack.SpherePack([((0,0,0),.033),((.033,0,0),.033)]) # clump de dos esferas triesfera=pack.SpherePack([((0,0,0),.04),((.02,.02,0),.04),((0,.02,.02),.04)]) # clump de tres esferas tetraesfera=pack.SpherePack([((0,0,0),.033),((.033,0,0),.033),((0,.033,0),.033),((0,0,.033),.033)]) # clump decuatro esferas sp.makeClumpCloud(Vol4,[biesfera,triesfera, tetraesfera]) # Empaquetamiento suelto y aleatoriode clumps sp.toSimulation() # llama las esferas a la simulacion #Imprime la cantidad total de cuerpos Q=len(O.bodies) print 'numero de cuerpos', Q # Definicion del piso piso=utils.wall(position=Vector3(0,0,0), axis=2, sense=0, material=mat1, color=(0.5,0.5,0.5)) O.bodies.append(piso) # Definicion de las fuerzas en aplicacion 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+0.1 and x2.state.pos[2]>H/2.000+0.1: listSphereB.append(x2.id) for x3 in O.bodies: if isinstance(x3.shape,Sphere) and x3.state.pos[1]<3*r and x3.state.pos[2]<=H/2.000+0.1 and x3.state.pos[2]>H/4.000+0.1: 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+0.1 and x4.state.pos[2]>0.000+0.1: 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) print 'Fzas en c/ esfera (N)', Fa, Fb, Fc, Fd 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)) # CONDICIONES DE BORDE ********************************************************* Base=utils.wall(0.1,axis=2,sense=0) Borde1=utils.wall(0.22, axis=0, sense=0) Borde2=utils.wall(L-.22, axis=0, sense=0) Borde3=utils.wall(-.0, axis=1, sense=0) O.bodies.append(Base) O.bodies.append(Borde1) O.bodies.append(Borde2) O.bodies.append(Borde3) listSphereBase=[] for y1 in O.bodies: if isinstance(y1.shape,Sphere) and y1.state.pos[2]<=3*r+0.1: listSphereBase.append(y1.id) #for y2 in O.bodies: # if isinstance(y2.shape,Sphere) and y2.state.pos[0]<=2.5*r: # listSphereBorde1.append(y2.id) #for y3 in O.bodies: # if isinstance(y3.shape,Sphere) and y3.state.pos[0]>L-2.5*r: # listSphereBorde2.append(y3.id) for i in listSphereBase: O.bodies[i].state.blockedDOFs='xyzXYZ' #for i in listSphereBorde1: # O.bodies[i].state.blockedDOFs='xyzXYZ' #for i in listSphereBorde2: # O.bodies[i].state.blockedDOFs='xyzXYZ' def checkUnbalanced(): print ('iter %d, unbalanced forces = %f'%(O.iter, utils.unbalancedForce())) # if unbalancedForce()<.00001: listTopBodies=[] for z1 in O.bodies: if isinstance(z1.shape,Sphere) and z1.state.pos[2]<=H+0.1 and z1.state.pos[2]>H-3.000*r and z1.state.pos[0]>(L/2.000)-2.000*r and z1.state.pos[0]<(L/2.000)+2.000*r and z1.state.pos[1]>B/2.000: listTopBodies.append(z1.id) i=max(listTopBodies) def addPlotData(): plot.addData(Ypos=O.bodies[i].state.pos[1],Zpos=O.bodies[i].state.pos[2],iteracion=O.iter,unbal=unbalancedForce()) plot.plots={'Ypos':('Zpos'),'iteracion':('unbal')} # MOTOR DEL ALGORITMO DEM ******************************* 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()') ] # 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,500) # tamano de pantalla

When	Who	Action	New State
2019-10-25	Jérôme Duriez	Comment	Open
2019-10-25	Grace Mejico	Comment	Open
2019-10-30	Grace Mejico	Comment	Open
2019-11-09	Launchpad Janitor	Expire	Expired

Yade

History of question #685401