yade-users team mailing list archive

[weijie <question692957@xxxxxxxxxxxxxxxxxxxxx>]

Question #692957 on Yade changed:
https://answers.launchpad.net/yade/+question/692957

    Status: Needs information => Open

weijie gave more information on the question:
Hi Jan,

I simplified my script and deleted many unnecessary parts, leaving only
a few polyhedra. But I found that the more polyhedrons, the more serious
the problem.

Best regards
Jie

Below is my script:
########
from yade import polyhedra_utils,pack,plot,utils,export,qt
import numpy as np
import math
import random
#-------------------------------------------
# Engines
Kn=1e8
Ks=Kn/10
O.engines=[
 ForceResetter(),
 InsertionSortCollider([PotentialBlock2AABB()],verletDist=0.005, avoidSelfInteractionMask=2),
 InteractionLoop(
  [Ig2_PB_PB_ScGeom(calContactArea=True)],
  [Ip2_FrictMat_FrictMat_KnKsPBPhys(kn_i=Kn, ks_i=Ks, Knormal=Kn, Kshear=Ks, viscousDamping=0.2)],
  [Law2_SCG_KnKsPBPhys_KnKsPBLaw(allowViscousAttraction=True, traceEnergy=False)]
 ),
 NewtonIntegrator(damping=0.0,exactAsphericalRot=True,gravity=[0,-9.81,0]),
]
O.materials.append(FrictMat(young=-1,poisson=-1,frictionAngle=radians(40.0),density=2500,label='frictionless'))
#-------------------------------------------
# Dimensions
lengthOfBase = 0.250
heightOfBase = 0.600
#-------------------------------------------
# Make cloud
sp=pack.SpherePack()
mn,mx=Vector3(-0.5*lengthOfBase,0.4*heightOfBase ,-0.5*lengthOfBase),Vector3(0.5*lengthOfBase,0.6*heightOfBase,0.5*lengthOfBase)
sp.makeCloud(mn,mx,psdSizes=[0.0035,0.005,0.01,0.02,0.04],psdCumm=(0,0.16,0.32,0.84,1),num=-1,distributeMass=False)
#Generate polyhedrons of different sizes
for center,radius in sp:   
    if radius<0.0025 :
        r=0.01*0.0025
        b=Body()
        b.mask=1
        b.aspherical=True
        color=[0.3,0.66,0.78]
        distanceToCentre=0.00175
        b.shape=PotentialBlock(k=0.0, r=r, R=0.0,a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1],d=[distanceToCentre-r,distanceToCentre-r,distanceToCentre-r,distanceToCentre-r,distanceToCentre-r,distanceToCentre-r], isBoundary=False, color=color, wire=False, AabbMinMax=True, fixedNormal=False)
        utils._commonBodySetup(b, b.shape.volume, b.shape.inertia, material='frictionless',pos=center,fixed=False)
        b.state.pos =center #s[0] stores center
        b.state.ori = Quaternion((random.random(),random.random(),random.random()),random.random()) #s[2]
        O.bodies.append(b)       
    if radius>0.0025 and radius<0.005 :
        r=0.01*0.0025
        b=Body()
        b.mask=1
        b.aspherical=True
        color=[0.1,0.2,0.3]
        distanceToCentre=0.00375
        b.shape=PotentialBlock(k=0.0, r=r, R=0.0,a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1],d=[distanceToCentre-r,distanceToCentre-r,distanceToCentre-r,distanceToCentre-r,distanceToCentre-r,distanceToCentre-r], isBoundary=False, color=color, wire=False, AabbMinMax=True, fixedNormal=False)
        utils._commonBodySetup(b, b.shape.volume, b.shape.inertia, material='frictionless',pos=center,fixed=False)
        b.state.pos =center #s[0] stores center
        b.state.ori = Quaternion((random.random(),random.random(),random.random()),random.random()) 
        O.bodies.append(b)
    if radius>0.01 and radius<0.02 :
        aa= [0.08376984881458294, -0.47781815584193316, -0.6272525512350887, 0.619009994702633, 0.03464399632141613, -0.9832172384126057, -0.2783328604495522, 0.764868485410703, -0.07639578575476973, -0.5670674455594352, 0.621538863800698, -0.7670822373612994, 0.766206355713753, -0.579033908027536, 0.14707077495765938]
        bb= [-0.8413821074662718, -0.8774642832791362, -0.34697163535772957, 0.10830905836646608, 0.9770674225401074, -0.13517146281570994, 0.6998398693326321, 0.6393484283187056, -0.7572080749137317, -0.19334270015810612, -0.5555042677443186, 0.6224892308049091, -0.516288957458027, 0.5202427575278886, 0.25146364935432863]
        cc= [0.533908945106932, -0.04178805471957184, -0.6972552769440008, 0.7778790229425239, -0.21009294450245802, 0.12252566151035581, 0.6578411480642227, 0.07880220321753764, -0.6486906930150549, 0.8006579247607992, 0.5523626067206088, 0.15521597423169936, -0.3825879413556488, -0.6277477252799387, -0.956627524278262]
        dd= [0.012044189,0.008907045,0.012620139,0.012629451,0.012161151,0.005857936,0.00904349,0.007525031,0.012808652,0.014389324,0.006000359,0.01150185,	0.009394797,0.0084963,0.008248162]
        r=min(dd)/2 #Suggested value                                                                                                                                                                   
        b=Body()
        b.mask=1
        b.aspherical=True
        color=[0.5,0.4,0.3]
        b.shape=PotentialBlock(k=0.0, r=r, R=0.0, a=aa, b=bb, c=cc, d=dd, isBoundary=False, color=color, wire=False, AabbMinMax=True, fixedNormal=False)
        utils._commonBodySetup(b, b.shape.volume, b.shape.inertia, material='frictionless',pos=center,fixed=False)
        b.state.pos =center #s[0] stores center
        b.state.ori = Quaternion((random.random(),random.random(),random.random()),random.random()) 
        O.bodies.append(b)   
#------------------------------------------
# Bottom plate
color=[0,0.5,1]
wallThickness=0.025
r=0.15*wallThickness
bbb=Body()
bbb.mask=3
wire=False
bbb.shape=PotentialBlock(k=0.0, r=r, R=0.0, a=[1,-1,0,0,0,0], b=[0,0,1,-1,0,0], c=[0,0,0,0,1,-1], d=[2.5*lengthOfBase-r,2.5*lengthOfBase-r,0.5*wallThickness-r,0.5*wallThickness-r,lengthOfBase/2.0-r,lengthOfBase/2.0-r], isBoundary=False, color=color, wire=wire, AabbMinMax=True, fixedNormal=False)
utils._commonBodySetup(bbb, bbb.shape.volume, bbb.shape.inertia, material='frictionless', pos=[0,0,0], fixed=True)
bbb.state.pos = [1*lengthOfBase,0,0]
O.bodies.append(bbb)
#-------------------------------------------
# Plot results
def myAddPlotData():
    KE = utils.kineticEnergy()
    plot.addData(timeStep1=O.iter,time=O.time,kineticEn=KE)
from yade import plot
plot.plots={'timeStep1':('kineticEn'),}
plot.plot() 
O.engines=O.engines+[PyRunner(iterPeriod=20,command='myAddPlotData()')]
#Timestep
O.dt = 1e-5
#-------------------------------------------
from yade import qt
v=qt.View()
v.sceneRadius=3.0
O.saveTmp()

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