on modify interactionLoop

As I understand, changing material for active contacts (interactions)
will not have any effect, only for new contacts.

You have to change those parameters directly in interactions (kn, ks etc.).
I am not sure, whether I am right,

Anton

But changing material parameters will perform, say changing Friction angle will have effect during the deviatoric loading

To "localize" the problem, you can simulate just 2 spheres, applying to them
different conditions and measure the response.

Anton

O.load('consolidated_state_'+key+'_'+str(phi)+'.xml') # we add new
materials to our material container
O.materials.append(MomentMat(eta=eyta,young=young,poisson=0.1,frictionAngle=radians(fricFinal),density=2600,label='spheres-moment'))
O.materials.append(MomentMat(eta=0.0,young=young,poisson=0.1,frictionAngle=0,density=0,label='walls-moment'))
# here we modify the material type for i_body in O.bodies: if i_body.id
< 6 : i_body.material = O.materials[3] else: i_body.material =
O.materials[2] #here we move to modify the interaction loop Interaction
= O.engines[2] Interaction = InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
[Ip2_MomentMat_MomentMat_MomentPhys(Knormal=knormal,Alpha =
alpha,Beta=beta,useAlphaBeta=True,userInputStiffness=True)],
[Law2_SCG_MomentPhys_CohesionlessMomentRotation()] )

Please send whole script.. what do you do with Interaction instance?
O.engines is processed differently then other python objects and must be
assigned at once, one item change has no effect [1]

Jan

[1] https://www.yade-dem.org/doc/user.html#base-engines

I would say that Anton is right, you need to change contact properties for every existing interactions cause changing bodies properties won't be effective for existing interactions (interactions are not updated), only new interactions will benefit from the new bodies properties.

so you should do something like

for b in O.bodies:
  b.mat.frictionAngle=radians(simulationFrictDeg)

AND

for i in O.interactions:
  i.phys.frictionAngle=radians(simulationFrictDeg)

Luc

Or just:

tt=TriaxialCompressionEngine()
tt.setContactProperties(30)

but it only assign friction. For stiffness and/or cohesion, you need Luc's loops.

Hi Luc,

 I know change friction angle using either your technique or Chareyre's are ok, but it is for FrictMat.

Now I am using Moment Material and Moment Law, changing kn ks eta and beta result same response using your technique,

I don't know why. My code is below

        O.load('consolidated_state_'+key+'_'+str(phi)+'.xml')

# we add new materials to our material container
        O.materials.append(MomentMat(eta=eyta,young=young,poisson=0.1,frictionAngle=radians(fricFinal),density=2600,label='spheres-moment'))
        O.materials.append(MomentMat(eta=0.0,young=young,poisson=0.1,frictionAngle=0,density=0,label='walls-moment'))

        # here we modify the material and interaction physics
        for i_body in O.bodies:
                if i_body.id < 6 :
                        i_body.material = O.materials[3]
                else:
                        i_body.material = O.materials[2]
        for i_interactions in O.interactions:
                if i_interactions.id1 > 5 and i_interactions.id2 > 5:
                        i_interactions.Eta = eyta # only interaction between particles are changed to eyta > 0

        #here we move to modify the interaction loop
        Interaction = O.engines[2]
        Interaction = InteractionLoop(
                [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
                [Ip2_MomentMat_MomentMat_MomentPhys(Knormal=knormal,Alpha = alpha,Beta=beta,useAlphaBeta=True,userInputStiffness=True)],
                [Law2_SCG_MomentPhys_CohesionlessMomentRotation()]
        )

        # here we redifine the triaxial engine
        triax=O.engines[4] #redifine 3d triaxial engine
        ##We move to deviatoric loading, let us turn internal compaction off to keep particles sizes constant
        triax.internalCompaction=False

        ## Change contact friction (remember that decreasing it would generate instantaneous instabilities)
        triax.setContactProperties(fricFinal)

   this statement "i_interactions.Eta = eyta" doesn't change eta in interaction physics, but I think it should change.

     Here we have parameters in this function below, doesn't this function change kn, ks beta in interaction physics?

             [Ip2_MomentMat_MomentMat_MomentPhys(Knormal=knormal,Alpha = alpha,Beta=beta,useAlphaBeta=True,userInputStiffness=True)],

Hello,

a full script would be really helpful in finding a mistake :-) as now we do
not know exactly what some details are..

       # here we modify the material and interaction physics
> for i_body in O.bodies:
> if i_body.id < 6 :
> i_body.material = O.materials[3]
> else:
> i_body.material = O.materials[2]
> for i_interactions in O.interactions:
> if i_interactions.id1 > 5 and i_interactions.id2 > 5:
> i_interactions.Eta = eyta # only interaction
> between particles are changed to eyta > 0
>
> Does these interactions exists at all? anyway, you should assign Eta to
i_interaction.phys, not i_interaction itself..
i_interaction.phys.Eta = eyta # [1]

>
> #here we move to modify the interaction loop
> Interaction = O.engines[2]
> Interaction = InteractionLoop(
> [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
> [Ip2_MomentMat_MomentMat_MomentPhys(Knormal=knormal,Alpha =
> alpha,Beta=beta,useAlphaBeta=True,userInputStiffness=True)],
> [Law2_SCG_MomentPhys_CohesionlessMomentRotation()]
> )
>
> actually this modifies nothing at all.. The fist line (Interaction =
O.engines[2]) creates temporary variable Interaction, which is on the next
line (Interaction = ..) overwritten with newly created InteractionLoop,
which is not used (according to the part of the script) anywhere

> # here we redifine the triaxial engine
> triax=O.engines[4] #redifine 3d triaxial engine
> ##We move to deviatoric loading, let us turn internal compaction
> off to keep particles sizes constant
> triax.internalCompaction=False
>
> ## Change contact friction (remember that decreasing it would
> generate instantaneous instabilities)
> triax.setContactProperties(fricFinal)
>
>
> this statement "i_interactions.Eta = eyta" doesn't change eta in
> interaction physics, but I think it should change.
>

explained above [1]..

> Here we have parameters in this function below, doesn't this function
> change kn, ks beta in interaction physics?
>
> [Ip2_MomentMat_MomentMat_MomentPhys(Knormal=knormal,Alpha =
> alpha,Beta=beta,useAlphaBeta=True,userInputStiffness=True)],
>
> no, the InteractionLoop is not used enywhere.. if you want to redefined
O.engines, you must reassigned it, see [2]. Example:

newIntLoop = InteractionLoop(...)
newEngines = [
   ..., # no change from previous steps
   newIntLoop, # only this engine is new
   ..., # no changes
]
O.engines = newEngines # to make new engines work, they *must* be
explicitely assigned like here, see [2]

Jan

[2] https://www.yade-dem.org/doc/user.html#base-engines

well, I have three parts,

in part 1 and part 2, I generate a sample with specific relative density and then compress to 100kPa isotropicly, which generated a xml file to be loaded in my triaxial test.

In part 3, I start my deviatoric load process to simulate the drained triaxial test process. My code is below

nRead=utils.readParamsFromTable(
                knormal=170e6,
                alpha = 0.1,
                fricFinal = 30,
                beta = 1.8,
                eyta = 1.8,
                unknow=True
)
from yade.params import table

knormal = table.knormal
alpha = table.alpha
fricFinal = table.fricFinal
beta = table.beta
eyta = table.eyta

young = 6e6
poisson = 0.1
num_spheres = 10000
key = 'Ec'+str(young)+'poisson'+str(poisson)+'num'+str(num_spheres)
key_output = 'knormal_'+str(knormal)+'alpha_'+str(alpha)+'fric_'+str(fricFinal)+'beta_'+str(beta)+'eta_'+str(eyta)
rate = 0.02
damp=0.2 # damping coefficient
stabilityThreshold=0.01

 O.load('consolidated_state_'+key+'_'+str(phi)+'.xml')

# we add new materials to our material container
        O.materials.append(MomentMat(eta=eyta,young=young,poisson=0.1,frictionAngle=radians(fricFinal),density=2600,label='spheres-moment'))
        O.materials.append(MomentMat(eta=0.0,young=young,poisson=0.1,frictionAngle=0,density=0,label='walls-moment'))

        # here we modify the material and interaction physics
        for i_body in O.bodies:
                if i_body.id < 6 :
                        i_body.material = O.materials[3]
                else:
                        i_body.material = O.materials[2]
        for i_interactions in O.interactions:
                if i_interactions.id1 > 5 and i_interactions.id2 > 5:
                        i_interactions.phys.Eta = eyta # only interaction between particles are changed to eyta > 0
#here we move to modify the interaction loop
        Interaction = O.engines[2]
        Interaction = InteractionLoop(
                [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
                [Ip2_MomentMat_MomentMat_MomentPhys(Knormal=knormal,Alpha = alpha,Beta=beta,useAlphaBeta=True,userInputStiffness=True)],
                [Law2_SCG_MomentPhys_CohesionlessMomentRotation()]
        )

        # here we redifine the triaxial engine
        triax=O.engines[4] #redifine 3d triaxial engine
        ##We move to deviatoric loading, let us turn internal compaction off to keep particles sizes constant
        triax.internalCompaction=False

        ## Change contact friction (remember that decreasing it would generate instantaneous instabilities)
        triax.setContactProperties(fricFinal)

        ##set independant stress control on each axis
        triax.stressControl_1=triax.stressControl_2=triax.stressControl_3=True
        ## We turn all these flags true, else boundaries will be fixed
        triax.wall_bottom_activated=True
        triax.wall_top_activated=True
        triax.wall_left_activated=True
        triax.wall_right_activated=True
        triax.wall_back_activated=True
        triax.wall_front_activated=True

        ##If we want a triaxial loading at imposed strain rate, let's assign srain rate instead of stress
        triax.stressControl_2=0 #we are tired of typing "True" and "False", we use implicit conversion from integer to boolean
        triax.strainRate2=rate
        triax.strainRate1=100*rate
        triax.strainRate3=100*rate

        ##we can change damping here. What is the effect in your opinion?
        #newton.damping=0.1

        ##Save temporary state in live memory. This state will be reloaded from the interface with the "reload" button.
        O.saveTmp()

#here we redefine the history recorder
        from yade import plot

### a function saving variables
        def history():
                plot.addData(e11=triax.strain[0], e22=triax.strain[1], e33=triax.strain[2],
                    s11=triax.stress(triax.wall_right_id)[0],
                    s22=triax.stress(triax.wall_top_id)[1],
                    s33=triax.stress(triax.wall_front_id)[2],
                    i=O.iter,CN = utils.avgNumInteractions())
        O.engines[4]=PyRunner(iterPeriod=1500,command='history()',label='recorder')

        while 1:
                O.run(1000,True)
                if O.iter%100 == 0:
                        print "iter = ", O.iter
                        print "Syy = ",triax.strain[1]
                if triax.strain[1] > 0.3:
                        break

        ### declare what is to plot. "None" is for separating y and y2 axis
        plot.plots={'i':('e11','e22','e33',None,'s11','s22','s33','CN')}

##display on the screen (doesn't work on VMware image it seems)
##plot.plot()

## In that case we can still save the data to a text file at the the end of the simulation, with:
        O.save('Final'+key_output+'.xml')
        plot.saveDataTxt('shear'+key_output+'_'+str(phi)+'.txt')

For such kind of problems I would _strongly_ recommend to use just 2-3
spheres (if it is possible) to understand, where the problem comes
from.

In these complex scripts there are a lot of "opportunities" to loose
something and low chance to find the problem place.

Anton

Anyway I would start with

i_interaction.phys.Eta = ...

instead of

i_interaction.Eta = ...

Jan
Dne 19.5.2012 18:26 "Anton Gladky" <email address hidden>
napsal(a):

> Question #197661 on Yade changed:
> https://answers.launchpad.net/yade/+question/197661
>
> Status: Open => Answered
>
> Anton Gladky proposed the following answer:
> For such kind of problems I would _strongly_ recommend to use just 2-3
> spheres (if it is possible) to understand, where the problem comes
> from.
>
> In these complex scripts there are a lot of "opportunities" to loose
> something and low chance to find the problem place.
>
> Anton
>
> --
> 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
>