Commit ba5647f7 authored by François Kneib's avatar François Kneib Committed by Janek Kozicki

[Python3] Clean examples as python3 doesn't deal with mix between spaces and tabs in indent lines.

parent 3385d1b9
......@@ -25,9 +25,9 @@ O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom_FrictPhys_CundallStrack()]
[Ig2_Sphere_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom_FrictPhys_CundallStrack()]
),
NewtonIntegrator(damping=0.1, gravity=(0.,0.,-10)),
]
......
......@@ -116,7 +116,7 @@ b5=utils.box(center=[center[0],center[1],center[2]],extents=[halfSize[0],halfSiz
O.bodies.append(b5)
#--
center=((lowerCornerW[0]+upperCornerW[0])/2,(lowerCornerW[1]+upperCornerW[1])/2,upperCornerW[2]+thickness/2.0)
halfSize=(wallOversizeFactor*fabs(lowerCornerW[0]-upperCornerW[0])/2+thickness,wallOversizeFactor*fabs(lowerCornerW[1]-upperCornerW[1])/2+thickness,thickness/2.0);
halfSize=(wallOversizeFactor*fabs(lowerCornerW[0]-upperCornerW[0])/2+thickness,wallOversizeFactor*fabs(lowerCornerW[1]-upperCornerW[1])/2+thickness,thickness/2.0);
b6=utils.box(center=[center[0],center[1],center[2]],extents=[halfSize[0],halfSize[1],halfSize[2]],color=[0,1,0],fixed=True,wire=True,material='walls')
O.bodies.append(b6)
......
......@@ -92,14 +92,14 @@ from yade import plot
O.engines=O.engines+[PyRunner(iterPeriod=20,command='history()',dead=1,label='recorder')]
def history():
plot.addData(e11=-triax.strain[0]-ei0, e22=-triax.strain[1]-ei1, e33=-triax.strain[2]-ei2,
s11=-triax.stress(0)[0]-si0,
s22=-triax.stress(2)[1]-si1,
s33=-triax.stress(4)[2]-si2,
pc=-unsat.bndCondValue[2],
sw=unsat.getSaturation(False),
i=O.iter
)
plot.addData(e11=-triax.strain[0]-ei0, e22=-triax.strain[1]-ei1, e33=-triax.strain[2]-ei2,
s11=-triax.stress(0)[0]-si0,
s22=-triax.stress(2)[1]-si1,
s33=-triax.stress(4)[2]-si2,
pc=-unsat.bndCondValue[2],
sw=unsat.getSaturation(False),
i=O.iter
)
plot.plots={'pc':('sw',None,'e22')}
plot.plot()
......
......@@ -148,8 +148,8 @@ from yade import plot
## a function saving variables
def history():
plot.addData(e22=-triax.strain[1]-zeroe22,e22_theory=drye22+(1-dryFraction)*consolidation((O.time-zeroTime)*Cv/wetHeight**2)*1000./modulus,t=O.time,p=flow.getPorePressure((0.5,0.1,0.5)),s22=-triax.stress(3)[1]-10000)
#plot.addData(e22=-triax.strain[1],t=O.time,s22=-triax.stress(2)[1],p=flow.MeasurePorePressure((0.5,0.5,0.5)))
plot.addData(e22=-triax.strain[1]-zeroe22,e22_theory=drye22+(1-dryFraction)*consolidation((O.time-zeroTime)*Cv/wetHeight**2)*1000./modulus,t=O.time,p=flow.getPorePressure((0.5,0.1,0.5)),s22=-triax.stress(3)[1]-10000)
#plot.addData(e22=-triax.strain[1],t=O.time,s22=-triax.stress(2)[1],p=flow.MeasurePorePressure((0.5,0.5,0.5)))
O.engines=O.engines+[PyRunner(iterPeriod=200,command='history()',label='recorder')]
##make nice animations:
......
......@@ -78,9 +78,9 @@ O.engines = [
InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Wall_Aabb(),Bo1_Facet_Aabb(),Bo1_Box_Aabb()],label='contactDetection',allowBiggerThanPeriod = True),
# Calculate the different interactions
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(), Ig2_Box_Sphere_ScGeom()],
[Ip2_ViscElMat_ViscElMat_ViscElPhys()],
[Law2_ScGeom_ViscElPhys_Basic()]
[Ig2_Sphere_Sphere_ScGeom(), Ig2_Box_Sphere_ScGeom()],
[Ip2_ViscElMat_ViscElMat_ViscElPhys()],
[Law2_ScGeom_ViscElPhys_Basic()]
,label = 'interactionLoop'),
#Apply an hydrodynamic force to the particles
HydroForceEngine(densFluid = densFluidPY,viscoDyn = kinematicViscoFluid*densFluidPY,zRef = groundPosition,gravity = gravityVector,deltaZ = fluidHeight/ndimz,expoRZ = 0.,lift = False,nCell = ndimz,vCell = 1.,vxFluid = np.zeros(ndimz),phiPart = np.zeros(ndimz),vxPart = np.zeros(ndimz),vFluctX = np.zeros(len(O.bodies)),vFluctY = np.zeros(len(O.bodies)),vFluctZ = np.zeros(len(O.bodies)),ids = idApplyForce, label = 'hydroEngine'),
......
......@@ -117,9 +117,9 @@ O.engines = [
InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Wall_Aabb(),Bo1_Facet_Aabb(),Bo1_Box_Aabb()],label='contactDetection',allowBiggerThanPeriod = True),
# Calculate the different interactions
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(), Ig2_Box_Sphere_ScGeom()],
[Ip2_ViscElMat_ViscElMat_ViscElPhys()],
[Law2_ScGeom_ViscElPhys_Basic()]
[Ig2_Sphere_Sphere_ScGeom(), Ig2_Box_Sphere_ScGeom()],
[Ip2_ViscElMat_ViscElMat_ViscElPhys()],
[Law2_ScGeom_ViscElPhys_Basic()]
,label = 'interactionLoop'),
#Apply an hydrodynamic force to the particles
HydroForceEngine(densFluid = densFluidPY,viscoDyn = kinematicViscoFluid*densFluidPY,zRef = groundPosition,gravity = gravityVector,deltaZ = dz,expoRZ = expoDrag_PY,lift = False,nCell = ndimz,vCell = length*width*dz ,vxFluid = vxFluidPY,phiPart = phiPartPY,vxPart = vxPartPY,ids = idApplyForce,vFluctX = np.zeros(len(O.bodies)),vFluctY = np.zeros(len(O.bodies)),vFluctZ = np.zeros(len(O.bodies)), label = 'hydroEngine', dead = True),
......@@ -152,7 +152,7 @@ def gravityDeposition(lim):
else :
print('\n Gravity deposition finished, apply fluid forces !\n')
newtonIntegr.damping = 0.0 # Set the artificial numerical damping to zero
gravDepo.dead = True # Remove the present engine for the following
gravDepo.dead = True # Remove the present engine for the following
hydroEngine.dead = False # Activate the HydroForceEngine
turbFluct.dead = False # Activate the turbulent fluctuations model
return
......
......@@ -135,10 +135,10 @@ O.engines = [
InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Wall_Aabb(),Bo1_Facet_Aabb(),Bo1_Box_Aabb()],label='contactDetection',allowBiggerThanPeriod = True),
# Calculate the different interactions
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(), Ig2_Box_Sphere_ScGeom()],
[Ip2_ViscElMat_ViscElMat_ViscElPhys()],
[Law2_ScGeom_ViscElPhys_Basic()]
,label = 'interactionLoop'),
[Ig2_Sphere_Sphere_ScGeom(), Ig2_Box_Sphere_ScGeom()],
[Ip2_ViscElMat_ViscElMat_ViscElPhys()],
[Law2_ScGeom_ViscElPhys_Basic()]
,label = 'interactionLoop'),
#Apply an hydrodynamic force to the particles
HydroForceEngine(densFluid = densFluidPY,viscoDyn = kinematicViscoFluid*densFluidPY,zRef = groundPosition,gravity = gravityVector,deltaZ = dz,expoRZ = expoDrag_PY,nCell = ndimz,vCell = length*width*dz ,vxFluid = vxFluidPY,phiPart = phiPartPY,vxPart = vxPartPY,ids = idApplyForce, label = 'hydroEngine', dead = True),
#Measurement, output files
......@@ -170,7 +170,7 @@ def gravityDeposition(lim):
else :
print('\n Gravity deposition finished, apply fluid forces !\n')
newtonIntegr.damping = 0.0 # Set the artificial numerical damping to zero
gravDepo.dead = True # Remove the present engine for the following
gravDepo.dead = True # Remove the present engine for the following
hydroEngine.dead = False # Activate the HydroForceEngine
hydroEngine.vxFluid = vxFluidPY # Send the fluid velocity vector used to apply the drag fluid force on particles in HydroForceEngine (see c++ code)
......
......@@ -145,9 +145,9 @@ O.engines = [
InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Wall_Aabb(),Bo1_Facet_Aabb(),Bo1_Box_Aabb()],label='contactDetection',allowBiggerThanPeriod = True),
# Calculate the different interactions
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(), Ig2_Box_Sphere_ScGeom()],
[Ip2_ViscElMat_ViscElMat_ViscElPhys()],
[Law2_ScGeom_ViscElPhys_Basic()]
[Ig2_Sphere_Sphere_ScGeom(), Ig2_Box_Sphere_ScGeom()],
[Ip2_ViscElMat_ViscElMat_ViscElPhys()],
[Law2_ScGeom_ViscElPhys_Basic()]
,label = 'interactionLoop'),
#Apply an hydrodynamic force to the particles
HydroForceEngine(densFluid = densFluidPY,viscoDyn = dynamicVisco,zRef = groundPosition,gravity = gravityVector,deltaZ = dz,expoRZ = expoDrag_PY,lift = False,nCell = ndimz,vCell = length*width*dz,radiusPart=diameterPart/2.,vxFluid = initVxFluid,phiPart = phiPartPY,vxPart = vxPartPY,irheolf=0, iusl = 0, uTop = uTop,iturbu = 0,ids = idApplyForce, label = 'hydroEngine', dead = True),
......@@ -181,7 +181,7 @@ def gravityDeposition(lim):
else :
print('\n Gravity deposition finished, apply fluid forces !\n')
newtonIntegr.damping = 0.0 # Set the artificial numerical damping to zero
gravDepo.dead = True # Remove the present engine for the following
gravDepo.dead = True # Remove the present engine for the following
hydroEngine.dead = False # Activate the HydroForceEngine
hydroEngine.vxFluid = vxFluidPY # Send the fluid velocity vector used to apply the drag fluid force on particles in HydroForceEngine (see c++ code)
hydroEngine.ReynoldStresses = np.zeros(ndimz) # Send the simplified fluid Reynolds stresses Rxz/\rho^f used to account for the fluid velocity fluctuations in HydroForceEngine (see c++ code)
......
......@@ -131,9 +131,9 @@ O.engines = [
InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Wall_Aabb(),Bo1_Facet_Aabb(),Bo1_Box_Aabb()],label='contactDetection',allowBiggerThanPeriod = True),
# Calculate the different interactions
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(), Ig2_Box_Sphere_ScGeom()],
[Ip2_ViscElMat_ViscElMat_ViscElPhys()],
[Law2_ScGeom_ViscElPhys_Basic()]
[Ig2_Sphere_Sphere_ScGeom(), Ig2_Box_Sphere_ScGeom()],
[Ip2_ViscElMat_ViscElMat_ViscElPhys()],
[Law2_ScGeom_ViscElPhys_Basic()]
,label = 'interactionLoop'),
#Apply an hydrodynamic force to the particles
HydroForceEngine(densFluid = densFluidPY,viscoDyn = kinematicViscoFluid*densFluidPY,zRef = groundPosition,gravity = gravityVector,deltaZ = dz,expoRZ = expoDrag_PY,lift = False,nCell = ndimz,vCell = length*width*dz,radiusPart=diameterPart/2.,vxFluid = np.array(vxFluidPY),phiPart = phiPartPY,vxPart = vxPartPY,ids = idApplyForce, label = 'hydroEngine', dead = True),
......@@ -167,7 +167,7 @@ def gravityDeposition(lim):
else :
print('\n Gravity deposition finished, apply fluid forces !\n')
newtonIntegr.damping = 0.0 # Set the artificial numerical damping to zero
gravDepo.dead = True # Remove the present engine for the following
gravDepo.dead = True # Remove the present engine for the following
hydroEngine.dead = False # Activate the HydroForceEngine
hydroEngine.vxFluid = vxFluidPY # Send the fluid velocity vector used to apply the drag fluid force on particles in HydroForceEngine (see c++ code)
hydroEngine.ReynoldStresses = np.ones(ndimz)*1e-4 # Send the simplified fluid Reynolds stresses Rxz/\rho^f used to account for the fluid velocity fluctuations in HydroForceEngine (see c++ code)
......
......@@ -160,9 +160,9 @@ O.engines = [
InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Wall_Aabb(),Bo1_Facet_Aabb(),Bo1_Box_Aabb()],label='contactDetection',allowBiggerThanPeriod = True),
# Calculate the different interactions
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(), Ig2_Box_Sphere_ScGeom()],
[Ip2_ViscElMat_ViscElMat_ViscElPhys()],
[Law2_ScGeom_ViscElPhys_Basic()]
[Ig2_Sphere_Sphere_ScGeom(), Ig2_Box_Sphere_ScGeom()],
[Ip2_ViscElMat_ViscElMat_ViscElPhys()],
[Law2_ScGeom_ViscElPhys_Basic()]
,label = 'interactionLoop'),
#Apply an hydrodynamic force to the particles
HydroForceEngine(densFluid = densFluidPY,viscoDyn = kinematicViscoFluid*densFluidPY,zRef = groundPosition,gravity = gravityVector,deltaZ = dz,expoRZ = expoDrag_PY,lift = False,nCell = ndimz,vCell = length*width*dz,radiusPart=diameterPart/2.,vxFluid = np.array(vxFluidPY),phiPart = phiPartPY,vxPart = vxPartPY,ids = idApplyForce, label = 'hydroEngine', dead = True,fluidWallFriction=True,channelWidth=width,phiMax = phiPartMax,iturbu = 1,ilm=2,iusl=1,irheolf=0),
......@@ -198,7 +198,7 @@ def gravityDeposition(lim):
else :
print('\n Gravity deposition finished, apply fluid forces !\n')
newtonIntegr.damping = 0.0 # Set the artificial numerical damping to zero
gravDepo.dead = True # Remove the present engine for the following
gravDepo.dead = True # Remove the present engine for the following
hydroEngine.dead = False # Activate the HydroForceEngine
hydroEngine.vxFluid = vxFluidPY # Send the fluid velocity vector used to apply the drag fluid force on particles in HydroForceEngine (see c++ code)
hydroEngine.ReynoldStresses = np.ones(ndimz)*1e-4 # Send the simplified fluid Reynolds stresses Rxz/\rho^f used to account for the fluid velocity fluctuations in HydroForceEngine (see c++ code)
......
......@@ -134,9 +134,9 @@ O.engines = [
InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Wall_Aabb(),Bo1_Facet_Aabb(),Bo1_Box_Aabb()],label='contactDetection',allowBiggerThanPeriod = True),
# Calculate the different interactions
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(), Ig2_Box_Sphere_ScGeom()],
[Ip2_ViscElMat_ViscElMat_ViscElPhys()],
[Law2_ScGeom_ViscElPhys_Basic()]
[Ig2_Sphere_Sphere_ScGeom(), Ig2_Box_Sphere_ScGeom()],
[Ip2_ViscElMat_ViscElMat_ViscElPhys()],
[Law2_ScGeom_ViscElPhys_Basic()]
,label = 'interactionLoop'),
#Apply an hydrodynamic force to the particles
HydroForceEngine(densFluid = densFluidPY,viscoDyn = kinematicViscoFluid*densFluidPY,zRef = groundPosition,gravity = gravityVector,deltaZ = dz,expoRZ = expoDrag_PY,lift = False,nCell = ndimz,vCell = length*width*dz,radiusPart=diameterPart/2.,vxFluid = np.array(vxFluidPY),phiPart = phiPartPY,vxPart = vxPartPY,ids = idApplyForce, label = 'hydroEngine', dead = True),
......@@ -170,7 +170,7 @@ def gravityDeposition(lim):
else :
print('\n Gravity deposition finished, apply fluid forces !\n')
newtonIntegr.damping = 0.0 # Set the artificial numerical damping to zero
gravDepo.dead = True # Remove the present engine for the following
gravDepo.dead = True # Remove the present engine for the following
hydroEngine.dead = False # Activate the HydroForceEngine
hydroEngine.vxFluid = vxFluidPY # Send the fluid velocity vector used to apply the drag fluid force on particles in HydroForceEngine (see c++ code)
hydroEngine.ReynoldStresses = np.ones(ndimz)*1e-4 # Send the simplified fluid Reynolds stresses Rxz/\rho^f used to account for the fluid velocity fluctuations in HydroForceEngine (see c++ code)
......
......@@ -6,8 +6,8 @@ from yade import pack,qt
O.periodic=True
O.cell.hSize=Matrix3(1.0, -0.15, -0.10,
-0.2 ,1.5, 0.3,
0.3, -0.3, 1.0)
-0.2 ,1.5, 0.3,
0.3, -0.3, 1.0)
sp=pack.SpherePack()
num=sp.makeCloud(hSize=O.cell.hSize, rMean=-0.01,rRelFuzz=.2, num=500,periodic=True, porosity=0.52,distributeMass=False)
O.bodies.append([sphere(s[0],s[1]) for s in sp])
......
......@@ -6,9 +6,9 @@ from yade import pack,qt
O.periodic=True
O.cell.hSize=Matrix3(0.1, 0, 0,
0 ,0.1, 0,
0, 0, 0.1)
0 ,0.1, 0,
0, 0, 0.1)
sp=pack.SpherePack()
radius=5e-3
num=sp.makeCloud(Vector3().Zero,O.cell.refSize,radius,.2,500,periodic=True) # min,max,radius,rRelFuzz,spheresInCell,periodic
......
......@@ -15,8 +15,8 @@ width=0.2
thickness=0.01
O.cell.hSize=Matrix3(length, 0, 0,
0 ,3.*height, 0,
0, 0, width)
0 ,3.*height, 0,
0, 0, width)
O.materials.append(FrictMat(density=1,young=1e5,poisson=0.3,frictionAngle=radians(30),label='boxMat'))
lowBox = box( center=(length/2.0,height-thickness/2.0,width/2.0), extents=(length*1000.0,thickness/2.0,width*1000.0) ,fixed=True,wire=False)
......
......@@ -71,8 +71,8 @@ for s in sp:
highlight=False
b.shape=PotentialBlock(k=0.2, r=0.05*meanSize, R=1.02*sphereRad, a=[1.0,-1.0,0.0,0.0,0.0,0.0], b=[0.0,0.0,1.0,-1.0,0.0,0.0], c=[0.0,0.0,0.0,0.0,1.0,-1.0], d=[distanceToCentre-rPP,distanceToCentre-rPP,distanceToCentre-rPP,distanceToCentre-rPP,distanceToCentre-rPP,distanceToCentre-rPP],isBoundary=False,color=color,wire=wire,highlight=highlight,minAabb=Vector3(3.0*sphereRad,3.0*sphereRad,3.0*sphereRad),maxAabb=Vector3(3.0*sphereRad,3.0*sphereRad,3.0*sphereRad),maxAabbRotated=Vector3(3.0*sphereRad,3.0*sphereRad,3.0*sphereRad),minAabbRotated=Vector3(3.0*sphereRad,3.0*sphereRad,3.0*sphereRad),AabbMinMax=True,fixedNormal=False,id=count)
length=meanSize
V= 1.0
geomInert=(2./5.)*powderDensity*V*sphereRad**2
V= 1.0
geomInert=(2./5.)*powderDensity*V*sphereRad**2
utils._commonBodySetup(b,V,Vector3(geomInert,geomInert,geomInert), material='frictionless',pos=s[0], dynamic=True, fixed=False)
b.state.pos = s[0] #s[0] stores center
b.state.ori = Quaternion((random.random(),random.random(),random.random()),random.random()) #s[2]
......
......@@ -53,8 +53,8 @@ for s in sp:
highlight=False
b.shape=PotentialParticle(k=0.2, r=0.05*meanSize, R=1.02*sphereRad, a=[1.0,-1.0,0.0,0.0,0.0,0.0], b=[0.0,0.0,1.0,-1.0,0.0,0.0], c=[0.0,0.0,0.0,0.0,1.0,-1.0], d=[distanceToCentre-rPP,distanceToCentre-rPP,distanceToCentre-rPP,distanceToCentre-rPP,distanceToCentre-rPP,distanceToCentre-rPP],isBoundary=False,color=color,wire=wire,highlight=highlight,minAabb=Vector3(1.2*sphereRad,1.2*sphereRad,1.2*sphereRad),maxAabb=Vector3(1.2*sphereRad,1.2*sphereRad,1.2*sphereRad),maxAabbRotated=Vector3(1.2*sphereRad,1.2*sphereRad,1.2*sphereRad),minAabbRotated=Vector3(1.2*sphereRad,1.2*sphereRad,1.2**sphereRad),AabbMinMax=True, id=count)
length=meanSize
V= 1.0
geomInert=(2./5.)*powderDensity*V*sphereRad**2
V= 1.0
geomInert=(2./5.)*powderDensity*V*sphereRad**2
utils._commonBodySetup(b,V,Vector3(geomInert,geomInert,geomInert), material='frictionless',pos=[0.0,0,0], dynamic=True, fixed=False)
b.state.pos = s[0] #s[0] stores center
b.state.ori = Quaternion((random.random(),random.random(),random.random()),random.random()) #s[2]
......
......@@ -55,9 +55,9 @@ yade.qt.View();
from yade import plot
plot.plots={'t':('pos1',None,'vel1')}
def history():
plot.addData(pos1=O.bodies[0].state.pos[1], # potential elastic energy
vel1=O.bodies[0].state.vel[1],
t=O.time)
plot.addData(pos1=O.bodies[0].state.pos[1], # potential elastic energy
vel1=O.bodies[0].state.vel[1],
t=O.time)
#yade.qt.Renderer().bound=True
plot.plot(subPlots=False)
......
......@@ -39,17 +39,17 @@ O.materials.append( FrictMat( young=E,poisson=0.3,density=1000,frictionAngle=rad
nodesIds=[]
cylIds=[]
cylinder((0,0,0),(L,0,0),radius=r,nodesIds=nodesIds,cylIds=cylIds,
fixed=True,color=[1,0,0],intMaterial='cMat',extMaterial='fMat')
fixed=True,color=[1,0,0],intMaterial='cMat',extMaterial='fMat')
cylinder((L/4,2*L/3,L),(L/4,-L/3,L),radius=r,nodesIds=nodesIds,cylIds=cylIds,
fixed=False,color=[0,1,0],intMaterial='cMat',extMaterial='fMat')
fixed=False,color=[0,1,0],intMaterial='cMat',extMaterial='fMat')
cylinder((0,2*L/3,L),(0,-L/3,L),radius=r,nodesIds=nodesIds,cylIds=cylIds,
fixed=False,color=[0,1,0],intMaterial='cMat',extMaterial='fMat')
fixed=False,color=[0,1,0],intMaterial='cMat',extMaterial='fMat')
cylinder((L/2,L/2,L),(L/2,-L/2,L),radius=r,nodesIds=nodesIds,cylIds=cylIds,
fixed=False,color=[0,1,1],intMaterial='cMat',extMaterial='fMat')
fixed=False,color=[0,1,1],intMaterial='cMat',extMaterial='fMat')
cylinder((3*L/4,L/3,L),(3*L/4,-2*L/3,L),radius=r,nodesIds=nodesIds,cylIds=cylIds,
fixed=False,color=[0,0,1],intMaterial='cMat',extMaterial='fMat')
fixed=False,color=[0,0,1],intMaterial='cMat',extMaterial='fMat')
cylinder((L,L/3,L),(L,-2*L/3,L),radius=r,nodesIds=nodesIds,cylIds=cylIds,
fixed=False,color=[0,0,1],intMaterial='cMat',extMaterial='fMat')
fixed=False,color=[0,0,1],intMaterial='cMat',extMaterial='fMat')
#### For viewing ####
from yade import qt
......
......@@ -58,7 +58,7 @@ for i in range(0,n):
L=Lmin+dL*random.random()
color=[random.random(),random.random(),random.random()]
cylinder((x,y,2*r),(x,y,L+2*r),radius=r,nodesIds=nodesIds,cylIds=cylIds,
fixed=False,color=color,intMaterial='cMat',extMaterial='fMat')
fixed=False,color=color,intMaterial='cMat',extMaterial='fMat')
#### Creat ground with pFacets ####
color=[255./255.,102./255.,0./255.]
......
......@@ -76,14 +76,14 @@ fixed_tail_body_upperlimit=1.01;
O.dt=1e-9;
O.engines=[
ForceResetter(),
## Apply internal force to the deformable elements and internal force of the interaction element
FEInternalForceEngine([If2_Lin4NodeTetra_LinIsoRayleighDampElast(),If2_2xLin4NodeTetra_LinCohesiveStiffPropDampElastMat()]),
NewtonIntegrator(damping=0,gravity=[0,0,0]),
#PyRunner(virtPeriod=1e-99,command='applyforcetoelements()'),
# ## Plotting data: adds plots after one step of the integrator engine
#PyRunner(virtPeriod=1e-99,command='addplot()')
]
ForceResetter(),
## Apply internal force to the deformable elements and internal force of the interaction element
FEInternalForceEngine([If2_Lin4NodeTetra_LinIsoRayleighDampElast(),If2_2xLin4NodeTetra_LinCohesiveStiffPropDampElastMat()]),
NewtonIntegrator(damping=0,gravity=[0,0,0]),
#PyRunner(virtPeriod=1e-99,command='applyforcetoelements()'),
# ## Plotting data: adds plots after one step of the integrator engine
#PyRunner(virtPeriod=1e-99,command='addplot()')
]
from yade import plot
......
......@@ -152,8 +152,8 @@ fixboundryelements();
def addplot():
vecpos=(O.bodies[nodebodies1[0].id].state.pos-initpos0)
plot.addData(displacement=abs(vecpos[2]),time=O.time)
vecpos=(O.bodies[nodebodies1[0].id].state.pos-initpos0)
plot.addData(displacement=abs(vecpos[2]),time=O.time)
......@@ -178,13 +178,11 @@ def addplot():
##integratoreng.abs_err=1e-3;
O.engines=[ForceResetter(),
PyRunner(virtPeriod=1e-99,command='applyforcetoelements()'),
## Apply internal force to the deformable elements and internal force of the interaction element
FEInternalForceEngine([If2_Lin4NodeTetra_LinIsoRayleighDampElast(),If2_2xLin4NodeTetra_LinCohesiveStiffPropDampElastMat()]),
NewtonIntegrator(damping=0,gravity=[0,0,0]),
PyRunner(virtPeriod=1e-99,command='addplot()'),
PyRunner(virtPeriod=1e-99,command='applyforcetoelements()'),
## Apply internal force to the deformable elements and internal force of the interaction element
FEInternalForceEngine([If2_Lin4NodeTetra_LinIsoRayleighDampElast(),If2_2xLin4NodeTetra_LinCohesiveStiffPropDampElastMat()]),
NewtonIntegrator(damping=0,gravity=[0,0,0]),
PyRunner(virtPeriod=1e-99,command='addplot()'),
];
......
......@@ -98,7 +98,7 @@ def applyforcetoelements():
# plot.addData(force=O.forces.f(i));
def addplot():
plot.addData(force=O.forces.f(forcebodies[1])[0],pos=(O.bodies[forcebodies[1]].state.pos[0]-initialpositions[forcebodies[1]][0]),vel=O.bodies[forcebodies[1]].state.vel[0],t=O.time,time=O.time,tm=O.time)
plot.addData(force=O.forces.f(forcebodies[1])[0],pos=(O.bodies[forcebodies[1]].state.pos[0]-initialpositions[forcebodies[1]][0]),vel=O.bodies[forcebodies[1]].state.vel[0],t=O.time,time=O.time,tm=O.time)
for i in forcebodies:
O.bodies[i].state.pos+=Vector3(0,0,0.05);
......@@ -125,13 +125,13 @@ for i in forcebodies:
O.dt=1e-8;
O.engines=[
ForceResetter(),
## Apply internal force to the deformable elements and internal force of the interaction element
FEInternalForceEngine([If2_Lin4NodeTetra_LinIsoRayleighDampElast(),If2_2xLin4NodeTetra_LinCohesiveStiffPropDampElastMat()]),
PyRunner(iterPeriod=1,command='applyforcetoelements()'),
NewtonIntegrator(damping=0,gravity=[0,0,0]),
ForceResetter(),
## Apply internal force to the deformable elements and internal force of the interaction element
FEInternalForceEngine([If2_Lin4NodeTetra_LinIsoRayleighDampElast(),If2_2xLin4NodeTetra_LinCohesiveStiffPropDampElastMat()]),
PyRunner(iterPeriod=1,command='applyforcetoelements()'),
NewtonIntegrator(damping=0,gravity=[0,0,0]),
# ## Plotting data: adds plots after one step of the integrator engine
PyRunner(iterPeriod=1,command='addplot()')
PyRunner(iterPeriod=1,command='addplot()')
];
......
......@@ -177,14 +177,14 @@ O.dt=1e-3;
O.engines=[
# integratoreng,
ForceResetter(),
# ## Apply internal force to the deformable elements and internal force of the interaction element
FEInternalForceEngine([If2_Lin4NodeTetra_LinIsoRayleighDampElast(),If2_2xLin4NodeTetra_LinCohesiveStiffPropDampElastMat()]),
PyRunner(iterPeriod=1,command='applyforcetoelements()'),
NewtonIntegrator(damping=0,gravity=[0,0,0]),
# ## Plotting data: adds plots after one step of the integrator engine
PyRunner(iterPeriod=1,command='addplot()')
]
ForceResetter(),
# Apply internal force to the deformable elements and internal force of the interaction element
FEInternalForceEngine([If2_Lin4NodeTetra_LinIsoRayleighDampElast(),If2_2xLin4NodeTetra_LinCohesiveStiffPropDampElastMat()]),
PyRunner(iterPeriod=1,command='applyforcetoelements()'),
NewtonIntegrator(damping=0,gravity=[0,0,0]),
# Plotting data: adds plots after one step of the integrator engine
PyRunner(iterPeriod=1,command='addplot()')
]
from yade import plot
......
......@@ -75,9 +75,9 @@ yade.qt.View();
from yade import plot
plot.plots={'t':('pos1',None,'vel1')}
def history():
plot.addData(pos1=O.bodies[0].state.pos[1], # potential elastic energy
vel1=O.bodies[0].state.vel[1],
t=O.time)
plot.addData(pos1=O.bodies[0].state.pos[1], # potential elastic energy
vel1=O.bodies[0].state.vel[1],
t=O.time)
#yade.qt.Renderer().bound=True
plot.plot(subPlots=False)
......
......@@ -73,12 +73,12 @@ for o in O.bodies:
o.shape.color=(0.9,0.8,0.6)
## to fix boundary particles on ground
if o.state.pos[2]<(zinf+2*e) :
o.state.blockedDOFs+='xyz'
o.shape.color=(1,1,1)
o.state.blockedDOFs+='xyz'
o.shape.color=(1,1,1)
## to identify indicator on top
if o.state.pos[2]>(zsup-e) and o.state.pos[0]>(xsup-e) and o.state.pos[1]>((yinf+ysup-e)/2.0) and o.state.pos[1]<((yinf+ysup+e)/2) :
refPoint=o.id
refPoint=o.id
O.bodies[refPoint].shape.highlight=True
......@@ -112,16 +112,16 @@ stableIter=1000
stableVel=0.001
degrade=True
def jointStrengthDegradation():
global degrade
if degrade and O.iter>=stableIter and abs(O.bodies[refPoint].state.vel[2])<stableVel :
global degrade
if degrade and O.iter>=stableIter and abs(O.bodies[refPoint].state.vel[2])<stableVel :
print('Equilibrium reached \nJoint cohesion canceled now !', ' | iteration=', O.iter)
degrade=False
for i in O.interactions:
if i.phys.isOnJoint :
if i.phys.isOnJoint :
if i.phys.isCohesive:
i.phys.isCohesive=False
i.phys.FnMax=0.
i.phys.FsMax=0.
i.phys.isCohesive=False
i.phys.FnMax=0.
i.phys.FsMax=0.
print('Seeking after an initial equilibrium state')
print('')
......
......@@ -76,14 +76,14 @@ for o in O.bodies:
o.shape.color=(0.9,0.8,0.6)
## to fix boundary particles on ground
if o.state.pos[1]<(yinf+2*e) :
o.state.blockedDOFs+='xyz'
baseBodies.append(o.id)
o.shape.color=(1,1,1)
o.state.blockedDOFs+='xyz'
baseBodies.append(o.id)
o.shape.color=(1,1,1)
## to identify indicator on top
if o.state.pos[1]>(ysup-e) and o.state.pos[0]>(xsup-e) and o.state.pos[2]>(zinf+(Z-e)/2) and o.state.pos[2]<(zsup-(Z-e)/2) :
refPoint=o.id
p0=o.state.pos[1]
refPoint=o.id
p0=o.state.pos[1]
baseBodies=tuple(baseBodies)
O.bodies[refPoint].shape.color=(1,0,0)
......@@ -124,12 +124,12 @@ def jointStrengthDegradation():
print('!joint cohesion total degradation!', ' | iteration=', O.iter)
degrade=False
for i in O.interactions:
if i.phys.isOnJoint :
if i.phys.isOnJoint :
if i.phys.isCohesive:
i.phys.isCohesive=False
i.phys.FnMax=0.
i.phys.FsMax=0.
i.phys.isCohesive=False
i.phys.FnMax=0.
i.phys.FsMax=0.
#### YADE windows
from yade import qt
v=qt.Controller()
......