← Back to team overview

yade-users team mailing list archive

Re: [Question #238289]: Use Law2_ScGeom_CpmPhys_Cpm() with two spheres

 

Question #238289 on Yade changed:
https://answers.launchpad.net/yade/+question/238289

Hicham BENNIOU posted a new comment:
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 !

-- 
You received this question notification because you are a member of
yade-users, which is an answer contact for Yade.