clumps generating randomly in 2D simulation

Hi,

> The first one is that making a generation using PFC2D and then importing the
> information.txt to YADE. But I failed to get a correct simulation in this way and
> gave the reason yesterday.

solution for this is here:

https://answers.launchpad.net/yade/+question/229898

For the second option I think actually there is no tool in yade, that can do "random generation of clumps in 2D".
BTW, the clump sizes should be equal and positions randomly, right?
I think you have to write your own python code, that creates your packing for 2D.

Hope it helps,

christian

Dear Christian Jakob:
     Thanks a lot for your perfect suggestions yesterday. You are exactly right. The only purpose for me is to get clumps in 2D simulation with the equal sizes but the random positions.
     I think that task is much more complex than "spherepack" but feasible in YADE. Because they have something in common. In the first step work, the relationship between "unit cells" and "clumps" is the same as that between "spheres" and "spherepack". And the unit cells defination is just a very easy task in YADE.
     The second step work, which can be called clumps generation according to unit cells should be difficult for me for that I now don't know either the algorithm of "spherepack" or the c++ language. If now the only way to solve my problem is to write my own code for use, what can I do now? Seeking your help!

Dear Christian Jakob and Bruno Chareyre and all users:
     I have consulted some handbooks about the PFC-2D clumps generation and think that it is almost no use for me to employ PFC-2D for generating clumps. Because its generation pattern is totally different from YADE's. As I know, in PFC-2D the first step working for generating clumps is to firstly create all the balls in specific region. Then clumps should be determined according to the "spheres ids" as the second step working. Sometimes the determination condition is the desirable shape of the clumps, which can model many different types of particles.

I would do somethink like this:

- generate random points in your desired area with random.seed() [1]
- use this points as balance points for your clumps and generate them, use [2] or [3]
- after one clump is generated, check if there are overlaps with other particles -> if there is an overlap, do not generate and try another random point ... something like:

for b in O.bodies:
 if isinstance(b.shape,Sphere):
  if (abs(b.state.pos - randomClumpMemberPos) > (b.shape.radius + ClumpMemberRadius)):
   #create the sphere
  else:
   #try another random point

- then you will get a loose clump packing
- switch on gravity and let your packing settle down

[1] http://docs.python.org/2/library/random.html#random.seed
[2] https://www.yade-dem.org/doc/yade.wrapper.html?highlight=clump#yade.wrapper.BodyContainer.clump
[3] https://www.yade-dem.org/doc/yade.wrapper.html?highlight=clump#yade.wrapper.BodyContainer.appendClumped

hope it helps,

christian

Dear Christian Jakob:
     Thanks a lot for your perfect suggestions last time. I have done some test following your 3 steps working and of course there are some questions in this in this process, which can be described as that:
    (1).Firstly, I wanna employ the "Pseudo Random Generation Method"(can be called PRGM here) to define the coordinates of the balance points of the clumps and the workflow is like this:
     1. Do mesh for the generation area. The generation area is set up as a rectangle and can be averagely divided into N copies on the X direction as well as N copies on the Y direction. So the number of all points in this generation area which can be picked out as the balance points of clumps is N*N.
     2. Determination of this points' coordinations. As you know, the generation area in YADE can be defined as "mn,mx=Vector3(,,),Vector3(,,)"(Z-direction can be assigned as the same value in my 2D simulation) and the balance points' coordinations can only sucessfully determined by this algorithm "(mx-mn)/N*random.randint(0,N)"
    (2).Secondly, generate only spheres with this points only for testing the feasibility of this method. All need to do is just only to input the number of shperes desired and,then make a "O.bodies.append". All the scripts for these first two steps working are as follow:

import random

i=int(input("i=")) # mesh numbers in one direction
N=int(input("N=")) # sphere numbers desired

for i in range (1,i):
    mn,mx=Vector3(0,0,0),Vector3(0.07,0.07,0)
    p=(mx[0]-mn[0])/N*random.randint(0,N)
    q=(mx[1]-mn[1])/N*random.randint(0,N)
    O.bodies.append(sphere((p,q,0),radius=0.001))

    (3).There may be of course some overlaps occuring between the generated spheres and the next step working for me is nothing but detecting and removing the overlaps. But lack of experience in Python makes me only use the commands below:

    for b in O.bodies:
        if isinstance(b.shape,Sphere):
            if (abs(b.state.pos[0] - p) > radius):
                O.bodies.append(sphere((p,q,0),radius=radius))
            else:
                continue

There must be some errors in this scripts for that YADE's core dumped without any error noticed.
What is wrong with it and how to correct it? Seeking your help!(all the scripts are that)

import random

i=int(input("i=")) # mesh numbers in one direction
N=int(input("N=")) # sphere numbers desired
radius=0.001

for i in range (1,i):
    mn,mx=Vector3(0,0,0),Vector3(0.07,0.07,0)
    p=(mx[0]-mn[0])/N*random.randint(0,N)
    q=(mx[1]-mn[1])/N*random.randint(0,N)
    O.bodies.append(sphere((p,q,0),radius=radius))
    for b in O.bodies:
        if isinstance(b.shape,Sphere):
            if (abs(b.state.pos[0] - p) > radius):
                O.bodies.append(sphere((p,q,0),radius=radius))
            else:
                continue

Sorry for this very late reply. I do not completely understand your script. There are several mistakes in this section:

    for b in O.bodies:
        if isinstance(b.shape,Sphere):
            if (abs(b.state.pos[0] - p) > radius):
                O.bodies.append(sphere((p,q,0),radius=radius))
            else:
                continue

Did you get it fixed meanwhile?

Christian

Is this what you need?

#!/usr/bin/python
# -*- coding: utf-8 -*-

''' script for generating clumps in 2D'''

### user input:
numberOfClumps = 100
radius = 0.001
#boundaries:
mn,mx=Vector3(0,0,0),Vector3(0.07,0.07,0)

### definition of material and engines:

#define material for all bodies:
id_Mat=O.materials.append(FrictMat(young=1e7,poisson=0.3,density=1000,frictionAngle=1))
Mat=O.materials[id_Mat]

#define engines:
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 NewtonIntegrator(damping=0.7,gravity=[0,-10,0]) #gravity in y-direction
]

### generation of a 2D clump assembly
notCreated = 0
for ii in range(0,numberOfClumps):
 #get a random point within boundaries:
 rx = mn[0] + mx[0]*random.random()
 ry = mn[1] + mx[1]*random.random()
 point = Vector3(rx,ry,0)
 print point
 #now we check if clump members, that should be generated from this point would have overlaps with other clump members:
 everythingIsFine = True
 memberPointList = []
 for x in [radius,-radius]:
  for y in [radius,-radius]:
   xmember = rx + x
   ymember = ry + y
   memberPoint = Vector3(xmember,ymember,0)
   memberPointList.append(memberPoint)
   for b in O.bodies:
    if isinstance(b.shape,Sphere):
     distance = memberPoint - b.state.pos
     if abs(distance.norm()) < 2*radius:
      everythingIsFine = False #overlapping spheres
 if everythingIsFine:
  memberList = []
  for memberPoint in memberPointList:
   memberList.append(O.bodies.append(sphere(memberPoint,radius=radius,material=Mat)))
  #clump them together and adapt masses:
  idClump=O.bodies.clump(memberList)
  massInfo = O.bodies.adaptClumpMasses([],10000)
  #make it 2D:
  O.bodies[idClump].blockedDOFs='zXY'
 else:
  notCreated +=1

print '%i of specified %i clumps created'%(numberOfClumps-notCreated,numberOfClumps)

#create a box:
O.bodies.append(box(((mx[0]-mn[0])/2,-2*radius,0),((mx[0]-mn[0])/2,0.001,0.001),fixed=True,material=Mat))

O.dt = 1e-5

Dear Christian Jakob and all producers:
     Working for a long time in these days made me forget to check emails. Thanks for provide me such a perfect script and sorry for this very late responding.
     I have taken a correct running in YADE with your script and there is nothing wrong in it. I think that the core conception of getting clumps generation is according to the distances between any two points. If the value of this distance is larger than (2*radius), allow to generate the clumps, else allow to generate no. I am afraid that that method can only be employed in small-scale simulations, but in large-scale it can not be used any more for that it need to spend a long time in determinative calculations. For instance, I wanna get 1000 clumps, the times of determinative calculations is approximately (1000-1)!. It is a terrible thing when 10000 clumps should be needed.
     As a student pythonner, I have finished my own codes for solution with very simple and basic commands in python, that can be described as that:

#################################################################################################
from yade import pack,os,utils
import random
from numpy import arange

########################################## clumps generation in 2D-simulation ###################
i = (input( "the number of clumps " ) ) # sphere numbers desired
N = int (pow (2*i,0.5) ) # assigning the number of meshes and nodes of pack area, which is 2*i.
mn = Vector3 ( 0,0,0 ) ; mx = Vector3 ( 0.07,0.07,0 )
areaX = mx[0] - mn[0] ; areaY = mx[1] - mn[1]
coorint = [] ; coordes = [] ; sphereList = []

if areaX > areaY :
   radius = (mx[1] - mn[1]) / N * 0.2

else:
   radius = (mx[0] - mn[0]) / N * 0.2

for i in xrange(i):
    coorX = (mx[0]-mn[0]) / N * random.randint(0,N)
    coorY = (mx[1]-mn[1]) / N * random.randint(0,N)
    coorint.append([coorX,coorY])

    if not coorint[i] in coordes:
       coordes.append(coorint[i])

print len(coorint) ; print len(coordes)

for b in coordes:
    coorb1 = [b[0] - radius , b[1] - radius , 0 ]
    coorb2 = [b[0] + radius , b[1] - radius , 0 ]
    coorb3 = [b[0] + radius , b[1] + radius , 0 ]
    coorb4 = [b[0] - radius , b[1] + radius , 0 ]

    sphere1 = O.bodies.append ( sphere ( coorb1 , radius = radius , color=(1,1,1) ) )
    sphere2 = O.bodies.append ( sphere ( coorb2 , radius = radius , color=(1,1,1) ) )
    sphere3 = O.bodies.append ( sphere ( coorb3 , radius = radius , color=(1,1,1) ) )
    sphere4 = O.bodies.append ( sphere ( coorb4 , radius = radius , color=(1,1,1) ) )
    sphereList.append ([sphere1 , sphere2 , sphere3 , sphere4])

print len(sphereList)

for b in sphereList:
    clump_ids = O.bodies.clump (b)

for b in O.bodies:
    O.bodies[clump_ids].blockedDOFs='zXY'

print ( O.bodies[clump_ids].shape.members.keys()[0] )

##################################### some codes only for testing ############
id_Mat=O.materials.append(FrictMat(young=1e7,poisson=0.3,density=1000,frictionAngle=1))
Mat=O.materials[id_Mat]

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 NewtonIntegrator(damping=0.7,gravity=[0,0,-10]) #gravity in z-direction
]

O.bodies.append(box(((mx[0]-mn[0])/2,-2*radius,0),((mx[0]-mn[0])/2,0.001,0.001),fixed=True,material=Mat))

O.dt = 0.2e-4
##########################################################################################

My conception of these codes can be described step by step as follows:
(1).generate meshes and nodes in pack area, which are related to the number of clumps desired. After it is assigned(the number is i), the whole pack area should be divided equally into (2*i) and the number of nodes in this area is also equal to it and it is obvious that the coordinates of all these nodes can be calculated by:

          N = √(2*i)
          coorX = (mx[0]-mn[0]) / N * random.randint(0,N)
          coorY = (mx[1]-mn[1]) / N * random.randint(0,N)

All these things done here is to assure these two things: one is that the number of nodes is enough for generating clumps and the other is that the minimum distance of any two nodes can be easily predicted.

(2).After determining the types and shapes of desired clumps, the radius of spheres for creating clumps can be confirmed according to the minimum distance of any two nodes so that there is no overlaps occuring.

(3).make clumps with nodes and input some necessary testing codes in script in order to essure that all steps of generation can be sucessfully done one by one.

Though there is nothing wrong after testing, the srcipt is not flexible enough for adapting complex requirements for that it can only be used in unique type and size clumps generation. And there are some problems occuring in these codes. So my questions today are that:

(1).What can I do to optimize this srcipt for adapting complex requirements(such as random orientation and size of clumps)

(2).How to get the information of all clumps. I have used this command "print ( O.bodies[clump_ids].shape.members.keys()[0] )" but it does not work any more.

(3).How to define DOF condition to all clumps. I have used this command "for b in O.bodies:
O.bodies[clump_ids].blockedDOFs='zXY' " and then have done "gravity in z-direction" with codes "NewtonIntegrator(damping=0.7,gravity=[0,0,-10])" but there is still a movement in z-direction. What can I do to correct this error.

Seeking your help!

to (1): well, maybe include another random parameter for radii ?!

to (2): have a look at releaseFromClump-example.py ... there is a definition you can use:

def getClumpInfo():
 for b in O.bodies:
  if b.isClump:
   print 'Clump ',b.id,' has following members:'
   keys = b.shape.members.keys()
   for ii in range(0,len(keys)):
    print '- Body ',keys[ii]

Inside the loop you can access information about the clump members by:

O.bodies[keys[ii]].<Info>

e.g. get positions and density:

print O.bodies[keys[ii]].state.pos
print O.bodies[keys[ii]].mat.density

to (3): if clump_ids is a list of clump ids, then you may try:

for ii in clump_ids:
  O.bodies[ii].blockedDOFs='zXY' #for 2D simulation in x-y-plane

OR

for ii in clump_ids:
  O.bodies[ii].blockedDOFs='yXZ' #for 2D simulation in x-z-plane

OR

for ii in clump_ids:
  O.bodies[ii].blockedDOFs='xYZ' #for 2D simulation in y-z-plane

or you loop over all bodies and make something like this:

for b in O.bodies:
  b.blockedDOFs='zXY' #for 2D simulation in x-y-plane

#note: if you have 2D simulation in x-y direction, gravity should act either in x or in y direction ... ^^

Hope it helps,

Christian

Dear Christian Jakob and all producers:
     Thanks for helping me last time and now I can sucessfully get 2D clumps using these codes:

#############################################################################################
from yade import pack,os,utils
import random
from numpy import arange

def Parameters(): # input some basic parameters for clumps generation

    global q,mn,mx,lengX,lengY
    q = (input( "the number of clumps " ) )
    mn,mx=Vector3(0,0,0),Vector3(0.07,0.07,0.1)
    [lengX,lengY] = [(mx[0] - mn[0]),(mx[1] - mn[1])]

def Radius(): # assign radius according to the size of pack area

    global NX,NY,radius

    if lengX > lengY :

       NX = int ( pow ( 2 * q , 0.5 ) )
       NY = int ( ( lengY / lengX ) * NX )
       rx = lengX / NX ; ry = lengY / NY

       if rx > ry :
          radius = ry * 0.22
       else :
          radius = rx * 0.22

    else:

       NY = int ( pow ( 2 * q , 0.5) )
       NX = int ( ( lengX / lengY ) * NY )
       rx = lengX / NX ; ry = lengY / NY

       if rx > ry :
          radius = ry * 0.22
       else :
          radius = rx * 0.22

    print ('The X-direction is equally divided into {0}'.format(NX) )
    print ('The Y-direction is equally divided into {0}'.format(NY) )
    print ('The radius of clumps desired is {0}'.format(radius) )

def Coor(): # determine the coordinations of each point for generating

    global coorint,coordes
    coorint = []
    coordes = []

    for i in xrange(q):
        coorX = lengX / NX * random.randint( 0, NX )
        coorY = lengY / NY * random.randint( 0, NY )
        coorint.append( [coorX , coorY] )

        if not coorint[i] in coordes:
           coordes.append (coorint[i])

    a = len(coorint); b = len(coordes)
    print ('There are {0} clumps desired'.format(a) )
    print ('There are {0} clumps generated'.format(b) )

def Unitcells(): # define unitcell for clumps

    global sphereList
    sphereList = []

    for b in coordes:
        coors1 = [ b[0] - radius , b[1] - radius , 0 ]
        coors2 = [ b[0] + radius , b[1] - radius , 0 ]
        coors3 = [ b[0] + radius , b[1] + radius , 0 ]
        coors4 = [ b[0] - radius , b[1] + radius , 0 ]
        spherepara1 = sphere ( coors1 , radius )
        spherepara2 = sphere ( coors2 , radius )
        spherepara3 = sphere ( coors3 , radius )
        spherepara4 = sphere ( coors4 , radius )
        sphere1 = O.bodies.append ( spherepara1 )
        sphere2 = O.bodies.append ( spherepara2 )
        sphere3 = O.bodies.append ( spherepara3 )
        sphere4 = O.bodies.append ( spherepara4 )
        sphereList.append ([sphere1 , sphere2 , sphere3 , sphere4])

def Clumps(): # generate the clumps desired

    for b in sphereList:
        clump_ids = O.bodies.clump (b)

def material():

    for b in O.bodies:
        if b.isClumps:
           O.bodies[b].material = O.material[clumpsmat]

def DOF(): # of course, DOF condition only for 2D simulation

    for b in O.bodies:
        if b.isClump:
           b.state.blockedDOFs = 'zXY'

# control flows
Parameters()
Radius()
Coor()
Unitcells()
Clumps()
DOF()

###########################################################################################

And I have taken a correct testing with it in gravity condition simulation. But when it is employed in 2D biaxial compression simulation, it does not work any more. The python shell gave me these warning:

###########################################################################################
File "new-test.py", line 140, in <module>
    Clumps()
  File "new-test.py", line 121, in Clumps
    clump_ids = O.bodies.clump (b)
RuntimeError: Clump::updateProperties: NaNs in eigen-decomposition of inertia matrix?!
###########################################################################################

When I forced to run this simulation, the shell noticed that "ValueError: cannot convert float NaN to integer" I do not understand why and wanna know how to solve this problem. Seeking your help!

hmm, i think your script is not working, as you would expect it. some small modifications may help:

#...

def Unitcells(): # define unitcell for clumps

    global sphereList

    for b in coordes:
     sphereList = []
        coors1 = [ b[0] - radius , b[1] - radius , 0 ]
        coors2 = [ b[0] + radius , b[1] - radius , 0 ]
        coors3 = [ b[0] + radius , b[1] + radius , 0 ]
        coors4 = [ b[0] - radius , b[1] + radius , 0 ]
        spherepara1 = sphere ( coors1 , radius )
        spherepara2 = sphere ( coors2 , radius )
        spherepara3 = sphere ( coors3 , radius )
        spherepara4 = sphere ( coors4 , radius )
        sphere1 = O.bodies.append ( spherepara1 )
        sphere2 = O.bodies.append ( spherepara2 )
        sphere3 = O.bodies.append ( spherepara3 )
        sphere4 = O.bodies.append ( spherepara4 )
        sphereList.append ([sphere1 , sphere2 , sphere3 , sphere4])
        Clumps()

def Clumps(): # generate the clumps desired

    for b in sphereList:
        clump_ids = O.bodies.clump (b)

def DOF(): # of course, DOF condition only for 2D simulation

    for b in O.bodies:
        if b.isClump:
           b.state.blockedDOFs = 'zXY'

# control flows
Parameters()
Radius()
Coor()
Unitcells()
#Clumps()
DOF()

# ...

for your problem with biaxial compression:

clump_ids = O.bodies.clump (b)

will not work since b is a body pointer, try b.id

if it is still not working, please append "new-test.py"

christian

Dear Christian Jakob:
     My whole script is as follow:

##/home/fzg/fu/2D-simulation.py### fundamental details of application ###
# unicode: UTF-8
Filename='2D-simulation-vertex'
from yade import pack,os,utils
import random
from numpy import arange

################################# preprocessing for simulation ##########################################
### prescribing variables and functions for simulation controller ###
# material defination
clumpsmat = O.materials.append(FrictMat(poisson=0.5,density=2500,young=1e10,frictionAngle=0.5))
wallmat = O.materials.append(FrictMat(poisson=0.5,density=0,young=1e10,frictionAngle=0))
# walls defination
mn,mx=Vector3(0,0,0),Vector3(0.07,0.07,0.1)
wallids=O.bodies.append(utils.aabbWalls([mn,mx],thickness=.00001,material=wallmat))
# ThreeDTriaxialEngine defination for initial-state determination(the first calculation step)
triax01=TriaxialStressController(
 wall_bottom_id=wallids[2],wall_top_id=wallids[3],
 wall_left_id=wallids[0],wall_right_id=wallids[1],
 wall_back_id=wallids[4],wall_front_id=wallids[5],
 wall_front_activated = False,wall_back_activated = False,
 internalCompaction=False,
 stressMask=7,
        goal1=30000,
        goal2=30000,
        goal3=30000,
        max_vel=5,
)
# ThreeDTriaxialEngine defination for triaxial compression(the second calculation step)
triax02=TriaxialStressController(
 wall_bottom_id=wallids[2],wall_top_id=wallids[3],
 wall_left_id=wallids[0],wall_right_id=wallids[1],
 wall_back_id=wallids[4],wall_front_id=wallids[5],
 wall_front_activated = False,wall_back_activated = False,
 internalCompaction=False,
 stressMask = 5,
        goal2=-100,
        goal1=30000,
        goal3=30000,
)
################################# control flow for simulation ##########################################
# particles generation
def Parameters(): # input some basic parameters for clumps generation

    global q,mn,mx,lengX,lengY
    q = (input( "the number of clumps " ) )

    [lengX,lengY] = [(mx[0] - mn[0]),(mx[1] - mn[1])]

def Radius(): # assign radius according to the size of pack area

    global NX,NY,radius

    if lengX > lengY :

       NX = int ( pow ( 2 * q , 0.5 ) )
       NY = int ( ( lengY / lengX ) * NX )
       rx = lengX / NX ; ry = lengY / NY

       if rx > ry :
          radius = ry * 0.22
       else :
          radius = rx * 0.22

    else:

       NY = int ( pow ( 2 * q , 0.5) )
       NX = int ( ( lengX / lengY ) * NY )
       rx = lengX / NX ; ry = lengY / NY

       if rx > ry :
          radius = ry * 0.22
       else :
          radius = rx * 0.22

    print ('The X-direction is equally divided into {0}'.format(NX) )
    print ('The Y-direction is equally divided into {0}'.format(NY) )
    print ('The radius of clumps desired is {0}'.format(radius) )

def Coor(): # determine the coordinations of each point for generating

    global coorint,coordes
    coorint = []
    coordes = []

    for i in xrange(q):
        coorX = lengX / NX * random.randint( 0, NX )
        coorY = lengY / NY * random.randint( 0, NY )
        coorint.append( [coorX , coorY] )

        if not coorint[i] in coordes:
           coordes.append (coorint[i])

    a = len(coorint); b = len(coordes)
    print ('There are {0} clumps desired'.format(a) )
    print ('There are {0} clumps generated'.format(b) )

def Unitcells(): # define unitcell for clumps

    global sphereList
    sphereList = []

    for b in coordes:
        coors1 = [ b[0] - radius , b[1] - radius , 0 ]
        coors2 = [ b[0] + radius , b[1] - radius , 0 ]
        coors3 = [ b[0] + radius , b[1] + radius , 0 ]
        coors4 = [ b[0] - radius , b[1] + radius , 0 ]
        spherepara1 = sphere ( coors1 , radius )
        spherepara2 = sphere ( coors2 , radius )
        spherepara3 = sphere ( coors3 , radius )
        spherepara4 = sphere ( coors4 , radius )
        sphere1 = O.bodies.append ( spherepara1 )
        sphere2 = O.bodies.append ( spherepara2 )
        sphere3 = O.bodies.append ( spherepara3 )
        sphere4 = O.bodies.append ( spherepara4 )
        sphereList.append ([sphere1 , sphere2 , sphere3 , sphere4])

def Clumps(): # generate the clumps desired

    for b in sphereList:
        clump_ids = O.bodies.clump (b)

def material():

    for b in O.bodies:
        if b.isClumps:
           O.bodies[b].material = O.material[clumpsmat]

def DOF(): # of course, DOF condition only for 2D simulation

    for b in O.bodies:
        if b.isClump:
           b.state.blockedDOFs = 'zXY'

# control flows
Parameters()
Radius()
Coor()
Unitcells()
Clumps()
DOF()

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_Wall_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax01,
 NewtonIntegrator(damping=.1),
 PyRunner(command='resultfirst()',iterPeriod=1000)
]
# first step of simulation startting with a correct inheriting for the next step
O.dt = 5e-6

O.run(30000,True)
################################################################################

Now I still have no idea about what you just said some small modifications because there is nothing
wrong when only the gravity condition is supplied. Seeking your help!

> Now I still have no idea about what you just said some small modifications because there is nothing

sorry, my fault, everthing looks fine.

your problem comes from density = 0.0 for your spheres. somewhere in the code inertia is devided by density, so you get nans for inertia of clumps ...

please define a material like:

id_Mat=O.materials.append(FrictMat(young=1e7,poisson=0.3,density=1000,frictionAngle=1))
Mat=O.materials[id_Mat]

and set it to your spheres:

        spherepara1 = sphere ( coors1 , radius , material=Mat )
        spherepara2 = sphere ( coors2 , radius , material=Mat )
        spherepara3 = sphere ( coors3 , radius , material=Mat )
        spherepara4 = sphere ( coors4 , radius , material=Mat )

then everthing works fine ;)

Dear Christian Jakob:
     Thannks for giving me that perfect advice. Now I have taken a correct running with modified script in which I input some codes about defining the spheres material. So everything looks fine.
But I still have some problems about this update which can be described as that:
(1). In material defination, YADE requires me to assign such parameters for "O.bodies" twice. For the first time, these parameters should be determined when these spheres are generated but not clumps, it just like that:
##############################################################################################
coors1 = [ b[0] - radius , b[1] - radius , 0 ]
coors2 = [ b[0] + radius , b[1] - radius , 0 ]
coors3 = [ b[0] + radius , b[1] + radius , 0 ]
coors4 = [ b[0] - radius , b[1] + radius , 0 ]
spherepara1 = sphere ( coors1 , radius , material=clumpsmat, color = [1,1,1] )
spherepara2 = sphere ( coors2 , radius , material=clumpsmat, color = [1,1,1] )
spherepara3 = sphere ( coors3 , radius , material=clumpsmat, color = [1,1,1] )
spherepara4 = sphere ( coors4 , radius , material=clumpsmat, color = [1,1,1] )
sphere1 = O.bodies.append ( spherepara1 )
sphere2 = O.bodies.append ( spherepara2 )
sphere3 = O.bodies.append ( spherepara3 )
sphere4 = O.bodies.append ( spherepara4 )
################################################################################################
And for the second time, they should be determined when all clumps generating is finished. it just like that:
################################################################################################
def material():
    for b in O.bodies:
        if b.isClumps:
           O.bodies[b].material = O.material[clumpsmat]
################################################################################################
I now understand that the absence of the first time working may cause something wrong in running this simulation and I am not sure that absence of the second time working may cause what. I think there also be something wrong of course.

(2).If they are both necessary for running this simulation. It would make the script complex because there still exists some repetitive tasks. What can I do to make the whole script more concise.

(3).I find out that some other attributes of clumps defination is not the same as "material-def", such as DOF condition. I also use these two methods for supplying the DOF condition, one is that:
before the clumps generating, I employ these codes for my purpose:
#################################################################################################
for b in O.bodies:
 if isinstance(b.shape,Sphere):
   b.state.blockedDOFs='zXY'
#################################################################################################
This is alao doing directly for spheses. But after the clumps generation, these can not work any more when the simulation starts.

And the other is that:
#################################################################################################
def DOF(): # of course, DOF condition only for 2D simulation
    for b in O.bodies:
        if b.isClump:
           b.state.blockedDOFs = 'zXY'
################################################################################################
There is nothing wrong with it when simulation starts.

Why these differences occur? Seeking your help!

to (1) You just need to define material for spheres. The volume/masses/inertia of clumps is calculated automatically in Yade. I think clumps itself can not have a material.

to (2) concise? in your case there is no need for any definition ... also this can be shortened:

#old:
coors1 = [ b[0] - radius , b[1] - radius , 0 ]
spherepara1 = sphere ( coors1 , radius , material=clumpsmat, color = [1,1,1] )
sphere1 = O.bodies.append ( spherepara1 )

#new:
sphere1 = O.bodies.append ( sphere ( [ b[0] - radius , b[1] - radius , 0 ] , radius , material=clumpsmat, color = [1,1,1] ) )

to (3): you just have to block DOFs for clumps, since you have only clumps in your sim. blocking DOFs of spheres makes no sense in your case.

    for b in O.bodies:
        if b.isClump:
           b.state.blockedDOFs = 'zXY'

is enough

Dear Christian Jakob and all users:
          Long time no see and thanks for helping last time, I can do some analysis with 2D clumps now buy at the same time find some problems with DOF for those. Following my comprehension, all the centroids of spheres(particles) should be at the same geometric surface during the entire time of simulation but this is no in fact. So my all problems can be shown as that:

(1).After the clumps generation finishing, the state is just like that(the beginning of my simulation):
http://i.imgur.com/ZgTCuhW.png

(2).My test simulation is compose of a large number of clumps and 4 fixed wall and all clumps should getting a free falling since
they have contacts with the bottom wall or some other ones below them. If use no function for DOF, the state of my finished tast is of course as that(this is 3D simulation):
http://i.imgur.com/ElBAWMJ.png

(3).After that above, I have another try for employing a function your told me last time for DOF , which can be described as:

def DOF(): # of course, DOF condition only for 2D simulation

    for b in O.bodies:
        if b.isClump:
           b.state.blockedDOFs = 'zXY'

and the state of my finished tast this time is of course as that:
http://i.imgur.com/h0LuITR.png

This can demonstrate obviously that this function does the positive effect for DOF. But I also find that all the centroids of spheres are not in the same plane by checking the status displayed by YADE just in Y-Z plane. It is as that:
http://i.imgur.com/q11ISgH.png

I think there is nothing wrong in what I just said above. But I am not sure the reason of it and it could do harm for my simulation.

Seeking for your help!

Hello Fu,

could you create a minimal working example (i.e. a short script without
unnecessary parts, that reproduces your problem)?

Just an idea: do you have also some standalone (non-clump) particles in the
simulation? if yes, this is the reason (because you prescribe blockedDOFs
only to clumps). If no, the problem is somewhere else :-)
cheers
Jan

2013/8/5 Fu zuoguang <email address hidden>

> Question #230139 on Yade changed:
> https://answers.launchpad.net/yade/+question/230139
>
> Status: Answered => Open
>
> Fu zuoguang is still having a problem:
> Dear Christian Jakob and all users:
> Long time no see and thanks for helping last time, I can do some
> analysis with 2D clumps now buy at the same time find some problems with
> DOF for those. Following my comprehension, all the centroids of
> spheres(particles) should be at the same geometric surface during the
> entire time of simulation but this is no in fact. So my all problems can
> be shown as that:
>
> (1).After the clumps generation finishing, the state is just like that(the
> beginning of my simulation):
> http://i.imgur.com/ZgTCuhW.png
>
> (2).My test simulation is compose of a large number of clumps and 4 fixed
> wall and all clumps should getting a free falling since
> they have contacts with the bottom wall or some other ones below them. If
> use no function for DOF, the state of my finished tast is of course as
> that(this is 3D simulation):
> http://i.imgur.com/ElBAWMJ.png
>
> (3).After that above, I have another try for employing a function your
> told me last time for DOF , which can be described as:
>
> def DOF(): # of course, DOF condition only for 2D simulation
>
> for b in O.bodies:
> if b.isClump:
> b.state.blockedDOFs = 'zXY'
>
> and the state of my finished tast this time is of course as that:
> http://i.imgur.com/h0LuITR.png
>
> This can demonstrate obviously that this function does the positive effect
> for DOF. But I also find that all the centroids of spheres are not in the
> same plane by checking the status displayed by YADE just in Y-Z plane. It
> is as that:
> http://i.imgur.com/q11ISgH.png
>
> I think there is nothing wrong in what I just said above. But I am not
> sure the reason of it and it could do harm for my simulation.
>
> Seeking for your help!
>
> --
> You received this question notification because you are a member of
> yade-users, which is an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

Dear Jan Stránský:
          I am sorry for forgetting to describe these very important detials with a intact script, which can be shown at the end of this conmmunication. There are 2 types of bodies in my simulatiion, one is the wall and the other is of course the clumps. In my own script you can obviously see only 2 functions for Bodies defination. Seeking for you help!

#######################################################################################

from yade import pack,os,utils
import random
from numpy import arange

Mat = O.materials.append(FrictMat(poisson=0.5,density=1500,young=1e7,frictionAngle=0.5))
clumpsmat = O.materials[Mat]
wallmat = O.materials.append(FrictMat(poisson=0.5,density=1500,young=1e7,frictionAngle=0))

def Parameters(): # input some basic parameters for clumps generation

    global sphnums,radius,clunums,callenY
    clunums = (input( "Enter the number of clumps: " ) )
    sphnums = (input( "Enter the number of spheres employed for generating only one clump: " ) )
    radius = (input( "Enter the radius of master particle: " ) )
    callenY = (input( "Enter the ratio of box-lengthY/lengthX: " ) )

    global ratios,radii,overlaps,distances,unitcellrads
    ratios = []; radii = []; overlaps = []; distances = []; unitcellrads=[]

    for i in xrange(sphnums-1):
        Enti01 = 'Enter the radius-ratio of number {0} slave sphere of clump: '.format(i+2)
        rat = input (Enti01)
        ratios.append(rat)

    for i in xrange(sphnums-1):
        Enti02 = 'Enter the overlap of number {0} slave sphere with the master: '.format(i+2)
        ola = input (Enti02)
        overlaps.append(ola)

    unitcellrads.append( radius )
    for i in xrange(sphnums-1):
        radii.append( radius * rat)
        Priti01 = 'the radius of number {0} slave sphere: '.format(i+2)
# print Priti01, radii[i]
        unitcellrads.append( radii[i] )

    for i in xrange(sphnums-1):
        distances.append( (1+rat)*(1-ola)*radius )
        Priti02 = 'distance of number {0} slave sphere with the master: '.format(i+2)
# print Priti02, distances[i]

    global NX,rx,lengX,lengY,NY
    NX = int ( pow ( 0.55*clunums, 0.5 ) ) + 1
    rx = 1.1 * (2 * distances[0] + 2 * radii[0])
    lengX = NX * rx
    lengY = lengX * callenY
    NY = int( lengY / rx ) - 1

def Coor(): # determine the coordinations of each point for generating

    global mastcoorinp,mastcoordes,slavecoordes,coordes
    mastcoorinp = []; mastcoordes = []
    slavecoordes = []; coordes = []

    for i in xrange(clunums):
        coorX = lengX / NX * random.randint( 0, NX )
        coorY = lengY / NY * random.randint( 0, NY )
        mastcoorinp.append( [coorX, coorY, 0] )

        if not mastcoorinp[i] in mastcoordes:
           mastcoordes.append (mastcoorinp[i])

    for i in xrange( len(mastcoordes) ):
        slacoor = []
        for j in xrange( len(radii) ):
            ramx = (-1)**random.randint(0,1)
            ramy = (-1)**random.randint(0,1)
            sina = ramy*random.random()
            cosa = ramx*pow( (1-sina*sina), 0.5 )
            slavecoorX = mastcoordes[i][0] + distances[j] * cosa
            slavecoorY = mastcoordes[i][1] + distances[j] * sina
            sphcoor = []
            sphcoor = ( [slavecoorX, slavecoorY, 0] )
            slacoor.append(sphcoor)
        slavecoordes.append(slacoor)

    for i in xrange( len(mastcoordes) ):
        clucoor = []
        clucoor.append(mastcoordes[i])
        for j in xrange( len(slavecoordes[i]) ):
            clucoor.append(slavecoordes[i][j])
        coordes.append(clucoor)
# print coordes

def spheres(): # define unitcell for clumps

    global spherespara,sphereslist,spheres_con
    spherespara = []; sphereslist = [];spheres_con = []

    for i in xrange( len(mastcoordes) ):
        cluparas = []
        for j in xrange( len(unitcellrads) ):
            clupara = []
            clupara.append([coordes[i][j],unitcellrads[j]])
            cluparas.append(clupara)
        spherespara.append(cluparas)

    for i in xrange( len(spherespara) ):
        spheres_con = []
        for j in xrange( len(spherespara[i]) ):
            spherepar = []
            spheres = []
            if j == 0 :
               spherepar = sphere (spherespara[i][j][0][0], radius=spherespara[i][j][0][1], color = [1,1,1])
               spheres = O.bodies.append(spherepar)
            else :
               spherepar = sphere (spherespara[i][j][0][0], radius=spherespara[i][j][0][1], color = [1,0,0])
               spheres = O.bodies.append(spherepar)
            spheres_con.append(spheres)
        sphereslist.append(spheres_con)

def clumps():
    for b in sphereslist:
        O.bodies.clump (b)

def material():

    for b in O.bodies:
        if b.isClump:
           O.bodies.material = clumpsmat

def DOF(): # of course, DOF condition only for 2D simulation

    for b in O.bodies:
        if b.isClump:
           b.state.blockedDOFs = 'zXY'

def walls():
    lwall = 0 - (max(distances)+radius);
    rwall = lengX + (max(distances)+radius);
    bwall = 0 - (max(distances)+radius);
    twall = lengY + (max(distances)+radius);
    mn,mx=Vector3(lwall,bwall,-0.05),Vector3(rwall,twall,0.05)
    wallids=O.bodies.append(utils.aabbWalls([mn,mx],thickness=.00001,material=wallmat))

def dynamiccontrol():
    a = len( mastcoorinp ); b = len( mastcoordes )
    print ( 'There are {0} clumps desired'.format(a) )
    print ('But only {} clumps generated'.format(b) )

# control flows
Parameters()
Coor()
walls()
spheres()
clumps()
material()
DOF()
dynamiccontrol()

#define engines:
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 NewtonIntegrator(damping=0.7,gravity=[0,-10,0]) #gravity in y-direction
]

This question was expired because it remained in the 'Open' state without activity for the last 15 days.