Use Law2_ScGeom_CpmPhys_Cpm() with two spheres

Hello Hicham,

then defined the engines :
>
> UX=uniaxialTestFeatures()
>
> negIds,posIds,axis,crossSectionArea=UX['negIds'],UX['posIds'],UX['axis'],UX['area']
>

Could you try to "debug" the four values values from UX (just 'print axis'
is enough)? This is the first suspicious place..

>
>
> Question : The data plot stays empty even if the simulation shows that the
> two spheres are interacting (tension and compression tests are running
> well) and the strengths fc and and ft printed on screen show "NAN".
>
> - Is using two spheres only can show results similar to a test on a normal
> specimen (of 2000 spheres) ?
>

It might be the reason. The UniaxialStrainer is not used much AFAIK and
probably nobody used it for two spheres before :-)

> - Am I using the engines right ? I am new to YADE and that may be the
> problem
>

The engines seems to be ok..

>
> Another question not related to this topic : Are all the parameters of a
> Law hard coded or is it possible to modify them when writing the python
> script ?
>

It depends on what you mean by "Law". Som eparts of the law are determined
by the material parameters, which are defined on CpmMat material. There are
a few parameters on Law2_ScGeom_CpmPhys_Cpm that might be assigned, see [1]

In similar cases (when the script has reasonable length, i.e. < 100 lines),
please include all the script to the mail, we can then try it ourselves and
give some suggestions based on personal testing. On the other hand (for the
future), do not include too long scripts :-)

cheers
Jan

[1]
https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.Law2_ScGeom_CpmPhys_Cpm

Thanks for the quick answer Jan !

->Could you try to "debug" the four values values from UX (just 'print axis'
is enough)? This is the first suspicious place..

negIds=[0], posIds=[1], axis=2, area=0.0

for further info, here is the complete script (might be familiar to you as it is a modified version of uniax.py :-) )

from yade import plot,pack,timing,utils
import time, sys, os, copy

readParamsFromTable(
        noTableOk=True,
        young=24e9,
        poisson=.2,
        sigmaT=3.5e6,
        frictionAngle=0,
        epsCrackOnset=1e-4,
        relDuctility=30,
        intRadius=1,
        dtSafety=1,
        damping=0.4,
        strainRateTension=.05,
        strainRateCompression=.05,
        setSpeeds=True,
        doModes=3, # 1=tension, 2=compression , 3=both
        sphereRadius=3.5e-3,
        isoPrestress=0, # isotropic confinement
)

from yade.params.table import *

if 'description' in O.tags.keys(): O.tags['id']=O.tags['id']+O.tags['description']

concreteId=O.materials.append(CpmMat(young=young,frictionAngle=frictionAngle,poisson=poisson,density=4800,sigmaT=sigmaT,relDuctility=relDuctility,epsCrackOnset=epsCrackOnset,isoPrestress=isoPrestress))

O.bodies.append([
   utils.sphere(center=(0,0,-1.5*sphereRadius),radius=sphereRadius,material=concreteId,fixed=True),
   utils.sphere((0,0,1.5*sphereRadius),radius=sphereRadius,material=concreteId)
])

yade.qt.View()

UX=uniaxialTestFeatures()
negIds,posIds,axis,crossSectionArea=UX['negIds'],UX['posIds'],UX['axis'],UX['area']
O.dt=dtSafety*PWaveTimeStep()
print 'Timestep',O.dt

O.engines=[
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intRadius,label='is2aabb')],verletDist=.5*sphereRadius),
        InteractionLoop(
                [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intRadius,label='ss2sc')],
                [Ip2_CpmMat_CpmMat_CpmPhys()],
                [Law2_ScGeom_CpmPhys_Cpm()],
        ),
        NewtonIntegrator(damping=damping,label='damper'),
        CpmStateUpdater(realPeriod=.5),

        UniaxialStrainer(strainRate=strainRateTension,axis=axis,asymmetry=0,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=False,blockRotations=False,setSpeeds=setSpeeds,label='strainer'),

        PyRunner(virtPeriod=1e-6/strainRateTension,realPeriod=1,command= 'addPlotData()',label='plotDataCollector',initRun=True),
        PyRunner(realPeriod=4,command='stopIfDamaged()',label='damageChecker'),
]

plot.plots={'eps':('sigma',)}

O.saveTmp('initial');

O.timingEnabled=False

global mode
mode='tension' if doModes & 1 else 'compression'

def initTest():
        global mode
        print "init"
        if O.iter>0:
                O.wait();
                O.loadTmp('initial')
                print "Reversing plot data"; plot.reverseData()
        else: plot.plot(subPlots=False)
        strainer.strainRate=abs(strainRateTension) if mode=='tension' else -abs(strainRateCompression)
        try:
                from yade import qt
                renderer=qt.Renderer()
                renderer.dispScale=(1000,1000,1000) if mode=='tension' else (100,100,100)
        except ImportError: pass
        print "init done, will now run."
        O.step();
        ss2sc.interactionDetectionFactor=1.
        is2aabb.aabbEnlargeFactor=1.

        O.run()

def stopIfDamaged():
        global mode
        if O.iter<2 or not plot.data.has_key('sigma'): return
        sigma,eps=plot.data['sigma'],plot.data['eps']
        extremum=max(sigma) if (strainer.strainRate>0) else min(sigma)
        minMaxRatio=0.5 if mode=='tension' else 1
        if extremum==0: return
        import sys; sys.stdout.flush()
        if abs(sigma[-1]/extremum)<minMaxRatio or abs(strainer.strain)>(5e-3 if isoPrestress==0 else 5e-2):
                if mode=='tension' and doModes & 2:
                        mode='compression'
                        print "Damaged, switching to compression... "; O.pause()
                        import thread; thread.start_new_thread(initTest,())
                        return
                else:
                        print "Damaged, stopping."
                        ft,fc=max(sigma),min(sigma)
                        print 'Strengths fc=%g, ft=%g, |fc/ft|=%g'%(fc,ft,abs(fc/ft))
                        title=O.tags['description'] if 'description' in O.tags.keys() else O.tags['params']
                        print 'gnuplot',plot.saveGnuplot(O.tags['id'],title=title)
                        print 'Bye.'
                        O.pause()

def addPlotData():
        yade.plot.addData({'t':O.time,'i':O.iter,'eps':strainer.strain,'sigma':strainer.avgStress+isoPrestress,
                })

plot.plot(subPlots=False)
initTest()
waitIfBatch()

Cheers !

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