← Back to team overview

yade-users team mailing list archive

[Question #687565]: problems when using templates to generate a clumps model

 

New question #687565 on Yade:
https://answers.launchpad.net/yade/+question/687565

 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

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.