no interactions for JCFpMat() under JCFpmPhys() engines

> no interations were generated

please provide a way how you obtained this.

I tested your code and have ~6000 interactions.
print(len([i for i in O.interactions]))

unbalancedForce() is nan, because it is computed as "ratio of ... force on bodies and mean force magnitude on interactions." [1].
However, the interaction forces are all zero (thus their mean, too) and nan is result of zero division, not not generated interactions...

cheers
Jan

[1] https://yade-dem.org/doc/yade.utils.html#yade._utils.unbalancedForce

Dear Jan,
Thanks for your reply.

I used to check the number of the interactions by command line and also found its nonzero.

However, when I click ‘inspect’ in GUI invoked by yade.qt.controller(), and select ‘interactions’, errors will appear in the command line and indicate no interactions exist.

If only the force is zero, the unbalanced force should be zero, but not NaN. The infinite value may come from the zero divisor, i.e. the number of the interactions.

Basically, in fact since there are already some overlaps between the spheres in the above model, the force should also not be zero. When I change the O.engines to be the second group as above, both the force and number of interactions are nonzero.

Would you please help me continue to identify this problem? Thanks a lot.
Yang

> in GUI

I don't use Yade GUI, so I cannot help here.. just have tried and indeed it behaves strangely.. You can open an issue [2] if you like

> If only the force is zero, the unbalanced force should be zero, but not NaN. The infinite value may come from the zero divisor, i.e. the number of the interactions.

unbalancedForce() is defined as bodyForce / interactionForce. If interactionForce is zero, then I see no reason why the result should be zero.. NaN makes total sense to me, even in the case 0/0.
You can consider NaN as a special case - zero interaction force, no interactions (pointed by you if I got it correctly) etc.
Of course, you can make up different definitions of unbalanced force more suitable for a specific situation.

> Basically, in fact since there are already some overlaps between the spheres in the above model, the force should also not be zero.

Some models are designed to specific situation. I have experience with CpmMat being used to model a continuum, where initial "overlaps" (in that case rather gaps) are considered as stress-free / force-free state.
I guess JCFpmMat does the same.

> When I change the O.engines to be the second group as above, both the force and number of interactions are nonzero.

because Frict* model just takes penetration depth and compute forces. Now they are nonzero

cheers
Jan

[2] https://gitlab.com/yade-dev/trunk/issues

Hi Yang,

What do you want to achieve?

It seems that there is no loading applied to your sample. No loading means that the forces are all null...

Please try the following script and tell me if it works (uniaxial compression with the JCFPM).

Luc

###

from yade import pack, plot

################# SIMULATIONS DEFINED HERO.bodies.append(pack.randomDensePack(pred,radius=D/30.,rRelFuzz=0.4,spheresInCell=1000,memoizeDb='/tmp/gts-triax-packings.sqlite',returnSpherePack=False,color=(0.9,0.8,0.6),material=sphereMat)) # RK: radius will determine the nb of particles!
E

#### packing (previously constructed)
OUT='compressionTest'

#### Simulation Control
rate=-0.05 #deformation rate
iterMax=200000 # maximum number of iterations
saveVTK=int(iterMax/5.) # saving output files for paraview

#### Microproperties (interparticle parameters)
intR=1.4 # allows near neighbour interaction (can be adjusted for every packing / the bigger -> the more brittle / careful when intR is too large -> bonds can be created "over" particles) -> intR can be calibrated to reach a certain coordination number K (see calculation on line 115)
DENS=3000 # this one can be adjusted for different reasons (porosity of packing vs porosity of material / increase time step (no gravity -> no real effect on the result)
YOUNG=10e9 # this one controls the Young's modulus of the material
ALPHA=0.15 # this one controls the material Poisson's ratio of the material
TENS=3e6 # this one controls the tensile strength UTS of the material
COH=30e6 # this one controls the compressive strength UCS of the material, more precisely, the ratio UCS/UTS (from my experience: COH should be >= to TENS, >= 10*TENS for competent materials like concrete)
FRICT=30 # this one controls the slope of the failure envelope (effect mainly visible on triaxial compression tests)

#### material definition
def sphereMat(): return JCFpmMat(type=1,density=DENS,young=YOUNG,poisson=ALPHA,tensileStrength=TENS,cohesion=COH,frictionAngle=radians(FRICT))

#### create the specimen
L=0.10
D=0.05
pred=pack.inCylinder((0,0,0),(0,0,L),D/2.)
#O.bodies.append(pack.regularHexa(pred,radius=D/20.,gap=0.,material=sphereMat)) # regular packings should be avoided as failure is controlled by particle arrangement
O.bodies.append(pack.randomDensePack(pred,radius=D/30.,rRelFuzz=0.4,spheresInCell=1000,memoizeDb='/tmp/gts-triax-packings.sqlite',returnSpherePack=False,color=(0.9,0.8,0.6),material=sphereMat)) # RK: radius will determine the nb of particles!

R=0
Rmax=0
Rmin=1e6
nbSpheres=0.
for o in O.bodies:
 if isinstance(o.shape,Sphere):
   nbSpheres+=1
   R+=o.shape.radius
   if o.shape.radius>Rmax:
     Rmax=o.shape.radius
   if o.shape.radius<Rmin:
     Rmin=o.shape.radius
Rmean=R/nbSpheres

print('nbSpheres=',nbSpheres,' | Rmean=',Rmean, ' | Rmax/Rmin=', Rmax/Rmin)

#### help define boundary conditions (see utils.uniaxialTestFeatures)
bb=utils.uniaxialTestFeatures()
negIds,posIds,longerAxis,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area']

################# DEM loop + ENGINES DEFINED HERE

O.engines=[
 ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intR,label='Saabb')]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intR,label='SSgeom')],
  [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,label='interactionPhys')],
  [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=True,Key=OUT,label='interactionLaw')]
 ),
 UniaxialStrainer(strainRate=rate,axis=longerAxis,asymmetry=0,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=1,blockRotations=1,setSpeeds=0,stopStrain=0.1,dead=1,label='strainer'),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.8,defaultDt=utils.PWaveTimeStep()),
 NewtonIntegrator(damping=0.4,label='newton'),
 PyRunner(iterPeriod=int(100),initRun=True,command='recorder()',label='data'),
        VTKRecorder(iterPeriod=1,initRun=True,fileName=OUT+'-',recorders=['spheres','jcfpm','cracks','bstresses'],Key=OUT,dead=1,label='vtk')
]

################# RECORDER DEFINED HERE

def recorder():
    yade.plot.addData({'i':O.iter,
         'eps':strainer.strain,
         'sigma':strainer.avgStress,
        'tc':interactionLaw.nbTensCracks,
        'sc':interactionLaw.nbShearCracks,
        #'te':interactionLaw.totalTensCracksE,
        #'se':interactionLaw.totalShearCracksE,
         'unbF':utils.unbalancedForce()})
    plot.saveDataTxt(OUT)

################# PREPROCESSING

#### manage interaction detection factor during the first timestep and then set default interaction range
O.step();
### initializes the interaction detection factor
SSgeom.interactionDetectionFactor=-1.
Saabb.aabbEnlargeFactor=-1.

#### remove particles with less than Kmin contacts (avoid "floating particles") -> Kmin can be adjusted
Kmin=3
erase=0
for o in O.bodies:
    if (not o) or (o.id in negIds) or (o.id in posIds) : continue
    cont=0
    for i in O.interactions.withBody(o.id) :
 if not i.isReal : continue
 if i.phys.isCohesive :
     cont+=1
    if cont<Kmin :
 erase+=1
 O.bodies.erase(o.id)
 #print erase, ' | id=', o.id, ' | cont=', cont
print("we removed ", erase, " particles with less than ", Kmin, " contacts")

#### coordination number calculation
numSSlinks=0
numCohesivelinks=0
for i in O.interactions:
    if not i.isReal : continue
    if isinstance(O.bodies[i.id1].shape,Sphere) and isinstance(O.bodies[i.id2].shape,Sphere):
      numSSlinks+=1
    if i.phys.isCohesive :
      numCohesivelinks+=1
print ("Kmean=", 2.0*numCohesivelinks/nbSpheres)

vtk.dead=0
O.step()

################# SIMULATION REALLY STARTS HERE
strainer.dead=0
vtk.iterPeriod=saveVTK
O.run(iterMax)

Dear Jan and luc,

 Thank you for your information. The problem has been solved. i.e. as you notified:

>'No loading means that the forces are all null' and

>'Some models are designed to specific situation. I have experience with CpmMat being used to model a continuum, where initial "overlaps" (in that case rather gaps) are considered as stress-free / force-free state.
I guess JCFpmMat does the same.'

I wants to establish an underground model with faults and the first step is to establish an initial pack with homogeneous stress state within a confining box. To achieve this,

Firstly, I use the makeCloud command to generate loose spheres and extend their radius to make the pack be dense and have some original overlaps between the spheres.

Then, I plan to make calculation without external loading to adjust the original overlaps automatically and reach an equilibrium state.

During this step, I found that no force will be found with the JCFpm engine group, but with the Frictm engine group, the spheres can adjust their positions according to the unbalanced forces acting on them during the calculation.

So, in the future work, I will invoke the JCFpm engine group only after some loading is applied, and use the Frictm engine group to adjust the initial pack.

However, even the loading is applied, errors will still appear when check the interactions in GUI. I think this error can be put aside and will not affect the calculation.

Best regards
Yang