Engine Error (JCFpmMat)

Hi Sam,
do you know at which iteration number the problem occur?
Jan

Dear Jan,

To be honest, no I do not know.

However, I am sure this error because of my wall material. This is FricMat and I need to consider it in Ip2 and Law2. But I do not know how can I bring it to Law 2.

Regards
Sam

I had a look to the source code and could not imagine a situation to trigger this situation.. Could you post a full code producing the error?
thanks
Jan

Thanks for the code in personal communication. However, next time please post it here, too.
I suspect O.bodies.replaceByClumps..

What happen if you O.materials.append first FrictMat and then JCFpmMat (updating material of walls and spheres accordingly)? Do you still have the same problem?

For the replaceByClumps issue, could you please try to test a very simple example (I have never used the function myself and don't know to use it) with:
1) O.materials.append(JCFpmMat(...))
2) O.materials.append(FrictMat(...))
3) create a few spheres
4) replaceByClumps
5) O.step()

the code crashes at the very first step.

thanks
Jan

Hi,

Did you consider also using JCFpmMat (another one, with eg zero friction) for the walls ?

related to https://answers.launchpad.net/yade/+question/676746

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

Daer Jan and Jerome,

Thanks for your guides but I think it should be a bug in this contact model.

So, I think it should be modified from the source.

Regards
Sam

Hello,
Please fill a bug report if relevant. It has more chance to trigger someone than calling it bug in a comment.
It could help to answer the questions asked in #4 and #5.
Bruno

Dear Bruno,

I tried to run a triaxial simulation based on the JCFpmMat with a simple script but still is not working properly.

Please find here in my script:

from yade import pack,plot
from yade import export, ymport
import random
from yade import utils
import numpy as np

############################################
### DEFINING VARIABLES AND MATERIALS ###
############################################

# The following 5 lines will be used later for batch execution
nRead=readParamsFromTable(
 num_spheres=1000,# number of spheres
 compFricDegree =30, # contact friction during the confining phase
 key='_TriaxJ2_', # put you simulation's name here
 unknownOk=True
)
from yade.params import table

num_spheres=table.num_spheres# number of spheres
key=table.key
targetPorosity = 0.5 #the porosity we want for the packing
frictionAngle = 30# contact friction during the deviatoric loading
finalFricDegree = 30
rate=-0.02 # loading rate (strain rate)
damp=0.1 # Newton damping
stabilityThreshold=0.01 # we test unbalancedForce against this value in different loops (see below)
iterper=1000
young=4.5e5# contact stiffness
poisson=0.25
cohesion=20e6
tensileStrength = 7e6
density = 1950
intRadius= 1

mn,mx=Vector3(0,0,0),Vector3(.1,.1,.1) # corners of the initial packing

O.materials.append(JCFpmMat(young=young, cohesion=cohesion, density=density, frictionAngle=radians(frictionAngle), tensileStrength=tensileStrength, poisson=poisson,label='spheres'))

#frictionless walls
O.materials.append(FrictMat(young=80e9,poisson=.45,frictionAngle=radians(frictionAngle),density=7000,label='walls'))

## create walls around the packing
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

## use a SpherePack object to generate a random loose particles packing
sp=pack.SpherePack()

###########################
## DEFINING ENGINES ###
###########################

triax=TriaxialStressController(
 ## TriaxialStressController will be used to control stress and strain. It controls particles size and plates positions.
 ## this control of boundary conditions was used for instance in http://dx.doi.org/10.1016/j.ijengsci.2008.07.002
 maxMultiplier=1.+2e4/young, # spheres growing factor (fast growth)
 finalMaxMultiplier=1.+2e3/young, # spheres growing factor (slow growth)
 thickness = 0,
 stressMask = 7,
 internalCompaction=False, # If true the confining pressure is generated by growing particles
)

newton=NewtonIntegrator(damping=damp, gravity=[0,0,0])
########################################
#Modified engine
##################################
O.engines=[
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intRadius),Bo1_Box_Aabb()]),
        InteractionLoop(
                [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intRadius),Ig2_Box_Sphere_ScGeom()],
                [Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1, label='jcf')],
  #[Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,label="cohesiveIp")],
                [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(smoothJoint=True,label='interactionLaw', recordCracks=True),Law2_ScGeom_FrictPhys_CundallStrack()]
  #[Law2_ScGeom_CohFrictPhys_CohesionMoment(
   #useIncrementalForm=True, #useIncrementalForm is turned on as we want plasticity on the contact moments
   #always_use_moment_law=True, #if we want "rolling" friction even if the contact is not cohesive (or cohesion is broken), we will have to turn this true somewhere
   #label='cohesiveLaw')]
        ),
        ## We will use the global stiffness of each body to determine an optimal timestep (see https://yade-dem.org/w/images/1/1b/Chareyre&Villard2005_licensed.pdf)
        GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.05),
        triax,
  #VTKRecorder(iterPeriod=100,initRun=True,fileName='VTKFiles/'+'-',recorders=['all']),
        TriaxialStateRecorder(iterPeriod=100,file='Wall stress,'+key),
        newton
]

##########################################################
#Display spheres with 2 colors for seeing rotations better
Gl1_Sphere.stripes=0
if nRead==0: yade.qt.Controller(), yade.qt.View()

## UNCOMMENT THE FOLLOWING SECTIONS ONE BY ONE
## DEPENDING ON YOUR EDITOR, IT COULD BE DONE
## BY SELECTING THE CODE BLOCKS BETWEEN THE SUBTITLES
## AND PRESSING CTRL+SHIFT+D
#if nRead==0: yade.qt.Controller(), yade.qt.View()
print 'Number of elements: ', len(O.bodies)
print 'Box Volume: ', triax.boxVolume
#######################################
### APPLYING CONFINING PRESSURE ###
#######################################

#the value of (isotropic) confining stress defines the target stress to be applied in all three directions
triax.goal1=triax.goal2=triax.goal3=-150000

###################################################
### REACHING A SPECIFIED POROSITY PRECISELY ###
###################################################

### We will reach a prescribed value of porosity with the REFD algorithm
### (see http://dx.doi.org/10.2516/ogst/2012032 and
### http://www.geosyntheticssociety.org/Resources/Archive/GI/src/V9I2/GI-V9-N2-Paper1.pdf)

while triax.porosity>targetPorosity:
 ## we decrease friction value and apply it to all the bodies and contacts
 frictionAngle = 0.95*frictionAngle
 setContactFriction(radians(frictionAngle))
 print "\r Friction: ",frictionAngle," porosity:",triax.porosity,
 sys.stdout.flush()
 ## while we run steps, triax will tend to grow particles as the packing
 ## keeps shrinking as a consequence of decreasing friction. Consequently
 ## porosity will decrease
 O.run(500,1)

O.save('compactedStateJ2'+key+'.yade.gz')
print "### Compacted state saved ###"

##############################
### DEVIATORIC LOADING ###
##############################

##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)
setContactFriction(radians(finalFricDegree))

##set stress control on x and z, we will impose strain rate on y
triax.stressMask = 5
##now goal2 is the target strain rate
triax.goal2=rate
## we define the lateral stresses during the test, here the same 10kPa as for the initial confinement.
triax.goal1=-150000
triax.goal3=-150000

##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()

#####################################################
### Example of how to record and plot data ###
#####################################################

#from yade import plot

### a function saving variables
def history():
   plot.addData(e11=-triax.strain[0], e22=-triax.strain[1], e33=-triax.strain[2],
        ev=-triax.strain[0]-triax.strain[1]-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)

if 1:
  ## include a periodic engine calling that function in the simulation loop
  O.engines=O.engines[0:5]+[PyRunner(iterPeriod=20,command='history()',label='recorder')]+O.engines[5:7]
  ##O.engines.insert(4,PyRunner(iterPeriod=20,command='history()',label='recorder'))
else:
  ## With the line above, we are recording some variables twice. We could in fact replace the previous
  ## TriaxialRecorder
  ## by our periodic engine. Uncomment the following line:
  O.engines[4]=PyRunner(iterPeriod=20,command='history()',label='recorder')

  #NOTE, that after replacing some overlaps may occur.
#So after replacing calm() function may be helpful:
#if 0:
 #print '\nPlease wait a minute ...\n'
 #O.engines=O.engines+[PyRunner(iterPeriod=10000,command='calm()',label='calmRunner')]
 #O.run(1000000,True)
 #calmRunner.dead=True

O.run(100,True)

### declare what is to plot. "None" is for separating y and y2 axis
#plot.plots={'i':('e11','e22','e33',None,'s11','s22','s33')}
### the traditional triaxial curves would be more like this:
#plot.plots={'e22':('s11','s22','s33',None,'ev')}
plot.plots={'e22':('s11','s22')}

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

##### PLAY THE SIMULATION HERE WITH "PLAY" BUTTON OR WITH THE COMMAND O.run(N) #####

## In that case we can still save the data to a text file at the the end of the simulation, with:
plot.saveDataTxt('resultsJ2'+key)
##or even generate a script for gnuplot. Open another terminal and type "gnuplot plotScriptKEY.gnuplot:
plot.saveGnuplot('plotScriptJ2'+key)
rr=yade.qt.Renderer()
rr.shape=True
rr.intrPhys=False

Regards
Hessam

Hi, sorry to hear that you still have a problem. I don't understand why
you are sending another script.
Fixing the issue would need your cooperation and consistent replies to our
questions, ideally in a single thread.
At the moment it looks more like spamming - sorry.
Bruno

Le mar. 25 déc. 2018 03:08, SayedHessam <
<email address hidden>> a écrit :

> Question #676382 on Yade changed:
> https://answers.launchpad.net/yade/+question/676382
>
> SayedHessam posted a new comment:
> Dear Bruno,
>
> I tried to run a triaxial simulation based on the JCFpmMat with a simple
> script but still is not working properly.
>
> Please find here in my script:
>
>
> from yade import pack,plot
> from yade import export, ymport
> import random
> from yade import utils
> import numpy as np
>
> ############################################
> ### DEFINING VARIABLES AND MATERIALS ###
> ############################################
>
> # The following 5 lines will be used later for batch execution
> nRead=readParamsFromTable(
> num_spheres=1000,# number of spheres
> compFricDegree =30, # contact friction during the confining phase
> key='_TriaxJ2_', # put you simulation's name here
> unknownOk=True
> )
> from yade.params import table
>
> num_spheres=table.num_spheres# number of spheres
> key=table.key
> targetPorosity = 0.5 #the porosity we want for the packing
> frictionAngle = 30# contact friction during the deviatoric loading
> finalFricDegree = 30
> rate=-0.02 # loading rate (strain rate)
> damp=0.1 # Newton damping
> stabilityThreshold=0.01 # we test unbalancedForce against this value in
> different loops (see below)
> iterper=1000
> young=4.5e5# contact stiffness
> poisson=0.25
> cohesion=20e6
> tensileStrength = 7e6
> density = 1950
> intRadius= 1
>
>
> mn,mx=Vector3(0,0,0),Vector3(.1,.1,.1) # corners of the initial packing
>
>
> O.materials.append(JCFpmMat(young=young, cohesion=cohesion,
> density=density, frictionAngle=radians(frictionAngle),
> tensileStrength=tensileStrength, poisson=poisson,label='spheres'))
>
> #frictionless walls
>
> O.materials.append(FrictMat(young=80e9,poisson=.45,frictionAngle=radians(frictionAngle),density=7000,label='walls'))
>
> ## create walls around the packing
> walls=aabbWalls([mn,mx],thickness=0,material='walls')
> wallIds=O.bodies.append(walls)
>
> ## use a SpherePack object to generate a random loose particles packing
> sp=pack.SpherePack()
>
>
> ###########################
> ## DEFINING ENGINES ###
> ###########################
>
> triax=TriaxialStressController(
> ## TriaxialStressController will be used to control stress and strain. It
> controls particles size and plates positions.
> ## this control of boundary conditions was used for instance in
> http://dx.doi.org/10.1016/j.ijengsci.2008.07.002
> maxMultiplier=1.+2e4/young, # spheres growing factor (fast growth)
> finalMaxMultiplier=1.+2e3/young, # spheres growing factor (slow growth)
> thickness = 0,
> stressMask = 7,
> internalCompaction=False, # If true the confining pressure is generated
> by growing particles
> )
>
> newton=NewtonIntegrator(damping=damp, gravity=[0,0,0])
> ########################################
> #Modified engine
> ##################################
> O.engines=[
> ForceResetter(),
>
> InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intRadius),Bo1_Box_Aabb()]),
> InteractionLoop(
>
> [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intRadius),Ig2_Box_Sphere_ScGeom()],
>
> [Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,
> label='jcf')],
>
> #[Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,label="cohesiveIp")],
>
> [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(smoothJoint=True,label='interactionLaw',
> recordCracks=True),Law2_ScGeom_FrictPhys_CundallStrack()]
> #[Law2_ScGeom_CohFrictPhys_CohesionMoment(
> #useIncrementalForm=True, #useIncrementalForm is turned on as we want
> plasticity on the contact moments
> #always_use_moment_law=True, #if we want "rolling" friction even if the
> contact is not cohesive (or cohesion is broken), we will have to turn this
> true somewhere
> #label='cohesiveLaw')]
> ),
> ## We will use the global stiffness of each body to determine an
> optimal timestep (see
> https://yade-dem.org/w/images/1/1b/Chareyre&Villard2005_licensed.pdf)
>
> GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.05),
> triax,
>
> #VTKRecorder(iterPeriod=100,initRun=True,fileName='VTKFiles/'+'-',recorders=['all']),
> TriaxialStateRecorder(iterPeriod=100,file='Wall stress,'+key),
> newton
> ]
>
> ##########################################################
> #Display spheres with 2 colors for seeing rotations better
> Gl1_Sphere.stripes=0
> if nRead==0: yade.qt.Controller(), yade.qt.View()
>
> ## UNCOMMENT THE FOLLOWING SECTIONS ONE BY ONE
> ## DEPENDING ON YOUR EDITOR, IT COULD BE DONE
> ## BY SELECTING THE CODE BLOCKS BETWEEN THE SUBTITLES
> ## AND PRESSING CTRL+SHIFT+D
> #if nRead==0: yade.qt.Controller(), yade.qt.View()
> print 'Number of elements: ', len(O.bodies)
> print 'Box Volume: ', triax.boxVolume
> #######################################
> ### APPLYING CONFINING PRESSURE ###
> #######################################
>
> #the value of (isotropic) confining stress defines the target stress to be
> applied in all three directions
> triax.goal1=triax.goal2=triax.goal3=-150000
>
> ###################################################
> ### REACHING A SPECIFIED POROSITY PRECISELY ###
> ###################################################
>
> ### We will reach a prescribed value of porosity with the REFD algorithm
> ### (see http://dx.doi.org/10.2516/ogst/2012032 and
> ###
> http://www.geosyntheticssociety.org/Resources/Archive/GI/src/V9I2/GI-V9-N2-Paper1.pdf
> )
>
>
> while triax.porosity>targetPorosity:
> ## we decrease friction value and apply it to all the bodies and contacts
> frictionAngle = 0.95*frictionAngle
> setContactFriction(radians(frictionAngle))
> print "\r Friction: ",frictionAngle," porosity:",triax.porosity,
> sys.stdout.flush()
> ## while we run steps, triax will tend to grow particles as the packing
> ## keeps shrinking as a consequence of decreasing friction. Consequently
> ## porosity will decrease
> O.run(500,1)
>
> O.save('compactedStateJ2'+key+'.yade.gz')
> print "### Compacted state saved ###"
>
> ##############################
> ### DEVIATORIC LOADING ###
> ##############################
>
> ##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)
> setContactFriction(radians(finalFricDegree))
>
> ##set stress control on x and z, we will impose strain rate on y
> triax.stressMask = 5
> ##now goal2 is the target strain rate
> triax.goal2=rate
> ## we define the lateral stresses during the test, here the same 10kPa as
> for the initial confinement.
> triax.goal1=-150000
> triax.goal3=-150000
>
> ##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()
>
> #####################################################
> ### Example of how to record and plot data ###
> #####################################################
>
> #from yade import plot
>
> ### a function saving variables
> def history():
> plot.addData(e11=-triax.strain[0], e22=-triax.strain[1],
> e33=-triax.strain[2],
> ev=-triax.strain[0]-triax.strain[1]-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)
>
> if 1:
> ## include a periodic engine calling that function in the simulation loop
>
> O.engines=O.engines[0:5]+[PyRunner(iterPeriod=20,command='history()',label='recorder')]+O.engines[5:7]
>
> ##O.engines.insert(4,PyRunner(iterPeriod=20,command='history()',label='recorder'))
> else:
> ## With the line above, we are recording some variables twice. We could
> in fact replace the previous
> ## TriaxialRecorder
> ## by our periodic engine. Uncomment the following line:
> O.engines[4]=PyRunner(iterPeriod=20,command='history()',label='recorder')
>
> #NOTE, that after replacing some overlaps may occur.
> #So after replacing calm() function may be helpful:
> #if 0:
> #print '\nPlease wait a minute ...\n'
>
> #O.engines=O.engines+[PyRunner(iterPeriod=10000,command='calm()',label='calmRunner')]
> #O.run(1000000,True)
> #calmRunner.dead=True
>
> O.run(100,True)
>
> ### declare what is to plot. "None" is for separating y and y2 axis
> #plot.plots={'i':('e11','e22','e33',None,'s11','s22','s33')}
> ### the traditional triaxial curves would be more like this:
> #plot.plots={'e22':('s11','s22','s33',None,'ev')}
> plot.plots={'e22':('s11','s22')}
>
> ## display on the screen (doesn't work on VMware image it seems)
> plot.plot()
>
> ##### PLAY THE SIMULATION HERE WITH "PLAY" BUTTON OR WITH THE COMMAND
> O.run(N) #####
>
> ## In that case we can still save the data to a text file at the the end
> of the simulation, with:
> plot.saveDataTxt('resultsJ2'+key)
> ##or even generate a script for gnuplot. Open another terminal and type
> "gnuplot plotScriptKEY.gnuplot:
> plot.saveGnuplot('plotScriptJ2'+key)
> rr=yade.qt.Renderer()
> rr.shape=True
> rr.intrPhys=False
>
>
> Regards
> Hessam
>
> --
> You received this question notification because your team yade-users 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 Bruno,

You right but this not another script. I've just tried to run a simple triaxial based on JCFpmMat contact model but I could not run the model properly. Anyway, I think it is a bug from the source.

Regards
Hessam