Yade

Avoid interaction between two particles in a four particles experiment

Asked by ZHU Huaxiang

Hello all,

I am doing a special biaxial test which only content Four particles laying to be a rhombus, like this

                         1
                         o
                 2 o --- o 0
                         o
                         3

how can I avoid interaction between axial line particles 1 and 3 if they contact?

I already tried the 'mask' and 'avoidSelfInteractionMask' to cooperate, as

O.bodies[0].groupMask=1;O.bodies[2].groupMask=1
O.bodies[1].groupMask=2;O.bodies[3].groupMask=2
[engines]collider.avoidSelfInteractionMask=2

but when the 1,3 interaction is avoided it is also impossible to establish interactions between two different mask's particles, then all interaction will lose by using 'avoidSelfInteractionMask'.

Do you have any idea to achieve that? thank you.

Huaxiang

here is my code
-------------------------------------------------------------------------------------------------------

import math
from yade import plot

conf=100000
poisson=0.5
kn=1.0e8
E=kn/1.0
ks=poisson*kn
rate=5e-4
status=[0.0,0.0,0.0,0.0]

O.materials.append(FrictMat(young=E,poisson=poisson,frictionAngle=radians(99.0),density=2600,label='spheres'))
O.materials.append(FrictMat(young=1.0e3*kn,poisson=0.3,frictionAngle=radians(99.0),density=0,label='walls'))

mn,mx=Vector3(-math.sqrt(2)-1,-math.sqrt(2)-1,-1.1),Vector3(math.sqrt(2)+1,math.sqrt(2)+1,1.1)
walls=utils.aabbWalls([mn,mx],thickness=0.001,material='walls')
wallIds=O.bodies.append(walls)

O.bodies.append([
   utils.sphere(center=(math.sqrt(2),0.0,0.0),radius=1.0,material='spheres'),
   utils.sphere(center=(0.0,math.sqrt(2),0.0),radius=1.0,material='spheres'),
   utils.sphere(center=(-math.sqrt(2),0.0,0.0),radius=1.0,material='spheres'),
   utils.sphere((0.0,-math.sqrt(2),0.0),1.0,material='spheres')
])

for k in O.bodies:
   if isinstance(k.shape, Sphere): k.state.blockedDOFs='zXY'

O.dt=.5*utils.PWaveTimeStep()
O.usesTimeStepper=True

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()]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8, defaultDt=4*utils.PWaveTimeStep()),
 NewtonIntegrator(damping=0.01),
]

O.dt=.5e-3*utils.PWaveTimeStep()

def setVel(velX,velY):
  O.bodies[0].state.vel=Vector3(velX,0.0,0.0)
  O.bodies[1].state.vel=Vector3(-velX,0.0,0.0)
  O.bodies[2].state.vel=Vector3(0.0,velY,0.0)
  O.bodies[3].state.vel=Vector3(0.0,-velY,0.0)

rateTmp=rate
while 1:
  meanS=(abs(O.forces.f(1)[0])+abs(O.forces.f(3)[1]))
  unb=unbalancedForce()
  if rateTmp>0.0 and meanS>conf:
     print 'convert'
     rateTmp=-0.5*rateTmp
  elif rateTmp<0.0 and meanS<conf:
     print 'convert'
     rateTmp=-0.5*rateTmp
  setVel(rateTmp,rateTmp)
  O.run(1, True)
  unb=unbalancedForce()
  meanS=(abs(O.forces.f(1)[0])+abs(O.forces.f(3)[1]))*0.5
  print 'unbalanced force:',unb,' mean stress: ',meanS
  if unb<0.001 and abs(meanS-conf)/conf<0.005:
    break

plot.plots={'e2':'q'}
plot.plot()
raw_input('confining finish')

'''
ini=O.bodies[7].state.se3[0][1]
e22=0.0
e11=0.0
#rate=1e2*rate
while 1:
   setVel(0.0,rate);e22+=rate
   O.run(1,True)
   while unbalancedForce()>0.001:
       setVel(0.0,0.0)
       O.run(1,True)
   rateTmp=-0.5*rate
   lateralS=(abs(O.forces.f(0)[0])+abs(O.forces.f(1)[0]))*0.5
   while abs(lateralS-conf)/conf>0.001:
       lateralS=(abs(O.forces.f(0)[0])+abs(O.forces.f(1)[0]))*0.5
       if rateTmp>0.0 and lateralS>conf:
          rateTmp=-0.5*rateTmp
       elif rateTmp<0.0 and lateralS<conf:
          rateTmp=-0.5*rateTmp
       setVel(rateTmp,0.0)
       O.run(1, True);e11+=rateTmp
       lateralS=(abs(O.forces.f(0)[0])+abs(O.forces.f(1)[0]))*0.5
   print O.forces.f(3)[1]-O.forces.f(1)[0],O.forces.f(1)[0],(ini-O.bodies[7].state.se3[0][1])/ini
   plot.addData(e2=e22,e1=e11,q=O.forces.f(3)[1]-O.forces.f(1)[0],i=O.iter)
'''

rate=10.0*rate
adjust=math.cos(75.0*pi/180.0)/math.sin(75.0*pi/180.0)*rate*1e-1
ini=O.bodies[7].state.se3[0][1]
e22=0.0
e11=0.0
rateX=0.0
while 1:
   lateralS=(abs(O.forces.f(0)[0])+abs(O.forces.f(1)[0]))*0.5
   if lateralS>conf:
       adjust=-abs(adjust)
   elif lateralS<conf:
       adjust=abs(adjust)
   rateX+=adjust
   setVel(rateX,rate);e22+=rate
   O.run(1,True)
   while unbalancedForce()>0.01:
       setVel(0.0,0.0)
       O.run(1,True)
   print O.forces.f(3)[1]-O.forces.f(1)[0],O.forces.f(1)[0],(ini-O.bodies[7].state.se3[0][1])/ini,unbalancedForce()
   plot.addData(e2=e22,e1=e11,q=O.forces.f(3)[1]-O.forces.f(1)[0],i=O.iter)

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
Last query:
Last reply:
Revision history for this message
ZHU Huaxiang (zhuhuaxiang1985) said : 		#1

By the way, I am playing a middle school physics teaching software named 'Algodoo', run by physical engine.
Each body can be set in two mask or more, if explained by yade code, it is like

O.bodies[0].groupMask=[1];O.bodies[2].groupMask=[1]
O.bodies[1].groupMask=[1,2];O.bodies[3].groupMask=[1,2]
[engines]collider.avoidSelfInteractionMask=[2]

in the end, interactions 0,1 ; 1,2 ; 2,3 ; 3,4 can function as usual, moreover, interaction 1,3 which I d like to avoid is really ignored.

I have already solved this problem by Algodoo, I think it will be interesting that Yade have this configuration.

by the way, this is a very funny software for all ages, enjoy it if you want.
http://www.algodoo.com/

Revision history for this message
Jérôme Duriez (jduriez) said : 		#2

Hi,

Use mask =3 (=11 in binary) for your bodies 1 and 3.

In this case it still has one bit in common with other bodies 0 and 2 (where mask=1 = 01 in binary), so bodies 1 or 3 may interact with bodies 0/2. Masks do not have to be strictly equal so that bodies interact.

And, if I'm right, collider.avoidSelfInteractionMask can in fact in this case still be equal to 2 (=10 in binary) because it has one bit in common with 3

Revision history for this message
Best Jan Stránský (honzik) said : 		#3

Hello Zhu,

you can achive it (as you mentioned) by setting groupMask. Let me explain
on your example:

O.bodies[0].mask = 1 # 0b0001
O.bodies[1].mask = 2 # 0b0010
O.bodies[2].mask = 1 # 0b0001
O.bodies[3].mask = 2 # 0b0010

where 0bXXXX stands for number binary representation. The interaction of
two bodies is avoided, if mask1&mask2==0, see [1], section "AND". You can
play e.g. in IPython:

a = 0b0001
b = 0b0010
a&b
# etc..

so in your case, 1&1 = 1, 2&2 = 2, 1&2 = 0, so body pairs 0-2 and 1-3 will
interact and no other pair. To achieve what you want, do something like

O.bodies[0].mask = 0b11 # 3
O.bodies[1].mask = 0b01 # 1
O.bodies[2].mask = 0b11 # 3
O.bodies[3].mask = 0b10 # 2

in this case 3&1=1, 3&2=2, 3&3=3, 2&1=0, 2&2=2, 1&1=1. So all particles
will interact witch each other, except the pairs with mask1=1 and mask2=2,
as 1&2=0. I.e. particles 1 and 3 in your example will not interact with
each other.

Maybe we could improve the documentation about mask attribute a little..

cheers
Jan

[1] http://en.wikipedia.org/wiki/Bitwise_operation

2013/12/4 ZHU Huaxiang <email address hidden>

> Question #240275 on Yade changed:
> https://answers.launchpad.net/yade/+question/240275
>
> Description changed to:
> Hello all,
>
> I am doing a special biaxial test which only content Four particles laying
> to be a rhombus, like this
>
> 1
> o
> 2 o o 0
> o
> 3
>
> how can I avoid interaction between axial line particles 1 and 3 if they
> contact?
>
> I already tried the 'mask' and 'avoidSelfInteractionMask' to cooperate,
> as
>
> O.bodies[0].groupMask=1;O.bodies[2].groupMask=1
> O.bodies[1].groupMask=2;O.bodies[3].groupMask=2
> [engines]collider.avoidSelfInteractionMask=2
>
> but when the 1,3 interaction is avoided it is also impossible to
> establish interactions between two different mask's particles, then all
> interaction will lose by using 'avoidSelfInteractionMask'.
>
> Do you have any idea to achieve that? thank you.
>
>
> Huaxiang
>
>
> here is my code
>
> -------------------------------------------------------------------------------------------------------
>
> import math
> from yade import plot
>
> conf=100000
> poisson=0.5
> kn=1.0e8
> E=kn/1.0
> ks=poisson*kn
> rate=5e-4
> status=[0.0,0.0,0.0,0.0]
>
>
> O.materials.append(FrictMat(young=E,poisson=poisson,frictionAngle=radians(99.0),density=2600,label='spheres'))
>
> O.materials.append(FrictMat(young=1.0e3*kn,poisson=0.3,frictionAngle=radians(99.0),density=0,label='walls'))
>
>
> mn,mx=Vector3(-math.sqrt(2)-1,-math.sqrt(2)-1,-1.1),Vector3(math.sqrt(2)+1,math.sqrt(2)+1,1.1)
> walls=utils.aabbWalls([mn,mx],thickness=0.001,material='walls')
> wallIds=O.bodies.append(walls)
>
> O.bodies.append([
>
> utils.sphere(center=(math.sqrt(2),0.0,0.0),radius=1.0,material='spheres'),
>
> utils.sphere(center=(0.0,math.sqrt(2),0.0),radius=1.0,material='spheres'),
>
> utils.sphere(center=(-math.sqrt(2),0.0,0.0),radius=1.0,material='spheres'),
> utils.sphere((0.0,-math.sqrt(2),0.0),1.0,material='spheres')
> ])
>
>
> for k in O.bodies:
> if isinstance(k.shape, Sphere): k.state.blockedDOFs='zXY'
>
> O.dt=.5*utils.PWaveTimeStep()
> O.usesTimeStepper=True
>
>
> 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()]
> ),
> GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8,
> defaultDt=4*utils.PWaveTimeStep()),
> NewtonIntegrator(damping=0.01),
> ]
>
> O.dt=.5e-3*utils.PWaveTimeStep()
>
> def setVel(velX,velY):
> O.bodies[0].state.vel=Vector3(velX,0.0,0.0)
> O.bodies[1].state.vel=Vector3(-velX,0.0,0.0)
> O.bodies[2].state.vel=Vector3(0.0,velY,0.0)
> O.bodies[3].state.vel=Vector3(0.0,-velY,0.0)
>
> rateTmp=rate
> while 1:
> meanS=(abs(O.forces.f(1)[0])+abs(O.forces.f(3)[1]))
> unb=unbalancedForce()
> if rateTmp>0.0 and meanS>conf:
> print 'convert'
> rateTmp=-0.5*rateTmp
> elif rateTmp<0.0 and meanS<conf:
> print 'convert'
> rateTmp=-0.5*rateTmp
> setVel(rateTmp,rateTmp)
> O.run(1, True)
> unb=unbalancedForce()
> meanS=(abs(O.forces.f(1)[0])+abs(O.forces.f(3)[1]))*0.5
> print 'unbalanced force:',unb,' mean stress: ',meanS
> if unb<0.001 and abs(meanS-conf)/conf<0.005:
> break
>
> plot.plots={'e2':'q'}
> plot.plot()
> raw_input('confining finish')
>
> '''
> ini=O.bodies[7].state.se3[0][1]
> e22=0.0
> e11=0.0
> #rate=1e2*rate
> while 1:
> setVel(0.0,rate);e22+=rate
> O.run(1,True)
> while unbalancedForce()>0.001:
> setVel(0.0,0.0)
> O.run(1,True)
> rateTmp=-0.5*rate
> lateralS=(abs(O.forces.f(0)[0])+abs(O.forces.f(1)[0]))*0.5
> while abs(lateralS-conf)/conf>0.001:
> lateralS=(abs(O.forces.f(0)[0])+abs(O.forces.f(1)[0]))*0.5
> if rateTmp>0.0 and lateralS>conf:
> rateTmp=-0.5*rateTmp
> elif rateTmp<0.0 and lateralS<conf:
> rateTmp=-0.5*rateTmp
> setVel(rateTmp,0.0)
> O.run(1, True);e11+=rateTmp
> lateralS=(abs(O.forces.f(0)[0])+abs(O.forces.f(1)[0]))*0.5
> print
> O.forces.f(3)[1]-O.forces.f(1)[0],O.forces.f(1)[0],(ini-O.bodies[7].state.se3[0][1])/ini
> plot.addData(e2=e22,e1=e11,q=O.forces.f(3)[1]-O.forces.f(1)[0],i=O.iter)
> '''
>
> rate=10.0*rate
> adjust=math.cos(75.0*pi/180.0)/math.sin(75.0*pi/180.0)*rate*1e-1
> ini=O.bodies[7].state.se3[0][1]
> e22=0.0
> e11=0.0
> rateX=0.0
> while 1:
> lateralS=(abs(O.forces.f(0)[0])+abs(O.forces.f(1)[0]))*0.5
> if lateralS>conf:
> adjust=-abs(adjust)
> elif lateralS<conf:
> adjust=abs(adjust)
> rateX+=adjust
> setVel(rateX,rate);e22+=rate
> O.run(1,True)
> while unbalancedForce()>0.01:
> setVel(0.0,0.0)
> O.run(1,True)
> print
> O.forces.f(3)[1]-O.forces.f(1)[0],O.forces.f(1)[0],(ini-O.bodies[7].state.se3[0][1])/ini,unbalancedForce()
> plot.addData(e2=e22,e1=e11,q=O.forces.f(3)[1]-O.forces.f(1)[0],i=O.iter)
>
> --
> 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
>

Revision history for this message
ZHU Huaxiang (zhuhuaxiang1985) said : 		#4

Thank jduriez very much!!

You are right, 1 and 3 interact avoid, other contact normally function.

I am still wondering in what rule the mask determine interaction or not.
if there is a bit series in common at same place in binary, two apparently unequal mask can contact?

Revision history for this message
Bruno Chareyre (bruno-chareyre) said : 		#5

>if explained by yade code, it is like
>O.bodies[1].groupMask=[1,2]

Jerome's suggestion is:
> O.bodies[1].groupMask=3

Which is indeed the correct translation of [1,2] (Algodoo) into a bitmask (Yade). It seems to me that any combination that you can define in the Algodoo way can similarly be reproduced with bitmasks.

Note that there is no absolute need to use avoidSelfInteractionMask if you stick to four particles only.
This should be ok:
mask0 = 11 (0 interacts with 1 and 3 => 1+2+8=11)
mask1 = 7 (1 interacts with 0 and 2 => 2+1+4=7)
mask2 = 14 (etc.)
mask3 = 13

Since this technique can be painfull to generalize to N particles, the purpose of avoidSelfInteractionMask was to make some logical operations easier in the general case.

Revision history for this message
Bruno Chareyre (bruno-chareyre) said : 		#6

Forget about my last note. It will absolutely not work. Bitmasks are sometimes confusing for me as well...
If you understand why 3 is the translation of [1,2] (Jan's post should help) then you should find your way.

Revision history for this message
ZHU Huaxiang (zhuhuaxiang1985) said : 		#7

Thanks Jan Stránský, that solved my question.

Revision history for this message
ZHU Huaxiang (zhuhuaxiang1985) said : 		#8

Thank you all for the very detailed explanation, problem solved!