yade-users team mailing list archive

[Mahdeyeh <question686970@xxxxxxxxxxxxxxxxxxxxx>]

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

    Status: Answered => Open

Mahdeyeh is still having a problem:
Than you Jan
I think so.

I ran the codes with :  yade j-8 code.py    , but it shows this error:

"ERROR /build/yade-
KKgSmd/yade-1.20.0/pkg/common/InsertionSortCollider.cpp:65
insertionSortParallel: Parallel insertion sort needs verletDist>0"

I think it because of my particles are polyhedra not sphere, yes? so is
this error is important? How I can solve it?

Here is my code:


from yade import export,polyhedra_utils
import os
from yade import plot
import math
from yade import utils
import pylab
import matplotlib; matplotlib.rc('axes',grid=True)
from matplotlib import pyplot
import numpy as np
from numpy import *
from yade import export as expt


# Input data:

RawVer=np.genfromtxt('ring.txt',names=True,dtype=None)
# ListVer is list of all the vertices of ring`s polygons
Ver=()
ListVer=[]
for b in RawVer:
    if b[0]=='*LWPOLYLINE':
        ListVer.append(Ver)
        Ver=()
        continue
    Cordn=b[0]
    Cordn=np.fromstring(Cordn, sep=',')
    Cordn=tuple(Cordn.tolist())
    Cordn1=Cordn+(0.1,)   # add z vertex to coordinates 
    Cordn2=Cordn+(0.2,)   # add z vertex to coordinates
    if not Cordn1 in Ver:
        Ver=Ver+(Cordn1,Cordn2)
ListVer.append(Ver)

RawVer1=np.genfromtxt('boundary.txt',names=True,dtype=None)
# ListVer1 is list of all the vertices of boundary`s polygons
Ver1=()
ListVer1=[]
for b in RawVer1:
    if b[0]=='*LWPOLYLINE':
        ListVer1.append(Ver1)
        Ver1=()
        continue
    Cordn=b[0]
    Cordn=np.fromstring(Cordn, sep=',')
    Cordn=tuple(Cordn.tolist())
    Cordn1=Cordn+(-2,)   # add z vertex to coordinates 
    Cordn2=Cordn+(2,)   # add z vertex to coordinates
    if not Cordn1 in Ver1:
        Ver1=Ver1+(Cordn1,Cordn2)
ListVer1.append(Ver1)

RawVer2=np.genfromtxt('waste.txt',names=True,dtype=None)
# ListVer2 is list of all the vertices of caved waste rock`s polygons
Ver2=()
ListVer2=[]
for b in RawVer2:
    if b[0]=='*LWPOLYLINE':
        ListVer2.append(Ver2)
        Ver2=()
        continue
    Cordn=b[0]
    Cordn=np.fromstring(Cordn, sep=',')
    Cordn=tuple(Cordn.tolist())
    Cordn1=Cordn+(0.1,)   # add z vertex to coordinates 
    Cordn2=Cordn+(0.2,)   # add z vertex to coordinates
    if not Cordn1 in Ver2:
        Ver2=Ver2+(Cordn1,Cordn2)
ListVer2.append(Ver2)


# Materials type:

Dolomite = PolyhedraMat()
Dolomite.density = 2870 #kg/m^3 
Dolomite.young = 24.36e9 #Pa 
Dolomite.poisson = 0.2
Dolomite.frictionAngle = radians(55.12) #rad

Shale = PolyhedraMat()
Shale.density = 2750 #kg/m^3 
Shale.young = 6e9 #Pa 
Shale.poisson = 0.23
Shale.frictionAngle = radians(42) #rad

for ii in ListVer:
    O.bodies.append(polyhedra_utils.polyhedra(Dolomite,v=ii,fixed=False, color=(0.1,0.5,0.2), mask=3))

for iii in ListVer1:
    O.bodies.append(polyhedra_utils.polyhedra(Dolomite,v=iii,fixed=True, color=(1,0,1), mask=4))

for iiii in ListVer2:
    O.bodies.append(polyhedra_utils.polyhedra(Shale,v=iiii,fixed=False, color=(0.9,0.81,0.45), mask=5))


O.bodies.erase(340)   # delete wall under ring: id: 340

# block z translation and block x, y, z rotations
for b in O.bodies:
    if b.mask is 3:
        b.state.blockedDOFs='zXY'
for b in O.bodies:
    if b.mask is 5:
        b.state.blockedDOFs='zXY'

# function for calming particles
def calm():
    for c in O.bodies:
        c.state.vel=Vector3(0,0,0)    
        c.state.angVel=Vector3(0,0,0)

# returns a value that can be useful for evaluating the stability of the packing. It is defined as (mean force on particles)/(mean contact force), so that it tends to 0 in a stable packing.
def checkUnbalanced():   
    print "iter %d , unbalanced forces %f"%(O.iter, utils.unbalancedForce())  # %[(keyname)][flags][width][.precision]typecode : String Formatting		
    iter00=O.iter
        
    Unbalanced=open("Unbalanced iter Unbalanced forces.txt","a")
    Unbalanced.write(repr(iter00)+'	'+repr(utils.unbalancedForce())+'	'+"\n")
    Unbalanced.close()

# Engines:

#os.mkdir(O.tags['id'])   # Created separate directory and put all files
of VTK in it, till external files be seperate because of large
simulation.


O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Polyhedra_Aabb(),]), # We can set collider's verletDist to a fraction of the polyhedra minimum size, since it determines how much is each body enlarged to avoid collision detection at every step.
   InteractionLoop(
      [Ig2_Polyhedra_Polyhedra_PolyhedraGeom(),], 
      [Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys()], # collision "physics"
      [Law2_PolyhedraGeom_PolyhedraPhys_Volumetric()]   # contact law -- apply forces
   ),
   NewtonIntegrator(gravity=(0,-9.81,0),damping=0.2),
]

utils.calm()

# the model has to calm, because there are some overlaps in particles. 
O.engines=O.engines+[PyRunner(iterPeriod=20,command='calm()',label="calmRunner")]  #  because we need to calm only on the first few steps in our model.
O.engines=O.engines+[PyRunner(command='checkUnbalanced()',iterPeriod=10000,label="checker")]  # call our function defined above 1000 cycles


O.dt=10e-6


O.saveTmp('model1')
O.save('model1.bz2')


# first run of model
O.run(50000,True)

O.save('model2.bz2')


# Outputs:


clrOre=[0.1,0.5,0.2]
clrWaste=[0.9,0.81,0.45]


def positions():
    for b in O.bodies:
        if b.shape.color==clrOre:
            time1=O.time
        
            Positions=open("Positions time x y z.txt","a")
            Positions.write(repr(time1)+'	'+repr(b.state.pos[0])+'	'+repr(b.state.pos[1])+'	'+repr(b.state.pos[2])+'	'"\n")
            Positions.close() 



def oreAmount():
    clrOre=[0.1,0.5,0.2]
    M_o=0
    for b in O.bodies:
        if b.shape.color==clrOre:
                M_o+=b.state.mass
    return M_o


os.mkdir(O.tags['id'])   # Created separate directory and put all files of VTK in it, till external files be seperate because of large simulation.
O.engines=O.engines+[PyRunner(iterPeriod=20000,command='VTKview()',label="VTKview")]          # every 20000 cycles we can see picture of model in Paraview 
PolVtkData = expt.VTKExporter(O.tags['id']+'/'+'polData')            # save animation and images of simulation in VTK format (paraview) in defined cycle
def VTKview():
    PolVtkData.exportPolyhedra()


m_o  =  0.00005         # mass of ore from the block
m_w  =  0.00005         # mass of waste

R    =  0              # R total ore utilization

D    =  0
D_m  =  0

E    =  0             # E extraction

M_o  =  oreAmount()               # M_o total ore mass


totalEMass = 1.25 * M_o         # Convergence criteria for the models are set as 125% extraction of the ring mass.

def excavation():
        m_o  =  0.00005         
        m_w  =  0.00005
        for b in O.bodies:
            if b.state.pos[0]>(-2.4) and b.state.pos[0]<1 and b.state.pos[1]>0 and b.state.pos[1]<1:     # -2.4<x<1 because we don`t know about repose angle and distance for bucket yet 
                
                if b.shape.color==clrOre:
                    m_o  +=  b.state.mass      #total ore mass extracted
                    print 'ore' 
                if b.shape.color==clrWaste:
                    m_w  +=  b.state.mass    #total waste mass
                    print 'waste'

        print 'm_o' , m_o , 'm_w' , m_w
        iter1=O.iter
        
        Mass=open("Mass iter massore masswaste.txt","a")
        Mass.write(repr(iter1)+'	'+repr(m_o)+'	'+repr(m_w)+'	'"\n")
        Mass.close() 


        E       = ( m_o + m_w) / M_o * 100              ##extraction definition

        D       = (m_w / (m_w + m_o)) * 100       ##ore dilution

        D_m     = (1- ( m_o / (m_w + m_o))) * 100   ##metal dilution

        R       = (m_o / M_o) * 100                             ## Total ore recovery
        print 'E' , E , 'D' , D , 'D_m' , D_m , 'R' , R
        iter2=O.iter
        
        DataAnalysis=open("DataAnalysis iter E D DM R.txt","a")
        DataAnalysis.write(repr(iter2)+'	'+repr(E)+'	'+repr(D)+'	'+repr(D_m)+ '	'+repr(R)+"\n")
        DataAnalysis.close()


# Extra Engines:

O.engines=O.engines+[PyRunner(command='positions()',realPeriod=0.5,label="flow")]
# call our function defined above every 60 seconds

O.engines=O.engines+[PyRunner(iterPeriod=200,command='excavation()',label="excavate")]
# call our function defined above every 200 cycles


polyhedra_utils.SizeRatio()    # save sizes of polyhedra into file
export.textPolyhedra('poly.txt')   # Save polyhedra into a text file, id, number of vertices and number of surfaces, particleVelocity and ... .

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