Question #703169 on Yade changed: https://answers.launchpad.net/yade/+question/703169
Xu gave more information on the question: Hi Karol, Sorry, the format of the code changes when i copy and paste it. I have already corrected it shown as below. There are two scripts, first is the consolidation and the second is the undrained test. After you run through the first script, a file of consolidated sample will be get, which can be directly loaded by the second script. > - minimal = short. It is difficult to find mistakes in long pieces of code > and it is usually much more annoying. Sorry for that, I have done my best to shorten the code. This is the corrected script ######Consolidation####### from yade import pack,plot,qt,export,ymport,utils import numpy as np IsoSigma = -1.e5 O.periodic=True O.trackEnergy=True idSand=O.materials.append(FrictMat(young=1e8,poisson=.8,frictionAngle=0.2,density=26500000,label='spheres')) sp1=pack.SpherePack() sp1.makeCloud(maxCorner=(0.1, 0.1, 0.1), psdSizes=[0.0017, 0.00191, 0.002285, 0.0026, 0.00292, 0.00325, 0.0035], psdCumm=[0.0, 0.1, 0.3, 0.5, 0.6, 0.9, 1], periodic=True,num=2000,seed=1) sp1.toSimulation(color=(0,0,1),material=idSand) v=utils.getSpheresVolume() print utils.getSpheresVolume() fout = file('Packing_volume.txt','w') fout.write(str(v)) fout.close() relRadList2 = [1,1] relPosList2 = [[0.4,0,0],[-0.4,0,0]] templates= [] templates.append(clumpTemplate(relRadii=relRadList2,relPositions=relPosList2)) O.bodies.replaceByClumps(templates,[1.],discretization=10) law=Law2_ScGeom_FrictPhys_CundallStrack(traceEnergy=True) O.engines=[ ForceResetter(), InsertionSortCollider([Bo1_Sphere_Aabb()]), InteractionLoop( [Ig2_Sphere_Sphere_ScGeom()], [Ip2_FrictMat_FrictMat_FrictPhys()], [law]), GlobalStiffnessTimeStepper(), NewtonIntegrator(damping=0.20), PeriTriaxController( goal=(IsoSigma,IsoSigma,IsoSigma), # Vector6 of prescribed final values stressMask=7, dynCell=True, maxStrainRate=(5.e+0,5.e+0,5.e+0), maxUnbalanced=0.001, relStressTol=1.e-3, doneHook='Finished()', label='p3d' ), ] def Finished(): O.save('post_iso.xml') print 'Finished' print (O.cell.size[0]*O.cell.size[1]*O.cell.size[2]-v)/v print getStress() print unbalancedForce() O.pause() O.run() ########Undrained triaxial test######### from yade import pack,plot,qt,export from math import fabs import numpy as np O.load('post_iso.xml') O.trackEnergy=True setContactFriction(0.5) filename1='Packing_volume.txt' a=np.loadtxt(filename1) Ep_ini=O.energy['elastPotential'] Edamp_ini=O.energy['nonviscDamp'] Ekin_ini=utils.kineticEnergy() work=0 law=Law2_ScGeom_FrictPhys_CundallStrack(traceEnergy=True) O.engines=[ ForceResetter(), InsertionSortCollider([Bo1_Sphere_Aabb()]), InteractionLoop( [Ig2_Sphere_Sphere_ScGeom()], [Ip2_FrictMat_FrictMat_FrictPhys()], [law]), NewtonIntegrator(damping=0.2), GlobalStiffnessTimeStepper(), PyRunner(command='inputwork()',iterPeriod=1), PyRunner(command='plotAddData()',iterPeriod=10), ] def inputwork(): global work dwork=(0.5*(O.cell.prevVelGrad+O.cell.velGrad)*utils.getStress()).trace()*O.dt*(O.cell.hSize).determinant() work=work+dwork def plotAddData(): global Ep_ini global Edamp_ini global Ekin_ini global work plot.addData( iter=O.iter,iter_=O.iter, sxx=utils.getStress()[0,0], syy=utils.getStress()[1,1], szz=utils.getStress()[2,2], exx=O.cell.getSmallStrain()[0,0], eyy=O.cell.getSmallStrain()[1,1], ezz=O.cell.getSmallStrain()[2,2], Z=avgNumInteractions(), Zm=avgNumInteractions(skipFree=True), voidratio=(O.cell.size[0]*O.cell.size[1]*O.cell.size[2]-a)/a, unbalanced=utils.unbalancedForce(), t=O.time, inwork=work, gWork=O.energy['gravWork'], # Ep=O.energy['elastPotential'], Ep=law.elasticEnergy()-Ep_ini, Edamp=O.energy['nonviscDamp']-Edamp_ini, Ediss=law.plasticDissipation(), # Ediss=O.energy['plastDissip'], Ekin=utils.kineticEnergy(), Etot=law.elasticEnergy()-Ep_ini+O.energy['nonviscDamp']-Edamp_ini+law.plasticDissipation()+utils.kineticEnergy()-Ekin_ini, # Etot=O.energy.total(),**O.energy ) plot.saveDataTxt('macroFile',vars=('t','exx','eyy','ezz','sxx','syy','szz','Z','Zm','voidratio')) plot.saveDataTxt('energyFile',vars=('t','unbalanced','gWork','Edamp','Ekin','Ep','Ediss','Etot','inwork')) O.cell.trsf=Matrix3.Identity O.cell.velGrad=Matrix3(-0.1,0,0,0,0.05,0,0,0,0.05) O.run() ########################################### My problom is the difference of "inwork" and "Etot" in the second script when using Clump particles. I really appreciate if you can confirm the issue or help me to solve the problem. Thanks a lot, Xu -- 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 : yade-users@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp