I try to use uniaxialstrainer to model a uniaxial compression test. And my script is written referred to https://github.com/yade/trunk/blob/master/examples/concrete/uniax.py.
My specific script is this:
----------------------------------------------------------------------------------------------------------------------
O.reset()
from __future__ import division
from yade import plot, pack, timing
import time, sys, os, copy
#matplotlib.rc('text',usetex=True)
#matplotlib.rc('text.latex',preamble=r'\usepackage{concrete}\usepackage{euler}')
runGunplot=False
runInGui=True
globals()['runGunplot']=locals()['runGunplot']
globals()['runInGui']=locals()['runInGui']
density=2350
young=1e8
poisson=.31
sigmaT=5e3
frictionAngle=radians(41.7)
epsCrackOnset=5e-1
relDuctility=30e0
rRelFuzz=0.1
globals()['density']=locals()['density']
globals()['young']=locals()['young']
globals()['poisson']=locals()['poisson']
globals()['sigmaT']=locals()['sigmaT']
globals()['frictionAngle']=locals()['frictionAngle']
globals()['epsCrackOnset']=locals()['epsCrackOnset']
globals()['relDuctility']=locals()['relDuctility']
globals()['rRelFuzz']=locals()['rRelFuzz']
intRadius=1.5
dtSafety=.8
damping=0.4
strainRateTension=.08
strainRateCompression=.5
setSpeeds=True
globals()['intRadius']=locals()['intRadius']
globals()['dtSafety']=locals()['dtSafety']
globals()['damping']=locals()['damping']
globals()['strainRateTension']=locals()['strainRateTension']
globals()['strainRateCompression']=locals()['strainRateCompression']
globals()['setSpeeds']=locals()['setSpeeds']
doModes=3
globals()['doModes']=locals()['doModes']
specimenLength=.1
specimenwidth=.05
specimenwidth2=specimenwidth/2.0
sphereRadius=1.7e-3
globals()['specimenLength']=locals()['specimenLength']
globals()['specimenwidth']=locals()['specimenwidth']
globals()['specimenwidth2']=locals()['specimenwidth2']
globals()['sphereRadius']=locals()['sphereRadius']
isoPrestress=0
globals()['isoPrestress']=locals()['isoPrestress']
concreteId=O.materials.append(CpmMat(density=density,young=young,frictionAngle=frictionAngle,poisson=poisson,sigmaT=sigmaT,relDuctility=relDuctility,epsCrackOnset=epsCrackOnset,isoPrestress=isoPrestress))
globals()['concreteId']=locals()['concreteId']
sps=SpherePack()
pred=pack.inCylinder((0,0,0),(0,0,specimenLength),specimenwidth2)
sp=pack.randomDensePack(pred,spheresInCell=2000,rRelFuzz=rRelFuzz,radius=sphereRadius,memoizeDb='/home/dj/yade/tmp/Uniaxialyuan.sqlite',returnSpherePack=True)
sp.toSimulation(material=concreteId)
bb=uniaxialTestFeatures()
negIds,posIds,axis,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area']
O.dt=dtSafety*PWaveTimeStep()
print 'Timestep',O.dt
globals()['negIds']=locals()['negIds']
globals()['posIds']=locals()['posIds']
globals()['axis']=locals()['axis']
globals()['crossSectionArea']=locals()['crossSectionArea']
mm,mx=[pt[axis] for pt in aabbExtrema()]
globals()['mm']=locals()['mm']
globals()['mx']=locals()['mx']
coord_00,coord_25,coord_50,coord_75,coord_100=mm+.0*(mx-mm),mm+.25*(mx-mm),mm+.5*(mx-mm),mm+.75*(mx-mm),mm+1.00*(mx-mm)
area_00,area_25,area_50,area_75,area_100=approxSectionArea(coord_00,axis),approxSectionArea(coord_25,axis),approxSectionArea(coord_50,axis),approxSectionArea(coord_75,axis),approxSectionArea(coord_100,axis)
globals()['coord_00']=locals()['coord_00']
globals()['coord_25']=locals()['coord_25']
globals()['coord_50']=locals()['coord_50']
globals()['coord_75']=locals()['coord_75']
globals()['coord_100']=locals()['coord_100']
globals()['area_00']=locals()['area_00']
globals()['area_25']=locals()['area_25']
globals()['area_50']=locals()['area_50']
globals()['area_75']=locals()['area_75']
globals()['area_100']=locals()['area_100']
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intRadius,label='is2aabb'),],verletDist=.05*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),
]
plot.plots={'eps':('sigma',)}
globals()['plot']=locals()['plot']
def addPlotData():
coord=coord_50
area=area_50
yade.plot.addData({'t':O.time,'i':O.iter,'eps':strainer.strain,'sigma':strainer.avgStress+isoPrestress,
'sigma':abs(forcesOnCoordPlane(coord,axis)[axis]/area+isoPrestress), })
globals()['addPlotData']=locals()['addPlotData']
O.step()
ss2sc.interactionDetectionFactor=1.
is2aabb.aabbEnlargeFactor=1.
if runInGui:
plot.plot()
try:
from yade import qt
renderer=qt.Renderer()
except:
pass
---------------------------------------------------------------------------------------------------
I type O.run(), and after some time, the simulation break and I get this information:
----------------------------------------------------------------------------------------------------
Yade [102]: FATAL /home/dj/yade/trunk/pkg/dem/ConcretePM.cpp:410 go: Verification `!std::isnan(sigmaN)' failed!
FATAL /home/dj/yade/trunk/pkg/dem/ConcretePM.cpp:410 go: in interaction #4176+#4343
terminate called without an active exception
----------------------------------------------------------------------------------------------------
Can someone help me? Thanks very much!
